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

    }

  }
    

  // 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 (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 tb2;
  double eta2;  
  double E_NPT;
  double U;
  double TS;
  double PV;
  double extra;

  U = tStats->getTotalE();

  TS = fkBT * 
    (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;

  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;

  tb2 = tauBarostat * tauBarostat;
  eta2 = eta * eta;


  extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / eConvert;

  cout.width(8);
  cout.precision(8);

  
  cout << info->getTime() << "\t" 
       << chi << "\t"
       << eta << "\t"
       << U << "\t" 
       << TS << "\t"
       << PV << "\t"
       << extra << "\t"
       << U+TS+PV+extra << endl;

  conservedQuantity = U+TS+PV+extra;
  return conservedQuantity; 
}
Revision:	767
Committed:	Tue Sep 16 20:02:11 2003 UTC (20 years, 9 months ago) by tim
File size:	9539 byte(s)
Log Message:	fixed ecr grow in SimInfo fixed conserved quantity in NPT (Still some small bug) NPTi appears very stable.
#	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			}
178
179			}
180
181
182	gezelter	577	// Scale the box after all the positions have been moved:
183	gezelter	600
184	gezelter	611	scaleFactor = exp(dt*eta);
185
186	mmeineke	614	if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
187	gezelter	611	sprintf( painCave.errMsg,
188			"NPTi error: Attempting a Box scaling of more than 10 percent"
189			" check your tauBarostat, as it is probably too small!\n"
190			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
191			);
192			painCave.isFatal = 1;
193			simError();
194			} else {
195	tim	763	info->scaleBox(scaleFactor);
196			}
197	mmeineke	614
198	gezelter	574	}
199
200	tim	645	template<typename T> void NPTi<T>::moveB( void ){
201	gezelter	574
202	tim	763	//new version of NPTi
203			int i, j, k;
204			DirectionalAtom* dAtom;
205			double Tb[3], ji[3];
206			double vel[3], frc[3];
207			double mass;
208
209			double instTemp, instPress, instVol;
210			double tt2, tb2;
211			double oldChi, prevChi;
212			double oldEta, preEta;
213
214			tt2 = tauThermostat * tauThermostat;
215			tb2 = tauBarostat * tauBarostat;
216
217
218			// Set things up for the iteration:
219
220			oldChi = chi;
221			oldEta = eta;
222
223			for( i=0; i<nAtoms; i++ ){
224
225			atoms[i]->getVel( vel );
226
227			for (j=0; j < 3; j++)
228			oldVel[3*i + j] = vel[j];
229
230			if( atoms[i]->isDirectional() ){
231
232			dAtom = (DirectionalAtom *)atoms[i];
233
234			dAtom->getJ( ji );
235
236			for (j=0; j < 3; j++)
237			oldJi[3*i + j] = ji[j];
238
239			}
240			}
241
242			// do the iteration:
243
244			instVol = tStats->getVolume();
245
246			for (k=0; k < 4; k++) {
247
248			instTemp = tStats->getTemperature();
249			instPress = tStats->getPressure();
250
251			// evolve chi another half step using the temperature at t + dt/2
252
253			prevChi = chi;
254			chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
255			(tauThermostat*tauThermostat);
256
257			preEta = eta;
258			eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
259			(p_convertNkBTtb2));
260
261
262			for( i=0; i<nAtoms; i++ ){
263
264			atoms[i]->getFrc( frc );
265			atoms[i]->getVel(vel);
266
267			mass = atoms[i]->getMass();
268
269			// velocity half step
270			for (j=0; j < 3; j++)
271			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j](chi + eta));
272
273			atoms[i]->setVel( vel );
274
275			if( atoms[i]->isDirectional() ){
276
277			dAtom = (DirectionalAtom *)atoms[i];
278
279			// get and convert the torque to body frame
280
281			dAtom->getTrq( Tb );
282			dAtom->lab2Body( Tb );
283
284			for (j=0; j < 3; j++)
285			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
286
287			dAtom->setJ( ji );
288			}
289			}
290
291			if (fabs(prevChi - chi) <= chiTolerance && fabs(preEta -eta) <= etaTolerance)
292			break;
293			}
294
295			//calculate integral of chida
296			integralOfChidt += dt2*chi;
297
298
299	gezelter	574	}
300
301	mmeineke	746	template<typename T> void NPTi<T>::resetIntegrator() {
302			chi = 0.0;
303			eta = 0.0;
304			}
305
306	tim	645	template<typename T> int NPTi<T>::readyCheck() {
307	tim	658
308			//check parent's readyCheck() first
309			if (T::readyCheck() == -1)
310			return -1;
311	gezelter	574
312			// First check to see if we have a target temperature.
313			// Not having one is fatal.
314
315			if (!have_target_temp) {
316			sprintf( painCave.errMsg,
317			"NPTi error: You can't use the NPTi integrator\n"
318			" without a targetTemp!\n"
319			);
320			painCave.isFatal = 1;
321			simError();
322			return -1;
323			}
324
325			if (!have_target_pressure) {
326			sprintf( painCave.errMsg,
327			"NPTi error: You can't use the NPTi integrator\n"
328			" without a targetPressure!\n"
329			);
330			painCave.isFatal = 1;
331			simError();
332			return -1;
333			}
334
335			// We must set tauThermostat.
336
337			if (!have_tau_thermostat) {
338			sprintf( painCave.errMsg,
339			"NPTi error: If you use the NPTi\n"
340			" integrator, you must set tauThermostat.\n");
341			painCave.isFatal = 1;
342			simError();
343			return -1;
344			}
345
346			// We must set tauBarostat.
347
348			if (!have_tau_barostat) {
349			sprintf( painCave.errMsg,
350			"NPTi error: If you use the NPTi\n"
351			" integrator, you must set tauBarostat.\n");
352			painCave.isFatal = 1;
353			simError();
354			return -1;
355			}
356
357	tim	763	if (!have_chi_tolerance) {
358			sprintf( painCave.errMsg,
359			"NPTi warning: setting chi tolerance to 1e-6\n");
360			chiTolerance = 1e-6;
361			have_chi_tolerance = 1;
362			painCave.isFatal = 0;
363			simError();
364			}
365
366			if (!have_eta_tolerance) {
367			sprintf( painCave.errMsg,
368			"NPTi warning: setting eta tolerance to 1e-6\n");
369			etaTolerance = 1e-6;
370			have_eta_tolerance = 1;
371			painCave.isFatal = 0;
372			simError();
373			}
374	gezelter	574	// We need NkBT a lot, so just set it here:
375
376	tim	763	NkBT = (double)Nparticles * kB * targetTemp;
377			fkBT = (double)info->ndf * kB * targetTemp;
378	gezelter	574
379			return 1;
380			}
381	tim	763
382			template<typename T> double NPTi<T>::getConservedQuantity(void){
383
384			double conservedQuantity;
385			double tb2;
386			double eta2;
387			double E_NPT;
388			double U;
389			double TS;
390			double PV;
391			double extra;
392
393			U = tStats->getTotalE();
394
395			TS = fkBT *
396			(integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
397
398			PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
399
400			tb2 = tauBarostat * tauBarostat;
401			eta2 = eta * eta;
402
403
404	tim	767	extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / eConvert;
405
406	tim	763	cout.width(8);
407			cout.precision(8);
408
409
410			cout << info->getTime() << "\t"
411			<< chi << "\t"
412			<< eta << "\t"
413			<< U << "\t"
414			<< TS << "\t"
415			<< PV << "\t"
416			<< extra << "\t"
417			<< U+TS+PV+extra << endl;
418
419			conservedQuantity = U+TS+PV+extra;
420			return conservedQuantity;
421			}