--- trunk/src/integrators/NPT.cpp 2005/01/12 22:41:40 246 +++ trunk/src/integrators/NPT.cpp 2015/03/07 21:41:51 2071 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -6,19 +6,10 @@ * 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 + * 1. 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 + * 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. @@ -37,6 +28,16 @@ * 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 @@ -46,7 +47,7 @@ #include "integrators/NPT.hpp" #include "math/SquareMatrix3.hpp" #include "primitives/Molecule.hpp" -#include "utils/OOPSEConstant.hpp" +#include "utils/PhysicalConstants.hpp" #include "utils/simError.h" // Basic isotropic thermostating and barostating via the Melchionna @@ -59,283 +60,309 @@ // // Hoover, W. G., 1986, Phys. Rev. A, 34, 2499. -namespace oopse { +namespace OpenMD { -NPT::NPT(SimInfo* info) : - VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) { + NPT::NPT(SimInfo* info) : + VelocityVerletIntegrator(info), etaTolerance(1e-6), chiTolerance(1e-6), + maxIterNum_(4) { - Globals* simParams = info_->getSimParams(); + Globals* simParams = info_->getSimParams(); - if (!simParams->getUseInitXSstate()) { + if (!simParams->getUseIntialExtendedSystemState()) { Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); - currSnapshot->setChi(0.0); - currSnapshot->setIntegralOfChiDt(0.0); - currSnapshot->setEta(Mat3x3d(0.0)); - } + currSnapshot->setThermostat(make_pair(0.0, 0.0)); + currSnapshot->setBarostat(Mat3x3d(0.0)); + } - if (!simParams->haveTargetTemp()) { + if (!simParams->haveTargetTemp()) { sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n"); painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; simError(); - } else { + } else { targetTemp = simParams->getTargetTemp(); - } + } - // We must set tauThermostat - if (!simParams->haveTauThermostat()) { + // We must set tauThermostat + if (!simParams->haveTauThermostat()) { sprintf(painCave.errMsg, "If you use the constant temperature\n" - "\tintegrator, you must set tauThermostat_.\n"); + "\tintegrator, you must set tauThermostat.\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); - } else { + } else { tauThermostat = simParams->getTauThermostat(); - } + } - if (!simParams->haveTargetPressure()) { + if (!simParams->haveTargetPressure()) { sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n" - " without a targetPressure!\n"); + " without a targetPressure!\n"); painCave.isFatal = 1; simError(); - } else { + } else { targetPressure = simParams->getTargetPressure(); - } + } - if (!simParams->haveTauBarostat()) { + if (!simParams->haveTauBarostat()) { sprintf(painCave.errMsg, "If you use the NPT integrator, you must set tauBarostat.\n"); - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); - } else { + } else { tauBarostat = simParams->getTauBarostat(); - } + } - tt2 = tauThermostat * tauThermostat; - tb2 = tauBarostat * tauBarostat; + tt2 = tauThermostat * tauThermostat; + tb2 = tauBarostat * tauBarostat; - update(); -} + updateSizes(); + } -NPT::~NPT() { -} + NPT::~NPT() { + } -void NPT::doUpdate() { + void NPT::doUpdateSizes() { oldPos.resize(info_->getNIntegrableObjects()); oldVel.resize(info_->getNIntegrableObjects()); oldJi.resize(info_->getNIntegrableObjects()); -} + } -void NPT::moveA() { + void NPT::moveA() { SimInfo::MoleculeIterator i; Molecule::IntegrableObjectIterator j; Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; Vector3d Tb, ji; - double mass; + RealType mass; Vector3d vel; Vector3d pos; Vector3d frc; Vector3d sc; int index; - chi= currentSnapshot_->getChi(); - integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); + thermostat = snap->getThermostat(); loadEta(); instaTemp =thermo.getTemperature(); press = thermo.getPressureTensor(); - instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0; + instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0; instaVol =thermo.getVolume(); - Vector3d COM = info_->getCom(); + Vector3d COM = thermo.getCom(); //evolve velocity half step calcVelScale(); - for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { - vel = integrableObject->getVel(); - frc = integrableObject->getFrc(); + vel = sd->getVel(); + frc = sd->getFrc(); - mass = integrableObject->getMass(); + mass = sd->getMass(); - getVelScaleA(sc, vel); + getVelScaleA(sc, vel); - // velocity half step (use chi from previous step here): - //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]); - vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc; - integrableObject->setVel(vel); + // velocity half step (use chi from previous step here): - if (integrableObject->isDirectional()) { + vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc; + sd->setVel(vel); - // get and convert the torque to body frame + if (sd->isDirectional()) { - Tb = integrableObject->lab2Body(integrableObject->getTrq()); + // get and convert the torque to body frame - // get the angular momentum, and propagate a half step + Tb = sd->lab2Body(sd->getTrq()); - ji = integrableObject->getJ(); + // get the angular momentum, and propagate a half step - //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi); - ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji; + ji = sd->getJ(); + + ji += dt2*PhysicalConstants::energyConvert * Tb + - dt2*thermostat.first* ji; - rotAlgo->rotate(integrableObject, ji, dt); + rotAlgo_->rotate(sd, ji, dt); - integrableObject->setJ(ji); - } + sd->setJ(ji); + } - } + } } // evolve chi and eta half step - chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2; + thermostat.first += dt2 * (instaTemp / targetTemp - 1.0) / tt2; evolveEtaA(); //calculate the integral of chidt - integralOfChidt += dt2 * chi; + thermostat.second += dt2 * thermostat.first; + flucQ_->moveA(); + + index = 0; - for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { - oldPos[index++] = integrableObject->getPos(); - } + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { + + oldPos[index++] = sd->getPos(); + + } } //the first estimation of r(t+dt) is equal to r(t) for(int k = 0; k < maxIterNum_; k++) { - index = 0; - for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + index = 0; + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { - vel = integrableObject->getVel(); - pos = integrableObject->getPos(); + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { - this->getPosScale(pos, COM, index, sc); + vel = sd->getVel(); + pos = sd->getPos(); - pos = oldPos[index] + dt * (vel + sc); - integrableObject->setPos(pos); + this->getPosScale(pos, COM, index, sc); - ++index; - } - } + pos = oldPos[index] + dt * (vel + sc); + sd->setPos(pos); - rattle->constraintA(); + ++index; + } + } + + rattle_->constraintA(); } // Scale the box after all the positions have been moved: this->scaleSimBox(); - currentSnapshot_->setChi(chi); - currentSnapshot_->setIntegralOfChiDt(integralOfChidt); + snap->setThermostat(thermostat); saveEta(); -} + } -void NPT::moveB(void) { + void NPT::moveB(void) { SimInfo::MoleculeIterator i; Molecule::IntegrableObjectIterator j; Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; int index; Vector3d Tb; Vector3d ji; Vector3d sc; Vector3d vel; Vector3d frc; - double mass; + RealType mass; + thermostat = snap->getThermostat(); + RealType oldChi = thermostat.first; + RealType prevChi; - chi= currentSnapshot_->getChi(); - integralOfChidt = currentSnapshot_->getIntegralOfChiDt(); - double oldChi = chi; - double prevChi; - loadEta(); //save velocity and angular momentum index = 0; - for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { - oldVel[index] = integrableObject->getVel(); - oldJi[index] = integrableObject->getJ(); - ++index; - } + oldVel[index] = sd->getVel(); + + if (sd->isDirectional()) + oldJi[index] = sd->getJ(); + + ++index; + } } // do the iteration: instaVol =thermo.getVolume(); for(int k = 0; k < maxIterNum_; k++) { - instaTemp =thermo.getTemperature(); - instaPress =thermo.getPressure(); + instaTemp =thermo.getTemperature(); + instaPress =thermo.getPressure(); - // evolve chi another half step using the temperature at t + dt/2 - prevChi = chi; - chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2; + // evolve chi another half step using the temperature at t + dt/2 + prevChi = thermostat.first; + thermostat.first = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2; - //evolve eta - this->evolveEtaB(); - this->calcVelScale(); + //evolve eta + this->evolveEtaB(); + this->calcVelScale(); - index = 0; - for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + index = 0; + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { - frc = integrableObject->getFrc(); - vel = integrableObject->getVel(); + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { - mass = integrableObject->getMass(); + frc = sd->getFrc(); + mass = sd->getMass(); - getVelScaleB(sc, index); + getVelScaleB(sc, index); - // velocity half step - //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]); - vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc; - integrableObject->setVel(vel); + // velocity half step + vel = oldVel[index] + + dt2*PhysicalConstants::energyConvert/mass* frc + - dt2*sc; - if (integrableObject->isDirectional()) { - // get and convert the torque to body frame - Tb = integrableObject->lab2Body(integrableObject->getTrq()); + sd->setVel(vel); - //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); - } + if (sd->isDirectional()) { + // get and convert the torque to body frame + Tb = sd->lab2Body(sd->getTrq()); - ++index; - } - } + ji = oldJi[index] + + dt2*PhysicalConstants::energyConvert*Tb + - dt2*thermostat.first*oldJi[index]; + + sd->setJ(ji); + } + + ++index; + } + } - rattle->constraintB(); + rattle_->constraintB(); - if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged()) - break; + if ((fabs(prevChi - thermostat.first) <= chiTolerance) && + this->etaConverged()) + break; } //calculate integral of chidt - integralOfChidt += dt2 * chi; + thermostat.second += dt2 * thermostat.first; - currentSnapshot_->setChi(chi); - currentSnapshot_->setIntegralOfChiDt(integralOfChidt); + snap->setThermostat(thermostat); + flucQ_->moveB(); saveEta(); -} + } + void NPT::resetIntegrator(){ + snap->setThermostat(make_pair(0.0, 0.0)); + resetEta(); + } + + void NPT::resetEta() { + Mat3x3d etaMat(0.0); + snap->setBarostat(etaMat); + } }