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