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