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
76 #include "compute_io.h"
78 #include "checkpoint.h"
79 #include "mtop_util.h"
80 #include "sighandler.h"
96 double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
97 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
99 gmx_vsite_t *vsite,gmx_constr_t constr,
100 int stepout,t_inputrec *ir,
101 gmx_mtop_t *top_global,
103 t_state *state_global,
105 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
106 gmx_edsam_t ed,t_forcerec *fr,
107 int repl_ex_nst,int repl_ex_nex,int repl_ex_seed,gmx_membed_t membed,
108 real cpt_period,real max_hours,
109 const char *deviceOptions,
111 gmx_runtime_t *runtime)
114 gmx_large_int_t step,step_rel;
116 double t,t0,lam0[efptNR];
117 gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres,bEnergyHere;
118 gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
119 bFirstStep,bStateFromCP,bStateFromTPX,bInitStep,bLastStep,
120 bBornRadii,bStartingFromCpt;
121 gmx_bool bDoDHDL=FALSE,bDoFEP=FALSE,bDoExpanded=FALSE;
122 gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
123 bForceUpdate=FALSE,bCPT;
125 gmx_bool bMasterState;
126 int force_flags,cglo_flags;
127 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
132 t_state *bufstate=NULL;
133 matrix *scale_tot,pcoupl_mu,M,ebox;
136 gmx_repl_ex_t repl_ex=NULL;
139 t_mdebin *mdebin=NULL;
140 df_history_t df_history;
145 gmx_enerdata_t *enerd;
147 gmx_global_stat_t gstat;
148 gmx_update_t upd=NULL;
151 gmx_rng_t mcrng=NULL;
153 gmx_groups_t *groups;
154 gmx_ekindata_t *ekind, *ekind_save;
155 gmx_shellfc_t shellfc;
156 int count,nconverged=0;
159 gmx_bool bIonize=FALSE;
160 gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
162 gmx_bool bResetCountersHalfMaxH=FALSE;
163 gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
166 atom_id *grpindex=NULL;
168 t_coupl_rec *tcr=NULL;
169 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
170 matrix boxcopy={{0}},lastbox;
172 real fom,oldfom,veta_save,pcurr,scalevir,tracevir;
180 real saved_conserved_quantity = 0;
185 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
186 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
187 gmx_iterate_t iterate;
188 gmx_large_int_t multisim_nsteps=-1; /* number of steps to do before first multisim
189 simulation stops. If equal to zero, don't
190 communicate any more between multisims.*/
195 "\n* WARNING * WARNING * WARNING * WARNING * WARNING * WARNING *\n"
196 "We have just committed the new CPU detection code in this branch,\n"
197 "and will commit new SSE/AVX kernels in a few days. However, this\n"
198 "means that currently only the NxN kernels are accelerated!\n"
199 "In the mean time, you might want to avoid production runs in 4.6.\n\n");
203 /* Temporary addition for FAHCORE checkpointing */
207 /* Check for special mdrun options */
208 bRerunMD = (Flags & MD_RERUN);
209 bIonize = (Flags & MD_IONIZE);
210 bFFscan = (Flags & MD_FFSCAN);
211 bAppend = (Flags & MD_APPENDFILES);
212 if (Flags & MD_RESETCOUNTERSHALFWAY)
216 /* Signal to reset the counters half the simulation steps. */
217 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
219 /* Signal to reset the counters halfway the simulation time. */
220 bResetCountersHalfMaxH = (max_hours > 0);
223 /* md-vv uses averaged full step velocities for T-control
224 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
225 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
227 if (bVV) /* to store the initial velocities while computing virial */
229 snew(cbuf,top_global->natoms);
231 /* all the iteratative cases - only if there are constraints */
232 bIterations = ((IR_NPH_TROTTER(ir) || IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
233 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir)));
237 /* Since we don't know if the frames read are related in any way,
238 * rebuild the neighborlist at every step.
241 ir->nstcalcenergy = 1;
245 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
247 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
248 bGStatEveryStep = (nstglobalcomm == 1);
250 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
253 "To reduce the energy communication with nstlist = -1\n"
254 "the neighbor list validity should not be checked at every step,\n"
255 "this means that exact integration is not guaranteed.\n"
256 "The neighbor list validity is checked after:\n"
257 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
258 "In most cases this will result in exact integration.\n"
259 "This reduces the energy communication by a factor of 2 to 3.\n"
260 "If you want less energy communication, set nstlist > 3.\n\n");
263 if (bRerunMD || bFFscan)
267 groups = &top_global->groups;
270 init_md(fplog,cr,ir,oenv,&t,&t0,state_global->lambda,
271 &(state_global->fep_state),lam0,
272 nrnb,top_global,&upd,
273 nfile,fnm,&outf,&mdebin,
274 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
276 clear_mat(total_vir);
278 /* Energy terms and groups */
280 init_enerdata(top_global->groups.grps[egcENER].nr,ir->fepvals->n_lambda,
282 if (DOMAINDECOMP(cr))
288 snew(f,top_global->natoms);
291 /* lambda Monte carlo random number generator */
294 mcrng = gmx_rng_init(ir->expandedvals->lmc_seed);
296 /* copy the state into df_history */
297 copy_df_history(&df_history,&state_global->dfhist);
299 /* Kinetic energy data */
301 init_ekindata(fplog,top_global,&(ir->opts),ekind);
302 /* needed for iteration of constraints */
304 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
305 /* Copy the cos acceleration to the groups struct */
306 ekind->cosacc.cos_accel = ir->cos_accel;
308 gstat = global_stat_init(ir);
311 /* Check for polarizable models and flexible constraints */
312 shellfc = init_shell_flexcon(fplog,
313 top_global,n_flexible_constraints(constr),
314 (ir->bContinuation ||
315 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
316 NULL : state_global->x);
320 #ifdef GMX_THREAD_MPI
321 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
323 set_deform_reference_box(upd,
324 deform_init_init_step_tpx,
325 deform_init_box_tpx);
326 #ifdef GMX_THREAD_MPI
327 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
332 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
333 if ((io > 2000) && MASTER(cr))
335 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
339 if (DOMAINDECOMP(cr)) {
340 top = dd_init_local_top(top_global);
343 dd_init_local_state(cr->dd,state_global,state);
345 if (DDMASTER(cr->dd) && ir->nstfout) {
346 snew(f_global,state_global->natoms);
350 /* Initialize the particle decomposition and split the topology */
351 top = split_system(fplog,top_global,ir,cr);
353 pd_cg_range(cr,&fr->cg0,&fr->hcg);
354 pd_at_range(cr,&a0,&a1);
356 top = gmx_mtop_generate_local_top(top_global,ir);
359 a1 = top_global->natoms;
362 state = partdec_init_local_state(cr,state_global);
365 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
368 set_vsite_top(vsite,top,mdatoms,cr);
371 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
372 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
376 make_local_shells(cr,mdatoms,shellfc);
379 if (ir->pull && PAR(cr)) {
380 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
384 if (DOMAINDECOMP(cr))
386 /* Distribute the charge groups over the nodes from the master node */
387 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
388 state_global,top_global,ir,
389 state,&f,mdatoms,top,fr,
390 vsite,shellfc,constr,
394 update_mdatoms(mdatoms,state->lambda[efptMASS]);
396 if (opt2bSet("-cpi",nfile,fnm))
398 bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr);
402 bStateFromCP = FALSE;
409 /* Update mdebin with energy history if appending to output files */
410 if ( Flags & MD_APPENDFILES )
412 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
416 /* We might have read an energy history from checkpoint,
417 * free the allocated memory and reset the counts.
419 done_energyhistory(&state_global->enerhist);
420 init_energyhistory(&state_global->enerhist);
423 /* Set the initial energy history in state by updating once */
424 update_energyhistory(&state_global->enerhist,mdebin);
427 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG))
429 /* Set the random state if we read a checkpoint file */
430 set_stochd_state(upd,state);
433 if (state->flags & (1<<estMC_RNG))
435 set_mc_state(mcrng,state);
438 /* Initialize constraints */
440 if (!DOMAINDECOMP(cr))
441 set_constraints(constr,top,ir,mdatoms,cr);
444 /* Check whether we have to GCT stuff */
445 bTCR = ftp2bSet(efGCT,nfile,fnm);
448 fprintf(stderr,"Will do General Coupling Theory!\n");
450 gnx = top_global->mols.nr;
452 for(i=0; (i<gnx); i++) {
459 /* We need to be sure replica exchange can only occur
460 * when the energies are current */
461 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
462 "repl_ex_nst",&repl_ex_nst);
463 /* This check needs to happen before inter-simulation
464 * signals are initialized, too */
466 if (repl_ex_nst > 0 && MASTER(cr))
468 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
469 repl_ex_nst,repl_ex_nex,repl_ex_seed);
471 if (!ir->bContinuation && !bRerunMD)
473 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
475 /* Set the velocities of frozen particles to zero */
476 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
480 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
490 /* Constrain the initial coordinates and velocities */
491 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
492 graph,cr,nrnb,fr,top,shake_vir);
496 /* Construct the virtual sites for the initial configuration */
497 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
498 top->idef.iparams,top->idef.il,
499 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
505 /* set free energy calculation frequency as the minimum of nstdhdl, nstexpanded, and nstrepl_ex_nst*/
506 nstfep = ir->fepvals->nstdhdl;
507 if (ir->bExpanded && (nstfep > ir->expandedvals->nstexpanded))
509 nstfep = ir->expandedvals->nstexpanded;
511 if (repl_ex_nst > 0 && repl_ex_nst > nstfep)
513 nstfep = repl_ex_nst;
516 /* I'm assuming we need global communication the first time! MRS */
517 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
518 | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM:0)
519 | (bVV ? CGLO_PRESSURE:0)
520 | (bVV ? CGLO_CONSTRAINT:0)
521 | (bRerunMD ? CGLO_RERUNMD:0)
522 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
524 bSumEkinhOld = FALSE;
525 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
526 NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
527 constr,NULL,FALSE,state->box,
528 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
529 if (ir->eI == eiVVAK) {
530 /* a second call to get the half step temperature initialized as well */
531 /* we do the same call as above, but turn the pressure off -- internally to
532 compute_globals, this is recognized as a velocity verlet half-step
533 kinetic energy calculation. This minimized excess variables, but
534 perhaps loses some logic?*/
536 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
537 NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
538 constr,NULL,FALSE,state->box,
539 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
540 cglo_flags &~ (CGLO_STOPCM | CGLO_PRESSURE));
543 /* Calculate the initial half step temperature, and save the ekinh_old */
544 if (!(Flags & MD_STARTFROMCPT))
546 for(i=0; (i<ir->opts.ngtc); i++)
548 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
553 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
554 and there is no previous step */
557 /* if using an iterative algorithm, we need to create a working directory for the state. */
560 bufstate = init_bufstate(state);
564 snew(xcopy,state->natoms);
565 snew(vcopy,state->natoms);
566 copy_rvecn(state->x,xcopy,0,state->natoms);
567 copy_rvecn(state->v,vcopy,0,state->natoms);
568 copy_mat(state->box,boxcopy);
571 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
572 temperature control */
573 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
577 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
580 "RMS relative constraint deviation after constraining: %.2e\n",
581 constr_rmsd(constr,FALSE));
583 if (EI_STATE_VELOCITY(ir->eI))
585 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
589 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
590 " input trajectory '%s'\n\n",
591 *(top_global->name),opt2fn("-rerun",nfile,fnm));
594 fprintf(stderr,"Calculated time to finish depends on nsteps from "
595 "run input file,\nwhich may not correspond to the time "
596 "needed to process input trajectory.\n\n");
602 fprintf(stderr,"starting mdrun '%s'\n",
603 *(top_global->name));
606 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
610 sprintf(tbuf,"%s","infinite");
612 if (ir->init_step > 0)
614 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
615 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
616 gmx_step_str(ir->init_step,sbuf2),
617 ir->init_step*ir->delta_t);
621 fprintf(stderr,"%s steps, %s ps.\n",
622 gmx_step_str(ir->nsteps,sbuf),tbuf);
628 /* Set and write start time */
629 runtime_start(runtime);
630 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
631 wallcycle_start(wcycle,ewcRUN);
637 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
639 chkpt_ret=fcCheckPointParallel( cr->nodeid,
641 if ( chkpt_ret == 0 )
642 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
646 /***********************************************************
650 ************************************************************/
652 /* if rerunMD then read coordinates and velocities from input trajectory */
655 if (getenv("GMX_FORCE_UPDATE"))
663 bNotLastFrame = read_first_frame(oenv,&status,
664 opt2fn("-rerun",nfile,fnm),
665 &rerun_fr,TRX_NEED_X | TRX_READ_V);
666 if (rerun_fr.natoms != top_global->natoms)
669 "Number of atoms in trajectory (%d) does not match the "
670 "run input file (%d)\n",
671 rerun_fr.natoms,top_global->natoms);
673 if (ir->ePBC != epbcNONE)
677 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);
679 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
681 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
688 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
691 if (ir->ePBC != epbcNONE)
693 /* Set the shift vectors.
694 * Necessary here when have a static box different from the tpr box.
696 calc_shifts(rerun_fr.box,fr->shift_vec);
700 /* loop over MD steps or if rerunMD to end of input trajectory */
702 /* Skip the first Nose-Hoover integration when we get the state from tpx */
703 bStateFromTPX = !bStateFromCP;
704 bInitStep = bFirstStep && (bStateFromTPX || bVV);
705 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
707 bSumEkinhOld = FALSE;
710 init_global_signals(&gs,cr,ir,repl_ex_nst);
712 step = ir->init_step;
715 if (ir->nstlist == -1)
717 init_nlistheuristics(&nlh,bGStatEveryStep,step);
720 if (MULTISIM(cr) && (repl_ex_nst <=0 ))
722 /* check how many steps are left in other sims */
723 multisim_nsteps=get_multisim_nsteps(cr, ir->nsteps);
727 /* and stop now if we should */
728 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
729 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
730 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
732 wallcycle_start(wcycle,ewcSTEP);
735 if (rerun_fr.bStep) {
736 step = rerun_fr.step;
737 step_rel = step - ir->init_step;
739 if (rerun_fr.bTime) {
749 bLastStep = (step_rel == ir->nsteps);
750 t = t0 + step*ir->delta_t;
753 if (ir->efep != efepNO || ir->bSimTemp)
755 /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
756 requiring different logic. */
758 set_current_lambdas(step,ir->fepvals,bRerunMD,&rerun_fr,state_global,state,lam0);
759 bDoDHDL = do_per_step(step,ir->fepvals->nstdhdl);
760 bDoFEP = (do_per_step(step,nstfep) && (ir->efep != efepNO));
761 bDoExpanded = (do_per_step(step,ir->expandedvals->nstexpanded) && (ir->bExpanded) && (step > 0));
766 update_annealing_target_temp(&(ir->opts),t);
771 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
773 for(i=0; i<state_global->natoms; i++)
775 copy_rvec(rerun_fr.x[i],state_global->x[i]);
779 for(i=0; i<state_global->natoms; i++)
781 copy_rvec(rerun_fr.v[i],state_global->v[i]);
786 for(i=0; i<state_global->natoms; i++)
788 clear_rvec(state_global->v[i]);
792 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
793 " Ekin, temperature and pressure are incorrect,\n"
794 " the virial will be incorrect when constraints are present.\n"
796 bRerunWarnNoV = FALSE;
800 copy_mat(rerun_fr.box,state_global->box);
801 copy_mat(state_global->box,state->box);
803 if (vsite && (Flags & MD_RERUN_VSITE))
805 if (DOMAINDECOMP(cr))
807 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
811 /* Following is necessary because the graph may get out of sync
812 * with the coordinates if we only have every N'th coordinate set
814 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
815 shift_self(graph,state->box,state->x);
817 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
818 top->idef.iparams,top->idef.il,
819 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
822 unshift_self(graph,state->box,state->x);
827 /* Stop Center of Mass motion */
828 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
830 /* Copy back starting coordinates in case we're doing a forcefield scan */
833 for(ii=0; (ii<state->natoms); ii++)
835 copy_rvec(xcopy[ii],state->x[ii]);
836 copy_rvec(vcopy[ii],state->v[ii]);
838 copy_mat(boxcopy,state->box);
843 /* for rerun MD always do Neighbour Searching */
844 bNS = (bFirstStep || ir->nstlist != 0);
849 /* Determine whether or not to do Neighbour Searching and LR */
850 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
852 bNS = (bFirstStep || bExchanged || bNStList || bDoFEP ||
853 (ir->nstlist == -1 && nlh.nabnsb > 0));
855 if (bNS && ir->nstlist == -1)
857 set_nlistheuristics(&nlh,bFirstStep || bExchanged || bDoFEP, step);
861 /* check whether we should stop because another simulation has
865 if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
866 (multisim_nsteps != ir->nsteps) )
873 "Stopping simulation %d because another one has finished\n",
877 gs.sig[eglsCHKPT] = 1;
882 /* < 0 means stop at next step, > 0 means stop at next NS step */
883 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
884 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
889 /* Determine whether or not to update the Born radii if doing GB */
890 bBornRadii=bFirstStep;
891 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
896 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
897 do_verbose = bVerbose &&
898 (step % stepout == 0 || bFirstStep || bLastStep);
900 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
908 bMasterState = FALSE;
909 /* Correct the new box if it is too skewed */
910 if (DYNAMIC_BOX(*ir))
912 if (correct_box(fplog,step,state->box,graph))
917 if (DOMAINDECOMP(cr) && bMasterState)
919 dd_collect_state(cr->dd,state,state_global);
923 if (DOMAINDECOMP(cr))
925 /* Repartition the domain decomposition */
926 wallcycle_start(wcycle,ewcDOMDEC);
927 dd_partition_system(fplog,step,cr,
928 bMasterState,nstglobalcomm,
929 state_global,top_global,ir,
930 state,&f,mdatoms,top,fr,
931 vsite,shellfc,constr,
932 nrnb,wcycle,do_verbose);
933 wallcycle_stop(wcycle,ewcDOMDEC);
934 /* If using an iterative integrator, reallocate space to match the decomposition */
938 if (MASTER(cr) && do_log && !bFFscan)
940 print_ebin_header(fplog,step,t,state->lambda[efptFEP]); /* can we improve the information printed here? */
943 if (ir->efep != efepNO)
945 update_mdatoms(mdatoms,state->lambda[efptMASS]);
948 if (bRerunMD && rerun_fr.bV)
951 /* We need the kinetic energy at minus the half step for determining
952 * the full step kinetic energy and possibly for T-coupling.*/
953 /* This may not be quite working correctly yet . . . . */
954 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
955 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
956 constr,NULL,FALSE,state->box,
957 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
958 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
960 clear_mat(force_vir);
962 /* Ionize the atoms if necessary */
965 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
966 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
969 /* Update force field in ffscan program */
972 if (update_forcefield(fplog,
974 mdatoms->nr,state->x,state->box))
982 /* We write a checkpoint at this MD step when:
983 * either at an NS step when we signalled through gs,
984 * or at the last step (but not when we do not want confout),
985 * but never at the first step or with rerun.
987 bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
988 (bLastStep && (Flags & MD_CONFOUT))) &&
989 step > ir->init_step && !bRerunMD);
992 gs.set[eglsCHKPT] = 0;
995 /* Determine the energy and pressure:
996 * at nstcalcenergy steps and at energy output steps (set below).
999 if (EI_VV(ir->eI) && (!bInitStep)) { /* for vv, the first half actually corresponds to the last step */
1000 bNstEner = do_per_step(step-1,ir->nstcalcenergy);
1002 bNstEner = do_per_step(step,ir->nstcalcenergy);
1006 (ir->epc > epcNO && do_per_step(step,ir->nstpcouple)));
1008 /* Do we need global communication ? */
1009 bGStat = (bCalcEnerPres || bStopCM ||
1010 do_per_step(step,nstglobalcomm) ||
1011 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1013 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1015 if (do_ene || do_log)
1017 bCalcEnerPres = TRUE;
1021 /* these CGLO_ options remain the same throughout the iteration */
1022 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1023 (bGStat ? CGLO_GSTAT : 0)
1026 force_flags = (GMX_FORCE_STATECHANGED |
1027 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1028 GMX_FORCE_ALLFORCES |
1029 (bNStList ? GMX_FORCE_DOLR : 0) |
1031 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
1032 (bDoFEP ? GMX_FORCE_DHDL : 0)
1037 /* Now is the time to relax the shells */
1038 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
1040 bStopCM,top,top_global,
1042 state,f,force_vir,mdatoms,
1043 nrnb,wcycle,graph,groups,
1044 shellfc,fr,bBornRadii,t,mu_tot,
1045 state->natoms,&bConverged,vsite,
1056 /* The coordinates (x) are shifted (to get whole molecules)
1058 * This is parallellized as well, and does communication too.
1059 * Check comments in sim_util.c
1061 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1062 state->box,state->x,&state->hist,
1063 f,force_vir,mdatoms,enerd,fcd,
1064 state->lambda,graph,
1065 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1066 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1071 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1072 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1075 if (bTCR && bFirstStep)
1077 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1078 fprintf(fplog,"Done init_coupling\n");
1082 if (bVV && !bStartingFromCpt && !bRerunMD)
1083 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1085 if (ir->eI==eiVV && bInitStep)
1087 /* if using velocity verlet with full time step Ekin,
1088 * take the first half step only to compute the
1089 * virial for the first step. From there,
1090 * revert back to the initial coordinates
1091 * so that the input is actually the initial step.
1093 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1095 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1096 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1099 update_coords(fplog,step,ir,mdatoms,state,
1100 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1101 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1102 cr,nrnb,constr,&top->idef);
1106 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1108 /* for iterations, we save these vectors, as we will be self-consistently iterating
1111 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1113 /* save the state */
1114 if (bIterations && iterate.bIterate) {
1115 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1118 bFirstIterate = TRUE;
1119 while (bFirstIterate || (bIterations && iterate.bIterate))
1121 if (bIterations && iterate.bIterate)
1123 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1124 if (bFirstIterate && bTrotter)
1126 /* The first time through, we need a decent first estimate
1127 of veta(t+dt) to compute the constraints. Do
1128 this by computing the box volume part of the
1129 trotter integration at this time. Nothing else
1130 should be changed by this routine here. If
1131 !(first time), we start with the previous value
1134 veta_save = state->veta;
1135 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1136 vetanew = state->veta;
1137 state->veta = veta_save;
1142 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
1145 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1146 &top->idef,shake_vir,NULL,
1147 cr,nrnb,wcycle,upd,constr,
1148 bInitStep,TRUE,bCalcEnerPres,vetanew);
1150 if (!bOK && !bFFscan)
1152 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1157 { /* Need to unshift here if a do_force has been
1158 called in the previous step */
1159 unshift_self(graph,state->box,state->x);
1163 /* if VV, compute the pressure and constraints */
1164 /* For VV2, we strictly only need this if using pressure
1165 * control, but we really would like to have accurate pressures
1167 * Think about ways around this in the future?
1168 * For now, keep this choice in comments.
1170 /* bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1171 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1173 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
1174 if (bNstEner && ir->eI==eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/
1176 bSumEkinhOld = TRUE;
1178 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1179 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1180 constr,NULL,FALSE,state->box,
1181 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1184 | (bStopCM ? CGLO_STOPCM : 0)
1185 | (bTemp ? CGLO_TEMPERATURE:0)
1186 | (bPres ? CGLO_PRESSURE : 0)
1187 | (bPres ? CGLO_CONSTRAINT : 0)
1188 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
1189 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1192 /* explanation of above:
1193 a) We compute Ekin at the full time step
1194 if 1) we are using the AveVel Ekin, and it's not the
1195 initial step, or 2) if we are using AveEkin, but need the full
1196 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1197 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1198 EkinAveVel because it's needed for the pressure */
1200 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1205 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
1209 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1214 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1215 state->veta,&vetanew))
1219 bFirstIterate = FALSE;
1222 if (bTrotter && !bInitStep) {
1223 enerd->term[F_DVDL_BONDED] += dvdl; /* only add after iterations */
1224 copy_mat(shake_vir,state->svir_prev);
1225 copy_mat(force_vir,state->fvir_prev);
1226 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
1227 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1228 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
1229 enerd->term[F_EKIN] = trace(ekind->ekin);
1232 /* if it's the initial step, we performed this first step just to get the constraint virial */
1233 if (bInitStep && ir->eI==eiVV) {
1234 copy_rvecn(cbuf,state->v,0,state->natoms);
1237 if (fr->bSepDVDL && fplog && do_log)
1239 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1241 enerd->term[F_DVDL_BONDED] += dvdl;
1244 /* MRS -- now done iterating -- compute the conserved quantity */
1246 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
1249 last_ekin = enerd->term[F_EKIN];
1251 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1253 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1255 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1256 sum_dhdl(enerd,state->lambda,ir->fepvals);
1259 /* ######## END FIRST UPDATE STEP ############## */
1260 /* ######## If doing VV, we now have v(dt) ###### */
1262 /* perform extended ensemble sampling in lambda - we don't
1263 actually move to the new state before outputting
1264 statistics, but if performing simulated tempering, we
1265 do update the velocities and the tau_t. */
1267 lamnew = ExpandedEnsembleDynamics(fplog,ir,enerd,state,&MassQ,&df_history,step,mcrng,state->v,mdatoms);
1269 /* ################## START TRAJECTORY OUTPUT ################# */
1271 /* Now we have the energies and forces corresponding to the
1272 * coordinates at time t. We must output all of this before
1274 * for RerunMD t is read from input trajectory
1277 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
1278 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
1279 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
1280 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
1281 if (bCPT) { mdof_flags |= MDOF_CPT; };
1283 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1286 /* Enforce writing positions and velocities at end of run */
1287 mdof_flags |= (MDOF_X | MDOF_V);
1292 fcReportProgress( ir->nsteps, step );
1294 /* sync bCPT and fc record-keeping */
1295 if (bCPT && MASTER(cr))
1296 fcRequestCheckPoint();
1299 if (mdof_flags != 0)
1301 wallcycle_start(wcycle,ewcTRAJ);
1304 if (state->flags & (1<<estLD_RNG))
1306 get_stochd_state(upd,state);
1308 if (state->flags & (1<<estMC_RNG))
1310 get_mc_state(mcrng,state);
1316 state_global->ekinstate.bUpToDate = FALSE;
1320 update_ekinstate(&state_global->ekinstate,ekind);
1321 state_global->ekinstate.bUpToDate = TRUE;
1323 update_energyhistory(&state_global->enerhist,mdebin);
1324 if (ir->efep!=efepNO || ir->bSimTemp)
1326 state_global->fep_state = state->fep_state; /* MRS: seems kludgy. The code should be
1327 structured so this isn't necessary.
1328 Note this reassignment is only necessary
1329 for single threads.*/
1330 copy_df_history(&state_global->dfhist,&df_history);
1334 write_traj(fplog,cr,outf,mdof_flags,top_global,
1335 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
1342 if (bLastStep && step_rel == ir->nsteps &&
1343 (Flags & MD_CONFOUT) && MASTER(cr) &&
1344 !bRerunMD && !bFFscan)
1346 /* x and v have been collected in write_traj,
1347 * because a checkpoint file will always be written
1350 fprintf(stderr,"\nWriting final coordinates.\n");
1351 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
1354 /* Make molecules whole only for confout writing */
1355 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
1357 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
1358 *top_global->name,top_global,
1359 state_global->x,state_global->v,
1360 ir->ePBC,state->box);
1363 wallcycle_stop(wcycle,ewcTRAJ);
1366 /* kludge -- virial is lost with restart for NPT control. Must restart */
1367 if (bStartingFromCpt && bVV)
1369 copy_mat(state->svir_prev,shake_vir);
1370 copy_mat(state->fvir_prev,force_vir);
1372 /* ################## END TRAJECTORY OUTPUT ################ */
1374 /* Determine the pressure:
1375 * always when we want exact averages in the energy file,
1376 * at ns steps when we have pressure coupling,
1377 * otherwise only at energy output steps (set below).
1381 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
1382 bCalcEnerPres = bNstEner;
1384 /* Do we need global communication ? */
1385 bGStat = (bGStatEveryStep || bStopCM || bNS ||
1386 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1388 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1390 if (do_ene || do_log)
1392 bCalcEnerPres = TRUE;
1396 /* Determine the wallclock run time up till now */
1397 run_time = gmx_gettime() - (double)runtime->real;
1398 /* Check whether everything is still allright */
1399 if (((int)gmx_get_stop_condition() > handled_stop_condition)
1400 #ifdef GMX_THREAD_MPI
1405 /* this is just make gs.sig compatible with the hack
1406 of sending signals around by MPI_Reduce with together with
1408 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
1409 gs.sig[eglsSTOPCOND]=1;
1410 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
1411 gs.sig[eglsSTOPCOND]=-1;
1412 /* < 0 means stop at next step, > 0 means stop at next NS step */
1416 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1417 gmx_get_signal_name(),
1418 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1422 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1423 gmx_get_signal_name(),
1424 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1426 handled_stop_condition=(int)gmx_get_stop_condition();
1428 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1429 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
1430 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
1432 /* Signal to terminate the run */
1433 gs.sig[eglsSTOPCOND] = 1;
1436 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1438 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1441 if (bResetCountersHalfMaxH && MASTER(cr) &&
1442 run_time > max_hours*60.0*60.0*0.495)
1444 gs.sig[eglsRESETCOUNTERS] = 1;
1447 if (ir->nstlist == -1 && !bRerunMD)
1449 /* When bGStatEveryStep=FALSE, global_stat is only called
1450 * when we check the atom displacements, not at NS steps.
1451 * This means that also the bonded interaction count check is not
1452 * performed immediately after NS. Therefore a few MD steps could
1453 * be performed with missing interactions.
1454 * But wrong energies are never written to file,
1455 * since energies are only written after global_stat
1458 if (step >= nlh.step_nscheck)
1460 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
1461 nlh.scale_tot,state->x);
1465 /* This is not necessarily true,
1466 * but step_nscheck is determined quite conservatively.
1472 /* In parallel we only have to check for checkpointing in steps
1473 * where we do global communication,
1474 * otherwise the other nodes don't know.
1476 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1479 run_time >= nchkpt*cpt_period*60.0)) &&
1480 gs.set[eglsCHKPT] == 0)
1482 gs.sig[eglsCHKPT] = 1;
1486 /* at the start of step, randomize the velocities */
1487 if (ETC_ANDERSEN(ir->etc) && EI_VV(ir->eI))
1489 gmx_bool bDoAndersenConstr;
1490 bDoAndersenConstr = (constr && update_randomize_velocities(ir,step,mdatoms,state,upd,&top->idef,constr));
1491 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1492 if (bDoAndersenConstr)
1494 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1495 &top->idef,tmp_vir,NULL,
1496 cr,nrnb,wcycle,upd,constr,
1497 bInitStep,TRUE,FALSE,vetanew);
1503 gmx_iterate_init(&iterate,bIterations);
1506 /* for iterations, we save these vectors, as we will be redoing the calculations */
1507 if (bIterations && iterate.bIterate)
1509 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1511 bFirstIterate = TRUE;
1512 while (bFirstIterate || (bIterations && iterate.bIterate))
1514 /* We now restore these vectors to redo the calculation with improved extended variables */
1517 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1520 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1521 so scroll down for that logic */
1523 /* ######### START SECOND UPDATE STEP ################# */
1524 /* Box is changed in update() when we do pressure coupling,
1525 * but we should still use the old box for energy corrections and when
1526 * writing it to the energy file, so it matches the trajectory files for
1527 * the same timestep above. Make a copy in a separate array.
1529 copy_mat(state->box,lastbox);
1532 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
1534 wallcycle_start(wcycle,ewcUPDATE);
1536 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1539 if (bIterations && iterate.bIterate)
1547 /* we use a new value of scalevir to converge the iterations faster */
1548 scalevir = tracevir/trace(shake_vir);
1550 msmul(shake_vir,scalevir,shake_vir);
1551 m_add(force_vir,shake_vir,total_vir);
1552 clear_mat(shake_vir);
1554 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
1555 /* We can only do Berendsen coupling after we have summed
1556 * the kinetic energy or virial. Since the happens
1557 * in global_state after update, we should only do it at
1558 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1563 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1564 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
1570 /* velocity half-step update */
1571 update_coords(fplog,step,ir,mdatoms,state,f,
1572 fr->bTwinRange && bNStList,fr->f_twin,fcd,
1573 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
1574 cr,nrnb,constr,&top->idef);
1577 /* Above, initialize just copies ekinh into ekin,
1578 * it doesn't copy position (for VV),
1579 * and entire integrator for MD.
1584 copy_rvecn(state->x,cbuf,0,state->natoms);
1587 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1588 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1589 wallcycle_stop(wcycle,ewcUPDATE);
1591 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1592 &top->idef,shake_vir,force_vir,
1593 cr,nrnb,wcycle,upd,constr,
1594 bInitStep,FALSE,bCalcEnerPres,state->veta);
1598 /* erase F_EKIN and F_TEMP here? */
1599 /* just compute the kinetic energy at the half step to perform a trotter step */
1600 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1601 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1602 constr,NULL,FALSE,lastbox,
1603 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1604 cglo_flags | CGLO_TEMPERATURE
1606 wallcycle_start(wcycle,ewcUPDATE);
1607 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
1608 /* now we know the scaling, we can compute the positions again again */
1609 copy_rvecn(cbuf,state->x,0,state->natoms);
1611 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1612 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1613 wallcycle_stop(wcycle,ewcUPDATE);
1615 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1616 /* are the small terms in the shake_vir here due
1617 * to numerical errors, or are they important
1618 * physically? I'm thinking they are just errors, but not completely sure.
1619 * For now, will call without actually constraining, constr=NULL*/
1620 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1621 &top->idef,tmp_vir,force_vir,
1622 cr,nrnb,wcycle,upd,NULL,
1623 bInitStep,FALSE,bCalcEnerPres,
1626 if (!bOK && !bFFscan)
1628 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1631 if (fr->bSepDVDL && fplog && do_log)
1633 fprintf(fplog,sepdvdlformat,"Constraint dV/dl",0.0,dvdl);
1635 enerd->term[F_DVDL_BONDED] += dvdl;
1639 /* Need to unshift here */
1640 unshift_self(graph,state->box,state->x);
1645 wallcycle_start(wcycle,ewcVSITECONSTR);
1648 shift_self(graph,state->box,state->x);
1650 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1651 top->idef.iparams,top->idef.il,
1652 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1656 unshift_self(graph,state->box,state->x);
1658 wallcycle_stop(wcycle,ewcVSITECONSTR);
1661 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1662 if (ir->nstlist == -1 && bFirstIterate)
1664 gs.sig[eglsNABNSB] = nlh.nabnsb;
1666 bEnergyHere = (!EI_VV(ir->eI) || (EI_VV(ir->eI) && bRerunMD)); /* this is not quite working for vv and rerun! fails for running rerun on multiple threads. This is caught in runner.c. */
1667 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1668 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1670 bFirstIterate ? &gs : NULL,
1671 (step_rel % gs.nstms == 0) &&
1672 (multisim_nsteps<0 || (step_rel<multisim_nsteps)),
1674 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1676 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
1677 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1678 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1679 | (bEnergyHere || bRerunMD ? CGLO_PRESSURE : 0)
1680 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
1681 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1684 if (ir->nstlist == -1 && bFirstIterate)
1686 nlh.nabnsb = gs.set[eglsNABNSB];
1687 gs.set[eglsNABNSB] = 0;
1689 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1690 /* ############# END CALC EKIN AND PRESSURE ################# */
1692 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1693 the virial that should probably be addressed eventually. state->veta has better properies,
1694 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1695 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1698 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1699 trace(shake_vir),&tracevir))
1703 bFirstIterate = FALSE;
1706 /* only add constraint dvdl after constraints */
1707 enerd->term[F_DVDL_BONDED] += dvdl;
1710 /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */
1711 sum_dhdl(enerd,state->lambda,ir->fepvals);
1713 update_box(fplog,step,ir,mdatoms,state,graph,f,
1714 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
1716 /* ################# END UPDATE STEP 2 ################# */
1717 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1719 /* The coordinates (x) were unshifted in update */
1720 if (bFFscan && (shellfc==NULL || bConverged))
1722 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
1724 &(top_global->mols),mdatoms->massT,pres))
1728 fprintf(stderr,"\n");
1734 /* We will not sum ekinh_old,
1735 * so signal that we still have to do it.
1737 bSumEkinhOld = TRUE;
1742 /* Only do GCT when the relaxation of shells (minimization) has converged,
1743 * otherwise we might be coupling to bogus energies.
1744 * In parallel we must always do this, because the other sims might
1748 /* Since this is called with the new coordinates state->x, I assume
1749 * we want the new box state->box too. / EL 20040121
1751 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
1753 mdatoms,&(top->idef),mu_aver,
1754 top_global->mols.nr,cr,
1755 state->box,total_vir,pres,
1756 mu_tot,state->x,f,bConverged);
1760 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1762 /* use the directly determined last velocity, not actually the averaged half steps */
1763 if (bTrotter && ir->eI==eiVV)
1765 enerd->term[F_EKIN] = last_ekin;
1767 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1771 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1775 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
1777 /* Check for excessively large energies */
1781 real etot_max = 1e200;
1783 real etot_max = 1e30;
1785 if (fabs(enerd->term[F_ETOT]) > etot_max)
1787 fprintf(stderr,"Energy too large (%g), giving up\n",
1788 enerd->term[F_ETOT]);
1791 /* ######### END PREPARING EDR OUTPUT ########### */
1793 /* Time for performance */
1794 if (((step % stepout) == 0) || bLastStep)
1796 runtime_upd_proc(runtime);
1802 gmx_bool do_dr,do_or;
1804 if (fplog && do_log && bDoExpanded)
1806 /* only needed if doing expanded ensemble */
1807 PrintFreeEnergyInfoToFile(fplog,ir->fepvals,ir->expandedvals,ir->bSimTemp?ir->simtempvals:NULL,
1808 &df_history,state->fep_state,ir->nstlog,step);
1810 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
1814 upd_mdebin(mdebin,bDoDHDL,TRUE,
1815 t,mdatoms->tmass,enerd,state,
1816 ir->fepvals,ir->expandedvals,lastbox,
1817 shake_vir,force_vir,total_vir,pres,
1818 ekind,mu_tot,constr);
1822 upd_mdebin_step(mdebin);
1825 do_dr = do_per_step(step,ir->nstdisreout);
1826 do_or = do_per_step(step,ir->nstorireout);
1828 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
1830 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
1832 if (ir->ePull != epullNO)
1834 pull_print_output(ir->pull,step,t);
1837 if (do_per_step(step,ir->nstlog))
1839 if(fflush(fplog) != 0)
1841 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
1847 /* Have to do this part after outputting the logfile and the edr file */
1848 state->fep_state = lamnew;
1849 for (i=0;i<efptNR;i++)
1851 state->lambda[i] = ir->fepvals->all_lambda[i][lamnew];
1854 /* Remaining runtime */
1855 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
1859 fprintf(stderr,"\n");
1861 print_time(stderr,runtime,step,ir,cr);
1864 /* Replica exchange */
1866 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
1867 do_per_step(step,repl_ex_nst))
1869 bExchanged = replica_exchange(fplog,cr,repl_ex,
1873 if (bExchanged && DOMAINDECOMP(cr))
1875 dd_partition_system(fplog,step,cr,TRUE,1,
1876 state_global,top_global,ir,
1877 state,&f,mdatoms,top,fr,
1878 vsite,shellfc,constr,
1885 bStartingFromCpt = FALSE;
1887 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1888 /* With all integrators, except VV, we need to retain the pressure
1889 * at the current step for coupling at the next step.
1891 if ((state->flags & (1<<estPRES_PREV)) &&
1893 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1895 /* Store the pressure in t_state for pressure coupling
1896 * at the next MD step.
1898 copy_mat(pres,state->pres_prev);
1901 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1903 if ( (membed!=NULL) && (!bLastStep) )
1905 rescale_membed(step_rel,membed,state_global->x);
1912 /* read next frame from input trajectory */
1913 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
1918 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
1922 if (!bRerunMD || !rerun_fr.bStep)
1924 /* increase the MD step number */
1929 cycles = wallcycle_stop(wcycle,ewcSTEP);
1930 if (DOMAINDECOMP(cr) && wcycle)
1932 dd_cycles_add(cr->dd,cycles,ddCyclStep);
1935 if (step_rel == wcycle_get_reset_counters(wcycle) ||
1936 gs.set[eglsRESETCOUNTERS] != 0)
1938 /* Reset all the counters related to performance over the run */
1939 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
1940 wcycle_set_reset_counters(wcycle,-1);
1941 /* Correct max_hours for the elapsed time */
1942 max_hours -= run_time/(60.0*60.0);
1943 bResetCountersHalfMaxH = FALSE;
1944 gs.set[eglsRESETCOUNTERS] = 0;
1948 /* End of main MD loop */
1952 runtime_end(runtime);
1954 if (bRerunMD && MASTER(cr))
1959 if (!(cr->duty & DUTY_PME))
1961 /* Tell the PME only node to finish */
1967 if (ir->nstcalcenergy > 0 && !bRerunMD)
1969 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
1970 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
1978 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
1980 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)));
1981 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
1984 if (shellfc && fplog)
1986 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
1987 (nconverged*100.0)/step_rel);
1988 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
1992 if (repl_ex_nst > 0 && MASTER(cr))
1994 print_replica_exchange_statistics(fplog,repl_ex);
1997 runtime->nsteps_done = step_rel;