OOPSE/libmdtools/NPTi.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> NPTi<T>::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  GenericData* data;
  DoubleArrayData * 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<DoubleArrayData*>(data);
      
      if(etaValue){
        etaArray = etaValue->getData();
        eta = etaArray[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>::getVelScaleA(double sc[3], double vel[3]) {
  int i;

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

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] * ( chi + eta );
}


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