1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
71 #include "mpelogging.h"
78 #include "compute_io.h"
80 #include "checkpoint.h"
81 #include "mtop_util.h"
82 #include "sighandler.h"
97 double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
98 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
100 gmx_vsite_t *vsite,gmx_constr_t constr,
101 int stepout,t_inputrec *ir,
102 gmx_mtop_t *top_global,
104 t_state *state_global,
106 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
107 gmx_edsam_t ed,t_forcerec *fr,
108 int repl_ex_nst,int repl_ex_seed,
109 real cpt_period,real max_hours,
110 const char *deviceOptions,
112 gmx_runtime_t *runtime)
115 gmx_large_int_t step,step_rel;
118 gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
119 gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
120 bFirstStep,bStateFromTPX,bInitStep,bLastStep,
121 bBornRadii,bStartingFromCpt;
122 gmx_bool bDoDHDL=FALSE;
123 gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
124 bForceUpdate=FALSE,bCPT;
126 gmx_bool bMasterState;
127 int force_flags,cglo_flags;
128 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
133 t_state *bufstate=NULL;
134 matrix *scale_tot,pcoupl_mu,M,ebox;
137 gmx_repl_ex_t repl_ex=NULL;
141 t_mdebin *mdebin=NULL;
146 gmx_enerdata_t *enerd;
148 gmx_global_stat_t gstat;
149 gmx_update_t upd=NULL;
154 gmx_groups_t *groups;
155 gmx_ekindata_t *ekind, *ekind_save;
156 gmx_shellfc_t shellfc;
157 int count,nconverged=0;
160 gmx_bool bIonize=FALSE;
161 gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
163 gmx_bool bResetCountersHalfMaxH=FALSE;
164 gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
167 atom_id *grpindex=NULL;
169 t_coupl_rec *tcr=NULL;
170 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
171 matrix boxcopy={{0}},lastbox;
173 real fom,oldfom,veta_save,pcurr,scalevir,tracevir;
176 real saved_conserved_quantity = 0;
181 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
182 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
183 gmx_iterate_t iterate;
184 gmx_large_int_t multisim_nsteps=-1; /* number of steps to do before first multisim
185 simulation stops. If equal to zero, don't
186 communicate any more between multisims.*/
188 /* Temporary addition for FAHCORE checkpointing */
192 /* Check for special mdrun options */
193 bRerunMD = (Flags & MD_RERUN);
194 bIonize = (Flags & MD_IONIZE);
195 bFFscan = (Flags & MD_FFSCAN);
196 bAppend = (Flags & MD_APPENDFILES);
197 if (Flags & MD_RESETCOUNTERSHALFWAY)
201 /* Signal to reset the counters half the simulation steps. */
202 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
204 /* Signal to reset the counters halfway the simulation time. */
205 bResetCountersHalfMaxH = (max_hours > 0);
208 /* md-vv uses averaged full step velocities for T-control
209 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
210 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
212 if (bVV) /* to store the initial velocities while computing virial */
214 snew(cbuf,top_global->natoms);
216 /* all the iteratative cases - only if there are constraints */
217 bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
218 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
222 /* Since we don't know if the frames read are related in any way,
223 * rebuild the neighborlist at every step.
226 ir->nstcalcenergy = 1;
230 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
232 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
233 bGStatEveryStep = (nstglobalcomm == 1);
235 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
238 "To reduce the energy communication with nstlist = -1\n"
239 "the neighbor list validity should not be checked at every step,\n"
240 "this means that exact integration is not guaranteed.\n"
241 "The neighbor list validity is checked after:\n"
242 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
243 "In most cases this will result in exact integration.\n"
244 "This reduces the energy communication by a factor of 2 to 3.\n"
245 "If you want less energy communication, set nstlist > 3.\n\n");
248 if (bRerunMD || bFFscan)
252 groups = &top_global->groups;
255 init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
256 nrnb,top_global,&upd,
257 nfile,fnm,&outf,&mdebin,
258 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
260 clear_mat(total_vir);
262 /* Energy terms and groups */
264 init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
265 if (DOMAINDECOMP(cr))
271 snew(f,top_global->natoms);
274 /* Kinetic energy data */
276 init_ekindata(fplog,top_global,&(ir->opts),ekind);
277 /* needed for iteration of constraints */
279 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
280 /* Copy the cos acceleration to the groups struct */
281 ekind->cosacc.cos_accel = ir->cos_accel;
283 gstat = global_stat_init(ir);
286 /* Check for polarizable models and flexible constraints */
287 shellfc = init_shell_flexcon(fplog,
288 top_global,n_flexible_constraints(constr),
289 (ir->bContinuation ||
290 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
291 NULL : state_global->x);
295 #ifdef GMX_THREAD_MPI
296 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
298 set_deform_reference_box(upd,
299 deform_init_init_step_tpx,
300 deform_init_box_tpx);
301 #ifdef GMX_THREAD_MPI
302 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
307 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
308 if ((io > 2000) && MASTER(cr))
310 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
314 if (DOMAINDECOMP(cr)) {
315 top = dd_init_local_top(top_global);
318 dd_init_local_state(cr->dd,state_global,state);
320 if (DDMASTER(cr->dd) && ir->nstfout) {
321 snew(f_global,state_global->natoms);
325 /* Initialize the particle decomposition and split the topology */
326 top = split_system(fplog,top_global,ir,cr);
328 pd_cg_range(cr,&fr->cg0,&fr->hcg);
329 pd_at_range(cr,&a0,&a1);
331 top = gmx_mtop_generate_local_top(top_global,ir);
334 a1 = top_global->natoms;
337 state = partdec_init_local_state(cr,state_global);
340 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
343 set_vsite_top(vsite,top,mdatoms,cr);
346 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
347 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
351 make_local_shells(cr,mdatoms,shellfc);
354 if (ir->pull && PAR(cr)) {
355 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
359 if (DOMAINDECOMP(cr))
361 /* Distribute the charge groups over the nodes from the master node */
362 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
363 state_global,top_global,ir,
364 state,&f,mdatoms,top,fr,
365 vsite,shellfc,constr,
369 update_mdatoms(mdatoms,state->lambda);
373 if (opt2bSet("-cpi",nfile,fnm))
375 /* Update mdebin with energy history if appending to output files */
376 if ( Flags & MD_APPENDFILES )
378 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
382 /* We might have read an energy history from checkpoint,
383 * free the allocated memory and reset the counts.
385 done_energyhistory(&state_global->enerhist);
386 init_energyhistory(&state_global->enerhist);
389 /* Set the initial energy history in state by updating once */
390 update_energyhistory(&state_global->enerhist,mdebin);
393 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
394 /* Set the random state if we read a checkpoint file */
395 set_stochd_state(upd,state);
398 /* Initialize constraints */
400 if (!DOMAINDECOMP(cr))
401 set_constraints(constr,top,ir,mdatoms,cr);
404 /* Check whether we have to GCT stuff */
405 bTCR = ftp2bSet(efGCT,nfile,fnm);
408 fprintf(stderr,"Will do General Coupling Theory!\n");
410 gnx = top_global->mols.nr;
412 for(i=0; (i<gnx); i++) {
419 /* We need to be sure replica exchange can only occur
420 * when the energies are current */
421 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
422 "repl_ex_nst",&repl_ex_nst);
423 /* This check needs to happen before inter-simulation
424 * signals are initialized, too */
426 if (repl_ex_nst > 0 && MASTER(cr))
427 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
428 repl_ex_nst,repl_ex_seed);
430 if (!ir->bContinuation && !bRerunMD)
432 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
434 /* Set the velocities of frozen particles to zero */
435 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
439 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
449 /* Constrain the initial coordinates and velocities */
450 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
451 graph,cr,nrnb,fr,top,shake_vir);
455 /* Construct the virtual sites for the initial configuration */
456 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
457 top->idef.iparams,top->idef.il,
458 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
464 /* I'm assuming we need global communication the first time! MRS */
465 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
466 | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM:0)
467 | (bVV ? CGLO_PRESSURE:0)
468 | (bVV ? CGLO_CONSTRAINT:0)
469 | (bRerunMD ? CGLO_RERUNMD:0)
470 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
472 bSumEkinhOld = FALSE;
473 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
474 NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
475 constr,NULL,FALSE,state->box,
476 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
477 if (ir->eI == eiVVAK) {
478 /* a second call to get the half step temperature initialized as well */
479 /* we do the same call as above, but turn the pressure off -- internally to
480 compute_globals, this is recognized as a velocity verlet half-step
481 kinetic energy calculation. This minimized excess variables, but
482 perhaps loses some logic?*/
484 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
485 NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
486 constr,NULL,FALSE,state->box,
487 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
488 cglo_flags &~ (CGLO_STOPCM | CGLO_PRESSURE));
491 /* Calculate the initial half step temperature, and save the ekinh_old */
492 if (!(Flags & MD_STARTFROMCPT))
494 for(i=0; (i<ir->opts.ngtc); i++)
496 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
501 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
502 and there is no previous step */
505 /* if using an iterative algorithm, we need to create a working directory for the state. */
508 bufstate = init_bufstate(state);
512 snew(xcopy,state->natoms);
513 snew(vcopy,state->natoms);
514 copy_rvecn(state->x,xcopy,0,state->natoms);
515 copy_rvecn(state->v,vcopy,0,state->natoms);
516 copy_mat(state->box,boxcopy);
519 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
520 temperature control */
521 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
525 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
528 "RMS relative constraint deviation after constraining: %.2e\n",
529 constr_rmsd(constr,FALSE));
531 if (EI_STATE_VELOCITY(ir->eI))
533 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
537 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
538 " input trajectory '%s'\n\n",
539 *(top_global->name),opt2fn("-rerun",nfile,fnm));
542 fprintf(stderr,"Calculated time to finish depends on nsteps from "
543 "run input file,\nwhich may not correspond to the time "
544 "needed to process input trajectory.\n\n");
550 fprintf(stderr,"starting mdrun '%s'\n",
551 *(top_global->name));
554 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
558 sprintf(tbuf,"%s","infinite");
560 if (ir->init_step > 0)
562 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
563 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
564 gmx_step_str(ir->init_step,sbuf2),
565 ir->init_step*ir->delta_t);
569 fprintf(stderr,"%s steps, %s ps.\n",
570 gmx_step_str(ir->nsteps,sbuf),tbuf);
576 /* Set and write start time */
577 runtime_start(runtime);
578 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
579 wallcycle_start(wcycle,ewcRUN);
583 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
585 chkpt_ret=fcCheckPointParallel( cr->nodeid,
587 if ( chkpt_ret == 0 )
588 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
592 /***********************************************************
596 ************************************************************/
598 /* if rerunMD then read coordinates and velocities from input trajectory */
601 if (getenv("GMX_FORCE_UPDATE"))
609 bNotLastFrame = read_first_frame(oenv,&status,
610 opt2fn("-rerun",nfile,fnm),
611 &rerun_fr,TRX_NEED_X | TRX_READ_V);
612 if (rerun_fr.natoms != top_global->natoms)
615 "Number of atoms in trajectory (%d) does not match the "
616 "run input file (%d)\n",
617 rerun_fr.natoms,top_global->natoms);
619 if (ir->ePBC != epbcNONE)
623 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);
625 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
627 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
634 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
637 if (ir->ePBC != epbcNONE)
639 /* Set the shift vectors.
640 * Necessary here when have a static box different from the tpr box.
642 calc_shifts(rerun_fr.box,fr->shift_vec);
646 /* loop over MD steps or if rerunMD to end of input trajectory */
648 /* Skip the first Nose-Hoover integration when we get the state from tpx */
649 bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
650 bInitStep = bFirstStep && (bStateFromTPX || bVV);
651 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
653 bSumEkinhOld = FALSE;
656 init_global_signals(&gs,cr,ir,repl_ex_nst);
658 step = ir->init_step;
661 if (ir->nstlist == -1)
663 init_nlistheuristics(&nlh,bGStatEveryStep,step);
666 if (MULTISIM(cr) && (repl_ex_nst <=0 ))
668 /* check how many steps are left in other sims */
669 multisim_nsteps=get_multisim_nsteps(cr, ir->nsteps);
673 /* and stop now if we should */
674 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
675 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
676 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
678 wallcycle_start(wcycle,ewcSTEP);
680 GMX_MPE_LOG(ev_timestep1);
683 if (rerun_fr.bStep) {
684 step = rerun_fr.step;
685 step_rel = step - ir->init_step;
687 if (rerun_fr.bTime) {
697 bLastStep = (step_rel == ir->nsteps);
698 t = t0 + step*ir->delta_t;
701 if (ir->efep != efepNO)
703 if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
705 state_global->lambda = rerun_fr.lambda;
709 state_global->lambda = lam0 + step*ir->delta_lambda;
711 state->lambda = state_global->lambda;
712 bDoDHDL = do_per_step(step,ir->nstdhdl);
717 update_annealing_target_temp(&(ir->opts),t);
722 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
724 for(i=0; i<state_global->natoms; i++)
726 copy_rvec(rerun_fr.x[i],state_global->x[i]);
730 for(i=0; i<state_global->natoms; i++)
732 copy_rvec(rerun_fr.v[i],state_global->v[i]);
737 for(i=0; i<state_global->natoms; i++)
739 clear_rvec(state_global->v[i]);
743 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
744 " Ekin, temperature and pressure are incorrect,\n"
745 " the virial will be incorrect when constraints are present.\n"
747 bRerunWarnNoV = FALSE;
751 copy_mat(rerun_fr.box,state_global->box);
752 copy_mat(state_global->box,state->box);
754 if (vsite && (Flags & MD_RERUN_VSITE))
756 if (DOMAINDECOMP(cr))
758 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
762 /* Following is necessary because the graph may get out of sync
763 * with the coordinates if we only have every N'th coordinate set
765 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
766 shift_self(graph,state->box,state->x);
768 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
769 top->idef.iparams,top->idef.il,
770 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
773 unshift_self(graph,state->box,state->x);
778 /* Stop Center of Mass motion */
779 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
781 /* Copy back starting coordinates in case we're doing a forcefield scan */
784 for(ii=0; (ii<state->natoms); ii++)
786 copy_rvec(xcopy[ii],state->x[ii]);
787 copy_rvec(vcopy[ii],state->v[ii]);
789 copy_mat(boxcopy,state->box);
794 /* for rerun MD always do Neighbour Searching */
795 bNS = (bFirstStep || ir->nstlist != 0);
800 /* Determine whether or not to do Neighbour Searching and LR */
801 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
803 bNS = (bFirstStep || bExchanged || bNStList ||
804 (ir->nstlist == -1 && nlh.nabnsb > 0));
806 if (bNS && ir->nstlist == -1)
808 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
812 /* check whether we should stop because another simulation has
816 if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
817 (multisim_nsteps != ir->nsteps) )
824 "Stopping simulation %d because another one has finished\n",
828 gs.sig[eglsCHKPT] = 1;
833 /* < 0 means stop at next step, > 0 means stop at next NS step */
834 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
835 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
840 /* Determine whether or not to update the Born radii if doing GB */
841 bBornRadii=bFirstStep;
842 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
847 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
848 do_verbose = bVerbose &&
849 (step % stepout == 0 || bFirstStep || bLastStep);
851 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
859 bMasterState = FALSE;
860 /* Correct the new box if it is too skewed */
861 if (DYNAMIC_BOX(*ir))
863 if (correct_box(fplog,step,state->box,graph))
868 if (DOMAINDECOMP(cr) && bMasterState)
870 dd_collect_state(cr->dd,state,state_global);
874 if (DOMAINDECOMP(cr))
876 /* Repartition the domain decomposition */
877 wallcycle_start(wcycle,ewcDOMDEC);
878 dd_partition_system(fplog,step,cr,
879 bMasterState,nstglobalcomm,
880 state_global,top_global,ir,
881 state,&f,mdatoms,top,fr,
882 vsite,shellfc,constr,
883 nrnb,wcycle,do_verbose);
884 wallcycle_stop(wcycle,ewcDOMDEC);
885 /* If using an iterative integrator, reallocate space to match the decomposition */
889 if (MASTER(cr) && do_log && !bFFscan)
891 print_ebin_header(fplog,step,t,state->lambda);
894 if (ir->efep != efepNO)
896 update_mdatoms(mdatoms,state->lambda);
899 if (bRerunMD && rerun_fr.bV)
902 /* We need the kinetic energy at minus the half step for determining
903 * the full step kinetic energy and possibly for T-coupling.*/
904 /* This may not be quite working correctly yet . . . . */
905 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
906 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
907 constr,NULL,FALSE,state->box,
908 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
909 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
911 clear_mat(force_vir);
913 /* Ionize the atoms if necessary */
916 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
917 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
920 /* Update force field in ffscan program */
923 if (update_forcefield(fplog,
925 mdatoms->nr,state->x,state->box)) {
926 if (gmx_parallel_env_initialized())
934 GMX_MPE_LOG(ev_timestep2);
936 /* We write a checkpoint at this MD step when:
937 * either at an NS step when we signalled through gs,
938 * or at the last step (but not when we do not want confout),
939 * but never at the first step or with rerun.
941 bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
942 (bLastStep && (Flags & MD_CONFOUT))) &&
943 step > ir->init_step && !bRerunMD);
946 gs.set[eglsCHKPT] = 0;
949 /* Determine the energy and pressure:
950 * at nstcalcenergy steps and at energy output steps (set below).
952 bNstEner = do_per_step(step,ir->nstcalcenergy);
955 (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)));
957 /* Do we need global communication ? */
958 bGStat = (bCalcEnerPres || bStopCM ||
959 do_per_step(step,nstglobalcomm) ||
960 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
962 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
964 if (do_ene || do_log)
966 bCalcEnerPres = TRUE;
970 /* these CGLO_ options remain the same throughout the iteration */
971 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
972 (bGStat ? CGLO_GSTAT : 0)
975 force_flags = (GMX_FORCE_STATECHANGED |
976 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
977 GMX_FORCE_ALLFORCES |
978 (bNStList ? GMX_FORCE_DOLR : 0) |
980 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
981 (bDoDHDL ? GMX_FORCE_DHDL : 0)
986 /* Now is the time to relax the shells */
987 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
989 bStopCM,top,top_global,
991 state,f,force_vir,mdatoms,
992 nrnb,wcycle,graph,groups,
993 shellfc,fr,bBornRadii,t,mu_tot,
994 state->natoms,&bConverged,vsite,
1005 /* The coordinates (x) are shifted (to get whole molecules)
1007 * This is parallellized as well, and does communication too.
1008 * Check comments in sim_util.c
1011 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1012 state->box,state->x,&state->hist,
1013 f,force_vir,mdatoms,enerd,fcd,
1014 state->lambda,graph,
1015 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1016 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1019 GMX_BARRIER(cr->mpi_comm_mygroup);
1023 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1024 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1027 if (bTCR && bFirstStep)
1029 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1030 fprintf(fplog,"Done init_coupling\n");
1034 if (bVV && !bStartingFromCpt && !bRerunMD)
1035 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1037 if (ir->eI==eiVV && bInitStep)
1039 /* if using velocity verlet with full time step Ekin,
1040 * take the first half step only to compute the
1041 * virial for the first step. From there,
1042 * revert back to the initial coordinates
1043 * so that the input is actually the initial step.
1045 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1047 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1048 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1051 update_coords(fplog,step,ir,mdatoms,state,
1052 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1053 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1054 cr,nrnb,constr,&top->idef);
1058 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1060 /* for iterations, we save these vectors, as we will be self-consistently iterating
1063 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1065 /* save the state */
1066 if (bIterations && iterate.bIterate) {
1067 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1070 bFirstIterate = TRUE;
1071 while (bFirstIterate || (bIterations && iterate.bIterate))
1073 if (bIterations && iterate.bIterate)
1075 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1076 if (bFirstIterate && bTrotter)
1078 /* The first time through, we need a decent first estimate
1079 of veta(t+dt) to compute the constraints. Do
1080 this by computing the box volume part of the
1081 trotter integration at this time. Nothing else
1082 should be changed by this routine here. If
1083 !(first time), we start with the previous value
1086 veta_save = state->veta;
1087 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1088 vetanew = state->veta;
1089 state->veta = veta_save;
1094 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
1097 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1098 &top->idef,shake_vir,NULL,
1099 cr,nrnb,wcycle,upd,constr,
1100 bInitStep,TRUE,bCalcEnerPres,vetanew);
1102 if (!bOK && !bFFscan)
1104 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1109 { /* Need to unshift here if a do_force has been
1110 called in the previous step */
1111 unshift_self(graph,state->box,state->x);
1115 /* if VV, compute the pressure and constraints */
1116 /* For VV2, we strictly only need this if using pressure
1117 * control, but we really would like to have accurate pressures
1119 * Think about ways around this in the future?
1120 * For now, keep this choice in comments.
1122 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1123 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1125 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
1126 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1127 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1128 constr,NULL,FALSE,state->box,
1129 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1132 | (bStopCM ? CGLO_STOPCM : 0)
1133 | (bTemp ? CGLO_TEMPERATURE:0)
1134 | (bPres ? CGLO_PRESSURE : 0)
1135 | (bPres ? CGLO_CONSTRAINT : 0)
1136 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
1137 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1140 /* explanation of above:
1141 a) We compute Ekin at the full time step
1142 if 1) we are using the AveVel Ekin, and it's not the
1143 initial step, or 2) if we are using AveEkin, but need the full
1144 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1145 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1146 EkinAveVel because it's needed for the pressure */
1148 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1153 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
1157 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1162 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1163 state->veta,&vetanew))
1167 bFirstIterate = FALSE;
1170 if (bTrotter && !bInitStep) {
1171 copy_mat(shake_vir,state->svir_prev);
1172 copy_mat(force_vir,state->fvir_prev);
1173 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
1174 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1175 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
1176 enerd->term[F_EKIN] = trace(ekind->ekin);
1179 /* if it's the initial step, we performed this first step just to get the constraint virial */
1180 if (bInitStep && ir->eI==eiVV) {
1181 copy_rvecn(cbuf,state->v,0,state->natoms);
1184 if (fr->bSepDVDL && fplog && do_log)
1186 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1188 enerd->term[F_DHDL_CON] += dvdl;
1190 GMX_MPE_LOG(ev_timestep1);
1193 /* MRS -- now done iterating -- compute the conserved quantity */
1195 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
1198 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
1200 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1202 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1206 /* ######## END FIRST UPDATE STEP ############## */
1207 /* ######## If doing VV, we now have v(dt) ###### */
1209 /* ################## START TRAJECTORY OUTPUT ################# */
1211 /* Now we have the energies and forces corresponding to the
1212 * coordinates at time t. We must output all of this before
1214 * for RerunMD t is read from input trajectory
1216 GMX_MPE_LOG(ev_output_start);
1219 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
1220 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
1221 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
1222 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
1223 if (bCPT) { mdof_flags |= MDOF_CPT; };
1225 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1228 /* Enforce writing positions and velocities at end of run */
1229 mdof_flags |= (MDOF_X | MDOF_V);
1234 fcReportProgress( ir->nsteps, step );
1236 /* sync bCPT and fc record-keeping */
1237 if (bCPT && MASTER(cr))
1238 fcRequestCheckPoint();
1241 if (mdof_flags != 0)
1243 wallcycle_start(wcycle,ewcTRAJ);
1246 if (state->flags & (1<<estLD_RNG))
1248 get_stochd_state(upd,state);
1254 state_global->ekinstate.bUpToDate = FALSE;
1258 update_ekinstate(&state_global->ekinstate,ekind);
1259 state_global->ekinstate.bUpToDate = TRUE;
1261 update_energyhistory(&state_global->enerhist,mdebin);
1264 write_traj(fplog,cr,outf,mdof_flags,top_global,
1265 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
1272 if (bLastStep && step_rel == ir->nsteps &&
1273 (Flags & MD_CONFOUT) && MASTER(cr) &&
1274 !bRerunMD && !bFFscan)
1276 /* x and v have been collected in write_traj,
1277 * because a checkpoint file will always be written
1280 fprintf(stderr,"\nWriting final coordinates.\n");
1281 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
1284 /* Make molecules whole only for confout writing */
1285 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
1287 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
1288 *top_global->name,top_global,
1289 state_global->x,state_global->v,
1290 ir->ePBC,state->box);
1293 wallcycle_stop(wcycle,ewcTRAJ);
1295 GMX_MPE_LOG(ev_output_finish);
1297 /* kludge -- virial is lost with restart for NPT control. Must restart */
1298 if (bStartingFromCpt && bVV)
1300 copy_mat(state->svir_prev,shake_vir);
1301 copy_mat(state->fvir_prev,force_vir);
1303 /* ################## END TRAJECTORY OUTPUT ################ */
1305 /* Determine the wallclock run time up till now */
1306 run_time = gmx_gettime() - (double)runtime->real;
1308 /* Check whether everything is still allright */
1309 if (((int)gmx_get_stop_condition() > handled_stop_condition)
1310 #ifdef GMX_THREAD_MPI
1315 /* this is just make gs.sig compatible with the hack
1316 of sending signals around by MPI_Reduce with together with
1318 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
1319 gs.sig[eglsSTOPCOND]=1;
1320 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
1321 gs.sig[eglsSTOPCOND]=-1;
1322 /* < 0 means stop at next step, > 0 means stop at next NS step */
1326 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1327 gmx_get_signal_name(),
1328 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1332 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1333 gmx_get_signal_name(),
1334 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1336 handled_stop_condition=(int)gmx_get_stop_condition();
1338 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1339 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
1340 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
1342 /* Signal to terminate the run */
1343 gs.sig[eglsSTOPCOND] = 1;
1346 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1348 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1351 if (bResetCountersHalfMaxH && MASTER(cr) &&
1352 run_time > max_hours*60.0*60.0*0.495)
1354 gs.sig[eglsRESETCOUNTERS] = 1;
1357 if (ir->nstlist == -1 && !bRerunMD)
1359 /* When bGStatEveryStep=FALSE, global_stat is only called
1360 * when we check the atom displacements, not at NS steps.
1361 * This means that also the bonded interaction count check is not
1362 * performed immediately after NS. Therefore a few MD steps could
1363 * be performed with missing interactions.
1364 * But wrong energies are never written to file,
1365 * since energies are only written after global_stat
1368 if (step >= nlh.step_nscheck)
1370 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
1371 nlh.scale_tot,state->x);
1375 /* This is not necessarily true,
1376 * but step_nscheck is determined quite conservatively.
1382 /* In parallel we only have to check for checkpointing in steps
1383 * where we do global communication,
1384 * otherwise the other nodes don't know.
1386 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1389 run_time >= nchkpt*cpt_period*60.0)) &&
1390 gs.set[eglsCHKPT] == 0)
1392 gs.sig[eglsCHKPT] = 1;
1397 gmx_iterate_init(&iterate,bIterations);
1400 /* for iterations, we save these vectors, as we will be redoing the calculations */
1401 if (bIterations && iterate.bIterate)
1403 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1405 bFirstIterate = TRUE;
1406 while (bFirstIterate || (bIterations && iterate.bIterate))
1408 /* We now restore these vectors to redo the calculation with improved extended variables */
1411 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1414 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1415 so scroll down for that logic */
1417 /* ######### START SECOND UPDATE STEP ################# */
1418 GMX_MPE_LOG(ev_update_start);
1419 /* Box is changed in update() when we do pressure coupling,
1420 * but we should still use the old box for energy corrections and when
1421 * writing it to the energy file, so it matches the trajectory files for
1422 * the same timestep above. Make a copy in a separate array.
1424 copy_mat(state->box,lastbox);
1427 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
1429 wallcycle_start(wcycle,ewcUPDATE);
1431 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1434 if (bIterations && iterate.bIterate)
1442 /* we use a new value of scalevir to converge the iterations faster */
1443 scalevir = tracevir/trace(shake_vir);
1445 msmul(shake_vir,scalevir,shake_vir);
1446 m_add(force_vir,shake_vir,total_vir);
1447 clear_mat(shake_vir);
1449 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
1450 /* We can only do Berendsen coupling after we have summed
1451 * the kinetic energy or virial. Since the happens
1452 * in global_state after update, we should only do it at
1453 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1458 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1459 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
1465 /* velocity half-step update */
1466 update_coords(fplog,step,ir,mdatoms,state,f,
1467 fr->bTwinRange && bNStList,fr->f_twin,fcd,
1468 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
1469 cr,nrnb,constr,&top->idef);
1472 /* Above, initialize just copies ekinh into ekin,
1473 * it doesn't copy position (for VV),
1474 * and entire integrator for MD.
1479 copy_rvecn(state->x,cbuf,0,state->natoms);
1482 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1483 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1484 wallcycle_stop(wcycle,ewcUPDATE);
1486 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1487 &top->idef,shake_vir,force_vir,
1488 cr,nrnb,wcycle,upd,constr,
1489 bInitStep,FALSE,bCalcEnerPres,state->veta);
1493 /* erase F_EKIN and F_TEMP here? */
1494 /* just compute the kinetic energy at the half step to perform a trotter step */
1495 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1496 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1497 constr,NULL,FALSE,lastbox,
1498 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1499 cglo_flags | CGLO_TEMPERATURE
1501 wallcycle_start(wcycle,ewcUPDATE);
1502 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
1503 /* now we know the scaling, we can compute the positions again again */
1504 copy_rvecn(cbuf,state->x,0,state->natoms);
1506 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1507 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1508 wallcycle_stop(wcycle,ewcUPDATE);
1510 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1511 /* are the small terms in the shake_vir here due
1512 * to numerical errors, or are they important
1513 * physically? I'm thinking they are just errors, but not completely sure.
1514 * For now, will call without actually constraining, constr=NULL*/
1515 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1516 &top->idef,tmp_vir,force_vir,
1517 cr,nrnb,wcycle,upd,NULL,
1518 bInitStep,FALSE,bCalcEnerPres,
1521 if (!bOK && !bFFscan)
1523 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1526 if (fr->bSepDVDL && fplog && do_log)
1528 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1530 enerd->term[F_DHDL_CON] += dvdl;
1534 /* Need to unshift here */
1535 unshift_self(graph,state->box,state->x);
1538 GMX_BARRIER(cr->mpi_comm_mygroup);
1539 GMX_MPE_LOG(ev_update_finish);
1543 wallcycle_start(wcycle,ewcVSITECONSTR);
1546 shift_self(graph,state->box,state->x);
1548 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1549 top->idef.iparams,top->idef.il,
1550 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1554 unshift_self(graph,state->box,state->x);
1556 wallcycle_stop(wcycle,ewcVSITECONSTR);
1559 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1560 if (ir->nstlist == -1 && bFirstIterate)
1562 gs.sig[eglsNABNSB] = nlh.nabnsb;
1564 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1565 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1567 bFirstIterate ? &gs : NULL,
1568 (step_rel % gs.nstms == 0) &&
1569 (multisim_nsteps<0 || (step_rel<multisim_nsteps)),
1571 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1573 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
1574 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1575 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1576 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1577 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
1578 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1581 if (ir->nstlist == -1 && bFirstIterate)
1583 nlh.nabnsb = gs.set[eglsNABNSB];
1584 gs.set[eglsNABNSB] = 0;
1586 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1587 /* ############# END CALC EKIN AND PRESSURE ################# */
1589 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1590 the virial that should probably be addressed eventually. state->veta has better properies,
1591 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1592 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1595 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1596 trace(shake_vir),&tracevir))
1600 bFirstIterate = FALSE;
1603 update_box(fplog,step,ir,mdatoms,state,graph,f,
1604 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
1606 /* ################# END UPDATE STEP 2 ################# */
1607 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1609 /* The coordinates (x) were unshifted in update */
1610 if (bFFscan && (shellfc==NULL || bConverged))
1612 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
1614 &(top_global->mols),mdatoms->massT,pres))
1616 if (gmx_parallel_env_initialized())
1620 fprintf(stderr,"\n");
1626 /* We will not sum ekinh_old,
1627 * so signal that we still have to do it.
1629 bSumEkinhOld = TRUE;
1634 /* Only do GCT when the relaxation of shells (minimization) has converged,
1635 * otherwise we might be coupling to bogus energies.
1636 * In parallel we must always do this, because the other sims might
1640 /* Since this is called with the new coordinates state->x, I assume
1641 * we want the new box state->box too. / EL 20040121
1643 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
1645 mdatoms,&(top->idef),mu_aver,
1646 top_global->mols.nr,cr,
1647 state->box,total_vir,pres,
1648 mu_tot,state->x,f,bConverged);
1652 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1654 /* sum up the foreign energy and dhdl terms */
1655 sum_dhdl(enerd,state->lambda,ir);
1657 /* use the directly determined last velocity, not actually the averaged half steps */
1658 if (bTrotter && ir->eI==eiVV)
1660 enerd->term[F_EKIN] = last_ekin;
1662 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1666 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1670 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
1672 /* Check for excessively large energies */
1676 real etot_max = 1e200;
1678 real etot_max = 1e30;
1680 if (fabs(enerd->term[F_ETOT]) > etot_max)
1682 fprintf(stderr,"Energy too large (%g), giving up\n",
1683 enerd->term[F_ETOT]);
1686 /* ######### END PREPARING EDR OUTPUT ########### */
1688 /* Time for performance */
1689 if (((step % stepout) == 0) || bLastStep)
1691 runtime_upd_proc(runtime);
1697 gmx_bool do_dr,do_or;
1699 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
1703 upd_mdebin(mdebin,bDoDHDL, TRUE,
1704 t,mdatoms->tmass,enerd,state,lastbox,
1705 shake_vir,force_vir,total_vir,pres,
1706 ekind,mu_tot,constr);
1710 upd_mdebin_step(mdebin);
1713 do_dr = do_per_step(step,ir->nstdisreout);
1714 do_or = do_per_step(step,ir->nstorireout);
1716 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
1718 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
1720 if (ir->ePull != epullNO)
1722 pull_print_output(ir->pull,step,t);
1725 if (do_per_step(step,ir->nstlog))
1727 if(fflush(fplog) != 0)
1729 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
1735 /* Remaining runtime */
1736 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
1740 fprintf(stderr,"\n");
1742 print_time(stderr,runtime,step,ir,cr);
1745 /* Replica exchange */
1747 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
1748 do_per_step(step,repl_ex_nst))
1750 bExchanged = replica_exchange(fplog,cr,repl_ex,
1751 state_global,enerd->term,
1754 if (bExchanged && DOMAINDECOMP(cr))
1756 dd_partition_system(fplog,step,cr,TRUE,1,
1757 state_global,top_global,ir,
1758 state,&f,mdatoms,top,fr,
1759 vsite,shellfc,constr,
1766 bStartingFromCpt = FALSE;
1768 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1769 /* With all integrators, except VV, we need to retain the pressure
1770 * at the current step for coupling at the next step.
1772 if ((state->flags & (1<<estPRES_PREV)) &&
1774 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1776 /* Store the pressure in t_state for pressure coupling
1777 * at the next MD step.
1779 copy_mat(pres,state->pres_prev);
1782 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1788 /* read next frame from input trajectory */
1789 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
1794 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
1798 if (!bRerunMD || !rerun_fr.bStep)
1800 /* increase the MD step number */
1805 cycles = wallcycle_stop(wcycle,ewcSTEP);
1806 if (DOMAINDECOMP(cr) && wcycle)
1808 dd_cycles_add(cr->dd,cycles,ddCyclStep);
1811 if (step_rel == wcycle_get_reset_counters(wcycle) ||
1812 gs.set[eglsRESETCOUNTERS] != 0)
1814 /* Reset all the counters related to performance over the run */
1815 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
1816 wcycle_set_reset_counters(wcycle,-1);
1817 /* Correct max_hours for the elapsed time */
1818 max_hours -= run_time/(60.0*60.0);
1819 bResetCountersHalfMaxH = FALSE;
1820 gs.set[eglsRESETCOUNTERS] = 0;
1824 /* End of main MD loop */
1828 runtime_end(runtime);
1830 if (bRerunMD && MASTER(cr))
1835 if (!(cr->duty & DUTY_PME))
1837 /* Tell the PME only node to finish */
1843 if (ir->nstcalcenergy > 0 && !bRerunMD)
1845 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
1846 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
1854 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
1856 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)));
1857 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
1860 if (shellfc && fplog)
1862 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
1863 (nconverged*100.0)/step_rel);
1864 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
1868 if (repl_ex_nst > 0 && MASTER(cr))
1870 print_replica_exchange_statistics(fplog,repl_ex);
1873 runtime->nsteps_done = step_rel;