src/integrators/NPTi.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> NPTi<T>::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  GenericData* data;
  DoubleVectorGenericData * etaValue;
  vector<double> etaArray;

  eta = 0.0;
  oldEta = 0.0;

  if( theInfo->useInitXSstate ){
    // retrieve eta from simInfo if
    data = info->getProperty(ETAVALUE_ID);
    if(data){
      etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
      
      if(etaValue){
        eta = (*etaValue)[0];
        oldEta = eta;
      }
    }
  }
}

template<typename T> NPTi<T>::~NPTi() {
  //nothing for now
}

template<typename T> void NPTi<T>::resetIntegrator() {
  eta = 0.0;
  T::resetIntegrator();
}

template<typename T> void NPTi<T>::evolveEtaA() {
  eta += dt2 * ( instaVol * (instaPress - targetPressure) /
                 (p_convert*NkBT*tb2));
  oldEta = eta;
}

template<typename T> void NPTi<T>::evolveEtaB() {

  prevEta = eta;
  eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
                 (p_convert*NkBT*tb2));
}

template<typename T> void NPTi<T>::calcVelScale(void) {
  vScale = chi + eta;
}

template<typename T> void NPTi<T>::getVelScaleA(double sc[3], double vel[3]) {
  int i;

  for(i=0; i<3; i++) sc[i] = vel[i] * vScale;
}

template<typename T> void NPTi<T>::getVelScaleB(double sc[3], int index ){
  int i;

  for(i=0; i<3; i++) sc[i] = oldVel[index*3 + i] * vScale;
}


template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
                                               int index, double sc[3]){
  int j;

  for(j=0; j<3; j++)
    sc[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];

  for(j=0; j<3; j++)
    sc[j] *= eta;
}

template<typename T> void NPTi<T>::scaleSimBox( void ){

  double scaleFactor;

  scaleFactor = exp(dt*eta);

  if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
    sprintf( painCave.errMsg,
             "NPTi error: Attempting a Box scaling of more than 10 percent"
             " check your tauBarostat, as it is probably too small!\n"
             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
             );
    painCave.isFatal = 1;
    simError();
  } else {
    info->scaleBox(scaleFactor);
  }

}

template<typename T> bool NPTi<T>::etaConverged() {

  return ( fabs(prevEta - eta) <= etaTolerance );
}

template<typename T> double NPTi<T>::getConservedQuantity(void){

  double conservedQuantity;
  double Energy;
  double thermostat_kinetic;
  double thermostat_potential;
  double barostat_kinetic;
  double barostat_potential;

  Energy = tStats->getTotalE();

  thermostat_kinetic = fkBT* tt2 * chi * chi /
    (2.0 * eConvert);

  thermostat_potential = fkBT* integralOfChidt / eConvert;


  barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
    (2.0 * eConvert);

  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
    eConvert;

  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
    barostat_kinetic + barostat_potential;

//   cout.width(8);
//   cout.precision(8);

//   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
//       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
//       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
  return conservedQuantity;
}

template<typename T> string NPTi<T>::getAdditionalParameters(void){
  string parameters;
  const int BUFFERSIZE = 2000; // size of the read buffer
  char buffer[BUFFERSIZE];

  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
  parameters += buffer;

  sprintf(buffer,"\t%G\t0\t0;", eta);
  parameters += buffer;

  sprintf(buffer,"\t0\t%G\t0;", eta);
  parameters += buffer;

  sprintf(buffer,"\t0\t0\t%G;", eta);
  parameters += buffer;

  return parameters;

}
Revision:	1625
Committed:	Thu Oct 21 16:22:01 2004 UTC (19 years, 8 months ago) by tim
File size:	4300 byte(s)
Log Message:	replace old GebericData with new GenericData
#	User	Rev	Content
1	gezelter	1490	#include <math.h>
2	tim	1492	#include "primitives/Atom.hpp"
3			#include "primitives/SRI.hpp"
4			#include "primitives/AbstractClasses.hpp"
5			#include "brains/SimInfo.hpp"
6			#include "UseTheForce/ForceFields.hpp"
7			#include "brains/Thermo.hpp"
8			#include "io/ReadWrite.hpp"
9			#include "integrators/Integrator.hpp"
10			#include "utils/simError.h"
11	gezelter	1490
12			#ifdef IS_MPI
13	tim	1492	#include "brains/mpiSimulation.hpp"
14	gezelter	1490	#endif
15
16			// Basic isotropic thermostating and barostating via the Melchionna
17			// modification of the Hoover algorithm:
18			//
19			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
20			// Molec. Phys., 78, 533.
21			//
22			// and
23			//
24			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
25
26			template<typename T> NPTi<T>::NPTi ( SimInfo theInfo, ForceFields the_ff):
27			T( theInfo, the_ff )
28			{
29			GenericData* data;
30	tim	1625	DoubleVectorGenericData * etaValue;
31	gezelter	1490	vector<double> etaArray;
32
33			eta = 0.0;
34			oldEta = 0.0;
35
36			if( theInfo->useInitXSstate ){
37			// retrieve eta from simInfo if
38			data = info->getProperty(ETAVALUE_ID);
39			if(data){
40	tim	1625	etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
41	gezelter	1490
42			if(etaValue){
43	tim	1625	eta = (*etaValue)[0];
44	gezelter	1490	oldEta = eta;
45			}
46			}
47			}
48			}
49
50			template<typename T> NPTi<T>::~NPTi() {
51			//nothing for now
52			}
53
54			template<typename T> void NPTi<T>::resetIntegrator() {
55			eta = 0.0;
56			T::resetIntegrator();
57			}
58
59			template<typename T> void NPTi<T>::evolveEtaA() {
60			eta += dt2 * ( instaVol * (instaPress - targetPressure) /
61			(p_convertNkBTtb2));
62			oldEta = eta;
63			}
64
65			template<typename T> void NPTi<T>::evolveEtaB() {
66
67			prevEta = eta;
68			eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
69			(p_convertNkBTtb2));
70			}
71
72			template<typename T> void NPTi<T>::calcVelScale(void) {
73			vScale = chi + eta;
74			}
75
76			template<typename T> void NPTi<T>::getVelScaleA(double sc[3], double vel[3]) {
77			int i;
78
79			for(i=0; i<3; i++) sc[i] = vel[i] * vScale;
80			}
81
82			template<typename T> void NPTi<T>::getVelScaleB(double sc[3], int index ){
83			int i;
84
85			for(i=0; i<3; i++) sc[i] = oldVel[index3 + i] vScale;
86			}
87
88
89			template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
90			int index, double sc[3]){
91			int j;
92
93			for(j=0; j<3; j++)
94			sc[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
95
96			for(j=0; j<3; j++)
97			sc[j] *= eta;
98			}
99
100			template<typename T> void NPTi<T>::scaleSimBox( void ){
101
102			double scaleFactor;
103
104			scaleFactor = exp(dt*eta);
105
106			if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
107			sprintf( painCave.errMsg,
108			"NPTi error: Attempting a Box scaling of more than 10 percent"
109			" check your tauBarostat, as it is probably too small!\n"
110			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
111			);
112			painCave.isFatal = 1;
113			simError();
114			} else {
115			info->scaleBox(scaleFactor);
116			}
117
118			}
119
120			template<typename T> bool NPTi<T>::etaConverged() {
121
122			return ( fabs(prevEta - eta) <= etaTolerance );
123			}
124
125			template<typename T> double NPTi<T>::getConservedQuantity(void){
126
127			double conservedQuantity;
128			double Energy;
129			double thermostat_kinetic;
130			double thermostat_potential;
131			double barostat_kinetic;
132			double barostat_potential;
133
134			Energy = tStats->getTotalE();
135
136			thermostat_kinetic = fkBT* tt2 * chi * chi /
137			(2.0 * eConvert);
138
139			thermostat_potential = fkBT* integralOfChidt / eConvert;
140
141
142			barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
143			(2.0 * eConvert);
144
145			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
146			eConvert;
147
148			conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
149			barostat_kinetic + barostat_potential;
150
151			// cout.width(8);
152			// cout.precision(8);
153
154			// cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
155			// "\t" << thermostat_potential << "\t" << barostat_kinetic <<
156			// "\t" << barostat_potential << "\t" << conservedQuantity << endl;
157			return conservedQuantity;
158			}
159
160			template<typename T> string NPTi<T>::getAdditionalParameters(void){
161			string parameters;
162			const int BUFFERSIZE = 2000; // size of the read buffer
163			char buffer[BUFFERSIZE];
164
165			sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
166			parameters += buffer;
167
168			sprintf(buffer,"\t%G\t0\t0;", eta);
169			parameters += buffer;
170
171			sprintf(buffer,"\t0\t%G\t0;", eta);
172			parameters += buffer;
173
174			sprintf(buffer,"\t0\t0\t%G;", eta);
175			parameters += buffer;
176
177			return parameters;
178
179			}