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, 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/mdrun.h"
52 #include "gromacs/mdlib/sim_util.h"
53 #include "gromacs/mdlib/simulationsignal.h"
54 #include "gromacs/mdlib/tgroup.h"
55 #include "gromacs/mdlib/update.h"
56 #include "gromacs/mdlib/vcm.h"
57 #include "gromacs/mdtypes/commrec.h"
58 #include "gromacs/mdtypes/df_history.h"
59 #include "gromacs/mdtypes/energyhistory.h"
60 #include "gromacs/mdtypes/forcerec.h"
61 #include "gromacs/mdtypes/group.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/mdtypes/state.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/pulling/pull.h"
67 #include "gromacs/timing/wallcycle.h"
68 #include "gromacs/topology/mtop_util.h"
69 #include "gromacs/trajectory/trajectoryframe.h"
70 #include "gromacs/utility/arrayref.h"
71 #include "gromacs/utility/cstringutil.h"
72 #include "gromacs/utility/fatalerror.h"
73 #include "gromacs/utility/gmxassert.h"
74 #include "gromacs/utility/logger.h"
75 #include "gromacs/utility/smalloc.h"
76 #include "gromacs/utility/snprintf.h"
78 // TODO move this to multi-sim module
79 bool multisim_int_all_are_equal(const gmx_multisim_t
*ms
,
82 bool allValuesAreEqual
= true;
85 GMX_RELEASE_ASSERT(ms
, "Invalid use of multi-simulation pointer");
88 /* send our value to all other master ranks, receive all of theirs */
90 gmx_sumli_sim(ms
->nsim
, buf
, ms
);
92 for (int s
= 0; s
< ms
->nsim
; s
++)
96 allValuesAreEqual
= false;
103 return allValuesAreEqual
;
106 int multisim_min(const gmx_multisim_t
*ms
, int nmin
, int n
)
109 gmx_bool bPos
, bEqual
;
114 gmx_sumi_sim(ms
->nsim
, buf
, ms
);
117 for (s
= 0; s
< ms
->nsim
; s
++)
119 bPos
= bPos
&& (buf
[s
] > 0);
120 bEqual
= bEqual
&& (buf
[s
] == buf
[0]);
126 nmin
= std::min(nmin
, buf
[0]);
130 /* Find the least common multiple */
131 for (d
= 2; d
< nmin
; d
++)
134 while (s
< ms
->nsim
&& d
% buf
[s
] == 0)
140 /* We found the LCM and it is less than nmin */
152 /* TODO Specialize this routine into init-time and loop-time versions?
153 e.g. bReadEkin is only true when restoring from checkpoint */
154 void compute_globals(FILE *fplog
, gmx_global_stat
*gstat
, t_commrec
*cr
, t_inputrec
*ir
,
155 t_forcerec
*fr
, gmx_ekindata_t
*ekind
,
156 t_state
*state
, t_mdatoms
*mdatoms
,
157 t_nrnb
*nrnb
, t_vcm
*vcm
, gmx_wallcycle_t wcycle
,
158 gmx_enerdata_t
*enerd
, tensor force_vir
, tensor shake_vir
, tensor total_vir
,
159 tensor pres
, rvec mu_tot
, gmx::Constraints
*constr
,
160 gmx::SimulationSignaller
*signalCoordinator
,
161 matrix box
, int *totalNumberOfBondedInteractions
,
162 gmx_bool
*bSumEkinhOld
, int flags
)
164 tensor corr_vir
, corr_pres
;
165 gmx_bool bEner
, bPres
, bTemp
;
166 gmx_bool bStopCM
, bGStat
,
167 bReadEkin
, bEkinAveVel
, bScaleEkin
, bConstrain
;
168 gmx_bool bCheckNumberOfBondedInteractions
;
169 real prescorr
, enercorr
, dvdlcorr
, dvdl_ekin
;
171 /* translate CGLO flags to gmx_booleans */
172 bStopCM
= ((flags
& CGLO_STOPCM
) != 0);
173 bGStat
= ((flags
& CGLO_GSTAT
) != 0);
174 bReadEkin
= ((flags
& CGLO_READEKIN
) != 0);
175 bScaleEkin
= ((flags
& CGLO_SCALEEKIN
) != 0);
176 bEner
= ((flags
& CGLO_ENERGY
) != 0);
177 bTemp
= ((flags
& CGLO_TEMPERATURE
) != 0);
178 bPres
= ((flags
& CGLO_PRESSURE
) != 0);
179 bConstrain
= ((flags
& CGLO_CONSTRAINT
) != 0);
180 bCheckNumberOfBondedInteractions
= ((flags
& CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
) != 0);
182 /* we calculate a full state kinetic energy either with full-step velocity verlet
183 or half step where we need the pressure */
185 bEkinAveVel
= (ir
->eI
== eiVV
|| (ir
->eI
== eiVVAK
&& bPres
) || bReadEkin
);
187 /* in initalization, it sums the shake virial in vv, and to
188 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
190 /* ########## Kinetic energy ############## */
194 /* Non-equilibrium MD: this is parallellized, but only does communication
195 * when there really is NEMD.
198 if (PAR(cr
) && (ekind
->bNEMD
))
200 accumulate_u(cr
, &(ir
->opts
), ekind
);
204 calc_ke_part(state
, &(ir
->opts
), mdatoms
, ekind
, nrnb
, bEkinAveVel
);
208 /* Calculate center of mass velocity if necessary, also parallellized */
211 calc_vcm_grp(0, mdatoms
->homenr
, mdatoms
,
212 state
->x
.rvec_array(), state
->v
.rvec_array(), vcm
);
215 if (bTemp
|| bStopCM
|| bPres
|| bEner
|| bConstrain
|| bCheckNumberOfBondedInteractions
)
219 /* We will not sum ekinh_old,
220 * so signal that we still have to do it.
222 *bSumEkinhOld
= TRUE
;
227 gmx::ArrayRef
<real
> signalBuffer
= signalCoordinator
->getCommunicationBuffer();
230 wallcycle_start(wcycle
, ewcMoveE
);
231 global_stat(gstat
, cr
, enerd
, force_vir
, shake_vir
, mu_tot
,
232 ir
, ekind
, constr
, bStopCM
? vcm
: nullptr,
233 signalBuffer
.size(), signalBuffer
.data(),
234 totalNumberOfBondedInteractions
,
235 *bSumEkinhOld
, flags
);
236 wallcycle_stop(wcycle
, ewcMoveE
);
238 signalCoordinator
->finalizeSignals();
239 *bSumEkinhOld
= FALSE
;
243 /* Do center of mass motion removal */
246 check_cm_grp(fplog
, vcm
, ir
, 1);
247 /* At initialization, do not pass x with acceleration-correction mode
248 * to avoid (incorrect) correction of the initial coordinates.
250 rvec
*xPtr
= nullptr;
251 if (vcm
->mode
== ecmANGULAR
|| (vcm
->mode
== ecmLINEAR_ACCELERATION_CORRECTION
&& !(flags
& CGLO_INITIALIZATION
)))
253 xPtr
= state
->x
.rvec_array();
255 do_stopcm_grp(*mdatoms
,
256 xPtr
, state
->v
.rvec_array(), *vcm
);
257 inc_nrnb(nrnb
, eNR_STOPCM
, mdatoms
->homenr
);
262 /* Calculate the amplitude of the cosine velocity profile */
263 ekind
->cosacc
.vcos
= ekind
->cosacc
.mvcos
/mdatoms
->tmass
;
268 /* Sum the kinetic energies of the groups & calc temp */
269 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with inputrecNptTrotter */
270 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
271 Leap with AveVel is not supported; it's not clear that it will actually work.
272 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
273 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
275 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
), ekind
, &dvdl_ekin
,
276 bEkinAveVel
, bScaleEkin
);
277 enerd
->dvdl_lin
[efptMASS
] = static_cast<double>(dvdl_ekin
);
279 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
282 /* ########## Long range energy information ###### */
284 if (bEner
|| bPres
|| bConstrain
)
286 calc_dispcorr(ir
, fr
, box
, state
->lambda
[efptVDW
],
287 corr_pres
, corr_vir
, &prescorr
, &enercorr
, &dvdlcorr
);
292 enerd
->term
[F_DISPCORR
] = enercorr
;
293 enerd
->term
[F_EPOT
] += enercorr
;
294 enerd
->term
[F_DVDL_VDW
] += dvdlcorr
;
297 /* ########## Now pressure ############## */
298 if (bPres
|| bConstrain
)
301 m_add(force_vir
, shake_vir
, total_vir
);
303 /* Calculate pressure and apply LR correction if PPPM is used.
304 * Use the box from last timestep since we already called update().
307 enerd
->term
[F_PRES
] = calc_pres(fr
->ePBC
, ir
->nwall
, box
, ekind
->ekin
, total_vir
, pres
);
309 /* Calculate long range corrections to pressure and energy */
310 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
311 and computes enerd->term[F_DISPCORR]. Also modifies the
312 total_vir and pres tesors */
314 m_add(total_vir
, corr_vir
, total_vir
);
315 m_add(pres
, corr_pres
, pres
);
316 enerd
->term
[F_PDISPCORR
] = prescorr
;
317 enerd
->term
[F_PRES
] += prescorr
;
321 /* check whether an 'nst'-style parameter p is a multiple of nst, and
322 set it to be one if not, with a warning. */
323 static void check_nst_param(const gmx::MDLogger
&mdlog
,
324 const char *desc_nst
, int nst
,
325 const char *desc_p
, int *p
)
327 if (*p
> 0 && *p
% nst
!= 0)
329 /* Round up to the next multiple of nst */
330 *p
= ((*p
)/nst
+ 1)*nst
;
331 GMX_LOG(mdlog
.warning
).asParagraph().appendTextFormatted(
332 "NOTE: %s changes %s to %d", desc_nst
, desc_p
, *p
);
336 void setCurrentLambdasRerun(int64_t step
, const t_lambda
*fepvals
,
337 const t_trxframe
*rerun_fr
, const double *lam0
,
338 t_state
*globalState
)
340 GMX_RELEASE_ASSERT(globalState
!= nullptr, "setCurrentLambdasGlobalRerun should be called with a valid state object");
342 if (rerun_fr
->bLambda
)
344 if (fepvals
->delta_lambda
== 0)
346 globalState
->lambda
[efptFEP
] = rerun_fr
->lambda
;
350 /* find out between which two value of lambda we should be */
351 real frac
= step
*fepvals
->delta_lambda
;
352 int fep_state
= static_cast<int>(std::floor(frac
*fepvals
->n_lambda
));
353 /* interpolate between this state and the next */
354 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
355 frac
= frac
*fepvals
->n_lambda
- fep_state
;
356 for (int i
= 0; i
< efptNR
; i
++)
358 globalState
->lambda
[i
] = lam0
[i
] + (fepvals
->all_lambda
[i
][fep_state
]) +
359 frac
*(fepvals
->all_lambda
[i
][fep_state
+1] - fepvals
->all_lambda
[i
][fep_state
]);
363 else if (rerun_fr
->bFepState
)
365 globalState
->fep_state
= rerun_fr
->fep_state
;
366 for (int i
= 0; i
< efptNR
; i
++)
368 globalState
->lambda
[i
] = fepvals
->all_lambda
[i
][globalState
->fep_state
];
373 void setCurrentLambdasLocal(int64_t step
, const t_lambda
*fepvals
,
374 const double *lam0
, t_state
*state
)
375 /* find the current lambdas. If rerunning, we either read in a state, or a lambda value,
376 requiring different logic. */
378 if (fepvals
->delta_lambda
!= 0)
380 /* find out between which two value of lambda we should be */
381 real frac
= step
*fepvals
->delta_lambda
;
382 if (fepvals
->n_lambda
> 0)
384 int fep_state
= static_cast<int>(std::floor(frac
*fepvals
->n_lambda
));
385 /* interpolate between this state and the next */
386 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
387 frac
= frac
*fepvals
->n_lambda
- fep_state
;
388 for (int i
= 0; i
< efptNR
; i
++)
390 state
->lambda
[i
] = lam0
[i
] + (fepvals
->all_lambda
[i
][fep_state
]) +
391 frac
*(fepvals
->all_lambda
[i
][fep_state
+ 1] - fepvals
->all_lambda
[i
][fep_state
]);
396 for (int i
= 0; i
< efptNR
; i
++)
398 state
->lambda
[i
] = lam0
[i
] + frac
;
404 /* if < 0, fep_state was never defined, and we should not set lambda from the state */
405 if (state
->fep_state
> -1)
407 for (int i
= 0; i
< efptNR
; i
++)
409 state
->lambda
[i
] = fepvals
->all_lambda
[i
][state
->fep_state
];
415 static void min_zero(int *n
, int i
)
417 if (i
> 0 && (*n
== 0 || i
< *n
))
423 static int lcd4(int i1
, int i2
, int i3
, int i4
)
434 gmx_incons("All 4 inputs for determining nstglobalcomm are <= 0");
437 while (nst
> 1 && ((i1
> 0 && i1
% nst
!= 0) ||
438 (i2
> 0 && i2
% nst
!= 0) ||
439 (i3
> 0 && i3
% nst
!= 0) ||
440 (i4
> 0 && i4
% nst
!= 0)))
448 int check_nstglobalcomm(const gmx::MDLogger
&mdlog
, int nstglobalcomm
, t_inputrec
*ir
)
450 if (!EI_DYNAMICS(ir
->eI
))
455 if (nstglobalcomm
== -1)
457 // Set up the default behaviour
458 if (!(ir
->nstcalcenergy
> 0 ||
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. */
467 if (ir
->nstenergy
> 0 && ir
->nstenergy
< nstglobalcomm
)
469 nstglobalcomm
= ir
->nstenergy
;
474 /* The user has made a choice (perhaps implicitly), so we
475 * ensure that we do timely intra-simulation communication
476 * for (possibly) each of the four parts that care.
478 * TODO Does the Verlet scheme (+ DD) need any
479 * communication at nstlist steps? Is the use of nstlist
480 * here a leftover of the twin-range scheme? Can we remove
481 * nstlist when we remove the group scheme?
483 nstglobalcomm
= lcd4(ir
->nstcalcenergy
,
485 ir
->etc
!= etcNO
? ir
->nsttcouple
: 0,
486 ir
->epc
!= epcNO
? ir
->nstpcouple
: 0);
491 // Check that the user's choice of mdrun -gcom will work
492 if (ir
->nstlist
> 0 &&
493 nstglobalcomm
> ir
->nstlist
&& nstglobalcomm
% ir
->nstlist
!= 0)
495 nstglobalcomm
= (nstglobalcomm
/ ir
->nstlist
)*ir
->nstlist
;
496 GMX_LOG(mdlog
.warning
).asParagraph().appendTextFormatted(
497 "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d",
500 if (ir
->nstcalcenergy
> 0)
502 check_nst_param(mdlog
, "-gcom", nstglobalcomm
,
503 "nstcalcenergy", &ir
->nstcalcenergy
);
505 if (ir
->etc
!= etcNO
&& ir
->nsttcouple
> 0)
507 check_nst_param(mdlog
, "-gcom", nstglobalcomm
,
508 "nsttcouple", &ir
->nsttcouple
);
510 if (ir
->epc
!= epcNO
&& ir
->nstpcouple
> 0)
512 check_nst_param(mdlog
, "-gcom", nstglobalcomm
,
513 "nstpcouple", &ir
->nstpcouple
);
516 check_nst_param(mdlog
, "-gcom", nstglobalcomm
,
517 "nstenergy", &ir
->nstenergy
);
519 check_nst_param(mdlog
, "-gcom", nstglobalcomm
,
520 "nstlog", &ir
->nstlog
);
523 if (ir
->comm_mode
!= ecmNO
&& ir
->nstcomm
< nstglobalcomm
)
525 GMX_LOG(mdlog
.warning
).asParagraph().appendTextFormatted(
526 "WARNING: Changing nstcomm from %d to %d",
527 ir
->nstcomm
, nstglobalcomm
);
528 ir
->nstcomm
= nstglobalcomm
;
531 GMX_LOG(mdlog
.info
).appendTextFormatted(
532 "Intra-simulation communication will occur every %d steps.\n", nstglobalcomm
);
533 return nstglobalcomm
;
536 void rerun_parallel_comm(t_commrec
*cr
, t_trxframe
*fr
,
541 if (MASTER(cr
) && *bLastStep
)
547 gmx_bcast(sizeof(*fr
), fr
, cr
);
551 *bLastStep
= (fr
->natoms
< 0);
555 // TODO Most of this logic seems to belong in the respective modules
556 void set_state_entries(t_state
*state
, const t_inputrec
*ir
)
558 /* The entries in the state in the tpx file might not correspond
559 * with what is needed, so we correct this here.
562 if (ir
->efep
!= efepNO
|| ir
->bExpanded
)
564 state
->flags
|= (1<<estLAMBDA
);
565 state
->flags
|= (1<<estFEPSTATE
);
567 state
->flags
|= (1<<estX
);
568 GMX_RELEASE_ASSERT(state
->x
.size() == state
->natoms
, "We should start a run with an initialized state->x");
569 if (EI_DYNAMICS(ir
->eI
))
571 state
->flags
|= (1<<estV
);
575 if (ir
->ePBC
!= epbcNONE
)
577 state
->flags
|= (1<<estBOX
);
578 if (inputrecPreserveShape(ir
))
580 state
->flags
|= (1<<estBOX_REL
);
582 if ((ir
->epc
== epcPARRINELLORAHMAN
) || (ir
->epc
== epcMTTK
))
584 state
->flags
|= (1<<estBOXV
);
585 state
->flags
|= (1<<estPRES_PREV
);
587 if (inputrecNptTrotter(ir
) || (inputrecNphTrotter(ir
)))
590 state
->flags
|= (1<<estNHPRES_XI
);
591 state
->flags
|= (1<<estNHPRES_VXI
);
592 state
->flags
|= (1<<estSVIR_PREV
);
593 state
->flags
|= (1<<estFVIR_PREV
);
594 state
->flags
|= (1<<estVETA
);
595 state
->flags
|= (1<<estVOL0
);
597 if (ir
->epc
== epcBERENDSEN
)
599 state
->flags
|= (1<<estBAROS_INT
);
603 if (ir
->etc
== etcNOSEHOOVER
)
605 state
->flags
|= (1<<estNH_XI
);
606 state
->flags
|= (1<<estNH_VXI
);
609 if (ir
->etc
== etcVRESCALE
|| ir
->etc
== etcBERENDSEN
)
611 state
->flags
|= (1<<estTHERM_INT
);
614 init_gtc_state(state
, state
->ngtc
, state
->nnhpres
, ir
->opts
.nhchainlength
); /* allocate the space for nose-hoover chains */
615 init_ekinstate(&state
->ekinstate
, ir
);
619 snew(state
->dfhist
, 1);
620 init_df_history(state
->dfhist
, ir
->fepvals
->n_lambda
);
623 if (ir
->pull
&& ir
->pull
->bSetPbcRefToPrevStepCOM
)
625 state
->flags
|= (1<<estPULLCOMPREVSTEP
);