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

Comparing trunk/src/applications/hydrodynamics/ApproximationModel.cpp (file contents):
Revision 954 by tim, Wed May 10 01:44:48 2006 UTC vs.
Revision 1208 by xsun, Wed Jan 16 20:19:28 2008 UTC

# Line 44 | Line 44
44   #include "math/DynamicRectMatrix.hpp"
45   #include "math/SquareMatrix3.hpp"
46   #include "utils/OOPSEConstant.hpp"
47 < #include "applications/hydrodynamics/Spheric.hpp"
48 < #include "applications/hydrodynamics/Ellipsoid.hpp"
47 > #include "hydrodynamics/Sphere.hpp"
48 > #include "hydrodynamics/Ellipsoid.hpp"
49   #include "applications/hydrodynamics/CompositeShape.hpp"
50   #include "math/LU.hpp"
51   #include "utils/simError.h"
# Line 57 | Line 57 | namespace oopse {
57   * Biophysical Journal, 75(6), 3044, 1999
58   */
59  
60 < ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){
61 <
62 < }
63 <
64 < bool ApproximationModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) {
65 <    return internalCalcHydroProps(static_cast<Shape*>(spheric), viscosity, temperature);
66 < }
67 <
68 < bool ApproximationModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
69 <    return internalCalcHydroProps(static_cast<Shape*>(ellipsoid), viscosity, temperature);
70 < }
71 < bool ApproximationModel::calcHydroProps(CompositeShape* compositeShape, double viscosity, double temperature) {
72 <    return internalCalcHydroProps(static_cast<Shape*>(compositeShape), viscosity, temperature);
73 < }
74 <
75 < void ApproximationModel::init() {
60 >  ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){    
61 >  }
62 >  
63 >  void ApproximationModel::init() {
64      if (!createBeads(beads_)) {
65        sprintf(painCave.errMsg, "ApproximationModel::init() : Can not create beads\n");
66        painCave.isFatal = 1;
67        simError();        
68      }
69 <
70 < }
71 <
72 < bool ApproximationModel::internalCalcHydroProps(Shape* shape, double viscosity, double temperature) {
73 <
69 >    
70 >  }
71 >  
72 >  bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
73 >    
74      bool ret = true;
75 <    HydroProps cr;
76 <    HydroProps cd;
75 >    HydroProp* cr = new HydroProp();
76 >    HydroProp* cd = new HydroProp();
77      calcHydroPropsAtCR(beads_, viscosity, temperature, cr);
78 <    //calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
78 >    calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
79      setCR(cr);
80      setCD(cd);
93 –    
81      return true;    
82 < }
83 <
84 < bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
85 <
82 >  }
83 >  
84 >  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cr) {
85 >    
86      int nbeads = beads.size();
87 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
88 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
87 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
88 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
89      Mat3x3d I;
90      I(0, 0) = 1.0;
91      I(1, 1) = 1.0;
92      I(2, 2) = 1.0;
93      
94      for (std::size_t i = 0; i < nbeads; ++i) {
95 <        for (std::size_t j = 0; j < nbeads; ++j) {
96 <            Mat3x3d Tij;
95 >      for (std::size_t j = 0; j < nbeads; ++j) {
96 >        Mat3x3d Tij;
97              if (i != j ) {
98 <                Vector3d Rij = beads[i].pos - beads[j].pos;
99 <                double rij = Rij.length();
100 <                double rij2 = rij * rij;
101 <                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
102 <                Mat3x3d tmpMat;
103 <                tmpMat = outProduct(Rij, Rij) / rij2;
104 <                double constant = 8.0 * NumericConstant::PI * viscosity * rij;
105 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
98 >              Vector3d Rij = beads[i].pos - beads[j].pos;
99 >              RealType rij = Rij.length();
100 >              RealType rij2 = rij * rij;
101 >              RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
102 >              Mat3x3d tmpMat;
103 >              tmpMat = outProduct(Rij, Rij) / rij2;
104 >              RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
105 >              RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
106 >              RealType tmp2 = 1.0 - sumSigma2OverRij2;
107 >              Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
108              }else {
109 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
110 <                Tij(0, 0) = constant;
111 <                Tij(1, 1) = constant;
112 <                Tij(2, 2) = constant;
109 >              RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
110 >              Tij(0, 0) = constant;
111 >              Tij(1, 1) = constant;
112 >              Tij(2, 2) = constant;
113              }
114              B.setSubMatrix(i*3, j*3, Tij);
115 <        }
115 >      }
116      }
117 <
117 >    
118      //invert B Matrix
119      invertMatrix(B, C);
120      
121      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
122      std::vector<Mat3x3d> U;
123      for (int i = 0; i < nbeads; ++i) {
124 <        Mat3x3d currU;
125 <        currU.setupSkewMat(beads[i].pos);
126 <        U.push_back(currU);
124 >      Mat3x3d currU;
125 >      currU.setupSkewMat(beads[i].pos);
126 >      U.push_back(currU);
127      }
128      
129      //calculate Xi matrix at arbitrary origin O
130      Mat3x3d Xiott;
131      Mat3x3d Xiorr;
132      Mat3x3d Xiotr;
133 <
133 >    
134      //calculate the total volume
135 <
136 <    double volume = 0.0;
135 >    
136 >    RealType volume = 0.0;
137      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
138 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
138 >      volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
139      }
140 <        
140 >    
141      for (std::size_t i = 0; i < nbeads; ++i) {
142 <        for (std::size_t j = 0; j < nbeads; ++j) {
143 <            Mat3x3d Cij;
144 <            C.getSubMatrix(i*3, j*3, Cij);
145 <            
146 <            Xiott += Cij;
147 <            Xiotr += U[i] * Cij;
148 <            //Xiorr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;    
149 <            Xiorr += -U[i] * Cij * U[j];
150 <        }
142 >      for (std::size_t j = 0; j < nbeads; ++j) {
143 >        Mat3x3d Cij;
144 >        C.getSubMatrix(i*3, j*3, Cij);
145 >        
146 >        Xiott += Cij;
147 >        Xiotr += U[i] * Cij;
148 >        // uncorrected here.  Volume correction is added after we assemble Xiorr
149 >        Xiorr += -U[i] * Cij * U[j];
150 >      }
151      }
152  
153 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
153 >    // add the volume correction
154 >    Xiorr += (6.0 * viscosity * volume) * I;    
155 >    
156 >    const RealType convertConstant = 1.439326479e4; //converts Poise angstroms
157 >                                                    // to kcal fs mol^-1 Angstrom^-1
158 >
159      Xiott *= convertConstant;
160      Xiotr *= convertConstant;
161      Xiorr *= convertConstant;
162      
169 –
170 –    
163      Mat3x3d tmp;
164      Mat3x3d tmpInv;
165      Vector3d tmpVec;
# Line 197 | Line 189 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
189      Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose();
190      
191  
192 <    SquareMatrix<double,6> Xir6x6;
193 <    SquareMatrix<double,6> Dr6x6;
192 >    SquareMatrix<RealType,6> Xir6x6;
193 >    SquareMatrix<RealType,6> Dr6x6;
194  
195      Xir6x6.setSubMatrix(0, 0, Xirtt);
196      Xir6x6.setSubMatrix(0, 3, Xirtr.transpose());
# Line 214 | Line 206 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
206      Dr6x6.getSubMatrix(0, 3, Drrt);
207      Dr6x6.getSubMatrix(3, 0, Drtr);
208      Dr6x6.getSubMatrix(3, 3, Drrr);
209 <    double kt = OOPSEConstant::kB * temperature ;
209 >    RealType kt = OOPSEConstant::kb * temperature ; // in kcal mol^-1
210      Drtt *= kt;
211      Drrt *= kt;
212      Drtr *= kt;
213      Drrr *= kt;
214 <    Xirtt *= OOPSEConstant::kb * temperature;
215 <    Xirtr *= OOPSEConstant::kb * temperature;
216 <    Xirrr *= OOPSEConstant::kb * temperature;
214 >    //Xirtt *= OOPSEConstant::kb * temperature;
215 >    //Xirtr *= OOPSEConstant::kb * temperature;
216 >    //Xirrr *= OOPSEConstant::kb * temperature;
217      
218 +    Mat6x6d Xi, D;
219  
220 <    cr.center = ror;
221 <    cr.Xi.setSubMatrix(0, 0, Xirtt);
222 <    cr.Xi.setSubMatrix(0, 3, Xirtr);
223 <    cr.Xi.setSubMatrix(3, 0, Xirtr);
224 <    cr.Xi.setSubMatrix(3, 3, Xirrr);
225 <    cr.D.setSubMatrix(0, 0, Drtt);
226 <    cr.D.setSubMatrix(0, 3, Drrt);
227 <    cr.D.setSubMatrix(3, 0, Drtr);
228 <    cr.D.setSubMatrix(3, 3, Drrr);    
220 >    cr->setCOR(ror);
221 >
222 >    Xi.setSubMatrix(0, 0, Xirtt);
223 >    Xi.setSubMatrix(0, 3, Xirtr);
224 >    Xi.setSubMatrix(3, 0, Xirtr);
225 >    Xi.setSubMatrix(3, 3, Xirrr);
226 >
227 >    cr->setXi(Xi);
228 >
229 >    D.setSubMatrix(0, 0, Drtt);
230 >    D.setSubMatrix(0, 3, Drrt);
231 >    D.setSubMatrix(3, 0, Drtr);
232 >    D.setSubMatrix(3, 3, Drrr);    
233 >
234 >    cr->setD(D);
235      
236      std::cout << "-----------------------------------------\n";
237      std::cout << "center of resistance :" << std::endl;
# Line 257 | Line 256 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
256  
257      return true;
258   }
259 <
260 < bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
261 <
259 >  
260 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProp* cd) {
261 >    
262      int nbeads = beads.size();
263 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
264 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
263 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
264 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
265      Mat3x3d I;
266      I(0, 0) = 1.0;
267      I(1, 1) = 1.0;
268      I(2, 2) = 1.0;
269      
270      for (std::size_t i = 0; i < nbeads; ++i) {
271 <        for (std::size_t j = 0; j < nbeads; ++j) {
272 <            Mat3x3d Tij;
273 <            if (i != j ) {
274 <                Vector3d Rij = beads[i].pos - beads[j].pos;
275 <                double rij = Rij.length();
276 <                double rij2 = rij * rij;
277 <                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
278 <                Mat3x3d tmpMat;
279 <                tmpMat = outProduct(Rij, Rij) / rij2;
280 <                double constant = 8.0 * NumericConstant::PI * viscosity * rij;
281 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
282 <            }else {
283 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
284 <                Tij(0, 0) = constant;
285 <                Tij(1, 1) = constant;
286 <                Tij(2, 2) = constant;
287 <            }
288 <            B.setSubMatrix(i*3, j*3, Tij);
271 >      for (std::size_t j = 0; j < nbeads; ++j) {
272 >        Mat3x3d Tij;
273 >        if (i != j ) {
274 >          Vector3d Rij = beads[i].pos - beads[j].pos;
275 >          RealType rij = Rij.length();
276 >          RealType rij2 = rij * rij;
277 >          RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
278 >          Mat3x3d tmpMat;
279 >          tmpMat = outProduct(Rij, Rij) / rij2;
280 >          RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
281 >          RealType tmp1 = 1.0 + sumSigma2OverRij2/3.0;
282 >          RealType tmp2 = 1.0 - sumSigma2OverRij2;
283 >          Tij = (tmp1 * I + tmp2 * tmpMat ) / constant;
284 >        }else {
285 >          RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
286 >          Tij(0, 0) = constant;
287 >          Tij(1, 1) = constant;
288 >          Tij(2, 2) = constant;
289          }
290 +        B.setSubMatrix(i*3, j*3, Tij);
291 +      }
292      }
293 <
293 >    
294      //invert B Matrix
295      invertMatrix(B, C);
296 <
296 >    
297      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
298      std::vector<Mat3x3d> U;
299      for (int i = 0; i < nbeads; ++i) {
300 <        Mat3x3d currU;
301 <        currU.setupSkewMat(beads[i].pos);
302 <        U.push_back(currU);
300 >      Mat3x3d currU;
301 >      currU.setupSkewMat(beads[i].pos);
302 >      U.push_back(currU);
303      }
304      
305      //calculate Xi matrix at arbitrary origin O
# Line 308 | Line 309 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
309  
310      //calculate the total volume
311  
312 <    double volume = 0.0;
312 >    RealType volume = 0.0;
313      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
314 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
314 >      volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
315      }
316 <        
316 >    
317      for (std::size_t i = 0; i < nbeads; ++i) {
318 <        for (std::size_t j = 0; j < nbeads; ++j) {
319 <            Mat3x3d Cij;
320 <            C.getSubMatrix(i*3, j*3, Cij);
318 >      for (std::size_t j = 0; j < nbeads; ++j) {
319 >        Mat3x3d Cij;
320 >        C.getSubMatrix(i*3, j*3, Cij);
321              
322 <            Xitt += Cij;
323 <            Xitr += U[i] * Cij;
324 <            //Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
325 <            Xirr += -U[i] * Cij * U[j];
326 <        }
322 >        Xitt += Cij;
323 >        Xitr += U[i] * Cij;
324 >        // uncorrected here.  Volume correction is added after we assemble Xiorr
325 >        Xirr += -U[i] * Cij * U[j];
326 >      }
327      }
328 <
329 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
328 >    // add the volume correction here:
329 >    Xirr += (6.0 * viscosity * volume) * I;    
330 >    
331 >    const RealType convertConstant = 1.439326479e4; //converts Poise angstroms
332 >                                                    // to kcal fs mol^-1 Angstrom^-1
333      Xitt *= convertConstant;
334      Xitr *= convertConstant;
335      Xirr *= convertConstant;
336 <
337 <    double kt = OOPSEConstant::kB * temperature;
338 <
336 >    
337 >    RealType kt = OOPSEConstant::kb * temperature; // in kcal mol^-1
338 >    
339      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
340      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
341      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
342 <
342 >    
343      const static Mat3x3d zeroMat(0.0);
344      
345      Mat3x3d XittInv(0.0);
# Line 389 | Line 393 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
393      Ddrr = Dorr;
394      Ddtr = Dotr + Dorr * Uod;
395  
396 <    SquareMatrix<double, 6> Dd;
396 >    SquareMatrix<RealType, 6> Dd;
397      Dd.setSubMatrix(0, 0, Ddtt);
398      Dd.setSubMatrix(0, 3, Ddtr.transpose());
399      Dd.setSubMatrix(3, 0, Ddtr);
400      Dd.setSubMatrix(3, 3, Ddrr);    
401 <    SquareMatrix<double, 6> Xid;
401 >    SquareMatrix<RealType, 6> Xid;
402      Ddtt *= kt;
403      Ddtr *=kt;
404      Ddrr *= kt;
# Line 406 | Line 410 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
410      //Xid /= OOPSEConstant::energyConvert;
411      Xid *= OOPSEConstant::kb * temperature;
412  
413 <    cr.center = rod;
410 <    cr.D.setSubMatrix(0, 0, Ddtt);
411 <    cr.D.setSubMatrix(0, 3, Ddtr);
412 <    cr.D.setSubMatrix(3, 0, Ddtr);
413 <    cr.D.setSubMatrix(3, 3, Ddrr);
414 <    cr.Xi = Xid;
413 >    Mat6x6d Xi, D;
414  
415 +    cd->setCOR(rod);
416 +
417 +    cd->setXi(Xid);
418 +
419 +    D.setSubMatrix(0, 0, Ddtt);
420 +    D.setSubMatrix(0, 3, Ddtr);
421 +    D.setSubMatrix(3, 0, Ddtr);
422 +    D.setSubMatrix(3, 3, Ddrr);
423 +
424 +    cd->setD(D);
425 +
426      std::cout << "viscosity = " << viscosity << std::endl;
427      std::cout << "temperature = " << temperature << std::endl;
428      std::cout << "center of diffusion :" << std::endl;
429      std::cout << rod << std::endl;
430      std::cout << "diffusion tensor at center of diffusion " << std::endl;
431 <    std::cout << "translation(A^2/fs) :" << std::endl;
431 >    std::cout << "translation(A^2 / fs) :" << std::endl;
432      std::cout << Ddtt << std::endl;
433 <    std::cout << "translation-rotation(A^3/fs):" << std::endl;
433 >    std::cout << "translation-rotation(A / fs):" << std::endl;
434      std::cout << Ddtr << std::endl;
435 <    std::cout << "rotation(A^4/fs):" << std::endl;
435 >    std::cout << "rotation(fs^-1):" << std::endl;
436      std::cout << Ddrr << std::endl;
437  
438      std::cout << "resistance tensor at center of diffusion " << std::endl;
# Line 438 | Line 448 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
448      Xid.getSubMatrix(3, 3, Xidrr);
449  
450      std::cout << Xidtt << std::endl;
451 <    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-3):" << std::endl;
451 >    std::cout << "rotation-translation (kcal*fs*mol^-1*Ang^-1):" << std::endl;
452      std::cout << Xidrt << std::endl;
453 <    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-3):" << std::endl;
453 >    std::cout << "translation-rotation(kcal*fs*mol^-1*Ang^-1):" << std::endl;
454      std::cout << Xidtr << std::endl;
455 <    std::cout << "rotation(kcal*fs*mol^-1*Ang^-4):" << std::endl;
455 >    std::cout << "rotation(kcal*fs*mol^-1):" << std::endl;
456      std::cout << Xidrr << std::endl;
457  
458      return true;
459 <      
460 < }
459 >    
460 >  }
461  
462 <
453 < void ApproximationModel::writeBeads(std::ostream& os) {
462 >  void ApproximationModel::writeBeads(std::ostream& os) {
463      std::vector<BeadParam>::iterator iter;
464      os << beads_.size() << std::endl;
465      os << "Generated by Hydro" << std::endl;
466      for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
467 <        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
467 >      os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
468      }
469 <
469 >    
470 >  }    
471   }
462 –
463 –
464 –
465 – }

Diff Legend

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