src/gromacs/mdlib/md_support.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, 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
  39 #include "gmxpre.h"
  40
  41 #include "md_support.h"
  42
  43 #include <climits>
  44 #include <cmath>
  45
  46 #include <algorithm>
  47
  48 #include "gromacs/domdec/domdec.h"
  49 #include "gromacs/gmxlib/network.h"
  50 #include "gromacs/gmxlib/nrnb.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/mdlib/dispersioncorrection.h"
  53 #include "gromacs/mdlib/simulationsignal.h"
  54 #include "gromacs/mdlib/stat.h"
  55 #include "gromacs/mdlib/tgroup.h"
  56 #include "gromacs/mdlib/update.h"
  57 #include "gromacs/mdlib/vcm.h"
  58 #include "gromacs/mdrunutility/multisim.h"
  59 #include "gromacs/mdtypes/commrec.h"
  60 #include "gromacs/mdtypes/df_history.h"
  61 #include "gromacs/mdtypes/enerdata.h"
  62 #include "gromacs/mdtypes/energyhistory.h"
  63 #include "gromacs/mdtypes/forcerec.h"
  64 #include "gromacs/mdtypes/group.h"
  65 #include "gromacs/mdtypes/inputrec.h"
  66 #include "gromacs/mdtypes/md_enums.h"
  67 #include "gromacs/mdtypes/state.h"
  68 #include "gromacs/pbcutil/pbc.h"
  69 #include "gromacs/pulling/pull.h"
  70 #include "gromacs/timing/wallcycle.h"
  71 #include "gromacs/topology/mtop_util.h"
  72 #include "gromacs/trajectory/trajectoryframe.h"
  73 #include "gromacs/utility/arrayref.h"
  74 #include "gromacs/utility/cstringutil.h"
  75 #include "gromacs/utility/fatalerror.h"
  76 #include "gromacs/utility/gmxassert.h"
  77 #include "gromacs/utility/logger.h"
  78 #include "gromacs/utility/smalloc.h"
  79 #include "gromacs/utility/snprintf.h"
  80
  81 // TODO move this to multi-sim module
  82 bool multisim_int_all_are_equal(const gmx_multisim_t* ms, int64_t value)
  83 {
  84     bool     allValuesAreEqual = true;
  85     int64_t* buf;
  86
  87     GMX_RELEASE_ASSERT(ms, "Invalid use of multi-simulation pointer");
  88
  89     snew(buf, ms->nsim);
  90     /* send our value to all other master ranks, receive all of theirs */
  91     buf[ms->sim] = value;
  92     gmx_sumli_sim(ms->nsim, buf, ms);
  93
  94     for (int s = 0; s < ms->nsim; s++)
  95     {
  96         if (buf[s] != value)
  97         {
  98             allValuesAreEqual = false;
  99             break;
 100         }
 101     }
 102
 103     sfree(buf);
 104
 105     return allValuesAreEqual;
 106 }
 107
 108 int multisim_min(const gmx_multisim_t* ms, int nmin, int n)
 109 {
 110     int*     buf;
 111     gmx_bool bPos, bEqual;
 112     int      s, d;
 113
 114     snew(buf, ms->nsim);
 115     buf[ms->sim] = n;
 116     gmx_sumi_sim(ms->nsim, buf, ms);
 117     bPos   = TRUE;
 118     bEqual = TRUE;
 119     for (s = 0; s < ms->nsim; s++)
 120     {
 121         bPos   = bPos && (buf[s] > 0);
 122         bEqual = bEqual && (buf[s] == buf[0]);
 123     }
 124     if (bPos)
 125     {
 126         if (bEqual)
 127         {
 128             nmin = std::min(nmin, buf[0]);
 129         }
 130         else
 131         {
 132             /* Find the least common multiple */
 133             for (d = 2; d < nmin; d++)
 134             {
 135                 s = 0;
 136                 while (s < ms->nsim && d % buf[s] == 0)
 137                 {
 138                     s++;
 139                 }
 140                 if (s == ms->nsim)
 141                 {
 142                     /* We found the LCM and it is less than nmin */
 143                     nmin = d;
 144                     break;
 145                 }
 146             }
 147         }
 148     }
 149     sfree(buf);
 150
 151     return nmin;
 152 }
 153
 154 /* TODO Specialize this routine into init-time and loop-time versions?
 155    e.g. bReadEkin is only true when restoring from checkpoint */
 156 void compute_globals(gmx_global_stat*          gstat,
 157                      t_commrec*                cr,
 158                      const t_inputrec*         ir,
 159                      t_forcerec*               fr,
 160                      gmx_ekindata_t*           ekind,
 161                      const rvec*               x,
 162                      const rvec*               v,
 163                      const matrix              box,
 164                      real                      vdwLambda,
 165                      const t_mdatoms*          mdatoms,
 166                      t_nrnb*                   nrnb,
 167                      t_vcm*                    vcm,
 168                      gmx_wallcycle_t           wcycle,
 169                      gmx_enerdata_t*           enerd,
 170                      tensor                    force_vir,
 171                      tensor                    shake_vir,
 172                      tensor                    total_vir,
 173                      tensor                    pres,
 174                      rvec                      mu_tot,
 175                      gmx::Constraints*         constr,
 176                      gmx::SimulationSignaller* signalCoordinator,
 177                      const matrix              lastbox,
 178                      int*                      totalNumberOfBondedInteractions,
 179                      gmx_bool*                 bSumEkinhOld,
 180                      const int                 flags)
 181 {
 182     gmx_bool bEner, bPres, bTemp;
 183     gmx_bool bStopCM, bGStat, bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
 184     gmx_bool bCheckNumberOfBondedInteractions;
 185     real     dvdl_ekin;
 186
 187     /* translate CGLO flags to gmx_booleans */
 188     bStopCM                          = ((flags & CGLO_STOPCM) != 0);
 189     bGStat                           = ((flags & CGLO_GSTAT) != 0);
 190     bReadEkin                        = ((flags & CGLO_READEKIN) != 0);
 191     bScaleEkin                       = ((flags & CGLO_SCALEEKIN) != 0);
 192     bEner                            = ((flags & CGLO_ENERGY) != 0);
 193     bTemp                            = ((flags & CGLO_TEMPERATURE) != 0);
 194     bPres                            = ((flags & CGLO_PRESSURE) != 0);
 195     bConstrain                       = ((flags & CGLO_CONSTRAINT) != 0);
 196     bCheckNumberOfBondedInteractions = ((flags & CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS) != 0);
 197
 198     /* we calculate a full state kinetic energy either with full-step velocity verlet
 199        or half step where we need the pressure */
 200
 201     bEkinAveVel = (ir->eI == eiVV || (ir->eI == eiVVAK && bPres) || bReadEkin);
 202
 203     /* in initalization, it sums the shake virial in vv, and to
 204        sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
 205
 206     /* ########## Kinetic energy  ############## */
 207
 208     if (bTemp)
 209     {
 210         /* Non-equilibrium MD: this is parallellized, but only does communication
 211          * when there really is NEMD.
 212          */
 213
 214         if (PAR(cr) && (ekind->bNEMD))
 215         {
 216             accumulate_u(cr, &(ir->opts), ekind);
 217         }
 218         if (!bReadEkin)
 219         {
 220             calc_ke_part(x, v, box, &(ir->opts), mdatoms, ekind, nrnb, bEkinAveVel);
 221         }
 222     }
 223
 224     /* Calculate center of mass velocity if necessary, also parallellized */
 225     if (bStopCM)
 226     {
 227         calc_vcm_grp(*mdatoms, x, v, vcm);
 228     }
 229
 230     if (bTemp || bStopCM || bPres || bEner || bConstrain || bCheckNumberOfBondedInteractions)
 231     {
 232         if (!bGStat)
 233         {
 234             /* We will not sum ekinh_old,
 235              * so signal that we still have to do it.
 236              */
 237             *bSumEkinhOld = TRUE;
 238         }
 239         else
 240         {
 241             gmx::ArrayRef<real> signalBuffer = signalCoordinator->getCommunicationBuffer();
 242             if (PAR(cr))
 243             {
 244                 wallcycle_start(wcycle, ewcMoveE);
 245                 global_stat(gstat, cr, enerd, force_vir, shake_vir, mu_tot, ir, ekind, constr,
 246                             bStopCM ? vcm : nullptr, signalBuffer.size(), signalBuffer.data(),
 247                             totalNumberOfBondedInteractions, *bSumEkinhOld, flags);
 248                 wallcycle_stop(wcycle, ewcMoveE);
 249             }
 250             signalCoordinator->finalizeSignals();
 251             *bSumEkinhOld = FALSE;
 252         }
 253     }
 254
 255     if (bEner)
 256     {
 257         /* Calculate the amplitude of the cosine velocity profile */
 258         ekind->cosacc.vcos = ekind->cosacc.mvcos / mdatoms->tmass;
 259     }
 260
 261     if (bTemp)
 262     {
 263         /* Sum the kinetic energies of the groups & calc temp */
 264         /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with inputrecNptTrotter */
 265         /* three maincase:  VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
 266            Leap with AveVel is not supported; it's not clear that it will actually work.
 267            bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
 268            If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
 269          */
 270         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin, bEkinAveVel, bScaleEkin);
 271         enerd->dvdl_lin[efptMASS] = static_cast<double>(dvdl_ekin);
 272
 273         enerd->term[F_EKIN] = trace(ekind->ekin);
 274     }
 275
 276     /* ########## Now pressure ############## */
 277     // TODO: For the VV integrator bConstrain is needed in the conditional. This is confusing, so get rid of this.
 278     if (bPres || bConstrain)
 279     {
 280         m_add(force_vir, shake_vir, total_vir);
 281
 282         /* Calculate pressure and apply LR correction if PPPM is used.
 283          * Use the box from last timestep since we already called update().
 284          */
 285
 286         enerd->term[F_PRES] = calc_pres(fr->ePBC, ir->nwall, lastbox, ekind->ekin, total_vir, pres);
 287     }
 288
 289     /* ##########  Long range energy information ###### */
 290     if ((bEner || bPres) && fr->dispersionCorrection)
 291     {
 292         /* Calculate long range corrections to pressure and energy */
 293         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
 294            and computes enerd->term[F_DISPCORR].  Also modifies the
 295            total_vir and pres tensors */
 296
 297         const DispersionCorrection::Correction correction =
 298                 fr->dispersionCorrection->calculate(lastbox, vdwLambda);
 299
 300         if (bEner)
 301         {
 302             enerd->term[F_DISPCORR] = correction.energy;
 303             enerd->term[F_EPOT] += correction.energy;
 304             enerd->term[F_DVDL_VDW] += correction.dvdl;
 305         }
 306         if (bPres)
 307         {
 308             correction.correctVirial(total_vir);
 309             correction.correctPressure(pres);
 310             enerd->term[F_PDISPCORR] = correction.pressure;
 311             enerd->term[F_PRES] += correction.pressure;
 312         }
 313     }
 314 }
 315
 316 void setCurrentLambdasRerun(int64_t           step,
 317                             const t_lambda*   fepvals,
 318                             const t_trxframe* rerun_fr,
 319                             const double*     lam0,
 320                             t_state*          globalState)
 321 {
 322     GMX_RELEASE_ASSERT(globalState != nullptr,
 323                        "setCurrentLambdasGlobalRerun should be called with a valid state object");
 324
 325     if (rerun_fr->bLambda)
 326     {
 327         if (fepvals->delta_lambda == 0)
 328         {
 329             globalState->lambda[efptFEP] = rerun_fr->lambda;
 330         }
 331         else
 332         {
 333             /* find out between which two value of lambda we should be */
 334             real frac      = step * fepvals->delta_lambda;
 335             int  fep_state = static_cast<int>(std::floor(frac * fepvals->n_lambda));
 336             /* interpolate between this state and the next */
 337             /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 338             frac = frac * fepvals->n_lambda - fep_state;
 339             for (int i = 0; i < efptNR; i++)
 340             {
 341                 globalState->lambda[i] =
 342                         lam0[i] + (fepvals->all_lambda[i][fep_state])
 343                         + frac * (fepvals->all_lambda[i][fep_state + 1] - fepvals->all_lambda[i][fep_state]);
 344             }
 345         }
 346     }
 347     else if (rerun_fr->bFepState)
 348     {
 349         globalState->fep_state = rerun_fr->fep_state;
 350         for (int i = 0; i < efptNR; i++)
 351         {
 352             globalState->lambda[i] = fepvals->all_lambda[i][globalState->fep_state];
 353         }
 354     }
 355 }
 356
 357 void setCurrentLambdasLocal(const int64_t       step,
 358                             const t_lambda*     fepvals,
 359                             const double*       lam0,
 360                             gmx::ArrayRef<real> lambda,
 361                             const int           currentFEPState)
 362 /* find the current lambdas.  If rerunning, we either read in a state, or a lambda value,
 363    requiring different logic. */
 364 {
 365     if (fepvals->delta_lambda != 0)
 366     {
 367         /* find out between which two value of lambda we should be */
 368         real frac = step * fepvals->delta_lambda;
 369         if (fepvals->n_lambda > 0)
 370         {
 371             int fep_state = static_cast<int>(std::floor(frac * fepvals->n_lambda));
 372             /* interpolate between this state and the next */
 373             /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 374             frac = frac * fepvals->n_lambda - fep_state;
 375             for (int i = 0; i < efptNR; i++)
 376             {
 377                 lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state])
 378                             + frac * (fepvals->all_lambda[i][fep_state + 1] - fepvals->all_lambda[i][fep_state]);
 379             }
 380         }
 381         else
 382         {
 383             for (int i = 0; i < efptNR; i++)
 384             {
 385                 lambda[i] = lam0[i] + frac;
 386             }
 387         }
 388     }
 389     else
 390     {
 391         /* if < 0, fep_state was never defined, and we should not set lambda from the state */
 392         if (currentFEPState > -1)
 393         {
 394             for (int i = 0; i < efptNR; i++)
 395             {
 396                 lambda[i] = fepvals->all_lambda[i][currentFEPState];
 397             }
 398         }
 399     }
 400 }
 401
 402 static void min_zero(int* n, int i)
 403 {
 404     if (i > 0 && (*n == 0 || i < *n))
 405     {
 406         *n = i;
 407     }
 408 }
 409
 410 static int lcd4(int i1, int i2, int i3, int i4)
 411 {
 412     int nst;
 413
 414     nst = 0;
 415     min_zero(&nst, i1);
 416     min_zero(&nst, i2);
 417     min_zero(&nst, i3);
 418     min_zero(&nst, i4);
 419     if (nst == 0)
 420     {
 421         gmx_incons("All 4 inputs for determining nstglobalcomm are <= 0");
 422     }
 423
 424     while (nst > 1
 425            && ((i1 > 0 && i1 % nst != 0) || (i2 > 0 && i2 % nst != 0) || (i3 > 0 && i3 % nst != 0)
 426                || (i4 > 0 && i4 % nst != 0)))
 427     {
 428         nst--;
 429     }
 430
 431     return nst;
 432 }
 433
 434 int computeGlobalCommunicationPeriod(const gmx::MDLogger& mdlog, t_inputrec* ir, const t_commrec* cr)
 435 {
 436     int nstglobalcomm;
 437     {
 438         // Set up the default behaviour
 439         if (!(ir->nstcalcenergy > 0 || ir->nstlist > 0 || ir->etc != etcNO || ir->epc != epcNO))
 440         {
 441             /* The user didn't choose the period for anything
 442                important, so we just make sure we can send signals and
 443                write output suitably. */
 444             nstglobalcomm = 10;
 445             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
 446             {
 447                 nstglobalcomm = ir->nstenergy;
 448             }
 449         }
 450         else
 451         {
 452             /* The user has made a choice (perhaps implicitly), so we
 453              * ensure that we do timely intra-simulation communication
 454              * for (possibly) each of the four parts that care.
 455              *
 456              * TODO Does the Verlet scheme (+ DD) need any
 457              * communication at nstlist steps? Is the use of nstlist
 458              * here a leftover of the twin-range scheme? Can we remove
 459              * nstlist when we remove the group scheme?
 460              */
 461             nstglobalcomm = lcd4(ir->nstcalcenergy, ir->nstlist, ir->etc != etcNO ? ir->nsttcouple : 0,
 462                                  ir->epc != epcNO ? ir->nstpcouple : 0);
 463         }
 464     }
 465
 466     // TODO change this behaviour. Instead grompp should print
 467     // a (performance) note and mdrun should not change ir.
 468     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
 469     {
 470         GMX_LOG(mdlog.warning)
 471                 .asParagraph()
 472                 .appendTextFormatted("WARNING: Changing nstcomm from %d to %d", ir->nstcomm, nstglobalcomm);
 473         ir->nstcomm = nstglobalcomm;
 474     }
 475
 476     if (cr->nnodes > 1)
 477     {
 478         GMX_LOG(mdlog.info)
 479                 .appendTextFormatted("Intra-simulation communication will occur every %d steps.\n",
 480                                      nstglobalcomm);
 481     }
 482     return nstglobalcomm;
 483 }
 484
 485 void rerun_parallel_comm(t_commrec* cr, t_trxframe* fr, gmx_bool* bLastStep)
 486 {
 487     rvec *xp, *vp;
 488
 489     if (MASTER(cr) && *bLastStep)
 490     {
 491         fr->natoms = -1;
 492     }
 493     xp = fr->x;
 494     vp = fr->v;
 495     gmx_bcast(sizeof(*fr), fr, cr);
 496     fr->x = xp;
 497     fr->v = vp;
 498
 499     *bLastStep = (fr->natoms < 0);
 500 }
 501
 502 // TODO Most of this logic seems to belong in the respective modules
 503 void set_state_entries(t_state* state, const t_inputrec* ir)
 504 {
 505     /* The entries in the state in the tpx file might not correspond
 506      * with what is needed, so we correct this here.
 507      */
 508     state->flags = 0;
 509     if (ir->efep != efepNO || ir->bExpanded)
 510     {
 511         state->flags |= (1 << estLAMBDA);
 512         state->flags |= (1 << estFEPSTATE);
 513     }
 514     state->flags |= (1 << estX);
 515     GMX_RELEASE_ASSERT(state->x.size() == state->natoms,
 516                        "We should start a run with an initialized state->x");
 517     if (EI_DYNAMICS(ir->eI))
 518     {
 519         state->flags |= (1 << estV);
 520     }
 521
 522     state->nnhpres = 0;
 523     if (ir->ePBC != epbcNONE)
 524     {
 525         state->flags |= (1 << estBOX);
 526         if (inputrecPreserveShape(ir))
 527         {
 528             state->flags |= (1 << estBOX_REL);
 529         }
 530         if ((ir->epc == epcPARRINELLORAHMAN) || (ir->epc == epcMTTK))
 531         {
 532             state->flags |= (1 << estBOXV);
 533             state->flags |= (1 << estPRES_PREV);
 534         }
 535         if (inputrecNptTrotter(ir) || (inputrecNphTrotter(ir)))
 536         {
 537             state->nnhpres = 1;
 538             state->flags |= (1 << estNHPRES_XI);
 539             state->flags |= (1 << estNHPRES_VXI);
 540             state->flags |= (1 << estSVIR_PREV);
 541             state->flags |= (1 << estFVIR_PREV);
 542             state->flags |= (1 << estVETA);
 543             state->flags |= (1 << estVOL0);
 544         }
 545         if (ir->epc == epcBERENDSEN)
 546         {
 547             state->flags |= (1 << estBAROS_INT);
 548         }
 549     }
 550
 551     if (ir->etc == etcNOSEHOOVER)
 552     {
 553         state->flags |= (1 << estNH_XI);
 554         state->flags |= (1 << estNH_VXI);
 555     }
 556
 557     if (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)
 558     {
 559         state->flags |= (1 << estTHERM_INT);
 560     }
 561
 562     init_gtc_state(state, state->ngtc, state->nnhpres,
 563                    ir->opts.nhchainlength); /* allocate the space for nose-hoover chains */
 564     init_ekinstate(&state->ekinstate, ir);
 565
 566     if (ir->bExpanded)
 567     {
 568         snew(state->dfhist, 1);
 569         init_df_history(state->dfhist, ir->fepvals->n_lambda);
 570     }
 571
 572     if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
 573     {
 574         state->flags |= (1 << estPULLCOMPREVSTEP);
 575     }
 576 }