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

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

# Line 77 | Line 77 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
77          std::cout << "can not create beads" << std::endl;
78          return false;
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);
# Line 95 | Line 99 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
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(beads_[i].pos, beads_[j].pos) / rij2;
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 {
# Line 105 | Line 109 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
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) {
# Line 122 | Line 133 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
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) {
# Line 129 | Line 147 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
147              C.getSubMatrix(i*3, j*3, Cij);
148              
149              Xitt += Cij;
150 <            Xirr += U[i] * Cij;
151 <            Xitr += U[i] * Cij * U[j];            
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  
# Line 150 | Line 169 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
169      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
170      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
171  
172 <    Mat3x3d XirrInv(0.0);
154 <    Mat3x3d XirrCopy;
155 <    XirrCopy = Xirr;
172 >    const static Mat3x3d zeroMat(0.0);
173      
174      Mat3x3d XittInv(0.0);
175 <    Mat3x3d XittCopy;
176 <    XittCopy = Xitt;
177 <    invertMatrix(XittCopy, XittInv);
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  
166 –    const static Mat3x3d zeroMat(0.0);
167 –    if (!invertMatrix(tmp, tmpInv)) {
168 –        tmpInv = zeroMat;
169 –    }
170 –
189      Dott = kt * tmpInv;
190 <    Dotr = -kt*XirrInv * Xitr * tmpInv;
190 >    Dotr = -kt*XirrInv * Xitr * tmpInv* 1.0E8;
191  
192 <    tmp = Xirr - Xitr * XittInv * Xitr.transpose();
192 >    tmp = Xirr - Xitr * XittInv * Xitr.transpose();    
193 >    tmpInv = tmp.inverse();
194      
195 <    if(!invertMatrix(tmp, tmpInv)) {
177 <        tmpInv = zeroMat;
178 <    }
179 <    Dorr = kt * tmpInv;
195 >    Dorr = kt * tmpInv*1.0E16;
196  
197      //Do.getSubMatrix(0, 0 , Dott);
198      //Do.getSubMatrix(3, 0, Dotr);
# Line 198 | Line 214 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
214      tmpVec[1] = Dotr(2, 0) - Dotr(0, 2);
215      tmpVec[2] = Dotr(0, 1) - Dotr(1, 0);
216  
217 <    if(!invertMatrix(tmp, tmpInv)) {
202 <        tmpInv = zeroMat;
203 <    }
217 >    tmpInv = tmp.inverse();
218      
219      Vector3d rod = tmpInv * tmpVec;
220  
# Line 220 | Line 234 | bool HydrodynamicsModel::calcHydrodyanmicsProps() {
234      props_.transDiff = Ddtt;
235      props_.transRotDiff = Ddtr;
236      props_.rotDiff = Ddrr;
237 <
237 > */
238      return true;    
239   }
240  
241 < void HydrodynamicsModel::writeBeads(std::ostream& os) {
242 <    std::vector<BeadParam>::iterator iter;
243 <    os << beads_.size() << std::endl;
244 <    os << "Generated by Hydro" << std::endl;
245 <    for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
246 <        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
241 > 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 < }
273 >  
274 >    //invert B Matrix
275 >    invertMatrix(B, C);
276  
277 < void HydrodynamicsModel::writeDiffCenterAndDiffTensor(std::ostream& os) {
278 <    os << "//viscosity = " << viscosity_ << std::endl;
279 <    os << "//temperature = " << temperature_<< std::endl;
280 <    std::vector<BeadParam>::iterator iter;
281 <    os << sd_->getType() << "\n";
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 <    os << "//diffusion center" << std::endl;
244 <    os << props_.diffCenter << std::endl;
245 <
246 <    os << "//translational diffusion tensor" << std::endl;
247 <    os << props_.transDiff << std::endl;
290 >    //calculate the total volume
291  
292 <    os << "//translational diffusion tensor" << std::endl;
293 <    os << props_.transRotDiff << std::endl;
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 <    os << "//rotational diffusion tensor" << std::endl;
310 <    os << props_.rotDiff << std::endl;
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 <    /*
330 <    os << props_.diffCenter[0] << "\t" << props_.diffCenter[1] << "\t" << props_.diffCenter[2] << "\n"
329 >    Mat3x3d Xirtt;
330 >    Mat3x3d Xirrr;
331 >    Mat3x3d Xirtr;
332  
333 <    os << props_.transDiff(0, 0) << "\t" << props_.transDiff(0, 1) << "\t" << props_.transDiff(0, 2) << "\t"
334 <        << props_.transDiff(1, 0) << "\t" << props_.transDiff(1, 1) << "\t" << props_.transDiff(1, 2) << "\t"
335 <        << props_.transDiff(2, 0) << "\t" << props_.transDiff(2, 1) << "\t" << props_.transDiff(2, 2) << "\n";
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 <    os << props_.transRotDiff(0, 0) << "\t" << props_.transRotDiff(0, 1) << "\t" << props_.transRotDiff(0, 2) << "\t"
363 <        << props_.transRotDiff(1, 0) << "\t" << props_.transRotDiff(1, 1) << "\t" << props_.transRotDiff(1, 2) << "\t"
364 <        << props_.transRotDiff(2, 0) << "\t" << props_.transRotDiff(2, 1) << "\t" << props_.transRotDiff(2, 2) << "\t"
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 <    os << props_.rotDiff(0, 0) << "\t" << props_.rotDiff(0, 1) << "\t" << props_.rotDiff(0, 2) << "\t"
375 <        << props_.rotDiff(1, 0) << "\t" << props_.rotDiff(1, 1) << "\t" << props_.rotDiff(1, 2) << "\t"
376 <        << props_.rotDiff(2, 0) << "\t" << props_.rotDiff(2, 1) << "\t" << props_.rotDiff(2, 2) << ";"
377 <        << std::endl;
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 + 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