54#include "brains/Thermo.hpp"
55#include "integrators/IntegratorCreator.hpp"
56#include "integrators/LangevinPiston.hpp"
58#include "utils/Constants.hpp"
59#include "utils/simError.h"
63 LangevinPiston::LangevinPiston(
SimInfo* info) :
NPT(info) {
65 NkBT = info_->getNGlobalIntegrableObjects() * Constants::kB * targetTemp;
69 W_ = 3.0 * NkBT * tb2;
72 if (!simParams->haveLangevinPistonDrag()) {
73 sprintf(painCave.errMsg,
"To use the LangevinPiston integrator, you must "
74 "set langevinPistonDrag "
76 painCave.severity = OPENMD_ERROR;
80 gamma_ = simParams->getLangevinPistonDrag();
86 randNumGen_ = info->getRandomNumberGenerator();
90 std::sqrt(2.0 * W_ * gamma_ * Constants::kB * targetTemp / dt);
92 forceDistribution_ = std::normal_distribution<RealType>(0.0, stdDev);
98 genRandomForce(randomForce_);
101 void LangevinPiston::moveA() {
102 SimInfo::MoleculeIterator i;
103 Molecule::IntegrableObjectIterator j;
116 instaTemp = thermo.getTemperature();
118 instaPress = Constants::pressureConvert *
119 (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
120 instaVol = thermo.getVolume();
130 for (sd = mol->beginIntegrableObject(j); sd != NULL;
131 sd = mol->nextIntegrableObject(j)) {
137 getVelScaleA(sc, vel);
141 vel += dt2 * Constants::energyConvert / mass * frc - dt2 * sc;
153 ji += dt2 * Constants::energyConvert * Tb;
155 rotAlgo_->rotate(sd, ji, dt);
169 for (sd = mol->beginIntegrableObject(j); sd != NULL;
170 sd = mol->nextIntegrableObject(j)) {
171 oldPos[index++] = sd->
getPos();
177 for (
int k = 0; k < maxIterNum_; k++) {
181 for (sd = mol->beginIntegrableObject(j); sd != NULL;
182 sd = mol->nextIntegrableObject(j)) {
186 this->getPosScale(pos, COM, index, sc);
188 pos = oldPos[index] + dt * (vel + sc);
195 rattle_->constraintA();
205 void LangevinPiston::moveB(
void) {
206 SimInfo::MoleculeIterator i;
207 Molecule::IntegrableObjectIterator j;
224 for (sd = mol->beginIntegrableObject(j); sd != NULL;
225 sd = mol->nextIntegrableObject(j)) {
226 oldVel[index] = sd->
getVel();
234 instaVol = thermo.getVolume();
236 for (
int k = 0; k < maxIterNum_; k++) {
237 instaTemp = thermo.getTemperature();
238 instaPress = thermo.getPressure();
242 this->calcVelScale();
247 for (sd = mol->beginIntegrableObject(j); sd != NULL;
248 sd = mol->nextIntegrableObject(j)) {
252 getVelScaleB(sc, index);
255 vel = oldVel[index] + dt2 * Constants::energyConvert / mass * frc -
264 ji = oldJi[index] + dt2 * Constants::energyConvert * Tb;
273 rattle_->constraintB();
274 if (this->etaConverged())
break;
281 void LangevinPiston::evolveEtaA() {
285 eta += dt2 * (instaVol * (instaPress - targetPressure) /
286 (Constants::pressureConvert * W_) -
287 gamma_ * eta + randomForce_ / W_);
291 void LangevinPiston::evolveEtaB() {
294 genRandomForce(randomForce_);
296 eta = oldEta + dt2 * (instaVol * (instaPress - targetPressure) /
297 (Constants::pressureConvert * W_) -
298 gamma_ * eta + randomForce_ / W_);
301 void LangevinPiston::calcVelScale() { vScale = eta; }
307 void LangevinPiston::getVelScaleB(
Vector3d& sc,
int index) {
308 sc = vScale * oldVel[index];
313 Vector3d rj = (oldPos[index] + pos) / (RealType)2.0 - COM;
317 void LangevinPiston::scaleSimBox() {
318 RealType scaleFactor;
321 scaleFactor = exp(dt * eta);
323 if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
324 sprintf(painCave.errMsg,
325 "LangevinPiston error: Attempting a Box scaling of more than 10 "
327 " check your tauBarostat, as it is probably too small!\n"
328 " eta = %lf, scaleFactor = %lf\n",
330 painCave.isFatal = 1;
339 bool LangevinPiston::etaConverged() {
340 return (fabs(prevEta - eta) <= etaTolerance);
343 void LangevinPiston::loadEta() {
344 Mat3x3d etaMat = snap->getBarostat();
348 void LangevinPiston::saveEta() {
353 snap->setBarostat(etaMat);
356 void LangevinPiston::genRandomForce(RealType& randomForce) {
358 if (worldRank == 0) {
360 randomForce = forceDistribution_(*randNumGen_);
365 MPI_Bcast(&randomForce, 1, MPI_REALTYPE, 0, MPI_COMM_WORLD);
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.
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
Mat3x3d getHmat()
Returns the H-Matrix.
void setHmat(const Mat3x3d &m)
Sets the H-Matrix.
"Don't move, or you're dead! Stand up! Captain, we've got them!"
Vector3d getTrq()
Returns the current torque of this stuntDouble.
Vector3d lab2Body(const Vector3d &v)
Converts a lab fixed vector to a body fixed vector.
Vector3d getVel()
Returns the current velocity of this stuntDouble.
RealType getMass()
Returns the mass of this stuntDouble.
Vector3d getPos()
Returns the current position of this stuntDouble.
void setPos(const Vector3d &pos)
Sets the current position of this stuntDouble.
void setVel(const Vector3d &vel)
Sets the current velocity of this stuntDouble.
Vector3d getJ()
Returns the current angular momentum of this stuntDouble (body -fixed).
bool isDirectional()
Tests if this stuntDouble is a directional one.
Vector3d getFrc()
Returns the current force of this stuntDouble.
void setJ(const Vector3d &angMom)
Sets the current angular momentum of this stuntDouble (body-fixed).
Mat3x3d getPressureTensor()
gives the pressure tensor in amu*fs^-2*Ang^-1
Vector3d getCom()
Returns the center of the mass of the whole system.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.