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 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/development/src/integrators/NPT.cpp (file contents), Revision 1759 by gezelter, Thu Jun 21 18:55:33 2012 UTC

#	Line 1 | Line 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
5	+	* non-exclusive, royalty free, license to use, modify and
6	+	* redistribute this software in source and binary code form, provided
7	+	* that the following conditions are met:
8	+	*
9	+	* 1. Redistributions of source code must retain the above copyright
10	+	*    notice, this list of conditions and the following disclaimer.
11	+	*
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.
16	+	*
17	+	* This software is provided "AS IS," without a warranty of any
18	+	* kind. All express or implied conditions, representations and
19	+	* warranties, including any implied warranty of merchantability,
20	+	* fitness for a particular purpose or non-infringement, are hereby
21	+	* excluded.  The University of Notre Dame and its licensors shall not
22	+	* be liable for any damages suffered by licensee as a result of
23	+	* using, modifying or distributing the software or its
24	+	* derivatives. In no event will the University of Notre Dame or its
25	+	* licensors be liable for any lost revenue, profit or data, or for
26	+	* direct, indirect, special, consequential, incidental or punitive
27	+	* damages, however caused and regardless of the theory of liability,
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>
44
45	<	#include "Atom.hpp"
46	<	#include "SRI.hpp"
47	<	#include "AbstractClasses.hpp"
48	<	#include "SimInfo.hpp"
49	<	#include "ForceFields.hpp"
50	<	#include "Thermo.hpp"
51	<	#include "ReadWrite.hpp"
10	<	#include "Integrator.hpp"
11	<	#include "simError.h"
45	>	#include "brains/SimInfo.hpp"
46	>	#include "brains/Thermo.hpp"
47	>	#include "integrators/NPT.hpp"
48	>	#include "math/SquareMatrix3.hpp"
49	>	#include "primitives/Molecule.hpp"
50	>	#include "utils/PhysicalConstants.hpp"
51	>	#include "utils/simError.h"
52
13	–	#ifdef IS_MPI
14	–	#include "mpiSimulation.hpp"
15	–	#endif
16	–
17	–
53		// Basic isotropic thermostating and barostating via the Melchionna
54		// modification of the Hoover algorithm:
55		//
#	Line 25 | Line 60
60		//
61		//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
62
63	<	template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff):
29	<	T( theInfo, the_ff )
30	<	{
31	<	GenericData* data;
32	<	DoubleData * chiValue;
33	<	DoubleData * integralOfChidtValue;
63	>	namespace OpenMD {
64
65	<	chiValue = NULL;
66	<	integralOfChidtValue = NULL;
65	>	NPT::NPT(SimInfo* info) :
66	>	VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
67
68	<	chi = 0.0;
69	<	integralOfChidt = 0.0;
70	<	have_tau_thermostat = 0;
71	<	have_tau_barostat = 0;
72	<	have_target_temp = 0;
73	<	have_target_pressure = 0;
74	<	have_chi_tolerance = 0;
75	<	have_eta_tolerance = 0;
76	<	have_pos_iter_tolerance = 0;
68	>	Globals* simParams = info_->getSimParams();
69	>
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	>	}
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 = OPENMD_ERROR;
81	>	simError();
82	>	} else {
83	>	targetTemp = simParams->getTargetTemp();
84	>	}
85
86	<	// retrieve chi and integralOfChidt from simInfo
87	<	data = info->getProperty(CHIVALUE_ID);
88	<	if(data){
89	<	chiValue = dynamic_cast<DoubleData*>(data);
52	<	}
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");
90
91	<	data = info->getProperty(INTEGRALOFCHIDT_ID);
92	<	if(data){
93	<	integralOfChidtValue = dynamic_cast<DoubleData*>(data);
94	<	}
91	>	painCave.severity = OPENMD_ERROR;
92	>	painCave.isFatal = 1;
93	>	simError();
94	>	} else {
95	>	tauThermostat = simParams->getTauThermostat();
96	>	}
97
98	<	// chi and integralOfChidt should appear by pair
99	<	if(chiValue && integralOfChidtValue){
100	<	chi = chiValue->getData();
62	<	integralOfChidt = integralOfChidtValue->getData();
63	<	}
98	>	if (!simParams->haveTargetPressure()) {
99	>	sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
100	>	"   without a targetPressure!\n");
101
102	<	oldPos = new double[3*integrableObjects.size()];
103	<	oldVel = new double[3*integrableObjects.size()];
104	<	oldJi = new double[3*integrableObjects.size()];
102	>	painCave.isFatal = 1;
103	>	simError();
104	>	} else {
105	>	targetPressure = simParams->getTargetPressure();
106	>	}
107	>
108	>	if (!simParams->haveTauBarostat()) {
109	>	sprintf(painCave.errMsg,
110	>	"If you use the NPT integrator, you must set tauBarostat.\n");
111	>	painCave.severity = OPENMD_ERROR;
112	>	painCave.isFatal = 1;
113	>	simError();
114	>	} else {
115	>	tauBarostat = simParams->getTauBarostat();
116	>	}
117	>
118	>	tt2 = tauThermostat * tauThermostat;
119	>	tb2 = tauBarostat * tauBarostat;
120
121	<	}
121	>	updateSizes();
122	>	}
123
124	<	template<typename T> NPT<T>::~NPT() {
125	<	delete[] oldPos;
73	<	delete[] oldVel;
74	<	delete[] oldJi;
75	<	}
124	>	NPT::~NPT() {
125	>	}
126
127	<	template<typename T> void NPT<T>::moveA() {
127	>	void NPT::doUpdateSizes() {
128
129	<	//new version of NPT
130	<	int i, j, k;
131	<	double Tb[3], ji[3];
82	<	double mass;
83	<	double vel[3], pos[3], frc[3];
84	<	double sc[3];
85	<	double COM[3];
129	>	oldPos.resize(info_->getNIntegrableObjects());
130	>	oldVel.resize(info_->getNIntegrableObjects());
131	>	oldJi.resize(info_->getNIntegrableObjects());
132
133	<	instaTemp = tStats->getTemperature();
88	<	tStats->getPressureTensor( press );
89	<	instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
90	<	instaVol = tStats->getVolume();
133	>	}
134
135	<	tStats->getCOM(COM);
135	>	void NPT::moveA() {
136	>	SimInfo::MoleculeIterator i;
137	>	Molecule::IntegrableObjectIterator  j;
138	>	Molecule* mol;
139	>	StuntDouble* integrableObject;
140	>	Vector3d Tb, ji;
141	>	RealType mass;
142	>	Vector3d vel;
143	>	Vector3d pos;
144	>	Vector3d frc;
145	>	Vector3d sc;
146	>	int index;
147
148	<	//evolve velocity half step
148	>	chi= currentSnapshot_->getChi();
149	>	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
150	>	loadEta();
151	>
152	>	instaTemp =thermo.getTemperature();
153	>	press = thermo.getPressureTensor();
154	>	instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
155	>	instaVol =thermo.getVolume();
156
157	<	calcVelScale();
97	<
98	<	for( i=0; i<integrableObjects.size(); i++ ){
157	>	Vector3d  COM = info_->getCom();
158
159	<	integrableObjects[i]->getVel( vel );
101	<	integrableObjects[i]->getFrc( frc );
159	>	//evolve velocity half step
160
161	<	mass = integrableObjects[i]->getMass();
161	>	calcVelScale();
162
163	<	getVelScaleA( sc, vel );
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)) {
166	>
167	>	vel = integrableObject->getVel();
168	>	frc = integrableObject->getFrc();
169
170	<	for (j=0; j < 3; j++) {
170	>	mass = integrableObject->getMass();
171
172	<	// velocity half step  (use chi from previous step here):
110	<	vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
172	>	getVelScaleA(sc, vel);
173
174	<	}
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	<	integrableObjects[i]->setVel( vel );
179	>	if (integrableObject->isDirectional()) {
180
181	<	if( integrableObjects[i]->isDirectional() ){
181	>	// get and convert the torque to body frame
182
183	<	// get and convert the torque to body frame
183	>	Tb = integrableObject->lab2Body(integrableObject->getTrq());
184
185	<	integrableObjects[i]->getTrq( Tb );
121	<	integrableObjects[i]->lab2Body( Tb );
185	>	// get the angular momentum, and propagate a half step
186
187	<	// get the angular momentum, and propagate a half step
187	>	ji = integrableObject->getJ();
188
189	<	integrableObjects[i]->getJ( 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);
193
194	<	for (j=0; j < 3; j++)
195	<	ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
194	>	integrableObject->setJ(ji);
195	>	}
196	>
197	>	}
198	>	}
199	>	// evolve chi and eta  half step
200
201	<	this->rotationPropagation( integrableObjects[i], ji );
201	>	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
202	>
203	>	evolveEtaA();
204
205	<	integrableObjects[i]->setJ( ji );
206	<	}
207	<	}
205	>	//calculate the integral of chidt
206	>	integralOfChidt += dt2 * chi;
207	>
208	>	flucQ_->moveA();
209
136	–	// evolve chi and eta  half step
210
211	<	evolveChiA();
212	<	evolveEtaA();
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	>	}
217	>	}
218	>
219	>	//the first estimation of r(t+dt) is equal to  r(t)
220
221	<	//calculate the integral of chidt
222	<	integralOfChidt += dt2*chi;
221	>	for(int k = 0; k < maxIterNum_; k++) {
222	>	index = 0;
223	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
224	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
225	>	integrableObject = mol->nextIntegrableObject(j)) {
226
227	<	//save the old positions
228	<	for(i = 0; i < integrableObjects.size(); i++){
146	<	integrableObjects[i]->getPos(pos);
147	<	for(j = 0; j < 3; j++)
148	<	oldPos[i*3 + j] = pos[j];
149	<	}
227	>	vel = integrableObject->getVel();
228	>	pos = integrableObject->getPos();
229
230	<	//the first estimation of r(t+dt) is equal to  r(t)
230	>	this->getPosScale(pos, COM, index, sc);
231
232	<	for(k = 0; k < 5; k ++){
232	>	pos = oldPos[index] + dt * (vel + sc);
233	>	integrableObject->setPos(pos);
234
235	<	for(i =0 ; i < integrableObjects.size(); i++){
235	>	++index;
236	>	}
237	>	}
238
239	<	integrableObjects[i]->getVel(vel);
240	<	integrableObjects[i]->getPos(pos);
239	>	rattle_->constraintA();
240	>	}
241
242	<	this->getPosScale( pos, COM, i, sc );
242	>	// Scale the box after all the positions have been moved:
243
244	<	for(j = 0; j < 3; j++)
245	<	pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
244	>	this->scaleSimBox();
245	>
246	>	currentSnapshot_->setChi(chi);
247	>	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
248
249	<	integrableObjects[i]->setPos( pos );
166	<	}
167	<
168	<	if(nConstrained)
169	<	constrainA();
249	>	saveEta();
250		}
251
252	+	void NPT::moveB(void) {
253	+	SimInfo::MoleculeIterator i;
254	+	Molecule::IntegrableObjectIterator  j;
255	+	Molecule* mol;
256	+	StuntDouble* integrableObject;
257	+	int index;
258	+	Vector3d Tb;
259	+	Vector3d ji;
260	+	Vector3d sc;
261	+	Vector3d vel;
262	+	Vector3d frc;
263	+	RealType mass;
264
173	–	// Scale the box after all the positions have been moved:
265
266	<	this->scaleSimBox();
267	<	}
266	>	chi= currentSnapshot_->getChi();
267	>	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
268	>	RealType oldChi  = chi;
269	>	RealType prevChi;
270
271	<	template<typename T> void NPT<T>::moveB( void ){
271	>	loadEta();
272	>
273	>	//save velocity and angular momentum
274	>	index = 0;
275	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
276	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
277	>	integrableObject = mol->nextIntegrableObject(j)) {
278	>
279	>	oldVel[index] = integrableObject->getVel();
280
281	<	//new version of NPT
282	<	int i, j, k;
182	<	double Tb[3], ji[3], sc[3];
183	<	double vel[3], frc[3];
184	<	double mass;
281	>	if (integrableObject->isDirectional())
282	>	oldJi[index] = integrableObject->getJ();
283
284	<	// Set things up for the iteration:
285	<
188	<	for( i=0; i<integrableObjects.size(); i++ ){
189	<
190	<	integrableObjects[i]->getVel( vel );
191	<
192	<	for (j=0; j < 3; j++)
193	<	oldVel[3*i + j]  = vel[j];
194	<
195	<	if( integrableObjects[i]->isDirectional() ){
196	<
197	<	integrableObjects[i]->getJ( ji );
198	<
199	<	for (j=0; j < 3; j++)
200	<	oldJi[3*i + j] = ji[j];
201	<
284	>	++index;
285	>	}
286		}
203	–	}
287
288	<	// do the iteration:
288	>	// do the iteration:
289	>	instaVol =thermo.getVolume();
290
291	<	instaVol = tStats->getVolume();
291	>	for(int k = 0; k < maxIterNum_; k++) {
292	>	instaTemp =thermo.getTemperature();
293	>	instaPress =thermo.getPressure();
294
295	<	for (k=0; k < 4; k++) {
295	>	// evolve chi another half step using the temperature at t + dt/2
296	>	prevChi = chi;
297	>	chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
298
299	<	instaTemp = tStats->getTemperature();
300	<	instaPress = tStats->getPressure();
299	>	//evolve eta
300	>	this->evolveEtaB();
301	>	this->calcVelScale();
302
303	<	// evolve chi another half step using the temperature at t + dt/2
303	>	index = 0;
304	>	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
305	>	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
306	>	integrableObject = mol->nextIntegrableObject(j)) {
307
308	<	this->evolveChiB();
309	<	this->evolveEtaB();
218	<	this->calcVelScale();
308	>	frc = integrableObject->getFrc();
309	>	vel = integrableObject->getVel();
310
311	<	for( i=0; i<integrableObjects.size(); i++ ){
311	>	mass = integrableObject->getMass();
312
313	<	integrableObjects[i]->getFrc( frc );
223	<	integrableObjects[i]->getVel(vel);
313	>	getVelScaleB(sc, index);
314
315	<	mass = integrableObjects[i]->getMass();
315	>	// velocity half step
316	>	//vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
317	>	vel = oldVel[index] + dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
318	>	integrableObject->setVel(vel);
319
320	<	getVelScaleB( sc, i );
320	>	if (integrableObject->isDirectional()) {
321	>	// get and convert the torque to body frame
322	>	Tb = integrableObject->lab2Body(integrableObject->getTrq());
323
324	<	// velocity half step
325	<	for (j=0; j < 3; j++)
326	<	vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
324	>	//ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi[3*i+j]*chi);
325	>	ji = oldJi[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi*oldJi[index];
326	>	integrableObject->setJ(ji);
327	>	}
328
329	<	integrableObjects[i]->setVel( vel );
330	<
235	<	if( integrableObjects[i]->isDirectional() ){
236	<
237	<	// get and convert the torque to body frame
238	<
239	<	integrableObjects[i]->getTrq( Tb );
240	<	integrableObjects[i]->lab2Body( Tb );
241	<
242	<	for (j=0; j < 3; j++)
243	<	ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
244	<
245	<	integrableObjects[i]->setJ( ji );
329	>	++index;
330	>	}
331		}
332	+
333	+	rattle_->constraintB();
334	+
335	+	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
336	+	break;
337		}
338
339	<	if(nConstrained)
340	<	constrainB();
339	>	//calculate integral of chidt
340	>	integralOfChidt += dt2 * chi;
341
342	<	if ( this->chiConverged() && this->etaConverged() ) break;
343	<	}
342	>	currentSnapshot_->setChi(chi);
343	>	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
344
345	<	//calculate integral of chida
346	<	integralOfChidt += dt2*chi;
257	<
258	<
259	<	}
260	<
261	<	template<typename T> void NPT<T>::resetIntegrator() {
262	<	chi = 0.0;
263	<	T::resetIntegrator();
264	<	}
265	<
266	<	template<typename T> void NPT<T>::evolveChiA() {
267	<	chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
268	<	oldChi = chi;
269	<	}
270	<
271	<	template<typename T> void NPT<T>::evolveChiB() {
272	<
273	<	prevChi = chi;
274	<	chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
275	<	}
276	<
277	<	template<typename T> bool NPT<T>::chiConverged() {
278	<
279	<	return ( fabs( prevChi - chi ) <= chiTolerance );
280	<	}
281	<
282	<	template<typename T> int NPT<T>::readyCheck() {
283	<
284	<	//check parent's readyCheck() first
285	<	if (T::readyCheck() == -1)
286	<	return -1;
287	<
288	<	// First check to see if we have a target temperature.
289	<	// Not having one is fatal.
290	<
291	<	if (!have_target_temp) {
292	<	sprintf( painCave.errMsg,
293	<	"NPT error: You can't use the NPT integrator\n"
294	<	"   without a targetTemp!\n"
295	<	);
296	<	painCave.isFatal = 1;
297	<	simError();
298	<	return -1;
345	>	flucQ_->moveB();
346	>	saveEta();
347		}
348
349	<	if (!have_target_pressure) {
350	<	sprintf( painCave.errMsg,
351	<	"NPT error: You can't use the NPT integrator\n"
352	<	"   without a targetPressure!\n"
305	<	);
306	<	painCave.isFatal = 1;
307	<	simError();
308	<	return -1;
349	>	void NPT::resetIntegrator(){
350	>	currentSnapshot_->setChi(0.0);
351	>	currentSnapshot_->setIntegralOfChiDt(0.0);
352	>	resetEta();
353		}
354
311	–	// We must set tauThermostat.
355
356	<	if (!have_tau_thermostat) {
357	<	sprintf( painCave.errMsg,
358	<	"NPT error: If you use the NPT\n"
359	<	"   integrator, you must set tauThermostat.\n");
360	<	painCave.isFatal = 1;
318	<	simError();
319	<	return -1;
320	<	}
321	<
322	<	// We must set tauBarostat.
323	<
324	<	if (!have_tau_barostat) {
325	<	sprintf( painCave.errMsg,
326	<	"If you use the NPT integrator, you must set tauBarostat.\n");
327	<	painCave.severity = OOPSE_ERROR;
328	<	painCave.isFatal = 1;
329	<	simError();
330	<	return -1;
331	<	}
332	<
333	<	if (!have_chi_tolerance) {
334	<	sprintf( painCave.errMsg,
335	<	"Setting chi tolerance to 1e-6 in NPT integrator\n");
336	<	chiTolerance = 1e-6;
337	<	have_chi_tolerance = 1;
338	<	painCave.severity = OOPSE_INFO;
339	<	painCave.isFatal = 0;
340	<	simError();
341	<	}
342	<
343	<	if (!have_eta_tolerance) {
344	<	sprintf( painCave.errMsg,
345	<	"Setting eta tolerance to 1e-6 in NPT integrator");
346	<	etaTolerance = 1e-6;
347	<	have_eta_tolerance = 1;
348	<	painCave.severity = OOPSE_INFO;
349	<	painCave.isFatal = 0;
350	<	simError();
351	<	}
352	<
353	<	// We need NkBT a lot, so just set it here: This is the RAW number
354	<	// of integrableObjects, so no subtraction or addition of constraints or
355	<	// orientational degrees of freedom:
356	<
357	<	NkBT = (double)(info->getTotIntegrableObjects()) * kB * targetTemp;
358	<
359	<	// fkBT is used because the thermostat operates on more degrees of freedom
360	<	// than the barostat (when there are particles with orientational degrees
361	<	// of freedom).
362	<
363	<	fkBT = (double)(info->getNDF()) * kB * targetTemp;
364	<
365	<	tt2 = tauThermostat * tauThermostat;
366	<	tb2 = tauBarostat * tauBarostat;
367	<
368	<	return 1;
356	>	void NPT::resetEta() {
357	>	Mat3x3d etaMat(0.0);
358	>	currentSnapshot_->setEta(etaMat);
359	>	}
360	>
361		}

Comparing:
trunk/src/integrators/NPT.cpp (property svn:keywords), Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/development/src/integrators/NPT.cpp (property svn:keywords), Revision 1759 by gezelter, Thu Jun 21 18:55:33 2012 UTC

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

Diff Legend

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