OOPSE/libmdtools/NPTf.cpp

#include <cmath>
#include "Atom.hpp"
#include "SRI.hpp"
#include "AbstractClasses.hpp"
#include "SimInfo.hpp"
#include "ForceFields.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"
#include "Integrator.hpp"
#include "simError.h" 


// Basic non-isotropic thermostating and barostating via the Melchionna
// modification of the Hoover algorithm:
//
//    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
//       Molec. Phys., 78, 533. 
//
//           and
// 
//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.

template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  int i, j;
  chi = 0.0;
  integralOfChidt = 0.0;

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

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

  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() {
  
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double A[3][3], I[3][3];
  double angle, mass;
  double vel[3], pos[3], frc[3];

  double rj[3];
  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double sc[3];
  double eta2ij;
  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
  double bigScale, smallScale, offDiagMax;
  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  (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);
      
      // 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  );    
    }    
  }

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

    }

  }  

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

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

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

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

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

  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double sc[3];
  double press[3][3], vScale[3][3];
  double oldChi, prevChi;
  double oldEta[3][3], preEta[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++)
        preEta[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];
        }
    }      

    //advance velocity half step
    for( i=0; i<nAtoms; i++ ){

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

      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 );
      }
    }

    
    diffEta = 0;
    for(i = 0; i < 3; i++)
      diffEta += pow(preEta[i][i] - eta[i][i], 2);    
    
    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance) 
      break;
  }

  //calculate integral of chida
  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:

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

  return 1;
}

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

  double conservedQuantity;
  double tb2;
  double trEta;  
  double U;
  double thermo;
  double integral;
  double baro;
  double PV;

  U = tStats->getTotalE();
  thermo = (fkBT * tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;

  tb2 = tauBarostat * tauBarostat;
  trEta = info->matTrace3(eta);
  baro = ((double)info->ndfTrans * kB * targetTemp * tb2 * trEta * trEta / 2.0) / eConvert;

  integral = ((double)(info->ndf + 1) * kB * targetTemp * integralOfChidt) /eConvert;

  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;


  cout.width(8);
  cout.precision(8);
  
  cout << info->getTime() << "\t" 
       << chi << "\t"
       << trEta << "\t"
       << U << "\t" 
       << thermo << "\t"
       << baro << "\t"
       << integral << "\t"
       << PV << "\t"
       << U+thermo+integral+PV+baro << endl;

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