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

Comparing branches/development/src/integrators/LDForceManager.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1769 by gezelter, Mon Jul 9 14:15:52 2012 UTC

# Line 36 | Line 36
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).                        
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  
# Line 92 | Line 95 | namespace OpenMD {
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 127 | Line 130 | namespace OpenMD {
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());
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();  
# Line 148 | Line 151 | namespace OpenMD {
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 = 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 <              }
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 = OPENMD_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 228 | Line 200 | namespace OpenMD {
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          }
# Line 263 | Line 235 | namespace OpenMD {
235      SimInfo::MoleculeIterator i;
236      Molecule::IntegrableObjectIterator  j;
237      Molecule* mol;
238 <    StuntDouble* integrableObject;
238 >    StuntDouble* sd;
239      RealType mass;
240      Vector3d pos;
241      Vector3d frc;
# Line 300 | Line 272 | namespace OpenMD {
272          }
273        }
274        
275 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
276 <           integrableObject = mol->nextIntegrableObject(j)) {
275 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
276 >           sd = mol->nextIntegrableObject(j)) {
277            
278          if (freezeMolecule)
279 <          fdf += integrableObject->freeze();
279 >          fdf += sd->freeze();
280          
281          if (doLangevinForces) {  
282 <          mass = integrableObject->getMass();
283 <          if (integrableObject->isDirectional()){
282 >          mass = sd->getMass();
283 >          if (sd->isDirectional()){
284  
285              // preliminaries for directional objects:
286  
287 <            A = integrableObject->getA();
287 >            A = sd->getA();
288              Atrans = A.transpose();
289              Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
290  
# Line 323 | Line 295 | namespace OpenMD {
295              genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
296              Vector3d randomForceLab = Atrans * randomForceBody;
297              Vector3d randomTorqueLab = Atrans * randomTorqueBody;
298 <            integrableObject->addFrc(randomForceLab);            
299 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
298 >            sd->addFrc(randomForceLab);            
299 >            sd->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
300  
301 <            Mat3x3d I = integrableObject->getI();
301 >            Mat3x3d I = sd->getI();
302              Vector3d omegaBody;
303  
304              // What remains contains velocity explicitly, but the velocity required
# Line 336 | Line 308 | namespace OpenMD {
308  
309              // this is the velocity at the half-step:
310              
311 <            Vector3d vel =integrableObject->getVel();
312 <            Vector3d angMom = integrableObject->getJ();
311 >            Vector3d vel =sd->getVel();
312 >            Vector3d angMom = sd->getJ();
313  
314              //estimate velocity at full-step using everything but friction forces:          
315  
316 <            frc = integrableObject->getFrc();
316 >            frc = sd->getFrc();
317              Vector3d velStep = vel + (dt2_ /mass * PhysicalConstants::energyConvert) * frc;
318  
319 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
319 >            Tb = sd->lab2Body(sd->getTrq());
320              Vector3d angMomStep = angMom + (dt2_ * PhysicalConstants::energyConvert) * Tb;                            
321  
322              Vector3d omegaLab;
# Line 363 | Line 335 | namespace OpenMD {
335  
336              for (int k = 0; k < maxIterNum_; k++) {
337                              
338 <              if (integrableObject->isLinear()) {
339 <                int linearAxis = integrableObject->linearAxis();
338 >              if (sd->isLinear()) {
339 >                int linearAxis = sd->linearAxis();
340                  int l = (linearAxis +1 )%3;
341                  int m = (linearAxis +2 )%3;
342                  omegaBody[l] = angMomStep[l] /I(l, l);
# Line 403 | Line 375 | namespace OpenMD {
375                  break; // iteration ends here
376              }
377  
378 <            integrableObject->addFrc(frictionForceLab);
379 <            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
378 >            sd->addFrc(frictionForceLab);
379 >            sd->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
380  
381              
382            } else {
# Line 413 | Line 385 | namespace OpenMD {
385              Vector3d randomForce;
386              Vector3d randomTorque;
387              genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
388 <            integrableObject->addFrc(randomForce);            
388 >            sd->addFrc(randomForce);            
389  
390              // What remains contains velocity explicitly, but the velocity required
391              // is at the full step: v(t + h), while we have initially the velocity
# Line 422 | Line 394 | namespace OpenMD {
394  
395              // this is the velocity at the half-step:
396              
397 <            Vector3d vel =integrableObject->getVel();
397 >            Vector3d vel =sd->getVel();
398  
399              //estimate velocity at full-step using everything but friction forces:          
400  
401 <            frc = integrableObject->getFrc();
401 >            frc = sd->getFrc();
402              Vector3d velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * frc;
403  
404              Vector3d frictionForce(0.0);
# Line 452 | Line 424 | namespace OpenMD {
424                  break; // iteration ends here
425              }
426  
427 <            integrableObject->addFrc(frictionForce);
427 >            sd->addFrc(frictionForce);
428  
429            }
430          }

Diff Legend

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