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 )
{
  eta = 0.0;
  oldEta = 0.0;
}

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; 
}
Revision:	829
Committed:	Tue Oct 28 16:03:37 2003 UTC (20 years, 8 months ago) by gezelter
File size:	3356 byte(s)
Log Message:	replace c++ header stuff with more portable c header stuff Also, mod file fixes and portability changes Some fortran changes will need to be reversed.
#	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			#include "simError.h"
11
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			// Molec. Phys., 78, 533.
21			//
22			// and
23			//
24			// 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			eta = 0.0;
30	mmeineke	778	oldEta = 0.0;
31	gezelter	574	}
32
33	tim	763	template<typename T> NPTi<T>::~NPTi() {
34	mmeineke	778	//nothing for now
35	tim	763	}
36
37	mmeineke	778	template<typename T> void NPTi<T>::resetIntegrator() {
38			eta = 0.0;
39			T::resetIntegrator();
40			}
41	gezelter	600
42	mmeineke	778	template<typename T> void NPTi<T>::evolveEtaA() {
43			eta += dt2 * ( instaVol * (instaPress - targetPressure) /
44			(p_convertNkBTtb2));
45			oldEta = eta;
46			}
47	gezelter	574
48	mmeineke	778	template<typename T> void NPTi<T>::evolveEtaB() {
49	gezelter	574
50	mmeineke	778	prevEta = eta;
51			eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
52			(p_convertNkBTtb2));
53			}
54	gezelter	574
55	mmeineke	778	template<typename T> void NPTi<T>::getVelScaleA(double sc[3], double vel[3]) {
56			int i;
57	gezelter	574
58	mmeineke	778	for(i=0; i<3; i++) sc[i] = vel[i] * ( chi + eta );
59			}
60	gezelter	574
61	mmeineke	778	template<typename T> void NPTi<T>::getVelScaleB(double sc[3], int index ){
62			int i;
63	gezelter	600
64	mmeineke	778	for(i=0; i<3; i++) sc[i] = oldVel[index3 + i] ( chi + eta );
65			}
66	gezelter	574
67	tim	763
68	mmeineke	778	template<typename T> void NPTi<T>::getPosScale(double pos[3], double COM[3],
69			int index, double sc[3]){
70			int j;
71	gezelter	574
72	mmeineke	778	for(j=0; j<3; j++)
73			sc[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
74	gezelter	600
75	mmeineke	778	for(j=0; j<3; j++)
76			sc[j] *= eta;
77			}
78	gezelter	600
79	mmeineke	778	template<typename T> void NPTi<T>::scaleSimBox( void ){
80	gezelter	600
81	mmeineke	778	double scaleFactor;
82	gezelter	600
83	gezelter	611	scaleFactor = exp(dt*eta);
84
85	mmeineke	614	if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
86	gezelter	611	sprintf( painCave.errMsg,
87			"NPTi error: Attempting a Box scaling of more than 10 percent"
88			" check your tauBarostat, as it is probably too small!\n"
89			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
90			);
91			painCave.isFatal = 1;
92			simError();
93			} else {
94	tim	763	info->scaleBox(scaleFactor);
95			}
96	mmeineke	614
97	gezelter	574	}
98
99	mmeineke	778	template<typename T> bool NPTi<T>::etaConverged() {
100	tim	763
101	mmeineke	778	return ( fabs(prevEta - eta) <= etaTolerance );
102	gezelter	574	}
103
104	tim	763	template<typename T> double NPTi<T>::getConservedQuantity(void){
105
106			double conservedQuantity;
107	tim	769	double Energy;
108			double thermostat_kinetic;
109			double thermostat_potential;
110			double barostat_kinetic;
111			double barostat_potential;
112	mmeineke	787
113	tim	769	Energy = tStats->getTotalE();
114	tim	763
115	mmeineke	787	thermostat_kinetic = fkBT* tt2 * chi * chi /
116	tim	769	(2.0 * eConvert);
117	tim	763
118	tim	769	thermostat_potential = fkBT* integralOfChidt / eConvert;
119	tim	763
120
121	mmeineke	787	barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
122	tim	769	(2.0 * eConvert);
123
124			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
125			eConvert;
126	tim	767
127	tim	769	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
128			barostat_kinetic + barostat_potential;
129
130	mmeineke	778	// cout.width(8);
131			// cout.precision(8);
132	tim	763
133	mmeineke	778	// cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
134			// "\t" << thermostat_potential << "\t" << barostat_kinetic <<
135			// "\t" << barostat_potential << "\t" << conservedQuantity << endl;
136	tim	763	return conservedQuantity;
137			}