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 3459 by chuckv, Tue Oct 7 17:12:48 2008 UTC vs.
Revision 3464 by chuckv, Mon Oct 20 19:36:32 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 195 | 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 207 | 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 233 | 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 241 | Line 256 | namespace oopse {
256        }
257      }
258  
259 <    surfaceMesh_->computeHull(localSites_);
245 <    Area0_ = surfaceMesh_->getArea();
246 <    //variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
247 <    
259 >
260    }  
261  
262    std::map<std::string, HydroProp*> SMIPDForceManager::parseFrictionFile(const std::string& filename) {
# Line 288 | Line 300 | namespace oopse {
300       RealType area = surfaceMesh_->getArea();
301       int nSurfaceSDs = surfaceMesh_->getNs();
302  
291 –     /* Compute variance for random forces */
303      
293 –     variance_ = sqrt(2.0*NumericConstant::PI)*((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs);
294 –    
304      std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
305 <    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
306 <    
305 >    int nTriangles = sMesh.size();
306 >
307 >
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      
# Line 310 | Line 331 | namespace oopse {
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 <
340 <         // mass = integrableObject->getMass();
341 <
342 <                    (*vertex)->addFrc(randomForce/3.0);          
339 >        if ((*vertex) != NULL){
340 >          // mass = integrableObject->getMass();
341 >          Vector3d vertexForce = langevinForce/3.0;
342 >          (*vertex)->addFrc(vertexForce);
343 >          if (integrableObject->isDirectional()){
344 >            Vector3d vertexPos = (*vertex)->getPos();
345 >            Vector3d vertexCentroidVector = vertexPos - centroid;
346 >            (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
347 >          }
348 >        }  
349        }
350      }
351  
352      /* Now loop over all surface particles and apply the drag*/
353 <
353 >    /*
354      std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
355      for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
356        integrableObject = *vertex;
# Line 332 | Line 362 | namespace oopse {
362          A = integrableObject->getA();
363          Atrans = A.transpose();
364          Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
365 <
336 <        
365 >        //apply random force and torque at center of resistance
366          Mat3x3d I = integrableObject->getI();
367          Vector3d omegaBody;
368          
# Line 416 | Line 445 | namespace oopse {
445  
446              
447        } else {
448 <        //spherical atom
448 >        //spherical atom
449  
421 –        
450          // What remains contains velocity explicitly, but the velocity required
451          // is at the full step: v(t + h), while we have initially the velocity
452          // at the half step: v(t + h/2).  We need to iterate to converge the
# Line 443 | Line 471 | namespace oopse {
471            
472            oldFF = frictionForce;                            
473            frictionForce = -hydroProps_[index]->getXitt() * velStep;
474 <          
474 >          //frictionForce = -gamma_t*velStep;
475            // re-estimate velocities at full-step using friction forces:
476            
477            velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
# Line 460 | Line 488 | namespace oopse {
488          
489          
490        }
491 +  
492        
493 <      
494 <    }
493 >  }
494 >    */
495      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
496      currSnapshot->setVolume(surfaceMesh_->getVolume());
497      
# Line 528 | Line 557 | namespace oopse {
557      return gaussRand;
558    }
559  
560 +
561 +
562 +
563 +
564   }

Diff Legend

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