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 (19 years, 7 months ago) by tim
File size:	8103 byte(s)
Log Message:	adding section parsers
#	Content
1	#inlcude "integrators/NVT.hpp"
2
3	namespace oopse {
4
5	NVT::NVT(SimInfo* Info) : VelocityVerletIntegrator(info){
6	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	oldVel_.resize(info_->getNIntegrableObjects());
42	oldJi_.resize(info_->getNIntegrableObjects());
43	}
44
45	NVT::~NVT() {
46	}
47
48	void NVT::moveA() {
49	typename SimInfo::MoleculeIterator i;
50	typename 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 instTemp;
61
62	// We need the temperature at time = t for the chi update below:
63
64	instTemp = tStats->getTemperature();
65
66	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
70	vel = integrableObject->getVel();
71	pos = integrableObject->getPos();
72	frc = integrableObject->getFrc();
73
74	mass = integrableObject->getMass();
75
76	// 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
84	integrableObject->setVel(vel);
85	integrableObject->setPos(pos);
86
87	if (integrableObject->isDirectional()) {
88
89	// get and convert the torque to body frame
90
91	Tb = integrableObject->getTrq();
92	integrableObject->lab2Body(Tb);
93
94	// get the angular momentum, and propagate a half step
95
96	ji = integrableObject->getJ();
97
98	//ji[j] += halfStep * (Tb[j] * eConvert - ji[j]*chi);
99	ji += halfStepeConvertTb - halfStepch ji;
100	this->rotationPropagation(integrableObject, ji);
101
102	integrableObject->setJ(ji);
103	}
104	}
105
106	}
107
108	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	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
121	currSnapshot->setChi(chi);
122	currSnapshot->setIntegralOfChiDt(integralOfChidt);
123	}
124
125	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	double mass;
136	double instTemp;
137	double oldChi, prevChi;
138	int index;
139	// Set things up for the iteration:
140
141	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
142	double chi = currSnapshot->getChi();
143	double integralOfChidt = currSnapshot->getIntegralOfChiDt();
144	oldChi = chi;
145
146	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	}
153	++index;
154	}
155
156	// do the iteration:
157
158	for(int k = 0; k < maxIterNum_; k++) {
159	index = 0;
160	instTemp = tStats->getTemperature();
161
162	// evolve chi another half step using the temperature at t + dt/2
163
164	prevChi = chi;
165	chi = oldChi + halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
166
167	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
171	frc = integrableObject->getFrc();
172	vel = integrableObject->getVel();
173
174	mass = integrableObject->getMass();
175
176	// 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
183	if (integrableObject->isDirectional()) {
184
185	// get and convert the torque to body frame
186
187	Tb = integrableObject->getTrq();
188	integrableObject->lab2Body(Tb);
189
190	//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
194	integrableObject->setJ(ji);
195	}
196	}
197	}
198
199
200	if (nConstrained)
201	constrainB();
202
203	if (fabs(prevChi - chi) <= chiTolerance)
204	break;
205
206	++index;
207	}
208
209	integralOfChidt += halfStep * chi;
210
211	currSnapshot->setChi(chi);
212	currSnapshot->setIntegralOfChiDt(integralOfChidt);
213	}
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	sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
227	painCave.isFatal = 1;
228	painCave.severity = OOPSE_ERROR;
229	simError();
230	return -1;
231	}
232
233	// We must set tauThermostat_.
234
235	if (!have_tau_thermostat) {
236	sprintf(painCave.errMsg, "If you use the constant temperature\n"
237	"\tintegrator, you must set tauThermostat_.\n");
238
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	fkBT = info_->getNdf() kB targetTemp_;
266
267	Energy = tStats->getTotalE();
268
269	thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert);
270
271	thermostat_potential = fkBT * integralOfChidt / eConvert;
272
273	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
274
275	return conservedQuantity;
276	}
277
278
279	}//end namespace oopse