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 {

// 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){

}

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(){

    chi= currentSnapshot_->getChi();
    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
    loadEta();
    // We need NkBT a lot, so just set it here: This is the RAW number
    // of integrableObjects, so no subtraction or addition of constraints or
    // orientational degrees of freedom:
    NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp;

    // fkBT is used because the thermostat operates on more degrees of freedom
    // than the barostat (when there are particles with orientational degrees
    // of freedom).  
    fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
    
    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;

    std::cout << "--------------------------------------------------------------" << std::endl;
    std::cout << "time: " << currentSnapshot_->getTime() << std:: endl;
    std::cout << "chi : " << chi << std::endl;
    std::cout << "integralOfChidt : " << integralOfChidt << std::endl;
    std::cout << "eta : " << eta << std::endl;    
    std::cout << "NkBT: " << NkBT << std::endl;
    std::cout << "fkBT: " << fkBT << std::endl;    
    std::cout << "thermostat_kinetic : " << thermostat_kinetic<< std::endl;
    std::cout << "thermostat_potential : " << thermostat_potential << std::endl;
    std::cout << "barostat_kinetic : " << barostat_kinetic << std::endl;
    std::cout << "barostat_potential : " << barostat_potential << std::endl;
    std::cout << "Total Energy: " << Energy << std::endl;
    std::cout << "Conserved Quantity: " << conservedQuantity <<std::endl;
    std::cout << "--------------------------------------------------------------" << std::endl;
    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  eta matrix are not the 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:	1897
Committed:	Mon Dec 20 19:49:44 2004 UTC (21 years, 6 months ago) by tim
File size:	5050 byte(s)
Log Message:	NPT is working now
#	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			// Basic isotropic thermostating and barostating via the Melchionna
13			// modification of the Hoover algorithm:
14			//
15			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
16			// Molec. Phys., 78, 533.
17			//
18			// and
19			//
20			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
21
22	tim	1822	NPTi::NPTi ( SimInfo *info) : NPT(info){
23	gezelter	1490
24			}
25
26	tim	1774	void NPTi::evolveEtaA() {
27	tim	1822	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
28			(OOPSEConstant::pressureConvertNkBTtb2));
29			oldEta = eta;
30	gezelter	1490	}
31
32	tim	1774	void NPTi::evolveEtaB() {
33	gezelter	1490
34	tim	1822	prevEta = eta;
35			eta = oldEta + dt2 * ( instaVol * (instaPress - targetPressure) /
36			(OOPSEConstant::pressureConvertNkBTtb2));
37	gezelter	1490	}
38
39	tim	1774	void NPTi::calcVelScale() {
40	tim	1822	vScale = chi + eta;
41	gezelter	1490	}
42
43	tim	1774	void NPTi::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
44			sc = vel * vScale;
45	gezelter	1490	}
46
47	tim	1774	void NPTi::getVelScaleB(Vector3d& sc, int index ){
48			sc = oldVel[index] * vScale;
49	gezelter	1490	}
50
51
52	tim	1774	void NPTi::getPosScale(const Vector3d& pos, const Vector3d& COM,
53	tim	1822	int index, Vector3d& sc){
54	tim	1774	/*@todo/
55	tim	1897	sc = (oldPos[index] + pos)/2.0 -COM;
56	tim	1774	sc *= eta;
57	gezelter	1490	}
58
59	tim	1774	void NPTi::scaleSimBox(){
60	gezelter	1490
61	tim	1822	double scaleFactor;
62	gezelter	1490
63	tim	1822	scaleFactor = exp(dt*eta);
64	gezelter	1490
65	tim	1822	if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
66			sprintf( painCave.errMsg,
67	gezelter	1490	"NPTi error: Attempting a Box scaling of more than 10 percent"
68			" check your tauBarostat, as it is probably too small!\n"
69			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
70			);
71	tim	1822	painCave.isFatal = 1;
72			simError();
73			} else {
74			Mat3x3d hmat = currentSnapshot_->getHmat();
75			hmat *= scaleFactor;
76			currentSnapshot_->setHmat(hmat);
77			}
78	gezelter	1490
79			}
80
81	tim	1774	bool NPTi::etaConverged() {
82	gezelter	1490
83	tim	1822	return ( fabs(prevEta - eta) <= etaTolerance );
84	gezelter	1490	}
85
86	tim	1774	double NPTi::calcConservedQuantity(){
87	gezelter	1490
88	tim	1867	chi= currentSnapshot_->getChi();
89			integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
90			loadEta();
91			// We need NkBT a lot, so just set it here: This is the RAW number
92			// of integrableObjects, so no subtraction or addition of constraints or
93			// orientational degrees of freedom:
94			NkBT = info_->getNGlobalIntegrableObjects()OOPSEConstant::kB targetTemp;
95
96			// fkBT is used because the thermostat operates on more degrees of freedom
97			// than the barostat (when there are particles with orientational degrees
98			// of freedom).
99			fkBT = info_->getNdf()OOPSEConstant::kB targetTemp;
100
101	tim	1822	double conservedQuantity;
102			double Energy;
103			double thermostat_kinetic;
104			double thermostat_potential;
105			double barostat_kinetic;
106			double barostat_potential;
107	gezelter	1490
108	tim	1822	Energy =thermo.getTotalE();
109	gezelter	1490
110	tim	1822	thermostat_kinetic = fkBT* tt2 * chi * chi / (2.0 * OOPSEConstant::energyConvert);
111	gezelter	1490
112	tim	1822	thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
113	gezelter	1490
114
115	tim	1822	barostat_kinetic = 3.0 * NkBT * tb2 * eta * eta /(2.0 * OOPSEConstant::energyConvert);
116	gezelter	1490
117	tim	1822	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /
118			OOPSEConstant::energyConvert;
119	gezelter	1490
120	tim	1822	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential +
121			barostat_kinetic + barostat_potential;
122	gezelter	1490
123	tim	1867	std::cout << "--------------------------------------------------------------" << std::endl;
124			std::cout << "time: " << currentSnapshot_->getTime() << std:: endl;
125			std::cout << "chi : " << chi << std::endl;
126			std::cout << "integralOfChidt : " << integralOfChidt << std::endl;
127			std::cout << "eta : " << eta << std::endl;
128			std::cout << "NkBT: " << NkBT << std::endl;
129			std::cout << "fkBT: " << fkBT << std::endl;
130			std::cout << "thermostat_kinetic : " << thermostat_kinetic<< std::endl;
131			std::cout << "thermostat_potential : " << thermostat_potential << std::endl;
132			std::cout << "barostat_kinetic : " << barostat_kinetic << std::endl;
133			std::cout << "barostat_potential : " << barostat_potential << std::endl;
134			std::cout << "Total Energy: " << Energy << std::endl;
135			std::cout << "Conserved Quantity: " << conservedQuantity <<std::endl;
136			std::cout << "--------------------------------------------------------------" << std::endl;
137	tim	1822	return conservedQuantity;
138	gezelter	1490	}
139
140	tim	1822	void NPTi::loadEta() {
141			Mat3x3d etaMat = currentSnapshot_->getEta();
142			eta = etaMat(0,0);
143	tim	1883	//if (fabs(etaMat(1,1) - eta) >= oopse::epsilon \|\| fabs(etaMat(1,1) - eta) >= oopse::epsilon \|\| !etaMat.isDiagonal()) {
144			// sprintf( painCave.errMsg,
145			// "NPTi error: the diagonal elements of eta matrix are not the same or etaMat is not a diagonal matrix");
146			// painCave.isFatal = 1;
147			// simError();
148			//}
149	tim	1822	}
150
151			void NPTi::saveEta() {
152			Mat3x3d etaMat(0.0);
153			etaMat(0, 0) = eta;
154			etaMat(1, 1) = eta;
155			etaMat(2, 2) = eta;
156			currentSnapshot_->setEta(etaMat);
157			}
158
159			}