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