src/integrators/NVT.cpp

#include "integrators/IntegratorCreator.hpp"
#include "integrators/NVT.hpp"
#include "primitives/Molecule.hpp"
#include "utils/simError.h"
#include "utils/OOPSEConstant.hpp"

namespace oopse {

NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {

    Globals* simParams = info_->getSimParams();

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

    update();
}

void NVT::doUpdate() {
    oldVel_.resize(info_->getNIntegrableObjects());
    oldJi_.resize(info_->getNIntegrableObjects());    
}
void NVT::moveA() {
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    Vector3d Tb;
    Vector3d ji;
    double mass;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;

    double chi = currentSnapshot_->getChi();
    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
    
    // We need the temperature at time = t for the chi update below:

    double instTemp = thermo.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] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
        vel += dt2 *OOPSEConstant::energyConvert/mass*frc - dt2*chi*vel;
        
        // position whole step
        //pos[j] += dt * vel[j];
        pos += dt * vel;

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

        if (integrableObject->isDirectional()) {

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

    }
    
    rattle->constraintA();

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

    
    chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
    integralOfChidt += chi * dt2;

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

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

    double chi = currentSnapshot_->getChi();
    double oldChi = chi;
    double  prevChi;
    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();

    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 = thermo.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 (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] + dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - oldVel_[3*i + j]*chi);
                vel = oldVel_[index] + dt2/mass*OOPSEConstant::energyConvert * frc - dt2*chi*oldVel_[index];
            
                integrableObject->setVel(vel);

                if (integrableObject->isDirectional()) {

                    // get and convert the torque to body frame

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

                    //for(j = 0; j < 3; j++)
                    //    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_)
            break;

    }

    integralOfChidt += dt2 * chi;

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


double NVT::calcConservedQuantity() {

    double chi = currentSnapshot_->getChi();
    double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
    double conservedQuantity;
    double fkBT;
    double Energy;
    double thermostat_kinetic;
    double thermostat_potential;
    
    fkBT = info_->getNdf() *OOPSEConstant::kB *targetTemp_;

    Energy = thermo.getTotalE();

    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);

    thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;

    conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;

    return conservedQuantity;
}


}//end namespace oopse
Revision:	1901
Committed:	Tue Jan 4 22:18:36 2005 UTC (19 years, 6 months ago) by tim
File size:	7231 byte(s)
Log Message:	constraints in progress
#	Content
1	#include "integrators/IntegratorCreator.hpp"
2	#include "integrators/NVT.hpp"
3	#include "primitives/Molecule.hpp"
4	#include "utils/simError.h"
5	#include "utils/OOPSEConstant.hpp"
6
7	namespace oopse {
8
9	NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
10
11	Globals* simParams = info_->getSimParams();
12
13	if (simParams->getUseInitXSstate()) {
14	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
15	currSnapshot->setChi(0.0);
16	currSnapshot->setIntegralOfChiDt(0.0);
17	}
18
19	if (!simParams->haveTargetTemp()) {
20	sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
21	painCave.isFatal = 1;
22	painCave.severity = OOPSE_ERROR;
23	simError();
24	} else {
25	targetTemp_ = simParams->getTargetTemp();
26	}
27
28	// We must set tauThermostat_.
29
30	if (!simParams->haveTauThermostat()) {
31	sprintf(painCave.errMsg, "If you use the constant temperature\n"
32	"\tintegrator, you must set tauThermostat_.\n");
33
34	painCave.severity = OOPSE_ERROR;
35	painCave.isFatal = 1;
36	simError();
37	} else {
38	tauThermostat_ = simParams->getTauThermostat();
39	}
40
41	update();
42	}
43
44	void NVT::doUpdate() {
45	oldVel_.resize(info_->getNIntegrableObjects());
46	oldJi_.resize(info_->getNIntegrableObjects());
47	}
48	void NVT::moveA() {
49	SimInfo::MoleculeIterator i;
50	Molecule::IntegrableObjectIterator j;
51	Molecule* mol;
52	StuntDouble* integrableObject;
53	Vector3d Tb;
54	Vector3d ji;
55	double mass;
56	Vector3d vel;
57	Vector3d pos;
58	Vector3d frc;
59
60	double chi = currentSnapshot_->getChi();
61	double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
62
63	// We need the temperature at time = t for the chi update below:
64
65	double instTemp = thermo.getTemperature();
66
67	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
68	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
69	integrableObject = mol->nextIntegrableObject(j)) {
70
71	vel = integrableObject->getVel();
72	pos = integrableObject->getPos();
73	frc = integrableObject->getFrc();
74
75	mass = integrableObject->getMass();
76
77	// velocity half step (use chi from previous step here):
78	//vel[j] += dt2 * ((frc[j] / mass ) * OOPSEConstant::energyConvert - vel[j]*chi);
79	vel += dt2 OOPSEConstant::energyConvert/massfrc - dt2chivel;
80
81	// position whole step
82	//pos[j] += dt * vel[j];
83	pos += dt * vel;
84
85	integrableObject->setVel(vel);
86	integrableObject->setPos(pos);
87
88	if (integrableObject->isDirectional()) {
89
90	//convert the torque to body frame
91	Tb = integrableObject->lab2Body(integrableObject->getTrq());
92
93	// get the angular momentum, and propagate a half step
94
95	ji = integrableObject->getJ();
96
97	//ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
98	ji += dt2OOPSEConstant::energyConvertTb - dt2chi ji;
99	rotAlgo->rotate(integrableObject, ji, dt);
100
101	integrableObject->setJ(ji);
102	}
103	}
104
105	}
106
107	rattle->constraintA();
108
109	// Finally, evolve chi a half step (just like a velocity) using
110	// temperature at time t, not time t+dt/2
111
112
113	chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
114	integralOfChidt += chi * dt2;
115
116	currentSnapshot_->setChi(chi);
117	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
118	}
119
120	void NVT::moveB() {
121	SimInfo::MoleculeIterator i;
122	Molecule::IntegrableObjectIterator j;
123	Molecule* mol;
124	StuntDouble* integrableObject;
125
126	Vector3d Tb;
127	Vector3d ji;
128	Vector3d vel;
129	Vector3d frc;
130	double mass;
131	double instTemp;
132	int index;
133	// Set things up for the iteration:
134
135	double chi = currentSnapshot_->getChi();
136	double oldChi = chi;
137	double prevChi;
138	double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
139
140	index = 0;
141	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
142	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
143	integrableObject = mol->nextIntegrableObject(j)) {
144	oldVel_[index] = integrableObject->getVel();
145	oldJi_[index] = integrableObject->getJ();
146
147	++index;
148	}
149
150	}
151
152	// do the iteration:
153
154	for(int k = 0; k < maxIterNum_; k++) {
155	index = 0;
156	instTemp = thermo.getTemperature();
157
158	// evolve chi another half step using the temperature at t + dt/2
159
160	prevChi = chi;
161	chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
162
163	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
164	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
165	integrableObject = mol->nextIntegrableObject(j)) {
166
167	frc = integrableObject->getFrc();
168	vel = integrableObject->getVel();
169
170	mass = integrableObject->getMass();
171
172	// velocity half step
173	//for(j = 0; j < 3; j++)
174	// vel[j] = oldVel_[3i+j] + dt2 ((frc[j] / mass ) * OOPSEConstant::energyConvert - oldVel_[3i + j]chi);
175	vel = oldVel_[index] + dt2/massOOPSEConstant::energyConvert frc - dt2chioldVel_[index];
176
177	integrableObject->setVel(vel);
178
179	if (integrableObject->isDirectional()) {
180
181	// get and convert the torque to body frame
182
183	Tb = integrableObject->lab2Body(integrableObject->getTrq());
184
185	//for(j = 0; j < 3; j++)
186	// ji[j] = oldJi_[3i + j] + dt2 (Tb[j] * OOPSEConstant::energyConvert - oldJi_[3i+j]chi);
187	ji = oldJi_[index] + dt2OOPSEConstant::energyConvertTb - dt2chi oldJi_[index];
188
189	integrableObject->setJ(ji);
190	}
191
192
193	++index;
194	}
195	}
196
197
198	rattle->constraintB();
199
200	if (fabs(prevChi - chi) <= chiTolerance_)
201	break;
202
203	}
204
205	integralOfChidt += dt2 * chi;
206
207	currentSnapshot_->setChi(chi);
208	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
209	}
210
211
212	double NVT::calcConservedQuantity() {
213
214	double chi = currentSnapshot_->getChi();
215	double integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
216	double conservedQuantity;
217	double fkBT;
218	double Energy;
219	double thermostat_kinetic;
220	double thermostat_potential;
221
222	fkBT = info_->getNdf() OOPSEConstant::kB targetTemp_;
223
224	Energy = thermo.getTotalE();
225
226	thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert);
227
228	thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert;
229
230	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
231
232	return conservedQuantity;
233	}
234
235
236	}//end namespace oopse