ViewVC Help

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

root/OpenMD/branches/development/src/integrators/NVT.cpp

Comparing:
trunk/src/integrators/NVT.cpp (file contents), Revision 271 by tim, Mon Jan 17 16:58:32 2005 UTC vs.
branches/development/src/integrators/NVT.cpp (file contents), Revision 1715 by gezelter, Tue May 22 21:55:31 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
43		#include "integrators/NVT.hpp"
44		#include "primitives/Molecule.hpp"
45		#include "utils/simError.h"
46	<	#include "utils/OOPSEConstant.hpp"
46	>	#include "utils/PhysicalConstants.hpp"
47
48	<	namespace oopse {
48	>	namespace OpenMD {
49
50	<	NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6), maxIterNum_(4) {
50	>	NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6), maxIterNum_(4) {
51
52		Globals* simParams = info_->getSimParams();
53
54	<	if (!simParams->getUseInitXSstate()) {
55	<	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
56	<	currSnapshot->setChi(0.0);
57	<	currSnapshot->setIntegralOfChiDt(0.0);
54	>	if (!simParams->getUseIntialExtendedSystemState()) {
55	>	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
56	>	currSnapshot->setChi(0.0);
57	>	currSnapshot->setIntegralOfChiDt(0.0);
58		}
59
60		if (!simParams->haveTargetTemp()) {
61	<	sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
62	<	painCave.isFatal = 1;
63	<	painCave.severity = OOPSE_ERROR;
64	<	simError();
61	>	sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
62	>	painCave.isFatal = 1;
63	>	painCave.severity = OPENMD_ERROR;
64	>	simError();
65		} else {
66	<	targetTemp_ = simParams->getTargetTemp();
66	>	targetTemp_ = simParams->getTargetTemp();
67		}
68
69	<	// We must set tauThermostat_.
69	>	// We must set tauThermostat.
70
71		if (!simParams->haveTauThermostat()) {
72	<	sprintf(painCave.errMsg, "If you use the constant temperature\n"
73	<	"\tintegrator, you must set tauThermostat_.\n");
72	>	sprintf(painCave.errMsg, "If you use the constant temperature\n"
73	>	"\tintegrator, you must set tauThermostat.\n");
74
75	<	painCave.severity = OOPSE_ERROR;
76	<	painCave.isFatal = 1;
77	<	simError();
75	>	painCave.severity = OPENMD_ERROR;
76	>	painCave.isFatal = 1;
77	>	simError();
78		} else {
79	<	tauThermostat_ = simParams->getTauThermostat();
79	>	tauThermostat_ = simParams->getTauThermostat();
80		}
81
82	<	update();
83	<	}
82	>	updateSizes();
83	>	}
84
85	<	void NVT::doUpdate() {
85	>	void NVT::doUpdateSizes() {
86		oldVel_.resize(info_->getNIntegrableObjects());
87	<	oldJi_.resize(info_->getNIntegrableObjects());
88	<	}
89	<	void NVT::moveA() {
87	>	oldJi_.resize(info_->getNIntegrableObjects());
88	>	}
89	>	void NVT::moveA() {
90		SimInfo::MoleculeIterator i;
91		Molecule::IntegrableObjectIterator  j;
92		Molecule* mol;
93		StuntDouble* integrableObject;
94		Vector3d Tb;
95		Vector3d ji;
96	<	double mass;
96	>	RealType mass;
97		Vector3d vel;
98		Vector3d pos;
99		Vector3d frc;
100
101	<	double chi = currentSnapshot_->getChi();
102	<	double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
101	>	RealType chi = currentSnapshot_->getChi();
102	>	RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
103
104		// We need the temperature at time = t for the chi update below:
105
106	<	double instTemp = thermo.getTemperature();
106	>	RealType instTemp = thermo.getTemperature();
107
108		for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
109	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
110	<	integrableObject = mol->nextIntegrableObject(j)) {
109	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
110	>	integrableObject = mol->nextIntegrableObject(j)) {
111
112		vel = integrableObject->getVel();
113		pos = integrableObject->getPos();
#	Line 115 | Line 116 | void NVT::moveA() {
116		mass = integrableObject->getMass();
117
118		// velocity half step  (use chi from previous step here):
119	<	//vel[j] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
120	<	vel += dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
119	>	//vel[j] += dt2 * ((frc[j] / mass ) * PhysicalConstants::energyConvert - vel[j]*chi);
120	>	vel += dt2 *PhysicalConstants::energyConvert/mass*frc - dt2*chi*vel;
121
122		// position whole step
123		//pos[j] += dt * vel[j];
#	Line 127 | Line 128 | void NVT::moveA() {
128
129		if (integrableObject->isDirectional()) {
130
131	<	//convert the torque to body frame
132	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
131	>	//convert the torque to body frame
132	>	Tb = integrableObject->lab2Body(integrableObject->getTrq());
133
134	<	// get the angular momentum, and propagate a half step
134	>	// get the angular momentum, and propagate a half step
135
136	<	ji = integrableObject->getJ();
136	>	ji = integrableObject->getJ();
137
138	<	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
139	<	ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji;
140	<	rotAlgo->rotate(integrableObject, ji, dt);
138	>	//ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
139	>	ji += dt2*PhysicalConstants::energyConvert*Tb - dt2*chi *ji;
140	>	rotAlgo_->rotate(integrableObject, ji, dt);
141
142	<	integrableObject->setJ(ji);
142	>	integrableObject->setJ(ji);
143		}
144	<	}
144	>	}
145
146		}
147
148	<	rattle->constraintA();
148	>	flucQ_->moveA();
149	>	rattle_->constraintA();
150
151		// Finally, evolve chi a half step (just like a velocity) using
152		// temperature at time t, not time t+dt/2
#	Line 155 | Line 157 | void NVT::moveA() {
157
158		currentSnapshot_->setChi(chi);
159		currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
160	<	}
160	>	}
161
162	<	void NVT::moveB() {
162	>	void NVT::moveB() {
163		SimInfo::MoleculeIterator i;
164		Molecule::IntegrableObjectIterator  j;
165		Molecule* mol;
#	Line 167 | Line 169 | void NVT::moveB() {
169		Vector3d ji;
170		Vector3d vel;
171		Vector3d frc;
172	<	double mass;
173	<	double instTemp;
172	>	RealType mass;
173	>	RealType instTemp;
174		int index;
175		// Set things up for the iteration:
176
177	<	double chi = currentSnapshot_->getChi();
178	<	double oldChi = chi;
179	<	double  prevChi;
180	<	double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
177	>	RealType chi = currentSnapshot_->getChi();
178	>	RealType oldChi = chi;
179	>	RealType  prevChi;
180	>	RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
181
182		index = 0;
183		for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
184	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
185	<	integrableObject = mol->nextIntegrableObject(j)) {
184	<	oldVel_[index] = integrableObject->getVel();
185	<	oldJi_[index] = integrableObject->getJ();
184	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
185	>	integrableObject = mol->nextIntegrableObject(j)) {
186
187	<	++index;
188	<	}
189	<
187	>	oldVel_[index] = integrableObject->getVel();
188	>
189	>	if (integrableObject->isDirectional())
190	>	oldJi_[index] = integrableObject->getJ();
191	>
192	>	++index;
193	>	}
194		}
195
196		// do the iteration:
197
198		for(int k = 0; k < maxIterNum_; k++) {
199	<	index = 0;
200	<	instTemp = thermo.getTemperature();
199	>	index = 0;
200	>	instTemp = thermo.getTemperature();
201
202	<	// evolve chi another half step using the temperature at t + dt/2
202	>	// evolve chi another half step using the temperature at t + dt/2
203
204	<	prevChi = chi;
205	<	chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
204	>	prevChi = chi;
205	>	chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
206
207	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
208	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
209	<	integrableObject = mol->nextIntegrableObject(j)) {
207	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
208	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
209	>	integrableObject = mol->nextIntegrableObject(j)) {
210
211	<	frc = integrableObject->getFrc();
212	<	vel = integrableObject->getVel();
211	>	frc = integrableObject->getFrc();
212	>	vel = integrableObject->getVel();
213
214	<	mass = integrableObject->getMass();
214	>	mass = integrableObject->getMass();
215
216	<	// velocity half step
217	<	//for(j = 0; j < 3; j++)
218	<	//    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - oldVel_[3*i + j]*chi);
219	<	vel = oldVel_[index] + dt2/mass*OOPSEConstant::energyConvert * frc - dt2*chi*oldVel_[index];
216	>	// velocity half step
217	>	//for(j = 0; j < 3; j++)
218	>	//    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * PhysicalConstants::energyConvert - oldVel_[3*i + j]*chi);
219	>	vel = oldVel_[index] + dt2/mass*PhysicalConstants::energyConvert * frc - dt2*chi*oldVel_[index];
220
221	<	integrableObject->setVel(vel);
221	>	integrableObject->setVel(vel);
222
223	<	if (integrableObject->isDirectional()) {
223	>	if (integrableObject->isDirectional()) {
224
225	<	// get and convert the torque to body frame
225	>	// get and convert the torque to body frame
226
227	<	Tb =  integrableObject->lab2Body(integrableObject->getTrq());
227	>	Tb =  integrableObject->lab2Body(integrableObject->getTrq());
228
229	<	//for(j = 0; j < 3; j++)
230	<	//    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3*i+j]*chi);
231	<	ji = oldJi_[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *oldJi_[index];
229	>	//for(j = 0; j < 3; j++)
230	>	//    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi_[3*i+j]*chi);
231	>	ji = oldJi_[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi *oldJi_[index];
232
233	<	integrableObject->setJ(ji);
234	<	}
233	>	integrableObject->setJ(ji);
234	>	}
235
236
237	<	++index;
238	<	}
239	<	}
237	>	++index;
238	>	}
239	>	}
240
241	+	rattle_->constraintB();
242
243	<	rattle->constraintB();
243	>	if (fabs(prevChi - chi) <= chiTolerance_)
244	>	break;
245
240	–	if (fabs(prevChi - chi) <= chiTolerance_)
241	–	break;
242	–
246		}
247
248	<	integralOfChidt += dt2 * chi;
248	>	flucQ_->moveB();
249
250	+	integralOfChidt += dt2 * chi;
251		currentSnapshot_->setChi(chi);
252		currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
253	<	}
253	>	}
254
255	+	void NVT::resetIntegrator() {
256	+	currentSnapshot_->setChi(0.0);
257	+	currentSnapshot_->setIntegralOfChiDt(0.0);
258	+	}
259	+
260	+	RealType NVT::calcConservedQuantity() {
261
262	<	double NVT::calcConservedQuantity() {
263	<
264	<	double chi = currentSnapshot_->getChi();
265	<	double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
266	<	double conservedQuantity;
267	<	double fkBT;
268	<	double Energy;
259	<	double thermostat_kinetic;
260	<	double thermostat_potential;
262	>	RealType chi = currentSnapshot_->getChi();
263	>	RealType integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
264	>	RealType conservedQuantity;
265	>	RealType fkBT;
266	>	RealType Energy;
267	>	RealType thermostat_kinetic;
268	>	RealType thermostat_potential;
269
270	<	fkBT = info_->getNdf() *OOPSEConstant::kB *targetTemp_;
270	>	fkBT = info_->getNdf() *PhysicalConstants::kB *targetTemp_;
271
272		Energy = thermo.getTotalE();
273
274	<	thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);
274	>	thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * PhysicalConstants::energyConvert);
275
276	<	thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;
276	>	thermostat_potential = fkBT * integralOfChidt / PhysicalConstants::energyConvert;
277
278		conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
279
280		return conservedQuantity;
281	<	}
281	>	}
282
283
284	<	}//end namespace oopse
284	>	}//end namespace OpenMD

Comparing:
trunk/src/integrators/NVT.cpp (property svn:keywords), Revision 271 by tim, Mon Jan 17 16:58:32 2005 UTC vs.
branches/development/src/integrators/NVT.cpp (property svn:keywords), Revision 1715 by gezelter, Tue May 22 21:55:31 2012 UTC

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

Diff Legend

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