ViewVC Help

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

root/OpenMD/trunk/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.
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 <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));
74	<	}
72	>	currSnapshot->setThermostat(make_pair(0.0, 0.0));
73	>	currSnapshot->setBarostat(Mat3x3d(0.0));
74	>	}
75
76	<	if (!simParams->haveTargetTemp()) {
76	>	if (!simParams->haveTargetTemp()) {
77		sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
78		painCave.isFatal = 1;
79	<	painCave.severity = OOPSE_ERROR;
79	>	painCave.severity = OPENMD_ERROR;
80		simError();
81	<	} else {
81	>	} else {
82		targetTemp = simParams->getTargetTemp();
83	<	}
83	>	}
84
85	<	// We must set tauThermostat
86	<	if (!simParams->haveTauThermostat()) {
85	>	// We must set tauThermostat
86	>	if (!simParams->haveTauThermostat()) {
87		sprintf(painCave.errMsg, "If you use the constant temperature\n"
88	<	"\tintegrator, you must set tauThermostat_.\n");
88	>	"\tintegrator, you must set tauThermostat.\n");
89
90	<	painCave.severity = OOPSE_ERROR;
90	>	painCave.severity = OPENMD_ERROR;
91		painCave.isFatal = 1;
92		simError();
93	<	} else {
93	>	} else {
94		tauThermostat = simParams->getTauThermostat();
95	<	}
95	>	}
96
97	<	if (!simParams->haveTargetPressure()) {
97	>	if (!simParams->haveTargetPressure()) {
98		sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
99	<	"   without a targetPressure!\n");
99	>	"   without a targetPressure!\n");
100
101		painCave.isFatal = 1;
102		simError();
103	<	} else {
103	>	} else {
104		targetPressure = simParams->getTargetPressure();
105	<	}
105	>	}
106
107	<	if (!simParams->haveTauBarostat()) {
107	>	if (!simParams->haveTauBarostat()) {
108		sprintf(painCave.errMsg,
109		"If you use the NPT integrator, you must set tauBarostat.\n");
110	<	painCave.severity = OOPSE_ERROR;
110	>	painCave.severity = OPENMD_ERROR;
111		painCave.isFatal = 1;
112		simError();
113	<	} else {
113	>	} else {
114		tauBarostat = simParams->getTauBarostat();
115	<	}
115	>	}
116
117	<	tt2 = tauThermostat * tauThermostat;
118	<	tb2 = tauBarostat * tauBarostat;
117	>	tt2 = tauThermostat * tauThermostat;
118	>	tb2 = tauBarostat * tauBarostat;
119
120	<	update();
121	<	}
120	>	updateSizes();
121	>	}
122
123	<	NPT::~NPT() {
124	<	}
123	>	NPT::~NPT() {
124	>	}
125
126	<	void NPT::doUpdate() {
126	>	void NPT::doUpdateSizes() {
127
128		oldPos.resize(info_->getNIntegrableObjects());
129		oldVel.resize(info_->getNIntegrableObjects());
130		oldJi.resize(info_->getNIntegrableObjects());
131
132	<	}
132	>	}
133
134	<	void NPT::moveA() {
134	>	void NPT::moveA() {
135		SimInfo::MoleculeIterator i;
136		Molecule::IntegrableObjectIterator  j;
137		Molecule* mol;
138	<	StuntDouble* integrableObject;
138	>	StuntDouble* sd;
139		Vector3d Tb, ji;
140	<	double mass;
140	>	RealType mass;
141		Vector3d vel;
142		Vector3d pos;
143		Vector3d frc;
144		Vector3d sc;
145		int index;
146
147	<	chi= currentSnapshot_->getChi();
148	<	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
147	>	thermostat = snap->getThermostat();
148		loadEta();
149
150		instaTemp =thermo.getTemperature();
151		press = thermo.getPressureTensor();
152	<	instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
152	>	instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
153		instaVol =thermo.getVolume();
154
155	<	Vector3d  COM = info_->getCom();
155	>	Vector3d  COM = thermo.getCom();
156
157		//evolve velocity half step
158
159		calcVelScale();
160
161	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
162	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
163	<	integrableObject = mol->nextIntegrableObject(j)) {
161	>	for (mol = info_->beginMolecule(i); mol != NULL;
162	>	mol = info_->nextMolecule(i)) {
163	>
164	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
165	>	sd = mol->nextIntegrableObject(j)) {
166
167	<	vel = integrableObject->getVel();
168	<	frc = integrableObject->getFrc();
167	>	vel = sd->getVel();
168	>	frc = sd->getFrc();
169
170	<	mass = integrableObject->getMass();
170	>	mass = sd->getMass();
171
172	<	getVelScaleA(sc, vel);
172	>	getVelScaleA(sc, vel);
173
174	<	// velocity half step  (use chi from previous step here):
174	<	//vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
175	<	vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
176	<	integrableObject->setVel(vel);
174	>	// velocity half step  (use chi from previous step here):
175
176	<	if (integrableObject->isDirectional()) {
176	>	vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
177	>	sd->setVel(vel);
178
179	<	// get and convert the torque to body frame
179	>	if (sd->isDirectional()) {
180
181	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
181	>	// get and convert the torque to body frame
182
183	<	// get the angular momentum, and propagate a half step
183	>	Tb = sd->lab2Body(sd->getTrq());
184
185	<	ji = integrableObject->getJ();
185	>	// get the angular momentum, and propagate a half step
186
187	<	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
188	<	ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
187	>	ji = sd->getJ();
188	>
189	>	ji += dt2*PhysicalConstants::energyConvert * Tb
190	>	- dt2*thermostat.first* ji;
191
192	<	rotAlgo->rotate(integrableObject, ji, dt);
192	>	rotAlgo_->rotate(sd, ji, dt);
193
194	<	integrableObject->setJ(ji);
195	<	}
194	>	sd->setJ(ji);
195	>	}
196
197	<	}
197	>	}
198		}
199		// evolve chi and eta  half step
200
201	<	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
201	>	thermostat.first += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
202
203		evolveEtaA();
204
205		//calculate the integral of chidt
206	<	integralOfChidt += dt2 * chi;
206	>	thermostat.second += dt2 * thermostat.first;
207
208	+	flucQ_->moveA();
209	+
210	+
211		index = 0;
212	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
213	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
214	<	integrableObject = mol->nextIntegrableObject(j)) {
215	<	oldPos[index++] = integrableObject->getPos();
216	<	}
212	>	for (mol = info_->beginMolecule(i); mol != NULL;
213	>	mol = info_->nextMolecule(i)) {
214	>
215	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
216	>	sd = mol->nextIntegrableObject(j)) {
217	>
218	>	oldPos[index++] = sd->getPos();
219	>
220	>	}
221		}
222
223		//the first estimation of r(t+dt) is equal to  r(t)
224
225		for(int k = 0; k < maxIterNum_; k++) {
226	<	index = 0;
227	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
228	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
221	<	integrableObject = mol->nextIntegrableObject(j)) {
226	>	index = 0;
227	>	for (mol = info_->beginMolecule(i); mol != NULL;
228	>	mol = info_->nextMolecule(i)) {
229
230	<	vel = integrableObject->getVel();
231	<	pos = integrableObject->getPos();
230	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
231	>	sd = mol->nextIntegrableObject(j)) {
232
233	<	this->getPosScale(pos, COM, index, sc);
233	>	vel = sd->getVel();
234	>	pos = sd->getPos();
235	>
236	>	this->getPosScale(pos, COM, index, sc);
237
238	<	pos = oldPos[index] + dt * (vel + sc);
239	<	integrableObject->setPos(pos);
238	>	pos = oldPos[index] + dt * (vel + sc);
239	>	sd->setPos(pos);
240
241	<	++index;
242	<	}
243	<	}
241	>	++index;
242	>	}
243	>	}
244
245	<	rattle->constraintA();
245	>	rattle_->constraintA();
246		}
247
248		// Scale the box after all the positions have been moved:
249
250		this->scaleSimBox();
251
252	<	currentSnapshot_->setChi(chi);
243	<	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
252	>	snap->setThermostat(thermostat);
253
254		saveEta();
255	<	}
255	>	}
256
257	<	void NPT::moveB(void) {
257	>	void NPT::moveB(void) {
258		SimInfo::MoleculeIterator i;
259		Molecule::IntegrableObjectIterator  j;
260		Molecule* mol;
261	<	StuntDouble* integrableObject;
261	>	StuntDouble* sd;
262		int index;
263		Vector3d Tb;
264		Vector3d ji;
265		Vector3d sc;
266		Vector3d vel;
267		Vector3d frc;
268	<	double mass;
268	>	RealType mass;
269
270	+	thermostat = snap->getThermostat();
271	+	RealType oldChi  = thermostat.first;
272	+	RealType prevChi;
273
262	–	chi= currentSnapshot_->getChi();
263	–	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
264	–	double oldChi  = chi;
265	–	double prevChi;
266	–
274		loadEta();
275
276		//save velocity and angular momentum
277		index = 0;
278	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
279	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
280	<	integrableObject = mol->nextIntegrableObject(j)) {
278	>	for (mol = info_->beginMolecule(i); mol != NULL;
279	>	mol = info_->nextMolecule(i)) {
280	>
281	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
282	>	sd = mol->nextIntegrableObject(j)) {
283
284	<	oldVel[index] = integrableObject->getVel();
285	<	oldJi[index] = integrableObject->getJ();
286	<	++index;
287	<	}
284	>	oldVel[index] = sd->getVel();
285	>
286	>	if (sd->isDirectional())
287	>	oldJi[index] = sd->getJ();
288	>
289	>	++index;
290	>	}
291		}
292
293		// do the iteration:
294		instaVol =thermo.getVolume();
295
296		for(int k = 0; k < maxIterNum_; k++) {
297	<	instaTemp =thermo.getTemperature();
298	<	instaPress =thermo.getPressure();
297	>	instaTemp =thermo.getTemperature();
298	>	instaPress =thermo.getPressure();
299
300	<	// evolve chi another half step using the temperature at t + dt/2
301	<	prevChi = chi;
302	<	chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
300	>	// evolve chi another half step using the temperature at t + dt/2
301	>	prevChi = thermostat.first;
302	>	thermostat.first = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
303
304	<	//evolve eta
305	<	this->evolveEtaB();
306	<	this->calcVelScale();
304	>	//evolve eta
305	>	this->evolveEtaB();
306	>	this->calcVelScale();
307
308	<	index = 0;
309	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
310	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
299	<	integrableObject = mol->nextIntegrableObject(j)) {
308	>	index = 0;
309	>	for (mol = info_->beginMolecule(i); mol != NULL;
310	>	mol = info_->nextMolecule(i)) {
311
312	<	frc = integrableObject->getFrc();
313	<	vel = integrableObject->getVel();
312	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
313	>	sd = mol->nextIntegrableObject(j)) {
314
315	<	mass = integrableObject->getMass();
315	>	frc = sd->getFrc();
316	>	vel = sd->getVel();
317
318	<	getVelScaleB(sc, index);
318	>	mass = sd->getMass();
319
320	<	// velocity half step
309	<	//vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
310	<	vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
311	<	integrableObject->setVel(vel);
320	>	getVelScaleB(sc, index);
321
322	<	if (integrableObject->isDirectional()) {
323	<	// get and convert the torque to body frame
324	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
322	>	// velocity half step
323	>	vel = oldVel[index]
324	>	+ dt2*PhysicalConstants::energyConvert/mass* frc
325	>	- dt2*sc;
326
327	<	//ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
318	<	ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
319	<	integrableObject->setJ(ji);
320	<	}
327	>	sd->setVel(vel);
328
329	<	++index;
330	<	}
331	<	}
329	>	if (sd->isDirectional()) {
330	>	// get and convert the torque to body frame
331	>	Tb = sd->lab2Body(sd->getTrq());
332	>
333	>	ji = oldJi[index]
334	>	+ dt2*PhysicalConstants::energyConvert*Tb
335	>	- dt2*thermostat.first*oldJi[index];
336	>
337	>	sd->setJ(ji);
338	>	}
339	>
340	>	++index;
341	>	}
342	>	}
343
344	<	rattle->constraintB();
344	>	rattle_->constraintB();
345
346	<	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
347	<	break;
346	>	if ((fabs(prevChi - thermostat.first) <= chiTolerance) &&
347	>	this->etaConverged())
348	>	break;
349		}
350
351		//calculate integral of chidt
352	<	integralOfChidt += dt2 * chi;
352	>	thermostat.second += dt2 * thermostat.first;
353
354	<	currentSnapshot_->setChi(chi);
336	<	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
354	>	snap->setThermostat(thermostat);
355
356	+	flucQ_->moveB();
357		saveEta();
358	<	}
358	>	}
359
360	+	void NPT::resetIntegrator(){
361	+	snap->setThermostat(make_pair(0.0, 0.0));
362	+	resetEta();
363	+	}
364	+
365	+	void NPT::resetEta() {
366	+	Mat3x3d etaMat(0.0);
367	+	snap->setBarostat(etaMat);
368	+	}
369		}

Comparing trunk/src/integrators/NPT.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