ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/applications/hydrodynamics/HydrodynamicsModel.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/applications/hydrodynamics/HydrodynamicsModel.cpp (file contents):
Revision 2597 by tim, Thu Feb 23 23:16:43 2006 UTC vs.
Revision 2632 by tim, Thu Mar 16 22:50:48 2006 UTC

# Line 77 | Line 77 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
77          std::cout << "can not create beads" << std::endl;
78          return false;
79      }
80 <    
80 >
81 >    //calcResistanceTensor();
82 >    calcDiffusionTensor();    
83 >    return true;    
84 > }
85 >
86 > void HydrodynamicsModel::calcResistanceTensor() {
87 > }
88 >
89 > void HydrodynamicsModel::calcDiffusionTensor() {
90      int nbeads = beads_.size();
91      DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
92      DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
# Line 95 | Line 104 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
104                  double rij2 = rij * rij;
105                  double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
106                  Mat3x3d tmpMat;
107 <                tmpMat = outProduct(beads_[i].pos, beads_[j].pos) / rij2;
107 >                tmpMat = outProduct(Rij, Rij) / rij2;
108                  double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
109                  Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
110              }else {
# Line 110 | Line 119 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
119  
120      //invert B Matrix
121      invertMatrix(B, C);
122 +
123      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
124      std::vector<Mat3x3d> U;
125      for (int i = 0; i < nbeads; ++i) {
# Line 122 | Line 132 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
132      Mat3x3d Xitt;
133      Mat3x3d Xirr;
134      Mat3x3d Xitr;
135 +
136 +    //calculate the total volume
137 +
138 +    double volume = 0.0;
139 +    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
140 +        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
141 +    }
142          
143      for (std::size_t i = 0; i < nbeads; ++i) {
144          for (std::size_t j = 0; j < nbeads; ++j) {
# Line 129 | Line 146 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
146              C.getSubMatrix(i*3, j*3, Cij);
147              
148              Xitt += Cij;
149 <            Xirr += U[i] * Cij;
150 <            Xitr += U[i] * Cij * U[j];            
149 >            Xitr += U[i] * Cij;
150 >            Xirr += -U[i] * Cij * U[j] + (6 * viscosity_ * volume) * I;            
151          }
152      }
153  
154 <    //invert Xi to get Diffusion Tensor at arbitrary origin O
155 <    RectMatrix<double, 6, 6> Xi;    
156 <    RectMatrix<double, 6, 6> Do;
157 <    Xi.setSubMatrix(0, 0, Xitt);
141 <    Xi.setSubMatrix(0, 3, Xitr.transpose());
142 <    Xi.setSubMatrix(3, 0, Xitr);
143 <    Xi.setSubMatrix(3, 3, Xirr);
144 <    //invertMatrix(Xi, Do);
145 <    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
146 <    //Do *= kt;    
154 >    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
155 >    Xitt *= convertConstant;
156 >    Xitr *= convertConstant;
157 >    Xirr *= convertConstant;
158  
159 +    double kt = OOPSEConstant::kB * temperature_;
160  
161      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
162      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
163      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
164  
165 <    Mat3x3d XirrInv(0.0);
154 <    Mat3x3d XirrCopy;
155 <    XirrCopy = Xirr;
165 >    const static Mat3x3d zeroMat(0.0);
166      
167      Mat3x3d XittInv(0.0);
168 <    Mat3x3d XittCopy;
169 <    XittCopy = Xitt;
170 <    invertMatrix(XittCopy, XittInv);
168 >    XittInv = Xitt.inverse();
169 >    
170 >    Mat3x3d XirrInv;
171 >    XirrInv = Xirr.inverse();
172  
173      Mat3x3d tmp;
174      Mat3x3d tmpInv;
175      tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
176 +    tmpInv = tmp.inverse();
177  
178 <    const static Mat3x3d zeroMat(0.0);
179 <    if (!invertMatrix(tmp, tmpInv)) {
168 <        tmpInv = zeroMat;
169 <    }
170 <
171 <    Dott = kt * tmpInv;
172 <    Dotr = -kt*XirrInv * Xitr * tmpInv;
173 <
174 <    tmp = Xirr - Xitr * XittInv * Xitr.transpose();
178 >    Dott = tmpInv;
179 >    Dotr = -XirrInv * Xitr * tmpInv;
180      
181 <    if(!invertMatrix(tmp, tmpInv)) {
182 <        tmpInv = zeroMat;
183 <    }
184 <    Dorr = kt * tmpInv;
181 >    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
182 >    tmpInv = tmp.inverse();
183 >    
184 >    Dorr = tmpInv;
185  
181 –    //Do.getSubMatrix(0, 0 , Dott);
182 –    //Do.getSubMatrix(3, 0, Dotr);
183 –    //Do.getSubMatrix(3, 3, Dorr);
184 –
186      //calculate center of diffusion
187      tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
188      tmp(0, 1) = - Dorr(0, 1);
# Line 198 | Line 199 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
199      tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
200      tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
201  
202 <    if(!invertMatrix(tmp, tmpInv)) {
202 <        tmpInv = zeroMat;
203 <    }
202 >    tmpInv = tmp.inverse();
203      
204      Vector3d rod = tmpInv * tmpVec;
205  
# Line 215 | Line 214 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
214      Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
215      Ddrr = Dorr;
216      Ddtr = Dotr + Dorr * Uod;
218 –    
219 –    props_.diffCenter = rod;
220 –    props_.transDiff = Ddtt;
221 –    props_.transRotDiff = Ddtr;
222 –    props_.rotDiff = Ddrr;
217  
218 <    return true;    
218 >    SquareMatrix<double, 6> Dd;
219 >    Dd.setSubMatrix(0, 0, Ddtt);
220 >    Dd.setSubMatrix(0, 3, Ddtr.transpose());
221 >    Dd.setSubMatrix(3, 0, Ddtr);
222 >    Dd.setSubMatrix(3, 3, Ddrr);    
223 >    SquareMatrix<double, 6> Xid;
224 >    Ddtt *= kt;
225 >    Ddtr *=kt;
226 >    Ddrr *= kt;
227 >    invertMatrix(Dd, Xid);
228 >
229 >
230 >
231 >    //Xidtt in units of kcal*fs*mol^-1*Ang^-2
232 >    //Xid /= OOPSEConstant::energyConvert;
233 >    Xid *= OOPSEConstant::kb * temperature_;
234 >    props_.diffCenter = rod;
235 >    props_.Ddtt = Ddtt;
236 >    props_.Ddtr = Ddtr;
237 >    props_.Ddrr = Ddrr;
238 >    Xid.getSubMatrix(0, 0, props_.Xidtt);
239 >    Xid.getSubMatrix(0, 3, props_.Xidrt);
240 >    Xid.getSubMatrix(3, 0, props_.Xidtr);
241 >    Xid.getSubMatrix(3, 3, props_.Xidrr);
242 >
243 >
244 >    std::cout << "viscosity = " << viscosity_ << std::endl;
245 >    std::cout << "temperature = " << temperature_ << std::endl;
246 >    std::cout << "center of diffusion :" << std::endl;
247 >    std::cout << rod << std::endl;
248 >    std::cout << "diffusion tensor at center of diffusion " << std::endl;
249 >    std::cout << "translation(A^2/fs) :" << std::endl;
250 >    std::cout << Ddtt << std::endl;
251 >    std::cout << "translation-rotation(A^3/fs):" << std::endl;
252 >    std::cout << Ddtr << std::endl;
253 >    std::cout << "rotation(A^4/fs):" << std::endl;
254 >    std::cout << Ddrr << std::endl;
255 >
256 >    std::cout << "resistance tensor at center of diffusion " << std::endl;
257 >    std::cout << "translation(kcal*fs*mol^-1*Ang^-2) :" << std::endl;
258 >    std::cout << props_.Xidtt << std::endl;
259 >    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl;
260 >    std::cout << props_.Xidrt << std::endl;
261 >    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl;
262 >    std::cout << props_.Xidtr << std::endl;
263 >    std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl;
264 >    std::cout << props_.Xidrr << std::endl;
265 >    
266 >      
267   }
268  
269   void HydrodynamicsModel::writeBeads(std::ostream& os) {
# Line 235 | Line 277 | void HydrodynamicsModel::writeDiffCenterAndDiffTensor(
277   }
278  
279   void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
238 –    os << "//viscosity = " << viscosity_ << std::endl;
239 –    os << "//temperature = " << temperature_<< std::endl;
240 –    std::vector<BeadParam>::iterator iter;
241 –    os << sd_->getType() << "\n";
280  
281 <    os << "//diffusion center" << std::endl;
282 <    os << props_.diffCenter << std::endl;
281 >    os << sd_->getType() << "\t";
282 >    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t";
283  
284 <    os << "//translational diffusion tensor" << std::endl;
285 <    os << props_.transDiff << std::endl;
284 >    os << props_.Ddtt(0, 0) << "\t" << props_.Ddtt(0, 1) << "\t" << props_.Ddtt(0, 2) << "\t"
285 >        << props_.Ddtt(1, 0) << "\t" << props_.Ddtt(1, 1) << "\t" << props_.Ddtt(1, 2) << "\t"
286 >        << props_.Ddtt(2, 0) << "\t" << props_.Ddtt(2, 1) << "\t" << props_.Ddtt(2, 2) << "\t";
287 >    
288 >    os << props_.Ddtr(0, 0) << "\t" << props_.Ddtr(0, 1) << "\t" << props_.Ddtr(0, 2) << "\t"
289 >        << props_.Ddtr(1, 0) << "\t" << props_.Ddtr(1, 1) << "\t" << props_.Ddtr(1, 2) << "\t"
290 >        << props_.Ddtr(2, 0) << "\t" << props_.Ddtr(2, 1) << "\t" << props_.Ddtr(2, 2) << "\t";
291  
292 <    os << "//translational diffusion tensor" << std::endl;
293 <    os << props_.transRotDiff << std::endl;
292 >    os << props_.Ddrr(0, 0) << "\t" << props_.Ddrr(0, 1) << "\t" << props_.Ddrr(0, 2) << "\t"
293 >        << props_.Ddrr(1, 0) << "\t" << props_.Ddrr(1, 1) << "\t" << props_.Ddrr(1, 2) << "\t"
294 >        << props_.Ddrr(2, 0) << "\t" << props_.Ddrr(2, 1) << "\t" << props_.Ddrr(2, 2) <<"\t";        
295  
296 <    os << "//rotational diffusion tensor" << std::endl;
297 <    os << props_.rotDiff << std::endl;
298 <    
255 <    /*
256 <    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
296 >    os << props_.Xidtt(0, 0) << "\t" << props_.Xidtt(0, 1) << "\t" << props_.Xidtt(0, 2) << "\t"
297 >        << props_.Xidtt(1, 0) << "\t" << props_.Xidtt(1, 1) << "\t" << props_.Xidtt(1, 2) << "\t"
298 >        << props_.Xidtt(2, 0) << "\t" << props_.Xidtt(2, 1) << "\t" << props_.Xidtt(2, 2) << "\t";
299  
300 <    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
301 <        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
302 <        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
300 >    os << props_.Xidrt(0, 0) << "\t" << props_.Xidrt(0, 1) << "\t" << props_.Xidrt(0, 2) << "\t"
301 >        << props_.Xidrt(1, 0) << "\t" << props_.Xidrt(1, 1) << "\t" << props_.Xidrt(1, 2) << "\t"
302 >        << props_.Xidrt(2, 0) << "\t" << props_.Xidrt(2, 1) << "\t" << props_.Xidrt(2, 2) << "\t";
303      
304 <    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
305 <        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
306 <        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
304 >    os << props_.Xidtr(0, 0) << "\t" << props_.Xidtr(0, 1) << "\t" << props_.Xidtr(0, 2) << "\t"
305 >        << props_.Xidtr(1, 0) << "\t" << props_.Xidtr(1, 1) << "\t" << props_.Xidtr(1, 2) << "\t"
306 >        << props_.Xidtr(2, 0) << "\t" << props_.Xidtr(2, 1) << "\t" << props_.Xidtr(2, 2) << "\t";
307  
308 <    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
309 <        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
310 <        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
311 <        << std::endl;
270 <    */
308 >    os << props_.Xidrr(0, 0) << "\t" << props_.Xidrr(0, 1) << "\t" << props_.Xidrr(0, 2) << "\t"
309 >        << props_.Xidrr(1, 0) << "\t" << props_.Xidrr(1, 1) << "\t" << props_.Xidrr(1, 2) << "\t"
310 >        << props_.Xidrr(2, 0) << "\t" << props_.Xidrr(2, 1) << "\t" << props_.Xidrr(2, 2) << std::endl;
311 >    
312   }
313  
314   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines