ViewVC Help

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

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