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*integrableObjects.size()];
  oldVel = new double[3*integrableObjects.size()];
  oldJi = new double[3*integrableObjects.size()];

}

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;
  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

  calcVelScale();
  
  for( i=0; i<integrableObjects.size(); i++ ){

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

    mass = integrableObjects[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]);

    }

    integrableObjects[i]->setVel( vel );

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

      // get and convert the torque to body frame

      integrableObjects[i]->getTrq( Tb );
      integrableObjects[i]->lab2Body( Tb );

      // get the angular momentum, and propagate a half step

      integrableObjects[i]->getJ( ji );

      for (j=0; j < 3; j++)
        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);

      this->rotationPropagation( integrableObjects[i], ji );

      integrableObjects[i]->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 < integrableObjects.size(); i++){
    integrableObjects[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 < integrableObjects.size(); i++){

      integrableObjects[i]->getVel(vel);
      integrableObjects[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]);

      integrableObjects[i]->setPos( pos );
    }

    rattle->doRattleA();
  }


  // 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;
  double Tb[3], ji[3], sc[3];
  double vel[3], frc[3];
  double mass;

  // Set things up for the iteration:

  for( i=0; i<integrableObjects.size(); i++ ){

    integrableObjects[i]->getVel( vel );

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

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

      integrableObjects[i]->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();
    this->calcVelScale();

    for( i=0; i<integrableObjects.size(); i++ ){

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

      mass = integrableObjects[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]);

      integrableObjects[i]->setVel( vel );

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

        // get and convert the torque to body frame

        integrableObjects[i]->getTrq( Tb );
        integrableObjects[i]->lab2Body( Tb );

        for (j=0; j < 3; j++)
          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);

          integrableObjects[i]->setJ( ji );
      }
    }

    rattle->doRattleB();

    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,
             "If you use the NPT integrator, you must set tauBarostat.\n");
    painCave.severity = OOPSE_ERROR;
    painCave.isFatal = 1;
    simError();
    return -1;
  }

  if (!have_chi_tolerance) {
    sprintf( painCave.errMsg,
             "Setting chi tolerance to 1e-6 in NPT integrator\n");
    chiTolerance = 1e-6;
    have_chi_tolerance = 1;
    painCave.severity = OOPSE_INFO;
    painCave.isFatal = 0;
    simError();
  }

  if (!have_eta_tolerance) {
    sprintf( painCave.errMsg,
             "Setting eta tolerance to 1e-6 in NPT integrator");
    etaTolerance = 1e-6;
    have_eta_tolerance = 1;
    painCave.severity = OOPSE_INFO;
    painCave.isFatal = 0;
    simError();
  }

  // We need NkBT a lot, so just set it here: This is the RAW number
  // of integrableObjects, so no subtraction or addition of constraints or
  // orientational degrees of freedom:

  NkBT = (double)(info->getTotIntegrableObjects()) * 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).  

  fkBT = (double)(info->getNDF()) * kB * targetTemp;

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

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