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/mdatoms.h"
  49 #include "gromacs/mdlib/stat.h"
  50 #include "gromacs/mdlib/update.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     isInitStep_(true),
  82     scalingTensor_(scalingTensor),
  83     propagatorCallback_(std::move(propagatorCallback)),
  84     statePropagatorData_(statePropagatorData),
  85     energyElement_(energyElement),
  86     fplog_(fplog),
  87     inputrec_(inputrec),
  88     mdAtoms_(mdAtoms)
  89 {
  90     clear_mat(mu_);
  91     clear_mat(boxRel_);
  92     clear_mat(boxVelocity_);
  93
  94     // TODO: This is only needed to restore the thermostatIntegral_ from cpt. Remove this when
  95     //       switching to purely client-based checkpointing.
  96     if (isRestart)
  97     {
  98         if (MASTER(cr))
  99         {
 100             copy_mat(globalState->boxv, boxVelocity_);
 101             copy_mat(globalState->box_rel, boxRel_);
 102         }
 103         if (DOMAINDECOMP(cr))
 104         {
 105             dd_bcast(cr->dd, sizeof(boxVelocity_), boxVelocity_);
 106             dd_bcast(cr->dd, sizeof(boxRel_), boxRel_);
 107         }
 108     }
 109 }
 110
 111 void ParrinelloRahmanBarostat::scheduleTask(gmx::Step step,
 112                                             gmx::Time gmx_unused               time,
 113                                             const gmx::RegisterRunFunctionPtr& registerRunFunction)
 114 {
 115     const bool scaleOnNextStep = do_per_step(step + nstpcouple_ + offset_ + 1, nstpcouple_);
 116     const bool scaleOnThisStep = do_per_step(step + nstpcouple_ + offset_, nstpcouple_);
 117
 118     if (scaleOnThisStep)
 119     {
 120         (*registerRunFunction)(
 121                 std::make_unique<SimulatorRunFunction>([this]() { scaleBoxAndPositions(); }));
 122     }
 123     if (scaleOnNextStep)
 124     {
 125         (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
 126                 [this, step]() { integrateBoxVelocityEquations(step); }));
 127         // let propagator know that it will have to scale on next step
 128         (*propagatorCallback_)(step + 1);
 129     }
 130 }
 131
 132 void ParrinelloRahmanBarostat::integrateBoxVelocityEquations(Step step)
 133 {
 134     auto box = statePropagatorData_->constBox();
 135     parrinellorahman_pcoupl(fplog_, step, inputrec_, couplingTimeStep_, energyElement_->pressure(step),
 136                             box, boxRel_, boxVelocity_, scalingTensor_.data(), mu_, isInitStep_);
 137     // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
 138     msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
 139 }
 140
 141 void ParrinelloRahmanBarostat::scaleBoxAndPositions()
 142 {
 143     // Propagate the box by the box velocities
 144     auto box = statePropagatorData_->box();
 145     for (int i = 0; i < DIM; i++)
 146     {
 147         for (int m = 0; m <= i; m++)
 148         {
 149             box[i][m] += couplingTimeStep_ * boxVelocity_[i][m];
 150         }
 151     }
 152     preserve_box_shape(inputrec_, boxRel_, box);
 153
 154     // Scale the coordinates
 155     const int start  = 0;
 156     const int homenr = mdAtoms_->mdatoms()->homenr;
 157     auto      x      = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 158     for (int n = start; n < start + homenr; n++)
 159     {
 160         tmvmul_ur0(mu_, x[n], x[n]);
 161     }
 162 }
 163
 164 void ParrinelloRahmanBarostat::elementSetup()
 165 {
 166     if (inputrecPreserveShape(inputrec_))
 167     {
 168         auto      box  = statePropagatorData_->box();
 169         const int ndim = inputrec_->epct == epctSEMIISOTROPIC ? 2 : 3;
 170         do_box_rel(ndim, inputrec_->deform, boxRel_, box, true);
 171     }
 172
 173     const bool scaleOnInitStep = do_per_step(initStep_ + nstpcouple_ + offset_, nstpcouple_);
 174     if (scaleOnInitStep)
 175     {
 176         integrateBoxVelocityEquations(initStep_);
 177         (*propagatorCallback_)(initStep_ + 1);
 178     }
 179     isInitStep_ = false;
 180 }
 181
 182 const rvec* ParrinelloRahmanBarostat::boxVelocities() const
 183 {
 184     return boxVelocity_;
 185 }
 186
 187 void ParrinelloRahmanBarostat::writeCheckpoint(t_state* localState, t_state gmx_unused* globalState)
 188 {
 189     copy_mat(boxVelocity_, localState->boxv);
 190     copy_mat(boxRel_, localState->box_rel);
 191     localState->flags |= (1U << estBOXV) | (1U << estBOX_REL);
 192 }
 193
 194
 195 } // namespace gmx