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 1874 by gezelter, Wed May 15 15:09:35 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 <  
179 <  void LangevinHullForceManager::postCalculation(){
180 <    SimInfo::MoleculeIterator i;
181 <    Molecule::IntegrableObjectIterator  j;
182 <    Molecule* mol;
159 <    StuntDouble* integrableObject;
178 >
179 >  LangevinHullForceManager::~LangevinHullForceManager() {
180 >    delete surfaceMesh_;
181 >    delete veloMunge;
182 >  }
183    
184 +  void LangevinHullForceManager::postCalculation(){
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_);
163 –
203      // Get total area and number of surface stunt doubles
204 <    RealType area = surfaceMesh_->getArea();
205 <    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
206 <    int nTriangles = sMesh.size();
204 >    area = surfaceMesh_->getArea();
205 >    sMesh = surfaceMesh_->getMesh();
206 >    nTriangles = sMesh.size();
207  
208 <    // Generate all of the necessary random forces
209 <    std::vector<Vector3d>  randNums = genTriangleForces(nTriangles, variance_);
210 <
211 <    // Loop over the mesh faces and apply external pressure to each
212 <    // of the faces
213 <    std::vector<Triangle>::iterator face;
214 <    std::vector<StuntDouble*>::iterator vertex;
176 <    int thisFacet = 0;
208 >    if (doThermalCoupling_) {
209 >      // Generate all of the necessary random forces
210 >      randNums = genTriangleForces(nTriangles, variance_);
211 >    }
212 >    
213 >    // Loop over the mesh faces
214 >    thisFacet = 0;
215      for (face = sMesh.begin(); face != sMesh.end(); ++face){
216 <      Triangle thisTriangle = *face;
217 <      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
218 <      RealType thisArea = thisTriangle.getArea();
219 <      Vector3d unitNormal = thisTriangle.getUnitNormal();
220 <      //unitNormal.normalize();
221 <      Vector3d centroid = thisTriangle.getCentroid();
222 <      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;
216 >      thisTriangle = *face;
217 >      vertexSDs = thisTriangle.getVertices();
218 >      thisArea = thisTriangle.getArea();
219 >      unitNormal = thisTriangle.getUnitNormal();
220 >      centroid = thisTriangle.getCentroid();
221 >      facetVel = thisTriangle.getFacetVelocity();
222 >      thisMass = thisTriangle.getFacetMass();
223  
224 <      Vector3d randomForce = S * randNums[thisFacet++];
196 <      Vector3d dragForce = -hydroTensor * facetVel;
224 >      langevinForce = V3Zero;
225  
226 <      Vector3d langevinForce = (extPressure + randomForce + dragForce);
226 >      if (doPressureCoupling_) {
227 >        extPressure = -unitNormal * (targetPressure_ * thisArea) /
228 >          PhysicalConstants::energyConvert;
229 >        langevinForce += extPressure;
230 >      }
231 >
232 >      if (doThermalCoupling_) {
233 >        hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);      
234 >        hydroTensor *= PhysicalConstants::viscoConvert;
235 >        CholeskyDecomposition(hydroTensor, S);
236 >        randomForce = S * randNums[thisFacet++];
237 >        dragForce = -hydroTensor * facetVel;
238 >        langevinForce += randomForce + dragForce;
239 >      }
240        
241        // Apply triangle force to stuntdouble vertices
242        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
243          if ((*vertex) != NULL){
244 <          Vector3d vertexForce = langevinForce / RealType(3.0);
244 >          vertexForce = langevinForce / RealType(3.0);
245            (*vertex)->addFrc(vertexForce);          
246          }  
247        }
# Line 210 | Line 251 | namespace OpenMD {
251      veloMunge->removeAngularDrift();
252      
253      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
254 <    currSnapshot->setVolume(surfaceMesh_->getVolume());    
254 >    currSnapshot->setVolume(surfaceMesh_->getVolume());  
255 >    currSnapshot->setHullVolume(surfaceMesh_->getVolume());
256      ForceManager::postCalculation();  
257    }
216 –  
217 –  
218 –  std::vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
219 –                                                             RealType variance)
220 –  {
258      
259 +  vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
260 +                                                               RealType var) {
261      // zero fill the random vector before starting:
262 <    std::vector<Vector3d> gaussRand;
262 >    vector<Vector3d> gaussRand;
263      gaussRand.resize(nTriangles);
264      std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
265      
# Line 228 | Line 267 | namespace OpenMD {
267      if (worldRank == 0) {
268   #endif
269        for (int i = 0; i < nTriangles; i++) {
270 <        gaussRand[i][0] = randNumGen_.randNorm(0.0, variance);
271 <        gaussRand[i][1] = randNumGen_.randNorm(0.0, variance);
272 <        gaussRand[i][2] = randNumGen_.randNorm(0.0, variance);
270 >        gaussRand[i][0] = randNumGen_.randNorm(0.0, var);
271 >        gaussRand[i][1] = randNumGen_.randNorm(0.0, var);
272 >        gaussRand[i][2] = randNumGen_.randNorm(0.0, var);
273        }
274   #ifdef IS_MPI
275      }
# Line 239 | Line 278 | namespace OpenMD {
278      // push these out to the other processors
279      
280   #ifdef IS_MPI
281 <    if (worldRank == 0) {
282 <      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 <    }
281 >    // Same command on all nodes:
282 >    MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
283   #endif
284      
285      return gaussRand;

Diff Legend

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