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
40 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
77 #include "mpelogging.h"
84 #include "compute_io.h"
86 #include "checkpoint.h"
87 #include "mtop_util.h"
88 #include "sighandler.h"
100 #include "corewrap.h"
104 double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
105 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
107 gmx_vsite_t *vsite,gmx_constr_t constr,
108 int stepout,t_inputrec *ir,
109 gmx_mtop_t *top_global,
111 t_state *state_global,
113 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
114 gmx_edsam_t ed,t_forcerec *fr,
115 int repl_ex_nst,int repl_ex_seed,gmx_membed_t *membed,
116 real cpt_period,real max_hours,
117 const char *deviceOptions,
119 gmx_runtime_t *runtime)
122 gmx_large_int_t step,step_rel;
125 gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
126 gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
127 bFirstStep,bStateFromTPX,bInitStep,bLastStep,
128 bBornRadii,bStartingFromCpt;
129 gmx_bool bDoDHDL=FALSE;
130 gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
131 bForceUpdate=FALSE,bCPT;
133 gmx_bool bMasterState;
134 int force_flags,cglo_flags;
135 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
140 t_state *bufstate=NULL;
141 matrix *scale_tot,pcoupl_mu,M,ebox;
144 gmx_repl_ex_t repl_ex=NULL;
148 t_mdebin *mdebin=NULL;
153 gmx_enerdata_t *enerd;
155 gmx_global_stat_t gstat;
156 gmx_update_t upd=NULL;
161 gmx_groups_t *groups;
162 gmx_ekindata_t *ekind, *ekind_save;
163 gmx_shellfc_t shellfc;
164 int count,nconverged=0;
167 gmx_bool bIonize=FALSE;
168 gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
170 gmx_bool bResetCountersHalfMaxH=FALSE;
171 gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
172 real temp0,mu_aver=0,dvdl;
174 atom_id *grpindex=NULL;
176 t_coupl_rec *tcr=NULL;
177 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
178 matrix boxcopy={{0}},lastbox;
180 real fom,oldfom,veta_save,pcurr,scalevir,tracevir;
183 real saved_conserved_quantity = 0;
188 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
189 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
190 gmx_iterate_t iterate;
191 gmx_large_int_t multisim_nsteps=-1; /* number of steps to do before first multisim
192 simulation stops. If equal to zero, don't
193 communicate any more between multisims.*/
195 /* Temporary addition for FAHCORE checkpointing */
199 /* Check for special mdrun options */
200 bRerunMD = (Flags & MD_RERUN);
201 bIonize = (Flags & MD_IONIZE);
202 bFFscan = (Flags & MD_FFSCAN);
203 bAppend = (Flags & MD_APPENDFILES);
204 if (Flags & MD_RESETCOUNTERSHALFWAY)
208 /* Signal to reset the counters half the simulation steps. */
209 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
211 /* Signal to reset the counters halfway the simulation time. */
212 bResetCountersHalfMaxH = (max_hours > 0);
215 /* md-vv uses averaged full step velocities for T-control
216 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
217 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
219 if (bVV) /* to store the initial velocities while computing virial */
221 snew(cbuf,top_global->natoms);
223 /* all the iteratative cases - only if there are constraints */
224 bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
225 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
229 /* Since we don't know if the frames read are related in any way,
230 * rebuild the neighborlist at every step.
233 ir->nstcalcenergy = 1;
237 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
239 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
240 bGStatEveryStep = (nstglobalcomm == 1);
242 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
245 "To reduce the energy communication with nstlist = -1\n"
246 "the neighbor list validity should not be checked at every step,\n"
247 "this means that exact integration is not guaranteed.\n"
248 "The neighbor list validity is checked after:\n"
249 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
250 "In most cases this will result in exact integration.\n"
251 "This reduces the energy communication by a factor of 2 to 3.\n"
252 "If you want less energy communication, set nstlist > 3.\n\n");
255 if (bRerunMD || bFFscan)
259 groups = &top_global->groups;
262 init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
263 nrnb,top_global,&upd,
264 nfile,fnm,&outf,&mdebin,
265 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
267 clear_mat(total_vir);
269 /* Energy terms and groups */
271 init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
272 if (DOMAINDECOMP(cr))
278 snew(f,top_global->natoms);
281 /* Kinetic energy data */
283 init_ekindata(fplog,top_global,&(ir->opts),ekind);
284 /* needed for iteration of constraints */
286 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
287 /* Copy the cos acceleration to the groups struct */
288 ekind->cosacc.cos_accel = ir->cos_accel;
290 gstat = global_stat_init(ir);
293 /* Check for polarizable models and flexible constraints */
294 shellfc = init_shell_flexcon(fplog,
295 top_global,n_flexible_constraints(constr),
296 (ir->bContinuation ||
297 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
298 NULL : state_global->x);
303 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
305 set_deform_reference_box(upd,
306 deform_init_init_step_tpx,
307 deform_init_box_tpx);
309 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
314 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
315 if ((io > 2000) && MASTER(cr))
317 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
321 if (DOMAINDECOMP(cr)) {
322 top = dd_init_local_top(top_global);
325 dd_init_local_state(cr->dd,state_global,state);
327 if (DDMASTER(cr->dd) && ir->nstfout) {
328 snew(f_global,state_global->natoms);
332 /* Initialize the particle decomposition and split the topology */
333 top = split_system(fplog,top_global,ir,cr);
335 pd_cg_range(cr,&fr->cg0,&fr->hcg);
336 pd_at_range(cr,&a0,&a1);
338 top = gmx_mtop_generate_local_top(top_global,ir);
341 a1 = top_global->natoms;
344 state = partdec_init_local_state(cr,state_global);
347 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
350 set_vsite_top(vsite,top,mdatoms,cr);
353 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
354 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
358 make_local_shells(cr,mdatoms,shellfc);
361 if (ir->pull && PAR(cr)) {
362 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
366 if (DOMAINDECOMP(cr))
368 /* Distribute the charge groups over the nodes from the master node */
369 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
370 state_global,top_global,ir,
371 state,&f,mdatoms,top,fr,
372 vsite,shellfc,constr,
376 update_mdatoms(mdatoms,state->lambda);
380 if (opt2bSet("-cpi",nfile,fnm))
382 /* Update mdebin with energy history if appending to output files */
383 if ( Flags & MD_APPENDFILES )
385 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
389 /* We might have read an energy history from checkpoint,
390 * free the allocated memory and reset the counts.
392 done_energyhistory(&state_global->enerhist);
393 init_energyhistory(&state_global->enerhist);
396 /* Set the initial energy history in state by updating once */
397 update_energyhistory(&state_global->enerhist,mdebin);
400 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
401 /* Set the random state if we read a checkpoint file */
402 set_stochd_state(upd,state);
405 /* Initialize constraints */
407 if (!DOMAINDECOMP(cr))
408 set_constraints(constr,top,ir,mdatoms,cr);
411 /* Check whether we have to GCT stuff */
412 bTCR = ftp2bSet(efGCT,nfile,fnm);
415 fprintf(stderr,"Will do General Coupling Theory!\n");
417 gnx = top_global->mols.nr;
419 for(i=0; (i<gnx); i++) {
426 /* We need to be sure replica exchange can only occur
427 * when the energies are current */
428 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
429 "repl_ex_nst",&repl_ex_nst);
430 /* This check needs to happen before inter-simulation
431 * signals are initialized, too */
433 if (repl_ex_nst > 0 && MASTER(cr))
434 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
435 repl_ex_nst,repl_ex_seed);
437 if (!ir->bContinuation && !bRerunMD)
439 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
441 /* Set the velocities of frozen particles to zero */
442 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
446 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
456 /* Constrain the initial coordinates and velocities */
457 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
458 graph,cr,nrnb,fr,top,shake_vir);
462 /* Construct the virtual sites for the initial configuration */
463 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
464 top->idef.iparams,top->idef.il,
465 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
471 /* I'm assuming we need global communication the first time! MRS */
472 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
473 | (bVV ? CGLO_PRESSURE:0)
474 | (bVV ? CGLO_CONSTRAINT:0)
475 | (bRerunMD ? CGLO_RERUNMD:0)
476 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
478 bSumEkinhOld = FALSE;
479 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
480 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
481 constr,NULL,FALSE,state->box,
482 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
483 if (ir->eI == eiVVAK) {
484 /* a second call to get the half step temperature initialized as well */
485 /* we do the same call as above, but turn the pressure off -- internally to
486 compute_globals, this is recognized as a velocity verlet half-step
487 kinetic energy calculation. This minimized excess variables, but
488 perhaps loses some logic?*/
490 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
491 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
492 constr,NULL,FALSE,state->box,
493 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
494 cglo_flags &~ CGLO_PRESSURE);
497 /* Calculate the initial half step temperature, and save the ekinh_old */
498 if (!(Flags & MD_STARTFROMCPT))
500 for(i=0; (i<ir->opts.ngtc); i++)
502 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
507 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
508 and there is no previous step */
510 temp0 = enerd->term[F_TEMP];
512 /* if using an iterative algorithm, we need to create a working directory for the state. */
515 bufstate = init_bufstate(state);
519 snew(xcopy,state->natoms);
520 snew(vcopy,state->natoms);
521 copy_rvecn(state->x,xcopy,0,state->natoms);
522 copy_rvecn(state->v,vcopy,0,state->natoms);
523 copy_mat(state->box,boxcopy);
526 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
527 temperature control */
528 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
532 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
535 "RMS relative constraint deviation after constraining: %.2e\n",
536 constr_rmsd(constr,FALSE));
538 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
541 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
542 " input trajectory '%s'\n\n",
543 *(top_global->name),opt2fn("-rerun",nfile,fnm));
546 fprintf(stderr,"Calculated time to finish depends on nsteps from "
547 "run input file,\nwhich may not correspond to the time "
548 "needed to process input trajectory.\n\n");
554 fprintf(stderr,"starting mdrun '%s'\n",
555 *(top_global->name));
558 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
562 sprintf(tbuf,"%s","infinite");
564 if (ir->init_step > 0)
566 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
567 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
568 gmx_step_str(ir->init_step,sbuf2),
569 ir->init_step*ir->delta_t);
573 fprintf(stderr,"%s steps, %s ps.\n",
574 gmx_step_str(ir->nsteps,sbuf),tbuf);
580 /* Set and write start time */
581 runtime_start(runtime);
582 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
583 wallcycle_start(wcycle,ewcRUN);
587 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
589 chkpt_ret=fcCheckPointParallel( cr->nodeid,
591 if ( chkpt_ret == 0 )
592 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
596 /***********************************************************
600 ************************************************************/
602 /* if rerunMD then read coordinates and velocities from input trajectory */
605 if (getenv("GMX_FORCE_UPDATE"))
613 bNotLastFrame = read_first_frame(oenv,&status,
614 opt2fn("-rerun",nfile,fnm),
615 &rerun_fr,TRX_NEED_X | TRX_READ_V);
616 if (rerun_fr.natoms != top_global->natoms)
619 "Number of atoms in trajectory (%d) does not match the "
620 "run input file (%d)\n",
621 rerun_fr.natoms,top_global->natoms);
623 if (ir->ePBC != epbcNONE)
627 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);
629 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
631 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
638 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
641 if (ir->ePBC != epbcNONE)
643 /* Set the shift vectors.
644 * Necessary here when have a static box different from the tpr box.
646 calc_shifts(rerun_fr.box,fr->shift_vec);
650 /* loop over MD steps or if rerunMD to end of input trajectory */
652 /* Skip the first Nose-Hoover integration when we get the state from tpx */
653 bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
654 bInitStep = bFirstStep && (bStateFromTPX || bVV);
655 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
657 bSumEkinhOld = FALSE;
660 init_global_signals(&gs,cr,ir,repl_ex_nst);
662 step = ir->init_step;
665 if (ir->nstlist == -1)
667 init_nlistheuristics(&nlh,bGStatEveryStep,step);
670 if (MULTISIM(cr) && (repl_ex_nst <=0 ))
672 /* check how many steps are left in other sims */
673 multisim_nsteps=get_multisim_nsteps(cr, ir->nsteps);
677 /* and stop now if we should */
678 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
679 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
680 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
682 wallcycle_start(wcycle,ewcSTEP);
684 GMX_MPE_LOG(ev_timestep1);
687 if (rerun_fr.bStep) {
688 step = rerun_fr.step;
689 step_rel = step - ir->init_step;
691 if (rerun_fr.bTime) {
701 bLastStep = (step_rel == ir->nsteps);
702 t = t0 + step*ir->delta_t;
705 if (ir->efep != efepNO)
707 if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
709 state_global->lambda = rerun_fr.lambda;
713 state_global->lambda = lam0 + step*ir->delta_lambda;
715 state->lambda = state_global->lambda;
716 bDoDHDL = do_per_step(step,ir->nstdhdl);
721 update_annealing_target_temp(&(ir->opts),t);
726 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
728 for(i=0; i<state_global->natoms; i++)
730 copy_rvec(rerun_fr.x[i],state_global->x[i]);
734 for(i=0; i<state_global->natoms; i++)
736 copy_rvec(rerun_fr.v[i],state_global->v[i]);
741 for(i=0; i<state_global->natoms; i++)
743 clear_rvec(state_global->v[i]);
747 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
748 " Ekin, temperature and pressure are incorrect,\n"
749 " the virial will be incorrect when constraints are present.\n"
751 bRerunWarnNoV = FALSE;
755 copy_mat(rerun_fr.box,state_global->box);
756 copy_mat(state_global->box,state->box);
758 if (vsite && (Flags & MD_RERUN_VSITE))
760 if (DOMAINDECOMP(cr))
762 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
766 /* Following is necessary because the graph may get out of sync
767 * with the coordinates if we only have every N'th coordinate set
769 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
770 shift_self(graph,state->box,state->x);
772 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
773 top->idef.iparams,top->idef.il,
774 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
777 unshift_self(graph,state->box,state->x);
782 /* Stop Center of Mass motion */
783 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
785 /* Copy back starting coordinates in case we're doing a forcefield scan */
788 for(ii=0; (ii<state->natoms); ii++)
790 copy_rvec(xcopy[ii],state->x[ii]);
791 copy_rvec(vcopy[ii],state->v[ii]);
793 copy_mat(boxcopy,state->box);
798 /* for rerun MD always do Neighbour Searching */
799 bNS = (bFirstStep || ir->nstlist != 0);
804 /* Determine whether or not to do Neighbour Searching and LR */
805 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
807 bNS = (bFirstStep || bExchanged || bNStList ||
808 (ir->nstlist == -1 && nlh.nabnsb > 0));
810 if (bNS && ir->nstlist == -1)
812 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
816 /* check whether we should stop because another simulation has
820 if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
821 (multisim_nsteps != ir->nsteps) )
828 "Stopping simulation %d because another one has finished\n",
832 gs.sig[eglsCHKPT] = 1;
837 /* < 0 means stop at next step, > 0 means stop at next NS step */
838 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
839 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
844 /* Determine whether or not to update the Born radii if doing GB */
845 bBornRadii=bFirstStep;
846 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
851 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
852 do_verbose = bVerbose &&
853 (step % stepout == 0 || bFirstStep || bLastStep);
855 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
863 bMasterState = FALSE;
864 /* Correct the new box if it is too skewed */
865 if (DYNAMIC_BOX(*ir))
867 if (correct_box(fplog,step,state->box,graph))
872 if (DOMAINDECOMP(cr) && bMasterState)
874 dd_collect_state(cr->dd,state,state_global);
878 if (DOMAINDECOMP(cr))
880 /* Repartition the domain decomposition */
881 wallcycle_start(wcycle,ewcDOMDEC);
882 dd_partition_system(fplog,step,cr,
883 bMasterState,nstglobalcomm,
884 state_global,top_global,ir,
885 state,&f,mdatoms,top,fr,
886 vsite,shellfc,constr,
887 nrnb,wcycle,do_verbose);
888 wallcycle_stop(wcycle,ewcDOMDEC);
889 /* If using an iterative integrator, reallocate space to match the decomposition */
893 if (MASTER(cr) && do_log && !bFFscan)
895 print_ebin_header(fplog,step,t,state->lambda);
898 if (ir->efep != efepNO)
900 update_mdatoms(mdatoms,state->lambda);
903 if (bRerunMD && rerun_fr.bV)
906 /* We need the kinetic energy at minus the half step for determining
907 * the full step kinetic energy and possibly for T-coupling.*/
908 /* This may not be quite working correctly yet . . . . */
909 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
910 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
911 constr,NULL,FALSE,state->box,
912 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
913 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
915 clear_mat(force_vir);
917 /* Ionize the atoms if necessary */
920 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
921 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
924 /* Update force field in ffscan program */
927 if (update_forcefield(fplog,
929 mdatoms->nr,state->x,state->box)) {
930 if (gmx_parallel_env_initialized())
938 GMX_MPE_LOG(ev_timestep2);
940 /* We write a checkpoint at this MD step when:
941 * either at an NS step when we signalled through gs,
942 * or at the last step (but not when we do not want confout),
943 * but never at the first step or with rerun.
945 bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
946 (bLastStep && (Flags & MD_CONFOUT))) &&
947 step > ir->init_step && !bRerunMD);
950 gs.set[eglsCHKPT] = 0;
953 /* Determine the energy and pressure:
954 * at nstcalcenergy steps and at energy output steps (set below).
956 bNstEner = do_per_step(step,ir->nstcalcenergy);
959 (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)));
961 /* Do we need global communication ? */
962 bGStat = (bCalcEnerPres || bStopCM ||
963 do_per_step(step,nstglobalcomm) ||
964 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
966 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
968 if (do_ene || do_log)
970 bCalcEnerPres = TRUE;
974 /* these CGLO_ options remain the same throughout the iteration */
975 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
976 (bStopCM ? CGLO_STOPCM : 0) |
977 (bGStat ? CGLO_GSTAT : 0)
980 force_flags = (GMX_FORCE_STATECHANGED |
981 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
982 GMX_FORCE_ALLFORCES |
983 (bNStList ? GMX_FORCE_DOLR : 0) |
985 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
986 (bDoDHDL ? GMX_FORCE_DHDL : 0)
991 /* Now is the time to relax the shells */
992 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
994 bStopCM,top,top_global,
996 state,f,force_vir,mdatoms,
997 nrnb,wcycle,graph,groups,
998 shellfc,fr,bBornRadii,t,mu_tot,
999 state->natoms,&bConverged,vsite,
1010 /* The coordinates (x) are shifted (to get whole molecules)
1012 * This is parallellized as well, and does communication too.
1013 * Check comments in sim_util.c
1016 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1017 state->box,state->x,&state->hist,
1018 f,force_vir,mdatoms,enerd,fcd,
1019 state->lambda,graph,
1020 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1021 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1024 GMX_BARRIER(cr->mpi_comm_mygroup);
1028 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1029 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1032 if (bTCR && bFirstStep)
1034 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1035 fprintf(fplog,"Done init_coupling\n");
1039 if (bVV && !bStartingFromCpt && !bRerunMD)
1040 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1042 if (ir->eI==eiVV && bInitStep)
1044 /* if using velocity verlet with full time step Ekin,
1045 * take the first half step only to compute the
1046 * virial for the first step. From there,
1047 * revert back to the initial coordinates
1048 * so that the input is actually the initial step.
1050 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1052 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1053 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1056 update_coords(fplog,step,ir,mdatoms,state,
1057 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1058 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1059 cr,nrnb,constr,&top->idef);
1063 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1065 /* for iterations, we save these vectors, as we will be self-consistently iterating
1068 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1070 /* save the state */
1071 if (bIterations && iterate.bIterate) {
1072 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1075 bFirstIterate = TRUE;
1076 while (bFirstIterate || (bIterations && iterate.bIterate))
1078 if (bIterations && iterate.bIterate)
1080 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1081 if (bFirstIterate && bTrotter)
1083 /* The first time through, we need a decent first estimate
1084 of veta(t+dt) to compute the constraints. Do
1085 this by computing the box volume part of the
1086 trotter integration at this time. Nothing else
1087 should be changed by this routine here. If
1088 !(first time), we start with the previous value
1091 veta_save = state->veta;
1092 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1093 vetanew = state->veta;
1094 state->veta = veta_save;
1099 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
1102 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1103 &top->idef,shake_vir,NULL,
1104 cr,nrnb,wcycle,upd,constr,
1105 bInitStep,TRUE,bCalcEnerPres,vetanew);
1107 if (!bOK && !bFFscan)
1109 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1114 { /* Need to unshift here if a do_force has been
1115 called in the previous step */
1116 unshift_self(graph,state->box,state->x);
1120 /* if VV, compute the pressure and constraints */
1121 /* For VV2, we strictly only need this if using pressure
1122 * control, but we really would like to have accurate pressures
1124 * Think about ways around this in the future?
1125 * For now, keep this choice in comments.
1127 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1128 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1130 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
1131 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1132 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1133 constr,NULL,FALSE,state->box,
1134 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1137 | (bTemp ? CGLO_TEMPERATURE:0)
1138 | (bPres ? CGLO_PRESSURE : 0)
1139 | (bPres ? CGLO_CONSTRAINT : 0)
1140 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
1141 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1144 /* explanation of above:
1145 a) We compute Ekin at the full time step
1146 if 1) we are using the AveVel Ekin, and it's not the
1147 initial step, or 2) if we are using AveEkin, but need the full
1148 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1149 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1150 EkinAveVel because it's needed for the pressure */
1152 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1157 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
1161 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1166 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1167 state->veta,&vetanew))
1171 bFirstIterate = FALSE;
1174 if (bTrotter && !bInitStep) {
1175 copy_mat(shake_vir,state->svir_prev);
1176 copy_mat(force_vir,state->fvir_prev);
1177 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
1178 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1179 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
1180 enerd->term[F_EKIN] = trace(ekind->ekin);
1183 /* if it's the initial step, we performed this first step just to get the constraint virial */
1184 if (bInitStep && ir->eI==eiVV) {
1185 copy_rvecn(cbuf,state->v,0,state->natoms);
1188 if (fr->bSepDVDL && fplog && do_log)
1190 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1192 enerd->term[F_DHDL_CON] += dvdl;
1194 GMX_MPE_LOG(ev_timestep1);
1197 /* MRS -- now done iterating -- compute the conserved quantity */
1199 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
1202 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
1204 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1206 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1210 /* ######## END FIRST UPDATE STEP ############## */
1211 /* ######## If doing VV, we now have v(dt) ###### */
1213 /* ################## START TRAJECTORY OUTPUT ################# */
1215 /* Now we have the energies and forces corresponding to the
1216 * coordinates at time t. We must output all of this before
1218 * for RerunMD t is read from input trajectory
1220 GMX_MPE_LOG(ev_output_start);
1223 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
1224 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
1225 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
1226 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
1227 if (bCPT) { mdof_flags |= MDOF_CPT; };
1229 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1232 /* Enforce writing positions and velocities at end of run */
1233 mdof_flags |= (MDOF_X | MDOF_V);
1238 fcReportProgress( ir->nsteps, step );
1240 /* sync bCPT and fc record-keeping */
1241 if (bCPT && MASTER(cr))
1242 fcRequestCheckPoint();
1245 if (mdof_flags != 0)
1247 wallcycle_start(wcycle,ewcTRAJ);
1250 if (state->flags & (1<<estLD_RNG))
1252 get_stochd_state(upd,state);
1258 state_global->ekinstate.bUpToDate = FALSE;
1262 update_ekinstate(&state_global->ekinstate,ekind);
1263 state_global->ekinstate.bUpToDate = TRUE;
1265 update_energyhistory(&state_global->enerhist,mdebin);
1268 write_traj(fplog,cr,outf,mdof_flags,top_global,
1269 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
1276 if (bLastStep && step_rel == ir->nsteps &&
1277 (Flags & MD_CONFOUT) && MASTER(cr) &&
1278 !bRerunMD && !bFFscan)
1280 /* x and v have been collected in write_traj,
1281 * because a checkpoint file will always be written
1284 fprintf(stderr,"\nWriting final coordinates.\n");
1285 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
1288 /* Make molecules whole only for confout writing */
1289 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
1291 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
1292 *top_global->name,top_global,
1293 state_global->x,state_global->v,
1294 ir->ePBC,state->box);
1297 wallcycle_stop(wcycle,ewcTRAJ);
1299 GMX_MPE_LOG(ev_output_finish);
1301 /* kludge -- virial is lost with restart for NPT control. Must restart */
1302 if (bStartingFromCpt && bVV)
1304 copy_mat(state->svir_prev,shake_vir);
1305 copy_mat(state->fvir_prev,force_vir);
1307 /* ################## END TRAJECTORY OUTPUT ################ */
1309 /* Determine the wallclock run time up till now */
1310 run_time = gmx_gettime() - (double)runtime->real;
1312 /* Check whether everything is still allright */
1313 if (((int)gmx_get_stop_condition() > handled_stop_condition)
1319 /* this is just make gs.sig compatible with the hack
1320 of sending signals around by MPI_Reduce with together with
1322 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
1323 gs.sig[eglsSTOPCOND]=1;
1324 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
1325 gs.sig[eglsSTOPCOND]=-1;
1326 /* < 0 means stop at next step, > 0 means stop at next NS step */
1330 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1331 gmx_get_signal_name(),
1332 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1336 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1337 gmx_get_signal_name(),
1338 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1340 handled_stop_condition=(int)gmx_get_stop_condition();
1342 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1343 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
1344 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
1346 /* Signal to terminate the run */
1347 gs.sig[eglsSTOPCOND] = 1;
1350 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1352 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1355 if (bResetCountersHalfMaxH && MASTER(cr) &&
1356 run_time > max_hours*60.0*60.0*0.495)
1358 gs.sig[eglsRESETCOUNTERS] = 1;
1361 if (ir->nstlist == -1 && !bRerunMD)
1363 /* When bGStatEveryStep=FALSE, global_stat is only called
1364 * when we check the atom displacements, not at NS steps.
1365 * This means that also the bonded interaction count check is not
1366 * performed immediately after NS. Therefore a few MD steps could
1367 * be performed with missing interactions.
1368 * But wrong energies are never written to file,
1369 * since energies are only written after global_stat
1372 if (step >= nlh.step_nscheck)
1374 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
1375 nlh.scale_tot,state->x);
1379 /* This is not necessarily true,
1380 * but step_nscheck is determined quite conservatively.
1386 /* In parallel we only have to check for checkpointing in steps
1387 * where we do global communication,
1388 * otherwise the other nodes don't know.
1390 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1393 run_time >= nchkpt*cpt_period*60.0)) &&
1394 gs.set[eglsCHKPT] == 0)
1396 gs.sig[eglsCHKPT] = 1;
1401 gmx_iterate_init(&iterate,bIterations);
1404 /* for iterations, we save these vectors, as we will be redoing the calculations */
1405 if (bIterations && iterate.bIterate)
1407 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1409 bFirstIterate = TRUE;
1410 while (bFirstIterate || (bIterations && iterate.bIterate))
1412 /* We now restore these vectors to redo the calculation with improved extended variables */
1415 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1418 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1419 so scroll down for that logic */
1421 /* ######### START SECOND UPDATE STEP ################# */
1422 GMX_MPE_LOG(ev_update_start);
1423 /* Box is changed in update() when we do pressure coupling,
1424 * but we should still use the old box for energy corrections and when
1425 * writing it to the energy file, so it matches the trajectory files for
1426 * the same timestep above. Make a copy in a separate array.
1428 copy_mat(state->box,lastbox);
1431 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
1433 wallcycle_start(wcycle,ewcUPDATE);
1435 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1438 if (bIterations && iterate.bIterate)
1446 /* we use a new value of scalevir to converge the iterations faster */
1447 scalevir = tracevir/trace(shake_vir);
1449 msmul(shake_vir,scalevir,shake_vir);
1450 m_add(force_vir,shake_vir,total_vir);
1451 clear_mat(shake_vir);
1453 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
1454 /* We can only do Berendsen coupling after we have summed
1455 * the kinetic energy or virial. Since the happens
1456 * in global_state after update, we should only do it at
1457 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1462 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1463 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
1469 /* velocity half-step update */
1470 update_coords(fplog,step,ir,mdatoms,state,f,
1471 fr->bTwinRange && bNStList,fr->f_twin,fcd,
1472 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
1473 cr,nrnb,constr,&top->idef);
1476 /* Above, initialize just copies ekinh into ekin,
1477 * it doesn't copy position (for VV),
1478 * and entire integrator for MD.
1483 copy_rvecn(state->x,cbuf,0,state->natoms);
1486 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1487 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1488 wallcycle_stop(wcycle,ewcUPDATE);
1490 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1491 &top->idef,shake_vir,force_vir,
1492 cr,nrnb,wcycle,upd,constr,
1493 bInitStep,FALSE,bCalcEnerPres,state->veta);
1497 /* erase F_EKIN and F_TEMP here? */
1498 /* just compute the kinetic energy at the half step to perform a trotter step */
1499 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1500 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1501 constr,NULL,FALSE,lastbox,
1502 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1503 cglo_flags | CGLO_TEMPERATURE
1505 wallcycle_start(wcycle,ewcUPDATE);
1506 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
1507 /* now we know the scaling, we can compute the positions again again */
1508 copy_rvecn(cbuf,state->x,0,state->natoms);
1510 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1511 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1512 wallcycle_stop(wcycle,ewcUPDATE);
1514 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1515 /* are the small terms in the shake_vir here due
1516 * to numerical errors, or are they important
1517 * physically? I'm thinking they are just errors, but not completely sure.
1518 * For now, will call without actually constraining, constr=NULL*/
1519 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1520 &top->idef,tmp_vir,force_vir,
1521 cr,nrnb,wcycle,upd,NULL,
1522 bInitStep,FALSE,bCalcEnerPres,
1525 if (!bOK && !bFFscan)
1527 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1530 if (fr->bSepDVDL && fplog && do_log)
1532 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1534 enerd->term[F_DHDL_CON] += dvdl;
1538 /* Need to unshift here */
1539 unshift_self(graph,state->box,state->x);
1542 GMX_BARRIER(cr->mpi_comm_mygroup);
1543 GMX_MPE_LOG(ev_update_finish);
1547 wallcycle_start(wcycle,ewcVSITECONSTR);
1550 shift_self(graph,state->box,state->x);
1552 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1553 top->idef.iparams,top->idef.il,
1554 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1558 unshift_self(graph,state->box,state->x);
1560 wallcycle_stop(wcycle,ewcVSITECONSTR);
1563 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1564 if (ir->nstlist == -1 && bFirstIterate)
1566 gs.sig[eglsNABNSB] = nlh.nabnsb;
1568 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1569 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1571 bFirstIterate ? &gs : NULL,
1572 (step_rel % gs.nstms == 0) &&
1573 (multisim_nsteps<0 || (step_rel<multisim_nsteps)),
1575 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1577 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
1578 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1579 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1580 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
1581 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1584 if (ir->nstlist == -1 && bFirstIterate)
1586 nlh.nabnsb = gs.set[eglsNABNSB];
1587 gs.set[eglsNABNSB] = 0;
1589 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1590 /* ############# END CALC EKIN AND PRESSURE ################# */
1592 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1593 the virial that should probably be addressed eventually. state->veta has better properies,
1594 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1595 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1598 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1599 trace(shake_vir),&tracevir))
1603 bFirstIterate = FALSE;
1606 update_box(fplog,step,ir,mdatoms,state,graph,f,
1607 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
1609 /* ################# END UPDATE STEP 2 ################# */
1610 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1612 /* The coordinates (x) were unshifted in update */
1613 if (bFFscan && (shellfc==NULL || bConverged))
1615 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
1617 &(top_global->mols),mdatoms->massT,pres))
1619 if (gmx_parallel_env_initialized())
1623 fprintf(stderr,"\n");
1629 /* We will not sum ekinh_old,
1630 * so signal that we still have to do it.
1632 bSumEkinhOld = TRUE;
1637 /* Only do GCT when the relaxation of shells (minimization) has converged,
1638 * otherwise we might be coupling to bogus energies.
1639 * In parallel we must always do this, because the other sims might
1643 /* Since this is called with the new coordinates state->x, I assume
1644 * we want the new box state->box too. / EL 20040121
1646 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
1648 mdatoms,&(top->idef),mu_aver,
1649 top_global->mols.nr,cr,
1650 state->box,total_vir,pres,
1651 mu_tot,state->x,f,bConverged);
1655 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1657 /* sum up the foreign energy and dhdl terms */
1658 sum_dhdl(enerd,state->lambda,ir);
1660 /* use the directly determined last velocity, not actually the averaged half steps */
1661 if (bTrotter && ir->eI==eiVV)
1663 enerd->term[F_EKIN] = last_ekin;
1665 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1669 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1673 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
1675 /* Check for excessively large energies */
1679 real etot_max = 1e200;
1681 real etot_max = 1e30;
1683 if (fabs(enerd->term[F_ETOT]) > etot_max)
1685 fprintf(stderr,"Energy too large (%g), giving up\n",
1686 enerd->term[F_ETOT]);
1689 /* ######### END PREPARING EDR OUTPUT ########### */
1691 /* Time for performance */
1692 if (((step % stepout) == 0) || bLastStep)
1694 runtime_upd_proc(runtime);
1700 gmx_bool do_dr,do_or;
1702 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
1706 upd_mdebin(mdebin,bDoDHDL, TRUE,
1707 t,mdatoms->tmass,enerd,state,lastbox,
1708 shake_vir,force_vir,total_vir,pres,
1709 ekind,mu_tot,constr);
1713 upd_mdebin_step(mdebin);
1716 do_dr = do_per_step(step,ir->nstdisreout);
1717 do_or = do_per_step(step,ir->nstorireout);
1719 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
1721 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
1723 if (ir->ePull != epullNO)
1725 pull_print_output(ir->pull,step,t);
1728 if (do_per_step(step,ir->nstlog))
1730 if(fflush(fplog) != 0)
1732 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
1738 /* Remaining runtime */
1739 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
1743 fprintf(stderr,"\n");
1745 print_time(stderr,runtime,step,ir,cr);
1748 /* Replica exchange */
1750 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
1751 do_per_step(step,repl_ex_nst))
1753 bExchanged = replica_exchange(fplog,cr,repl_ex,
1754 state_global,enerd->term,
1757 if (bExchanged && DOMAINDECOMP(cr))
1759 dd_partition_system(fplog,step,cr,TRUE,1,
1760 state_global,top_global,ir,
1761 state,&f,mdatoms,top,fr,
1762 vsite,shellfc,constr,
1769 bStartingFromCpt = FALSE;
1771 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1772 /* With all integrators, except VV, we need to retain the pressure
1773 * at the current step for coupling at the next step.
1775 if ((state->flags & (1<<estPRES_PREV)) &&
1777 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1779 /* Store the pressure in t_state for pressure coupling
1780 * at the next MD step.
1782 copy_mat(pres,state->pres_prev);
1785 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1787 if ( (membed!=NULL) && (!bLastStep) )
1788 rescale_membed(step_rel,membed,state_global->x);
1794 /* read next frame from input trajectory */
1795 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
1800 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
1804 if (!bRerunMD || !rerun_fr.bStep)
1806 /* increase the MD step number */
1811 cycles = wallcycle_stop(wcycle,ewcSTEP);
1812 if (DOMAINDECOMP(cr) && wcycle)
1814 dd_cycles_add(cr->dd,cycles,ddCyclStep);
1817 if (step_rel == wcycle_get_reset_counters(wcycle) ||
1818 gs.set[eglsRESETCOUNTERS] != 0)
1820 /* Reset all the counters related to performance over the run */
1821 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
1822 wcycle_set_reset_counters(wcycle,-1);
1823 /* Correct max_hours for the elapsed time */
1824 max_hours -= run_time/(60.0*60.0);
1825 bResetCountersHalfMaxH = FALSE;
1826 gs.set[eglsRESETCOUNTERS] = 0;
1830 /* End of main MD loop */
1834 runtime_end(runtime);
1836 if (bRerunMD && MASTER(cr))
1841 if (!(cr->duty & DUTY_PME))
1843 /* Tell the PME only node to finish */
1849 if (ir->nstcalcenergy > 0 && !bRerunMD)
1851 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
1852 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
1860 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
1862 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)));
1863 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
1866 if (shellfc && fplog)
1868 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
1869 (nconverged*100.0)/step_rel);
1870 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
1874 if (repl_ex_nst > 0 && MASTER(cr))
1876 print_replica_exchange_statistics(fplog,repl_ex);
1879 runtime->nsteps_done = step_rel;