*/
#include "gmxpre.h"
+#include "md.h"
+
#include "config.h"
#include <math.h>
#include "thread_mpi/threads.h"
+#include "gromacs/commandline/filenm.h"
#include "gromacs/domdec/domdec.h"
#include "gromacs/domdec/domdec_network.h"
-#include "gromacs/ewald/pme-load-balancing.h"
+#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/ewald/pme.h"
-#include "gromacs/fileio/filenm.h"
-#include "gromacs/fileio/mdoutf.h"
-#include "gromacs/fileio/trajectory_writing.h"
-#include "gromacs/fileio/trx.h"
+#include "gromacs/ewald/pme-load-balancing.h"
#include "gromacs/fileio/trxio.h"
+#include "gromacs/gmxlib/md_logging.h"
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/imd/imd.h"
-#include "gromacs/legacyheaders/constr.h"
-#include "gromacs/legacyheaders/ebin.h"
-#include "gromacs/legacyheaders/force.h"
-#include "gromacs/legacyheaders/md_logging.h"
-#include "gromacs/legacyheaders/md_support.h"
-#include "gromacs/legacyheaders/mdatoms.h"
-#include "gromacs/legacyheaders/mdebin.h"
-#include "gromacs/legacyheaders/mdrun.h"
-#include "gromacs/legacyheaders/network.h"
-#include "gromacs/legacyheaders/nrnb.h"
-#include "gromacs/legacyheaders/ns.h"
-#include "gromacs/legacyheaders/shellfc.h"
-#include "gromacs/legacyheaders/sighandler.h"
-#include "gromacs/legacyheaders/sim_util.h"
-#include "gromacs/legacyheaders/tgroup.h"
-#include "gromacs/legacyheaders/typedefs.h"
-#include "gromacs/legacyheaders/update.h"
-#include "gromacs/legacyheaders/vcm.h"
-#include "gromacs/legacyheaders/vsite.h"
-#include "gromacs/legacyheaders/types/commrec.h"
-#include "gromacs/legacyheaders/types/constr.h"
-#include "gromacs/legacyheaders/types/enums.h"
-#include "gromacs/legacyheaders/types/fcdata.h"
-#include "gromacs/legacyheaders/types/force_flags.h"
-#include "gromacs/legacyheaders/types/forcerec.h"
-#include "gromacs/legacyheaders/types/group.h"
-#include "gromacs/legacyheaders/types/inputrec.h"
-#include "gromacs/legacyheaders/types/interaction_const.h"
-#include "gromacs/legacyheaders/types/mdatom.h"
-#include "gromacs/legacyheaders/types/membedt.h"
-#include "gromacs/legacyheaders/types/nrnb.h"
-#include "gromacs/legacyheaders/types/oenv.h"
-#include "gromacs/legacyheaders/types/shellfc.h"
-#include "gromacs/legacyheaders/types/state.h"
#include "gromacs/listed-forces/manage-threading.h"
+#include "gromacs/math/functions.h"
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/mdlib/compute_io.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/ebin.h"
+#include "gromacs/mdlib/force.h"
+#include "gromacs/mdlib/force_flags.h"
+#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/md_support.h"
+#include "gromacs/mdlib/mdatoms.h"
+#include "gromacs/mdlib/mdebin.h"
+#include "gromacs/mdlib/mdoutf.h"
+#include "gromacs/mdlib/mdrun.h"
#include "gromacs/mdlib/mdrun_signalling.h"
#include "gromacs/mdlib/nb_verlet.h"
#include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
+#include "gromacs/mdlib/ns.h"
+#include "gromacs/mdlib/shellfc.h"
+#include "gromacs/mdlib/sighandler.h"
+#include "gromacs/mdlib/sim_util.h"
+#include "gromacs/mdlib/tgroup.h"
+#include "gromacs/mdlib/trajectory_writing.h"
+#include "gromacs/mdlib/update.h"
+#include "gromacs/mdlib/vcm.h"
+#include "gromacs/mdlib/vsite.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/df_history.h"
+#include "gromacs/mdtypes/energyhistory.h"
+#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/group.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/interaction_const.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/trajectory/trajectoryframe.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "corewrap.h"
#endif
+/*! \brief Check whether bonded interactions are missing, if appropriate
+ *
+ * \param[in] fplog Log file pointer
+ * \param[in] cr Communication object
+ * \param[in] totalNumberOfBondedInteractions Result of the global reduction over the number of bonds treated in each domain
+ * \param[in] top_global Global topology for the error message
+ * \param[in] top_local Local topology for the error message
+ * \param[in] state Global state for the error message
+ * \param[inout] shouldCheckNumberOfBondedInteractions Whether we should do the check.
+ *
+ * \return Nothing, except that shouldCheckNumberOfBondedInteractions
+ * is always set to false after exit.
+ */
+static void checkNumberOfBondedInteractions(FILE *fplog, t_commrec *cr, int totalNumberOfBondedInteractions,
+ gmx_mtop_t *top_global, gmx_localtop_t *top_local, t_state *state,
+ bool *shouldCheckNumberOfBondedInteractions)
+{
+ if (*shouldCheckNumberOfBondedInteractions)
+ {
+ if (totalNumberOfBondedInteractions != cr->dd->nbonded_global)
+ {
+ dd_print_missing_interactions(fplog, cr, totalNumberOfBondedInteractions, top_global, top_local, state); // Does not return
+ }
+ *shouldCheckNumberOfBondedInteractions = false;
+ }
+}
+
static void reset_all_counters(FILE *fplog, t_commrec *cr,
gmx_int64_t step,
gmx_int64_t *step_rel, t_inputrec *ir,
print_date_and_time(fplog, cr->nodeid, "Restarted time", gmx_gettime());
}
-double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
- const output_env_t oenv, gmx_bool bVerbose, gmx_bool bCompact,
- int nstglobalcomm,
- gmx_vsite_t *vsite, gmx_constr_t constr,
- int stepout, t_inputrec *ir,
- gmx_mtop_t *top_global,
- t_fcdata *fcd,
- t_state *state_global,
- t_mdatoms *mdatoms,
- t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- gmx_edsam_t ed, t_forcerec *fr,
- int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t membed,
- real cpt_period, real max_hours,
- int imdport,
- unsigned long Flags,
- gmx_walltime_accounting_t walltime_accounting)
+/*! \libinternal
+ \copydoc integrator_t (FILE *fplog, t_commrec *cr,
+ int nfile, const t_filenm fnm[],
+ const gmx_output_env_t *oenv, gmx_bool bVerbose,
+ int nstglobalcomm,
+ gmx_vsite_t *vsite, gmx_constr_t constr,
+ int stepout,
+ t_inputrec *inputrec,
+ gmx_mtop_t *top_global, t_fcdata *fcd,
+ t_state *state_global,
+ t_mdatoms *mdatoms,
+ t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ gmx_edsam_t ed,
+ t_forcerec *fr,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ gmx_membed_t *membed,
+ real cpt_period, real max_hours,
+ int imdport,
+ unsigned long Flags,
+ gmx_walltime_accounting_t walltime_accounting)
+ */
+double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
+ const gmx_output_env_t *oenv, gmx_bool bVerbose,
+ int nstglobalcomm,
+ gmx_vsite_t *vsite, gmx_constr_t constr,
+ int stepout, t_inputrec *ir,
+ gmx_mtop_t *top_global,
+ t_fcdata *fcd,
+ t_state *state_global,
+ t_mdatoms *mdatoms,
+ t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ gmx_edsam_t ed, t_forcerec *fr,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t *membed,
+ real cpt_period, real max_hours,
+ int imdport,
+ unsigned long Flags,
+ gmx_walltime_accounting_t walltime_accounting)
{
gmx_mdoutf_t outf = NULL;
gmx_int64_t step, step_rel;
double elapsed_time;
double t, t0, lam0[efptNR];
- gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEner;
+ gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
gmx_bool bNS, bNStList, bSimAnn, bStopCM, bRerunMD, bNotLastFrame = FALSE,
- bFirstStep, bStateFromCP, bStateFromTPX, bInitStep, bLastStep,
- bBornRadii, bStartingFromCpt;
+ bFirstStep, startingFromCheckpoint, bInitStep, bLastStep,
+ bBornRadii;
gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
gmx_bool do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE,
bForceUpdate = FALSE, bCPT;
gmx_localtop_t *top;
t_mdebin *mdebin = NULL;
t_state *state = NULL;
- rvec *f_global = NULL;
gmx_enerdata_t *enerd;
rvec *f = NULL;
gmx_global_stat_t gstat;
- gmx_update_t upd = NULL;
+ gmx_update_t *upd = NULL;
t_graph *graph = NULL;
gmx_signalling_t gs;
gmx_groups_t *groups;
gmx_ekindata_t *ekind;
- gmx_shellfc_t shellfc;
- int count, nconverged = 0;
- double tcount = 0;
- gmx_bool bConverged = TRUE, bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
+ gmx_shellfc_t *shellfc;
+ gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
gmx_bool bResetCountersHalfMaxH = FALSE;
- gmx_bool bVV, bTemp, bPres, bTrotter;
- gmx_bool bUpdateDoLR;
+ gmx_bool bTemp, bPres, bTrotter;
real dvdl_constr;
rvec *cbuf = NULL;
int cbuf_nalloc = 0;
/* Temporary addition for FAHCORE checkpointing */
int chkpt_ret;
#endif
+ /* Domain decomposition could incorrectly miss a bonded
+ interaction, but checking for that requires a global
+ communication stage, which does not otherwise happen in DD
+ code. So we do that alongside the first global energy reduction
+ after a new DD is made. These variables handle whether the
+ check happens, and the result it returns. */
+ bool shouldCheckNumberOfBondedInteractions = false;
+ int totalNumberOfBondedInteractions = -1;
/* Check for special mdrun options */
bRerunMD = (Flags & MD_RERUN);
/* md-vv uses averaged full step velocities for T-control
md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
- bVV = EI_VV(ir->eI);
- bTrotter = (bVV && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir)));
+ bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
if (bRerunMD)
{
nstglobalcomm = 1;
}
- check_ir_old_tpx_versions(cr, fplog, ir, top_global);
-
nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
bGStatEveryStep = (nstglobalcomm == 1);
ekind->cosacc.cos_accel = ir->cos_accel;
gstat = global_stat_init(ir);
- debug_gmx();
/* Check for polarizable models and flexible constraints */
shellfc = init_shell_flexcon(fplog,
top_global, n_flexible_constraints(constr),
- (ir->bContinuation ||
- (DOMAINDECOMP(cr) && !MASTER(cr))) ?
- NULL : state_global->x);
+ ir->nstcalcenergy, DOMAINDECOMP(cr));
+
if (shellfc && ir->nstcalcenergy != 1)
{
gmx_fatal(FARGS, "You have nstcalcenergy set to a value (%d) that is different from 1.\nThis is not supported in combinations with shell particles.\nPlease make a new tpr file.", ir->nstcalcenergy);
{
gmx_fatal(FARGS, "Shell particles are not implemented with domain decomposition, use a single rank");
}
- if (shellfc && ir->eI == eiNM)
- {
- /* Currently shells don't work with Normal Modes */
- gmx_fatal(FARGS, "Normal Mode analysis is not supported with shells.\nIf you'd like to help with adding support, we have an open discussion at http://redmine.gromacs.org/issues/879\n");
- }
- if (vsite && ir->eI == eiNM)
- {
- /* Currently virtual sites don't work with Normal Modes */
- gmx_fatal(FARGS, "Normal Mode analysis is not supported with virtual sites.\nIf you'd like to help with adding support, we have an open discussion at http://redmine.gromacs.org/issues/879\n");
- }
-
- if (DEFORM(*ir))
+ if (inputrecDeform(ir))
{
tMPI_Thread_mutex_lock(&deform_init_box_mutex);
set_deform_reference_box(upd,
snew(state, 1);
dd_init_local_state(cr->dd, state_global, state);
-
- if (DDMASTER(cr->dd) && ir->nstfout)
- {
- snew(f_global, state_global->natoms);
- }
}
else
{
- top = gmx_mtop_generate_local_top(top_global, ir);
+ top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
forcerec_set_excl_load(fr, top);
state = serial_init_local_state(state_global);
- f_global = f;
atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
dd_partition_system(fplog, ir->init_step, cr, TRUE, 1,
state_global, top_global, ir,
state, &f, mdatoms, top, fr,
- vsite, shellfc, constr,
+ vsite, constr,
nrnb, NULL, FALSE);
-
+ shouldCheckNumberOfBondedInteractions = true;
}
update_mdatoms(mdatoms, state->lambda[efptMASS]);
- if (opt2bSet("-cpi", nfile, fnm))
- {
- bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi", nfile, fnm, cr), cr);
- }
- else
- {
- bStateFromCP = FALSE;
- }
+ startingFromCheckpoint = Flags & MD_STARTFROMCPT;
if (ir->bExpanded)
{
- init_expanded_ensemble(bStateFromCP, ir, &state->dfhist);
+ init_expanded_ensemble(startingFromCheckpoint, ir, &state->dfhist);
}
if (MASTER(cr))
{
- if (bStateFromCP)
+ if (startingFromCheckpoint)
{
/* Update mdebin with energy history if appending to output files */
if (Flags & MD_APPENDFILES)
{
- restore_energyhistory_from_state(mdebin, &state_global->enerhist);
+ restore_energyhistory_from_state(mdebin, state_global->enerhist);
}
else
{
/* We might have read an energy history from checkpoint,
* free the allocated memory and reset the counts.
*/
- done_energyhistory(&state_global->enerhist);
- init_energyhistory(&state_global->enerhist);
+ done_energyhistory(state_global->enerhist);
+ init_energyhistory(state_global->enerhist);
}
}
/* Set the initial energy history in state by updating once */
- update_energyhistory(&state_global->enerhist, mdebin);
+ update_energyhistory(state_global->enerhist, mdebin);
}
/* Initialize constraints */
}
}
- debug_gmx();
-
- if (IR_TWINRANGE(*ir) && repl_ex_nst % ir->nstcalclr != 0)
- {
- /* We should exchange at nstcalclr steps to get correct integration */
- gmx_fatal(FARGS, "The replica exchange period (%d) is not divisible by nstcalclr (%d)", repl_ex_nst, ir->nstcalclr);
- }
-
if (ir->efep != efepNO)
{
/* Set free energy calculation frequency as the greatest common
- * denominator of nstdhdl and repl_ex_nst.
- * Check for nstcalclr with twin-range, since we need the long-range
- * contribution to the free-energy at the correct (nstcalclr) steps.
- */
+ * denominator of nstdhdl and repl_ex_nst. */
nstfep = ir->fepvals->nstdhdl;
if (ir->bExpanded)
{
- if (IR_TWINRANGE(*ir) &&
- ir->expandedvals->nstexpanded % ir->nstcalclr != 0)
- {
- gmx_fatal(FARGS, "nstexpanded should be divisible by nstcalclr");
- }
nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
}
if (repl_ex_nst > 0)
{
nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
}
- /* We checked divisibility of repl_ex_nst and nstcalclr above */
- if (IR_TWINRANGE(*ir) && nstfep % ir->nstcalclr != 0)
- {
- gmx_incons("nstfep not divisible by nstcalclr");
- }
}
/* Be REALLY careful about what flags you set here. You CANNOT assume
*/
bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
+ if (Flags & MD_READ_EKIN)
+ {
+ restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
+ }
+
cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
| (bStopCM ? CGLO_STOPCM : 0)
- | (bVV ? CGLO_PRESSURE : 0)
- | (bVV ? CGLO_CONSTRAINT : 0)
- | (bRerunMD ? CGLO_RERUNMD : 0)
+ | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
+ | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
| ((Flags & MD_READ_EKIN) ? CGLO_READEKIN : 0));
bSumEkinhOld = FALSE;
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, NULL, FALSE, state->box,
- top_global, &bSumEkinhOld, cglo_flags);
+ &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags
+ | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
+ checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
+ top_global, top, state,
+ &shouldCheckNumberOfBondedInteractions);
if (ir->eI == eiVVAK)
{
/* a second call to get the half step temperature initialized as well */
kinetic energy calculation. This minimized excess variables, but
perhaps loses some logic?*/
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, NULL, FALSE, state->box,
- top_global, &bSumEkinhOld,
+ NULL, &bSumEkinhOld,
cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
}
{
fprintf(fplog,
"RMS relative constraint deviation after constraining: %.2e\n",
- constr_rmsd(constr, FALSE));
+ constr_rmsd(constr));
}
if (EI_STATE_VELOCITY(ir->eI))
{
}
#endif
- debug_gmx();
/***********************************************************
*
* Loop over MD steps
{
gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f does not contain a box, while pbc is used", rerun_fr.step, rerun_fr.time);
}
- if (max_cutoff2(ir->ePBC, rerun_fr.box) < sqr(fr->rlistlong))
+ if (max_cutoff2(ir->ePBC, rerun_fr.box) < gmx::square(fr->rlist))
{
gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr.step, rerun_fr.time);
}
/* loop over MD steps or if rerunMD to end of input trajectory */
bFirstStep = TRUE;
/* Skip the first Nose-Hoover integration when we get the state from tpx */
- bStateFromTPX = !bStateFromCP;
- bInitStep = bFirstStep && (bStateFromTPX || bVV);
- bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
+ bInitStep = !startingFromCheckpoint || EI_VV(ir->eI);
bSumEkinhOld = FALSE;
bExchanged = FALSE;
bNeedRepartition = FALSE;
pme_loadbal_do(pme_loadbal, cr,
(bVerbose && MASTER(cr)) ? stderr : NULL,
fplog,
- ir, fr, state, wcycle,
+ ir, fr, state,
+ wcycle,
step, step_rel,
&bPMETunePrinting);
}
bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
- && (ir->bExpanded) && (step > 0) && (!bStartingFromCpt));
+ && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
}
bDoReplEx = ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
if (bSimAnn)
{
- update_annealing_target_temp(&(ir->opts), t);
+ update_annealing_target_temp(ir, t, upd);
}
if (bRerunMD)
}
else
{
- /* Determine whether or not to do Neighbour Searching and LR */
+ /* Determine whether or not to do Neighbour Searching */
bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
}
bBornRadii = TRUE;
}
- do_log = do_per_step(step, ir->nstlog) || bFirstStep || bLastStep;
+ /* do_log triggers energy and virial calculation. Because this leads
+ * to different code paths, forces can be different. Thus for exact
+ * continuation we should avoid extra log output.
+ * Note that the || bLastStep can result in non-exact continuation
+ * beyond the last step. But we don't consider that to be an issue.
+ */
- do_log = do_per_step(step, ir->nstlog) || (bFirstStep && !bStateFromCP) || bLastStep;
++ do_log = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep;
do_verbose = bVerbose &&
(step % stepout == 0 || bFirstStep || bLastStep);
{
bMasterState = FALSE;
/* Correct the new box if it is too skewed */
- if (DYNAMIC_BOX(*ir))
+ if (inputrecDynamicBox(ir))
{
if (correct_box(fplog, step, state->box, graph))
{
bMasterState, nstglobalcomm,
state_global, top_global, ir,
state, &f, mdatoms, top, fr,
- vsite, shellfc, constr,
+ vsite, constr,
nrnb, wcycle,
do_verbose && !bPMETunePrinting);
+ shouldCheckNumberOfBondedInteractions = true;
}
}
if (MASTER(cr) && do_log)
{
- print_ebin_header(fplog, step, t, state->lambda[efptFEP]); /* can we improve the information printed here? */
+ print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
}
if (ir->efep != efepNO)
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
constr, NULL, FALSE, state->box,
- top_global, &bSumEkinhOld,
- CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
+ &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
+ checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
+ top_global, top, state,
+ &shouldCheckNumberOfBondedInteractions);
}
clear_mat(force_vir);
but the virial needs to be calculated on both the current step and the 'next' step. Future
reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
- bCalcEner = do_per_step(step-1, ir->nstcalcenergy);
- bCalcVir = bCalcEner ||
+ /* TODO: This is probably not what we want, we will write to energy file one step after nstcalcenergy steps. */
+ bCalcEnerStep = do_per_step(step - 1, ir->nstcalcenergy);
+ bCalcVir = bCalcEnerStep ||
(ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
}
else
{
- bCalcEner = do_per_step(step, ir->nstcalcenergy);
- bCalcVir = bCalcEner ||
+ bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
+ bCalcVir = bCalcEnerStep ||
(ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
}
-
- /* Do we need global communication ? */
- bGStat = (bCalcVir || bCalcEner || bStopCM ||
- do_per_step(step, nstglobalcomm) ||
- (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
+ bCalcEner = bCalcEnerStep;
do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
{
bCalcVir = TRUE;
bCalcEner = TRUE;
- bGStat = TRUE;
}
- (bVV && IR_NVT_TROTTER(ir) && do_per_step(step-1, nstglobalcomm)));
-
- /* these CGLO_ options remain the same throughout the iteration */
- cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
- (bGStat ? CGLO_GSTAT : 0)
- );
+ /* Do we need global communication ? */
+ bGStat = (bCalcVir || bCalcEner || bStopCM ||
+ do_per_step(step, nstglobalcomm) ||
++ (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
+
force_flags = (GMX_FORCE_STATECHANGED |
- ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
+ ((inputrecDynamicBox(ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
GMX_FORCE_ALLFORCES |
- GMX_FORCE_SEPLRF |
(bCalcVir ? GMX_FORCE_VIRIAL : 0) |
(bCalcEner ? GMX_FORCE_ENERGY : 0) |
(bDoFEP ? GMX_FORCE_DHDL : 0)
);
- if (fr->bTwinRange)
- {
- if (do_per_step(step, ir->nstcalclr))
- {
- force_flags |= GMX_FORCE_DO_LR;
- }
- }
-
if (shellfc)
{
/* Now is the time to relax the shells */
- count = relax_shell_flexcon(fplog, cr, bVerbose, step,
- ir, bNS, force_flags,
- top,
- constr, enerd, fcd,
- state, f, force_vir, mdatoms,
- nrnb, wcycle, graph, groups,
- shellfc, fr, bBornRadii, t, mu_tot,
- &bConverged, vsite,
- mdoutf_get_fp_field(outf));
- tcount += count;
-
- if (bConverged)
- {
- nconverged++;
- }
+ relax_shell_flexcon(fplog, cr, bVerbose, step,
+ ir, bNS, force_flags, top,
+ constr, enerd, fcd,
+ state, f, force_vir, mdatoms,
+ nrnb, wcycle, graph, groups,
+ shellfc, fr, bBornRadii, t, mu_tot,
+ vsite, mdoutf_get_fp_field(outf));
}
else
{
(bNS ? GMX_FORCE_NS : 0) | force_flags);
}
- if (bVV && !bStartingFromCpt && !bRerunMD)
+ if (EI_VV(ir->eI) && !startingFromCheckpoint && !bRerunMD)
/* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
{
rvec *vbuf = NULL;
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
}
- /* If we are using twin-range interactions where the long-range component
- * is only evaluated every nstcalclr>1 steps, we should do a special update
- * step to combine the long-range forces on these steps.
- * For nstcalclr=1 this is not done, since the forces would have been added
- * directly to the short-range forces already.
- *
- * TODO Remove various aspects of VV+twin-range in master
- * branch, because VV integrators did not ever support
- * twin-range multiple time stepping with constraints.
- */
- bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
-
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC,
- f, bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
- ekind, M, upd, bInitStep, etrtVELOCITY1,
- cr, nrnb, constr, &top->idef);
+ update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ ekind, M, upd, etrtVELOCITY1,
+ cr, constr);
if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
{
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir);
wallcycle_start(wcycle, ewcUPDATE);
- if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
- {
- /* Correct the virial for multiple time stepping */
- m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
- }
}
else if (graph)
{
* Think about ways around this in the future?
* For now, keep this choice in comments.
*/
- /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
- /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
+ /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
+ /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
bPres = TRUE;
bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
if (bCalcEner && ir->eI == eiVVAK)
if (bGStat || do_per_step(step-1, nstglobalcomm))
{
wallcycle_stop(wcycle, ewcUPDATE);
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, NULL, FALSE, state->box,
- top_global, &bSumEkinhOld,
- cglo_flags
+ &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ (bGStat ? CGLO_GSTAT : 0)
| CGLO_ENERGY
| (bTemp ? CGLO_TEMPERATURE : 0)
| (bPres ? CGLO_PRESSURE : 0)
| (bPres ? CGLO_CONSTRAINT : 0)
| (bStopCM ? CGLO_STOPCM : 0)
+ | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
| CGLO_SCALEEKIN
);
/* explanation of above:
time step kinetic energy for the pressure (always true now, since we want accurate statistics).
b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
EkinAveVel because it's needed for the pressure */
+ checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
+ top_global, top, state,
+ &shouldCheckNumberOfBondedInteractions);
wallcycle_start(wcycle, ewcUPDATE);
}
/* temperature scaling and pressure scaling to produce the extended variables at t+dt */
{
m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
- }
- else
- {
- if (bExchanged)
+
+ copy_mat(shake_vir, state->svir_prev);
+ copy_mat(force_vir, state->fvir_prev);
+ if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
{
- wallcycle_stop(wcycle, ewcUPDATE);
- /* We need the kinetic energy at minus the half step for determining
- * the full step kinetic energy and possibly for T-coupling.*/
- /* This may not be quite working correctly yet . . . . */
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
- wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
- constr, NULL, FALSE, state->box,
- top_global, &bSumEkinhOld,
- CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
- wallcycle_start(wcycle, ewcUPDATE);
+ /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
+ enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE);
+ enerd->term[F_EKIN] = trace(ekind->ekin);
}
}
- }
- if (bTrotter && !bInitStep)
- {
- copy_mat(shake_vir, state->svir_prev);
- copy_mat(force_vir, state->fvir_prev);
- if (IR_NVT_TROTTER(ir) && ir->eI == eiVV)
+ else if (bExchanged)
{
- /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
- enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE);
- enerd->term[F_EKIN] = trace(ekind->ekin);
+ wallcycle_stop(wcycle, ewcUPDATE);
+ /* We need the kinetic energy at minus the half step for determining
+ * the full step kinetic energy and possibly for T-coupling.*/
+ /* This may not be quite working correctly yet . . . . */
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+ constr, NULL, FALSE, state->box,
+ NULL, &bSumEkinhOld,
+ CGLO_GSTAT | CGLO_TEMPERATURE);
+ wallcycle_start(wcycle, ewcUPDATE);
}
}
/* if it's the initial step, we performed this first step just to get the constraint virial */
}
/* compute the conserved quantity */
- if (bVV)
+ if (EI_VV(ir->eI))
{
saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ);
if (ir->eI == eiVV)
*/
do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
ir, state, state_global, top_global, fr,
- outf, mdebin, ekind, f, f_global,
+ outf, mdebin, ekind, f,
&nchkpt,
bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
bSumEkinhOld);
bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
/* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
- if (bStartingFromCpt && bTrotter)
+ if (startingFromCheckpoint && bTrotter)
{
copy_mat(state->svir_prev, shake_vir);
copy_mat(state->fvir_prev, force_vir);
/* Check whether everything is still allright */
if (((int)gmx_get_stop_condition() > handled_stop_condition)
-#ifdef GMX_THREAD_MPI
+#if GMX_THREAD_MPI
&& MASTER(cr)
#endif
)
gs.sig[eglsCHKPT] = 1;
}
+ /* ######### START SECOND UPDATE STEP ################# */
+
/* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
in preprocessing */
TRUE, bCalcVir);
}
}
- /* ######### START SECOND UPDATE STEP ################# */
/* Box is changed in update() when we do pressure coupling,
* but we should still use the old box for energy corrections and when
* writing it to the energy file, so it matches the trajectory files for
update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep);
}
- if (bVV)
+ if (EI_VV(ir->eI))
{
- bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
-
/* velocity half-step update */
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
- ekind, M, upd, FALSE, etrtVELOCITY2,
- cr, nrnb, constr, &top->idef);
+ update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ ekind, M, upd, etrtVELOCITY2,
+ cr, constr);
}
/* Above, initialize just copies ekinh into ekin,
}
copy_rvecn(state->x, cbuf, 0, state->natoms);
}
- bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
- ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
+ update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
cr, nrnb, wcycle, upd, constr,
FALSE, bCalcVir);
- if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
- {
- /* Correct the virial for multiple time stepping */
- m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
- }
-
if (ir->eI == eiVVAK)
{
/* erase F_EKIN and F_TEMP here? */
/* just compute the kinetic energy at the half step to perform a trotter step */
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, NULL, FALSE, lastbox,
- top_global, &bSumEkinhOld,
- cglo_flags | CGLO_TEMPERATURE
+ NULL, &bSumEkinhOld,
+ (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
);
wallcycle_start(wcycle, ewcUPDATE);
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
/* now we know the scaling, we can compute the positions again again */
copy_rvecn(cbuf, state->x, 0, state->natoms);
- bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
-
- update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
- ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
+ update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
/* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
cr, nrnb, wcycle, upd, NULL,
FALSE, bCalcVir);
}
- if (bVV)
+ if (EI_VV(ir->eI))
{
/* this factor or 2 correction is necessary
because half of the constraint force is removed
*/
if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)))
{
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &gs,
(step_rel % gs.nstms == 0) &&
(multisim_nsteps < 0 || (step_rel < multisim_nsteps)),
lastbox,
- top_global, &bSumEkinhOld,
- cglo_flags
+ &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ (bGStat ? CGLO_GSTAT : 0)
| (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
| (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
| (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
| (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
| CGLO_CONSTRAINT
+ | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
);
+ checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
+ top_global, top, state,
+ &shouldCheckNumberOfBondedInteractions);
}
/* ############# END CALC EKIN AND PRESSURE ################# */
but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
- if (ir->efep != efepNO && (!bVV || bRerunMD))
+ if (ir->efep != efepNO && (!EI_VV(ir->eI) || bRerunMD))
{
/* Sum up the foreign energy and dhdl terms for md and sd.
Currently done every step so that dhdl is correct in the .edr */
}
enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
- if (bVV)
+ if (EI_VV(ir->eI))
{
enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
}
/* Output stuff */
if (MASTER(cr))
{
- gmx_bool do_dr, do_or;
-
if (fplog && do_log && bDoExpanded)
{
/* only needed if doing expanded ensemble */
PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
&state_global->dfhist, state->fep_state, ir->nstlog, step);
}
- if (!(startingFromCheckpoint && (EI_VV(ir->eI))))
+ if (bCalcEner)
{
- if (bCalcEner)
- {
- upd_mdebin(mdebin, bDoDHDL, TRUE,
- t, mdatoms->tmass, enerd, state,
- ir->fepvals, ir->expandedvals, lastbox,
- shake_vir, force_vir, total_vir, pres,
- ekind, mu_tot, constr);
- }
- else
- {
- upd_mdebin_step(mdebin);
- }
+ upd_mdebin(mdebin, bDoDHDL, bCalcEnerStep,
+ t, mdatoms->tmass, enerd, state,
+ ir->fepvals, ir->expandedvals, lastbox,
+ shake_vir, force_vir, total_vir, pres,
+ ekind, mu_tot, constr);
+ }
+ else
+ {
+ upd_mdebin_step(mdebin);
+ }
- do_dr = do_per_step(step, ir->nstdisreout);
- do_or = do_per_step(step, ir->nstorireout);
+ gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
+ gmx_bool do_or = do_per_step(step, ir->nstorireout);
+
+ print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : NULL,
+ step, t,
- eprNORMAL, bCompact, mdebin, fcd, groups, &(ir->opts));
++ eprNORMAL, mdebin, fcd, groups, &(ir->opts));
- print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : NULL,
- step, t,
- eprNORMAL, mdebin, fcd, groups, &(ir->opts));
- }
if (ir->bPull)
{
pull_print_output(ir->pull_work, step, t);
dd_partition_system(fplog, step, cr, TRUE, 1,
state_global, top_global, ir,
state, &f, mdatoms, top, fr,
- vsite, shellfc, constr,
+ vsite, constr,
nrnb, wcycle, FALSE);
+ shouldCheckNumberOfBondedInteractions = true;
}
- bFirstStep = FALSE;
- bInitStep = FALSE;
- bStartingFromCpt = FALSE;
+ bFirstStep = FALSE;
+ bInitStep = FALSE;
+ startingFromCheckpoint = FALSE;
/* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
/* With all integrators, except VV, we need to retain the pressure
}
/* End of main MD loop */
- debug_gmx();
/* Closing TNG files can include compressing data. Therefore it is good to do that
* before stopping the time measurements. */
if (ir->nstcalcenergy > 0 && !bRerunMD)
{
print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
- eprAVER, FALSE, mdebin, fcd, groups, &(ir->opts));
+ eprAVER, mdebin, fcd, groups, &(ir->opts));
}
}
done_mdoutf(outf);
- debug_gmx();
if (bPMETune)
{
pme_loadbal_done(pme_loadbal, cr, fplog, use_GPU(fr->nbv));
}
- if (shellfc && fplog)
- {
- fprintf(fplog, "Fraction of iterations that converged: %.2f %%\n",
- (nconverged*100.0)/step_rel);
- fprintf(fplog, "Average number of force evaluations per MD step: %.2f\n\n",
- tcount/step_rel);
- }
+ done_shellfc(fplog, shellfc, step_rel);
if (repl_ex_nst > 0 && MASTER(cr))
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff