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

  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]->getPos( pos );
    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]);
      rj[j] = pos[j];
    }

    atoms[i]->setVel( vel );

    // position whole step    

    info->wrapVector(rj);

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

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