--- trunk/src/integrators/NVT.cpp 2004/09/24 16:27:58 3 +++ trunk/src/integrators/NVT.cpp 2013/06/16 15:15:42 1879 @@ -1,286 +1,283 @@ -#include - -#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" +/* + * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * + * The University of Notre Dame grants you ("Licensee") a + * non-exclusive, royalty free, license to use, modify and + * redistribute this software in source and binary code form, provided + * that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the + * distribution. + * + * This software is provided "AS IS," without a warranty of any + * kind. All express or implied conditions, representations and + * warranties, including any implied warranty of merchantability, + * fitness for a particular purpose or non-infringement, are hereby + * excluded. The University of Notre Dame and its licensors shall not + * be liable for any damages suffered by licensee as a result of + * using, modifying or distributing the software or its + * derivatives. In no event will the University of Notre Dame or its + * licensors be liable for any lost revenue, profit or data, or for + * direct, indirect, special, consequential, incidental or punitive + * damages, however caused and regardless of the theory of liability, + * arising out of the use of or inability to use software, even if the + * University of Notre Dame has been advised of the possibility of + * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). + */ + +#include "integrators/NVT.hpp" +#include "primitives/Molecule.hpp" #include "utils/simError.h" +#include "utils/PhysicalConstants.hpp" +namespace OpenMD { -// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 + NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6), maxIterNum_(4) { -template NVT::NVT ( SimInfo *theInfo, ForceFields* the_ff): - T( theInfo, the_ff ) -{ - GenericData* data; - DoubleData * chiValue; - DoubleData * integralOfChidtValue; + Globals* simParams = info_->getSimParams(); - 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->getProperty(CHIVALUE_ID); - if(data){ - chiValue = dynamic_cast(data); + if (!simParams->getUseIntialExtendedSystemState()) { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); + snap->setThermostat(make_pair(0.0, 0.0)); } - data = info->getProperty(INTEGRALOFCHIDT_ID); - if(data){ - integralOfChidtValue = dynamic_cast(data); + if (!simParams->haveTargetTemp()) { + sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n"); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } else { + targetTemp_ = simParams->getTargetTemp(); } - - // 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()]; -} + // We must set tauThermostat. -template NVT::~NVT() { - delete[] oldVel; - delete[] oldJi; -} + if (!simParams->haveTauThermostat()) { + sprintf(painCave.errMsg, "If you use the constant temperature\n" + "\tintegrator, you must set tauThermostat.\n"); -template void NVT::moveA() { - - int i, j; - DirectionalAtom* dAtom; - double Tb[3], ji[3]; - double mass; - double vel[3], pos[3], frc[3]; - - double instTemp; - - // We need the temperature at time = t for the chi update below: - - instTemp = tStats->getTemperature(); - - for( i=0; i < integrableObjects.size(); i++ ){ - - integrableObjects[i]->getVel( vel ); - integrableObjects[i]->getPos( pos ); - integrableObjects[i]->getFrc( frc ); - - mass = integrableObjects[i]->getMass(); - - 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]; + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + } else { + tauThermostat_ = simParams->getTauThermostat(); } - integrableObjects[i]->setVel( vel ); - integrableObjects[i]->setPos( pos ); + updateSizes(); + } - if( integrableObjects[i]->isDirectional() ){ + void NVT::doUpdateSizes() { + oldVel_.resize(info_->getNIntegrableObjects()); + oldJi_.resize(info_->getNIntegrableObjects()); + } - // get and convert the torque to body frame + void NVT::moveA() { + SimInfo::MoleculeIterator i; + Molecule::IntegrableObjectIterator j; + Molecule* mol; + StuntDouble* sd; + Vector3d Tb; + Vector3d ji; + RealType mass; + Vector3d vel; + Vector3d pos; + Vector3d frc; - integrableObjects[i]->getTrq( Tb ); - integrableObjects[i]->lab2Body( Tb ); + pair thermostat = snap->getThermostat(); - // get the angular momentum, and propagate a half step + // We need the temperature at time = t for the chi update below: - integrableObjects[i]->getJ( ji ); + RealType instTemp = thermo.getTemperature(); - for (j=0; j < 3; j++) - ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { - this->rotationPropagation( integrableObjects[i], ji ); + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { - integrableObjects[i]->setJ( ji ); - } - } - - if(nConstrained) - constrainA(); + vel = sd->getVel(); + pos = sd->getPos(); + frc = sd->getFrc(); - // Finally, evolve chi a half step (just like a velocity) using - // temperature at time t, not time t+dt/2 + mass = sd->getMass(); - //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n"; - - chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat); - integralOfChidt += chi*dt2; + // velocity half step (use chi from previous step here): + vel += dt2 *PhysicalConstants::energyConvert/mass*frc + - dt2*thermostat.first*vel; + + // position whole step + pos += dt * vel; -} + sd->setVel(vel); + sd->setPos(pos); -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; + if (sd->isDirectional()) { - // Set things up for the iteration: + //convert the torque to body frame + Tb = sd->lab2Body(sd->getTrq()); - oldChi = chi; + // get the angular momentum, and propagate a half step - for( i=0; i < integrableObjects.size(); i++ ){ + ji = sd->getJ(); - integrableObjects[i]->getVel( vel ); + ji += dt2*PhysicalConstants::energyConvert*Tb + - dt2*thermostat.first *ji; - for (j=0; j < 3; j++) - oldVel[3*i + j] = vel[j]; + rotAlgo_->rotate(sd, ji, dt); - if( integrableObjects[i]->isDirectional() ){ + sd->setJ(ji); + } + } - integrableObjects[i]->getJ( ji ); + } + + flucQ_->moveA(); + rattle_->constraintA(); - for (j=0; j < 3; j++) - oldJi[3*i + j] = ji[j]; + // Finally, evolve chi a half step (just like a velocity) using + // temperature at time t, not time t+dt/2 - } + thermostat.first += dt2 * (instTemp / targetTemp_ - 1.0) + / (tauThermostat_ * tauThermostat_); + thermostat.second += thermostat.first * dt2; + + snap->setThermostat(thermostat); } - // do the iteration: + void NVT::moveB() { + SimInfo::MoleculeIterator i; + Molecule::IntegrableObjectIterator j; + Molecule* mol; + StuntDouble* sd; + + Vector3d Tb; + Vector3d ji; + Vector3d vel; + Vector3d frc; + RealType mass; + RealType instTemp; + int index; + // Set things up for the iteration: - for (k=0; k < 4; k++) { + pair thermostat = snap->getThermostat(); + RealType oldChi = thermostat.first; + RealType prevChi; - instTemp = tStats->getTemperature(); + index = 0; + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { - // evolve chi another half step using the temperature at t + dt/2 + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { - prevChi = chi; - chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / - (tauThermostat*tauThermostat); + oldVel_[index] = sd->getVel(); + + if (sd->isDirectional()) + oldJi_[index] = sd->getJ(); + + ++index; + } + } - for( i=0; i < integrableObjects.size(); i++ ){ + // do the iteration: - integrableObjects[i]->getFrc( frc ); - integrableObjects[i]->getVel(vel); + for(int k = 0; k < maxIterNum_; k++) { + index = 0; + instTemp = thermo.getTemperature(); - mass = integrableObjects[i]->getMass(); + // evolve chi another half step using the temperature at t + dt/2 - // 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); + prevChi = thermostat.first; + thermostat.first = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) + / (tauThermostat_ * tauThermostat_); - integrableObjects[i]->setVel( vel ); + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { - if( integrableObjects[i]->isDirectional() ){ + frc = sd->getFrc(); + mass = sd->getMass(); - // get and convert the torque to body frame + // velocity half step - integrableObjects[i]->getTrq( Tb ); - integrableObjects[i]->lab2Body( Tb ); + vel = oldVel_[index] + + dt2/mass*PhysicalConstants::energyConvert * frc + - dt2*thermostat.first*oldVel_[index]; + + sd->setVel(vel); - for (j=0; j < 3; j++) - ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi); + if (sd->isDirectional()) { - integrableObjects[i]->setJ( ji ); - } - } - - if(nConstrained) - constrainB(); + // get and convert the torque to body frame - if (fabs(prevChi - chi) <= chiTolerance) break; - } + Tb = sd->lab2Body(sd->getTrq()); - integralOfChidt += dt2*chi; -} + ji = oldJi_[index] + dt2*PhysicalConstants::energyConvert*Tb + - dt2*thermostat.first *oldJi_[index]; -template void NVT::resetIntegrator( void ){ + sd->setJ(ji); + } - chi = 0.0; - integralOfChidt = 0.0; -} -template int NVT::readyCheck() { + ++index; + } + } + + rattle_->constraintB(); - //check parent's readyCheck() first - if (T::readyCheck() == -1) - return -1; + if (fabs(prevChi - thermostat.first) <= chiTolerance_) + break; - // 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; - } + flucQ_->moveB(); - // 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; + thermostat.second += dt2 * thermostat.first; + snap->setThermostat(thermostat); } - 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(); + void NVT::resetIntegrator() { + snap->setThermostat(make_pair(0.0, 0.0)); } + + RealType NVT::calcConservedQuantity() { - return 1; + pair thermostat = snap->getThermostat(); + RealType conservedQuantity; + RealType fkBT; + RealType Energy; + RealType thermostat_kinetic; + RealType thermostat_potential; + + fkBT = info_->getNdf() *PhysicalConstants::kB *targetTemp_; -} + Energy = thermo.getTotalEnergy(); -template double NVT::getConservedQuantity(void){ + thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * thermostat.first * thermostat.first / (2.0 * PhysicalConstants::energyConvert); - double conservedQuantity; - double fkBT; - double Energy; - double thermostat_kinetic; - double thermostat_potential; + thermostat_potential = fkBT * thermostat.second / PhysicalConstants::energyConvert; - fkBT = (double)(info->ndf) * kB * targetTemp; + conservedQuantity = Energy + thermostat_kinetic + thermostat_potential; - Energy = tStats->getTotalE(); + return conservedQuantity; + } - thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / - (2.0 * eConvert); - thermostat_potential = fkBT * integralOfChidt / eConvert; - - conservedQuantity = Energy + thermostat_kinetic + thermostat_potential; - - 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 OpenMD