src/integrators/NPT.cpp

#include <math.h>

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

// 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.

namespace oopse {

NPT::NPT(SimInfo* info) :
    VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6) {

    Globals* globals = info_->getGlobals();
    
    if (globals->getUseInitXSstate()) {
        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
        currSnapshot->setChi(0.0);
        currSnapshot->setIntegralOfChiDt(0.0);
        currSnapshot->setEta(Mat3x3d(0.0));
    }
    
    if (!globals->haveTargetTemp()) {
        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
        painCave.isFatal = 1;
        painCave.severity = OOPSE_ERROR;
        simError();
    } else {
        targetTemp = globals->getTargetTemp();
    }

    // We must set tauThermostat
    if (!globals->haveTauThermostat()) {
        sprintf(painCave.errMsg, "If you use the constant temperature\n"
                                     "\tintegrator, you must set tauThermostat_.\n");

        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();
    } else {
        tauThermostat = globals->getTauThermostat();
    }

    if (!globals->haveTargetPressure()) {
        sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
                                     "   without a targetPressure!\n");

        painCave.isFatal = 1;
        simError();
    } else {
        targetPressure = globals->getTargetPressure();
    }
    
    if (!globals->haveTauBarostat()) {
        sprintf(painCave.errMsg,
                "If you use the NPT integrator, you must set tauBarostat.\n");
        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();
    } else {
        tauBarostat = globals->getTauBarostat();
    }
    
    tt2 = tauThermostat * tauThermostat;
    tb2 = tauBarostat * tauBarostat;

    update();
}

NPT::~NPT() {
}

void NPT::doUpdate() {
    VelocityVerletIntegrator::update();
    oldPos.resize(info_->getNIntegrableObjects());
    oldVel.resize(info_->getNIntegrableObjects());
    oldJi.resize(info_->getNIntegrableObjects());
    // 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;    
}

void NPT::moveA() {
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    Vector3d Tb, ji;
    double mass;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;
    Vector3d sc;
    int index;

    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    chi= currSnapshot->getChi();
    integralOfChidt = currSnapshot->getIntegralOfChiDt();
    loadEta();
    
    instaTemp =thermo.getTemperature();
    press = thermo.getPressureTensor();
    instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
    instaVol =thermo.getVolume();

    Vector3d  COM = info_->getCom();

    //evolve velocity half step

    calcVelScale();

    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
               integrableObject = mol->nextIntegrableObject(j)) {
                
            vel = integrableObject->getVel();
            frc = integrableObject->getFrc();

            mass = integrableObject->getMass();

            getVelScaleA(sc, vel);

            // velocity half step  (use chi from previous step here):
            //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
            vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
            integrableObject->setVel(vel);

            if (integrableObject->isDirectional()) {

                // get and convert the torque to body frame

                Tb = integrableObject->getTrq();
                integrableObject->lab2Body(Tb);

                // get the angular momentum, and propagate a half step

                ji = integrableObject->getJ();

                //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
                ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
                
                rotAlgo->rotate(integrableObject, ji, dt);

                integrableObject->setJ(ji);
            }
            
        }
    }
    // evolve chi and eta  half step

    chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
    
    evolveEtaA();

    //calculate the integral of chidt
    integralOfChidt += dt2 * chi;
    
    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
               integrableObject = mol->nextIntegrableObject(j)) {
            oldPos[index++] = integrableObject->getPos();            
        }
    }
    
    //the first estimation of r(t+dt) is equal to  r(t)

    for(int k = 0; k < maxIterNum_; k++) {
        index = 0;
        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
                   integrableObject = mol->nextIntegrableObject(j)) {

                vel = integrableObject->getVel();
                pos = integrableObject->getPos();

                this->getPosScale(pos, COM, index, sc);

                pos = oldPos[index] + dt * (vel + sc);
                integrableObject->setPos(pos);     

                ++index;
           }
        }

        //constraintAlgorithm->doConstrainA();
    }

    // Scale the box after all the positions have been moved:

    this->scaleSimBox();

    currSnapshot->setChi(chi);
    currSnapshot->setIntegralOfChiDt(integralOfChidt);

    saveEta();
}

void NPT::moveB(void) {
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    int index;
    Vector3d Tb;
    Vector3d ji;
    Vector3d sc;
    Vector3d vel;
    Vector3d frc;
    double mass;


    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    chi= currSnapshot->getChi();
    integralOfChidt = currSnapshot->getIntegralOfChiDt();
    double oldChi  = chi;
    double prevChi;

    loadEta();
    
    //save velocity and angular momentum
    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
               integrableObject = mol->nextIntegrableObject(j)) {
                
            oldVel[index] = integrableObject->getVel();
            oldJi[index] = integrableObject->getJ();
            ++index;
        }
    }

    // do the iteration:
    instaVol =thermo.getVolume();

    for(int k = 0; k < maxIterNum_; k++) {
        instaTemp =thermo.getTemperature();
        instaPress =thermo.getPressure();

        // evolve chi another half step using the temperature at t + dt/2
        prevChi = chi;
        chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;

        //evolve eta
        this->evolveEtaB();
        this->calcVelScale();

        index = 0;
        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
                   integrableObject = mol->nextIntegrableObject(j)) {            

                frc = integrableObject->getFrc();
                vel = integrableObject->getVel();

                mass = integrableObject->getMass();

                getVelScaleB(sc, index);

                // velocity half step
                //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
                vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
                integrableObject->setVel(vel);

                if (integrableObject->isDirectional()) {
                    // get and convert the torque to body frame
                    Tb = integrableObject->getTrq();
                    integrableObject->lab2Body(Tb);

                    //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
                    ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
                    integrableObject->setJ(ji);
                }

                ++index;
            }
        }
        
        //constraintAlgorithm->doConstrainB();

        if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
            break;
    }

    //calculate integral of chidt
    integralOfChidt += dt2 * chi;

    currSnapshot->setChi(chi);
    currSnapshot->setIntegralOfChiDt(integralOfChidt);    

    saveEta();
}

}
Revision:	1822
Committed:	Thu Dec 2 02:08:29 2004 UTC (19 years, 9 months ago) by tim
File size:	9804 byte(s)
Log Message:	oopse get compiled, still has some linking problem
#	User	Rev	Content
1	tim	1774	#include <math.h>
2
3			#include "brains/SimInfo.hpp"
4			#include "brains/Thermo.hpp"
5			#include "integrators/NPT.hpp"
6	tim	1822	#include "math/SquareMatrix3.hpp"
7			#include "primitives/Molecule.hpp"
8			#include "utils/OOPSEConstant.hpp"
9	tim	1774	#include "utils/simError.h"
10
11			// Basic isotropic thermostating and barostating via the Melchionna
12			// modification of the Hoover algorithm:
13			//
14			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
15			// Molec. Phys., 78, 533.
16			//
17			// and
18			//
19			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
20
21			namespace oopse {
22
23			NPT::NPT(SimInfo* info) :
24			VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6) {
25
26			Globals* globals = info_->getGlobals();
27
28			if (globals->getUseInitXSstate()) {
29			Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
30			currSnapshot->setChi(0.0);
31			currSnapshot->setIntegralOfChiDt(0.0);
32	tim	1822	currSnapshot->setEta(Mat3x3d(0.0));
33	tim	1774	}
34
35			if (!globals->haveTargetTemp()) {
36			sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
37			painCave.isFatal = 1;
38			painCave.severity = OOPSE_ERROR;
39			simError();
40			} else {
41			targetTemp = globals->getTargetTemp();
42			}
43
44			// We must set tauThermostat
45			if (!globals->haveTauThermostat()) {
46			sprintf(painCave.errMsg, "If you use the constant temperature\n"
47			"\tintegrator, you must set tauThermostat_.\n");
48
49			painCave.severity = OOPSE_ERROR;
50			painCave.isFatal = 1;
51			simError();
52			} else {
53			tauThermostat = globals->getTauThermostat();
54			}
55
56			if (!globals->haveTargetPressure()) {
57			sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
58			" without a targetPressure!\n");
59
60			painCave.isFatal = 1;
61			simError();
62			} else {
63			targetPressure = globals->getTargetPressure();
64			}
65
66			if (!globals->haveTauBarostat()) {
67			sprintf(painCave.errMsg,
68			"If you use the NPT integrator, you must set tauBarostat.\n");
69			painCave.severity = OOPSE_ERROR;
70			painCave.isFatal = 1;
71			simError();
72			} else {
73			tauBarostat = globals->getTauBarostat();
74			}
75
76			tt2 = tauThermostat * tauThermostat;
77			tb2 = tauBarostat * tauBarostat;
78
79			update();
80			}
81
82			NPT::~NPT() {
83			}
84
85	tim	1822	void NPT::doUpdate() {
86			VelocityVerletIntegrator::update();
87	tim	1774	oldPos.resize(info_->getNIntegrableObjects());
88			oldVel.resize(info_->getNIntegrableObjects());
89			oldJi.resize(info_->getNIntegrableObjects());
90			// We need NkBT a lot, so just set it here: This is the RAW number
91			// of integrableObjects, so no subtraction or addition of constraints or
92			// orientational degrees of freedom:
93	tim	1822	NkBT = info_->getNGlobalIntegrableObjects()OOPSEConstant::kB targetTemp;
94	tim	1774
95			// fkBT is used because the thermostat operates on more degrees of freedom
96			// than the barostat (when there are particles with orientational degrees
97			// of freedom).
98	tim	1822	fkBT = info_->getNdf()OOPSEConstant::kB targetTemp;
99	tim	1774	}
100
101			void NPT::moveA() {
102	tim	1822	SimInfo::MoleculeIterator i;
103			Molecule::IntegrableObjectIterator j;
104	tim	1774	Molecule* mol;
105			StuntDouble* integrableObject;
106			Vector3d Tb, ji;
107			double mass;
108			Vector3d vel;
109			Vector3d pos;
110			Vector3d frc;
111			Vector3d sc;
112			int index;
113	tim	1822
114			Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
115			chi= currSnapshot->getChi();
116			integralOfChidt = currSnapshot->getIntegralOfChiDt();
117			loadEta();
118	tim	1774
119	tim	1822	instaTemp =thermo.getTemperature();
120			press = thermo.getPressureTensor();
121			instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
122			instaVol =thermo.getVolume();
123	tim	1774
124	tim	1822	Vector3d COM = info_->getCom();
125	tim	1774
126			//evolve velocity half step
127
128			calcVelScale();
129
130			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
131			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
132			integrableObject = mol->nextIntegrableObject(j)) {
133
134			vel = integrableObject->getVel();
135			frc = integrableObject->getFrc();
136
137			mass = integrableObject->getMass();
138
139			getVelScaleA(sc, vel);
140
141			// velocity half step (use chi from previous step here):
142	tim	1822	//vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
143			vel += dt2OOPSEConstant::energyConvert/mass frc - dt2*sc;
144	tim	1774	integrableObject->setVel(vel);
145
146			if (integrableObject->isDirectional()) {
147
148			// get and convert the torque to body frame
149
150			Tb = integrableObject->getTrq();
151			integrableObject->lab2Body(Tb);
152
153			// get the angular momentum, and propagate a half step
154
155			ji = integrableObject->getJ();
156
157	tim	1822	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
158			ji += dt2OOPSEConstant::energyConvert Tb - dt2chi ji;
159	tim	1774
160	tim	1822	rotAlgo->rotate(integrableObject, ji, dt);
161	tim	1774
162			integrableObject->setJ(ji);
163			}
164
165			}
166			}
167			// evolve chi and eta half step
168
169	tim	1822	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
170
171	tim	1774	evolveEtaA();
172
173			//calculate the integral of chidt
174			integralOfChidt += dt2 * chi;
175
176			index = 0;
177			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
178			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
179			integrableObject = mol->nextIntegrableObject(j)) {
180			oldPos[index++] = integrableObject->getPos();
181			}
182			}
183
184			//the first estimation of r(t+dt) is equal to r(t)
185
186			for(int k = 0; k < maxIterNum_; k++) {
187			index = 0;
188			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
189			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
190			integrableObject = mol->nextIntegrableObject(j)) {
191
192			vel = integrableObject->getVel();
193			pos = integrableObject->getPos();
194
195			this->getPosScale(pos, COM, index, sc);
196
197			pos = oldPos[index] + dt * (vel + sc);
198			integrableObject->setPos(pos);
199
200			++index;
201			}
202			}
203
204			//constraintAlgorithm->doConstrainA();
205			}
206
207			// Scale the box after all the positions have been moved:
208
209			this->scaleSimBox();
210	tim	1822
211			currSnapshot->setChi(chi);
212			currSnapshot->setIntegralOfChiDt(integralOfChidt);
213
214			saveEta();
215	tim	1774	}
216
217			void NPT::moveB(void) {
218	tim	1822	SimInfo::MoleculeIterator i;
219			Molecule::IntegrableObjectIterator j;
220	tim	1774	Molecule* mol;
221			StuntDouble* integrableObject;
222			int index;
223			Vector3d Tb;
224			Vector3d ji;
225			Vector3d sc;
226			Vector3d vel;
227			Vector3d frc;
228			double mass;
229
230	tim	1822
231			Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
232			chi= currSnapshot->getChi();
233			integralOfChidt = currSnapshot->getIntegralOfChiDt();
234			double oldChi = chi;
235			double prevChi;
236
237			loadEta();
238
239	tim	1774	//save velocity and angular momentum
240			index = 0;
241			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
242			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
243			integrableObject = mol->nextIntegrableObject(j)) {
244
245			oldVel[index] = integrableObject->getVel();
246	tim	1822	oldJi[index] = integrableObject->getJ();
247	tim	1774	++index;
248			}
249			}
250
251			// do the iteration:
252	tim	1822	instaVol =thermo.getVolume();
253	tim	1774
254			for(int k = 0; k < maxIterNum_; k++) {
255	tim	1822	instaTemp =thermo.getTemperature();
256			instaPress =thermo.getPressure();
257	tim	1774
258			// evolve chi another half step using the temperature at t + dt/2
259	tim	1822	prevChi = chi;
260			chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
261
262			//evolve eta
263	tim	1774	this->evolveEtaB();
264			this->calcVelScale();
265
266			index = 0;
267			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
268			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
269			integrableObject = mol->nextIntegrableObject(j)) {
270	tim	1822
271			frc = integrableObject->getFrc();
272	tim	1774	vel = integrableObject->getVel();
273
274			mass = integrableObject->getMass();
275
276	tim	1822	getVelScaleB(sc, index);
277	tim	1774
278			// velocity half step
279	tim	1822	//vel[j] = oldVel[3 * i + j] + dt2 ((frc[j] / mass) OOPSEConstant::energyConvert - sc[j]);
280			vel = oldVel[index] + dt2OOPSEConstant::energyConvert/mass frc - dt2*sc;
281	tim	1774	integrableObject->setVel(vel);
282
283			if (integrableObject->isDirectional()) {
284			// get and convert the torque to body frame
285			Tb = integrableObject->getTrq();
286			integrableObject->lab2Body(Tb);
287
288	tim	1822	//ji[j] = oldJi[3i + j] + dt2 (Tb[j] * OOPSEConstant::energyConvert - oldJi[3i+j]chi);
289			ji = oldJi[index] + dt2OOPSEConstant::energyConvertTb - dt2chioldJi[index];
290	tim	1774	integrableObject->setJ(ji);
291			}
292
293			++index;
294			}
295			}
296
297			//constraintAlgorithm->doConstrainB();
298
299	tim	1822	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
300	tim	1774	break;
301			}
302
303			//calculate integral of chidt
304			integralOfChidt += dt2 * chi;
305
306	tim	1822	currSnapshot->setChi(chi);
307			currSnapshot->setIntegralOfChiDt(integralOfChidt);
308	tim	1774
309	tim	1822	saveEta();
310	tim	1774	}
311
312			}