return couplingFrequency_;
}
+void AndersenTemperatureCoupling::updateReferenceTemperature(ArrayRef<const real> gmx_unused temperatures,
+ ReferenceTemperatureChangeAlgorithm gmx_unused algorithm)
+{
+ // Currently, we don't know about any temperature change algorithms, so we assert this never gets called
+ GMX_ASSERT(false, "AndersenTemperatureCoupling: Unknown ReferenceTemperatureChangeAlgorithm.");
+}
+
void AndersenTemperatureCoupling::elementSetup() {}
ISimulatorElement* AndersenTemperatureCoupling::getElementPointerImpl(
LegacySimulatorData* legacySimulatorData,
statePropagatorData,
legacySimulatorData->mdAtoms,
legacySimulatorData->cr);
+ auto* andersenThermostatPtr = andersenThermostat.get();
+ builderHelper->registerReferenceTemperatureUpdate(
+ [andersenThermostatPtr](ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm) {
+ andersenThermostatPtr->updateReferenceTemperature(temperatures, algorithm);
+ });
// T-coupling frequency will be composite element frequency
const auto frequency = andersenThermostat->frequency();
[[nodiscard]] int frequency() const;
private:
+ //! Update the reference temperature
+ static void updateReferenceTemperature(ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm);
+
//! Whether we're doing massive Andersen thermostatting
const bool doMassive_;
//! The rate
class LastStepSignaller;
class LoggingSignaller;
class NeighborSearchSignaller;
+enum class ReferenceTemperatureChangeAlgorithm;
enum class ScaleVelocities;
template<class Signaller>
class SignallerBuilder;
virtual DomDecCallback registerDomDecCallback() = 0;
};
+//! Callback updating the reference temperature
+using ReferenceTemperatureCallback =
+ std::function<void(ArrayRef<const real>, ReferenceTemperatureChangeAlgorithm algorithm)>;
+
//! /}
} // namespace gmx
initialVolume,
legacySimulatorData->inputrec->compress,
statePropagatorData));
- const auto* ptrToDataObject = builderHelper->simulationData<MttkData>(MttkData::dataID()).value();
+ auto* ptrToDataObject = builderHelper->simulationData<MttkData>(MttkData::dataID()).value();
energyData->addConservedEnergyContribution([ptrToDataObject](Step /*unused*/, Time time) {
return ptrToDataObject->temperatureCouplingIntegral(time);
});
energyData->setParrinelloRahmanBoxVelocities(
[ptrToDataObject]() { return ptrToDataObject->boxVelocity_; });
+ builderHelper->registerReferenceTemperatureUpdate(
+ [ptrToDataObject](ArrayRef<const real> temperatures, ReferenceTemperatureChangeAlgorithm algorithm) {
+ ptrToDataObject->updateReferenceTemperature(temperatures[0], algorithm);
+ });
}
std::string MttkData::dataID()
integralTime_(0.0),
referencePressure_(referencePressure),
boxVelocity_{ { 0 } },
+ referenceTemperature_(referenceTemperature),
statePropagatorData_(statePropagatorData)
{
}
return boxVelocity_;
}
+void MttkData::updateReferenceTemperature(real temperature,
+ ReferenceTemperatureChangeAlgorithm gmx_unused algorithm)
+{
+ // Currently, we don't know about any temperature change algorithms, so we assert this never gets called
+ GMX_ASSERT(false, "MttkData: Unknown ReferenceTemperatureChangeAlgorithm.");
+ invMass_ *= temperature / referenceTemperature_;
+ referenceTemperature_ = temperature;
+}
+
namespace
{
/*!
double temperatureCouplingIntegral(Time time) const;
//! Calculate the current value of the MTTK conserved energy if it is needed
void calculateIntegralIfNeeded();
+ //! Update the reference temperature
+ void updateReferenceTemperature(real temperature, ReferenceTemperatureChangeAlgorithm algorithm);
//! The coupling time step
const real couplingTimeStep_;
const real referencePressure_;
//! The current box velocities (used for reporting only)
tensor boxVelocity_;
+ //! The reference temperature the mass is based on
+ real referenceTemperature_;
// TODO: Clarify relationship to data objects and find a more robust alternative to raw pointers (#3583)
//! Pointer to the micro state data (access to the current box)
//! Return the current time of the NHC integral for the group
real integralTime() const;
+ //! Set the reference temperature
+ void updateReferenceTemperature(real temperature);
+
private:
//! Increment coordinate time and update integral if applicable
void finalizeUpdate(real couplingTimeStep);
//! The reference temperature of this group
- const real referenceTemperature_;
+ real referenceTemperature_;
+ //! The coupling time of this group
+ const real couplingTime_;
//! The number of degrees of freedom in this group
const real numDegreesOfFreedom_;
//! The chain length of this group
real couplingTimeStep,
NhcUsage nhcUsage) :
referenceTemperature_(referenceTemperature),
+ couplingTime_(couplingTime),
numDegreesOfFreedom_(numDegreesOfFreedom),
chainLength_(chainLength),
couplingTimeStep_(couplingTimeStep),
ArrayRef<real>(),
nhcUsage));
}
+ auto* nhcDataPtr =
+ builderHelper
+ ->simulationData<NoseHooverChainsData>(NoseHooverChainsData::dataID(nhcUsage))
+ .value();
+ builderHelper->registerReferenceTemperatureUpdate(
+ [nhcDataPtr](ArrayRef<const real> temperatures, ReferenceTemperatureChangeAlgorithm algorithm) {
+ nhcDataPtr->updateReferenceTemperature(temperatures, algorithm);
+ });
const auto* ptrToDataObject =
builderHelper
return timesClose(std::lround(coordinateTime_ / couplingTimeStep_) * couplingTimeStep_, coordinateTime_);
}
+void NoseHooverChainsData::updateReferenceTemperature(ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm gmx_unused algorithm)
+{
+ // Currently, we don't know about any temperature change algorithms, so we assert this never gets called
+ GMX_ASSERT(false, "NoseHooverChainsData: Unknown ReferenceTemperatureChangeAlgorithm.");
+ for (auto temperatureGroup = 0; temperatureGroup < numTemperatureGroups_; ++temperatureGroup)
+ {
+ noseHooverGroups_[temperatureGroup].updateReferenceTemperature(temperatures[temperatureGroup]);
+ if (noseHooverGroups_[temperatureGroup].isAtFullCouplingTimeStep())
+ {
+ noseHooverGroups_[temperatureGroup].calculateIntegral();
+ }
+ }
+}
+
+void NoseHooverGroup::updateReferenceTemperature(real temperature)
+{
+ const bool newTemperatureIsValid = (temperature > 0 && couplingTime_ > 0 && numDegreesOfFreedom_ > 0);
+ const bool oldTemperatureIsValid =
+ (referenceTemperature_ > 0 && couplingTime_ > 0 && numDegreesOfFreedom_ > 0);
+ GMX_RELEASE_ASSERT(newTemperatureIsValid == oldTemperatureIsValid,
+ "Cannot turn temperature coupling on / off during simulation run.");
+ if (oldTemperatureIsValid && newTemperatureIsValid)
+ {
+ const real velocityFactor = std::sqrt(temperature / referenceTemperature_);
+ for (auto chainPosition = 0; chainPosition < chainLength_; ++chainPosition)
+ {
+ invXiMass_[chainPosition] *= (referenceTemperature_ / temperature);
+ xiVelocities_[chainPosition] *= velocityFactor;
+ }
+ }
+ referenceTemperature_ = temperature;
+ if (isAtFullCouplingTimeStep())
+ {
+ calculateIntegral();
+ }
+}
+
namespace
{
/*!
private:
//! Return the value of the coupling integral at a specific time
double temperatureCouplingIntegral(Time time) const;
+ //! Update the reference temperature
+ void updateReferenceTemperature(ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm);
//! CheckpointHelper identifier
const std::string identifier_;
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief Defines the a helper struct managing reference temperature changes
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \ingroup module_modularsimulator
+ */
+
+#include "gmxpre.h"
+
+#include "referencetemperaturemanager.h"
+
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/utility/gmxassert.h"
+
+namespace gmx
+{
+
+ReferenceTemperatureManager::ReferenceTemperatureManager(t_inputrec* inputrec) : inputrec_(inputrec)
+{
+}
+
+void ReferenceTemperatureManager::registerUpdateCallback(ReferenceTemperatureCallback referenceTemperatureCallback)
+{
+ callbacks_.emplace_back(std::move(referenceTemperatureCallback));
+}
+
+void ReferenceTemperatureManager::setReferenceTemperature(ArrayRef<const real> newReferenceTemperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm)
+{
+ GMX_RELEASE_ASSERT(newReferenceTemperatures.ssize() == inputrec_->opts.ngtc,
+ "Expected one new reference temperature per temperature group.");
+
+ std::copy(newReferenceTemperatures.begin(), newReferenceTemperatures.end(), inputrec_->opts.ref_t);
+ for (const auto& callback : callbacks_)
+ {
+ callback(constArrayRefFromArray(inputrec_->opts.ref_t, inputrec_->opts.ngtc), algorithm);
+ }
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief Declares the a helper struct managing reference temperature changes
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \ingroup module_modularsimulator
+ *
+ * This header is only used within the modular simulator module
+ */
+
+#ifndef GMX_MODULARSIMULATOR_REFERENCETEMPERATUREMANAGER_H
+#define GMX_MODULARSIMULATOR_REFERENCETEMPERATUREMANAGER_H
+
+#include <functional>
+
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/real.h"
+
+#include "modularsimulatorinterfaces.h"
+
+struct t_inputrec;
+
+namespace gmx
+{
+
+/*! \internal
+ * \brief The algorithm changing the reference temperature
+ *
+ * In the legacy implementation, reference temperature changes by
+ * different algorithms are not handled identically. This enum is
+ * used to inform clients what algorithm caused the temperature
+ * change, allowing them to customize their response.
+ */
+enum class ReferenceTemperatureChangeAlgorithm
+{
+};
+
+/*! \internal
+ * \brief Object managing reference temperature changes
+ *
+ * The ReferenceTemperatureManager allows to change the reference
+ * temperatures of the temperature groups. Elements can register a callback
+ * if they need to be informed about changes.
+ *
+ * The ReferenceTemperatureManager updates the inputrec. Elements
+ * might, however, have a copy of the reference temperature they
+ * need updated, or perform another action upon change of the
+ * reference temperature (e.g. velocity scaling or recalculating
+ * a temperature coupling integral).
+ */
+class ReferenceTemperatureManager final
+{
+public:
+ //! Constructor
+ ReferenceTemperatureManager(t_inputrec* inputrec);
+ //! Register a callback for reference temperature update
+ void registerUpdateCallback(ReferenceTemperatureCallback referenceTemperatureCallback);
+ //! Set reference temperatures (one per temperature group)
+ void setReferenceTemperature(ArrayRef<const real> newReferenceTemperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm);
+
+private:
+ //! List of callbacks
+ std::vector<ReferenceTemperatureCallback> callbacks_;
+ //! Pointer to the input record
+ t_inputrec* inputrec_;
+};
+
+} // namespace gmx
+
+#endif // GMX_MODULARSIMULATOR_REFERENCETEMPERATUREMANAGER_H
#include "modularsimulator.h"
#include "pmeloadbalancehelper.h"
#include "propagator.h"
+#include "referencetemperaturemanager.h"
#include "statepropagatordata.h"
namespace gmx
legacySimulatorData->startingBehavior,
simulationsShareState);
registerExistingElement(energyData_->element());
+
+ // This is the only modular simulator object which changes the inputrec
+ // TODO: Avoid changing inputrec (#3854)
+ storeSimulationData(
+ "ReferenceTemperatureManager",
+ ReferenceTemperatureManager(const_cast<t_inputrec*>(legacySimulatorData->inputrec)));
+ auto* referenceTemperatureManager =
+ simulationData<ReferenceTemperatureManager>("ReferenceTemperatureManager").value();
+
+ // State propagator data is scaling velocities if reference temperature is updated
+ auto* statePropagatorDataPtr = statePropagatorData_.get();
+ referenceTemperatureManager->registerUpdateCallback(
+ [statePropagatorDataPtr](ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm) {
+ statePropagatorDataPtr->updateReferenceTemperature(temperatures, algorithm);
+ });
}
ModularSimulatorAlgorithm ModularSimulatorAlgorithmBuilder::build()
builder_->propagatorConnections_.emplace_back(std::move(connectionData));
}
+void ModularSimulatorAlgorithmBuilderHelper::registerReferenceTemperatureUpdate(
+ ReferenceTemperatureCallback referenceTemperatureCallback)
+{
+ auto* referenceTemperatureManager =
+ simulationData<ReferenceTemperatureManager>("ReferenceTemperatureManager").value();
+ referenceTemperatureManager->registerUpdateCallback(std::move(referenceTemperatureCallback));
+}
+
+ReferenceTemperatureCallback ModularSimulatorAlgorithmBuilderHelper::changeReferenceTemperatureCallback()
+{
+ // Capture is safe because SimulatorAlgorithm will manage life time of both the
+ // recipient of the callback and the reference temperature manager
+ auto* referenceTemperatureManager =
+ simulationData<ReferenceTemperatureManager>("ReferenceTemperatureManager").value();
+ return [referenceTemperatureManager](ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm) {
+ referenceTemperatureManager->setReferenceTemperature(temperatures, algorithm);
+ };
+}
+
} // namespace gmx
void registerTemperaturePressureControl(std::function<void(const PropagatorConnection&)> registrationFunction);
//! Register a propagator to be used with a temperature / pressure control algorithm
void registerPropagator(PropagatorConnection connectionData);
+ //! Register for callback after an update to the reference temperature
+ void registerReferenceTemperatureUpdate(ReferenceTemperatureCallback referenceTemperatureCallback);
+ //! Get a callback to change reference temperature
+ ReferenceTemperatureCallback changeReferenceTemperatureCallback();
private:
//! Pointer to the associated ModularSimulatorAlgorithmBuilder
namespace gmx
{
+/*! \internal
+ * \brief Helper object to scale velocities according to reference temperature change
+ */
+class StatePropagatorData::ReferenceTemperatureHelper
+{
+public:
+ //! Constructor
+ ReferenceTemperatureHelper(const t_inputrec* inputrec,
+ StatePropagatorData* statePropagatorData,
+ const t_mdatoms* mdatoms) :
+ numTemperatureGroups_(inputrec->opts.ngtc),
+ referenceTemperature_(inputrec->opts.ref_t, inputrec->opts.ref_t + inputrec->opts.ngtc),
+ velocityScalingFactors_(numTemperatureGroups_),
+ statePropagatorData_(statePropagatorData),
+ mdatoms_(mdatoms)
+ {
+ }
+
+ /*! \brief Update the reference temperature
+ *
+ * Changing the reference temperature requires scaling the velocities, which
+ * is done here.
+ *
+ * \param temperatures New reference temperatures
+ * \param algorithm The algorithm which initiated the temperature update
+ */
+ void updateReferenceTemperature(ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm gmx_unused algorithm)
+ {
+ // Currently, we don't know about any temperature change algorithms, so we assert this never gets called
+ GMX_ASSERT(false, "StatePropagatorData: Unknown ReferenceTemperatureChangeAlgorithm.");
+ for (int temperatureGroup = 0; temperatureGroup < numTemperatureGroups_; ++temperatureGroup)
+ {
+ velocityScalingFactors_[temperatureGroup] =
+ std::sqrt(temperatures[temperatureGroup] / referenceTemperature_[temperatureGroup]);
+ }
+
+ auto velocities = statePropagatorData_->velocitiesView().unpaddedArrayRef();
+ int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
+#pragma omp parallel for num_threads(nth) schedule(static) default(none) \
+ shared(nth, velocityScalingFactors_, velocities)
+ for (int threadIndex = 0; threadIndex < nth; threadIndex++)
+ {
+ int startAtom = 0;
+ int endAtom = 0;
+ getThreadAtomRange(nth, threadIndex, statePropagatorData_->localNAtoms_, &startAtom, &endAtom);
+ for (int atomIdx = startAtom; atomIdx < endAtom; ++atomIdx)
+ {
+ const int temperatureGroup = (mdatoms_->cTC != nullptr) ? mdatoms_->cTC[atomIdx] : 0;
+ velocities[atomIdx] *= velocityScalingFactors_[temperatureGroup];
+ }
+ }
+ std::copy(temperatures.begin(), temperatures.end(), referenceTemperature_.begin());
+ }
+
+private:
+ //! The number of temperature groups
+ const int numTemperatureGroups_;
+ //! Coupling temperature per group
+ std::vector<real> referenceTemperature_;
+ //! Working array used at every temperature update
+ std::vector<real> velocityScalingFactors_;
+
+ //! Pointer to StatePropagatorData to scale velocities
+ StatePropagatorData* statePropagatorData_;
+ //! Atom parameters for this domain (temperature group information)
+ const t_mdatoms* mdatoms_;
+};
+
StatePropagatorData::StatePropagatorData(int numAtoms,
FILE* fplog,
const t_commrec* cr,
finalConfigurationFilename,
inputrec,
globalTop)),
+ referenceTemperatureHelper_(std::make_unique<ReferenceTemperatureHelper>(inputrec, this, mdatoms)),
vvResetVelocities_(false),
isRegularSimulationEnd_(false),
lastStep_(-1),
}
}
+StatePropagatorData::~StatePropagatorData() = default;
+
StatePropagatorData::Element* StatePropagatorData::element()
{
return element_.get();
}
}
+void StatePropagatorData::updateReferenceTemperature(ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm)
+{
+ referenceTemperatureHelper_->updateReferenceTemperature(temperatures, algorithm);
+}
+
void StatePropagatorData::Element::scheduleTask(Step step,
Time gmx_unused time,
const RegisterRunFunction& registerRunFunction)
const t_mdatoms* mdatoms,
const gmx_mtop_t& globalTop);
+ //! Destructor (allows forward declaration of internal type)
+ ~StatePropagatorData();
+
// Allow access to state
//! Get write access to position vector
ArrayRefWithPadding<RVec> positionsView();
//! Initial set up for the associated element
void setup();
+ //! Update the reference temperature
+ void updateReferenceTemperature(ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm);
+ //! Helper class handling reference temperature change
+ class ReferenceTemperatureHelper;
+
//! Read everything that can be stored in t_trxframe from a checkpoint file
static void readCheckpointToTrxFrame(t_trxframe* trxFrame, ReadCheckpointData readCheckpointData);
//! CheckpointHelper identifier
//! The element
std::unique_ptr<Element> element_;
+ //! Instance of reference temperature helper
+ std::unique_ptr<ReferenceTemperatureHelper> referenceTemperatureHelper_;
//! Move x_ to previousX_
void copyPosition();
//! Read private data from checkpoint
virtual void readCheckpoint(std::optional<ReadCheckpointData> checkpointData, const t_commrec* cr) = 0;
+ //! Update the reference temperature and update and return the temperature coupling integral
+ virtual real updateReferenceTemperatureAndIntegral(int temperatureGroup,
+ real newTemperature,
+ ReferenceTemperatureChangeAlgorithm algorithm,
+ const TemperatureCouplingData& temperatureCouplingData) = 0;
+
//! Standard virtual destructor
virtual ~ITemperatureCouplingImpl() = default;
};
{
}
+ //! No changes needed
+ real updateReferenceTemperatureAndIntegral(int temperatureGroup,
+ real gmx_unused newTemperature,
+ ReferenceTemperatureChangeAlgorithm gmx_unused algorithm,
+ const TemperatureCouplingData& temperatureCouplingData) override
+ {
+ return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup];
+ }
+
//! Constructor
VRescaleTemperatureCoupling(int64_t seed) : seed_(seed) {}
{
}
+ //! No changes needed
+ real updateReferenceTemperatureAndIntegral(int temperatureGroup,
+ real gmx_unused newTemperature,
+ ReferenceTemperatureChangeAlgorithm gmx_unused algorithm,
+ const TemperatureCouplingData& temperatureCouplingData) override
+ {
+ return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup];
+ }
+
private:
//! View on the scaling factor of the propagator (pre-step velocities)
ArrayRef<real> lambdaStartVelocities_;
class NoseHooverTemperatureCoupling final : public ITemperatureCouplingImpl
{
public:
+ //! Calculate the current value of the temperature coupling integral
+ real integral(int temperatureGroup, real numDegreesOfFreedom, real referenceTemperature)
+ {
+ return 0.5 * c_boltz * numDegreesOfFreedom
+ * (xiVelocities_[temperatureGroup] * xiVelocities_[temperatureGroup])
+ / invXiMass_[temperatureGroup]
+ + numDegreesOfFreedom * xi_[temperatureGroup] * c_boltz * referenceTemperature;
+ }
+
//! Apply the Nose-Hoover temperature control
real apply(Step gmx_unused step,
int temperatureGroup,
1. / (1 + 0.5 * thermostatData.couplingTimeStep * xiVelocities_[temperatureGroup]);
// Current value of the thermostat integral
- return 0.5 * c_boltz * thermostatData.numDegreesOfFreedom[temperatureGroup]
- * (xiVelocities_[temperatureGroup] * xiVelocities_[temperatureGroup])
- / invXiMass_[temperatureGroup]
- + thermostatData.numDegreesOfFreedom[temperatureGroup] * xi_[temperatureGroup]
- * c_boltz * thermostatData.referenceTemperature[temperatureGroup];
+ return integral(temperatureGroup,
+ thermostatData.numDegreesOfFreedom[temperatureGroup],
+ thermostatData.referenceTemperature[temperatureGroup]);
}
//! Connect with propagator - Nose-Hoover scales start and end step velocities
}
}
+ //! Adapt masses
+ real updateReferenceTemperatureAndIntegral(int temperatureGroup,
+ real gmx_unused newTemperature,
+ ReferenceTemperatureChangeAlgorithm gmx_unused algorithm,
+ const TemperatureCouplingData& temperatureCouplingData) override
+ {
+ // Currently, we don't know about any temperature change algorithms, so we assert this never gets called
+ GMX_ASSERT(false,
+ "NoseHooverTemperatureCoupling: Unknown ReferenceTemperatureChangeAlgorithm.");
+ const bool newTemperatureIsValid =
+ (newTemperature > 0 && temperatureCouplingData.couplingTime[temperatureGroup] > 0
+ && temperatureCouplingData.numDegreesOfFreedom[temperatureGroup] > 0);
+ const bool oldTemperatureIsValid =
+ (temperatureCouplingData.referenceTemperature[temperatureGroup] > 0
+ && temperatureCouplingData.couplingTime[temperatureGroup] > 0
+ && temperatureCouplingData.numDegreesOfFreedom[temperatureGroup] > 0);
+ GMX_RELEASE_ASSERT(newTemperatureIsValid == oldTemperatureIsValid,
+ "Cannot turn temperature coupling on / off during simulation run.");
+ if (oldTemperatureIsValid && newTemperatureIsValid)
+ {
+ invXiMass_[temperatureGroup] *=
+ (temperatureCouplingData.referenceTemperature[temperatureGroup] / newTemperature);
+ xiVelocities_[temperatureGroup] *= std::sqrt(
+ newTemperature / temperatureCouplingData.referenceTemperature[temperatureGroup]);
+ }
+ return integral(temperatureGroup,
+ temperatureCouplingData.numDegreesOfFreedom[temperatureGroup],
+ newTemperature);
+ }
+
private:
//! The thermostat degree of freedom
std::vector<real> xi_;
}
}
+void VelocityScalingTemperatureCoupling::updateReferenceTemperature(ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm)
+{
+ TemperatureCouplingData thermostatData = {
+ couplingTimeStep_, referenceTemperature_, couplingTime_, numDegreesOfFreedom_, temperatureCouplingIntegral_
+ };
+ for (int temperatureGroup = 0; (temperatureGroup < numTemperatureGroups_); temperatureGroup++)
+ {
+ temperatureCouplingIntegral_[temperatureGroup] =
+ temperatureCouplingImpl_->updateReferenceTemperatureAndIntegral(
+ temperatureGroup, temperatures[temperatureGroup], algorithm, thermostatData);
+ }
+ // Currently, we don't know about any temperature change algorithms, so we assert this never gets called
+ GMX_ASSERT(false,
+ "VelocityScalingTemperatureCoupling: Unknown ReferenceTemperatureChangeAlgorithm.");
+ std::copy(temperatures.begin(), temperatures.end(), referenceTemperature_.begin());
+}
+
namespace
{
/*!
[thermostat, propagatorTag](const PropagatorConnection& connection) {
thermostat->connectWithMatchingPropagator(connection, propagatorTag);
});
+ builderHelper->registerReferenceTemperatureUpdate(
+ [thermostat](ArrayRef<const real> temperatures, ReferenceTemperatureChangeAlgorithm algorithm) {
+ thermostat->updateReferenceTemperature(temperatures, algorithm);
+ });
return element;
}
const PropagatorTag& propagatorTag);
private:
+ //! Update the reference temperature
+ void updateReferenceTemperature(ArrayRef<const real> temperatures,
+ ReferenceTemperatureChangeAlgorithm algorithm);
+
//! The frequency at which the thermostat is applied
const int nstcouple_;
//! If != 0, offset the step at which the thermostat is applied
//! The coupling time step - simulation time step x nstcouple_
const double couplingTimeStep_;
//! Coupling temperature per group
- const std::vector<real> referenceTemperature_;
+ std::vector<real> referenceTemperature_;
//! Coupling time per group
const std::vector<real> couplingTime_;
//! Number of degrees of freedom per group
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff