src/gromacs/modularsimulator/andersentemperaturecoupling.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief Defines Andersen temperature coupling for the modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "andersentemperaturecoupling.h"
  43
  44 #include <vector>
  45
  46 #include "gromacs/domdec/domdec_struct.h"
  47 #include "gromacs/math/functions.h"
  48 #include "gromacs/math/units.h"
  49 #include "gromacs/mdlib/constr.h"
  50 #include "gromacs/mdlib/mdatoms.h"
  51 #include "gromacs/mdlib/stat.h"
  52 #include "gromacs/mdtypes/commrec.h"
  53 #include "gromacs/mdtypes/inputrec.h"
  54 #include "gromacs/mdtypes/mdatom.h"
  55 #include "gromacs/mdrun/isimulator.h"
  56 #include "gromacs/random/tabulatednormaldistribution.h"
  57 #include "gromacs/random/threefry.h"
  58 #include "gromacs/random/uniformrealdistribution.h"
  59
  60 #include "compositesimulatorelement.h"
  61 #include "constraintelement.h"
  62 #include "simulatoralgorithm.h"
  63 #include "statepropagatordata.h"
  64
  65 namespace gmx
  66 {
  67 AndersenTemperatureCoupling::AndersenTemperatureCoupling(double               simulationTimestep,
  68                                                          bool                 doMassive,
  69                                                          int64_t              seed,
  70                                                          ArrayRef<const real> referenceTemperature,
  71                                                          ArrayRef<const real> couplingTime,
  72                                                          StatePropagatorData* statePropagatorData,
  73                                                          const MDAtoms*       mdAtoms,
  74                                                          const t_commrec*     cr) :
  75     doMassive_(doMassive),
  76     randomizationRate_(simulationTimestep / couplingTime[0]),
  77     couplingFrequency_(doMassive ? roundToInt(1. / randomizationRate_) : 1),
  78     seed_(seed),
  79     referenceTemperature_(referenceTemperature),
  80     couplingTime_(couplingTime),
  81     statePropagatorData_(statePropagatorData),
  82     mdAtoms_(mdAtoms->mdatoms()),
  83     cr_(cr)
  84 {
  85 }
  86
  87 void AndersenTemperatureCoupling::scheduleTask(Step                       step,
  88                                                Time gmx_unused            time,
  89                                                const RegisterRunFunction& registerRunFunction)
  90 {
  91     if (do_per_step(step, couplingFrequency_))
  92     {
  93         registerRunFunction([this, step]() { apply(step); });
  94     }
  95 }
  96
  97 void AndersenTemperatureCoupling::apply(Step step)
  98 {
  99     ThreeFry2x64<0>                       rng(seed_, RandomDomain::Thermostat);
 100     UniformRealDistribution<real>         uniformDist;
 101     TabulatedNormalDistribution<real, 14> normalDist;
 102
 103     const bool doDomainDecomposition = DOMAINDECOMP(cr_);
 104
 105     auto velocities = statePropagatorData_->velocitiesView().unpaddedArrayRef();
 106
 107     for (int atomIdx = 0; atomIdx < mdAtoms_->homenr; ++atomIdx)
 108     {
 109         const int temperatureGroup = mdAtoms_->cTC ? mdAtoms_->cTC[atomIdx] : 0;
 110         if (referenceTemperature_[temperatureGroup] <= 0 || couplingTime_[temperatureGroup] <= 0)
 111         {
 112             continue;
 113         }
 114
 115         const int globalAtomIdx = doDomainDecomposition ? cr_->dd->globalAtomIndices[atomIdx] : atomIdx;
 116         rng.restart(step, globalAtomIdx);
 117
 118         // For massive Andersen, this function is only called periodically, but we apply each time
 119         // Otherwise, this function is called every step, but we randomize atoms probabilistically
 120         if (!doMassive_)
 121         {
 122             uniformDist.reset();
 123         }
 124         if (doMassive_ || (uniformDist(rng) < randomizationRate_))
 125         {
 126             const real scalingFactor = std::sqrt(c_boltz * referenceTemperature_[temperatureGroup]
 127                                                  * mdAtoms_->invmass[atomIdx]);
 128             normalDist.reset();
 129             for (int d = 0; d < DIM; d++)
 130             {
 131                 velocities[atomIdx][d] = scalingFactor * normalDist(rng);
 132             }
 133         }
 134     }
 135 }
 136
 137 int AndersenTemperatureCoupling::frequency() const
 138 {
 139     return couplingFrequency_;
 140 }
 141
 142 void AndersenTemperatureCoupling::updateReferenceTemperature(ArrayRef<const real> gmx_unused temperatures,
 143                                                              ReferenceTemperatureChangeAlgorithm gmx_unused algorithm)
 144 {
 145     // Currently, we don't know about any temperature change algorithms, so we assert this never gets called
 146     GMX_ASSERT(false, "AndersenTemperatureCoupling: Unknown ReferenceTemperatureChangeAlgorithm.");
 147 }
 148
 149 void               AndersenTemperatureCoupling::elementSetup() {}
 150 ISimulatorElement* AndersenTemperatureCoupling::getElementPointerImpl(
 151         LegacySimulatorData*                    legacySimulatorData,
 152         ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 153         StatePropagatorData*                    statePropagatorData,
 154         EnergyData*                             energyData,
 155         FreeEnergyPerturbationData*             freeEnergyPerturbationData,
 156         GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
 157         ObservablesReducer gmx_unused* observablesReducer)
 158 {
 159     GMX_RELEASE_ASSERT(legacySimulatorData->inputrec->etc == TemperatureCoupling::Andersen
 160                                || legacySimulatorData->inputrec->etc == TemperatureCoupling::AndersenMassive,
 161                        "Expected the thermostat type to be andersen or andersen-massive.");
 162     auto andersenThermostat = std::make_unique<AndersenTemperatureCoupling>(
 163             legacySimulatorData->inputrec->delta_t,
 164             legacySimulatorData->inputrec->etc == TemperatureCoupling::AndersenMassive,
 165             legacySimulatorData->inputrec->andersen_seed,
 166             constArrayRefFromArray(legacySimulatorData->inputrec->opts.ref_t,
 167                                    legacySimulatorData->inputrec->opts.ngtc),
 168             constArrayRefFromArray(legacySimulatorData->inputrec->opts.tau_t,
 169                                    legacySimulatorData->inputrec->opts.ngtc),
 170             statePropagatorData,
 171             legacySimulatorData->mdAtoms,
 172             legacySimulatorData->cr);
 173     auto* andersenThermostatPtr = andersenThermostat.get();
 174     builderHelper->registerReferenceTemperatureUpdate(
 175             [andersenThermostatPtr](ArrayRef<const real>                temperatures,
 176                                     ReferenceTemperatureChangeAlgorithm algorithm) {
 177                 andersenThermostatPtr->updateReferenceTemperature(temperatures, algorithm);
 178             });
 179
 180     // T-coupling frequency will be composite element frequency
 181     const auto frequency = andersenThermostat->frequency();
 182     // Set up call list for composite element
 183     std::vector<compat::not_null<ISimulatorElement*>> elementCallList = { compat::make_not_null(
 184             andersenThermostat.get()) };
 185     // Set up element list for composite element
 186     std::vector<std::unique_ptr<gmx::ISimulatorElement>> elements;
 187     elements.emplace_back(std::move(andersenThermostat));
 188
 189     // If there are constraints, add constraint element after Andersen element
 190     if (legacySimulatorData->constr)
 191     {
 192         // This is excluded in preprocessing -
 193         // asserted here to make sure things don't get out of sync
 194         GMX_RELEASE_ASSERT(
 195                 legacySimulatorData->inputrec->etc == TemperatureCoupling::AndersenMassive,
 196                 "Per-particle Andersen thermostat is not implemented for systems with constrains.");
 197         // Build constraint element
 198         auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Velocities>>(
 199                 legacySimulatorData->constr,
 200                 statePropagatorData,
 201                 energyData,
 202                 freeEnergyPerturbationData,
 203                 MASTER(legacySimulatorData->cr),
 204                 legacySimulatorData->fplog,
 205                 legacySimulatorData->inputrec,
 206                 legacySimulatorData->mdAtoms->mdatoms());
 207         // Add call to composite element call list
 208         elementCallList.emplace_back(compat::make_not_null(constraintElement.get()));
 209         // Move ownership of constraint element to composite element
 210         elements.emplace_back(std::move(constraintElement));
 211     }
 212
 213     // Store composite element in builder helper and return pointer
 214     return builderHelper->storeElement(std::make_unique<CompositeSimulatorElement>(
 215             std::move(elementCallList), std::move(elements), frequency));
 216 }
 217
 218 } // namespace gmx