OOPSE/libmdtools/NPTi.cpp

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

NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
  Integrator( theInfo, the_ff )
{
  chi = 0.0;
  eta = 0.0;
  have_tau_thermostat = 0;
  have_tau_barostat = 0;
  have_target_temp = 0;
  have_target_pressure = 0;
}

void NPTi::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;


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

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  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));

  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;
    aMatIndex = i * 9;
    
    // velocity half step
    for( j=atomIndex; j<(atomIndex+3); j++ )
      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert 
                       - vel[j]*(chi+eta));

    // position whole step    

    rj[0] = pos[atomIndex];
    rj[1] = pos[atomIndex+1];
    rj[2] = pos[atomIndex+2];
    
    info->wrapVector(rj);

    pos[atomIndex] += dt * (vel[atomIndex] + eta*rj[0]);
    pos[atomIndex+1] += dt * (vel[atomIndex+1] + eta*rj[1]);
    pos[atomIndex+2] += dt * (vel[atomIndex+2] + eta*rj[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:

  cerr << "eta = " << eta 
       << "; exp(dt*eta) = " << exp(eta*dt) << "\n";

  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) / 
                 (p_convert*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:	586
Committed:	Wed Jul 9 22:14:06 2003 UTC (21 years ago) by mmeineke
File size:	6410 byte(s)
Log Message:	Bug fixing NPTi and NPTf. there is some error in the caclulation of HmatInverse.
#	User	Rev	Content
1	gezelter	578	#include <cmath>
2	gezelter	574	#include "Atom.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
13			// Basic isotropic thermostating and barostating via the Melchionna
14			// modification of the Hoover algorithm:
15			//
16			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
17			// Molec. Phys., 78, 533.
18			//
19			// and
20			//
21			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
22
23			NPTi::NPTi ( SimInfo theInfo, ForceFields the_ff):
24			Integrator( theInfo, the_ff )
25			{
26			chi = 0.0;
27			eta = 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			void NPTi::moveA() {
35
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
46	mmeineke	586
47	gezelter	574	tt2 = tauThermostat * tauThermostat;
48			tb2 = tauBarostat * tauBarostat;
49
50			instaTemp = tStats->getTemperature();
51			instaPress = tStats->getPressure();
52			instaVol = tStats->getVolume();
53
54	mmeineke	586	// first evolve chi a half step
55	gezelter	574
56			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
57	mmeineke	586	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
58			(p_convertNkBTtb2));
59	gezelter	574
60			for( i=0; i<nAtoms; i++ ){
61			atomIndex = i * 3;
62			aMatIndex = i * 9;
63
64			// velocity half step
65			for( j=atomIndex; j<(atomIndex+3); j++ )
66			vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
67			- vel[j]*(chi+eta));
68
69			// position whole step
70
71	gezelter	577	rj[0] = pos[atomIndex];
72			rj[1] = pos[atomIndex+1];
73			rj[2] = pos[atomIndex+2];
74
75			info->wrapVector(rj);
76	gezelter	574
77	gezelter	577	pos[atomIndex] += dt * (vel[atomIndex] + eta*rj[0]);
78			pos[atomIndex+1] += dt * (vel[atomIndex+1] + eta*rj[1]);
79			pos[atomIndex+2] += dt * (vel[atomIndex+2] + eta*rj[2]);
80	gezelter	574
81			if( atoms[i]->isDirectional() ){
82
83			dAtom = (DirectionalAtom *)atoms[i];
84
85			// get and convert the torque to body frame
86
87			Tb[0] = dAtom->getTx();
88			Tb[1] = dAtom->getTy();
89			Tb[2] = dAtom->getTz();
90
91			dAtom->lab2Body( Tb );
92
93			// get the angular momentum, and propagate a half step
94
95			ji[0] = dAtom->getJx();
96			ji[1] = dAtom->getJy();
97			ji[2] = dAtom->getJz();
98
99			ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
100			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
101			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
102
103			// use the angular velocities to propagate the rotation matrix a
104			// full time step
105
106			// rotate about the x-axis
107			angle = dt2 * ji[0] / dAtom->getIxx();
108			this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
109
110			// rotate about the y-axis
111			angle = dt2 * ji[1] / dAtom->getIyy();
112			this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
113
114			// rotate about the z-axis
115			angle = dt * ji[2] / dAtom->getIzz();
116			this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
117
118			// rotate about the y-axis
119			angle = dt2 * ji[1] / dAtom->getIyy();
120			this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
121
122			// rotate about the x-axis
123			angle = dt2 * ji[0] / dAtom->getIxx();
124			this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
125
126			dAtom->setJx( ji[0] );
127			dAtom->setJy( ji[1] );
128			dAtom->setJz( ji[2] );
129			}
130
131			}
132	gezelter	577	// Scale the box after all the positions have been moved:
133
134	mmeineke	586	cerr << "eta = " << eta
135			<< "; exp(dteta) = " << exp(etadt) << "\n";
136
137	gezelter	577	info->scaleBox(exp(dt*eta));
138
139	gezelter	574	}
140
141			void NPTi::moveB( void ){
142			int i,j,k;
143			int atomIndex;
144			DirectionalAtom* dAtom;
145			double Tb[3];
146			double ji[3];
147			double instaTemp, instaPress, instaVol;
148			double tt2, tb2;
149
150			tt2 = tauThermostat * tauThermostat;
151			tb2 = tauBarostat * tauBarostat;
152
153			instaTemp = tStats->getTemperature();
154			instaPress = tStats->getPressure();
155			instaVol = tStats->getVolume();
156
157			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
158	mmeineke	586	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
159			(p_convertNkBTtb2));
160	gezelter	574
161			for( i=0; i<nAtoms; i++ ){
162			atomIndex = i * 3;
163
164			// velocity half step
165			for( j=atomIndex; j<(atomIndex+3); j++ )
166			for( j=atomIndex; j<(atomIndex+3); j++ )
167			vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert
168			- vel[j]*(chi+eta));
169
170			if( atoms[i]->isDirectional() ){
171
172			dAtom = (DirectionalAtom *)atoms[i];
173
174			// get and convert the torque to body frame
175
176			Tb[0] = dAtom->getTx();
177			Tb[1] = dAtom->getTy();
178			Tb[2] = dAtom->getTz();
179
180			dAtom->lab2Body( Tb );
181
182			// get the angular momentum, and complete the angular momentum
183			// half step
184
185			ji[0] = dAtom->getJx();
186			ji[1] = dAtom->getJy();
187			ji[2] = dAtom->getJz();
188
189			ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
190			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
191			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
192
193			dAtom->setJx( ji[0] );
194			dAtom->setJy( ji[1] );
195			dAtom->setJz( ji[2] );
196			}
197			}
198			}
199
200			int NPTi::readyCheck() {
201
202			// First check to see if we have a target temperature.
203			// Not having one is fatal.
204
205			if (!have_target_temp) {
206			sprintf( painCave.errMsg,
207			"NPTi error: You can't use the NPTi integrator\n"
208			" without a targetTemp!\n"
209			);
210			painCave.isFatal = 1;
211			simError();
212			return -1;
213			}
214
215			if (!have_target_pressure) {
216			sprintf( painCave.errMsg,
217			"NPTi error: You can't use the NPTi integrator\n"
218			" without a targetPressure!\n"
219			);
220			painCave.isFatal = 1;
221			simError();
222			return -1;
223			}
224
225			// We must set tauThermostat.
226
227			if (!have_tau_thermostat) {
228			sprintf( painCave.errMsg,
229			"NPTi error: If you use the NPTi\n"
230			" integrator, you must set tauThermostat.\n");
231			painCave.isFatal = 1;
232			simError();
233			return -1;
234			}
235
236			// We must set tauBarostat.
237
238			if (!have_tau_barostat) {
239			sprintf( painCave.errMsg,
240			"NPTi error: If you use the NPTi\n"
241			" integrator, you must set tauBarostat.\n");
242			painCave.isFatal = 1;
243			simError();
244			return -1;
245			}
246
247			// We need NkBT a lot, so just set it here:
248
249			NkBT = (double)info->ndf * kB * targetTemp;
250
251			return 1;
252			}