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 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

#	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, 234107 (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) :
66	<	VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
66	>	VelocityVerletIntegrator(info), etaTolerance(1e-6), chiTolerance(1e-6),
67	>	maxIterNum_(4) {
68
69		Globals* simParams = info_->getSimParams();
70
71		if (!simParams->getUseIntialExtendedSystemState()) {
72		Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
73	<	currSnapshot->setChi(0.0);
74	<	currSnapshot->setIntegralOfChiDt(0.0);
73	<	currSnapshot->setEta(Mat3x3d(0.0));
73	>	currSnapshot->setThermostat(make_pair(0.0, 0.0));
74	>	currSnapshot->setBarostat(Mat3x3d(0.0));
75		}
76
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 {
83		targetTemp = simParams->getTargetTemp();
#	Line 87 | Line 88 | namespace oopse {
88		sprintf(painCave.errMsg, "If you use the constant temperature\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 {
#	Line 107 | Line 108 | namespace oopse {
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 {
#	Line 117 | Line 118 | namespace oopse {
118		tt2 = tauThermostat * tauThermostat;
119		tb2 = tauBarostat * tauBarostat;
120
121	<	update();
121	>	updateSizes();
122		}
123
124		NPT::~NPT() {
125		}
126
127	<	void NPT::doUpdate() {
127	>	void NPT::doUpdateSizes() {
128
129		oldPos.resize(info_->getNIntegrableObjects());
130		oldVel.resize(info_->getNIntegrableObjects());
#	Line 135 | Line 136 | namespace oopse {
136		SimInfo::MoleculeIterator i;
137		Molecule::IntegrableObjectIterator  j;
138		Molecule* mol;
139	<	StuntDouble* integrableObject;
139	>	StuntDouble* sd;
140		Vector3d Tb, ji;
141		RealType mass;
142		Vector3d vel;
#	Line 144 | Line 145 | namespace oopse {
145		Vector3d sc;
146		int index;
147
148	<	chi= currentSnapshot_->getChi();
148	<	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
148	>	thermostat = snap->getThermostat();
149		loadEta();
150
151		instaTemp =thermo.getTemperature();
152		press = thermo.getPressureTensor();
153	<	instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
153	>	instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
154		instaVol =thermo.getVolume();
155
156	<	Vector3d  COM = info_->getCom();
156	>	Vector3d  COM = thermo.getCom();
157
158		//evolve velocity half step
159
160		calcVelScale();
161
162	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
163	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
164	<	integrableObject = mol->nextIntegrableObject(j)) {
162	>	for (mol = info_->beginMolecule(i); mol != NULL;
163	>	mol = info_->nextMolecule(i)) {
164	>
165	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
166	>	sd = mol->nextIntegrableObject(j)) {
167
168	<	vel = integrableObject->getVel();
169	<	frc = integrableObject->getFrc();
168	>	vel = sd->getVel();
169	>	frc = sd->getFrc();
170
171	<	mass = integrableObject->getMass();
171	>	mass = sd->getMass();
172
173		getVelScaleA(sc, vel);
174
175		// 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);
176
177	<	if (integrableObject->isDirectional()) {
177	>	vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
178	>	sd->setVel(vel);
179
180	+	if (sd->isDirectional()) {
181	+
182		// get and convert the torque to body frame
183
184	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
184	>	Tb = sd->lab2Body(sd->getTrq());
185
186		// get the angular momentum, and propagate a half step
187
188	<	ji = integrableObject->getJ();
188	>	ji = sd->getJ();
189
190	<	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
191	<	ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
190	>	ji += dt2*PhysicalConstants::energyConvert * Tb
191	>	- dt2*thermostat.first* ji;
192
193	<	rotAlgo->rotate(integrableObject, ji, dt);
193	>	rotAlgo_->rotate(sd, ji, dt);
194
195	<	integrableObject->setJ(ji);
195	>	sd->setJ(ji);
196		}
197
198		}
199		}
200		// evolve chi and eta  half step
201
202	<	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
202	>	thermostat.first += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
203
204		evolveEtaA();
205
206		//calculate the integral of chidt
207	<	integralOfChidt += dt2 * chi;
207	>	thermostat.second += dt2 * thermostat.first;
208
209	+	flucQ_->moveA();
210	+
211	+
212		index = 0;
213	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
214	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
215	<	integrableObject = mol->nextIntegrableObject(j)) {
216	<	oldPos[index++] = integrableObject->getPos();
213	>	for (mol = info_->beginMolecule(i); mol != NULL;
214	>	mol = info_->nextMolecule(i)) {
215	>
216	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
217	>	sd = mol->nextIntegrableObject(j)) {
218	>
219	>	oldPos[index++] = sd->getPos();
220	>
221		}
222		}
223
#	Line 216 | Line 225 | namespace oopse {
225
226		for(int k = 0; k < maxIterNum_; k++) {
227		index = 0;
228	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
229	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
221	<	integrableObject = mol->nextIntegrableObject(j)) {
228	>	for (mol = info_->beginMolecule(i); mol != NULL;
229	>	mol = info_->nextMolecule(i)) {
230
231	<	vel = integrableObject->getVel();
232	<	pos = integrableObject->getPos();
231	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
232	>	sd = mol->nextIntegrableObject(j)) {
233
234	+	vel = sd->getVel();
235	+	pos = sd->getPos();
236	+
237		this->getPosScale(pos, COM, index, sc);
238
239		pos = oldPos[index] + dt * (vel + sc);
240	<	integrableObject->setPos(pos);
240	>	sd->setPos(pos);
241
242		++index;
243		}
244		}
245
246	<	rattle->constraintA();
246	>	rattle_->constraintA();
247		}
248
249		// Scale the box after all the positions have been moved:
250
251		this->scaleSimBox();
252
253	<	currentSnapshot_->setChi(chi);
243	<	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
253	>	snap->setThermostat(thermostat);
254
255		saveEta();
256		}
#	Line 249 | Line 259 | namespace oopse {
259		SimInfo::MoleculeIterator i;
260		Molecule::IntegrableObjectIterator  j;
261		Molecule* mol;
262	<	StuntDouble* integrableObject;
262	>	StuntDouble* sd;
263		int index;
264		Vector3d Tb;
265		Vector3d ji;
#	Line 258 | Line 268 | namespace oopse {
268		Vector3d frc;
269		RealType mass;
270
271	<
272	<	chi= currentSnapshot_->getChi();
263	<	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
264	<	RealType oldChi  = chi;
271	>	thermostat = snap->getThermostat();
272	>	RealType oldChi  = thermostat.first;
273		RealType prevChi;
274
275		loadEta();
276
277		//save velocity and angular momentum
278		index = 0;
279	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
280	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
281	<	integrableObject = mol->nextIntegrableObject(j)) {
279	>	for (mol = info_->beginMolecule(i); mol != NULL;
280	>	mol = info_->nextMolecule(i)) {
281	>
282	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
283	>	sd = mol->nextIntegrableObject(j)) {
284
285	<	oldVel[index] = integrableObject->getVel();
286	<	oldJi[index] = integrableObject->getJ();
285	>	oldVel[index] = sd->getVel();
286	>
287	>	if (sd->isDirectional())
288	>	oldJi[index] = sd->getJ();
289	>
290		++index;
291		}
292		}
#	Line 286 | Line 299 | namespace oopse {
299		instaPress =thermo.getPressure();
300
301		// evolve chi another half step using the temperature at t + dt/2
302	<	prevChi = chi;
303	<	chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
302	>	prevChi = thermostat.first;
303	>	thermostat.first = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
304
305		//evolve eta
306		this->evolveEtaB();
307		this->calcVelScale();
308
309		index = 0;
310	<	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
311	<	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
299	<	integrableObject = mol->nextIntegrableObject(j)) {
310	>	for (mol = info_->beginMolecule(i); mol != NULL;
311	>	mol = info_->nextMolecule(i)) {
312
313	<	frc = integrableObject->getFrc();
314	<	vel = integrableObject->getVel();
313	>	for (sd = mol->beginIntegrableObject(j); sd != NULL;
314	>	sd = mol->nextIntegrableObject(j)) {
315
316	<	mass = integrableObject->getMass();
316	>	frc = sd->getFrc();
317	>	mass = sd->getMass();
318
319		getVelScaleB(sc, index);
320
321		// velocity half step
322	<	//vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
323	<	vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
324	<	integrableObject->setVel(vel);
322	>	vel = oldVel[index]
323	>	+ dt2*PhysicalConstants::energyConvert/mass* frc
324	>	- dt2*sc;
325
326	<	if (integrableObject->isDirectional()) {
326	>	sd->setVel(vel);
327	>
328	>	if (sd->isDirectional()) {
329		// get and convert the torque to body frame
330	<	Tb = integrableObject->lab2Body(integrableObject->getTrq());
330	>	Tb = sd->lab2Body(sd->getTrq());
331
332	<	//ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
333	<	ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
334	<	integrableObject->setJ(ji);
332	>	ji = oldJi[index]
333	>	+ dt2*PhysicalConstants::energyConvert*Tb
334	>	- dt2*thermostat.first*oldJi[index];
335	>
336	>	sd->setJ(ji);
337		}
338
339		++index;
340		}
341		}
342
343	<	rattle->constraintB();
343	>	rattle_->constraintB();
344
345	<	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
345	>	if ((fabs(prevChi - thermostat.first) <= chiTolerance) &&
346	>	this->etaConverged())
347		break;
348		}
349
350		//calculate integral of chidt
351	<	integralOfChidt += dt2 * chi;
351	>	thermostat.second += dt2 * thermostat.first;
352
353	<	currentSnapshot_->setChi(chi);
336	<	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
353	>	snap->setThermostat(thermostat);
354
355	+	flucQ_->moveB();
356		saveEta();
357		}
358
359		void NPT::resetIntegrator(){
360	<	currentSnapshot_->setChi(0.0);
361	<	currentSnapshot_->setIntegralOfChiDt(0.0);
344	<	resetEta();
360	>	snap->setThermostat(make_pair(0.0, 0.0));
361	>	resetEta();
362		}
363
364	<
365	<	void NPT::resetEta() {
366	<	Mat3x3d etaMat(0.0);
367	<	currentSnapshot_->setEta(etaMat);
351	<	}
352	<
364	>	void NPT::resetEta() {
365	>	Mat3x3d etaMat(0.0);
366	>	snap->setBarostat(etaMat);
367	>	}
368		}

Comparing trunk/src/integrators/NPT.cpp (property svn:keywords):
Revision 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

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

Diff Legend

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