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: 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:	770
Committed:	Fri Sep 19 14:55:41 2003 UTC (20 years, 9 months ago) by gezelter
File size:	10250 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			// Set things up for the iteration:
220
221			oldChi = chi;
222			oldEta = eta;
223
224			for( i=0; i<nAtoms; i++ ){
225
226			atoms[i]->getVel( vel );
227
228			for (j=0; j < 3; j++)
229			oldVel[3*i + j] = vel[j];
230
231			if( atoms[i]->isDirectional() ){
232
233			dAtom = (DirectionalAtom *)atoms[i];
234
235			dAtom->getJ( ji );
236
237			for (j=0; j < 3; j++)
238			oldJi[3*i + j] = ji[j];
239
240			}
241			}
242
243			// do the iteration:
244
245			instVol = tStats->getVolume();
246
247			for (k=0; k < 4; k++) {
248
249			instTemp = tStats->getTemperature();
250			instPress = tStats->getPressure();
251
252			// evolve chi another half step using the temperature at t + dt/2
253
254			prevChi = chi;
255			chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
256			(tauThermostat*tauThermostat);
257
258			preEta = eta;
259			eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
260			(p_convertNkBTtb2));
261
262
263			for( i=0; i<nAtoms; i++ ){
264
265			atoms[i]->getFrc( frc );
266			atoms[i]->getVel(vel);
267
268			mass = atoms[i]->getMass();
269
270			// velocity half step
271			for (j=0; j < 3; j++)
272			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j](chi + eta));
273
274			atoms[i]->setVel( vel );
275
276			if( atoms[i]->isDirectional() ){
277
278			dAtom = (DirectionalAtom *)atoms[i];
279
280			// get and convert the torque to body frame
281
282			dAtom->getTrq( Tb );
283			dAtom->lab2Body( Tb );
284
285			for (j=0; j < 3; j++)
286			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
287
288			dAtom->setJ( ji );
289			}
290			}
291	mmeineke	768
292			if (nConstrained){
293			constrainB();
294			}
295
296			if (fabs(prevChi - chi) <=
297			chiTolerance && fabs(preEta -eta) <= etaTolerance)
298	tim	763	break;
299			}
300
301			//calculate integral of chida
302			integralOfChidt += dt2*chi;
303
304
305	gezelter	574	}
306
307	mmeineke	746	template<typename T> void NPTi<T>::resetIntegrator() {
308			chi = 0.0;
309			eta = 0.0;
310			}
311
312	tim	645	template<typename T> int NPTi<T>::readyCheck() {
313	tim	658
314			//check parent's readyCheck() first
315			if (T::readyCheck() == -1)
316			return -1;
317	gezelter	574
318			// First check to see if we have a target temperature.
319			// Not having one is fatal.
320
321			if (!have_target_temp) {
322			sprintf( painCave.errMsg,
323			"NPTi error: You can't use the NPTi integrator\n"
324			" without a targetTemp!\n"
325			);
326			painCave.isFatal = 1;
327			simError();
328			return -1;
329			}
330
331			if (!have_target_pressure) {
332			sprintf( painCave.errMsg,
333			"NPTi error: You can't use the NPTi integrator\n"
334			" without a targetPressure!\n"
335			);
336			painCave.isFatal = 1;
337			simError();
338			return -1;
339			}
340
341			// We must set tauThermostat.
342
343			if (!have_tau_thermostat) {
344			sprintf( painCave.errMsg,
345			"NPTi error: If you use the NPTi\n"
346			" integrator, you must set tauThermostat.\n");
347			painCave.isFatal = 1;
348			simError();
349			return -1;
350			}
351
352			// We must set tauBarostat.
353
354			if (!have_tau_barostat) {
355			sprintf( painCave.errMsg,
356			"NPTi error: If you use the NPTi\n"
357			" integrator, you must set tauBarostat.\n");
358			painCave.isFatal = 1;
359			simError();
360			return -1;
361			}
362
363	tim	763	if (!have_chi_tolerance) {
364			sprintf( painCave.errMsg,
365			"NPTi warning: setting chi tolerance to 1e-6\n");
366			chiTolerance = 1e-6;
367			have_chi_tolerance = 1;
368			painCave.isFatal = 0;
369			simError();
370			}
371
372	gezelter	770	if (!have_eta_tolerance) {
373	tim	763	sprintf( painCave.errMsg,
374			"NPTi warning: setting eta tolerance to 1e-6\n");
375			etaTolerance = 1e-6;
376			have_eta_tolerance = 1;
377			painCave.isFatal = 0;
378			simError();
379			}
380	gezelter	770
381
382			// We need NkBT a lot, so just set it here: This is the RAW number
383			// of particles, so no subtraction or addition of constraints or
384			// orientational degrees of freedom:
385
386	tim	763	NkBT = (double)Nparticles * kB * targetTemp;
387	gezelter	770
388			// fkBT is used because the thermostat operates on more degrees of freedom
389			// than the barostat (when there are particles with orientational degrees
390			// of freedom). ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
391
392	tim	763	fkBT = (double)info->ndf * kB * targetTemp;
393	gezelter	574
394			return 1;
395			}
396	tim	763
397			template<typename T> double NPTi<T>::getConservedQuantity(void){
398
399			double conservedQuantity;
400	gezelter	770	double Three_NkBT;
401	tim	769	double Energy;
402			double thermostat_kinetic;
403			double thermostat_potential;
404			double barostat_kinetic;
405			double barostat_potential;
406	tim	763	double tb2;
407	gezelter	770	double eta2;
408	tim	763
409	tim	769	Energy = tStats->getTotalE();
410	tim	763
411	tim	769	thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
412			(2.0 * eConvert);
413	tim	763
414	tim	769	thermostat_potential = fkBT* integralOfChidt / eConvert;
415	tim	763
416
417	gezelter	770	barostat_kinetic = 3.0 * NkBT * tauBarostat * tauBarostat * eta * eta /
418	tim	769	(2.0 * eConvert);
419
420			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
421			eConvert;
422	tim	767
423	tim	769	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
424			barostat_kinetic + barostat_potential;
425
426	tim	763	cout.width(8);
427			cout.precision(8);
428
429	tim	769	cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
430			"\t" << thermostat_potential << "\t" << barostat_kinetic <<
431			"\t" << barostat_potential << "\t" << conservedQuantity << endl;
432	tim	763
433			return conservedQuantity;
434			}