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 1237 by gezelter, Fri Apr 18 16:55:15 2008 UTC vs.
branches/development/src/integrators/LDForceManager.cpp (file contents), Revision 1874 by gezelter, Wed May 15 15:09:35 2013 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
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, 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>
44   #include "integrators/LDForceManager.hpp"
45   #include "math/CholeskyDecomposition.hpp"
46 < #include "utils/OOPSEConstant.hpp"
46 > #include "utils/PhysicalConstants.hpp"
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 oopse {
53 > namespace OpenMD {
54  
55    LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
56      simParams = info->getSimParams();
# Line 62 | Line 65 | namespace oopse {
65          sprintf( painCave.errMsg,
66                   "langevinBufferRadius must be specified "
67                   "when useSphericalBoundaryConditions is turned on.\n");
68 <        painCave.severity = OOPSE_ERROR;
68 >        painCave.severity = OPENMD_ERROR;
69          painCave.isFatal = 1;
70          simError();  
71        }
# Line 73 | Line 76 | namespace oopse {
76          sprintf( painCave.errMsg,
77                   "frozenBufferRadius must be specified "
78                   "when useSphericalBoundaryConditions is turned on.\n");
79 <        painCave.severity = OOPSE_ERROR;
79 >        painCave.severity = OPENMD_ERROR;
80          painCave.isFatal = 1;
81          simError();  
82        }
# Line 82 | Line 85 | namespace oopse {
85          sprintf( painCave.errMsg,
86                   "frozenBufferRadius has been set smaller than the "
87                   "langevinBufferRadius.  This is probably an error.\n");
88 <        painCave.severity = OOPSE_WARNING;
88 >        painCave.severity = OPENMD_WARNING;
89          painCave.isFatal = 0;
90          simError();  
91        }
# Line 92 | Line 95 | namespace oopse {
95      std::map<std::string, HydroProp*> hydroPropMap;
96  
97      Molecule* mol;
98 <    StuntDouble* integrableObject;
98 >    StuntDouble* sd;
99      SimInfo::MoleculeIterator i;
100      Molecule::IntegrableObjectIterator  j;              
101      bool needHydroPropFile = false;
102      
103      for (mol = info->beginMolecule(i); mol != NULL;
104           mol = info->nextMolecule(i)) {
105 <      for (integrableObject = mol->beginIntegrableObject(j);
106 <           integrableObject != NULL;
107 <           integrableObject = mol->nextIntegrableObject(j)) {
105 >
106 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
107 >           sd = mol->nextIntegrableObject(j)) {
108          
109 <        if (integrableObject->isRigidBody()) {
110 <          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
109 >        if (sd->isRigidBody()) {
110 >          RigidBody* rb = static_cast<RigidBody*>(sd);
111            if (rb->getNumAtoms() > 1) needHydroPropFile = true;
112          }
113          
# Line 120 | Line 123 | namespace oopse {
123                   "HydroPropFile must be set to a file name if Langevin Dynamics\n"
124                   "\tis specified for rigidBodies which contain more than one atom\n"
125                   "\tTo create a HydroPropFile, run the \"Hydro\" program.\n");
126 <        painCave.severity = OOPSE_ERROR;
126 >        painCave.severity = OPENMD_ERROR;
127          painCave.isFatal = 1;
128          simError();  
129        }      
130  
131        for (mol = info->beginMolecule(i); mol != NULL;
132             mol = info->nextMolecule(i)) {
130        for (integrableObject = mol->beginIntegrableObject(j);
131             integrableObject != NULL;
132             integrableObject = mol->nextIntegrableObject(j)) {
133  
134 <          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
134 >        for (sd = mol->beginIntegrableObject(j);  sd != NULL;
135 >             sd = mol->nextIntegrableObject(j)) {
136 >
137 >          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(sd->getType());
138            if (iter != hydroPropMap.end()) {
139              hydroProps_.push_back(iter->second);
140            } else {
141              sprintf( painCave.errMsg,
142 <                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
143 <            painCave.severity = OOPSE_ERROR;
142 >                     "Can not find resistance tensor for atom [%s]\n", sd->getType().c_str());
143 >            painCave.severity = OPENMD_ERROR;
144              painCave.isFatal = 1;
145              simError();  
146            }        
# Line 148 | Line 151 | namespace oopse {
151        std::map<std::string, HydroProp*> hydroPropMap;
152        for (mol = info->beginMolecule(i); mol != NULL;
153             mol = info->nextMolecule(i)) {
154 <        for (integrableObject = mol->beginIntegrableObject(j);
155 <             integrableObject != NULL;
156 <             integrableObject = mol->nextIntegrableObject(j)) {
154 >
155 >        for (sd = mol->beginIntegrableObject(j); sd != NULL;
156 >             sd = mol->nextIntegrableObject(j)) {
157 >
158            Shape* currShape = NULL;
159  
160 <          if (integrableObject->isAtom()){
161 <            Atom* atom = static_cast<Atom*>(integrableObject);
160 >          if (sd->isAtom()){
161 >            Atom* atom = static_cast<Atom*>(sd);
162              AtomType* atomType = atom->getAtomType();
163 <            if (atomType->isGayBerne()) {
164 <              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
165 <              GenericData* data = dAtomType->getPropertyByName("GayBerne");
166 <              if (data != NULL) {
167 <                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 = OOPSE_ERROR;
175 <                  painCave.isFatal = 1;
176 <                  simError();  
177 <                }
178 <              } else {
179 <                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
180 <                painCave.severity = OOPSE_ERROR;
181 <                painCave.isFatal = 1;
182 <                simError();    
183 <              }
163 >            GayBerneAdapter gba = GayBerneAdapter(atomType);
164 >            if (gba.isGayBerne()) {
165 >              currShape = new Ellipsoid(V3Zero, gba.getL() / 2.0,
166 >                                        gba.getD() / 2.0,
167 >                                        Mat3x3d::identity());
168              } else {
169 <              if (atomType->isLennardJones()){
170 <                GenericData* data = atomType->getPropertyByName("LennardJones");
171 <                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 = OOPSE_ERROR;
196 <                    painCave.isFatal = 1;
197 <                    simError();          
198 <                  }      
199 <                }
169 >              LennardJonesAdapter lja = LennardJonesAdapter(atomType);
170 >              if (lja.isLennardJones()){
171 >                currShape = new Sphere(atom->getPos(), lja.getSigma()/2.0);
172                } else {
173                  int aNum = etab.GetAtomicNum((atom->getType()).c_str());
174                  if (aNum != 0) {
# Line 204 | Line 176 | namespace oopse {
176                  } else {
177                    sprintf( painCave.errMsg,
178                             "Could not find atom type in default element.txt\n");
179 <                  painCave.severity = OOPSE_ERROR;
179 >                  painCave.severity = OPENMD_ERROR;
180                    painCave.isFatal = 1;
181                    simError();          
182                  }
183                }
184              }
185            }
186 +
187 +          if (!simParams->haveTargetTemp()) {
188 +            sprintf(painCave.errMsg, "You can't use LangevinDynamics without a targetTemp!\n");
189 +            painCave.isFatal = 1;
190 +            painCave.severity = OPENMD_ERROR;
191 +            simError();
192 +          }
193 +
194 +          if (!simParams->haveViscosity()) {
195 +            sprintf(painCave.errMsg, "You can't use LangevinDynamics without a viscosity!\n");
196 +            painCave.isFatal = 1;
197 +            painCave.severity = OPENMD_ERROR;
198 +            simError();
199 +          }
200 +
201 +
202            HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
203 <          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
203 >          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(sd->getType());
204            if (iter != hydroPropMap.end())
205              hydroProps_.push_back(iter->second);
206            else {
207              currHydroProp->complete();
208 <            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
208 >            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(sd->getType(), currHydroProp));
209              hydroProps_.push_back(currHydroProp);
210            }
211 +          delete currShape;
212          }
213        }
214      }
215 <    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
215 >    variance_ = 2.0 * PhysicalConstants::kb*simParams->getTargetTemp()/simParams->getDt();
216    }  
217  
218    std::map<std::string, HydroProp*> LDForceManager::parseFrictionFile(const std::string& filename) {
# Line 243 | Line 232 | namespace oopse {
232      return props;
233    }
234    
235 <  void LDForceManager::postCalculation(bool needStress){
235 >  void LDForceManager::postCalculation(){
236      SimInfo::MoleculeIterator i;
237      Molecule::IntegrableObjectIterator  j;
238      Molecule* mol;
239 <    StuntDouble* integrableObject;
239 >    StuntDouble* sd;
240      RealType mass;
241      Vector3d pos;
242      Vector3d frc;
# Line 284 | Line 273 | namespace oopse {
273          }
274        }
275        
276 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
277 <           integrableObject = mol->nextIntegrableObject(j)) {
276 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
277 >           sd = mol->nextIntegrableObject(j)) {
278            
279          if (freezeMolecule)
280 <          fdf += integrableObject->freeze();
280 >          fdf += sd->freeze();
281          
282          if (doLangevinForces) {  
283 <          mass = integrableObject->getMass();
284 <          if (integrableObject->isDirectional()){
283 >          mass = sd->getMass();
284 >          if (sd->isDirectional()){
285  
286              // preliminaries for directional objects:
287  
288 <            A = integrableObject->getA();
288 >            A = sd->getA();
289              Atrans = A.transpose();
290              Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
291  
# Line 307 | Line 296 | namespace oopse {
296              genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
297              Vector3d randomForceLab = Atrans * randomForceBody;
298              Vector3d randomTorqueLab = Atrans * randomTorqueBody;
299 <            integrableObject->addFrc(randomForceLab);            
300 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
299 >            sd->addFrc(randomForceLab);            
300 >            sd->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
301  
302 <            Mat3x3d I = integrableObject->getI();
302 >            Mat3x3d I = sd->getI();
303              Vector3d omegaBody;
304  
305              // What remains contains velocity explicitly, but the velocity required
# Line 320 | Line 309 | namespace oopse {
309  
310              // this is the velocity at the half-step:
311              
312 <            Vector3d vel =integrableObject->getVel();
313 <            Vector3d angMom = integrableObject->getJ();
312 >            Vector3d vel =sd->getVel();
313 >            Vector3d angMom = sd->getJ();
314  
315              //estimate velocity at full-step using everything but friction forces:          
316  
317 <            frc = integrableObject->getFrc();
318 <            Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
317 >            frc = sd->getFrc();
318 >            Vector3d velStep = vel + (dt2_ /mass * PhysicalConstants::energyConvert) * frc;
319  
320 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
321 <            Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
320 >            Tb = sd->lab2Body(sd->getTrq());
321 >            Vector3d angMomStep = angMom + (dt2_ * PhysicalConstants::energyConvert) * Tb;                            
322  
323              Vector3d omegaLab;
324              Vector3d vcdLab;
# Line 347 | Line 336 | namespace oopse {
336  
337              for (int k = 0; k < maxIterNum_; k++) {
338                              
339 <              if (integrableObject->isLinear()) {
340 <                int linearAxis = integrableObject->linearAxis();
339 >              if (sd->isLinear()) {
340 >                int linearAxis = sd->linearAxis();
341                  int l = (linearAxis +1 )%3;
342                  int m = (linearAxis +2 )%3;
343                  omegaBody[l] = angMomStep[l] /I(l, l);
# Line 375 | Line 364 | namespace oopse {
364                
365                // re-estimate velocities at full-step using friction forces:
366                
367 <              velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
368 <              angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
367 >              velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * (frc + frictionForceLab);
368 >              angMomStep = angMom + (dt2_ * PhysicalConstants::energyConvert) * (Tb + frictionTorqueBody);
369  
370                // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
371                
# Line 387 | Line 376 | namespace oopse {
376                  break; // iteration ends here
377              }
378  
379 <            integrableObject->addFrc(frictionForceLab);
380 <            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
379 >            sd->addFrc(frictionForceLab);
380 >            sd->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
381  
382              
383            } else {
# Line 397 | Line 386 | namespace oopse {
386              Vector3d randomForce;
387              Vector3d randomTorque;
388              genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
389 <            integrableObject->addFrc(randomForce);            
389 >            sd->addFrc(randomForce);            
390  
391              // What remains contains velocity explicitly, but the velocity required
392              // is at the full step: v(t + h), while we have initially the velocity
# Line 406 | Line 395 | namespace oopse {
395  
396              // this is the velocity at the half-step:
397              
398 <            Vector3d vel =integrableObject->getVel();
398 >            Vector3d vel =sd->getVel();
399  
400              //estimate velocity at full-step using everything but friction forces:          
401  
402 <            frc = integrableObject->getFrc();
403 <            Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
402 >            frc = sd->getFrc();
403 >            Vector3d velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * frc;
404  
405              Vector3d frictionForce(0.0);
406              Vector3d oldFF;  // used to test for convergence
# Line 426 | Line 415 | namespace oopse {
415  
416                // re-estimate velocities at full-step using friction forces:
417                
418 <              velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
418 >              velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * (frc + frictionForce);
419  
420                // check for convergence (if the vector has converged, fdot will be 1.0):
421                
# Line 436 | Line 425 | namespace oopse {
425                  break; // iteration ends here
426              }
427  
428 <            integrableObject->addFrc(frictionForce);
428 >            sd->addFrc(frictionForce);
429  
430            }
431          }
# Line 452 | Line 441 | namespace oopse {
441      if(!simParams->getUsePeriodicBoundaryConditions())
442        veloMunge->removeAngularDrift();
443  
444 <    ForceManager::postCalculation(needStress);  
444 >    ForceManager::postCalculation();  
445    }
446  
447   void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {

Comparing:
trunk/src/integrators/LDForceManager.cpp (property svn:keywords), Revision 1237 by gezelter, Fri Apr 18 16:55:15 2008 UTC vs.
branches/development/src/integrators/LDForceManager.cpp (property svn:keywords), Revision 1874 by gezelter, Wed May 15 15:09:35 2013 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines