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

Comparing trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp (file contents):
Revision 895 by tim, Mon Mar 13 22:42:40 2006 UTC vs.
Revision 900 by tim, Wed Mar 15 20:58:47 2006 UTC

# Line 79 | Line 79 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
79      }
80  
81      //calcResistanceTensor();
82 <    calcDiffusionTensor();
83 <    
84 < /*    
85 <    int nbeads = beads_.size();
86 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
87 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
88 <    Mat3x3d I;
89 <    I(0, 0) = 1.0;
90 <    I(1, 1) = 1.0;
91 <    I(2, 2) = 1.0;
92 <    
93 <    for (std::size_t i = 0; i < nbeads; ++i) {
94 <        for (std::size_t j = 0; j < nbeads; ++j) {
95 <            Mat3x3d Tij;
96 <            if (i != j ) {
97 <                Vector3d Rij = beads_[i].pos - beads_[j].pos;
98 <                double rij = Rij.length();
99 <                double rij2 = rij * rij;
100 <                double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
101 <                Mat3x3d tmpMat;
102 <                tmpMat = outProduct(Rij, Rij) / rij2;
103 <                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
104 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
105 <            }else {
106 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
107 <                Tij(0, 0) = constant;
108 <                Tij(1, 1) = constant;
109 <                Tij(2, 2) = constant;
110 <            }
111 <            B.setSubMatrix(i*3, j*3, Tij);
112 <            std::cout << Tij << std::endl;
113 <        }
114 <    }
115 <
116 <    std::cout << "B=\n"
117 <                  << B << std::endl;
118 <    //invert B Matrix
119 <    invertMatrix(B, C);
120 <
121 <    std::cout << "C=\n"
122 <                  << C << std::endl;
123 <
124 <    //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
125 <    std::vector<Mat3x3d> U;
126 <    for (int i = 0; i < nbeads; ++i) {
127 <        Mat3x3d currU;
128 <        currU.setupSkewMat(beads_[i].pos);
129 <        U.push_back(currU);
130 <    }
131 <    
132 <    //calculate Xi matrix at arbitrary origin O
133 <    Mat3x3d Xitt;
134 <    Mat3x3d Xirr;
135 <    Mat3x3d Xitr;
136 <
137 <    //calculate the total volume
138 <
139 <    double volume = 0.0;
140 <    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
141 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
142 <    }
143 <        
144 <    for (std::size_t i = 0; i < nbeads; ++i) {
145 <        for (std::size_t j = 0; j < nbeads; ++j) {
146 <            Mat3x3d Cij;
147 <            C.getSubMatrix(i*3, j*3, Cij);
148 <            
149 <            Xitt += Cij;
150 <            Xitr += U[i] * Cij;
151 <            //Xirr += -U[i] * Cij * U[j];            
152 <            Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I;            
153 <        }
154 <    }
155 <
156 <    //invert Xi to get Diffusion Tensor at arbitrary origin O
157 <    RectMatrix<double, 6, 6> Xi;    
158 <    RectMatrix<double, 6, 6> Do;
159 <    Xi.setSubMatrix(0, 0, Xitt);
160 <    Xi.setSubMatrix(0, 3, Xitr.transpose());
161 <    Xi.setSubMatrix(3, 0, Xitr);
162 <    Xi.setSubMatrix(3, 3, Xirr);
163 <    //invertMatrix(Xi, Do);
164 <    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
165 <    //Do *= kt;    
166 <
167 <
168 <    Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
169 <    Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
170 <    Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
171 <
172 <    const static Mat3x3d zeroMat(0.0);
173 <    
174 <    Mat3x3d XittInv(0.0);
175 <    XittInv = Xitt.inverse();
176 <    
177 <    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
178 <    //( a11 a12 0)
179 <    //( a21 a22 0)
180 <    //( 0    0    0)
181 <    Mat3x3d XirrInv;
182 <    XirrInv = Xirr.inverse();
183 <
184 <    Mat3x3d tmp;
185 <    Mat3x3d tmpInv;
186 <    tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
187 <    tmpInv = tmp.inverse();
188 <
189 <    Dott = kt * tmpInv;
190 <    Dotr = -kt*XirrInv * Xitr * tmpInv* 1.0E8;
191 <
192 <    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
193 <    tmpInv = tmp.inverse();
194 <    
195 <    Dorr = kt * tmpInv*1.0E16;
196 <
197 <    //Do.getSubMatrix(0, 0 , Dott);
198 <    //Do.getSubMatrix(3, 0, Dotr);
199 <    //Do.getSubMatrix(3, 3, Dorr);
200 <
201 <    //calculate center of diffusion
202 <    tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
203 <    tmp(0, 1) = - Dorr(0, 1);
204 <    tmp(0, 2) = -Dorr(0, 2);
205 <    tmp(1, 0) = -Dorr(0, 1);
206 <    tmp(1, 1) = Dorr(0, 0)  + Dorr(2, 2);
207 <    tmp(1, 2) = -Dorr(1, 2);
208 <    tmp(2, 0) = -Dorr(0, 2);
209 <    tmp(2, 1) = -Dorr(1, 2);
210 <    tmp(2, 2) = Dorr(1, 1) + Dorr(0, 0);
211 <
212 <    Vector3d tmpVec;
213 <    tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
214 <    tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
215 <    tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
216 <
217 <    tmpInv = tmp.inverse();
218 <    
219 <    Vector3d rod = tmpInv * tmpVec;
220 <
221 <    //calculate Diffusion Tensor at center of diffusion
222 <    Mat3x3d Uod;
223 <    Uod.setupSkewMat(rod);
224 <    
225 <    Mat3x3d Ddtt; //translational diffusion tensor at diffusion center
226 <    Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center
227 <    Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor
228 <    
229 <    Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
230 <    Ddrr = Dorr;
231 <    Ddtr = Dotr + Dorr * Uod;
232 <    
233 <    props_.diffCenter = rod;
234 <    props_.transDiff = Ddtt;
235 <    props_.transRotDiff = Ddtr;
236 <    props_.rotDiff = Ddrr;
237 < */
82 >    calcDiffusionTensor();    
83      return true;    
84   }
85  
86   void HydrodynamicsModel::calcResistanceTensor() {
242 –    
243 –    int nbeads = beads_.size();
244 –    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
245 –    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
246 –    Mat3x3d I;
247 –    I(0, 0) = 1.0;
248 –    I(1, 1) = 1.0;
249 –    I(2, 2) = 1.0;
250 –    
251 –    for (std::size_t i = 0; i < nbeads; ++i) {
252 –        for (std::size_t j = 0; j < nbeads; ++j) {
253 –            Mat3x3d Tij;
254 –            if (i != j ) {
255 –                Vector3d Rij = beads_[i].pos - beads_[j].pos;
256 –                double rij = Rij.length();
257 –                double rij2 = rij * rij;
258 –                double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;                
259 –                Mat3x3d tmpMat;
260 –                tmpMat = outProduct(Rij, Rij) / rij2;
261 –                double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
262 –                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
263 –            }else {
264 –                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
265 –                Tij(0, 0) = constant;
266 –                Tij(1, 1) = constant;
267 –                Tij(2, 2) = constant;
268 –            }
269 –            B.setSubMatrix(i*3, j*3, Tij);
270 –        }
271 –    }
272 –
273 –  
274 –    //invert B Matrix
275 –    invertMatrix(B, C);
276 –
277 –    //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
278 –    std::vector<Mat3x3d> U;
279 –    for (int i = 0; i < nbeads; ++i) {
280 –        Mat3x3d currU;
281 –        currU.setupSkewMat(beads_[i].pos);
282 –        U.push_back(currU);
283 –    }
284 –    
285 –    //calculate Xi matrix at arbitrary origin O
286 –    Mat3x3d Xiott;
287 –    Mat3x3d Xiorr;
288 –    Mat3x3d Xiotr;
289 –
290 –    //calculate the total volume
291 –
292 –    double volume = 0.0;
293 –    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
294 –        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
295 –    }
296 –        
297 –    for (std::size_t i = 0; i < nbeads; ++i) {
298 –        for (std::size_t j = 0; j < nbeads; ++j) {
299 –            Mat3x3d Cij;
300 –            C.getSubMatrix(i*3, j*3, Cij);
301 –            
302 –            Xiott += Cij;
303 –            Xiotr += U[i] * Cij;
304 –            //Xiorr += -U[i] * Cij * U[j];            
305 –            Xiorr += -U[i] * Cij * U[j] + (6 * viscosity_ * volume) * I;            
306 –        }
307 –    }
308 –
309 –    Mat3x3d tmp;
310 –    Mat3x3d tmpInv;
311 –    Vector3d tmpVec;
312 –    tmp(0, 0) = Xiott(1, 1) + Xiott(2, 2);
313 –    tmp(0, 1) = - Xiott(0, 1);
314 –    tmp(0, 2) = -Xiott(0, 2);
315 –    tmp(1, 0) = -Xiott(0, 1);
316 –    tmp(1, 1) = Xiott(0, 0)  + Xiott(2, 2);
317 –    tmp(1, 2) = -Xiott(1, 2);
318 –    tmp(2, 0) = -Xiott(0, 2);
319 –    tmp(2, 1) = -Xiott(1, 2);
320 –    tmp(2, 2) = Xiott(1, 1) + Xiott(0, 0);
321 –    tmpVec[0] = Xiotr(2, 1) - Xiotr(1, 2);
322 –    tmpVec[1] = Xiotr(0, 2) - Xiotr(2, 0);
323 –    tmpVec[2] = Xiotr(1, 0) - Xiotr(0, 1);
324 –    tmpInv = tmp.inverse();    
325 –    Vector3d ror = tmpInv * tmpVec; //center of resistance
326 –    Mat3x3d Uor;
327 –    Uor.setupSkewMat(ror);
328 –    
329 –    Mat3x3d Xirtt;
330 –    Mat3x3d Xirrr;
331 –    Mat3x3d Xirtr;
332 –
333 –    Xirtt = Xiott;
334 –    Xirtr = (Xiotr - Uor * Xiott) * 1E-8;
335 –    Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose() * 1E-16;
336 – /*
337 –    SquareMatrix<double,6> Xir6x6;
338 –    SquareMatrix<double,6> Dr6x6;
339 –
340 –    Xir6x6.setSubMatrix(0, 0, Xirtt);
341 –    Xir6x6.setSubMatrix(0, 3, Xirtr.transpose());
342 –    Xir6x6.setSubMatrix(3, 0, Xirtr);
343 –    Xir6x6.setSubMatrix(3, 3, Xirrr);
344 –
345 –    invertMatrix(Xir6x6, Dr6x6);
346 –    Mat3x3d Drtt;
347 –    Mat3x3d Drtr;
348 –    Mat3x3d Drrr;
349 –    Dr6x6.getSubMatrix(0, 0, Drtt);
350 –    Dr6x6.getSubMatrix(3, 0, Drtr);
351 –    Dr6x6.getSubMatrix(3, 3, Drrr);
352 –    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
353 –    Drtt *= kt;
354 –    Drtr *= kt*1E8;
355 –    Drrr *= kt*1E16;
356 – */  
357 –
358 –    const static Mat3x3d zeroMat(0.0);
359 –
360 –
361 –    
362 –    Mat3x3d XirttInv(0.0);
363 –    XirttInv = Xirtt.inverse();
364 –    
365 –    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
366 –    //( a11 a12 0)
367 –    //( a21 a22 0)
368 –    //( 0    0    0)
369 –    Mat3x3d XirrrInv;
370 –    XirrrInv = Xirrr.inverse();
371 –    tmp = Xirtt - Xirtr.transpose() * XirrrInv * Xirtr;
372 –    tmpInv = tmp.inverse();
373 –
374 –    Mat3x3d Drtt;
375 –    Mat3x3d Drtr;
376 –    Mat3x3d Drrr;
377 –    double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
378 –    Drtt = kt * tmpInv;
379 –    Drtr = -kt*XirrrInv * Xirtr * tmpInv* 1.0E8;
380 –
381 –    tmp = Xirrr - Xirtr * XirttInv * Xirtr.transpose();    
382 –    tmpInv = tmp.inverse();
383 –    
384 –    Drrr = kt * tmpInv*1.0E16;
385 –
386 –    std::cout << "-----------------------------------------\n";
387 –    std::cout << "center of resistance :" << std::endl;
388 –    std::cout << ror << std::endl;
389 –    std::cout << "resistant tensor at center of resistance" << std::endl;
390 –    std::cout << "translation:" << std::endl;
391 –    std::cout << Xirtt << std::endl;
392 –    std::cout << "translation-rotation:" << std::endl;
393 –    std::cout << Xirtr << std::endl;
394 –    std::cout << "rotation:" << std::endl;
395 –    std::cout << Xirrr << std::endl;
396 –    std::cout << "diffusion tensor at center of resistance" << std::endl;
397 –    std::cout << "translation:" << std::endl;
398 –    std::cout << Drtt << std::endl;
399 –    std::cout << "translation-rotation:" << std::endl;
400 –    std::cout << Drtr << std::endl;
401 –    std::cout << "rotation:" << std::endl;
402 –    std::cout << Drrr << std::endl;    
403 –    std::cout << "-----------------------------------------\n";
404 –
87   }
88  
89   void HydrodynamicsModel::calcDiffusionTensor() {
# Line 465 | Line 147 | void HydrodynamicsModel::calcDiffusionTensor() {
147              
148              Xitt += Cij;
149              Xitr += U[i] * Cij;
468 –            //Xirr += -U[i] * Cij * U[j];            
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);
477 <    Xi.setSubMatrix(0, 3, Xitr.transpose());
478 <    Xi.setSubMatrix(3, 0, Xitr);
479 <    Xi.setSubMatrix(3, 3, Xirr);
480 <    //invertMatrix(Xi, Do);
481 <    //double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
154 >    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
155 >    Xitt *= convertConstant;
156 >    Xitr *= convertConstant;
157 >    Xirr *= convertConstant;
158  
159 <    //1 poise = 0.1 N.S/m^2 = 1.661E-3 amu/ (Angstrom*fs)
484 <    double kt = OOPSEConstant::kB * temperature_ * 1.66E-3;
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
# Line 492 | Line 167 | void HydrodynamicsModel::calcDiffusionTensor() {
167      Mat3x3d XittInv(0.0);
168      XittInv = Xitt.inverse();
169      
495 –    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
496 –    //( a11 a12 0)
497 –    //( a21 a22 0)
498 –    //( 0    0    0)
170      Mat3x3d XirrInv;
171      XirrInv = Xirr.inverse();
172  
# Line 504 | Line 175 | void HydrodynamicsModel::calcDiffusionTensor() {
175      tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
176      tmpInv = tmp.inverse();
177  
178 <    //Dott = kt * tmpInv; //unit in A^2/fs
179 <    Dott = tmpInv;
509 <    //Dotr = -kt*XirrInv * Xitr * tmpInv*1E8;
510 <    //Dotr = -kt*XirrInv * Xitr * tmpInv;
511 <    Dotr = -XirrInv* Xitr * tmpInv;
178 >    Dott = kt*tmpInv;
179 >    Dotr = -kt*XirrInv * Xitr * tmpInv;
180      
181      tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
182      tmpInv = tmp.inverse();
183      
184 <    //Dorr = kt * tmpInv*1E16;
185 <    //Dorr = kt * tmpInv;
518 <    Dorr = tmpInv;
184 >    Dorr = kt * tmpInv;
185 >
186      //calculate center of diffusion
187      tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
188      tmp(0, 1) = - Dorr(0, 1);
# Line 561 | Line 228 | void HydrodynamicsModel::calcDiffusionTensor() {
228      SquareMatrix<double, 6> Xid;
229      invertMatrix(Dd, Xid);
230  
231 <    Ddtt *= kt;
232 <    Ddtr *= kt;
566 <    Ddrr *= kt;
567 <    Xid /= 1.66E-3;
231 >    //Xidtt in units of kcal*fs*mol^-1*Ang^-2
232 >    Xid *= OOPSEConstant::kb*temperature_;
233  
234      Xid.getSubMatrix(0, 0, props_.Xidtt);
235      Xid.getSubMatrix(0, 3, props_.Xidrt);
236      Xid.getSubMatrix(3, 0, props_.Xidtr);
237      Xid.getSubMatrix(3, 3, props_.Xidrr);
238  
239 <    /*        
239 >
240      std::cout << "center of diffusion :" << std::endl;
241      std::cout << rod << std::endl;
242 <    std::cout << "diffusion tensor at center of diffusion" << std::endl;
243 <    std::cout << "translation:" << std::endl;
242 >    std::cout << "diffusion tensor at center of diffusion " << std::endl;
243 >    std::cout << "translation(A^2/fs) :" << std::endl;
244      std::cout << Ddtt << std::endl;
245 <    std::cout << "translation-rotation:" << std::endl;
245 >    std::cout << "translation-rotation(A^3/fs):" << std::endl;
246      std::cout << Ddtr << std::endl;
247 <    std::cout << "rotation:" << std::endl;
247 >    std::cout << "rotation(A^4/fs):" << std::endl;
248      std::cout << Ddrr << std::endl;
249 <    */
249 >
250 >    std::cout << "resistance tensor at center of diffusion " << std::endl;
251 >    std::cout << "translation(kcal*fs*mol^-1*Ang^-2) :" << std::endl;
252 >    std::cout << props_.Xidtt << std::endl;
253 >    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl;
254 >    std::cout << props_.Xidrt << std::endl;
255 >    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl;
256 >    std::cout << props_.Xidtr << std::endl;
257 >    std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl;
258 >    std::cout << props_.Xidrr << std::endl;
259 >    
260        
261   }
262  

Diff Legend

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