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,2021, 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                           const 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 == IntegrationAlgorithm::VV && 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, state->natoms); /* should make this better for parallelizing? */
 134         }
 135         else
 136         {
 137             /* this is for NHC in the Ekin(t+dt/2) version of vv */
 138             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, MassQ, trotter_seq, ettTSEQ1);
 139         }
 140
 141         upd->update_coords(
 142                 *ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind, M, etrtVELOCITY1, cr, constr != nullptr);
 143
 144         wallcycle_stop(wcycle, ewcUPDATE);
 145         constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir);
 146         wallcycle_start(wcycle, ewcUPDATE);
 147         /* if VV, compute the pressure and constraints */
 148         /* For VV2, we strictly only need this if using pressure
 149          * control, but we really would like to have accurate pressures
 150          * printed out.
 151          * Think about ways around this in the future?
 152          * For now, keep this choice in comments.
 153          */
 154         /*bPres = (ir->eI==IntegrationAlgorithm::VV || inputrecNptTrotter(ir)); */
 155         /*bTemp = ((ir->eI==IntegrationAlgorithm::VV &&(!bInitStep)) || (ir->eI==IntegrationAlgorithm::VVAK && inputrecNptTrotter(ir)));*/
 156         bool bPres = TRUE;
 157         bool bTemp = ((ir->eI == IntegrationAlgorithm::VV && (!bInitStep))
 158                       || (ir->eI == IntegrationAlgorithm::VVAK));
 159         if (bCalcEner && ir->eI == IntegrationAlgorithm::VVAK)
 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,
 170                             cr,
 171                             ir,
 172                             fr,
 173                             ekind,
 174                             makeConstArrayRef(state->x),
 175                             makeConstArrayRef(state->v),
 176                             state->box,
 177                             mdatoms,
 178                             nrnb,
 179                             vcm,
 180                             wcycle,
 181                             enerd,
 182                             force_vir,
 183                             shake_vir,
 184                             total_vir,
 185                             pres,
 186                             (bCalcEner && constr != nullptr) ? constr->rmsdData() : gmx::ArrayRef<real>{},
 187                             nullSignaller,
 188                             state->box,
 189                             &totalNumberOfBondedInteractions,
 190                             bSumEkinhOld,
 191                             (bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
 192                                     | (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
 193                                     | (bPres ? CGLO_CONSTRAINT : 0) | (bStopCM ? CGLO_STOPCM : 0)
 194                                     | (*shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
 195                                                                               : 0)
 196                                     | CGLO_SCALEEKIN);
 197             /* explanation of above:
 198                 a) We compute Ekin at the full time step
 199                 if 1) we are using the AveVel Ekin, and it's not the
 200                 initial step, or 2) if we are using AveEkin, but need the full
 201                 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
 202                 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
 203                 EkinAveVel because it's needed for the pressure */
 204             checkNumberOfBondedInteractions(mdlog,
 205                                             cr,
 206                                             totalNumberOfBondedInteractions,
 207                                             *top_global,
 208                                             &top,
 209                                             makeConstArrayRef(state->x),
 210                                             state->box,
 211                                             shouldCheckNumberOfBondedInteractions);
 212             if (bStopCM)
 213             {
 214                 process_and_stopcm_grp(
 215                         fplog, vcm, *mdatoms, makeArrayRef(state->x), makeArrayRef(state->v));
 216                 inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 217             }
 218             wallcycle_start(wcycle, ewcUPDATE);
 219         }
 220         /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
 221         if (!bInitStep)
 222         {
 223             if (bTrotter)
 224             {
 225                 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
 226                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, MassQ, trotter_seq, ettTSEQ2);
 227
 228                 /* TODO This is only needed when we're about to write
 229                  * a checkpoint, because we use it after the restart
 230                  * (in a kludge?). But what should we be doing if
 231                  * the startingBehavior is NewSimulation or bInitStep are true? */
 232                 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
 233                 {
 234                     copy_mat(shake_vir, state->svir_prev);
 235                     copy_mat(force_vir, state->fvir_prev);
 236                 }
 237                 if ((inputrecNptTrotter(ir) || inputrecNvtTrotter(ir)) && ir->eI == IntegrationAlgorithm::VV)
 238                 {
 239                     /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
 240                     enerd->term[F_TEMP] = sum_ekin(
 241                             &(ir->opts), ekind, nullptr, (ir->eI == IntegrationAlgorithm::VV), FALSE);
 242                     enerd->term[F_EKIN] = trace(ekind->ekin);
 243                 }
 244             }
 245             else if (bExchanged)
 246             {
 247                 wallcycle_stop(wcycle, ewcUPDATE);
 248                 /* We need the kinetic energy at minus the half step for determining
 249                  * the full step kinetic energy and possibly for T-coupling.*/
 250                 /* This may not be quite working correctly yet . . . . */
 251                 compute_globals(gstat,
 252                                 cr,
 253                                 ir,
 254                                 fr,
 255                                 ekind,
 256                                 makeConstArrayRef(state->x),
 257                                 makeConstArrayRef(state->v),
 258                                 state->box,
 259                                 mdatoms,
 260                                 nrnb,
 261                                 vcm,
 262                                 wcycle,
 263                                 enerd,
 264                                 nullptr,
 265                                 nullptr,
 266                                 nullptr,
 267                                 nullptr,
 268                                 gmx::ArrayRef<real>{},
 269                                 nullSignaller,
 270                                 state->box,
 271                                 nullptr,
 272                                 bSumEkinhOld,
 273                                 CGLO_GSTAT | CGLO_TEMPERATURE);
 274                 wallcycle_start(wcycle, ewcUPDATE);
 275             }
 276         }
 277         /* if it's the initial step, we performed this first step just to get the constraint virial */
 278         if (ir->eI == IntegrationAlgorithm::VV && bInitStep)
 279         {
 280             copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
 281             sfree(vbuf);
 282         }
 283         wallcycle_stop(wcycle, ewcUPDATE);
 284     }
 285
 286     /* compute the conserved quantity */
 287     *saved_conserved_quantity = NPT_energy(ir, state, MassQ);
 288     if (ir->eI == IntegrationAlgorithm::VV)
 289     {
 290         *last_ekin = enerd->term[F_EKIN];
 291     }
 292     if ((ir->eDispCorr != DispersionCorrectionType::EnerPres)
 293         && (ir->eDispCorr != DispersionCorrectionType::AllEnerPres))
 294     {
 295         *saved_conserved_quantity -= enerd->term[F_DISPCORR];
 296     }
 297     /* sum up the foreign kinetic energy and dK/dl terms for vv.  currently done every step so that dhdl is correct in the .edr */
 298     if (ir->efep != FreeEnergyPerturbationType::No)
 299     {
 300         accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
 301     }
 302 }
 303
 304 void integrateVVSecondStep(int64_t                                  step,
 305                            const t_inputrec*                        ir,
 306                            t_forcerec*                              fr,
 307                            t_commrec*                               cr,
 308                            t_state*                                 state,
 309                            t_mdatoms*                               mdatoms,
 310                            const t_fcdata&                          fcdata,
 311                            t_extmass*                               MassQ,
 312                            t_vcm*                                   vcm,
 313                            pull_t*                                  pull_work,
 314                            gmx_enerdata_t*                          enerd,
 315                            gmx_ekindata_t*                          ekind,
 316                            gmx_global_stat*                         gstat,
 317                            real*                                    dvdl_constr,
 318                            bool                                     bCalcVir,
 319                            tensor                                   total_vir,
 320                            tensor                                   shake_vir,
 321                            tensor                                   force_vir,
 322                            tensor                                   pres,
 323                            matrix                                   M,
 324                            matrix                                   lastbox,
 325                            bool                                     do_log,
 326                            bool                                     do_ene,
 327                            bool                                     bGStat,
 328                            bool*                                    bSumEkinhOld,
 329                            gmx::ForceBuffers*                       f,
 330                            std::vector<gmx::RVec>*                  cbuf,
 331                            gmx::Update*                             upd,
 332                            gmx::Constraints*                        constr,
 333                            gmx::SimulationSignaller*                nullSignaller,
 334                            std::array<std::vector<int>, ettTSEQMAX> trotter_seq,
 335                            t_nrnb*                                  nrnb,
 336                            gmx_wallcycle*                           wcycle)
 337 {
 338     /* velocity half-step update */
 339     upd->update_coords(
 340             *ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind, M, etrtVELOCITY2, cr, constr != nullptr);
 341
 342
 343     /* Above, initialize just copies ekinh into ekin,
 344      * it doesn't copy position (for VV),
 345      * and entire integrator for MD.
 346      */
 347
 348     if (ir->eI == IntegrationAlgorithm::VVAK)
 349     {
 350         cbuf->resize(state->x.size());
 351         std::copy(state->x.begin(), state->x.end(), cbuf->begin());
 352     }
 353
 354     if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
 355     {
 356         updatePrevStepPullCom(pull_work, state);
 357     }
 358
 359     upd->update_coords(
 360             *ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind, M, etrtPOSITION, cr, constr != nullptr);
 361
 362     wallcycle_stop(wcycle, ewcUPDATE);
 363
 364     constrain_coordinates(
 365             constr, do_log, do_ene, step, state, upd->xp()->arrayRefWithPadding(), dvdl_constr, bCalcVir, shake_vir);
 366
 367     upd->update_sd_second_half(
 368             *ir, step, dvdl_constr, mdatoms, state, cr, nrnb, wcycle, constr, do_log, do_ene);
 369     upd->finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
 370
 371     if (ir->eI == IntegrationAlgorithm::VVAK)
 372     {
 373         /* erase F_EKIN and F_TEMP here? */
 374         /* just compute the kinetic energy at the half step to perform a trotter step */
 375         compute_globals(gstat,
 376                         cr,
 377                         ir,
 378                         fr,
 379                         ekind,
 380                         makeConstArrayRef(state->x),
 381                         makeConstArrayRef(state->v),
 382                         state->box,
 383                         mdatoms,
 384                         nrnb,
 385                         vcm,
 386                         wcycle,
 387                         enerd,
 388                         force_vir,
 389                         shake_vir,
 390                         total_vir,
 391                         pres,
 392                         gmx::ArrayRef<real>{},
 393                         nullSignaller,
 394                         lastbox,
 395                         nullptr,
 396                         bSumEkinhOld,
 397                         (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE);
 398         wallcycle_start(wcycle, ewcUPDATE);
 399         trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, MassQ, trotter_seq, ettTSEQ4);
 400         /* now we know the scaling, we can compute the positions again */
 401         std::copy(cbuf->begin(), cbuf->end(), state->x.begin());
 402
 403         upd->update_coords(
 404                 *ir, step, mdatoms, state, f->view().forceWithPadding(), fcdata, ekind, M, etrtPOSITION, cr, constr != nullptr);
 405         wallcycle_stop(wcycle, ewcUPDATE);
 406
 407         /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
 408         /* are the small terms in the shake_vir here due
 409          * to numerical errors, or are they important
 410          * physically? I'm thinking they are just errors, but not completely sure.
 411          * For now, will call without actually constraining, constr=NULL*/
 412         upd->finish_update(*ir, mdatoms, state, wcycle, false);
 413     }
 414     /* this factor or 2 correction is necessary
 415         because half of the constraint force is removed
 416         in the vv step, so we have to double it.  See
 417         the Issue #1255.  It is not yet clear
 418         if the factor of 2 is exact, or just a very
 419         good approximation, and this will be
 420         investigated.  The next step is to see if this
 421         can be done adding a dhdl contribution from the
 422         rattle step, but this is somewhat more
 423         complicated with the current code. Will be
 424         investigated, hopefully for 4.6.3. However,
 425         this current solution is much better than
 426         having it completely wrong.
 427         */
 428     enerd->term[F_DVDL_CONSTR] += 2 * *dvdl_constr;
 429 }