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
  44 #include <algorithm>
  45
  46 #include "gromacs/domdec/domdec.h"
  47 #include "gromacs/gmxlib/nrnb.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/mdlib/constr.h"
  50 #include "gromacs/mdlib/coupling.h"
  51 #include "gromacs/mdlib/enerdata_utils.h"
  52 #include "gromacs/mdlib/mdatoms.h"
  53 #include "gromacs/mdlib/md_support.h"
  54 #include "gromacs/mdlib/stat.h"
  55 #include "gromacs/mdlib/tgroup.h"
  56 #include "gromacs/mdlib/update.h"
  57 #include "gromacs/mdrunutility/handlerestart.h"
  58 #include "gromacs/mdtypes/commrec.h"
  59 #include "gromacs/mdtypes/enerdata.h"
  60 #include "gromacs/mdtypes/fcdata.h"
  61 #include "gromacs/mdtypes/forcebuffers.h"
  62 #include "gromacs/mdtypes/forcerec.h"
  63 #include "gromacs/mdtypes/group.h"
  64 #include "gromacs/mdtypes/inputrec.h"
  65 #include "gromacs/mdtypes/mdatom.h"
  66 #include "gromacs/mdtypes/state.h"
  67 #include "gromacs/pulling/pull.h"
  68 #include "gromacs/timing/wallcycle.h"
  69 #include "gromacs/topology/topology.h"
  70
  71 void integrateVVFirstStep(int64_t                   step,
  72                           bool                      bFirstStep,
  73                           bool                      bInitStep,
  74                           gmx::StartingBehavior     startingBehavior,
  75                           int                       nstglobalcomm,
  76                           const t_inputrec*         ir,
  77                           t_forcerec*               fr,
  78                           t_commrec*                cr,
  79                           t_state*                  state,
  80                           t_mdatoms*                mdatoms,
  81                           const t_fcdata&           fcdata,
  82                           t_extmass*                MassQ,
  83                           t_vcm*                    vcm,
  84                           const gmx_mtop_t&         top_global,
  85                           const gmx_localtop_t&     top,
  86                           gmx_enerdata_t*           enerd,
  87                           gmx_ekindata_t*           ekind,
  88                           gmx_global_stat*          gstat,
  89                           real*                     last_ekin,
  90                           bool                      bCalcVir,
  91                           tensor                    total_vir,
  92                           tensor                    shake_vir,
  93                           tensor                    force_vir,
  94                           tensor                    pres,
  95                           matrix                    M,
  96                           bool                      do_log,
  97                           bool                      do_ene,
  98                           bool                      bCalcEner,
  99                           bool                      bGStat,
 100                           bool                      bStopCM,
 101                           bool                      bTrotter,
 102                           bool                      bExchanged,
 103                           bool*                     bSumEkinhOld,
 104                           real*                     saved_conserved_quantity,
 105                           gmx::ForceBuffers*        f,
 106                           gmx::Update*              upd,
 107                           gmx::Constraints*         constr,
 108                           gmx::SimulationSignaller* nullSignaller,
 109                           gmx::EnumerationArray<TrotterSequence, std::vector<int>> trotter_seq,
 110                           t_nrnb*                                                  nrnb,
 111                           const gmx::MDLogger&                                     mdlog,
 112                           FILE*                                                    fplog,
 113                           gmx_wallcycle*                                           wcycle)
 114 {
 115     if (!bFirstStep || startingBehavior == gmx::StartingBehavior::NewSimulation)
 116     {
 117         /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
 118         rvec* vbuf = nullptr;
 119
 120         wallcycle_start(wcycle, WallCycleCounter::Update);
 121         if (ir->eI == IntegrationAlgorithm::VV && bInitStep)
 122         {
 123             /* if using velocity verlet with full time step Ekin,
 124              * take the first half step only to compute the
 125              * virial for the first step. From there,
 126              * revert back to the initial coordinates
 127              * so that the input is actually the initial step.
 128              */
 129             snew(vbuf, state->natoms);
 130             copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
 131         }
 132         else
 133         {
 134             /* this is for NHC in the Ekin(t+dt/2) version of vv */
 135             trotter_update(ir,
 136                            step,
 137                            ekind,
 138                            enerd,
 139                            state,
 140                            total_vir,
 141                            mdatoms->homenr,
 142                            mdatoms->cTC ? gmx::arrayRefFromArray(mdatoms->cTC, mdatoms->nr)
 143                                         : gmx::ArrayRef<const unsigned short>(),
 144                            gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
 145                            MassQ,
 146                            trotter_seq,
 147                            TrotterSequence::One);
 148         }
 149
 150         upd->update_coords(*ir,
 151                            step,
 152                            mdatoms->homenr,
 153                            mdatoms->havePartiallyFrozenAtoms,
 154                            gmx::arrayRefFromArray(mdatoms->ptype, mdatoms->nr),
 155                            gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
 156                            gmx::arrayRefFromArray(mdatoms->invMassPerDim, mdatoms->nr),
 157                            state,
 158                            f->view().forceWithPadding(),
 159                            fcdata,
 160                            ekind,
 161                            M,
 162                            etrtVELOCITY1,
 163                            cr,
 164                            constr != nullptr);
 165
 166         wallcycle_stop(wcycle, WallCycleCounter::Update);
 167         constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir);
 168         wallcycle_start(wcycle, WallCycleCounter::Update);
 169         /* if VV, compute the pressure and constraints */
 170         /* For VV2, we strictly only need this if using pressure
 171          * control, but we really would like to have accurate pressures
 172          * printed out.
 173          * Think about ways around this in the future?
 174          * For now, keep this choice in comments.
 175          */
 176         /*bPres = (ir->eI==IntegrationAlgorithm::VV || inputrecNptTrotter(ir)); */
 177         /*bTemp = ((ir->eI==IntegrationAlgorithm::VV &&(!bInitStep)) || (ir->eI==IntegrationAlgorithm::VVAK && inputrecNptTrotter(ir)));*/
 178         bool bPres = TRUE;
 179         bool bTemp = ((ir->eI == IntegrationAlgorithm::VV && (!bInitStep))
 180                       || (ir->eI == IntegrationAlgorithm::VVAK));
 181         if (bCalcEner && ir->eI == IntegrationAlgorithm::VVAK)
 182         {
 183             *bSumEkinhOld = TRUE;
 184         }
 185         /* for vv, the first half of the integration actually corresponds to the previous step.
 186             So we need information from the last step in the first half of the integration */
 187         if (bGStat || do_per_step(step - 1, nstglobalcomm))
 188         {
 189             wallcycle_stop(wcycle, WallCycleCounter::Update);
 190             int cglo_flags =
 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) | CGLO_SCALEEKIN);
 194             if (DOMAINDECOMP(cr) && shouldCheckNumberOfBondedInteractions(*cr->dd))
 195             {
 196                 cglo_flags |= CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS;
 197             }
 198             compute_globals(gstat,
 199                             cr,
 200                             ir,
 201                             fr,
 202                             ekind,
 203                             makeConstArrayRef(state->x),
 204                             makeConstArrayRef(state->v),
 205                             state->box,
 206                             mdatoms,
 207                             nrnb,
 208                             vcm,
 209                             wcycle,
 210                             enerd,
 211                             force_vir,
 212                             shake_vir,
 213                             total_vir,
 214                             pres,
 215                             (bCalcEner && constr != nullptr) ? constr->rmsdData() : gmx::ArrayRef<real>{},
 216                             nullSignaller,
 217                             state->box,
 218                             bSumEkinhOld,
 219                             cglo_flags);
 220             /* explanation of above:
 221                 a) We compute Ekin at the full time step
 222                 if 1) we are using the AveVel Ekin, and it's not the
 223                 initial step, or 2) if we are using AveEkin, but need the full
 224                 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
 225                 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
 226                 EkinAveVel because it's needed for the pressure */
 227             if (DOMAINDECOMP(cr))
 228             {
 229                 checkNumberOfBondedInteractions(
 230                         mdlog, cr, top_global, &top, makeConstArrayRef(state->x), state->box);
 231             }
 232             if (bStopCM)
 233             {
 234                 process_and_stopcm_grp(
 235                         fplog, vcm, *mdatoms, makeArrayRef(state->x), makeArrayRef(state->v));
 236                 inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 237             }
 238             wallcycle_start(wcycle, WallCycleCounter::Update);
 239         }
 240         /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
 241         if (!bInitStep)
 242         {
 243             if (bTrotter)
 244             {
 245                 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
 246                 trotter_update(ir,
 247                                step,
 248                                ekind,
 249                                enerd,
 250                                state,
 251                                total_vir,
 252                                mdatoms->homenr,
 253                                mdatoms->cTC ? gmx::arrayRefFromArray(mdatoms->cTC, mdatoms->nr)
 254                                             : gmx::ArrayRef<const unsigned short>(),
 255                                gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
 256                                MassQ,
 257                                trotter_seq,
 258                                TrotterSequence::Two);
 259
 260                 /* TODO This is only needed when we're about to write
 261                  * a checkpoint, because we use it after the restart
 262                  * (in a kludge?). But what should we be doing if
 263                  * the startingBehavior is NewSimulation or bInitStep are true? */
 264                 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
 265                 {
 266                     copy_mat(shake_vir, state->svir_prev);
 267                     copy_mat(force_vir, state->fvir_prev);
 268                 }
 269                 if ((inputrecNptTrotter(ir) || inputrecNvtTrotter(ir)) && ir->eI == IntegrationAlgorithm::VV)
 270                 {
 271                     /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
 272                     enerd->term[F_TEMP] = sum_ekin(
 273                             &(ir->opts), ekind, nullptr, (ir->eI == IntegrationAlgorithm::VV), FALSE);
 274                     enerd->term[F_EKIN] = trace(ekind->ekin);
 275                 }
 276             }
 277             else if (bExchanged)
 278             {
 279                 wallcycle_stop(wcycle, WallCycleCounter::Update);
 280                 /* We need the kinetic energy at minus the half step for determining
 281                  * the full step kinetic energy and possibly for T-coupling.*/
 282                 /* This may not be quite working correctly yet . . . . */
 283                 compute_globals(gstat,
 284                                 cr,
 285                                 ir,
 286                                 fr,
 287                                 ekind,
 288                                 makeConstArrayRef(state->x),
 289                                 makeConstArrayRef(state->v),
 290                                 state->box,
 291                                 mdatoms,
 292                                 nrnb,
 293                                 vcm,
 294                                 wcycle,
 295                                 enerd,
 296                                 nullptr,
 297                                 nullptr,
 298                                 nullptr,
 299                                 nullptr,
 300                                 gmx::ArrayRef<real>{},
 301                                 nullSignaller,
 302                                 state->box,
 303                                 bSumEkinhOld,
 304                                 CGLO_GSTAT | CGLO_TEMPERATURE);
 305                 wallcycle_start(wcycle, WallCycleCounter::Update);
 306             }
 307         }
 308         /* if it's the initial step, we performed this first step just to get the constraint virial */
 309         if (ir->eI == IntegrationAlgorithm::VV && bInitStep)
 310         {
 311             copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
 312             sfree(vbuf);
 313         }
 314         wallcycle_stop(wcycle, WallCycleCounter::Update);
 315     }
 316
 317     /* compute the conserved quantity */
 318     *saved_conserved_quantity = NPT_energy(ir, state, MassQ);
 319     if (ir->eI == IntegrationAlgorithm::VV)
 320     {
 321         *last_ekin = enerd->term[F_EKIN];
 322     }
 323     if ((ir->eDispCorr != DispersionCorrectionType::EnerPres)
 324         && (ir->eDispCorr != DispersionCorrectionType::AllEnerPres))
 325     {
 326         *saved_conserved_quantity -= enerd->term[F_DISPCORR];
 327     }
 328     /* sum up the foreign kinetic energy and dK/dl terms for vv.  currently done every step so that dhdl is correct in the .edr */
 329     if (ir->efep != FreeEnergyPerturbationType::No)
 330     {
 331         accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
 332     }
 333 }
 334
 335 void integrateVVSecondStep(int64_t                                                  step,
 336                            const t_inputrec*                                        ir,
 337                            t_forcerec*                                              fr,
 338                            t_commrec*                                               cr,
 339                            t_state*                                                 state,
 340                            t_mdatoms*                                               mdatoms,
 341                            const t_fcdata&                                          fcdata,
 342                            t_extmass*                                               MassQ,
 343                            t_vcm*                                                   vcm,
 344                            pull_t*                                                  pull_work,
 345                            gmx_enerdata_t*                                          enerd,
 346                            gmx_ekindata_t*                                          ekind,
 347                            gmx_global_stat*                                         gstat,
 348                            real*                                                    dvdl_constr,
 349                            bool                                                     bCalcVir,
 350                            tensor                                                   total_vir,
 351                            tensor                                                   shake_vir,
 352                            tensor                                                   force_vir,
 353                            tensor                                                   pres,
 354                            matrix                                                   M,
 355                            matrix                                                   lastbox,
 356                            bool                                                     do_log,
 357                            bool                                                     do_ene,
 358                            bool                                                     bGStat,
 359                            bool*                                                    bSumEkinhOld,
 360                            gmx::ForceBuffers*                                       f,
 361                            std::vector<gmx::RVec>*                                  cbuf,
 362                            gmx::Update*                                             upd,
 363                            gmx::Constraints*                                        constr,
 364                            gmx::SimulationSignaller*                                nullSignaller,
 365                            gmx::EnumerationArray<TrotterSequence, std::vector<int>> trotter_seq,
 366                            t_nrnb*                                                  nrnb,
 367                            gmx_wallcycle*                                           wcycle)
 368 {
 369     /* velocity half-step update */
 370     upd->update_coords(*ir,
 371                        step,
 372                        mdatoms->homenr,
 373                        mdatoms->havePartiallyFrozenAtoms,
 374                        gmx::arrayRefFromArray(mdatoms->ptype, mdatoms->nr),
 375                        gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
 376                        gmx::arrayRefFromArray(mdatoms->invMassPerDim, mdatoms->nr),
 377                        state,
 378                        f->view().forceWithPadding(),
 379                        fcdata,
 380                        ekind,
 381                        M,
 382                        etrtVELOCITY2,
 383                        cr,
 384                        constr != nullptr);
 385
 386
 387     /* Above, initialize just copies ekinh into ekin,
 388      * it doesn't copy position (for VV),
 389      * and entire integrator for MD.
 390      */
 391
 392     if (ir->eI == IntegrationAlgorithm::VVAK)
 393     {
 394         cbuf->resize(state->x.size());
 395         std::copy(state->x.begin(), state->x.end(), cbuf->begin());
 396     }
 397
 398     if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
 399     {
 400         updatePrevStepPullCom(pull_work, state);
 401     }
 402
 403     upd->update_coords(*ir,
 404                        step,
 405                        mdatoms->homenr,
 406                        mdatoms->havePartiallyFrozenAtoms,
 407                        gmx::arrayRefFromArray(mdatoms->ptype, mdatoms->nr),
 408                        gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
 409                        gmx::arrayRefFromArray(mdatoms->invMassPerDim, mdatoms->nr),
 410                        state,
 411                        f->view().forceWithPadding(),
 412                        fcdata,
 413                        ekind,
 414                        M,
 415                        etrtPOSITION,
 416                        cr,
 417                        constr != nullptr);
 418
 419     wallcycle_stop(wcycle, WallCycleCounter::Update);
 420
 421     constrain_coordinates(
 422             constr, do_log, do_ene, step, state, upd->xp()->arrayRefWithPadding(), dvdl_constr, bCalcVir, shake_vir);
 423
 424     upd->update_sd_second_half(*ir,
 425                                step,
 426                                dvdl_constr,
 427                                mdatoms->homenr,
 428                                gmx::arrayRefFromArray(mdatoms->ptype, mdatoms->nr),
 429                                gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
 430                                state,
 431                                cr,
 432                                nrnb,
 433                                wcycle,
 434                                constr,
 435                                do_log,
 436                                do_ene);
 437     upd->finish_update(
 438             *ir, mdatoms->havePartiallyFrozenAtoms, mdatoms->homenr, state, wcycle, constr != nullptr);
 439
 440     if (ir->eI == IntegrationAlgorithm::VVAK)
 441     {
 442         /* erase F_EKIN and F_TEMP here? */
 443         /* just compute the kinetic energy at the half step to perform a trotter step */
 444         compute_globals(gstat,
 445                         cr,
 446                         ir,
 447                         fr,
 448                         ekind,
 449                         makeConstArrayRef(state->x),
 450                         makeConstArrayRef(state->v),
 451                         state->box,
 452                         mdatoms,
 453                         nrnb,
 454                         vcm,
 455                         wcycle,
 456                         enerd,
 457                         force_vir,
 458                         shake_vir,
 459                         total_vir,
 460                         pres,
 461                         gmx::ArrayRef<real>{},
 462                         nullSignaller,
 463                         lastbox,
 464                         bSumEkinhOld,
 465                         (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE);
 466         wallcycle_start(wcycle, WallCycleCounter::Update);
 467         trotter_update(ir,
 468                        step,
 469                        ekind,
 470                        enerd,
 471                        state,
 472                        total_vir,
 473                        mdatoms->homenr,
 474                        mdatoms->cTC ? gmx::arrayRefFromArray(mdatoms->cTC, mdatoms->nr)
 475                                     : gmx::ArrayRef<const unsigned short>(),
 476                        gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
 477                        MassQ,
 478                        trotter_seq,
 479                        TrotterSequence::Four);
 480         /* now we know the scaling, we can compute the positions again */
 481         std::copy(cbuf->begin(), cbuf->end(), state->x.begin());
 482
 483         upd->update_coords(*ir,
 484                            step,
 485                            mdatoms->homenr,
 486                            mdatoms->havePartiallyFrozenAtoms,
 487                            gmx::arrayRefFromArray(mdatoms->ptype, mdatoms->nr),
 488                            gmx::arrayRefFromArray(mdatoms->invmass, mdatoms->nr),
 489                            gmx::arrayRefFromArray(mdatoms->invMassPerDim, mdatoms->nr),
 490                            state,
 491                            f->view().forceWithPadding(),
 492                            fcdata,
 493                            ekind,
 494                            M,
 495                            etrtPOSITION,
 496                            cr,
 497                            constr != nullptr);
 498         wallcycle_stop(wcycle, WallCycleCounter::Update);
 499
 500         /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
 501         /* are the small terms in the shake_vir here due
 502          * to numerical errors, or are they important
 503          * physically? I'm thinking they are just errors, but not completely sure.
 504          * For now, will call without actually constraining, constr=NULL*/
 505         upd->finish_update(*ir, mdatoms->havePartiallyFrozenAtoms, mdatoms->homenr, state, wcycle, false);
 506     }
 507     /* this factor or 2 correction is necessary
 508         because half of the constraint force is removed
 509         in the vv step, so we have to double it.  See
 510         the Issue #1255.  It is not yet clear
 511         if the factor of 2 is exact, or just a very
 512         good approximation, and this will be
 513         investigated.  The next step is to see if this
 514         can be done adding a dhdl contribution from the
 515         rattle step, but this is somewhat more
 516         complicated with the current code. Will be
 517         investigated, hopefully for 4.6.3. However,
 518         this current solution is much better than
 519         having it completely wrong.
 520         */
 521     enerd->term[F_DVDL_CONSTR] += 2 * *dvdl_constr;
 522 }