src/integrators/NVT.cpp

#inlcude "integrators/NVT.hpp"

namespace oopse {

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

    Globals* globals = info_->getGlobals();

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

    // We must set tauThermostat_.

    if (!globals->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_ = globals->getTauThermostat();
    }

    update();
}

NVT::~NVT() {
}

void NVT::update() {
    oldVel_.resize(info_->getNIntegrableObjects());
    oldJi_.resize(info_->getNIntegrableObjects());    
}
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] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
        vel = dt2 *eConvert/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()) {

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

            integrableObject->setJ(ji);
        }
    }

    }
    
    //constraintAlgorithm->doConstrainA();

    // 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 += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
    integralOfChidt += chi * dt2;

    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 + 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 ) * eConvert - oldVel_[3*i + j]*chi);
                vel = oldVel_ + dt2/mass*eConvert * frc - dt2*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] + dt2 * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
                    ji += dt2*eConvert*Tb - dt2*ch *oldJi_[index];

                    integrableObject->setJ(ji);
                }
            }
        }
    

        //constraintAlgorithm->doConstrainB();

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

        ++index;
    }

    integralOfChidt += dt2 * chi;

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


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