--- trunk/src/integrators/NPT.cpp 2005/10/13 22:26:47 665 +++ trunk/src/integrators/NPT.cpp 2013/06/16 15:15:42 1879 @@ -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,7 +60,7 @@ // // 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) { @@ -68,15 +69,14 @@ namespace oopse { 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()) { 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 { targetTemp = simParams->getTargetTemp(); @@ -85,9 +85,9 @@ namespace oopse { // 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 { @@ -107,7 +107,7 @@ namespace oopse { 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 { @@ -117,13 +117,13 @@ namespace oopse { tt2 = tauThermostat * tauThermostat; tb2 = tauBarostat * tauBarostat; - update(); + updateSizes(); } NPT::~NPT() { } - void NPT::doUpdate() { + void NPT::doUpdateSizes() { oldPos.resize(info_->getNIntegrableObjects()); oldVel.resize(info_->getNIntegrableObjects()); @@ -135,80 +135,88 @@ namespace oopse { 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); // 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); - if (integrableObject->isDirectional()) { + vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc; + sd->setVel(vel); + if (sd->isDirectional()) { + // get and convert the torque to body frame - Tb = integrableObject->lab2Body(integrableObject->getTrq()); + Tb = sd->lab2Body(sd->getTrq()); // get the angular momentum, and propagate a half step - ji = integrableObject->getJ(); + ji = sd->getJ(); - //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi); - ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji; + 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(); + } } @@ -216,31 +224,32 @@ namespace oopse { 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)) { + 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)) { + vel = sd->getVel(); + pos = sd->getPos(); + this->getPosScale(pos, COM, index, sc); pos = oldPos[index] + dt * (vel + sc); - integrableObject->setPos(pos); + sd->setPos(pos); ++index; } } - rattle->constraintA(); + rattle_->constraintA(); } // Scale the box after all the positions have been moved: this->scaleSimBox(); - currentSnapshot_->setChi(chi); - currentSnapshot_->setIntegralOfChiDt(integralOfChidt); + snap->setThermostat(thermostat); saveEta(); } @@ -249,31 +258,34 @@ namespace oopse { 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(); + oldVel[index] = sd->getVel(); + + if (sd->isDirectional()) + oldJi[index] = sd->getJ(); + ++index; } } @@ -286,68 +298,70 @@ namespace oopse { 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; + prevChi = thermostat.first; + thermostat.first = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2; //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)) { + 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); // 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); + vel = oldVel[index] + + dt2*PhysicalConstants::energyConvert/mass* frc + - dt2*sc; - if (integrableObject->isDirectional()) { + sd->setVel(vel); + + if (sd->isDirectional()) { // get and convert the torque to body frame - Tb = integrableObject->lab2Body(integrableObject->getTrq()); + Tb = sd->lab2Body(sd->getTrq()); - //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); + 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()) + 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(){ - currentSnapshot_->setChi(0.0); - currentSnapshot_->setIntegralOfChiDt(0.0); - resetEta(); + snap->setThermostat(make_pair(0.0, 0.0)); + resetEta(); } - - void NPT::resetEta() { - Mat3x3d etaMat(0.0); - currentSnapshot_->setEta(etaMat); - } - + void NPT::resetEta() { + Mat3x3d etaMat(0.0); + snap->setBarostat(etaMat); + } }