src/integrators/NVT.cpp

#inlcude "integrators/NVT.hpp"

namespace oopse {

NVT::NVT(SimInfo* Info) : VelocityVerletIntegrator(info){
    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_.resize(info_->getNIntegrableObjects());
    oldJi_.resize(info_->getNIntegrableObjects());
}

NVT::~NVT() {
}

void NVT::moveA() {
    typename SimInfo::MoleculeIterator i;
    typename Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    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 (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();
        frc = integrableObject->getFrc();

        mass = integrableObject->getMass();

        // velocity half step  (use chi from previous step here):
        //vel[j] += halfStep * ((frc[j] / mass ) * eConvert - vel[j]*chi);
        vel = halfStep *eConvert/mass*frc - halfStep*chi*vel;
        
        // position whole step
        //pos[j] += dt * vel[j];
        pos += fullStep * vel;

        integrableObject->setVel(vel);
        integrableObject->setPos(pos);

        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] += halfStep * (Tb[j] * eConvert - ji[j]*chi);
            ji += halfStep*eConvert*Tb - halfStep*ch *ji;
            this->rotationPropagation(integrableObject, ji);

            integrableObject->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

    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    double chi = currSnapshot->getChi();
    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
    
    chi += halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
    integralOfChidt += chi * halfStep;

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

void NVT::moveB() {
    typename SimInfo::MoleculeIterator i;
    typename Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    
    Vector3d Tb;
    Vector3d ji;    
    Vector3d vel;
    Vector3d frc;
    double mass;
    double instTemp;
    double oldChi, prevChi;
    int index;
    // Set things up for the iteration:

    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
    double chi = currSnapshot->getChi();
    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
    oldChi = 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)) {
                oldVel_[index] = integrableObject->getVel();
                oldJi_[index] = integrableObject->getJ();                
        }
        ++index;               
    }

    // do the iteration:

    for(int k = 0; k < maxIterNum_; k++) {
        index = 0;
        instTemp = tStats->getTemperature();

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

        prevChi = chi;
        chi = oldChi + halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);

        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();

                // velocity half step
                //for(j = 0; j < 3; j++)
                //    vel[j] = oldVel_[3*i+j] + halfStep * ((frc[j] / mass ) * eConvert - oldVel_[3*i + j]*chi);
                vel = oldVel_ + halfStep/mass*eConvert * frc - halfStep*chi*oldVel_[index];
            
                integrableObject->setVel(vel);

                if (integrableObject->isDirectional()) {

                    // get and convert the torque to body frame

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

                    //for(j = 0; j < 3; j++)
                    //    ji[j] = oldJi_[3*i + j] + halfStep * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
                    ji += halfStep*eConvert*Tb - halfStep*ch *oldJi_[index];

                    integrableObject->setJ(ji);
                }
            }
        }
    

        if (nConstrained)
            constrainB();

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

        ++index;
    }

    integralOfChidt += halfStep * chi;

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

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 = info_->getNdf() *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;
}


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