ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/integrators/LDForceManager.cpp
(Generate patch)

Comparing trunk/src/integrators/LDForceManager.cpp (file contents):
Revision 945 by gezelter, Tue Apr 25 02:09:01 2006 UTC vs.
Revision 1120 by chuckv, Fri Feb 2 18:55:21 2007 UTC

# Line 39 | Line 39
39   * such damages.
40   */
41   #include <fstream>
42 + #include <iostream>
43   #include "integrators/LDForceManager.hpp"
44   #include "math/CholeskyDecomposition.hpp"
45   #include "utils/OOPSEConstant.hpp"
46 + #include "hydrodynamics/Sphere.hpp"
47 + #include "hydrodynamics/Ellipsoid.hpp"
48 + #include "openbabel/mol.hpp"
49 +
50 + using namespace OpenBabel;
51   namespace oopse {
52  
53    LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){
54 <    Globals* simParams = info->getSimParams();
55 <    
56 <    std::map<std::string, HydroProp> hydroPropMap;
51 <    if (simParams->haveHydroPropFile()) {
52 <      hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
53 <    } else {
54 <      sprintf( painCave.errMsg,
55 <               "HydroPropFile must be set if Langevin Dynamics is specified.\n");
56 <      painCave.severity = OOPSE_ERROR;
57 <      painCave.isFatal = 1;
58 <      simError();  
59 <    }
60 <    
54 >    simParams = info->getSimParams();
55 >    veloMunge = new Velocitizer(info);
56 >
57      sphericalBoundaryConditions_ = false;
58      if (simParams->getUseSphericalBoundaryConditions()) {
59        sphericalBoundaryConditions_ = true;
# Line 92 | Line 88 | namespace oopse {
88          simError();  
89        }
90      }
91 <      
92 <    SimInfo::MoleculeIterator i;
93 <    Molecule::IntegrableObjectIterator  j;
91 >
92 >    // Build the hydroProp map:
93 >    std::map<std::string, HydroProp*> hydroPropMap;
94 >
95      Molecule* mol;
96      StuntDouble* integrableObject;
97 <    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
98 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
97 >    SimInfo::MoleculeIterator i;
98 >    Molecule::IntegrableObjectIterator  j;              
99 >    bool needHydroPropFile = false;
100 >    
101 >    for (mol = info->beginMolecule(i); mol != NULL;
102 >         mol = info->nextMolecule(i)) {
103 >      for (integrableObject = mol->beginIntegrableObject(j);
104 >           integrableObject != NULL;
105             integrableObject = mol->nextIntegrableObject(j)) {
106 <        std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType());
107 <        if (iter != hydroPropMap.end()) {
108 <          hydroProps_.push_back(iter->second);
109 <        } else {
107 <          sprintf( painCave.errMsg,
108 <                   "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
109 <          painCave.severity = OOPSE_ERROR;
110 <          painCave.isFatal = 1;
111 <          simError();  
106 >        
107 >        if (integrableObject->isRigidBody()) {
108 >          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
109 >          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
110          }
111          
112        }
113      }
114 +        
115 +
116 +    if (needHydroPropFile) {              
117 +      if (simParams->haveHydroPropFile()) {
118 +        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
119 +      } else {              
120 +        sprintf( painCave.errMsg,
121 +                 "HydroPropFile must be set to a file name if Langevin\n"
122 +                 "\tDynamics is specified for rigidBodies which contain more\n"
123 +                 "\tthan one atom.  To create a HydroPropFile, run \"Hydro\".\n");
124 +        painCave.severity = OOPSE_ERROR;
125 +        painCave.isFatal = 1;
126 +        simError();  
127 +      }      
128 +
129 +      for (mol = info->beginMolecule(i); mol != NULL;
130 +           mol = info->nextMolecule(i)) {
131 +        for (integrableObject = mol->beginIntegrableObject(j);
132 +             integrableObject != NULL;
133 +             integrableObject = mol->nextIntegrableObject(j)) {
134 +
135 +          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
136 +          if (iter != hydroPropMap.end()) {
137 +            hydroProps_.push_back(iter->second);
138 +          } else {
139 +            sprintf( painCave.errMsg,
140 +                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
141 +            painCave.severity = OOPSE_ERROR;
142 +            painCave.isFatal = 1;
143 +            simError();  
144 +          }        
145 +        }
146 +      }
147 +    } else {
148 +      
149 +      std::map<std::string, HydroProp*> hydroPropMap;
150 +      for (mol = info->beginMolecule(i); mol != NULL;
151 +           mol = info->nextMolecule(i)) {
152 +        for (integrableObject = mol->beginIntegrableObject(j);
153 +             integrableObject != NULL;
154 +             integrableObject = mol->nextIntegrableObject(j)) {
155 +          Shape* currShape = NULL;
156 +          if (integrableObject->isDirectionalAtom()) {
157 +            DirectionalAtom* dAtom = static_cast<DirectionalAtom*>(integrableObject);
158 +            AtomType* atomType = dAtom->getAtomType();
159 +            if (atomType->isGayBerne()) {
160 +              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
161 +              
162 +              GenericData* data = dAtomType->getPropertyByName("GayBerne");
163 +              if (data != NULL) {
164 +                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
165 +                
166 +                if (gayBerneData != NULL) {  
167 +                  GayBerneParam gayBerneParam = gayBerneData->getData();
168 +                  currShape = new Ellipsoid(V3Zero,
169 +                                            gayBerneParam.GB_d / 2.0,
170 +                                            gayBerneParam.GB_l / 2.0,
171 +                                            Mat3x3d::identity());
172 +                } else {
173 +                  sprintf( painCave.errMsg,
174 +                           "Can not cast GenericData to GayBerneParam\n");
175 +                  painCave.severity = OOPSE_ERROR;
176 +                  painCave.isFatal = 1;
177 +                  simError();  
178 +                }
179 +              } else {
180 +                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
181 +                painCave.severity = OOPSE_ERROR;
182 +                painCave.isFatal = 1;
183 +                simError();    
184 +              }
185 +            }
186 +          } else {
187 +            Atom* atom = static_cast<Atom*>(integrableObject);
188 +            AtomType* atomType = atom->getAtomType();
189 +            if (atomType->isLennardJones()){
190 +              GenericData* data = atomType->getPropertyByName("LennardJones");
191 +              if (data != NULL) {
192 +                LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
193 +                
194 +                if (ljData != NULL) {
195 +                  LJParam ljParam = ljData->getData();
196 +                  currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
197 +                } else {
198 +                  sprintf( painCave.errMsg,
199 +                           "Can not cast GenericData to LJParam\n");
200 +                  painCave.severity = OOPSE_ERROR;
201 +                  painCave.isFatal = 1;
202 +                  simError();          
203 +                }      
204 +              }
205 +            } else {
206 +              int obanum = etab.GetAtomicNum((atom->getType()).c_str());
207 +              if (obanum != 0) {
208 +                currShape = new Sphere(atom->getPos(), etab.GetVdwRad(obanum));
209 +              } else {
210 +                sprintf( painCave.errMsg,
211 +                         "Could not find atom type in default element.txt\n");
212 +                painCave.severity = OOPSE_ERROR;
213 +                painCave.isFatal = 1;
214 +                simError();          
215 +              }
216 +            }
217 +          }
218 +          HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
219 +          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
220 +          if (iter != hydroPropMap.end())
221 +            hydroProps_.push_back(iter->second);
222 +          else {
223 +            currHydroProp->complete();
224 +            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
225 +            hydroProps_.push_back(currHydroProp);
226 +          }
227 +        }
228 +      }
229 +    }
230      variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
231 <  }
232 <  std::map<std::string, HydroProp> LDForceManager::parseFrictionFile(const std::string& filename) {
233 <    std::map<std::string, HydroProp> props;
231 >  }  
232 >
233 >  std::map<std::string, HydroProp*> LDForceManager::parseFrictionFile(const std::string& filename) {
234 >    std::map<std::string, HydroProp*> props;
235      std::ifstream ifs(filename.c_str());
236      if (ifs.is_open()) {
237        
# Line 125 | Line 240 | namespace oopse {
240      const unsigned int BufferSize = 65535;
241      char buffer[BufferSize];  
242      while (ifs.getline(buffer, BufferSize)) {
243 <      StringTokenizer tokenizer(buffer);
244 <      HydroProp currProp;
130 <      if (tokenizer.countTokens() >= 40) {
131 <        std::string atomName = tokenizer.nextToken();
132 <        currProp.cor[0] = tokenizer.nextTokenAsDouble();
133 <        currProp.cor[1] = tokenizer.nextTokenAsDouble();
134 <        currProp.cor[2] = tokenizer.nextTokenAsDouble();
135 <        
136 <        currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble();
137 <        currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble();
138 <        currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble();
139 <        currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble();
140 <        currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble();
141 <        currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble();
142 <        currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble();
143 <        currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble();
144 <        currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble();
145 <        
146 <        currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble();
147 <        currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble();
148 <        currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble();
149 <        currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble();
150 <        currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble();
151 <        currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble();
152 <        currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble();
153 <        currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble();
154 <        currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble();
155 <        
156 <        currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble();
157 <        currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble();
158 <        currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble();
159 <        currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble();
160 <        currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble();
161 <        currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble();
162 <        currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble();
163 <        currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble();
164 <        currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble();
165 <        
166 <        currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble();
167 <        currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble();
168 <        currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble();
169 <        currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble();
170 <        currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble();
171 <        currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble();
172 <        currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble();
173 <        currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
174 <        currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
175 <        
176 <        SquareMatrix<double, 6> Xir;
177 <        Xir.setSubMatrix(0, 0, currProp.Xirtt);
178 <        Xir.setSubMatrix(0, 3, currProp.Xirrt);
179 <        Xir.setSubMatrix(3, 0, currProp.Xirtr);
180 <        Xir.setSubMatrix(3, 3, currProp.Xirrr);
181 <        CholeskyDecomposition(Xir, currProp.S);            
182 <        
183 <        props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
184 <      }
243 >      HydroProp* currProp = new HydroProp(buffer);
244 >      props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
245      }
246 <    
246 >
247      return props;
248    }
249 <  
249 >  
250    void LDForceManager::postCalculation() {
251      SimInfo::MoleculeIterator i;
252      Molecule::IntegrableObjectIterator  j;
# Line 199 | Line 259 | namespace oopse {
259      Mat3x3d Atrans;
260      Vector3d Tb;
261      Vector3d ji;
262 <    double mass;
262 >    RealType mass;
263      unsigned int index = 0;
264      bool doLangevinForces;
265      bool freezeMolecule;
266      int fdf;
267 <    
267 >    int nIntegrated;
268 >    int nFrozen;
269 >
270      fdf = 0;
271 +
272      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
273 <      
273 >
274 >      doLangevinForces = true;          
275 >      freezeMolecule = false;
276 >
277        if (sphericalBoundaryConditions_) {
278          
279          Vector3d molPos = mol->getCom();
280 <        double molRad = molPos.length();
281 <        
280 >        RealType molRad = molPos.length();
281 >
282          doLangevinForces = false;
217 –        freezeMolecule = false;
283          
284          if (molRad > langevinBufferRadius_) {
285            doLangevinForces = true;
# Line 226 | Line 291 | namespace oopse {
291          }
292        }
293        
294 <      if (doLangevinForces) {
295 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
231 <             integrableObject = mol->nextIntegrableObject(j)) {
294 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
295 >           integrableObject = mol->nextIntegrableObject(j)) {
296            
297 +        if (freezeMolecule)
298 +          fdf += integrableObject->freeze();
299 +        
300 +        if (doLangevinForces) {  
301            vel =integrableObject->getVel();
302            if (integrableObject->isDirectional()){
303              //calculate angular velocity in lab frame
# Line 253 | Line 321 | namespace oopse {
321              //apply friction force and torque at center of resistance
322              A = integrableObject->getA();
323              Atrans = A.transpose();
324 <            Vector3d rcr = Atrans * hydroProps_[index].cor;  
324 >            Vector3d rcr = Atrans * hydroProps_[index]->getCOR();  
325              Vector3d vcdLab = vel + cross(omega, rcr);
326              Vector3d vcdBody = A* vcdLab;
327 <            Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega);
327 >            Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omega);
328              Vector3d frictionForceLab = Atrans*frictionForceBody;
329              integrableObject->addFrc(frictionForceLab);
330 <            Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega);
330 >            Vector3d frictionTorqueBody = - (hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omega);
331              Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
332              integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
333              
# Line 274 | Line 342 | namespace oopse {
342              
343            } else {
344              //spherical atom
345 <            Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel);    
345 >            Vector3d frictionForce = -(hydroProps_[index]->getXitt() * vel);
346              Vector3d randomForce;
347              Vector3d randomTorque;
348              genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
349              
350              integrableObject->addFrc(frictionForce+randomForce);            
351            }
352 +        }
353            
354 <          ++index;
354 >        ++index;
355      
287 –        }
356        }
357 <      if (freezeMolecule)
358 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
291 <             integrableObject = mol->nextIntegrableObject(j)) {          
292 <          fdf += integrableObject->freeze();
293 <        }
294 <    }
295 <    
357 >    }    
358 >
359      info_->setFdf(fdf);
360 <    
360 >    veloMunge->removeComDrift();
361 >    // Remove angular drift if we are not using periodic boundary conditions.
362 >    if(!simParams->getUsePeriodicBoundaryConditions())
363 >      veloMunge->removeAngularDrift();
364 >
365      ForceManager::postCalculation();  
366    }
367  
368 < void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, double variance) {
368 > void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
369  
370  
371 <    Vector<double, 6> Z;
372 <    Vector<double, 6> generalForce;
371 >    Vector<RealType, 6> Z;
372 >    Vector<RealType, 6> generalForce;
373  
374          
375      Z[0] = randNumGen_.randNorm(0, variance);
# Line 313 | Line 380 | void LDForceManager::genRandomForceAndTorque(Vector3d&
380      Z[5] = randNumGen_.randNorm(0, variance);
381      
382  
383 <    generalForce = hydroProps_[index].S*Z;
383 >    generalForce = hydroProps_[index]->getS()*Z;
384      
385      force[0] = generalForce[0];
386      force[1] = generalForce[1];

Diff Legend

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