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), maxIterNum_(4) {

    Globals* simParams = info_->getSimParams();
    
    if (!simParams->getUseInitXSstate()) {
        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
        currSnapshot->setChi(0.0);
        currSnapshot->setIntegralOfChiDt(0.0);
        currSnapshot->setEta(Mat3x3d(0.0));
    }
    
    if (!simParams->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 = simParams->getTargetTemp();
    }

    // We must set tauThermostat
    if (!simParams->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 = simParams->getTauThermostat();
    }

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

        painCave.isFatal = 1;
        simError();
    } else {
        targetPressure = simParams->getTargetPressure();
    }
    
    if (!simParams->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 = simParams->getTauBarostat();
    }
    
    tt2 = tauThermostat * tauThermostat;
    tb2 = tauBarostat * tauBarostat;

    update();
}

NPT::~NPT() {
}

void NPT::doUpdate() {

    oldPos.resize(info_->getNIntegrableObjects());
    oldVel.resize(info_->getNIntegrableObjects());
    oldJi.resize(info_->getNIntegrableObjects());

}

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;

    chi= currentSnapshot_->getChi();
    integralOfChidt = currentSnapshot_->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->lab2Body(integrableObject->getTrq());

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

        rattle->constraintA();
    }

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

    this->scaleSimBox();

    currentSnapshot_->setChi(chi);
    currentSnapshot_->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;


    chi= currentSnapshot_->getChi();
    integralOfChidt = currentSnapshot_->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->lab2Body(integrableObject->getTrq());

                    //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;
            }
        }
        
        rattle->constraintB();

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

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

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

    saveEta();
}

}
Revision:	1901
Committed:	Tue Jan 4 22:18:36 2005 UTC (19 years, 6 months ago) by tim
File size:	9091 byte(s)
Log Message:	constraints in progress
#	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	tim	1883	VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
25	tim	1774
26	tim	1841	Globals* simParams = info_->getSimParams();
27	tim	1774
28	tim	1865	if (!simParams->getUseInitXSstate()) {
29	tim	1774	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	tim	1841	if (!simParams->haveTargetTemp()) {
36	tim	1774	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	tim	1841	targetTemp = simParams->getTargetTemp();
42	tim	1774	}
43
44			// We must set tauThermostat
45	tim	1841	if (!simParams->haveTauThermostat()) {
46	tim	1774	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	tim	1841	tauThermostat = simParams->getTauThermostat();
54	tim	1774	}
55
56	tim	1841	if (!simParams->haveTargetPressure()) {
57	tim	1774	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	tim	1841	targetPressure = simParams->getTargetPressure();
64	tim	1774	}
65
66	tim	1841	if (!simParams->haveTauBarostat()) {
67	tim	1774	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	tim	1841	tauBarostat = simParams->getTauBarostat();
74	tim	1774	}
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	tim	1867
87	tim	1774	oldPos.resize(info_->getNIntegrableObjects());
88			oldVel.resize(info_->getNIntegrableObjects());
89			oldJi.resize(info_->getNIntegrableObjects());
90
91			}
92
93			void NPT::moveA() {
94	tim	1822	SimInfo::MoleculeIterator i;
95			Molecule::IntegrableObjectIterator j;
96	tim	1774	Molecule* mol;
97			StuntDouble* integrableObject;
98			Vector3d Tb, ji;
99			double mass;
100			Vector3d vel;
101			Vector3d pos;
102			Vector3d frc;
103			Vector3d sc;
104			int index;
105	tim	1822
106	tim	1867	chi= currentSnapshot_->getChi();
107			integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
108	tim	1822	loadEta();
109	tim	1774
110	tim	1822	instaTemp =thermo.getTemperature();
111			press = thermo.getPressureTensor();
112			instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
113			instaVol =thermo.getVolume();
114	tim	1774
115	tim	1822	Vector3d COM = info_->getCom();
116	tim	1774
117			//evolve velocity half step
118
119			calcVelScale();
120
121			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
122			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
123			integrableObject = mol->nextIntegrableObject(j)) {
124
125			vel = integrableObject->getVel();
126			frc = integrableObject->getFrc();
127
128			mass = integrableObject->getMass();
129
130			getVelScaleA(sc, vel);
131
132			// velocity half step (use chi from previous step here):
133	tim	1822	//vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
134			vel += dt2OOPSEConstant::energyConvert/mass frc - dt2*sc;
135	tim	1774	integrableObject->setVel(vel);
136
137			if (integrableObject->isDirectional()) {
138
139			// get and convert the torque to body frame
140
141	tim	1871	Tb = integrableObject->lab2Body(integrableObject->getTrq());
142	tim	1774
143			// get the angular momentum, and propagate a half step
144
145			ji = integrableObject->getJ();
146
147	tim	1822	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
148			ji += dt2OOPSEConstant::energyConvert Tb - dt2chi ji;
149	tim	1774
150	tim	1822	rotAlgo->rotate(integrableObject, ji, dt);
151	tim	1774
152			integrableObject->setJ(ji);
153			}
154
155			}
156			}
157			// evolve chi and eta half step
158
159	tim	1822	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
160
161	tim	1774	evolveEtaA();
162
163			//calculate the integral of chidt
164			integralOfChidt += dt2 * chi;
165
166			index = 0;
167			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
168			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
169			integrableObject = mol->nextIntegrableObject(j)) {
170			oldPos[index++] = integrableObject->getPos();
171			}
172			}
173
174			//the first estimation of r(t+dt) is equal to r(t)
175
176			for(int k = 0; k < maxIterNum_; k++) {
177			index = 0;
178			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
179			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
180			integrableObject = mol->nextIntegrableObject(j)) {
181
182			vel = integrableObject->getVel();
183			pos = integrableObject->getPos();
184
185			this->getPosScale(pos, COM, index, sc);
186
187			pos = oldPos[index] + dt * (vel + sc);
188			integrableObject->setPos(pos);
189
190			++index;
191			}
192			}
193
194	tim	1901	rattle->constraintA();
195	tim	1774	}
196
197			// Scale the box after all the positions have been moved:
198
199			this->scaleSimBox();
200	tim	1822
201	tim	1867	currentSnapshot_->setChi(chi);
202			currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
203	tim	1822
204			saveEta();
205	tim	1774	}
206
207			void NPT::moveB(void) {
208	tim	1822	SimInfo::MoleculeIterator i;
209			Molecule::IntegrableObjectIterator j;
210	tim	1774	Molecule* mol;
211			StuntDouble* integrableObject;
212			int index;
213			Vector3d Tb;
214			Vector3d ji;
215			Vector3d sc;
216			Vector3d vel;
217			Vector3d frc;
218			double mass;
219
220	tim	1822
221	tim	1867	chi= currentSnapshot_->getChi();
222			integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
223	tim	1822	double oldChi = chi;
224			double prevChi;
225
226			loadEta();
227
228	tim	1774	//save velocity and angular momentum
229			index = 0;
230			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
231			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
232			integrableObject = mol->nextIntegrableObject(j)) {
233
234			oldVel[index] = integrableObject->getVel();
235	tim	1822	oldJi[index] = integrableObject->getJ();
236	tim	1774	++index;
237			}
238			}
239
240			// do the iteration:
241	tim	1822	instaVol =thermo.getVolume();
242	tim	1774
243			for(int k = 0; k < maxIterNum_; k++) {
244	tim	1822	instaTemp =thermo.getTemperature();
245			instaPress =thermo.getPressure();
246	tim	1774
247			// evolve chi another half step using the temperature at t + dt/2
248	tim	1822	prevChi = chi;
249			chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
250
251			//evolve eta
252	tim	1774	this->evolveEtaB();
253			this->calcVelScale();
254
255			index = 0;
256			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
257			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
258			integrableObject = mol->nextIntegrableObject(j)) {
259	tim	1822
260			frc = integrableObject->getFrc();
261	tim	1774	vel = integrableObject->getVel();
262
263			mass = integrableObject->getMass();
264
265	tim	1822	getVelScaleB(sc, index);
266	tim	1774
267			// velocity half step
268	tim	1822	//vel[j] = oldVel[3 * i + j] + dt2 ((frc[j] / mass) OOPSEConstant::energyConvert - sc[j]);
269			vel = oldVel[index] + dt2OOPSEConstant::energyConvert/mass frc - dt2*sc;
270	tim	1774	integrableObject->setVel(vel);
271
272			if (integrableObject->isDirectional()) {
273			// get and convert the torque to body frame
274	tim	1871	Tb = integrableObject->lab2Body(integrableObject->getTrq());
275	tim	1774
276	tim	1822	//ji[j] = oldJi[3i + j] + dt2 (Tb[j] * OOPSEConstant::energyConvert - oldJi[3i+j]chi);
277			ji = oldJi[index] + dt2OOPSEConstant::energyConvertTb - dt2chioldJi[index];
278	tim	1774	integrableObject->setJ(ji);
279			}
280
281			++index;
282			}
283			}
284
285	tim	1901	rattle->constraintB();
286	tim	1774
287	tim	1822	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
288	tim	1774	break;
289			}
290
291			//calculate integral of chidt
292			integralOfChidt += dt2 * chi;
293
294	tim	1867	currentSnapshot_->setChi(chi);
295			currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
296	tim	1774
297	tim	1822	saveEta();
298	tim	1774	}
299
300			}