src/integrators/NPTi.cpp

#include "NPTi.hpp"
#include "brains/SimInfo.hpp"
#include "brains/Thermo.hpp"
#include "integrators/IntegratorCreator.hpp"
#include "integrators/NPT.hpp"
#include "primitives/Molecule.hpp"
#include "utils/OOPSEConstant.hpp"
#include "utils/simError.h"

namespace oopse {

static IntegratorBuilder<NPTi>* NPTiCreator = new IntegratorBuilder<NPTi>("NPTi");

// 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 *info) : NPT(info){

    Globals* simParams = info_->getSimParams();

}

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

void NPTi::evolveEtaB() {

    prevEta = eta;
    eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
         (OOPSEConstant::pressureConvert*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 {
        Mat3x3d hmat = currentSnapshot_->getHmat();
        hmat *= scaleFactor;
        currentSnapshot_->setHmat(hmat);
    }

}

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 =thermo.getTotalE();

    thermostat_kinetic = fkBT* tt2 * chi * chi / (2.0 * OOPSEConstant::energyConvert);

    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;


    barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /(2.0 * OOPSEConstant::energyConvert);

    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /
        OOPSEConstant::energyConvert;

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

    return conservedQuantity;
}

void NPTi::loadEta() {
    Mat3x3d etaMat = currentSnapshot_->getEta();
    eta = etaMat(0,0);
    if (fabs(etaMat(1,1) - eta) >= oopse::epsilon || fabs(etaMat(1,1) - eta) >= oopse::epsilon || etaMat.isDiagonal()) {
        sprintf( painCave.errMsg,
                 "NPTi error: the diagonal elements of are eta matrix is not same or etaMat is not a diagonal matrix");
        painCave.isFatal = 1;
        simError();
    }
}

void NPTi::saveEta() {
    Mat3x3d etaMat(0.0);
    etaMat(0, 0) = eta;
    etaMat(1, 1) = eta;
    etaMat(2, 2) = eta;
    currentSnapshot_->setEta(etaMat);
}

}
Revision:	1841
Committed:	Fri Dec 3 17:59:45 2004 UTC (19 years, 9 months ago) by tim
File size:	3565 byte(s)
Log Message:	begin to fix bugs
#	User	Rev	Content
1	tim	1822	#include "NPTi.hpp"
2	tim	1492	#include "brains/SimInfo.hpp"
3			#include "brains/Thermo.hpp"
4	tim	1837	#include "integrators/IntegratorCreator.hpp"
5	tim	1822	#include "integrators/NPT.hpp"
6			#include "primitives/Molecule.hpp"
7			#include "utils/OOPSEConstant.hpp"
8	tim	1492	#include "utils/simError.h"
9	gezelter	1490
10	tim	1822	namespace oopse {
11	gezelter	1490
12	tim	1837	static IntegratorBuilder<NPTi>* NPTiCreator = new IntegratorBuilder<NPTi>("NPTi");
13
14	gezelter	1490	// Basic isotropic thermostating and barostating via the Melchionna
15			// modification of the Hoover algorithm:
16			//
17			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
18			// Molec. Phys., 78, 533.
19			//
20			// and
21			//
22			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
23
24	tim	1822	NPTi::NPTi ( SimInfo *info) : NPT(info){
25	gezelter	1490
26	tim	1841	Globals* simParams = info_->getSimParams();
27	gezelter	1490
28			}
29
30	tim	1774	void NPTi::evolveEtaA() {
31	tim	1822	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
32			(OOPSEConstant::pressureConvertNkBTtb2));
33			oldEta = eta;
34	gezelter	1490	}
35
36	tim	1774	void NPTi::evolveEtaB() {
37	gezelter	1490
38	tim	1822	prevEta = eta;
39			eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
40			(OOPSEConstant::pressureConvertNkBTtb2));
41	gezelter	1490	}
42
43	tim	1774	void NPTi::calcVelScale() {
44	tim	1822	vScale = chi + eta;
45	gezelter	1490	}
46
47	tim	1774	void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
48			sc = vel * vScale;
49	gezelter	1490	}
50
51	tim	1774	void NPTi::getVelScaleB(Vector3d& sc, int index ){
52			sc = oldVel[index] * vScale;
53	gezelter	1490	}
54
55
56	tim	1774	void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
57	tim	1822	int index, Vector3d& sc){
58	tim	1774	/*@todo/
59			sc = oldPos[index] + pos/2.0 -COM;
60			sc *= eta;
61	gezelter	1490	}
62
63	tim	1774	void NPTi::scaleSimBox(){
64	gezelter	1490
65	tim	1822	double scaleFactor;
66	gezelter	1490
67	tim	1822	scaleFactor = exp(dt*eta);
68	gezelter	1490
69	tim	1822	if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
70			sprintf( painCave.errMsg,
71	gezelter	1490	"NPTi error: Attempting a Box scaling of more than 10 percent"
72			" check your tauBarostat, as it is probably too small!\n"
73			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
74			);
75	tim	1822	painCave.isFatal = 1;
76			simError();
77			} else {
78			Mat3x3d hmat = currentSnapshot_->getHmat();
79			hmat *= scaleFactor;
80			currentSnapshot_->setHmat(hmat);
81			}
82	gezelter	1490
83			}
84
85	tim	1774	bool NPTi::etaConverged() {
86	gezelter	1490
87	tim	1822	return ( fabs(prevEta - eta) <= etaTolerance );
88	gezelter	1490	}
89
90	tim	1774	double NPTi::calcConservedQuantity(){
91	gezelter	1490
92	tim	1822	double conservedQuantity;
93			double Energy;
94			double thermostat_kinetic;
95			double thermostat_potential;
96			double barostat_kinetic;
97			double barostat_potential;
98	gezelter	1490
99	tim	1822	Energy =thermo.getTotalE();
100	gezelter	1490
101	tim	1822	thermostat_kinetic = fkBT* tt2 * chi * chi / (2.0 * OOPSEConstant::energyConvert);
102	gezelter	1490
103	tim	1822	thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
104	gezelter	1490
105
106	tim	1822	barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /(2.0 * OOPSEConstant::energyConvert);
107	gezelter	1490
108	tim	1822	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /
109			OOPSEConstant::energyConvert;
110	gezelter	1490
111	tim	1822	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
112			barostat_kinetic + barostat_potential;
113	gezelter	1490
114	tim	1822	return conservedQuantity;
115	gezelter	1490	}
116
117	tim	1822	void NPTi::loadEta() {
118			Mat3x3d etaMat = currentSnapshot_->getEta();
119			eta = etaMat(0,0);
120			if (fabs(etaMat(1,1) - eta) >= oopse::epsilon \|\| fabs(etaMat(1,1) - eta) >= oopse::epsilon \|\| etaMat.isDiagonal()) {
121			sprintf( painCave.errMsg,
122			"NPTi error: the diagonal elements of are eta matrix is not same or etaMat is not a diagonal matrix");
123			painCave.isFatal = 1;
124			simError();
125			}
126			}
127
128			void NPTi::saveEta() {
129			Mat3x3d etaMat(0.0);
130			etaMat(0, 0) = eta;
131			etaMat(1, 1) = eta;
132			etaMat(2, 2) = eta;
133			currentSnapshot_->setEta(etaMat);
134			}
135
136			}