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 2611 by tim, Mon Mar 13 22:42:40 2006 UTC vs.
Revision 2768 by gezelter, Wed May 24 18:31:12 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 < */
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();
48 >  
49 >  bool HydrodynamicsModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
50 >    return false;
51 >  }
52 >  
53 >  void HydrodynamicsModel::writeHydroProps(std::ostream& os) {
54      
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;
55      
56 <    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 <    }
56 >    os << sd_->getType() << "\t";
57      
58 <    //calculate Xi matrix at arbitrary origin O
59 <    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);
58 >    //center of resistance
59 >    os << cr_.center[0] <<  "\t" << cr_.center[1] <<  "\t" << cr_.center[2] <<  "\t";
60      
61 <    Mat3x3d XittInv(0.0);
62 <    XittInv = Xitt.inverse();
61 >    //resistance tensor at center of resistance
62 >    //translation
63 >    os << cr_.Xi(0, 0) <<  "\t" << cr_.Xi(0, 1) <<  "\t" << cr_.Xi(0, 2) <<  "\t"
64 >       << cr_.Xi(1, 0) <<  "\t" << cr_.Xi(1, 1) <<  "\t" << cr_.Xi(1, 2) <<  "\t"
65 >       << cr_.Xi(2, 0) <<  "\t" << cr_.Xi(2, 1) <<  "\t" << cr_.Xi(2, 2) <<  "\t";
66      
67 <    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
68 <    //( a11 a12 0)
69 <    //( a21 a22 0)
70 <    //( 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();
67 >    //rotation-translation
68 >    os << cr_.Xi(0, 3) <<  "\t" << cr_.Xi(0, 4) <<  "\t" << cr_.Xi(0, 5) <<  "\t"
69 >       << cr_.Xi(1, 3) <<  "\t" << cr_.Xi(1, 4) <<  "\t" << cr_.Xi(1, 5) <<  "\t"
70 >       << cr_.Xi(2, 3) <<  "\t" << cr_.Xi(2, 4) <<  "\t" << cr_.Xi(2, 5) <<  "\t";
71      
72 <    Dorr = kt * tmpInv*1.0E16;
73 <
74 <    //Do.getSubMatrix(0, 0 , Dott);
75 <    //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();
72 >    //translation-rotation
73 >    os << cr_.Xi(3, 0) <<  "\t" << cr_.Xi(3, 1) <<  "\t" << cr_.Xi(3, 2) <<  "\t"
74 >       << cr_.Xi(4, 0) <<  "\t" << cr_.Xi(4, 1) <<  "\t" << cr_.Xi(4, 2) <<  "\t"
75 >       << cr_.Xi(5, 0) <<  "\t" << cr_.Xi(5, 1) <<  "\t" << cr_.Xi(5, 2) <<  "\t";
76      
77 <    Vector3d rod = tmpInv * tmpVec;
78 <
79 <    //calculate Diffusion Tensor at center of diffusion
80 <    Mat3x3d Uod;
223 <    Uod.setupSkewMat(rod);
77 >    //rotation
78 >    os << cr_.Xi(3, 3) <<  "\t" << cr_.Xi(3, 4) <<  "\t" << cr_.Xi(3, 5) <<  "\t"
79 >       << cr_.Xi(4, 3) <<  "\t" << cr_.Xi(4, 4) <<  "\t" << cr_.Xi(4, 5) <<  "\t"
80 >       << cr_.Xi(5, 3) <<  "\t" << cr_.Xi(5, 4) <<  "\t" << cr_.Xi(5, 5) <<  "\t";
81      
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
82      
83 <    Ddtt = Dott - Uod * Dorr * Uod + Dotr.transpose() * Uod - Uod * Dotr;
84 <    Ddrr = Dorr;
85 <    Ddtr = Dotr + Dorr * Uod;
83 >    //diffusion tensor at center of resistance
84 >    //translation
85 >    os << cr_.D(0, 0) <<  "\t" << cr_.D(0, 1) <<  "\t" << cr_.D(0, 2) <<  "\t"
86 >       << cr_.D(1, 0) <<  "\t" << cr_.D(1, 1) <<  "\t" << cr_.D(1, 2) <<  "\t"
87 >       << cr_.D(2, 0) <<  "\t" << cr_.D(2, 1) <<  "\t" << cr_.D(2, 2) <<  "\t";
88      
89 <    props_.diffCenter = rod;
90 <    props_.transDiff = Ddtt;
91 <    props_.transRotDiff = Ddtr;
92 <    props_.rotDiff = Ddrr;
237 < */
238 <    return true;    
239 < }
240 <
241 < void HydrodynamicsModel::calcResistanceTensor() {
89 >    //rotation-translation
90 >    os << cr_.D(0, 3) <<  "\t" << cr_.D(0, 4) <<  "\t" << cr_.D(0, 5) <<  "\t"
91 >       << cr_.D(1, 3) <<  "\t" << cr_.D(1, 4) <<  "\t" << cr_.D(1, 5) <<  "\t"
92 >       << cr_.D(2, 3) <<  "\t" << cr_.D(2, 4) <<  "\t" << cr_.D(2, 5) <<  "\t";
93      
94 <    int nbeads = beads_.size();
95 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
96 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
97 <    Mat3x3d I;
247 <    I(0, 0) = 1.0;
248 <    I(1, 1) = 1.0;
249 <    I(2, 2) = 1.0;
94 >    //translation-rotation
95 >    os << cr_.D(3, 0) <<  "\t" << cr_.D(3, 1) <<  "\t" << cr_.D(3, 2) <<  "\t"
96 >       << cr_.D(4, 0) <<  "\t" << cr_.D(4, 1) <<  "\t" << cr_.D(4, 2) <<  "\t"
97 >       << cr_.D(5, 0) <<  "\t" << cr_.D(5, 1) <<  "\t" << cr_.D(5, 2) <<  "\t";
98      
99 <    for (std::size_t i = 0; i < nbeads; ++i) {
100 <        for (std::size_t j = 0; j < nbeads; ++j) {
101 <            Mat3x3d Tij;
102 <            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 <    }
99 >    //rotation
100 >    os << cr_.D(3, 3) <<  "\t" << cr_.D(3, 4) <<  "\t" << cr_.D(3, 5) <<  "\t"
101 >       << cr_.D(4, 3) <<  "\t" << cr_.D(4, 4) <<  "\t" << cr_.D(4, 5) <<  "\t"
102 >       << cr_.D(5, 3) <<  "\t" << cr_.D(5, 4) <<  "\t" << cr_.D(5, 5) <<  "\t";
103      
104 <    //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);
104 >    //---------------------------------------------------------------------
105      
106 <    Mat3x3d Xirtt;
107 <    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 <
106 >    //center of diffusion
107 >    os << cd_.center[0] <<  "\t" << cd_.center[1] <<  "\t" << cd_.center[2] <<  "\t";
108      
109 <    Mat3x3d XirttInv(0.0);
110 <    XirttInv = Xirtt.inverse();
109 >    //resistance tensor at center of diffusion
110 >    //translation
111 >    os << cd_.Xi(0, 0) <<  "\t" << cd_.Xi(0, 1) <<  "\t" << cd_.Xi(0, 2) <<  "\t"
112 >       << cd_.Xi(1, 0) <<  "\t" << cd_.Xi(1, 1) <<  "\t" << cd_.Xi(1, 2) <<  "\t"
113 >       << cd_.Xi(2, 0) <<  "\t" << cd_.Xi(2, 1) <<  "\t" << cd_.Xi(2, 2) <<  "\t";
114      
115 <    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
116 <    //( a11 a12 0)
117 <    //( a21 a22 0)
118 <    //( 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();
115 >    //rotation-translation
116 >    os << cd_.Xi(0, 3) <<  "\t" << cd_.Xi(0, 4) <<  "\t" << cd_.Xi(0, 5) <<  "\t"
117 >       << cd_.Xi(1, 3) <<  "\t" << cd_.Xi(1, 4) <<  "\t" << cd_.Xi(1, 5) <<  "\t"
118 >       << cd_.Xi(2, 3) <<  "\t" << cd_.Xi(2, 4) <<  "\t" << cd_.Xi(2, 5) <<  "\t";
119      
120 <    Drrr = kt * tmpInv*1.0E16;
121 <
122 <    std::cout << "-----------------------------------------\n";
123 <    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;
120 >    //translation-rotation
121 >    os << cd_.Xi(3, 0) <<  "\t" << cd_.Xi(3, 1) <<  "\t" << cd_.Xi(3, 2) <<  "\t"
122 >       << cd_.Xi(4, 0) <<  "\t" << cd_.Xi(4, 1) <<  "\t" << cd_.Xi(4, 2) <<  "\t"
123 >       << cd_.Xi(5, 0) <<  "\t" << cd_.Xi(5, 1) <<  "\t" << cd_.Xi(5, 2) <<  "\t";
124      
125 <    for (std::size_t i = 0; i < nbeads; ++i) {
126 <        for (std::size_t j = 0; j < nbeads; ++j) {
127 <            Mat3x3d Tij;
128 <            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 <    }
125 >    //rotation
126 >    os << cd_.Xi(3, 3) <<  "\t" << cd_.Xi(3, 4) <<  "\t" << cd_.Xi(3, 5) <<  "\t"
127 >       << cd_.Xi(4, 3) <<  "\t" << cd_.Xi(4, 4) <<  "\t" << cd_.Xi(4, 5) <<  "\t"
128 >       << cd_.Xi(5, 3) <<  "\t" << cd_.Xi(5, 4) <<  "\t" << cd_.Xi(5, 5) <<  "\t";
129      
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);
130      
131 <    Mat3x3d XittInv(0.0);
132 <    XittInv = Xitt.inverse();
131 >    //diffusion tensor at center of diffusion
132 >    //translation
133 >    os << cd_.D(0, 0) <<  "\t" << cd_.D(0, 1) <<  "\t" << cd_.D(0, 2) <<  "\t"
134 >       << cd_.D(1, 0) <<  "\t" << cd_.D(1, 1) <<  "\t" << cd_.D(1, 2) <<  "\t"
135 >       << cd_.D(2, 0) <<  "\t" << cd_.D(2, 1) <<  "\t" << cd_.D(2, 2) <<  "\t";
136      
137 <    //Xirr may not be inverted,if it one of the diagonal element is zero, for example
138 <    //( a11 a12 0)
139 <    //( a21 a22 0)
140 <    //( 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;
137 >    //rotation-translation
138 >    os << cd_.D(0, 3) <<  "\t" << cd_.D(0, 4) <<  "\t" << cd_.D(0, 5) <<  "\t"
139 >       << cd_.D(1, 3) <<  "\t" << cd_.D(1, 4) <<  "\t" << cd_.D(1, 5) <<  "\t"
140 >       << cd_.D(2, 3) <<  "\t" << cd_.D(2, 4) <<  "\t" << cd_.D(2, 5) <<  "\t";
141      
142 <    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
143 <    tmpInv = tmp.inverse();
142 >    //translation-rotation
143 >    os << cd_.D(3, 0) <<  "\t" << cd_.D(3, 1) <<  "\t" << cd_.D(3, 2) <<  "\t"
144 >       << cd_.D(4, 0) <<  "\t" << cd_.D(4, 1) <<  "\t" << cd_.D(4, 2) <<  "\t"
145 >       << cd_.D(5, 0) <<  "\t" << cd_.D(5, 1) <<  "\t" << cd_.D(5, 2) <<  "\t";
146      
147 <    //Dorr = kt * tmpInv*1E16;
148 <    //Dorr = kt * tmpInv;
149 <    Dorr = tmpInv;
150 <    //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();
147 >    //rotation
148 >    os << cd_.D(3, 3) <<  "\t" << cd_.D(3, 4) <<  "\t" << cd_.D(3, 5) <<  "\t"
149 >       << cd_.D(4, 3) <<  "\t" << cd_.D(4, 4) <<  "\t" << cd_.D(4, 5) <<  "\t"
150 >       << cd_.D(5, 3) <<  "\t" << cd_.D(5, 4) <<  "\t" << cd_.D(5, 5) <<  "\n";    
151      
152 <    Vector3d rod = tmpInv * tmpVec;
153 <
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 <      
152 >  }
153 >  
154   }
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