OOPSE/libmdtools/NPT.cpp

#include <stdlib.h>
#include <math.h>

#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.

NPT::NPT ( SimInfo *theInfo, ForceFields* the_ff):
  Integrator( theInfo, the_ff )
{
  GenericData* data;
  
  chi = 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;

  // retrieve chi and integralOfChidt from simInfo
  data = info->getProperty(CHIVALUE_ID);
  if(data != NULL ){
    chi = data->getDval();
  }

  data = info->getProperty(INTEGRALOFCHIDT_ID);
  if(data != NULL ){
    integralOfChidt = data->getDval();
  }

  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

}

NPT::~NPT() {
  delete[] oldPos;
  delete[] oldVel;
  delete[] oldJi;
}

void NPT::moveA() {

  //new version of NPT
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double mass;
  double vel[3], pos[3], frc[3];
  double sc[3];
  double COM[3];

  instaTemp = tStats->getTemperature();
  tStats->getPressureTensor( press );
  instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
  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();

    getVelScaleA( sc, vel );

    for (j=0; j < 3; j++) {

      // velocity half step  (use chi from previous step here):
      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);

      this->rotationPropagation( dAtom, ji );

      dAtom->setJ( ji );
    }
  }

  // evolve chi and eta  half step

  evolveChiA();
  evolveEtaA();

  //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 < 5; k ++){

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

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

      this->getPosScale( pos, COM, i, sc );

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

      atoms[i]->setPos( pos );
    }

    if (nConstrained){
      constrainA();
    }
  }


  // Scale the box after all the positions have been moved:

  this->scaleSimBox();
}

void NPT::moveB( void ){

  //new version of NPT
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3], sc[3];
  double vel[3], frc[3], qVel[3];
  double mass;

  // Set things up for the iteration:

  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++) {

    atoms[0]->getVel(vel);

 //    std::cerr << "vel atom0 = "
//            << vel[0] << ", " << vel[1] << ", " << vel[2] << "\n";
    
    instaTemp = tStats->getTemperature();
    
    std::cerr << "instaTemp MoveB = " << instaTemp << "\n";
    
    instaPress = tStats->getPressure();

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

    this->evolveChiB();
    this->evolveEtaB();

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

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

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

      getVelScaleB( sc, i );

      for(j=0;j<3;j++)
        qVel[j] = vel[j];

      // velocity half step
      for (j=0; j < 3; j++)
        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);

//       std::cerr << " vel of " << i << " => " 
//              << vel[0] << ", " << vel[1] << ", " << vel[2] << "\n";

      double delVel = 0.0;
      for(j=0;j<3;j++)
        delVel += pow( (qVel[j] - vel[j]), 2 );
      
      delVel = sqrt(delVel);
      
      std::cerr << "delVel[" << i << "] = " << delVel << "\n";

      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 ( this->chiConverged() && this->etaConverged() ) break;
  }

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


}

void NPT::resetIntegrator() {
  chi = 0.0;
  Integrator::resetIntegrator();
}

void NPT::evolveChiA() {

  std::cerr << "\n"
            << "evolveChiA before:\n"
            << "chi = " << chi << "\n"
            << "oldChi = " << oldChi << "\n"
            << "\n"
            << "instaTemp = " << instaTemp << "\n"
            << "targetTemp = " << targetTemp << "\n"
            << "dt2 = " << dt2 << "\n"
            << "tt2 = " << tt2 << "\n";

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

  std::cerr << "\n"
            << "evolveChiA after:\n"
            << "chi = " << chi << "\n"
            << "oldChi = " << oldChi << "\n"
            << "\n"
            << "instaTemp = " << instaTemp << "\n"
            << "targetTemp = " << targetTemp << "\n"
            << "dt2 = " << dt2 << "\n"
            << "tt2 = " << tt2 << "\n";

}

void NPT::evolveChiB() {

  std::cerr << "\n"
            << "evolveChiB before:\n"
            << "chi = " << chi << "\n"
            << "oldChi = " << oldChi << "\n"
            << "\n"
            << "instaTemp = " << instaTemp << "\n"
            << "targetTemp = " << targetTemp << "\n"
            << "dt2 = " << dt2 << "\n"
            << "tt2 = " << tt2 << "\n";

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

  std::cerr << "\n"
            << "evolveChiB after:\n"
            << "chi = " << chi << "\n"
            << "oldChi = " << oldChi << "\n"
            << "\n"
            << "instaTemp = " << instaTemp << "\n"
            << "targetTemp = " << targetTemp << "\n"
            << "dt2 = " << dt2 << "\n"
            << "tt2 = " << tt2 << "\n";
}

bool NPT::chiConverged() {

  return ( fabs( prevChi - chi ) <= chiTolerance );
}

int NPT::readyCheck() {

  //check parent's readyCheck() first
  if (Integrator::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,
             "NPT error: You can't use the NPT integrator\n"
             "   without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

  if (!have_target_pressure) {
    sprintf( painCave.errMsg,
             "NPT error: You can't use the NPT integrator\n"
             "   without a targetPressure!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

  // We must set tauThermostat.

  if (!have_tau_thermostat) {
    sprintf( painCave.errMsg,
             "NPT error: If you use the NPT\n"
             "   integrator, you must set tauThermostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }

  // We must set tauBarostat.

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

  if (!have_chi_tolerance) {
    sprintf( painCave.errMsg,
             "NPT 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,
             "NPT 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;

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

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