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 t_state *state = NULL;
185 rvec *f_global = NULL;
188 gmx_enerdata_t *enerd;
190 gmx_global_stat_t gstat;
191 gmx_update_t upd = NULL;
192 t_graph *graph = NULL;
194 gmx_rng_t mcrng = NULL;
196 gmx_groups_t *groups;
197 gmx_ekindata_t *ekind, *ekind_save;
198 gmx_shellfc_t shellfc;
199 int count, nconverged = 0;
202 gmx_bool bIonize = FALSE;
203 gmx_bool bTCR = FALSE, bConverged = TRUE, bOK, bSumEkinhOld, bExchanged;
205 gmx_bool bResetCountersHalfMaxH = FALSE;
206 gmx_bool bVV, bIterativeCase, bFirstIterate, bTemp, bPres, bTrotter;
207 gmx_bool bUpdateDoLR;
208 real mu_aver = 0, dvdl_constr;
209 int a0, a1, gnx = 0, ii;
210 atom_id *grpindex = NULL;
212 t_coupl_rec *tcr = NULL;
213 rvec *xcopy = NULL, *vcopy = NULL, *cbuf = NULL;
214 matrix boxcopy = {{0}}, lastbox;
216 real fom, oldfom, veta_save, pcurr, scalevir, tracevir;
223 real saved_conserved_quantity = 0;
228 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
229 int handled_stop_condition = gmx_stop_cond_none; /* compare to get_stop_condition*/
230 gmx_iterate_t iterate;
231 gmx_large_int_t multisim_nsteps = -1; /* number of steps to do before first multisim
232 simulation stops. If equal to zero, don't
233 communicate any more between multisims.*/
234 /* PME load balancing data for GPU kernels */
235 pme_load_balancing_t pme_loadbal = NULL;
237 gmx_bool bPMETuneTry = FALSE, bPMETuneRunning = FALSE;
240 /* Temporary addition for FAHCORE checkpointing */
244 /* Check for special mdrun options */
245 bRerunMD = (Flags & MD_RERUN);
246 bIonize = (Flags & MD_IONIZE);
247 bFFscan = (Flags & MD_FFSCAN);
248 bAppend = (Flags & MD_APPENDFILES);
249 if (Flags & MD_RESETCOUNTERSHALFWAY)
253 /* Signal to reset the counters half the simulation steps. */
254 wcycle_set_reset_counters(wcycle, ir->nsteps/2);
256 /* Signal to reset the counters halfway the simulation time. */
257 bResetCountersHalfMaxH = (max_hours > 0);
260 /* md-vv uses averaged full step velocities for T-control
261 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
262 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
264 if (bVV) /* to store the initial velocities while computing virial */
266 snew(cbuf, top_global->natoms);
268 /* all the iteratative cases - only if there are constraints */
269 bIterativeCase = ((IR_NPH_TROTTER(ir) || IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
270 gmx_iterate_init(&iterate, FALSE); /* The default value of iterate->bIterationActive is set to
271 false in this step. The correct value, true or false,
272 is set at each step, as it depends on the frequency of temperature
273 and pressure control.*/
274 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir)));
278 /* Since we don't know if the frames read are related in any way,
279 * rebuild the neighborlist at every step.
282 ir->nstcalcenergy = 1;
286 check_ir_old_tpx_versions(cr, fplog, ir, top_global);
288 nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
289 bGStatEveryStep = (nstglobalcomm == 1);
291 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
294 "To reduce the energy communication with nstlist = -1\n"
295 "the neighbor list validity should not be checked at every step,\n"
296 "this means that exact integration is not guaranteed.\n"
297 "The neighbor list validity is checked after:\n"
298 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
299 "In most cases this will result in exact integration.\n"
300 "This reduces the energy communication by a factor of 2 to 3.\n"
301 "If you want less energy communication, set nstlist > 3.\n\n");
304 if (bRerunMD || bFFscan)
308 groups = &top_global->groups;
311 init_md(fplog, cr, ir, oenv, &t, &t0, state_global->lambda,
312 &(state_global->fep_state), lam0,
313 nrnb, top_global, &upd,
314 nfile, fnm, &outf, &mdebin,
315 force_vir, shake_vir, mu_tot, &bSimAnn, &vcm, state_global, Flags);
317 clear_mat(total_vir);
319 /* Energy terms and groups */
321 init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
323 if (DOMAINDECOMP(cr))
329 snew(f, top_global->natoms);
332 /* Kinetic energy data */
334 init_ekindata(fplog, top_global, &(ir->opts), ekind);
335 /* needed for iteration of constraints */
337 init_ekindata(fplog, top_global, &(ir->opts), ekind_save);
338 /* Copy the cos acceleration to the groups struct */
339 ekind->cosacc.cos_accel = ir->cos_accel;
341 gstat = global_stat_init(ir);
344 /* Check for polarizable models and flexible constraints */
345 shellfc = init_shell_flexcon(fplog,
346 top_global, n_flexible_constraints(constr),
347 (ir->bContinuation ||
348 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
349 NULL : state_global->x);
350 if (shellfc && ir->nstcalcenergy != 1)
352 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);
354 if (shellfc && ir->eI == eiNM)
356 /* Currently shells don't work with Normal Modes */
357 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");
360 if (vsite && ir->eI == eiNM)
362 /* Currently virtual sites don't work with Normal Modes */
363 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");
368 #ifdef GMX_THREAD_MPI
369 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
371 set_deform_reference_box(upd,
372 deform_init_init_step_tpx,
373 deform_init_box_tpx);
374 #ifdef GMX_THREAD_MPI
375 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
380 double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1);
381 if ((io > 2000) && MASTER(cr))
384 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
389 if (DOMAINDECOMP(cr))
391 top = dd_init_local_top(top_global);
394 dd_init_local_state(cr->dd, state_global, state);
396 if (DDMASTER(cr->dd) && ir->nstfout)
398 snew(f_global, state_global->natoms);
405 /* Initialize the particle decomposition and split the topology */
406 top = split_system(fplog, top_global, ir, cr);
408 pd_cg_range(cr, &fr->cg0, &fr->hcg);
409 pd_at_range(cr, &a0, &a1);
413 top = gmx_mtop_generate_local_top(top_global, ir);
416 a1 = top_global->natoms;
419 forcerec_set_excl_load(fr, top, cr);
421 state = partdec_init_local_state(cr, state_global);
424 atoms2md(top_global, ir, 0, NULL, a0, a1-a0, mdatoms);
428 set_vsite_top(vsite, top, mdatoms, cr);
431 if (ir->ePBC != epbcNONE && !fr->bMolPBC)
433 graph = mk_graph(fplog, &(top->idef), 0, top_global->natoms, FALSE, FALSE);
438 make_local_shells(cr, mdatoms, shellfc);
441 setup_bonded_threading(fr, &top->idef);
443 if (ir->pull && PAR(cr))
445 dd_make_local_pull_groups(NULL, ir->pull, mdatoms);
449 if (DOMAINDECOMP(cr))
451 /* Distribute the charge groups over the nodes from the master node */
452 dd_partition_system(fplog, ir->init_step, cr, TRUE, 1,
453 state_global, top_global, ir,
454 state, &f, mdatoms, top, fr,
455 vsite, shellfc, constr,
456 nrnb, wcycle, FALSE);
460 update_mdatoms(mdatoms, state->lambda[efptMASS]);
462 if (opt2bSet("-cpi", nfile, fnm))
464 bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi", nfile, fnm, cr), cr);
468 bStateFromCP = FALSE;
473 init_expanded_ensemble(bStateFromCP,ir,&mcrng,&state->dfhist);
480 /* Update mdebin with energy history if appending to output files */
481 if (Flags & MD_APPENDFILES)
483 restore_energyhistory_from_state(mdebin, &state_global->enerhist);
487 /* We might have read an energy history from checkpoint,
488 * free the allocated memory and reset the counts.
490 done_energyhistory(&state_global->enerhist);
491 init_energyhistory(&state_global->enerhist);
494 /* Set the initial energy history in state by updating once */
495 update_energyhistory(&state_global->enerhist, mdebin);
498 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG))
500 /* Set the random state if we read a checkpoint file */
501 set_stochd_state(upd, state);
504 if (state->flags & (1<<estMC_RNG))
506 set_mc_state(mcrng, state);
509 /* Initialize constraints */
512 if (!DOMAINDECOMP(cr))
514 set_constraints(constr, top, ir, mdatoms, cr);
518 /* Check whether we have to GCT stuff */
519 bTCR = ftp2bSet(efGCT, nfile, fnm);
524 fprintf(stderr, "Will do General Coupling Theory!\n");
526 gnx = top_global->mols.nr;
528 for (i = 0; (i < gnx); i++)
536 /* We need to be sure replica exchange can only occur
537 * when the energies are current */
538 check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
539 "repl_ex_nst", &repl_ex_nst);
540 /* This check needs to happen before inter-simulation
541 * signals are initialized, too */
543 if (repl_ex_nst > 0 && MASTER(cr))
545 repl_ex = init_replica_exchange(fplog, cr->ms, state_global, ir,
546 repl_ex_nst, repl_ex_nex, repl_ex_seed);
549 /* PME tuning is only supported with GPUs or PME nodes and not with rerun.
550 * With perturbed charges with soft-core we should not change the cut-off.
552 if ((Flags & MD_TUNEPME) &&
553 EEL_PME(fr->eeltype) &&
554 ( (fr->cutoff_scheme == ecutsVERLET && fr->nbv->bUseGPU) || !(cr->duty & DUTY_PME)) &&
555 !(ir->efep != efepNO && mdatoms->nChargePerturbed > 0 && ir->fepvals->bScCoul) &&
558 pme_loadbal_init(&pme_loadbal, ir, state->box, fr->ic, fr->pmedata);
560 if (cr->duty & DUTY_PME)
562 /* Start tuning right away, as we can't measure the load */
563 bPMETuneRunning = TRUE;
567 /* Separate PME nodes, we can measure the PP/PME load balance */
572 if (!ir->bContinuation && !bRerunMD)
574 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
576 /* Set the velocities of frozen particles to zero */
577 for (i = mdatoms->start; i < mdatoms->start+mdatoms->homenr; i++)
579 for (m = 0; m < DIM; m++)
581 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
591 /* Constrain the initial coordinates and velocities */
592 do_constrain_first(fplog, constr, ir, mdatoms, state, f,
593 graph, cr, nrnb, fr, top, shake_vir);
597 /* Construct the virtual sites for the initial configuration */
598 construct_vsites(fplog, vsite, state->x, nrnb, ir->delta_t, NULL,
599 top->idef.iparams, top->idef.il,
600 fr->ePBC, fr->bMolPBC, graph, cr, state->box);
606 /* set free energy calculation frequency as the minimum
607 greatest common denominator of nstdhdl, nstexpanded, and repl_ex_nst*/
608 nstfep = ir->fepvals->nstdhdl;
611 nstfep = gmx_greatest_common_divisor(ir->fepvals->nstdhdl,nstfep);
615 nstfep = gmx_greatest_common_divisor(repl_ex_nst,nstfep);
618 /* I'm assuming we need global communication the first time! MRS */
619 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
620 | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM : 0)
621 | (bVV ? CGLO_PRESSURE : 0)
622 | (bVV ? CGLO_CONSTRAINT : 0)
623 | (bRerunMD ? CGLO_RERUNMD : 0)
624 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN : 0));
626 bSumEkinhOld = FALSE;
627 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
628 NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
629 constr, NULL, FALSE, state->box,
630 top_global, &pcurr, top_global->natoms, &bSumEkinhOld, cglo_flags);
631 if (ir->eI == eiVVAK)
633 /* a second call to get the half step temperature initialized as well */
634 /* we do the same call as above, but turn the pressure off -- internally to
635 compute_globals, this is recognized as a velocity verlet half-step
636 kinetic energy calculation. This minimized excess variables, but
637 perhaps loses some logic?*/
639 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
640 NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
641 constr, NULL, FALSE, state->box,
642 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
643 cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
646 /* Calculate the initial half step temperature, and save the ekinh_old */
647 if (!(Flags & MD_STARTFROMCPT))
649 for (i = 0; (i < ir->opts.ngtc); i++)
651 copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
656 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
657 and there is no previous step */
660 /* if using an iterative algorithm, we need to create a working directory for the state. */
663 bufstate = init_bufstate(state);
667 snew(xcopy, state->natoms);
668 snew(vcopy, state->natoms);
669 copy_rvecn(state->x, xcopy, 0, state->natoms);
670 copy_rvecn(state->v, vcopy, 0, state->natoms);
671 copy_mat(state->box, boxcopy);
674 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
675 temperature control */
676 trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
680 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
683 "RMS relative constraint deviation after constraining: %.2e\n",
684 constr_rmsd(constr, FALSE));
686 if (EI_STATE_VELOCITY(ir->eI))
688 fprintf(fplog, "Initial temperature: %g K\n", enerd->term[F_TEMP]);
692 fprintf(stderr, "starting md rerun '%s', reading coordinates from"
693 " input trajectory '%s'\n\n",
694 *(top_global->name), opt2fn("-rerun", nfile, fnm));
697 fprintf(stderr, "Calculated time to finish depends on nsteps from "
698 "run input file,\nwhich may not correspond to the time "
699 "needed to process input trajectory.\n\n");
705 fprintf(stderr, "starting mdrun '%s'\n",
706 *(top_global->name));
709 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
713 sprintf(tbuf, "%s", "infinite");
715 if (ir->init_step > 0)
717 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
718 gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
719 gmx_step_str(ir->init_step, sbuf2),
720 ir->init_step*ir->delta_t);
724 fprintf(stderr, "%s steps, %s ps.\n",
725 gmx_step_str(ir->nsteps, sbuf), tbuf);
728 fprintf(fplog, "\n");
731 print_start(fplog, cr, runtime, "mdrun");
732 runtime_start(runtime);
733 wallcycle_start(wcycle, ewcRUN);
735 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
737 chkpt_ret = fcCheckPointParallel( cr->nodeid,
741 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
746 /***********************************************************
750 ************************************************************/
752 /* if rerunMD then read coordinates and velocities from input trajectory */
755 if (getenv("GMX_FORCE_UPDATE"))
763 bNotLastFrame = read_first_frame(oenv, &status,
764 opt2fn("-rerun", nfile, fnm),
765 &rerun_fr, TRX_NEED_X | TRX_READ_V);
766 if (rerun_fr.natoms != top_global->natoms)
769 "Number of atoms in trajectory (%d) does not match the "
770 "run input file (%d)\n",
771 rerun_fr.natoms, top_global->natoms);
773 if (ir->ePBC != epbcNONE)
777 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);
779 if (max_cutoff2(ir->ePBC, rerun_fr.box) < sqr(fr->rlistlong))
781 gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr.step, rerun_fr.time);
788 rerun_parallel_comm(cr, &rerun_fr, &bNotLastFrame);
791 if (ir->ePBC != epbcNONE)
793 /* Set the shift vectors.
794 * Necessary here when have a static box different from the tpr box.
796 calc_shifts(rerun_fr.box, fr->shift_vec);
800 /* loop over MD steps or if rerunMD to end of input trajectory */
802 /* Skip the first Nose-Hoover integration when we get the state from tpx */
803 bStateFromTPX = !bStateFromCP;
804 bInitStep = bFirstStep && (bStateFromTPX || bVV);
805 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
807 bSumEkinhOld = FALSE;
810 init_global_signals(&gs, cr, ir, repl_ex_nst);
812 step = ir->init_step;
815 if (ir->nstlist == -1)
817 init_nlistheuristics(&nlh, bGStatEveryStep, step);
820 if (MULTISIM(cr) && (repl_ex_nst <= 0 ))
822 /* check how many steps are left in other sims */
823 multisim_nsteps = get_multisim_nsteps(cr, ir->nsteps);
827 /* and stop now if we should */
828 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
829 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
830 while (!bLastStep || (bRerunMD && bNotLastFrame))
833 wallcycle_start(wcycle, ewcSTEP);
835 GMX_MPE_LOG(ev_timestep1);
841 step = rerun_fr.step;
842 step_rel = step - ir->init_step;
855 bLastStep = (step_rel == ir->nsteps);
856 t = t0 + step*ir->delta_t;
859 if (ir->efep != efepNO || ir->bSimTemp)
861 /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
862 requiring different logic. */
864 set_current_lambdas(step, ir->fepvals, bRerunMD, &rerun_fr, state_global, state, lam0);
865 bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
866 bDoFEP = (do_per_step(step, nstfep) && (ir->efep != efepNO));
867 bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
868 && (ir->bExpanded) && (step > 0) && (!bStartingFromCpt));
873 update_annealing_target_temp(&(ir->opts), t);
878 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
880 for (i = 0; i < state_global->natoms; i++)
882 copy_rvec(rerun_fr.x[i], state_global->x[i]);
886 for (i = 0; i < state_global->natoms; i++)
888 copy_rvec(rerun_fr.v[i], state_global->v[i]);
893 for (i = 0; i < state_global->natoms; i++)
895 clear_rvec(state_global->v[i]);
899 fprintf(stderr, "\nWARNING: Some frames do not contain velocities.\n"
900 " Ekin, temperature and pressure are incorrect,\n"
901 " the virial will be incorrect when constraints are present.\n"
903 bRerunWarnNoV = FALSE;
907 copy_mat(rerun_fr.box, state_global->box);
908 copy_mat(state_global->box, state->box);
910 if (vsite && (Flags & MD_RERUN_VSITE))
912 if (DOMAINDECOMP(cr))
914 gmx_fatal(FARGS, "Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
918 /* Following is necessary because the graph may get out of sync
919 * with the coordinates if we only have every N'th coordinate set
921 mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
922 shift_self(graph, state->box, state->x);
924 construct_vsites(fplog, vsite, state->x, nrnb, ir->delta_t, state->v,
925 top->idef.iparams, top->idef.il,
926 fr->ePBC, fr->bMolPBC, graph, cr, state->box);
929 unshift_self(graph, state->box, state->x);
934 /* Stop Center of Mass motion */
935 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
937 /* Copy back starting coordinates in case we're doing a forcefield scan */
940 for (ii = 0; (ii < state->natoms); ii++)
942 copy_rvec(xcopy[ii], state->x[ii]);
943 copy_rvec(vcopy[ii], state->v[ii]);
945 copy_mat(boxcopy, state->box);
950 /* for rerun MD always do Neighbour Searching */
951 bNS = (bFirstStep || ir->nstlist != 0);
956 /* Determine whether or not to do Neighbour Searching and LR */
957 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
959 bNS = (bFirstStep || bExchanged || bNStList || bDoFEP ||
960 (ir->nstlist == -1 && nlh.nabnsb > 0));
962 if (bNS && ir->nstlist == -1)
964 set_nlistheuristics(&nlh, bFirstStep || bExchanged || bDoFEP, step);
968 /* check whether we should stop because another simulation has
972 if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
973 (multisim_nsteps != ir->nsteps) )
980 "Stopping simulation %d because another one has finished\n",
984 gs.sig[eglsCHKPT] = 1;
989 /* < 0 means stop at next step, > 0 means stop at next NS step */
990 if ( (gs.set[eglsSTOPCOND] < 0) ||
991 ( (gs.set[eglsSTOPCOND] > 0) && (bNStList || ir->nstlist == 0) ) )
996 /* Determine whether or not to update the Born radii if doing GB */
997 bBornRadii = bFirstStep;
998 if (ir->implicit_solvent && (step % ir->nstgbradii == 0))
1003 do_log = do_per_step(step, ir->nstlog) || bFirstStep || bLastStep;
1004 do_verbose = bVerbose &&
1005 (step % stepout == 0 || bFirstStep || bLastStep);
1007 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
1011 bMasterState = TRUE;
1015 bMasterState = FALSE;
1016 /* Correct the new box if it is too skewed */
1017 if (DYNAMIC_BOX(*ir))
1019 if (correct_box(fplog, step, state->box, graph))
1021 bMasterState = TRUE;
1024 if (DOMAINDECOMP(cr) && bMasterState)
1026 dd_collect_state(cr->dd, state, state_global);
1030 if (DOMAINDECOMP(cr))
1032 /* Repartition the domain decomposition */
1033 wallcycle_start(wcycle, ewcDOMDEC);
1034 dd_partition_system(fplog, step, cr,
1035 bMasterState, nstglobalcomm,
1036 state_global, top_global, ir,
1037 state, &f, mdatoms, top, fr,
1038 vsite, shellfc, constr,
1040 do_verbose && !bPMETuneRunning);
1041 wallcycle_stop(wcycle, ewcDOMDEC);
1042 /* If using an iterative integrator, reallocate space to match the decomposition */
1046 if (MASTER(cr) && do_log && !bFFscan)
1048 print_ebin_header(fplog, step, t, state->lambda[efptFEP]); /* can we improve the information printed here? */
1051 if (ir->efep != efepNO)
1053 update_mdatoms(mdatoms, state->lambda[efptMASS]);
1056 if ((bRerunMD && rerun_fr.bV) || bExchanged)
1059 /* We need the kinetic energy at minus the half step for determining
1060 * the full step kinetic energy and possibly for T-coupling.*/
1061 /* This may not be quite working correctly yet . . . . */
1062 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1063 wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
1064 constr, NULL, FALSE, state->box,
1065 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1066 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1068 clear_mat(force_vir);
1070 /* Ionize the atoms if necessary */
1073 ionize(fplog, oenv, mdatoms, top_global, t, ir, state->x, state->v,
1074 mdatoms->start, mdatoms->start+mdatoms->homenr, state->box, cr);
1077 /* Update force field in ffscan program */
1080 if (update_forcefield(fplog,
1082 mdatoms->nr, state->x, state->box))
1090 GMX_MPE_LOG(ev_timestep2);
1092 /* We write a checkpoint at this MD step when:
1093 * either at an NS step when we signalled through gs,
1094 * or at the last step (but not when we do not want confout),
1095 * but never at the first step or with rerun.
1097 bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
1098 (bLastStep && (Flags & MD_CONFOUT))) &&
1099 step > ir->init_step && !bRerunMD);
1102 gs.set[eglsCHKPT] = 0;
1105 /* Determine the energy and pressure:
1106 * at nstcalcenergy steps and at energy output steps (set below).
1108 if (EI_VV(ir->eI) && (!bInitStep))
1110 /* for vv, the first half of the integration actually corresponds
1111 to the previous step. bCalcEner is only required to be evaluated on the 'next' step,
1112 but the virial needs to be calculated on both the current step and the 'next' step. Future
1113 reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
1115 bCalcEner = do_per_step(step-1, ir->nstcalcenergy);
1116 bCalcVir = bCalcEner ||
1117 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
1121 bCalcEner = do_per_step(step, ir->nstcalcenergy);
1122 bCalcVir = bCalcEner ||
1123 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
1126 /* Do we need global communication ? */
1127 bGStat = (bCalcVir || bCalcEner || bStopCM ||
1128 do_per_step(step, nstglobalcomm) || (bVV && IR_NVT_TROTTER(ir) && do_per_step(step-1, nstglobalcomm)) ||
1129 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1131 do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
1133 if (do_ene || do_log)
1140 /* these CGLO_ options remain the same throughout the iteration */
1141 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1142 (bGStat ? CGLO_GSTAT : 0)
1145 force_flags = (GMX_FORCE_STATECHANGED |
1146 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1147 GMX_FORCE_ALLFORCES |
1149 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
1150 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
1151 (bDoFEP ? GMX_FORCE_DHDL : 0)
1156 if (do_per_step(step, ir->nstcalclr))
1158 force_flags |= GMX_FORCE_DO_LR;
1164 /* Now is the time to relax the shells */
1165 count = relax_shell_flexcon(fplog, cr, bVerbose, bFFscan ? step+1 : step,
1166 ir, bNS, force_flags,
1167 bStopCM, top, top_global,
1169 state, f, force_vir, mdatoms,
1170 nrnb, wcycle, graph, groups,
1171 shellfc, fr, bBornRadii, t, mu_tot,
1172 state->natoms, &bConverged, vsite,
1183 /* The coordinates (x) are shifted (to get whole molecules)
1185 * This is parallellized as well, and does communication too.
1186 * Check comments in sim_util.c
1188 do_force(fplog, cr, ir, step, nrnb, wcycle, top, top_global, groups,
1189 state->box, state->x, &state->hist,
1190 f, force_vir, mdatoms, enerd, fcd,
1191 state->lambda, graph,
1192 fr, vsite, mu_tot, t, outf->fp_field, ed, bBornRadii,
1193 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1196 GMX_BARRIER(cr->mpi_comm_mygroup);
1200 mu_aver = calc_mu_aver(cr, state->x, mdatoms->chargeA,
1201 mu_tot, &top_global->mols, mdatoms, gnx, grpindex);
1204 if (bTCR && bFirstStep)
1206 tcr = init_coupling(fplog, nfile, fnm, cr, fr, mdatoms, &(top->idef));
1207 fprintf(fplog, "Done init_coupling\n");
1211 if (bVV && !bStartingFromCpt && !bRerunMD)
1212 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1214 if (ir->eI == eiVV && bInitStep)
1216 /* if using velocity verlet with full time step Ekin,
1217 * take the first half step only to compute the
1218 * virial for the first step. From there,
1219 * revert back to the initial coordinates
1220 * so that the input is actually the initial step.
1222 copy_rvecn(state->v, cbuf, 0, state->natoms); /* should make this better for parallelizing? */
1226 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1227 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
1230 /* If we are using twin-range interactions where the long-range component
1231 * is only evaluated every nstcalclr>1 steps, we should do a special update
1232 * step to combine the long-range forces on these steps.
1233 * For nstcalclr=1 this is not done, since the forces would have been added
1234 * directly to the short-range forces already.
1236 * TODO Remove various aspects of VV+twin-range in master
1237 * branch, because VV integrators did not ever support
1238 * twin-range multiple time stepping with constraints.
1240 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
1242 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC,
1243 f, bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
1244 ekind, M, wcycle, upd, bInitStep, etrtVELOCITY1,
1245 cr, nrnb, constr, &top->idef);
1247 if (bIterativeCase && do_per_step(step-1, ir->nstpcouple) && !bInitStep)
1249 gmx_iterate_init(&iterate, TRUE);
1251 /* for iterations, we save these vectors, as we will be self-consistently iterating
1254 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1256 /* save the state */
1257 if (iterate.bIterationActive)
1259 copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts));
1262 bFirstIterate = TRUE;
1263 while (bFirstIterate || iterate.bIterationActive)
1265 if (iterate.bIterationActive)
1267 copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts));
1268 if (bFirstIterate && bTrotter)
1270 /* The first time through, we need a decent first estimate
1271 of veta(t+dt) to compute the constraints. Do
1272 this by computing the box volume part of the
1273 trotter integration at this time. Nothing else
1274 should be changed by this routine here. If
1275 !(first time), we start with the previous value
1278 veta_save = state->veta;
1279 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ0);
1280 vetanew = state->veta;
1281 state->veta = veta_save;
1286 if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
1288 update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
1289 state, fr->bMolPBC, graph, f,
1290 &top->idef, shake_vir, NULL,
1291 cr, nrnb, wcycle, upd, constr,
1292 bInitStep, TRUE, bCalcVir, vetanew);
1294 if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
1296 /* Correct the virial for multiple time stepping */
1297 m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
1300 if (!bOK && !bFFscan)
1302 gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
1308 /* Need to unshift here if a do_force has been
1309 called in the previous step */
1310 unshift_self(graph, state->box, state->x);
1313 /* if VV, compute the pressure and constraints */
1314 /* For VV2, we strictly only need this if using pressure
1315 * control, but we really would like to have accurate pressures
1317 * Think about ways around this in the future?
1318 * For now, keep this choice in comments.
1320 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1321 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1323 bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
1324 if (bCalcEner && ir->eI == eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/
1326 bSumEkinhOld = TRUE;
1328 /* for vv, the first half of the integration actually corresponds to the previous step.
1329 So we need information from the last step in the first half of the integration */
1330 if (bGStat || do_per_step(step-1, nstglobalcomm))
1332 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1333 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1334 constr, NULL, FALSE, state->box,
1335 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1338 | (bTemp ? CGLO_TEMPERATURE : 0)
1339 | (bPres ? CGLO_PRESSURE : 0)
1340 | (bPres ? CGLO_CONSTRAINT : 0)
1341 | ((iterate.bIterationActive) ? CGLO_ITERATE : 0)
1342 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1345 /* explanation of above:
1346 a) We compute Ekin at the full time step
1347 if 1) we are using the AveVel Ekin, and it's not the
1348 initial step, or 2) if we are using AveEkin, but need the full
1349 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1350 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1351 EkinAveVel because it's needed for the pressure */
1353 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1358 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
1359 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
1366 /* We need the kinetic energy at minus the half step for determining
1367 * the full step kinetic energy and possibly for T-coupling.*/
1368 /* This may not be quite working correctly yet . . . . */
1369 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1370 wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
1371 constr, NULL, FALSE, state->box,
1372 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1373 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1378 if (iterate.bIterationActive &&
1379 done_iterating(cr, fplog, step, &iterate, bFirstIterate,
1380 state->veta, &vetanew))
1384 bFirstIterate = FALSE;
1387 if (bTrotter && !bInitStep)
1389 copy_mat(shake_vir, state->svir_prev);
1390 copy_mat(force_vir, state->fvir_prev);
1391 if (IR_NVT_TROTTER(ir) && ir->eI == eiVV)
1393 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1394 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE, FALSE);
1395 enerd->term[F_EKIN] = trace(ekind->ekin);
1398 /* if it's the initial step, we performed this first step just to get the constraint virial */
1399 if (bInitStep && ir->eI == eiVV)
1401 copy_rvecn(cbuf, state->v, 0, state->natoms);
1404 GMX_MPE_LOG(ev_timestep1);
1407 /* MRS -- now done iterating -- compute the conserved quantity */
1410 saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ);
1413 last_ekin = enerd->term[F_EKIN];
1415 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1417 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1419 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1422 sum_dhdl(enerd, state->lambda, ir->fepvals);
1426 /* ######## END FIRST UPDATE STEP ############## */
1427 /* ######## If doing VV, we now have v(dt) ###### */
1430 /* perform extended ensemble sampling in lambda - we don't
1431 actually move to the new state before outputting
1432 statistics, but if performing simulated tempering, we
1433 do update the velocities and the tau_t. */
1435 lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, &state->dfhist, step, mcrng, state->v, mdatoms);
1436 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1437 copy_df_history(&state_global->dfhist,&state->dfhist);
1439 /* ################## START TRAJECTORY OUTPUT ################# */
1441 /* Now we have the energies and forces corresponding to the
1442 * coordinates at time t. We must output all of this before
1444 * for RerunMD t is read from input trajectory
1446 GMX_MPE_LOG(ev_output_start);
1449 if (do_per_step(step, ir->nstxout))
1451 mdof_flags |= MDOF_X;
1453 if (do_per_step(step, ir->nstvout))
1455 mdof_flags |= MDOF_V;
1457 if (do_per_step(step, ir->nstfout))
1459 mdof_flags |= MDOF_F;
1461 if (do_per_step(step, ir->nstxtcout))
1463 mdof_flags |= MDOF_XTC;
1467 mdof_flags |= MDOF_CPT;
1471 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1474 /* Enforce writing positions and velocities at end of run */
1475 mdof_flags |= (MDOF_X | MDOF_V);
1481 fcReportProgress( ir->nsteps, step );
1484 #if defined(__native_client__)
1485 fcCheckin(MASTER(cr));
1488 /* sync bCPT and fc record-keeping */
1489 if (bCPT && MASTER(cr))
1491 fcRequestCheckPoint();
1495 if (mdof_flags != 0)
1497 wallcycle_start(wcycle, ewcTRAJ);
1500 if (state->flags & (1<<estLD_RNG))
1502 get_stochd_state(upd, state);
1504 if (state->flags & (1<<estMC_RNG))
1506 get_mc_state(mcrng, state);
1512 state_global->ekinstate.bUpToDate = FALSE;
1516 update_ekinstate(&state_global->ekinstate, ekind);
1517 state_global->ekinstate.bUpToDate = TRUE;
1519 update_energyhistory(&state_global->enerhist, mdebin);
1522 write_traj(fplog, cr, outf, mdof_flags, top_global,
1523 step, t, state, state_global, f, f_global, &n_xtc, &x_xtc);
1530 if (bLastStep && step_rel == ir->nsteps &&
1531 (Flags & MD_CONFOUT) && MASTER(cr) &&
1532 !bRerunMD && !bFFscan)
1534 /* x and v have been collected in write_traj,
1535 * because a checkpoint file will always be written
1538 fprintf(stderr, "\nWriting final coordinates.\n");
1541 /* Make molecules whole only for confout writing */
1542 do_pbc_mtop(fplog, ir->ePBC, state->box, top_global, state_global->x);
1544 write_sto_conf_mtop(ftp2fn(efSTO, nfile, fnm),
1545 *top_global->name, top_global,
1546 state_global->x, state_global->v,
1547 ir->ePBC, state->box);
1550 wallcycle_stop(wcycle, ewcTRAJ);
1552 GMX_MPE_LOG(ev_output_finish);
1554 /* kludge -- virial is lost with restart for NPT control. Must restart */
1555 if (bStartingFromCpt && bVV)
1557 copy_mat(state->svir_prev, shake_vir);
1558 copy_mat(state->fvir_prev, force_vir);
1560 /* ################## END TRAJECTORY OUTPUT ################ */
1562 /* Determine the wallclock run time up till now */
1563 run_time = gmx_gettime() - (double)runtime->real;
1565 /* Check whether everything is still allright */
1566 if (((int)gmx_get_stop_condition() > handled_stop_condition)
1567 #ifdef GMX_THREAD_MPI
1572 /* this is just make gs.sig compatible with the hack
1573 of sending signals around by MPI_Reduce with together with
1575 if (gmx_get_stop_condition() == gmx_stop_cond_next_ns)
1577 gs.sig[eglsSTOPCOND] = 1;
1579 if (gmx_get_stop_condition() == gmx_stop_cond_next)
1581 gs.sig[eglsSTOPCOND] = -1;
1583 /* < 0 means stop at next step, > 0 means stop at next NS step */
1587 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1588 gmx_get_signal_name(),
1589 gs.sig[eglsSTOPCOND] == 1 ? "NS " : "");
1593 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1594 gmx_get_signal_name(),
1595 gs.sig[eglsSTOPCOND] == 1 ? "NS " : "");
1597 handled_stop_condition = (int)gmx_get_stop_condition();
1599 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1600 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
1601 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
1603 /* Signal to terminate the run */
1604 gs.sig[eglsSTOPCOND] = 1;
1607 fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
1609 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
1612 if (bResetCountersHalfMaxH && MASTER(cr) &&
1613 run_time > max_hours*60.0*60.0*0.495)
1615 gs.sig[eglsRESETCOUNTERS] = 1;
1618 if (ir->nstlist == -1 && !bRerunMD)
1620 /* When bGStatEveryStep=FALSE, global_stat is only called
1621 * when we check the atom displacements, not at NS steps.
1622 * This means that also the bonded interaction count check is not
1623 * performed immediately after NS. Therefore a few MD steps could
1624 * be performed with missing interactions.
1625 * But wrong energies are never written to file,
1626 * since energies are only written after global_stat
1629 if (step >= nlh.step_nscheck)
1631 nlh.nabnsb = natoms_beyond_ns_buffer(ir, fr, &top->cgs,
1632 nlh.scale_tot, state->x);
1636 /* This is not necessarily true,
1637 * but step_nscheck is determined quite conservatively.
1643 /* In parallel we only have to check for checkpointing in steps
1644 * where we do global communication,
1645 * otherwise the other nodes don't know.
1647 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1650 run_time >= nchkpt*cpt_period*60.0)) &&
1651 gs.set[eglsCHKPT] == 0)
1653 gs.sig[eglsCHKPT] = 1;
1656 /* at the start of step, randomize or scale the velocities (trotter done elsewhere) */
1661 update_tcouple(fplog, step, ir, state, ekind, wcycle, upd, &MassQ, mdatoms);
1663 if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1665 gmx_bool bIfRandomize;
1666 bIfRandomize = update_randomize_velocities(ir, step, mdatoms, state, upd, &top->idef, constr, DOMAINDECOMP(cr));
1667 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1668 if (constr && bIfRandomize)
1670 update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
1671 state, fr->bMolPBC, graph, f,
1672 &top->idef, tmp_vir, NULL,
1673 cr, nrnb, wcycle, upd, constr,
1674 bInitStep, TRUE, bCalcVir, vetanew);
1679 if (bIterativeCase && do_per_step(step, ir->nstpcouple))
1681 gmx_iterate_init(&iterate, TRUE);
1682 /* for iterations, we save these vectors, as we will be redoing the calculations */
1683 copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts));
1686 bFirstIterate = TRUE;
1687 while (bFirstIterate || iterate.bIterationActive)
1689 /* We now restore these vectors to redo the calculation with improved extended variables */
1690 if (iterate.bIterationActive)
1692 copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts));
1695 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1696 so scroll down for that logic */
1698 /* ######### START SECOND UPDATE STEP ################# */
1699 GMX_MPE_LOG(ev_update_start);
1700 /* Box is changed in update() when we do pressure coupling,
1701 * but we should still use the old box for energy corrections and when
1702 * writing it to the energy file, so it matches the trajectory files for
1703 * the same timestep above. Make a copy in a separate array.
1705 copy_mat(state->box, lastbox);
1710 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
1712 wallcycle_start(wcycle, ewcUPDATE);
1713 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1716 if (iterate.bIterationActive)
1724 /* we use a new value of scalevir to converge the iterations faster */
1725 scalevir = tracevir/trace(shake_vir);
1727 msmul(shake_vir, scalevir, shake_vir);
1728 m_add(force_vir, shake_vir, total_vir);
1729 clear_mat(shake_vir);
1731 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1732 /* We can only do Berendsen coupling after we have summed
1733 * the kinetic energy or virial. Since the happens
1734 * in global_state after update, we should only do it at
1735 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1740 update_tcouple(fplog, step, ir, state, ekind, wcycle, upd, &MassQ, mdatoms);
1741 update_pcouple(fplog, step, ir, state, pcoupl_mu, M, wcycle,
1747 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
1749 /* velocity half-step update */
1750 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
1751 bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
1752 ekind, M, wcycle, upd, FALSE, etrtVELOCITY2,
1753 cr, nrnb, constr, &top->idef);
1756 /* Above, initialize just copies ekinh into ekin,
1757 * it doesn't copy position (for VV),
1758 * and entire integrator for MD.
1761 if (ir->eI == eiVVAK)
1763 copy_rvecn(state->x, cbuf, 0, state->natoms);
1765 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
1767 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
1768 bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
1769 ekind, M, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
1770 wallcycle_stop(wcycle, ewcUPDATE);
1772 update_constraints(fplog, step, &dvdl_constr, ir, ekind, mdatoms, state,
1773 fr->bMolPBC, graph, f,
1774 &top->idef, shake_vir, force_vir,
1775 cr, nrnb, wcycle, upd, constr,
1776 bInitStep, FALSE, bCalcVir, state->veta);
1778 if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
1780 /* Correct the virial for multiple time stepping */
1781 m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
1784 if (ir->eI == eiVVAK)
1786 /* erase F_EKIN and F_TEMP here? */
1787 /* just compute the kinetic energy at the half step to perform a trotter step */
1788 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1789 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1790 constr, NULL, FALSE, lastbox,
1791 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1792 cglo_flags | CGLO_TEMPERATURE
1794 wallcycle_start(wcycle, ewcUPDATE);
1795 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1796 /* now we know the scaling, we can compute the positions again again */
1797 copy_rvecn(cbuf, state->x, 0, state->natoms);
1799 bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
1801 update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
1802 bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
1803 ekind, M, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
1804 wallcycle_stop(wcycle, ewcUPDATE);
1806 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1807 /* are the small terms in the shake_vir here due
1808 * to numerical errors, or are they important
1809 * physically? I'm thinking they are just errors, but not completely sure.
1810 * For now, will call without actually constraining, constr=NULL*/
1811 update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
1812 state, fr->bMolPBC, graph, f,
1813 &top->idef, tmp_vir, force_vir,
1814 cr, nrnb, wcycle, upd, NULL,
1815 bInitStep, FALSE, bCalcVir,
1818 if (!bOK && !bFFscan)
1820 gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
1823 if (fr->bSepDVDL && fplog && do_log)
1825 fprintf(fplog, sepdvdlformat, "Constraint dV/dl", 0.0, dvdl_constr);
1829 /* this factor or 2 correction is necessary
1830 because half of the constraint force is removed
1831 in the vv step, so we have to double it. See
1832 the Redmine issue #1255. It is not yet clear
1833 if the factor of 2 is exact, or just a very
1834 good approximation, and this will be
1835 investigated. The next step is to see if this
1836 can be done adding a dhdl contribution from the
1837 rattle step, but this is somewhat more
1838 complicated with the current code. Will be
1839 investigated, hopefully for 4.6.3. However,
1840 this current solution is much better than
1841 having it completely wrong.
1843 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1847 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1852 /* Need to unshift here */
1853 unshift_self(graph, state->box, state->x);
1856 GMX_BARRIER(cr->mpi_comm_mygroup);
1857 GMX_MPE_LOG(ev_update_finish);
1861 wallcycle_start(wcycle, ewcVSITECONSTR);
1864 shift_self(graph, state->box, state->x);
1866 construct_vsites(fplog, vsite, state->x, nrnb, ir->delta_t, state->v,
1867 top->idef.iparams, top->idef.il,
1868 fr->ePBC, fr->bMolPBC, graph, cr, state->box);
1872 unshift_self(graph, state->box, state->x);
1874 wallcycle_stop(wcycle, ewcVSITECONSTR);
1877 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1878 /* With Leap-Frog we can skip compute_globals at
1879 * non-communication steps, but we need to calculate
1880 * the kinetic energy one step before communication.
1882 if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)))
1884 if (ir->nstlist == -1 && bFirstIterate)
1886 gs.sig[eglsNABNSB] = nlh.nabnsb;
1888 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm,
1889 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1891 bFirstIterate ? &gs : NULL,
1892 (step_rel % gs.nstms == 0) &&
1893 (multisim_nsteps < 0 || (step_rel < multisim_nsteps)),
1895 top_global, &pcurr, top_global->natoms, &bSumEkinhOld,
1897 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
1898 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1899 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1900 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1901 | (iterate.bIterationActive ? CGLO_ITERATE : 0)
1902 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1905 if (ir->nstlist == -1 && bFirstIterate)
1907 nlh.nabnsb = gs.set[eglsNABNSB];
1908 gs.set[eglsNABNSB] = 0;
1911 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1912 /* ############# END CALC EKIN AND PRESSURE ################# */
1914 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1915 the virial that should probably be addressed eventually. state->veta has better properies,
1916 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1917 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1919 if (iterate.bIterationActive &&
1920 done_iterating(cr, fplog, step, &iterate, bFirstIterate,
1921 trace(shake_vir), &tracevir))
1925 bFirstIterate = FALSE;
1928 if (!bVV || bRerunMD)
1930 /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */
1931 sum_dhdl(enerd, state->lambda, ir->fepvals);
1933 update_box(fplog, step, ir, mdatoms, state, graph, f,
1934 ir->nstlist == -1 ? &nlh.scale_tot : NULL, pcoupl_mu, nrnb, wcycle, upd, bInitStep, FALSE);
1936 /* ################# END UPDATE STEP 2 ################# */
1937 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1939 /* The coordinates (x) were unshifted in update */
1940 if (bFFscan && (shellfc == NULL || bConverged))
1942 if (print_forcefield(fplog, enerd->term, mdatoms->homenr,
1944 &(top_global->mols), mdatoms->massT, pres))
1948 fprintf(stderr, "\n");
1954 /* We will not sum ekinh_old,
1955 * so signal that we still have to do it.
1957 bSumEkinhOld = TRUE;
1962 /* Only do GCT when the relaxation of shells (minimization) has converged,
1963 * otherwise we might be coupling to bogus energies.
1964 * In parallel we must always do this, because the other sims might
1968 /* Since this is called with the new coordinates state->x, I assume
1969 * we want the new box state->box too. / EL 20040121
1971 do_coupling(fplog, oenv, nfile, fnm, tcr, t, step, enerd->term, fr,
1973 mdatoms, &(top->idef), mu_aver,
1974 top_global->mols.nr, cr,
1975 state->box, total_vir, pres,
1976 mu_tot, state->x, f, bConverged);
1980 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1982 /* use the directly determined last velocity, not actually the averaged half steps */
1983 if (bTrotter && ir->eI == eiVV)
1985 enerd->term[F_EKIN] = last_ekin;
1987 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1991 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1995 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir, state, &MassQ);
1997 /* Check for excessively large energies */
2001 real etot_max = 1e200;
2003 real etot_max = 1e30;
2005 if (fabs(enerd->term[F_ETOT]) > etot_max)
2007 fprintf(stderr, "Energy too large (%g), giving up\n",
2008 enerd->term[F_ETOT]);
2011 /* ######### END PREPARING EDR OUTPUT ########### */
2013 /* Time for performance */
2014 if (((step % stepout) == 0) || bLastStep)
2016 runtime_upd_proc(runtime);
2022 gmx_bool do_dr, do_or;
2024 if (fplog && do_log && bDoExpanded)
2026 /* only needed if doing expanded ensemble */
2027 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
2028 &state_global->dfhist, state->fep_state, ir->nstlog, step);
2030 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
2034 upd_mdebin(mdebin, bDoDHDL, TRUE,
2035 t, mdatoms->tmass, enerd, state,
2036 ir->fepvals, ir->expandedvals, lastbox,
2037 shake_vir, force_vir, total_vir, pres,
2038 ekind, mu_tot, constr);
2042 upd_mdebin_step(mdebin);
2045 do_dr = do_per_step(step, ir->nstdisreout);
2046 do_or = do_per_step(step, ir->nstorireout);
2048 print_ebin(outf->fp_ene, do_ene, do_dr, do_or, do_log ? fplog : NULL,
2050 eprNORMAL, bCompact, mdebin, fcd, groups, &(ir->opts));
2052 if (ir->ePull != epullNO)
2054 pull_print_output(ir->pull, step, t);
2057 if (do_per_step(step, ir->nstlog))
2059 if (fflush(fplog) != 0)
2061 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
2067 /* Have to do this part _after_ outputting the logfile and the edr file */
2068 /* Gets written into the state at the beginning of next loop*/
2069 state->fep_state = lamnew;
2072 /* Remaining runtime */
2073 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal()) && !bPMETuneRunning)
2077 fprintf(stderr, "\n");
2079 print_time(stderr, runtime, step, ir, cr);
2082 /* Replica exchange */
2084 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
2085 do_per_step(step, repl_ex_nst))
2087 bExchanged = replica_exchange(fplog, cr, repl_ex,
2088 state_global, enerd,
2091 if (bExchanged && DOMAINDECOMP(cr))
2093 dd_partition_system(fplog, step, cr, TRUE, 1,
2094 state_global, top_global, ir,
2095 state, &f, mdatoms, top, fr,
2096 vsite, shellfc, constr,
2097 nrnb, wcycle, FALSE);
2103 bStartingFromCpt = FALSE;
2105 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2106 /* With all integrators, except VV, we need to retain the pressure
2107 * at the current step for coupling at the next step.
2109 if ((state->flags & (1<<estPRES_PREV)) &&
2111 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
2113 /* Store the pressure in t_state for pressure coupling
2114 * at the next MD step.
2116 copy_mat(pres, state->pres_prev);
2119 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2121 if ( (membed != NULL) && (!bLastStep) )
2123 rescale_membed(step_rel, membed, state_global->x);
2130 /* read next frame from input trajectory */
2131 bNotLastFrame = read_next_frame(oenv, status, &rerun_fr);
2136 rerun_parallel_comm(cr, &rerun_fr, &bNotLastFrame);
2140 if (!bRerunMD || !rerun_fr.bStep)
2142 /* increase the MD step number */
2147 cycles = wallcycle_stop(wcycle, ewcSTEP);
2148 if (DOMAINDECOMP(cr) && wcycle)
2150 dd_cycles_add(cr->dd, cycles, ddCyclStep);
2153 if (bPMETuneRunning || bPMETuneTry)
2155 /* PME grid + cut-off optimization with GPUs or PME nodes */
2157 /* Count the total cycles over the last steps */
2158 cycles_pmes += cycles;
2160 /* We can only switch cut-off at NS steps */
2161 if (step % ir->nstlist == 0)
2163 /* PME grid + cut-off optimization with GPUs or PME nodes */
2166 if (DDMASTER(cr->dd))
2168 /* PME node load is too high, start tuning */
2169 bPMETuneRunning = (dd_pme_f_ratio(cr->dd) >= 1.05);
2171 dd_bcast(cr->dd, sizeof(gmx_bool), &bPMETuneRunning);
2173 if (bPMETuneRunning || step_rel > ir->nstlist*50)
2175 bPMETuneTry = FALSE;
2178 if (bPMETuneRunning)
2180 /* init_step might not be a multiple of nstlist,
2181 * but the first cycle is always skipped anyhow.
2184 pme_load_balance(pme_loadbal, cr,
2185 (bVerbose && MASTER(cr)) ? stderr : NULL,
2187 ir, state, cycles_pmes,
2188 fr->ic, fr->nbv, &fr->pmedata,
2191 /* Update constants in forcerec/inputrec to keep them in sync with fr->ic */
2192 fr->ewaldcoeff = fr->ic->ewaldcoeff;
2193 fr->rlist = fr->ic->rlist;
2194 fr->rlistlong = fr->ic->rlistlong;
2195 fr->rcoulomb = fr->ic->rcoulomb;
2196 fr->rvdw = fr->ic->rvdw;
2202 if (step_rel == wcycle_get_reset_counters(wcycle) ||
2203 gs.set[eglsRESETCOUNTERS] != 0)
2205 /* Reset all the counters related to performance over the run */
2206 reset_all_counters(fplog, cr, step, &step_rel, ir, wcycle, nrnb, runtime,
2207 fr->nbv != NULL && fr->nbv->bUseGPU ? fr->nbv->cu_nbv : NULL);
2208 wcycle_set_reset_counters(wcycle, -1);
2209 if (!(cr->duty & DUTY_PME))
2211 /* Tell our PME node to reset its counters */
2212 gmx_pme_send_resetcounters(cr, step);
2214 /* Correct max_hours for the elapsed time */
2215 max_hours -= run_time/(60.0*60.0);
2216 bResetCountersHalfMaxH = FALSE;
2217 gs.set[eglsRESETCOUNTERS] = 0;
2221 /* End of main MD loop */
2225 runtime_end(runtime);
2227 if (bRerunMD && MASTER(cr))
2232 if (!(cr->duty & DUTY_PME))
2234 /* Tell the PME only node to finish */
2235 gmx_pme_send_finish(cr);
2240 if (ir->nstcalcenergy > 0 && !bRerunMD)
2242 print_ebin(outf->fp_ene, FALSE, FALSE, FALSE, fplog, step, t,
2243 eprAVER, FALSE, mdebin, fcd, groups, &(ir->opts));
2251 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
2253 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)));
2254 fprintf(fplog, "Average number of atoms that crossed the half buffer length: %.1f\n\n", nlh.ab/nlh.nns);
2257 if (pme_loadbal != NULL)
2259 pme_loadbal_done(pme_loadbal, cr, fplog,
2260 fr->nbv != NULL && fr->nbv->bUseGPU);
2263 if (shellfc && fplog)
2265 fprintf(fplog, "Fraction of iterations that converged: %.2f %%\n",
2266 (nconverged*100.0)/step_rel);
2267 fprintf(fplog, "Average number of force evaluations per MD step: %.2f\n\n",
2271 if (repl_ex_nst > 0 && MASTER(cr))
2273 print_replica_exchange_statistics(fplog, repl_ex);
2276 runtime->nsteps_done = step_rel;