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 1868 by gezelter, Tue Apr 30 15:56:54 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 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
178    }  
179 <  
179 >
180 >  LangevinHullForceManager::~LangevinHullForceManager() {
181 >    delete surfaceMesh_;
182 >    delete veloMunge;
183 >  }
184 >  
185    void LangevinHullForceManager::postCalculation(){
156 –    SimInfo::MoleculeIterator i;
157 –    Molecule::IntegrableObjectIterator  j;
158 –    Molecule* mol;
159 –    StuntDouble* integrableObject;
186    
187 +    int nTriangles, thisFacet;
188 +    RealType area, thisArea, thisMass;
189 +    vector<Triangle> sMesh;
190 +    Triangle thisTriangle;
191 +    vector<Triangle>::iterator face;
192 +    vector<StuntDouble*> vertexSDs;
193 +    vector<StuntDouble*>::iterator vertex;
194 +
195 +    Vector3d unitNormal, centroid, facetVel;
196 +    Vector3d langevinForce, vertexForce;
197 +    Vector3d extPressure, randomForce, dragForce;
198 +
199 +    Mat3x3d hydroTensor, S;
200 +    vector<Vector3d> randNums;
201 +
202      // Compute surface Mesh
203      surfaceMesh_->computeHull(localSites_);
163 –
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      ForceManager::postCalculation();  
258    }
216 –  
217 –  
218 –  std::vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
219 –                                                             RealType variance)
220 –  {
259      
260 +  vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
261 +                                                               RealType var) {
262      // zero fill the random vector before starting:
263 <    std::vector<Vector3d> gaussRand;
263 >    vector<Vector3d> gaussRand;
264      gaussRand.resize(nTriangles);
265      std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
266      
# Line 228 | Line 268 | namespace OpenMD {
268      if (worldRank == 0) {
269   #endif
270        for (int i = 0; i < nTriangles; i++) {
271 <        gaussRand[i][0] = randNumGen_.randNorm(0.0, variance);
272 <        gaussRand[i][1] = randNumGen_.randNorm(0.0, variance);
273 <        gaussRand[i][2] = randNumGen_.randNorm(0.0, variance);
271 >        gaussRand[i][0] = randNumGen_.randNorm(0.0, var);
272 >        gaussRand[i][1] = randNumGen_.randNorm(0.0, var);
273 >        gaussRand[i][2] = randNumGen_.randNorm(0.0, var);
274        }
275   #ifdef IS_MPI
276      }

Diff Legend

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