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 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 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      
304 <     variance_ = sqrt(2.0*NumericConstant::PI)*((targetPressure_/OOPSEConstant::pressureConvert)*area/nSurfaceSDs);
305 <    
306 <    std::vector<Triangle*> sMesh = surfaceMesh_->getMesh();
307 <    std::vector<RealType>  randNums = genTriangleForces(sMesh.size(),variance_);
308 <    
304 >    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
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      
322 <    std::vector<Triangle*>::iterator face;
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();
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 +        if ((*vertex) != NULL){
340 +          // mass = integrableObject->getMass();
341 +          Vector3d vertexForce = langevinForce/3.0;
342 +          (*vertex)->addFrc(vertexForce);
343  
344 <         // mass = integrableObject->getMass();
345 <
346 <                    (*vertex)->addFrc(randomForce/3.0);          
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      /* Now loop over all surface particles and apply the drag*/
354 <
354 >    /*
355      std::vector<StuntDouble*> surfaceSDs = surfaceMesh_->getSurfaceAtoms();
356      for (vertex = surfaceSDs.begin(); vertex != surfaceSDs.end(); ++vertex){
357        integrableObject = *vertex;
# Line 332 | Line 363 | namespace oopse {
363          A = integrableObject->getA();
364          Atrans = A.transpose();
365          Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();  
366 <
336 <        
366 >        //apply random force and torque at center of resistance
367          Mat3x3d I = integrableObject->getI();
368          Vector3d omegaBody;
369          
# Line 416 | Line 446 | namespace oopse {
446  
447              
448        } else {
449 <        //spherical atom
449 >        //spherical atom
450  
421 –        
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
# Line 443 | Line 472 | namespace oopse {
472            
473            oldFF = frictionForce;                            
474            frictionForce = -hydroProps_[index]->getXitt() * velStep;
475 <          
475 >          //frictionForce = -gamma_t*velStep;
476            // re-estimate velocities at full-step using friction forces:
477            
478            velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
# Line 460 | Line 489 | namespace oopse {
489          
490          
491        }
492 +  
493        
494 <      
495 <    }
494 >  }
495 >    */
496      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
497      currSnapshot->setVolume(surfaceMesh_->getVolume());
498      
# Line 528 | Line 558 | namespace oopse {
558      return gaussRand;
559    }
560  
561 +
562 +
563 +
564 +
565   }

Diff Legend

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