src/gromacs/modularsimulator/parrinellorahmanbarostat.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 Parrinello-Rahman barostat 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 "parrinellorahmanbarostat.h"
  45
  46 #include "gromacs/domdec/domdec_network.h"
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/mdlib/coupling.h"
  49 #include "gromacs/mdlib/mdatoms.h"
  50 #include "gromacs/mdlib/stat.h"
  51 #include "gromacs/mdtypes/commrec.h"
  52 #include "gromacs/mdtypes/inputrec.h"
  53 #include "gromacs/mdtypes/mdatom.h"
  54 #include "gromacs/mdtypes/state.h"
  55 #include "gromacs/pbcutil/boxutilities.h"
  56
  57 #include "energyelement.h"
  58 #include "statepropagatordata.h"
  59
  60 namespace gmx
  61 {
  62
  63 ParrinelloRahmanBarostat::ParrinelloRahmanBarostat(int                   nstpcouple,
  64                                                    int                   offset,
  65                                                    real                  couplingTimeStep,
  66                                                    Step                  initStep,
  67                                                    ArrayRef<rvec>        scalingTensor,
  68                                                    PropagatorCallbackPtr propagatorCallback,
  69                                                    StatePropagatorData*  statePropagatorData,
  70                                                    EnergyElement*        energyElement,
  71                                                    FILE*                 fplog,
  72                                                    const t_inputrec*     inputrec,
  73                                                    const MDAtoms*        mdAtoms,
  74                                                    const t_state*        globalState,
  75                                                    t_commrec*            cr,
  76                                                    bool                  isRestart) :
  77     nstpcouple_(nstpcouple),
  78     offset_(offset),
  79     couplingTimeStep_(couplingTimeStep),
  80     initStep_(initStep),
  81     scalingTensor_(scalingTensor),
  82     propagatorCallback_(std::move(propagatorCallback)),
  83     statePropagatorData_(statePropagatorData),
  84     energyElement_(energyElement),
  85     fplog_(fplog),
  86     inputrec_(inputrec),
  87     mdAtoms_(mdAtoms)
  88 {
  89     clear_mat(mu_);
  90     clear_mat(boxRel_);
  91     clear_mat(boxVelocity_);
  92
  93     // TODO: This is only needed to restore the thermostatIntegral_ from cpt. Remove this when
  94     //       switching to purely client-based checkpointing.
  95     if (isRestart)
  96     {
  97         if (MASTER(cr))
  98         {
  99             copy_mat(globalState->boxv, boxVelocity_);
 100             copy_mat(globalState->box_rel, boxRel_);
 101         }
 102         if (DOMAINDECOMP(cr))
 103         {
 104             dd_bcast(cr->dd, sizeof(boxVelocity_), boxVelocity_);
 105             dd_bcast(cr->dd, sizeof(boxRel_), boxRel_);
 106         }
 107     }
 108 }
 109
 110 void ParrinelloRahmanBarostat::scheduleTask(gmx::Step step,
 111                                             gmx::Time gmx_unused               time,
 112                                             const gmx::RegisterRunFunctionPtr& registerRunFunction)
 113 {
 114     const bool scaleOnNextStep = do_per_step(step + nstpcouple_ + offset_ + 1, nstpcouple_);
 115     const bool scaleOnThisStep = do_per_step(step + nstpcouple_ + offset_, nstpcouple_);
 116
 117     if (scaleOnThisStep)
 118     {
 119         (*registerRunFunction)(
 120                 std::make_unique<SimulatorRunFunction>([this]() { scaleBoxAndPositions(); }));
 121     }
 122     if (scaleOnNextStep)
 123     {
 124         (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
 125                 [this, step]() { integrateBoxVelocityEquations(step); }));
 126         // let propagator know that it will have to scale on next step
 127         (*propagatorCallback_)(step + 1);
 128     }
 129 }
 130
 131 void ParrinelloRahmanBarostat::integrateBoxVelocityEquations(Step step)
 132 {
 133     auto box = statePropagatorData_->constBox();
 134     parrinellorahman_pcoupl(fplog_, step, inputrec_, couplingTimeStep_, energyElement_->pressure(step),
 135                             box, boxRel_, boxVelocity_, scalingTensor_.data(), mu_, false);
 136     // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
 137     msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
 138 }
 139
 140 void ParrinelloRahmanBarostat::scaleBoxAndPositions()
 141 {
 142     // Propagate the box by the box velocities
 143     auto box = statePropagatorData_->box();
 144     for (int i = 0; i < DIM; i++)
 145     {
 146         for (int m = 0; m <= i; m++)
 147         {
 148             box[i][m] += couplingTimeStep_ * boxVelocity_[i][m];
 149         }
 150     }
 151     preserve_box_shape(inputrec_, boxRel_, box);
 152
 153     // Scale the coordinates
 154     const int start  = 0;
 155     const int homenr = mdAtoms_->mdatoms()->homenr;
 156     auto      x      = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 157     for (int n = start; n < start + homenr; n++)
 158     {
 159         tmvmul_ur0(mu_, x[n], x[n]);
 160     }
 161 }
 162
 163 void ParrinelloRahmanBarostat::elementSetup()
 164 {
 165     if (inputrecPreserveShape(inputrec_))
 166     {
 167         auto      box  = statePropagatorData_->box();
 168         const int ndim = inputrec_->epct == epctSEMIISOTROPIC ? 2 : 3;
 169         do_box_rel(ndim, inputrec_->deform, boxRel_, box, true);
 170     }
 171
 172     const bool scaleOnInitStep = do_per_step(initStep_ + nstpcouple_ + offset_, nstpcouple_);
 173     if (scaleOnInitStep)
 174     {
 175         // If we need to scale on the first step, we need to set the scaling matrix using the current
 176         // box velocity. If this is a fresh start, we will hence not move the box (this does currently
 177         // never happen as the offset is set to -1 in all cases). If this is a restart, we will use
 178         // the saved box velocity which we would have updated right before checkpointing.
 179         // Setting bFirstStep = true in parrinellorahman_pcoupl (last argument) makes sure that only
 180         // the scaling matrix is calculated, without updating the box velocities.
 181         // The call to parrinellorahman_pcoupl is using nullptr for fplog (since we don't expect any
 182         // output here) and for the pressure (since it might not be calculated yet, and we don't need it).
 183         auto box = statePropagatorData_->constBox();
 184         parrinellorahman_pcoupl(nullptr, initStep_, inputrec_, couplingTimeStep_, nullptr, box,
 185                                 boxRel_, boxVelocity_, scalingTensor_.data(), mu_, true);
 186         // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
 187         msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
 188
 189         (*propagatorCallback_)(initStep_);
 190     }
 191 }
 192
 193 const rvec* ParrinelloRahmanBarostat::boxVelocities() const
 194 {
 195     return boxVelocity_;
 196 }
 197
 198 void ParrinelloRahmanBarostat::writeCheckpoint(t_state* localState, t_state gmx_unused* globalState)
 199 {
 200     copy_mat(boxVelocity_, localState->boxv);
 201     copy_mat(boxRel_, localState->box_rel);
 202     localState->flags |= (1U << estBOXV) | (1U << estBOX_REL);
 203 }
 204
 205
 206 } // namespace gmx