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 1610 by gezelter, Fri Aug 12 14:37:25 2011 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42   #include <fstream>
43   #include <iostream>
# Line 46 | Line 47
47   #include "hydrodynamics/Sphere.hpp"
48   #include "hydrodynamics/Ellipsoid.hpp"
49   #include "utils/ElementsTable.hpp"
50 + #include "types/LennardJonesAdapter.hpp"
51 + #include "types/GayBerneAdapter.hpp"
52  
53   namespace OpenMD {
54  
55 <  LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
55 >  LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info),
56 >                                                  maxIterNum_(4),
57 >                                                  forceTolerance_(1e-6) {
58      simParams = info->getSimParams();
59      veloMunge = new Velocitizer(info);
60  
# Line 92 | Line 97 | namespace OpenMD {
97      std::map<std::string, HydroProp*> hydroPropMap;
98  
99      Molecule* mol;
100 <    StuntDouble* integrableObject;
100 >    StuntDouble* sd;
101      SimInfo::MoleculeIterator i;
102      Molecule::IntegrableObjectIterator  j;              
103      bool needHydroPropFile = false;
104      
105      for (mol = info->beginMolecule(i); mol != NULL;
106           mol = info->nextMolecule(i)) {
107 <      for (integrableObject = mol->beginIntegrableObject(j);
108 <           integrableObject != NULL;
109 <           integrableObject = mol->nextIntegrableObject(j)) {
107 >
108 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
109 >           sd = mol->nextIntegrableObject(j)) {
110          
111 <        if (integrableObject->isRigidBody()) {
112 <          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
111 >        if (sd->isRigidBody()) {
112 >          RigidBody* rb = static_cast<RigidBody*>(sd);
113            if (rb->getNumAtoms() > 1) needHydroPropFile = true;
114          }
115          
# Line 127 | Line 132 | namespace OpenMD {
132  
133        for (mol = info->beginMolecule(i); mol != NULL;
134             mol = info->nextMolecule(i)) {
130 –        for (integrableObject = mol->beginIntegrableObject(j);
131 –             integrableObject != NULL;
132 –             integrableObject = mol->nextIntegrableObject(j)) {
135  
136 <          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
136 >        for (sd = mol->beginIntegrableObject(j);  sd != NULL;
137 >             sd = mol->nextIntegrableObject(j)) {
138 >
139 >          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(sd->getType());
140            if (iter != hydroPropMap.end()) {
141              hydroProps_.push_back(iter->second);
142            } else {
143              sprintf( painCave.errMsg,
144 <                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
144 >                     "Can not find resistance tensor for atom [%s]\n", sd->getType().c_str());
145              painCave.severity = OPENMD_ERROR;
146              painCave.isFatal = 1;
147              simError();  
# Line 148 | Line 153 | namespace OpenMD {
153        std::map<std::string, HydroProp*> hydroPropMap;
154        for (mol = info->beginMolecule(i); mol != NULL;
155             mol = info->nextMolecule(i)) {
156 <        for (integrableObject = mol->beginIntegrableObject(j);
157 <             integrableObject != NULL;
158 <             integrableObject = mol->nextIntegrableObject(j)) {
156 >
157 >        for (sd = mol->beginIntegrableObject(j); sd != NULL;
158 >             sd = mol->nextIntegrableObject(j)) {
159 >
160            Shape* currShape = NULL;
161  
162 <          if (integrableObject->isAtom()){
163 <            Atom* atom = static_cast<Atom*>(integrableObject);
162 >          if (sd->isAtom()){
163 >            Atom* atom = static_cast<Atom*>(sd);
164              AtomType* atomType = atom->getAtomType();
165 <            if (atomType->isGayBerne()) {
166 <              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
167 <              GenericData* data = dAtomType->getPropertyByName("GayBerne");
168 <              if (data != NULL) {
169 <                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
164 <                
165 <                if (gayBerneData != NULL) {  
166 <                  GayBerneParam gayBerneParam = gayBerneData->getData();
167 <                  currShape = new Ellipsoid(V3Zero,
168 <                                            gayBerneParam.GB_l / 2.0,
169 <                                            gayBerneParam.GB_d / 2.0,
170 <                                            Mat3x3d::identity());
171 <                } else {
172 <                  sprintf( painCave.errMsg,
173 <                           "Can not cast GenericData to GayBerneParam\n");
174 <                  painCave.severity = OPENMD_ERROR;
175 <                  painCave.isFatal = 1;
176 <                  simError();  
177 <                }
178 <              } else {
179 <                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
180 <                painCave.severity = OPENMD_ERROR;
181 <                painCave.isFatal = 1;
182 <                simError();    
183 <              }
165 >            GayBerneAdapter gba = GayBerneAdapter(atomType);
166 >            if (gba.isGayBerne()) {
167 >              currShape = new Ellipsoid(V3Zero, gba.getL() / 2.0,
168 >                                        gba.getD() / 2.0,
169 >                                        Mat3x3d::identity());
170              } else {
171 <              if (atomType->isLennardJones()){
172 <                GenericData* data = atomType->getPropertyByName("LennardJones");
173 <                if (data != NULL) {
188 <                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
189 <                  if (ljData != NULL) {
190 <                    LJParam ljParam = ljData->getData();
191 <                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
192 <                  } else {
193 <                    sprintf( painCave.errMsg,
194 <                             "Can not cast GenericData to LJParam\n");
195 <                    painCave.severity = OPENMD_ERROR;
196 <                    painCave.isFatal = 1;
197 <                    simError();          
198 <                  }      
199 <                }
171 >              LennardJonesAdapter lja = LennardJonesAdapter(atomType);
172 >              if (lja.isLennardJones()){
173 >                currShape = new Sphere(atom->getPos(), lja.getSigma()/2.0);
174                } else {
175                  int aNum = etab.GetAtomicNum((atom->getType()).c_str());
176                  if (aNum != 0) {
# Line 226 | Line 200 | namespace OpenMD {
200              simError();
201            }
202  
229 –          HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
203  
204 <          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
204 >          HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
205 >          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(sd->getType());
206            if (iter != hydroPropMap.end())
207              hydroProps_.push_back(iter->second);
208            else {
209              currHydroProp->complete();
210 <            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
210 >            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(sd->getType(), currHydroProp));
211              hydroProps_.push_back(currHydroProp);
212            }
213 +          delete currShape;
214          }
215        }
216      }
# Line 263 | Line 238 | namespace OpenMD {
238      SimInfo::MoleculeIterator i;
239      Molecule::IntegrableObjectIterator  j;
240      Molecule* mol;
241 <    StuntDouble* integrableObject;
241 >    StuntDouble* sd;
242      RealType mass;
243      Vector3d pos;
244      Vector3d frc;
# Line 300 | Line 275 | namespace OpenMD {
275          }
276        }
277        
278 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
279 <           integrableObject = mol->nextIntegrableObject(j)) {
278 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
279 >           sd = mol->nextIntegrableObject(j)) {
280            
281          if (freezeMolecule)
282 <          fdf += integrableObject->freeze();
282 >          fdf += sd->freeze();
283          
284          if (doLangevinForces) {  
285 <          mass = integrableObject->getMass();
286 <          if (integrableObject->isDirectional()){
285 >          mass = sd->getMass();
286 >          if (sd->isDirectional()){
287  
288              // preliminaries for directional objects:
289  
290 <            A = integrableObject->getA();
290 >            A = sd->getA();
291              Atrans = A.transpose();
292              Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
293  
# Line 323 | Line 298 | namespace OpenMD {
298              genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
299              Vector3d randomForceLab = Atrans * randomForceBody;
300              Vector3d randomTorqueLab = Atrans * randomTorqueBody;
301 +            sd->addFrc(randomForceLab);            
302 +            sd->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
303  
304 <            integrableObject->addFrc(randomForceLab);            
328 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
329 <
330 <            Mat3x3d I = integrableObject->getI();
304 >            Mat3x3d I = sd->getI();
305              Vector3d omegaBody;
306  
307              // What remains contains velocity explicitly, but the velocity required
# Line 337 | Line 311 | namespace OpenMD {
311  
312              // this is the velocity at the half-step:
313              
314 <            Vector3d vel =integrableObject->getVel();
315 <            Vector3d angMom = integrableObject->getJ();
314 >            Vector3d vel =sd->getVel();
315 >            Vector3d angMom = sd->getJ();
316  
317              //estimate velocity at full-step using everything but friction forces:          
318  
319 <            frc = integrableObject->getFrc();
319 >            frc = sd->getFrc();
320              Vector3d velStep = vel + (dt2_ /mass * PhysicalConstants::energyConvert) * frc;
321  
322 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
322 >            Tb = sd->lab2Body(sd->getTrq());
323              Vector3d angMomStep = angMom + (dt2_ * PhysicalConstants::energyConvert) * Tb;                            
324  
325              Vector3d omegaLab;
# Line 364 | Line 338 | namespace OpenMD {
338  
339              for (int k = 0; k < maxIterNum_; k++) {
340                              
341 <              if (integrableObject->isLinear()) {
342 <                int linearAxis = integrableObject->linearAxis();
341 >              if (sd->isLinear()) {
342 >                int linearAxis = sd->linearAxis();
343                  int l = (linearAxis +1 )%3;
344                  int m = (linearAxis +2 )%3;
345                  omegaBody[l] = angMomStep[l] /I(l, l);
# Line 404 | Line 378 | namespace OpenMD {
378                  break; // iteration ends here
379              }
380  
381 <            integrableObject->addFrc(frictionForceLab);
382 <            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
381 >            sd->addFrc(frictionForceLab);
382 >            sd->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
383  
384              
385            } else {
# Line 414 | Line 388 | namespace OpenMD {
388              Vector3d randomForce;
389              Vector3d randomTorque;
390              genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
391 <            integrableObject->addFrc(randomForce);            
391 >            sd->addFrc(randomForce);            
392  
393              // What remains contains velocity explicitly, but the velocity required
394              // is at the full step: v(t + h), while we have initially the velocity
# Line 423 | Line 397 | namespace OpenMD {
397  
398              // this is the velocity at the half-step:
399              
400 <            Vector3d vel =integrableObject->getVel();
400 >            Vector3d vel =sd->getVel();
401  
402              //estimate velocity at full-step using everything but friction forces:          
403  
404 <            frc = integrableObject->getFrc();
404 >            frc = sd->getFrc();
405              Vector3d velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * frc;
406  
407              Vector3d frictionForce(0.0);
# Line 453 | Line 427 | namespace OpenMD {
427                  break; // iteration ends here
428              }
429  
430 <            integrableObject->addFrc(frictionForce);
430 >            sd->addFrc(frictionForce);
431  
432            }
433          }

Diff Legend

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