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