ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/integrators/LangevinHullForceManager.cpp
(Generate patch)

Comparing branches/development/src/integrators/LangevinHullForceManager.cpp (file contents):
Revision 1668 by gezelter, Fri Jan 6 19:03:05 2012 UTC vs.
Revision 1858 by gezelter, Wed Apr 3 21:32:13 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# Line 51 | Line 51
51   #include <mpi.h>
52   #endif
53  
54 + using namespace std;
55   namespace OpenMD {
56 <
57 <  LangevinHullForceManager::LangevinHullForceManager(SimInfo* info) : ForceManager(info) {
56 >  
57 >  LangevinHullForceManager::LangevinHullForceManager(SimInfo* info) :
58 >    ForceManager(info) {
59      
60      simParams = info->getSimParams();
61      veloMunge = new Velocitizer(info);
# Line 66 | Line 68 | namespace OpenMD {
68      
69      const std::string ht = simParams->getHULL_Method();
70      
69 –    
70 –    
71      std::map<std::string, HullTypeEnum>::iterator iter;
72      iter = stringToEnumMap_.find(ht);
73 <    hullType_ = (iter == stringToEnumMap_.end()) ? LangevinHullForceManager::hullUnknown : iter->second;
74 <    if (hullType_ == hullUnknown) {
75 <      std::cerr << "WARNING! Hull Type Unknown!\n";
76 <    }
77 <    
73 >    hullType_ = (iter == stringToEnumMap_.end()) ?
74 >      LangevinHullForceManager::hullUnknown : iter->second;
75 >
76      switch(hullType_) {
77      case hullConvex :
78        surfaceMesh_ = new ConvexHull();
# Line 84 | Line 82 | namespace OpenMD {
82        break;
83      case hullUnknown :
84      default :
85 +      sprintf(painCave.errMsg,
86 +              "LangevinHallForceManager: Unknown Hull_Method was requested!\n");
87 +      painCave.isFatal = 1;
88 +      simError();      
89        break;
90      }
91 +
92 +    doThermalCoupling_ = true;
93 +    doPressureCoupling_ = true;
94 +
95      /* Check that the simulation has target pressure and target
96 <       temperature set */
91 <    
96 >       temperature set */    
97      if (!simParams->haveTargetTemp()) {
98        sprintf(painCave.errMsg,
99 <              "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
100 <              "\twithout a targetTemp (K)!\n");      
101 <      painCave.isFatal = 1;
102 <      painCave.severity = OPENMD_ERROR;
99 >              "LangevinHullForceManager: no targetTemp (K) was set.\n"
100 >              "\tOpenMD is turning off the thermal coupling to the bath.\n");
101 >      painCave.isFatal = 0;
102 >      painCave.severity = OPENMD_INFO;
103        simError();
104 +      doThermalCoupling_ = false;
105      } else {
106        targetTemp_ = simParams->getTargetTemp();
107      }
108      
109      if (!simParams->haveTargetPressure()) {
110        sprintf(painCave.errMsg,
111 <              "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
112 <              "\twithout a targetPressure (atm)!\n");      
113 <      painCave.isFatal = 1;
111 >              "LangevinHullForceManager: no targetPressure (atm) was set.\n"
112 >              "\tOpenMD is turning off the pressure coupling to the bath.\n");
113 >      painCave.isFatal = 0;
114 >      painCave.severity = OPENMD_INFO;
115        simError();
116 +      doPressureCoupling_ = false;
117      } else {
118        // Convert pressure from atm -> amu/(fs^2*Ang)
119        targetPressure_ = simParams->getTargetPressure() /
# Line 114 | Line 122 | namespace OpenMD {
122    
123      if (simParams->getUsePeriodicBoundaryConditions()) {
124        sprintf(painCave.errMsg,
125 <              "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
126 <              "\twith periodic boundary conditions!\n");    
125 >              "LangevinHallForceManager: You can't use the Langevin Hull\n"
126 >              "\tintegrator with periodic boundary conditions turned on!\n");
127        painCave.isFatal = 1;
128        simError();
129      }
130      
131      if (!simParams->haveViscosity()) {
132        sprintf(painCave.errMsg,
133 <              "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
134 <              "\twithout a viscosity!\n");
135 <      painCave.isFatal = 1;
136 <      painCave.severity = OPENMD_ERROR;
133 >              "LangevinHullForceManager: no viscosity was set.\n"
134 >              "\tOpenMD is turning off the thermal coupling to the bath.\n");
135 >      painCave.isFatal = 0;
136 >      painCave.severity = OPENMD_INFO;
137        simError();
138 +      doThermalCoupling_ = false;
139      }else{
140        viscosity_ = simParams->getViscosity();
141      }
142      
143 +    if ( fabs(viscosity_) < 1e-6 ) {
144 +      sprintf(painCave.errMsg,
145 +              "LangevinHullDynamics: The bath viscosity was set lower than\n"
146 +              "\t1e-6 poise.  OpenMD is turning off the thermal coupling to\n"
147 +              "\tthe bath.\n");
148 +      painCave.isFatal = 0;
149 +      painCave.severity = OPENMD_INFO;
150 +      simError();
151 +      doThermalCoupling_ = false;
152 +    }
153 +
154      dt_ = simParams->getDt();
155    
156 <    variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;
156 >    if (doThermalCoupling_)
157 >      variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;
158  
159      // Build a vector of integrable objects to determine if the are
160      // surface atoms
161      Molecule* mol;
162 <    StuntDouble* integrableObject;
162 >    StuntDouble* sd;
163      SimInfo::MoleculeIterator i;
164      Molecule::IntegrableObjectIterator  j;
165 <
165 >    
166      for (mol = info_->beginMolecule(i); mol != NULL;
167           mol = info_->nextMolecule(i)) {          
168 <      for (integrableObject = mol->beginIntegrableObject(j);
169 <           integrableObject != NULL;
170 <           integrableObject = mol->nextIntegrableObject(j)) {  
171 <        localSites_.push_back(integrableObject);
168 >      for (sd = mol->beginIntegrableObject(j);
169 >           sd != NULL;
170 >           sd = mol->nextIntegrableObject(j)) {
171 >        localSites_.push_back(sd);
172        }
173 <    }  
173 >    }
174 >    
175 >    // We need to make an initial guess at the bounding box in order
176 >    // to compute long range forces in ForceMatrixDecomposition:
177 >
178 >    // Compute surface Mesh
179 >    surfaceMesh_->computeHull(localSites_);
180 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
181 >    currSnapshot->setBoundingBox(surfaceMesh_->getBoundingBox());
182    }  
183 <  
183 >  
184    void LangevinHullForceManager::postCalculation(){
156 –    SimInfo::MoleculeIterator i;
157 –    Molecule::IntegrableObjectIterator  j;
158 –    Molecule* mol;
159 –    StuntDouble* integrableObject;
185    
186 +    int nTriangles, thisFacet;
187 +    RealType area, thisArea, thisMass;
188 +    vector<Triangle> sMesh;
189 +    Triangle thisTriangle;
190 +    vector<Triangle>::iterator face;
191 +    vector<StuntDouble*> vertexSDs;
192 +    vector<StuntDouble*>::iterator vertex;
193 +
194 +    Vector3d unitNormal, centroid, facetVel;
195 +    Vector3d langevinForce, vertexForce;
196 +    Vector3d extPressure, randomForce, dragForce;
197 +
198 +    Mat3x3d hydroTensor, S;
199 +    vector<Vector3d> randNums;
200 +
201      // Compute surface Mesh
202      surfaceMesh_->computeHull(localSites_);
203  
204      // Get total area and number of surface stunt doubles
205 <    RealType area = surfaceMesh_->getArea();
206 <    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
207 <    int nTriangles = sMesh.size();
205 >    area = surfaceMesh_->getArea();
206 >    sMesh = surfaceMesh_->getMesh();
207 >    nTriangles = sMesh.size();
208  
209 <    // Generate all of the necessary random forces
210 <    std::vector<Vector3d>  randNums = genTriangleForces(nTriangles, variance_);
211 <
212 <    // Loop over the mesh faces and apply external pressure to each
213 <    // of the faces
214 <    std::vector<Triangle>::iterator face;
215 <    std::vector<StuntDouble*>::iterator vertex;
176 <    int thisFacet = 0;
209 >    if (doThermalCoupling_) {
210 >      // Generate all of the necessary random forces
211 >      randNums = genTriangleForces(nTriangles, variance_);
212 >    }
213 >    
214 >    // Loop over the mesh faces
215 >    thisFacet = 0;
216      for (face = sMesh.begin(); face != sMesh.end(); ++face){
217 <      Triangle thisTriangle = *face;
218 <      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
219 <      RealType thisArea = thisTriangle.getArea();
220 <      Vector3d unitNormal = thisTriangle.getUnitNormal();
221 <      //unitNormal.normalize();
222 <      Vector3d centroid = thisTriangle.getCentroid();
223 <      Vector3d facetVel = thisTriangle.getFacetVelocity();
185 <      RealType thisMass = thisTriangle.getFacetMass();
186 <      Mat3x3d hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);
187 <      
188 <      hydroTensor *= PhysicalConstants::viscoConvert;
189 <      Mat3x3d S;
190 <      CholeskyDecomposition(hydroTensor, S);
191 <      
192 <      Vector3d extPressure = -unitNormal * (targetPressure_ * thisArea) /
193 <        PhysicalConstants::energyConvert;
217 >      thisTriangle = *face;
218 >      vertexSDs = thisTriangle.getVertices();
219 >      thisArea = thisTriangle.getArea();
220 >      unitNormal = thisTriangle.getUnitNormal();
221 >      centroid = thisTriangle.getCentroid();
222 >      facetVel = thisTriangle.getFacetVelocity();
223 >      thisMass = thisTriangle.getFacetMass();
224  
225 <      Vector3d randomForce = S * randNums[thisFacet++];
196 <      Vector3d dragForce = -hydroTensor * facetVel;
225 >      langevinForce = V3Zero;
226  
227 <      Vector3d langevinForce = (extPressure + randomForce + dragForce);
227 >      if (doPressureCoupling_) {
228 >        extPressure = -unitNormal * (targetPressure_ * thisArea) /
229 >          PhysicalConstants::energyConvert;
230 >        langevinForce += extPressure;
231 >      }
232 >
233 >      if (doThermalCoupling_) {
234 >        hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);      
235 >        hydroTensor *= PhysicalConstants::viscoConvert;
236 >        CholeskyDecomposition(hydroTensor, S);
237 >        randomForce = S * randNums[thisFacet++];
238 >        dragForce = -hydroTensor * facetVel;
239 >        langevinForce += randomForce + dragForce;
240 >      }
241        
242        // Apply triangle force to stuntdouble vertices
243        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
244          if ((*vertex) != NULL){
245 <          Vector3d vertexForce = langevinForce / RealType(3.0);
245 >          vertexForce = langevinForce / RealType(3.0);
246            (*vertex)->addFrc(vertexForce);          
247          }  
248        }
# Line 210 | Line 252 | namespace OpenMD {
252      veloMunge->removeAngularDrift();
253      
254      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
255 <    currSnapshot->setVolume(surfaceMesh_->getVolume());    
255 >    currSnapshot->setVolume(surfaceMesh_->getVolume());  
256 >    currSnapshot->setHullVolume(surfaceMesh_->getVolume());
257 >    // update the bounding box for use by ForceMatrixDecomposition:
258 >    currSnapshot->setBoundingBox(surfaceMesh_->getBoundingBox());
259      ForceManager::postCalculation();  
260    }
216 –  
217 –  
218 –  std::vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
219 –                                                             RealType variance)
220 –  {
261      
262 +  vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
263 +                                                               RealType var) {
264      // zero fill the random vector before starting:
265 <    std::vector<Vector3d> gaussRand;
265 >    vector<Vector3d> gaussRand;
266      gaussRand.resize(nTriangles);
267      std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
268      
# Line 228 | Line 270 | namespace OpenMD {
270      if (worldRank == 0) {
271   #endif
272        for (int i = 0; i < nTriangles; i++) {
273 <        gaussRand[i][0] = randNumGen_.randNorm(0.0, variance);
274 <        gaussRand[i][1] = randNumGen_.randNorm(0.0, variance);
275 <        gaussRand[i][2] = randNumGen_.randNorm(0.0, variance);
273 >        gaussRand[i][0] = randNumGen_.randNorm(0.0, var);
274 >        gaussRand[i][1] = randNumGen_.randNorm(0.0, var);
275 >        gaussRand[i][2] = randNumGen_.randNorm(0.0, var);
276        }
277   #ifdef IS_MPI
278      }

Diff Legend

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