--- branches/development/src/integrators/FluctuatingChargeNVT.cpp 2012/05/23 01:26:15 1716 +++ trunk/src/flucq/FluctuatingChargeNVT.cpp 2015/03/07 21:41:51 2071 @@ -35,7 +35,7 @@ * * [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, 24107 (2008). + * [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). */ @@ -45,82 +45,58 @@ #include "utils/simError.h" #include "utils/PhysicalConstants.hpp" + namespace OpenMD { FluctuatingChargeNVT::FluctuatingChargeNVT(SimInfo* info) : - FluctuatingChargePropagator(info), chiTolerance_ (1e-6), maxIterNum_(4), - thermo(info), - currentSnapshot_(info->getSnapshotManager()->getCurrentSnapshot()) { - - if (info_->usesFluctuatingCharges()) { - if (info_->getNFluctuatingCharges() > 0) { - - hasFlucQ_ = true; - Globals* simParams = info_->getSimParams(); - - if (simParams->haveDt()) { - dt_ = simParams->getDt(); - dt2_ = dt_ * 0.5; - } else { - sprintf(painCave.errMsg, - "FluctuatingChargeNVT Error: dt is not set\n"); - painCave.isFatal = 1; - simError(); - } - - if (!simParams->getUseIntialExtendedSystemState()) { - currentSnapshot_->setChiElectronic(0.0); - currentSnapshot_->setIntegralOfChiElectronicDt(0.0); - } - - if (!simParams->haveFlucQTargetTemp()) { - sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT " - "propagator without a flucQ.targetTemp!\n"); - painCave.isFatal = 1; - painCave.severity = OPENMD_ERROR; - simError(); - } else { - targetTemp_ = simParams->getFlucQTargetTemp(); - } - - // We must set tauThermostat. - - if (!simParams->haveFlucQtauThermostat()) { - sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n" - "\tpropagator, you must set flucQ.tauThermostat .\n"); - - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } else { - tauThermostat_ = simParams->getFlucQtauThermostat(); - } - updateSizes(); - } - } + FluctuatingChargePropagator(info), maxIterNum_(4), chiTolerance_ (1e-6), + snap(info->getSnapshotManager()->getCurrentSnapshot()), thermo(info) { } void FluctuatingChargeNVT::initialize() { - - if (!hasFlucQ_) return; - - SimInfo::MoleculeIterator i; - Molecule::FluctuatingChargeIterator j; - Molecule* mol; - Atom* atom; - - for (mol = info_->beginMolecule(i); mol != NULL; - mol = info_->nextMolecule(i)) { - for (atom = mol->beginFluctuatingCharge(j); atom != NULL; - atom = mol->nextFluctuatingCharge(j)) { - atom->setFlucQPos(0.0); - atom->setFlucQVel(0.0); + FluctuatingChargePropagator::initialize(); + if (hasFlucQ_) { + if (info_->getSimParams()->haveDt()) { + dt_ = info_->getSimParams()->getDt(); + dt2_ = dt_ * 0.5; + } else { + sprintf(painCave.errMsg, + "FluctuatingChargeNVT Error: dt is not set\n"); + painCave.isFatal = 1; + simError(); } + + if (!info_->getSimParams()->getUseIntialExtendedSystemState()) { + snap->setElectronicThermostat(make_pair(0.0, 0.0)); + } + + if (!fqParams_->haveTargetTemp()) { + sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT " + "propagator without a flucQ.targetTemp!\n"); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } else { + targetTemp_ = fqParams_->getTargetTemp(); + } + + // We must set tauThermostat. + + if (!fqParams_->haveTauThermostat()) { + sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n" + "\tpropagator, you must set flucQ.tauThermostat .\n"); + + painCave.severity = OPENMD_ERROR; + painCave.isFatal = 1; + simError(); + } else { + tauThermostat_ = fqParams_->getTauThermostat(); + } + updateSizes(); } - - cerr << "Yeah, you should probably implement this\n"; } + void FluctuatingChargeNVT::moveA() { if (!hasFlucQ_) return; @@ -131,12 +107,11 @@ namespace OpenMD { Atom* atom; RealType cvel, cpos, cfrc, cmass; - RealType chi = currentSnapshot_->getChiElectronic(); - RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt(); + pair thermostat = snap->getElectronicThermostat(); + RealType chi = thermostat.first; + RealType integralOfChidt = thermostat.second; RealType instTemp = thermo.getElectronicTemperature(); - cerr << "why are we here?\n"; - for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { for (atom = mol->beginFluctuatingCharge(j); atom != NULL; @@ -145,13 +120,13 @@ namespace OpenMD { cvel = atom->getFlucQVel(); cpos = atom->getFlucQPos(); cfrc = atom->getFlucQFrc(); - cmass = atom->getChargeMass(); - - // velocity half step - cvel += dt2_ *PhysicalConstants::energyConvert/cmass*cfrc - dt2_*chi*cvel; + cmass = atom->getChargeMass(); + + // velocity half step + cvel += dt2_ * cfrc / cmass - dt2_*chi*cvel; // position whole step cpos += dt_ * cvel; - + atom->setFlucQVel(cvel); atom->setFlucQPos(cpos); } @@ -161,13 +136,10 @@ namespace OpenMD { (tauThermostat_ * tauThermostat_); integralOfChidt += chi * dt2_; - currentSnapshot_->setChiElectronic(chi); - currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt); - + snap->setElectronicThermostat(make_pair(chi, integralOfChidt)); } void FluctuatingChargeNVT::updateSizes() { - if (!hasFlucQ_) return; oldVel_.resize(info_->getNFluctuatingCharges()); } @@ -178,10 +150,11 @@ namespace OpenMD { Molecule* mol; Atom* atom; RealType instTemp; - RealType chi = currentSnapshot_->getChiElectronic(); + pair thermostat = snap->getElectronicThermostat(); + RealType chi = thermostat.first; RealType oldChi = chi; RealType prevChi; - RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt(); + RealType integralOfChidt = thermostat.second; int index; RealType cfrc, cvel, cmass; @@ -201,7 +174,6 @@ namespace OpenMD { for(int k = 0; k < maxIterNum_; k++) { index = 0; instTemp = thermo.getElectronicTemperature(); - // evolve chi another half step using the temperature at t + dt/2 prevChi = chi; chi = oldChi + dt2_ * (instTemp / targetTemp_ - 1.0) / @@ -213,12 +185,10 @@ namespace OpenMD { atom = mol->nextFluctuatingCharge(j)) { cfrc = atom->getFlucQFrc(); - cvel =atom->getFlucQVel(); cmass = atom->getChargeMass(); // velocity half step - cvel = oldVel_[index] + dt2_/cmass*PhysicalConstants::energyConvert * cfrc - dt2_*chi*oldVel_[index]; - + cvel = oldVel_[index] + dt2_ * cfrc / cmass - dt2_*chi*oldVel_[index]; atom->setFlucQVel(cvel); ++index; } @@ -227,20 +197,19 @@ namespace OpenMD { break; } integralOfChidt += dt2_ * chi; - currentSnapshot_->setChiElectronic(chi); - currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt); + snap->setElectronicThermostat(make_pair(chi, integralOfChidt)); } void FluctuatingChargeNVT::resetPropagator() { if (!hasFlucQ_) return; - currentSnapshot_->setChiElectronic(0.0); - currentSnapshot_->setIntegralOfChiElectronicDt(0.0); + snap->setElectronicThermostat(make_pair(0.0, 0.0)); } RealType FluctuatingChargeNVT::calcConservedQuantity() { if (!hasFlucQ_) return 0.0; - RealType chi = currentSnapshot_->getChiElectronic(); - RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt(); + pair thermostat = snap->getElectronicThermostat(); + RealType chi = thermostat.first; + RealType integralOfChidt = thermostat.second; RealType fkBT = info_->getNFluctuatingCharges() * PhysicalConstants::kB *targetTemp_;