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
#	Content
1	#include <math.h>
2
3	#include "brains/SimInfo.hpp"
4	#include "brains/Thermo.hpp"
5	#include "integrators/NPT.hpp"
6	#include "math/SquareMatrix3.hpp"
7	#include "primitives/Molecule.hpp"
8	#include "utils/OOPSEConstant.hpp"
9	#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), maxIterNum_(4) {
25
26	Globals* simParams = info_->getSimParams();
27
28	if (!simParams->getUseInitXSstate()) {
29	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
30	currSnapshot->setChi(0.0);
31	currSnapshot->setIntegralOfChiDt(0.0);
32	currSnapshot->setEta(Mat3x3d(0.0));
33	}
34
35	if (!simParams->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 = simParams->getTargetTemp();
42	}
43
44	// We must set tauThermostat
45	if (!simParams->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 = simParams->getTauThermostat();
54	}
55
56	if (!simParams->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 = simParams->getTargetPressure();
64	}
65
66	if (!simParams->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 = simParams->getTauBarostat();
74	}
75
76	tt2 = tauThermostat * tauThermostat;
77	tb2 = tauBarostat * tauBarostat;
78
79	update();
80	}
81
82	NPT::~NPT() {
83	}
84
85	void NPT::doUpdate() {
86
87	oldPos.resize(info_->getNIntegrableObjects());
88	oldVel.resize(info_->getNIntegrableObjects());
89	oldJi.resize(info_->getNIntegrableObjects());
90
91	}
92
93	void NPT::moveA() {
94	SimInfo::MoleculeIterator i;
95	Molecule::IntegrableObjectIterator j;
96	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
106	chi= currentSnapshot_->getChi();
107	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
108	loadEta();
109
110	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
115	Vector3d COM = info_->getCom();
116
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	//vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
134	vel += dt2OOPSEConstant::energyConvert/mass frc - dt2*sc;
135	integrableObject->setVel(vel);
136
137	if (integrableObject->isDirectional()) {
138
139	// get and convert the torque to body frame
140
141	Tb = integrableObject->lab2Body(integrableObject->getTrq());
142
143	// get the angular momentum, and propagate a half step
144
145	ji = integrableObject->getJ();
146
147	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
148	ji += dt2OOPSEConstant::energyConvert Tb - dt2chi ji;
149
150	rotAlgo->rotate(integrableObject, ji, dt);
151
152	integrableObject->setJ(ji);
153	}
154
155	}
156	}
157	// evolve chi and eta half step
158
159	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
160
161	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	rattle->constraintA();
195	}
196
197	// Scale the box after all the positions have been moved:
198
199	this->scaleSimBox();
200
201	currentSnapshot_->setChi(chi);
202	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
203
204	saveEta();
205	}
206
207	void NPT::moveB(void) {
208	SimInfo::MoleculeIterator i;
209	Molecule::IntegrableObjectIterator j;
210	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
221	chi= currentSnapshot_->getChi();
222	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
223	double oldChi = chi;
224	double prevChi;
225
226	loadEta();
227
228	//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	oldJi[index] = integrableObject->getJ();
236	++index;
237	}
238	}
239
240	// do the iteration:
241	instaVol =thermo.getVolume();
242
243	for(int k = 0; k < maxIterNum_; k++) {
244	instaTemp =thermo.getTemperature();
245	instaPress =thermo.getPressure();
246
247	// evolve chi another half step using the temperature at t + dt/2
248	prevChi = chi;
249	chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
250
251	//evolve eta
252	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
260	frc = integrableObject->getFrc();
261	vel = integrableObject->getVel();
262
263	mass = integrableObject->getMass();
264
265	getVelScaleB(sc, index);
266
267	// velocity half step
268	//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	integrableObject->setVel(vel);
271
272	if (integrableObject->isDirectional()) {
273	// get and convert the torque to body frame
274	Tb = integrableObject->lab2Body(integrableObject->getTrq());
275
276	//ji[j] = oldJi[3i + j] + dt2 (Tb[j] * OOPSEConstant::energyConvert - oldJi[3i+j]chi);
277	ji = oldJi[index] + dt2OOPSEConstant::energyConvertTb - dt2chioldJi[index];
278	integrableObject->setJ(ji);
279	}
280
281	++index;
282	}
283	}
284
285	rattle->constraintB();
286
287	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
288	break;
289	}
290
291	//calculate integral of chidt
292	integralOfChidt += dt2 * chi;
293
294	currentSnapshot_->setChi(chi);
295	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
296
297	saveEta();
298	}
299
300	}