ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/integrators/SMLDForceManager.cpp
Revision: 1145
Committed: Thu May 31 17:46:49 2007 UTC (17 years, 11 months ago) by chuckv
File size: 15803 byte(s)
Log Message:
Starting NPT Langevin. Added new force manager files.

File Contents

# User Rev Content
1 chuckv 1145 /*
2     * Copyright (c) 2005 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. Acknowledgement of the program authors must be made in any
10     * publication of scientific results based in part on use of the
11     * program. An acceptable form of acknowledgement is citation of
12     * the article in which the program was described (Matthew
13     * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14     * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15     * Parallel Simulation Engine for Molecular Dynamics,"
16     * J. Comput. Chem. 26, pp. 252-271 (2005))
17     *
18     * 2. Redistributions of source code must retain the above copyright
19     * notice, this list of conditions and the following disclaimer.
20     *
21     * 3. Redistributions in binary form must reproduce the above copyright
22     * notice, this list of conditions and the following disclaimer in the
23     * documentation and/or other materials provided with the
24     * distribution.
25     *
26     * This software is provided "AS IS," without a warranty of any
27     * kind. All express or implied conditions, representations and
28     * warranties, including any implied warranty of merchantability,
29     * fitness for a particular purpose or non-infringement, are hereby
30     * excluded. The University of Notre Dame and its licensors shall not
31     * be liable for any damages suffered by licensee as a result of
32     * using, modifying or distributing the software or its
33     * derivatives. In no event will the University of Notre Dame or its
34     * licensors be liable for any lost revenue, profit or data, or for
35     * direct, indirect, special, consequential, incidental or punitive
36     * damages, however caused and regardless of the theory of liability,
37     * arising out of the use of or inability to use software, even if the
38     * University of Notre Dame has been advised of the possibility of
39     * such damages.
40     */
41    
42     /*
43     Provides the force manager for Isothermal-Isobaric Langevin Dynamics where the stochastic force
44     is applied to the surface atoms anisotropically so as to provide a constant pressure. The
45     surface atoms are determined by computing the convex hull and then triangulating that hull. The force
46     applied to the facets of the triangulation and mapped back onto the surface atoms.
47     See: Kohanoff et.al. CHEMPHYSCHEM,2005,6,1848-1852.
48     */
49    
50     #include <fstream>
51     #include <iostream>
52     #include "integrators/SMLDForceManager.hpp"
53     #include "math/CholeskyDecomposition.hpp"
54     #include "utils/OOPSEConstant.hpp"
55     #include "hydrodynamics/Sphere.hpp"
56     #include "hydrodynamics/Ellipsoid.hpp"
57     #include "math/ConvexHull.hpp"
58     #include "openbabel/mol.hpp"
59    
60     using namespace OpenBabel;
61     namespace oopse {
62    
63     SMLDForceManager::SMLDForceManager(SimInfo* info) : ForceManager(info){
64     simParams = info->getSimParams();
65     veloMunge = new Velocitizer(info);
66    
67     sphericalBoundaryConditions_ = false;
68     if (simParams->getUseSphericalBoundaryConditions()) {
69     sphericalBoundaryConditions_ = true;
70     if (simParams->haveLangevinBufferRadius()) {
71     langevinBufferRadius_ = simParams->getLangevinBufferRadius();
72     } else {
73     sprintf( painCave.errMsg,
74     "langevinBufferRadius must be specified "
75     "when useSphericalBoundaryConditions is turned on.\n");
76     painCave.severity = OOPSE_ERROR;
77     painCave.isFatal = 1;
78     simError();
79     }
80    
81     if (simParams->haveFrozenBufferRadius()) {
82     frozenBufferRadius_ = simParams->getFrozenBufferRadius();
83     } else {
84     sprintf( painCave.errMsg,
85     "frozenBufferRadius must be specified "
86     "when useSphericalBoundaryConditions is turned on.\n");
87     painCave.severity = OOPSE_ERROR;
88     painCave.isFatal = 1;
89     simError();
90     }
91    
92     if (frozenBufferRadius_ < langevinBufferRadius_) {
93     sprintf( painCave.errMsg,
94     "frozenBufferRadius has been set smaller than the "
95     "langevinBufferRadius. This is probably an error.\n");
96     painCave.severity = OOPSE_WARNING;
97     painCave.isFatal = 0;
98     simError();
99     }
100     }
101    
102     // Build the hydroProp map:
103     std::map<std::string, HydroProp*> hydroPropMap;
104    
105     Molecule* mol;
106     StuntDouble* integrableObject;
107     SimInfo::MoleculeIterator i;
108     Molecule::IntegrableObjectIterator j;
109     bool needHydroPropFile = false;
110    
111     for (mol = info->beginMolecule(i); mol != NULL;
112     mol = info->nextMolecule(i)) {
113     for (integrableObject = mol->beginIntegrableObject(j);
114     integrableObject != NULL;
115     integrableObject = mol->nextIntegrableObject(j)) {
116    
117     if (integrableObject->isRigidBody()) {
118     RigidBody* rb = static_cast<RigidBody*>(integrableObject);
119     if (rb->getNumAtoms() > 1) needHydroPropFile = true;
120     }
121    
122     }
123     }
124    
125    
126     if (needHydroPropFile) {
127     if (simParams->haveHydroPropFile()) {
128     hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
129     } else {
130     sprintf( painCave.errMsg,
131     "HydroPropFile must be set to a file name if Langevin\n"
132     "\tDynamics is specified for rigidBodies which contain more\n"
133     "\tthan one atom. To create a HydroPropFile, run \"Hydro\".\n");
134     painCave.severity = OOPSE_ERROR;
135     painCave.isFatal = 1;
136     simError();
137     }
138    
139     for (mol = info->beginMolecule(i); mol != NULL;
140     mol = info->nextMolecule(i)) {
141     for (integrableObject = mol->beginIntegrableObject(j);
142     integrableObject != NULL;
143     integrableObject = mol->nextIntegrableObject(j)) {
144    
145     std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
146     if (iter != hydroPropMap.end()) {
147     hydroProps_.push_back(iter->second);
148     } else {
149     sprintf( painCave.errMsg,
150     "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
151     painCave.severity = OOPSE_ERROR;
152     painCave.isFatal = 1;
153     simError();
154     }
155     }
156     }
157     } else {
158    
159     std::map<std::string, HydroProp*> hydroPropMap;
160     for (mol = info->beginMolecule(i); mol != NULL;
161     mol = info->nextMolecule(i)) {
162     for (integrableObject = mol->beginIntegrableObject(j);
163     integrableObject != NULL;
164     integrableObject = mol->nextIntegrableObject(j)) {
165     Shape* currShape = NULL;
166     if (integrableObject->isDirectionalAtom()) {
167     DirectionalAtom* dAtom = static_cast<DirectionalAtom*>(integrableObject);
168     AtomType* atomType = dAtom->getAtomType();
169     if (atomType->isGayBerne()) {
170     DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
171    
172     GenericData* data = dAtomType->getPropertyByName("GayBerne");
173     if (data != NULL) {
174     GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
175    
176     if (gayBerneData != NULL) {
177     GayBerneParam gayBerneParam = gayBerneData->getData();
178     currShape = new Ellipsoid(V3Zero,
179     gayBerneParam.GB_d / 2.0,
180     gayBerneParam.GB_l / 2.0,
181     Mat3x3d::identity());
182     } else {
183     sprintf( painCave.errMsg,
184     "Can not cast GenericData to GayBerneParam\n");
185     painCave.severity = OOPSE_ERROR;
186     painCave.isFatal = 1;
187     simError();
188     }
189     } else {
190     sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
191     painCave.severity = OOPSE_ERROR;
192     painCave.isFatal = 1;
193     simError();
194     }
195     }
196     } else {
197     Atom* atom = static_cast<Atom*>(integrableObject);
198     AtomType* atomType = atom->getAtomType();
199     if (atomType->isLennardJones()){
200     GenericData* data = atomType->getPropertyByName("LennardJones");
201     if (data != NULL) {
202     LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
203    
204     if (ljData != NULL) {
205     LJParam ljParam = ljData->getData();
206     currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
207     } else {
208     sprintf( painCave.errMsg,
209     "Can not cast GenericData to LJParam\n");
210     painCave.severity = OOPSE_ERROR;
211     painCave.isFatal = 1;
212     simError();
213     }
214     }
215     } else {
216     int obanum = etab.GetAtomicNum((atom->getType()).c_str());
217     if (obanum != 0) {
218     currShape = new Sphere(atom->getPos(), etab.GetVdwRad(obanum));
219     } else {
220     sprintf( painCave.errMsg,
221     "Could not find atom type in default element.txt\n");
222     painCave.severity = OOPSE_ERROR;
223     painCave.isFatal = 1;
224     simError();
225     }
226     }
227     }
228     HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
229     std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
230     if (iter != hydroPropMap.end())
231     hydroProps_.push_back(iter->second);
232     else {
233     currHydroProp->complete();
234     hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
235     hydroProps_.push_back(currHydroProp);
236     }
237     }
238     }
239     }
240     variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
241     }
242    
243     std::map<std::string, HydroProp*> SMLDForceManager::parseFrictionFile(const std::string& filename) {
244     std::map<std::string, HydroProp*> props;
245     std::ifstream ifs(filename.c_str());
246     if (ifs.is_open()) {
247    
248     }
249    
250     const unsigned int BufferSize = 65535;
251     char buffer[BufferSize];
252     while (ifs.getline(buffer, BufferSize)) {
253     HydroProp* currProp = new HydroProp(buffer);
254     props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp));
255     }
256    
257     return props;
258     }
259    
260     void SMLDForceManager::postCalculation(bool needStress){
261     SimInfo::MoleculeIterator i;
262     Molecule::IntegrableObjectIterator j;
263     Molecule* mol;
264     StuntDouble* integrableObject;
265     Vector3d vel;
266     Vector3d pos;
267     Vector3d frc;
268     Mat3x3d A;
269     Mat3x3d Atrans;
270     Vector3d Tb;
271     Vector3d ji;
272     unsigned int index = 0;
273     bool doLangevinForces;
274     bool freezeMolecule;
275     int fdf;
276    
277     fdf = 0;
278    
279     for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
280    
281     doLangevinForces = true;
282     freezeMolecule = false;
283    
284     if (sphericalBoundaryConditions_) {
285    
286     Vector3d molPos = mol->getCom();
287     RealType molRad = molPos.length();
288    
289     doLangevinForces = false;
290    
291     if (molRad > langevinBufferRadius_) {
292     doLangevinForces = true;
293     freezeMolecule = false;
294     }
295     if (molRad > frozenBufferRadius_) {
296     doLangevinForces = false;
297     freezeMolecule = true;
298     }
299     }
300    
301     for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
302     integrableObject = mol->nextIntegrableObject(j)) {
303    
304     if (freezeMolecule)
305     fdf += integrableObject->freeze();
306    
307     if (doLangevinForces) {
308     vel =integrableObject->getVel();
309     if (integrableObject->isDirectional()){
310     //calculate angular velocity in lab frame
311     Mat3x3d I = integrableObject->getI();
312     Vector3d angMom = integrableObject->getJ();
313     Vector3d omega;
314    
315     if (integrableObject->isLinear()) {
316     int linearAxis = integrableObject->linearAxis();
317     int l = (linearAxis +1 )%3;
318     int m = (linearAxis +2 )%3;
319     omega[l] = angMom[l] /I(l, l);
320     omega[m] = angMom[m] /I(m, m);
321    
322     } else {
323     omega[0] = angMom[0] /I(0, 0);
324     omega[1] = angMom[1] /I(1, 1);
325     omega[2] = angMom[2] /I(2, 2);
326     }
327    
328     //apply friction force and torque at center of resistance
329     A = integrableObject->getA();
330     Atrans = A.transpose();
331     Vector3d rcr = Atrans * hydroProps_[index]->getCOR();
332     Vector3d vcdLab = vel + cross(omega, rcr);
333     Vector3d vcdBody = A* vcdLab;
334     Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omega);
335     Vector3d frictionForceLab = Atrans*frictionForceBody;
336     integrableObject->addFrc(frictionForceLab);
337     Vector3d frictionTorqueBody = - (hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omega);
338     Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
339     integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
340    
341     //apply random force and torque at center of resistance
342     Vector3d randomForceBody;
343     Vector3d randomTorqueBody;
344     genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
345     Vector3d randomForceLab = Atrans*randomForceBody;
346     Vector3d randomTorqueLab = Atrans* randomTorqueBody;
347     integrableObject->addFrc(randomForceLab);
348     integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));
349    
350     } else {
351     //spherical atom
352     Vector3d frictionForce = -(hydroProps_[index]->getXitt() * vel);
353     Vector3d randomForce;
354     Vector3d randomTorque;
355     genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
356    
357     integrableObject->addFrc(frictionForce+randomForce);
358     }
359     }
360    
361     ++index;
362    
363     }
364     }
365    
366     info_->setFdf(fdf);
367     veloMunge->removeComDrift();
368     // Remove angular drift if we are not using periodic boundary conditions.
369     if(!simParams->getUsePeriodicBoundaryConditions())
370     veloMunge->removeAngularDrift();
371    
372     ForceManager::postCalculation(needStress);
373     }
374    
375     void SMLDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
376    
377    
378     Vector<RealType, 6> Z;
379     Vector<RealType, 6> generalForce;
380    
381     Z[0] = randNumGen_.randNorm(0, variance);
382     Z[1] = randNumGen_.randNorm(0, variance);
383     Z[2] = randNumGen_.randNorm(0, variance);
384     Z[3] = randNumGen_.randNorm(0, variance);
385     Z[4] = randNumGen_.randNorm(0, variance);
386     Z[5] = randNumGen_.randNorm(0, variance);
387    
388    
389     generalForce = hydroProps_[index]->getS()*Z;
390    
391     force[0] = generalForce[0];
392     force[1] = generalForce[1];
393     force[2] = generalForce[2];
394     torque[0] = generalForce[3];
395     torque[1] = generalForce[4];
396     torque[2] = generalForce[5];
397    
398     }
399    
400     }

Properties

Name Value
svn:executable *