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.*/
192 /* Temporary addition for FAHCORE checkpointing */
196 /* Check for special mdrun options */
197 bRerunMD = (Flags & MD_RERUN);
198 bIonize = (Flags & MD_IONIZE);
199 bFFscan = (Flags & MD_FFSCAN);
200 bAppend = (Flags & MD_APPENDFILES);
201 if (Flags & MD_RESETCOUNTERSHALFWAY)
205 /* Signal to reset the counters half the simulation steps. */
206 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
208 /* Signal to reset the counters halfway the simulation time. */
209 bResetCountersHalfMaxH = (max_hours > 0);
212 /* md-vv uses averaged full step velocities for T-control
213 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
214 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
216 if (bVV) /* to store the initial velocities while computing virial */
218 snew(cbuf,top_global->natoms);
220 /* all the iteratative cases - only if there are constraints */
221 bIterations = ((IR_NPH_TROTTER(ir) || IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
222 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir)));
226 /* Since we don't know if the frames read are related in any way,
227 * rebuild the neighborlist at every step.
230 ir->nstcalcenergy = 1;
234 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
236 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
237 bGStatEveryStep = (nstglobalcomm == 1);
239 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
242 "To reduce the energy communication with nstlist = -1\n"
243 "the neighbor list validity should not be checked at every step,\n"
244 "this means that exact integration is not guaranteed.\n"
245 "The neighbor list validity is checked after:\n"
246 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
247 "In most cases this will result in exact integration.\n"
248 "This reduces the energy communication by a factor of 2 to 3.\n"
249 "If you want less energy communication, set nstlist > 3.\n\n");
252 if (bRerunMD || bFFscan)
256 groups = &top_global->groups;
259 init_md(fplog,cr,ir,oenv,&t,&t0,state_global->lambda,
260 &(state_global->fep_state),lam0,
261 nrnb,top_global,&upd,
262 nfile,fnm,&outf,&mdebin,
263 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
265 clear_mat(total_vir);
267 /* Energy terms and groups */
269 init_enerdata(top_global->groups.grps[egcENER].nr,ir->fepvals->n_lambda,
271 if (DOMAINDECOMP(cr))
277 snew(f,top_global->natoms);
280 /* lambda Monte carlo random number generator */
283 mcrng = gmx_rng_init(ir->expandedvals->lmc_seed);
285 /* copy the state into df_history */
286 copy_df_history(&df_history,&state_global->dfhist);
288 /* Kinetic energy data */
290 init_ekindata(fplog,top_global,&(ir->opts),ekind);
291 /* needed for iteration of constraints */
293 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
294 /* Copy the cos acceleration to the groups struct */
295 ekind->cosacc.cos_accel = ir->cos_accel;
297 gstat = global_stat_init(ir);
300 /* Check for polarizable models and flexible constraints */
301 shellfc = init_shell_flexcon(fplog,
302 top_global,n_flexible_constraints(constr),
303 (ir->bContinuation ||
304 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
305 NULL : state_global->x);
309 #ifdef GMX_THREAD_MPI
310 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
312 set_deform_reference_box(upd,
313 deform_init_init_step_tpx,
314 deform_init_box_tpx);
315 #ifdef GMX_THREAD_MPI
316 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
321 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
322 if ((io > 2000) && MASTER(cr))
324 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
328 if (DOMAINDECOMP(cr)) {
329 top = dd_init_local_top(top_global);
332 dd_init_local_state(cr->dd,state_global,state);
334 if (DDMASTER(cr->dd) && ir->nstfout) {
335 snew(f_global,state_global->natoms);
339 /* Initialize the particle decomposition and split the topology */
340 top = split_system(fplog,top_global,ir,cr);
342 pd_cg_range(cr,&fr->cg0,&fr->hcg);
343 pd_at_range(cr,&a0,&a1);
345 top = gmx_mtop_generate_local_top(top_global,ir);
348 a1 = top_global->natoms;
351 state = partdec_init_local_state(cr,state_global);
354 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
357 set_vsite_top(vsite,top,mdatoms,cr);
360 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
361 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
365 make_local_shells(cr,mdatoms,shellfc);
368 if (ir->pull && PAR(cr)) {
369 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
373 if (DOMAINDECOMP(cr))
375 /* Distribute the charge groups over the nodes from the master node */
376 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
377 state_global,top_global,ir,
378 state,&f,mdatoms,top,fr,
379 vsite,shellfc,constr,
383 update_mdatoms(mdatoms,state->lambda[efptMASS]);
385 if (opt2bSet("-cpi",nfile,fnm))
387 bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr);
391 bStateFromCP = FALSE;
398 /* Update mdebin with energy history if appending to output files */
399 if ( Flags & MD_APPENDFILES )
401 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
405 /* We might have read an energy history from checkpoint,
406 * free the allocated memory and reset the counts.
408 done_energyhistory(&state_global->enerhist);
409 init_energyhistory(&state_global->enerhist);
412 /* Set the initial energy history in state by updating once */
413 update_energyhistory(&state_global->enerhist,mdebin);
416 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG))
418 /* Set the random state if we read a checkpoint file */
419 set_stochd_state(upd,state);
422 if (state->flags & (1<<estMC_RNG))
424 set_mc_state(mcrng,state);
427 /* Initialize constraints */
429 if (!DOMAINDECOMP(cr))
430 set_constraints(constr,top,ir,mdatoms,cr);
433 /* Check whether we have to GCT stuff */
434 bTCR = ftp2bSet(efGCT,nfile,fnm);
437 fprintf(stderr,"Will do General Coupling Theory!\n");
439 gnx = top_global->mols.nr;
441 for(i=0; (i<gnx); i++) {
448 /* We need to be sure replica exchange can only occur
449 * when the energies are current */
450 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
451 "repl_ex_nst",&repl_ex_nst);
452 /* This check needs to happen before inter-simulation
453 * signals are initialized, too */
455 if (repl_ex_nst > 0 && MASTER(cr))
457 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
458 repl_ex_nst,repl_ex_nex,repl_ex_seed);
460 if (!ir->bContinuation && !bRerunMD)
462 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
464 /* Set the velocities of frozen particles to zero */
465 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
469 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
479 /* Constrain the initial coordinates and velocities */
480 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
481 graph,cr,nrnb,fr,top,shake_vir);
485 /* Construct the virtual sites for the initial configuration */
486 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
487 top->idef.iparams,top->idef.il,
488 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
494 /* set free energy calculation frequency as the minimum of nstdhdl, nstexpanded, and nstrepl_ex_nst*/
495 nstfep = ir->fepvals->nstdhdl;
496 if (ir->bExpanded && (nstfep > ir->expandedvals->nstexpanded))
498 nstfep = ir->expandedvals->nstexpanded;
500 if (repl_ex_nst > 0 && repl_ex_nst > nstfep)
502 nstfep = repl_ex_nst;
505 /* I'm assuming we need global communication the first time! MRS */
506 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
507 | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM:0)
508 | (bVV ? CGLO_PRESSURE:0)
509 | (bVV ? CGLO_CONSTRAINT:0)
510 | (bRerunMD ? CGLO_RERUNMD:0)
511 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
513 bSumEkinhOld = FALSE;
514 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
515 NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
516 constr,NULL,FALSE,state->box,
517 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
518 if (ir->eI == eiVVAK) {
519 /* a second call to get the half step temperature initialized as well */
520 /* we do the same call as above, but turn the pressure off -- internally to
521 compute_globals, this is recognized as a velocity verlet half-step
522 kinetic energy calculation. This minimized excess variables, but
523 perhaps loses some logic?*/
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,
529 cglo_flags &~ (CGLO_STOPCM | CGLO_PRESSURE));
532 /* Calculate the initial half step temperature, and save the ekinh_old */
533 if (!(Flags & MD_STARTFROMCPT))
535 for(i=0; (i<ir->opts.ngtc); i++)
537 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
542 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
543 and there is no previous step */
546 /* if using an iterative algorithm, we need to create a working directory for the state. */
549 bufstate = init_bufstate(state);
553 snew(xcopy,state->natoms);
554 snew(vcopy,state->natoms);
555 copy_rvecn(state->x,xcopy,0,state->natoms);
556 copy_rvecn(state->v,vcopy,0,state->natoms);
557 copy_mat(state->box,boxcopy);
560 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
561 temperature control */
562 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
566 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
569 "RMS relative constraint deviation after constraining: %.2e\n",
570 constr_rmsd(constr,FALSE));
572 if (EI_STATE_VELOCITY(ir->eI))
574 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
578 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
579 " input trajectory '%s'\n\n",
580 *(top_global->name),opt2fn("-rerun",nfile,fnm));
583 fprintf(stderr,"Calculated time to finish depends on nsteps from "
584 "run input file,\nwhich may not correspond to the time "
585 "needed to process input trajectory.\n\n");
591 fprintf(stderr,"starting mdrun '%s'\n",
592 *(top_global->name));
595 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
599 sprintf(tbuf,"%s","infinite");
601 if (ir->init_step > 0)
603 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
604 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
605 gmx_step_str(ir->init_step,sbuf2),
606 ir->init_step*ir->delta_t);
610 fprintf(stderr,"%s steps, %s ps.\n",
611 gmx_step_str(ir->nsteps,sbuf),tbuf);
617 /* Set and write start time */
618 runtime_start(runtime);
619 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
620 wallcycle_start(wcycle,ewcRUN);
626 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
628 chkpt_ret=fcCheckPointParallel( cr->nodeid,
630 if ( chkpt_ret == 0 )
631 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
635 /***********************************************************
639 ************************************************************/
641 /* if rerunMD then read coordinates and velocities from input trajectory */
644 if (getenv("GMX_FORCE_UPDATE"))
652 bNotLastFrame = read_first_frame(oenv,&status,
653 opt2fn("-rerun",nfile,fnm),
654 &rerun_fr,TRX_NEED_X | TRX_READ_V);
655 if (rerun_fr.natoms != top_global->natoms)
658 "Number of atoms in trajectory (%d) does not match the "
659 "run input file (%d)\n",
660 rerun_fr.natoms,top_global->natoms);
662 if (ir->ePBC != epbcNONE)
666 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);
668 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
670 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
677 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
680 if (ir->ePBC != epbcNONE)
682 /* Set the shift vectors.
683 * Necessary here when have a static box different from the tpr box.
685 calc_shifts(rerun_fr.box,fr->shift_vec);
689 /* loop over MD steps or if rerunMD to end of input trajectory */
691 /* Skip the first Nose-Hoover integration when we get the state from tpx */
692 bStateFromTPX = !bStateFromCP;
693 bInitStep = bFirstStep && (bStateFromTPX || bVV);
694 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
696 bSumEkinhOld = FALSE;
699 init_global_signals(&gs,cr,ir,repl_ex_nst);
701 step = ir->init_step;
704 if (ir->nstlist == -1)
706 init_nlistheuristics(&nlh,bGStatEveryStep,step);
709 if (MULTISIM(cr) && (repl_ex_nst <=0 ))
711 /* check how many steps are left in other sims */
712 multisim_nsteps=get_multisim_nsteps(cr, ir->nsteps);
716 /* and stop now if we should */
717 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
718 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
719 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
721 wallcycle_start(wcycle,ewcSTEP);
724 if (rerun_fr.bStep) {
725 step = rerun_fr.step;
726 step_rel = step - ir->init_step;
728 if (rerun_fr.bTime) {
738 bLastStep = (step_rel == ir->nsteps);
739 t = t0 + step*ir->delta_t;
742 if (ir->efep != efepNO || ir->bSimTemp)
744 /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
745 requiring different logic. */
747 set_current_lambdas(step,ir->fepvals,bRerunMD,&rerun_fr,state_global,state,lam0);
748 bDoDHDL = do_per_step(step,ir->fepvals->nstdhdl);
749 bDoFEP = (do_per_step(step,nstfep) && (ir->efep != efepNO));
750 bDoExpanded = (do_per_step(step,ir->expandedvals->nstexpanded) && (ir->bExpanded) && (step > 0));
755 update_annealing_target_temp(&(ir->opts),t);
760 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
762 for(i=0; i<state_global->natoms; i++)
764 copy_rvec(rerun_fr.x[i],state_global->x[i]);
768 for(i=0; i<state_global->natoms; i++)
770 copy_rvec(rerun_fr.v[i],state_global->v[i]);
775 for(i=0; i<state_global->natoms; i++)
777 clear_rvec(state_global->v[i]);
781 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
782 " Ekin, temperature and pressure are incorrect,\n"
783 " the virial will be incorrect when constraints are present.\n"
785 bRerunWarnNoV = FALSE;
789 copy_mat(rerun_fr.box,state_global->box);
790 copy_mat(state_global->box,state->box);
792 if (vsite && (Flags & MD_RERUN_VSITE))
794 if (DOMAINDECOMP(cr))
796 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
800 /* Following is necessary because the graph may get out of sync
801 * with the coordinates if we only have every N'th coordinate set
803 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
804 shift_self(graph,state->box,state->x);
806 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
807 top->idef.iparams,top->idef.il,
808 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
811 unshift_self(graph,state->box,state->x);
816 /* Stop Center of Mass motion */
817 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
819 /* Copy back starting coordinates in case we're doing a forcefield scan */
822 for(ii=0; (ii<state->natoms); ii++)
824 copy_rvec(xcopy[ii],state->x[ii]);
825 copy_rvec(vcopy[ii],state->v[ii]);
827 copy_mat(boxcopy,state->box);
832 /* for rerun MD always do Neighbour Searching */
833 bNS = (bFirstStep || ir->nstlist != 0);
838 /* Determine whether or not to do Neighbour Searching and LR */
839 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
841 bNS = (bFirstStep || bExchanged || bNStList || bDoFEP ||
842 (ir->nstlist == -1 && nlh.nabnsb > 0));
844 if (bNS && ir->nstlist == -1)
846 set_nlistheuristics(&nlh,bFirstStep || bExchanged || bDoFEP, step);
850 /* check whether we should stop because another simulation has
854 if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
855 (multisim_nsteps != ir->nsteps) )
862 "Stopping simulation %d because another one has finished\n",
866 gs.sig[eglsCHKPT] = 1;
871 /* < 0 means stop at next step, > 0 means stop at next NS step */
872 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
873 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
878 /* Determine whether or not to update the Born radii if doing GB */
879 bBornRadii=bFirstStep;
880 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
885 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
886 do_verbose = bVerbose &&
887 (step % stepout == 0 || bFirstStep || bLastStep);
889 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
897 bMasterState = FALSE;
898 /* Correct the new box if it is too skewed */
899 if (DYNAMIC_BOX(*ir))
901 if (correct_box(fplog,step,state->box,graph))
906 if (DOMAINDECOMP(cr) && bMasterState)
908 dd_collect_state(cr->dd,state,state_global);
912 if (DOMAINDECOMP(cr))
914 /* Repartition the domain decomposition */
915 wallcycle_start(wcycle,ewcDOMDEC);
916 dd_partition_system(fplog,step,cr,
917 bMasterState,nstglobalcomm,
918 state_global,top_global,ir,
919 state,&f,mdatoms,top,fr,
920 vsite,shellfc,constr,
921 nrnb,wcycle,do_verbose);
922 wallcycle_stop(wcycle,ewcDOMDEC);
923 /* If using an iterative integrator, reallocate space to match the decomposition */
927 if (MASTER(cr) && do_log && !bFFscan)
929 print_ebin_header(fplog,step,t,state->lambda[efptFEP]); /* can we improve the information printed here? */
932 if (ir->efep != efepNO)
934 update_mdatoms(mdatoms,state->lambda[efptMASS]);
937 if (bRerunMD && rerun_fr.bV)
940 /* We need the kinetic energy at minus the half step for determining
941 * the full step kinetic energy and possibly for T-coupling.*/
942 /* This may not be quite working correctly yet . . . . */
943 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
944 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
945 constr,NULL,FALSE,state->box,
946 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
947 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
949 clear_mat(force_vir);
951 /* Ionize the atoms if necessary */
954 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
955 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
958 /* Update force field in ffscan program */
961 if (update_forcefield(fplog,
963 mdatoms->nr,state->x,state->box))
971 /* We write a checkpoint at this MD step when:
972 * either at an NS step when we signalled through gs,
973 * or at the last step (but not when we do not want confout),
974 * but never at the first step or with rerun.
976 bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
977 (bLastStep && (Flags & MD_CONFOUT))) &&
978 step > ir->init_step && !bRerunMD);
981 gs.set[eglsCHKPT] = 0;
984 /* Determine the energy and pressure:
985 * at nstcalcenergy steps and at energy output steps (set below).
988 if (EI_VV(ir->eI) && (!bInitStep)) { /* for vv, the first half actually corresponds to the last step */
989 bNstEner = do_per_step(step-1,ir->nstcalcenergy);
991 bNstEner = do_per_step(step,ir->nstcalcenergy);
995 (ir->epc > epcNO && do_per_step(step,ir->nstpcouple)));
997 /* Do we need global communication ? */
998 bGStat = (bCalcEnerPres || bStopCM ||
999 do_per_step(step,nstglobalcomm) ||
1000 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1002 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1004 if (do_ene || do_log)
1006 bCalcEnerPres = TRUE;
1010 /* these CGLO_ options remain the same throughout the iteration */
1011 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1012 (bGStat ? CGLO_GSTAT : 0)
1015 force_flags = (GMX_FORCE_STATECHANGED |
1016 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1017 GMX_FORCE_ALLFORCES |
1018 (bNStList ? GMX_FORCE_DOLR : 0) |
1020 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
1021 (bDoFEP ? GMX_FORCE_DHDL : 0)
1026 /* Now is the time to relax the shells */
1027 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
1029 bStopCM,top,top_global,
1031 state,f,force_vir,mdatoms,
1032 nrnb,wcycle,graph,groups,
1033 shellfc,fr,bBornRadii,t,mu_tot,
1034 state->natoms,&bConverged,vsite,
1045 /* The coordinates (x) are shifted (to get whole molecules)
1047 * This is parallellized as well, and does communication too.
1048 * Check comments in sim_util.c
1050 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1051 state->box,state->x,&state->hist,
1052 f,force_vir,mdatoms,enerd,fcd,
1053 state->lambda,graph,
1054 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1055 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1060 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1061 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1064 if (bTCR && bFirstStep)
1066 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1067 fprintf(fplog,"Done init_coupling\n");
1071 if (bVV && !bStartingFromCpt && !bRerunMD)
1072 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1074 if (ir->eI==eiVV && bInitStep)
1076 /* if using velocity verlet with full time step Ekin,
1077 * take the first half step only to compute the
1078 * virial for the first step. From there,
1079 * revert back to the initial coordinates
1080 * so that the input is actually the initial step.
1082 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1084 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1085 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1088 update_coords(fplog,step,ir,mdatoms,state,
1089 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1090 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1091 cr,nrnb,constr,&top->idef);
1095 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1097 /* for iterations, we save these vectors, as we will be self-consistently iterating
1100 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1102 /* save the state */
1103 if (bIterations && iterate.bIterate) {
1104 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1107 bFirstIterate = TRUE;
1108 while (bFirstIterate || (bIterations && iterate.bIterate))
1110 if (bIterations && iterate.bIterate)
1112 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1113 if (bFirstIterate && bTrotter)
1115 /* The first time through, we need a decent first estimate
1116 of veta(t+dt) to compute the constraints. Do
1117 this by computing the box volume part of the
1118 trotter integration at this time. Nothing else
1119 should be changed by this routine here. If
1120 !(first time), we start with the previous value
1123 veta_save = state->veta;
1124 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1125 vetanew = state->veta;
1126 state->veta = veta_save;
1131 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
1134 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1135 &top->idef,shake_vir,NULL,
1136 cr,nrnb,wcycle,upd,constr,
1137 bInitStep,TRUE,bCalcEnerPres,vetanew);
1139 if (!bOK && !bFFscan)
1141 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1146 { /* Need to unshift here if a do_force has been
1147 called in the previous step */
1148 unshift_self(graph,state->box,state->x);
1152 /* if VV, compute the pressure and constraints */
1153 /* For VV2, we strictly only need this if using pressure
1154 * control, but we really would like to have accurate pressures
1156 * Think about ways around this in the future?
1157 * For now, keep this choice in comments.
1159 /* bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1160 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1162 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
1163 if (bNstEner && ir->eI==eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/
1165 bSumEkinhOld = TRUE;
1167 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1168 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1169 constr,NULL,FALSE,state->box,
1170 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1173 | (bStopCM ? CGLO_STOPCM : 0)
1174 | (bTemp ? CGLO_TEMPERATURE:0)
1175 | (bPres ? CGLO_PRESSURE : 0)
1176 | (bPres ? CGLO_CONSTRAINT : 0)
1177 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
1178 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1181 /* explanation of above:
1182 a) We compute Ekin at the full time step
1183 if 1) we are using the AveVel Ekin, and it's not the
1184 initial step, or 2) if we are using AveEkin, but need the full
1185 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1186 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1187 EkinAveVel because it's needed for the pressure */
1189 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1194 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
1198 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1203 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1204 state->veta,&vetanew))
1208 bFirstIterate = FALSE;
1211 if (bTrotter && !bInitStep) {
1212 enerd->term[F_DVDL_BONDED] += dvdl; /* only add after iterations */
1213 copy_mat(shake_vir,state->svir_prev);
1214 copy_mat(force_vir,state->fvir_prev);
1215 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
1216 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1217 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
1218 enerd->term[F_EKIN] = trace(ekind->ekin);
1221 /* if it's the initial step, we performed this first step just to get the constraint virial */
1222 if (bInitStep && ir->eI==eiVV) {
1223 copy_rvecn(cbuf,state->v,0,state->natoms);
1226 if (fr->bSepDVDL && fplog && do_log)
1228 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1230 enerd->term[F_DVDL_BONDED] += dvdl;
1233 /* MRS -- now done iterating -- compute the conserved quantity */
1235 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
1238 last_ekin = enerd->term[F_EKIN];
1240 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1242 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1244 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1245 sum_dhdl(enerd,state->lambda,ir->fepvals);
1248 /* ######## END FIRST UPDATE STEP ############## */
1249 /* ######## If doing VV, we now have v(dt) ###### */
1251 /* perform extended ensemble sampling in lambda - we don't
1252 actually move to the new state before outputting
1253 statistics, but if performing simulated tempering, we
1254 do update the velocities and the tau_t. */
1256 lamnew = ExpandedEnsembleDynamics(fplog,ir,enerd,state,&MassQ,&df_history,step,mcrng,state->v,mdatoms);
1258 /* ################## START TRAJECTORY OUTPUT ################# */
1260 /* Now we have the energies and forces corresponding to the
1261 * coordinates at time t. We must output all of this before
1263 * for RerunMD t is read from input trajectory
1266 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
1267 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
1268 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
1269 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
1270 if (bCPT) { mdof_flags |= MDOF_CPT; };
1272 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1275 /* Enforce writing positions and velocities at end of run */
1276 mdof_flags |= (MDOF_X | MDOF_V);
1281 fcReportProgress( ir->nsteps, step );
1283 /* sync bCPT and fc record-keeping */
1284 if (bCPT && MASTER(cr))
1285 fcRequestCheckPoint();
1288 if (mdof_flags != 0)
1290 wallcycle_start(wcycle,ewcTRAJ);
1293 if (state->flags & (1<<estLD_RNG))
1295 get_stochd_state(upd,state);
1297 if (state->flags & (1<<estMC_RNG))
1299 get_mc_state(mcrng,state);
1305 state_global->ekinstate.bUpToDate = FALSE;
1309 update_ekinstate(&state_global->ekinstate,ekind);
1310 state_global->ekinstate.bUpToDate = TRUE;
1312 update_energyhistory(&state_global->enerhist,mdebin);
1313 if (ir->efep!=efepNO || ir->bSimTemp)
1315 state_global->fep_state = state->fep_state; /* MRS: seems kludgy. The code should be
1316 structured so this isn't necessary.
1317 Note this reassignment is only necessary
1318 for single threads.*/
1319 copy_df_history(&state_global->dfhist,&df_history);
1323 write_traj(fplog,cr,outf,mdof_flags,top_global,
1324 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
1331 if (bLastStep && step_rel == ir->nsteps &&
1332 (Flags & MD_CONFOUT) && MASTER(cr) &&
1333 !bRerunMD && !bFFscan)
1335 /* x and v have been collected in write_traj,
1336 * because a checkpoint file will always be written
1339 fprintf(stderr,"\nWriting final coordinates.\n");
1340 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
1343 /* Make molecules whole only for confout writing */
1344 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
1346 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
1347 *top_global->name,top_global,
1348 state_global->x,state_global->v,
1349 ir->ePBC,state->box);
1352 wallcycle_stop(wcycle,ewcTRAJ);
1355 /* kludge -- virial is lost with restart for NPT control. Must restart */
1356 if (bStartingFromCpt && bVV)
1358 copy_mat(state->svir_prev,shake_vir);
1359 copy_mat(state->fvir_prev,force_vir);
1361 /* ################## END TRAJECTORY OUTPUT ################ */
1363 /* Determine the pressure:
1364 * always when we want exact averages in the energy file,
1365 * at ns steps when we have pressure coupling,
1366 * otherwise only at energy output steps (set below).
1370 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
1371 bCalcEnerPres = bNstEner;
1373 /* Do we need global communication ? */
1374 bGStat = (bGStatEveryStep || bStopCM || bNS ||
1375 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1377 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1379 if (do_ene || do_log)
1381 bCalcEnerPres = TRUE;
1385 /* Determine the wallclock run time up till now */
1386 run_time = gmx_gettime() - (double)runtime->real;
1387 /* Check whether everything is still allright */
1388 if (((int)gmx_get_stop_condition() > handled_stop_condition)
1389 #ifdef GMX_THREAD_MPI
1394 /* this is just make gs.sig compatible with the hack
1395 of sending signals around by MPI_Reduce with together with
1397 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
1398 gs.sig[eglsSTOPCOND]=1;
1399 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
1400 gs.sig[eglsSTOPCOND]=-1;
1401 /* < 0 means stop at next step, > 0 means stop at next NS step */
1405 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1406 gmx_get_signal_name(),
1407 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1411 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1412 gmx_get_signal_name(),
1413 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1415 handled_stop_condition=(int)gmx_get_stop_condition();
1417 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1418 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
1419 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
1421 /* Signal to terminate the run */
1422 gs.sig[eglsSTOPCOND] = 1;
1425 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1427 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1430 if (bResetCountersHalfMaxH && MASTER(cr) &&
1431 run_time > max_hours*60.0*60.0*0.495)
1433 gs.sig[eglsRESETCOUNTERS] = 1;
1436 if (ir->nstlist == -1 && !bRerunMD)
1438 /* When bGStatEveryStep=FALSE, global_stat is only called
1439 * when we check the atom displacements, not at NS steps.
1440 * This means that also the bonded interaction count check is not
1441 * performed immediately after NS. Therefore a few MD steps could
1442 * be performed with missing interactions.
1443 * But wrong energies are never written to file,
1444 * since energies are only written after global_stat
1447 if (step >= nlh.step_nscheck)
1449 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
1450 nlh.scale_tot,state->x);
1454 /* This is not necessarily true,
1455 * but step_nscheck is determined quite conservatively.
1461 /* In parallel we only have to check for checkpointing in steps
1462 * where we do global communication,
1463 * otherwise the other nodes don't know.
1465 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1468 run_time >= nchkpt*cpt_period*60.0)) &&
1469 gs.set[eglsCHKPT] == 0)
1471 gs.sig[eglsCHKPT] = 1;
1475 /* at the start of step, randomize the velocities */
1476 if (ETC_ANDERSEN(ir->etc) && EI_VV(ir->eI))
1478 gmx_bool bDoAndersenConstr;
1479 bDoAndersenConstr = update_randomize_velocities(ir,step,mdatoms,state,upd,&top->idef,constr);
1480 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1481 if (bDoAndersenConstr)
1483 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1484 &top->idef,tmp_vir,NULL,
1485 cr,nrnb,wcycle,upd,constr,
1486 bInitStep,TRUE,FALSE,vetanew);
1492 gmx_iterate_init(&iterate,bIterations);
1495 /* for iterations, we save these vectors, as we will be redoing the calculations */
1496 if (bIterations && iterate.bIterate)
1498 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1500 bFirstIterate = TRUE;
1501 while (bFirstIterate || (bIterations && iterate.bIterate))
1503 /* We now restore these vectors to redo the calculation with improved extended variables */
1506 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1509 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1510 so scroll down for that logic */
1512 /* ######### START SECOND UPDATE STEP ################# */
1513 /* Box is changed in update() when we do pressure coupling,
1514 * but we should still use the old box for energy corrections and when
1515 * writing it to the energy file, so it matches the trajectory files for
1516 * the same timestep above. Make a copy in a separate array.
1518 copy_mat(state->box,lastbox);
1521 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
1523 wallcycle_start(wcycle,ewcUPDATE);
1525 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1528 if (bIterations && iterate.bIterate)
1536 /* we use a new value of scalevir to converge the iterations faster */
1537 scalevir = tracevir/trace(shake_vir);
1539 msmul(shake_vir,scalevir,shake_vir);
1540 m_add(force_vir,shake_vir,total_vir);
1541 clear_mat(shake_vir);
1543 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
1544 /* We can only do Berendsen coupling after we have summed
1545 * the kinetic energy or virial. Since the happens
1546 * in global_state after update, we should only do it at
1547 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1552 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1553 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
1559 /* velocity half-step update */
1560 update_coords(fplog,step,ir,mdatoms,state,f,
1561 fr->bTwinRange && bNStList,fr->f_twin,fcd,
1562 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
1563 cr,nrnb,constr,&top->idef);
1566 /* Above, initialize just copies ekinh into ekin,
1567 * it doesn't copy position (for VV),
1568 * and entire integrator for MD.
1573 copy_rvecn(state->x,cbuf,0,state->natoms);
1576 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1577 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1578 wallcycle_stop(wcycle,ewcUPDATE);
1580 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1581 &top->idef,shake_vir,force_vir,
1582 cr,nrnb,wcycle,upd,constr,
1583 bInitStep,FALSE,bCalcEnerPres,state->veta);
1587 /* erase F_EKIN and F_TEMP here? */
1588 /* just compute the kinetic energy at the half step to perform a trotter step */
1589 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1590 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1591 constr,NULL,FALSE,lastbox,
1592 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1593 cglo_flags | CGLO_TEMPERATURE
1595 wallcycle_start(wcycle,ewcUPDATE);
1596 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
1597 /* now we know the scaling, we can compute the positions again again */
1598 copy_rvecn(cbuf,state->x,0,state->natoms);
1600 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1601 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1602 wallcycle_stop(wcycle,ewcUPDATE);
1604 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1605 /* are the small terms in the shake_vir here due
1606 * to numerical errors, or are they important
1607 * physically? I'm thinking they are just errors, but not completely sure.
1608 * For now, will call without actually constraining, constr=NULL*/
1609 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
1610 &top->idef,tmp_vir,force_vir,
1611 cr,nrnb,wcycle,upd,NULL,
1612 bInitStep,FALSE,bCalcEnerPres,
1615 if (!bOK && !bFFscan)
1617 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1620 if (fr->bSepDVDL && fplog && do_log)
1622 fprintf(fplog,sepdvdlformat,"Constraint dV/dl",0.0,dvdl);
1624 enerd->term[F_DVDL_BONDED] += dvdl;
1628 /* Need to unshift here */
1629 unshift_self(graph,state->box,state->x);
1634 wallcycle_start(wcycle,ewcVSITECONSTR);
1637 shift_self(graph,state->box,state->x);
1639 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1640 top->idef.iparams,top->idef.il,
1641 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1645 unshift_self(graph,state->box,state->x);
1647 wallcycle_stop(wcycle,ewcVSITECONSTR);
1650 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1651 if (ir->nstlist == -1 && bFirstIterate)
1653 gs.sig[eglsNABNSB] = nlh.nabnsb;
1655 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. */
1656 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1657 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1659 bFirstIterate ? &gs : NULL,
1660 (step_rel % gs.nstms == 0) &&
1661 (multisim_nsteps<0 || (step_rel<multisim_nsteps)),
1663 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1665 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
1666 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1667 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1668 | (bEnergyHere || bRerunMD ? CGLO_PRESSURE : 0)
1669 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
1670 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1673 if (ir->nstlist == -1 && bFirstIterate)
1675 nlh.nabnsb = gs.set[eglsNABNSB];
1676 gs.set[eglsNABNSB] = 0;
1678 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1679 /* ############# END CALC EKIN AND PRESSURE ################# */
1681 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1682 the virial that should probably be addressed eventually. state->veta has better properies,
1683 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1684 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1687 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1688 trace(shake_vir),&tracevir))
1692 bFirstIterate = FALSE;
1695 /* only add constraint dvdl after constraints */
1696 enerd->term[F_DVDL_BONDED] += dvdl;
1699 /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */
1700 sum_dhdl(enerd,state->lambda,ir->fepvals);
1702 update_box(fplog,step,ir,mdatoms,state,graph,f,
1703 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
1705 /* ################# END UPDATE STEP 2 ################# */
1706 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1708 /* The coordinates (x) were unshifted in update */
1709 if (bFFscan && (shellfc==NULL || bConverged))
1711 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
1713 &(top_global->mols),mdatoms->massT,pres))
1717 fprintf(stderr,"\n");
1723 /* We will not sum ekinh_old,
1724 * so signal that we still have to do it.
1726 bSumEkinhOld = TRUE;
1731 /* Only do GCT when the relaxation of shells (minimization) has converged,
1732 * otherwise we might be coupling to bogus energies.
1733 * In parallel we must always do this, because the other sims might
1737 /* Since this is called with the new coordinates state->x, I assume
1738 * we want the new box state->box too. / EL 20040121
1740 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
1742 mdatoms,&(top->idef),mu_aver,
1743 top_global->mols.nr,cr,
1744 state->box,total_vir,pres,
1745 mu_tot,state->x,f,bConverged);
1749 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1751 /* use the directly determined last velocity, not actually the averaged half steps */
1752 if (bTrotter && ir->eI==eiVV)
1754 enerd->term[F_EKIN] = last_ekin;
1756 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1760 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1764 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
1766 /* Check for excessively large energies */
1770 real etot_max = 1e200;
1772 real etot_max = 1e30;
1774 if (fabs(enerd->term[F_ETOT]) > etot_max)
1776 fprintf(stderr,"Energy too large (%g), giving up\n",
1777 enerd->term[F_ETOT]);
1780 /* ######### END PREPARING EDR OUTPUT ########### */
1782 /* Time for performance */
1783 if (((step % stepout) == 0) || bLastStep)
1785 runtime_upd_proc(runtime);
1791 gmx_bool do_dr,do_or;
1793 if (fplog && do_log && bDoExpanded)
1795 /* only needed if doing expanded ensemble */
1796 PrintFreeEnergyInfoToFile(fplog,ir->fepvals,ir->expandedvals,ir->bSimTemp?ir->simtempvals:NULL,
1797 &df_history,state->fep_state,ir->nstlog,step);
1799 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
1803 upd_mdebin(mdebin,bDoDHDL,TRUE,
1804 t,mdatoms->tmass,enerd,state,
1805 ir->fepvals,ir->expandedvals,lastbox,
1806 shake_vir,force_vir,total_vir,pres,
1807 ekind,mu_tot,constr);
1811 upd_mdebin_step(mdebin);
1814 do_dr = do_per_step(step,ir->nstdisreout);
1815 do_or = do_per_step(step,ir->nstorireout);
1817 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
1819 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
1821 if (ir->ePull != epullNO)
1823 pull_print_output(ir->pull,step,t);
1826 if (do_per_step(step,ir->nstlog))
1828 if(fflush(fplog) != 0)
1830 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
1836 /* Have to do this part after outputting the logfile and the edr file */
1837 state->fep_state = lamnew;
1838 for (i=0;i<efptNR;i++)
1840 state->lambda[i] = ir->fepvals->all_lambda[i][lamnew];
1843 /* Remaining runtime */
1844 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
1848 fprintf(stderr,"\n");
1850 print_time(stderr,runtime,step,ir,cr);
1853 /* Replica exchange */
1855 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
1856 do_per_step(step,repl_ex_nst))
1858 bExchanged = replica_exchange(fplog,cr,repl_ex,
1862 if (bExchanged && DOMAINDECOMP(cr))
1864 dd_partition_system(fplog,step,cr,TRUE,1,
1865 state_global,top_global,ir,
1866 state,&f,mdatoms,top,fr,
1867 vsite,shellfc,constr,
1874 bStartingFromCpt = FALSE;
1876 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1877 /* With all integrators, except VV, we need to retain the pressure
1878 * at the current step for coupling at the next step.
1880 if ((state->flags & (1<<estPRES_PREV)) &&
1882 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1884 /* Store the pressure in t_state for pressure coupling
1885 * at the next MD step.
1887 copy_mat(pres,state->pres_prev);
1890 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1892 if ( (membed!=NULL) && (!bLastStep) )
1894 rescale_membed(step_rel,membed,state_global->x);
1901 /* read next frame from input trajectory */
1902 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
1907 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
1911 if (!bRerunMD || !rerun_fr.bStep)
1913 /* increase the MD step number */
1918 cycles = wallcycle_stop(wcycle,ewcSTEP);
1919 if (DOMAINDECOMP(cr) && wcycle)
1921 dd_cycles_add(cr->dd,cycles,ddCyclStep);
1924 if (step_rel == wcycle_get_reset_counters(wcycle) ||
1925 gs.set[eglsRESETCOUNTERS] != 0)
1927 /* Reset all the counters related to performance over the run */
1928 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
1929 wcycle_set_reset_counters(wcycle,-1);
1930 /* Correct max_hours for the elapsed time */
1931 max_hours -= run_time/(60.0*60.0);
1932 bResetCountersHalfMaxH = FALSE;
1933 gs.set[eglsRESETCOUNTERS] = 0;
1937 /* End of main MD loop */
1941 runtime_end(runtime);
1943 if (bRerunMD && MASTER(cr))
1948 if (!(cr->duty & DUTY_PME))
1950 /* Tell the PME only node to finish */
1956 if (ir->nstcalcenergy > 0 && !bRerunMD)
1958 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
1959 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
1967 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
1969 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)));
1970 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
1973 if (shellfc && fplog)
1975 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
1976 (nconverged*100.0)/step_rel);
1977 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
1981 if (repl_ex_nst > 0 && MASTER(cr))
1983 print_replica_exchange_statistics(fplog,repl_ex);
1986 runtime->nsteps_done = step_rel;