OOPSE/libmdtools/NPTf.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" 


// Basic non-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> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  int i, j;
  chi = 0.0;
  integralOfChidt = 0.0;

  for(i = 0; i < 3; i++) 
    for (j = 0; j < 3; j++) 
      eta[i][j] = 0.0;

  have_tau_thermostat = 0;
  have_tau_barostat = 0;
  have_target_temp = 0;
  have_target_pressure = 0;
}

template<typename T> void NPTf<T>::moveA() {
  
  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;
  double sc[3];
  double eta2ij;
  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
  double bigScale, smallScale, offDiagMax;

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

  instaTemp = tStats->getTemperature();
  tStats->getPressureTensor(press);
  instaVol = tStats->getVolume();
   
  // first evolve chi a half step
  
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;

  for (i = 0; i < 3; i++ ) {
    for (j = 0; j < 3; j++ ) {
      if (i == j) {
        
        eta[i][j] += dt2 * instaVol * 
          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
        
        vScale[i][j] = eta[i][j] + chi;
          
      } else {
        
        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);

        vScale[i][j] = eta[i][j];
        
      }
    }
  }

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

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

    mass = atoms[i]->getMass();
    
    // velocity half step
        
    info->matVecMul3( vScale, vel, sc );
    
    for (j = 0; j < 3; j++) {
      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
      rj[j] = pos[j];
    }

    atoms[i]->setVel( vel );

    // position whole step    

    info->wrapVector(rj);

    info->matVecMul3( eta, rj, sc );

    for (j = 0; j < 3; j++ ) 
      pos[j] += dt * (vel[j] + sc[j]);

    atoms[i]->setPos( pos );
   
    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  );    
    }                    
  }
  
  // Scale the box after all the positions have been moved:
  
  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
  
  bigScale = 1.0;
  smallScale = 1.0;
  offDiagMax = 0.0;
  
  for(i=0; i<3; i++){
    for(j=0; j<3; j++){
      
      // Calculate the matrix Product of the eta array (we only need
      // the ij element right now):
      
      eta2ij = 0.0;
      for(k=0; k<3; k++){
        eta2ij += eta[i][k] * eta[k][j];
      }
      
      scaleMat[i][j] = 0.0;
      // identity matrix (see above):
      if (i == j) scaleMat[i][j] = 1.0;
      // Taylor expansion for the exponential truncated at second order:
      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;

      if (i != j)
        if (fabs(scaleMat[i][j]) > offDiagMax) 
          offDiagMax = fabs(scaleMat[i][j]);
    }

    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
  }
  
  if ((bigScale > 1.1) || (smallScale < 0.9)) {
    sprintf( painCave.errMsg,
             "NPTf error: Attempting a Box scaling of more than 10 percent.\n"
             " Check your tauBarostat, as it is probably too small!\n\n"
             " scaleMat = [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n",
             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
    painCave.isFatal = 1;
    simError();
  } else if (offDiagMax > 0.1) {
    sprintf( painCave.errMsg,
             "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
             " Check your tauBarostat, as it is probably too small!\n\n"
             " scaleMat = [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n",
             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
    painCave.isFatal = 1;
    simError();
  } else {
    info->getBoxM(hm);
    info->matMul3(hm, scaleMat, hmnew);
    info->setBoxM(hmnew);
  }
  
}

template<typename T> void NPTf<T>::moveB( void ){

  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double vel[3], frc[3];
  double mass;

  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double sc[3];
  double press[3][3], vScale[3][3];
  
  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  instaTemp = tStats->getTemperature();
  tStats->getPressureTensor(press);
  instaVol = tStats->getVolume();
   
  // first evolve chi a half step
  
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  
  for (i = 0; i < 3; i++ ) {
    for (j = 0; j < 3; j++ ) {
      if (i == j) {

        eta[i][j] += dt2 * instaVol * 
          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);

        vScale[i][j] = eta[i][j] + chi;
        
      } else {
        
        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);

        vScale[i][j] = eta[i][j];
        
      }
    }
  }

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

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

    mass = atoms[i]->getMass();
    
    // velocity half step
        
    info->matVecMul3( vScale, vel, sc );
    
    for (j = 0; j < 3; j++) {
      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
    }

    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);
      
      dAtom->setJ( ji );

    }                    
  }
}

template<typename T> void NPTf<T>::resetIntegrator() {
  int i,j;
  
  chi = 0.0;

  for(i = 0; i < 3; i++) 
    for (j = 0; j < 3; j++) 
      eta[i][j] = 0.0;

}

template<typename T> int NPTf<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,
             "NPTf error: You can't use the NPTf integrator\n"
             "   without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

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

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

  // We need NkBT a lot, so just set it here:

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

  return 1;
}

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

  double conservedQuantity;
  double tb2;
  double eta2[3][3];  
  double trEta;

  //HNVE
  conservedQuantity = tStats->getTotalE();

  //HNVT
  conservedQuantity += (info->getNDF() * kB * targetTemp * 
    (integralOfChidt + tauThermostat * tauThermostat * chi * chi /2)) / eConvert;

  //HNPT
  tb2 = tauBarostat *tauBarostat;

  trEta = info->matTrace3(eta);
  
  conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
                        3*NkBT/2 * tb2 * trEta * trEta) / eConvert;
  
  return conservedQuantity; 
}
Revision:	763
Committed:	Mon Sep 15 16:52:02 2003 UTC (20 years, 9 months ago) by tim
File size:	9931 byte(s)
Log Message:	add conserved quantity to statWriter fix bug of vector wrapping at NPTi
#	User	Rev	Content
1	gezelter	617	#include <cmath>
2	gezelter	576	#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
13	gezelter	578	// Basic non-isotropic thermostating and barostating via the Melchionna
14	gezelter	576	// modification of the Hoover algorithm:
15			//
16			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
17			// Molec. Phys., 78, 533.
18			//
19			// and
20			//
21			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
22
23	tim	645	template<typename T> NPTf<T>::NPTf ( SimInfo theInfo, ForceFields the_ff):
24			T( theInfo, the_ff )
25	gezelter	576	{
26	gezelter	588	int i, j;
27	gezelter	576	chi = 0.0;
28	tim	763	integralOfChidt = 0.0;
29	gezelter	588
30			for(i = 0; i < 3; i++)
31	mmeineke	590	for (j = 0; j < 3; j++)
32	gezelter	588	eta[i][j] = 0.0;
33
34	gezelter	576	have_tau_thermostat = 0;
35			have_tau_barostat = 0;
36			have_target_temp = 0;
37			have_target_pressure = 0;
38			}
39
40	tim	645	template<typename T> void NPTf<T>::moveA() {
41	gezelter	576
42	gezelter	600	int i, j, k;
43	gezelter	576	DirectionalAtom* dAtom;
44	gezelter	600	double Tb[3], ji[3];
45			double A[3][3], I[3][3];
46			double angle, mass;
47			double vel[3], pos[3], frc[3];
48
49			double rj[3];
50			double instaTemp, instaPress, instaVol;
51			double tt2, tb2;
52			double sc[3];
53			double eta2ij;
54	gezelter	588	double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
55	gezelter	617	double bigScale, smallScale, offDiagMax;
56	gezelter	576
57			tt2 = tauThermostat * tauThermostat;
58			tb2 = tauBarostat * tauBarostat;
59
60			instaTemp = tStats->getTemperature();
61	gezelter	577	tStats->getPressureTensor(press);
62	gezelter	576	instaVol = tStats->getVolume();
63
64			// first evolve chi a half step
65
66			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
67	gezelter	588
68			for (i = 0; i < 3; i++ ) {
69			for (j = 0; j < 3; j++ ) {
70			if (i == j) {
71	gezelter	600
72	gezelter	588	eta[i][j] += dt2 * instaVol *
73			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
74	gezelter	600
75	gezelter	588	vScale[i][j] = eta[i][j] + chi;
76	tim	763
77	gezelter	588	} else {
78
79			eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
80
81			vScale[i][j] = eta[i][j];
82
83			}
84			}
85			}
86
87	gezelter	576	for( i=0; i<nAtoms; i++ ){
88	gezelter	600
89			atoms[i]->getVel( vel );
90			atoms[i]->getPos( pos );
91			atoms[i]->getFrc( frc );
92
93			mass = atoms[i]->getMass();
94	gezelter	576
95			// velocity half step
96	gezelter	600
97			info->matVecMul3( vScale, vel, sc );
98	gezelter	577
99	gezelter	600	for (j = 0; j < 3; j++) {
100			vel[j] += dt2 * ((frc[j] / mass) * eConvert - sc[j]);
101			rj[j] = pos[j];
102			}
103	gezelter	576
104	gezelter	600	atoms[i]->setVel( vel );
105	gezelter	577
106	gezelter	576	// position whole step
107
108	gezelter	600	info->wrapVector(rj);
109	gezelter	576
110	gezelter	600	info->matVecMul3( eta, rj, sc );
111	gezelter	576
112	gezelter	600	for (j = 0; j < 3; j++ )
113			pos[j] += dt * (vel[j] + sc[j]);
114	gezelter	617
115			atoms[i]->setPos( pos );
116	gezelter	576
117			if( atoms[i]->isDirectional() ){
118
119			dAtom = (DirectionalAtom *)atoms[i];
120
121			// get and convert the torque to body frame
122
123	gezelter	600	dAtom->getTrq( Tb );
124	gezelter	576	dAtom->lab2Body( Tb );
125
126			// get the angular momentum, and propagate a half step
127
128	gezelter	600	dAtom->getJ( ji );
129
130			for (j=0; j < 3; j++)
131			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
132	gezelter	576
133			// use the angular velocities to propagate the rotation matrix a
134			// full time step
135	gezelter	600
136			dAtom->getA(A);
137			dAtom->getI(I);
138
139	gezelter	576	// rotate about the x-axis
140	gezelter	600	angle = dt2 * ji[0] / I[0][0];
141			this->rotate( 1, 2, angle, ji, A );
142
143	gezelter	576	// rotate about the y-axis
144	gezelter	600	angle = dt2 * ji[1] / I[1][1];
145			this->rotate( 2, 0, angle, ji, A );
146	gezelter	576
147			// rotate about the z-axis
148	gezelter	600	angle = dt * ji[2] / I[2][2];
149			this->rotate( 0, 1, angle, ji, A);
150	gezelter	576
151			// rotate about the y-axis
152	gezelter	600	angle = dt2 * ji[1] / I[1][1];
153			this->rotate( 2, 0, angle, ji, A );
154	gezelter	576
155			// rotate about the x-axis
156	gezelter	600	angle = dt2 * ji[0] / I[0][0];
157			this->rotate( 1, 2, angle, ji, A );
158	gezelter	576
159	gezelter	600	dAtom->setJ( ji );
160			dAtom->setA( A );
161			}
162	gezelter	576	}
163	gezelter	600
164	gezelter	577	// Scale the box after all the positions have been moved:
165	gezelter	600
166	gezelter	578	// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
167			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
168	gezelter	600
169	gezelter	617	bigScale = 1.0;
170			smallScale = 1.0;
171			offDiagMax = 0.0;
172	gezelter	600
173	gezelter	578	for(i=0; i<3; i++){
174			for(j=0; j<3; j++){
175	gezelter	600
176	gezelter	588	// Calculate the matrix Product of the eta array (we only need
177			// the ij element right now):
178	gezelter	600
179	gezelter	588	eta2ij = 0.0;
180	gezelter	578	for(k=0; k<3; k++){
181	gezelter	588	eta2ij += eta[i][k] * eta[k][j];
182	gezelter	578	}
183	gezelter	588
184			scaleMat[i][j] = 0.0;
185			// identity matrix (see above):
186			if (i == j) scaleMat[i][j] = 1.0;
187			// Taylor expansion for the exponential truncated at second order:
188			scaleMat[i][j] += dteta[i][j] + 0.5dtdteta2ij;
189	gezelter	617
190			if (i != j)
191			if (fabs(scaleMat[i][j]) > offDiagMax)
192			offDiagMax = fabs(scaleMat[i][j]);
193	gezelter	578	}
194	gezelter	617
195			if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
196			if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
197	gezelter	578	}
198	gezelter	600
199	gezelter	617	if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
200			sprintf( painCave.errMsg,
201			"NPTf error: Attempting a Box scaling of more than 10 percent.\n"
202			" Check your tauBarostat, as it is probably too small!\n\n"
203			" scaleMat = [%lf\t%lf\t%lf]\n"
204			" [%lf\t%lf\t%lf]\n"
205			" [%lf\t%lf\t%lf]\n",
206			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
207			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
208			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
209			painCave.isFatal = 1;
210			simError();
211			} else if (offDiagMax > 0.1) {
212			sprintf( painCave.errMsg,
213			"NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
214			" Check your tauBarostat, as it is probably too small!\n\n"
215			" scaleMat = [%lf\t%lf\t%lf]\n"
216			" [%lf\t%lf\t%lf]\n"
217			" [%lf\t%lf\t%lf]\n",
218			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
219			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
220			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
221			painCave.isFatal = 1;
222			simError();
223			} else {
224			info->getBoxM(hm);
225			info->matMul3(hm, scaleMat, hmnew);
226			info->setBoxM(hmnew);
227			}
228	gezelter	577
229	gezelter	576	}
230
231	tim	645	template<typename T> void NPTf<T>::moveB( void ){
232	gezelter	600
233			int i, j;
234	gezelter	576	DirectionalAtom* dAtom;
235	gezelter	600	double Tb[3], ji[3];
236			double vel[3], frc[3];
237			double mass;
238
239			double instaTemp, instaPress, instaVol;
240	gezelter	576	double tt2, tb2;
241	gezelter	600	double sc[3];
242	gezelter	588	double press[3][3], vScale[3][3];
243	gezelter	576
244			tt2 = tauThermostat * tauThermostat;
245			tb2 = tauBarostat * tauBarostat;
246
247			instaTemp = tStats->getTemperature();
248	gezelter	578	tStats->getPressureTensor(press);
249	gezelter	576	instaVol = tStats->getVolume();
250	gezelter	578
251			// first evolve chi a half step
252
253	gezelter	576	chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
254
255	gezelter	588	for (i = 0; i < 3; i++ ) {
256			for (j = 0; j < 3; j++ ) {
257			if (i == j) {
258	gezelter	578
259	gezelter	588	eta[i][j] += dt2 * instaVol *
260			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
261
262			vScale[i][j] = eta[i][j] + chi;
263
264			} else {
265
266			eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
267
268			vScale[i][j] = eta[i][j];
269
270			}
271			}
272			}
273
274	gezelter	576	for( i=0; i<nAtoms; i++ ){
275	gezelter	578
276	gezelter	600	atoms[i]->getVel( vel );
277			atoms[i]->getFrc( frc );
278
279			mass = atoms[i]->getMass();
280
281	gezelter	576	// velocity half step
282	gezelter	600
283			info->matVecMul3( vScale, vel, sc );
284	gezelter	576
285	gezelter	600	for (j = 0; j < 3; j++) {
286			vel[j] += dt2 * ((frc[j] / mass) * eConvert - sc[j]);
287			}
288	gezelter	578
289	gezelter	600	atoms[i]->setVel( vel );
290	gezelter	578
291	gezelter	576	if( atoms[i]->isDirectional() ){
292	gezelter	600
293	gezelter	576	dAtom = (DirectionalAtom *)atoms[i];
294	gezelter	600
295	gezelter	576	// get and convert the torque to body frame
296
297	gezelter	600	dAtom->getTrq( Tb );
298	gezelter	576	dAtom->lab2Body( Tb );
299
300	gezelter	600	// get the angular momentum, and propagate a half step
301	gezelter	576
302	gezelter	600	dAtom->getJ( ji );
303	gezelter	576
304	gezelter	600	for (j=0; j < 3; j++)
305			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
306	gezelter	576
307	gezelter	600	dAtom->setJ( ji );
308
309			}
310	gezelter	576	}
311			}
312
313	mmeineke	746	template<typename T> void NPTf<T>::resetIntegrator() {
314			int i,j;
315
316			chi = 0.0;
317
318			for(i = 0; i < 3; i++)
319			for (j = 0; j < 3; j++)
320			eta[i][j] = 0.0;
321
322			}
323
324	tim	645	template<typename T> int NPTf<T>::readyCheck() {
325	tim	658
326			//check parent's readyCheck() first
327			if (T::readyCheck() == -1)
328			return -1;
329	gezelter	576
330			// First check to see if we have a target temperature.
331			// Not having one is fatal.
332
333			if (!have_target_temp) {
334			sprintf( painCave.errMsg,
335	gezelter	580	"NPTf error: You can't use the NPTf integrator\n"
336	gezelter	576	" without a targetTemp!\n"
337			);
338			painCave.isFatal = 1;
339			simError();
340			return -1;
341			}
342
343			if (!have_target_pressure) {
344			sprintf( painCave.errMsg,
345	gezelter	580	"NPTf error: You can't use the NPTf integrator\n"
346	gezelter	576	" without a targetPressure!\n"
347			);
348			painCave.isFatal = 1;
349			simError();
350			return -1;
351			}
352
353			// We must set tauThermostat.
354
355			if (!have_tau_thermostat) {
356			sprintf( painCave.errMsg,
357	gezelter	580	"NPTf error: If you use the NPTf\n"
358	gezelter	576	" integrator, you must set tauThermostat.\n");
359			painCave.isFatal = 1;
360			simError();
361			return -1;
362			}
363
364			// We must set tauBarostat.
365
366			if (!have_tau_barostat) {
367			sprintf( painCave.errMsg,
368	gezelter	580	"NPTf error: If you use the NPTf\n"
369	gezelter	576	" integrator, you must set tauBarostat.\n");
370			painCave.isFatal = 1;
371			simError();
372			return -1;
373			}
374
375			// We need NkBT a lot, so just set it here:
376
377			NkBT = (double)info->ndf * kB * targetTemp;
378
379			return 1;
380			}
381	tim	763
382			template<typename T> double NPTf<T>::getConservedQuantity(void){
383
384			double conservedQuantity;
385			double tb2;
386			double eta2[3][3];
387			double trEta;
388
389			//HNVE
390			conservedQuantity = tStats->getTotalE();
391
392			//HNVT
393			conservedQuantity += (info->getNDF() * kB * targetTemp *
394			(integralOfChidt + tauThermostat * tauThermostat * chi * chi /2)) / eConvert;
395
396			//HNPT
397			tb2 = tauBarostat *tauBarostat;
398
399			trEta = info->matTrace3(eta);
400
401			conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
402			3NkBT/2 tb2 * trEta * trEta) / eConvert;
403
404			return conservedQuantity;
405			}