src/integrators/NVT.cpp

#include <math.h>

#include "primitives/Atom.hpp"
#include "primitives/SRI.hpp"
#include "primitives/AbstractClasses.hpp"
#include "brains/SimInfo.hpp"
#include "UseTheForce/ForceFields.hpp"
#include "brains/Thermo.hpp"
#include "io/ReadWrite.hpp"
#include "integrators/Integrator.hpp"
#include "utils/simError.h"

// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697

template <typename T>NVT<T>::NVT(SimInfo *theInfo, ForceFields *the_ff) :
    T(theInfo, the_ff) {
    GenericData *data;
    DoubleGenericData *chiValue;
    DoubleGenericData *integralOfChidtValue;

    chiValue = NULL;
    integralOfChidtValue = NULL;

    chi = 0.0;
    have_tau_thermostat = 0;
    have_target_temp = 0;
    have_chi_tolerance = 0;
    integralOfChidt = 0.0;

    if (theInfo->useInitXSstate) {

        // retrieve chi and integralOfChidt from simInfo
        data = info->getPropertyByName(CHIVALUE_ID);

        if (data) {
            chiValue = dynamic_cast<DoubleGenericData *>(data);
        }

        data = info->getPropertyByName(INTEGRALOFCHIDT_ID);

        if (data) {
            integralOfChidtValue = dynamic_cast<DoubleGenericData *>(data);
        }

        // chi and integralOfChidt should appear by pair
        if (chiValue && integralOfChidtValue) {
            chi = chiValue->getData();
            integralOfChidt = integralOfChidtValue->getData();
        }
    }

    oldVel = new double[3 * integrableObjects.size()];
    oldJi = new double[3 * integrableObjects.size()];
}

template <typename T>NVT<T>::~NVT() {
    delete [] oldVel;
    delete [] oldJi;
}

template <typename T>
void NVT<T>::moveA() {
    int i, j;
    DirectionalAtom *dAtom;
    Vector3d Tb;
    Vector3d ji;
    double mass;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;

    double instTemp;

    // We need the temperature at time = t for the chi update below:

    instTemp = tStats->getTemperature();

    for(i = 0; i < integrableObjects.size(); i++) {
        vel = integrableObjects[i]->getVel();
        pos = integrableObjects[i]->getPos();
        integrableObjects[i]->getFrc(frc);

        mass = integrableObjects[i]->getMass();

        for(j = 0; j < 3; j++) {
            // velocity half step  (use chi from previous step here):
            vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);

            // position whole step
            pos[j] += dt * vel[j];
        }

        integrableObjects[i]->setVel(vel);
        integrableObjects[i]->setPos(pos);

        if (integrableObjects[i]->isDirectional()) {

            // get and convert the torque to body frame

            Tb = integrableObjects[i]->getTrq();
            integrableObjects[i]->lab2Body(Tb);

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

            ji = integrableObjects[i]->getJ();

            for(j = 0; j < 3; j++)
        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);

            this->rotationPropagation(integrableObjects[i], ji);

            integrableObjects[i]->setJ(ji);
        }
    }

    if (nConstrained)
        constrainA();

    // Finally, evolve chi a half step (just like a velocity) using
    // temperature at time t, not time t+dt/2

    //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";

    chi += dt2 * (instTemp / targetTemp - 1.0) / (tauThermostat * tauThermostat);
    integralOfChidt += chi * dt2;
}

template <typename T>
void NVT<T>::moveB(void) {
    int i, j, k;
    double Tb[3], ji[3];
    double vel[3], frc[3];
    double mass;
    double instTemp;
    double oldChi, prevChi;

    // Set things up for the iteration:

    oldChi = chi;

    for(i = 0; i < integrableObjects.size(); i++) {
        vel = integrableObjects[i]->getVel();

        for(j = 0; j < 3; j++)
    oldVel[3 * i + j] = vel[j];

        if (integrableObjects[i]->isDirectional()) {
            ji = integrableObjects[i]->getJ();

            for(j = 0; j < 3; j++)
        oldJi[3 * i + j] = ji[j];
        }
    }

    // do the iteration:

    for(k = 0; k < 4; k++) {
        instTemp = tStats->getTemperature();

        // evolve chi another half step using the temperature at t + dt/2

        prevChi = chi;
        chi = oldChi + dt2 * (instTemp / targetTemp - 1.0) / (tauThermostat * tauThermostat);

        for(i = 0; i < integrableObjects.size(); i++) {
            integrableObjects[i]->getFrc(frc);
            vel = integrableObjects[i]->getVel();

            mass = integrableObjects[i]->getMass();

            // velocity half step
            for(j = 0; j < 3; j++)
        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*chi);

            integrableObjects[i]->setVel(vel);

            if (integrableObjects[i]->isDirectional()) {

                // get and convert the torque to body frame

                Tb = integrableObjects[i]->getTrq();
                integrableObjects[i]->lab2Body(Tb);

                for(j = 0; j < 3; j++)
            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);

                integrableObjects[i]->setJ(ji);
            }
        }

        if (nConstrained)
            constrainB();

        if (fabs(prevChi - chi) <= chiTolerance)
            break;
    }

    integralOfChidt += dt2 * chi;
}

template <typename T>
void NVT<T>::resetIntegrator(void) {
    chi = 0.0;
    integralOfChidt = 0.0;
}

template <typename T>
int NVT<T>::readyCheck() {

    //check parent's readyCheck() first
    if (T::readyCheck() == -1)
        return -1;

    // First check to see if we have a target temperature.
    // Not having one is fatal.

    if (!have_target_temp) {
        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
        painCave.isFatal = 1;
        painCave.severity = OOPSE_ERROR;
        simError();
        return -1;
    }

    // We must set tauThermostat.

    if (!have_tau_thermostat) {
        sprintf(painCave.errMsg, "If you use the constant temperature\n"
                                     "\tintegrator, you must set tauThermostat.\n");

        painCave.severity = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();
        return -1;
    }

    if (!have_chi_tolerance) {
        sprintf(painCave.errMsg, "In NVT integrator: setting chi tolerance to 1e-6\n");
        chiTolerance = 1e - 6;
        have_chi_tolerance = 1;
        painCave.severity = OOPSE_INFO;
        painCave.isFatal = 0;
        simError();
    }

    return 1;
}

template <typename T>
double NVT<T>::getConservedQuantity(void) {
    double conservedQuantity;
    double fkBT;
    double Energy;
    double thermostat_kinetic;
    double thermostat_potential;

    fkBT = (double)(info->ndf) *kB *targetTemp;

    Energy = tStats->getTotalE();

    thermostat_kinetic = fkBT * tauThermostat * tauThermostat * chi * chi / (2.0 * eConvert);

    thermostat_potential = fkBT * integralOfChidt / eConvert;

    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;

    return conservedQuantity;
}

template <typename T>
string NVT<T>::getAdditionalParameters(void) {
    string parameters;
    const int BUFFERSIZE = 2000; // size of the read buffer
    char buffer[BUFFERSIZE];

    sprintf(buffer, "\t%G\t%G;", chi, integralOfChidt);
    parameters += buffer;

    return parameters;
}
Revision:	1762
Committed:	Fri Nov 19 21:38:22 2004 UTC (19 years, 9 months ago) by tim
File size:	7452 byte(s)
Log Message:	refactory integrator
#	User	Rev	Content
1	tim	1762	#include <math.h>
2
3			#include "primitives/Atom.hpp"
4			#include "primitives/SRI.hpp"
5			#include "primitives/AbstractClasses.hpp"
6			#include "brains/SimInfo.hpp"
7			#include "UseTheForce/ForceFields.hpp"
8			#include "brains/Thermo.hpp"
9			#include "io/ReadWrite.hpp"
10			#include "integrators/Integrator.hpp"
11			#include "utils/simError.h"
12
13			// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
14
15			template <typename T>NVT<T>::NVT(SimInfo theInfo, ForceFields the_ff) :
16			T(theInfo, the_ff) {
17			GenericData *data;
18			DoubleGenericData *chiValue;
19			DoubleGenericData *integralOfChidtValue;
20
21			chiValue = NULL;
22			integralOfChidtValue = NULL;
23
24			chi = 0.0;
25			have_tau_thermostat = 0;
26			have_target_temp = 0;
27			have_chi_tolerance = 0;
28			integralOfChidt = 0.0;
29
30			if (theInfo->useInitXSstate) {
31
32			// retrieve chi and integralOfChidt from simInfo
33			data = info->getPropertyByName(CHIVALUE_ID);
34
35			if (data) {
36			chiValue = dynamic_cast<DoubleGenericData *>(data);
37			}
38
39			data = info->getPropertyByName(INTEGRALOFCHIDT_ID);
40
41			if (data) {
42			integralOfChidtValue = dynamic_cast<DoubleGenericData *>(data);
43			}
44
45			// chi and integralOfChidt should appear by pair
46			if (chiValue && integralOfChidtValue) {
47			chi = chiValue->getData();
48			integralOfChidt = integralOfChidtValue->getData();
49			}
50			}
51
52			oldVel = new double[3 * integrableObjects.size()];
53			oldJi = new double[3 * integrableObjects.size()];
54			}
55
56			template <typename T>NVT<T>::~NVT() {
57			delete [] oldVel;
58			delete [] oldJi;
59			}
60
61			template <typename T>
62			void NVT<T>::moveA() {
63			int i, j;
64			DirectionalAtom *dAtom;
65			Vector3d Tb;
66			Vector3d ji;
67			double mass;
68			Vector3d vel;
69			Vector3d pos;
70			Vector3d frc;
71
72			double instTemp;
73
74			// We need the temperature at time = t for the chi update below:
75
76			instTemp = tStats->getTemperature();
77
78			for(i = 0; i < integrableObjects.size(); i++) {
79			vel = integrableObjects[i]->getVel();
80			pos = integrableObjects[i]->getPos();
81			integrableObjects[i]->getFrc(frc);
82
83			mass = integrableObjects[i]->getMass();
84
85			for(j = 0; j < 3; j++) {
86			// velocity half step (use chi from previous step here):
87			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
88
89			// position whole step
90			pos[j] += dt * vel[j];
91			}
92
93			integrableObjects[i]->setVel(vel);
94			integrableObjects[i]->setPos(pos);
95
96			if (integrableObjects[i]->isDirectional()) {
97
98			// get and convert the torque to body frame
99
100			Tb = integrableObjects[i]->getTrq();
101			integrableObjects[i]->lab2Body(Tb);
102
103			// get the angular momentum, and propagate a half step
104
105			ji = integrableObjects[i]->getJ();
106
107			for(j = 0; j < 3; j++)
108			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
109
110			this->rotationPropagation(integrableObjects[i], ji);
111
112			integrableObjects[i]->setJ(ji);
113			}
114			}
115
116			if (nConstrained)
117			constrainA();
118
119			// Finally, evolve chi a half step (just like a velocity) using
120			// temperature at time t, not time t+dt/2
121
122			//std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
123
124			chi += dt2 * (instTemp / targetTemp - 1.0) / (tauThermostat * tauThermostat);
125			integralOfChidt += chi * dt2;
126			}
127
128			template <typename T>
129			void NVT<T>::moveB(void) {
130			int i, j, k;
131			double Tb[3], ji[3];
132			double vel[3], frc[3];
133			double mass;
134			double instTemp;
135			double oldChi, prevChi;
136
137			// Set things up for the iteration:
138
139			oldChi = chi;
140
141			for(i = 0; i < integrableObjects.size(); i++) {
142			vel = integrableObjects[i]->getVel();
143
144			for(j = 0; j < 3; j++)
145			oldVel[3 * i + j] = vel[j];
146
147			if (integrableObjects[i]->isDirectional()) {
148			ji = integrableObjects[i]->getJ();
149
150			for(j = 0; j < 3; j++)
151			oldJi[3 * i + j] = ji[j];
152			}
153			}
154
155			// do the iteration:
156
157			for(k = 0; k < 4; k++) {
158			instTemp = tStats->getTemperature();
159
160			// evolve chi another half step using the temperature at t + dt/2
161
162			prevChi = chi;
163			chi = oldChi + dt2 * (instTemp / targetTemp - 1.0) / (tauThermostat * tauThermostat);
164
165			for(i = 0; i < integrableObjects.size(); i++) {
166			integrableObjects[i]->getFrc(frc);
167			vel = integrableObjects[i]->getVel();
168
169			mass = integrableObjects[i]->getMass();
170
171			// velocity half step
172			for(j = 0; j < 3; j++)
173			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j]chi);
174
175			integrableObjects[i]->setVel(vel);
176
177			if (integrableObjects[i]->isDirectional()) {
178
179			// get and convert the torque to body frame
180
181			Tb = integrableObjects[i]->getTrq();
182			integrableObjects[i]->lab2Body(Tb);
183
184			for(j = 0; j < 3; j++)
185			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
186
187			integrableObjects[i]->setJ(ji);
188			}
189			}
190
191			if (nConstrained)
192			constrainB();
193
194			if (fabs(prevChi - chi) <= chiTolerance)
195			break;
196			}
197
198			integralOfChidt += dt2 * chi;
199			}
200
201			template <typename T>
202			void NVT<T>::resetIntegrator(void) {
203			chi = 0.0;
204			integralOfChidt = 0.0;
205			}
206
207			template <typename T>
208			int NVT<T>::readyCheck() {
209
210			//check parent's readyCheck() first
211			if (T::readyCheck() == -1)
212			return -1;
213
214			// First check to see if we have a target temperature.
215			// Not having one is fatal.
216
217			if (!have_target_temp) {
218			sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
219			painCave.isFatal = 1;
220			painCave.severity = OOPSE_ERROR;
221			simError();
222			return -1;
223			}
224
225			// We must set tauThermostat.
226
227			if (!have_tau_thermostat) {
228			sprintf(painCave.errMsg, "If you use the constant temperature\n"
229			"\tintegrator, you must set tauThermostat.\n");
230
231			painCave.severity = OOPSE_ERROR;
232			painCave.isFatal = 1;
233			simError();
234			return -1;
235			}
236
237			if (!have_chi_tolerance) {
238			sprintf(painCave.errMsg, "In NVT integrator: setting chi tolerance to 1e-6\n");
239			chiTolerance = 1e - 6;
240			have_chi_tolerance = 1;
241			painCave.severity = OOPSE_INFO;
242			painCave.isFatal = 0;
243			simError();
244			}
245
246			return 1;
247			}
248
249			template <typename T>
250			double NVT<T>::getConservedQuantity(void) {
251			double conservedQuantity;
252			double fkBT;
253			double Energy;
254			double thermostat_kinetic;
255			double thermostat_potential;
256
257			fkBT = (double)(info->ndf) kB targetTemp;
258
259			Energy = tStats->getTotalE();
260
261			thermostat_kinetic = fkBT * tauThermostat * tauThermostat * chi * chi / (2.0 * eConvert);
262
263			thermostat_potential = fkBT * integralOfChidt / eConvert;
264
265			conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
266
267			return conservedQuantity;
268			}
269
270			template <typename T>
271			string NVT<T>::getAdditionalParameters(void) {
272			string parameters;
273			const int BUFFERSIZE = 2000; // size of the read buffer
274			char buffer[BUFFERSIZE];
275
276			sprintf(buffer, "\t%G\t%G;", chi, integralOfChidt);
277			parameters += buffer;
278
279			return parameters;
280			}