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 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.
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.
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.
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.
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.
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.
41 #include "md_support.h"
48 #include "gromacs/domdec/domdec.h"
49 #include "gromacs/gmxlib/network.h"
50 #include "gromacs/gmxlib/nrnb.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/math/utilities.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdlib/coupling.h"
55 #include "gromacs/mdlib/gmx_omp_nthreads.h"
56 #include "gromacs/mdlib/simulationsignal.h"
57 #include "gromacs/mdlib/stat.h"
58 #include "gromacs/mdlib/tgroup.h"
59 #include "gromacs/mdlib/update.h"
60 #include "gromacs/mdlib/vcm.h"
61 #include "gromacs/mdtypes/commrec.h"
62 #include "gromacs/mdtypes/df_history.h"
63 #include "gromacs/mdtypes/enerdata.h"
64 #include "gromacs/mdtypes/energyhistory.h"
65 #include "gromacs/mdtypes/forcerec.h"
66 #include "gromacs/mdtypes/group.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/md_enums.h"
69 #include "gromacs/mdtypes/mdatom.h"
70 #include "gromacs/mdtypes/observablesreducer.h"
71 #include "gromacs/mdtypes/state.h"
72 #include "gromacs/pbcutil/pbc.h"
73 #include "gromacs/pulling/pull.h"
74 #include "gromacs/timing/wallcycle.h"
75 #include "gromacs/topology/mtop_util.h"
76 #include "gromacs/trajectory/trajectoryframe.h"
77 #include "gromacs/utility/arrayref.h"
78 #include "gromacs/utility/cstringutil.h"
79 #include "gromacs/utility/fatalerror.h"
80 #include "gromacs/utility/gmxassert.h"
81 #include "gromacs/utility/logger.h"
82 #include "gromacs/utility/smalloc.h"
83 #include "gromacs/utility/snprintf.h"
85 static void calc_ke_part_normal(gmx::ArrayRef<const gmx::RVec> v,
86 const t_grpopts* opts,
88 gmx_ekindata_t* ekind,
93 gmx::ArrayRef<t_grp_tcstat> tcstat = ekind->tcstat;
95 /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin. Leap with AveVel is also
96 an option, but not supported now.
97 bEkinAveVel: If TRUE, we sum into ekin, if FALSE, into ekinh.
100 // Now accumulate the partial global and groups ekin.
101 for (g = 0; (g < opts->ngtc); g++)
103 copy_mat(tcstat[g].ekinh, tcstat[g].ekinh_old);
106 clear_mat(tcstat[g].ekinf);
107 tcstat[g].ekinscalef_nhc = 1.0; /* need to clear this -- logic is complicated! */
111 clear_mat(tcstat[g].ekinh);
114 ekind->dekindl_old = ekind->dekindl;
115 int nthread = gmx_omp_nthreads_get(ModuleMultiThread::Update);
117 #pragma omp parallel for num_threads(nthread) schedule(static)
118 for (int thread = 0; thread < nthread; thread++)
120 // This OpenMP only loops over arrays and does not call any functions
121 // or memory allocation. It should not be able to throw, so for now
122 // we do not need a try/catch wrapper.
123 int start_t, end_t, n;
130 start_t = ((thread + 0) * md->homenr) / nthread;
131 end_t = ((thread + 1) * md->homenr) / nthread;
133 ekin_sum = ekind->ekin_work[thread];
134 dekindl_sum = ekind->dekindl_work[thread];
136 for (gt = 0; gt < opts->ngtc; gt++)
138 clear_mat(ekin_sum[gt]);
143 for (n = start_t; n < end_t; n++)
149 hm = 0.5 * md->massT[n];
151 for (d = 0; (d < DIM); d++)
153 for (m = 0; (m < DIM); m++)
155 /* if we're computing a full step velocity, v[d] has v(t). Otherwise, v(t+dt/2) */
156 ekin_sum[gt][m][d] += hm * v[n][m] * v[n][d];
159 if (md->nMassPerturbed && md->bPerturbed[n])
161 *dekindl_sum += 0.5 * (md->massB[n] - md->massA[n]) * iprod(v[n], v[n]);
167 for (int thread = 0; thread < nthread; thread++)
169 for (g = 0; g < opts->ngtc; g++)
173 m_add(tcstat[g].ekinf, ekind->ekin_work[thread][g], tcstat[g].ekinf);
177 m_add(tcstat[g].ekinh, ekind->ekin_work[thread][g], tcstat[g].ekinh);
181 ekind->dekindl += *ekind->dekindl_work[thread];
184 inc_nrnb(nrnb, eNR_EKIN, md->homenr);
187 static void calc_ke_part_visc(const matrix box,
188 gmx::ArrayRef<const gmx::RVec> x,
189 gmx::ArrayRef<const gmx::RVec> v,
190 const t_grpopts* opts,
192 gmx_ekindata_t* ekind,
194 gmx_bool bEkinAveVel)
196 int start = 0, homenr = md->homenr;
197 int g, d, n, m, gt = 0;
200 gmx::ArrayRef<t_grp_tcstat> tcstat = ekind->tcstat;
201 t_cos_acc* cosacc = &(ekind->cosacc);
206 for (g = 0; g < opts->ngtc; g++)
208 copy_mat(ekind->tcstat[g].ekinh, ekind->tcstat[g].ekinh_old);
209 clear_mat(ekind->tcstat[g].ekinh);
211 ekind->dekindl_old = ekind->dekindl;
213 fac = 2 * M_PI / box[ZZ][ZZ];
216 for (n = start; n < start + homenr; n++)
222 hm = 0.5 * md->massT[n];
224 /* Note that the times of x and v differ by half a step */
225 /* MRS -- would have to be changed for VV */
226 cosz = std::cos(fac * x[n][ZZ]);
227 /* Calculate the amplitude of the new velocity profile */
228 mvcos += 2 * cosz * md->massT[n] * v[n][XX];
230 copy_rvec(v[n], v_corrt);
231 /* Subtract the profile for the kinetic energy */
232 v_corrt[XX] -= cosz * cosacc->vcos;
233 for (d = 0; (d < DIM); d++)
235 for (m = 0; (m < DIM); m++)
237 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
240 tcstat[gt].ekinf[m][d] += hm * v_corrt[m] * v_corrt[d];
244 tcstat[gt].ekinh[m][d] += hm * v_corrt[m] * v_corrt[d];
248 if (md->nPerturbed && md->bPerturbed[n])
250 /* The minus sign here might be confusing.
251 * The kinetic contribution from dH/dl doesn't come from
252 * d m(l)/2 v^2 / dl, but rather from d p^2/2m(l) / dl,
253 * where p are the momenta. The difference is only a minus sign.
255 dekindl -= 0.5 * (md->massB[n] - md->massA[n]) * iprod(v_corrt, v_corrt);
258 ekind->dekindl = dekindl;
259 cosacc->mvcos = mvcos;
261 inc_nrnb(nrnb, eNR_EKIN, homenr);
264 static void calc_ke_part(gmx::ArrayRef<const gmx::RVec> x,
265 gmx::ArrayRef<const gmx::RVec> v,
267 const t_grpopts* opts,
269 gmx_ekindata_t* ekind,
271 gmx_bool bEkinAveVel)
273 if (ekind->cosacc.cos_accel == 0)
275 calc_ke_part_normal(v, opts, md, ekind, nrnb, bEkinAveVel);
279 calc_ke_part_visc(box, x, v, opts, md, ekind, nrnb, bEkinAveVel);
283 /* TODO Specialize this routine into init-time and loop-time versions?
284 e.g. bReadEkin is only true when restoring from checkpoint */
285 void compute_globals(gmx_global_stat* gstat,
287 const t_inputrec* ir,
289 gmx_ekindata_t* ekind,
290 gmx::ArrayRef<const gmx::RVec> x,
291 gmx::ArrayRef<const gmx::RVec> v,
293 const t_mdatoms* mdatoms,
296 gmx_wallcycle* wcycle,
297 gmx_enerdata_t* enerd,
302 gmx::ArrayRef<real> constraintsRmsdData,
303 gmx::SimulationSignaller* signalCoordinator,
304 const matrix lastbox,
305 gmx_bool* bSumEkinhOld,
308 gmx::ObservablesReducer* observablesReducer)
310 gmx_bool bEner, bPres, bTemp;
311 gmx_bool bStopCM, bGStat, bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
312 gmx_bool bCheckNumberOfBondedInteractions;
315 /* translate CGLO flags to gmx_booleans */
316 bStopCM = ((flags & CGLO_STOPCM) != 0);
317 bGStat = ((flags & CGLO_GSTAT) != 0);
318 bReadEkin = ((flags & CGLO_READEKIN) != 0);
319 bScaleEkin = ((flags & CGLO_SCALEEKIN) != 0);
320 bEner = ((flags & CGLO_ENERGY) != 0);
321 bTemp = ((flags & CGLO_TEMPERATURE) != 0);
322 bPres = ((flags & CGLO_PRESSURE) != 0);
323 bConstrain = ((flags & CGLO_CONSTRAINT) != 0);
324 bCheckNumberOfBondedInteractions = ((flags & CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS) != 0);
326 /* we calculate a full state kinetic energy either with full-step velocity verlet
327 or half step where we need the pressure */
329 bEkinAveVel = (ir->eI == IntegrationAlgorithm::VV
330 || (ir->eI == IntegrationAlgorithm::VVAK && bPres) || bReadEkin);
332 /* in initalization, it sums the shake virial in vv, and to
333 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
335 /* ########## Kinetic energy ############## */
341 calc_ke_part(x, v, box, &(ir->opts), mdatoms, ekind, nrnb, bEkinAveVel);
345 /* Calculate center of mass velocity if necessary, also parallellized */
348 calc_vcm_grp(*mdatoms, x, v, vcm);
351 if (bTemp || bStopCM || bPres || bEner || bConstrain || bCheckNumberOfBondedInteractions
352 || !observablesReducer->communicationBuffer().empty())
356 /* We will not sum ekinh_old,
357 * so signal that we still have to do it.
359 *bSumEkinhOld = TRUE;
363 gmx::ArrayRef<real> signalBuffer = signalCoordinator->getCommunicationBuffer();
366 wallcycle_start(wcycle, WallCycleCounter::MoveE);
375 bStopCM ? vcm : nullptr,
381 wallcycle_stop(wcycle, WallCycleCounter::MoveE);
383 signalCoordinator->finalizeSignals();
384 *bSumEkinhOld = FALSE;
390 /* Calculate the amplitude of the cosine velocity profile */
391 ekind->cosacc.vcos = ekind->cosacc.mvcos / mdatoms->tmass;
396 /* Sum the kinetic energies of the groups & calc temp */
397 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with inputrecNptTrotter */
398 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
399 Leap with AveVel is not supported; it's not clear that it will actually work.
400 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
401 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
403 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin, bEkinAveVel, bScaleEkin);
404 enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Mass] = static_cast<double>(dvdl_ekin);
406 enerd->term[F_EKIN] = trace(ekind->ekin);
409 /* ########## Now pressure ############## */
410 // TODO: For the VV integrator bConstrain is needed in the conditional. This is confusing, so get rid of this.
411 if (bPres || bConstrain)
413 m_add(force_vir, shake_vir, total_vir);
415 /* Calculate pressure and apply LR correction if PPPM is used.
416 * Use the box from last timestep since we already called update().
419 enerd->term[F_PRES] = calc_pres(fr->pbcType, ir->nwall, lastbox, ekind->ekin, total_vir, pres);
423 static void min_zero(int* n, int i)
425 if (i > 0 && (*n == 0 || i < *n))
431 static int lcd4(int i1, int i2, int i3, int i4)
442 gmx_incons("All 4 inputs for determining nstglobalcomm are <= 0");
446 && ((i1 > 0 && i1 % nst != 0) || (i2 > 0 && i2 % nst != 0) || (i3 > 0 && i3 % nst != 0)
447 || (i4 > 0 && i4 % nst != 0)))
455 int computeGlobalCommunicationPeriod(const t_inputrec* ir)
457 int nstglobalcomm = 10;
459 // Set up the default behaviour
460 if (!(ir->nstcalcenergy > 0 || ir->nstlist > 0 || ir->etc != TemperatureCoupling::No
461 || ir->epc != PressureCoupling::No))
463 /* The user didn't choose the period for anything
464 important, so we just make sure we can send signals and
465 write output suitably. */
466 if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
468 nstglobalcomm = ir->nstenergy;
473 /* The user has made a choice (perhaps implicitly), so we
474 * ensure that we do timely intra-simulation communication
475 * for (possibly) each of the four parts that care.
477 * TODO Does the Verlet scheme (+ DD) need any
478 * communication at nstlist steps? Is the use of nstlist
479 * here a leftover of the twin-range scheme? Can we remove
480 * nstlist when we remove the group scheme?
482 nstglobalcomm = lcd4(ir->nstcalcenergy,
484 ir->etc != TemperatureCoupling::No ? ir->nsttcouple : 0,
485 ir->epc != PressureCoupling::No ? ir->nstpcouple : 0);
488 return nstglobalcomm;
491 int computeGlobalCommunicationPeriod(const gmx::MDLogger& mdlog, const t_inputrec* ir, const t_commrec* cr)
493 const int nstglobalcomm = computeGlobalCommunicationPeriod(ir);
498 .appendTextFormatted("Intra-simulation communication will occur every %d steps.\n",
501 return nstglobalcomm;
504 void rerun_parallel_comm(t_commrec* cr, t_trxframe* fr, gmx_bool* bLastStep)
508 if (MASTER(cr) && *bLastStep)
514 gmx_bcast(sizeof(*fr), fr, cr->mpi_comm_mygroup);
518 *bLastStep = (fr->natoms < 0);
521 // TODO Most of this logic seems to belong in the respective modules
522 void set_state_entries(t_state* state, const t_inputrec* ir, bool useModularSimulator)
524 /* The entries in the state in the tpx file might not correspond
525 * with what is needed, so we correct this here.
528 if (ir->efep != FreeEnergyPerturbationType::No || ir->bExpanded)
530 state->flags |= enumValueToBitMask(StateEntry::Lambda);
531 state->flags |= enumValueToBitMask(StateEntry::FepState);
533 state->flags |= enumValueToBitMask(StateEntry::X);
534 GMX_RELEASE_ASSERT(state->x.size() == state->natoms,
535 "We should start a run with an initialized state->x");
536 if (EI_DYNAMICS(ir->eI))
538 state->flags |= enumValueToBitMask(StateEntry::V);
542 if (ir->pbcType != PbcType::No)
544 state->flags |= enumValueToBitMask(StateEntry::Box);
545 if (inputrecPreserveShape(ir))
547 state->flags |= enumValueToBitMask(StateEntry::BoxRel);
549 if ((ir->epc == PressureCoupling::ParrinelloRahman) || (ir->epc == PressureCoupling::Mttk))
551 state->flags |= enumValueToBitMask(StateEntry::BoxV);
552 if (!useModularSimulator)
554 state->flags |= enumValueToBitMask(StateEntry::PressurePrevious);
557 if (inputrecNptTrotter(ir) || (inputrecNphTrotter(ir)))
560 state->flags |= enumValueToBitMask(StateEntry::Nhpresxi);
561 state->flags |= enumValueToBitMask(StateEntry::Nhpresvxi);
562 state->flags |= enumValueToBitMask(StateEntry::SVirPrev);
563 state->flags |= enumValueToBitMask(StateEntry::FVirPrev);
564 state->flags |= enumValueToBitMask(StateEntry::Veta);
565 state->flags |= enumValueToBitMask(StateEntry::Vol0);
567 if (ir->epc == PressureCoupling::Berendsen || ir->epc == PressureCoupling::CRescale)
569 state->flags |= enumValueToBitMask(StateEntry::BarosInt);
573 if (ir->etc == TemperatureCoupling::NoseHoover)
575 state->flags |= enumValueToBitMask(StateEntry::Nhxi);
576 state->flags |= enumValueToBitMask(StateEntry::Nhvxi);
579 if (ir->etc == TemperatureCoupling::VRescale || ir->etc == TemperatureCoupling::Berendsen)
581 state->flags |= enumValueToBitMask(StateEntry::ThermInt);
584 init_gtc_state(state, state->ngtc, state->nnhpres, ir->opts.nhchainlength); /* allocate the space for nose-hoover chains */
585 init_ekinstate(&state->ekinstate, ir);
587 if (ir->bExpanded && !useModularSimulator)
589 snew(state->dfhist, 1);
590 init_df_history(state->dfhist, ir->fepvals->n_lambda);
593 if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
595 state->flags |= enumValueToBitMask(StateEntry::PullComPrevStep);