OOPSE/libmdtools/NPTi.cpp

#include <cmath>
#include "Atom.hpp"
#include "SRI.hpp"
#include "AbstractClasses.hpp"
#include "SimInfo.hpp"
#include "ForceFields.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"
#include "Integrator.hpp"
#include "simError.h" 

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#endif


// Basic isotropic thermostating and barostating via the Melchionna
// modification of the Hoover algorithm:
//
//    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
//       Molec. Phys., 78, 533. 
//
//           and
// 
//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.

template<typename T> NPTi<T>::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  chi = 0.0;
  eta = 0.0;
  integralOfChidt = 0.0;
  have_tau_thermostat = 0;
  have_tau_barostat = 0;
  have_target_temp = 0;
  have_target_pressure = 0;
  have_chi_tolerance = 0;
  have_eta_tolerance = 0;
  have_pos_iter_tolerance = 0;

  oldPos = new double[3*nAtoms];
  oldVel = new double[3*nAtoms];
  oldJi = new double[3*nAtoms];
#ifdef IS_MPI
  Nparticles = mpiSim->getTotAtoms();
#else
  Nparticles = theInfo->n_atoms;
#endif

}

template<typename T> NPTi<T>::~NPTi() {
  delete[] oldPos;
  delete[] oldVel;
  delete[] oldJi;
}

template<typename T> void NPTi<T>::moveA() {
 
  //new version of NPTi
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double A[3][3], I[3][3];
  double angle, mass;
  double vel[3], pos[3], frc[3];

  double rj[3];
  double instaTemp, instaPress, instaVol;
  double tt2, tb2, scaleFactor;
  double COM[3];

  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  instaVol = tStats->getVolume();
  
  tStats->getCOM(COM);
  
  //evolve velocity half step
  for( i=0; i<nAtoms; i++ ){

    atoms[i]->getVel( vel );
    atoms[i]->getFrc( frc );

    mass = atoms[i]->getMass();

    for (j=0; j < 3; j++) {
      // velocity half step  (use chi from previous step here):
      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta));
  
    }

    atoms[i]->setVel( vel );
   
    if( atoms[i]->isDirectional() ){

      dAtom = (DirectionalAtom *)atoms[i];

      // get and convert the torque to body frame
      
      dAtom->getTrq( Tb );
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and propagate a half step

      dAtom->getJ( ji );

      for (j=0; j < 3; j++) 
        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
      
      // use the angular velocities to propagate the rotation matrix a
      // full time step

      dAtom->getA(A);
      dAtom->getI(I);
    
      // rotate about the x-axis      
      angle = dt2 * ji[0] / I[0][0];
      this->rotate( 1, 2, angle, ji, A ); 

      // rotate about the y-axis
      angle = dt2 * ji[1] / I[1][1];
      this->rotate( 2, 0, angle, ji, A );
      
      // rotate about the z-axis
      angle = dt * ji[2] / I[2][2];
      this->rotate( 0, 1, angle, ji, A);
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / I[1][1];
      this->rotate( 2, 0, angle, ji, A );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / I[0][0];
      this->rotate( 1, 2, angle, ji, A );
      
      dAtom->setJ( ji );
      dAtom->setA( A  );    
    }    
  }

  // evolve chi and eta  half step
  
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convert*NkBT*tb2));

  //calculate the integral of chidt
  integralOfChidt += dt2*chi;

  //save the old positions
  for(i = 0; i < nAtoms; i++){
    atoms[i]->getPos(pos);
    for(j = 0; j < 3; j++)
      oldPos[i*3 + j] = pos[j];
  }
  
  //the first estimation of r(t+dt) is equal to  r(t) 
    
  for(k = 0; k < 4; k ++){

    for(i =0 ; i < nAtoms; i++){

      atoms[i]->getVel(vel);
      atoms[i]->getPos(pos);

      for(j = 0; j < 3; j++)
        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];     
      
      for(j = 0; j < 3; j++)
        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]);

      atoms[i]->setPos( pos );
    }
    
    if (nConstrained){
      constrainA();
    }
  }
    

  // Scale the box after all the positions have been moved:
  
  scaleFactor = exp(dt*eta);

  if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
    sprintf( painCave.errMsg,
             "NPTi error: Attempting a Box scaling of more than 10 percent"
             " check your tauBarostat, as it is probably too small!\n"
             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
             );
    painCave.isFatal = 1;
    simError();
  } else {         
    info->scaleBox(scaleFactor);      
  }  

}

template<typename T> void NPTi<T>::moveB( void ){
  
  //new version of NPTi
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double vel[3], frc[3];
  double mass;

  double instTemp, instPress, instVol;
  double tt2, tb2;
  double oldChi, prevChi;
  double oldEta, preEta;
  
  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;


  // Set things up for the iteration:

  oldChi = chi;
  oldEta = eta;

  for( i=0; i<nAtoms; i++ ){

    atoms[i]->getVel( vel );

    for (j=0; j < 3; j++)
      oldVel[3*i + j]  = vel[j];

    if( atoms[i]->isDirectional() ){

      dAtom = (DirectionalAtom *)atoms[i];

      dAtom->getJ( ji );

      for (j=0; j < 3; j++)
        oldJi[3*i + j] = ji[j];

    }
  }

  // do the iteration:

  instVol = tStats->getVolume();
  
  for (k=0; k < 4; k++) {
    
    instTemp = tStats->getTemperature();
    instPress = tStats->getPressure();

    // evolve chi another half step using the temperature at t + dt/2

    prevChi = chi;
    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / 
      (tauThermostat*tauThermostat);

    preEta = eta;
    eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) / 
       (p_convert*NkBT*tb2));

  
    for( i=0; i<nAtoms; i++ ){

      atoms[i]->getFrc( frc );
      atoms[i]->getVel(vel);
      
      mass = atoms[i]->getMass();
      
      // velocity half step
      for (j=0; j < 3; j++) 
        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*(chi + eta));
      
      atoms[i]->setVel( vel );
      
      if( atoms[i]->isDirectional() ){

        dAtom = (DirectionalAtom *)atoms[i];
  
        // get and convert the torque to body frame      
  
        dAtom->getTrq( Tb );
        dAtom->lab2Body( Tb );       
             
        for (j=0; j < 3; j++) 
          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
      
          dAtom->setJ( ji );
      }
    }
    
    if (nConstrained){
      constrainB();
    }    
    
    if (fabs(prevChi - chi) <= 
        chiTolerance && fabs(preEta -eta) <= etaTolerance) 
      break;
  }

  //calculate integral of chida
  integralOfChidt += dt2*chi;


}

template<typename T> void NPTi<T>::resetIntegrator() {
  chi = 0.0;
  eta = 0.0;
}

template<typename T> int NPTi<T>::readyCheck() {

  //check parent's readyCheck() first
  if (T::readyCheck() == -1)
    return -1;
 
  // First check to see if we have a target temperature. 
  // Not having one is fatal. 
  
  if (!have_target_temp) {
    sprintf( painCave.errMsg,
             "NPTi error: You can't use the NPTi integrator\n"
             "   without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

  if (!have_target_pressure) {
    sprintf( painCave.errMsg,
             "NPTi error: You can't use the NPTi integrator\n"
             "   without a targetPressure!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }
  
  // We must set tauThermostat.
   
  if (!have_tau_thermostat) {
    sprintf( painCave.errMsg,
             "NPTi error: If you use the NPTi\n"
             "   integrator, you must set tauThermostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }    

  // We must set tauBarostat.
   
  if (!have_tau_barostat) {
    sprintf( painCave.errMsg,
             "NPTi error: If you use the NPTi\n"
             "   integrator, you must set tauBarostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }    

  if (!have_chi_tolerance) {
    sprintf( painCave.errMsg,
             "NPTi warning: setting chi tolerance to 1e-6\n");
    chiTolerance = 1e-6;
    have_chi_tolerance = 1;
    painCave.isFatal = 0;
    simError();
  } 

    if (!have_eta_tolerance) {
    sprintf( painCave.errMsg,
             "NPTi warning: setting eta tolerance to 1e-6\n");
    etaTolerance = 1e-6;
    have_eta_tolerance = 1;
    painCave.isFatal = 0;
    simError();
  } 
  // We need NkBT a lot, so just set it here:

  NkBT = (double)Nparticles * kB * targetTemp;
  fkBT = (double)info->ndf * kB * targetTemp;

  return 1;
}

template<typename T> double NPTi<T>::getConservedQuantity(void){

  double conservedQuantity;
  double LkBT;
  double fkBT;
  double f1kBT;
  double f2kBT;
  double NkBT;
  double Energy;
  double thermostat_kinetic;
  double thermostat_potential;
  double barostat_kinetic;
  double barostat_potential;
  double tb2;
  double eta2;  

  LkBT = (double)(info->getNDF() + 4) * kB * targetTemp;  // 3N + 1
  fkBT = (double)(info->getNDF()    ) * kB * targetTemp;  // 3N - 3
  f1kBT = (double)(info->getNDF()+ 1) * kB * targetTemp;  // 3N - 3 + 1
  NkBT = (double)(info->getNDF() + 3) * kB * targetTemp;  // 3N
  f2kBT = (double)(info->getNDF()+ 2) * kB * targetTemp;  // 3N - 3 + 1

  Energy = tStats->getTotalE();

  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / 
    (2.0 * eConvert);

  thermostat_potential = fkBT* integralOfChidt / eConvert;


  barostat_kinetic = fkBT * tauBarostat * tauBarostat * eta * eta / 
    (2.0 * eConvert);
  
  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) / 
    eConvert;

  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
    barostat_kinetic + barostat_potential;
  
  cout.width(8);
  cout.precision(8);

  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << 
      "\t" << thermostat_potential << "\t" << barostat_kinetic << 
      "\t" << barostat_potential << "\t" << conservedQuantity << endl;

  return conservedQuantity; 
}
Revision:	769
Committed:	Fri Sep 19 14:22:29 2003 UTC (20 years, 9 months ago) by tim
File size:	10274 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	gezelter	578	#include <cmath>
2	gezelter	574	#include "Atom.hpp"
3			#include "SRI.hpp"
4			#include "AbstractClasses.hpp"
5			#include "SimInfo.hpp"
6			#include "ForceFields.hpp"
7			#include "Thermo.hpp"
8			#include "ReadWrite.hpp"
9			#include "Integrator.hpp"
10			#include "simError.h"
11
12	tim	763	#ifdef IS_MPI
13			#include "mpiSimulation.hpp"
14			#endif
15	gezelter	574
16	tim	763
17	gezelter	574	// Basic isotropic thermostating and barostating via the Melchionna
18			// modification of the Hoover algorithm:
19			//
20			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
21			// Molec. Phys., 78, 533.
22			//
23			// and
24			//
25			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
26
27	tim	645	template<typename T> NPTi<T>::NPTi ( SimInfo theInfo, ForceFields the_ff):
28			T( theInfo, the_ff )
29	gezelter	574	{
30			chi = 0.0;
31			eta = 0.0;
32	tim	763	integralOfChidt = 0.0;
33	gezelter	574	have_tau_thermostat = 0;
34			have_tau_barostat = 0;
35			have_target_temp = 0;
36			have_target_pressure = 0;
37	tim	763	have_chi_tolerance = 0;
38			have_eta_tolerance = 0;
39			have_pos_iter_tolerance = 0;
40
41			oldPos = new double[3*nAtoms];
42			oldVel = new double[3*nAtoms];
43			oldJi = new double[3*nAtoms];
44			#ifdef IS_MPI
45			Nparticles = mpiSim->getTotAtoms();
46			#else
47			Nparticles = theInfo->n_atoms;
48			#endif
49
50	gezelter	574	}
51
52	tim	763	template<typename T> NPTi<T>::~NPTi() {
53			delete[] oldPos;
54			delete[] oldVel;
55			delete[] oldJi;
56			}
57
58	tim	645	template<typename T> void NPTi<T>::moveA() {
59	tim	763
60			//new version of NPTi
61			int i, j, k;
62	gezelter	574	DirectionalAtom* dAtom;
63	gezelter	600	double Tb[3], ji[3];
64			double A[3][3], I[3][3];
65			double angle, mass;
66			double vel[3], pos[3], frc[3];
67
68	gezelter	574	double rj[3];
69			double instaTemp, instaPress, instaVol;
70	gezelter	611	double tt2, tb2, scaleFactor;
71	tim	763	double COM[3];
72	gezelter	574
73			tt2 = tauThermostat * tauThermostat;
74			tb2 = tauBarostat * tauBarostat;
75
76			instaTemp = tStats->getTemperature();
77			instaPress = tStats->getPressure();
78			instaVol = tStats->getVolume();
79
80	tim	763	tStats->getCOM(COM);
81
82			//evolve velocity half step
83			for( i=0; i<nAtoms; i++ ){
84	gezelter	574
85	gezelter	600	atoms[i]->getVel( vel );
86			atoms[i]->getFrc( frc );
87	gezelter	574
88	gezelter	600	mass = atoms[i]->getMass();
89	gezelter	574
90	gezelter	600	for (j=0; j < 3; j++) {
91	tim	763	// velocity half step (use chi from previous step here):
92			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta));
93
94	gezelter	600	}
95
96			atoms[i]->setVel( vel );
97	tim	763
98	gezelter	574	if( atoms[i]->isDirectional() ){
99
100			dAtom = (DirectionalAtom *)atoms[i];
101	tim	763
102	gezelter	574	// get and convert the torque to body frame
103
104	gezelter	600	dAtom->getTrq( Tb );
105	gezelter	574	dAtom->lab2Body( Tb );
106
107			// get the angular momentum, and propagate a half step
108
109	gezelter	600	dAtom->getJ( ji );
110
111			for (j=0; j < 3; j++)
112			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
113	gezelter	574
114			// use the angular velocities to propagate the rotation matrix a
115			// full time step
116	gezelter	600
117			dAtom->getA(A);
118			dAtom->getI(I);
119
120	gezelter	574	// rotate about the x-axis
121	gezelter	600	angle = dt2 * ji[0] / I[0][0];
122			this->rotate( 1, 2, angle, ji, A );
123
124	gezelter	574	// rotate about the y-axis
125	gezelter	600	angle = dt2 * ji[1] / I[1][1];
126			this->rotate( 2, 0, angle, ji, A );
127	gezelter	574
128			// rotate about the z-axis
129	gezelter	600	angle = dt * ji[2] / I[2][2];
130			this->rotate( 0, 1, angle, ji, A);
131	gezelter	574
132			// rotate about the y-axis
133	gezelter	600	angle = dt2 * ji[1] / I[1][1];
134			this->rotate( 2, 0, angle, ji, A );
135	gezelter	574
136			// rotate about the x-axis
137	gezelter	600	angle = dt2 * ji[0] / I[0][0];
138			this->rotate( 1, 2, angle, ji, A );
139	gezelter	574
140	gezelter	600	dAtom->setJ( ji );
141			dAtom->setA( A );
142	tim	763	}
143			}
144	gezelter	600
145	tim	763	// evolve chi and eta half step
146
147			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
148			eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convertNkBTtb2));
149
150			//calculate the integral of chidt
151			integralOfChidt += dt2*chi;
152
153			//save the old positions
154			for(i = 0; i < nAtoms; i++){
155			atoms[i]->getPos(pos);
156			for(j = 0; j < 3; j++)
157			oldPos[i*3 + j] = pos[j];
158	gezelter	574	}
159	tim	763
160			//the first estimation of r(t+dt) is equal to r(t)
161
162			for(k = 0; k < 4; k ++){
163	gezelter	611
164	tim	763	for(i =0 ; i < nAtoms; i++){
165
166			atoms[i]->getVel(vel);
167			atoms[i]->getPos(pos);
168
169			for(j = 0; j < 3; j++)
170	tim	767	rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
171	tim	763
172			for(j = 0; j < 3; j++)
173			pos[j] = oldPos[i3 + j] + dt(vel[j] + eta*rj[j]);
174
175			atoms[i]->setPos( pos );
176			}
177	mmeineke	768
178			if (nConstrained){
179			constrainA();
180			}
181	tim	763	}
182
183
184	gezelter	577	// Scale the box after all the positions have been moved:
185	gezelter	600
186	gezelter	611	scaleFactor = exp(dt*eta);
187
188	mmeineke	614	if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
189	gezelter	611	sprintf( painCave.errMsg,
190			"NPTi error: Attempting a Box scaling of more than 10 percent"
191			" check your tauBarostat, as it is probably too small!\n"
192			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
193			);
194			painCave.isFatal = 1;
195			simError();
196			} else {
197	tim	763	info->scaleBox(scaleFactor);
198			}
199	mmeineke	614
200	gezelter	574	}
201
202	tim	645	template<typename T> void NPTi<T>::moveB( void ){
203	gezelter	574
204	tim	763	//new version of NPTi
205			int i, j, k;
206			DirectionalAtom* dAtom;
207			double Tb[3], ji[3];
208			double vel[3], frc[3];
209			double mass;
210
211			double instTemp, instPress, instVol;
212			double tt2, tb2;
213			double oldChi, prevChi;
214			double oldEta, preEta;
215
216			tt2 = tauThermostat * tauThermostat;
217			tb2 = tauBarostat * tauBarostat;
218
219
220			// Set things up for the iteration:
221
222			oldChi = chi;
223			oldEta = eta;
224
225			for( i=0; i<nAtoms; i++ ){
226
227			atoms[i]->getVel( vel );
228
229			for (j=0; j < 3; j++)
230			oldVel[3*i + j] = vel[j];
231
232			if( atoms[i]->isDirectional() ){
233
234			dAtom = (DirectionalAtom *)atoms[i];
235
236			dAtom->getJ( ji );
237
238			for (j=0; j < 3; j++)
239			oldJi[3*i + j] = ji[j];
240
241			}
242			}
243
244			// do the iteration:
245
246			instVol = tStats->getVolume();
247
248			for (k=0; k < 4; k++) {
249
250			instTemp = tStats->getTemperature();
251			instPress = tStats->getPressure();
252
253			// evolve chi another half step using the temperature at t + dt/2
254
255			prevChi = chi;
256			chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
257			(tauThermostat*tauThermostat);
258
259			preEta = eta;
260			eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
261			(p_convertNkBTtb2));
262
263
264			for( i=0; i<nAtoms; i++ ){
265
266			atoms[i]->getFrc( frc );
267			atoms[i]->getVel(vel);
268
269			mass = atoms[i]->getMass();
270
271			// velocity half step
272			for (j=0; j < 3; j++)
273			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j](chi + eta));
274
275			atoms[i]->setVel( vel );
276
277			if( atoms[i]->isDirectional() ){
278
279			dAtom = (DirectionalAtom *)atoms[i];
280
281			// get and convert the torque to body frame
282
283			dAtom->getTrq( Tb );
284			dAtom->lab2Body( Tb );
285
286			for (j=0; j < 3; j++)
287			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
288
289			dAtom->setJ( ji );
290			}
291			}
292	mmeineke	768
293			if (nConstrained){
294			constrainB();
295			}
296
297			if (fabs(prevChi - chi) <=
298			chiTolerance && fabs(preEta -eta) <= etaTolerance)
299	tim	763	break;
300			}
301
302			//calculate integral of chida
303			integralOfChidt += dt2*chi;
304
305
306	gezelter	574	}
307
308	mmeineke	746	template<typename T> void NPTi<T>::resetIntegrator() {
309			chi = 0.0;
310			eta = 0.0;
311			}
312
313	tim	645	template<typename T> int NPTi<T>::readyCheck() {
314	tim	658
315			//check parent's readyCheck() first
316			if (T::readyCheck() == -1)
317			return -1;
318	gezelter	574
319			// First check to see if we have a target temperature.
320			// Not having one is fatal.
321
322			if (!have_target_temp) {
323			sprintf( painCave.errMsg,
324			"NPTi error: You can't use the NPTi integrator\n"
325			" without a targetTemp!\n"
326			);
327			painCave.isFatal = 1;
328			simError();
329			return -1;
330			}
331
332			if (!have_target_pressure) {
333			sprintf( painCave.errMsg,
334			"NPTi error: You can't use the NPTi integrator\n"
335			" without a targetPressure!\n"
336			);
337			painCave.isFatal = 1;
338			simError();
339			return -1;
340			}
341
342			// We must set tauThermostat.
343
344			if (!have_tau_thermostat) {
345			sprintf( painCave.errMsg,
346			"NPTi error: If you use the NPTi\n"
347			" integrator, you must set tauThermostat.\n");
348			painCave.isFatal = 1;
349			simError();
350			return -1;
351			}
352
353			// We must set tauBarostat.
354
355			if (!have_tau_barostat) {
356			sprintf( painCave.errMsg,
357			"NPTi error: If you use the NPTi\n"
358			" integrator, you must set tauBarostat.\n");
359			painCave.isFatal = 1;
360			simError();
361			return -1;
362			}
363
364	tim	763	if (!have_chi_tolerance) {
365			sprintf( painCave.errMsg,
366			"NPTi warning: setting chi tolerance to 1e-6\n");
367			chiTolerance = 1e-6;
368			have_chi_tolerance = 1;
369			painCave.isFatal = 0;
370			simError();
371			}
372
373			if (!have_eta_tolerance) {
374			sprintf( painCave.errMsg,
375			"NPTi warning: setting eta tolerance to 1e-6\n");
376			etaTolerance = 1e-6;
377			have_eta_tolerance = 1;
378			painCave.isFatal = 0;
379			simError();
380			}
381	gezelter	574	// We need NkBT a lot, so just set it here:
382
383	tim	763	NkBT = (double)Nparticles * kB * targetTemp;
384			fkBT = (double)info->ndf * kB * targetTemp;
385	gezelter	574
386			return 1;
387			}
388	tim	763
389			template<typename T> double NPTi<T>::getConservedQuantity(void){
390
391			double conservedQuantity;
392	tim	769	double LkBT;
393			double fkBT;
394			double f1kBT;
395			double f2kBT;
396			double NkBT;
397			double Energy;
398			double thermostat_kinetic;
399			double thermostat_potential;
400			double barostat_kinetic;
401			double barostat_potential;
402	tim	763	double tb2;
403			double eta2;
404
405	tim	769	LkBT = (double)(info->getNDF() + 4) * kB * targetTemp; // 3N + 1
406			fkBT = (double)(info->getNDF() ) * kB * targetTemp; // 3N - 3
407			f1kBT = (double)(info->getNDF()+ 1) * kB * targetTemp; // 3N - 3 + 1
408			NkBT = (double)(info->getNDF() + 3) * kB * targetTemp; // 3N
409			f2kBT = (double)(info->getNDF()+ 2) * kB * targetTemp; // 3N - 3 + 1
410	tim	763
411	tim	769	Energy = tStats->getTotalE();
412	tim	763
413	tim	769	thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
414			(2.0 * eConvert);
415	tim	763
416	tim	769	thermostat_potential = fkBT* integralOfChidt / eConvert;
417	tim	763
418
419	tim	769	barostat_kinetic = fkBT * tauBarostat * tauBarostat * eta * eta /
420			(2.0 * eConvert);
421
422			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
423			eConvert;
424	tim	767
425	tim	769	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
426			barostat_kinetic + barostat_potential;
427
428	tim	763	cout.width(8);
429			cout.precision(8);
430
431	tim	769	cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
432			"\t" << thermostat_potential << "\t" << barostat_kinetic <<
433			"\t" << barostat_potential << "\t" << conservedQuantity << endl;
434	tim	763
435			return conservedQuantity;
436			}