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

    }

    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:	772
Committed:	Fri Sep 19 16:01:07 2003 UTC (22 years, 1 month ago) by gezelter
Original Path:	*trunk/OOPSE/libmdtools/NPTf.cpp*
File size:	13271 byte(s)
Log Message:	fixed bugs in NPTf, found (nearly) conserved quantities for both NPTi and NPTf
#	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			double A[3][3], I[3][3];
70			double angle, mass;
71			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			// use the angular velocities to propagate the rotation matrix a
136			// full time step
137	gezelter	600
138			dAtom->getA(A);
139			dAtom->getI(I);
140
141	gezelter	576	// rotate about the x-axis
142	gezelter	600	angle = dt2 * ji[0] / I[0][0];
143			this->rotate( 1, 2, angle, ji, A );
144
145	gezelter	576	// rotate about the y-axis
146	gezelter	600	angle = dt2 * ji[1] / I[1][1];
147			this->rotate( 2, 0, angle, ji, A );
148	gezelter	576
149			// rotate about the z-axis
150	gezelter	600	angle = dt * ji[2] / I[2][2];
151			this->rotate( 0, 1, angle, ji, A);
152	gezelter	576
153			// rotate about the y-axis
154	gezelter	600	angle = dt2 * ji[1] / I[1][1];
155			this->rotate( 2, 0, angle, ji, A );
156	gezelter	576
157			// rotate about the x-axis
158	gezelter	600	angle = dt2 * ji[0] / I[0][0];
159			this->rotate( 1, 2, angle, ji, A );
160	gezelter	576
161	gezelter	600	dAtom->setJ( ji );
162			dAtom->setA( A );
163	tim	767	}
164	gezelter	576	}
165	tim	767
166			// advance chi half step
167			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
168
169	gezelter	772	// calculate the integral of chidt
170	tim	767	integralOfChidt += dt2*chi;
171
172	gezelter	772	// advance eta half step
173
174	tim	767	for(i = 0; i < 3; i ++)
175			for(j = 0; j < 3; j++){
176			if( i == j)
177			eta[i][j] += dt2 * instaVol *
178			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
179			else
180	gezelter	772	eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
181	tim	767	}
182
183			//save the old positions
184			for(i = 0; i < nAtoms; i++){
185			atoms[i]->getPos(pos);
186			for(j = 0; j < 3; j++)
187			oldPos[i*3 + j] = pos[j];
188			}
189	gezelter	600
190	tim	767	//the first estimation of r(t+dt) is equal to r(t)
191
192			for(k = 0; k < 4; k ++){
193
194			for(i =0 ; i < nAtoms; i++){
195
196			atoms[i]->getVel(vel);
197			atoms[i]->getPos(pos);
198
199			for(j = 0; j < 3; j++)
200			rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
201
202			info->matVecMul3( eta, rj, sc );
203
204			for(j = 0; j < 3; j++)
205			pos[j] = oldPos[i3 + j] + dt(vel[j] + sc[j]);
206
207			atoms[i]->setPos( pos );
208
209			}
210
211	gezelter	772	if (nConstrained) {
212			constrainA();
213			}
214	tim	767	}
215
216
217	gezelter	577	// Scale the box after all the positions have been moved:
218	gezelter	600
219	gezelter	578	// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
220			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
221	gezelter	600
222	gezelter	617	bigScale = 1.0;
223			smallScale = 1.0;
224			offDiagMax = 0.0;
225	gezelter	600
226	gezelter	578	for(i=0; i<3; i++){
227			for(j=0; j<3; j++){
228	gezelter	600
229	gezelter	588	// Calculate the matrix Product of the eta array (we only need
230			// the ij element right now):
231	gezelter	600
232	gezelter	588	eta2ij = 0.0;
233	gezelter	578	for(k=0; k<3; k++){
234	gezelter	588	eta2ij += eta[i][k] * eta[k][j];
235	gezelter	578	}
236	gezelter	588
237			scaleMat[i][j] = 0.0;
238			// identity matrix (see above):
239			if (i == j) scaleMat[i][j] = 1.0;
240			// Taylor expansion for the exponential truncated at second order:
241			scaleMat[i][j] += dteta[i][j] + 0.5dtdteta2ij;
242	gezelter	617
243			if (i != j)
244			if (fabs(scaleMat[i][j]) > offDiagMax)
245			offDiagMax = fabs(scaleMat[i][j]);
246	gezelter	578	}
247	gezelter	617
248			if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
249			if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
250	gezelter	578	}
251	gezelter	600
252	gezelter	617	if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
253			sprintf( painCave.errMsg,
254			"NPTf error: Attempting a Box scaling of more than 10 percent.\n"
255			" Check your tauBarostat, as it is probably too small!\n\n"
256			" scaleMat = [%lf\t%lf\t%lf]\n"
257			" [%lf\t%lf\t%lf]\n"
258			" [%lf\t%lf\t%lf]\n",
259			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
260			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
261			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
262			painCave.isFatal = 1;
263			simError();
264			} else if (offDiagMax > 0.1) {
265			sprintf( painCave.errMsg,
266			"NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
267			" Check your tauBarostat, as it is probably too small!\n\n"
268			" scaleMat = [%lf\t%lf\t%lf]\n"
269			" [%lf\t%lf\t%lf]\n"
270			" [%lf\t%lf\t%lf]\n",
271			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
272			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
273			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
274			painCave.isFatal = 1;
275			simError();
276			} else {
277			info->getBoxM(hm);
278			info->matMul3(hm, scaleMat, hmnew);
279			info->setBoxM(hmnew);
280			}
281	gezelter	577
282	gezelter	576	}
283
284	tim	645	template<typename T> void NPTf<T>::moveB( void ){
285	gezelter	600
286	gezelter	772	//new version of NPTf
287	tim	767	int i, j, k;
288	gezelter	576	DirectionalAtom* dAtom;
289	gezelter	600	double Tb[3], ji[3];
290	gezelter	772	double vel[3], myVel[3], frc[3];
291	gezelter	600	double mass;
292
293			double instaTemp, instaPress, instaVol;
294	gezelter	576	double tt2, tb2;
295	gezelter	600	double sc[3];
296	gezelter	588	double press[3][3], vScale[3][3];
297	tim	767	double oldChi, prevChi;
298	gezelter	772	double oldEta[3][3], prevEta[3][3], diffEta;
299	gezelter	576
300			tt2 = tauThermostat * tauThermostat;
301			tb2 = tauBarostat * tauBarostat;
302
303	tim	767	// Set things up for the iteration:
304
305			oldChi = chi;
306	gezelter	578
307	tim	767	for(i = 0; i < 3; i++)
308			for(j = 0; j < 3; j++)
309			oldEta[i][j] = eta[i][j];
310	gezelter	578
311	tim	767	for( i=0; i<nAtoms; i++ ){
312	gezelter	588
313	tim	767	atoms[i]->getVel( vel );
314	gezelter	588
315	tim	767	for (j=0; j < 3; j++)
316			oldVel[3*i + j] = vel[j];
317
318			if( atoms[i]->isDirectional() ){
319
320			dAtom = (DirectionalAtom *)atoms[i];
321
322			dAtom->getJ( ji );
323
324			for (j=0; j < 3; j++)
325			oldJi[3*i + j] = ji[j];
326
327	gezelter	588	}
328			}
329
330	tim	767	// do the iteration:
331	gezelter	578
332	tim	767	instaVol = tStats->getVolume();
333
334			for (k=0; k < 4; k++) {
335	gezelter	600
336	tim	767	instaTemp = tStats->getTemperature();
337			tStats->getPressureTensor(press);
338	gezelter	578
339	tim	767	// evolve chi another half step using the temperature at t + dt/2
340
341			prevChi = chi;
342			chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
343	gezelter	578
344	tim	767	for(i = 0; i < 3; i++)
345			for(j = 0; j < 3; j++)
346	gezelter	772	prevEta[i][j] = eta[i][j];
347	gezelter	600
348	tim	767	//advance eta half step and calculate scale factor for velocity
349	gezelter	772
350	tim	767	for(i = 0; i < 3; i ++)
351			for(j = 0; j < 3; j++){
352	gezelter	772	if( i == j) {
353	tim	767	eta[i][j] = oldEta[i][j] + dt2 * instaVol *
354	gezelter	772	(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
355	tim	767	vScale[i][j] = eta[i][j] + chi;
356	gezelter	772	} else {
357	tim	767	eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
358			vScale[i][j] = eta[i][j];
359			}
360	gezelter	772	}
361
362	tim	767	for( i=0; i<nAtoms; i++ ){
363
364			atoms[i]->getFrc( frc );
365			atoms[i]->getVel(vel);
366	gezelter	576
367	tim	767	mass = atoms[i]->getMass();
368	gezelter	772
369			for (j = 0; j < 3; j++)
370			myVel[j] = oldVel[3*i + j];
371	gezelter	576
372	gezelter	772	info->matVecMul3( vScale, myVel, sc );
373
374			// velocity half step
375	tim	767	for (j=0; j < 3; j++) {
376			// velocity half step (use chi from previous step here):
377			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass) * eConvert - sc[j]);
378			}
379	gezelter	576
380	tim	767	atoms[i]->setVel( vel );
381	gezelter	576
382	tim	767	if( atoms[i]->isDirectional() ){
383
384			dAtom = (DirectionalAtom *)atoms[i];
385
386			// get and convert the torque to body frame
387
388			dAtom->getTrq( Tb );
389			dAtom->lab2Body( Tb );
390
391			for (j=0; j < 3; j++)
392			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
393	gezelter	576
394	tim	767	dAtom->setJ( ji );
395			}
396			}
397	gezelter	600
398	gezelter	772	if (nConstrained) {
399			constrainB();
400			}
401	tim	767
402			diffEta = 0;
403			for(i = 0; i < 3; i++)
404	gezelter	772	diffEta += pow(prevEta[i][i] - eta[i][i], 2);
405	tim	767
406			if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance)
407			break;
408	gezelter	576	}
409	tim	767
410	gezelter	772	//calculate integral of chidt
411	tim	767	integralOfChidt += dt2*chi;
412
413	gezelter	576	}
414
415	mmeineke	746	template<typename T> void NPTf<T>::resetIntegrator() {
416			int i,j;
417
418			chi = 0.0;
419
420			for(i = 0; i < 3; i++)
421			for (j = 0; j < 3; j++)
422			eta[i][j] = 0.0;
423
424			}
425
426	tim	645	template<typename T> int NPTf<T>::readyCheck() {
427	tim	658
428			//check parent's readyCheck() first
429			if (T::readyCheck() == -1)
430			return -1;
431	gezelter	576
432			// First check to see if we have a target temperature.
433			// Not having one is fatal.
434
435			if (!have_target_temp) {
436			sprintf( painCave.errMsg,
437	gezelter	580	"NPTf error: You can't use the NPTf integrator\n"
438	gezelter	576	" without a targetTemp!\n"
439			);
440			painCave.isFatal = 1;
441			simError();
442			return -1;
443			}
444
445			if (!have_target_pressure) {
446			sprintf( painCave.errMsg,
447	gezelter	580	"NPTf error: You can't use the NPTf integrator\n"
448	gezelter	576	" without a targetPressure!\n"
449			);
450			painCave.isFatal = 1;
451			simError();
452			return -1;
453			}
454
455			// We must set tauThermostat.
456
457			if (!have_tau_thermostat) {
458			sprintf( painCave.errMsg,
459	gezelter	580	"NPTf error: If you use the NPTf\n"
460	gezelter	576	" integrator, you must set tauThermostat.\n");
461			painCave.isFatal = 1;
462			simError();
463			return -1;
464			}
465
466			// We must set tauBarostat.
467
468			if (!have_tau_barostat) {
469			sprintf( painCave.errMsg,
470	gezelter	580	"NPTf error: If you use the NPTf\n"
471	gezelter	576	" integrator, you must set tauBarostat.\n");
472			painCave.isFatal = 1;
473			simError();
474			return -1;
475			}
476
477	gezelter	772
478			// We need NkBT a lot, so just set it here: This is the RAW number
479			// of particles, so no subtraction or addition of constraints or
480			// orientational degrees of freedom:
481
482	tim	767	NkBT = (double)Nparticles * kB * targetTemp;
483	gezelter	772
484			// fkBT is used because the thermostat operates on more degrees of freedom
485			// than the barostat (when there are particles with orientational degrees
486			// of freedom). ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
487
488	tim	767	fkBT = (double)info->ndf * kB * targetTemp;
489	gezelter	576
490			return 1;
491			}
492	tim	763
493			template<typename T> double NPTf<T>::getConservedQuantity(void){
494
495			double conservedQuantity;
496	gezelter	772	double Energy;
497			double thermostat_kinetic;
498			double thermostat_potential;
499			double barostat_kinetic;
500			double barostat_potential;
501			double trEta;
502			double a[3][3], b[3][3];
503	tim	763
504	gezelter	772	Energy = tStats->getTotalE();
505	tim	763
506	gezelter	772	thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
507			(2.0 * eConvert);
508	tim	763
509	gezelter	772	thermostat_potential = fkBT* integralOfChidt / eConvert;
510	tim	763
511	gezelter	772	info->transposeMat3(eta, a);
512			info->matMul3(a, eta, b);
513			trEta = info->matTrace3(b);
514	tim	767
515	gezelter	772	barostat_kinetic = NkBT * tauBarostat * tauBarostat * trEta /
516			(2.0 * eConvert);
517
518			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
519			eConvert;
520	tim	767
521	gezelter	772	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
522			barostat_kinetic + barostat_potential;
523
524	tim	767	cout.width(8);
525			cout.precision(8);
526
527	gezelter	772	cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
528			"\t" << thermostat_potential << "\t" << barostat_kinetic <<
529			"\t" << barostat_potential << "\t" << conservedQuantity << endl;
530
531			return conservedQuantity;
532	tim	763	}