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 972 by gezelter, Wed May 24 18:31:12 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, double viscosity, double 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, double viscosity, double 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);
# Line 118 | Line 93 | bool ApproximationModel::calcHydroPropsAtCR(std::vecto
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 >              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;
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 >              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;
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 <
134 >    
135      double 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 <
151 >    
152      const double 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 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, double viscosity, double 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);
# Line 281 | Line 257 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
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 >          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;
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 322 | Line 298 | bool ApproximationModel::calcHydroPropsAtCD(std::vecto
298  
299      double 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 <
315 >    
316      const double convertConstant = 6.023; //convert poise.angstrom to amu/fs
317      Xitt *= convertConstant;
318      Xitr *= convertConstant;
319      Xirr *= convertConstant;
320 <
320 >    
321      double kt = OOPSEConstant::kB * temperature;
322 <
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 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