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 3480 by chuckv, Tue Nov 25 16:32:47 2008 UTC

# Line 51 | Line 51 | namespace oopse {
51  
52  
53   namespace oopse {
54 <
54 >  
55    SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
56      simParams = info->getSimParams();
57      veloMunge = new Velocitizer(info);
58      
59      // Create Hull, Convex Hull for now, other options later.
60      surfaceMesh_ = new ConvexHull();
61 +
62      
63      
64      /* Check that the simulation has target pressure and target
65         temperature set*/
66 <
66 >    
67      if (!simParams->haveTargetTemp()) {
68        sprintf(painCave.errMsg, "You can't use the SMIPDynamics integrator without a targetTemp!\n");
69        painCave.isFatal = 1;
# Line 71 | Line 72 | namespace oopse {
72      } else {
73        targetTemp_ = simParams->getTargetTemp();
74      }
75 <
75 >    
76      if (!simParams->haveTargetPressure()) {
77        sprintf(painCave.errMsg, "SMIPDynamics error: You can't use the SMIPD integrator\n"
78                "   without a targetPressure!\n");
# Line 79 | Line 80 | namespace oopse {
80        painCave.isFatal = 1;
81        simError();
82      } else {
83 <      targetPressure_ = simParams->getTargetPressure();
83 >      /* Convert pressure from atm -> amu/(fs^2*Ang)*/
84 >      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
85      }
86  
87    
# Line 91 | Line 93 | namespace oopse {
93        simError();
94      }
95  
96 +    if (!simParams->haveViscosity()) {
97 +      sprintf(painCave.errMsg, "You can't use SMIPDynamics without a viscosity!\n");
98 +      painCave.isFatal = 1;
99 +      painCave.severity = OOPSE_ERROR;
100 +      simError();
101 +    }else{
102 +      viscosity_ = simParams->getViscosity();
103 +    }
104  
105 <    // Build the hydroProp map:
96 <    std::map<std::string, HydroProp*> hydroPropMap;
105 >    dt_ = simParams->getDt();
106  
107 +    
108 +
109 +    //Compute initial hull
110 +    /*
111 +    surfaceMesh_->computeHull(localSites_);
112 +    Area0_ = surfaceMesh_->getArea();
113 +    int nTriangles = surfaceMesh_->getNMeshElements();
114 +    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
115 +    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
116 +      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
117 +    //RealType eta0 = gamma_0/
118 +    */
119 +
120 +
121      Molecule* mol;
122      StuntDouble* integrableObject;
123      SimInfo::MoleculeIterator i;
124      Molecule::IntegrableObjectIterator  j;              
125 <    bool needHydroPropFile = false;
125 >
126      
104 –    for (mol = info->beginMolecule(i); mol != NULL;
105 –         mol = info->nextMolecule(i)) {
106 –      for (integrableObject = mol->beginIntegrableObject(j);
107 –           integrableObject != NULL;
108 –           integrableObject = mol->nextIntegrableObject(j)) {
109 –        
110 –        if (integrableObject->isRigidBody()) {
111 –          RigidBody* rb = static_cast<RigidBody*>(integrableObject);
112 –          if (rb->getNumAtoms() > 1) needHydroPropFile = true;
113 –        }
114 –        
115 –      }
116 –    }
117 –        
127  
119 –    if (needHydroPropFile) {              
120 –      if (simParams->haveHydroPropFile()) {
121 –        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
122 –      } else {              
123 –        sprintf( painCave.errMsg,
124 –                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
125 –                 "\tis specified for rigidBodies which contain more than one atom\n"
126 –                 "\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");
130 –        painCave.severity = OOPSE_ERROR;
131 –        painCave.isFatal = 1;
132 –        simError();  
133 –      }      
134 –
135 –      for (mol = info->beginMolecule(i); mol != NULL;
136 –           mol = info->nextMolecule(i)) {
137 –        for (integrableObject = mol->beginIntegrableObject(j);
138 –             integrableObject != NULL;
139 –             integrableObject = mol->nextIntegrableObject(j)) {
140 –
141 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
142 –          if (iter != hydroPropMap.end()) {
143 –            hydroProps_.push_back(iter->second);
144 –          } else {
145 –            sprintf( painCave.errMsg,
146 –                     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
147 –            painCave.severity = OOPSE_ERROR;
148 –            painCave.isFatal = 1;
149 –            simError();  
150 –          }        
151 –        }
152 –      }
153 –    } else {
154 –      
155 –      std::map<std::string, HydroProp*> hydroPropMap;
156 –      for (mol = info->beginMolecule(i); mol != NULL;
157 –           mol = info->nextMolecule(i)) {
158 –        for (integrableObject = mol->beginIntegrableObject(j);
159 –             integrableObject != NULL;
160 –             integrableObject = mol->nextIntegrableObject(j)) {
161 –          Shape* currShape = NULL;
162 –
163 –          if (integrableObject->isAtom()){
164 –            Atom* atom = static_cast<Atom*>(integrableObject);
165 –            AtomType* atomType = atom->getAtomType();
166 –            if (atomType->isGayBerne()) {
167 –              DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);              
168 –              GenericData* data = dAtomType->getPropertyByName("GayBerne");
169 –              if (data != NULL) {
170 –                GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
171 –                
172 –                if (gayBerneData != NULL) {  
173 –                  GayBerneParam gayBerneParam = gayBerneData->getData();
174 –                  currShape = new Ellipsoid(V3Zero,
175 –                                            gayBerneParam.GB_l / 2.0,
176 –                                            gayBerneParam.GB_d / 2.0,
177 –                                            Mat3x3d::identity());
178 –                } else {
179 –                  sprintf( painCave.errMsg,
180 –                           "Can not cast GenericData to GayBerneParam\n");
181 –                  painCave.severity = OOPSE_ERROR;
182 –                  painCave.isFatal = 1;
183 –                  simError();  
184 –                }
185 –              } else {
186 –                sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
187 –                painCave.severity = OOPSE_ERROR;
188 –                painCave.isFatal = 1;
189 –                simError();    
190 –              }
191 –            } else {
192 –              if (atomType->isLennardJones()){
193 –                GenericData* data = atomType->getPropertyByName("LennardJones");
194 –                if (data != NULL) {
195 –                  LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
196 –                  if (ljData != NULL) {
197 –                    LJParam ljParam = ljData->getData();
198 –                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
199 –                  } else {
200 –                    sprintf( painCave.errMsg,
201 –                             "Can not cast GenericData to LJParam\n");
202 –                    painCave.severity = OOPSE_ERROR;
203 –                    painCave.isFatal = 1;
204 –                    simError();          
205 –                  }      
206 –                }
207 –              } else {
208 –                int aNum = etab.GetAtomicNum((atom->getType()).c_str());
209 –                if (aNum != 0) {
210 –                  currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum));
211 –                } else {
212 –                  sprintf( painCave.errMsg,
213 –                           "Could not find atom type in default element.txt\n");
214 –                  painCave.severity = OOPSE_ERROR;
215 –                  painCave.isFatal = 1;
216 –                  simError();          
217 –                }
218 –              }
219 –            }
220 –          }
221 –          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
222 –          std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
223 –          if (iter != hydroPropMap.end())
224 –            hydroProps_.push_back(iter->second);
225 –          else {
226 –            currHydroProp->complete();
227 –            hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
228 –            hydroProps_.push_back(currHydroProp);
229 –          }
230 –        }
231 –      }
232 –    }
233 –
128      /* Compute hull first time through to get the area of t=0*/
129  
130 <    /* Build a vector of integrable objects to determine if the are surface atoms */
130 >    //Build a vector of integrable objects to determine if the are surface atoms
131      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
132        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
133             integrableObject = mol->nextIntegrableObject(j)) {  
# Line 241 | Line 135 | namespace oopse {
135        }
136      }
137  
138 <    surfaceMesh_->computeHull(localSites_);
245 <    Area0_ = surfaceMesh_->getArea();
246 <    //variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
247 <    
138 >
139    }  
140  
141    std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
# Line 288 | Line 179 | namespace oopse {
179       RealType area = surfaceMesh_->getArea();
180       int nSurfaceSDs = surfaceMesh_->getNs();
181  
291 –     /* Compute variance for random forces */
182      
183 <     variance_ = sqrt(2.0*NumericConstant::PI)*((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs)
184 <       /OOPSEConstant::energyConvert;
185 <    
186 <    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
187 <    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
188 <    
189 <    /* Loop over the mesh faces and apply random force to each of the faces*/
190 <    
191 <    
192 <    std::vector<Triangle*>::iterator face;
183 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
184 >    int nTriangles = sMesh.size();
185 >
186 >
187 >
188 >     /* Compute variance for random forces */
189 >    std::vector<RealType>  randNums = genTriangleForces(nTriangles, 1.0);
190 >
191 >    /* Loop over the mesh faces and apply random force to each of the faces*/    
192 >    std::vector<Triangle>::iterator face;
193      std::vector<StuntDouble*>::iterator vertex;
194      int thisNumber = 0;
195      for (face = sMesh.begin(); face != sMesh.end(); ++face){
196      
197 <      Triangle* thisTriangle = *face;
198 <      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
199 <      
200 <      /* Get Random Force */
201 <      Vector3d unitNormal = thisTriangle->getNormal();
312 <      unitNormal.normalize();
313 <      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
314 <      Vector3d centroid = thisTriangle->getCentroid();
197 >      Triangle thisTriangle = *face;
198 >      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
199 >      RealType thisArea = thisTriangle.getArea();
200 >      // RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*thisArea) /OOPSEConstant::energyConvert;
201 >      // gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * thisArea * thisArea * simParams->getDt()) /(4.0 * OOPSEConstant::kB*simParams->getTargetTemp());
202  
203 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
203 >      /* Get Random Force */
204 >      Vector3d unitNormal = thisTriangle.getNormal();
205 >      unitNormal.normalize();
206 >      //Vector3d randomForce = -randNums[thisNumber] * sigma_t * unitNormal;
207 >      Vector3d centroid = thisTriangle.getCentroid();
208 >      Vector3d facetVel = thisTriangle.getFacetVelocity();
209 >      RealType hydroLength = thisTriangle.getIncircleRadius()*2.0/NumericConstant::PI;
210  
211 <         // mass = integrableObject->getMass();
212 <        Vector3d vertexForce = randomForce/3.0;
213 <        (*vertex)->addFrc(vertexForce);
321 <        if (integrableObject->isDirectional()){
322 <          Vector3d vertexPos = (*vertex)->getPos();
323 <          Vector3d vertexCentroidVector = vertexPos - centroid;
324 <          (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
325 <        }
326 <          
327 <      }
328 <    }
211 >      RealType f_normal = viscosity_*hydroLength*OOPSEConstant::viscoConvert;
212 >      RealType extPressure = -(targetPressure_ * thisArea)/OOPSEConstant::energyConvert;
213 >      RealType randomForce = randNums[thisNumber++] * sqrt(2.0 * f_normal * OOPSEConstant::kb*targetTemp_/dt_);
214  
215 <    /* Now loop over all surface particles and apply the drag*/
215 >      RealType dragForce = -f_normal * dot(facetVel, unitNormal);
216  
332 –    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
333 –    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
334 –      integrableObject = *vertex;
335 –      mass = integrableObject->getMass();
336 –      if (integrableObject->isDirectional()){
337 –        
338 –        // preliminaries for directional objects:
339 –        
340 –        A = integrableObject->getA();
341 –        Atrans = A.transpose();
342 –        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
217  
218 <        
219 <        Mat3x3d I = integrableObject->getI();
220 <        Vector3d omegaBody;
221 <        
222 <        // What remains contains velocity explicitly, but the velocity required
349 <        // is at the full step: v(t + h), while we have initially the velocity
350 <        // at the half step: v(t + h/2).  We need to iterate to converge the
351 <        // friction force and friction torque vectors.
352 <        
353 <        // this is the velocity at the half-step:
354 <        
355 <        Vector3d vel =integrableObject->getVel();
356 <        Vector3d angMom = integrableObject->getJ();
357 <        
358 <        //estimate velocity at full-step using everything but friction forces:          
359 <        
360 <        frc = integrableObject->getFrc();
361 <        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
362 <        
363 <        Tb = integrableObject->lab2Body(integrableObject->getTrq());
364 <        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
365 <        
366 <        Vector3d omegaLab;
367 <        Vector3d vcdLab;
368 <        Vector3d vcdBody;
369 <        Vector3d frictionForceBody;
370 <        Vector3d frictionForceLab(0.0);
371 <        Vector3d oldFFL;  // used to test for convergence
372 <        Vector3d frictionTorqueBody(0.0);
373 <        Vector3d oldFTB;  // used to test for convergence
374 <        Vector3d frictionTorqueLab;
375 <        RealType fdot;
376 <        RealType tdot;
218 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
219 >      
220 >      //      Vector3d dragForce = - gamma_t_ * dot(facetVel, unitNormal) * unitNormal / OOPSEConstant::energyConvert;
221 >      
222 >      // std::cout << " " << randomForce << " " << f_normal <<   std::endl;
223  
224 <        //iteration starts here:
225 <        
226 <        for (int k = 0; k < maxIterNum_; k++) {
227 <          
228 <          if (integrableObject->isLinear()) {
229 <            int linearAxis = integrableObject->linearAxis();
384 <            int l = (linearAxis +1 )%3;
385 <            int m = (linearAxis +2 )%3;
386 <            omegaBody[l] = angMomStep[l] /I(l, l);
387 <            omegaBody[m] = angMomStep[m] /I(m, m);
388 <            
389 <          } else {
390 <            omegaBody[0] = angMomStep[0] /I(0, 0);
391 <            omegaBody[1] = angMomStep[1] /I(1, 1);
392 <            omegaBody[2] = angMomStep[2] /I(2, 2);
393 <          }
394 <          
395 <          omegaLab = Atrans * omegaBody;
396 <          
397 <          // apply friction force and torque at center of resistance
398 <          
399 <          vcdLab = velStep + cross(omegaLab, rcrLab);      
400 <          vcdBody = A * vcdLab;
401 <          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
402 <          oldFFL = frictionForceLab;
403 <          frictionForceLab = Atrans * frictionForceBody;
404 <          oldFTB = frictionTorqueBody;
405 <          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
406 <          frictionTorqueLab = Atrans * frictionTorqueBody;
407 <          
408 <          // re-estimate velocities at full-step using friction forces:
409 <              
410 <          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
411 <          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
224 >      /* Apply triangle force to stuntdouble vertices */
225 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
226 >        if ((*vertex) != NULL){
227 >          // mass = integrableObject->getMass();
228 >          Vector3d vertexForce = langevinForce/3.0;
229 >          (*vertex)->addFrc(vertexForce);
230  
231 <          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
414 <              
415 <          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
416 <          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
417 <          
418 <          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
419 <            break; // iteration ends here
420 <        }
421 <        
422 <        integrableObject->addFrc(frictionForceLab);
423 <        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
231 >          if ((*vertex)->isDirectional()){
232  
233 <            
234 <      } else {
235 <        //spherical atom
236 <
237 <        
430 <        // What remains contains velocity explicitly, but the velocity required
431 <        // is at the full step: v(t + h), while we have initially the velocity
432 <        // at the half step: v(t + h/2).  We need to iterate to converge the
433 <        // friction force vector.
434 <        
435 <        // this is the velocity at the half-step:
436 <        
437 <        Vector3d vel =integrableObject->getVel();
438 <        
439 <        //estimate velocity at full-step using everything but friction forces:          
440 <        
441 <        frc = integrableObject->getFrc();
442 <        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
443 <        
444 <        Vector3d frictionForce(0.0);
445 <        Vector3d oldFF;  // used to test for convergence
446 <        RealType fdot;
447 <        
448 <        //iteration starts here:
449 <        
450 <        for (int k = 0; k < maxIterNum_; k++) {
451 <          
452 <          oldFF = frictionForce;                            
453 <          frictionForce = -hydroProps_[index]->getXitt() * velStep;
454 <          
455 <          // re-estimate velocities at full-step using friction forces:
456 <          
457 <          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
458 <          
459 <          // check for convergence (if the vector has converged, fdot will be 1.0):
460 <          
461 <          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
462 <          
463 <          if (fabs(1.0 - fdot) <= forceTolerance_)
464 <            break; // iteration ends here
465 <        }
466 <        
467 <        integrableObject->addFrc(frictionForce);
468 <        
469 <        
233 >            Vector3d vertexPos = (*vertex)->getPos();
234 >            Vector3d vertexCentroidVector = vertexPos - centroid;
235 >            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
236 >          }
237 >        }  
238        }
239 <      
240 <      
473 <    }
239 >    }
240 >
241      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
242 <    currSnapshot->setVolume(surfaceMesh_->getVolume());
476 <    
242 >    currSnapshot->setVolume(surfaceMesh_->getVolume());    
243      ForceManager::postCalculation(needStress);  
244    }
245  
# Line 508 | Line 274 | namespace oopse {
274      gaussRand.resize(nTriangles);
275      std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
276  
277 +    
278  
279   #ifdef IS_MPI
280      if (worldRank == 0) {
281   #endif
282        for (int i = 0; i < nTriangles; i++) {
283 <        gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
283 >        //gaussRand[i] = fabs(randNumGen_.randNorm(0.0, 1.0));    
284 >        gaussRand[i] = randNumGen_.randNorm(0.0, 1.0);    
285        }
286   #ifdef IS_MPI
287      }
# Line 528 | Line 296 | namespace oopse {
296        MPI_Bcast(&gaussRand[0], nTriangles, MPI_REAL, 0, MPI_COMM_WORLD);
297      }
298   #endif
299 <
299 >  
300      for (int i = 0; i < nTriangles; i++) {
301        gaussRand[i] = gaussRand[i] * variance;
302      }
303 <
303 >  
304      return gaussRand;
305    }
306  
307 +
308   }

Diff Legend

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