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" 

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#endif

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

  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> NPTf<T>::~NPTf() {
  delete[] oldPos;
  delete[] oldVel;
  delete[] oldJi;
}

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

  // new version of NPTf
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];

  double 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;
  double COM[3];

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

  instaTemp = tStats->getTemperature();
  tStats->getPressureTensor(press);
  instaVol = tStats->getVolume();
  
  tStats->getCOM(COM);

  //calculate scale factor of veloity
  for (i = 0; i < 3; i++ ) {
    for (j = 0; j < 3; j++ ) {
      vScale[i][j] = eta[i][j];
      
      if (i == j) {
        vScale[i][j] += chi;          
      }               
    }
  }
  
  //evolve velocity half step
  for( i=0; i<nAtoms; i++ ){

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

    mass = atoms[i]->getMass();
    
    info->matVecMul3( vScale, vel, sc );

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

  // advance chi half step
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;

  // calculate the integral of chidt
  integralOfChidt += dt2*chi;

  // advance eta half step

  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);
      else
        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
    }
    
  //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]; 
      
      info->matVecMul3( eta, rj, 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:
  
  // 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 ){

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

  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double sc[3];
  double press[3][3], vScale[3][3];
  double oldChi, prevChi;
  double oldEta[3][3], prevEta[3][3], diffEta;
  
  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  // Set things up for the iteration:

  oldChi = chi;
  
  for(i = 0; i < 3; i++)
    for(j = 0; j < 3; j++)
      oldEta[i][j] = eta[i][j];

  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++) {
    
    instaTemp = tStats->getTemperature();
    tStats->getPressureTensor(press);

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

    prevChi = chi;
    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
    
    for(i = 0; i < 3; i++)
      for(j = 0; j < 3; j++)
        prevEta[i][j] = eta[i][j];

    //advance eta half step and calculate scale factor for velocity

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

      atoms[i]->getFrc( frc );
      atoms[i]->getVel(vel);
      
      mass = atoms[i]->getMass();
     
      for (j = 0; j < 3; j++)
        myVel[j] = oldVel[3*i + j];
      
      info->matVecMul3( vScale, myVel, sc );
      
      // velocity half step
      for (j=0; j < 3; j++) {
        // velocity half step  (use chi from previous step here):
        vel[j] = oldVel[3*i+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 );       
             
        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();
    }
    
    diffEta = 0;
    for(i = 0; i < 3; i++)
      diffEta += pow(prevEta[i][i] - eta[i][i], 2);    
    
    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance) 
      break;
  }

  //calculate integral of chidt
  integralOfChidt += dt2*chi;
  
}

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: This is the RAW number
  // of particles, so no subtraction or addition of constraints or
  // orientational degrees of freedom:
  
  NkBT = (double)Nparticles * kB * targetTemp;
  
  // fkBT is used because the thermostat operates on more degrees of freedom
  // than the barostat (when there are particles with orientational degrees
  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
  
  fkBT = (double)info->ndf * kB * targetTemp;

  return 1;
}

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

  double conservedQuantity;
  double Energy;
  double thermostat_kinetic;
  double thermostat_potential;
  double barostat_kinetic;
  double barostat_potential;
  double trEta;
  double a[3][3], b[3][3];

  Energy = tStats->getTotalE();

  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / 
    (2.0 * eConvert);

  thermostat_potential = fkBT* integralOfChidt / eConvert;

  info->transposeMat3(eta, a);
  info->matMul3(a, eta, b);
  trEta = info->matTrace3(b);

  barostat_kinetic = NkBT * tauBarostat * tauBarostat * trEta / 
    (2.0 * eConvert);
  
  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) / 
    eConvert;

  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
    barostat_kinetic + barostat_potential;
  
  cout.width(8);
  cout.precision(8);

  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << 
      "\t" << thermostat_potential << "\t" << barostat_kinetic << 
      "\t" << barostat_potential << "\t" << conservedQuantity << endl;

  return conservedQuantity; 
}
Revision:	778
Committed:	Fri Sep 19 20:00:27 2003 UTC (20 years, 9 months ago) by mmeineke
File size:	12519 byte(s)
Log Message:	added NPT base class. NPTi is up to date. NPTf is not.
#	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	gezelter	772	#ifdef IS_MPI
13			#include "mpiSimulation.hpp"
14			#endif
15	gezelter	576
16	gezelter	578	// Basic non-isotropic thermostating and barostating via the Melchionna
17	gezelter	576	// modification of the Hoover algorithm:
18			//
19			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
20			// Molec. Phys., 78, 533.
21			//
22			// and
23			//
24			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
25
26	tim	645	template<typename T> NPTf<T>::NPTf ( SimInfo theInfo, ForceFields the_ff):
27			T( theInfo, the_ff )
28	gezelter	576	{
29	gezelter	588	int i, j;
30	gezelter	576	chi = 0.0;
31	tim	763	integralOfChidt = 0.0;
32	gezelter	588
33			for(i = 0; i < 3; i++)
34	mmeineke	590	for (j = 0; j < 3; j++)
35	gezelter	588	eta[i][j] = 0.0;
36
37	gezelter	576	have_tau_thermostat = 0;
38			have_tau_barostat = 0;
39			have_target_temp = 0;
40			have_target_pressure = 0;
41	tim	767
42			have_chi_tolerance = 0;
43			have_eta_tolerance = 0;
44			have_pos_iter_tolerance = 0;
45
46			oldPos = new double[3*nAtoms];
47			oldVel = new double[3*nAtoms];
48			oldJi = new double[3*nAtoms];
49			#ifdef IS_MPI
50			Nparticles = mpiSim->getTotAtoms();
51			#else
52			Nparticles = theInfo->n_atoms;
53			#endif
54	gezelter	772
55	gezelter	576	}
56
57	tim	767	template<typename T> NPTf<T>::~NPTf() {
58			delete[] oldPos;
59			delete[] oldVel;
60			delete[] oldJi;
61			}
62
63	tim	645	template<typename T> void NPTf<T>::moveA() {
64	gezelter	772
65			// new version of NPTf
66	gezelter	600	int i, j, k;
67	gezelter	576	DirectionalAtom* dAtom;
68	gezelter	600	double Tb[3], ji[3];
69	mmeineke	778
70			double mass;
71	gezelter	600	double vel[3], pos[3], frc[3];
72
73			double rj[3];
74			double instaTemp, instaPress, instaVol;
75			double tt2, tb2;
76			double sc[3];
77			double eta2ij;
78	gezelter	588	double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
79	gezelter	617	double bigScale, smallScale, offDiagMax;
80	tim	767	double COM[3];
81	gezelter	576
82			tt2 = tauThermostat * tauThermostat;
83			tb2 = tauBarostat * tauBarostat;
84
85			instaTemp = tStats->getTemperature();
86	gezelter	577	tStats->getPressureTensor(press);
87	gezelter	576	instaVol = tStats->getVolume();
88
89	tim	767	tStats->getCOM(COM);
90	gezelter	588
91	tim	767	//calculate scale factor of veloity
92	gezelter	588	for (i = 0; i < 3; i++ ) {
93			for (j = 0; j < 3; j++ ) {
94	tim	767	vScale[i][j] = eta[i][j];
95
96	gezelter	588	if (i == j) {
97	tim	767	vScale[i][j] += chi;
98			}
99	gezelter	588	}
100			}
101	tim	767
102			//evolve velocity half step
103	gezelter	576	for( i=0; i<nAtoms; i++ ){
104	gezelter	600
105			atoms[i]->getVel( vel );
106			atoms[i]->getFrc( frc );
107
108			mass = atoms[i]->getMass();
109	gezelter	576
110	gezelter	600	info->matVecMul3( vScale, vel, sc );
111	tim	767
112			for (j=0; j < 3; j++) {
113	gezelter	772	// velocity half step
114	gezelter	600	vel[j] += dt2 * ((frc[j] / mass) * eConvert - sc[j]);
115			}
116	gezelter	576
117	gezelter	600	atoms[i]->setVel( vel );
118	gezelter	576
119			if( atoms[i]->isDirectional() ){
120
121			dAtom = (DirectionalAtom *)atoms[i];
122	tim	767
123	gezelter	576	// get and convert the torque to body frame
124
125	gezelter	600	dAtom->getTrq( Tb );
126	gezelter	576	dAtom->lab2Body( Tb );
127
128			// get the angular momentum, and propagate a half step
129
130	gezelter	600	dAtom->getJ( ji );
131
132			for (j=0; j < 3; j++)
133			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
134	gezelter	576
135	mmeineke	778	this->rotationPropagation( dAtom, ji );
136
137	gezelter	600	dAtom->setJ( ji );
138	tim	767	}
139	gezelter	576	}
140	tim	767
141			// advance chi half step
142			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
143
144	gezelter	772	// calculate the integral of chidt
145	tim	767	integralOfChidt += dt2*chi;
146
147	gezelter	772	// advance eta half step
148
149	tim	767	for(i = 0; i < 3; i ++)
150			for(j = 0; j < 3; j++){
151			if( i == j)
152			eta[i][j] += dt2 * instaVol *
153			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
154			else
155	gezelter	772	eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
156	tim	767	}
157
158			//save the old positions
159			for(i = 0; i < nAtoms; i++){
160			atoms[i]->getPos(pos);
161			for(j = 0; j < 3; j++)
162			oldPos[i*3 + j] = pos[j];
163			}
164	gezelter	600
165	tim	767	//the first estimation of r(t+dt) is equal to r(t)
166
167			for(k = 0; k < 4; k ++){
168
169			for(i =0 ; i < nAtoms; i++){
170
171			atoms[i]->getVel(vel);
172			atoms[i]->getPos(pos);
173
174			for(j = 0; j < 3; j++)
175			rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
176
177			info->matVecMul3( eta, rj, sc );
178
179			for(j = 0; j < 3; j++)
180			pos[j] = oldPos[i3 + j] + dt(vel[j] + sc[j]);
181
182			atoms[i]->setPos( pos );
183
184			}
185
186	gezelter	772	if (nConstrained) {
187			constrainA();
188			}
189	tim	767	}
190
191
192	gezelter	577	// Scale the box after all the positions have been moved:
193	gezelter	600
194	gezelter	578	// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
195			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
196	gezelter	600
197	gezelter	617	bigScale = 1.0;
198			smallScale = 1.0;
199			offDiagMax = 0.0;
200	gezelter	600
201	gezelter	578	for(i=0; i<3; i++){
202			for(j=0; j<3; j++){
203	gezelter	600
204	gezelter	588	// Calculate the matrix Product of the eta array (we only need
205			// the ij element right now):
206	gezelter	600
207	gezelter	588	eta2ij = 0.0;
208	gezelter	578	for(k=0; k<3; k++){
209	gezelter	588	eta2ij += eta[i][k] * eta[k][j];
210	gezelter	578	}
211	gezelter	588
212			scaleMat[i][j] = 0.0;
213			// identity matrix (see above):
214			if (i == j) scaleMat[i][j] = 1.0;
215			// Taylor expansion for the exponential truncated at second order:
216			scaleMat[i][j] += dteta[i][j] + 0.5dtdteta2ij;
217	gezelter	617
218			if (i != j)
219			if (fabs(scaleMat[i][j]) > offDiagMax)
220			offDiagMax = fabs(scaleMat[i][j]);
221	gezelter	578	}
222	gezelter	617
223			if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
224			if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
225	gezelter	578	}
226	gezelter	600
227	gezelter	617	if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
228			sprintf( painCave.errMsg,
229			"NPTf error: Attempting a Box scaling of more than 10 percent.\n"
230			" Check your tauBarostat, as it is probably too small!\n\n"
231			" scaleMat = [%lf\t%lf\t%lf]\n"
232			" [%lf\t%lf\t%lf]\n"
233			" [%lf\t%lf\t%lf]\n",
234			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
235			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
236			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
237			painCave.isFatal = 1;
238			simError();
239			} else if (offDiagMax > 0.1) {
240			sprintf( painCave.errMsg,
241			"NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
242			" Check your tauBarostat, as it is probably too small!\n\n"
243			" scaleMat = [%lf\t%lf\t%lf]\n"
244			" [%lf\t%lf\t%lf]\n"
245			" [%lf\t%lf\t%lf]\n",
246			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
247			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
248			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
249			painCave.isFatal = 1;
250			simError();
251			} else {
252			info->getBoxM(hm);
253			info->matMul3(hm, scaleMat, hmnew);
254			info->setBoxM(hmnew);
255			}
256	gezelter	577
257	gezelter	576	}
258
259	tim	645	template<typename T> void NPTf<T>::moveB( void ){
260	gezelter	600
261	gezelter	772	//new version of NPTf
262	tim	767	int i, j, k;
263	gezelter	576	DirectionalAtom* dAtom;
264	gezelter	600	double Tb[3], ji[3];
265	gezelter	772	double vel[3], myVel[3], frc[3];
266	gezelter	600	double mass;
267
268			double instaTemp, instaPress, instaVol;
269	gezelter	576	double tt2, tb2;
270	gezelter	600	double sc[3];
271	gezelter	588	double press[3][3], vScale[3][3];
272	tim	767	double oldChi, prevChi;
273	gezelter	772	double oldEta[3][3], prevEta[3][3], diffEta;
274	gezelter	576
275			tt2 = tauThermostat * tauThermostat;
276			tb2 = tauBarostat * tauBarostat;
277
278	tim	767	// Set things up for the iteration:
279
280			oldChi = chi;
281	gezelter	578
282	tim	767	for(i = 0; i < 3; i++)
283			for(j = 0; j < 3; j++)
284			oldEta[i][j] = eta[i][j];
285	gezelter	578
286	tim	767	for( i=0; i<nAtoms; i++ ){
287	gezelter	588
288	tim	767	atoms[i]->getVel( vel );
289	gezelter	588
290	tim	767	for (j=0; j < 3; j++)
291			oldVel[3*i + j] = vel[j];
292
293			if( atoms[i]->isDirectional() ){
294
295			dAtom = (DirectionalAtom *)atoms[i];
296
297			dAtom->getJ( ji );
298
299			for (j=0; j < 3; j++)
300			oldJi[3*i + j] = ji[j];
301
302	gezelter	588	}
303			}
304
305	tim	767	// do the iteration:
306	gezelter	578
307	tim	767	instaVol = tStats->getVolume();
308
309			for (k=0; k < 4; k++) {
310	gezelter	600
311	tim	767	instaTemp = tStats->getTemperature();
312			tStats->getPressureTensor(press);
313	gezelter	578
314	tim	767	// evolve chi another half step using the temperature at t + dt/2
315
316			prevChi = chi;
317			chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
318	gezelter	578
319	tim	767	for(i = 0; i < 3; i++)
320			for(j = 0; j < 3; j++)
321	gezelter	772	prevEta[i][j] = eta[i][j];
322	gezelter	600
323	tim	767	//advance eta half step and calculate scale factor for velocity
324	gezelter	772
325	tim	767	for(i = 0; i < 3; i ++)
326			for(j = 0; j < 3; j++){
327	gezelter	772	if( i == j) {
328	tim	767	eta[i][j] = oldEta[i][j] + dt2 * instaVol *
329	gezelter	772	(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
330	tim	767	vScale[i][j] = eta[i][j] + chi;
331	gezelter	772	} else {
332	tim	767	eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
333			vScale[i][j] = eta[i][j];
334			}
335	gezelter	772	}
336
337	tim	767	for( i=0; i<nAtoms; i++ ){
338
339			atoms[i]->getFrc( frc );
340			atoms[i]->getVel(vel);
341	gezelter	576
342	tim	767	mass = atoms[i]->getMass();
343	gezelter	772
344			for (j = 0; j < 3; j++)
345			myVel[j] = oldVel[3*i + j];
346	gezelter	576
347	gezelter	772	info->matVecMul3( vScale, myVel, sc );
348
349			// velocity half step
350	tim	767	for (j=0; j < 3; j++) {
351			// velocity half step (use chi from previous step here):
352			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass) * eConvert - sc[j]);
353			}
354	gezelter	576
355	tim	767	atoms[i]->setVel( vel );
356	gezelter	576
357	tim	767	if( atoms[i]->isDirectional() ){
358
359			dAtom = (DirectionalAtom *)atoms[i];
360
361			// get and convert the torque to body frame
362
363			dAtom->getTrq( Tb );
364			dAtom->lab2Body( Tb );
365
366			for (j=0; j < 3; j++)
367			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
368	gezelter	576
369	tim	767	dAtom->setJ( ji );
370			}
371			}
372	gezelter	600
373	gezelter	772	if (nConstrained) {
374			constrainB();
375			}
376	tim	767
377			diffEta = 0;
378			for(i = 0; i < 3; i++)
379	gezelter	772	diffEta += pow(prevEta[i][i] - eta[i][i], 2);
380	tim	767
381			if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
382			break;
383	gezelter	576	}
384	tim	767
385	gezelter	772	//calculate integral of chidt
386	tim	767	integralOfChidt += dt2*chi;
387
388	gezelter	576	}
389
390	mmeineke	746	template<typename T> void NPTf<T>::resetIntegrator() {
391			int i,j;
392
393			chi = 0.0;
394
395			for(i = 0; i < 3; i++)
396			for (j = 0; j < 3; j++)
397			eta[i][j] = 0.0;
398
399			}
400
401	tim	645	template<typename T> int NPTf<T>::readyCheck() {
402	tim	658
403			//check parent's readyCheck() first
404			if (T::readyCheck() == -1)
405			return -1;
406	gezelter	576
407			// First check to see if we have a target temperature.
408			// Not having one is fatal.
409
410			if (!have_target_temp) {
411			sprintf( painCave.errMsg,
412	gezelter	580	"NPTf error: You can't use the NPTf integrator\n"
413	gezelter	576	" without a targetTemp!\n"
414			);
415			painCave.isFatal = 1;
416			simError();
417			return -1;
418			}
419
420			if (!have_target_pressure) {
421			sprintf( painCave.errMsg,
422	gezelter	580	"NPTf error: You can't use the NPTf integrator\n"
423	gezelter	576	" without a targetPressure!\n"
424			);
425			painCave.isFatal = 1;
426			simError();
427			return -1;
428			}
429
430			// We must set tauThermostat.
431
432			if (!have_tau_thermostat) {
433			sprintf( painCave.errMsg,
434	gezelter	580	"NPTf error: If you use the NPTf\n"
435	gezelter	576	" integrator, you must set tauThermostat.\n");
436			painCave.isFatal = 1;
437			simError();
438			return -1;
439			}
440
441			// We must set tauBarostat.
442
443			if (!have_tau_barostat) {
444			sprintf( painCave.errMsg,
445	gezelter	580	"NPTf error: If you use the NPTf\n"
446	gezelter	576	" integrator, you must set tauBarostat.\n");
447			painCave.isFatal = 1;
448			simError();
449			return -1;
450			}
451
452	gezelter	772
453			// We need NkBT a lot, so just set it here: This is the RAW number
454			// of particles, so no subtraction or addition of constraints or
455			// orientational degrees of freedom:
456
457	tim	767	NkBT = (double)Nparticles * kB * targetTemp;
458	gezelter	772
459			// fkBT is used because the thermostat operates on more degrees of freedom
460			// than the barostat (when there are particles with orientational degrees
461			// of freedom). ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
462
463	tim	767	fkBT = (double)info->ndf * kB * targetTemp;
464	gezelter	576
465			return 1;
466			}
467	tim	763
468			template<typename T> double NPTf<T>::getConservedQuantity(void){
469
470			double conservedQuantity;
471	gezelter	772	double Energy;
472			double thermostat_kinetic;
473			double thermostat_potential;
474			double barostat_kinetic;
475			double barostat_potential;
476			double trEta;
477			double a[3][3], b[3][3];
478	tim	763
479	gezelter	772	Energy = tStats->getTotalE();
480	tim	763
481	gezelter	772	thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
482			(2.0 * eConvert);
483	tim	763
484	gezelter	772	thermostat_potential = fkBT* integralOfChidt / eConvert;
485	tim	763
486	gezelter	772	info->transposeMat3(eta, a);
487			info->matMul3(a, eta, b);
488			trEta = info->matTrace3(b);
489	tim	767
490	gezelter	772	barostat_kinetic = NkBT * tauBarostat * tauBarostat * trEta /
491			(2.0 * eConvert);
492
493			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
494			eConvert;
495	tim	767
496	gezelter	772	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
497			barostat_kinetic + barostat_potential;
498
499	tim	767	cout.width(8);
500			cout.precision(8);
501
502	gezelter	772	cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
503			"\t" << thermostat_potential << "\t" << barostat_kinetic <<
504			"\t" << barostat_potential << "\t" << conservedQuantity << endl;
505
506			return conservedQuantity;
507	tim	763	}