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

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#endif


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

template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  chi = 0.0;
  integralOfChidt = 0.0;
  have_tau_thermostat = 0;
  have_tau_barostat = 0;
  have_target_temp = 0;
  have_target_pressure = 0;
  have_chi_tolerance = 0;
  have_eta_tolerance = 0;
  have_pos_iter_tolerance = 0;

  oldPos = new double[3*nAtoms];
  oldVel = new double[3*nAtoms];
  oldJi = new double[3*nAtoms];
#ifdef IS_MPI
  Nparticles = mpiSim->getTotAtoms();
#else
  Nparticles = theInfo->n_atoms;
#endif

}

template<typename T> NPT<T>::~NPT() {
  delete[] oldPos;
  delete[] oldVel;
  delete[] oldJi;
}

template<typename T> void NPT<T>::moveA() {
 
  //new version of NPT
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double mass;
  double vel[3], pos[3], frc[3];
  double sc[3];
  double COM[3];

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  instaVol = tStats->getVolume();
  
  tStats->getCOM(COM);
  
  //evolve velocity half step
  for( i=0; i<nAtoms; i++ ){

    atoms[i]->getVel( vel );
    atoms[i]->getFrc( frc );

    mass = atoms[i]->getMass();

    getVelScaleA( sc, vel );

    for (j=0; j < 3; j++) {
      
      // velocity half step  (use chi from previous step here):
      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);
      
      this->rotationPropagation( dAtom, ji );
      
      dAtom->setJ( ji );
    }    
  }

  // evolve chi and eta  half step
  
  evolveChiA();
  evolveEtaA();

  //calculate the integral of chidt
  integralOfChidt += dt2*chi;

  //save the old positions
  for(i = 0; i < nAtoms; i++){
    atoms[i]->getPos(pos);
    for(j = 0; j < 3; j++)
      oldPos[i*3 + j] = pos[j];
  }
  
  //the first estimation of r(t+dt) is equal to  r(t) 
    
  for(k = 0; k < 5; k ++){

    for(i =0 ; i < nAtoms; i++){

      atoms[i]->getVel(vel);
      atoms[i]->getPos(pos);

      this->getPosScale( pos, COM, i, sc );
   
      for(j = 0; j < 3; j++)
        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);

      atoms[i]->setPos( pos );
    }
    
    if (nConstrained){
      constrainA();
    }
  }
    

  // Scale the box after all the positions have been moved:
  
  this->scaleSimBox();
}

template<typename T> void NPT<T>::moveB( void ){
  
  //new version of NPT
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3], sc[3];
  double vel[3], frc[3];
  double mass;
  
  // Set things up for the iteration:

  for( i=0; i<nAtoms; i++ ){

    atoms[i]->getVel( vel );

    for (j=0; j < 3; j++)
      oldVel[3*i + j]  = vel[j];

    if( atoms[i]->isDirectional() ){

      dAtom = (DirectionalAtom *)atoms[i];

      dAtom->getJ( ji );

      for (j=0; j < 3; j++)
        oldJi[3*i + j] = ji[j];

    }
  }

  // do the iteration:

  instaVol = tStats->getVolume();
  
  for (k=0; k < 4; k++) {
    
    instaTemp = tStats->getTemperature();
    instaPress = tStats->getPressure();

    // evolve chi another half step using the temperature at t + dt/2

    this->evolveChiB();
    this->evolveEtaB();
    
    for( i=0; i<nAtoms; i++ ){

      atoms[i]->getFrc( frc );
      atoms[i]->getVel(vel);
      
      mass = atoms[i]->getMass();
      
      getVelScaleB( sc, i );

      // velocity half step
      for (j=0; j < 3; j++) 
        vel[j] = oldVel[3*i+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 );       
             
        for (j=0; j < 3; j++) 
          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
      
          dAtom->setJ( ji );
      }
    }
    
    if (nConstrained){
      constrainB();
    }    
    
    if ( this->chiConverged() && this->etaConverged() ) break;
  }

  //calculate integral of chida
  integralOfChidt += dt2*chi;


}

template<typename T> void NPT<T>::resetIntegrator() {
  chi = 0.0;
  T::resetIntegrator();
}

template<typename T> void NPT<T>::evolveChiA() {
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  oldChi = chi;
}

template<typename T> void NPT<T>::evolveChiB() {
  
  prevChi = chi;
  chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
}

template<typename T> bool NPT<T>::chiConverged() {
  
  return ( fabs( prevChi - chi ) <= chiTolerance );  
}

template<typename T> int NPT<T>::readyCheck() {

  //check parent's readyCheck() first
  if (T::readyCheck() == -1)
    return -1;
 
  // First check to see if we have a target temperature. 
  // Not having one is fatal. 
  
  if (!have_target_temp) {
    sprintf( painCave.errMsg,
             "NPT error: You can't use the NPT integrator\n"
             "   without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

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

  // We must set tauBarostat.
   
  if (!have_tau_barostat) {
    sprintf( painCave.errMsg,
             "NPT error: If you use the NPT\n"
             "   integrator, you must set tauBarostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }    

  if (!have_chi_tolerance) {
    sprintf( painCave.errMsg,
             "NPT warning: setting chi tolerance to 1e-6\n");
    chiTolerance = 1e-6;
    have_chi_tolerance = 1;
    painCave.isFatal = 0;
    simError();
  } 

  if (!have_eta_tolerance) {
    sprintf( painCave.errMsg,
             "NPT warning: setting eta tolerance to 1e-6\n");
    etaTolerance = 1e-6;
    have_eta_tolerance = 1;
    painCave.isFatal = 0;
    simError();
  } 
  
  // We need NkBT a lot, so just set it here: This is the RAW number
  // of particles, so no subtraction or addition of constraints or
  // orientational degrees of freedom:
  
  NkBT = (double)Nparticles * kB * targetTemp;
  
  // fkBT is used because the thermostat operates on more degrees of freedom
  // than the barostat (when there are particles with orientational degrees
  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
  
  fkBT = (double)info->ndf * kB * targetTemp;

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

  return 1;
}
Revision:	778
Committed:	Fri Sep 19 20:00:27 2003 UTC (20 years, 9 months ago) by mmeineke
File size:	7719 byte(s)
Log Message:	added NPT base class. NPTi is up to date. NPTf is not.
#	User	Rev	Content
1	mmeineke	778	#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			#ifdef IS_MPI
13			#include "mpiSimulation.hpp"
14			#endif
15
16
17			// Basic isotropic thermostating and barostating via the Melchionna
18			// modification of the Hoover algorithm:
19			//
20			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
21			// Molec. Phys., 78, 533.
22			//
23			// and
24			//
25			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
26
27			template<typename T> NPT<T>::NPT ( SimInfo theInfo, ForceFields the_ff):
28			T( theInfo, the_ff )
29			{
30			chi = 0.0;
31			integralOfChidt = 0.0;
32			have_tau_thermostat = 0;
33			have_tau_barostat = 0;
34			have_target_temp = 0;
35			have_target_pressure = 0;
36			have_chi_tolerance = 0;
37			have_eta_tolerance = 0;
38			have_pos_iter_tolerance = 0;
39
40			oldPos = new double[3*nAtoms];
41			oldVel = new double[3*nAtoms];
42			oldJi = new double[3*nAtoms];
43			#ifdef IS_MPI
44			Nparticles = mpiSim->getTotAtoms();
45			#else
46			Nparticles = theInfo->n_atoms;
47			#endif
48
49			}
50
51			template<typename T> NPT<T>::~NPT() {
52			delete[] oldPos;
53			delete[] oldVel;
54			delete[] oldJi;
55			}
56
57			template<typename T> void NPT<T>::moveA() {
58
59			//new version of NPT
60			int i, j, k;
61			DirectionalAtom* dAtom;
62			double Tb[3], ji[3];
63			double mass;
64			double vel[3], pos[3], frc[3];
65			double sc[3];
66			double COM[3];
67
68			instaTemp = tStats->getTemperature();
69			instaPress = tStats->getPressure();
70			instaVol = tStats->getVolume();
71
72			tStats->getCOM(COM);
73
74			//evolve velocity half step
75			for( i=0; i<nAtoms; i++ ){
76
77			atoms[i]->getVel( vel );
78			atoms[i]->getFrc( frc );
79
80			mass = atoms[i]->getMass();
81
82			getVelScaleA( sc, vel );
83
84			for (j=0; j < 3; j++) {
85
86			// velocity half step (use chi from previous step here):
87			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
88
89			}
90
91			atoms[i]->setVel( vel );
92
93			if( atoms[i]->isDirectional() ){
94
95			dAtom = (DirectionalAtom *)atoms[i];
96
97			// get and convert the torque to body frame
98
99			dAtom->getTrq( Tb );
100			dAtom->lab2Body( Tb );
101
102			// get the angular momentum, and propagate a half step
103
104			dAtom->getJ( ji );
105
106			for (j=0; j < 3; j++)
107			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
108
109			this->rotationPropagation( dAtom, ji );
110
111			dAtom->setJ( ji );
112			}
113			}
114
115			// evolve chi and eta half step
116
117			evolveChiA();
118			evolveEtaA();
119
120			//calculate the integral of chidt
121			integralOfChidt += dt2*chi;
122
123			//save the old positions
124			for(i = 0; i < nAtoms; i++){
125			atoms[i]->getPos(pos);
126			for(j = 0; j < 3; j++)
127			oldPos[i*3 + j] = pos[j];
128			}
129
130			//the first estimation of r(t+dt) is equal to r(t)
131
132			for(k = 0; k < 5; k ++){
133
134			for(i =0 ; i < nAtoms; i++){
135
136			atoms[i]->getVel(vel);
137			atoms[i]->getPos(pos);
138
139			this->getPosScale( pos, COM, i, sc );
140
141			for(j = 0; j < 3; j++)
142			pos[j] = oldPos[i3 + j] + dt(vel[j] + sc[j]);
143
144			atoms[i]->setPos( pos );
145			}
146
147			if (nConstrained){
148			constrainA();
149			}
150			}
151
152
153			// Scale the box after all the positions have been moved:
154
155			this->scaleSimBox();
156			}
157
158			template<typename T> void NPT<T>::moveB( void ){
159
160			//new version of NPT
161			int i, j, k;
162			DirectionalAtom* dAtom;
163			double Tb[3], ji[3], sc[3];
164			double vel[3], frc[3];
165			double mass;
166
167			// Set things up for the iteration:
168
169			for( i=0; i<nAtoms; i++ ){
170
171			atoms[i]->getVel( vel );
172
173			for (j=0; j < 3; j++)
174			oldVel[3*i + j] = vel[j];
175
176			if( atoms[i]->isDirectional() ){
177
178			dAtom = (DirectionalAtom *)atoms[i];
179
180			dAtom->getJ( ji );
181
182			for (j=0; j < 3; j++)
183			oldJi[3*i + j] = ji[j];
184
185			}
186			}
187
188			// do the iteration:
189
190			instaVol = tStats->getVolume();
191
192			for (k=0; k < 4; k++) {
193
194			instaTemp = tStats->getTemperature();
195			instaPress = tStats->getPressure();
196
197			// evolve chi another half step using the temperature at t + dt/2
198
199			this->evolveChiB();
200			this->evolveEtaB();
201
202			for( i=0; i<nAtoms; i++ ){
203
204			atoms[i]->getFrc( frc );
205			atoms[i]->getVel(vel);
206
207			mass = atoms[i]->getMass();
208
209			getVelScaleB( sc, i );
210
211			// velocity half step
212			for (j=0; j < 3; j++)
213			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - sc[j]);
214
215			atoms[i]->setVel( vel );
216
217			if( atoms[i]->isDirectional() ){
218
219			dAtom = (DirectionalAtom *)atoms[i];
220
221			// get and convert the torque to body frame
222
223			dAtom->getTrq( Tb );
224			dAtom->lab2Body( Tb );
225
226			for (j=0; j < 3; j++)
227			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
228
229			dAtom->setJ( ji );
230			}
231			}
232
233			if (nConstrained){
234			constrainB();
235			}
236
237			if ( this->chiConverged() && this->etaConverged() ) break;
238			}
239
240			//calculate integral of chida
241			integralOfChidt += dt2*chi;
242
243
244			}
245
246			template<typename T> void NPT<T>::resetIntegrator() {
247			chi = 0.0;
248			T::resetIntegrator();
249			}
250
251			template<typename T> void NPT<T>::evolveChiA() {
252			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
253			oldChi = chi;
254			}
255
256			template<typename T> void NPT<T>::evolveChiB() {
257
258			prevChi = chi;
259			chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
260			}
261
262			template<typename T> bool NPT<T>::chiConverged() {
263
264			return ( fabs( prevChi - chi ) <= chiTolerance );
265			}
266
267			template<typename T> int NPT<T>::readyCheck() {
268
269			//check parent's readyCheck() first
270			if (T::readyCheck() == -1)
271			return -1;
272
273			// First check to see if we have a target temperature.
274			// Not having one is fatal.
275
276			if (!have_target_temp) {
277			sprintf( painCave.errMsg,
278			"NPT error: You can't use the NPT integrator\n"
279			" without a targetTemp!\n"
280			);
281			painCave.isFatal = 1;
282			simError();
283			return -1;
284			}
285
286			if (!have_target_pressure) {
287			sprintf( painCave.errMsg,
288			"NPT error: You can't use the NPT integrator\n"
289			" without a targetPressure!\n"
290			);
291			painCave.isFatal = 1;
292			simError();
293			return -1;
294			}
295
296			// We must set tauThermostat.
297
298			if (!have_tau_thermostat) {
299			sprintf( painCave.errMsg,
300			"NPT error: If you use the NPT\n"
301			" integrator, you must set tauThermostat.\n");
302			painCave.isFatal = 1;
303			simError();
304			return -1;
305			}
306
307			// We must set tauBarostat.
308
309			if (!have_tau_barostat) {
310			sprintf( painCave.errMsg,
311			"NPT error: If you use the NPT\n"
312			" integrator, you must set tauBarostat.\n");
313			painCave.isFatal = 1;
314			simError();
315			return -1;
316			}
317
318			if (!have_chi_tolerance) {
319			sprintf( painCave.errMsg,
320			"NPT warning: setting chi tolerance to 1e-6\n");
321			chiTolerance = 1e-6;
322			have_chi_tolerance = 1;
323			painCave.isFatal = 0;
324			simError();
325			}
326
327			if (!have_eta_tolerance) {
328			sprintf( painCave.errMsg,
329			"NPT warning: setting eta tolerance to 1e-6\n");
330			etaTolerance = 1e-6;
331			have_eta_tolerance = 1;
332			painCave.isFatal = 0;
333			simError();
334			}
335
336			// We need NkBT a lot, so just set it here: This is the RAW number
337			// of particles, so no subtraction or addition of constraints or
338			// orientational degrees of freedom:
339
340			NkBT = (double)Nparticles * kB * targetTemp;
341
342			// fkBT is used because the thermostat operates on more degrees of freedom
343			// than the barostat (when there are particles with orientational degrees
344			// of freedom). ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
345
346			fkBT = (double)info->ndf * kB * targetTemp;
347
348			tt2 = tauThermostat * tauThermostat;
349			tb2 = tauBarostat * tauBarostat;
350
351			return 1;
352			}