OOPSE/libmdtools/NPTfm.cpp

#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 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.

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

NPTfm::NPTfm ( SimInfo *theInfo, ForceFields* the_ff):
  Integrator( theInfo, the_ff )
{
  int i, j;
  chi = 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;
}

void NPTfm::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];

  int nInMol;
  double rc[3];

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

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

  instaTemp = tStats->getTemperature();
  tStats->getPressureTensor(press);
  instaVol = tStats->getVolume();
   
  // first evolve chi a half step
  
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;

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

        vScale[i][j] = eta[i][j];
        
      }
    }
  }


  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();
    
        // velocity half step
        
        info->matVecMul3( vScale, vel, sc );
    
        for (k = 0; k < 3; k++) 
          vel[k] += dt2 * ((frc[k]  / mass) * eConvert - sc[k]);

        myAtoms[j]->setVel( vel );

        // position whole step    
        
        info->matVecMul3( eta, rc, sc );

        for (k = 0; k < 3; k++ ) 
          pos[k] += dt * (vel[k] + sc[k]);
   
        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:
  
  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
  
  
  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;
      
    }
  }
  
  info->getBoxM(hm);
  info->matMul3(hm, scaleMat, hmnew);
  info->setBoxM(hmnew);
  
}

void NPTfm::moveB( void ){

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

  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double sc[3];
  double press[3][3], vScale[3][3];
  
  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  instaTemp = tStats->getTemperature();
  tStats->getPressureTensor(press);
  instaVol = tStats->getVolume();
   
  // first evolve chi a half step
  
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  
  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);

        vScale[i][j] = eta[i][j] + chi;
        
      } else {
        
        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);

        vScale[i][j] = eta[i][j];
        
      }
    }
  }

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

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

    mass = atoms[i]->getMass();
    
    // velocity half step
        
    info->matVecMul3( vScale, vel, sc );
    
    for (j = 0; j < 3; j++) {
      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);
      
      dAtom->setJ( ji );

    }                    
  }
}

int NPTfm::readyCheck() {
 
  // First check to see if we have a target temperature. 
  // Not having one is fatal. 
  
  if (!have_target_temp) {
    sprintf( painCave.errMsg,
             "NPTfm error: You can't use the NPTfm integrator\n"
             "   without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

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

  // We must set tauBarostat.
   
  if (!have_tau_barostat) {
    sprintf( painCave.errMsg,
             "NPTfm error: If you use the NPTfm\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:	606
Committed:	Tue Jul 15 03:28:05 2003 UTC (20 years, 11 months ago) by gezelter
File size:	8152 byte(s)
Log Message:	Added NPTfm
#	User	Rev	Content
1	gezelter	606	#include "Atom.hpp"
2			#include "Molecule.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			// Basic non-isotropic thermostating and barostating via the Melchionna
13			// modification of the Hoover algorithm:
14			//
15			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
16			// Molec. Phys., 78, 533.
17			//
18			// and
19			//
20			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
21
22			// The NPTfm variant scales the molecular center-of-mass coordinates
23			// instead of the atomic coordinates
24
25			NPTfm::NPTfm ( SimInfo theInfo, ForceFields the_ff):
26			Integrator( theInfo, the_ff )
27			{
28			int i, j;
29			chi = 0.0;
30
31			for(i = 0; i < 3; i++)
32			for (j = 0; j < 3; j++)
33			eta[i][j] = 0.0;
34
35			have_tau_thermostat = 0;
36			have_tau_barostat = 0;
37			have_target_temp = 0;
38			have_target_pressure = 0;
39			}
40
41			void NPTfm::moveA() {
42
43			int i, j, k;
44			DirectionalAtom* dAtom;
45			double Tb[3], ji[3];
46			double A[3][3], I[3][3];
47			double angle, mass;
48			double vel[3], pos[3], frc[3];
49
50			double rj[3];
51			double instaTemp, instaPress, instaVol;
52			double tt2, tb2;
53			double sc[3];
54			double eta2ij;
55			double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
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			tStats->getPressureTensor(press);
68			instaVol = tStats->getVolume();
69
70			// first evolve chi a half step
71
72			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
73
74			for (i = 0; i < 3; i++ ) {
75			for (j = 0; j < 3; j++ ) {
76			if (i == j) {
77
78			eta[i][j] += dt2 * instaVol *
79			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
80
81			vScale[i][j] = eta[i][j] + chi;
82
83			} else {
84
85			eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
86
87			vScale[i][j] = eta[i][j];
88
89			}
90			}
91			}
92
93
94			for (i = 0; i < nMols; i++) {
95
96			myMolecules[i].getCOM(rc);
97
98			nInMol = myMolecules[i].getNAtoms();
99			myAtoms = myMolecules[i].getMyAtoms();
100
101			// find the minimum image coordinates of the molecular centers of mass:
102
103			info->wrapVector(rc);
104
105			for( j=0; j< nInMol; j++ ){
106
107			if(myAtoms[j] != NULL) {
108
109			myAtoms[j]->getVel( vel );
110			myAtoms[j]->getPos( pos );
111			myAtoms[j]->getFrc( frc );
112
113			mass = myAtoms[j]->getMass();
114
115			// velocity half step
116
117			info->matVecMul3( vScale, vel, sc );
118
119			for (k = 0; k < 3; k++)
120			vel[k] += dt2 * ((frc[k] / mass) * eConvert - sc[k]);
121
122			myAtoms[j]->setVel( vel );
123
124			// position whole step
125
126			info->matVecMul3( eta, rc, sc );
127
128			for (k = 0; k < 3; k++ )
129			pos[k] += dt * (vel[k] + sc[k]);
130
131			if( myAtoms[j]->isDirectional() ){
132
133			dAtom = (DirectionalAtom *)myAtoms[j];
134
135			// get and convert the torque to body frame
136
137			dAtom->getTrq( Tb );
138			dAtom->lab2Body( Tb );
139
140			// get the angular momentum, and propagate a half step
141
142			dAtom->getJ( ji );
143
144			for (k=0; k < 3; k++)
145			ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
146
147			// use the angular velocities to propagate the rotation matrix a
148			// full time step
149
150			dAtom->getA(A);
151			dAtom->getI(I);
152
153			// rotate about the x-axis
154			angle = dt2 * ji[0] / I[0][0];
155			this->rotate( 1, 2, angle, ji, A );
156
157			// rotate about the y-axis
158			angle = dt2 * ji[1] / I[1][1];
159			this->rotate( 2, 0, angle, ji, A );
160
161			// rotate about the z-axis
162			angle = dt * ji[2] / I[2][2];
163			this->rotate( 0, 1, angle, ji, A);
164
165			// rotate about the y-axis
166			angle = dt2 * ji[1] / I[1][1];
167			this->rotate( 2, 0, angle, ji, A );
168
169			// rotate about the x-axis
170			angle = dt2 * ji[0] / I[0][0];
171			this->rotate( 1, 2, angle, ji, A );
172
173			dAtom->setJ( ji );
174			dAtom->setA( A );
175			}
176			}
177			}
178			}
179
180			// Scale the box after all the positions have been moved:
181
182			// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
183			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
184
185
186			for(i=0; i<3; i++){
187			for(j=0; j<3; j++){
188
189			// Calculate the matrix Product of the eta array (we only need
190			// the ij element right now):
191
192			eta2ij = 0.0;
193			for(k=0; k<3; k++){
194			eta2ij += eta[i][k] * eta[k][j];
195			}
196
197			scaleMat[i][j] = 0.0;
198			// identity matrix (see above):
199			if (i == j) scaleMat[i][j] = 1.0;
200			// Taylor expansion for the exponential truncated at second order:
201			scaleMat[i][j] += dteta[i][j] + 0.5dtdteta2ij;
202
203			}
204			}
205
206			info->getBoxM(hm);
207			info->matMul3(hm, scaleMat, hmnew);
208			info->setBoxM(hmnew);
209
210			}
211
212			void NPTfm::moveB( void ){
213
214			int i, j;
215			DirectionalAtom* dAtom;
216			double Tb[3], ji[3];
217			double vel[3], frc[3];
218			double mass;
219
220			double instaTemp, instaPress, instaVol;
221			double tt2, tb2;
222			double sc[3];
223			double press[3][3], vScale[3][3];
224
225			tt2 = tauThermostat * tauThermostat;
226			tb2 = tauBarostat * tauBarostat;
227
228			instaTemp = tStats->getTemperature();
229			tStats->getPressureTensor(press);
230			instaVol = tStats->getVolume();
231
232			// first evolve chi a half step
233
234			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
235
236			for (i = 0; i < 3; i++ ) {
237			for (j = 0; j < 3; j++ ) {
238			if (i == j) {
239
240			eta[i][j] += dt2 * instaVol *
241			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
242
243			vScale[i][j] = eta[i][j] + chi;
244
245			} else {
246
247			eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
248
249			vScale[i][j] = eta[i][j];
250
251			}
252			}
253			}
254
255			for( i=0; i<nAtoms; i++ ){
256
257			atoms[i]->getVel( vel );
258			atoms[i]->getFrc( frc );
259
260			mass = atoms[i]->getMass();
261
262			// velocity half step
263
264			info->matVecMul3( vScale, vel, sc );
265
266			for (j = 0; j < 3; j++) {
267			vel[j] += dt2 * ((frc[j] / mass) * eConvert - sc[j]);
268			}
269
270			atoms[i]->setVel( vel );
271
272			if( atoms[i]->isDirectional() ){
273
274			dAtom = (DirectionalAtom *)atoms[i];
275
276			// get and convert the torque to body frame
277
278			dAtom->getTrq( Tb );
279			dAtom->lab2Body( Tb );
280
281			// get the angular momentum, and propagate a half step
282
283			dAtom->getJ( ji );
284
285			for (j=0; j < 3; j++)
286			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
287
288			dAtom->setJ( ji );
289
290			}
291			}
292			}
293
294			int NPTfm::readyCheck() {
295
296			// First check to see if we have a target temperature.
297			// Not having one is fatal.
298
299			if (!have_target_temp) {
300			sprintf( painCave.errMsg,
301			"NPTfm error: You can't use the NPTfm integrator\n"
302			" without a targetTemp!\n"
303			);
304			painCave.isFatal = 1;
305			simError();
306			return -1;
307			}
308
309			if (!have_target_pressure) {
310			sprintf( painCave.errMsg,
311			"NPTfm error: You can't use the NPTfm integrator\n"
312			" without a targetPressure!\n"
313			);
314			painCave.isFatal = 1;
315			simError();
316			return -1;
317			}
318
319			// We must set tauThermostat.
320
321			if (!have_tau_thermostat) {
322			sprintf( painCave.errMsg,
323			"NPTfm error: If you use the NPTfm\n"
324			" integrator, you must set tauThermostat.\n");
325			painCave.isFatal = 1;
326			simError();
327			return -1;
328			}
329
330			// We must set tauBarostat.
331
332			if (!have_tau_barostat) {
333			sprintf( painCave.errMsg,
334			"NPTfm error: If you use the NPTfm\n"
335			" integrator, you must set tauBarostat.\n");
336			painCave.isFatal = 1;
337			simError();
338			return -1;
339			}
340
341			// We need NkBT a lot, so just set it here:
342
343			NkBT = (double)info->ndf * kB * targetTemp;
344
345			return 1;
346			}