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