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