OOPSE/libmdtools/NPTf.cpp

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

NPTf::NPTf ( 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 NPTf::moveA() {
  
  int i,j,k;
  int atomIndex, aMatIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];
  double ri[3], vi[3], sc[3];
  double instaTemp, instaVol;
  double tt2, tb2, eta2ij;
  double angle;
  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[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++ ){
    atomIndex = i * 3;
    aMatIndex = i * 9;
    
    // velocity half step
    
    vi[0] = vel[atomIndex];
    vi[1] = vel[atomIndex+1];
    vi[2] = vel[atomIndex+2];
    
    info->matVecMul3( vScale, vi, sc );
    
    vi[0] += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - sc[0]);
    vi[1] += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - sc[1]);
    vi[2] += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - sc[2]);

    vel[atomIndex]   = vi[0];
    vel[atomIndex+1] = vi[1];
    vel[atomIndex+2] = vi[2];

    // position whole step    

    ri[0] = pos[atomIndex];
    ri[1] = pos[atomIndex+1];
    ri[2] = pos[atomIndex+2];

    info->wrapVector(ri);

    info->matVecMul3( eta, ri, sc );

    pos[atomIndex]   += dt * (vel[atomIndex] + sc[0]);
    pos[atomIndex+1] += dt * (vel[atomIndex+1] + sc[1]);
    pos[atomIndex+2] += dt * (vel[atomIndex+2] + sc[2]);
   
    if( atoms[i]->isDirectional() ){

      dAtom = (DirectionalAtom *)atoms[i];
          
      // get and convert the torque to body frame
      
      Tb[0] = dAtom->getTx();
      Tb[1] = dAtom->getTy();
      Tb[2] = dAtom->getTz();
      
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and propagate a half step

      ji[0] = dAtom->getJx();
      ji[1] = dAtom->getJy();
      ji[2] = dAtom->getJz();
      
      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
      
      // use the angular velocities to propagate the rotation matrix a
      // full time step
      
      // rotate about the x-axis      
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); 
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
      
      // rotate about the z-axis
      angle = dt * ji[2] / dAtom->getIzz();
      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
    
  }

  // 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 NPTf::moveB( void ){
  int i,j, k;
  int atomIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];
  double vi[3], sc[3];
  double instaTemp, instaVol;
  double tt2, tb2;
  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++ ){
    atomIndex = i * 3;

    // velocity half step
    
    vi[0] = vel[atomIndex];
    vi[1] = vel[atomIndex+1];
    vi[2] = vel[atomIndex+2];
    
    info->matVecMul3( vScale, vi, sc );
    
    vi[0] += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - sc[0]);
    vi[1] += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - sc[1]);
    vi[2] += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - sc[2]);

    vel[atomIndex]   = vi[0];
    vel[atomIndex+1] = vi[1];
    vel[atomIndex+2] = vi[2];
    
    if( atoms[i]->isDirectional() ){
      
      dAtom = (DirectionalAtom *)atoms[i];
      
      // get and convert the torque to body frame
      
      Tb[0] = dAtom->getTx();
      Tb[1] = dAtom->getTy();
      Tb[2] = dAtom->getTz();
      
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and complete the angular momentum
      // half step
      
      ji[0] = dAtom->getJx();
      ji[1] = dAtom->getJy();
      ji[2] = dAtom->getJz();
      
      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
  }
}

int NPTf::readyCheck() {
 
  // 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)info->ndf * kB * targetTemp;

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