src/gromacs/modularsimulator/constraintelement.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020, 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/mdtypes/enerdata.h"
  48 #include "gromacs/mdtypes/inputrec.h"
  49 #include "gromacs/mdtypes/state.h"
  50 #include "gromacs/utility/fatalerror.h"
  51
  52 #include "energyelement.h"
  53 #include "freeenergyperturbationelement.h"
  54 #include "statepropagatordata.h"
  55
  56 namespace gmx
  57 {
  58 template<ConstraintVariable variable>
  59 ConstraintsElement<variable>::ConstraintsElement(Constraints*                   constr,
  60                                                  StatePropagatorData*           statePropagatorData,
  61                                                  EnergyElement*                 energyElement,
  62                                                  FreeEnergyPerturbationElement* freeEnergyPerturbationElement,
  63                                                  bool                           isMaster,
  64                                                  FILE*                          fplog,
  65                                                  const t_inputrec*              inputrec,
  66                                                  const t_mdatoms*               mdAtoms) :
  67     nextVirialCalculationStep_(-1),
  68     nextEnergyWritingStep_(-1),
  69     nextLogWritingStep_(-1),
  70     isMasterRank_(isMaster),
  71     statePropagatorData_(statePropagatorData),
  72     energyElement_(energyElement),
  73     freeEnergyPerturbationElement_(freeEnergyPerturbationElement),
  74     constr_(constr),
  75     fplog_(fplog),
  76     inputrec_(inputrec),
  77     mdAtoms_(mdAtoms)
  78 {
  79     GMX_ASSERT(constr_, "Constraint element created but constr == nullptr");
  80 }
  81
  82 template<ConstraintVariable variable>
  83 void ConstraintsElement<variable>::elementSetup()
  84 {
  85     if (!inputrec_->bContinuation
  86         && ((variable == ConstraintVariable::Positions && inputrec_->eI == eiMD)
  87             || (variable == ConstraintVariable::Velocities && inputrec_->eI == eiVV)))
  88     {
  89         const real lambdaBonded = freeEnergyPerturbationElement_
  90                                           ? freeEnergyPerturbationElement_->constLambdaView()[efptBONDED]
  91                                           : 0;
  92         // Constrain the initial coordinates and velocities
  93         do_constrain_first(fplog_, constr_, inputrec_, mdAtoms_, statePropagatorData_->localNumAtoms(),
  94                            statePropagatorData_->positionsView(), statePropagatorData_->velocitiesView(),
  95                            statePropagatorData_->box(), lambdaBonded);
  96
  97         if (isMasterRank_)
  98         {
  99             if (inputrec_->eConstrAlg == econtLINCS)
 100             {
 101                 fprintf(fplog_, "RMS relative constraint deviation after constraining: %.2e\n",
 102                         constr_->rmsd());
 103             }
 104         }
 105     }
 106 }
 107
 108 template<ConstraintVariable variable>
 109 void ConstraintsElement<variable>::scheduleTask(Step step,
 110                                                 Time gmx_unused               time,
 111                                                 const RegisterRunFunctionPtr& registerRunFunction)
 112 {
 113     bool calculateVirial = (step == nextVirialCalculationStep_);
 114     bool writeLog        = (step == nextLogWritingStep_);
 115     bool writeEnergy     = (step == nextEnergyWritingStep_);
 116
 117     // register constraining
 118     (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
 119             [this, step, calculateVirial, writeLog, writeEnergy]() {
 120                 apply(step, calculateVirial, writeLog, writeEnergy);
 121             }));
 122 }
 123
 124 template<ConstraintVariable variable>
 125 void ConstraintsElement<variable>::apply(Step step, bool calculateVirial, bool writeLog, bool writeEnergy)
 126 {
 127     tensor vir_con;
 128
 129     ArrayRefWithPadding<RVec> x;
 130     ArrayRefWithPadding<RVec> xprime;
 131     ArrayRef<RVec>            min_proj;
 132     ArrayRefWithPadding<RVec> v;
 133
 134     const real lambdaBonded = freeEnergyPerturbationElement_
 135                                       ? freeEnergyPerturbationElement_->constLambdaView()[efptBONDED]
 136                                       : 0;
 137     real dvdlambda = 0;
 138
 139     switch (variable)
 140     {
 141         case ConstraintVariable::Positions:
 142             x      = statePropagatorData_->previousPositionsView();
 143             xprime = statePropagatorData_->positionsView();
 144             v      = statePropagatorData_->velocitiesView();
 145             break;
 146         case ConstraintVariable::Velocities:
 147             x        = statePropagatorData_->positionsView();
 148             xprime   = statePropagatorData_->velocitiesView();
 149             min_proj = statePropagatorData_->velocitiesView().unpaddedArrayRef();
 150             break;
 151         default: gmx_fatal(FARGS, "Constraint algorithm not implemented for modular simulator.");
 152     }
 153
 154     constr_->apply(writeLog, writeEnergy, step, 1, 1.0, x, xprime, min_proj, statePropagatorData_->box(),
 155                    lambdaBonded, &dvdlambda, v, calculateVirial ? &vir_con : nullptr, variable);
 156
 157     if (calculateVirial)
 158     {
 159         if (inputrec_->eI == eiVV)
 160         {
 161             // For some reason, the shake virial in VV is reset twice a step.
 162             // Energy element will only do this once per step.
 163             // TODO: Investigate this
 164             clear_mat(energyElement_->constraintVirial(step));
 165         }
 166         energyElement_->addToConstraintVirial(vir_con, step);
 167     }
 168
 169     /* The factor of 2 correction is necessary because half of the constraint
 170      * force is removed in the VV step. This factor is either exact or a very
 171      * good approximation, statistically insignificant in any real free energy
 172      * calculation. Any possible error is not a simulation propagation error,
 173      * but a potential reporting error in the data that goes to dh/dlambda.
 174      * Cf. Issue #1255
 175      */
 176     const real c_dvdlConstraintCorrectionFactor = EI_VV(inputrec_->eI) ? 2.0 : 1.0;
 177     energyElement_->enerdata()->term[F_DVDL_CONSTR] += c_dvdlConstraintCorrectionFactor * dvdlambda;
 178 }
 179
 180 template<ConstraintVariable variable>
 181 SignallerCallbackPtr ConstraintsElement<variable>::registerEnergyCallback(EnergySignallerEvent event)
 182 {
 183     if (event == EnergySignallerEvent::VirialCalculationStep)
 184     {
 185         return std::make_unique<SignallerCallback>(
 186                 [this](Step step, Time /*unused*/) { nextVirialCalculationStep_ = step; });
 187     }
 188     return nullptr;
 189 }
 190
 191 template<ConstraintVariable variable>
 192 SignallerCallbackPtr ConstraintsElement<variable>::registerTrajectorySignallerCallback(TrajectoryEvent event)
 193 {
 194     if (event == TrajectoryEvent::EnergyWritingStep)
 195     {
 196         return std::make_unique<SignallerCallback>(
 197                 [this](Step step, Time /*unused*/) { nextEnergyWritingStep_ = step; });
 198     }
 199     return nullptr;
 200 }
 201
 202 template<ConstraintVariable variable>
 203 SignallerCallbackPtr ConstraintsElement<variable>::registerLoggingCallback()
 204 {
 205     return std::make_unique<SignallerCallback>(
 206             [this](Step step, Time /*unused*/) { nextLogWritingStep_ = step; });
 207 }
 208
 209 //! Explicit template initialization
 210 //! @{
 211 template class ConstraintsElement<ConstraintVariable::Positions>;
 212 template class ConstraintsElement<ConstraintVariable::Velocities>;
 213 //! @}
 214
 215 } // namespace gmx