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 1874 by gezelter, Wed May 15 15:09:35 2013 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  
# 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 +          delete currShape;
212          }
213        }
214      }
# Line 263 | Line 236 | namespace OpenMD {
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 300 | Line 273 | namespace OpenMD {
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 323 | Line 296 | namespace OpenMD {
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 336 | Line 309 | namespace OpenMD {
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();
317 >            frc = sd->getFrc();
318              Vector3d velStep = vel + (dt2_ /mass * PhysicalConstants::energyConvert) * frc;
319  
320 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
320 >            Tb = sd->lab2Body(sd->getTrq());
321              Vector3d angMomStep = angMom + (dt2_ * PhysicalConstants::energyConvert) * Tb;                            
322  
323              Vector3d omegaLab;
# Line 363 | Line 336 | namespace OpenMD {
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 403 | Line 376 | namespace OpenMD {
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 413 | Line 386 | namespace OpenMD {
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 422 | Line 395 | namespace OpenMD {
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();
402 >            frc = sd->getFrc();
403              Vector3d velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * frc;
404  
405              Vector3d frictionForce(0.0);
# Line 452 | Line 425 | namespace OpenMD {
425                  break; // iteration ends here
426              }
427  
428 <            integrableObject->addFrc(frictionForce);
428 >            sd->addFrc(frictionForce);
429  
430            }
431          }

Diff Legend

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