ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/integrators/LDForceManager.cpp
Revision: 2759
Committed: Wed May 17 21:51:42 2006 UTC (18 years, 2 months ago) by tim
File size: 18004 byte(s)
Log Message:
Adding single precision capabilities to c++ side

File Contents

# Content
1 /*
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 #include <fstream>
42 #include "integrators/LDForceManager.hpp"
43 #include "math/CholeskyDecomposition.hpp"
44 #include "utils/OOPSEConstant.hpp"
45 #include "hydrodynamics/Sphere.hpp"
46 #include "hydrodynamics/Ellipsoid.hpp"
47 #include "openbabel/mol.hpp"
48
49 using namespace OpenBabel;
50 namespace oopse {
51
52 LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){
53 Globals* simParams = info->getSimParams();
54
55 sphericalBoundaryConditions_ = false;
56 if (simParams->getUseSphericalBoundaryConditions()) {
57 sphericalBoundaryConditions_ = true;
58 if (simParams->haveLangevinBufferRadius()) {
59 langevinBufferRadius_ = simParams->getLangevinBufferRadius();
60 } else {
61 sprintf( painCave.errMsg,
62 "langevinBufferRadius must be specified "
63 "when useSphericalBoundaryConditions is turned on.\n");
64 painCave.severity = OOPSE_ERROR;
65 painCave.isFatal = 1;
66 simError();
67 }
68
69 if (simParams->haveFrozenBufferRadius()) {
70 frozenBufferRadius_ = simParams->getFrozenBufferRadius();
71 } else {
72 sprintf( painCave.errMsg,
73 "frozenBufferRadius must be specified "
74 "when useSphericalBoundaryConditions is turned on.\n");
75 painCave.severity = OOPSE_ERROR;
76 painCave.isFatal = 1;
77 simError();
78 }
79
80 if (frozenBufferRadius_ < langevinBufferRadius_) {
81 sprintf( painCave.errMsg,
82 "frozenBufferRadius has been set smaller than the "
83 "langevinBufferRadius. This is probably an error.\n");
84 painCave.severity = OOPSE_WARNING;
85 painCave.isFatal = 0;
86 simError();
87 }
88 }
89
90 // Build the hydroProp map:
91 std::map<std::string, HydroProp> hydroPropMap;
92
93 Molecule* mol;
94 StuntDouble* integrableObject;
95 SimInfo::MoleculeIterator i;
96 Molecule::IntegrableObjectIterator j;
97 bool needHydroPropFile = false;
98
99 for (mol = info->beginMolecule(i); mol != NULL;
100 mol = info->nextMolecule(i)) {
101 for (integrableObject = mol->beginIntegrableObject(j);
102 integrableObject != NULL;
103 integrableObject = mol->nextIntegrableObject(j)) {
104
105 if (integrableObject->isRigidBody()) {
106 RigidBody* rb = static_cast<RigidBody*>(integrableObject);
107 if (rb->getNumAtoms() > 1) needHydroPropFile = true;
108 }
109
110 }
111 }
112
113
114 if (needHydroPropFile) {
115 if (simParams->haveHydroPropFile()) {
116 hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
117 } else {
118 sprintf( painCave.errMsg,
119 "HydroPropFile must be set to a file name if Langevin\n"
120 "\tDynamics is specified for rigidBodies which contain more\n"
121 "\tthan one atom. To create a HydroPropFile, run \"Hydro\".\n");
122 painCave.severity = OOPSE_ERROR;
123 painCave.isFatal = 1;
124 simError();
125 }
126 std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType());
127 if (iter != hydroPropMap.end()) {
128 hydroProps_.push_back(iter->second);
129 } else {
130 sprintf( painCave.errMsg,
131 "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
132 painCave.severity = OOPSE_ERROR;
133 painCave.isFatal = 1;
134 simError();
135 }
136 } else {
137
138 std::map<std::string, HydroProp> hydroPropMap;
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 Shape* currShape = NULL;
145 if (integrableObject->isDirectionalAtom()) {
146 DirectionalAtom* dAtom = static_cast<DirectionalAtom*>(integrableObject);
147 AtomType* atomType = dAtom->getAtomType();
148 if (atomType->isGayBerne()) {
149 DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
150
151 GenericData* data = dAtomType->getPropertyByName("GayBerne");
152 if (data != NULL) {
153 GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
154
155 if (gayBerneData != NULL) {
156 GayBerneParam gayBerneParam = gayBerneData->getData();
157 currShape = new Ellipsoid(V3Zero,
158 gayBerneParam.GB_sigma/2.0,
159 gayBerneParam.GB_l2b_ratio*gayBerneParam.GB_sigma/2.0,
160 Mat3x3d::identity());
161 } else {
162 sprintf( painCave.errMsg,
163 "Can not cast GenericData to GayBerneParam\n");
164 painCave.severity = OOPSE_ERROR;
165 painCave.isFatal = 1;
166 simError();
167 }
168 } else {
169 sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
170 painCave.severity = OOPSE_ERROR;
171 painCave.isFatal = 1;
172 simError();
173 }
174 }
175 } else {
176 Atom* atom = static_cast<Atom*>(integrableObject);
177 AtomType* atomType = atom->getAtomType();
178 if (atomType->isLennardJones()){
179 GenericData* data = atomType->getPropertyByName("LennardJones");
180 if (data != NULL) {
181 LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
182
183 if (ljData != NULL) {
184 LJParam ljParam = ljData->getData();
185 currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
186 } else {
187 sprintf( painCave.errMsg,
188 "Can not cast GenericData to LJParam\n");
189 painCave.severity = OOPSE_ERROR;
190 painCave.isFatal = 1;
191 simError();
192 }
193 }
194 } else {
195 int obanum = etab.GetAtomicNum((atom->getType()).c_str());
196 if (obanum != 0) {
197 currShape = new Sphere(atom->getPos(), etab.GetVdwRad(obanum));
198 } else {
199 sprintf( painCave.errMsg,
200 "Could not find atom type in default element.txt\n");
201 painCave.severity = OOPSE_ERROR;
202 painCave.isFatal = 1;
203 simError();
204 }
205 }
206 }
207 HydroProps currHydroProp = currShape->getHydroProps(simParams->getViscosity(),simParams->getTargetTemp());
208 std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType());
209 if (iter != hydroPropMap.end())
210 hydroProps_.push_back(iter->second);
211 else {
212 HydroProp myProp;
213 myProp.cor = V3Zero;
214 for (int i1 = 0; i1 < 3; i1++) {
215 for (int j1 = 0; j1 < 3; j1++) {
216 myProp.Xirtt(i1,j1) = currHydroProp.Xi(i1,j1);
217 myProp.Xirrt(i1,j1) = currHydroProp.Xi(i1,j1+3);
218 myProp.Xirtr(i1,j1) = currHydroProp.Xi(i1+3,j1);
219 myProp.Xirrr(i1,j1) = currHydroProp.Xi(i1+3,j1+3);
220 }
221 }
222 CholeskyDecomposition(currHydroProp.Xi, myProp.S);
223 hydroPropMap.insert(std::map<std::string, HydroProp>::value_type(integrableObject->getType(), myProp));
224 hydroProps_.push_back(myProp);
225 }
226 }
227 }
228 }
229 variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
230 }
231
232
233
234
235
236 std::map<std::string, HydroProp> LDForceManager::parseFrictionFile(const std::string& filename) {
237 std::map<std::string, HydroProp> props;
238 std::ifstream ifs(filename.c_str());
239 if (ifs.is_open()) {
240
241 }
242
243 const unsigned int BufferSize = 65535;
244 char buffer[BufferSize];
245 while (ifs.getline(buffer, BufferSize)) {
246 StringTokenizer tokenizer(buffer);
247 HydroProp currProp;
248 if (tokenizer.countTokens() >= 40) {
249 std::string atomName = tokenizer.nextToken();
250 currProp.cor[0] = tokenizer.nextTokenAsDouble();
251 currProp.cor[1] = tokenizer.nextTokenAsDouble();
252 currProp.cor[2] = tokenizer.nextTokenAsDouble();
253
254 currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble();
255 currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble();
256 currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble();
257 currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble();
258 currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble();
259 currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble();
260 currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble();
261 currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble();
262 currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble();
263
264 currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble();
265 currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble();
266 currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble();
267 currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble();
268 currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble();
269 currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble();
270 currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble();
271 currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble();
272 currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble();
273
274 currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble();
275 currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble();
276 currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble();
277 currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble();
278 currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble();
279 currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble();
280 currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble();
281 currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble();
282 currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble();
283
284 currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble();
285 currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble();
286 currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble();
287 currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble();
288 currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble();
289 currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble();
290 currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble();
291 currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
292 currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
293
294 SquareMatrix<RealType, 6> Xir;
295 Xir.setSubMatrix(0, 0, currProp.Xirtt);
296 Xir.setSubMatrix(0, 3, currProp.Xirrt);
297 Xir.setSubMatrix(3, 0, currProp.Xirtr);
298 Xir.setSubMatrix(3, 3, currProp.Xirrr);
299 CholeskyDecomposition(Xir, currProp.S);
300
301 props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
302 }
303 }
304
305 return props;
306 }
307
308 void LDForceManager::postCalculation() {
309 SimInfo::MoleculeIterator i;
310 Molecule::IntegrableObjectIterator j;
311 Molecule* mol;
312 StuntDouble* integrableObject;
313 Vector3d vel;
314 Vector3d pos;
315 Vector3d frc;
316 Mat3x3d A;
317 Mat3x3d Atrans;
318 Vector3d Tb;
319 Vector3d ji;
320 RealType mass;
321 unsigned int index = 0;
322 bool doLangevinForces;
323 bool freezeMolecule;
324 int fdf;
325
326 fdf = 0;
327 for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
328
329 if (sphericalBoundaryConditions_) {
330
331 Vector3d molPos = mol->getCom();
332 RealType molRad = molPos.length();
333
334 doLangevinForces = false;
335 freezeMolecule = false;
336
337 if (molRad > langevinBufferRadius_) {
338 doLangevinForces = true;
339 freezeMolecule = false;
340 }
341 if (molRad > frozenBufferRadius_) {
342 doLangevinForces = false;
343 freezeMolecule = true;
344 }
345 }
346
347 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
348 integrableObject = mol->nextIntegrableObject(j)) {
349
350 if (freezeMolecule)
351 fdf += integrableObject->freeze();
352
353 if (doLangevinForces) {
354 vel =integrableObject->getVel();
355 if (integrableObject->isDirectional()){
356 //calculate angular velocity in lab frame
357 Mat3x3d I = integrableObject->getI();
358 Vector3d angMom = integrableObject->getJ();
359 Vector3d omega;
360
361 if (integrableObject->isLinear()) {
362 int linearAxis = integrableObject->linearAxis();
363 int l = (linearAxis +1 )%3;
364 int m = (linearAxis +2 )%3;
365 omega[l] = angMom[l] /I(l, l);
366 omega[m] = angMom[m] /I(m, m);
367
368 } else {
369 omega[0] = angMom[0] /I(0, 0);
370 omega[1] = angMom[1] /I(1, 1);
371 omega[2] = angMom[2] /I(2, 2);
372 }
373
374 //apply friction force and torque at center of resistance
375 A = integrableObject->getA();
376 Atrans = A.transpose();
377 Vector3d rcr = Atrans * hydroProps_[index].cor;
378 Vector3d vcdLab = vel + cross(omega, rcr);
379 Vector3d vcdBody = A* vcdLab;
380 Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega);
381 Vector3d frictionForceLab = Atrans*frictionForceBody;
382 integrableObject->addFrc(frictionForceLab);
383 Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega);
384 Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
385 integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
386
387 //apply random force and torque at center of resistance
388 Vector3d randomForceBody;
389 Vector3d randomTorqueBody;
390 genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
391 Vector3d randomForceLab = Atrans*randomForceBody;
392 Vector3d randomTorqueLab = Atrans* randomTorqueBody;
393 integrableObject->addFrc(randomForceLab);
394 integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));
395
396 } else {
397 //spherical atom
398 Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel);
399 Vector3d randomForce;
400 Vector3d randomTorque;
401 genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
402
403 integrableObject->addFrc(frictionForce+randomForce);
404 }
405 }
406
407 ++index;
408
409 }
410 }
411 info_->setFdf(fdf);
412
413 ForceManager::postCalculation();
414 }
415
416 void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {
417
418
419 Vector<RealType, 6> Z;
420 Vector<RealType, 6> generalForce;
421
422
423 Z[0] = randNumGen_.randNorm(0, variance);
424 Z[1] = randNumGen_.randNorm(0, variance);
425 Z[2] = randNumGen_.randNorm(0, variance);
426 Z[3] = randNumGen_.randNorm(0, variance);
427 Z[4] = randNumGen_.randNorm(0, variance);
428 Z[5] = randNumGen_.randNorm(0, variance);
429
430
431 generalForce = hydroProps_[index].S*Z;
432
433 force[0] = generalForce[0];
434 force[1] = generalForce[1];
435 force[2] = generalForce[2];
436 torque[0] = generalForce[3];
437 torque[1] = generalForce[4];
438 torque[2] = generalForce[5];
439
440 }
441
442 }

Properties

Name Value
svn:executable *