src/gromacs/mdlib/update_vv.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "update_vv.h"
  41
  42 #include <cmath>
  43 #include <cstdio>
  44
  45 #include <algorithm>
  46 #include <memory>
  47
  48 #include "gromacs/domdec/partition.h"
  49 #include "gromacs/gmxlib/nrnb.h"
  50 #include "gromacs/math/units.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/mdlib/constr.h"
  53 #include "gromacs/mdlib/coupling.h"
  54 #include "gromacs/mdlib/enerdata_utils.h"
  55 #include "gromacs/mdlib/mdatoms.h"
  56 #include "gromacs/mdlib/md_support.h"
  57 #include "gromacs/mdlib/stat.h"
  58 #include "gromacs/mdlib/tgroup.h"
  59 #include "gromacs/mdlib/update.h"
  60 #include "gromacs/mdtypes/commrec.h"
  61 #include "gromacs/mdtypes/enerdata.h"
  62 #include "gromacs/mdtypes/fcdata.h"
  63 #include "gromacs/mdtypes/forcebuffers.h"
  64 #include "gromacs/mdtypes/forcerec.h"
  65 #include "gromacs/mdtypes/group.h"
  66 #include "gromacs/mdtypes/inputrec.h"
  67 #include "gromacs/mdtypes/mdatom.h"
  68 #include "gromacs/mdtypes/state.h"
  69 #include "gromacs/pulling/pull.h"
  70 #include "gromacs/timing/wallcycle.h"
  71 #include "gromacs/topology/topology.h"
  72
  73 void integrateVVFirstStep(int64_t                   step,
  74                           bool                      bFirstStep,
  75                           bool                      bInitStep,
  76                           gmx::StartingBehavior     startingBehavior,
  77                           int                       nstglobalcomm,
  78                           t_inputrec*               ir,
  79                           t_forcerec*               fr,
  80                           t_commrec*                cr,
  81                           t_state*                  state,
  82                           t_mdatoms*                mdatoms,
  83                           const t_fcdata&           fcdata,
  84                           t_extmass*                MassQ,
  85                           t_vcm*                    vcm,
  86                           const gmx_mtop_t*         top_global,
  87                           const gmx_localtop_t&     top,
  88                           gmx_enerdata_t*           enerd,
  89                           gmx_ekindata_t*           ekind,
  90                           gmx_global_stat*          gstat,
  91                           real*                     last_ekin,
  92                           bool                      bCalcVir,
  93                           tensor                    total_vir,
  94                           tensor                    shake_vir,
  95                           tensor                    force_vir,
  96                           tensor                    pres,
  97                           matrix                    M,
  98                           bool                      do_log,
  99                           bool                      do_ene,
 100                           bool                      bCalcEner,
 101                           bool                      bGStat,
 102                           bool                      bStopCM,
 103                           bool                      bTrotter,
 104                           bool                      bExchanged,
 105                           bool*                     bSumEkinhOld,
 106                           bool*                     shouldCheckNumberOfBondedInteractions,
 107                           real*                     saved_conserved_quantity,
 108                           gmx::ForceBuffers*        f,
 109                           gmx::Update*              upd,
 110                           gmx::Constraints*         constr,
 111                           gmx::SimulationSignaller* nullSignaller,
 112                           std::array<std::vector<int>, ettTSEQMAX> trotter_seq,
 113                           t_nrnb*                                  nrnb,
 114                           const gmx::MDLogger&                     mdlog,
 115                           FILE*                                    fplog,
 116                           gmx_wallcycle*                           wcycle)
 117 {
 118     if (!bFirstStep || startingBehavior == gmx::StartingBehavior::NewSimulation)
 119     {
 120         /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
 121         rvec* vbuf = nullptr;
 122
 123         wallcycle_start(wcycle, ewcUPDATE);
 124         if (ir->eI == eiVV && bInitStep)
 125         {
 126             /* if using velocity verlet with full time step Ekin,
 127              * take the first half step only to compute the
 128              * virial for the first step. From there,
 129              * revert back to the initial coordinates
 130              * so that the input is actually the initial step.
 131              */
 132             snew(vbuf, state->natoms);
 133             copy_rvecn(state->v.rvec_array(), vbuf, 0,
 134                        state->natoms); /* should make this better for parallelizing? */
 135         }
 136         else
 137         {
 138             /* this is for NHC in the Ekin(t+dt/2) version of vv */
 139             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, MassQ, trotter_seq, ettTSEQ1);
 140         }
 141
 142         upd->update_coords(*ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind,
 143                            M, etrtVELOCITY1, cr, constr != nullptr);
 144
 145         wallcycle_stop(wcycle, ewcUPDATE);
 146         constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir);
 147         wallcycle_start(wcycle, ewcUPDATE);
 148         /* if VV, compute the pressure and constraints */
 149         /* For VV2, we strictly only need this if using pressure
 150          * control, but we really would like to have accurate pressures
 151          * printed out.
 152          * Think about ways around this in the future?
 153          * For now, keep this choice in comments.
 154          */
 155         /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
 156         /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
 157         bool bPres = TRUE;
 158         bool bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
 159         if (bCalcEner && ir->eI == eiVVAK)
 160         {
 161             *bSumEkinhOld = TRUE;
 162         }
 163         /* for vv, the first half of the integration actually corresponds to the previous step.
 164             So we need information from the last step in the first half of the integration */
 165         if (bGStat || do_per_step(step - 1, nstglobalcomm))
 166         {
 167             wallcycle_stop(wcycle, ewcUPDATE);
 168             int totalNumberOfBondedInteractions = -1;
 169             compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
 170                             makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle,
 171                             enerd, force_vir, shake_vir, total_vir, pres, constr, nullSignaller,
 172                             state->box, &totalNumberOfBondedInteractions, bSumEkinhOld,
 173                             (bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
 174                                     | (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
 175                                     | (bPres ? CGLO_CONSTRAINT : 0) | (bStopCM ? CGLO_STOPCM : 0)
 176                                     | (*shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
 177                                                                               : 0)
 178                                     | CGLO_SCALEEKIN);
 179             /* explanation of above:
 180                 a) We compute Ekin at the full time step
 181                 if 1) we are using the AveVel Ekin, and it's not the
 182                 initial step, or 2) if we are using AveEkin, but need the full
 183                 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
 184                 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
 185                 EkinAveVel because it's needed for the pressure */
 186             checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global,
 187                                             &top, makeConstArrayRef(state->x), state->box,
 188                                             shouldCheckNumberOfBondedInteractions);
 189             if (bStopCM)
 190             {
 191                 process_and_stopcm_grp(fplog, vcm, *mdatoms, makeArrayRef(state->x),
 192                                        makeArrayRef(state->v));
 193                 inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 194             }
 195             wallcycle_start(wcycle, ewcUPDATE);
 196         }
 197         /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
 198         if (!bInitStep)
 199         {
 200             if (bTrotter)
 201             {
 202                 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
 203                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, MassQ,
 204                                trotter_seq, ettTSEQ2);
 205
 206                 /* TODO This is only needed when we're about to write
 207                  * a checkpoint, because we use it after the restart
 208                  * (in a kludge?). But what should we be doing if
 209                  * the startingBehavior is NewSimulation or bInitStep are true? */
 210                 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
 211                 {
 212                     copy_mat(shake_vir, state->svir_prev);
 213                     copy_mat(force_vir, state->fvir_prev);
 214                 }
 215                 if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
 216                 {
 217                     /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
 218                     enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
 219                     enerd->term[F_EKIN] = trace(ekind->ekin);
 220                 }
 221             }
 222             else if (bExchanged)
 223             {
 224                 wallcycle_stop(wcycle, ewcUPDATE);
 225                 /* We need the kinetic energy at minus the half step for determining
 226                  * the full step kinetic energy and possibly for T-coupling.*/
 227                 /* This may not be quite working correctly yet . . . . */
 228                 compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
 229                                 makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle,
 230                                 enerd, nullptr, nullptr, nullptr, nullptr, constr, nullSignaller,
 231                                 state->box, nullptr, bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
 232                 wallcycle_start(wcycle, ewcUPDATE);
 233             }
 234         }
 235         /* if it's the initial step, we performed this first step just to get the constraint virial */
 236         if (ir->eI == eiVV && bInitStep)
 237         {
 238             copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
 239             sfree(vbuf);
 240         }
 241         wallcycle_stop(wcycle, ewcUPDATE);
 242     }
 243
 244     /* compute the conserved quantity */
 245     *saved_conserved_quantity = NPT_energy(ir, state, MassQ);
 246     if (ir->eI == eiVV)
 247     {
 248         *last_ekin = enerd->term[F_EKIN];
 249     }
 250     if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
 251     {
 252         *saved_conserved_quantity -= enerd->term[F_DISPCORR];
 253     }
 254     /* sum up the foreign kinetic energy and dK/dl terms for vv.  currently done every step so that dhdl is correct in the .edr */
 255     if (ir->efep != efepNO)
 256     {
 257         accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
 258     }
 259 }
 260
 261 void integrateVVSecondStep(int64_t                                  step,
 262                            t_inputrec*                              ir,
 263                            t_forcerec*                              fr,
 264                            t_commrec*                               cr,
 265                            t_state*                                 state,
 266                            t_mdatoms*                               mdatoms,
 267                            const t_fcdata&                          fcdata,
 268                            t_extmass*                               MassQ,
 269                            t_vcm*                                   vcm,
 270                            pull_t*                                  pull_work,
 271                            gmx_enerdata_t*                          enerd,
 272                            gmx_ekindata_t*                          ekind,
 273                            gmx_global_stat*                         gstat,
 274                            real*                                    dvdl_constr,
 275                            bool                                     bCalcVir,
 276                            tensor                                   total_vir,
 277                            tensor                                   shake_vir,
 278                            tensor                                   force_vir,
 279                            tensor                                   pres,
 280                            matrix                                   M,
 281                            matrix                                   lastbox,
 282                            bool                                     do_log,
 283                            bool                                     do_ene,
 284                            bool                                     bGStat,
 285                            bool*                                    bSumEkinhOld,
 286                            gmx::ForceBuffers*                       f,
 287                            std::vector<gmx::RVec>*                  cbuf,
 288                            gmx::Update*                             upd,
 289                            gmx::Constraints*                        constr,
 290                            gmx::SimulationSignaller*                nullSignaller,
 291                            std::array<std::vector<int>, ettTSEQMAX> trotter_seq,
 292                            t_nrnb*                                  nrnb,
 293                            gmx_wallcycle*                           wcycle)
 294 {
 295     /* velocity half-step update */
 296     upd->update_coords(*ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind, M,
 297                        etrtVELOCITY2, cr, constr != nullptr);
 298
 299
 300     /* Above, initialize just copies ekinh into ekin,
 301      * it doesn't copy position (for VV),
 302      * and entire integrator for MD.
 303      */
 304
 305     if (ir->eI == eiVVAK)
 306     {
 307         cbuf->resize(state->x.size());
 308         std::copy(state->x.begin(), state->x.end(), cbuf->begin());
 309     }
 310
 311     if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
 312     {
 313         updatePrevStepPullCom(pull_work, state);
 314     }
 315
 316     upd->update_coords(*ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind, M,
 317                        etrtPOSITION, cr, constr != nullptr);
 318
 319     wallcycle_stop(wcycle, ewcUPDATE);
 320
 321     constrain_coordinates(constr, do_log, do_ene, step, state, upd->xp()->arrayRefWithPadding(),
 322                           dvdl_constr, bCalcVir, shake_vir);
 323
 324     upd->update_sd_second_half(*ir, step, dvdl_constr, mdatoms, state, cr, nrnb, wcycle, constr,
 325                                do_log, do_ene);
 326     upd->finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
 327
 328     if (ir->eI == eiVVAK)
 329     {
 330         /* erase F_EKIN and F_TEMP here? */
 331         /* just compute the kinetic energy at the half step to perform a trotter step */
 332         compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
 333                         makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle, enerd,
 334                         force_vir, shake_vir, total_vir, pres, constr, nullSignaller, lastbox,
 335                         nullptr, bSumEkinhOld, (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE);
 336         wallcycle_start(wcycle, ewcUPDATE);
 337         trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, MassQ, trotter_seq, ettTSEQ4);
 338         /* now we know the scaling, we can compute the positions again */
 339         std::copy(cbuf->begin(), cbuf->end(), state->x.begin());
 340
 341         upd->update_coords(*ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind,
 342                            M, etrtPOSITION, cr, constr != nullptr);
 343         wallcycle_stop(wcycle, ewcUPDATE);
 344
 345         /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
 346         /* are the small terms in the shake_vir here due
 347          * to numerical errors, or are they important
 348          * physically? I'm thinking they are just errors, but not completely sure.
 349          * For now, will call without actually constraining, constr=NULL*/
 350         upd->finish_update(*ir, mdatoms, state, wcycle, false);
 351     }
 352     /* this factor or 2 correction is necessary
 353         because half of the constraint force is removed
 354         in the vv step, so we have to double it.  See
 355         the Issue #1255.  It is not yet clear
 356         if the factor of 2 is exact, or just a very
 357         good approximation, and this will be
 358         investigated.  The next step is to see if this
 359         can be done adding a dhdl contribution from the
 360         rattle step, but this is somewhat more
 361         complicated with the current code. Will be
 362         investigated, hopefully for 4.6.3. However,
 363         this current solution is much better than
 364         having it completely wrong.
 365         */
 366     enerd->term[F_DVDL_CONSTR] += 2 * *dvdl_constr;
 367 }