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 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.
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/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"
82 // TODO move this to multi-sim module
83 bool multisim_int_all_are_equal(const gmx_multisim_t
* ms
, int64_t value
)
85 bool allValuesAreEqual
= true;
88 GMX_RELEASE_ASSERT(ms
, "Invalid use of multi-simulation pointer");
91 /* send our value to all other master ranks, receive all of theirs */
93 gmx_sumli_sim(ms
->nsim
, buf
, ms
);
95 for (int s
= 0; s
< ms
->nsim
; s
++)
99 allValuesAreEqual
= false;
106 return allValuesAreEqual
;
109 int multisim_min(const gmx_multisim_t
* ms
, int nmin
, int n
)
112 gmx_bool bPos
, bEqual
;
117 gmx_sumi_sim(ms
->nsim
, buf
, ms
);
120 for (s
= 0; s
< ms
->nsim
; s
++)
122 bPos
= bPos
&& (buf
[s
] > 0);
123 bEqual
= bEqual
&& (buf
[s
] == buf
[0]);
129 nmin
= std::min(nmin
, buf
[0]);
133 /* Find the least common multiple */
134 for (d
= 2; d
< nmin
; d
++)
137 while (s
< ms
->nsim
&& d
% buf
[s
] == 0)
143 /* We found the LCM and it is less than nmin */
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
,
159 const t_inputrec
* ir
,
161 gmx_ekindata_t
* ekind
,
162 gmx::ArrayRef
<const gmx::RVec
> x
,
163 gmx::ArrayRef
<const gmx::RVec
> v
,
166 const t_mdatoms
* mdatoms
,
169 gmx_wallcycle_t wcycle
,
170 gmx_enerdata_t
* enerd
,
175 gmx::Constraints
* constr
,
176 gmx::SimulationSignaller
* signalCoordinator
,
177 const matrix lastbox
,
178 int* totalNumberOfBondedInteractions
,
179 gmx_bool
* bSumEkinhOld
,
182 gmx_bool bEner
, bPres
, bTemp
;
183 gmx_bool bStopCM
, bGStat
, bReadEkin
, bEkinAveVel
, bScaleEkin
, bConstrain
;
184 gmx_bool bCheckNumberOfBondedInteractions
;
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);
198 /* we calculate a full state kinetic energy either with full-step velocity verlet
199 or half step where we need the pressure */
201 bEkinAveVel
= (ir
->eI
== eiVV
|| (ir
->eI
== eiVVAK
&& bPres
) || bReadEkin
);
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 */
206 /* ########## Kinetic energy ############## */
210 /* Non-equilibrium MD: this is parallellized, but only does communication
211 * when there really is NEMD.
214 if (PAR(cr
) && (ekind
->bNEMD
))
216 accumulate_u(cr
, &(ir
->opts
), ekind
);
220 calc_ke_part(x
, v
, box
, &(ir
->opts
), mdatoms
, ekind
, nrnb
, bEkinAveVel
);
224 /* Calculate center of mass velocity if necessary, also parallellized */
227 calc_vcm_grp(*mdatoms
, x
, v
, vcm
);
230 if (bTemp
|| bStopCM
|| bPres
|| bEner
|| bConstrain
|| bCheckNumberOfBondedInteractions
)
234 /* We will not sum ekinh_old,
235 * so signal that we still have to do it.
237 *bSumEkinhOld
= TRUE
;
241 gmx::ArrayRef
<real
> signalBuffer
= signalCoordinator
->getCommunicationBuffer();
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
);
250 signalCoordinator
->finalizeSignals();
251 *bSumEkinhOld
= FALSE
;
257 /* Calculate the amplitude of the cosine velocity profile */
258 ekind
->cosacc
.vcos
= ekind
->cosacc
.mvcos
/ mdatoms
->tmass
;
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.
270 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
), ekind
, &dvdl_ekin
, bEkinAveVel
, bScaleEkin
);
271 enerd
->dvdl_lin
[efptMASS
] = static_cast<double>(dvdl_ekin
);
273 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
275 for (auto& dhdl
: enerd
->dhdlLambda
)
277 dhdl
+= enerd
->dvdl_lin
[efptMASS
];
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
)
285 m_add(force_vir
, shake_vir
, total_vir
);
287 /* Calculate pressure and apply LR correction if PPPM is used.
288 * Use the box from last timestep since we already called update().
291 enerd
->term
[F_PRES
] = calc_pres(fr
->pbcType
, ir
->nwall
, lastbox
, ekind
->ekin
, total_vir
, pres
);
294 /* ########## Long range energy information ###### */
295 if ((bEner
|| bPres
) && fr
->dispersionCorrection
)
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 */
302 const DispersionCorrection::Correction correction
=
303 fr
->dispersionCorrection
->calculate(lastbox
, vdwLambda
);
307 enerd
->term
[F_DISPCORR
] = correction
.energy
;
308 enerd
->term
[F_EPOT
] += correction
.energy
;
309 enerd
->term
[F_DVDL_VDW
] += correction
.dvdl
;
313 correction
.correctVirial(total_vir
);
314 correction
.correctPressure(pres
);
315 enerd
->term
[F_PDISPCORR
] = correction
.pressure
;
316 enerd
->term
[F_PRES
] += correction
.pressure
;
321 void setCurrentLambdasRerun(int64_t step
,
322 const t_lambda
* fepvals
,
323 const t_trxframe
* rerun_fr
,
325 t_state
* globalState
)
327 GMX_RELEASE_ASSERT(globalState
!= nullptr,
328 "setCurrentLambdasGlobalRerun should be called with a valid state object");
330 if (rerun_fr
->bLambda
)
332 if (fepvals
->delta_lambda
== 0)
334 globalState
->lambda
[efptFEP
] = rerun_fr
->lambda
;
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
++)
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
]);
352 else if (rerun_fr
->bFepState
)
354 globalState
->fep_state
= rerun_fr
->fep_state
;
355 for (int i
= 0; i
< efptNR
; i
++)
357 globalState
->lambda
[i
] = fepvals
->all_lambda
[i
][globalState
->fep_state
];
362 void setCurrentLambdasLocal(const int64_t step
,
363 const t_lambda
* fepvals
,
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. */
370 if (fepvals
->delta_lambda
!= 0)
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)
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
++)
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
]);
388 for (int i
= 0; i
< efptNR
; i
++)
390 lambda
[i
] = lam0
[i
] + frac
;
396 /* if < 0, fep_state was never defined, and we should not set lambda from the state */
397 if (currentFEPState
> -1)
399 for (int i
= 0; i
< efptNR
; i
++)
401 lambda
[i
] = fepvals
->all_lambda
[i
][currentFEPState
];
407 static void min_zero(int* n
, int i
)
409 if (i
> 0 && (*n
== 0 || i
< *n
))
415 static int lcd4(int i1
, int i2
, int i3
, int i4
)
426 gmx_incons("All 4 inputs for determining nstglobalcomm are <= 0");
430 && ((i1
> 0 && i1
% nst
!= 0) || (i2
> 0 && i2
% nst
!= 0) || (i3
> 0 && i3
% nst
!= 0)
431 || (i4
> 0 && i4
% nst
!= 0)))
439 int computeGlobalCommunicationPeriod(const gmx::MDLogger
& mdlog
, t_inputrec
* ir
, const t_commrec
* cr
)
443 // Set up the default behaviour
444 if (!(ir
->nstcalcenergy
> 0 || ir
->nstlist
> 0 || ir
->etc
!= etcNO
|| ir
->epc
!= epcNO
))
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. */
450 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
452 nstglobalcomm
= ir
->nstenergy
;
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.
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?
466 nstglobalcomm
= lcd4(ir
->nstcalcenergy
, ir
->nstlist
, ir
->etc
!= etcNO
? ir
->nsttcouple
: 0,
467 ir
->epc
!= epcNO
? ir
->nstpcouple
: 0);
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
)
475 GMX_LOG(mdlog
.warning
)
477 .appendTextFormatted("WARNING: Changing nstcomm from %d to %d", ir
->nstcomm
, nstglobalcomm
);
478 ir
->nstcomm
= nstglobalcomm
;
484 .appendTextFormatted("Intra-simulation communication will occur every %d steps.\n",
487 return nstglobalcomm
;
490 void rerun_parallel_comm(t_commrec
* cr
, t_trxframe
* fr
, gmx_bool
* bLastStep
)
494 if (MASTER(cr
) && *bLastStep
)
500 gmx_bcast(sizeof(*fr
), fr
, cr
);
504 *bLastStep
= (fr
->natoms
< 0);
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
)
510 /* The entries in the state in the tpx file might not correspond
511 * with what is needed, so we correct this here.
514 if (ir
->efep
!= efepNO
|| ir
->bExpanded
)
516 state
->flags
|= (1 << estLAMBDA
);
517 state
->flags
|= (1 << estFEPSTATE
);
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
))
524 state
->flags
|= (1 << estV
);
528 if (ir
->pbcType
!= PbcType::No
)
530 state
->flags
|= (1 << estBOX
);
531 if (inputrecPreserveShape(ir
))
533 state
->flags
|= (1 << estBOX_REL
);
535 if ((ir
->epc
== epcPARRINELLORAHMAN
) || (ir
->epc
== epcMTTK
))
537 state
->flags
|= (1 << estBOXV
);
538 state
->flags
|= (1 << estPRES_PREV
);
540 if (inputrecNptTrotter(ir
) || (inputrecNphTrotter(ir
)))
543 state
->flags
|= (1 << estNHPRES_XI
);
544 state
->flags
|= (1 << estNHPRES_VXI
);
545 state
->flags
|= (1 << estSVIR_PREV
);
546 state
->flags
|= (1 << estFVIR_PREV
);
547 state
->flags
|= (1 << estVETA
);
548 state
->flags
|= (1 << estVOL0
);
550 if (ir
->epc
== epcBERENDSEN
)
552 state
->flags
|= (1 << estBAROS_INT
);
556 if (ir
->etc
== etcNOSEHOOVER
)
558 state
->flags
|= (1 << estNH_XI
);
559 state
->flags
|= (1 << estNH_VXI
);
562 if (ir
->etc
== etcVRESCALE
|| ir
->etc
== etcBERENDSEN
)
564 state
->flags
|= (1 << estTHERM_INT
);
567 init_gtc_state(state
, state
->ngtc
, state
->nnhpres
,
568 ir
->opts
.nhchainlength
); /* allocate the space for nose-hoover chains */
569 init_ekinstate(&state
->ekinstate
, ir
);
573 snew(state
->dfhist
, 1);
574 init_df_history(state
->dfhist
, ir
->fepvals
->n_lambda
);
577 if (ir
->pull
&& ir
->pull
->bSetPbcRefToPrevStepCOM
)
579 state
->flags
|= (1 << estPULLCOMPREVSTEP
);