2 * This file is part of the GROMACS molecular simulation package.
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,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
40 #include "md_support.h"
47 #include "gromacs/domdec/domdec.h"
48 #include "gromacs/gmxlib/network.h"
49 #include "gromacs/gmxlib/nrnb.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/mdlib/dispersioncorrection.h"
52 #include "gromacs/mdlib/mdrun.h"
53 #include "gromacs/mdlib/sim_util.h"
54 #include "gromacs/mdlib/simulationsignal.h"
55 #include "gromacs/mdlib/tgroup.h"
56 #include "gromacs/mdlib/update.h"
57 #include "gromacs/mdlib/vcm.h"
58 #include "gromacs/mdtypes/commrec.h"
59 #include "gromacs/mdtypes/df_history.h"
60 #include "gromacs/mdtypes/enerdata.h"
61 #include "gromacs/mdtypes/energyhistory.h"
62 #include "gromacs/mdtypes/forcerec.h"
63 #include "gromacs/mdtypes/group.h"
64 #include "gromacs/mdtypes/inputrec.h"
65 #include "gromacs/mdtypes/md_enums.h"
66 #include "gromacs/mdtypes/state.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/pulling/pull.h"
69 #include "gromacs/timing/wallcycle.h"
70 #include "gromacs/topology/mtop_util.h"
71 #include "gromacs/trajectory/trajectoryframe.h"
72 #include "gromacs/utility/arrayref.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/gmxassert.h"
76 #include "gromacs/utility/logger.h"
77 #include "gromacs/utility/smalloc.h"
78 #include "gromacs/utility/snprintf.h"
80 // TODO move this to multi-sim module
81 bool multisim_int_all_are_equal(const gmx_multisim_t *ms,
84 bool allValuesAreEqual = true;
87 GMX_RELEASE_ASSERT(ms, "Invalid use of multi-simulation pointer");
90 /* send our value to all other master ranks, receive all of theirs */
92 gmx_sumli_sim(ms->nsim, buf, ms);
94 for (int s = 0; s < ms->nsim; s++)
98 allValuesAreEqual = false;
105 return allValuesAreEqual;
108 int multisim_min(const gmx_multisim_t *ms, int nmin, int n)
111 gmx_bool bPos, bEqual;
116 gmx_sumi_sim(ms->nsim, buf, ms);
119 for (s = 0; s < ms->nsim; s++)
121 bPos = bPos && (buf[s] > 0);
122 bEqual = bEqual && (buf[s] == buf[0]);
128 nmin = std::min(nmin, buf[0]);
132 /* Find the least common multiple */
133 for (d = 2; d < nmin; d++)
136 while (s < ms->nsim && d % buf[s] == 0)
142 /* We found the LCM and it is less than nmin */
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(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_inputrec *ir,
157 t_forcerec *fr, gmx_ekindata_t *ekind,
158 t_state *state, t_mdatoms *mdatoms,
159 t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
160 gmx_enerdata_t *enerd, tensor force_vir, tensor shake_vir, tensor total_vir,
161 tensor pres, rvec mu_tot, gmx::Constraints *constr,
162 gmx::SimulationSignaller *signalCoordinator,
163 matrix box, int *totalNumberOfBondedInteractions,
164 gmx_bool *bSumEkinhOld, int flags)
166 tensor corr_vir, corr_pres;
167 gmx_bool bEner, bPres, bTemp;
168 gmx_bool bStopCM, bGStat,
169 bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
170 gmx_bool bCheckNumberOfBondedInteractions;
171 real prescorr, enercorr, dvdlcorr, dvdl_ekin;
173 /* translate CGLO flags to gmx_booleans */
174 bStopCM = ((flags & CGLO_STOPCM) != 0);
175 bGStat = ((flags & CGLO_GSTAT) != 0);
176 bReadEkin = ((flags & CGLO_READEKIN) != 0);
177 bScaleEkin = ((flags & CGLO_SCALEEKIN) != 0);
178 bEner = ((flags & CGLO_ENERGY) != 0);
179 bTemp = ((flags & CGLO_TEMPERATURE) != 0);
180 bPres = ((flags & CGLO_PRESSURE) != 0);
181 bConstrain = ((flags & CGLO_CONSTRAINT) != 0);
182 bCheckNumberOfBondedInteractions = ((flags & CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS) != 0);
184 /* we calculate a full state kinetic energy either with full-step velocity verlet
185 or half step where we need the pressure */
187 bEkinAveVel = (ir->eI == eiVV || (ir->eI == eiVVAK && bPres) || bReadEkin);
189 /* in initalization, it sums the shake virial in vv, and to
190 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
192 /* ########## Kinetic energy ############## */
196 /* Non-equilibrium MD: this is parallellized, but only does communication
197 * when there really is NEMD.
200 if (PAR(cr) && (ekind->bNEMD))
202 accumulate_u(cr, &(ir->opts), ekind);
206 calc_ke_part(state, &(ir->opts), mdatoms, ekind, nrnb, bEkinAveVel);
210 /* Calculate center of mass velocity if necessary, also parallellized */
213 calc_vcm_grp(0, mdatoms->homenr, mdatoms,
214 state->x.rvec_array(), state->v.rvec_array(), vcm);
217 if (bTemp || bStopCM || bPres || bEner || bConstrain || bCheckNumberOfBondedInteractions)
221 /* We will not sum ekinh_old,
222 * so signal that we still have to do it.
224 *bSumEkinhOld = TRUE;
229 gmx::ArrayRef<real> signalBuffer = signalCoordinator->getCommunicationBuffer();
232 wallcycle_start(wcycle, ewcMoveE);
233 global_stat(gstat, cr, enerd, force_vir, shake_vir, mu_tot,
234 ir, ekind, constr, bStopCM ? vcm : nullptr,
235 signalBuffer.size(), signalBuffer.data(),
236 totalNumberOfBondedInteractions,
237 *bSumEkinhOld, flags);
238 wallcycle_stop(wcycle, ewcMoveE);
240 signalCoordinator->finalizeSignals();
241 *bSumEkinhOld = FALSE;
245 /* Do center of mass motion removal */
248 check_cm_grp(fplog, vcm, ir, 1);
249 /* At initialization, do not pass x with acceleration-correction mode
250 * to avoid (incorrect) correction of the initial coordinates.
252 rvec *xPtr = nullptr;
253 if (vcm->mode == ecmANGULAR || (vcm->mode == ecmLINEAR_ACCELERATION_CORRECTION && !(flags & CGLO_INITIALIZATION)))
255 xPtr = state->x.rvec_array();
257 do_stopcm_grp(*mdatoms,
258 xPtr, state->v.rvec_array(), *vcm);
259 inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
264 /* Calculate the amplitude of the cosine velocity profile */
265 ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
270 /* Sum the kinetic energies of the groups & calc temp */
271 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with inputrecNptTrotter */
272 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
273 Leap with AveVel is not supported; it's not clear that it will actually work.
274 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
275 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
277 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin,
278 bEkinAveVel, bScaleEkin);
279 enerd->dvdl_lin[efptMASS] = static_cast<double>(dvdl_ekin);
281 enerd->term[F_EKIN] = trace(ekind->ekin);
284 /* ########## Long range energy information ###### */
286 if (bEner || bPres || bConstrain)
288 calc_dispcorr(ir, fr, box, state->lambda[efptVDW],
289 corr_pres, corr_vir, &prescorr, &enercorr, &dvdlcorr);
294 enerd->term[F_DISPCORR] = enercorr;
295 enerd->term[F_EPOT] += enercorr;
296 enerd->term[F_DVDL_VDW] += dvdlcorr;
299 /* ########## Now pressure ############## */
300 if (bPres || bConstrain)
303 m_add(force_vir, shake_vir, total_vir);
305 /* Calculate pressure and apply LR correction if PPPM is used.
306 * Use the box from last timestep since we already called update().
309 enerd->term[F_PRES] = calc_pres(fr->ePBC, ir->nwall, box, ekind->ekin, total_vir, pres);
311 /* Calculate long range corrections to pressure and energy */
312 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
313 and computes enerd->term[F_DISPCORR]. Also modifies the
314 total_vir and pres tesors */
316 m_add(total_vir, corr_vir, total_vir);
317 m_add(pres, corr_pres, pres);
318 enerd->term[F_PDISPCORR] = prescorr;
319 enerd->term[F_PRES] += prescorr;
323 /* check whether an 'nst'-style parameter p is a multiple of nst, and
324 set it to be one if not, with a warning. */
325 static void check_nst_param(const gmx::MDLogger &mdlog,
326 const char *desc_nst, int nst,
327 const char *desc_p, int *p)
329 if (*p > 0 && *p % nst != 0)
331 /* Round up to the next multiple of nst */
332 *p = ((*p)/nst + 1)*nst;
333 GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
334 "NOTE: %s changes %s to %d", desc_nst, desc_p, *p);
338 void setCurrentLambdasRerun(int64_t step, const t_lambda *fepvals,
339 const t_trxframe *rerun_fr, const double *lam0,
340 t_state *globalState)
342 GMX_RELEASE_ASSERT(globalState != nullptr, "setCurrentLambdasGlobalRerun should be called with a valid state object");
344 if (rerun_fr->bLambda)
346 if (fepvals->delta_lambda == 0)
348 globalState->lambda[efptFEP] = rerun_fr->lambda;
352 /* find out between which two value of lambda we should be */
353 real frac = step*fepvals->delta_lambda;
354 int fep_state = static_cast<int>(std::floor(frac*fepvals->n_lambda));
355 /* interpolate between this state and the next */
356 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
357 frac = frac*fepvals->n_lambda - fep_state;
358 for (int i = 0; i < efptNR; i++)
360 globalState->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
361 frac*(fepvals->all_lambda[i][fep_state+1] - fepvals->all_lambda[i][fep_state]);
365 else if (rerun_fr->bFepState)
367 globalState->fep_state = rerun_fr->fep_state;
368 for (int i = 0; i < efptNR; i++)
370 globalState->lambda[i] = fepvals->all_lambda[i][globalState->fep_state];
375 void setCurrentLambdasLocal(int64_t step, const t_lambda *fepvals,
376 const double *lam0, t_state *state)
377 /* find the current lambdas. If rerunning, we either read in a state, or a lambda value,
378 requiring different logic. */
380 if (fepvals->delta_lambda != 0)
382 /* find out between which two value of lambda we should be */
383 real frac = step*fepvals->delta_lambda;
384 if (fepvals->n_lambda > 0)
386 int fep_state = static_cast<int>(std::floor(frac*fepvals->n_lambda));
387 /* interpolate between this state and the next */
388 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
389 frac = frac*fepvals->n_lambda - fep_state;
390 for (int i = 0; i < efptNR; i++)
392 state->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
393 frac*(fepvals->all_lambda[i][fep_state + 1] - fepvals->all_lambda[i][fep_state]);
398 for (int i = 0; i < efptNR; i++)
400 state->lambda[i] = lam0[i] + frac;
406 /* if < 0, fep_state was never defined, and we should not set lambda from the state */
407 if (state->fep_state > -1)
409 for (int i = 0; i < efptNR; i++)
411 state->lambda[i] = fepvals->all_lambda[i][state->fep_state];
417 static void min_zero(int *n, int i)
419 if (i > 0 && (*n == 0 || i < *n))
425 static int lcd4(int i1, int i2, int i3, int i4)
436 gmx_incons("All 4 inputs for determining nstglobalcomm are <= 0");
439 while (nst > 1 && ((i1 > 0 && i1 % nst != 0) ||
440 (i2 > 0 && i2 % nst != 0) ||
441 (i3 > 0 && i3 % nst != 0) ||
442 (i4 > 0 && i4 % nst != 0)))
450 int check_nstglobalcomm(const gmx::MDLogger &mdlog, int nstglobalcomm, t_inputrec *ir, const t_commrec * cr)
452 if (!EI_DYNAMICS(ir->eI))
457 if (nstglobalcomm == -1)
459 // Set up the default behaviour
460 if (!(ir->nstcalcenergy > 0 ||
465 /* The user didn't choose the period for anything
466 important, so we just make sure we can send signals and
467 write output suitably. */
469 if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
471 nstglobalcomm = ir->nstenergy;
476 /* The user has made a choice (perhaps implicitly), so we
477 * ensure that we do timely intra-simulation communication
478 * for (possibly) each of the four parts that care.
480 * TODO Does the Verlet scheme (+ DD) need any
481 * communication at nstlist steps? Is the use of nstlist
482 * here a leftover of the twin-range scheme? Can we remove
483 * nstlist when we remove the group scheme?
485 nstglobalcomm = lcd4(ir->nstcalcenergy,
487 ir->etc != etcNO ? ir->nsttcouple : 0,
488 ir->epc != epcNO ? ir->nstpcouple : 0);
493 // Check that the user's choice of mdrun -gcom will work
494 if (ir->nstlist > 0 &&
495 nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
497 nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
498 GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
499 "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d",
502 if (ir->nstcalcenergy > 0)
504 check_nst_param(mdlog, "-gcom", nstglobalcomm,
505 "nstcalcenergy", &ir->nstcalcenergy);
507 if (ir->etc != etcNO && ir->nsttcouple > 0)
509 check_nst_param(mdlog, "-gcom", nstglobalcomm,
510 "nsttcouple", &ir->nsttcouple);
512 if (ir->epc != epcNO && ir->nstpcouple > 0)
514 check_nst_param(mdlog, "-gcom", nstglobalcomm,
515 "nstpcouple", &ir->nstpcouple);
518 check_nst_param(mdlog, "-gcom", nstglobalcomm,
519 "nstenergy", &ir->nstenergy);
521 check_nst_param(mdlog, "-gcom", nstglobalcomm,
522 "nstlog", &ir->nstlog);
525 if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
527 GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
528 "WARNING: Changing nstcomm from %d to %d",
529 ir->nstcomm, nstglobalcomm);
530 ir->nstcomm = nstglobalcomm;
535 GMX_LOG(mdlog.info).appendTextFormatted(
536 "Intra-simulation communication will occur every %d steps.\n", nstglobalcomm);
538 return nstglobalcomm;
542 void rerun_parallel_comm(t_commrec *cr, t_trxframe *fr,
547 if (MASTER(cr) && *bLastStep)
553 gmx_bcast(sizeof(*fr), fr, cr);
557 *bLastStep = (fr->natoms < 0);
561 // TODO Most of this logic seems to belong in the respective modules
562 void set_state_entries(t_state *state, const t_inputrec *ir)
564 /* The entries in the state in the tpx file might not correspond
565 * with what is needed, so we correct this here.
568 if (ir->efep != efepNO || ir->bExpanded)
570 state->flags |= (1<<estLAMBDA);
571 state->flags |= (1<<estFEPSTATE);
573 state->flags |= (1<<estX);
574 GMX_RELEASE_ASSERT(state->x.size() == state->natoms, "We should start a run with an initialized state->x");
575 if (EI_DYNAMICS(ir->eI))
577 state->flags |= (1<<estV);
581 if (ir->ePBC != epbcNONE)
583 state->flags |= (1<<estBOX);
584 if (inputrecPreserveShape(ir))
586 state->flags |= (1<<estBOX_REL);
588 if ((ir->epc == epcPARRINELLORAHMAN) || (ir->epc == epcMTTK))
590 state->flags |= (1<<estBOXV);
591 state->flags |= (1<<estPRES_PREV);
593 if (inputrecNptTrotter(ir) || (inputrecNphTrotter(ir)))
596 state->flags |= (1<<estNHPRES_XI);
597 state->flags |= (1<<estNHPRES_VXI);
598 state->flags |= (1<<estSVIR_PREV);
599 state->flags |= (1<<estFVIR_PREV);
600 state->flags |= (1<<estVETA);
601 state->flags |= (1<<estVOL0);
603 if (ir->epc == epcBERENDSEN)
605 state->flags |= (1<<estBAROS_INT);
609 if (ir->etc == etcNOSEHOOVER)
611 state->flags |= (1<<estNH_XI);
612 state->flags |= (1<<estNH_VXI);
615 if (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)
617 state->flags |= (1<<estTHERM_INT);
620 init_gtc_state(state, state->ngtc, state->nnhpres, ir->opts.nhchainlength); /* allocate the space for nose-hoover chains */
621 init_ekinstate(&state->ekinstate, ir);
625 snew(state->dfhist, 1);
626 init_df_history(state->dfhist, ir->fepvals->n_lambda);
629 if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
631 state->flags |= (1<<estPULLCOMPREVSTEP);