src/gromacs/modularsimulator/constraintelement.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,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 the constraint element 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 "constraintelement.h"
  45
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/mdlib/mdatoms.h"
  48 #include "gromacs/mdtypes/commrec.h"
  49 #include "gromacs/mdtypes/enerdata.h"
  50 #include "gromacs/mdtypes/inputrec.h"
  51 #include "gromacs/mdtypes/mdatom.h"
  52 #include "gromacs/mdtypes/state.h"
  53 #include "gromacs/utility/fatalerror.h"
  54
  55 #include "energydata.h"
  56 #include "freeenergyperturbationdata.h"
  57 #include "modularsimulator.h"
  58 #include "simulatoralgorithm.h"
  59 #include "statepropagatordata.h"
  60
  61 namespace gmx
  62 {
  63 template<ConstraintVariable variable>
  64 ConstraintsElement<variable>::ConstraintsElement(Constraints*                constr,
  65                                                  StatePropagatorData*        statePropagatorData,
  66                                                  EnergyData*                 energyData,
  67                                                  FreeEnergyPerturbationData* freeEnergyPerturbationData,
  68                                                  bool                        isMaster,
  69                                                  FILE*                       fplog,
  70                                                  const t_inputrec*           inputrec,
  71                                                  const t_mdatoms*            mdAtoms) :
  72     nextVirialCalculationStep_(-1),
  73     nextEnergyWritingStep_(-1),
  74     nextLogWritingStep_(-1),
  75     isMasterRank_(isMaster),
  76     statePropagatorData_(statePropagatorData),
  77     energyData_(energyData),
  78     freeEnergyPerturbationData_(freeEnergyPerturbationData),
  79     constr_(constr),
  80     fplog_(fplog),
  81     inputrec_(inputrec),
  82     mdAtoms_(mdAtoms)
  83 {
  84     GMX_ASSERT(constr_, "Constraint element created but constr == nullptr");
  85 }
  86
  87 template<ConstraintVariable variable>
  88 void ConstraintsElement<variable>::elementSetup()
  89 {
  90     if (!inputrec_->bContinuation
  91         && ((variable == ConstraintVariable::Positions && inputrec_->eI == IntegrationAlgorithm::MD)
  92             || (variable == ConstraintVariable::Velocities && inputrec_->eI == IntegrationAlgorithm::VV)))
  93     {
  94         const real lambdaBonded =
  95                 freeEnergyPerturbationData_
  96                         ? freeEnergyPerturbationData_->constLambdaView()[static_cast<int>(
  97                                   FreeEnergyPerturbationCouplingType::Bonded)]
  98                         : 0;
  99         // Constrain the initial coordinates and velocities
 100         do_constrain_first(fplog_,
 101                            constr_,
 102                            inputrec_,
 103                            statePropagatorData_->totalNumAtoms(),
 104                            statePropagatorData_->localNumAtoms(),
 105                            statePropagatorData_->positionsView(),
 106                            statePropagatorData_->velocitiesView(),
 107                            statePropagatorData_->box(),
 108                            lambdaBonded);
 109
 110         if (isMasterRank_)
 111         {
 112             if (inputrec_->eConstrAlg == ConstraintAlgorithm::Lincs)
 113             {
 114                 fprintf(fplog_,
 115                         "RMS relative constraint deviation after constraining: %.2e\n",
 116                         constr_->rmsd());
 117             }
 118         }
 119     }
 120 }
 121
 122 template<ConstraintVariable variable>
 123 void ConstraintsElement<variable>::scheduleTask(Step step,
 124                                                 Time gmx_unused            time,
 125                                                 const RegisterRunFunction& registerRunFunction)
 126 {
 127     bool calculateVirial = (step == nextVirialCalculationStep_);
 128     bool writeLog        = (step == nextLogWritingStep_);
 129     bool writeEnergy     = (step == nextEnergyWritingStep_);
 130
 131     // register constraining
 132     registerRunFunction([this, step, calculateVirial, writeLog, writeEnergy]() {
 133         apply(step, calculateVirial, writeLog, writeEnergy);
 134     });
 135 }
 136
 137 template<ConstraintVariable variable>
 138 void ConstraintsElement<variable>::apply(Step step, bool calculateVirial, bool writeLog, bool writeEnergy)
 139 {
 140     tensor vir_con;
 141
 142     ArrayRefWithPadding<RVec> x;
 143     ArrayRefWithPadding<RVec> xprime;
 144     ArrayRef<RVec>            min_proj;
 145     ArrayRefWithPadding<RVec> v;
 146
 147     const real lambdaBonded =
 148             freeEnergyPerturbationData_
 149                     ? freeEnergyPerturbationData_
 150                               ->constLambdaView()[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)]
 151                     : 0;
 152     real dvdlambda = 0;
 153
 154     switch (variable)
 155     {
 156         case ConstraintVariable::Positions:
 157             x      = statePropagatorData_->previousPositionsView();
 158             xprime = statePropagatorData_->positionsView();
 159             v      = statePropagatorData_->velocitiesView();
 160             break;
 161         case ConstraintVariable::Velocities:
 162             x        = statePropagatorData_->positionsView();
 163             xprime   = statePropagatorData_->velocitiesView();
 164             min_proj = statePropagatorData_->velocitiesView().unpaddedArrayRef();
 165             break;
 166         default: gmx_fatal(FARGS, "Constraint algorithm not implemented for modular simulator.");
 167     }
 168
 169     constr_->apply(writeLog,
 170                    writeEnergy,
 171                    step,
 172                    1,
 173                    1.0,
 174                    x,
 175                    xprime,
 176                    min_proj,
 177                    statePropagatorData_->box(),
 178                    lambdaBonded,
 179                    &dvdlambda,
 180                    v,
 181                    calculateVirial,
 182                    vir_con,
 183                    variable);
 184
 185     if (calculateVirial)
 186     {
 187         if (inputrec_->eI == IntegrationAlgorithm::VV)
 188         {
 189             // For some reason, the shake virial in VV is reset twice a step.
 190             // Energy element will only do this once per step.
 191             // TODO: Investigate this
 192             clear_mat(energyData_->constraintVirial(step));
 193         }
 194         energyData_->addToConstraintVirial(vir_con, step);
 195     }
 196
 197     /* The factor of 2 correction is necessary because half of the constraint
 198      * force is removed in the VV step. This factor is either exact or a very
 199      * good approximation, statistically insignificant in any real free energy
 200      * calculation. Any possible error is not a simulation propagation error,
 201      * but a potential reporting error in the data that goes to dh/dlambda.
 202      * Cf. Issue #1255
 203      */
 204     const real c_dvdlConstraintCorrectionFactor = EI_VV(inputrec_->eI) ? 2.0 : 1.0;
 205     energyData_->enerdata()->term[F_DVDL_CONSTR] += c_dvdlConstraintCorrectionFactor * dvdlambda;
 206 }
 207
 208 template<ConstraintVariable variable>
 209 std::optional<SignallerCallback> ConstraintsElement<variable>::registerEnergyCallback(EnergySignallerEvent event)
 210 {
 211     if (event == EnergySignallerEvent::VirialCalculationStep)
 212     {
 213         return [this](Step step, Time /*unused*/) { nextVirialCalculationStep_ = step; };
 214     }
 215     return std::nullopt;
 216 }
 217
 218 template<ConstraintVariable variable>
 219 std::optional<SignallerCallback> ConstraintsElement<variable>::registerTrajectorySignallerCallback(TrajectoryEvent event)
 220 {
 221     if (event == TrajectoryEvent::EnergyWritingStep)
 222     {
 223         return [this](Step step, Time /*unused*/) { nextEnergyWritingStep_ = step; };
 224     }
 225     return std::nullopt;
 226 }
 227
 228 template<ConstraintVariable variable>
 229 std::optional<SignallerCallback> ConstraintsElement<variable>::registerLoggingCallback()
 230 {
 231     return [this](Step step, Time /*unused*/) { nextLogWritingStep_ = step; };
 232 }
 233
 234 template<ConstraintVariable variable>
 235 ISimulatorElement* ConstraintsElement<variable>::getElementPointerImpl(
 236         LegacySimulatorData*                    legacySimulatorData,
 237         ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 238         StatePropagatorData*                    statePropagatorData,
 239         EnergyData*                             energyData,
 240         FreeEnergyPerturbationData*             freeEnergyPerturbationData,
 241         GlobalCommunicationHelper gmx_unused* globalCommunicationHelper)
 242 {
 243     return builderHelper->storeElement(
 244             std::make_unique<ConstraintsElement<variable>>(legacySimulatorData->constr,
 245                                                            statePropagatorData,
 246                                                            energyData,
 247                                                            freeEnergyPerturbationData,
 248                                                            MASTER(legacySimulatorData->cr),
 249                                                            legacySimulatorData->fplog,
 250                                                            legacySimulatorData->inputrec,
 251                                                            legacySimulatorData->mdAtoms->mdatoms()));
 252 }
 253
 254 // Explicit template initializations
 255 template class ConstraintsElement<ConstraintVariable::Positions>;
 256 template class ConstraintsElement<ConstraintVariable::Velocities>;
 257
 258 } // namespace gmx