src/gromacs/modularsimulator/vrescalethermostat.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, 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 the v-rescale thermostat for the modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "vrescalethermostat.h"
  45
  46 #include "gromacs/domdec/domdec_network.h"
  47 #include "gromacs/math/units.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/mdlib/stat.h"
  50 #include "gromacs/mdlib/update.h"
  51 #include "gromacs/mdtypes/commrec.h"
  52 #include "gromacs/mdtypes/group.h"
  53 #include "gromacs/mdtypes/state.h"
  54 #include "gromacs/utility/fatalerror.h"
  55
  56 namespace gmx
  57 {
  58
  59 VRescaleThermostat::VRescaleThermostat(int                   nstcouple,
  60                                        int                   offset,
  61                                        bool                  useFullStepKE,
  62                                        int64_t               seed,
  63                                        int                   numTemperatureGroups,
  64                                        double                couplingTimeStep,
  65                                        const real*           referenceTemperature,
  66                                        const real*           couplingTime,
  67                                        const real*           numDegreesOfFreedom,
  68                                        EnergyElement*        energyElement,
  69                                        ArrayRef<real>        lambdaView,
  70                                        PropagatorCallbackPtr propagatorCallback,
  71                                        const t_state*        globalState,
  72                                        t_commrec*            cr,
  73                                        bool                  isRestart) :
  74     nstcouple_(nstcouple),
  75     offset_(offset),
  76     useFullStepKE_(useFullStepKE),
  77     seed_(seed),
  78     numTemperatureGroups_(numTemperatureGroups),
  79     couplingTimeStep_(couplingTimeStep),
  80     referenceTemperature_(referenceTemperature, referenceTemperature + numTemperatureGroups),
  81     couplingTime_(couplingTime, couplingTime + numTemperatureGroups),
  82     numDegreesOfFreedom_(numDegreesOfFreedom, numDegreesOfFreedom + numTemperatureGroups),
  83     thermostatIntegral_(numTemperatureGroups, 0.0),
  84     energyElement_(energyElement),
  85     lambda_(lambdaView),
  86     propagatorCallback_(std::move(propagatorCallback))
  87 {
  88     // TODO: This is only needed to restore the thermostatIntegral_ from cpt. Remove this when
  89     //       switching to purely client-based checkpointing.
  90     if (isRestart)
  91     {
  92         if (MASTER(cr))
  93         {
  94             for (unsigned long i = 0; i < thermostatIntegral_.size(); ++i)
  95             {
  96                 thermostatIntegral_[i] = globalState->therm_integral[i];
  97             }
  98         }
  99         if (DOMAINDECOMP(cr))
 100         {
 101             dd_bcast(cr->dd, int(thermostatIntegral_.size() * sizeof(double)), thermostatIntegral_.data());
 102         }
 103     }
 104 }
 105
 106 void VRescaleThermostat::scheduleTask(Step step, Time gmx_unused time, const RegisterRunFunctionPtr& registerRunFunction)
 107 {
 108     /* The thermostat will need a valid kinetic energy when it is running.
 109      * Currently, computeGlobalCommunicationPeriod() is making sure this
 110      * happens on time.
 111      * TODO: Once we're switching to a new global communication scheme, we
 112      *       will want the thermostat to signal that global reduction
 113      *       of the kinetic energy is needed.
 114      *
 115      */
 116     if (do_per_step(step + nstcouple_ + offset_, nstcouple_))
 117     {
 118         // do T-coupling this step
 119         (*registerRunFunction)(
 120                 std::make_unique<SimulatorRunFunction>([this, step]() { setLambda(step); }));
 121
 122         // Let propagator know that we want to do T-coupling
 123         (*propagatorCallback_)(step);
 124     }
 125 }
 126
 127 void VRescaleThermostat::setLambda(Step step)
 128 {
 129     real currentKineticEnergy, referenceKineticEnergy, newKineticEnergy;
 130
 131     auto ekind = energyElement_->ekindata();
 132
 133     for (int i = 0; (i < numTemperatureGroups_); i++)
 134     {
 135         if (useFullStepKE_)
 136         {
 137             currentKineticEnergy = trace(ekind->tcstat[i].ekinf);
 138         }
 139         else
 140         {
 141             currentKineticEnergy = trace(ekind->tcstat[i].ekinh);
 142         }
 143
 144         if (couplingTime_[i] >= 0 && numDegreesOfFreedom_[i] > 0 && currentKineticEnergy > 0)
 145         {
 146             referenceKineticEnergy = 0.5 * referenceTemperature_[i] * BOLTZ * numDegreesOfFreedom_[i];
 147
 148             newKineticEnergy = vrescale_resamplekin(currentKineticEnergy, referenceKineticEnergy,
 149                                                     numDegreesOfFreedom_[i],
 150                                                     couplingTime_[i] / couplingTimeStep_, step, seed_);
 151
 152             // Analytically newKineticEnergy >= 0, but we check for rounding errors
 153             if (newKineticEnergy <= 0)
 154             {
 155                 lambda_[i] = 0.0;
 156             }
 157             else
 158             {
 159                 lambda_[i] = std::sqrt(newKineticEnergy / currentKineticEnergy);
 160             }
 161
 162             thermostatIntegral_[i] -= newKineticEnergy - currentKineticEnergy;
 163
 164             if (debug)
 165             {
 166                 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n", i,
 167                         referenceKineticEnergy, currentKineticEnergy, newKineticEnergy, lambda_[i]);
 168             }
 169         }
 170         else
 171         {
 172             lambda_[i] = 1.0;
 173         }
 174     }
 175 }
 176
 177 void VRescaleThermostat::writeCheckpoint(t_state* localState, t_state gmx_unused* globalState)
 178 {
 179     localState->therm_integral = thermostatIntegral_;
 180     localState->flags |= (1U << estTHERM_INT);
 181 }
 182
 183 const std::vector<double>& VRescaleThermostat::thermostatIntegral() const
 184 {
 185     return thermostatIntegral_;
 186 }
 187
 188 } // namespace gmx