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     rvec *x, *xprime, *min_proj, *v;
 130
 131     const real lambdaBonded = freeEnergyPerturbationElement_
 132                                       ? freeEnergyPerturbationElement_->constLambdaView()[efptBONDED]
 133                                       : 0;
 134     real dvdlambda = 0;
 135
 136     switch (variable)
 137     {
 138         case ConstraintVariable::Positions:
 139             x = as_rvec_array(statePropagatorData_->previousPositionsView().paddedArrayRef().data());
 140             xprime   = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 141             min_proj = nullptr;
 142             v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
 143             break;
 144         case ConstraintVariable::Velocities:
 145             x      = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 146             xprime = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
 147             min_proj = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
 148             v        = nullptr;
 149             break;
 150         default: gmx_fatal(FARGS, "Constraint algorithm not implemented for modular simulator.");
 151     }
 152
 153     constr_->apply(writeLog, writeEnergy, step, 1, 1.0, x, xprime, min_proj, statePropagatorData_->box(),
 154                    lambdaBonded, &dvdlambda, v, calculateVirial ? &vir_con : nullptr, variable);
 155
 156     if (calculateVirial)
 157     {
 158         if (inputrec_->eI == eiVV)
 159         {
 160             // For some reason, the shake virial in VV is reset twice a step.
 161             // Energy element will only do this once per step.
 162             // TODO: Investigate this
 163             clear_mat(energyElement_->constraintVirial(step));
 164         }
 165         energyElement_->addToConstraintVirial(vir_con, step);
 166     }
 167
 168     /* The factor of 2 correction is necessary because half of the constraint
 169      * force is removed in the VV step. This factor is either exact or a very
 170      * good approximation, statistically insignificant in any real free energy
 171      * calculation. Any possible error is not a simulation propagation error,
 172      * but a potential reporting error in the data that goes to dh/dlambda.
 173      * Cf. redmine issue #1255
 174      */
 175     const real c_dvdlConstraintCorrectionFactor = EI_VV(inputrec_->eI) ? 2.0 : 1.0;
 176     energyElement_->enerdata()->term[F_DVDL_CONSTR] += c_dvdlConstraintCorrectionFactor * dvdlambda;
 177 }
 178
 179 template<ConstraintVariable variable>
 180 SignallerCallbackPtr ConstraintsElement<variable>::registerEnergyCallback(EnergySignallerEvent event)
 181 {
 182     if (event == EnergySignallerEvent::VirialCalculationStep)
 183     {
 184         return std::make_unique<SignallerCallback>(
 185                 [this](Step step, Time /*unused*/) { nextVirialCalculationStep_ = step; });
 186     }
 187     return nullptr;
 188 }
 189
 190 template<ConstraintVariable variable>
 191 SignallerCallbackPtr ConstraintsElement<variable>::registerTrajectorySignallerCallback(TrajectoryEvent event)
 192 {
 193     if (event == TrajectoryEvent::EnergyWritingStep)
 194     {
 195         return std::make_unique<SignallerCallback>(
 196                 [this](Step step, Time /*unused*/) { nextEnergyWritingStep_ = step; });
 197     }
 198     return nullptr;
 199 }
 200
 201 template<ConstraintVariable variable>
 202 SignallerCallbackPtr ConstraintsElement<variable>::registerLoggingCallback()
 203 {
 204     return std::make_unique<SignallerCallback>(
 205             [this](Step step, Time /*unused*/) { nextLogWritingStep_ = step; });
 206 }
 207
 208 //! Explicit template initialization
 209 //! @{
 210 template class ConstraintsElement<ConstraintVariable::Positions>;
 211 template class ConstraintsElement<ConstraintVariable::Velocities>;
 212 //! @}
 213
 214 } // namespace gmx