ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/integrators/LangevinHullForceManager.cpp
Revision: 1767
Committed: Fri Jul 6 22:01:58 2012 UTC (12 years, 9 months ago) by gezelter
File size: 8840 byte(s)
Log Message:
Various fixes required to compile OpenMD with the MS Visual C++ compiler

File Contents

# Content
1 /*
2 * Copyright (c) 2008, 2009, 2010 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
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).
39 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 */
42 #include <fstream>
43 #include <iostream>
44 #include "integrators/LangevinHullForceManager.hpp"
45 #include "utils/PhysicalConstants.hpp"
46 #include "math/ConvexHull.hpp"
47 #include "math/AlphaHull.hpp"
48 #include "math/Triangle.hpp"
49 #include "math/CholeskyDecomposition.hpp"
50 #ifdef IS_MPI
51 #include <mpi.h>
52 #endif
53
54 namespace OpenMD {
55
56 LangevinHullForceManager::LangevinHullForceManager(SimInfo* info) : ForceManager(info) {
57
58 simParams = info->getSimParams();
59 veloMunge = new Velocitizer(info);
60
61 // Create Hull, Convex Hull for now, other options later.
62
63 stringToEnumMap_["Convex"] = hullConvex;
64 stringToEnumMap_["AlphaShape"] = hullAlphaShape;
65 stringToEnumMap_["Unknown"] = hullUnknown;
66
67 const std::string ht = simParams->getHULL_Method();
68
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
78 switch(hullType_) {
79 case hullConvex :
80 surfaceMesh_ = new ConvexHull();
81 break;
82 case hullAlphaShape :
83 surfaceMesh_ = new AlphaHull(simParams->getAlpha());
84 break;
85 case hullUnknown :
86 default :
87 break;
88 }
89 /* Check that the simulation has target pressure and target
90 temperature set */
91
92 if (!simParams->haveTargetTemp()) {
93 sprintf(painCave.errMsg,
94 "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
95 "\twithout a targetTemp (K)!\n");
96 painCave.isFatal = 1;
97 painCave.severity = OPENMD_ERROR;
98 simError();
99 } else {
100 targetTemp_ = simParams->getTargetTemp();
101 }
102
103 if (!simParams->haveTargetPressure()) {
104 sprintf(painCave.errMsg,
105 "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
106 "\twithout a targetPressure (atm)!\n");
107 painCave.isFatal = 1;
108 simError();
109 } else {
110 // Convert pressure from atm -> amu/(fs^2*Ang)
111 targetPressure_ = simParams->getTargetPressure() /
112 PhysicalConstants::pressureConvert;
113 }
114
115 if (simParams->getUsePeriodicBoundaryConditions()) {
116 sprintf(painCave.errMsg,
117 "LangevinHullDynamics error: You can't use the Langevin Hull integrator\n"
118 "\twith periodic boundary conditions!\n");
119 painCave.isFatal = 1;
120 simError();
121 }
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 dt_ = simParams->getDt();
135
136 variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;
137
138 // Build a vector of integrable objects to determine if the are
139 // surface atoms
140 Molecule* mol;
141 StuntDouble* sd;
142 SimInfo::MoleculeIterator i;
143 Molecule::IntegrableObjectIterator j;
144
145 for (mol = info_->beginMolecule(i); mol != NULL;
146 mol = info_->nextMolecule(i)) {
147 for (sd = mol->beginIntegrableObject(j);
148 sd != NULL;
149 sd = mol->nextIntegrableObject(j)) {
150 localSites_.push_back(sd);
151 }
152 }
153 }
154
155 void LangevinHullForceManager::postCalculation(){
156
157 // Compute surface Mesh
158 surfaceMesh_->computeHull(localSites_);
159
160 // Get total area and number of surface stunt doubles
161 RealType area = surfaceMesh_->getArea();
162 std::vector<Triangle> sMesh = surfaceMesh_->getMesh();
163 int nTriangles = sMesh.size();
164
165 // Generate all of the necessary random forces
166 std::vector<Vector3d> randNums = genTriangleForces(nTriangles, variance_);
167
168 // Loop over the mesh faces and apply external pressure to each
169 // of the faces
170 std::vector<Triangle>::iterator face;
171 std::vector<StuntDouble*>::iterator vertex;
172 int thisFacet = 0;
173 for (face = sMesh.begin(); face != sMesh.end(); ++face){
174 Triangle thisTriangle = *face;
175 std::vector<StuntDouble*> vertexSDs = thisTriangle.getVertices();
176 RealType thisArea = thisTriangle.getArea();
177 Vector3d unitNormal = thisTriangle.getUnitNormal();
178 //unitNormal.normalize();
179 Vector3d centroid = thisTriangle.getCentroid();
180 Vector3d facetVel = thisTriangle.getFacetVelocity();
181 RealType thisMass = thisTriangle.getFacetMass();
182 Mat3x3d hydroTensor = thisTriangle.computeHydrodynamicTensor(viscosity_);
183
184 hydroTensor *= PhysicalConstants::viscoConvert;
185 Mat3x3d S;
186 CholeskyDecomposition(hydroTensor, S);
187
188 Vector3d extPressure = -unitNormal * (targetPressure_ * thisArea) /
189 PhysicalConstants::energyConvert;
190
191 Vector3d randomForce = S * randNums[thisFacet++];
192 Vector3d dragForce = -hydroTensor * facetVel;
193
194 Vector3d langevinForce = (extPressure + randomForce + dragForce);
195
196 // Apply triangle force to stuntdouble vertices
197 for (vertex = vertexSDs.begin(); vertex != vertexSDs.end(); ++vertex){
198 if ((*vertex) != NULL){
199 Vector3d vertexForce = langevinForce / RealType(3.0);
200 (*vertex)->addFrc(vertexForce);
201 }
202 }
203 }
204
205 veloMunge->removeComDrift();
206 veloMunge->removeAngularDrift();
207
208 Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
209 currSnapshot->setVolume(surfaceMesh_->getVolume());
210 ForceManager::postCalculation();
211 }
212
213
214 std::vector<Vector3d> LangevinHullForceManager::genTriangleForces(int nTriangles,
215 RealType variance)
216 {
217
218 // zero fill the random vector before starting:
219 std::vector<Vector3d> gaussRand;
220 gaussRand.resize(nTriangles);
221 std::fill(gaussRand.begin(), gaussRand.end(), V3Zero);
222
223 #ifdef IS_MPI
224 if (worldRank == 0) {
225 #endif
226 for (int i = 0; i < nTriangles; i++) {
227 gaussRand[i][0] = randNumGen_.randNorm(0.0, variance);
228 gaussRand[i][1] = randNumGen_.randNorm(0.0, variance);
229 gaussRand[i][2] = randNumGen_.randNorm(0.0, variance);
230 }
231 #ifdef IS_MPI
232 }
233 #endif
234
235 // push these out to the other processors
236
237 #ifdef IS_MPI
238 if (worldRank == 0) {
239 MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
240 } else {
241 MPI::COMM_WORLD.Bcast(&gaussRand[0], nTriangles*3, MPI::REALTYPE, 0);
242 }
243 #endif
244
245 return gaussRand;
246 }
247 }

Properties

Name Value
svn:eol-style native
svn:executable *