OOPSE/libmdtools/NPTf.cpp

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

NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
  Integrator( theInfo, the_ff )
{
  int i;
  chi = 0.0;
  for(i = 0; i < 9; i++) eta[i] = 0.0;
  have_tau_thermostat = 0;
  have_tau_barostat = 0;
  have_target_temp = 0;
  have_target_pressure = 0;
}

void NPTf::moveA() {
  
  int i,j,k;
  int atomIndex, aMatIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];
  double rj[3];
  double ident[3][3], eta1[3][3], eta2[3][3], hmnew[3][3];
  double hm[9];
  double vx, vy, vz;
  double scx, scy, scz;
  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double angle;
  double press[9];

  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;
  
  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) / 
    (NkBT*tb2);
  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
  eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) / 
    (NkBT*tb2);
  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
  eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) / 
    (NkBT*tb2);
  
  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;
    aMatIndex = i * 9;
    
    // velocity half step
    
    vx = vel[atomIndex];
    vy = vel[atomIndex+1];
    vz = vel[atomIndex+2];
    
    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
    
    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);

    vel[atomIndex] = vx;
    vel[atomIndex+1] = vy;
    vel[atomIndex+2] = vz;

    // position whole step    

    rj[0] = pos[atomIndex];
    rj[1] = pos[atomIndex+1];
    rj[2] = pos[atomIndex+2];

    info->wrapVector(rj);

    scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2];
    scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2];
    scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2];

    pos[atomIndex] += dt * (vel[atomIndex] + scx);
    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
   
    if( atoms[i]->isDirectional() ){

      dAtom = (DirectionalAtom *)atoms[i];
          
      // get and convert the torque to body frame
      
      Tb[0] = dAtom->getTx();
      Tb[1] = dAtom->getTy();
      Tb[2] = dAtom->getTz();
      
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and propagate a half step

      ji[0] = dAtom->getJx();
      ji[1] = dAtom->getJy();
      ji[2] = dAtom->getJz();
      
      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
      
      // use the angular velocities to propagate the rotation matrix a
      // full time step
      
      // rotate about the x-axis      
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); 
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
      
      // rotate about the z-axis
      angle = dt * ji[2] / dAtom->getIzz();
      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
    
  }

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


  for(i=0; i<3; i++){
    for(j=0; j<3; j++){
      ident[i][j] = 0.0;
      eta1[i][j] = eta[3*i+j];
      eta2[i][j] = 0.0;
      for(k=0; k<3; k++){
        eta2[i][j] += eta[3*i+k] * eta[3*k+j];
      }
    }
    ident[i][i] = 1.0;
  }

   
  info->getBoxM(hm);
 
  for(i=0; i<3; i++){
    for(j=0; j<3; j++){      
      hmnew[i][j] = 0.0;
      for(k=0; k<3; k++){
        // remember that hmat has transpose ordering for Fortran compat:
        hmnew[i][j] += hm[3*k+i] * (ident[k][j] 
                                    + dt * eta1[k][j] 
                                    + 0.5 * dt * dt * eta2[k][j]);
      }
    }
  }
  
  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      // remember that hmat has transpose ordering for Fortran compat:
      hm[3*j + i] = hmnew[i][j];
    }
  }

  info->setBoxM(hm);
  
}

void NPTf::moveB( void ){
  int i,j,k;
  int atomIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];
  double press[9];
  double instaTemp, instaVol;
  double tt2, tb2;
  double vx, vy, vz;
  double scx, scy, scz;
  const double p_convert = 1.63882576e8;
  
  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;
  
  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) / 
    (NkBT*tb2);
  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
  eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) / 
    (NkBT*tb2);
  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
  eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) / 
    (NkBT*tb2);

  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;

    // velocity half step
    
    vx = vel[atomIndex];
    vy = vel[atomIndex+1];
    vz = vel[atomIndex+2];
    
    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
    
    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);

    vel[atomIndex] = vx;
    vel[atomIndex+1] = vy;
    vel[atomIndex+2] = vz;
    
    if( atoms[i]->isDirectional() ){
      
      dAtom = (DirectionalAtom *)atoms[i];
      
      // get and convert the torque to body frame
      
      Tb[0] = dAtom->getTx();
      Tb[1] = dAtom->getTy();
      Tb[2] = dAtom->getTz();
      
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and complete the angular momentum
      // half step
      
      ji[0] = dAtom->getJx();
      ji[1] = dAtom->getJy();
      ji[2] = dAtom->getJz();
      
      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
  }
}

int NPTf::readyCheck() {
 
  // 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;
}
Revision:	586
Committed:	Wed Jul 9 22:14:06 2003 UTC (21 years ago) by mmeineke
File size:	9436 byte(s)
Log Message:	Bug fixing NPTi and NPTf. there is some error in the caclulation of HmatInverse.
#	User	Rev	Content
1	gezelter	576	#include "Atom.hpp"
2			#include "SRI.hpp"
3			#include "AbstractClasses.hpp"
4			#include "SimInfo.hpp"
5			#include "ForceFields.hpp"
6			#include "Thermo.hpp"
7			#include "ReadWrite.hpp"
8			#include "Integrator.hpp"
9			#include "simError.h"
10
11
12	gezelter	578	// Basic non-isotropic thermostating and barostating via the Melchionna
13	gezelter	576	// modification of the Hoover algorithm:
14			//
15			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
16			// Molec. Phys., 78, 533.
17			//
18			// and
19			//
20			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
21
22	gezelter	577	NPTf::NPTf ( SimInfo theInfo, ForceFields the_ff):
23	gezelter	576	Integrator( theInfo, the_ff )
24			{
25			int i;
26			chi = 0.0;
27			for(i = 0; i < 9; i++) eta[i] = 0.0;
28			have_tau_thermostat = 0;
29			have_tau_barostat = 0;
30			have_target_temp = 0;
31			have_target_pressure = 0;
32			}
33
34	gezelter	577	void NPTf::moveA() {
35	gezelter	576
36			int i,j,k;
37			int atomIndex, aMatIndex;
38			DirectionalAtom* dAtom;
39			double Tb[3];
40			double ji[3];
41			double rj[3];
42	gezelter	578	double ident[3][3], eta1[3][3], eta2[3][3], hmnew[3][3];
43			double hm[9];
44			double vx, vy, vz;
45			double scx, scy, scz;
46	gezelter	576	double instaTemp, instaPress, instaVol;
47			double tt2, tb2;
48			double angle;
49	gezelter	577	double press[9];
50	gezelter	576
51			tt2 = tauThermostat * tauThermostat;
52			tb2 = tauBarostat * tauBarostat;
53
54			instaTemp = tStats->getTemperature();
55	gezelter	577	tStats->getPressureTensor(press);
56	gezelter	576	instaVol = tStats->getVolume();
57
58			// first evolve chi a half step
59
60			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
61
62	mmeineke	586	eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) /
63			(NkBT*tb2);
64	gezelter	577	eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
65			eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
66			eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
67	mmeineke	586	eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) /
68			(NkBT*tb2);
69	gezelter	577	eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
70			eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
71			eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
72	mmeineke	586	eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) /
73			(NkBT*tb2);
74	gezelter	577
75	gezelter	576	for( i=0; i<nAtoms; i++ ){
76			atomIndex = i * 3;
77			aMatIndex = i * 9;
78
79			// velocity half step
80	gezelter	577
81			vx = vel[atomIndex];
82			vy = vel[atomIndex+1];
83			vz = vel[atomIndex+2];
84
85			scx = (chi + eta[0])vx + eta[1]vy + eta[2]*vz;
86			scy = eta[3]vx + (chi + eta[4])vy + eta[5]*vz;
87			scz = eta[6]vx + eta[7]vy + (chi + eta[8])*vz;
88
89			vx += dt2 * ((frc[atomIndex] /atoms[i]->getMass())*eConvert - scx);
90			vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
91			vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
92	gezelter	576
93	gezelter	577	vel[atomIndex] = vx;
94			vel[atomIndex+1] = vy;
95			vel[atomIndex+2] = vz;
96
97	gezelter	576	// position whole step
98
99	gezelter	577	rj[0] = pos[atomIndex];
100			rj[1] = pos[atomIndex+1];
101			rj[2] = pos[atomIndex+2];
102	gezelter	576
103	gezelter	577	info->wrapVector(rj);
104	gezelter	576
105	gezelter	577	scx = eta[0]rj[0] + eta[1]rj[1] + eta[2]*rj[2];
106			scy = eta[3]rj[0] + eta[4]rj[1] + eta[5]*rj[2];
107			scz = eta[6]rj[0] + eta[7]rj[1] + eta[8]*rj[2];
108	gezelter	576
109	gezelter	577	pos[atomIndex] += dt * (vel[atomIndex] + scx);
110			pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
111			pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
112	gezelter	576
113			if( atoms[i]->isDirectional() ){
114
115			dAtom = (DirectionalAtom *)atoms[i];
116
117			// get and convert the torque to body frame
118
119			Tb[0] = dAtom->getTx();
120			Tb[1] = dAtom->getTy();
121			Tb[2] = dAtom->getTz();
122
123			dAtom->lab2Body( Tb );
124
125			// get the angular momentum, and propagate a half step
126
127			ji[0] = dAtom->getJx();
128			ji[1] = dAtom->getJy();
129			ji[2] = dAtom->getJz();
130
131			ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
132			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
133			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
134
135			// use the angular velocities to propagate the rotation matrix a
136			// full time step
137
138			// rotate about the x-axis
139			angle = dt2 * ji[0] / dAtom->getIxx();
140			this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
141
142			// rotate about the y-axis
143			angle = dt2 * ji[1] / dAtom->getIyy();
144			this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
145
146			// rotate about the z-axis
147			angle = dt * ji[2] / dAtom->getIzz();
148			this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
149
150			// rotate about the y-axis
151			angle = dt2 * ji[1] / dAtom->getIyy();
152			this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
153
154			// rotate about the x-axis
155			angle = dt2 * ji[0] / dAtom->getIxx();
156			this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
157
158			dAtom->setJx( ji[0] );
159			dAtom->setJy( ji[1] );
160			dAtom->setJz( ji[2] );
161			}
162
163			}
164	gezelter	577
165			// Scale the box after all the positions have been moved:
166
167	gezelter	578	// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
168			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
169	gezelter	577
170	gezelter	578
171			for(i=0; i<3; i++){
172			for(j=0; j<3; j++){
173			ident[i][j] = 0.0;
174			eta1[i][j] = eta[3*i+j];
175			eta2[i][j] = 0.0;
176			for(k=0; k<3; k++){
177			eta2[i][j] += eta[3i+k] eta[3*k+j];
178			}
179			}
180			ident[i][i] = 1.0;
181			}
182
183	gezelter	577
184			info->getBoxM(hm);
185	gezelter	578
186			for(i=0; i<3; i++){
187			for(j=0; j<3; j++){
188			hmnew[i][j] = 0.0;
189			for(k=0; k<3; k++){
190			// remember that hmat has transpose ordering for Fortran compat:
191			hmnew[i][j] += hm[3k+i] (ident[k][j]
192			+ dt * eta1[k][j]
193			+ 0.5 * dt * dt * eta2[k][j]);
194			}
195			}
196			}
197	gezelter	577
198	gezelter	578	for (i = 0; i < 3; i++) {
199			for (j = 0; j < 3; j++) {
200			// remember that hmat has transpose ordering for Fortran compat:
201	mmeineke	586	hm[3*j + i] = hmnew[i][j];
202	gezelter	578	}
203			}
204	gezelter	577
205	gezelter	578	info->setBoxM(hm);
206
207	gezelter	576	}
208
209	gezelter	578	void NPTf::moveB( void ){
210	gezelter	576	int i,j,k;
211			int atomIndex;
212			DirectionalAtom* dAtom;
213			double Tb[3];
214			double ji[3];
215	gezelter	578	double press[9];
216			double instaTemp, instaVol;
217	gezelter	576	double tt2, tb2;
218	gezelter	578	double vx, vy, vz;
219			double scx, scy, scz;
220			const double p_convert = 1.63882576e8;
221	gezelter	576
222			tt2 = tauThermostat * tauThermostat;
223			tb2 = tauBarostat * tauBarostat;
224
225			instaTemp = tStats->getTemperature();
226	gezelter	578	tStats->getPressureTensor(press);
227	gezelter	576	instaVol = tStats->getVolume();
228	gezelter	578
229			// first evolve chi a half step
230
231	gezelter	576	chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
232
233	mmeineke	586	eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) /
234			(NkBT*tb2);
235	gezelter	578	eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
236			eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
237			eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
238	mmeineke	586	eta[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) /
239			(NkBT*tb2);
240	gezelter	578	eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
241			eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
242			eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
243	mmeineke	586	eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) /
244			(NkBT*tb2);
245	gezelter	578
246	gezelter	576	for( i=0; i<nAtoms; i++ ){
247			atomIndex = i * 3;
248	gezelter	578
249	gezelter	576	// velocity half step
250
251	gezelter	578	vx = vel[atomIndex];
252			vy = vel[atomIndex+1];
253			vz = vel[atomIndex+2];
254
255			scx = (chi + eta[0])vx + eta[1]vy + eta[2]*vz;
256			scy = eta[3]vx + (chi + eta[4])vy + eta[5]*vz;
257			scz = eta[6]vx + eta[7]vy + (chi + eta[8])*vz;
258
259			vx += dt2 * ((frc[atomIndex] /atoms[i]->getMass())*eConvert - scx);
260			vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
261			vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
262
263			vel[atomIndex] = vx;
264			vel[atomIndex+1] = vy;
265			vel[atomIndex+2] = vz;
266
267	gezelter	576	if( atoms[i]->isDirectional() ){
268
269			dAtom = (DirectionalAtom *)atoms[i];
270
271			// get and convert the torque to body frame
272
273			Tb[0] = dAtom->getTx();
274			Tb[1] = dAtom->getTy();
275			Tb[2] = dAtom->getTz();
276
277			dAtom->lab2Body( Tb );
278
279			// get the angular momentum, and complete the angular momentum
280			// half step
281
282			ji[0] = dAtom->getJx();
283			ji[1] = dAtom->getJy();
284			ji[2] = dAtom->getJz();
285
286			ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
287			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
288			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
289
290			dAtom->setJx( ji[0] );
291			dAtom->setJy( ji[1] );
292			dAtom->setJz( ji[2] );
293			}
294			}
295			}
296
297	gezelter	580	int NPTf::readyCheck() {
298	gezelter	576
299			// First check to see if we have a target temperature.
300			// Not having one is fatal.
301
302			if (!have_target_temp) {
303			sprintf( painCave.errMsg,
304	gezelter	580	"NPTf error: You can't use the NPTf integrator\n"
305	gezelter	576	" without a targetTemp!\n"
306			);
307			painCave.isFatal = 1;
308			simError();
309			return -1;
310			}
311
312			if (!have_target_pressure) {
313			sprintf( painCave.errMsg,
314	gezelter	580	"NPTf error: You can't use the NPTf integrator\n"
315	gezelter	576	" without a targetPressure!\n"
316			);
317			painCave.isFatal = 1;
318			simError();
319			return -1;
320			}
321
322			// We must set tauThermostat.
323
324			if (!have_tau_thermostat) {
325			sprintf( painCave.errMsg,
326	gezelter	580	"NPTf error: If you use the NPTf\n"
327	gezelter	576	" integrator, you must set tauThermostat.\n");
328			painCave.isFatal = 1;
329			simError();
330			return -1;
331			}
332
333			// We must set tauBarostat.
334
335			if (!have_tau_barostat) {
336			sprintf( painCave.errMsg,
337	gezelter	580	"NPTf error: If you use the NPTf\n"
338	gezelter	576	" integrator, you must set tauBarostat.\n");
339			painCave.isFatal = 1;
340			simError();
341			return -1;
342			}
343
344			// We need NkBT a lot, so just set it here:
345
346			NkBT = (double)info->ndf * kB * targetTemp;
347
348			return 1;
349			}