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* 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() {
    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:	1865
Committed:	Tue Dec 7 05:12:41 2004 UTC (19 years, 6 months ago) by tim
File size:	9827 byte(s)
Log Message:	fix a bug in UseInitXSstate
#	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) {
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	VelocityVerletIntegrator::update();
87	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	NkBT = info_->getNGlobalIntegrableObjects()OOPSEConstant::kB targetTemp;
94
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	fkBT = info_->getNdf()OOPSEConstant::kB targetTemp;
99	}
100
101	void NPT::moveA() {
102	SimInfo::MoleculeIterator i;
103	Molecule::IntegrableObjectIterator j;
104	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
114	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
115	chi= currSnapshot->getChi();
116	integralOfChidt = currSnapshot->getIntegralOfChiDt();
117	loadEta();
118
119	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
124	Vector3d COM = info_->getCom();
125
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	//vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
143	vel += dt2OOPSEConstant::energyConvert/mass frc - dt2*sc;
144	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	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
158	ji += dt2OOPSEConstant::energyConvert Tb - dt2chi ji;
159
160	rotAlgo->rotate(integrableObject, ji, dt);
161
162	integrableObject->setJ(ji);
163	}
164
165	}
166	}
167	// evolve chi and eta half step
168
169	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
170
171	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
211	currSnapshot->setChi(chi);
212	currSnapshot->setIntegralOfChiDt(integralOfChidt);
213
214	saveEta();
215	}
216
217	void NPT::moveB(void) {
218	SimInfo::MoleculeIterator i;
219	Molecule::IntegrableObjectIterator j;
220	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
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	//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	oldJi[index] = integrableObject->getJ();
247	++index;
248	}
249	}
250
251	// do the iteration:
252	instaVol =thermo.getVolume();
253
254	for(int k = 0; k < maxIterNum_; k++) {
255	instaTemp =thermo.getTemperature();
256	instaPress =thermo.getPressure();
257
258	// evolve chi another half step using the temperature at t + dt/2
259	prevChi = chi;
260	chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
261
262	//evolve eta
263	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
271	frc = integrableObject->getFrc();
272	vel = integrableObject->getVel();
273
274	mass = integrableObject->getMass();
275
276	getVelScaleB(sc, index);
277
278	// velocity half step
279	//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	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	//ji[j] = oldJi[3i + j] + dt2 (Tb[j] * OOPSEConstant::energyConvert - oldJi[3i+j]chi);
289	ji = oldJi[index] + dt2OOPSEConstant::energyConvertTb - dt2chioldJi[index];
290	integrableObject->setJ(ji);
291	}
292
293	++index;
294	}
295	}
296
297	//constraintAlgorithm->doConstrainB();
298
299	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
300	break;
301	}
302
303	//calculate integral of chidt
304	integralOfChidt += dt2 * chi;
305
306	currSnapshot->setChi(chi);
307	currSnapshot->setIntegralOfChiDt(integralOfChidt);
308
309	saveEta();
310	}
311
312	}