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