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 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;
  chi = 0.0;
  for(i = 0; i < 9; i++) eta[i] = 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 rj[3];
  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double angle;
  double press[9];
  const double p_convert = 1.63882576e8;

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

  instaTemp = tStats->getTemperature();
  tStats->getPressureTensor(press);

  for (i=0; i < 9; i++) press[i] *= p_convert;

  instaVol = tStats->getVolume();
   
  // first evolve chi a half step
  
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  
  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
  eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
  eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
  
  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;
    aMatIndex = i * 9;
    
    // velocity half step
    
    vx = vel[atomIndex];
    vy = vel[atomIndex+1];
    vz = vel[atomIndex+2];
    
    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
    
    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);

    vel[atomIndex] = vx;
    vel[atomIndex+1] = vy;
    vel[atomIndex+2] = vz;

    // position whole step    

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

    info->wrapVector(rj);

    scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2];
    scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2];
    scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2];

    pos[atomIndex] += dt * (vel[atomIndex] + scx);
    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
   
    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
   
  info->getBoxM(hm);
  

   info->scaleBox(exp(dt*eta));


}

void NPTi::moveB( void ){
  int i,j,k;
  int atomIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];
  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  
  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  instaVol = tStats->getVolume();

  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
  
  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;
    
    // velocity half step
    for( j=atomIndex; j<(atomIndex+3); j++ )
    for( j=atomIndex; j<(atomIndex+3); j++ )
      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert 
                       - vel[j]*(chi+eta));
    
    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 NPTi::readyCheck() {
 
  // First check to see if we have a target temperature. 
  // Not having one is fatal. 
  
  if (!have_target_temp) {
    sprintf( painCave.errMsg,
             "NPTi error: You can't use the NPTi integrator\n"
             "   without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

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

  // We must set tauBarostat.
   
  if (!have_tau_barostat) {
    sprintf( painCave.errMsg,
             "NPTi error: If you use the NPTi\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:	577
Committed:	Wed Jul 9 01:41:11 2003 UTC (21 years ago) by gezelter
File size:	7508 byte(s)
Log Message:	Fixes in NPTi migrated into NPTf
#	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			// Basic 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	gezelter	577	NPTf::NPTf ( SimInfo theInfo, ForceFields the_ff):
23	gezelter	576	Integrator( theInfo, the_ff )
24			{
25			int i;
26			chi = 0.0;
27			for(i = 0; i < 9; i++) eta[i] = 0.0;
28			have_tau_thermostat = 0;
29			have_tau_barostat = 0;
30			have_target_temp = 0;
31			have_target_pressure = 0;
32			}
33
34	gezelter	577	void NPTf::moveA() {
35	gezelter	576
36			int i,j,k;
37			int atomIndex, aMatIndex;
38			DirectionalAtom* dAtom;
39			double Tb[3];
40			double ji[3];
41			double rj[3];
42			double instaTemp, instaPress, instaVol;
43			double tt2, tb2;
44			double angle;
45	gezelter	577	double press[9];
46			const double p_convert = 1.63882576e8;
47	gezelter	576
48			tt2 = tauThermostat * tauThermostat;
49			tb2 = tauBarostat * tauBarostat;
50
51			instaTemp = tStats->getTemperature();
52	gezelter	577	tStats->getPressureTensor(press);
53
54			for (i=0; i < 9; i++) press[i] *= p_convert;
55
56	gezelter	576	instaVol = tStats->getVolume();
57
58			// first evolve chi a half step
59
60			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
61
62	gezelter	577	eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
63			eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
64			eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
65			eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
66			eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
67			eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
68			eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
69			eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
70			eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
71
72	gezelter	576	for( i=0; i<nAtoms; i++ ){
73			atomIndex = i * 3;
74			aMatIndex = i * 9;
75
76			// velocity half step
77	gezelter	577
78			vx = vel[atomIndex];
79			vy = vel[atomIndex+1];
80			vz = vel[atomIndex+2];
81
82			scx = (chi + eta[0])vx + eta[1]vy + eta[2]*vz;
83			scy = eta[3]vx + (chi + eta[4])vy + eta[5]*vz;
84			scz = eta[6]vx + eta[7]vy + (chi + eta[8])*vz;
85
86			vx += dt2 * ((frc[atomIndex] /atoms[i]->getMass())*eConvert - scx);
87			vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
88			vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
89	gezelter	576
90	gezelter	577	vel[atomIndex] = vx;
91			vel[atomIndex+1] = vy;
92			vel[atomIndex+2] = vz;
93
94	gezelter	576	// position whole step
95
96	gezelter	577	rj[0] = pos[atomIndex];
97			rj[1] = pos[atomIndex+1];
98			rj[2] = pos[atomIndex+2];
99	gezelter	576
100	gezelter	577	info->wrapVector(rj);
101	gezelter	576
102	gezelter	577	scx = eta[0]rj[0] + eta[1]rj[1] + eta[2]*rj[2];
103			scy = eta[3]rj[0] + eta[4]rj[1] + eta[5]*rj[2];
104			scz = eta[6]rj[0] + eta[7]rj[1] + eta[8]*rj[2];
105	gezelter	576
106	gezelter	577	pos[atomIndex] += dt * (vel[atomIndex] + scx);
107			pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
108			pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
109	gezelter	576
110			if( atoms[i]->isDirectional() ){
111
112			dAtom = (DirectionalAtom *)atoms[i];
113
114			// get and convert the torque to body frame
115
116			Tb[0] = dAtom->getTx();
117			Tb[1] = dAtom->getTy();
118			Tb[2] = dAtom->getTz();
119
120			dAtom->lab2Body( Tb );
121
122			// get the angular momentum, and propagate a half step
123
124			ji[0] = dAtom->getJx();
125			ji[1] = dAtom->getJy();
126			ji[2] = dAtom->getJz();
127
128			ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
129			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
130			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
131
132			// use the angular velocities to propagate the rotation matrix a
133			// full time step
134
135			// rotate about the x-axis
136			angle = dt2 * ji[0] / dAtom->getIxx();
137			this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
138
139			// rotate about the y-axis
140			angle = dt2 * ji[1] / dAtom->getIyy();
141			this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
142
143			// rotate about the z-axis
144			angle = dt * ji[2] / dAtom->getIzz();
145			this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
146
147			// rotate about the y-axis
148			angle = dt2 * ji[1] / dAtom->getIyy();
149			this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
150
151			// rotate about the x-axis
152			angle = dt2 * ji[0] / dAtom->getIxx();
153			this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
154
155			dAtom->setJx( ji[0] );
156			dAtom->setJy( ji[1] );
157			dAtom->setJz( ji[2] );
158			}
159
160			}
161	gezelter	577
162			// Scale the box after all the positions have been moved:
163
164
165
166			// Use a taylor expansion for eta products
167
168			info->getBoxM(hm);
169
170
171
172
173
174
175			info->scaleBox(exp(dt*eta));
176
177
178	gezelter	576	}
179
180			void NPTi::moveB( void ){
181			int i,j,k;
182			int atomIndex;
183			DirectionalAtom* dAtom;
184			double Tb[3];
185			double ji[3];
186			double instaTemp, instaPress, instaVol;
187			double tt2, tb2;
188
189			tt2 = tauThermostat * tauThermostat;
190			tb2 = tauBarostat * tauBarostat;
191
192			instaTemp = tStats->getTemperature();
193			instaPress = tStats->getPressure();
194			instaVol = tStats->getVolume();
195
196			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
197			eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2));
198
199			for( i=0; i<nAtoms; i++ ){
200			atomIndex = i * 3;
201
202			// velocity half step
203			for( j=atomIndex; j<(atomIndex+3); j++ )
204			for( j=atomIndex; j<(atomIndex+3); j++ )
205			vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
206			- vel[j]*(chi+eta));
207
208			if( atoms[i]->isDirectional() ){
209
210			dAtom = (DirectionalAtom *)atoms[i];
211
212			// get and convert the torque to body frame
213
214			Tb[0] = dAtom->getTx();
215			Tb[1] = dAtom->getTy();
216			Tb[2] = dAtom->getTz();
217
218			dAtom->lab2Body( Tb );
219
220			// get the angular momentum, and complete the angular momentum
221			// half step
222
223			ji[0] = dAtom->getJx();
224			ji[1] = dAtom->getJy();
225			ji[2] = dAtom->getJz();
226
227			ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
228			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
229			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
230
231			dAtom->setJx( ji[0] );
232			dAtom->setJy( ji[1] );
233			dAtom->setJz( ji[2] );
234			}
235			}
236			}
237
238			int NPTi::readyCheck() {
239
240			// First check to see if we have a target temperature.
241			// Not having one is fatal.
242
243			if (!have_target_temp) {
244			sprintf( painCave.errMsg,
245			"NPTi error: You can't use the NPTi integrator\n"
246			" without a targetTemp!\n"
247			);
248			painCave.isFatal = 1;
249			simError();
250			return -1;
251			}
252
253			if (!have_target_pressure) {
254			sprintf( painCave.errMsg,
255			"NPTi error: You can't use the NPTi integrator\n"
256			" without a targetPressure!\n"
257			);
258			painCave.isFatal = 1;
259			simError();
260			return -1;
261			}
262
263			// We must set tauThermostat.
264
265			if (!have_tau_thermostat) {
266			sprintf( painCave.errMsg,
267			"NPTi error: If you use the NPTi\n"
268			" integrator, you must set tauThermostat.\n");
269			painCave.isFatal = 1;
270			simError();
271			return -1;
272			}
273
274			// We must set tauBarostat.
275
276			if (!have_tau_barostat) {
277			sprintf( painCave.errMsg,
278			"NPTi error: If you use the NPTi\n"
279			" integrator, you must set tauBarostat.\n");
280			painCave.isFatal = 1;
281			simError();
282			return -1;
283			}
284
285			// We need NkBT a lot, so just set it here:
286
287			NkBT = (double)info->ndf * kB * targetTemp;
288
289			return 1;
290			}