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.

NPTi::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->getPropertyByName(ETAVALUE_ID);
    if(data){
      etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
      
      if(etaValue){
        eta = (*etaValue)[0];
        oldEta = eta;
      }
    }
  }
}

NPTi::~NPTi() {
  //nothing for now
}

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

void NPTi::evolveEtaB() {

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

void NPTi::calcVelScale() {
  vScale = chi + eta;
}

void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
    sc = vel * vScale;
}

void NPTi::getVelScaleB(Vector3d& sc, int index ){
    sc = oldVel[index] * vScale;    
}


void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
                                               int index, Vector3d& sc){
    /**@todo*/
    sc  = oldPos[index] + pos/2.0 -COM;
    sc *= eta;
}

void NPTi::scaleSimBox(){

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

}

bool NPTi::etaConverged() {

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

double NPTi::calcConservedQuantity(){

  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;

  return conservedQuantity;
}

Revision:	1774
Committed:	Tue Nov 23 23:12:23 2004 UTC (19 years, 9 months ago) by tim
File size:	3094 byte(s)
Log Message:	refactory NPT integrator
#	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	tim	1774	NPTi::NPTi ( SimInfo theInfo, ForceFields the_ff):
27	gezelter	1490	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	tim	1701	data = info->getPropertyByName(ETAVALUE_ID);
39	gezelter	1490	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	tim	1774	NPTi::~NPTi() {
51	gezelter	1490	//nothing for now
52			}
53
54	tim	1774	void NPTi::evolveEtaA() {
55	gezelter	1490	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
56			(p_convertNkBTtb2));
57			oldEta = eta;
58			}
59
60	tim	1774	void NPTi::evolveEtaB() {
61	gezelter	1490
62			prevEta = eta;
63			eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
64			(p_convertNkBTtb2));
65			}
66
67	tim	1774	void NPTi::calcVelScale() {
68	gezelter	1490	vScale = chi + eta;
69			}
70
71	tim	1774	void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
72			sc = vel * vScale;
73	gezelter	1490	}
74
75	tim	1774	void NPTi::getVelScaleB(Vector3d& sc, int index ){
76			sc = oldVel[index] * vScale;
77	gezelter	1490	}
78
79
80	tim	1774	void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
81			int index, Vector3d& sc){
82			/*@todo/
83			sc = oldPos[index] + pos/2.0 -COM;
84			sc *= eta;
85	gezelter	1490	}
86
87	tim	1774	void NPTi::scaleSimBox(){
88	gezelter	1490
89			double scaleFactor;
90
91			scaleFactor = exp(dt*eta);
92
93			if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
94			sprintf( painCave.errMsg,
95			"NPTi error: Attempting a Box scaling of more than 10 percent"
96			" check your tauBarostat, as it is probably too small!\n"
97			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
98			);
99			painCave.isFatal = 1;
100			simError();
101			} else {
102			info->scaleBox(scaleFactor);
103			}
104
105			}
106
107	tim	1774	bool NPTi::etaConverged() {
108	gezelter	1490
109			return ( fabs(prevEta - eta) <= etaTolerance );
110			}
111
112	tim	1774	double NPTi::calcConservedQuantity(){
113	gezelter	1490
114			double conservedQuantity;
115			double Energy;
116			double thermostat_kinetic;
117			double thermostat_potential;
118			double barostat_kinetic;
119			double barostat_potential;
120
121			Energy = tStats->getTotalE();
122
123			thermostat_kinetic = fkBT* tt2 * chi * chi /
124			(2.0 * eConvert);
125
126			thermostat_potential = fkBT* integralOfChidt / eConvert;
127
128
129			barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /
130			(2.0 * eConvert);
131
132			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
133			eConvert;
134
135			conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
136			barostat_kinetic + barostat_potential;
137
138			return conservedQuantity;
139			}
140