src/integrators/NPT.cpp

#include <math.h>

#include "primitives/Atom.hpp"
#include "primitives/SRI.hpp"
#include "primitives/AbstractClasses.hpp"
#include "brains/SimInfo.hpp"
#include "UseTheForce/ForceFields.hpp"
#include "brains/Thermo.hpp"
#include "io/ReadWrite.hpp"
#include "integrators/Integrator.hpp"
#include "utils/simError.h"

#ifdef IS_MPI
#include "brains/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;
  DoubleGenericData * chiValue;
  DoubleGenericData * 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<DoubleGenericData*>(data);
  }

  data = info->getProperty(INTEGRALOFCHIDT_ID);
  if(data){
    integralOfChidtValue = dynamic_cast<DoubleGenericData*>(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;
  Vector3d Tb, ji;
  double mass;
  Vector3d vel;
  Vector3d pos;
  Vector3d frc;
  Vector3d sc;
  Vector3d COM;

  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++ ){

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

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

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

      ji = integrableObjects[i]->getJ();

      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++){
    pos = integrableObjects[i]->getPos();
    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++){

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

      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 );
    }
    
    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;
  Vector3d Tb;
  Vector3d ji;
  Vector3d sc;
  Vector3d vel;
  Vector3d frc;
  double mass;

  // Set things up for the iteration:

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

    vel = integrableObjects[i]->getVel();

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

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

      ji = integrableObjects[i]->getJ();

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

      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

        Tb = integrableObjects[i]->getTrq();
        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 );
      }
    }

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