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 3450 by chuckv, Sun Sep 14 01:32:26 2008 UTC vs.
Revision 3469 by chuckv, Wed Oct 22 17:52:40 2008 UTC

# Line 58 | Line 58 | namespace oopse {
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
# Line 79 | Line 80 | namespace oopse {
80        painCave.isFatal = 1;
81        simError();
82      } else {
83 <      targetPressure_ = simParams->getTargetPressure();
83 >      targetPressure_ = simParams->getTargetPressure()/OOPSEConstant::pressureConvert;
84      }
85  
86    
# Line 92 | Line 93 | namespace oopse {
93      }
94  
95  
96 +
97 +    
98 +
99 +    //Compute initial hull
100 +    /*
101 +    surfaceMesh_->computeHull(localSites_);
102 +    Area0_ = surfaceMesh_->getArea();
103 +    int nTriangles = surfaceMesh_->getNMeshElements();
104 +    //    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
105 +    gamma_0_ = (NumericConstant::PI * targetPressure_* targetPressure_ * Area0_ * Area0_ * simParams->getDt()) /
106 +      (4.0 * nTriangles * nTriangles* OOPSEConstant::kb*simParams->getTargetTemp());
107 +    //RealType eta0 = gamma_0/
108 +    */
109 +
110      // Build the hydroProp map:
111      std::map<std::string, HydroProp*> hydroPropMap;
112  
# Line 121 | Line 136 | namespace oopse {
136          hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
137        } else {              
138          sprintf( painCave.errMsg,
139 <                 "HydroPropFile must be set to a file name if Langevin Dynamics\n"
139 >                 "HydroPropFile must be set to a file name if SMIPDynamics\n"
140                   "\tis specified for rigidBodies which contain more than one atom\n"
141 <                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n");
141 >                 "\tTo create a HydroPropFile, run the \"Hydro\" program.\n\n"
142 >                 "\tFor use with SMIPD, the default viscosity in Hydro should be\n"
143 >                 "\tset to 1.0 because the friction and random forces will be\n"
144 >                 "\tdynamically re-set assuming this is true.\n");
145          painCave.severity = OOPSE_ERROR;
146          painCave.isFatal = 1;
147          simError();  
# Line 192 | Line 210 | namespace oopse {
210                    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
211                    if (ljData != NULL) {
212                      LJParam ljParam = ljData->getData();
213 <                    currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
213 >                    currShape = new Sphere(atom->getPos(), 2.0);
214                    } else {
215                      sprintf( painCave.errMsg,
216                               "Can not cast GenericData to LJParam\n");
# Line 204 | Line 222 | namespace oopse {
222                } else {
223                  int aNum = etab.GetAtomicNum((atom->getType()).c_str());
224                  if (aNum != 0) {
225 <                  currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum));
225 >                  currShape = new Sphere(atom->getPos(), 2.0);
226                  } else {
227                    sprintf( painCave.errMsg,
228                             "Could not find atom type in default element.txt\n");
# Line 215 | Line 233 | namespace oopse {
233                }
234              }
235            }
236 <          HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
236 >          HydroProp* currHydroProp = currShape->getHydroProp(1.0,simParams->getTargetTemp());
237            std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
238            if (iter != hydroPropMap.end())
239              hydroProps_.push_back(iter->second);
# Line 230 | Line 248 | namespace oopse {
248  
249      /* Compute hull first time through to get the area of t=0*/
250  
251 <    /* Build a vector of integrable objects to determine if the are surface atoms */
251 >    //Build a vector of integrable objects to determine if the are surface atoms
252      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {          
253        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
254             integrableObject = mol->nextIntegrableObject(j)) {  
# Line 238 | Line 256 | namespace oopse {
256        }
257      }
258  
241 –    surfaceMesh_->computeHull(localSites_);
242 –    Area0_ = surfaceMesh_->getArea();
243 –
259  
260    }  
261  
# Line 277 | Line 292 | namespace oopse {
292      int fdf;
293    
294      fdf = 0;
295 <
281 <  
295 >  
296      /*Compute surface Mesh*/
297      surfaceMesh_->computeHull(localSites_);
298  
299      /* Get area and number of surface stunt doubles and compute new variance */
300 <    RealType area = surfaceMesh_->getArea();
301 <    RealType nSurfaceSDs = surfaceMesh_->getNs();
288 <    
289 <    std::cerr << "Surface Area is: " << area << " nSurfaceSDs is: " << nSurfaceSDs << std::endl;
300 >     RealType area = surfaceMesh_->getArea();
301 >     int nSurfaceSDs = surfaceMesh_->getNs();
302  
291 –    /* Compute variance for random forces */
303      
304 <    variance_ = sqrt(2.0*NumericConstant::PI)*(targetPressure_*area/nSurfaceSDs);
304 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
305 >    int nTriangles = sMesh.size();
306  
295 –    /* Loop over the mesh faces and apply random force to each of the faces*/
307  
308 <    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
309 <    std::vector<Triangle*>::iterator face;
308 >
309 >     /* Compute variance for random forces */
310 >  
311 >    RealType sigma_t = sqrt(NumericConstant::PI/2.0)*((targetPressure_)*area/nTriangles)
312 >       /OOPSEConstant::energyConvert;
313 >
314 >    gamma_t_ = (NumericConstant::PI * targetPressure_* targetPressure_ * area * area * simParams->getDt()) /
315 >      (4.0 * nTriangles * nTriangles* OOPSEConstant::kB*simParams->getTargetTemp()) /OOPSEConstant::energyConvert;
316 >
317 >    std::vector<RealType>  randNums = genTriangleForces(nTriangles, sigma_t);
318 >
319 >    /* Loop over the mesh faces and apply random force to each of the faces*/
320 >    
321 >    
322 >    std::vector<Triangle>::iterator face;
323      std::vector<StuntDouble*>::iterator vertex;
324 +    int thisNumber = 0;
325      for (face = sMesh.begin(); face != sMesh.end(); ++face){
326      
327 <      Triangle* thisTriangle = *face;
328 <      std::vector<StuntDouble*> vertexSDs = thisTriangle->getVertices();
327 >      Triangle thisTriangle = *face;
328 >      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
329 >      
330 >      /* Get Random Force */
331 >      Vector3d unitNormal = thisTriangle.getNormal();
332 >      unitNormal.normalize();
333 >      Vector3d randomForce = -randNums[thisNumber] * unitNormal;
334 >      Vector3d centroid = thisTriangle.getCentroid();
335  
336 +      Vector3d langevinForce = randomForce - gamma_t_*thisTriangle.getFacetVelocity();
337 +      
338        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
339 <         Vector3d randomForce;
340 <         Vector3d randomTorque;
341 <         genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
342 <         mass = integrableObject->getMass();
339 >        if ((*vertex) != NULL){
340 >          // mass = integrableObject->getMass();
341 >          Vector3d vertexForce = langevinForce/3.0;
342 >          (*vertex)->addFrc(vertexForce);
343  
344 <         integrableObject->addFrc(randomForce);          
344 >          if ((*vertex)->isDirectional()){
345 >            Vector3d vertexPos = (*vertex)->getPos();
346 >            Vector3d vertexCentroidVector = vertexPos - centroid;
347 >            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
348 >          }
349 >        }  
350        }
351 +    }
352  
353 <
354 <    }
355 <
356 <
357 <
358 <    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
353 >    /* Now loop over all surface particles and apply the drag*/
354 >    /*
355 >    std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
356 >    for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
357 >      integrableObject = *vertex;
358 >      mass = integrableObject->getMass();
359 >      if (integrableObject->isDirectional()){
360 >        
361 >        // preliminaries for directional objects:
362 >        
363 >        A = integrableObject->getA();
364 >        Atrans = A.transpose();
365 >        Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
366 >        //apply random force and torque at center of resistance
367 >        Mat3x3d I = integrableObject->getI();
368 >        Vector3d omegaBody;
369 >        
370 >        // What remains contains velocity explicitly, but the velocity required
371 >        // is at the full step: v(t + h), while we have initially the velocity
372 >        // at the half step: v(t + h/2).  We need to iterate to converge the
373 >        // friction force and friction torque vectors.
374 >        
375 >        // this is the velocity at the half-step:
376 >        
377 >        Vector3d vel =integrableObject->getVel();
378 >        Vector3d angMom = integrableObject->getJ();
379 >        
380 >        //estimate velocity at full-step using everything but friction forces:          
381 >        
382 >        frc = integrableObject->getFrc();
383 >        Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
384 >        
385 >        Tb = integrableObject->lab2Body(integrableObject->getTrq());
386 >        Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
387 >        
388 >        Vector3d omegaLab;
389 >        Vector3d vcdLab;
390 >        Vector3d vcdBody;
391 >        Vector3d frictionForceBody;
392 >        Vector3d frictionForceLab(0.0);
393 >        Vector3d oldFFL;  // used to test for convergence
394 >        Vector3d frictionTorqueBody(0.0);
395 >        Vector3d oldFTB;  // used to test for convergence
396 >        Vector3d frictionTorqueLab;
397 >        RealType fdot;
398 >        RealType tdot;
399  
400 <      
401 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
402 <           integrableObject = mol->nextIntegrableObject(j)) {
403 <          
404 <          mass = integrableObject->getMass();
405 <          if (integrableObject->isDirectional()){
406 <
407 <            // preliminaries for directional objects:
408 <
409 <            A = integrableObject->getA();
331 <            Atrans = A.transpose();
332 <            Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
333 <
334 <            //apply random force and torque at center of resistance
335 <
336 <            Vector3d randomForceBody;
337 <            Vector3d randomTorqueBody;
338 <            genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
339 <            Vector3d randomForceLab = Atrans * randomForceBody;
340 <            Vector3d randomTorqueLab = Atrans * randomTorqueBody;
341 <            integrableObject->addFrc(randomForceLab);            
342 <            integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));            
343 <
344 <            Mat3x3d I = integrableObject->getI();
345 <            Vector3d omegaBody;
346 <
347 <            // What remains contains velocity explicitly, but the velocity required
348 <            // is at the full step: v(t + h), while we have initially the velocity
349 <            // at the half step: v(t + h/2).  We need to iterate to converge the
350 <            // friction force and friction torque vectors.
351 <
352 <            // this is the velocity at the half-step:
400 >        //iteration starts here:
401 >        
402 >        for (int k = 0; k < maxIterNum_; k++) {
403 >          
404 >          if (integrableObject->isLinear()) {
405 >            int linearAxis = integrableObject->linearAxis();
406 >            int l = (linearAxis +1 )%3;
407 >            int m = (linearAxis +2 )%3;
408 >            omegaBody[l] = angMomStep[l] /I(l, l);
409 >            omegaBody[m] = angMomStep[m] /I(m, m);
410              
411 <            Vector3d vel =integrableObject->getVel();
412 <            Vector3d angMom = integrableObject->getJ();
413 <
414 <            //estimate velocity at full-step using everything but friction forces:          
415 <
416 <            frc = integrableObject->getFrc();
417 <            Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
418 <
419 <            Tb = integrableObject->lab2Body(integrableObject->getTrq());
420 <            Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;                            
421 <
422 <            Vector3d omegaLab;
423 <            Vector3d vcdLab;
424 <            Vector3d vcdBody;
425 <            Vector3d frictionForceBody;
426 <            Vector3d frictionForceLab(0.0);
427 <            Vector3d oldFFL;  // used to test for convergence
428 <            Vector3d frictionTorqueBody(0.0);
429 <            Vector3d oldFTB;  // used to test for convergence
430 <            Vector3d frictionTorqueLab;
374 <            RealType fdot;
375 <            RealType tdot;
376 <
377 <            //iteration starts here:
378 <
379 <            for (int k = 0; k < maxIterNum_; k++) {
380 <                            
381 <              if (integrableObject->isLinear()) {
382 <                int linearAxis = integrableObject->linearAxis();
383 <                int l = (linearAxis +1 )%3;
384 <                int m = (linearAxis +2 )%3;
385 <                omegaBody[l] = angMomStep[l] /I(l, l);
386 <                omegaBody[m] = angMomStep[m] /I(m, m);
387 <                
388 <              } else {
389 <                omegaBody[0] = angMomStep[0] /I(0, 0);
390 <                omegaBody[1] = angMomStep[1] /I(1, 1);
391 <                omegaBody[2] = angMomStep[2] /I(2, 2);
392 <              }
411 >          } else {
412 >            omegaBody[0] = angMomStep[0] /I(0, 0);
413 >            omegaBody[1] = angMomStep[1] /I(1, 1);
414 >            omegaBody[2] = angMomStep[2] /I(2, 2);
415 >          }
416 >          
417 >          omegaLab = Atrans * omegaBody;
418 >          
419 >          // apply friction force and torque at center of resistance
420 >          
421 >          vcdLab = velStep + cross(omegaLab, rcrLab);      
422 >          vcdBody = A * vcdLab;
423 >          frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
424 >          oldFFL = frictionForceLab;
425 >          frictionForceLab = Atrans * frictionForceBody;
426 >          oldFTB = frictionTorqueBody;
427 >          frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
428 >          frictionTorqueLab = Atrans * frictionTorqueBody;
429 >          
430 >          // re-estimate velocities at full-step using friction forces:
431                
432 <              omegaLab = Atrans * omegaBody;
433 <              
396 <              // apply friction force and torque at center of resistance
397 <              
398 <              vcdLab = velStep + cross(omegaLab, rcrLab);      
399 <              vcdBody = A * vcdLab;
400 <              frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody);
401 <              oldFFL = frictionForceLab;
402 <              frictionForceLab = Atrans * frictionForceBody;
403 <              oldFTB = frictionTorqueBody;
404 <              frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody);
405 <              frictionTorqueLab = Atrans * frictionTorqueBody;
406 <              
407 <              // re-estimate velocities at full-step using friction forces:
408 <              
409 <              velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
410 <              angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
432 >          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
433 >          angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
434  
435 <              // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
435 >          // check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
436                
437 <              fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
438 <              tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
439 <              
440 <              if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
441 <                break; // iteration ends here
442 <            }
437 >          fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare();
438 >          tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare();
439 >          
440 >          if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_)
441 >            break; // iteration ends here
442 >        }
443 >        
444 >        integrableObject->addFrc(frictionForceLab);
445 >        integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
446  
421 –            integrableObject->addFrc(frictionForceLab);
422 –            integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
423 –
447              
448 <          } else {
449 <            //spherical atom
448 >      } else {
449 >        //spherical atom
450  
451 <            Vector3d randomForce;
452 <            Vector3d randomTorque;
453 <            genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
454 <            integrableObject->addFrc(randomForce);            
455 <
456 <            // What remains contains velocity explicitly, but the velocity required
457 <            // is at the full step: v(t + h), while we have initially the velocity
458 <            // at the half step: v(t + h/2).  We need to iterate to converge the
459 <            // friction force vector.
460 <
461 <            // this is the velocity at the half-step:
462 <            
463 <            Vector3d vel =integrableObject->getVel();
464 <
465 <            //estimate velocity at full-step using everything but friction forces:          
466 <
467 <            frc = integrableObject->getFrc();
468 <            Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
469 <
470 <            Vector3d frictionForce(0.0);
471 <            Vector3d oldFF;  // used to test for convergence
472 <            RealType fdot;
473 <
474 <            //iteration starts here:
475 <
476 <            for (int k = 0; k < maxIterNum_; k++) {
454 <
455 <              oldFF = frictionForce;                            
456 <              frictionForce = -hydroProps_[index]->getXitt() * velStep;
457 <
458 <              // re-estimate velocities at full-step using friction forces:
459 <              
460 <              velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
461 <
462 <              // check for convergence (if the vector has converged, fdot will be 1.0):
463 <              
464 <              fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
465 <              
466 <              if (fabs(1.0 - fdot) <= forceTolerance_)
467 <                break; // iteration ends here
468 <            }
469 <
470 <            integrableObject->addFrc(frictionForce);
471 <
451 >        // What remains contains velocity explicitly, but the velocity required
452 >        // is at the full step: v(t + h), while we have initially the velocity
453 >        // at the half step: v(t + h/2).  We need to iterate to converge the
454 >        // friction force vector.
455 >        
456 >        // this is the velocity at the half-step:
457 >        
458 >        Vector3d vel =integrableObject->getVel();
459 >        
460 >        //estimate velocity at full-step using everything but friction forces:          
461 >        
462 >        frc = integrableObject->getFrc();
463 >        Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
464 >        
465 >        Vector3d frictionForce(0.0);
466 >        Vector3d oldFF;  // used to test for convergence
467 >        RealType fdot;
468 >        
469 >        //iteration starts here:
470 >        
471 >        for (int k = 0; k < maxIterNum_; k++) {
472 >          
473 >          oldFF = frictionForce;                            
474 >          frictionForce = -hydroProps_[index]->getXitt() * velStep;
475 >          //frictionForce = -gamma_t*velStep;
476 >          // re-estimate velocities at full-step using friction forces:
477            
478 <          }
478 >          velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
479 >          
480 >          // check for convergence (if the vector has converged, fdot will be 1.0):
481            
482 <        ++index;
483 <    
482 >          fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare();
483 >          
484 >          if (fabs(1.0 - fdot) <= forceTolerance_)
485 >            break; // iteration ends here
486 >        }
487 >        
488 >        integrableObject->addFrc(frictionForce);
489 >        
490 >        
491        }
492 <    }    
493 <
494 <    info_->setFdf(fdf);
495 <    // veloMunge->removeComDrift();
496 <    // Remove angular drift if we are not using periodic boundary conditions.
497 <    //if(!simParams->getUsePeriodicBoundaryConditions())
498 <    //  veloMunge->removeAngularDrift();
485 <
492 >  
493 >      
494 >  }
495 >    */
496 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
497 >    currSnapshot->setVolume(surfaceMesh_->getVolume());
498 >    
499      ForceManager::postCalculation(needStress);  
500    }
501  
502 < void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
502 >  void SMIPDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
503  
504 <
504 >    
505      Vector<RealType, 6> Z;
506      Vector<RealType, 6> generalForce;
507 <        
507 >    
508 >
509      Z[0] = randNumGen_.randNorm(0, variance);
510      Z[1] = randNumGen_.randNorm(0, variance);
511      Z[2] = randNumGen_.randNorm(0, variance);
# Line 544 | Line 558 | void SMIPDForceManager::genRandomForceAndTorque(Vector
558      return gaussRand;
559    }
560  
561 +
562 +
563 +
564 +
565   }

Diff Legend

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