--- branches/new_design/OOPSE-4/src/integrators/NVT.cpp 2004/11/22 20:55:52 1765 +++ branches/new_design/OOPSE-4/src/integrators/NVT.cpp 2004/12/09 16:22:42 1871 @@ -1,53 +1,55 @@ -#inlcude "integrators/NVT.hpp" +#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){ - GenericData *data; - DoubleGenericData *chiValue; - DoubleGenericData *integralOfChidtValue; +static IntegratorBuilder* NVTCreator = new IntegratorBuilder("NVT"); - chiValue = NULL; - integralOfChidtValue = NULL; +NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) { - chi = 0.0; - have_tau_thermostat = 0; - have_target_temp = 0; - have_chi_tolerance = 0; - integralOfChidt = 0.0; + Globals* simParams = info_->getSimParams(); - if (theInfo->useInitXSstate) { + 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(); + } - // retrieve chi and integralOfChidt from simInfo - data = info->getPropertyByName(CHIVALUE_ID); + // We must set tauThermostat_. - if (data) { - chiValue = dynamic_cast(data); - } + if (!simParams->haveTauThermostat()) { + sprintf(painCave.errMsg, "If you use the constant temperature\n" + "\tintegrator, you must set tauThermostat_.\n"); - data = info->getPropertyByName(INTEGRALOFCHIDT_ID); - - if (data) { - integralOfChidtValue = dynamic_cast(data); - } - - // chi and integralOfChidt should appear by pair - if (chiValue && integralOfChidtValue) { - chi = chiValue->getData(); - integralOfChidt = integralOfChidtValue->getData(); - } + painCave.severity = OOPSE_ERROR; + painCave.isFatal = 1; + simError(); + } else { + tauThermostat_ = simParams->getTauThermostat(); } - oldVel_.resize(info_->getNIntegrableObjects()); - oldJi_.resize(info_->getNIntegrableObjects()); + update(); } -NVT::~NVT() { +void NVT::doUpdate() { + oldVel_.resize(info_->getNIntegrableObjects()); + oldJi_.resize(info_->getNIntegrableObjects()); } - void NVT::moveA() { - typename SimInfo::MoleculeIterator i; - typename Molecule::IntegrableObjectIterator j; + SimInfo::MoleculeIterator i; + Molecule::IntegrableObjectIterator j; Molecule* mol; StuntDouble* integrableObject; Vector3d Tb; @@ -57,11 +59,12 @@ void NVT::moveA() { Vector3d pos; Vector3d frc; - double instTemp; - + double chi = currentSnapshot_->getChi(); + double integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); + // We need the temperature at time = t for the chi update below: - instTemp = tStats->getTemperature(); + double instTemp = thermo.getTemperature(); for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; @@ -74,30 +77,28 @@ void NVT::moveA() { 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; + //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 += fullStep * vel; + pos += dt * vel; integrableObject->setVel(vel); integrableObject->setPos(pos); if (integrableObject->isDirectional()) { - // get and convert the torque to body frame + //convert the torque to body frame + Tb = integrableObject->lab2Body(integrableObject->getTrq()); - 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); + //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); } @@ -105,26 +106,22 @@ void NVT::moveA() { } - if (nConstrained) - constrainA(); + //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 += halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_); - integralOfChidt += chi * halfStep; + chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_); + integralOfChidt += chi * dt2; - currSnapshot->setChi(chi); - currSnapshot->setIntegralOfChiDt(integralOfChidt); + currentSnapshot_->setChi(chi); + currentSnapshot_->setIntegralOfChiDt(integralOfChidt); } void NVT::moveB() { - typename SimInfo::MoleculeIterator i; - typename Molecule::IntegrableObjectIterator j; + SimInfo::MoleculeIterator i; + Molecule::IntegrableObjectIterator j; Molecule* mol; StuntDouble* integrableObject; @@ -134,14 +131,13 @@ void NVT::moveB() { 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; + 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)) { @@ -149,20 +145,22 @@ void NVT::moveB() { integrableObject = mol->nextIntegrableObject(j)) { oldVel_[index] = integrableObject->getVel(); oldJi_[index] = integrableObject->getJ(); + + ++index; } - ++index; + } // do the iteration: for(int k = 0; k < maxIterNum_; k++) { index = 0; - instTemp = tStats->getTemperature(); + instTemp = thermo.getTemperature(); // evolve chi another half step using the temperature at t + dt/2 prevChi = chi; - chi = oldChi + halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_); + 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; @@ -175,8 +173,8 @@ void NVT::moveB() { // 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]; + // 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); @@ -184,92 +182,53 @@ void NVT::moveB() { // get and convert the torque to body frame - Tb = integrableObject->getTrq(); - integrableObject->lab2Body(Tb); + Tb = integrableObject->lab2Body(integrableObject->getTrq()); //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]; + // 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; } } - if (nConstrained) - constrainB(); + //constraintAlgorithm->doConstrainB(); - if (fabs(prevChi - chi) <= chiTolerance) + if (fabs(prevChi - chi) <= chiTolerance_) break; - ++index; } - integralOfChidt += halfStep * chi; + integralOfChidt += dt2 * chi; - currSnapshot->setChi(chi); - currSnapshot->setIntegralOfChiDt(integralOfChidt); + currentSnapshot_->setChi(chi); + currentSnapshot_->setIntegralOfChiDt(integralOfChidt); } -template -int NVT::readyCheck() { - //check parent's readyCheck() first - if (T::readyCheck() == -1) - return -1; +double NVT::calcConservedQuantity() { - // 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 -double NVT::getConservedQuantity(void) { + 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_; - fkBT = info_->getNdf() *kB *targetTemp_; + Energy = thermo.getTotalE(); - Energy = tStats->getTotalE(); + thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert); - thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert); + thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::energyConvert; - thermostat_potential = fkBT * integralOfChidt / eConvert; - conservedQuantity = Energy + thermostat_kinetic + thermostat_potential; return conservedQuantity;