OOPSE/libmdtools/NPTxym.cpp

#include <cmath>
#include "Atom.hpp"
#include "Molecule.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 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.

//  The NPTxym variant scales the molecular center-of-mass coordinates
//  instead of the atomic coordinates

template<typename T> NPTxym<T>::NPTxym ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  chi = 0.0;
  eta = 0.0;
  etaX = 0.0;
  etaY = 0.0;
  have_tau_thermostat = 0;
  have_tau_barostat = 0;
  have_target_temp = 0;
  have_target_pressure = 0;
}

template<typename T> void NPTxym<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, instaPressZ, instaVol;
  double tt2, tb2;
  double scaleFactor, scaleFactorZ, scaleFactorX, scaleFactorY,  bigFactor;
  double hm[3][3], hmnew[3][3], scaleMat[3][3];
  double scF3;


  int nInMol;
  double rc[3];

  nMols = info->n_mol;
  myMolecules = info->molecules;

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

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  instaPressX = tStats->getPressureX();
  instaPressY = tstats->getPressureY();
  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));
  etaX += dt2 * ( instaVol * (instaPressX - targetPressure) / 
                 (p_convert*NkBT*tb2));
  etaY += dt2 * ( instaVol * (instaPressY - targetPressure) / 
                 (p_convert*NkBT*tb2));

  for( i = 0; i < nMols; i++) {

    myMolecules[i].getCOM(rc);

    nInMol = myMolecules[i].getNAtoms();
    myAtoms = myMolecules[i].getMyAtoms();

    // find the minimum image coordinates of the molecular centers of mass:

    info->wrapVector(rc);

    for (j = 0; j < nInMol; j++) {

      if(myAtoms[j] != NULL) {

        myAtoms[j]->getVel( vel );
        myAtoms[j]->getPos( pos );
        myAtoms[j]->getFrc( frc );

        mass = myAtoms[j]->getMass();

        for (k=0; k < 3; k++)
          vel[k] += dt2 * ((frc[k] / mass ) * eConvert - vel[k]*(chi+eta));

        myAtoms[j]->setVel( vel );

        for (k = 0; k < 3; k++) 
          pos[k] += dt * (vel[k] + eta*rc[k]);

        myAtoms[j]->setPos( pos );

        if( myAtoms[j]->isDirectional() ){

          dAtom = (DirectionalAtom *)myAtoms[j];
          
          // 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 (k=0; k < 3; k++) 
            ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*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:
  
  scaleFactorX = exp(dt*etaX);
  scaleFactorY = exp(dt*etaY);
  scaleFactor = exp(dt*eta);
  
  bigFactor = abs( 1.0 - scaleFactorX );
  if( abs(1.0 - scaleFactorY) > bigFactor )
    bigFactor = abs(1.0 - scaleFactorY);
  if( abs(1.0 - scaleFactor) > bigFactor )
    bigFactor = abs(1.0 - scaleFactor);
  
  if (bigFactor > 0.1) {
    sprintf( painCave.errMsg,
             "NPTxym 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"
             " etaX = %lf, scaleFactorX = %lf\n", 
             " etaY = %lf, scaleFactorY = %lf\n", 
             eta, scaleFactor,
             etaX, scaleFactorX,
             etaY, scaleFactorY
             );
    painCave.isFatal = 1;
    simError();
  } else {         
    
    for(i=0;i<3;i++)
      for(j=0;j<3;j++)
        scaleBox[i][j] = 0.0;

    scF3 = scaleFactor * scaleFactor * scaleFactor;
    scaleFactorZ = scF3 / (scaleFactorX * scaleFactorY);

    scaleBox[0][0] = scaleFactorX;
    scaleBox[1][1] = scaleFactorY;
    scaleBox[2][2] = scaleFactorZ;

    info->getBoxM( hm );
    info->matMul3( hm, scaleBox, hmnew );

    info->setBoxM( hmnew );
  }
}

template<typename T> void NPTxym<T>::moveB( void ){
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double vel[3], frc[3];
  double mass;

  double instaTemp, instaPress, instaPressX, instaPressY, instaVol;
  double tt2, tb2;
 
  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  instaPressX = tStats->getPressureX();
  instaPressY = tstats->getPressureY();
  instaVol = tStats->getVolume();

  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  eta += dt2 * ( instaVol * (instaPress - targetPressure) / 
                 (p_convert*NkBT*tb2));
  etaX += dt2 * ( instaVol * (instaPressX - targetPressure) / 
                 (p_convert*NkBT*tb2));
  etaY += dt2 * ( instaVol * (instaPressY - targetPressure) / 
                 (p_convert*NkBT*tb2));

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

template<typename T> void NPTxym<T>::resetIntegrator() {
  chi = 0.0;
  eta = 0.0;
  etaX = 0.0;
  etaY = 0.0;  
}

template<typename T> int NPTxym<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,
             "NPTxym error: You can't use the NPTxym integrator\n"
             "   without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

  if (!have_target_pressure) {
    sprintf( painCave.errMsg,
             "NPTxym error: You can't use the NPTxym integrator\n"
             "   without a targetPressure!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }
  
  // We must set tauThermostat.
   
  if (!have_tau_thermostat) {
    sprintf( painCave.errMsg,
             "NPTxym error: If you use the NPTxym\n"
             "   integrator, you must set tauThermostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }    

  // We must set tauBarostat.
   
  if (!have_tau_barostat) {
    sprintf( painCave.errMsg,
             "NPTxym error: If you use the NPTxym\n"
             "   integrator, you must set tauBarostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }    

  // We need NkBT a lot, so just set it here:

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

  return 1;
}
Revision:	755
Committed:	Tue Sep 9 20:35:25 2003 UTC (20 years, 10 months ago) by mmeineke
File size:	8536 byte(s)
Log Message:	updated the ChangeLog added two new NPT integrators, they still need work.
#	User	Rev	Content
1	mmeineke	755	#include <cmath>
2			#include "Atom.hpp"
3			#include "Molecule.hpp"
4			#include "SRI.hpp"
5			#include "AbstractClasses.hpp"
6			#include "SimInfo.hpp"
7			#include "ForceFields.hpp"
8			#include "Thermo.hpp"
9			#include "ReadWrite.hpp"
10			#include "Integrator.hpp"
11			#include "simError.h"
12
13
14			// Basic isotropic thermostating and barostating via the Melchionna
15			// modification of the Hoover algorithm:
16			//
17			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
18			// Molec. Phys., 78, 533.
19			//
20			// and
21			//
22			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
23
24			// The NPTxym variant scales the molecular center-of-mass coordinates
25			// instead of the atomic coordinates
26
27			template<typename T> NPTxym<T>::NPTxym ( SimInfo theInfo, ForceFields the_ff):
28			T( theInfo, the_ff )
29			{
30			chi = 0.0;
31			eta = 0.0;
32			etaX = 0.0;
33			etaY = 0.0;
34			have_tau_thermostat = 0;
35			have_tau_barostat = 0;
36			have_target_temp = 0;
37			have_target_pressure = 0;
38			}
39
40			template<typename T> void NPTxym<T>::moveA() {
41
42			int i, j, k;
43			DirectionalAtom* dAtom;
44			double Tb[3], ji[3];
45			double A[3][3], I[3][3];
46			double angle, mass;
47			double vel[3], pos[3], frc[3];
48
49			double rj[3];
50			double instaTemp, instaPress, instaPressZ, instaVol;
51			double tt2, tb2;
52			double scaleFactor, scaleFactorZ, scaleFactorX, scaleFactorY, bigFactor;
53			double hm[3][3], hmnew[3][3], scaleMat[3][3];
54			double scF3;
55
56
57			int nInMol;
58			double rc[3];
59
60			nMols = info->n_mol;
61			myMolecules = info->molecules;
62
63			tt2 = tauThermostat * tauThermostat;
64			tb2 = tauBarostat * tauBarostat;
65
66			instaTemp = tStats->getTemperature();
67			instaPress = tStats->getPressure();
68			instaPressX = tStats->getPressureX();
69			instaPressY = tstats->getPressureY();
70			instaVol = tStats->getVolume();
71
72			// first evolve chi a half step
73
74			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
75			eta += dt2 * ( instaVol * (instaPress - targetPressure) /
76			(p_convertNkBTtb2));
77			etaX += dt2 * ( instaVol * (instaPressX - targetPressure) /
78			(p_convertNkBTtb2));
79			etaY += dt2 * ( instaVol * (instaPressY - targetPressure) /
80			(p_convertNkBTtb2));
81
82			for( i = 0; i < nMols; i++) {
83
84			myMolecules[i].getCOM(rc);
85
86			nInMol = myMolecules[i].getNAtoms();
87			myAtoms = myMolecules[i].getMyAtoms();
88
89			// find the minimum image coordinates of the molecular centers of mass:
90
91			info->wrapVector(rc);
92
93			for (j = 0; j < nInMol; j++) {
94
95			if(myAtoms[j] != NULL) {
96
97			myAtoms[j]->getVel( vel );
98			myAtoms[j]->getPos( pos );
99			myAtoms[j]->getFrc( frc );
100
101			mass = myAtoms[j]->getMass();
102
103			for (k=0; k < 3; k++)
104			vel[k] += dt2 * ((frc[k] / mass ) * eConvert - vel[k]*(chi+eta));
105
106			myAtoms[j]->setVel( vel );
107
108			for (k = 0; k < 3; k++)
109			pos[k] += dt * (vel[k] + eta*rc[k]);
110
111			myAtoms[j]->setPos( pos );
112
113			if( myAtoms[j]->isDirectional() ){
114
115			dAtom = (DirectionalAtom *)myAtoms[j];
116
117			// get and convert the torque to body frame
118
119			dAtom->getTrq( Tb );
120			dAtom->lab2Body( Tb );
121
122			// get the angular momentum, and propagate a half step
123
124			dAtom->getJ( ji );
125
126			for (k=0; k < 3; k++)
127			ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
128
129			// use the angular velocities to propagate the rotation matrix a
130			// full time step
131
132			dAtom->getA(A);
133			dAtom->getI(I);
134
135			// rotate about the x-axis
136			angle = dt2 * ji[0] / I[0][0];
137			this->rotate( 1, 2, angle, ji, A );
138
139			// rotate about the y-axis
140			angle = dt2 * ji[1] / I[1][1];
141			this->rotate( 2, 0, angle, ji, A );
142
143			// rotate about the z-axis
144			angle = dt * ji[2] / I[2][2];
145			this->rotate( 0, 1, angle, ji, A);
146
147			// rotate about the y-axis
148			angle = dt2 * ji[1] / I[1][1];
149			this->rotate( 2, 0, angle, ji, A );
150
151			// rotate about the x-axis
152			angle = dt2 * ji[0] / I[0][0];
153			this->rotate( 1, 2, angle, ji, A );
154
155			dAtom->setJ( ji );
156			dAtom->setA( A );
157			}
158			}
159			}
160			}
161
162			// Scale the box after all the positions have been moved:
163
164			scaleFactorX = exp(dt*etaX);
165			scaleFactorY = exp(dt*etaY);
166			scaleFactor = exp(dt*eta);
167
168			bigFactor = abs( 1.0 - scaleFactorX );
169			if( abs(1.0 - scaleFactorY) > bigFactor )
170			bigFactor = abs(1.0 - scaleFactorY);
171			if( abs(1.0 - scaleFactor) > bigFactor )
172			bigFactor = abs(1.0 - scaleFactor);
173
174			if (bigFactor > 0.1) {
175			sprintf( painCave.errMsg,
176			"NPTxym error: Attempting a Box scaling of more than 10 percent"
177			" check your tauBarostat, as it is probably too small!\n"
178			" eta = %lf, scaleFactor = %lf\n"
179			" etaX = %lf, scaleFactorX = %lf\n",
180			" etaY = %lf, scaleFactorY = %lf\n",
181			eta, scaleFactor,
182			etaX, scaleFactorX,
183			etaY, scaleFactorY
184			);
185			painCave.isFatal = 1;
186			simError();
187			} else {
188
189			for(i=0;i<3;i++)
190			for(j=0;j<3;j++)
191			scaleBox[i][j] = 0.0;
192
193			scF3 = scaleFactor * scaleFactor * scaleFactor;
194			scaleFactorZ = scF3 / (scaleFactorX * scaleFactorY);
195
196			scaleBox[0][0] = scaleFactorX;
197			scaleBox[1][1] = scaleFactorY;
198			scaleBox[2][2] = scaleFactorZ;
199
200			info->getBoxM( hm );
201			info->matMul3( hm, scaleBox, hmnew );
202
203			info->setBoxM( hmnew );
204			}
205			}
206
207			template<typename T> void NPTxym<T>::moveB( void ){
208			int i, j;
209			DirectionalAtom* dAtom;
210			double Tb[3], ji[3];
211			double vel[3], frc[3];
212			double mass;
213
214			double instaTemp, instaPress, instaPressX, instaPressY, instaVol;
215			double tt2, tb2;
216
217			tt2 = tauThermostat * tauThermostat;
218			tb2 = tauBarostat * tauBarostat;
219
220			instaTemp = tStats->getTemperature();
221			instaPress = tStats->getPressure();
222			instaPressX = tStats->getPressureX();
223			instaPressY = tstats->getPressureY();
224			instaVol = tStats->getVolume();
225
226			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
227			eta += dt2 * ( instaVol * (instaPress - targetPressure) /
228			(p_convertNkBTtb2));
229			etaX += dt2 * ( instaVol * (instaPressX - targetPressure) /
230			(p_convertNkBTtb2));
231			etaY += dt2 * ( instaVol * (instaPressY - targetPressure) /
232			(p_convertNkBTtb2));
233
234			for( i=0; i<nAtoms; i++ ){
235
236			atoms[i]->getVel( vel );
237			atoms[i]->getFrc( frc );
238
239			mass = atoms[i]->getMass();
240
241			// velocity half step
242			for (j=0; j < 3; j++)
243			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
244
245			atoms[i]->setVel( vel );
246
247			if( atoms[i]->isDirectional() ){
248
249			dAtom = (DirectionalAtom *)atoms[i];
250
251			// get and convert the torque to body frame
252
253			dAtom->getTrq( Tb );
254			dAtom->lab2Body( Tb );
255
256			// get the angular momentum, and propagate a half step
257
258			dAtom->getJ( ji );
259
260			for (j=0; j < 3; j++)
261			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
262
263			dAtom->setJ( ji );
264			}
265			}
266			}
267
268			template<typename T> void NPTxym<T>::resetIntegrator() {
269			chi = 0.0;
270			eta = 0.0;
271			etaX = 0.0;
272			etaY = 0.0;
273			}
274
275			template<typename T> int NPTxym<T>::readyCheck() {
276
277			//check parent's readyCheck() first
278			if (T::readyCheck() == -1)
279			return -1;
280
281			// First check to see if we have a target temperature.
282			// Not having one is fatal.
283
284			if (!have_target_temp) {
285			sprintf( painCave.errMsg,
286			"NPTxym error: You can't use the NPTxym integrator\n"
287			" without a targetTemp!\n"
288			);
289			painCave.isFatal = 1;
290			simError();
291			return -1;
292			}
293
294			if (!have_target_pressure) {
295			sprintf( painCave.errMsg,
296			"NPTxym error: You can't use the NPTxym integrator\n"
297			" without a targetPressure!\n"
298			);
299			painCave.isFatal = 1;
300			simError();
301			return -1;
302			}
303
304			// We must set tauThermostat.
305
306			if (!have_tau_thermostat) {
307			sprintf( painCave.errMsg,
308			"NPTxym error: If you use the NPTxym\n"
309			" integrator, you must set tauThermostat.\n");
310			painCave.isFatal = 1;
311			simError();
312			return -1;
313			}
314
315			// We must set tauBarostat.
316
317			if (!have_tau_barostat) {
318			sprintf( painCave.errMsg,
319			"NPTxym error: If you use the NPTxym\n"
320			" integrator, you must set tauBarostat.\n");
321			painCave.isFatal = 1;
322			simError();
323			return -1;
324			}
325
326			// We need NkBT a lot, so just set it here:
327
328			NkBT = (double)info->ndf * kB * targetTemp;
329
330			return 1;
331			}