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 3473 by chuckv, Fri Nov 14 15:44:34 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();
61 –
63      
63 –    
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()/OOPSEConstant::pressureConvert;
85 >      // Convert pressure from atm -> amu/(fs^2*Ang)
86 >      targetPressure_ = simParams->getTargetPressure() /
87 >        OOPSEConstant::pressureConvert;
88      }
85 –
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 <
97 >    
98      if (!simParams->haveViscosity()) {
99 <      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
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      }
101 –
102 –
108      
109 +    dt_ = simParams->getDt();
110 +  
111 +    variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
112  
105 –    //Compute initial hull
106 –    /*
107 –    surfaceMesh_->computeHull(localSites_);
108 –    Area0_ = surfaceMesh_->getArea();
109 –    int nTriangles = surfaceMesh_->getNMeshElements();
110 –    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
111 –    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
112 –      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
113 –    //RealType eta0 = gamma_0/
114 –    */
115 –
113      // Build the hydroProp map:
114      std::map<std::string, HydroProp*> hydroPropMap;
115  
# Line 135 | Line 132 | namespace oopse {
132          
133        }
134      }
138 –        
135  
136 +    hydroProps_.resize(info->getNIntegrableObjects());
137 +    
138      if (needHydroPropFile) {              
139        if (simParams->haveHydroPropFile()) {
140          hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
# Line 144 | 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"
148 <                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
149 <                 "\tset to 1.0 because the friction and random forces will be\n"
150 <                 "\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 161 | 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 216 | Line 217 | namespace oopse {
217                    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
218                    if (ljData != NULL) {
219                      LJParam ljParam = ljData->getData();
220 <                    currShape = new Sphere(atom->getPos(), 2.0);
220 >                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
221                    } else {
222                      sprintf( painCave.errMsg,
223                               "Can not cast GenericData to LJParam\n");
# Line 228 | Line 229 | namespace oopse {
229                } else {
230                  int aNum = etab.GetAtomicNum((atom->getType()).c_str());
231                  if (aNum != 0) {
232 <                  currShape = new Sphere(atom->getPos(), 2.0);
232 >                  currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum));
233                  } else {
234                    sprintf( painCave.errMsg,
235                             "Could not find atom type in default element.txt\n");
# Line 239 | 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 <    }
264 <
265 <
266 >    }  
267    }  
268  
269    std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
# Line 278 | 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 287 | 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 294 | Line 296 | namespace oopse {
296      Mat3x3d Atrans;
297      Vector3d Tb;
298      Vector3d ji;
299 <    unsigned int index = 0;
298 <    int fdf;
299 <  
300 <    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();
308 <
309 <    
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 <
311 <
315 <     /* Compute variance for random forces */
316 <  
317 <
318 <
319 <
320 <    std::vector<RealType>  randNums = genTriangleForces(nTriangles, 1.0);
321 <
322 <    /* Loop over the mesh faces and apply random force to each of the faces*/
323 <    
324 <    
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;
# Line 330 | Line 317 | namespace oopse {
317        Triangle thisTriangle = *face;
318        std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
319        RealType thisArea = thisTriangle.getArea();
333 –      // RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*thisArea) /OOPSEConstant::energyConvert;
334 –      // gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());
335 –
336 –      /* Get Random Force */
320        Vector3d unitNormal = thisTriangle.getNormal();
321        unitNormal.normalize();
339 –      //Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
322        Vector3d centroid = thisTriangle.getCentroid();
323 <      Vector3d facetVel = thisTriangle.getFacetVelocity();
324 <      RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/3.14;
343 <
344 <      RealType f_normal = simParams->getViscosity()*hydroLength*1.439326479e4*OOPSEConstant::energyConvert;
345 <      RealType extPressure = -(targetPressure_ * thisArea);
346 <      RealType randomForce = randNums[thisNumber] * f_normal * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
347 <      RealType dragForce = -f_normal * dot(facetVel, unitNormal);      
348 <      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
323 >      Vector3d extPressure = - unitNormal * (targetPressure_ * thisArea)
324 >        / OOPSEConstant::energyConvert;
325        
350 –      //      Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
351 –      
352 –      //std::cout << "randomForce " << randomForce << " dragForce " << dragForce <<  " hydro  " << hydroLength << std::endl;
353 –
354 –
326        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
327 <        if ((*vertex) != NULL){
328 <          // mass = integrableObject->getMass();
358 <          Vector3d vertexForce = langevinForce/3.0;
327 >        if ((*vertex) != NULL){
328 >          Vector3d vertexForce = extPressure/3.0;
329            (*vertex)->addFrc(vertexForce);
330 <
331 <          if ((*vertex)->isDirectional()){
362 <
330 >          
331 >          if ((*vertex)->isDirectional()){          
332              Vector3d vertexPos = (*vertex)->getPos();
333              Vector3d vertexCentroidVector = vertexPos - centroid;
334              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 368 | Line 337 | namespace oopse {
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 381 | Line 352 | namespace oopse {
352          A = integrableObject->getA();
353          Atrans = A.transpose();
354          Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
355 <        //apply random force and torque at center of resistance
355 >
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 395 | 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 466 | Line 449 | namespace oopse {
449        } else {
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
459          // at the half step: v(t + h/2).  We need to iterate to converge the
# Line 504 | Line 492 | namespace oopse {
492          }
493          
494          integrableObject->addFrc(frictionForce);
495 <        
495 >                
496 >      }      
497 >    }
498          
499 <      }
500 <  
501 <      
512 <  }
513 <    */
499 >    veloMunge->removeComDrift();
500 >    veloMunge->removeAngularDrift();
501 >
502      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
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      
525 –
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 537 | 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 <    
542 < }
543 <  std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles, RealType variance) {
544 <
545 <    // zero fill the random vector before starting:
546 <    std::vector<RealType> gaussRand;
547 <    gaussRand.resize(nTriangles);
548 <    std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
549 <
550 <
551 < #ifdef IS_MPI
552 <    if (worldRank == 0) {
553 < #endif
554 <      for (int i = 0; i < nTriangles; i++) {
555 <        //gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
556 <        gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);    
557 <      }
558 < #ifdef IS_MPI
559 <    }
560 < #endif
561 <
562 <    // push these out to the other processors
563 <
564 < #ifdef IS_MPI
565 <    if (worldRank == 0) {
566 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
567 <    } else {
568 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
569 <    }
570 < #endif
571 <  
572 <    for (int i = 0; i < nTriangles; i++) {
573 <      gaussRand[i] = gaussRand[i] * variance;
574 <    }
575 <  
576 <    return gaussRand;
577 <  }
578 <
579 <
580 <
581 <
582 <
530 >    torque[2] = generalForce[5];    
531 >  }
532   }

Diff Legend

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