--- trunk/OOPSE-4/src/integrators/NVT.cpp 2004/10/21 16:22:01 1625 +++ branches/new_design/OOPSE-4/src/integrators/NVT.cpp 2004/11/22 20:55:52 1765 @@ -1,286 +1,279 @@ -#include +#inlcude "integrators/NVT.hpp" -#include "primitives/Atom.hpp" -#include "primitives/SRI.hpp" -#include "primitives/AbstractClasses.hpp" -#include "brains/SimInfo.hpp" -#include "UseTheForce/ForceFields.hpp" -#include "brains/Thermo.hpp" -#include "io/ReadWrite.hpp" -#include "integrators/Integrator.hpp" -#include "utils/simError.h" +namespace oopse { +NVT::NVT(SimInfo* Info) : VelocityVerletIntegrator(info){ + GenericData *data; + DoubleGenericData *chiValue; + DoubleGenericData *integralOfChidtValue; -// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 + chiValue = NULL; + integralOfChidtValue = NULL; -template NVT::NVT ( SimInfo *theInfo, ForceFields* the_ff): - T( theInfo, the_ff ) -{ - GenericData* data; - DoubleGenericData * chiValue; - DoubleGenericData * integralOfChidtValue; + chi = 0.0; + have_tau_thermostat = 0; + have_target_temp = 0; + have_chi_tolerance = 0; + integralOfChidt = 0.0; - chiValue = NULL; - integralOfChidtValue = NULL; + if (theInfo->useInitXSstate) { - chi = 0.0; - have_tau_thermostat = 0; - have_target_temp = 0; - have_chi_tolerance = 0; - integralOfChidt = 0.0; + // retrieve chi and integralOfChidt from simInfo + data = info->getPropertyByName(CHIVALUE_ID); + if (data) { + chiValue = dynamic_cast(data); + } - if( theInfo->useInitXSstate ){ + data = info->getPropertyByName(INTEGRALOFCHIDT_ID); - // retrieve chi and integralOfChidt from simInfo - data = info->getProperty(CHIVALUE_ID); - if(data){ - chiValue = dynamic_cast(data); + if (data) { + integralOfChidtValue = dynamic_cast(data); + } + + // chi and integralOfChidt should appear by pair + if (chiValue && integralOfChidtValue) { + chi = chiValue->getData(); + integralOfChidt = integralOfChidtValue->getData(); + } } - - data = info->getProperty(INTEGRALOFCHIDT_ID); - if(data){ - integralOfChidtValue = dynamic_cast(data); - } - - // chi and integralOfChidt should appear by pair - if(chiValue && integralOfChidtValue){ - chi = chiValue->getData(); - integralOfChidt = integralOfChidtValue->getData(); - } - } - oldVel = new double[3*integrableObjects.size()]; - oldJi = new double[3*integrableObjects.size()]; + oldVel_.resize(info_->getNIntegrableObjects()); + oldJi_.resize(info_->getNIntegrableObjects()); } -template NVT::~NVT() { - delete[] oldVel; - delete[] oldJi; +NVT::~NVT() { } -template void NVT::moveA() { +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; - int i, j; - DirectionalAtom* dAtom; - double Tb[3], ji[3]; - double mass; - double vel[3], pos[3], frc[3]; + double instTemp; - double instTemp; + // We need the temperature at time = t for the chi update below: - // We need the temperature at time = t for the chi update below: + instTemp = tStats->getTemperature(); - 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)) { - for( i=0; i < integrableObjects.size(); i++ ){ + vel = integrableObject->getVel(); + pos = integrableObject->getPos(); + frc = integrableObject->getFrc(); - integrableObjects[i]->getVel( vel ); - integrableObjects[i]->getPos( pos ); - integrableObjects[i]->getFrc( frc ); + mass = integrableObject->getMass(); - mass = integrableObjects[i]->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; - for (j=0; j < 3; j++) { - // velocity half step (use chi from previous step here): - vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi); - // position whole step - pos[j] += dt * vel[j]; - } + integrableObject->setVel(vel); + integrableObject->setPos(pos); - integrableObjects[i]->setVel( vel ); - integrableObjects[i]->setPos( pos ); + if (integrableObject->isDirectional()) { - if( integrableObjects[i]->isDirectional() ){ + // get and convert the torque to body frame - // get and convert the torque to body frame + Tb = integrableObject->getTrq(); + integrableObject->lab2Body(Tb); - integrableObjects[i]->getTrq( Tb ); - integrableObjects[i]->lab2Body( Tb ); + // get the angular momentum, and propagate a half step - // get the angular momentum, and propagate a half step + ji = integrableObject->getJ(); - integrableObjects[i]->getJ( ji ); + //ji[j] += halfStep * (Tb[j] * eConvert - ji[j]*chi); + ji += halfStep*eConvert*Tb - halfStep*ch *ji; + this->rotationPropagation(integrableObject, ji); - for (j=0; j < 3; j++) - ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); + integrableObject->setJ(ji); + } + } - this->rotationPropagation( integrableObjects[i], ji ); - - integrableObjects[i]->setJ( ji ); } - } - - if(nConstrained) - constrainA(); + + if (nConstrained) + constrainA(); - // Finally, evolve chi a half step (just like a velocity) using - // temperature at time t, not time t+dt/2 + // Finally, evolve chi a half step (just like a velocity) using + // temperature at time t, not time t+dt/2 - //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n"; - - chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat); - integralOfChidt += chi*dt2; + 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); } -template void NVT::moveB( void ){ - int i, j, k; - double Tb[3], ji[3]; - double vel[3], frc[3]; - double mass; - double instTemp; - double oldChi, prevChi; +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: - // Set things up for the iteration: + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + double chi = currSnapshot->getChi(); + double integralOfChidt = currSnapshot->getIntegralOfChiDt(); + oldChi = chi; - 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; + } - for( i=0; i < integrableObjects.size(); i++ ){ + // do the iteration: - integrableObjects[i]->getVel( vel ); + for(int k = 0; k < maxIterNum_; k++) { + index = 0; + instTemp = tStats->getTemperature(); - for (j=0; j < 3; j++) - oldVel[3*i + j] = vel[j]; + // evolve chi another half step using the temperature at t + dt/2 - if( integrableObjects[i]->isDirectional() ){ + prevChi = chi; + chi = oldChi + halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_); - integrableObjects[i]->getJ( ji ); + for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - for (j=0; j < 3; j++) - oldJi[3*i + j] = ji[j]; + frc = integrableObject->getFrc(); + vel = integrableObject->getVel(); - } - } + mass = integrableObject->getMass(); - // do the iteration: + // 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); - for (k=0; k < 4; k++) { + if (integrableObject->isDirectional()) { - instTemp = tStats->getTemperature(); + // get and convert the torque to body frame - // evolve chi another half step using the temperature at t + dt/2 + Tb = integrableObject->getTrq(); + integrableObject->lab2Body(Tb); - prevChi = chi; - chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / - (tauThermostat*tauThermostat); + //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]; - for( i=0; i < integrableObjects.size(); i++ ){ + integrableObject->setJ(ji); + } + } + } + - integrableObjects[i]->getFrc( frc ); - integrableObjects[i]->getVel(vel); + if (nConstrained) + constrainB(); - mass = integrableObjects[i]->getMass(); + if (fabs(prevChi - chi) <= chiTolerance) + break; - // 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); - - integrableObjects[i]->setVel( vel ); - - if( integrableObjects[i]->isDirectional() ){ - - // get and convert the torque to body frame - - integrableObjects[i]->getTrq( Tb ); - integrableObjects[i]->lab2Body( Tb ); - - for (j=0; j < 3; j++) - ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi); - - integrableObjects[i]->setJ( ji ); - } + ++index; } - - if(nConstrained) - constrainB(); - if (fabs(prevChi - chi) <= chiTolerance) break; - } + integralOfChidt += halfStep * chi; - integralOfChidt += dt2*chi; + currSnapshot->setChi(chi); + currSnapshot->setIntegralOfChiDt(integralOfChidt); } -template void NVT::resetIntegrator( void ){ +template +int NVT::readyCheck() { - chi = 0.0; - integralOfChidt = 0.0; -} + //check parent's readyCheck() first + if (T::readyCheck() == -1) + return -1; -template int NVT::readyCheck() { + // First check to see if we have a target temperature. + // Not having one is fatal. - //check parent's readyCheck() first - if (T::readyCheck() == -1) - return -1; + 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; + } - // First check to see if we have a target temperature. - // Not having one is fatal. + // We must set tauThermostat_. - 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; - } + if (!have_tau_thermostat) { + sprintf(painCave.errMsg, "If you use the constant temperature\n" + "\tintegrator, you must set tauThermostat_.\n"); - // We must set tauThermostat. + painCave.severity = OOPSE_ERROR; + painCave.isFatal = 1; + simError(); + return -1; + } - 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(); + } - 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; - + return 1; } -template double NVT::getConservedQuantity(void){ +template +double NVT::getConservedQuantity(void) { + double conservedQuantity; + double fkBT; + double Energy; + double thermostat_kinetic; + double thermostat_potential; - double conservedQuantity; - double fkBT; - double Energy; - double thermostat_kinetic; - double thermostat_potential; + fkBT = info_->getNdf() *kB *targetTemp_; - fkBT = (double)(info->ndf) * kB * targetTemp; + Energy = tStats->getTotalE(); - Energy = tStats->getTotalE(); + thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert); - thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / - (2.0 * eConvert); + thermostat_potential = fkBT * integralOfChidt / eConvert; - thermostat_potential = fkBT * integralOfChidt / eConvert; + conservedQuantity = Energy + thermostat_kinetic + thermostat_potential; - conservedQuantity = Energy + thermostat_kinetic + thermostat_potential; - - return conservedQuantity; + return conservedQuantity; } -template string NVT::getAdditionalParameters(void){ - string parameters; - const int BUFFERSIZE = 2000; // size of the read buffer - char buffer[BUFFERSIZE]; - sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt); - parameters += buffer; - - return parameters; -} +}//end namespace oopse