--- trunk/src/integrators/NVT.cpp 2004/09/24 04:16:43 2 +++ trunk/src/integrators/NVT.cpp 2005/10/13 22:26:47 665 @@ -1,286 +1,280 @@ -#include +/* + * 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. Acknowledgement of the program authors must be made in any + * publication of scientific results based in part on use of the + * program. An acceptable form of acknowledgement is citation of + * the article in which the program was described (Matthew + * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher + * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented + * Parallel Simulation Engine for Molecular Dynamics," + * J. Comput. Chem. 26, pp. 252-271 (2005)) + * + * 2. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 3. 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. + */ + +#include "integrators/NVT.hpp" +#include "primitives/Molecule.hpp" +#include "utils/simError.h" +#include "utils/OOPSEConstant.hpp" -#include "Atom.hpp" -#include "SRI.hpp" -#include "AbstractClasses.hpp" -#include "SimInfo.hpp" -#include "ForceFields.hpp" -#include "Thermo.hpp" -#include "ReadWrite.hpp" -#include "Integrator.hpp" -#include "simError.h" +namespace oopse { + NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6), maxIterNum_(4) { -// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 + Globals* simParams = info_->getSimParams(); -template NVT::NVT ( SimInfo *theInfo, ForceFields* the_ff): - T( theInfo, the_ff ) -{ - GenericData* data; - DoubleData * chiValue; - DoubleData * 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->getProperty(CHIVALUE_ID); - if(data){ - chiValue = dynamic_cast(data); + if (!simParams->getUseIntialExtendedSystemState()) { + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + currSnapshot->setChi(0.0); + currSnapshot->setIntegralOfChiDt(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 = OOPSE_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() { + painCave.severity = OOPSE_ERROR; + painCave.isFatal = 1; + simError(); + } else { + tauThermostat_ = simParams->getTauThermostat(); + } - int i, j; - DirectionalAtom* dAtom; - double Tb[3], ji[3]; - double mass; - double vel[3], pos[3], frc[3]; + update(); + } - double instTemp; + void NVT::doUpdate() { + oldVel_.resize(info_->getNIntegrableObjects()); + oldJi_.resize(info_->getNIntegrableObjects()); + } + void NVT::moveA() { + SimInfo::MoleculeIterator i; + Molecule::IntegrableObjectIterator j; + Molecule* mol; + StuntDouble* integrableObject; + Vector3d Tb; + Vector3d ji; + double mass; + Vector3d vel; + Vector3d pos; + Vector3d frc; - // We need the temperature at time = t for the chi update below: + 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( i=0; i < integrableObjects.size(); i++ ){ + for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - integrableObjects[i]->getVel( vel ); - integrableObjects[i]->getPos( pos ); - integrableObjects[i]->getFrc( frc ); + vel = integrableObject->getVel(); + pos = integrableObject->getPos(); + frc = integrableObject->getFrc(); - mass = integrableObjects[i]->getMass(); + mass = integrableObject->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]; - } + // velocity half step (use chi from previous step here): + //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 += dt * vel; - integrableObjects[i]->setVel( vel ); - integrableObjects[i]->setPos( pos ); + integrableObject->setVel(vel); + integrableObject->setPos(pos); - if( integrableObjects[i]->isDirectional() ){ + if (integrableObject->isDirectional()) { - // get and convert the torque to body frame + //convert the torque to body frame + Tb = integrableObject->lab2Body(integrableObject->getTrq()); - 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] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi); + ji += dt2*OOPSEConstant::energyConvert*Tb - dt2*chi *ji; + rotAlgo->rotate(integrableObject, ji, dt); - 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(); + + rattle->constraintA(); - // 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; + + chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_); + integralOfChidt += chi * dt2; -} - -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; - - // Set things up for the iteration: - - oldChi = chi; - - for( i=0; i < integrableObjects.size(); i++ ){ - - integrableObjects[i]->getVel( vel ); - - for (j=0; j < 3; j++) - oldVel[3*i + j] = vel[j]; - - if( integrableObjects[i]->isDirectional() ){ - - integrableObjects[i]->getJ( ji ); - - for (j=0; j < 3; j++) - oldJi[3*i + j] = ji[j]; - - } + currentSnapshot_->setChi(chi); + currentSnapshot_->setIntegralOfChiDt(integralOfChidt); } - // do the iteration: + void NVT::moveB() { + SimInfo::MoleculeIterator i; + Molecule::IntegrableObjectIterator j; + Molecule* mol; + StuntDouble* integrableObject; + + Vector3d Tb; + Vector3d ji; + Vector3d vel; + Vector3d frc; + double mass; + double instTemp; + int index; + // Set things up for the iteration: - for (k=0; k < 4; k++) { + double chi = currentSnapshot_->getChi(); + double oldChi = chi; + double prevChi; + double integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); - instTemp = tStats->getTemperature(); + 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(); - // evolve chi another half step using the temperature at t + dt/2 + ++index; + } + + } - prevChi = chi; - chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / - (tauThermostat*tauThermostat); + // do the iteration: - for( i=0; i < integrableObjects.size(); i++ ){ + for(int k = 0; k < maxIterNum_; k++) { + index = 0; + instTemp = thermo.getTemperature(); - integrableObjects[i]->getFrc( frc ); - integrableObjects[i]->getVel(vel); + // evolve chi another half step using the temperature at t + dt/2 - mass = integrableObjects[i]->getMass(); + prevChi = chi; + chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_); - // 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); + for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - integrableObjects[i]->setVel( vel ); + frc = integrableObject->getFrc(); + vel = integrableObject->getVel(); - if( integrableObjects[i]->isDirectional() ){ + mass = integrableObject->getMass(); - // get and convert the torque to body frame + // velocity half step + //for(j = 0; j < 3; j++) + // 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); - integrableObjects[i]->getTrq( Tb ); - integrableObjects[i]->lab2Body( Tb ); + if (integrableObject->isDirectional()) { - for (j=0; j < 3; j++) - ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi); + // get and convert the torque to body frame - integrableObjects[i]->setJ( ji ); - } - } - - if(nConstrained) - constrainB(); + Tb = integrableObject->lab2Body(integrableObject->getTrq()); - if (fabs(prevChi - chi) <= chiTolerance) break; - } + //for(j = 0; j < 3; j++) + // 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]; - integralOfChidt += dt2*chi; -} + integrableObject->setJ(ji); + } -template void NVT::resetIntegrator( void ){ - chi = 0.0; - integralOfChidt = 0.0; -} + ++index; + } + } + -template int NVT::readyCheck() { + rattle->constraintB(); - //check parent's readyCheck() first - if (T::readyCheck() == -1) - return -1; + if (fabs(prevChi - chi) <= 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; - } + integralOfChidt += dt2 * chi; - // 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; + currentSnapshot_->setChi(chi); + currentSnapshot_->setIntegralOfChiDt(integralOfChidt); } - 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() { + currentSnapshot_->setChi(0.0); + currentSnapshot_->setIntegralOfChiDt(0.0); } + + double NVT::calcConservedQuantity() { - return 1; + 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_; -} + Energy = thermo.getTotalE(); -template double NVT::getConservedQuantity(void){ + thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * OOPSEConstant::energyConvert); - double conservedQuantity; - double fkBT; - double Energy; - double thermostat_kinetic; - double thermostat_potential; + thermostat_potential = fkBT * integralOfChidt / OOPSEConstant::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 oopse