ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/applications/hydrodynamics/HydrodynamicsModel.cpp
Revision: 895
Committed: Mon Mar 13 22:42:40 2006 UTC (19 years, 7 months ago) by tim
Original Path: trunk/src/applications/hydrodynamics/HydrodynamicsModel.cpp
File size: 22146 byte(s)
Log Message: 
LangevinDynamics in progress

File Contents

# User Rev Content
1 tim 891 /*
2     * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
6     * redistribute this software in source and binary code form, provided
7     * that the following conditions are met:
8     *
9     * 1. Acknowledgement of the program authors must be made in any
10     * publication of scientific results based in part on use of the
11     * program. An acceptable form of acknowledgement is citation of
12     * the article in which the program was described (Matthew
13     * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     * Parallel Simulation Engine for Molecular Dynamics,"
16     * J. Comput. Chem. 26, pp. 252-271 (2005))
17     *
18     * 2. Redistributions of source code must retain the above copyright
19     * notice, this list of conditions and the following disclaimer.
20     *
21     * 3. Redistributions in binary form must reproduce the above copyright
22     * notice, this list of conditions and the following disclaimer in the
23     * documentation and/or other materials provided with the
24     * distribution.
25     *
26     * This software is provided "AS IS," without a warranty of any
27     * kind. All express or implied conditions, representations and
28     * warranties, including any implied warranty of merchantability,
29     * fitness for a particular purpose or non-infringement, are hereby
30     * excluded. The University of Notre Dame and its licensors shall not
31     * be liable for any damages suffered by licensee as a result of
32     * using, modifying or distributing the software or its
33     * derivatives. In no event will the University of Notre Dame or its
34     * licensors be liable for any lost revenue, profit or data, or for
35     * direct, indirect, special, consequential, incidental or punitive
36     * damages, however caused and regardless of the theory of liability,
37     * arising out of the use of or inability to use software, even if the
38     * University of Notre Dame has been advised of the possibility of
39     * such damages.
40     */
41    
42     #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
43     #include "math/LU.hpp"
44     #include "math/DynamicRectMatrix.hpp"
45     #include "math/SquareMatrix3.hpp"
46 tim 892 #include "utils/OOPSEConstant.hpp"
47 tim 891 namespace oopse {
48     /**
49     * Reference:
50     * Beatriz Carrasco and Jose Gracia de la Torre, Hydrodynamic Properties of Rigid Particles:
51     * Comparison of Different Modeling and Computational Procedures.
52     * Biophysical Journal, 75(6), 3044, 1999
53     */
54 tim 892
55     HydrodynamicsModel::HydrodynamicsModel(StuntDouble* sd, const DynamicProperty& extraParams) : sd_(sd){
56     DynamicProperty::const_iterator iter;
57    
58     iter = extraParams.find("Viscosity");
59     if (iter != extraParams.end()) {
60     boost::any param = iter->second;
61     viscosity_ = boost::any_cast<double>(param);
62     }else {
63     std::cout << "HydrodynamicsModel Error\n" ;
64     }
65    
66     iter = extraParams.find("Temperature");
67     if (iter != extraParams.end()) {
68     boost::any param = iter->second;
69     temperature_ = boost::any_cast<double>(param);
70     }else {
71     std::cout << "HydrodynamicsModel Error\n" ;
72     }
73     }
74    
75     bool HydrodynamicsModel::calcHydrodyanmicsProps() {
76 tim 891 if (!createBeads(beads_)) {
77     std::cout << "can not create beads" << std::endl;
78     return false;
79     }
80 tim 895
81     //calcResistanceTensor();
82     calcDiffusionTensor();
83 tim 891
84 tim 895 /*
85 tim 891 int nbeads = beads_.size();
86     DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
87     DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
88     Mat3x3d I;
89 tim 892 I(0, 0) = 1.0;
90     I(1, 1) = 1.0;
91     I(2, 2) = 1.0;
92    
93 tim 891 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 tim 893 tmpMat = outProduct(Rij, Rij) / rij2;
103 tim 892 double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
104 tim 891 Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
105     }else {
106 tim 892 double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
107 tim 891 Tij(0, 0) = constant;
108     Tij(1, 1) = constant;
109     Tij(2, 2) = constant;
110     }
111     B.setSubMatrix(i*3, j*3, Tij);
112 tim 893 std::cout << Tij << std::endl;
113 tim 891 }
114     }
115    
116 tim 893 std::cout << "B=\n"
117     << B << std::endl;
118 tim 891 //invert B Matrix
119     invertMatrix(B, C);
120 tim 893
121     std::cout << "C=\n"
122     << C << std::endl;
123    
124 tim 891 //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 tim 893
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 tim 895 volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
142 tim 893 }
143 tim 891
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 tim 893 Xitr += U[i] * Cij;
151 tim 895 //Xirr += -U[i] * Cij * U[j];
152     Xirr += -U[i] * Cij * U[j] + (0.166*6 * viscosity_ * volume) * I;
153 tim 891 }
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 tim 892 Xi.setSubMatrix(3, 3, Xirr);
163     //invertMatrix(Xi, Do);
164     double kt = OOPSEConstant::kB * temperature_ * 1.66E-2;
165     //Do *= kt;
166 tim 891
167 tim 892
168 tim 891 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 tim 893 const static Mat3x3d zeroMat(0.0);
173 tim 892
174     Mat3x3d XittInv(0.0);
175 tim 893 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 tim 892
184     Mat3x3d tmp;
185     Mat3x3d tmpInv;
186     tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
187 tim 893 tmpInv = tmp.inverse();
188 tim 892
189     Dott = kt * tmpInv;
190 tim 893 Dotr = -kt*XirrInv * Xitr * tmpInv* 1.0E8;
191 tim 892
192 tim 893 tmp = Xirr - Xitr * XittInv * Xitr.transpose();
193     tmpInv = tmp.inverse();
194 tim 892
195 tim 893 Dorr = kt * tmpInv*1.0E16;
196 tim 892
197     //Do.getSubMatrix(0, 0 , Dott);
198     //Do.getSubMatrix(3, 0, Dotr);
199     //Do.getSubMatrix(3, 3, Dorr);
200    
201 tim 891 //calculate center of diffusion
202 tim 892 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 tim 891
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 tim 893 tmpInv = tmp.inverse();
218 tim 892
219     Vector3d rod = tmpInv * tmpVec;
220    
221 tim 891 //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 tim 895 */
238 tim 891 return true;
239     }
240    
241 tim 895 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    
405     }
406    
407     void HydrodynamicsModel::calcDiffusionTensor() {
408     int nbeads = beads_.size();
409     DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
410     DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
411     Mat3x3d I;
412     I(0, 0) = 1.0;
413     I(1, 1) = 1.0;
414     I(2, 2) = 1.0;
415    
416     for (std::size_t i = 0; i < nbeads; ++i) {
417     for (std::size_t j = 0; j < nbeads; ++j) {
418     Mat3x3d Tij;
419     if (i != j ) {
420     Vector3d Rij = beads_[i].pos - beads_[j].pos;
421     double rij = Rij.length();
422     double rij2 = rij * rij;
423     double sumSigma2OverRij2 = ((beads_[i].radius*beads_[i].radius) + (beads_[i].radius*beads_[i].radius)) / rij2;
424     Mat3x3d tmpMat;
425     tmpMat = outProduct(Rij, Rij) / rij2;
426     double constant = 8.0 * NumericConstant::PI * viscosity_ * rij;
427     Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
428     }else {
429     double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity_ * beads_[i].radius);
430     Tij(0, 0) = constant;
431     Tij(1, 1) = constant;
432     Tij(2, 2) = constant;
433     }
434     B.setSubMatrix(i*3, j*3, Tij);
435     }
436     }
437    
438     //invert B Matrix
439     invertMatrix(B, C);
440    
441     //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
442     std::vector<Mat3x3d> U;
443     for (int i = 0; i < nbeads; ++i) {
444     Mat3x3d currU;
445     currU.setupSkewMat(beads_[i].pos);
446     U.push_back(currU);
447     }
448    
449     //calculate Xi matrix at arbitrary origin O
450     Mat3x3d Xitt;
451     Mat3x3d Xirr;
452     Mat3x3d Xitr;
453    
454     //calculate the total volume
455    
456     double volume = 0.0;
457     for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
458     volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
459     }
460    
461     for (std::size_t i = 0; i < nbeads; ++i) {
462     for (std::size_t j = 0; j < nbeads; ++j) {
463     Mat3x3d Cij;
464     C.getSubMatrix(i*3, j*3, Cij);
465    
466     Xitt += Cij;
467     Xitr += U[i] * Cij;
468     //Xirr += -U[i] * Cij * U[j];
469     Xirr += -U[i] * Cij * U[j] + (6 * viscosity_ * volume) * I;
470     }
471     }
472    
473     //invert Xi to get Diffusion Tensor at arbitrary origin O
474     RectMatrix<double, 6, 6> Xi;
475     RectMatrix<double, 6, 6> Do;
476     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;
482    
483     //1 poise = 0.1 N.S/m^2 = 1.661E-3 amu/ (Angstrom*fs)
484     double kt = OOPSEConstant::kB * temperature_ * 1.66E-3;
485    
486     Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
487     Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
488     Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
489    
490     const static Mat3x3d zeroMat(0.0);
491    
492     Mat3x3d XittInv(0.0);
493     XittInv = Xitt.inverse();
494    
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)
499     Mat3x3d XirrInv;
500     XirrInv = Xirr.inverse();
501    
502     Mat3x3d tmp;
503     Mat3x3d tmpInv;
504     tmp = Xitt - Xitr.transpose() * XirrInv * Xitr;
505     tmpInv = tmp.inverse();
506    
507     //Dott = kt * tmpInv; //unit in A^2/fs
508     Dott = tmpInv;
509     //Dotr = -kt*XirrInv * Xitr * tmpInv*1E8;
510     //Dotr = -kt*XirrInv * Xitr * tmpInv;
511     Dotr = -XirrInv* Xitr * tmpInv;
512    
513     tmp = Xirr - Xitr * XittInv * Xitr.transpose();
514     tmpInv = tmp.inverse();
515    
516     //Dorr = kt * tmpInv*1E16;
517     //Dorr = kt * tmpInv;
518     Dorr = tmpInv;
519     //calculate center of diffusion
520     tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
521     tmp(0, 1) = - Dorr(0, 1);
522     tmp(0, 2) = -Dorr(0, 2);
523     tmp(1, 0) = -Dorr(0, 1);
524     tmp(1, 1) = Dorr(0, 0) + Dorr(2, 2);
525     tmp(1, 2) = -Dorr(1, 2);
526     tmp(2, 0) = -Dorr(0, 2);
527     tmp(2, 1) = -Dorr(1, 2);
528     tmp(2, 2) = Dorr(1, 1) + Dorr(0, 0);
529    
530     Vector3d tmpVec;
531     tmpVec[0] = Dotr(1, 2) - Dotr(2, 1);
532     tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
533     tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
534    
535     tmpInv = tmp.inverse();
536    
537     Vector3d rod = tmpInv * tmpVec;
538    
539     //calculate Diffusion Tensor at center of diffusion
540     Mat3x3d Uod;
541     Uod.setupSkewMat(rod);
542    
543     Mat3x3d Ddtt; //translational diffusion tensor at diffusion center
544     Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center
545     Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor
546    
547     Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
548     Ddrr = Dorr;
549     Ddtr = Dotr + Dorr * Uod;
550    
551     props_.diffCenter = rod;
552     props_.Ddtt = Ddtt;
553     props_.Ddtr = Ddtr;
554     props_.Ddrr = Ddrr;
555    
556     SquareMatrix<double, 6> Dd;
557     Dd.setSubMatrix(0, 0, Ddtt);
558     Dd.setSubMatrix(0, 3, Ddtr.transpose());
559     Dd.setSubMatrix(3, 0, Ddtr);
560     Dd.setSubMatrix(3, 3, Ddrr);
561     SquareMatrix<double, 6> Xid;
562     invertMatrix(Dd, Xid);
563    
564     Ddtt *= kt;
565     Ddtr *= kt;
566     Ddrr *= kt;
567     Xid /= 1.66E-3;
568    
569     Xid.getSubMatrix(0, 0, props_.Xidtt);
570     Xid.getSubMatrix(0, 3, props_.Xidrt);
571     Xid.getSubMatrix(3, 0, props_.Xidtr);
572     Xid.getSubMatrix(3, 3, props_.Xidrr);
573    
574     /*
575     std::cout << "center of diffusion :" << std::endl;
576     std::cout << rod << std::endl;
577     std::cout << "diffusion tensor at center of diffusion" << std::endl;
578     std::cout << "translation:" << std::endl;
579     std::cout << Ddtt << std::endl;
580     std::cout << "translation-rotation:" << std::endl;
581     std::cout << Ddtr << std::endl;
582     std::cout << "rotation:" << std::endl;
583     std::cout << Ddrr << std::endl;
584     */
585    
586     }
587    
588 tim 891 void HydrodynamicsModel::writeBeads(std::ostream& os) {
589 tim 892 std::vector<BeadParam>::iterator iter;
590     os << beads_.size() << std::endl;
591     os << "Generated by Hydro" << std::endl;
592     for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
593     os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
594     }
595 tim 891
596     }
597    
598     void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
599    
600 tim 895 os << sd_->getType() << "\t";
601     os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t";
602 tim 892
603 tim 895 os << props_.Ddtt(0, 0) << "\t" << props_.Ddtt(0, 1) << "\t" << props_.Ddtt(0, 2) << "\t"
604     << props_.Ddtt(1, 0) << "\t" << props_.Ddtt(1, 1) << "\t" << props_.Ddtt(1, 2) << "\t"
605     << props_.Ddtt(2, 0) << "\t" << props_.Ddtt(2, 1) << "\t" << props_.Ddtt(2, 2) << "\t";
606    
607     os << props_.Ddtr(0, 0) << "\t" << props_.Ddtr(0, 1) << "\t" << props_.Ddtr(0, 2) << "\t"
608     << props_.Ddtr(1, 0) << "\t" << props_.Ddtr(1, 1) << "\t" << props_.Ddtr(1, 2) << "\t"
609     << props_.Ddtr(2, 0) << "\t" << props_.Ddtr(2, 1) << "\t" << props_.Ddtr(2, 2) << "\t";
610 tim 892
611 tim 895 os << props_.Ddrr(0, 0) << "\t" << props_.Ddrr(0, 1) << "\t" << props_.Ddrr(0, 2) << "\t"
612     << props_.Ddrr(1, 0) << "\t" << props_.Ddrr(1, 1) << "\t" << props_.Ddrr(1, 2) << "\t"
613     << props_.Ddrr(2, 0) << "\t" << props_.Ddrr(2, 1) << "\t" << props_.Ddrr(2, 2) <<"\t";
614 tim 892
615 tim 895 os << props_.Xidtt(0, 0) << "\t" << props_.Xidtt(0, 1) << "\t" << props_.Xidtt(0, 2) << "\t"
616     << props_.Xidtt(1, 0) << "\t" << props_.Xidtt(1, 1) << "\t" << props_.Xidtt(1, 2) << "\t"
617     << props_.Xidtt(2, 0) << "\t" << props_.Xidtt(2, 1) << "\t" << props_.Xidtt(2, 2) << "\t";
618    
619     os << props_.Xidrt(0, 0) << "\t" << props_.Xidrt(0, 1) << "\t" << props_.Xidrt(0, 2) << "\t"
620     << props_.Xidrt(1, 0) << "\t" << props_.Xidrt(1, 1) << "\t" << props_.Xidrt(1, 2) << "\t"
621     << props_.Xidrt(2, 0) << "\t" << props_.Xidrt(2, 1) << "\t" << props_.Xidrt(2, 2) << "\t";
622 tim 892
623 tim 895 os << props_.Xidtr(0, 0) << "\t" << props_.Xidtr(0, 1) << "\t" << props_.Xidtr(0, 2) << "\t"
624     << props_.Xidtr(1, 0) << "\t" << props_.Xidtr(1, 1) << "\t" << props_.Xidtr(1, 2) << "\t"
625     << props_.Xidtr(2, 0) << "\t" << props_.Xidtr(2, 1) << "\t" << props_.Xidtr(2, 2) << "\t";
626 tim 892
627 tim 895 os << props_.Xidrr(0, 0) << "\t" << props_.Xidrr(0, 1) << "\t" << props_.Xidrr(0, 2) << "\t"
628     << props_.Xidrr(1, 0) << "\t" << props_.Xidrr(1, 1) << "\t" << props_.Xidrr(1, 2) << "\t"
629     << props_.Xidrr(2, 0) << "\t" << props_.Xidrr(2, 1) << "\t" << props_.Xidrr(2, 2) << std::endl;
630 tim 892
631 tim 891 }
632    
633     }

Properties

Name Value
svn:executable *