OOPSE/libmdtools/NPTim.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.

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

NPTim::NPTim ( 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 NPTim::moveA() {
  
  int i,j,k;
  int nInMol, aMatIndex;
  DirectionalAtom* dAtom;
  Atom** theAtoms;
  Molecule** myMols;
  double Tb[3];
  double ji[3];
  double rc[3];
  double mass;
  double rx, ry, rz, vx, vy, vz, fx, fy, fz;
  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
  double angle;

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

  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 < nMols; i++) {

    myMols[i].getCOM(rc);

    nInMol = myMols[i]->getNAtoms();
    theAtoms = myMols[i]->getMyAtoms();

    // find the minimum image coordinates of the molecular centers of mass:

    info->wrapVector(rc);

    for (j = 0; j < nInMol; j++) {

      if(theAtoms[j] != NULL) {

        aMatIndex = 9 * theAtoms[j]->getIndex();
        
        mass = theAtoms[j]->getMass();
        
        vx = theAtoms[j]->get_vx();
        vy = theAtoms[j]->get_vy();
        vz = theAtoms[j]->get_vz();
        
        fx = theAtoms[j]->getFx();
        fy = theAtoms[j]->getFy();
        fz = theAtoms[j]->getFz();

        rx = theAtoms[j]->getX();
        ry = theAtoms[j]->getY();
        rz = theAtoms[j]->getZ();

        // velocity half step

        vx += dt2 * ((fx / mass)*eConvert - vx*(chi+eta));
        vy += dt2 * ((fy / mass)*eConvert - vy*(chi+eta));
        vz += dt2 * ((fz / mass)*eConvert - vz*(chi+eta));

        // position whole step    

        rx += dt*(vx + eta*rc[0]);
        ry += dt*(vy + eta*rc[1]);
        rz += dt*(vz + eta*rc[2]);
        
        theAtoms[j]->set_vx(vx);        
        theAtoms[j]->set_vy(vy);
        theAtoms[j]->set_vz(vz);

        theAtoms[j]->setX(rx);        
        theAtoms[j]->setY(ry);
        theAtoms[j]->setZ(rz);

        if( theAtoms[j]->isDirectional() ){

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