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
77 #include "compute_io.h"
79 #include "checkpoint.h"
80 #include "mtop_util.h"
81 #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,gmx_membed_t membed,
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);
681 if (rerun_fr.bStep) {
682 step = rerun_fr.step;
683 step_rel = step - ir->init_step;
685 if (rerun_fr.bTime) {
695 bLastStep = (step_rel == ir->nsteps);
696 t = t0 + step*ir->delta_t;
699 if (ir->efep != efepNO)
701 if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
703 state_global->lambda = rerun_fr.lambda;
707 state_global->lambda = lam0 + step*ir->delta_lambda;
709 state->lambda = state_global->lambda;
710 bDoDHDL = do_per_step(step,ir->nstdhdl);
715 update_annealing_target_temp(&(ir->opts),t);
720 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
722 for(i=0; i<state_global->natoms; i++)
724 copy_rvec(rerun_fr.x[i],state_global->x[i]);
728 for(i=0; i<state_global->natoms; i++)
730 copy_rvec(rerun_fr.v[i],state_global->v[i]);
735 for(i=0; i<state_global->natoms; i++)
737 clear_rvec(state_global->v[i]);
741 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
742 " Ekin, temperature and pressure are incorrect,\n"
743 " the virial will be incorrect when constraints are present.\n"
745 bRerunWarnNoV = FALSE;
749 copy_mat(rerun_fr.box,state_global->box);
750 copy_mat(state_global->box,state->box);
752 if (vsite && (Flags & MD_RERUN_VSITE))
754 if (DOMAINDECOMP(cr))
756 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
760 /* Following is necessary because the graph may get out of sync
761 * with the coordinates if we only have every N'th coordinate set
763 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
764 shift_self(graph,state->box,state->x);
766 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
767 top->idef.iparams,top->idef.il,
768 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
771 unshift_self(graph,state->box,state->x);
776 /* Stop Center of Mass motion */
777 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
779 /* Copy back starting coordinates in case we're doing a forcefield scan */
782 for(ii=0; (ii<state->natoms); ii++)
784 copy_rvec(xcopy[ii],state->x[ii]);
785 copy_rvec(vcopy[ii],state->v[ii]);
787 copy_mat(boxcopy,state->box);
792 /* for rerun MD always do Neighbour Searching */
793 bNS = (bFirstStep || ir->nstlist != 0);
798 /* Determine whether or not to do Neighbour Searching and LR */
799 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
801 bNS = (bFirstStep || bExchanged || bNStList ||
802 (ir->nstlist == -1 && nlh.nabnsb > 0));
804 if (bNS && ir->nstlist == -1)
806 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
810 /* check whether we should stop because another simulation has
814 if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
815 (multisim_nsteps != ir->nsteps) )
822 "Stopping simulation %d because another one has finished\n",
826 gs.sig[eglsCHKPT] = 1;
831 /* < 0 means stop at next step, > 0 means stop at next NS step */
832 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
833 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
838 /* Determine whether or not to update the Born radii if doing GB */
839 bBornRadii=bFirstStep;
840 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
845 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
846 do_verbose = bVerbose &&
847 (step % stepout == 0 || bFirstStep || bLastStep);
849 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
857 bMasterState = FALSE;
858 /* Correct the new box if it is too skewed */
859 if (DYNAMIC_BOX(*ir))
861 if (correct_box(fplog,step,state->box,graph))
866 if (DOMAINDECOMP(cr) && bMasterState)
868 dd_collect_state(cr->dd,state,state_global);
872 if (DOMAINDECOMP(cr))
874 /* Repartition the domain decomposition */
875 wallcycle_start(wcycle,ewcDOMDEC);
876 dd_partition_system(fplog,step,cr,
877 bMasterState,nstglobalcomm,
878 state_global,top_global,ir,
879 state,&f,mdatoms,top,fr,
880 vsite,shellfc,constr,
881 nrnb,wcycle,do_verbose);
882 wallcycle_stop(wcycle,ewcDOMDEC);
883 /* If using an iterative integrator, reallocate space to match the decomposition */
887 if (MASTER(cr) && do_log && !bFFscan)
889 print_ebin_header(fplog,step,t,state->lambda);
892 if (ir->efep != efepNO)
894 update_mdatoms(mdatoms,state->lambda);
897 if (bRerunMD && rerun_fr.bV)
900 /* We need the kinetic energy at minus the half step for determining
901 * the full step kinetic energy and possibly for T-coupling.*/
902 /* This may not be quite working correctly yet . . . . */
903 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
904 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
905 constr,NULL,FALSE,state->box,
906 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
907 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
909 clear_mat(force_vir);
911 /* Ionize the atoms if necessary */
914 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
915 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
918 /* Update force field in ffscan program */
921 if (update_forcefield(fplog,
923 mdatoms->nr,state->x,state->box)) {
924 if (gmx_parallel_env_initialized())
932 /* We write a checkpoint at this MD step when:
933 * either at an NS step when we signalled through gs,
934 * or at the last step (but not when we do not want confout),
935 * but never at the first step or with rerun.
937 bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
938 (bLastStep && (Flags & MD_CONFOUT))) &&
939 step > ir->init_step && !bRerunMD);
942 gs.set[eglsCHKPT] = 0;
945 /* Determine the energy and pressure:
946 * at nstcalcenergy steps and at energy output steps (set below).
948 bNstEner = do_per_step(step,ir->nstcalcenergy);
951 (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)));
953 /* Do we need global communication ? */
954 bGStat = (bCalcEnerPres || bStopCM ||
955 do_per_step(step,nstglobalcomm) ||
956 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
958 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
960 if (do_ene || do_log)
962 bCalcEnerPres = TRUE;
966 /* these CGLO_ options remain the same throughout the iteration */
967 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
968 (bGStat ? CGLO_GSTAT : 0)
971 force_flags = (GMX_FORCE_STATECHANGED |
972 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
973 GMX_FORCE_ALLFORCES |
974 (bNStList ? GMX_FORCE_DOLR : 0) |
976 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
977 (bDoDHDL ? GMX_FORCE_DHDL : 0)
982 /* Now is the time to relax the shells */
983 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
985 bStopCM,top,top_global,
987 state,f,force_vir,mdatoms,
988 nrnb,wcycle,graph,groups,
989 shellfc,fr,bBornRadii,t,mu_tot,
990 state->natoms,&bConverged,vsite,
1001 /* The coordinates (x) are shifted (to get whole molecules)
1003 * This is parallellized as well, and does communication too.
1004 * Check comments in sim_util.c
1007 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1008 state->box,state->x,&state->hist,
1009 f,force_vir,mdatoms,enerd,fcd,
1010 state->lambda,graph,
1011 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1012 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1017 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1018 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1021 if (bTCR && bFirstStep)
1023 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1024 fprintf(fplog,"Done init_coupling\n");
1028 if (bVV && !bStartingFromCpt && !bRerunMD)
1029 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1031 if (ir->eI==eiVV && bInitStep)
1033 /* if using velocity verlet with full time step Ekin,
1034 * take the first half step only to compute the
1035 * virial for the first step. From there,
1036 * revert back to the initial coordinates
1037 * so that the input is actually the initial step.
1039 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1041 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1042 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1045 update_coords(fplog,step,ir,mdatoms,state,
1046 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1047 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1048 cr,nrnb,constr,&top->idef);
1052 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1054 /* for iterations, we save these vectors, as we will be self-consistently iterating
1057 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1059 /* save the state */
1060 if (bIterations && iterate.bIterate) {
1061 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1064 bFirstIterate = TRUE;
1065 while (bFirstIterate || (bIterations && iterate.bIterate))
1067 if (bIterations && iterate.bIterate)
1069 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1070 if (bFirstIterate && bTrotter)
1072 /* The first time through, we need a decent first estimate
1073 of veta(t+dt) to compute the constraints. Do
1074 this by computing the box volume part of the
1075 trotter integration at this time. Nothing else
1076 should be changed by this routine here. If
1077 !(first time), we start with the previous value
1080 veta_save = state->veta;
1081 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1082 vetanew = state->veta;
1083 state->veta = veta_save;
1088 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
1091 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1092 &top->idef,shake_vir,NULL,
1093 cr,nrnb,wcycle,upd,constr,
1094 bInitStep,TRUE,bCalcEnerPres,vetanew);
1096 if (!bOK && !bFFscan)
1098 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1103 { /* Need to unshift here if a do_force has been
1104 called in the previous step */
1105 unshift_self(graph,state->box,state->x);
1109 /* if VV, compute the pressure and constraints */
1110 /* For VV2, we strictly only need this if using pressure
1111 * control, but we really would like to have accurate pressures
1113 * Think about ways around this in the future?
1114 * For now, keep this choice in comments.
1116 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1117 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1119 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
1120 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1121 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1122 constr,NULL,FALSE,state->box,
1123 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1126 | (bStopCM ? CGLO_STOPCM : 0)
1127 | (bTemp ? CGLO_TEMPERATURE:0)
1128 | (bPres ? CGLO_PRESSURE : 0)
1129 | (bPres ? CGLO_CONSTRAINT : 0)
1130 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
1131 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1134 /* explanation of above:
1135 a) We compute Ekin at the full time step
1136 if 1) we are using the AveVel Ekin, and it's not the
1137 initial step, or 2) if we are using AveEkin, but need the full
1138 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1139 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1140 EkinAveVel because it's needed for the pressure */
1142 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1147 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
1151 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1156 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1157 state->veta,&vetanew))
1161 bFirstIterate = FALSE;
1164 if (bTrotter && !bInitStep) {
1165 copy_mat(shake_vir,state->svir_prev);
1166 copy_mat(force_vir,state->fvir_prev);
1167 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
1168 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1169 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
1170 enerd->term[F_EKIN] = trace(ekind->ekin);
1173 /* if it's the initial step, we performed this first step just to get the constraint virial */
1174 if (bInitStep && ir->eI==eiVV) {
1175 copy_rvecn(cbuf,state->v,0,state->natoms);
1178 if (fr->bSepDVDL && fplog && do_log)
1180 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1182 enerd->term[F_DHDL_CON] += dvdl;
1185 /* MRS -- now done iterating -- compute the conserved quantity */
1187 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
1190 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
1192 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1194 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1198 /* ######## END FIRST UPDATE STEP ############## */
1199 /* ######## If doing VV, we now have v(dt) ###### */
1201 /* ################## START TRAJECTORY OUTPUT ################# */
1203 /* Now we have the energies and forces corresponding to the
1204 * coordinates at time t. We must output all of this before
1206 * for RerunMD t is read from input trajectory
1209 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
1210 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
1211 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
1212 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
1213 if (bCPT) { mdof_flags |= MDOF_CPT; };
1215 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1218 /* Enforce writing positions and velocities at end of run */
1219 mdof_flags |= (MDOF_X | MDOF_V);
1224 fcReportProgress( ir->nsteps, step );
1226 /* sync bCPT and fc record-keeping */
1227 if (bCPT && MASTER(cr))
1228 fcRequestCheckPoint();
1231 if (mdof_flags != 0)
1233 wallcycle_start(wcycle,ewcTRAJ);
1236 if (state->flags & (1<<estLD_RNG))
1238 get_stochd_state(upd,state);
1244 state_global->ekinstate.bUpToDate = FALSE;
1248 update_ekinstate(&state_global->ekinstate,ekind);
1249 state_global->ekinstate.bUpToDate = TRUE;
1251 update_energyhistory(&state_global->enerhist,mdebin);
1254 write_traj(fplog,cr,outf,mdof_flags,top_global,
1255 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
1262 if (bLastStep && step_rel == ir->nsteps &&
1263 (Flags & MD_CONFOUT) && MASTER(cr) &&
1264 !bRerunMD && !bFFscan)
1266 /* x and v have been collected in write_traj,
1267 * because a checkpoint file will always be written
1270 fprintf(stderr,"\nWriting final coordinates.\n");
1271 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
1274 /* Make molecules whole only for confout writing */
1275 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
1277 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
1278 *top_global->name,top_global,
1279 state_global->x,state_global->v,
1280 ir->ePBC,state->box);
1283 wallcycle_stop(wcycle,ewcTRAJ);
1286 /* kludge -- virial is lost with restart for NPT control. Must restart */
1287 if (bStartingFromCpt && bVV)
1289 copy_mat(state->svir_prev,shake_vir);
1290 copy_mat(state->fvir_prev,force_vir);
1292 /* ################## END TRAJECTORY OUTPUT ################ */
1294 /* Determine the wallclock run time up till now */
1295 run_time = gmx_gettime() - (double)runtime->real;
1297 /* Check whether everything is still allright */
1298 if (((int)gmx_get_stop_condition() > handled_stop_condition)
1299 #ifdef GMX_THREAD_MPI
1304 /* this is just make gs.sig compatible with the hack
1305 of sending signals around by MPI_Reduce with together with
1307 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
1308 gs.sig[eglsSTOPCOND]=1;
1309 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
1310 gs.sig[eglsSTOPCOND]=-1;
1311 /* < 0 means stop at next step, > 0 means stop at next NS step */
1315 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1316 gmx_get_signal_name(),
1317 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1321 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1322 gmx_get_signal_name(),
1323 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1325 handled_stop_condition=(int)gmx_get_stop_condition();
1327 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1328 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
1329 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
1331 /* Signal to terminate the run */
1332 gs.sig[eglsSTOPCOND] = 1;
1335 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1337 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1340 if (bResetCountersHalfMaxH && MASTER(cr) &&
1341 run_time > max_hours*60.0*60.0*0.495)
1343 gs.sig[eglsRESETCOUNTERS] = 1;
1346 if (ir->nstlist == -1 && !bRerunMD)
1348 /* When bGStatEveryStep=FALSE, global_stat is only called
1349 * when we check the atom displacements, not at NS steps.
1350 * This means that also the bonded interaction count check is not
1351 * performed immediately after NS. Therefore a few MD steps could
1352 * be performed with missing interactions.
1353 * But wrong energies are never written to file,
1354 * since energies are only written after global_stat
1357 if (step >= nlh.step_nscheck)
1359 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
1360 nlh.scale_tot,state->x);
1364 /* This is not necessarily true,
1365 * but step_nscheck is determined quite conservatively.
1371 /* In parallel we only have to check for checkpointing in steps
1372 * where we do global communication,
1373 * otherwise the other nodes don't know.
1375 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1378 run_time >= nchkpt*cpt_period*60.0)) &&
1379 gs.set[eglsCHKPT] == 0)
1381 gs.sig[eglsCHKPT] = 1;
1386 gmx_iterate_init(&iterate,bIterations);
1389 /* for iterations, we save these vectors, as we will be redoing the calculations */
1390 if (bIterations && iterate.bIterate)
1392 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1394 bFirstIterate = TRUE;
1395 while (bFirstIterate || (bIterations && iterate.bIterate))
1397 /* We now restore these vectors to redo the calculation with improved extended variables */
1400 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1403 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1404 so scroll down for that logic */
1406 /* ######### START SECOND UPDATE STEP ################# */
1407 /* Box is changed in update() when we do pressure coupling,
1408 * but we should still use the old box for energy corrections and when
1409 * writing it to the energy file, so it matches the trajectory files for
1410 * the same timestep above. Make a copy in a separate array.
1412 copy_mat(state->box,lastbox);
1415 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
1417 wallcycle_start(wcycle,ewcUPDATE);
1419 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1422 if (bIterations && iterate.bIterate)
1430 /* we use a new value of scalevir to converge the iterations faster */
1431 scalevir = tracevir/trace(shake_vir);
1433 msmul(shake_vir,scalevir,shake_vir);
1434 m_add(force_vir,shake_vir,total_vir);
1435 clear_mat(shake_vir);
1437 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
1438 /* We can only do Berendsen coupling after we have summed
1439 * the kinetic energy or virial. Since the happens
1440 * in global_state after update, we should only do it at
1441 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1446 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1447 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
1453 /* velocity half-step update */
1454 update_coords(fplog,step,ir,mdatoms,state,f,
1455 fr->bTwinRange && bNStList,fr->f_twin,fcd,
1456 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
1457 cr,nrnb,constr,&top->idef);
1460 /* Above, initialize just copies ekinh into ekin,
1461 * it doesn't copy position (for VV),
1462 * and entire integrator for MD.
1467 copy_rvecn(state->x,cbuf,0,state->natoms);
1470 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1471 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1472 wallcycle_stop(wcycle,ewcUPDATE);
1474 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1475 &top->idef,shake_vir,force_vir,
1476 cr,nrnb,wcycle,upd,constr,
1477 bInitStep,FALSE,bCalcEnerPres,state->veta);
1481 /* erase F_EKIN and F_TEMP here? */
1482 /* just compute the kinetic energy at the half step to perform a trotter step */
1483 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1484 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1485 constr,NULL,FALSE,lastbox,
1486 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1487 cglo_flags | CGLO_TEMPERATURE
1489 wallcycle_start(wcycle,ewcUPDATE);
1490 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
1491 /* now we know the scaling, we can compute the positions again again */
1492 copy_rvecn(cbuf,state->x,0,state->natoms);
1494 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1495 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1496 wallcycle_stop(wcycle,ewcUPDATE);
1498 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1499 /* are the small terms in the shake_vir here due
1500 * to numerical errors, or are they important
1501 * physically? I'm thinking they are just errors, but not completely sure.
1502 * For now, will call without actually constraining, constr=NULL*/
1503 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1504 &top->idef,tmp_vir,force_vir,
1505 cr,nrnb,wcycle,upd,NULL,
1506 bInitStep,FALSE,bCalcEnerPres,
1509 if (!bOK && !bFFscan)
1511 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1514 if (fr->bSepDVDL && fplog && do_log)
1516 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1518 enerd->term[F_DHDL_CON] += dvdl;
1522 /* Need to unshift here */
1523 unshift_self(graph,state->box,state->x);
1528 wallcycle_start(wcycle,ewcVSITECONSTR);
1531 shift_self(graph,state->box,state->x);
1533 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1534 top->idef.iparams,top->idef.il,
1535 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1539 unshift_self(graph,state->box,state->x);
1541 wallcycle_stop(wcycle,ewcVSITECONSTR);
1544 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1545 if (ir->nstlist == -1 && bFirstIterate)
1547 gs.sig[eglsNABNSB] = nlh.nabnsb;
1549 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1550 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1552 bFirstIterate ? &gs : NULL,
1553 (step_rel % gs.nstms == 0) &&
1554 (multisim_nsteps<0 || (step_rel<multisim_nsteps)),
1556 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1558 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
1559 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1560 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1561 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1562 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
1563 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1566 if (ir->nstlist == -1 && bFirstIterate)
1568 nlh.nabnsb = gs.set[eglsNABNSB];
1569 gs.set[eglsNABNSB] = 0;
1571 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1572 /* ############# END CALC EKIN AND PRESSURE ################# */
1574 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1575 the virial that should probably be addressed eventually. state->veta has better properies,
1576 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1577 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1580 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1581 trace(shake_vir),&tracevir))
1585 bFirstIterate = FALSE;
1588 update_box(fplog,step,ir,mdatoms,state,graph,f,
1589 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
1591 /* ################# END UPDATE STEP 2 ################# */
1592 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1594 /* The coordinates (x) were unshifted in update */
1595 if (bFFscan && (shellfc==NULL || bConverged))
1597 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
1599 &(top_global->mols),mdatoms->massT,pres))
1601 if (gmx_parallel_env_initialized())
1605 fprintf(stderr,"\n");
1611 /* We will not sum ekinh_old,
1612 * so signal that we still have to do it.
1614 bSumEkinhOld = TRUE;
1619 /* Only do GCT when the relaxation of shells (minimization) has converged,
1620 * otherwise we might be coupling to bogus energies.
1621 * In parallel we must always do this, because the other sims might
1625 /* Since this is called with the new coordinates state->x, I assume
1626 * we want the new box state->box too. / EL 20040121
1628 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
1630 mdatoms,&(top->idef),mu_aver,
1631 top_global->mols.nr,cr,
1632 state->box,total_vir,pres,
1633 mu_tot,state->x,f,bConverged);
1637 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1639 /* sum up the foreign energy and dhdl terms */
1640 sum_dhdl(enerd,state->lambda,ir);
1642 /* use the directly determined last velocity, not actually the averaged half steps */
1643 if (bTrotter && ir->eI==eiVV)
1645 enerd->term[F_EKIN] = last_ekin;
1647 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1651 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1655 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
1657 /* Check for excessively large energies */
1661 real etot_max = 1e200;
1663 real etot_max = 1e30;
1665 if (fabs(enerd->term[F_ETOT]) > etot_max)
1667 fprintf(stderr,"Energy too large (%g), giving up\n",
1668 enerd->term[F_ETOT]);
1671 /* ######### END PREPARING EDR OUTPUT ########### */
1673 /* Time for performance */
1674 if (((step % stepout) == 0) || bLastStep)
1676 runtime_upd_proc(runtime);
1682 gmx_bool do_dr,do_or;
1684 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
1688 upd_mdebin(mdebin,bDoDHDL, TRUE,
1689 t,mdatoms->tmass,enerd,state,lastbox,
1690 shake_vir,force_vir,total_vir,pres,
1691 ekind,mu_tot,constr);
1695 upd_mdebin_step(mdebin);
1698 do_dr = do_per_step(step,ir->nstdisreout);
1699 do_or = do_per_step(step,ir->nstorireout);
1701 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
1703 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
1705 if (ir->ePull != epullNO)
1707 pull_print_output(ir->pull,step,t);
1710 if (do_per_step(step,ir->nstlog))
1712 if(fflush(fplog) != 0)
1714 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
1720 /* Remaining runtime */
1721 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
1725 fprintf(stderr,"\n");
1727 print_time(stderr,runtime,step,ir,cr);
1730 /* Replica exchange */
1732 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
1733 do_per_step(step,repl_ex_nst))
1735 bExchanged = replica_exchange(fplog,cr,repl_ex,
1736 state_global,enerd->term,
1739 if (bExchanged && DOMAINDECOMP(cr))
1741 dd_partition_system(fplog,step,cr,TRUE,1,
1742 state_global,top_global,ir,
1743 state,&f,mdatoms,top,fr,
1744 vsite,shellfc,constr,
1751 bStartingFromCpt = FALSE;
1753 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1754 /* With all integrators, except VV, we need to retain the pressure
1755 * at the current step for coupling at the next step.
1757 if ((state->flags & (1<<estPRES_PREV)) &&
1759 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1761 /* Store the pressure in t_state for pressure coupling
1762 * at the next MD step.
1764 copy_mat(pres,state->pres_prev);
1767 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1769 if ( (membed!=NULL) && (!bLastStep) )
1771 rescale_membed(step_rel,membed,state_global->x);
1778 /* read next frame from input trajectory */
1779 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
1784 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
1788 if (!bRerunMD || !rerun_fr.bStep)
1790 /* increase the MD step number */
1795 cycles = wallcycle_stop(wcycle,ewcSTEP);
1796 if (DOMAINDECOMP(cr) && wcycle)
1798 dd_cycles_add(cr->dd,cycles,ddCyclStep);
1801 if (step_rel == wcycle_get_reset_counters(wcycle) ||
1802 gs.set[eglsRESETCOUNTERS] != 0)
1804 /* Reset all the counters related to performance over the run */
1805 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
1806 wcycle_set_reset_counters(wcycle,-1);
1807 /* Correct max_hours for the elapsed time */
1808 max_hours -= run_time/(60.0*60.0);
1809 bResetCountersHalfMaxH = FALSE;
1810 gs.set[eglsRESETCOUNTERS] = 0;
1814 /* End of main MD loop */
1818 runtime_end(runtime);
1820 if (bRerunMD && MASTER(cr))
1825 if (!(cr->duty & DUTY_PME))
1827 /* Tell the PME only node to finish */
1833 if (ir->nstcalcenergy > 0 && !bRerunMD)
1835 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
1836 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
1844 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
1846 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)));
1847 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
1850 if (shellfc && fplog)
1852 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
1853 (nconverged*100.0)/step_rel);
1854 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
1858 if (repl_ex_nst > 0 && MASTER(cr))
1860 print_replica_exchange_statistics(fplog,repl_ex);
1863 runtime->nsteps_done = step_rel;