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 
      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  );    
    }    
  }

  // advance chi half step
  
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;

  // calculate the integral of chidt

  integralOfChidt += dt2*chi;

  // advance eta half step

  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convert*NkBT*tb2));

  //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 instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double oldChi, prevChi;
  double oldEta, prevEta;
  
  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:

  instaVol = tStats->getVolume();
  
  for (k=0; k < 4; k++) {
    
    instaTemp = tStats->getTemperature();
    instaPress = tStats->getPressure();

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

    prevChi = chi;
    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;

    prevEta = eta;

    // advance eta half step and calculate scale factor for velocity

    eta = oldEta + dt2 * ( instaVol * (instaPress - 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(prevEta -eta) <= etaTolerance) 
      break;
  }

  //calculate integral of chidt
  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: This is the RAW number
  // of particles, so no subtraction or addition of constraints or
  // orientational degrees of freedom:
  
  NkBT = (double)Nparticles * kB * targetTemp;
  
  // fkBT is used because the thermostat operates on more degrees of freedom
  // than the barostat (when there are particles with orientational degrees
  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
  
  fkBT = (double)info->ndf * kB * targetTemp;

  return 1;
}

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

  double conservedQuantity;
  double Three_NkBT;
  double Energy;
  double thermostat_kinetic;
  double thermostat_potential;
  double barostat_kinetic;
  double barostat_potential;
  double tb2;
  double eta2; 

  Energy = tStats->getTotalE();

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

  thermostat_potential = fkBT* integralOfChidt / eConvert;


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