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 3481 by gezelter, Wed Nov 26 14:26:17 2008 UTC vs.
Revision 3516 by gezelter, Wed Jul 22 15:00:21 2009 UTC

# Line 46 | Line 46 | namespace oopse {
46   #include "math/Triangle.hpp"
47  
48   namespace oopse {
49 <  
49 >
50    SMIPDForceManager::SMIPDForceManager(SimInfo* info) : ForceManager(info) {
51  
52      simParams = info->getSimParams();
53 +    thermo = new Thermo(info);
54 +    veloMunge = new Velocitizer(info);
55      
56      // Create Hull, Convex Hull for now, other options later.
57      
58      surfaceMesh_ = new ConvexHull();
59      
60      /* Check that the simulation has target pressure and target
61 <       temperature set*/
61 >       temperature set */
62      
63      if (!simParams->haveTargetTemp()) {
64        sprintf(painCave.errMsg,
65                "SMIPDynamics error: You can't use the SMIPD integrator\n"
66 <              "   without a targetTemp!\n");      
66 >              "   without a targetTemp (K)!\n");      
67        painCave.isFatal = 1;
68        painCave.severity = OOPSE_ERROR;
69        simError();
# Line 72 | Line 74 | namespace oopse {
74      if (!simParams->haveTargetPressure()) {
75        sprintf(painCave.errMsg,
76                "SMIPDynamics error: You can't use the SMIPD integrator\n"
77 <              "   without a targetPressure!\n");      
77 >              "   without a targetPressure (atm)!\n");      
78        painCave.isFatal = 1;
79        simError();
80      } else {
# Line 84 | Line 86 | namespace oopse {
86      if (simParams->getUsePeriodicBoundaryConditions()) {
87        sprintf(painCave.errMsg,
88                "SMIPDynamics error: You can't use the SMIPD integrator\n"
89 <              "   with periodic boundary conditions!\n");    
89 >              "   with periodic boundary conditions!\n");    
90        painCave.isFatal = 1;
91        simError();
92      }
93      
94 <    if (!simParams->haveViscosity()) {
94 >    if (!simParams->haveThermalConductivity()) {
95        sprintf(painCave.errMsg,
96                "SMIPDynamics error: You can't use the SMIPD integrator\n"
97 <              "   without a viscosity!\n");
97 >              "   without a thermalConductivity (W m^-1 K^-1)!\n");
98        painCave.isFatal = 1;
99        painCave.severity = OOPSE_ERROR;
100        simError();
101      }else{
102 <      viscosity_ = simParams->getViscosity();
102 >      thermalConductivity_ = simParams->getThermalConductivity() *
103 >        OOPSEConstant::thermalConductivityConvert;
104      }
105 +
106 +    if (!simParams->haveThermalLength()) {
107 +      sprintf(painCave.errMsg,
108 +              "SMIPDynamics error: You can't use the SMIPD integrator\n"
109 +              "   without a thermalLength (Angstroms)!\n");
110 +      painCave.isFatal = 1;
111 +      painCave.severity = OOPSE_ERROR;
112 +      simError();
113 +    }else{
114 +      thermalLength_ = simParams->getThermalLength();
115 +    }
116      
117      dt_ = simParams->getDt();
118    
119      variance_ = 2.0 * OOPSEConstant::kb * targetTemp_ / dt_;
120  
121 +    // Build a vector of integrable objects to determine if the are
122 +    // surface atoms
123      Molecule* mol;
124      StuntDouble* integrableObject;
125      SimInfo::MoleculeIterator i;
126 <    Molecule::IntegrableObjectIterator  j;              
126 >    Molecule::IntegrableObjectIterator  j;
127  
112 –    // Build a vector of integrable objects to determine if the are
113 –    // surface atoms
114 –
128      for (mol = info_->beginMolecule(i); mol != NULL;
129           mol = info_->nextMolecule(i)) {          
130        for (integrableObject = mol->beginIntegrableObject(j);
# Line 121 | Line 134 | namespace oopse {
134        }
135      }  
136    }  
137 <  
137 >  
138    void SMIPDForceManager::postCalculation(bool needStress){
139      SimInfo::MoleculeIterator i;
140      Molecule::IntegrableObjectIterator  j;
141      Molecule* mol;
142      StuntDouble* integrableObject;
143 <  
143 >  
144      // Compute surface Mesh
145      surfaceMesh_->computeHull(localSites_);
146 <    
146 >
147      // Get total area and number of surface stunt doubles
148      RealType area = surfaceMesh_->getArea();
149 <    int nSurfaceSDs = surfaceMesh_->getNs();        
149 >    int nSurfaceSDs = surfaceMesh_->getNs();
150      std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
151      int nTriangles = sMesh.size();
152 <        
152 >
153      // Generate all of the necessary random forces
154      std::vector<RealType>  randNums = genTriangleForces(nTriangles, variance_);
155  
156 <    // Loop over the mesh faces and apply random force to each of the faces
156 >    RealType instaTemp = thermo->getTemperature();
157 >
158 >    // Loop over the mesh faces and apply external pressure to each
159 >    // of the faces
160      std::vector<Triangle>::iterator face;
161      std::vector<StuntDouble*>::iterator vertex;
162 <    int thisNumber = 0;
162 >    int thisFacet = 0;
163      for (face = sMesh.begin(); face != sMesh.end(); ++face){
164      
165        Triangle thisTriangle = *face;
# Line 151 | Line 167 | namespace oopse {
167        RealType thisArea = thisTriangle.getArea();
168        Vector3d unitNormal = thisTriangle.getNormal();
169        unitNormal.normalize();
154 –
170        Vector3d centroid = thisTriangle.getCentroid();
171        Vector3d facetVel = thisTriangle.getFacetVelocity();
172 <      RealType hydroLength = thisTriangle.getIncircleRadius() * 2.0 /
158 <        NumericConstant::PI;
159 <      
160 <      // gamma is the drag coefficient normal to the face of the triangle      
161 <      RealType gamma = viscosity_ * hydroLength *
162 <        OOPSEConstant::viscoConvert;
172 >      RealType thisMass = thisTriangle.getFacetMass();
173  
174 +      // gamma is the drag coefficient normal to the face of the triangle      
175 +      RealType gamma = thermalConductivity_ * thisMass * thisArea  
176 +        / (2.0 * thermalLength_ * OOPSEConstant::kB);
177 +      
178 +      gamma *= fabs(1.0 - targetTemp_/instaTemp);      
179 +      
180        RealType extPressure = - (targetPressure_ * thisArea) /
181          OOPSEConstant::energyConvert;
182  
183 <      RealType randomForce = randNums[thisNumber++] * sqrt(gamma);
183 >      RealType randomForce = randNums[thisFacet++] * sqrt(gamma);
184        RealType dragForce = -gamma * dot(facetVel, unitNormal);
185  
186 <      Vector3d langevinForce = (extPressure + randomForce + dragForce) * unitNormal;
186 >      Vector3d langevinForce = (extPressure + randomForce + dragForce) *
187 >        unitNormal;
188        
189        // Apply triangle force to stuntdouble vertices
190 <      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex) {
191 <        if ((*vertex) != NULL) {
190 >      for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
191 >        if ((*vertex) != NULL){
192            Vector3d vertexForce = langevinForce / 3.0;
193 <          (*vertex)->addFrc(vertexForce);          
194 <          if ((*vertex)->isDirectional()) {
193 >          (*vertex)->addFrc(vertexForce);          
194 >          if ((*vertex)->isDirectional()){          
195              Vector3d vertexPos = (*vertex)->getPos();
196              Vector3d vertexCentroidVector = vertexPos - centroid;
197              (*vertex)->addTrq(cross(vertexCentroidVector,vertexForce));
# Line 183 | Line 200 | namespace oopse {
200        }
201      }
202      
203 +    veloMunge->removeComDrift();
204 +    veloMunge->removeAngularDrift();
205 +    
206      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
207      currSnapshot->setVolume(surfaceMesh_->getVolume());    
208      ForceManager::postCalculation(needStress);  
209    }
210 <
210 >  
211 >  
212    std::vector<RealType> SMIPDForceManager::genTriangleForces(int nTriangles,
213 <                                                             RealType variance) {
214 <
213 >                                                             RealType variance)
214 >  {
215 >    
216      // zero fill the random vector before starting:
217      std::vector<RealType> gaussRand;
218      gaussRand.resize(nTriangles);
219      std::fill(gaussRand.begin(), gaussRand.end(), 0.0);
220 <  
220 >    
221   #ifdef IS_MPI
222      if (worldRank == 0) {
223   #endif
224        for (int i = 0; i < nTriangles; i++) {
225 <        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
225 >        gaussRand[i] = randNumGen_.randNorm(0.0, variance);
226        }
227   #ifdef IS_MPI
228      }
229   #endif
230 <
230 >    
231      // push these out to the other processors
232 <
232 >    
233   #ifdef IS_MPI
234      if (worldRank == 0) {
235 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
235 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
236      } else {
237 <      MPI_Bcast(&gaussRand[0], nTriangles, MPI_REALTYPE, 0, MPI_COMM_WORLD);
237 >      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles, MPI::REALTYPE, 0);
238      }
239   #endif
240 <  
240 >    
241      return gaussRand;
242    }
243   }

Diff Legend

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