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 910 by tim, Tue Mar 21 15:33:04 2006 UTC vs.
Revision 977 by tim, Thu May 25 16:27:27 2006 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"
52   namespace oopse {
53   /**
54   * Reference:
# Line 56 | 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 <    DynamicProperty::const_iterator iter;
63 <
63 <    iter = extraParams.find("Viscosity");
64 <    if (iter != extraParams.end()) {
65 <        boost::any param = iter->second;
66 <        viscosity = boost::any_cast<double>(param);
67 <    }else {
68 <        std::cout << "ApproximationModel Error\n" ;
69 <    }
70 <
71 <    iter = extraParams.find("Temperature");
72 <    if (iter != extraParams.end()) {
73 <        boost::any param = iter->second;
74 <        temperature = boost::any_cast<double>(param);
75 <    }else {
76 <        std::cout << "ApproximationModel Error\n" ;
77 <    }    
78 < */
79 < }
80 <
81 < bool ApproximationModel::calcHydroProps(Spheric* spheric, double viscosity, double temperature) {
82 <    return internalCalcHydroProps(static_cast<Shape*>(spheric), viscosity, temperature);
83 < }
84 <
85 < bool ApproximationModel::calcHydroProps(Ellipsoid* ellipsoid, double viscosity, double temperature) {
86 <    return internalCalcHydroProps(static_cast<Shape*>(ellipsoid), viscosity, temperature);
87 < }
88 < bool ApproximationModel::calcHydroProps(CompositeShape* compositeShape, double viscosity, double temperature) {
89 <    return internalCalcHydroProps(static_cast<Shape*>(compositeShape), viscosity, temperature);
90 < }
91 <
92 <
93 < bool ApproximationModel::internalCalcHydroProps(Shape* shape, double viscosity, double temperature) {
60 >  ApproximationModel::ApproximationModel(StuntDouble* sd, SimInfo* info): HydrodynamicsModel(sd, info){    
61 >  }
62 >  
63 >  void ApproximationModel::init() {
64      if (!createBeads(beads_)) {
65 <        std::cout << "can not create beads" << std::endl;
66 <        return false;
65 >      sprintf(painCave.errMsg, "ApproximationModel::init() : Can not create beads\n");
66 >      painCave.isFatal = 1;
67 >      simError();        
68      }
69 <
69 >    
70 >  }
71 >  
72 >  bool ApproximationModel::calcHydroProps(Shape* shape, RealType viscosity, RealType temperature) {
73 >    
74      bool ret = true;
75      HydroProps cr;
76      HydroProps cd;
77      calcHydroPropsAtCR(beads_, viscosity, temperature, cr);
78 <    calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
78 >    //calcHydroPropsAtCD(beads_, viscosity, temperature, cd);
79      setCR(cr);
80      setCD(cd);
81      
82      return true;    
83 < }
84 <
85 < bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
86 <
83 >  }
84 >  
85 >  bool ApproximationModel::calcHydroPropsAtCR(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
86 >    
87      int nbeads = beads.size();
88 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
89 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
88 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
89 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
90      Mat3x3d I;
91      I(0, 0) = 1.0;
92      I(1, 1) = 1.0;
93      I(2, 2) = 1.0;
94      
95      for (std::size_t i = 0; i < nbeads; ++i) {
96 <        for (std::size_t j = 0; j < nbeads; ++j) {
97 <            Mat3x3d Tij;
96 >      for (std::size_t j = 0; j < nbeads; ++j) {
97 >        Mat3x3d Tij;
98              if (i != j ) {
99 <                Vector3d Rij = beads[i].pos - beads[j].pos;
100 <                double rij = Rij.length();
101 <                double rij2 = rij * rij;
102 <                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
103 <                Mat3x3d tmpMat;
104 <                tmpMat = outProduct(Rij, Rij) / rij2;
105 <                double constant = 8.0 * NumericConstant::PI * viscosity * rij;
106 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
99 >              Vector3d Rij = beads[i].pos - beads[j].pos;
100 >              RealType rij = Rij.length();
101 >              RealType rij2 = rij * rij;
102 >              RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
103 >              Mat3x3d tmpMat;
104 >              tmpMat = outProduct(Rij, Rij) / rij2;
105 >              RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
106 >              Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
107              }else {
108 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
109 <                Tij(0, 0) = constant;
110 <                Tij(1, 1) = constant;
111 <                Tij(2, 2) = constant;
108 >              RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
109 >              Tij(0, 0) = constant;
110 >              Tij(1, 1) = constant;
111 >              Tij(2, 2) = constant;
112              }
113              B.setSubMatrix(i*3, j*3, Tij);
114 <        }
114 >      }
115      }
116 <
116 >    
117      //invert B Matrix
118      invertMatrix(B, C);
119      
120      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
121      std::vector<Mat3x3d> U;
122      for (int i = 0; i < nbeads; ++i) {
123 <        Mat3x3d currU;
124 <        currU.setupSkewMat(beads[i].pos);
125 <        U.push_back(currU);
123 >      Mat3x3d currU;
124 >      currU.setupSkewMat(beads[i].pos);
125 >      U.push_back(currU);
126      }
127      
128      //calculate Xi matrix at arbitrary origin O
129      Mat3x3d Xiott;
130      Mat3x3d Xiorr;
131      Mat3x3d Xiotr;
132 <
132 >    
133      //calculate the total volume
134 <
135 <    double volume = 0.0;
134 >    
135 >    RealType volume = 0.0;
136      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
137 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
137 >      volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
138      }
139 <        
139 >    
140      for (std::size_t i = 0; i < nbeads; ++i) {
141 <        for (std::size_t j = 0; j < nbeads; ++j) {
142 <            Mat3x3d Cij;
143 <            C.getSubMatrix(i*3, j*3, Cij);
144 <            
145 <            Xiott += Cij;
146 <            Xiotr += U[i] * Cij;
147 <            Xiorr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
148 <        }
141 >      for (std::size_t j = 0; j < nbeads; ++j) {
142 >        Mat3x3d Cij;
143 >        C.getSubMatrix(i*3, j*3, Cij);
144 >        
145 >        Xiott += Cij;
146 >        Xiotr += U[i] * Cij;
147 >        //Xiorr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;    
148 >        Xiorr += -U[i] * Cij * U[j];
149 >      }
150      }
151 <
152 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
151 >    
152 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
153      Xiott *= convertConstant;
154      Xiotr *= convertConstant;
155      Xiorr *= convertConstant;
156      
181 –
157      
158 +    
159      Mat3x3d tmp;
160      Mat3x3d tmpInv;
161      Vector3d tmpVec;
# Line 209 | Line 185 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
185      Xirrr = Xiorr - Uor * Xiott * Uor + Xiotr * Uor - Uor * Xiotr.transpose();
186      
187  
188 <    SquareMatrix<double,6> Xir6x6;
189 <    SquareMatrix<double,6> Dr6x6;
188 >    SquareMatrix<RealType,6> Xir6x6;
189 >    SquareMatrix<RealType,6> Dr6x6;
190  
191      Xir6x6.setSubMatrix(0, 0, Xirtt);
192      Xir6x6.setSubMatrix(0, 3, Xirtr.transpose());
# Line 226 | Line 202 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
202      Dr6x6.getSubMatrix(0, 3, Drrt);
203      Dr6x6.getSubMatrix(3, 0, Drtr);
204      Dr6x6.getSubMatrix(3, 3, Drrr);
205 <    double kt = OOPSEConstant::kB * temperature ;
205 >    RealType kt = OOPSEConstant::kB * temperature ;
206      Drtt *= kt;
207      Drrt *= kt;
208      Drtr *= kt;
# Line 269 | Line 245 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
245  
246      return true;
247   }
248 <
249 < bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, double viscosity, double temperature, HydroProps& cr) {
250 <
248 >  
249 >  bool ApproximationModel::calcHydroPropsAtCD(std::vector<BeadParam>& beads, RealType viscosity, RealType temperature, HydroProps& cr) {
250 >    
251      int nbeads = beads.size();
252 <    DynamicRectMatrix<double> B(3*nbeads, 3*nbeads);
253 <    DynamicRectMatrix<double> C(3*nbeads, 3*nbeads);
252 >    DynamicRectMatrix<RealType> B(3*nbeads, 3*nbeads);
253 >    DynamicRectMatrix<RealType> C(3*nbeads, 3*nbeads);
254      Mat3x3d I;
255      I(0, 0) = 1.0;
256      I(1, 1) = 1.0;
257      I(2, 2) = 1.0;
258      
259      for (std::size_t i = 0; i < nbeads; ++i) {
260 <        for (std::size_t j = 0; j < nbeads; ++j) {
261 <            Mat3x3d Tij;
262 <            if (i != j ) {
263 <                Vector3d Rij = beads[i].pos - beads[j].pos;
264 <                double rij = Rij.length();
265 <                double rij2 = rij * rij;
266 <                double sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
267 <                Mat3x3d tmpMat;
268 <                tmpMat = outProduct(Rij, Rij) / rij2;
269 <                double constant = 8.0 * NumericConstant::PI * viscosity * rij;
270 <                Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
271 <            }else {
272 <                double constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
273 <                Tij(0, 0) = constant;
274 <                Tij(1, 1) = constant;
275 <                Tij(2, 2) = constant;
300 <            }
301 <            B.setSubMatrix(i*3, j*3, Tij);
260 >      for (std::size_t j = 0; j < nbeads; ++j) {
261 >        Mat3x3d Tij;
262 >        if (i != j ) {
263 >          Vector3d Rij = beads[i].pos - beads[j].pos;
264 >          RealType rij = Rij.length();
265 >          RealType rij2 = rij * rij;
266 >          RealType sumSigma2OverRij2 = ((beads[i].radius*beads[i].radius) + (beads[j].radius*beads[j].radius)) / rij2;                
267 >          Mat3x3d tmpMat;
268 >          tmpMat = outProduct(Rij, Rij) / rij2;
269 >          RealType constant = 8.0 * NumericConstant::PI * viscosity * rij;
270 >          Tij = ((1.0 + sumSigma2OverRij2/3.0) * I + (1.0 - sumSigma2OverRij2) * tmpMat ) / constant;
271 >        }else {
272 >          RealType constant = 1.0 / (6.0 * NumericConstant::PI * viscosity * beads[i].radius);
273 >          Tij(0, 0) = constant;
274 >          Tij(1, 1) = constant;
275 >          Tij(2, 2) = constant;
276          }
277 +        B.setSubMatrix(i*3, j*3, Tij);
278 +      }
279      }
280 <
280 >    
281      //invert B Matrix
282      invertMatrix(B, C);
283 <
283 >    
284      //prepare U Matrix relative to arbitrary origin O(0.0, 0.0, 0.0)
285      std::vector<Mat3x3d> U;
286      for (int i = 0; i < nbeads; ++i) {
287 <        Mat3x3d currU;
288 <        currU.setupSkewMat(beads[i].pos);
289 <        U.push_back(currU);
287 >      Mat3x3d currU;
288 >      currU.setupSkewMat(beads[i].pos);
289 >      U.push_back(currU);
290      }
291      
292      //calculate Xi matrix at arbitrary origin O
# Line 320 | Line 296 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
296  
297      //calculate the total volume
298  
299 <    double volume = 0.0;
299 >    RealType volume = 0.0;
300      for (std::vector<BeadParam>::iterator iter = beads.begin(); iter != beads.end(); ++iter) {
301 <        volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
301 >      volume += 4.0/3.0 * NumericConstant::PI * pow((*iter).radius,3);
302      }
303 <        
303 >    
304      for (std::size_t i = 0; i < nbeads; ++i) {
305 <        for (std::size_t j = 0; j < nbeads; ++j) {
306 <            Mat3x3d Cij;
307 <            C.getSubMatrix(i*3, j*3, Cij);
305 >      for (std::size_t j = 0; j < nbeads; ++j) {
306 >        Mat3x3d Cij;
307 >        C.getSubMatrix(i*3, j*3, Cij);
308              
309 <            Xitt += Cij;
310 <            Xitr += U[i] * Cij;
311 <            Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
312 <        }
309 >        Xitt += Cij;
310 >        Xitr += U[i] * Cij;
311 >            //Xirr += -U[i] * Cij * U[j] + (6 * viscosity * volume) * I;            
312 >        Xirr += -U[i] * Cij * U[j];
313 >      }
314      }
315 <
316 <    const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
315 >    
316 >    const RealType convertConstant = 6.023; //convert poise.angstrom to amu/fs
317      Xitt *= convertConstant;
318      Xitr *= convertConstant;
319      Xirr *= convertConstant;
320 <
321 <    double kt = OOPSEConstant::kB * temperature;
322 <
320 >    
321 >    RealType kt = OOPSEConstant::kB * temperature;
322 >    
323      Mat3x3d Dott; //translational diffusion tensor at arbitrary origin O
324      Mat3x3d Dorr; //rotational diffusion tensor at arbitrary origin O
325      Mat3x3d Dotr; //translation-rotation couplingl diffusion tensor at arbitrary origin O
326 <
326 >    
327      const static Mat3x3d zeroMat(0.0);
328      
329      Mat3x3d XittInv(0.0);
# Line 400 | Line 377 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
377      Ddrr = Dorr;
378      Ddtr = Dotr + Dorr * Uod;
379  
380 <    SquareMatrix<double, 6> Dd;
380 >    SquareMatrix<RealType, 6> Dd;
381      Dd.setSubMatrix(0, 0, Ddtt);
382      Dd.setSubMatrix(0, 3, Ddtr.transpose());
383      Dd.setSubMatrix(3, 0, Ddtr);
384      Dd.setSubMatrix(3, 3, Ddrr);    
385 <    SquareMatrix<double, 6> Xid;
385 >    SquareMatrix<RealType, 6> Xid;
386      Ddtt *= kt;
387      Ddtr *=kt;
388      Ddrr *= kt;
# Line 457 | Line 434 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
434      std::cout << Xidrr << std::endl;
435  
436      return true;
437 <      
438 < }
437 >    
438 >  }
439  
440 < /*
464 < void ApproximationModel::writeBeads(std::ostream& os) {
440 >  void ApproximationModel::writeBeads(std::ostream& os) {
441      std::vector<BeadParam>::iterator iter;
442      os << beads_.size() << std::endl;
443      os << "Generated by Hydro" << std::endl;
444      for (iter = beads_.begin(); iter != beads_.end(); ++iter) {
445 <        os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
445 >      os << iter->atomName << "\t" << iter->pos[0] << "\t" << iter->pos[1] << "\t" << iter->pos[2] << std::endl;
446      }
447 <
447 >    
448 >  }    
449   }
473 – */
474 –
475 –
476 – }

Diff Legend

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