ViewVC Help

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

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

Comparing:
trunk/src/integrators/NPT.cpp (file contents), Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
branches/development/src/integrators/NPT.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 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 <math.h>
#	Line 46 | Line 47
47		#include "integrators/NPT.hpp"
48		#include "math/SquareMatrix3.hpp"
49		#include "primitives/Molecule.hpp"
50	<	#include "utils/OOPSEConstant.hpp"
50	>	#include "utils/PhysicalConstants.hpp"
51		#include "utils/simError.h"
52
53		// Basic isotropic thermostating and barostating via the Melchionna
#	Line 59 | Line 60
60		//
61		//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
62
63	<	namespace oopse {
63	>	namespace OpenMD {
64
65	<	NPT::NPT(SimInfo* info) :
65	>	NPT::NPT(SimInfo* info) :
66		VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
67
68	<	Globals* simParams = info_->getSimParams();
68	>	Globals* simParams = info_->getSimParams();
69
70	<	if (!simParams->getUseInitXSstate()) {
70	>	if (!simParams->getUseIntialExtendedSystemState()) {
71		Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
72		currSnapshot->setChi(0.0);
73		currSnapshot->setIntegralOfChiDt(0.0);
74		currSnapshot->setEta(Mat3x3d(0.0));
75	<	}
75	>	}
76
77	<	if (!simParams->haveTargetTemp()) {
77	>	if (!simParams->haveTargetTemp()) {
78		sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
79		painCave.isFatal = 1;
80	<	painCave.severity = OOPSE_ERROR;
80	>	painCave.severity = OPENMD_ERROR;
81		simError();
82	<	} else {
82	>	} else {
83		targetTemp = simParams->getTargetTemp();
84	<	}
84	>	}
85
86	<	// We must set tauThermostat
87	<	if (!simParams->haveTauThermostat()) {
86	>	// We must set tauThermostat
87	>	if (!simParams->haveTauThermostat()) {
88		sprintf(painCave.errMsg, "If you use the constant temperature\n"
89	<	"\tintegrator, you must set tauThermostat_.\n");
89	>	"\tintegrator, you must set tauThermostat.\n");
90
91	<	painCave.severity = OOPSE_ERROR;
91	>	painCave.severity = OPENMD_ERROR;
92		painCave.isFatal = 1;
93		simError();
94	<	} else {
94	>	} else {
95		tauThermostat = simParams->getTauThermostat();
96	<	}
96	>	}
97
98	<	if (!simParams->haveTargetPressure()) {
98	>	if (!simParams->haveTargetPressure()) {
99		sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
100	<	"   without a targetPressure!\n");
100	>	"   without a targetPressure!\n");
101
102		painCave.isFatal = 1;
103		simError();
104	<	} else {
104	>	} else {
105		targetPressure = simParams->getTargetPressure();
106	<	}
106	>	}
107
108	<	if (!simParams->haveTauBarostat()) {
108	>	if (!simParams->haveTauBarostat()) {
109		sprintf(painCave.errMsg,
110		"If you use the NPT integrator, you must set tauBarostat.\n");
111	<	painCave.severity = OOPSE_ERROR;
111	>	painCave.severity = OPENMD_ERROR;
112		painCave.isFatal = 1;
113		simError();
114	<	} else {
114	>	} else {
115		tauBarostat = simParams->getTauBarostat();
116	<	}
116	>	}
117
118	<	tt2 = tauThermostat * tauThermostat;
119	<	tb2 = tauBarostat * tauBarostat;
118	>	tt2 = tauThermostat * tauThermostat;
119	>	tb2 = tauBarostat * tauBarostat;
120
121	<	update();
122	<	}
121	>	update();
122	>	}
123
124	<	NPT::~NPT() {
125	<	}
124	>	NPT::~NPT() {
125	>	}
126
127	<	void NPT::doUpdate() {
127	>	void NPT::doUpdate() {
128
129		oldPos.resize(info_->getNIntegrableObjects());
130		oldVel.resize(info_->getNIntegrableObjects());
131		oldJi.resize(info_->getNIntegrableObjects());
132
133	<	}
133	>	}
134
135	<	void NPT::moveA() {
135	>	void NPT::moveA() {
136		SimInfo::MoleculeIterator i;
137		Molecule::IntegrableObjectIterator  j;
138		Molecule* mol;
139		StuntDouble* integrableObject;
140		Vector3d Tb, ji;
141	<	double mass;
141	>	RealType mass;
142		Vector3d vel;
143		Vector3d pos;
144		Vector3d frc;
#	Line 150 | Line 151 | void NPT::moveA() {
151
152		instaTemp =thermo.getTemperature();
153		press = thermo.getPressureTensor();
154	<	instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
154	>	instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
155		instaVol =thermo.getVolume();
156
157		Vector3d  COM = info_->getCom();
#	Line 160 | Line 161 | void NPT::moveA() {
161		calcVelScale();
162
163		for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
164	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
165	<	integrableObject = mol->nextIntegrableObject(j)) {
164	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
165	>	integrableObject = mol->nextIntegrableObject(j)) {
166
167	<	vel = integrableObject->getVel();
168	<	frc = integrableObject->getFrc();
167	>	vel = integrableObject->getVel();
168	>	frc = integrableObject->getFrc();
169
170	<	mass = integrableObject->getMass();
170	>	mass = integrableObject->getMass();
171
172	<	getVelScaleA(sc, vel);
172	>	getVelScaleA(sc, vel);
173
174	<	// velocity half step  (use chi from previous step here):
175	<	//vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
176	<	vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
177	<	integrableObject->setVel(vel);
174	>	// velocity half step  (use chi from previous step here):
175	>	//vel[j] += dt2 * ((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
176	>	vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
177	>	integrableObject->setVel(vel);
178
179	<	if (integrableObject->isDirectional()) {
179	>	if (integrableObject->isDirectional()) {
180
181	<	// get and convert the torque to body frame
181	>	// get and convert the torque to body frame
182
183	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
183	>	Tb = integrableObject->lab2Body(integrableObject->getTrq());
184
185	<	// get the angular momentum, and propagate a half step
185	>	// get the angular momentum, and propagate a half step
186
187	<	ji = integrableObject->getJ();
187	>	ji = integrableObject->getJ();
188
189	<	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
190	<	ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
189	>	//ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
190	>	ji += dt2*PhysicalConstants::energyConvert * Tb - dt2*chi* ji;
191
192	<	rotAlgo->rotate(integrableObject, ji, dt);
192	>	rotAlgo->rotate(integrableObject, ji, dt);
193
194	<	integrableObject->setJ(ji);
195	<	}
194	>	integrableObject->setJ(ji);
195	>	}
196
197	<	}
197	>	}
198		}
199		// evolve chi and eta  half step
200
#	Line 206 | Line 207 | void NPT::moveA() {
207
208		index = 0;
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)) {
212	<	oldPos[index++] = integrableObject->getPos();
213	<	}
210	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
211	>	integrableObject = mol->nextIntegrableObject(j)) {
212	>	oldPos[index++] = integrableObject->getPos();
213	>	}
214		}
215
216		//the first estimation of r(t+dt) is equal to  r(t)
217
218		for(int k = 0; k < maxIterNum_; k++) {
219	<	index = 0;
220	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
221	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
222	<	integrableObject = mol->nextIntegrableObject(j)) {
219	>	index = 0;
220	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
221	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
222	>	integrableObject = mol->nextIntegrableObject(j)) {
223
224	<	vel = integrableObject->getVel();
225	<	pos = integrableObject->getPos();
224	>	vel = integrableObject->getVel();
225	>	pos = integrableObject->getPos();
226
227	<	this->getPosScale(pos, COM, index, sc);
227	>	this->getPosScale(pos, COM, index, sc);
228
229	<	pos = oldPos[index] + dt * (vel + sc);
230	<	integrableObject->setPos(pos);
229	>	pos = oldPos[index] + dt * (vel + sc);
230	>	integrableObject->setPos(pos);
231
232	<	++index;
233	<	}
234	<	}
232	>	++index;
233	>	}
234	>	}
235
236	<	rattle->constraintA();
236	>	rattle->constraintA();
237		}
238
239		// Scale the box after all the positions have been moved:
#	Line 243 | Line 244 | void NPT::moveA() {
244		currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
245
246		saveEta();
247	<	}
247	>	}
248
249	<	void NPT::moveB(void) {
249	>	void NPT::moveB(void) {
250		SimInfo::MoleculeIterator i;
251		Molecule::IntegrableObjectIterator  j;
252		Molecule* mol;
#	Line 256 | Line 257 | void NPT::moveB(void) {
257		Vector3d sc;
258		Vector3d vel;
259		Vector3d frc;
260	<	double mass;
260	>	RealType mass;
261
262
263		chi= currentSnapshot_->getChi();
264		integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
265	<	double oldChi  = chi;
266	<	double prevChi;
265	>	RealType oldChi  = chi;
266	>	RealType prevChi;
267
268		loadEta();
269
270		//save velocity and angular momentum
271		index = 0;
272		for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
273	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
274	<	integrableObject = mol->nextIntegrableObject(j)) {
273	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
274	>	integrableObject = mol->nextIntegrableObject(j)) {
275
276	<	oldVel[index] = integrableObject->getVel();
277	<	oldJi[index] = integrableObject->getJ();
278	<	++index;
279	<	}
276	>	oldVel[index] = integrableObject->getVel();
277	>	oldJi[index] = integrableObject->getJ();
278	>	++index;
279	>	}
280		}
281
282		// do the iteration:
283		instaVol =thermo.getVolume();
284
285		for(int k = 0; k < maxIterNum_; k++) {
286	<	instaTemp =thermo.getTemperature();
287	<	instaPress =thermo.getPressure();
286	>	instaTemp =thermo.getTemperature();
287	>	instaPress =thermo.getPressure();
288
289	<	// evolve chi another half step using the temperature at t + dt/2
290	<	prevChi = chi;
291	<	chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
289	>	// evolve chi another half step using the temperature at t + dt/2
290	>	prevChi = chi;
291	>	chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
292
293	<	//evolve eta
294	<	this->evolveEtaB();
295	<	this->calcVelScale();
293	>	//evolve eta
294	>	this->evolveEtaB();
295	>	this->calcVelScale();
296
297	<	index = 0;
298	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
299	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
300	<	integrableObject = mol->nextIntegrableObject(j)) {
297	>	index = 0;
298	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
299	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
300	>	integrableObject = mol->nextIntegrableObject(j)) {
301
302	<	frc = integrableObject->getFrc();
303	<	vel = integrableObject->getVel();
302	>	frc = integrableObject->getFrc();
303	>	vel = integrableObject->getVel();
304
305	<	mass = integrableObject->getMass();
305	>	mass = integrableObject->getMass();
306
307	<	getVelScaleB(sc, index);
307	>	getVelScaleB(sc, index);
308
309	<	// velocity half step
310	<	//vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
311	<	vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
312	<	integrableObject->setVel(vel);
309	>	// velocity half step
310	>	//vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
311	>	vel = oldVel[index] + dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
312	>	integrableObject->setVel(vel);
313
314	<	if (integrableObject->isDirectional()) {
315	<	// get and convert the torque to body frame
316	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
314	>	if (integrableObject->isDirectional()) {
315	>	// get and convert the torque to body frame
316	>	Tb = integrableObject->lab2Body(integrableObject->getTrq());
317
318	<	//ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
319	<	ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
320	<	integrableObject->setJ(ji);
321	<	}
318	>	//ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi[3*i+j]*chi);
319	>	ji = oldJi[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi*oldJi[index];
320	>	integrableObject->setJ(ji);
321	>	}
322
323	<	++index;
324	<	}
325	<	}
323	>	++index;
324	>	}
325	>	}
326
327	<	rattle->constraintB();
327	>	rattle->constraintB();
328
329	<	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
330	<	break;
329	>	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
330	>	break;
331		}
332
333		//calculate integral of chidt
#	Line 336 | Line 337 | void NPT::moveB(void) {
337		currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
338
339		saveEta();
340	<	}
340	>	}
341
342	+	void NPT::resetIntegrator(){
343	+	currentSnapshot_->setChi(0.0);
344	+	currentSnapshot_->setIntegralOfChiDt(0.0);
345	+	resetEta();
346	+	}
347	+
348	+
349	+	void NPT::resetEta() {
350	+	Mat3x3d etaMat(0.0);
351	+	currentSnapshot_->setEta(etaMat);
352	+	}
353	+
354		}

Comparing:
trunk/src/integrators/NPT.cpp (property svn:keywords), Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
branches/development/src/integrators/NPT.cpp (property svn:keywords), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

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

Diff Legend

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