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

 
//   int i, j;
//   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;

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

//   instaTemp = tStats->getTemperature();
//   instaPress = tStats->getPressure();
//   instaVol = tStats->getVolume();
   
//    // first evolve chi a half step
  
//   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
//   eta += dt2 * ( instaVol * (instaPress - targetPressure) / 
//               (p_convert*NkBT*tb2));

//   integralOfChidt += dt2* chi;

//   for( i=0; i<nAtoms; i++ ){
//     atoms[i]->getVel( vel );
//     atoms[i]->getPos( pos );
//     atoms[i]->getFrc( frc );

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

//     for (j=0; j < 3; j++) {
//       vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
//       rj[j] = pos[j];
//     }

//     atoms[i]->setVel( vel );

//     info->wrapVector(rj);

//     for (j = 0; j < 3; j++) 
//       pos[j] += dt * (vel[j] + eta*rj[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:
  
//   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(exp(dt*eta));      
//   }
  
 
  //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];

      
      //wrapVector(r(t)) = r(t)-R0
      //info->wrapVector(rj);
      
      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);      
  }  

  //advance volume;
  volume = volume * exp(dt*eta);
}

template<typename T> void NPTi<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;
  
  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  instaVol = tStats->getVolume();

  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  eta += dt2 * ( instaVol * (instaPress - targetPressure) / 
                 (p_convert*NkBT*tb2));
  integralOfChidt += dt2*chi;
  
  for( i=0; i<nAtoms; i++ ){

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

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

    // velocity half step
    for (j=0; j < 3; j++) 
      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);     

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

  static double pre_U;
  static double pre_TS;
  static double pre_PV;
  static double pre_extra;
  static int hackCount = 0;

  double delta_U;
  double delta_TS;
  double delta_PV;
  double delta_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 = (fkBT * tb2 * eta2 / 2.0 ) / eConvert;
  /*
  if(hackCount == 0){
    pre_U = U;
    pre_TS =TS;
    pre_PV = PV;
    pre_extra =extra;
    hackCount ++;
  }

  delta_U = U - pre_U;
  delta_TS = TS - pre_TS;
  delta_PV = PV - pre_PV;
  delta_extra = extra - pre_extra;
*/
  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;

/*
    pre_U = U;
    pre_TS =TS;
    pre_PV = PV;
    pre_extra =extra;


  cout << info->getTime() << "\t" 
       << U << "\t" 
       << U+TS << "\t"
       << U+TS+PV << "\t"
       << U+TS+PV+extra << endl;
*/
  conservedQuantity = U+TS+PV+extra;
  return conservedQuantity; 
}
Revision:	763
Committed:	Mon Sep 15 16:52:02 2003 UTC (20 years, 9 months ago) by tim
File size:	14537 byte(s)
Log Message:	add conserved quantity to statWriter fix bug of vector wrapping at NPTi
#	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
61			// int i, j;
62			// DirectionalAtom* dAtom;
63			// 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			// double rj[3];
69			// double instaTemp, instaPress, instaVol;
70			// double tt2, tb2, scaleFactor;
71
72			// tt2 = tauThermostat * tauThermostat;
73			// tb2 = tauBarostat * tauBarostat;
74
75			// instaTemp = tStats->getTemperature();
76			// instaPress = tStats->getPressure();
77			// instaVol = tStats->getVolume();
78
79			// // first evolve chi a half step
80	gezelter	574
81	tim	763	// chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
82			// eta += dt2 * ( instaVol * (instaPress - targetPressure) /
83			// (p_convertNkBTtb2));
84
85			// integralOfChidt += dt2* chi;
86
87			// for( i=0; i<nAtoms; i++ ){
88			// atoms[i]->getVel( vel );
89			// atoms[i]->getPos( pos );
90			// atoms[i]->getFrc( frc );
91
92			// mass = atoms[i]->getMass();
93
94			// for (j=0; j < 3; j++) {
95			// vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
96			// rj[j] = pos[j];
97			// }
98
99			// atoms[i]->setVel( vel );
100
101			// info->wrapVector(rj);
102
103			// for (j = 0; j < 3; j++)
104			// pos[j] += dt * (vel[j] + eta*rj[j]);
105
106			// atoms[i]->setPos( pos );
107
108			// if( atoms[i]->isDirectional() ){
109
110			// dAtom = (DirectionalAtom *)atoms[i];
111
112			// // get and convert the torque to body frame
113
114			// dAtom->getTrq( Tb );
115			// dAtom->lab2Body( Tb );
116
117			// // get the angular momentum, and propagate a half step
118
119			// dAtom->getJ( ji );
120
121			// for (j=0; j < 3; j++)
122			// ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
123
124			// // use the angular velocities to propagate the rotation matrix a
125			// // full time step
126
127			// dAtom->getA(A);
128			// dAtom->getI(I);
129
130			// // rotate about the x-axis
131			// angle = dt2 * ji[0] / I[0][0];
132			// this->rotate( 1, 2, angle, ji, A );
133
134			// // rotate about the y-axis
135			// angle = dt2 * ji[1] / I[1][1];
136			// this->rotate( 2, 0, angle, ji, A );
137
138			// // rotate about the z-axis
139			// angle = dt * ji[2] / I[2][2];
140			// this->rotate( 0, 1, angle, ji, A);
141
142			// // rotate about the y-axis
143			// angle = dt2 * ji[1] / I[1][1];
144			// this->rotate( 2, 0, angle, ji, A );
145
146			// // rotate about the x-axis
147			// angle = dt2 * ji[0] / I[0][0];
148			// this->rotate( 1, 2, angle, ji, A );
149
150			// dAtom->setJ( ji );
151			// dAtom->setA( A );
152			// }
153
154			// }
155
156			// // Scale the box after all the positions have been moved:
157
158			// scaleFactor = exp(dt*eta);
159
160			// if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
161			// sprintf( painCave.errMsg,
162			// "NPTi error: Attempting a Box scaling of more than 10 percent"
163			// " check your tauBarostat, as it is probably too small!\n"
164			// " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
165			// );
166			// painCave.isFatal = 1;
167			// simError();
168			// } else {
169			// info->scaleBox(exp(dt*eta));
170			// }
171
172
173			//new version of NPTi
174			int i, j, k;
175	gezelter	574	DirectionalAtom* dAtom;
176	gezelter	600	double Tb[3], ji[3];
177			double A[3][3], I[3][3];
178			double angle, mass;
179			double vel[3], pos[3], frc[3];
180
181	gezelter	574	double rj[3];
182			double instaTemp, instaPress, instaVol;
183	gezelter	611	double tt2, tb2, scaleFactor;
184	tim	763	double COM[3];
185	gezelter	574
186			tt2 = tauThermostat * tauThermostat;
187			tb2 = tauBarostat * tauBarostat;
188
189			instaTemp = tStats->getTemperature();
190			instaPress = tStats->getPressure();
191			instaVol = tStats->getVolume();
192
193	tim	763	tStats->getCOM(COM);
194
195			//evolve velocity half step
196			for( i=0; i<nAtoms; i++ ){
197	gezelter	574
198	gezelter	600	atoms[i]->getVel( vel );
199			atoms[i]->getFrc( frc );
200	gezelter	574
201	gezelter	600	mass = atoms[i]->getMass();
202	gezelter	574
203	gezelter	600	for (j=0; j < 3; j++) {
204	tim	763	// velocity half step (use chi from previous step here):
205			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta));
206
207	gezelter	600	}
208
209			atoms[i]->setVel( vel );
210	tim	763
211	gezelter	574	if( atoms[i]->isDirectional() ){
212
213			dAtom = (DirectionalAtom *)atoms[i];
214	tim	763
215	gezelter	574	// get and convert the torque to body frame
216
217	gezelter	600	dAtom->getTrq( Tb );
218	gezelter	574	dAtom->lab2Body( Tb );
219
220			// get the angular momentum, and propagate a half step
221
222	gezelter	600	dAtom->getJ( ji );
223
224			for (j=0; j < 3; j++)
225			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
226	gezelter	574
227			// use the angular velocities to propagate the rotation matrix a
228			// full time step
229	gezelter	600
230			dAtom->getA(A);
231			dAtom->getI(I);
232
233	gezelter	574	// rotate about the x-axis
234	gezelter	600	angle = dt2 * ji[0] / I[0][0];
235			this->rotate( 1, 2, angle, ji, A );
236
237	gezelter	574	// rotate about the y-axis
238	gezelter	600	angle = dt2 * ji[1] / I[1][1];
239			this->rotate( 2, 0, angle, ji, A );
240	gezelter	574
241			// rotate about the z-axis
242	gezelter	600	angle = dt * ji[2] / I[2][2];
243			this->rotate( 0, 1, angle, ji, A);
244	gezelter	574
245			// rotate about the y-axis
246	gezelter	600	angle = dt2 * ji[1] / I[1][1];
247			this->rotate( 2, 0, angle, ji, A );
248	gezelter	574
249			// rotate about the x-axis
250	gezelter	600	angle = dt2 * ji[0] / I[0][0];
251			this->rotate( 1, 2, angle, ji, A );
252	gezelter	574
253	gezelter	600	dAtom->setJ( ji );
254			dAtom->setA( A );
255	tim	763	}
256			}
257	gezelter	600
258	tim	763	// evolve chi and eta half step
259
260			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
261			eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convertNkBTtb2));
262
263			//calculate the integral of chidt
264			integralOfChidt += dt2*chi;
265
266			//save the old positions
267			for(i = 0; i < nAtoms; i++){
268			atoms[i]->getPos(pos);
269			for(j = 0; j < 3; j++)
270			oldPos[i*3 + j] = pos[j];
271	gezelter	574	}
272	tim	763
273			//the first estimation of r(t+dt) is equal to r(t)
274
275			for(k = 0; k < 4; k ++){
276	gezelter	611
277	tim	763	for(i =0 ; i < nAtoms; i++){
278
279			atoms[i]->getVel(vel);
280			atoms[i]->getPos(pos);
281
282			for(j = 0; j < 3; j++)
283			rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
284
285
286			//wrapVector(r(t)) = r(t)-R0
287			//info->wrapVector(rj);
288
289			for(j = 0; j < 3; j++)
290			pos[j] = oldPos[i3 + j] + dt(vel[j] + eta*rj[j]);
291
292			atoms[i]->setPos( pos );
293
294			}
295
296			}
297
298
299	gezelter	577	// Scale the box after all the positions have been moved:
300	gezelter	600
301	gezelter	611	scaleFactor = exp(dt*eta);
302
303	mmeineke	614	if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
304	gezelter	611	sprintf( painCave.errMsg,
305			"NPTi error: Attempting a Box scaling of more than 10 percent"
306			" check your tauBarostat, as it is probably too small!\n"
307			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
308			);
309			painCave.isFatal = 1;
310			simError();
311			} else {
312	tim	763	info->scaleBox(scaleFactor);
313			}
314	mmeineke	614
315	tim	763	//advance volume;
316			volume = volume * exp(dt*eta);
317	gezelter	574	}
318
319	tim	645	template<typename T> void NPTi<T>::moveB( void ){
320	gezelter	600
321	tim	763	/*
322	gezelter	600	int i, j;
323	gezelter	574	DirectionalAtom* dAtom;
324	gezelter	600	double Tb[3], ji[3];
325			double vel[3], frc[3];
326			double mass;
327
328	gezelter	574	double instaTemp, instaPress, instaVol;
329			double tt2, tb2;
330
331			tt2 = tauThermostat * tauThermostat;
332			tb2 = tauBarostat * tauBarostat;
333
334			instaTemp = tStats->getTemperature();
335			instaPress = tStats->getPressure();
336			instaVol = tStats->getVolume();
337
338			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
339	mmeineke	586	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
340			(p_convertNkBTtb2));
341	tim	763	integralOfChidt += dt2*chi;
342	gezelter	574
343			for( i=0; i<nAtoms; i++ ){
344	gezelter	600
345			atoms[i]->getVel( vel );
346			atoms[i]->getFrc( frc );
347
348			mass = atoms[i]->getMass();
349
350	gezelter	574	// velocity half step
351	gezelter	600	for (j=0; j < 3; j++)
352			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
353	gezelter	574
354	gezelter	600	atoms[i]->setVel( vel );
355
356	gezelter	574	if( atoms[i]->isDirectional() ){
357	gezelter	600
358	gezelter	574	dAtom = (DirectionalAtom *)atoms[i];
359	gezelter	600
360			// get and convert the torque to body frame
361
362			dAtom->getTrq( Tb );
363	gezelter	574	dAtom->lab2Body( Tb );
364	gezelter	600
365			// get the angular momentum, and propagate a half step
366
367			dAtom->getJ( ji );
368
369			for (j=0; j < 3; j++)
370			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
371
372			dAtom->setJ( ji );
373	gezelter	574	}
374			}
375	tim	763
376			*/
377
378			//new version of NPTi
379			int i, j, k;
380			DirectionalAtom* dAtom;
381			double Tb[3], ji[3];
382			double vel[3], frc[3];
383			double mass;
384
385			double instTemp, instPress, instVol;
386			double tt2, tb2;
387			double oldChi, prevChi;
388			double oldEta, preEta;
389
390			tt2 = tauThermostat * tauThermostat;
391			tb2 = tauBarostat * tauBarostat;
392
393
394			// Set things up for the iteration:
395
396			oldChi = chi;
397			oldEta = eta;
398
399			for( i=0; i<nAtoms; i++ ){
400
401			atoms[i]->getVel( vel );
402
403			for (j=0; j < 3; j++)
404			oldVel[3*i + j] = vel[j];
405
406			if( atoms[i]->isDirectional() ){
407
408			dAtom = (DirectionalAtom *)atoms[i];
409
410			dAtom->getJ( ji );
411
412			for (j=0; j < 3; j++)
413			oldJi[3*i + j] = ji[j];
414
415			}
416			}
417
418			// do the iteration:
419
420			instVol = tStats->getVolume();
421
422			for (k=0; k < 4; k++) {
423
424			instTemp = tStats->getTemperature();
425			instPress = tStats->getPressure();
426
427			// evolve chi another half step using the temperature at t + dt/2
428
429			prevChi = chi;
430			chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
431			(tauThermostat*tauThermostat);
432
433			preEta = eta;
434			eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) /
435			(p_convertNkBTtb2));
436
437
438			for( i=0; i<nAtoms; i++ ){
439
440			atoms[i]->getFrc( frc );
441			atoms[i]->getVel(vel);
442
443			mass = atoms[i]->getMass();
444
445			// velocity half step
446			for (j=0; j < 3; j++)
447			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j](chi + eta));
448
449			atoms[i]->setVel( vel );
450
451			if( atoms[i]->isDirectional() ){
452
453			dAtom = (DirectionalAtom *)atoms[i];
454
455			// get and convert the torque to body frame
456
457			dAtom->getTrq( Tb );
458			dAtom->lab2Body( Tb );
459
460			for (j=0; j < 3; j++)
461			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
462
463			dAtom->setJ( ji );
464			}
465			}
466
467			if (fabs(prevChi - chi) <= chiTolerance && fabs(preEta -eta) <= etaTolerance)
468			break;
469			}
470
471			//calculate integral of chida
472			integralOfChidt += dt2*chi;
473
474
475	gezelter	574	}
476
477	mmeineke	746	template<typename T> void NPTi<T>::resetIntegrator() {
478			chi = 0.0;
479			eta = 0.0;
480			}
481
482	tim	645	template<typename T> int NPTi<T>::readyCheck() {
483	tim	658
484			//check parent's readyCheck() first
485			if (T::readyCheck() == -1)
486			return -1;
487	gezelter	574
488			// First check to see if we have a target temperature.
489			// Not having one is fatal.
490
491			if (!have_target_temp) {
492			sprintf( painCave.errMsg,
493			"NPTi error: You can't use the NPTi integrator\n"
494			" without a targetTemp!\n"
495			);
496			painCave.isFatal = 1;
497			simError();
498			return -1;
499			}
500
501			if (!have_target_pressure) {
502			sprintf( painCave.errMsg,
503			"NPTi error: You can't use the NPTi integrator\n"
504			" without a targetPressure!\n"
505			);
506			painCave.isFatal = 1;
507			simError();
508			return -1;
509			}
510
511			// We must set tauThermostat.
512
513			if (!have_tau_thermostat) {
514			sprintf( painCave.errMsg,
515			"NPTi error: If you use the NPTi\n"
516			" integrator, you must set tauThermostat.\n");
517			painCave.isFatal = 1;
518			simError();
519			return -1;
520			}
521
522			// We must set tauBarostat.
523
524			if (!have_tau_barostat) {
525			sprintf( painCave.errMsg,
526			"NPTi error: If you use the NPTi\n"
527			" integrator, you must set tauBarostat.\n");
528			painCave.isFatal = 1;
529			simError();
530			return -1;
531			}
532
533	tim	763	if (!have_chi_tolerance) {
534			sprintf( painCave.errMsg,
535			"NPTi warning: setting chi tolerance to 1e-6\n");
536			chiTolerance = 1e-6;
537			have_chi_tolerance = 1;
538			painCave.isFatal = 0;
539			simError();
540			}
541
542			if (!have_eta_tolerance) {
543			sprintf( painCave.errMsg,
544			"NPTi warning: setting eta tolerance to 1e-6\n");
545			etaTolerance = 1e-6;
546			have_eta_tolerance = 1;
547			painCave.isFatal = 0;
548			simError();
549			}
550	gezelter	574	// We need NkBT a lot, so just set it here:
551
552	tim	763	NkBT = (double)Nparticles * kB * targetTemp;
553			fkBT = (double)info->ndf * kB * targetTemp;
554	gezelter	574
555			return 1;
556			}
557	tim	763
558			template<typename T> double NPTi<T>::getConservedQuantity(void){
559
560			double conservedQuantity;
561			double tb2;
562			double eta2;
563			double E_NPT;
564			double U;
565			double TS;
566			double PV;
567			double extra;
568
569			static double pre_U;
570			static double pre_TS;
571			static double pre_PV;
572			static double pre_extra;
573			static int hackCount = 0;
574
575			double delta_U;
576			double delta_TS;
577			double delta_PV;
578			double delta_extra;
579
580			U = tStats->getTotalE();
581
582			TS = fkBT *
583			(integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
584
585			PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
586
587			tb2 = tauBarostat * tauBarostat;
588			eta2 = eta * eta;
589
590			extra = (fkBT * tb2 * eta2 / 2.0 ) / eConvert;
591			/*
592			if(hackCount == 0){
593			pre_U = U;
594			pre_TS =TS;
595			pre_PV = PV;
596			pre_extra =extra;
597			hackCount ++;
598			}
599
600			delta_U = U - pre_U;
601			delta_TS = TS - pre_TS;
602			delta_PV = PV - pre_PV;
603			delta_extra = extra - pre_extra;
604			*/
605			cout.width(8);
606			cout.precision(8);
607
608
609			cout << info->getTime() << "\t"
610			<< chi << "\t"
611			<< eta << "\t"
612			<< U << "\t"
613			<< TS << "\t"
614			<< PV << "\t"
615			<< extra << "\t"
616			<< U+TS+PV+extra << endl;
617
618			/*
619			pre_U = U;
620			pre_TS =TS;
621			pre_PV = PV;
622			pre_extra =extra;
623
624
625			cout << info->getTime() << "\t"
626			<< U << "\t"
627			<< U+TS << "\t"
628			<< U+TS+PV << "\t"
629			<< U+TS+PV+extra << endl;
630			*/
631			conservedQuantity = U+TS+PV+extra;
632			return conservedQuantity;
633			}