ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/integrators/LDForceManager.cpp

Comparing trunk/src/integrators/LDForceManager.cpp (file contents):
Revision 1237 by gezelter, Fri Apr 18 16:55:15 2008 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

#	Line 6 | Line 6
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
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
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.
#	Line 37 | Line 28
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, 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		*/
42		#include <fstream>
43		#include <iostream>
44		#include "integrators/LDForceManager.hpp"
45		#include "math/CholeskyDecomposition.hpp"
46	<	#include "utils/OOPSEConstant.hpp"
46	>	#include "utils/PhysicalConstants.hpp"
47		#include "hydrodynamics/Sphere.hpp"
48		#include "hydrodynamics/Ellipsoid.hpp"
49		#include "utils/ElementsTable.hpp"
50	+	#include "types/LennardJonesAdapter.hpp"
51	+	#include "types/GayBerneAdapter.hpp"
52
53	<	namespace oopse {
53	>	namespace OpenMD {
54
55	<	LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) {
55	>	LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info),
56	>	maxIterNum_(4),
57	>	forceTolerance_(1e-6) {
58		simParams = info->getSimParams();
59		veloMunge = new Velocitizer(info);
60
#	Line 62 | Line 67 | namespace oopse {
67		sprintf( painCave.errMsg,
68		"langevinBufferRadius must be specified "
69		"when useSphericalBoundaryConditions is turned on.\n");
70	<	painCave.severity = OOPSE_ERROR;
70	>	painCave.severity = OPENMD_ERROR;
71		painCave.isFatal = 1;
72		simError();
73		}
#	Line 73 | Line 78 | namespace oopse {
78		sprintf( painCave.errMsg,
79		"frozenBufferRadius must be specified "
80		"when useSphericalBoundaryConditions is turned on.\n");
81	<	painCave.severity = OOPSE_ERROR;
81	>	painCave.severity = OPENMD_ERROR;
82		painCave.isFatal = 1;
83		simError();
84		}
#	Line 82 | Line 87 | namespace oopse {
87		sprintf( painCave.errMsg,
88		"frozenBufferRadius has been set smaller than the "
89		"langevinBufferRadius.  This is probably an error.\n");
90	<	painCave.severity = OOPSE_WARNING;
90	>	painCave.severity = OPENMD_WARNING;
91		painCave.isFatal = 0;
92		simError();
93		}
#	Line 92 | Line 97 | namespace oopse {
97		std::map<std::string, HydroProp*> hydroPropMap;
98
99		Molecule* mol;
100	<	StuntDouble* integrableObject;
100	>	StuntDouble* sd;
101		SimInfo::MoleculeIterator i;
102		Molecule::IntegrableObjectIterator  j;
103		bool needHydroPropFile = false;
104
105		for (mol = info->beginMolecule(i); mol != NULL;
106		mol = info->nextMolecule(i)) {
107	<	for (integrableObject = mol->beginIntegrableObject(j);
108	<	integrableObject != NULL;
109	<	integrableObject = mol->nextIntegrableObject(j)) {
107	>
108	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
109	>	sd = mol->nextIntegrableObject(j)) {
110
111	<	if (integrableObject->isRigidBody()) {
112	<	RigidBody* rb = static_cast<RigidBody*>(integrableObject);
111	>	if (sd->isRigidBody()) {
112	>	RigidBody* rb = static_cast<RigidBody*>(sd);
113		if (rb->getNumAtoms() > 1) needHydroPropFile = true;
114		}
115
#	Line 120 | Line 125 | namespace oopse {
125		"HydroPropFile must be set to a file name if Langevin Dynamics\n"
126		"\tis specified for rigidBodies which contain more than one atom\n"
127		"\tTo create a HydroPropFile, run the \"Hydro\" program.\n");
128	<	painCave.severity = OOPSE_ERROR;
128	>	painCave.severity = OPENMD_ERROR;
129		painCave.isFatal = 1;
130		simError();
131		}
132
133		for (mol = info->beginMolecule(i); mol != NULL;
134		mol = info->nextMolecule(i)) {
130	–	for (integrableObject = mol->beginIntegrableObject(j);
131	–	integrableObject != NULL;
132	–	integrableObject = mol->nextIntegrableObject(j)) {
135
136	<	std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
136	>	for (sd = mol->beginIntegrableObject(j);  sd != NULL;
137	>	sd = mol->nextIntegrableObject(j)) {
138	>
139	>	std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(sd->getType());
140		if (iter != hydroPropMap.end()) {
141		hydroProps_.push_back(iter->second);
142		} else {
143		sprintf( painCave.errMsg,
144	<	"Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
145	<	painCave.severity = OOPSE_ERROR;
144	>	"Can not find resistance tensor for atom [%s]\n", sd->getType().c_str());
145	>	painCave.severity = OPENMD_ERROR;
146		painCave.isFatal = 1;
147		simError();
148		}
#	Line 148 | Line 153 | namespace oopse {
153		std::map<std::string, HydroProp*> hydroPropMap;
154		for (mol = info->beginMolecule(i); mol != NULL;
155		mol = info->nextMolecule(i)) {
156	<	for (integrableObject = mol->beginIntegrableObject(j);
157	<	integrableObject != NULL;
158	<	integrableObject = mol->nextIntegrableObject(j)) {
156	>
157	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
158	>	sd = mol->nextIntegrableObject(j)) {
159	>
160		Shape* currShape = NULL;
161
162	<	if (integrableObject->isAtom()){
163	<	Atom* atom = static_cast<Atom*>(integrableObject);
162	>	if (sd->isAtom()){
163	>	Atom* atom = static_cast<Atom*>(sd);
164		AtomType* atomType = atom->getAtomType();
165	<	if (atomType->isGayBerne()) {
166	<	DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
167	<	GenericData* data = dAtomType->getPropertyByName("GayBerne");
168	<	if (data != NULL) {
169	<	GayBerneParamGenericData* gayBerneData = dynamic_cast<GayBerneParamGenericData*>(data);
164	<
165	<	if (gayBerneData != NULL) {
166	<	GayBerneParam gayBerneParam = gayBerneData->getData();
167	<	currShape = new Ellipsoid(V3Zero,
168	<	gayBerneParam.GB_l / 2.0,
169	<	gayBerneParam.GB_d / 2.0,
170	<	Mat3x3d::identity());
171	<	} else {
172	<	sprintf( painCave.errMsg,
173	<	"Can not cast GenericData to GayBerneParam\n");
174	<	painCave.severity = OOPSE_ERROR;
175	<	painCave.isFatal = 1;
176	<	simError();
177	<	}
178	<	} else {
179	<	sprintf( painCave.errMsg, "Can not find Parameters for GayBerne\n");
180	<	painCave.severity = OOPSE_ERROR;
181	<	painCave.isFatal = 1;
182	<	simError();
183	<	}
165	>	GayBerneAdapter gba = GayBerneAdapter(atomType);
166	>	if (gba.isGayBerne()) {
167	>	currShape = new Ellipsoid(V3Zero, gba.getL() / 2.0,
168	>	gba.getD() / 2.0,
169	>	Mat3x3d::identity());
170		} else {
171	<	if (atomType->isLennardJones()){
172	<	GenericData* data = atomType->getPropertyByName("LennardJones");
173	<	if (data != NULL) {
188	<	LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
189	<	if (ljData != NULL) {
190	<	LJParam ljParam = ljData->getData();
191	<	currShape = new Sphere(atom->getPos(), ljParam.sigma/2.0);
192	<	} else {
193	<	sprintf( painCave.errMsg,
194	<	"Can not cast GenericData to LJParam\n");
195	<	painCave.severity = OOPSE_ERROR;
196	<	painCave.isFatal = 1;
197	<	simError();
198	<	}
199	<	}
171	>	LennardJonesAdapter lja = LennardJonesAdapter(atomType);
172	>	if (lja.isLennardJones()){
173	>	currShape = new Sphere(atom->getPos(), lja.getSigma()/2.0);
174		} else {
175		int aNum = etab.GetAtomicNum((atom->getType()).c_str());
176		if (aNum != 0) {
#	Line 204 | Line 178 | namespace oopse {
178		} else {
179		sprintf( painCave.errMsg,
180		"Could not find atom type in default element.txt\n");
181	<	painCave.severity = OOPSE_ERROR;
181	>	painCave.severity = OPENMD_ERROR;
182		painCave.isFatal = 1;
183		simError();
184		}
185		}
186		}
187		}
188	+
189	+	if (!simParams->haveTargetTemp()) {
190	+	sprintf(painCave.errMsg, "You can't use LangevinDynamics without a targetTemp!\n");
191	+	painCave.isFatal = 1;
192	+	painCave.severity = OPENMD_ERROR;
193	+	simError();
194	+	}
195	+
196	+	if (!simParams->haveViscosity()) {
197	+	sprintf(painCave.errMsg, "You can't use LangevinDynamics without a viscosity!\n");
198	+	painCave.isFatal = 1;
199	+	painCave.severity = OPENMD_ERROR;
200	+	simError();
201	+	}
202	+
203	+
204		HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp());
205	<	std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType());
205	>	std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(sd->getType());
206		if (iter != hydroPropMap.end())
207		hydroProps_.push_back(iter->second);
208		else {
209		currHydroProp->complete();
210	<	hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp));
210	>	hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(sd->getType(), currHydroProp));
211		hydroProps_.push_back(currHydroProp);
212		}
213	+	delete currShape;
214		}
215		}
216		}
217	<	variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
217	>	variance_ = 2.0 * PhysicalConstants::kb*simParams->getTargetTemp()/simParams->getDt();
218		}
219
220		std::map<std::string, HydroProp*> LDForceManager::parseFrictionFile(const std::string& filename) {
#	Line 243 | Line 234 | namespace oopse {
234		return props;
235		}
236
237	<	void LDForceManager::postCalculation(bool needStress){
237	>	void LDForceManager::postCalculation(){
238		SimInfo::MoleculeIterator i;
239		Molecule::IntegrableObjectIterator  j;
240		Molecule* mol;
241	<	StuntDouble* integrableObject;
241	>	StuntDouble* sd;
242		RealType mass;
243		Vector3d pos;
244		Vector3d frc;
#	Line 284 | Line 275 | namespace oopse {
275		}
276		}
277
278	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
279	<	integrableObject = mol->nextIntegrableObject(j)) {
278	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
279	>	sd = mol->nextIntegrableObject(j)) {
280
281		if (freezeMolecule)
282	<	fdf += integrableObject->freeze();
282	>	fdf += sd->freeze();
283
284		if (doLangevinForces) {
285	<	mass = integrableObject->getMass();
286	<	if (integrableObject->isDirectional()){
285	>	mass = sd->getMass();
286	>	if (sd->isDirectional()){
287
288		// preliminaries for directional objects:
289
290	<	A = integrableObject->getA();
290	>	A = sd->getA();
291		Atrans = A.transpose();
292		Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR();
293
#	Line 307 | Line 298 | namespace oopse {
298		genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
299		Vector3d randomForceLab = Atrans * randomForceBody;
300		Vector3d randomTorqueLab = Atrans * randomTorqueBody;
301	<	integrableObject->addFrc(randomForceLab);
302	<	integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));
301	>	sd->addFrc(randomForceLab);
302	>	sd->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab ));
303
304	<	Mat3x3d I = integrableObject->getI();
304	>	Mat3x3d I = sd->getI();
305		Vector3d omegaBody;
306
307		// What remains contains velocity explicitly, but the velocity required
#	Line 320 | Line 311 | namespace oopse {
311
312		// this is the velocity at the half-step:
313
314	<	Vector3d vel =integrableObject->getVel();
315	<	Vector3d angMom = integrableObject->getJ();
314	>	Vector3d vel =sd->getVel();
315	>	Vector3d angMom = sd->getJ();
316
317		//estimate velocity at full-step using everything but friction forces:
318
319	<	frc = integrableObject->getFrc();
320	<	Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc;
319	>	frc = sd->getFrc();
320	>	Vector3d velStep = vel + (dt2_ /mass * PhysicalConstants::energyConvert) * frc;
321
322	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
323	<	Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb;
322	>	Tb = sd->lab2Body(sd->getTrq());
323	>	Vector3d angMomStep = angMom + (dt2_ * PhysicalConstants::energyConvert) * Tb;
324
325		Vector3d omegaLab;
326		Vector3d vcdLab;
#	Line 347 | Line 338 | namespace oopse {
338
339		for (int k = 0; k < maxIterNum_; k++) {
340
341	<	if (integrableObject->isLinear()) {
342	<	int linearAxis = integrableObject->linearAxis();
341	>	if (sd->isLinear()) {
342	>	int linearAxis = sd->linearAxis();
343		int l = (linearAxis +1 )%3;
344		int m = (linearAxis +2 )%3;
345		omegaBody[l] = angMomStep[l] /I(l, l);
#	Line 375 | Line 366 | namespace oopse {
366
367		// re-estimate velocities at full-step using friction forces:
368
369	<	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab);
370	<	angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody);
369	>	velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * (frc + frictionForceLab);
370	>	angMomStep = angMom + (dt2_ * PhysicalConstants::energyConvert) * (Tb + frictionTorqueBody);
371
372		// check for convergence (if the vectors have converged, fdot and tdot will both be 1.0):
373
#	Line 387 | Line 378 | namespace oopse {
378		break; // iteration ends here
379		}
380
381	<	integrableObject->addFrc(frictionForceLab);
382	<	integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
381	>	sd->addFrc(frictionForceLab);
382	>	sd->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab));
383
384
385		} else {
#	Line 397 | Line 388 | namespace oopse {
388		Vector3d randomForce;
389		Vector3d randomTorque;
390		genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
391	<	integrableObject->addFrc(randomForce);
391	>	sd->addFrc(randomForce);
392
393		// What remains contains velocity explicitly, but the velocity required
394		// is at the full step: v(t + h), while we have initially the velocity
#	Line 406 | Line 397 | namespace oopse {
397
398		// this is the velocity at the half-step:
399
400	<	Vector3d vel =integrableObject->getVel();
400	>	Vector3d vel =sd->getVel();
401
402		//estimate velocity at full-step using everything but friction forces:
403
404	<	frc = integrableObject->getFrc();
405	<	Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc;
404	>	frc = sd->getFrc();
405	>	Vector3d velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * frc;
406
407		Vector3d frictionForce(0.0);
408		Vector3d oldFF;  // used to test for convergence
#	Line 426 | Line 417 | namespace oopse {
417
418		// re-estimate velocities at full-step using friction forces:
419
420	<	velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce);
420	>	velStep = vel + (dt2_ / mass * PhysicalConstants::energyConvert) * (frc + frictionForce);
421
422		// check for convergence (if the vector has converged, fdot will be 1.0):
423
#	Line 436 | Line 427 | namespace oopse {
427		break; // iteration ends here
428		}
429
430	<	integrableObject->addFrc(frictionForce);
430	>	sd->addFrc(frictionForce);
431
432		}
433		}
#	Line 452 | Line 443 | namespace oopse {
443		if(!simParams->getUsePeriodicBoundaryConditions())
444		veloMunge->removeAngularDrift();
445
446	<	ForceManager::postCalculation(needStress);
446	>	ForceManager::postCalculation();
447		}
448
449		void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) {

Comparing trunk/src/integrators/LDForceManager.cpp (property svn:keywords):
Revision 1237 by gezelter, Fri Apr 18 16:55:15 2008 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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