OOPSE/libmdtools/NPT.cpp

#include <math.h>
#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 )
{
  GenericData* data;
  DoubleData * chiValue;
  DoubleData * integralOfChidtValue;

  chiValue = NULL;
  integralOfChidtValue = NULL;

  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;

  // retrieve chi and integralOfChidt from simInfo
  data = info->getProperty(CHIVALUE_ID);
  if(data){
    chiValue = dynamic_cast<DoubleData*>(data);
  }

  data = info->getProperty(INTEGRALOFCHIDT_ID);
  if(data){
    integralOfChidtValue = dynamic_cast<DoubleData*>(data);
  }

  // chi and integralOfChidt should appear by pair
  if(chiValue && integralOfChidtValue){
    chi = chiValue->getData();
    integralOfChidt = integralOfChidtValue->getData();
  }

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