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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * 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.
58 #include "md_support.h"
59 #include "md_logging.h"
74 #include "mpelogging.h"
76 #include "domdec_network.h"
82 #include "compute_io.h"
84 #include "checkpoint.h"
85 #include "mtop_util.h"
86 #include "sighandler.h"
89 #include "pme_loadbal.h"
92 #include "types/nlistheuristics.h"
93 #include "types/iteratedconstraints.h"
94 #include "nbnxn_cuda_data_mgmt.h"
104 #include "corewrap.h"
107 static void reset_all_counters(FILE *fplog, t_commrec *cr,
108 gmx_large_int_t step,
109 gmx_large_int_t *step_rel, t_inputrec *ir,
110 gmx_wallcycle_t wcycle, t_nrnb *nrnb,
111 gmx_runtime_t *runtime,
112 nbnxn_cuda_ptr_t cu_nbv)
114 char sbuf[STEPSTRSIZE];
116 /* Reset all the counters related to performance over the run */
117 md_print_warn(cr, fplog, "step %s: resetting all time and cycle counters\n",
118 gmx_step_str(step, sbuf));
122 nbnxn_cuda_reset_timings(cu_nbv);
125 wallcycle_stop(wcycle, ewcRUN);
126 wallcycle_reset_all(wcycle);
127 if (DOMAINDECOMP(cr))
129 reset_dd_statistics_counters(cr->dd);
132 ir->init_step += *step_rel;
133 ir->nsteps -= *step_rel;
135 wallcycle_start(wcycle, ewcRUN);
136 runtime_start(runtime);
137 print_date_and_time(fplog, cr->nodeid, "Restarted time", runtime);
140 double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
141 const output_env_t oenv, gmx_bool bVerbose, gmx_bool bCompact,
143 gmx_vsite_t *vsite, gmx_constr_t constr,
144 int stepout, t_inputrec *ir,
145 gmx_mtop_t *top_global,
147 t_state *state_global,
149 t_nrnb *nrnb, gmx_wallcycle_t wcycle,
150 gmx_edsam_t ed, t_forcerec *fr,
151 int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t membed,
152 real cpt_period, real max_hours,
153 const char *deviceOptions,
155 gmx_runtime_t *runtime)
158 gmx_large_int_t step, step_rel;
160 double t, t0, lam0[efptNR];
161 gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEner;
162 gmx_bool bNS, bNStList, bSimAnn, bStopCM, bRerunMD, bNotLastFrame = FALSE,
163 bFirstStep, bStateFromCP, bStateFromTPX, bInitStep, bLastStep,
164 bBornRadii, bStartingFromCpt;
165 gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
166 gmx_bool do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE,
167 bForceUpdate = FALSE, bCPT;
169 gmx_bool bMasterState;
170 int force_flags, cglo_flags;
171 tensor force_vir, shake_vir, total_vir, tmp_vir, pres;
176 t_state *bufstate = NULL;
177 matrix *scale_tot, pcoupl_mu, M, ebox;
180 gmx_repl_ex_t repl_ex = NULL;
183 t_mdebin *mdebin = NULL;
184 df_history_t df_history;
185 t_state *state = NULL;
186 rvec *f_global = NULL;
189 gmx_enerdata_t *enerd;
191 gmx_global_stat_t gstat;
192 gmx_update_t upd = NULL;
193 t_graph *graph = NULL;
195 gmx_rng_t mcrng = NULL;
197 gmx_groups_t *groups;
198 gmx_ekindata_t *ekind, *ekind_save;
199 gmx_shellfc_t shellfc;
200 int count, nconverged = 0;
203 gmx_bool bIonize = FALSE;
204 gmx_bool bTCR = FALSE, bConverged = TRUE, bOK, bSumEkinhOld, bExchanged;
206 gmx_bool bResetCountersHalfMaxH = FALSE;
207 gmx_bool bVV, bIterativeCase, bFirstIterate, bTemp, bPres, bTrotter;
208 gmx_bool bUpdateDoLR;
209 real mu_aver = 0, dvdl_constr;
210 int a0, a1, gnx = 0, ii;
211 atom_id *grpindex = NULL;
213 t_coupl_rec *tcr = NULL;
214 rvec *xcopy = NULL, *vcopy = NULL, *cbuf = NULL;
215 matrix boxcopy = {{0}}, lastbox;
217 real fom, oldfom, veta_save, pcurr, scalevir, tracevir;
224 real saved_conserved_quantity = 0;
229 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
230 int handled_stop_condition = gmx_stop_cond_none; /* compare to get_stop_condition*/
231 gmx_iterate_t iterate;
232 gmx_large_int_t multisim_nsteps = -1; /* number of steps to do before first multisim
233 simulation stops. If equal to zero, don't
234 communicate any more between multisims.*/
235 /* PME load balancing data for GPU kernels */
236 pme_load_balancing_t pme_loadbal = NULL;
238 gmx_bool bPMETuneTry = FALSE, bPMETuneRunning = FALSE;
241 /* Temporary addition for FAHCORE checkpointing */
245 /* Check for special mdrun options */
246 bRerunMD = (Flags & MD_RERUN);
247 bIonize = (Flags & MD_IONIZE);
248 bFFscan = (Flags & MD_FFSCAN);
249 bAppend = (Flags & MD_APPENDFILES);
250 if (Flags & MD_RESETCOUNTERSHALFWAY)
254 /* Signal to reset the counters half the simulation steps. */
255 wcycle_set_reset_counters(wcycle, ir->nsteps/2);
257 /* Signal to reset the counters halfway the simulation time. */
258 bResetCountersHalfMaxH = (max_hours > 0);
261 /* md-vv uses averaged full step velocities for T-control
262 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
263 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
265 if (bVV) /* to store the initial velocities while computing virial */
267 snew(cbuf, top_global->natoms);
269 /* all the iteratative cases - only if there are constraints */
270 bIterativeCase = ((IR_NPH_TROTTER(ir) || IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
271 gmx_iterate_init(&iterate, FALSE); /* The default value of iterate->bIterationActive is set to
272 false in this step. The correct value, true or false,
273 is set at each step, as it depends on the frequency of temperature
274 and pressure control.*/
275 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir)));
279 /* Since we don't know if the frames read are related in any way,
280 * rebuild the neighborlist at every step.
283 ir->nstcalcenergy = 1;
287 check_ir_old_tpx_versions(cr, fplog, ir, top_global);
289 nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
290 bGStatEveryStep = (nstglobalcomm == 1);
292 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
295 "To reduce the energy communication with nstlist = -1\n"
296 "the neighbor list validity should not be checked at every step,\n"
297 "this means that exact integration is not guaranteed.\n"
298 "The neighbor list validity is checked after:\n"
299 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
300 "In most cases this will result in exact integration.\n"
301 "This reduces the energy communication by a factor of 2 to 3.\n"
302 "If you want less energy communication, set nstlist > 3.\n\n");
305 if (bRerunMD || bFFscan)
309 groups = &top_global->groups;
312 init_md(fplog, cr, ir, oenv, &t, &t0, state_global->lambda,
313 &(state_global->fep_state), lam0,
314 nrnb, top_global, &upd,
315 nfile, fnm, &outf, &mdebin,
316 force_vir, shake_vir, mu_tot, &bSimAnn, &vcm, state_global, Flags);
318 clear_mat(total_vir);
320 /* Energy terms and groups */
322 init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
324 if (DOMAINDECOMP(cr))
330 snew(f, top_global->natoms);
333 /* lambda Monte carlo random number generator */
336 mcrng = gmx_rng_init(ir->expandedvals->lmc_seed);
338 /* copy the state into df_history */
339 copy_df_history(&df_history, &state_global->dfhist);
341 /* Kinetic energy data */
343 init_ekindata(fplog, top_global, &(ir->opts), ekind);
344 /* needed for iteration of constraints */
346 init_ekindata(fplog, top_global, &(ir->opts), ekind_save);
347 /* Copy the cos acceleration to the groups struct */
348 ekind->cosacc.cos_accel = ir->cos_accel;
350 gstat = global_stat_init(ir);
353 /* Check for polarizable models and flexible constraints */
354 shellfc = init_shell_flexcon(fplog,
355 top_global, n_flexible_constraints(constr),
356 (ir->bContinuation ||
357 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
358 NULL : state_global->x);
362 #ifdef GMX_THREAD_MPI
363 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
365 set_deform_reference_box(upd,
366 deform_init_init_step_tpx,
367 deform_init_box_tpx);
368 #ifdef GMX_THREAD_MPI
369 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
374 double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1);
375 if ((io > 2000) && MASTER(cr))
378 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
383 if (DOMAINDECOMP(cr))
385 top = dd_init_local_top(top_global);
388 dd_init_local_state(cr->dd, state_global, state);
390 if (DDMASTER(cr->dd) && ir->nstfout)
392 snew(f_global, state_global->natoms);
399 /* Initialize the particle decomposition and split the topology */
400 top = split_system(fplog, top_global, ir, cr);
402 pd_cg_range(cr, &fr->cg0, &fr->hcg);
403 pd_at_range(cr, &a0, &a1);
407 top = gmx_mtop_generate_local_top(top_global, ir);
410 a1 = top_global->natoms;
413 forcerec_set_excl_load(fr, top, cr);
415 state = partdec_init_local_state(cr, state_global);
418 atoms2md(top_global, ir, 0, NULL, a0, a1-a0, mdatoms);
422 set_vsite_top(vsite, top, mdatoms, cr);
425 if (ir->ePBC != epbcNONE && !fr->bMolPBC)
427 graph = mk_graph(fplog, &(top->idef), 0, top_global->natoms, FALSE, FALSE);
432 make_local_shells(cr, mdatoms, shellfc);
435 init_bonded_thread_force_reduction(fr, &top->idef);
437 if (ir->pull && PAR(cr))
439 dd_make_local_pull_groups(NULL, ir->pull, mdatoms);
443 if (DOMAINDECOMP(cr))
445 /* Distribute the charge groups over the nodes from the master node */
446 dd_partition_system(fplog, ir->init_step, cr, TRUE, 1,
447 state_global, top_global, ir,
448 state, &f, mdatoms, top, fr,
449 vsite, shellfc, constr,
450 nrnb, wcycle, FALSE);
454 update_mdatoms(mdatoms, state->lambda[efptMASS]);
456 if (opt2bSet("-cpi", nfile, fnm))
458 bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi", nfile, fnm, cr), cr);
462 bStateFromCP = FALSE;
469 /* Update mdebin with energy history if appending to output files */
470 if (Flags & MD_APPENDFILES)
472 restore_energyhistory_from_state(mdebin, &state_global->enerhist);
476 /* We might have read an energy history from checkpoint,
477 * free the allocated memory and reset the counts.
479 done_energyhistory(&state_global->enerhist);
480 init_energyhistory(&state_global->enerhist);
483 /* Set the initial energy history in state by updating once */
484 update_energyhistory(&state_global->enerhist, mdebin);
487 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG))
489 /* Set the random state if we read a checkpoint file */
490 set_stochd_state(upd, state);
493 if (state->flags & (1<<estMC_RNG))
495 set_mc_state(mcrng, state);
498 /* Initialize constraints */
501 if (!DOMAINDECOMP(cr))
503 set_constraints(constr, top, ir, mdatoms, cr);
507 /* Check whether we have to GCT stuff */
508 bTCR = ftp2bSet(efGCT, nfile, fnm);
513 fprintf(stderr, "Will do General Coupling Theory!\n");
515 gnx = top_global->mols.nr;
517 for (i = 0; (i < gnx); i++)
525 /* We need to be sure replica exchange can only occur
526 * when the energies are current */
527 check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
528 "repl_ex_nst", &repl_ex_nst);
529 /* This check needs to happen before inter-simulation
530 * signals are initialized, too */
532 if (repl_ex_nst > 0 && MASTER(cr))
534 repl_ex = init_replica_exchange(fplog, cr->ms, state_global, ir,
535 repl_ex_nst, repl_ex_nex, repl_ex_seed);
538 /* PME tuning is only supported with GPUs or PME nodes and not with rerun.
539 * With perturbed charges with soft-core we should not change the cut-off.
541 if ((Flags & MD_TUNEPME) &&
542 EEL_PME(fr->eeltype) &&
543 ( (fr->cutoff_scheme == ecutsVERLET && fr->nbv->bUseGPU) || !(cr->duty & DUTY_PME)) &&
544 !(ir->efep != efepNO && mdatoms->nChargePerturbed > 0 && ir->fepvals->bScCoul) &&
547 pme_loadbal_init(&pme_loadbal, ir, state->box, fr->ic, fr->pmedata);
549 if (cr->duty & DUTY_PME)
551 /* Start tuning right away, as we can't measure the load */
552 bPMETuneRunning = TRUE;
556 /* Separate PME nodes, we can measure the PP/PME load balance */
561 if (!ir->bContinuation && !bRerunMD)
563 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
565 /* Set the velocities of frozen particles to zero */
566 for (i = mdatoms->start; i < mdatoms->start+mdatoms->homenr; i++)
568 for (m = 0; m < DIM; m++)
570 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
580 /* Constrain the initial coordinates and velocities */
581 do_constrain_first(fplog, constr, ir, mdatoms, state, f,
582 graph, cr, nrnb, fr, top, shake_vir);
586 /* Construct the virtual sites for the initial configuration */
587 construct_vsites(fplog, vsite, state->x, nrnb, ir->delta_t, NULL,
588 top->idef.iparams, top->idef.il,
589 fr->ePBC, fr->bMolPBC, graph, cr, state->box);
595 /* set free energy calculation frequency as the minimum
596 greatest common denominator of nstdhdl, nstexpanded, and repl_ex_nst*/
597 nstfep = ir->fepvals->nstdhdl;
600 nstfep = gmx_greatest_common_divisor(ir->fepvals->nstdhdl,nstfep);
604 nstfep = gmx_greatest_common_divisor(repl_ex_nst,nstfep);
607 /* I'm assuming we need global communication the first time! MRS */
608 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
609 | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM : 0)
610 | (bVV ? CGLO_PRESSURE : 0)
611 | (bVV ? CGLO_CONSTRAINT : 0)
612 | (bRerunMD ? CGLO_RERUNMD : 0)
613 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN : 0));
615 bSumEkinhOld = FALSE;
616 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
617 NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
618 constr, NULL, FALSE, state->box,
619 top_global, &pcurr, top_global->natoms, &bSumEkinhOld, cglo_flags);
620 if (ir->eI == eiVVAK)
622 /* a second call to get the half step temperature initialized as well */
623 /* we do the same call as above, but turn the pressure off -- internally to
624 compute_globals, this is recognized as a velocity verlet half-step
625 kinetic energy calculation. This minimized excess variables, but
626 perhaps loses some logic?*/
628 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
629 NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
630 constr, NULL, FALSE, state->box,
631 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
632 cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
635 /* Calculate the initial half step temperature, and save the ekinh_old */
636 if (!(Flags & MD_STARTFROMCPT))
638 for (i = 0; (i < ir->opts.ngtc); i++)
640 copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
645 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
646 and there is no previous step */
649 /* if using an iterative algorithm, we need to create a working directory for the state. */
652 bufstate = init_bufstate(state);
656 snew(xcopy, state->natoms);
657 snew(vcopy, state->natoms);
658 copy_rvecn(state->x, xcopy, 0, state->natoms);
659 copy_rvecn(state->v, vcopy, 0, state->natoms);
660 copy_mat(state->box, boxcopy);
663 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
664 temperature control */
665 trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
669 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
672 "RMS relative constraint deviation after constraining: %.2e\n",
673 constr_rmsd(constr, FALSE));
675 if (EI_STATE_VELOCITY(ir->eI))
677 fprintf(fplog, "Initial temperature: %g K\n", enerd->term[F_TEMP]);
681 fprintf(stderr, "starting md rerun '%s', reading coordinates from"
682 " input trajectory '%s'\n\n",
683 *(top_global->name), opt2fn("-rerun", nfile, fnm));
686 fprintf(stderr, "Calculated time to finish depends on nsteps from "
687 "run input file,\nwhich may not correspond to the time "
688 "needed to process input trajectory.\n\n");
694 fprintf(stderr, "starting mdrun '%s'\n",
695 *(top_global->name));
698 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
702 sprintf(tbuf, "%s", "infinite");
704 if (ir->init_step > 0)
706 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
707 gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
708 gmx_step_str(ir->init_step, sbuf2),
709 ir->init_step*ir->delta_t);
713 fprintf(stderr, "%s steps, %s ps.\n",
714 gmx_step_str(ir->nsteps, sbuf), tbuf);
717 fprintf(fplog, "\n");
720 /* Set and write start time */
721 runtime_start(runtime);
722 print_date_and_time(fplog, cr->nodeid, "Started mdrun", runtime);
723 wallcycle_start(wcycle, ewcRUN);
726 fprintf(fplog, "\n");
729 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
731 chkpt_ret = fcCheckPointParallel( cr->nodeid,
735 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
740 /***********************************************************
744 ************************************************************/
746 /* if rerunMD then read coordinates and velocities from input trajectory */
749 if (getenv("GMX_FORCE_UPDATE"))
757 bNotLastFrame = read_first_frame(oenv, &status,
758 opt2fn("-rerun", nfile, fnm),
759 &rerun_fr, TRX_NEED_X | TRX_READ_V);
760 if (rerun_fr.natoms != top_global->natoms)
763 "Number of atoms in trajectory (%d) does not match the "
764 "run input file (%d)\n",
765 rerun_fr.natoms, top_global->natoms);
767 if (ir->ePBC != epbcNONE)
771 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);
773 if (max_cutoff2(ir->ePBC, rerun_fr.box) < sqr(fr->rlistlong))
775 gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr.step, rerun_fr.time);
782 rerun_parallel_comm(cr, &rerun_fr, &bNotLastFrame);
785 if (ir->ePBC != epbcNONE)
787 /* Set the shift vectors.
788 * Necessary here when have a static box different from the tpr box.
790 calc_shifts(rerun_fr.box, fr->shift_vec);
794 /* loop over MD steps or if rerunMD to end of input trajectory */
796 /* Skip the first Nose-Hoover integration when we get the state from tpx */
797 bStateFromTPX = !bStateFromCP;
798 bInitStep = bFirstStep && (bStateFromTPX || bVV);
799 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
801 bSumEkinhOld = FALSE;
804 init_global_signals(&gs, cr, ir, repl_ex_nst);
806 step = ir->init_step;
809 if (ir->nstlist == -1)
811 init_nlistheuristics(&nlh, bGStatEveryStep, step);
814 if (MULTISIM(cr) && (repl_ex_nst <= 0 ))
816 /* check how many steps are left in other sims */
817 multisim_nsteps = get_multisim_nsteps(cr, ir->nsteps);
821 /* and stop now if we should */
822 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
823 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
824 while (!bLastStep || (bRerunMD && bNotLastFrame))
827 wallcycle_start(wcycle, ewcSTEP);
829 GMX_MPE_LOG(ev_timestep1);
835 step = rerun_fr.step;
836 step_rel = step - ir->init_step;
849 bLastStep = (step_rel == ir->nsteps);
850 t = t0 + step*ir->delta_t;
853 if (ir->efep != efepNO || ir->bSimTemp)
855 /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
856 requiring different logic. */
858 set_current_lambdas(step, ir->fepvals, bRerunMD, &rerun_fr, state_global, state, lam0);
859 bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
860 bDoFEP = (do_per_step(step, nstfep) && (ir->efep != efepNO));
861 bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded) && (ir->bExpanded) && (step > 0));
866 update_annealing_target_temp(&(ir->opts), t);
871 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
873 for (i = 0; i < state_global->natoms; i++)
875 copy_rvec(rerun_fr.x[i], state_global->x[i]);
879 for (i = 0; i < state_global->natoms; i++)
881 copy_rvec(rerun_fr.v[i], state_global->v[i]);
886 for (i = 0; i < state_global->natoms; i++)
888 clear_rvec(state_global->v[i]);
892 fprintf(stderr, "\nWARNING: Some frames do not contain velocities.\n"
893 " Ekin, temperature and pressure are incorrect,\n"
894 " the virial will be incorrect when constraints are present.\n"
896 bRerunWarnNoV = FALSE;
900 copy_mat(rerun_fr.box, state_global->box);
901 copy_mat(state_global->box, state->box);
903 if (vsite && (Flags & MD_RERUN_VSITE))
905 if (DOMAINDECOMP(cr))
907 gmx_fatal(FARGS, "Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
911 /* Following is necessary because the graph may get out of sync
912 * with the coordinates if we only have every N'th coordinate set
914 mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
915 shift_self(graph, state->box, state->x);
917 construct_vsites(fplog, vsite, state->x, nrnb, ir->delta_t, state->v,
918 top->idef.iparams, top->idef.il,
919 fr->ePBC, fr->bMolPBC, graph, cr, state->box);
922 unshift_self(graph, state->box, state->x);
927 /* Stop Center of Mass motion */
928 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
930 /* Copy back starting coordinates in case we're doing a forcefield scan */
933 for (ii = 0; (ii < state->natoms); ii++)
935 copy_rvec(xcopy[ii], state->x[ii]);
936 copy_rvec(vcopy[ii], state->v[ii]);
938 copy_mat(boxcopy, state->box);
943 /* for rerun MD always do Neighbour Searching */
944 bNS = (bFirstStep || ir->nstlist != 0);
949 /* Determine whether or not to do Neighbour Searching and LR */
950 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
952 bNS = (bFirstStep || bExchanged || bNStList || bDoFEP ||
953 (ir->nstlist == -1 && nlh.nabnsb > 0));
955 if (bNS && ir->nstlist == -1)
957 set_nlistheuristics(&nlh, bFirstStep || bExchanged || bDoFEP, step);
961 /* check whether we should stop because another simulation has
965 if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
966 (multisim_nsteps != ir->nsteps) )
973 "Stopping simulation %d because another one has finished\n",
977 gs.sig[eglsCHKPT] = 1;
982 /* < 0 means stop at next step, > 0 means stop at next NS step */
983 if ( (gs.set[eglsSTOPCOND] < 0) ||
984 ( (gs.set[eglsSTOPCOND] > 0) && (bNStList || ir->nstlist == 0) ) )
989 /* Determine whether or not to update the Born radii if doing GB */
990 bBornRadii = bFirstStep;
991 if (ir->implicit_solvent && (step % ir->nstgbradii == 0))
996 do_log = do_per_step(step, ir->nstlog) || bFirstStep || bLastStep;
997 do_verbose = bVerbose &&
998 (step % stepout == 0 || bFirstStep || bLastStep);
1000 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
1004 bMasterState = TRUE;
1008 bMasterState = FALSE;
1009 /* Correct the new box if it is too skewed */
1010 if (DYNAMIC_BOX(*ir))
1012 if (correct_box(fplog, step, state->box, graph))
1014 bMasterState = TRUE;
1017 if (DOMAINDECOMP(cr) && bMasterState)
1019 dd_collect_state(cr->dd, state, state_global);
1023 if (DOMAINDECOMP(cr))
1025 /* Repartition the domain decomposition */
1026 wallcycle_start(wcycle, ewcDOMDEC);
1027 dd_partition_system(fplog, step, cr,
1028 bMasterState, nstglobalcomm,
1029 state_global, top_global, ir,
1030 state, &f, mdatoms, top, fr,
1031 vsite, shellfc, constr,
1033 do_verbose && !bPMETuneRunning);
1034 wallcycle_stop(wcycle, ewcDOMDEC);
1035 /* If using an iterative integrator, reallocate space to match the decomposition */
1039 if (MASTER(cr) && do_log && !bFFscan)
1041 print_ebin_header(fplog, step, t, state->lambda[efptFEP]); /* can we improve the information printed here? */
1044 if (ir->efep != efepNO)
1046 update_mdatoms(mdatoms, state->lambda[efptMASS]);
1049 if ((bRerunMD && rerun_fr.bV) || bExchanged)
1052 /* We need the kinetic energy at minus the half step for determining
1053 * the full step kinetic energy and possibly for T-coupling.*/
1054 /* This may not be quite working correctly yet . . . . */
1055 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1056 wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
1057 constr, NULL, FALSE, state->box,
1058 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1059 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1061 clear_mat(force_vir);
1063 /* Ionize the atoms if necessary */
1066 ionize(fplog, oenv, mdatoms, top_global, t, ir, state->x, state->v,
1067 mdatoms->start, mdatoms->start+mdatoms->homenr, state->box, cr);
1070 /* Update force field in ffscan program */
1073 if (update_forcefield(fplog,
1075 mdatoms->nr, state->x, state->box))
1083 GMX_MPE_LOG(ev_timestep2);
1085 /* We write a checkpoint at this MD step when:
1086 * either at an NS step when we signalled through gs,
1087 * or at the last step (but not when we do not want confout),
1088 * but never at the first step or with rerun.
1090 bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
1091 (bLastStep && (Flags & MD_CONFOUT))) &&
1092 step > ir->init_step && !bRerunMD);
1095 gs.set[eglsCHKPT] = 0;
1098 /* Determine the energy and pressure:
1099 * at nstcalcenergy steps and at energy output steps (set below).
1101 if (EI_VV(ir->eI) && (!bInitStep))
1103 /* for vv, the first half of the integration actually corresponds
1104 to the previous step. bCalcEner is only required to be evaluated on the 'next' step,
1105 but the virial needs to be calculated on both the current step and the 'next' step. Future
1106 reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
1108 bCalcEner = do_per_step(step-1, ir->nstcalcenergy);
1109 bCalcVir = bCalcEner ||
1110 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
1114 bCalcEner = do_per_step(step, ir->nstcalcenergy);
1115 bCalcVir = bCalcEner ||
1116 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
1119 /* Do we need global communication ? */
1120 bGStat = (bCalcVir || bCalcEner || bStopCM ||
1121 do_per_step(step, nstglobalcomm) || (bVV && IR_NVT_TROTTER(ir) && do_per_step(step-1, nstglobalcomm)) ||
1122 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1124 do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
1126 if (do_ene || do_log)
1133 /* these CGLO_ options remain the same throughout the iteration */
1134 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1135 (bGStat ? CGLO_GSTAT : 0)
1138 force_flags = (GMX_FORCE_STATECHANGED |
1139 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1140 GMX_FORCE_ALLFORCES |
1142 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
1143 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
1144 (bDoFEP ? GMX_FORCE_DHDL : 0)
1149 if (do_per_step(step, ir->nstcalclr))
1151 force_flags |= GMX_FORCE_DO_LR;
1157 /* Now is the time to relax the shells */
1158 count = relax_shell_flexcon(fplog, cr, bVerbose, bFFscan ? step+1 : step,
1159 ir, bNS, force_flags,
1160 bStopCM, top, top_global,
1162 state, f, force_vir, mdatoms,
1163 nrnb, wcycle, graph, groups,
1164 shellfc, fr, bBornRadii, t, mu_tot,
1165 state->natoms, &bConverged, vsite,
1176 /* The coordinates (x) are shifted (to get whole molecules)
1178 * This is parallellized as well, and does communication too.
1179 * Check comments in sim_util.c
1181 do_force(fplog, cr, ir, step, nrnb, wcycle, top, top_global, groups,
1182 state->box, state->x, &state->hist,
1183 f, force_vir, mdatoms, enerd, fcd,
1184 state->lambda, graph,
1185 fr, vsite, mu_tot, t, outf->fp_field, ed, bBornRadii,
1186 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1189 GMX_BARRIER(cr->mpi_comm_mygroup);
1193 mu_aver = calc_mu_aver(cr, state->x, mdatoms->chargeA,
1194 mu_tot, &top_global->mols, mdatoms, gnx, grpindex);
1197 if (bTCR && bFirstStep)
1199 tcr = init_coupling(fplog, nfile, fnm, cr, fr, mdatoms, &(top->idef));
1200 fprintf(fplog, "Done init_coupling\n");
1204 if (bVV && !bStartingFromCpt && !bRerunMD)
1205 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1207 if (ir->eI == eiVV && bInitStep)
1209 /* if using velocity verlet with full time step Ekin,
1210 * take the first half step only to compute the
1211 * virial for the first step. From there,
1212 * revert back to the initial coordinates
1213 * so that the input is actually the initial step.
1215 copy_rvecn(state->v, cbuf, 0, state->natoms); /* should make this better for parallelizing? */
1219 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1220 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
1223 /* If we are using twin-range interactions where the long-range component
1224 * is only evaluated every nstcalclr>1 steps, we should do a special update
1225 * step to combine the long-range forces on these steps.
1226 * For nstcalclr=1 this is not done, since the forces would have been added
1227 * directly to the short-range forces already.
1229 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
1231 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC,
1232 f, bUpdateDoLR, fr->f_twin, fcd,
1233 ekind, M, wcycle, upd, bInitStep, etrtVELOCITY1,
1234 cr, nrnb, constr, &top->idef);
1236 if (bIterativeCase && do_per_step(step-1, ir->nstpcouple) && !bInitStep)
1238 gmx_iterate_init(&iterate, TRUE);
1240 /* for iterations, we save these vectors, as we will be self-consistently iterating
1243 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1245 /* save the state */
1246 if (iterate.bIterationActive)
1248 copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts));
1251 bFirstIterate = TRUE;
1252 while (bFirstIterate || iterate.bIterationActive)
1254 if (iterate.bIterationActive)
1256 copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts));
1257 if (bFirstIterate && bTrotter)
1259 /* The first time through, we need a decent first estimate
1260 of veta(t+dt) to compute the constraints. Do
1261 this by computing the box volume part of the
1262 trotter integration at this time. Nothing else
1263 should be changed by this routine here. If
1264 !(first time), we start with the previous value
1267 veta_save = state->veta;
1268 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ0);
1269 vetanew = state->veta;
1270 state->veta = veta_save;
1275 if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
1277 update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
1278 state, fr->bMolPBC, graph, f,
1279 &top->idef, shake_vir, NULL,
1280 cr, nrnb, wcycle, upd, constr,
1281 bInitStep, TRUE, bCalcVir, vetanew);
1283 if (!bOK && !bFFscan)
1285 gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
1291 /* Need to unshift here if a do_force has been
1292 called in the previous step */
1293 unshift_self(graph, state->box, state->x);
1296 /* if VV, compute the pressure and constraints */
1297 /* For VV2, we strictly only need this if using pressure
1298 * control, but we really would like to have accurate pressures
1300 * Think about ways around this in the future?
1301 * For now, keep this choice in comments.
1303 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1304 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1306 bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
1307 if (bCalcEner && ir->eI == eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/
1309 bSumEkinhOld = TRUE;
1311 /* for vv, the first half of the integration actually corresponds to the previous step.
1312 So we need information from the last step in the first half of the integration */
1313 if (bGStat || do_per_step(step-1, nstglobalcomm))
1315 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1316 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1317 constr, NULL, FALSE, state->box,
1318 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1321 | (bTemp ? CGLO_TEMPERATURE : 0)
1322 | (bPres ? CGLO_PRESSURE : 0)
1323 | (bPres ? CGLO_CONSTRAINT : 0)
1324 | ((iterate.bIterationActive) ? CGLO_ITERATE : 0)
1325 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1328 /* explanation of above:
1329 a) We compute Ekin at the full time step
1330 if 1) we are using the AveVel Ekin, and it's not the
1331 initial step, or 2) if we are using AveEkin, but need the full
1332 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1333 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1334 EkinAveVel because it's needed for the pressure */
1336 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1341 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
1342 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
1349 /* We need the kinetic energy at minus the half step for determining
1350 * the full step kinetic energy and possibly for T-coupling.*/
1351 /* This may not be quite working correctly yet . . . . */
1352 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1353 wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
1354 constr, NULL, FALSE, state->box,
1355 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1356 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1361 if (iterate.bIterationActive &&
1362 done_iterating(cr, fplog, step, &iterate, bFirstIterate,
1363 state->veta, &vetanew))
1367 bFirstIterate = FALSE;
1370 if (bTrotter && !bInitStep)
1372 copy_mat(shake_vir, state->svir_prev);
1373 copy_mat(force_vir, state->fvir_prev);
1374 if (IR_NVT_TROTTER(ir) && ir->eI == eiVV)
1376 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1377 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE, FALSE);
1378 enerd->term[F_EKIN] = trace(ekind->ekin);
1381 /* if it's the initial step, we performed this first step just to get the constraint virial */
1382 if (bInitStep && ir->eI == eiVV)
1384 copy_rvecn(cbuf, state->v, 0, state->natoms);
1387 GMX_MPE_LOG(ev_timestep1);
1390 /* MRS -- now done iterating -- compute the conserved quantity */
1393 saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ);
1396 last_ekin = enerd->term[F_EKIN];
1398 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1400 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1402 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1405 sum_dhdl(enerd, state->lambda, ir->fepvals);
1409 /* ######## END FIRST UPDATE STEP ############## */
1410 /* ######## If doing VV, we now have v(dt) ###### */
1413 /* perform extended ensemble sampling in lambda - we don't
1414 actually move to the new state before outputting
1415 statistics, but if performing simulated tempering, we
1416 do update the velocities and the tau_t. */
1418 lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, &df_history, step, mcrng, state->v, mdatoms);
1420 /* ################## START TRAJECTORY OUTPUT ################# */
1422 /* Now we have the energies and forces corresponding to the
1423 * coordinates at time t. We must output all of this before
1425 * for RerunMD t is read from input trajectory
1427 GMX_MPE_LOG(ev_output_start);
1430 if (do_per_step(step, ir->nstxout))
1432 mdof_flags |= MDOF_X;
1434 if (do_per_step(step, ir->nstvout))
1436 mdof_flags |= MDOF_V;
1438 if (do_per_step(step, ir->nstfout))
1440 mdof_flags |= MDOF_F;
1442 if (do_per_step(step, ir->nstxtcout))
1444 mdof_flags |= MDOF_XTC;
1448 mdof_flags |= MDOF_CPT;
1452 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1455 /* Enforce writing positions and velocities at end of run */
1456 mdof_flags |= (MDOF_X | MDOF_V);
1462 fcReportProgress( ir->nsteps, step );
1465 /* sync bCPT and fc record-keeping */
1466 if (bCPT && MASTER(cr))
1468 fcRequestCheckPoint();
1472 if (mdof_flags != 0)
1474 wallcycle_start(wcycle, ewcTRAJ);
1477 if (state->flags & (1<<estLD_RNG))
1479 get_stochd_state(upd, state);
1481 if (state->flags & (1<<estMC_RNG))
1483 get_mc_state(mcrng, state);
1489 state_global->ekinstate.bUpToDate = FALSE;
1493 update_ekinstate(&state_global->ekinstate, ekind);
1494 state_global->ekinstate.bUpToDate = TRUE;
1496 update_energyhistory(&state_global->enerhist, mdebin);
1497 if (ir->efep != efepNO || ir->bSimTemp)
1499 state_global->fep_state = state->fep_state; /* MRS: seems kludgy. The code should be
1500 structured so this isn't necessary.
1501 Note this reassignment is only necessary
1502 for single threads.*/
1503 copy_df_history(&state_global->dfhist, &df_history);
1507 write_traj(fplog, cr, outf, mdof_flags, top_global,
1508 step, t, state, state_global, f, f_global, &n_xtc, &x_xtc);
1515 if (bLastStep && step_rel == ir->nsteps &&
1516 (Flags & MD_CONFOUT) && MASTER(cr) &&
1517 !bRerunMD && !bFFscan)
1519 /* x and v have been collected in write_traj,
1520 * because a checkpoint file will always be written
1523 fprintf(stderr, "\nWriting final coordinates.\n");
1526 /* Make molecules whole only for confout writing */
1527 do_pbc_mtop(fplog, ir->ePBC, state->box, top_global, state_global->x);
1529 write_sto_conf_mtop(ftp2fn(efSTO, nfile, fnm),
1530 *top_global->name, top_global,
1531 state_global->x, state_global->v,
1532 ir->ePBC, state->box);
1535 wallcycle_stop(wcycle, ewcTRAJ);
1537 GMX_MPE_LOG(ev_output_finish);
1539 /* kludge -- virial is lost with restart for NPT control. Must restart */
1540 if (bStartingFromCpt && bVV)
1542 copy_mat(state->svir_prev, shake_vir);
1543 copy_mat(state->fvir_prev, force_vir);
1545 /* ################## END TRAJECTORY OUTPUT ################ */
1547 /* Determine the wallclock run time up till now */
1548 run_time = gmx_gettime() - (double)runtime->real;
1550 /* Check whether everything is still allright */
1551 if (((int)gmx_get_stop_condition() > handled_stop_condition)
1552 #ifdef GMX_THREAD_MPI
1557 /* this is just make gs.sig compatible with the hack
1558 of sending signals around by MPI_Reduce with together with
1560 if (gmx_get_stop_condition() == gmx_stop_cond_next_ns)
1562 gs.sig[eglsSTOPCOND] = 1;
1564 if (gmx_get_stop_condition() == gmx_stop_cond_next)
1566 gs.sig[eglsSTOPCOND] = -1;
1568 /* < 0 means stop at next step, > 0 means stop at next NS step */
1572 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1573 gmx_get_signal_name(),
1574 gs.sig[eglsSTOPCOND] == 1 ? "NS " : "");
1578 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1579 gmx_get_signal_name(),
1580 gs.sig[eglsSTOPCOND] == 1 ? "NS " : "");
1582 handled_stop_condition = (int)gmx_get_stop_condition();
1584 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1585 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
1586 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
1588 /* Signal to terminate the run */
1589 gs.sig[eglsSTOPCOND] = 1;
1592 fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
1594 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
1597 if (bResetCountersHalfMaxH && MASTER(cr) &&
1598 run_time > max_hours*60.0*60.0*0.495)
1600 gs.sig[eglsRESETCOUNTERS] = 1;
1603 if (ir->nstlist == -1 && !bRerunMD)
1605 /* When bGStatEveryStep=FALSE, global_stat is only called
1606 * when we check the atom displacements, not at NS steps.
1607 * This means that also the bonded interaction count check is not
1608 * performed immediately after NS. Therefore a few MD steps could
1609 * be performed with missing interactions.
1610 * But wrong energies are never written to file,
1611 * since energies are only written after global_stat
1614 if (step >= nlh.step_nscheck)
1616 nlh.nabnsb = natoms_beyond_ns_buffer(ir, fr, &top->cgs,
1617 nlh.scale_tot, state->x);
1621 /* This is not necessarily true,
1622 * but step_nscheck is determined quite conservatively.
1628 /* In parallel we only have to check for checkpointing in steps
1629 * where we do global communication,
1630 * otherwise the other nodes don't know.
1632 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1635 run_time >= nchkpt*cpt_period*60.0)) &&
1636 gs.set[eglsCHKPT] == 0)
1638 gs.sig[eglsCHKPT] = 1;
1641 /* at the start of step, randomize or scale the velocities (trotter done elsewhere) */
1646 update_tcouple(fplog, step, ir, state, ekind, wcycle, upd, &MassQ, mdatoms);
1648 if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1650 gmx_bool bIfRandomize;
1651 bIfRandomize = update_randomize_velocities(ir, step, mdatoms, state, upd, &top->idef, constr);
1652 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1653 if (constr && bIfRandomize)
1655 update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
1656 state, fr->bMolPBC, graph, f,
1657 &top->idef, tmp_vir, NULL,
1658 cr, nrnb, wcycle, upd, constr,
1659 bInitStep, TRUE, bCalcVir, vetanew);
1664 if (bIterativeCase && do_per_step(step, ir->nstpcouple))
1666 gmx_iterate_init(&iterate, TRUE);
1667 /* for iterations, we save these vectors, as we will be redoing the calculations */
1668 copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts));
1671 bFirstIterate = TRUE;
1672 while (bFirstIterate || iterate.bIterationActive)
1674 /* We now restore these vectors to redo the calculation with improved extended variables */
1675 if (iterate.bIterationActive)
1677 copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts));
1680 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1681 so scroll down for that logic */
1683 /* ######### START SECOND UPDATE STEP ################# */
1684 GMX_MPE_LOG(ev_update_start);
1685 /* Box is changed in update() when we do pressure coupling,
1686 * but we should still use the old box for energy corrections and when
1687 * writing it to the energy file, so it matches the trajectory files for
1688 * the same timestep above. Make a copy in a separate array.
1690 copy_mat(state->box, lastbox);
1695 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
1697 wallcycle_start(wcycle, ewcUPDATE);
1698 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1701 if (iterate.bIterationActive)
1709 /* we use a new value of scalevir to converge the iterations faster */
1710 scalevir = tracevir/trace(shake_vir);
1712 msmul(shake_vir, scalevir, shake_vir);
1713 m_add(force_vir, shake_vir, total_vir);
1714 clear_mat(shake_vir);
1716 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1717 /* We can only do Berendsen coupling after we have summed
1718 * the kinetic energy or virial. Since the happens
1719 * in global_state after update, we should only do it at
1720 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1725 update_tcouple(fplog, step, ir, state, ekind, wcycle, upd, &MassQ, mdatoms);
1726 update_pcouple(fplog, step, ir, state, pcoupl_mu, M, wcycle,
1732 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
1734 /* velocity half-step update */
1735 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
1736 bUpdateDoLR, fr->f_twin, fcd,
1737 ekind, M, wcycle, upd, FALSE, etrtVELOCITY2,
1738 cr, nrnb, constr, &top->idef);
1741 /* Above, initialize just copies ekinh into ekin,
1742 * it doesn't copy position (for VV),
1743 * and entire integrator for MD.
1746 if (ir->eI == eiVVAK)
1748 copy_rvecn(state->x, cbuf, 0, state->natoms);
1750 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
1752 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
1753 bUpdateDoLR, fr->f_twin, fcd,
1754 ekind, M, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
1755 wallcycle_stop(wcycle, ewcUPDATE);
1757 update_constraints(fplog, step, &dvdl_constr, ir, ekind, mdatoms, state,
1758 fr->bMolPBC, graph, f,
1759 &top->idef, shake_vir, force_vir,
1760 cr, nrnb, wcycle, upd, constr,
1761 bInitStep, FALSE, bCalcVir, state->veta);
1763 if (ir->eI == eiVVAK)
1765 /* erase F_EKIN and F_TEMP here? */
1766 /* just compute the kinetic energy at the half step to perform a trotter step */
1767 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1768 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1769 constr, NULL, FALSE, lastbox,
1770 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1771 cglo_flags | CGLO_TEMPERATURE
1773 wallcycle_start(wcycle, ewcUPDATE);
1774 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1775 /* now we know the scaling, we can compute the positions again again */
1776 copy_rvecn(cbuf, state->x, 0, state->natoms);
1778 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
1780 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
1781 bUpdateDoLR, fr->f_twin, fcd,
1782 ekind, M, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
1783 wallcycle_stop(wcycle, ewcUPDATE);
1785 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1786 /* are the small terms in the shake_vir here due
1787 * to numerical errors, or are they important
1788 * physically? I'm thinking they are just errors, but not completely sure.
1789 * For now, will call without actually constraining, constr=NULL*/
1790 update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
1791 state, fr->bMolPBC, graph, f,
1792 &top->idef, tmp_vir, force_vir,
1793 cr, nrnb, wcycle, upd, NULL,
1794 bInitStep, FALSE, bCalcVir,
1797 if (!bOK && !bFFscan)
1799 gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
1802 if (fr->bSepDVDL && fplog && do_log)
1804 fprintf(fplog, sepdvdlformat, "Constraint dV/dl", 0.0, dvdl_constr);
1808 /* this factor or 2 correction is necessary
1809 because half of the constraint force is removed
1810 in the vv step, so we have to double it. See
1811 the Redmine issue #1255. It is not yet clear
1812 if the factor of 2 is exact, or just a very
1813 good approximation, and this will be
1814 investigated. The next step is to see if this
1815 can be done adding a dhdl contribution from the
1816 rattle step, but this is somewhat more
1817 complicated with the current code. Will be
1818 investigated, hopefully for 4.6.3. However,
1819 this current solution is much better than
1820 having it completely wrong.
1822 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1826 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1831 /* Need to unshift here */
1832 unshift_self(graph, state->box, state->x);
1835 GMX_BARRIER(cr->mpi_comm_mygroup);
1836 GMX_MPE_LOG(ev_update_finish);
1840 wallcycle_start(wcycle, ewcVSITECONSTR);
1843 shift_self(graph, state->box, state->x);
1845 construct_vsites(fplog, vsite, state->x, nrnb, ir->delta_t, state->v,
1846 top->idef.iparams, top->idef.il,
1847 fr->ePBC, fr->bMolPBC, graph, cr, state->box);
1851 unshift_self(graph, state->box, state->x);
1853 wallcycle_stop(wcycle, ewcVSITECONSTR);
1856 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1857 /* With Leap-Frog we can skip compute_globals at
1858 * non-communication steps, but we need to calculate
1859 * the kinetic energy one step before communication.
1861 if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)))
1863 if (ir->nstlist == -1 && bFirstIterate)
1865 gs.sig[eglsNABNSB] = nlh.nabnsb;
1867 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1868 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1870 bFirstIterate ? &gs : NULL,
1871 (step_rel % gs.nstms == 0) &&
1872 (multisim_nsteps < 0 || (step_rel < multisim_nsteps)),
1874 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1876 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
1877 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1878 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1879 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1880 | (iterate.bIterationActive ? CGLO_ITERATE : 0)
1881 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1884 if (ir->nstlist == -1 && bFirstIterate)
1886 nlh.nabnsb = gs.set[eglsNABNSB];
1887 gs.set[eglsNABNSB] = 0;
1890 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1891 /* ############# END CALC EKIN AND PRESSURE ################# */
1893 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1894 the virial that should probably be addressed eventually. state->veta has better properies,
1895 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1896 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1898 if (iterate.bIterationActive &&
1899 done_iterating(cr, fplog, step, &iterate, bFirstIterate,
1900 trace(shake_vir), &tracevir))
1904 bFirstIterate = FALSE;
1907 if (!bVV || bRerunMD)
1909 /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */
1910 sum_dhdl(enerd, state->lambda, ir->fepvals);
1912 update_box(fplog, step, ir, mdatoms, state, graph, f,
1913 ir->nstlist == -1 ? &nlh.scale_tot : NULL, pcoupl_mu, nrnb, wcycle, upd, bInitStep, FALSE);
1915 /* ################# END UPDATE STEP 2 ################# */
1916 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1918 /* The coordinates (x) were unshifted in update */
1919 if (bFFscan && (shellfc == NULL || bConverged))
1921 if (print_forcefield(fplog, enerd->term, mdatoms->homenr,
1923 &(top_global->mols), mdatoms->massT, pres))
1927 fprintf(stderr, "\n");
1933 /* We will not sum ekinh_old,
1934 * so signal that we still have to do it.
1936 bSumEkinhOld = TRUE;
1941 /* Only do GCT when the relaxation of shells (minimization) has converged,
1942 * otherwise we might be coupling to bogus energies.
1943 * In parallel we must always do this, because the other sims might
1947 /* Since this is called with the new coordinates state->x, I assume
1948 * we want the new box state->box too. / EL 20040121
1950 do_coupling(fplog, oenv, nfile, fnm, tcr, t, step, enerd->term, fr,
1952 mdatoms, &(top->idef), mu_aver,
1953 top_global->mols.nr, cr,
1954 state->box, total_vir, pres,
1955 mu_tot, state->x, f, bConverged);
1959 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1961 /* use the directly determined last velocity, not actually the averaged half steps */
1962 if (bTrotter && ir->eI == eiVV)
1964 enerd->term[F_EKIN] = last_ekin;
1966 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1970 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1974 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir, state, &MassQ);
1976 /* Check for excessively large energies */
1980 real etot_max = 1e200;
1982 real etot_max = 1e30;
1984 if (fabs(enerd->term[F_ETOT]) > etot_max)
1986 fprintf(stderr, "Energy too large (%g), giving up\n",
1987 enerd->term[F_ETOT]);
1990 /* ######### END PREPARING EDR OUTPUT ########### */
1992 /* Time for performance */
1993 if (((step % stepout) == 0) || bLastStep)
1995 runtime_upd_proc(runtime);
2001 gmx_bool do_dr, do_or;
2003 if (fplog && do_log && bDoExpanded)
2005 /* only needed if doing expanded ensemble */
2006 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
2007 &df_history, state->fep_state, ir->nstlog, step);
2009 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
2013 upd_mdebin(mdebin, bDoDHDL, TRUE,
2014 t, mdatoms->tmass, enerd, state,
2015 ir->fepvals, ir->expandedvals, lastbox,
2016 shake_vir, force_vir, total_vir, pres,
2017 ekind, mu_tot, constr);
2021 upd_mdebin_step(mdebin);
2024 do_dr = do_per_step(step, ir->nstdisreout);
2025 do_or = do_per_step(step, ir->nstorireout);
2027 print_ebin(outf->fp_ene, do_ene, do_dr, do_or, do_log ? fplog : NULL,
2029 eprNORMAL, bCompact, mdebin, fcd, groups, &(ir->opts));
2031 if (ir->ePull != epullNO)
2033 pull_print_output(ir->pull, step, t);
2036 if (do_per_step(step, ir->nstlog))
2038 if (fflush(fplog) != 0)
2040 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
2046 /* Have to do this part after outputting the logfile and the edr file */
2047 state->fep_state = lamnew;
2048 for (i = 0; i < efptNR; i++)
2050 state_global->lambda[i] = ir->fepvals->all_lambda[i][lamnew];
2053 /* Remaining runtime */
2054 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal()) && !bPMETuneRunning)
2058 fprintf(stderr, "\n");
2060 print_time(stderr, runtime, step, ir, cr);
2063 /* Replica exchange */
2065 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
2066 do_per_step(step, repl_ex_nst))
2068 bExchanged = replica_exchange(fplog, cr, repl_ex,
2069 state_global, enerd,
2072 if (bExchanged && DOMAINDECOMP(cr))
2074 dd_partition_system(fplog, step, cr, TRUE, 1,
2075 state_global, top_global, ir,
2076 state, &f, mdatoms, top, fr,
2077 vsite, shellfc, constr,
2078 nrnb, wcycle, FALSE);
2084 bStartingFromCpt = FALSE;
2086 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2087 /* With all integrators, except VV, we need to retain the pressure
2088 * at the current step for coupling at the next step.
2090 if ((state->flags & (1<<estPRES_PREV)) &&
2092 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
2094 /* Store the pressure in t_state for pressure coupling
2095 * at the next MD step.
2097 copy_mat(pres, state->pres_prev);
2100 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2102 if ( (membed != NULL) && (!bLastStep) )
2104 rescale_membed(step_rel, membed, state_global->x);
2111 /* read next frame from input trajectory */
2112 bNotLastFrame = read_next_frame(oenv, status, &rerun_fr);
2117 rerun_parallel_comm(cr, &rerun_fr, &bNotLastFrame);
2121 if (!bRerunMD || !rerun_fr.bStep)
2123 /* increase the MD step number */
2128 cycles = wallcycle_stop(wcycle, ewcSTEP);
2129 if (DOMAINDECOMP(cr) && wcycle)
2131 dd_cycles_add(cr->dd, cycles, ddCyclStep);
2134 if (bPMETuneRunning || bPMETuneTry)
2136 /* PME grid + cut-off optimization with GPUs or PME nodes */
2138 /* Count the total cycles over the last steps */
2139 cycles_pmes += cycles;
2141 /* We can only switch cut-off at NS steps */
2142 if (step % ir->nstlist == 0)
2144 /* PME grid + cut-off optimization with GPUs or PME nodes */
2147 if (DDMASTER(cr->dd))
2149 /* PME node load is too high, start tuning */
2150 bPMETuneRunning = (dd_pme_f_ratio(cr->dd) >= 1.05);
2152 dd_bcast(cr->dd, sizeof(gmx_bool), &bPMETuneRunning);
2154 if (bPMETuneRunning || step_rel > ir->nstlist*50)
2156 bPMETuneTry = FALSE;
2159 if (bPMETuneRunning)
2161 /* init_step might not be a multiple of nstlist,
2162 * but the first cycle is always skipped anyhow.
2165 pme_load_balance(pme_loadbal, cr,
2166 (bVerbose && MASTER(cr)) ? stderr : NULL,
2168 ir, state, cycles_pmes,
2169 fr->ic, fr->nbv, &fr->pmedata,
2172 /* Update constants in forcerec/inputrec to keep them in sync with fr->ic */
2173 fr->ewaldcoeff = fr->ic->ewaldcoeff;
2174 fr->rlist = fr->ic->rlist;
2175 fr->rlistlong = fr->ic->rlistlong;
2176 fr->rcoulomb = fr->ic->rcoulomb;
2177 fr->rvdw = fr->ic->rvdw;
2183 if (step_rel == wcycle_get_reset_counters(wcycle) ||
2184 gs.set[eglsRESETCOUNTERS] != 0)
2186 /* Reset all the counters related to performance over the run */
2187 reset_all_counters(fplog, cr, step, &step_rel, ir, wcycle, nrnb, runtime,
2188 fr->nbv != NULL && fr->nbv->bUseGPU ? fr->nbv->cu_nbv : NULL);
2189 wcycle_set_reset_counters(wcycle, -1);
2190 if (!(cr->duty & DUTY_PME))
2192 /* Tell our PME node to reset its counters */
2193 gmx_pme_send_resetcounters(cr, step);
2195 /* Correct max_hours for the elapsed time */
2196 max_hours -= run_time/(60.0*60.0);
2197 bResetCountersHalfMaxH = FALSE;
2198 gs.set[eglsRESETCOUNTERS] = 0;
2202 /* End of main MD loop */
2206 runtime_end(runtime);
2208 if (bRerunMD && MASTER(cr))
2213 if (!(cr->duty & DUTY_PME))
2215 /* Tell the PME only node to finish */
2216 gmx_pme_send_finish(cr);
2221 if (ir->nstcalcenergy > 0 && !bRerunMD)
2223 print_ebin(outf->fp_ene, FALSE, FALSE, FALSE, fplog, step, t,
2224 eprAVER, FALSE, mdebin, fcd, groups, &(ir->opts));
2232 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
2234 fprintf(fplog, "Average neighborlist lifetime: %.1f steps, std.dev.: %.1f steps\n", nlh.s1/nlh.nns, sqrt(nlh.s2/nlh.nns - sqr(nlh.s1/nlh.nns)));
2235 fprintf(fplog, "Average number of atoms that crossed the half buffer length: %.1f\n\n", nlh.ab/nlh.nns);
2238 if (pme_loadbal != NULL)
2240 pme_loadbal_done(pme_loadbal, cr, fplog,
2241 fr->nbv != NULL && fr->nbv->bUseGPU);
2244 if (shellfc && fplog)
2246 fprintf(fplog, "Fraction of iterations that converged: %.2f %%\n",
2247 (nconverged*100.0)/step_rel);
2248 fprintf(fplog, "Average number of force evaluations per MD step: %.2f\n\n",
2252 if (repl_ex_nst > 0 && MASTER(cr))
2254 print_replica_exchange_statistics(fplog, repl_ex);
2257 runtime->nsteps_done = step_rel;