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