ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/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 963 by tim, Wed May 17 21:51:42 2006 UTC

# Line 42 | Line 42
42   #include "integrators/LDForceManager.hpp"
43   #include "math/CholeskyDecomposition.hpp"
44   #include "utils/OOPSEConstant.hpp"
45 + #include "hydrodynamics/Sphere.hpp"
46 + #include "hydrodynamics/Ellipsoid.hpp"
47 + #include "openbabel/mol.hpp"
48 +
49 + using namespace OpenBabel;
50   namespace oopse {
51  
52    LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){
53      Globals* simParams = info->getSimParams();
54 <    
50 <    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 >        
55      sphericalBoundaryConditions_ = false;
56      if (simParams->getUseSphericalBoundaryConditions()) {
57        sphericalBoundaryConditions_ = true;
# Line 92 | Line 86 | namespace oopse {
86          simError();  
87        }
88      }
89 <      
90 <    SimInfo::MoleculeIterator i;
91 <    Molecule::IntegrableObjectIterator  j;
89 >
90 >    // Build the hydroProp map:
91 >    std::map<std::string, HydroProp> hydroPropMap;
92 >
93      Molecule* mol;
94      StuntDouble* integrableObject;
95 <    for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
96 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
95 >    SimInfo::MoleculeIterator i;
96 >    Molecule::IntegrableObjectIterator  j;              
97 >    bool needHydroPropFile = false;
98 >    
99 >    for (mol = info->beginMolecule(i); mol != NULL;
100 >         mol = info->nextMolecule(i)) {
101 >      for (integrableObject = mol->beginIntegrableObject(j);
102 >           integrableObject != NULL;
103             integrableObject = mol->nextIntegrableObject(j)) {
104 <        std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType());
105 <        if (iter != hydroPropMap.end()) {
106 <          hydroProps_.push_back(iter->second);
107 <        } 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();  
104 >        
105 >        if (integrableObject->isRigidBody()) {
106 >          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
107 >          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
108          }
109          
110        }
111      }
112 +        
113 +
114 +    if (needHydroPropFile) {              
115 +      if (simParams->haveHydroPropFile()) {
116 +        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
117 +      } else {              
118 +        sprintf( painCave.errMsg,
119 +                 "HydroPropFile must be set to a file name if Langevin\n"
120 +                 "\tDynamics is specified for rigidBodies which contain more\n"
121 +                 "\tthan one atom.  To create a HydroPropFile, run \"Hydro\".\n");
122 +        painCave.severity = OOPSE_ERROR;
123 +        painCave.isFatal = 1;
124 +        simError();  
125 +      }      
126 +      std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType());
127 +      if (iter != hydroPropMap.end()) {
128 +        hydroProps_.push_back(iter->second);
129 +      } else {
130 +        sprintf( painCave.errMsg,
131 +                 "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
132 +        painCave.severity = OOPSE_ERROR;
133 +        painCave.isFatal = 1;
134 +        simError();  
135 +      }
136 +    } else {
137 +
138 +      std::map<std::string, HydroProp> hydroPropMap;
139 +      for (mol = info->beginMolecule(i); mol != NULL;
140 +           mol = info->nextMolecule(i)) {
141 +        for (integrableObject = mol->beginIntegrableObject(j);
142 +             integrableObject != NULL;
143 +             integrableObject = mol->nextIntegrableObject(j)) {
144 +          Shape* currShape = NULL;
145 +          if (integrableObject->isDirectionalAtom()) {
146 +            DirectionalAtom* dAtom = static_cast<DirectionalAtom*>(integrableObject);
147 +            AtomType* atomType = dAtom->getAtomType();
148 +            if (atomType->isGayBerne()) {
149 +              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
150 +              
151 +              GenericData* data = dAtomType->getPropertyByName("GayBerne");
152 +              if (data != NULL) {
153 +                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
154 +                
155 +                if (gayBerneData != NULL) {  
156 +                  GayBerneParam gayBerneParam = gayBerneData->getData();
157 +                  currShape = new Ellipsoid(V3Zero,
158 +                                            gayBerneParam.GB_sigma/2.0,
159 +                                            gayBerneParam.GB_l2b_ratio*gayBerneParam.GB_sigma/2.0,
160 +                                            Mat3x3d::identity());
161 +                } else {
162 +                  sprintf( painCave.errMsg,
163 +                           "Can not cast GenericData to GayBerneParam\n");
164 +                  painCave.severity = OOPSE_ERROR;
165 +                  painCave.isFatal = 1;
166 +                  simError();  
167 +                }
168 +              } else {
169 +                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
170 +                painCave.severity = OOPSE_ERROR;
171 +                painCave.isFatal = 1;
172 +                simError();    
173 +              }
174 +            }
175 +          } else {
176 +            Atom* atom = static_cast<Atom*>(integrableObject);
177 +            AtomType* atomType = atom->getAtomType();
178 +            if (atomType->isLennardJones()){
179 +              GenericData* data = atomType->getPropertyByName("LennardJones");
180 +              if (data != NULL) {
181 +                LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
182 +                
183 +                if (ljData != NULL) {
184 +                  LJParam ljParam = ljData->getData();
185 +                  currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
186 +                } else {
187 +                  sprintf( painCave.errMsg,
188 +                           "Can not cast GenericData to LJParam\n");
189 +                  painCave.severity = OOPSE_ERROR;
190 +                  painCave.isFatal = 1;
191 +                  simError();          
192 +                }      
193 +              }
194 +            } else {
195 +              int obanum = etab.GetAtomicNum((atom->getType()).c_str());
196 +              if (obanum != 0) {
197 +                currShape = new Sphere(atom->getPos(), etab.GetVdwRad(obanum));
198 +              } else {
199 +                sprintf( painCave.errMsg,
200 +                         "Could not find atom type in default element.txt\n");
201 +                painCave.severity = OOPSE_ERROR;
202 +                painCave.isFatal = 1;
203 +                simError();          
204 +              }
205 +            }
206 +          }
207 +          HydroProps currHydroProp = currShape->getHydroProps(simParams->getViscosity(),simParams->getTargetTemp());
208 +          std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType());
209 +          if (iter != hydroPropMap.end())
210 +            hydroProps_.push_back(iter->second);
211 +          else {
212 +            HydroProp myProp;
213 +            myProp.cor = V3Zero;
214 +            for (int i1 = 0; i1 < 3; i1++) {
215 +              for (int j1 = 0; j1 < 3; j1++) {
216 +                myProp.Xirtt(i1,j1) = currHydroProp.Xi(i1,j1);
217 +                myProp.Xirrt(i1,j1) = currHydroProp.Xi(i1,j1+3);
218 +                myProp.Xirtr(i1,j1) = currHydroProp.Xi(i1+3,j1);
219 +                myProp.Xirrr(i1,j1) = currHydroProp.Xi(i1+3,j1+3);
220 +              }
221 +            }
222 +            CholeskyDecomposition(currHydroProp.Xi, myProp.S);
223 +            hydroPropMap.insert(std::map<std::string, HydroProp>::value_type(integrableObject->getType(), myProp));
224 +            hydroProps_.push_back(myProp);
225 +          }
226 +        }
227 +      }
228 +    }
229      variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
230    }
231 +  
232 +  
233 +
234 +
235 +
236    std::map<std::string, HydroProp> LDForceManager::parseFrictionFile(const std::string& filename) {
237      std::map<std::string, HydroProp> props;
238      std::ifstream ifs(filename.c_str());
# Line 173 | Line 291 | namespace oopse {
291          currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
292          currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
293          
294 <        SquareMatrix<double, 6> Xir;
294 >        SquareMatrix<RealType, 6> Xir;
295          Xir.setSubMatrix(0, 0, currProp.Xirtt);
296          Xir.setSubMatrix(0, 3, currProp.Xirrt);
297          Xir.setSubMatrix(3, 0, currProp.Xirtr);
# Line 199 | Line 317 | namespace oopse {
317      Mat3x3d Atrans;
318      Vector3d Tb;
319      Vector3d ji;
320 <    double mass;
320 >    RealType mass;
321      unsigned int index = 0;
322      bool doLangevinForces;
323      bool freezeMolecule;
# Line 211 | Line 329 | namespace oopse {
329        if (sphericalBoundaryConditions_) {
330          
331          Vector3d molPos = mol->getCom();
332 <        double molRad = molPos.length();
332 >        RealType molRad = molPos.length();
333          
334          doLangevinForces = false;
335          freezeMolecule = false;
# Line 226 | Line 344 | namespace oopse {
344          }
345        }
346        
347 <      if (doLangevinForces) {
348 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
231 <             integrableObject = mol->nextIntegrableObject(j)) {
347 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
348 >           integrableObject = mol->nextIntegrableObject(j)) {
349            
350 +        if (freezeMolecule)
351 +          fdf += integrableObject->freeze();
352 +        
353 +        if (doLangevinForces) {          
354            vel =integrableObject->getVel();
355            if (integrableObject->isDirectional()){
356              //calculate angular velocity in lab frame
# Line 281 | Line 402 | namespace oopse {
402              
403              integrableObject->addFrc(frictionForce+randomForce);            
404            }
405 +        }
406            
407 <          ++index;
407 >        ++index;
408      
287 –        }
409        }
410 <      if (freezeMolecule)
290 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
291 <             integrableObject = mol->nextIntegrableObject(j)) {          
292 <          fdf += integrableObject->freeze();
293 <        }
294 <    }
295 <    
410 >    }    
411      info_->setFdf(fdf);
412      
413      ForceManager::postCalculation();  
414    }
415  
416 < void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, double variance) {
416 > void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
417  
418  
419 <    Vector<double, 6> Z;
420 <    Vector<double, 6> generalForce;
419 >    Vector<RealType, 6> Z;
420 >    Vector<RealType, 6> generalForce;
421  
422          
423      Z[0] = randNumGen_.randNorm(0, variance);

Diff Legend

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