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