--- /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 Andersen temperature coupling for the modular simulator
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \ingroup module_modularsimulator
+ */
+
+#include "andersentemperaturecoupling.h"
+
+#include <vector>
+
+#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/math/functions.h"
+#include "gromacs/math/units.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/mdatoms.h"
+#include "gromacs/mdlib/stat.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/mdrun/isimulator.h"
+#include "gromacs/random/tabulatednormaldistribution.h"
+#include "gromacs/random/threefry.h"
+#include "gromacs/random/uniformrealdistribution.h"
+
+#include "compositesimulatorelement.h"
+#include "constraintelement.h"
+#include "simulatoralgorithm.h"
+#include "statepropagatordata.h"
+
+namespace gmx
+{
+AndersenTemperatureCoupling::AndersenTemperatureCoupling(double simulationTimestep,
+ bool doMassive,
+ int64_t seed,
+ ArrayRef<const real> referenceTemperature,
+ ArrayRef<const real> couplingTime,
+ StatePropagatorData* statePropagatorData,
+ const MDAtoms* mdAtoms,
+ const t_commrec* cr) :
+ doMassive_(doMassive),
+ randomizationRate_(simulationTimestep / couplingTime[0]),
+ couplingFrequency_(doMassive ? roundToInt(1. / randomizationRate_) : 1),
+ seed_(seed),
+ referenceTemperature_(referenceTemperature),
+ couplingTime_(couplingTime),
+ statePropagatorData_(statePropagatorData),
+ mdAtoms_(mdAtoms->mdatoms()),
+ cr_(cr)
+{
+}
+
+void AndersenTemperatureCoupling::scheduleTask(Step step,
+ Time gmx_unused time,
+ const RegisterRunFunction& registerRunFunction)
+{
+ if (do_per_step(step, couplingFrequency_))
+ {
+ registerRunFunction([this, step]() { apply(step); });
+ }
+}
+
+void AndersenTemperatureCoupling::apply(Step step)
+{
+ ThreeFry2x64<0> rng(seed_, RandomDomain::Thermostat);
+ UniformRealDistribution<real> uniformDist;
+ TabulatedNormalDistribution<real, 14> normalDist;
+
+ const bool doDomainDecomposition = DOMAINDECOMP(cr_);
+
+ auto velocities = statePropagatorData_->velocitiesView().unpaddedArrayRef();
+
+ for (int atomIdx = 0; atomIdx < mdAtoms_->homenr; ++atomIdx)
+ {
+ const int temperatureGroup = mdAtoms_->cTC ? mdAtoms_->cTC[atomIdx] : 0;
+ if (referenceTemperature_[temperatureGroup] <= 0 || couplingTime_[temperatureGroup] <= 0)
+ {
+ continue;
+ }
+
+ const int globalAtomIdx = doDomainDecomposition ? cr_->dd->globalAtomIndices[atomIdx] : atomIdx;
+ rng.restart(step, globalAtomIdx);
+
+ // For massive Andersen, this function is only called periodically, but we apply each time
+ // Otherwise, this function is called every step, but we randomize atoms probabilistically
+ if (!doMassive_)
+ {
+ uniformDist.reset();
+ }
+ if (doMassive_ || (uniformDist(rng) < randomizationRate_))
+ {
+ const real scalingFactor = std::sqrt(c_boltz * referenceTemperature_[temperatureGroup]
+ * mdAtoms_->invmass[atomIdx]);
+ normalDist.reset();
+ for (int d = 0; d < DIM; d++)
+ {
+ velocities[atomIdx][d] = scalingFactor * normalDist(rng);
+ }
+ }
+ }
+}
+
+int AndersenTemperatureCoupling::frequency() const
+{
+ return couplingFrequency_;
+}
+
+void AndersenTemperatureCoupling::elementSetup() {}
+ISimulatorElement* AndersenTemperatureCoupling::getElementPointerImpl(
+ LegacySimulatorData* legacySimulatorData,
+ ModularSimulatorAlgorithmBuilderHelper* builderHelper,
+ StatePropagatorData* statePropagatorData,
+ EnergyData* energyData,
+ FreeEnergyPerturbationData* freeEnergyPerturbationData,
+ GlobalCommunicationHelper gmx_unused* globalCommunicationHelper)
+{
+ GMX_RELEASE_ASSERT(legacySimulatorData->inputrec->etc == TemperatureCoupling::Andersen
+ || legacySimulatorData->inputrec->etc == TemperatureCoupling::AndersenMassive,
+ "Expected the thermostat type to be andersen or andersen-massive.");
+ auto andersenThermostat = std::make_unique<AndersenTemperatureCoupling>(
+ legacySimulatorData->inputrec->delta_t,
+ legacySimulatorData->inputrec->etc == TemperatureCoupling::AndersenMassive,
+ legacySimulatorData->inputrec->andersen_seed,
+ constArrayRefFromArray(legacySimulatorData->inputrec->opts.ref_t,
+ legacySimulatorData->inputrec->opts.ngtc),
+ constArrayRefFromArray(legacySimulatorData->inputrec->opts.tau_t,
+ legacySimulatorData->inputrec->opts.ngtc),
+ statePropagatorData,
+ legacySimulatorData->mdAtoms,
+ legacySimulatorData->cr);
+
+ // T-coupling frequency will be composite element frequency
+ const auto frequency = andersenThermostat->frequency();
+ // Set up call list for composite element
+ std::vector<compat::not_null<ISimulatorElement*>> elementCallList = { compat::make_not_null(
+ andersenThermostat.get()) };
+ // Set up element list for composite element
+ std::vector<std::unique_ptr<gmx::ISimulatorElement>> elements;
+ elements.emplace_back(std::move(andersenThermostat));
+
+ // If there are constraints, add constraint element after Andersen element
+ if (legacySimulatorData->constr)
+ {
+ // This is excluded in preprocessing -
+ // asserted here to make sure things don't get out of sync
+ GMX_RELEASE_ASSERT(
+ legacySimulatorData->inputrec->etc == TemperatureCoupling::AndersenMassive,
+ "Per-particle Andersen thermostat is not implemented for systems with constrains.");
+ // Build constraint element
+ auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Velocities>>(
+ legacySimulatorData->constr,
+ statePropagatorData,
+ energyData,
+ freeEnergyPerturbationData,
+ MASTER(legacySimulatorData->cr),
+ legacySimulatorData->fplog,
+ legacySimulatorData->inputrec,
+ legacySimulatorData->mdAtoms->mdatoms());
+ // Add call to composite element call list
+ elementCallList.emplace_back(compat::make_not_null(constraintElement.get()));
+ // Move ownership of constraint element to composite element
+ elements.emplace_back(std::move(constraintElement));
+ }
+
+ // Store composite element in builder helper and return pointer
+ return builderHelper->storeElement(std::make_unique<CompositeSimulatorElement>(
+ std::move(elementCallList), std::move(elements), frequency));
+}
+
+} // 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 Andersen temperature coupling for the modular simulator
+ *
+ * \author Pascal Merz <pascal.merz@me.com>
+ * \ingroup module_modularsimulator
+ */
+
+#ifndef GMX_MODULARSIMULATOR_ANDERSENTHERMOSTAT_H
+#define GMX_MODULARSIMULATOR_ANDERSENTHERMOSTAT_H
+
+#include "gromacs/utility/arrayref.h"
+
+#include "energydata.h"
+#include "modularsimulatorinterfaces.h"
+#include "propagator.h"
+
+struct t_commrec;
+struct t_mdatoms;
+
+namespace gmx
+{
+
+/*! \internal
+ * \ingroup module_modularsimulator
+ * \brief Element implementing the Andersen thermostat
+ *
+ */
+class AndersenTemperatureCoupling final : public ISimulatorElement
+{
+public:
+ //! Constructor
+ AndersenTemperatureCoupling(double simulationTimestep,
+ bool doMassive,
+ int64_t seed,
+ ArrayRef<const real> referenceTemperature,
+ ArrayRef<const real> couplingTime,
+ StatePropagatorData* statePropagatorData,
+ const MDAtoms* mdAtoms,
+ const t_commrec* cr);
+
+ /*! \brief Register run function for step / time
+ *
+ * \param step The step number
+ * \param time The time
+ * \param registerRunFunction Function allowing to register a run function
+ */
+ void scheduleTask(Step step, Time time, const RegisterRunFunction& registerRunFunction) override;
+
+ //! Sanity check at setup time
+ void elementSetup() override;
+ //! No element teardown needed
+ void elementTeardown() override {}
+
+ /*! \brief Factory method implementation
+ *
+ * \param legacySimulatorData Pointer allowing access to simulator level data
+ * \param builderHelper ModularSimulatorAlgorithmBuilder helper object
+ * \param statePropagatorData Pointer to the \c StatePropagatorData object
+ * \param energyData Pointer to the \c EnergyData object
+ * \param freeEnergyPerturbationData Pointer to the \c FreeEnergyPerturbationData object
+ * \param globalCommunicationHelper Pointer to the \c GlobalCommunicationHelper object
+ *
+ * \return Pointer to the element to be added. Element needs to have been stored using \c storeElement
+ */
+ static ISimulatorElement* getElementPointerImpl(LegacySimulatorData* legacySimulatorData,
+ ModularSimulatorAlgorithmBuilderHelper* builderHelper,
+ StatePropagatorData* statePropagatorData,
+ EnergyData* energyData,
+ FreeEnergyPerturbationData* freeEnergyPerturbationData,
+ GlobalCommunicationHelper* globalCommunicationHelper);
+
+ //! Returns the frequency at which temperature coupling is performed
+ [[nodiscard]] int frequency() const;
+
+private:
+ //! Whether we're doing massive Andersen thermostatting
+ const bool doMassive_;
+ //! The rate
+ const real randomizationRate_;
+ //! The frequency at which the thermostat is applied
+ const int couplingFrequency_;
+ //! The random seed
+ const int64_t seed_;
+ //! Coupling temperature per group
+ ArrayRef<const real> referenceTemperature_;
+ //! Coupling time per group
+ ArrayRef<const real> couplingTime_;
+
+ // TODO: Clarify relationship to data objects and find a more robust alternative to raw pointers (#3583)
+ //! Pointer to the micro state
+ StatePropagatorData* statePropagatorData_;
+ //! Atom parameters for this domain.
+ const t_mdatoms* mdAtoms_;
+ //! Handles communication.
+ const t_commrec* cr_;
+
+ //! Apply the thermostat at step
+ void apply(Step step);
+};
+
+} // namespace gmx
+
+#endif // GMX_MODULARSIMULATOR_ANDERSENTHERMOSTAT_H
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/int64_to_int.h"
+#include "andersentemperaturecoupling.h"
#include "computeglobalselement.h"
#include "constraintelement.h"
#include "firstorderpressurecoupling.h"
ReportPreviousStepConservedEnergy::Yes,
PropagatorTag("VelocityHalfAndPositionFullStep"));
}
+ else if (ETC_ANDERSEN(legacySimulatorData_->inputrec->etc))
+ {
+ builder->add<AndersenTemperatureCoupling>();
+ }
builder->add<Propagator<IntegrationStage::VelocityVerletPositionsAndVelocities>>(
PropagatorTag("VelocityHalfAndPositionFullStep"),
TimeStep(legacySimulatorData_->inputrec->delta_t));
|| inputrec->etc == TemperatureCoupling::VRescale
|| inputrec->etc == TemperatureCoupling::Berendsen
|| (inputrec->etc == TemperatureCoupling::NoseHoover
- && inputrec->eI == IntegrationAlgorithm::MD),
- "Only v-rescale, Berendsen and Nose-Hoover "
+ && inputrec->eI == IntegrationAlgorithm::MD)
+ || ETC_ANDERSEN(inputrec->etc),
+ "Only v-rescale, Berendsen, Nose-Hoover, "
+ "and Andersen / Andersen-massive "
"thermostats are supported by the modular simulator.");
isInputCompatible = isInputCompatible
&& conditionalAssert(inputrec->epc == PressureCoupling::No
int maxNumWarnings = 0;
+ const bool isAndersen = (tcoupling == "andersen-massive" || tcoupling == "andersen");
+ const bool hasConstraints = (simulationName != "argon12");
+
// TODO At some point we should also test PME-only ranks.
const int numRanksAvailable = getNumberOfTestMpiRanks();
if (!isNumberOfPpRanksSupported(simulationName, numRanksAvailable))
return;
}
+ if (isAndersen && pcoupling == "berendsen")
+ {
+ // "Using Berendsen pressure coupling invalidates the true ensemble for the thermostat"
+ maxNumWarnings++;
+ }
+ if (tcoupling == "andersen" && hasConstraints)
+ {
+ // Constraints are not allowed with non-massive Andersen
+ return;
+ }
+
if (tcoupling == "nose-hoover" && pcoupling == "berendsen")
{
if (integrator == "md-vv")
("Expected tested environment variable to be " + envVariableModSimOn + " or " + envVariableModSimOff)
.c_str());
- const auto hasConservedField = !(tcoupling == "no" && pcoupling == "no");
+ const auto hasConservedField = !(tcoupling == "no" && pcoupling == "no")
+ && !(tcoupling == "andersen-massive" || tcoupling == "andersen");
SCOPED_TRACE(formatString(
"Comparing two simulations of '%s' "
pcoupling.c_str(),
environmentVariable.c_str()));
- const auto mdpFieldValues = prepareMdpFieldValues(
+ auto mdpFieldValues = prepareMdpFieldValues(
simulationName.c_str(), integrator.c_str(), tcoupling.c_str(), pcoupling.c_str());
+ if (tcoupling == "andersen")
+ {
+ // Fixes error "nstcomm must be 1, not 4 for Andersen, as velocities of
+ // atoms in coupled groups are randomized every time step"
+ mdpFieldValues["nstcomm"] = "1";
+ mdpFieldValues["nstcalcenergy"] = "1";
+ }
EnergyTermsToCompare energyTermsToCompare{ {
{ interaction_function[F_EPOT].longname, relativeToleranceAsPrecisionDependentUlp(60.0, 200, 160) },
SimulatorComparisonTest,
::testing::Combine(::testing::Combine(::testing::Values("argon12", "tip3p5"),
::testing::Values("md-vv"),
- ::testing::Values("no", "v-rescale", "berendsen"),
+ ::testing::Values("no",
+ "v-rescale",
+ "berendsen",
+ "andersen-massive",
+ "andersen"),
::testing::Values("no", "berendsen", "c-rescale")),
::testing::Values("GMX_DISABLE_MODULAR_SIMULATOR")));
INSTANTIATE_TEST_CASE_P(
SimulatorComparisonTest,
::testing::Combine(::testing::Combine(::testing::Values("argon12", "tip3p5"),
::testing::Values("md-vv"),
- ::testing::Values("no", "v-rescale", "berendsen"),
+ ::testing::Values("no",
+ "v-rescale",
+ "berendsen",
+ "andersen-massive",
+ "andersen"),
::testing::Values("no", "berendsen", "c-rescale")),
::testing::Values("GMX_DISABLE_MODULAR_SIMULATOR")));
INSTANTIATE_TEST_CASE_P(
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff