ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/integrators/SMIPDForceManager.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/integrators/SMIPDForceManager.cpp (file contents):
Revision 3461 by chuckv, Wed Oct 15 18:26:01 2008 UTC vs.
Revision 3482 by chuckv, Fri Dec 5 16:20:39 2008 UTC

# Line 49 | Line 49
49   #include "math/ConvexHull.hpp"
50   #include "math/Triangle.hpp"
51  
52 –
52   namespace oopse {
53  
54 <  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
54 >  SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info),
55 >                                                        forceTolerance_(1e-6),
56 >                                                        maxIterNum_(4) {
57      simParams = info->getSimParams();
58      veloMunge = new Velocitizer(info);
59      
60      // Create Hull, Convex Hull for now, other options later.
61 +    
62      surfaceMesh_ = new ConvexHull();
63      
62 –    
64      /* Check that the simulation has target pressure and target
65 <       temperature set*/
66 <
65 >       temperature set */
66 >    
67      if (!simParams->haveTargetTemp()) {
68 <      sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
68 >      sprintf(painCave.errMsg,
69 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
70 >              "   without a targetTemp!\n");      
71        painCave.isFatal = 1;
72        painCave.severity = OOPSE_ERROR;
73        simError();
74      } else {
75        targetTemp_ = simParams->getTargetTemp();
76      }
77 <
77 >    
78      if (!simParams->haveTargetPressure()) {
79 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
80 <              "   without a targetPressure!\n");
81 <      
79 >      sprintf(painCave.errMsg,
80 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
81 >              "   without a targetPressure!\n");      
82        painCave.isFatal = 1;
83        simError();
84      } else {
85 <      targetPressure_ = simParams->getTargetPressure();
85 >      // Convert pressure from atm -> amu/(fs^2*Ang)
86 >      targetPressure_ = simParams->getTargetPressure() /
87 >        OOPSEConstant::pressureConvert;
88      }
84 –
89    
90      if (simParams->getUsePeriodicBoundaryConditions()) {
91 <      sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
92 <              "   with periodic boundary conditions !\n");
93 <      
91 >      sprintf(painCave.errMsg,
92 >              "SMIPDynamics error: You can't use the SMIPD integrator\n"
93 >              "   with periodic boundary conditions!\n");    
94        painCave.isFatal = 1;
95        simError();
96      }
97 +    
98 +    if (!simParams->haveViscosity()) {
99 +      sprintf(painCave.errMsg,
100 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
101 +              "   without a viscosity!\n");
102 +      painCave.isFatal = 1;
103 +      painCave.severity = OOPSE_ERROR;
104 +      simError();
105 +    }else{
106 +      viscosity_ = simParams->getViscosity();
107 +    }
108 +    
109 +    dt_ = simParams->getDt();
110 +  
111 +    variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
112  
94 –
113      // Build the hydroProp map:
114      std::map<std::string, HydroProp*> hydroPropMap;
115  
# Line 114 | Line 132 | namespace oopse {
132          
133        }
134      }
117 –        
135  
136 +    hydroProps_.resize(info->getNIntegrableObjects());
137 +    
138      if (needHydroPropFile) {              
139        if (simParams->haveHydroPropFile()) {
140          hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
# Line 123 | Line 142 | namespace oopse {
142          sprintf( painCave.errMsg,
143                   "HydroPropFile must be set to a file name if SMIPDynamics\n"
144                   "\tis specified for rigidBodies which contain more than one atom\n"
145 <                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
127 <                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
128 <                 "\tset to 1.0 because the friction and random forces will be\n"
129 <                 "\tdynamically re-set assuming this is true.\n");
145 >                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n");
146          painCave.severity = OOPSE_ERROR;
147          painCave.isFatal = 1;
148          simError();  
149        }      
150 +      
151  
152 +
153 +
154 +
155 +
156        for (mol = info->beginMolecule(i); mol != NULL;
157             mol = info->nextMolecule(i)) {
158          for (integrableObject = mol->beginIntegrableObject(j);
# Line 140 | Line 161 | namespace oopse {
161  
162            std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
163            if (iter != hydroPropMap.end()) {
164 <            hydroProps_.push_back(iter->second);
164 >            hydroProps_[integrableObject->getLocalIndex()] = iter->second;
165            } else {
166              sprintf( painCave.errMsg,
167 <                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
167 >                     "Can not find resistance tensor for atom [%s]\n",
168 >                     integrableObject->getType().c_str());
169              painCave.severity = OOPSE_ERROR;
170              painCave.isFatal = 1;
171              simError();  
# Line 218 | Line 240 | namespace oopse {
240                }
241              }
242            }
243 <          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
243 >          HydroProp* currHydroProp = currShape->getHydroProp(viscosity_, targetTemp_);
244            std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
245 <          if (iter != hydroPropMap.end())
246 <            hydroProps_.push_back(iter->second);
247 <          else {
245 >          if (iter != hydroPropMap.end()) {
246 >            hydroProps_[integrableObject->getLocalIndex()] = iter->second;
247 >          } else {
248              currHydroProp->complete();
249              hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
250 <            hydroProps_.push_back(currHydroProp);
250 >            hydroProps_[integrableObject->getLocalIndex()] = currHydroProp;
251            }
252          }
253        }
254      }
255  
256 <    /* Compute hull first time through to get the area of t=0*/
257 <
258 <    /* Build a vector of integrable objects to determine if the are surface atoms */
259 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
260 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
256 >    // Build a vector of integrable objects to determine if the are
257 >    // surface atoms
258 >    
259 >    for (mol = info_->beginMolecule(i); mol != NULL;
260 >         mol = info_->nextMolecule(i)) {          
261 >      for (integrableObject = mol->beginIntegrableObject(j);
262 >           integrableObject != NULL;
263             integrableObject = mol->nextIntegrableObject(j)) {  
264          localSites_.push_back(integrableObject);
265        }
266 <    }
243 <
244 <    surfaceMesh_->computeHull(localSites_);
245 <    Area0_ = surfaceMesh_->getArea();
246 <    //variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
247 <    
266 >    }  
267    }  
268  
269    std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
# Line 260 | Line 279 | namespace oopse {
279        HydroProp* currProp = new HydroProp(buffer);
280        props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
281      }
282 <
282 >    
283      return props;
284    }
285    
# Line 269 | Line 288 | namespace oopse {
288      Molecule::IntegrableObjectIterator  j;
289      Molecule* mol;
290      StuntDouble* integrableObject;
291 +
292      RealType mass;
293      Vector3d pos;
294      Vector3d frc;
# Line 276 | Line 296 | namespace oopse {
296      Mat3x3d Atrans;
297      Vector3d Tb;
298      Vector3d ji;
299 <    unsigned int index = 0;
280 <    int fdf;
281 <  
282 <    fdf = 0;
299 >    int index = 0;
300    
301 <    /*Compute surface Mesh*/
301 >    // Compute surface Mesh
302      surfaceMesh_->computeHull(localSites_);
303  
304 <    /* Get area and number of surface stunt doubles and compute new variance */
305 <     RealType area = surfaceMesh_->getArea();
306 <     int nSurfaceSDs = surfaceMesh_->getNs();
304 >    // Get total area and number of surface stunt doubles
305 >    RealType area = surfaceMesh_->getArea();
306 >    int nSurfaceSDs = surfaceMesh_->getNs();
307 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
308 >    int nTriangles = sMesh.size();
309  
310 <     /* Compute variance for random forces */
311 <    
312 <     variance_ = sqrt(2.0*NumericConstant::PI)*((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs)
294 <       /OOPSEConstant::energyConvert;
295 <    
296 <    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
297 <    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
298 <    
299 <    /* Loop over the mesh faces and apply random force to each of the faces*/
300 <    
301 <    
302 <    std::vector<Triangle*>::iterator face;
310 >    // Loop over the mesh faces and apply external pressure to each
311 >    // of the faces
312 >    std::vector<Triangle>::iterator face;
313      std::vector<StuntDouble*>::iterator vertex;
314      int thisNumber = 0;
315      for (face = sMesh.begin(); face != sMesh.end(); ++face){
316      
317 <      Triangle* thisTriangle = *face;
318 <      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
319 <      
320 <      /* Get Random Force */
311 <      Vector3d unitNormal = thisTriangle->getNormal();
317 >      Triangle thisTriangle = *face;
318 >      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
319 >      RealType thisArea = thisTriangle.getArea();
320 >      Vector3d unitNormal = thisTriangle.getNormal();
321        unitNormal.normalize();
322 <      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
323 <      Vector3d centroid = thisTriangle->getCentroid();
324 <
322 >      Vector3d centroid = thisTriangle.getCentroid();
323 >      Vector3d extPressure = - unitNormal * (targetPressure_ * thisArea)
324 >        / OOPSEConstant::energyConvert;
325 >      
326        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
327 <
328 <         // mass = integrableObject->getMass();
329 <        Vector3d vertexForce = randomForce/3.0;
330 <        (*vertex)->addFrc(vertexForce);
331 <        if (integrableObject->isDirectional()){
332 <          Vector3d vertexPos = (*vertex)->getPos();
333 <          Vector3d vertexCentroidVector = vertexPos - centroid;
334 <          (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
335 <        }
336 <          
327 >        if ((*vertex) != NULL){
328 >          Vector3d vertexForce = extPressure/3.0;
329 >          (*vertex)->addFrc(vertexForce);
330 >          
331 >          if ((*vertex)->isDirectional()){          
332 >            Vector3d vertexPos = (*vertex)->getPos();
333 >            Vector3d vertexCentroidVector = vertexPos - centroid;
334 >            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
335 >          }
336 >        }  
337        }
338      }
339  
340 <    /* Now loop over all surface particles and apply the drag*/
341 <
340 >    // Now loop over all surface particles and apply the drag
341 >    // and random forces
342 >    
343      std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
344      for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
345 <      integrableObject = *vertex;
345 >      integrableObject = *vertex;      
346 >      index = integrableObject->getLocalIndex();
347        mass = integrableObject->getMass();
348        if (integrableObject->isDirectional()){
349          
# Line 341 | Line 353 | namespace oopse {
353          Atrans = A.transpose();
354          Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
355  
356 <        
356 >        // apply random force and torque at center of resistance
357 >
358 >        Vector3d randomForceBody;
359 >        Vector3d randomTorqueBody;
360 >        genRandomForceAndTorque(randomForceBody, randomTorqueBody,
361 >                                index, variance_);
362 >        Vector3d randomForceLab = Atrans * randomForceBody;
363 >        Vector3d randomTorqueLab = Atrans * randomTorqueBody;
364 >        integrableObject->addFrc(randomForceLab);            
365 >        integrableObject->addTrq(randomTorqueLab + cross(rcrLab,
366 >                                                         randomForceLab ));
367 >
368          Mat3x3d I = integrableObject->getI();
369          Vector3d omegaBody;
370          
# Line 355 | Line 378 | namespace oopse {
378          Vector3d vel =integrableObject->getVel();
379          Vector3d angMom = integrableObject->getJ();
380          
381 <        //estimate velocity at full-step using everything but friction forces:          
381 >        // estimate velocity at full-step using everything but friction forces:
382          
383          frc = integrableObject->getFrc();
384          Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
# Line 424 | Line 447 | namespace oopse {
447  
448              
449        } else {
450 <        //spherical atom
450 >        //spherical atom
451  
452 +        Vector3d randomForce;
453 +        Vector3d randomTorque;
454 +        genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
455 +        integrableObject->addFrc(randomForce);            
456          
457          // What remains contains velocity explicitly, but the velocity required
458          // is at the full step: v(t + h), while we have initially the velocity
# Line 451 | Line 478 | namespace oopse {
478            
479            oldFF = frictionForce;                            
480            frictionForce = -hydroProps_[index]->getXitt() * velStep;
481 <          
481 >          //frictionForce = -gamma_t*velStep;
482            // re-estimate velocities at full-step using friction forces:
483            
484            velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
# Line 465 | Line 492 | namespace oopse {
492          }
493          
494          integrableObject->addFrc(frictionForce);
495 <        
496 <        
470 <      }
471 <      
472 <      
495 >                
496 >      }      
497      }
498 +        
499 +    veloMunge->removeComDrift();
500 +    veloMunge->removeAngularDrift();
501 +
502      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
503 <    currSnapshot->setVolume(surfaceMesh_->getVolume());
476 <    
503 >    currSnapshot->setVolume(surfaceMesh_->getVolume());    
504      ForceManager::postCalculation(needStress);  
505    }
506  
507 <  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
508 <
509 <    
507 >  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force,
508 >                                                  Vector3d& torque,
509 >                                                  unsigned int index,
510 >                                                  RealType variance) {
511 >        
512      Vector<RealType, 6> Z;
513      Vector<RealType, 6> generalForce;
514      
486 –
515      Z[0] = randNumGen_.randNorm(0, variance);
516      Z[1] = randNumGen_.randNorm(0, variance);
517      Z[2] = randNumGen_.randNorm(0, variance);
518      Z[3] = randNumGen_.randNorm(0, variance);
519      Z[4] = randNumGen_.randNorm(0, variance);
520      Z[5] = randNumGen_.randNorm(0, variance);
521 <    
521 >    
522 >
523      generalForce = hydroProps_[index]->getS()*Z;
524      
525      force[0] = generalForce[0];
# Line 498 | Line 527 | namespace oopse {
527      force[2] = generalForce[2];
528      torque[0] = generalForce[3];
529      torque[1] = generalForce[4];
530 <    torque[2] = generalForce[5];
531 <    
503 < }
504 <  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
505 <
506 <    // zero fill the random vector before starting:
507 <    std::vector<RealType> gaussRand;
508 <    gaussRand.resize(nTriangles);
509 <    std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
510 <
511 <
512 < #ifdef IS_MPI
513 <    if (worldRank == 0) {
514 < #endif
515 <      for (int i = 0; i < nTriangles; i++) {
516 <        gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
517 <      }
518 < #ifdef IS_MPI
519 <    }
520 < #endif
521 <
522 <    // push these out to the other processors
523 <
524 < #ifdef IS_MPI
525 <    if (worldRank == 0) {
526 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
527 <    } else {
528 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
529 <    }
530 < #endif
531 <
532 <    for (int i = 0; i < nTriangles; i++) {
533 <      gaussRand[i] = gaussRand[i] * variance;
534 <    }
535 <
536 <    return gaussRand;
537 <  }
538 <
530 >    torque[2] = generalForce[5];    
531 >  }
532   }

Diff Legend

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