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 1876 by gezelter, Fri May 17 17:10:11 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) {
58 <    
56 >  
57 >  LangevinHullForceManager::LangevinHullForceManager(SimInfo* info) :
58 >    ForceManager(info) {
59 >  
60      simParams = info->getSimParams();
61      veloMunge = new Velocitizer(info);
62      
# 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 +      if (!simParams->haveViscosity()) {
109 +        sprintf(painCave.errMsg,
110 +                "LangevinHullForceManager: no viscosity was set.\n"
111 +                "\tOpenMD is turning off the thermal coupling to the bath.\n");
112 +        painCave.isFatal = 0;
113 +        painCave.severity = OPENMD_INFO;
114 +        simError();      
115 +        doThermalCoupling_ = false;
116 +      }else{
117 +        viscosity_ = simParams->getViscosity();
118 +        if ( fabs(viscosity_) < 1e-6 ) {
119 +          sprintf(painCave.errMsg,
120 +                  "LangevinHullDynamics: The bath viscosity was set lower\n"
121 +                  "\tthan 1e-6 poise.  OpenMD is turning off the thermal\n"
122 +                  "\tcoupling to the bath.\n");
123 +          painCave.isFatal = 0;
124 +          painCave.severity = OPENMD_INFO;
125 +          simError();
126 +          doThermalCoupling_ = false;
127 +        }      
128 +      }      
129      }
102 –    
130      if (!simParams->haveTargetPressure()) {
131        sprintf(painCave.errMsg,
132 <              "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
133 <              "\twithout a targetPressure (atm)!\n");      
134 <      painCave.isFatal = 1;
132 >              "LangevinHullForceManager: no targetPressure (atm) was set.\n"
133 >              "\tOpenMD is turning off the pressure coupling to the bath.\n");
134 >      painCave.isFatal = 0;
135 >      painCave.severity = OPENMD_INFO;
136        simError();
137 +      doPressureCoupling_ = false;
138      } else {
139        // Convert pressure from atm -> amu/(fs^2*Ang)
140        targetPressure_ = simParams->getTargetPressure() /
141          PhysicalConstants::pressureConvert;
142      }
114 –  
143      if (simParams->getUsePeriodicBoundaryConditions()) {
144        sprintf(painCave.errMsg,
145 <              "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
146 <              "\twith periodic boundary conditions!\n");    
145 >              "LangevinHallForceManager: You can't use the Langevin Hull\n"
146 >              "\tintegrator with periodic boundary conditions turned on!\n");
147        painCave.isFatal = 1;
148        simError();
149      }
150 <    
123 <    if (!simParams->haveViscosity()) {
124 <      sprintf(painCave.errMsg,
125 <              "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
126 <              "\twithout a viscosity!\n");
127 <      painCave.isFatal = 1;
128 <      painCave.severity = OPENMD_ERROR;
129 <      simError();
130 <    }else{
131 <      viscosity_ = simParams->getViscosity();
132 <    }
133 <    
150 >
151      dt_ = simParams->getDt();
135 –  
136 –    variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;
152  
153 +    if (doThermalCoupling_)
154 +      variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;
155 +
156      // Build a vector of integrable objects to determine if the are
157      // surface atoms
158      Molecule* mol;
159 <    StuntDouble* integrableObject;
159 >    StuntDouble* sd;
160      SimInfo::MoleculeIterator i;
161      Molecule::IntegrableObjectIterator  j;
162 <
162 >    
163      for (mol = info_->beginMolecule(i); mol != NULL;
164           mol = info_->nextMolecule(i)) {          
165 <      for (integrableObject = mol->beginIntegrableObject(j);
166 <           integrableObject != NULL;
167 <           integrableObject = mol->nextIntegrableObject(j)) {  
168 <        localSites_.push_back(integrableObject);
165 >      for (sd = mol->beginIntegrableObject(j);
166 >           sd != NULL;
167 >           sd = mol->nextIntegrableObject(j)) {
168 >        localSites_.push_back(sd);
169        }
170 <    }  
170 >    }
171 >    
172 >    // We need to make an initial guess at the bounding box in order
173 >    // to compute long range forces in ForceMatrixDecomposition:
174 >
175 >    // Compute surface Mesh
176 >    surfaceMesh_->computeHull(localSites_);
177    }  
178 <  
178 >
179 >  LangevinHullForceManager::~LangevinHullForceManager() {
180 >    delete surfaceMesh_;
181 >    delete veloMunge;
182 >  }
183 >  
184    void LangevinHullForceManager::postCalculation(){
156 –    SimInfo::MoleculeIterator i;
157 –    Molecule::IntegrableObjectIterator  j;
158 –    Molecule* mol;
159 –    StuntDouble* integrableObject;
185    
186 <    // Compute surface Mesh
187 <    surfaceMesh_->computeHull(localSites_);
186 >    int nTriangles, thisFacet;
187 >    RealType thisArea;
188 >    vector<Triangle> sMesh;
189 >    Triangle thisTriangle;
190 >    vector<Triangle>::iterator face;
191 >    vector<StuntDouble*> vertexSDs;
192 >    vector<StuntDouble*>::iterator vertex;
193  
194 <    // Get total area and number of surface stunt doubles
195 <    RealType area = surfaceMesh_->getArea();
196 <    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
167 <    int nTriangles = sMesh.size();
194 >    Vector3d unitNormal, centroid, facetVel;
195 >    Vector3d langevinForce, vertexForce;
196 >    Vector3d extPressure, randomForce, dragForce;
197  
198 <    // Generate all of the necessary random forces
199 <    std::vector<Vector3d>  randNums = genTriangleForces(nTriangles, variance_);
198 >    Mat3x3d hydroTensor, S;
199 >    vector<Vector3d> randNums;
200  
201 <    // Loop over the mesh faces and apply external pressure to each
202 <    // of the faces
203 <    std::vector<Triangle>::iterator face;
204 <    std::vector<StuntDouble*>::iterator vertex;
205 <    int thisFacet = 0;
201 >    // Compute surface Mesh
202 >    surfaceMesh_->computeHull(localSites_);
203 >    // Get number of surface stunt doubles
204 >    sMesh = surfaceMesh_->getMesh();
205 >    nTriangles = sMesh.size();
206 >
207 >    if (doThermalCoupling_) {
208 >      // Generate all of the necessary random forces
209 >      randNums = genTriangleForces(nTriangles, variance_);
210 >    }
211 >    
212 >    // Loop over the mesh faces
213 >    thisFacet = 0;
214      for (face = sMesh.begin(); face != sMesh.end(); ++face){
215 <      Triangle thisTriangle = *face;
216 <      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
217 <      RealType thisArea = thisTriangle.getArea();
218 <      Vector3d unitNormal = thisTriangle.getUnitNormal();
219 <      //unitNormal.normalize();
220 <      Vector3d centroid = thisTriangle.getCentroid();
184 <      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;
215 >      thisTriangle = *face;
216 >      vertexSDs = thisTriangle.getVertices();
217 >      thisArea = thisTriangle.getArea();
218 >      unitNormal = thisTriangle.getUnitNormal();
219 >      centroid = thisTriangle.getCentroid();
220 >      facetVel = thisTriangle.getFacetVelocity();
221  
222 <      Vector3d randomForce = S * randNums[thisFacet++];
196 <      Vector3d dragForce = -hydroTensor * facetVel;
222 >      langevinForce = V3Zero;
223  
224 <      Vector3d langevinForce = (extPressure + randomForce + dragForce);
224 >      if (doPressureCoupling_) {
225 >        extPressure = -unitNormal * (targetPressure_ * thisArea) /
226 >          PhysicalConstants::energyConvert;
227 >        langevinForce += extPressure;
228 >      }
229 >
230 >      if (doThermalCoupling_) {
231 >        hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);      
232 >        hydroTensor *= PhysicalConstants::viscoConvert;
233 >        CholeskyDecomposition(hydroTensor, S);
234 >        randomForce = S * randNums[thisFacet++];
235 >        dragForce = -hydroTensor * facetVel;
236 >        langevinForce += randomForce + dragForce;
237 >      }
238        
239        // Apply triangle force to stuntdouble vertices
240        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
241          if ((*vertex) != NULL){
242 <          Vector3d vertexForce = langevinForce / RealType(3.0);
242 >          vertexForce = langevinForce / RealType(3.0);
243            (*vertex)->addFrc(vertexForce);          
244          }  
245        }
# Line 210 | Line 249 | namespace OpenMD {
249      veloMunge->removeAngularDrift();
250      
251      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
252 <    currSnapshot->setVolume(surfaceMesh_->getVolume());    
252 >    currSnapshot->setVolume(surfaceMesh_->getVolume());  
253 >    currSnapshot->setHullVolume(surfaceMesh_->getVolume());
254      ForceManager::postCalculation();  
255    }
216 –  
217 –  
218 –  std::vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
219 –                                                             RealType variance)
220 –  {
256      
257 +  vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
258 +                                                               RealType var) {
259      // zero fill the random vector before starting:
260 <    std::vector<Vector3d> gaussRand;
260 >    vector<Vector3d> gaussRand;
261      gaussRand.resize(nTriangles);
262      std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
263      
# Line 228 | Line 265 | namespace OpenMD {
265      if (worldRank == 0) {
266   #endif
267        for (int i = 0; i < nTriangles; i++) {
268 <        gaussRand[i][0] = randNumGen_.randNorm(0.0, variance);
269 <        gaussRand[i][1] = randNumGen_.randNorm(0.0, variance);
270 <        gaussRand[i][2] = randNumGen_.randNorm(0.0, variance);
268 >        gaussRand[i][0] = randNumGen_.randNorm(0.0, var);
269 >        gaussRand[i][1] = randNumGen_.randNorm(0.0, var);
270 >        gaussRand[i][2] = randNumGen_.randNorm(0.0, var);
271        }
272   #ifdef IS_MPI
273      }
# Line 239 | Line 276 | namespace OpenMD {
276      // push these out to the other processors
277      
278   #ifdef IS_MPI
279 <    if (worldRank == 0) {
280 <      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
244 <    } else {
245 <      MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
246 <    }
279 >    // Same command on all nodes:
280 >    MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
281   #endif
282      
283      return gaussRand;

Diff Legend

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