45#include "FluctuatingChargeNVT.hpp"
48#include "utils/Constants.hpp"
49#include "utils/simError.h"
53 FluctuatingChargeNVT::FluctuatingChargeNVT(
SimInfo* info) :
55 snap(info->getSnapshotManager()->getCurrentSnapshot()), thermo(info) {}
57 void FluctuatingChargeNVT::initialize() {
58 FluctuatingChargePropagator::initialize();
60 if (info_->getSimParams()->haveDt()) {
61 dt_ = info_->getSimParams()->getDt();
64 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
65 "FluctuatingChargeNVT Error: dt is not set\n");
70 if (!info_->getSimParams()->getUseIntialExtendedSystemState()) {
71 snap->setElectronicThermostat(make_pair(0.0, 0.0));
74 if (!fqParams_->haveTargetTemp()) {
75 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
76 "You can't use the FluctuatingChargeNVT "
77 "propagator without a flucQ.targetTemp!\n");
79 painCave.severity = OPENMD_ERROR;
82 targetTemp_ = fqParams_->getTargetTemp();
87 if (!fqParams_->haveTauThermostat()) {
88 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
89 "If you use the FluctuatingChargeNVT\n"
90 "\tpropagator, you must set flucQ.tauThermostat .\n");
92 painCave.severity = OPENMD_ERROR;
96 tauThermostat_ = fqParams_->getTauThermostat();
102 void FluctuatingChargeNVT::moveA() {
103 if (!hasFlucQ_)
return;
105 SimInfo::MoleculeIterator i;
106 Molecule::FluctuatingChargeIterator j;
109 RealType cvel, cpos, cfrc, cmass;
111 pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
112 RealType chi = thermostat.first;
113 RealType integralOfChidt = thermostat.second;
114 RealType instTemp = thermo.getElectronicTemperature();
118 for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
119 atom = mol->nextFluctuatingCharge(j)) {
120 cvel = atom->getFlucQVel();
121 cpos = atom->getFlucQPos();
122 cfrc = atom->getFlucQFrc();
123 cmass = atom->getChargeMass();
126 cvel += dt2_ * cfrc / cmass - dt2_ * chi * cvel;
130 atom->setFlucQVel(cvel);
131 atom->setFlucQPos(cpos);
135 chi += dt2_ * (instTemp / targetTemp_ - 1.0) /
136 (tauThermostat_ * tauThermostat_);
138 integralOfChidt += chi * dt2_;
139 snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
142 void FluctuatingChargeNVT::updateSizes() {
146 void FluctuatingChargeNVT::moveB() {
147 if (!hasFlucQ_)
return;
148 SimInfo::MoleculeIterator i;
149 Molecule::FluctuatingChargeIterator j;
153 pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
154 RealType chi = thermostat.first;
155 RealType oldChi = chi;
157 RealType integralOfChidt = thermostat.second;
159 RealType cfrc, cvel, cmass;
164 for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
165 atom = mol->nextFluctuatingCharge(j)) {
166 oldVel_[index] = atom->getFlucQVel();
173 for (
int k = 0; k < maxIterNum_; k++) {
175 instTemp = thermo.getElectronicTemperature();
178 chi = oldChi + dt2_ * (instTemp / targetTemp_ - 1.0) /
179 (tauThermostat_ * tauThermostat_);
183 for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
184 atom = mol->nextFluctuatingCharge(j)) {
185 cfrc = atom->getFlucQFrc();
186 cmass = atom->getChargeMass();
189 cvel = oldVel_[index] + dt2_ * cfrc / cmass -
190 dt2_ * chi * oldVel_[index];
191 atom->setFlucQVel(cvel);
195 if (fabs(prevChi - chi) <= chiTolerance_)
break;
197 integralOfChidt += dt2_ * chi;
198 snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
201 void FluctuatingChargeNVT::resetPropagator() {
202 if (!hasFlucQ_)
return;
203 snap->setElectronicThermostat(make_pair(0.0, 0.0));
206 RealType FluctuatingChargeNVT::calcConservedQuantity() {
207 if (!hasFlucQ_)
return 0.0;
208 pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
209 RealType chi = thermostat.first;
210 RealType integralOfChidt = thermostat.second;
214 RealType thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi *
215 chi / (2.0 * Constants::energyConvert);
217 RealType thermostat_potential =
218 fkBT * integralOfChidt / Constants::energyConvert;
220 return thermostat_kinetic + thermostat_potential;
abstract class for propagating fluctuating charge variables
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Molecule * beginMolecule(MoleculeIterator &i)
Returns the first molecule in this SimInfo and intialize the iterator.
int getNFluctuatingCharges()
Returns the total number of fluctuating charges that are present.
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.