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 1855 by gezelter, Tue Apr 2 18:31:51 2013 UTC vs.
Revision 1866 by gezelter, Thu Apr 25 14:32:56 2013 UTC

# 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      }
122 –    
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 –    
134 –    doThermalCoupling_ = true;
135 –    if ( fabs(viscosity_) < 1e-6 ) {
136 –      sprintf(painCave.errMsg,
137 –              "LangevinHullDynamics: The bath viscosity was set lower than\n"
138 –              "\t1e-6 poise.  OpenMD is turning off the thermal coupling to\n"
139 –              "\t the bath.\n");
140 –      painCave.isFatal = 0;
141 –      painCave.severity = OPENMD_INFO;
142 –      simError();
143 –      doThermalCoupling_ = false;
144 –    }
150  
151      dt_ = simParams->getDt();
147 –  
148 –    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* 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 (sd = mol->beginIntegrableObject(j);
# Line 161 | Line 167 | namespace OpenMD {
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    void LangevinHullForceManager::postCalculation(){
181    
182 +    int nTriangles, thisFacet;
183 +    RealType area, thisArea, thisMass;
184 +    vector<Triangle> sMesh;
185 +    Triangle thisTriangle;
186 +    vector<Triangle>::iterator face;
187 +    vector<StuntDouble*> vertexSDs;
188 +    vector<StuntDouble*>::iterator vertex;
189 +
190 +    Vector3d unitNormal, centroid, facetVel;
191 +    Vector3d langevinForce, vertexForce;
192 +    Vector3d extPressure, randomForce, dragForce;
193 +
194 +    Mat3x3d hydroTensor, S;
195 +    vector<Vector3d> randNums;
196 +
197      // Compute surface Mesh
198      surfaceMesh_->computeHull(localSites_);
171 –
199      // Get total area and number of surface stunt doubles
200 <    RealType area = surfaceMesh_->getArea();
201 <    std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
202 <    int nTriangles = sMesh.size();
200 >    area = surfaceMesh_->getArea();
201 >    sMesh = surfaceMesh_->getMesh();
202 >    nTriangles = sMesh.size();
203  
204 <    // Generate all of the necessary random forces
205 <    std::vector<Vector3d>  randNums = genTriangleForces(nTriangles, variance_);
206 <
207 <    // Loop over the mesh faces and apply external pressure to each
208 <    // of the faces
209 <    std::vector<Triangle>::iterator face;
210 <    std::vector<StuntDouble*>::iterator vertex;
184 <    int thisFacet = 0;
204 >    if (doThermalCoupling_) {
205 >      // Generate all of the necessary random forces
206 >      randNums = genTriangleForces(nTriangles, variance_);
207 >    }
208 >    
209 >    // Loop over the mesh faces
210 >    thisFacet = 0;
211      for (face = sMesh.begin(); face != sMesh.end(); ++face){
212 <      Triangle thisTriangle = *face;
213 <      std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
214 <      RealType thisArea = thisTriangle.getArea();
215 <      Vector3d unitNormal = thisTriangle.getUnitNormal();
216 <      //unitNormal.normalize();
217 <      Vector3d centroid = thisTriangle.getCentroid();
218 <      Vector3d facetVel = thisTriangle.getFacetVelocity();
193 <      RealType thisMass = thisTriangle.getFacetMass();
194 <      Mat3x3d hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);
195 <      
196 <      hydroTensor *= PhysicalConstants::viscoConvert;
197 <      Mat3x3d S;
198 <      CholeskyDecomposition(hydroTensor, S);
212 >      thisTriangle = *face;
213 >      vertexSDs = thisTriangle.getVertices();
214 >      thisArea = thisTriangle.getArea();
215 >      unitNormal = thisTriangle.getUnitNormal();
216 >      centroid = thisTriangle.getCentroid();
217 >      facetVel = thisTriangle.getFacetVelocity();
218 >      thisMass = thisTriangle.getFacetMass();
219  
220 <      Vector3d extPressure = -unitNormal * (targetPressure_ * thisArea) /
201 <        PhysicalConstants::energyConvert;
220 >      langevinForce = V3Zero;
221  
222 <      Vector3d randomForce(V3Zero);
223 <      Vector3d dragForce(V3Zero);
222 >      if (doPressureCoupling_) {
223 >        extPressure = -unitNormal * (targetPressure_ * thisArea) /
224 >          PhysicalConstants::energyConvert;
225 >        langevinForce += extPressure;
226 >      }
227 >
228        if (doThermalCoupling_) {
229 +        hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);      
230 +        hydroTensor *= PhysicalConstants::viscoConvert;
231 +        CholeskyDecomposition(hydroTensor, S);
232          randomForce = S * randNums[thisFacet++];
233          dragForce = -hydroTensor * facetVel;
234 +        langevinForce += randomForce + dragForce;
235        }
209 –
210 –      Vector3d langevinForce = (extPressure + randomForce + dragForce);
236        
237        // Apply triangle force to stuntdouble vertices
238        for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
239          if ((*vertex) != NULL){
240 <          Vector3d vertexForce = langevinForce / RealType(3.0);
240 >          vertexForce = langevinForce / RealType(3.0);
241            (*vertex)->addFrc(vertexForce);          
242          }  
243        }
# Line 222 | Line 247 | namespace OpenMD {
247      veloMunge->removeAngularDrift();
248      
249      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
250 <    currSnapshot->setVolume(surfaceMesh_->getVolume());    
250 >    currSnapshot->setVolume(surfaceMesh_->getVolume());  
251 >    currSnapshot->setHullVolume(surfaceMesh_->getVolume());
252      ForceManager::postCalculation();  
253    }
228 –  
229 –  
230 –  std::vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
231 –                                                             RealType variance)
232 –  {
254      
255 +  vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
256 +                                                               RealType var) {
257      // zero fill the random vector before starting:
258 <    std::vector<Vector3d> gaussRand;
258 >    vector<Vector3d> gaussRand;
259      gaussRand.resize(nTriangles);
260      std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
261      
# Line 240 | Line 263 | namespace OpenMD {
263      if (worldRank == 0) {
264   #endif
265        for (int i = 0; i < nTriangles; i++) {
266 <        gaussRand[i][0] = randNumGen_.randNorm(0.0, variance);
267 <        gaussRand[i][1] = randNumGen_.randNorm(0.0, variance);
268 <        gaussRand[i][2] = randNumGen_.randNorm(0.0, variance);
266 >        gaussRand[i][0] = randNumGen_.randNorm(0.0, var);
267 >        gaussRand[i][1] = randNumGen_.randNorm(0.0, var);
268 >        gaussRand[i][2] = randNumGen_.randNorm(0.0, var);
269        }
270   #ifdef IS_MPI
271      }

Diff Legend

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