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 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/mdatom.h"
  68 #include "gromacs/mdtypes/state.h"
  69 #include "gromacs/pbcutil/pbc.h"
  70 #include "gromacs/pulling/pull.h"
  71 #include "gromacs/timing/wallcycle.h"
  72 #include "gromacs/topology/mtop_util.h"
  73 #include "gromacs/trajectory/trajectoryframe.h"
  74 #include "gromacs/utility/arrayref.h"
  75 #include "gromacs/utility/cstringutil.h"
  76 #include "gromacs/utility/fatalerror.h"
  77 #include "gromacs/utility/gmxassert.h"
  78 #include "gromacs/utility/logger.h"
  79 #include "gromacs/utility/smalloc.h"
  80 #include "gromacs/utility/snprintf.h"
  81
  82 // TODO move this to multi-sim module
  83 bool multisim_int_all_are_equal(const gmx_multisim_t* ms, int64_t value)
  84 {
  85     bool     allValuesAreEqual = true;
  86     int64_t* buf;
  87
  88     GMX_RELEASE_ASSERT(ms, "Invalid use of multi-simulation pointer");
  89
  90     snew(buf, ms->nsim);
  91     /* send our value to all other master ranks, receive all of theirs */
  92     buf[ms->sim] = value;
  93     gmx_sumli_sim(ms->nsim, buf, ms);
  94
  95     for (int s = 0; s < ms->nsim; s++)
  96     {
  97         if (buf[s] != value)
  98         {
  99             allValuesAreEqual = false;
 100             break;
 101         }
 102     }
 103
 104     sfree(buf);
 105
 106     return allValuesAreEqual;
 107 }
 108
 109 int multisim_min(const gmx_multisim_t* ms, int nmin, int n)
 110 {
 111     int*     buf;
 112     gmx_bool bPos, bEqual;
 113     int      s, d;
 114
 115     snew(buf, ms->nsim);
 116     buf[ms->sim] = n;
 117     gmx_sumi_sim(ms->nsim, buf, ms);
 118     bPos   = TRUE;
 119     bEqual = TRUE;
 120     for (s = 0; s < ms->nsim; s++)
 121     {
 122         bPos   = bPos && (buf[s] > 0);
 123         bEqual = bEqual && (buf[s] == buf[0]);
 124     }
 125     if (bPos)
 126     {
 127         if (bEqual)
 128         {
 129             nmin = std::min(nmin, buf[0]);
 130         }
 131         else
 132         {
 133             /* Find the least common multiple */
 134             for (d = 2; d < nmin; d++)
 135             {
 136                 s = 0;
 137                 while (s < ms->nsim && d % buf[s] == 0)
 138                 {
 139                     s++;
 140                 }
 141                 if (s == ms->nsim)
 142                 {
 143                     /* We found the LCM and it is less than nmin */
 144                     nmin = d;
 145                     break;
 146                 }
 147             }
 148         }
 149     }
 150     sfree(buf);
 151
 152     return nmin;
 153 }
 154
 155 /* TODO Specialize this routine into init-time and loop-time versions?
 156    e.g. bReadEkin is only true when restoring from checkpoint */
 157 void compute_globals(gmx_global_stat*               gstat,
 158                      t_commrec*                     cr,
 159                      const t_inputrec*              ir,
 160                      t_forcerec*                    fr,
 161                      gmx_ekindata_t*                ekind,
 162                      gmx::ArrayRef<const gmx::RVec> x,
 163                      gmx::ArrayRef<const gmx::RVec> v,
 164                      const matrix                   box,
 165                      real                           vdwLambda,
 166                      const t_mdatoms*               mdatoms,
 167                      t_nrnb*                        nrnb,
 168                      t_vcm*                         vcm,
 169                      gmx_wallcycle_t                wcycle,
 170                      gmx_enerdata_t*                enerd,
 171                      tensor                         force_vir,
 172                      tensor                         shake_vir,
 173                      tensor                         total_vir,
 174                      tensor                         pres,
 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, 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         for (auto& dhdl : enerd->dhdlLambda)
 276         {
 277             dhdl += enerd->dvdl_lin[efptMASS];
 278         }
 279     }
 280
 281     /* ########## Now pressure ############## */
 282     // TODO: For the VV integrator bConstrain is needed in the conditional. This is confusing, so get rid of this.
 283     if (bPres || bConstrain)
 284     {
 285         m_add(force_vir, shake_vir, total_vir);
 286
 287         /* Calculate pressure and apply LR correction if PPPM is used.
 288          * Use the box from last timestep since we already called update().
 289          */
 290
 291         enerd->term[F_PRES] = calc_pres(fr->pbcType, ir->nwall, lastbox, ekind->ekin, total_vir, pres);
 292     }
 293
 294     /* ##########  Long range energy information ###### */
 295     if ((bEner || bPres) && fr->dispersionCorrection)
 296     {
 297         /* Calculate long range corrections to pressure and energy */
 298         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
 299            and computes enerd->term[F_DISPCORR].  Also modifies the
 300            total_vir and pres tensors */
 301
 302         const DispersionCorrection::Correction correction =
 303                 fr->dispersionCorrection->calculate(lastbox, vdwLambda);
 304
 305         if (bEner)
 306         {
 307             enerd->term[F_DISPCORR] = correction.energy;
 308             enerd->term[F_EPOT] += correction.energy;
 309             enerd->term[F_DVDL_VDW] += correction.dvdl;
 310         }
 311         if (bPres)
 312         {
 313             correction.correctVirial(total_vir);
 314             correction.correctPressure(pres);
 315             enerd->term[F_PDISPCORR] = correction.pressure;
 316             enerd->term[F_PRES] += correction.pressure;
 317         }
 318     }
 319 }
 320
 321 void setCurrentLambdasRerun(int64_t           step,
 322                             const t_lambda*   fepvals,
 323                             const t_trxframe* rerun_fr,
 324                             const double*     lam0,
 325                             t_state*          globalState)
 326 {
 327     GMX_RELEASE_ASSERT(globalState != nullptr,
 328                        "setCurrentLambdasGlobalRerun should be called with a valid state object");
 329
 330     if (rerun_fr->bLambda)
 331     {
 332         if (fepvals->delta_lambda == 0)
 333         {
 334             globalState->lambda[efptFEP] = rerun_fr->lambda;
 335         }
 336         else
 337         {
 338             /* find out between which two value of lambda we should be */
 339             real frac      = step * fepvals->delta_lambda;
 340             int  fep_state = static_cast<int>(std::floor(frac * fepvals->n_lambda));
 341             /* interpolate between this state and the next */
 342             /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 343             frac = frac * fepvals->n_lambda - fep_state;
 344             for (int i = 0; i < efptNR; i++)
 345             {
 346                 globalState->lambda[i] =
 347                         lam0[i] + (fepvals->all_lambda[i][fep_state])
 348                         + frac * (fepvals->all_lambda[i][fep_state + 1] - fepvals->all_lambda[i][fep_state]);
 349             }
 350         }
 351     }
 352     else if (rerun_fr->bFepState)
 353     {
 354         globalState->fep_state = rerun_fr->fep_state;
 355         for (int i = 0; i < efptNR; i++)
 356         {
 357             globalState->lambda[i] = fepvals->all_lambda[i][globalState->fep_state];
 358         }
 359     }
 360 }
 361
 362 void setCurrentLambdasLocal(const int64_t       step,
 363                             const t_lambda*     fepvals,
 364                             const double*       lam0,
 365                             gmx::ArrayRef<real> lambda,
 366                             const int           currentFEPState)
 367 /* find the current lambdas.  If rerunning, we either read in a state, or a lambda value,
 368    requiring different logic. */
 369 {
 370     if (fepvals->delta_lambda != 0)
 371     {
 372         /* find out between which two value of lambda we should be */
 373         real frac = step * fepvals->delta_lambda;
 374         if (fepvals->n_lambda > 0)
 375         {
 376             int fep_state = static_cast<int>(std::floor(frac * fepvals->n_lambda));
 377             /* interpolate between this state and the next */
 378             /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 379             frac = frac * fepvals->n_lambda - fep_state;
 380             for (int i = 0; i < efptNR; i++)
 381             {
 382                 lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state])
 383                             + frac * (fepvals->all_lambda[i][fep_state + 1] - fepvals->all_lambda[i][fep_state]);
 384             }
 385         }
 386         else
 387         {
 388             for (int i = 0; i < efptNR; i++)
 389             {
 390                 lambda[i] = lam0[i] + frac;
 391             }
 392         }
 393     }
 394     else
 395     {
 396         /* if < 0, fep_state was never defined, and we should not set lambda from the state */
 397         if (currentFEPState > -1)
 398         {
 399             for (int i = 0; i < efptNR; i++)
 400             {
 401                 lambda[i] = fepvals->all_lambda[i][currentFEPState];
 402             }
 403         }
 404     }
 405 }
 406
 407 static void min_zero(int* n, int i)
 408 {
 409     if (i > 0 && (*n == 0 || i < *n))
 410     {
 411         *n = i;
 412     }
 413 }
 414
 415 static int lcd4(int i1, int i2, int i3, int i4)
 416 {
 417     int nst;
 418
 419     nst = 0;
 420     min_zero(&nst, i1);
 421     min_zero(&nst, i2);
 422     min_zero(&nst, i3);
 423     min_zero(&nst, i4);
 424     if (nst == 0)
 425     {
 426         gmx_incons("All 4 inputs for determining nstglobalcomm are <= 0");
 427     }
 428
 429     while (nst > 1
 430            && ((i1 > 0 && i1 % nst != 0) || (i2 > 0 && i2 % nst != 0) || (i3 > 0 && i3 % nst != 0)
 431                || (i4 > 0 && i4 % nst != 0)))
 432     {
 433         nst--;
 434     }
 435
 436     return nst;
 437 }
 438
 439 int computeGlobalCommunicationPeriod(const gmx::MDLogger& mdlog, t_inputrec* ir, const t_commrec* cr)
 440 {
 441     int nstglobalcomm;
 442     {
 443         // Set up the default behaviour
 444         if (!(ir->nstcalcenergy > 0 || ir->nstlist > 0 || ir->etc != etcNO || ir->epc != epcNO))
 445         {
 446             /* The user didn't choose the period for anything
 447                important, so we just make sure we can send signals and
 448                write output suitably. */
 449             nstglobalcomm = 10;
 450             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
 451             {
 452                 nstglobalcomm = ir->nstenergy;
 453             }
 454         }
 455         else
 456         {
 457             /* The user has made a choice (perhaps implicitly), so we
 458              * ensure that we do timely intra-simulation communication
 459              * for (possibly) each of the four parts that care.
 460              *
 461              * TODO Does the Verlet scheme (+ DD) need any
 462              * communication at nstlist steps? Is the use of nstlist
 463              * here a leftover of the twin-range scheme? Can we remove
 464              * nstlist when we remove the group scheme?
 465              */
 466             nstglobalcomm = lcd4(ir->nstcalcenergy, ir->nstlist, ir->etc != etcNO ? ir->nsttcouple : 0,
 467                                  ir->epc != epcNO ? ir->nstpcouple : 0);
 468         }
 469     }
 470
 471     // TODO change this behaviour. Instead grompp should print
 472     // a (performance) note and mdrun should not change ir.
 473     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
 474     {
 475         GMX_LOG(mdlog.warning)
 476                 .asParagraph()
 477                 .appendTextFormatted("WARNING: Changing nstcomm from %d to %d", ir->nstcomm, nstglobalcomm);
 478         ir->nstcomm = nstglobalcomm;
 479     }
 480
 481     if (cr->nnodes > 1)
 482     {
 483         GMX_LOG(mdlog.info)
 484                 .appendTextFormatted("Intra-simulation communication will occur every %d steps.\n",
 485                                      nstglobalcomm);
 486     }
 487     return nstglobalcomm;
 488 }
 489
 490 void rerun_parallel_comm(t_commrec* cr, t_trxframe* fr, gmx_bool* bLastStep)
 491 {
 492     rvec *xp, *vp;
 493
 494     if (MASTER(cr) && *bLastStep)
 495     {
 496         fr->natoms = -1;
 497     }
 498     xp = fr->x;
 499     vp = fr->v;
 500     gmx_bcast(sizeof(*fr), fr, cr);
 501     fr->x = xp;
 502     fr->v = vp;
 503
 504     *bLastStep = (fr->natoms < 0);
 505 }
 506
 507 // TODO Most of this logic seems to belong in the respective modules
 508 void set_state_entries(t_state* state, const t_inputrec* ir, bool useModularSimulator)
 509 {
 510     /* The entries in the state in the tpx file might not correspond
 511      * with what is needed, so we correct this here.
 512      */
 513     state->flags = 0;
 514     if (ir->efep != efepNO || ir->bExpanded)
 515     {
 516         state->flags |= (1 << estLAMBDA);
 517         state->flags |= (1 << estFEPSTATE);
 518     }
 519     state->flags |= (1 << estX);
 520     GMX_RELEASE_ASSERT(state->x.size() == state->natoms,
 521                        "We should start a run with an initialized state->x");
 522     if (EI_DYNAMICS(ir->eI))
 523     {
 524         state->flags |= (1 << estV);
 525     }
 526
 527     state->nnhpres = 0;
 528     if (ir->pbcType != PbcType::No)
 529     {
 530         state->flags |= (1 << estBOX);
 531         if (inputrecPreserveShape(ir))
 532         {
 533             state->flags |= (1 << estBOX_REL);
 534         }
 535         if ((ir->epc == epcPARRINELLORAHMAN) || (ir->epc == epcMTTK))
 536         {
 537             state->flags |= (1 << estBOXV);
 538             if (!useModularSimulator)
 539             {
 540                 state->flags |= (1 << estPRES_PREV);
 541             }
 542         }
 543         if (inputrecNptTrotter(ir) || (inputrecNphTrotter(ir)))
 544         {
 545             state->nnhpres = 1;
 546             state->flags |= (1 << estNHPRES_XI);
 547             state->flags |= (1 << estNHPRES_VXI);
 548             state->flags |= (1 << estSVIR_PREV);
 549             state->flags |= (1 << estFVIR_PREV);
 550             state->flags |= (1 << estVETA);
 551             state->flags |= (1 << estVOL0);
 552         }
 553         if (ir->epc == epcBERENDSEN)
 554         {
 555             state->flags |= (1 << estBAROS_INT);
 556         }
 557     }
 558
 559     if (ir->etc == etcNOSEHOOVER)
 560     {
 561         state->flags |= (1 << estNH_XI);
 562         state->flags |= (1 << estNH_VXI);
 563     }
 564
 565     if (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)
 566     {
 567         state->flags |= (1 << estTHERM_INT);
 568     }
 569
 570     init_gtc_state(state, state->ngtc, state->nnhpres,
 571                    ir->opts.nhchainlength); /* allocate the space for nose-hoover chains */
 572     init_ekinstate(&state->ekinstate, ir);
 573
 574     if (ir->bExpanded)
 575     {
 576         snew(state->dfhist, 1);
 577         init_df_history(state->dfhist, ir->fepvals->n_lambda);
 578     }
 579
 580     if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
 581     {
 582         state->flags |= (1 << estPULLCOMPREVSTEP);
 583     }
 584 }