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