ViewVC Help

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

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

Comparing trunk/src/integrators/Velocitizer.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

#	Line 1 | Line 1
1	<	/*
1	>	/*
2		* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3		*
4		* The University of Notre Dame grants you ("Licensee") a
#	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, 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	<
42	>
43		#include "integrators/Velocitizer.hpp"
44		#include "math/SquareMatrix3.hpp"
45		#include "primitives/Molecule.hpp"
46		#include "primitives/StuntDouble.hpp"
46	–	#include "math/randomSPRNG.hpp"
47
48	<	namespace oopse {
48	>	#ifndef IS_MPI
49	>	#include "math/SeqRandNumGen.hpp"
50	>	#else
51	>	#include "math/ParallelRandNumGen.hpp"
52	>	#endif
53
54	<	void Velocitizer::velocitize(double temperature) {
54	>	namespace OpenMD {
55	>
56	>	Velocitizer::Velocitizer(SimInfo* info) : info_(info), thermo(info) {
57	>
58	>	int seedValue;
59	>	Globals * simParams = info->getSimParams();
60	>
61	>	#ifndef IS_MPI
62	>	if (simParams->haveSeed()) {
63	>	seedValue = simParams->getSeed();
64	>	randNumGen_ = new SeqRandNumGen(seedValue);
65	>	}else {
66	>	randNumGen_ = new SeqRandNumGen();
67	>	}
68	>	#else
69	>	if (simParams->haveSeed()) {
70	>	seedValue = simParams->getSeed();
71	>	randNumGen_ = new ParallelRandNumGen(seedValue);
72	>	}else {
73	>	randNumGen_ = new ParallelRandNumGen();
74	>	}
75	>	#endif
76	>	}
77	>
78	>	Velocitizer::~Velocitizer() {
79	>	delete randNumGen_;
80	>	}
81	>
82	>	void Velocitizer::velocitize(RealType temperature) {
83		Vector3d aVel;
84		Vector3d aJ;
85		Mat3x3d I;
86	<	int l;
55	<	int m;
56	<	int n;
86	>	int l, m, n;
87		Vector3d vdrift;
88	<	double vbar;
89	<	/**@todo refactory kb */
90	<	const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
91	<	double av2;
92	<	double kebar;
93	<
88	>	RealType vbar;
89	>	/**@todo refactor kb */
90	>	const RealType kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
91	>	RealType av2;
92	>	RealType kebar;
93	>
94	>	Globals * simParams = info_->getSimParams();
95	>
96		SimInfo::MoleculeIterator i;
97		Molecule::IntegrableObjectIterator j;
98		Molecule * mol;
99	<	StuntDouble * integrableObject;
68	<	gaussianSPRNG gaussStream(info_->getSeed());
99	>	StuntDouble * sd;
100
101		kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
102
103		for( mol = info_->beginMolecule(i); mol != NULL;
104	<	mol = info_->nextMolecule(i) ) {
74	<	for( integrableObject = mol->beginIntegrableObject(j);
75	<	integrableObject != NULL;
76	<	integrableObject = mol->nextIntegrableObject(j) ) {
104	>	mol = info_->nextMolecule(i) ) {
105
106	<	// uses equipartition theory to solve for vbar in angstrom/fs
107	<
108	<	av2 = 2.0 * kebar / integrableObject->getMass();
109	<	vbar = sqrt(av2);
110	<
111	<	// picks random velocities from a gaussian distribution
112	<	// centered on vbar
113	<
114	<	for( int k = 0; k < 3; k++ ) {
115	<	aVel[k] = vbar * gaussStream.getGaussian();
116	<	}
117	<
118	<	integrableObject->setVel(aVel);
119	<
120	<	if (integrableObject->isDirectional()) {
121	<	I = integrableObject->getI();
122	<
123	<	if (integrableObject->isLinear()) {
124	<	l = integrableObject->linearAxis();
125	<	m = (l + 1) % 3;
126	<	n = (l + 2) % 3;
127	<
128	<	aJ[l] = 0.0;
129	<	vbar = sqrt(2.0 * kebar * I(m, m));
130	<	aJ[m] = vbar * gaussStream.getGaussian();
131	<	vbar = sqrt(2.0 * kebar * I(n, n));
132	<	aJ[n] = vbar * gaussStream.getGaussian();
133	<	} else {
134	<	for( int k = 0; k < 3; k++ ) {
135	<	vbar = sqrt(2.0 * kebar * I(k, k));
136	<	aJ[k] = vbar * gaussStream.getGaussian();
137	<	}
138	<	} // else isLinear
139	<
140	<	integrableObject->setJ(aJ);
141	<	}     //isDirectional
142	<	}
143	<	}             //end for (mol = beginMolecule(i); ...)
144	<
145	<
146	<
106	>	for( sd = mol->beginIntegrableObject(j); sd != NULL;
107	>	sd = mol->nextIntegrableObject(j) ) {
108	>
109	>	// uses equipartition theory to solve for vbar in angstrom/fs
110	>
111	>	av2 = 2.0 * kebar / sd->getMass();
112	>	vbar = sqrt(av2);
113	>
114	>	// picks random velocities from a gaussian distribution
115	>	// centered on vbar
116	>
117	>	for( int k = 0; k < 3; k++ ) {
118	>	aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0);
119	>	}
120	>	sd->setVel(aVel);
121	>
122	>	if (sd->isDirectional()) {
123	>	I = sd->getI();
124	>
125	>	if (sd->isLinear()) {
126	>	l = sd->linearAxis();
127	>	m = (l + 1) % 3;
128	>	n = (l + 2) % 3;
129	>
130	>	aJ[l] = 0.0;
131	>	vbar = sqrt(2.0 * kebar * I(m, m));
132	>	aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0);
133	>	vbar = sqrt(2.0 * kebar * I(n, n));
134	>	aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0);
135	>	} else {
136	>	for( int k = 0; k < 3; k++ ) {
137	>	vbar = sqrt(2.0 * kebar * I(k, k));
138	>	aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0);
139	>	}
140	>	}
141	>
142	>	sd->setJ(aJ);
143	>	}
144	>	}
145	>	}
146	>
147		removeComDrift();
148
149	<	}
149	>	// Remove angular drift if we are not using periodic boundary
150	>	// conditions:
151
152	<
153	<
154	<	void Velocitizer::removeComDrift() {
152	>	if(!simParams->getUsePeriodicBoundaryConditions()) removeAngularDrift();
153	>	}
154	>
155	>	void Velocitizer::removeComDrift() {
156		// Get the Center of Mass drift velocity.
157	<	Vector3d vdrift = info_->getComVel();
157	>	Vector3d vdrift = thermo.getComVel();
158
159		SimInfo::MoleculeIterator i;
160		Molecule::IntegrableObjectIterator j;
161		Molecule * mol;
162	<	StuntDouble * integrableObject;
162	>	StuntDouble * sd;
163
164		//  Corrects for the center of mass drift.
165		// sums all the momentum and divides by total mass.
166		for( mol = info_->beginMolecule(i); mol != NULL;
167	<	mol = info_->nextMolecule(i) ) {
138	<	for( integrableObject = mol->beginIntegrableObject(j);
139	<	integrableObject != NULL;
140	<	integrableObject = mol->nextIntegrableObject(j) ) {
141	<	integrableObject->setVel(integrableObject->getVel() - vdrift);
142	<	}
143	<	}
167	>	mol = info_->nextMolecule(i) ) {
168
169	<	}
169	>	for( sd = mol->beginIntegrableObject(j); sd != NULL;
170	>	sd = mol->nextIntegrableObject(j) ) {
171
172	+	sd->setVel(sd->getVel() - vdrift);
173	+
174	+	}
175	+	}
176	+	}
177	+
178	+	void Velocitizer::removeAngularDrift() {
179	+	// Get the Center of Mass drift velocity.
180	+
181	+	Vector3d vdrift;
182	+	Vector3d com;
183	+
184	+	thermo.getComAll(com, vdrift);
185	+
186	+	Mat3x3d inertiaTensor;
187	+	Vector3d angularMomentum;
188	+	Vector3d omega;
189	+
190	+	thermo.getInertiaTensor(inertiaTensor, angularMomentum);
191	+
192	+	// We now need the inverse of the inertia tensor.
193	+	inertiaTensor = inertiaTensor.inverse();
194	+	omega = inertiaTensor * angularMomentum;
195	+
196	+	SimInfo::MoleculeIterator i;
197	+	Molecule::IntegrableObjectIterator j;
198	+	Molecule* mol;
199	+	StuntDouble* sd;
200	+	Vector3d tempComPos;
201	+
202	+	// Corrects for the center of mass angular drift by summing all
203	+	// the angular momentum and dividing by the total mass.
204	+
205	+	for( mol = info_->beginMolecule(i); mol != NULL;
206	+	mol = info_->nextMolecule(i) ) {
207	+
208	+	for( sd = mol->beginIntegrableObject(j); sd != NULL;
209	+	sd = mol->nextIntegrableObject(j) ) {
210	+
211	+	tempComPos = sd->getPos() - com;
212	+	sd->setVel((sd->getVel() - vdrift) - cross(omega, tempComPos));
213	+
214	+	}
215	+	}
216	+	}
217		}

Comparing trunk/src/integrators/Velocitizer.cpp (property svn:keywords):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC

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

Diff Legend

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