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 981 by gezelter, Mon Jun 5 18:24:45 2006 UTC

# Line 38 | Line 38
38   * University of Notre Dame has been advised of the possibility of
39   * such damages.
40   */
41 + #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
42 + #include "hydrodynamics/Shape.hpp"
43 + #include "hydrodynamics/Sphere.hpp"
44 + #include "hydrodynamics/Ellipsoid.hpp"
45 + #include "applications/hydrodynamics/CompositeShape.hpp"
46  
42 – #include "applications/hydrodynamics/HydrodynamicsModel.hpp"
43 – #include "math/LU.hpp"
44 – #include "math/DynamicRectMatrix.hpp"
45 – #include "math/SquareMatrix3.hpp"
46 – #include "utils/OOPSEConstant.hpp"
47   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 < */
48 >  
49 >  bool HydrodynamicsModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
50 >    return false;
51 >  }
52 >  
53 >  void HydrodynamicsModel::writeHydroProps(std::ostream& os) {
54  
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 <    if (!createBeads(beads_)) {
77 <        std::cout << "can not create beads" << std::endl;
78 <        return false;
79 <    }
80 <
81 <    //calcResistanceTensor();
82 <    calcDiffusionTensor();
55 >    Vector3d center;
56 >    Mat6x6d Xi, D;
57      
58 < /*    
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;
58 >    os << sd_->getType() << "\t";
59      
60 <    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 <    }
60 >    //center of resistance
61  
62 <    std::cout << "B=\n"
117 <                  << B << std::endl;
118 <    //invert B Matrix
119 <    invertMatrix(B, C);
62 >    center = cr_->getCOR();
63  
64 <    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 <    }
64 >    os << center[0] <<  "\t" << center[1] <<  "\t" << center[2] <<  "\t";
65      
66 <    //calculate Xi matrix at arbitrary origin O
67 <    Mat3x3d Xitt;
134 <    Mat3x3d Xirr;
135 <    Mat3x3d Xitr;
66 >    //resistance tensor at center of resistance
67 >    //translation
68  
69 <    //calculate the total volume
69 >    Xi = cr_->getXi();
70  
71 <    double volume = 0.0;
72 <    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
73 <        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);
71 >    os << Xi(0, 0) <<  "\t" << Xi(0, 1) <<  "\t" << Xi(0, 2) <<  "\t"
72 >       << Xi(1, 0) <<  "\t" << Xi(1, 1) <<  "\t" << Xi(1, 2) <<  "\t"
73 >       << Xi(2, 0) <<  "\t" << Xi(2, 1) <<  "\t" << Xi(2, 2) <<  "\t";
74      
75 <    Mat3x3d XittInv(0.0);
76 <    XittInv = Xitt.inverse();
75 >    //rotation-translation
76 >    os << Xi(0, 3) <<  "\t" << Xi(0, 4) <<  "\t" << Xi(0, 5) <<  "\t"
77 >       << Xi(1, 3) <<  "\t" << Xi(1, 4) <<  "\t" << Xi(1, 5) <<  "\t"
78 >       << Xi(2, 3) <<  "\t" << Xi(2, 4) <<  "\t" << Xi(2, 5) <<  "\t";
79      
80 <    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
81 <    //( a11 a12 0)
82 <    //( a21 a22 0)
83 <    //( 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();
80 >    //translation-rotation
81 >    os << Xi(3, 0) <<  "\t" << Xi(3, 1) <<  "\t" << Xi(3, 2) <<  "\t"
82 >       << Xi(4, 0) <<  "\t" << Xi(4, 1) <<  "\t" << Xi(4, 2) <<  "\t"
83 >       << Xi(5, 0) <<  "\t" << Xi(5, 1) <<  "\t" << Xi(5, 2) <<  "\t";
84      
85 <    Dorr = kt * tmpInv*1.0E16;
85 >    //rotation
86 >    os << Xi(3, 3) <<  "\t" << Xi(3, 4) <<  "\t" << Xi(3, 5) <<  "\t"
87 >       << Xi(4, 3) <<  "\t" << Xi(4, 4) <<  "\t" << Xi(4, 5) <<  "\t"
88 >       << Xi(5, 3) <<  "\t" << Xi(5, 4) <<  "\t" << Xi(5, 5) <<  "\t";
89 >    
90 >    
91 >    //diffusion tensor at center of resistance
92 >    //translation
93  
94 <    //Do.getSubMatrix(0, 0 , Dott);
198 <    //Do.getSubMatrix(3, 0, Dotr);
199 <    //Do.getSubMatrix(3, 3, Dorr);
94 >    D = cr_->getD();
95  
96 <    //calculate center of diffusion
97 <    tmp(0, 0) = Dorr(1, 1) + Dorr(2, 2);
98 <    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();
96 >    os << D(0, 0) <<  "\t" << D(0, 1) <<  "\t" << D(0, 2) <<  "\t"
97 >       << D(1, 0) <<  "\t" << D(1, 1) <<  "\t" << D(1, 2) <<  "\t"
98 >       << D(2, 0) <<  "\t" << D(2, 1) <<  "\t" << D(2, 2) <<  "\t";
99      
100 <    Vector3d rod = tmpInv * tmpVec;
101 <
102 <    //calculate Diffusion Tensor at center of diffusion
103 <    Mat3x3d Uod;
223 <    Uod.setupSkewMat(rod);
100 >    //rotation-translation
101 >    os << D(0, 3) <<  "\t" << D(0, 4) <<  "\t" << D(0, 5) <<  "\t"
102 >       << D(1, 3) <<  "\t" << D(1, 4) <<  "\t" << D(1, 5) <<  "\t"
103 >       << D(2, 3) <<  "\t" << D(2, 4) <<  "\t" << D(2, 5) <<  "\t";
104      
105 <    Mat3x3d Ddtt; //translational diffusion tensor at diffusion center
106 <    Mat3x3d Ddtr; //rotational diffusion tensor at diffusion center
107 <    Mat3x3d Ddrr; //translation-rotation couplingl diffusion tensor at diffusion tensor
105 >    //translation-rotation
106 >    os << D(3, 0) <<  "\t" << D(3, 1) <<  "\t" << D(3, 2) <<  "\t"
107 >       << D(4, 0) <<  "\t" << D(4, 1) <<  "\t" << D(4, 2) <<  "\t"
108 >       << D(5, 0) <<  "\t" << D(5, 1) <<  "\t" << D(5, 2) <<  "\t";
109      
110 <    Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
111 <    Ddrr = Dorr;
112 <    Ddtr = Dotr + Dorr * Uod;
110 >    //rotation
111 >    os << D(3, 3) <<  "\t" << D(3, 4) <<  "\t" << D(3, 5) <<  "\t"
112 >       << D(4, 3) <<  "\t" << D(4, 4) <<  "\t" << D(4, 5) <<  "\t"
113 >       << D(5, 3) <<  "\t" << D(5, 4) <<  "\t" << D(5, 5) <<  "\t";
114      
115 <    props_.diffCenter = rod;
234 <    props_.transDiff = Ddtt;
235 <    props_.transRotDiff = Ddtr;
236 <    props_.rotDiff = Ddrr;
237 < */
238 <    return true;    
239 < }
240 <
241 < void HydrodynamicsModel::calcResistanceTensor() {
115 >    //---------------------------------------------------------------------
116      
117 <    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 <    }
117 >    //center of diffusion
118  
119 <  
274 <    //invert B Matrix
275 <    invertMatrix(B, C);
119 >    center = cd_->getCOR();
120  
121 <    //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 <    }
121 >    os << center[0] <<  "\t" << center[1] <<  "\t" << center[2] <<  "\t";
122      
123 <    //calculate Xi matrix at arbitrary origin O
124 <    Mat3x3d Xiott;
287 <    Mat3x3d Xiorr;
288 <    Mat3x3d Xiotr;
123 >    //resistance tensor at center of diffusion
124 >    //translation
125  
126 <    //calculate the total volume
126 >    Xi = cd_->getXi();
127  
128 <    double volume = 0.0;
129 <    for (std::vector<BeadParam>::iterator iter = beads_.begin(); iter != beads_.end(); ++iter) {
130 <        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);
128 >    os << Xi(0, 0) <<  "\t" << Xi(0, 1) <<  "\t" << Xi(0, 2) <<  "\t"
129 >       << Xi(1, 0) <<  "\t" << Xi(1, 1) <<  "\t" << Xi(1, 2) <<  "\t"
130 >       << Xi(2, 0) <<  "\t" << Xi(2, 1) <<  "\t" << Xi(2, 2) <<  "\t";
131      
132 <    Mat3x3d Xirtt;
133 <    Mat3x3d Xirrr;
134 <    Mat3x3d Xirtr;
135 <
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 <
132 >    //rotation-translation
133 >    os << Xi(0, 3) <<  "\t" << Xi(0, 4) <<  "\t" << Xi(0, 5) <<  "\t"
134 >       << Xi(1, 3) <<  "\t" << Xi(1, 4) <<  "\t" << Xi(1, 5) <<  "\t"
135 >       << Xi(2, 3) <<  "\t" << Xi(2, 4) <<  "\t" << Xi(2, 5) <<  "\t";
136      
137 <    Mat3x3d XirttInv(0.0);
138 <    XirttInv = Xirtt.inverse();
137 >    //translation-rotation
138 >    os << Xi(3, 0) <<  "\t" << Xi(3, 1) <<  "\t" << Xi(3, 2) <<  "\t"
139 >       << Xi(4, 0) <<  "\t" << Xi(4, 1) <<  "\t" << Xi(4, 2) <<  "\t"
140 >       << Xi(5, 0) <<  "\t" << Xi(5, 1) <<  "\t" << Xi(5, 2) <<  "\t";
141      
142 <    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
143 <    //( a11 a12 0)
144 <    //( a21 a22 0)
145 <    //( 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();
142 >    //rotation
143 >    os << Xi(3, 3) <<  "\t" << Xi(3, 4) <<  "\t" << Xi(3, 5) <<  "\t"
144 >       << Xi(4, 3) <<  "\t" << Xi(4, 4) <<  "\t" << Xi(4, 5) <<  "\t"
145 >       << Xi(5, 3) <<  "\t" << Xi(5, 4) <<  "\t" << Xi(5, 5) <<  "\t";
146      
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;
147      
148 <    for (std::size_t i = 0; i < nbeads; ++i) {
149 <        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 <    }
148 >    //diffusion tensor at center of diffusion
149 >    //translation
150  
151 <    //invert B Matrix
439 <    invertMatrix(B, C);
151 >    D = cd_->getD();
152  
153 <    //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
154 <    std::vector<Mat3x3d> U;
155 <    for (int i = 0; i < nbeads; ++i) {
444 <        Mat3x3d currU;
445 <        currU.setupSkewMat(beads_[i].pos);
446 <        U.push_back(currU);
447 <    }
153 >    os << D(0, 0) <<  "\t" << D(0, 1) <<  "\t" << D(0, 2) <<  "\t"
154 >       << D(1, 0) <<  "\t" << D(1, 1) <<  "\t" << D(1, 2) <<  "\t"
155 >       << D(2, 0) <<  "\t" << D(2, 1) <<  "\t" << D(2, 2) <<  "\t";
156      
157 <    //calculate Xi matrix at arbitrary origin O
158 <    Mat3x3d Xitt;
159 <    Mat3x3d Xirr;
160 <    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);
157 >    //rotation-translation
158 >    os << D(0, 3) <<  "\t" << D(0, 4) <<  "\t" << D(0, 5) <<  "\t"
159 >       << D(1, 3) <<  "\t" << D(1, 4) <<  "\t" << D(1, 5) <<  "\t"
160 >       << D(2, 3) <<  "\t" << D(2, 4) <<  "\t" << D(2, 5) <<  "\t";
161      
162 <    Mat3x3d XittInv(0.0);
163 <    XittInv = Xitt.inverse();
162 >    //translation-rotation
163 >    os << D(3, 0) <<  "\t" << D(3, 1) <<  "\t" << D(3, 2) <<  "\t"
164 >       << D(4, 0) <<  "\t" << D(4, 1) <<  "\t" << D(4, 2) <<  "\t"
165 >       << D(5, 0) <<  "\t" << D(5, 1) <<  "\t" << D(5, 2) <<  "\t";
166      
167 <    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
168 <    //( a11 a12 0)
169 <    //( a21 a22 0)
170 <    //( 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;
167 >    //rotation
168 >    os << D(3, 3) <<  "\t" << D(3, 4) <<  "\t" << D(3, 5) <<  "\t"
169 >       << D(4, 3) <<  "\t" << D(4, 4) <<  "\t" << D(4, 5) <<  "\t"
170 >       << D(5, 3) <<  "\t" << D(5, 4) <<  "\t" << D(5, 5) <<  "\n";    
171      
172 <    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
173 <    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 <      
172 >  }
173 >  
174   }
587 –
588 – void HydrodynamicsModel::writeBeads(std::ostream& os) {
589 –    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 –
596 – }
597 –
598 – void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
599 –
600 –    os << sd_->getType() << "\t";
601 –    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\t";
602 –
603 –    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 –
611 –    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 –
615 –    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 –    
623 –    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 –
627 –    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 –    
631 – }
632 –
633 – }

Diff Legend

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