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
56 #include "md_support.h"
57 #include "md_logging.h"
72 #include "mpelogging.h"
74 #include "domdec_network.h"
80 #include "compute_io.h"
82 #include "checkpoint.h"
83 #include "mtop_util.h"
84 #include "sighandler.h"
87 #include "pme_loadbal.h"
90 #include "types/nlistheuristics.h"
91 #include "types/iteratedconstraints.h"
92 #include "nbnxn_cuda_data_mgmt.h"
102 #include "corewrap.h"
105 static void reset_all_counters(FILE *fplog,t_commrec *cr,
106 gmx_large_int_t step,
107 gmx_large_int_t *step_rel,t_inputrec *ir,
108 gmx_wallcycle_t wcycle,t_nrnb *nrnb,
109 gmx_runtime_t *runtime,
110 nbnxn_cuda_ptr_t cu_nbv)
112 char sbuf[STEPSTRSIZE];
114 /* Reset all the counters related to performance over the run */
115 md_print_warn(cr,fplog,"step %s: resetting all time and cycle counters\n",
116 gmx_step_str(step,sbuf));
120 nbnxn_cuda_reset_timings(cu_nbv);
123 wallcycle_stop(wcycle,ewcRUN);
124 wallcycle_reset_all(wcycle);
125 if (DOMAINDECOMP(cr))
127 reset_dd_statistics_counters(cr->dd);
130 ir->init_step += *step_rel;
131 ir->nsteps -= *step_rel;
133 wallcycle_start(wcycle,ewcRUN);
134 runtime_start(runtime);
135 print_date_and_time(fplog,cr->nodeid,"Restarted time",runtime);
138 double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
139 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
141 gmx_vsite_t *vsite,gmx_constr_t constr,
142 int stepout,t_inputrec *ir,
143 gmx_mtop_t *top_global,
145 t_state *state_global,
147 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
148 gmx_edsam_t ed,t_forcerec *fr,
149 int repl_ex_nst,int repl_ex_nex,int repl_ex_seed,gmx_membed_t membed,
150 real cpt_period,real max_hours,
151 const char *deviceOptions,
153 gmx_runtime_t *runtime)
156 gmx_large_int_t step,step_rel;
158 double t,t0,lam0[efptNR];
159 gmx_bool bGStatEveryStep,bGStat,bCalcVir,bCalcEner;
160 gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
161 bFirstStep,bStateFromCP,bStateFromTPX,bInitStep,bLastStep,
162 bBornRadii,bStartingFromCpt;
163 gmx_bool bDoDHDL=FALSE,bDoFEP=FALSE,bDoExpanded=FALSE;
164 gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
165 bForceUpdate=FALSE,bCPT;
167 gmx_bool bMasterState;
168 int force_flags,cglo_flags;
169 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
174 t_state *bufstate=NULL;
175 matrix *scale_tot,pcoupl_mu,M,ebox;
178 gmx_repl_ex_t repl_ex=NULL;
181 t_mdebin *mdebin=NULL;
182 df_history_t df_history;
187 gmx_enerdata_t *enerd;
189 gmx_global_stat_t gstat;
190 gmx_update_t upd=NULL;
193 gmx_rng_t mcrng=NULL;
195 gmx_groups_t *groups;
196 gmx_ekindata_t *ekind, *ekind_save;
197 gmx_shellfc_t shellfc;
198 int count,nconverged=0;
201 gmx_bool bIonize=FALSE;
202 gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
204 gmx_bool bResetCountersHalfMaxH=FALSE;
205 gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
206 gmx_bool bUpdateDoLR;
209 atom_id *grpindex=NULL;
211 t_coupl_rec *tcr=NULL;
212 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
213 matrix boxcopy={{0}},lastbox;
215 real fom,oldfom,veta_save,pcurr,scalevir,tracevir;
223 real saved_conserved_quantity = 0;
228 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
229 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
230 gmx_iterate_t iterate;
231 gmx_large_int_t multisim_nsteps=-1; /* number of steps to do before first multisim
232 simulation stops. If equal to zero, don't
233 communicate any more between multisims.*/
234 /* PME load balancing data for GPU kernels */
235 pme_load_balancing_t pme_loadbal=NULL;
237 gmx_bool bPMETuneTry=FALSE,bPMETuneRunning=FALSE;
240 /* Temporary addition for FAHCORE checkpointing */
244 /* Check for special mdrun options */
245 bRerunMD = (Flags & MD_RERUN);
246 bIonize = (Flags & MD_IONIZE);
247 bFFscan = (Flags & MD_FFSCAN);
248 bAppend = (Flags & MD_APPENDFILES);
249 if (Flags & MD_RESETCOUNTERSHALFWAY)
253 /* Signal to reset the counters half the simulation steps. */
254 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
256 /* Signal to reset the counters halfway the simulation time. */
257 bResetCountersHalfMaxH = (max_hours > 0);
260 /* md-vv uses averaged full step velocities for T-control
261 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
262 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
264 if (bVV) /* to store the initial velocities while computing virial */
266 snew(cbuf,top_global->natoms);
268 /* all the iteratative cases - only if there are constraints */
269 bIterations = ((IR_NPH_TROTTER(ir) || IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
270 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir) || IR_NVT_TROTTER(ir)));
274 /* Since we don't know if the frames read are related in any way,
275 * rebuild the neighborlist at every step.
278 ir->nstcalcenergy = 1;
282 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
284 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
285 bGStatEveryStep = (nstglobalcomm == 1);
287 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
290 "To reduce the energy communication with nstlist = -1\n"
291 "the neighbor list validity should not be checked at every step,\n"
292 "this means that exact integration is not guaranteed.\n"
293 "The neighbor list validity is checked after:\n"
294 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
295 "In most cases this will result in exact integration.\n"
296 "This reduces the energy communication by a factor of 2 to 3.\n"
297 "If you want less energy communication, set nstlist > 3.\n\n");
300 if (bRerunMD || bFFscan)
304 groups = &top_global->groups;
307 init_md(fplog,cr,ir,oenv,&t,&t0,state_global->lambda,
308 &(state_global->fep_state),lam0,
309 nrnb,top_global,&upd,
310 nfile,fnm,&outf,&mdebin,
311 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
313 clear_mat(total_vir);
315 /* Energy terms and groups */
317 init_enerdata(top_global->groups.grps[egcENER].nr,ir->fepvals->n_lambda,
319 if (DOMAINDECOMP(cr))
325 snew(f,top_global->natoms);
328 /* lambda Monte carlo random number generator */
331 mcrng = gmx_rng_init(ir->expandedvals->lmc_seed);
333 /* copy the state into df_history */
334 copy_df_history(&df_history,&state_global->dfhist);
336 /* Kinetic energy data */
338 init_ekindata(fplog,top_global,&(ir->opts),ekind);
339 /* needed for iteration of constraints */
341 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
342 /* Copy the cos acceleration to the groups struct */
343 ekind->cosacc.cos_accel = ir->cos_accel;
345 gstat = global_stat_init(ir);
348 /* Check for polarizable models and flexible constraints */
349 shellfc = init_shell_flexcon(fplog,
350 top_global,n_flexible_constraints(constr),
351 (ir->bContinuation ||
352 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
353 NULL : state_global->x);
357 #ifdef GMX_THREAD_MPI
358 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
360 set_deform_reference_box(upd,
361 deform_init_init_step_tpx,
362 deform_init_box_tpx);
363 #ifdef GMX_THREAD_MPI
364 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
369 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
370 if ((io > 2000) && MASTER(cr))
372 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
376 if (DOMAINDECOMP(cr)) {
377 top = dd_init_local_top(top_global);
380 dd_init_local_state(cr->dd,state_global,state);
382 if (DDMASTER(cr->dd) && ir->nstfout) {
383 snew(f_global,state_global->natoms);
387 /* Initialize the particle decomposition and split the topology */
388 top = split_system(fplog,top_global,ir,cr);
390 pd_cg_range(cr,&fr->cg0,&fr->hcg);
391 pd_at_range(cr,&a0,&a1);
393 top = gmx_mtop_generate_local_top(top_global,ir);
396 a1 = top_global->natoms;
399 forcerec_set_excl_load(fr,top,cr);
401 state = partdec_init_local_state(cr,state_global);
404 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
407 set_vsite_top(vsite,top,mdatoms,cr);
410 if (ir->ePBC != epbcNONE && !fr->bMolPBC) {
411 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
415 make_local_shells(cr,mdatoms,shellfc);
418 init_bonded_thread_force_reduction(fr,&top->idef);
420 if (ir->pull && PAR(cr)) {
421 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
425 if (DOMAINDECOMP(cr))
427 /* Distribute the charge groups over the nodes from the master node */
428 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
429 state_global,top_global,ir,
430 state,&f,mdatoms,top,fr,
431 vsite,shellfc,constr,
436 update_mdatoms(mdatoms,state->lambda[efptMASS]);
438 if (opt2bSet("-cpi",nfile,fnm))
440 bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr);
444 bStateFromCP = FALSE;
451 /* Update mdebin with energy history if appending to output files */
452 if ( Flags & MD_APPENDFILES )
454 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
458 /* We might have read an energy history from checkpoint,
459 * free the allocated memory and reset the counts.
461 done_energyhistory(&state_global->enerhist);
462 init_energyhistory(&state_global->enerhist);
465 /* Set the initial energy history in state by updating once */
466 update_energyhistory(&state_global->enerhist,mdebin);
469 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG))
471 /* Set the random state if we read a checkpoint file */
472 set_stochd_state(upd,state);
475 if (state->flags & (1<<estMC_RNG))
477 set_mc_state(mcrng,state);
480 /* Initialize constraints */
482 if (!DOMAINDECOMP(cr))
483 set_constraints(constr,top,ir,mdatoms,cr);
486 /* Check whether we have to GCT stuff */
487 bTCR = ftp2bSet(efGCT,nfile,fnm);
490 fprintf(stderr,"Will do General Coupling Theory!\n");
492 gnx = top_global->mols.nr;
494 for(i=0; (i<gnx); i++) {
501 /* We need to be sure replica exchange can only occur
502 * when the energies are current */
503 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
504 "repl_ex_nst",&repl_ex_nst);
505 /* This check needs to happen before inter-simulation
506 * signals are initialized, too */
508 if (repl_ex_nst > 0 && MASTER(cr))
510 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
511 repl_ex_nst,repl_ex_nex,repl_ex_seed);
514 /* PME tuning is only supported with GPUs or PME nodes and not with rerun */
515 if ((Flags & MD_TUNEPME) &&
516 EEL_PME(fr->eeltype) &&
517 ( (fr->cutoff_scheme == ecutsVERLET && fr->nbv->bUseGPU) || !(cr->duty & DUTY_PME)) &&
520 pme_loadbal_init(&pme_loadbal,ir,state->box,fr->ic,fr->pmedata);
522 if (cr->duty & DUTY_PME)
524 /* Start tuning right away, as we can't measure the load */
525 bPMETuneRunning = TRUE;
529 /* Separate PME nodes, we can measure the PP/PME load balance */
534 if (!ir->bContinuation && !bRerunMD)
536 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
538 /* Set the velocities of frozen particles to zero */
539 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
543 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
553 /* Constrain the initial coordinates and velocities */
554 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
555 graph,cr,nrnb,fr,top,shake_vir);
559 /* Construct the virtual sites for the initial configuration */
560 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
561 top->idef.iparams,top->idef.il,
562 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
568 /* set free energy calculation frequency as the minimum of nstdhdl, nstexpanded, and nstrepl_ex_nst*/
569 nstfep = ir->fepvals->nstdhdl;
570 if (ir->bExpanded && (nstfep > ir->expandedvals->nstexpanded))
572 nstfep = ir->expandedvals->nstexpanded;
574 if (repl_ex_nst > 0 && repl_ex_nst > nstfep)
576 nstfep = repl_ex_nst;
579 /* I'm assuming we need global communication the first time! MRS */
580 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
581 | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM:0)
582 | (bVV ? CGLO_PRESSURE:0)
583 | (bVV ? CGLO_CONSTRAINT:0)
584 | (bRerunMD ? CGLO_RERUNMD:0)
585 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
587 bSumEkinhOld = FALSE;
588 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
589 NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
590 constr,NULL,FALSE,state->box,
591 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
592 if (ir->eI == eiVVAK) {
593 /* a second call to get the half step temperature initialized as well */
594 /* we do the same call as above, but turn the pressure off -- internally to
595 compute_globals, this is recognized as a velocity verlet half-step
596 kinetic energy calculation. This minimized excess variables, but
597 perhaps loses some logic?*/
599 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
600 NULL,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
601 constr,NULL,FALSE,state->box,
602 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
603 cglo_flags &~ (CGLO_STOPCM | CGLO_PRESSURE));
606 /* Calculate the initial half step temperature, and save the ekinh_old */
607 if (!(Flags & MD_STARTFROMCPT))
609 for(i=0; (i<ir->opts.ngtc); i++)
611 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
616 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
617 and there is no previous step */
620 /* if using an iterative algorithm, we need to create a working directory for the state. */
623 bufstate = init_bufstate(state);
627 snew(xcopy,state->natoms);
628 snew(vcopy,state->natoms);
629 copy_rvecn(state->x,xcopy,0,state->natoms);
630 copy_rvecn(state->v,vcopy,0,state->natoms);
631 copy_mat(state->box,boxcopy);
634 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
635 temperature control */
636 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
640 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
643 "RMS relative constraint deviation after constraining: %.2e\n",
644 constr_rmsd(constr,FALSE));
646 if (EI_STATE_VELOCITY(ir->eI))
648 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
652 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
653 " input trajectory '%s'\n\n",
654 *(top_global->name),opt2fn("-rerun",nfile,fnm));
657 fprintf(stderr,"Calculated time to finish depends on nsteps from "
658 "run input file,\nwhich may not correspond to the time "
659 "needed to process input trajectory.\n\n");
665 fprintf(stderr,"starting mdrun '%s'\n",
666 *(top_global->name));
669 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
673 sprintf(tbuf,"%s","infinite");
675 if (ir->init_step > 0)
677 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
678 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
679 gmx_step_str(ir->init_step,sbuf2),
680 ir->init_step*ir->delta_t);
684 fprintf(stderr,"%s steps, %s ps.\n",
685 gmx_step_str(ir->nsteps,sbuf),tbuf);
691 /* Set and write start time */
692 runtime_start(runtime);
693 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
694 wallcycle_start(wcycle,ewcRUN);
700 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
702 chkpt_ret=fcCheckPointParallel( cr->nodeid,
704 if ( chkpt_ret == 0 )
705 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
709 /***********************************************************
713 ************************************************************/
715 /* if rerunMD then read coordinates and velocities from input trajectory */
718 if (getenv("GMX_FORCE_UPDATE"))
726 bNotLastFrame = read_first_frame(oenv,&status,
727 opt2fn("-rerun",nfile,fnm),
728 &rerun_fr,TRX_NEED_X | TRX_READ_V);
729 if (rerun_fr.natoms != top_global->natoms)
732 "Number of atoms in trajectory (%d) does not match the "
733 "run input file (%d)\n",
734 rerun_fr.natoms,top_global->natoms);
736 if (ir->ePBC != epbcNONE)
740 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);
742 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
744 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
751 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
754 if (ir->ePBC != epbcNONE)
756 /* Set the shift vectors.
757 * Necessary here when have a static box different from the tpr box.
759 calc_shifts(rerun_fr.box,fr->shift_vec);
763 /* loop over MD steps or if rerunMD to end of input trajectory */
765 /* Skip the first Nose-Hoover integration when we get the state from tpx */
766 bStateFromTPX = !bStateFromCP;
767 bInitStep = bFirstStep && (bStateFromTPX || bVV);
768 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
770 bSumEkinhOld = FALSE;
773 init_global_signals(&gs,cr,ir,repl_ex_nst);
775 step = ir->init_step;
778 if (ir->nstlist == -1)
780 init_nlistheuristics(&nlh,bGStatEveryStep,step);
783 if (MULTISIM(cr) && (repl_ex_nst <=0 ))
785 /* check how many steps are left in other sims */
786 multisim_nsteps=get_multisim_nsteps(cr, ir->nsteps);
790 /* and stop now if we should */
791 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
792 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
793 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
795 wallcycle_start(wcycle,ewcSTEP);
797 GMX_MPE_LOG(ev_timestep1);
800 if (rerun_fr.bStep) {
801 step = rerun_fr.step;
802 step_rel = step - ir->init_step;
804 if (rerun_fr.bTime) {
814 bLastStep = (step_rel == ir->nsteps);
815 t = t0 + step*ir->delta_t;
818 if (ir->efep != efepNO || ir->bSimTemp)
820 /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
821 requiring different logic. */
823 set_current_lambdas(step,ir->fepvals,bRerunMD,&rerun_fr,state_global,state,lam0);
824 bDoDHDL = do_per_step(step,ir->fepvals->nstdhdl);
825 bDoFEP = (do_per_step(step,nstfep) && (ir->efep != efepNO));
826 bDoExpanded = (do_per_step(step,ir->expandedvals->nstexpanded) && (ir->bExpanded) && (step > 0));
831 update_annealing_target_temp(&(ir->opts),t);
836 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
838 for(i=0; i<state_global->natoms; i++)
840 copy_rvec(rerun_fr.x[i],state_global->x[i]);
844 for(i=0; i<state_global->natoms; i++)
846 copy_rvec(rerun_fr.v[i],state_global->v[i]);
851 for(i=0; i<state_global->natoms; i++)
853 clear_rvec(state_global->v[i]);
857 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
858 " Ekin, temperature and pressure are incorrect,\n"
859 " the virial will be incorrect when constraints are present.\n"
861 bRerunWarnNoV = FALSE;
865 copy_mat(rerun_fr.box,state_global->box);
866 copy_mat(state_global->box,state->box);
868 if (vsite && (Flags & MD_RERUN_VSITE))
870 if (DOMAINDECOMP(cr))
872 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
876 /* Following is necessary because the graph may get out of sync
877 * with the coordinates if we only have every N'th coordinate set
879 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
880 shift_self(graph,state->box,state->x);
882 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
883 top->idef.iparams,top->idef.il,
884 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
887 unshift_self(graph,state->box,state->x);
892 /* Stop Center of Mass motion */
893 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
895 /* Copy back starting coordinates in case we're doing a forcefield scan */
898 for(ii=0; (ii<state->natoms); ii++)
900 copy_rvec(xcopy[ii],state->x[ii]);
901 copy_rvec(vcopy[ii],state->v[ii]);
903 copy_mat(boxcopy,state->box);
908 /* for rerun MD always do Neighbour Searching */
909 bNS = (bFirstStep || ir->nstlist != 0);
914 /* Determine whether or not to do Neighbour Searching and LR */
915 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
917 bNS = (bFirstStep || bExchanged || bNStList || bDoFEP ||
918 (ir->nstlist == -1 && nlh.nabnsb > 0));
920 if (bNS && ir->nstlist == -1)
922 set_nlistheuristics(&nlh,bFirstStep || bExchanged || bDoFEP, step);
926 /* check whether we should stop because another simulation has
930 if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
931 (multisim_nsteps != ir->nsteps) )
938 "Stopping simulation %d because another one has finished\n",
942 gs.sig[eglsCHKPT] = 1;
947 /* < 0 means stop at next step, > 0 means stop at next NS step */
948 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
949 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
954 /* Determine whether or not to update the Born radii if doing GB */
955 bBornRadii=bFirstStep;
956 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
961 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
962 do_verbose = bVerbose &&
963 (step % stepout == 0 || bFirstStep || bLastStep);
965 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
973 bMasterState = FALSE;
974 /* Correct the new box if it is too skewed */
975 if (DYNAMIC_BOX(*ir))
977 if (correct_box(fplog,step,state->box,graph))
982 if (DOMAINDECOMP(cr) && bMasterState)
984 dd_collect_state(cr->dd,state,state_global);
988 if (DOMAINDECOMP(cr))
990 /* Repartition the domain decomposition */
991 wallcycle_start(wcycle,ewcDOMDEC);
992 dd_partition_system(fplog,step,cr,
993 bMasterState,nstglobalcomm,
994 state_global,top_global,ir,
995 state,&f,mdatoms,top,fr,
996 vsite,shellfc,constr,
998 do_verbose && !bPMETuneRunning);
999 wallcycle_stop(wcycle,ewcDOMDEC);
1000 /* If using an iterative integrator, reallocate space to match the decomposition */
1004 if (MASTER(cr) && do_log && !bFFscan)
1006 print_ebin_header(fplog,step,t,state->lambda[efptFEP]); /* can we improve the information printed here? */
1009 if (ir->efep != efepNO)
1011 update_mdatoms(mdatoms,state->lambda[efptMASS]);
1014 if ((bRerunMD && rerun_fr.bV) || bExchanged)
1017 /* We need the kinetic energy at minus the half step for determining
1018 * the full step kinetic energy and possibly for T-coupling.*/
1019 /* This may not be quite working correctly yet . . . . */
1020 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1021 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
1022 constr,NULL,FALSE,state->box,
1023 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1024 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1026 clear_mat(force_vir);
1028 /* Ionize the atoms if necessary */
1031 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
1032 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
1035 /* Update force field in ffscan program */
1038 if (update_forcefield(fplog,
1040 mdatoms->nr,state->x,state->box))
1048 GMX_MPE_LOG(ev_timestep2);
1050 /* We write a checkpoint at this MD step when:
1051 * either at an NS step when we signalled through gs,
1052 * or at the last step (but not when we do not want confout),
1053 * but never at the first step or with rerun.
1055 bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
1056 (bLastStep && (Flags & MD_CONFOUT))) &&
1057 step > ir->init_step && !bRerunMD);
1060 gs.set[eglsCHKPT] = 0;
1063 /* Determine the energy and pressure:
1064 * at nstcalcenergy steps and at energy output steps (set below).
1066 if (EI_VV(ir->eI) && (!bInitStep))
1068 /* for vv, the first half actually corresponds to the last step */
1069 bCalcEner = do_per_step(step-1,ir->nstcalcenergy);
1073 bCalcEner = do_per_step(step,ir->nstcalcenergy);
1075 bCalcVir = bCalcEner ||
1076 (ir->epc != epcNO && do_per_step(step,ir->nstpcouple));
1078 /* Do we need global communication ? */
1079 bGStat = (bCalcVir || bCalcEner || bStopCM ||
1080 do_per_step(step,nstglobalcomm) ||
1081 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1083 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1085 if (do_ene || do_log)
1092 /* these CGLO_ options remain the same throughout the iteration */
1093 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1094 (bGStat ? CGLO_GSTAT : 0)
1097 force_flags = (GMX_FORCE_STATECHANGED |
1098 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1099 GMX_FORCE_ALLFORCES |
1101 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
1102 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
1103 (bDoFEP ? GMX_FORCE_DHDL : 0)
1108 if(do_per_step(step,ir->nstcalclr))
1110 force_flags |= GMX_FORCE_DO_LR;
1116 /* Now is the time to relax the shells */
1117 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
1119 bStopCM,top,top_global,
1121 state,f,force_vir,mdatoms,
1122 nrnb,wcycle,graph,groups,
1123 shellfc,fr,bBornRadii,t,mu_tot,
1124 state->natoms,&bConverged,vsite,
1135 /* The coordinates (x) are shifted (to get whole molecules)
1137 * This is parallellized as well, and does communication too.
1138 * Check comments in sim_util.c
1140 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1141 state->box,state->x,&state->hist,
1142 f,force_vir,mdatoms,enerd,fcd,
1143 state->lambda,graph,
1144 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1145 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1148 GMX_BARRIER(cr->mpi_comm_mygroup);
1152 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1153 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1156 if (bTCR && bFirstStep)
1158 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1159 fprintf(fplog,"Done init_coupling\n");
1163 if (bVV && !bStartingFromCpt && !bRerunMD)
1164 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1166 if (ir->eI==eiVV && bInitStep)
1168 /* if using velocity verlet with full time step Ekin,
1169 * take the first half step only to compute the
1170 * virial for the first step. From there,
1171 * revert back to the initial coordinates
1172 * so that the input is actually the initial step.
1174 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1176 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1177 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1180 /* If we are using twin-range interactions where the long-range component
1181 * is only evaluated every nstcalclr>1 steps, we should do a special update
1182 * step to combine the long-range forces on these steps.
1183 * For nstcalclr=1 this is not done, since the forces would have been added
1184 * directly to the short-range forces already.
1186 bUpdateDoLR = (fr->bTwinRange && do_per_step(step,ir->nstcalclr));
1188 update_coords(fplog,step,ir,mdatoms,state,fr->bMolPBC,
1189 f,bUpdateDoLR,fr->f_twin,fcd,
1190 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1191 cr,nrnb,constr,&top->idef);
1195 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1197 /* for iterations, we save these vectors, as we will be self-consistently iterating
1200 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1202 /* save the state */
1203 if (bIterations && iterate.bIterate) {
1204 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1207 bFirstIterate = TRUE;
1208 while (bFirstIterate || (bIterations && iterate.bIterate))
1210 if (bIterations && iterate.bIterate)
1212 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1213 if (bFirstIterate && bTrotter)
1215 /* The first time through, we need a decent first estimate
1216 of veta(t+dt) to compute the constraints. Do
1217 this by computing the box volume part of the
1218 trotter integration at this time. Nothing else
1219 should be changed by this routine here. If
1220 !(first time), we start with the previous value
1223 veta_save = state->veta;
1224 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1225 vetanew = state->veta;
1226 state->veta = veta_save;
1231 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
1234 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,
1235 state,fr->bMolPBC,graph,f,
1236 &top->idef,shake_vir,NULL,
1237 cr,nrnb,wcycle,upd,constr,
1238 bInitStep,TRUE,bCalcVir,vetanew);
1240 if (!bOK && !bFFscan)
1242 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1247 { /* Need to unshift here if a do_force has been
1248 called in the previous step */
1249 unshift_self(graph,state->box,state->x);
1253 /* if VV, compute the pressure and constraints */
1254 /* For VV2, we strictly only need this if using pressure
1255 * control, but we really would like to have accurate pressures
1257 * Think about ways around this in the future?
1258 * For now, keep this choice in comments.
1260 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
1261 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
1263 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
1264 if (bCalcEner && ir->eI==eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/
1266 bSumEkinhOld = TRUE;
1268 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1269 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1270 constr,NULL,FALSE,state->box,
1271 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1274 | (bStopCM ? CGLO_STOPCM : 0)
1275 | (bTemp ? CGLO_TEMPERATURE:0)
1276 | (bPres ? CGLO_PRESSURE : 0)
1277 | (bPres ? CGLO_CONSTRAINT : 0)
1278 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
1279 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1282 /* explanation of above:
1283 a) We compute Ekin at the full time step
1284 if 1) we are using the AveVel Ekin, and it's not the
1285 initial step, or 2) if we are using AveEkin, but need the full
1286 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1287 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
1288 EkinAveVel because it's needed for the pressure */
1290 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1295 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
1302 /* We need the kinetic energy at minus the half step for determining
1303 * the full step kinetic energy and possibly for T-coupling.*/
1304 /* This may not be quite working correctly yet . . . . */
1305 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1306 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
1307 constr,NULL,FALSE,state->box,
1308 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1309 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1313 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1318 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1319 state->veta,&vetanew))
1323 bFirstIterate = FALSE;
1326 if (bTrotter && !bInitStep) {
1327 enerd->term[F_DVDL_BONDED] += dvdl; /* only add after iterations */
1328 copy_mat(shake_vir,state->svir_prev);
1329 copy_mat(force_vir,state->fvir_prev);
1330 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
1331 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1332 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
1333 enerd->term[F_EKIN] = trace(ekind->ekin);
1336 /* if it's the initial step, we performed this first step just to get the constraint virial */
1337 if (bInitStep && ir->eI==eiVV) {
1338 copy_rvecn(cbuf,state->v,0,state->natoms);
1341 if (fr->bSepDVDL && fplog && do_log)
1343 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
1345 enerd->term[F_DVDL_BONDED] += dvdl;
1347 GMX_MPE_LOG(ev_timestep1);
1350 /* MRS -- now done iterating -- compute the conserved quantity */
1352 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
1355 last_ekin = enerd->term[F_EKIN];
1357 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1359 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1361 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1362 sum_dhdl(enerd,state->lambda,ir->fepvals);
1365 /* ######## END FIRST UPDATE STEP ############## */
1366 /* ######## If doing VV, we now have v(dt) ###### */
1368 /* perform extended ensemble sampling in lambda - we don't
1369 actually move to the new state before outputting
1370 statistics, but if performing simulated tempering, we
1371 do update the velocities and the tau_t. */
1373 lamnew = ExpandedEnsembleDynamics(fplog,ir,enerd,state,&MassQ,&df_history,step,mcrng,state->v,mdatoms);
1375 /* ################## START TRAJECTORY OUTPUT ################# */
1377 /* Now we have the energies and forces corresponding to the
1378 * coordinates at time t. We must output all of this before
1380 * for RerunMD t is read from input trajectory
1382 GMX_MPE_LOG(ev_output_start);
1385 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
1386 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
1387 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
1388 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
1389 if (bCPT) { mdof_flags |= MDOF_CPT; };
1391 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
1394 /* Enforce writing positions and velocities at end of run */
1395 mdof_flags |= (MDOF_X | MDOF_V);
1400 fcReportProgress( ir->nsteps, step );
1402 /* sync bCPT and fc record-keeping */
1403 if (bCPT && MASTER(cr))
1404 fcRequestCheckPoint();
1407 if (mdof_flags != 0)
1409 wallcycle_start(wcycle,ewcTRAJ);
1412 if (state->flags & (1<<estLD_RNG))
1414 get_stochd_state(upd,state);
1416 if (state->flags & (1<<estMC_RNG))
1418 get_mc_state(mcrng,state);
1424 state_global->ekinstate.bUpToDate = FALSE;
1428 update_ekinstate(&state_global->ekinstate,ekind);
1429 state_global->ekinstate.bUpToDate = TRUE;
1431 update_energyhistory(&state_global->enerhist,mdebin);
1432 if (ir->efep!=efepNO || ir->bSimTemp)
1434 state_global->fep_state = state->fep_state; /* MRS: seems kludgy. The code should be
1435 structured so this isn't necessary.
1436 Note this reassignment is only necessary
1437 for single threads.*/
1438 copy_df_history(&state_global->dfhist,&df_history);
1442 write_traj(fplog,cr,outf,mdof_flags,top_global,
1443 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
1450 if (bLastStep && step_rel == ir->nsteps &&
1451 (Flags & MD_CONFOUT) && MASTER(cr) &&
1452 !bRerunMD && !bFFscan)
1454 /* x and v have been collected in write_traj,
1455 * because a checkpoint file will always be written
1458 fprintf(stderr,"\nWriting final coordinates.\n");
1461 /* Make molecules whole only for confout writing */
1462 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
1464 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
1465 *top_global->name,top_global,
1466 state_global->x,state_global->v,
1467 ir->ePBC,state->box);
1470 wallcycle_stop(wcycle,ewcTRAJ);
1472 GMX_MPE_LOG(ev_output_finish);
1474 /* kludge -- virial is lost with restart for NPT control. Must restart */
1475 if (bStartingFromCpt && bVV)
1477 copy_mat(state->svir_prev,shake_vir);
1478 copy_mat(state->fvir_prev,force_vir);
1480 /* ################## END TRAJECTORY OUTPUT ################ */
1482 /* Determine the wallclock run time up till now */
1483 run_time = gmx_gettime() - (double)runtime->real;
1485 /* Check whether everything is still allright */
1486 if (((int)gmx_get_stop_condition() > handled_stop_condition)
1487 #ifdef GMX_THREAD_MPI
1492 /* this is just make gs.sig compatible with the hack
1493 of sending signals around by MPI_Reduce with together with
1495 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
1496 gs.sig[eglsSTOPCOND]=1;
1497 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
1498 gs.sig[eglsSTOPCOND]=-1;
1499 /* < 0 means stop at next step, > 0 means stop at next NS step */
1503 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1504 gmx_get_signal_name(),
1505 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1509 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
1510 gmx_get_signal_name(),
1511 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
1513 handled_stop_condition=(int)gmx_get_stop_condition();
1515 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1516 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
1517 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
1519 /* Signal to terminate the run */
1520 gs.sig[eglsSTOPCOND] = 1;
1523 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1525 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
1528 if (bResetCountersHalfMaxH && MASTER(cr) &&
1529 run_time > max_hours*60.0*60.0*0.495)
1531 gs.sig[eglsRESETCOUNTERS] = 1;
1534 if (ir->nstlist == -1 && !bRerunMD)
1536 /* When bGStatEveryStep=FALSE, global_stat is only called
1537 * when we check the atom displacements, not at NS steps.
1538 * This means that also the bonded interaction count check is not
1539 * performed immediately after NS. Therefore a few MD steps could
1540 * be performed with missing interactions.
1541 * But wrong energies are never written to file,
1542 * since energies are only written after global_stat
1545 if (step >= nlh.step_nscheck)
1547 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
1548 nlh.scale_tot,state->x);
1552 /* This is not necessarily true,
1553 * but step_nscheck is determined quite conservatively.
1559 /* In parallel we only have to check for checkpointing in steps
1560 * where we do global communication,
1561 * otherwise the other nodes don't know.
1563 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1566 run_time >= nchkpt*cpt_period*60.0)) &&
1567 gs.set[eglsCHKPT] == 0)
1569 gs.sig[eglsCHKPT] = 1;
1573 /* at the start of step, randomize the velocities */
1574 if (ETC_ANDERSEN(ir->etc) && EI_VV(ir->eI))
1576 gmx_bool bDoAndersenConstr;
1577 bDoAndersenConstr = (constr && update_randomize_velocities(ir,step,mdatoms,state,upd,&top->idef,constr));
1578 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1579 if (bDoAndersenConstr)
1581 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,
1582 state,fr->bMolPBC,graph,f,
1583 &top->idef,tmp_vir,NULL,
1584 cr,nrnb,wcycle,upd,constr,
1585 bInitStep,TRUE,bCalcVir,vetanew);
1591 gmx_iterate_init(&iterate,bIterations);
1594 /* for iterations, we save these vectors, as we will be redoing the calculations */
1595 if (bIterations && iterate.bIterate)
1597 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1599 bFirstIterate = TRUE;
1600 while (bFirstIterate || (bIterations && iterate.bIterate))
1602 /* We now restore these vectors to redo the calculation with improved extended variables */
1605 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1608 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
1609 so scroll down for that logic */
1611 /* ######### START SECOND UPDATE STEP ################# */
1612 GMX_MPE_LOG(ev_update_start);
1613 /* Box is changed in update() when we do pressure coupling,
1614 * but we should still use the old box for energy corrections and when
1615 * writing it to the energy file, so it matches the trajectory files for
1616 * the same timestep above. Make a copy in a separate array.
1618 copy_mat(state->box,lastbox);
1621 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
1623 wallcycle_start(wcycle,ewcUPDATE);
1625 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1628 if (bIterations && iterate.bIterate)
1636 /* we use a new value of scalevir to converge the iterations faster */
1637 scalevir = tracevir/trace(shake_vir);
1639 msmul(shake_vir,scalevir,shake_vir);
1640 m_add(force_vir,shake_vir,total_vir);
1641 clear_mat(shake_vir);
1643 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
1644 /* We can only do Berendsen coupling after we have summed
1645 * the kinetic energy or virial. Since the happens
1646 * in global_state after update, we should only do it at
1647 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1652 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
1653 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
1659 bUpdateDoLR = (fr->bTwinRange && do_per_step(step,ir->nstcalclr));
1661 /* velocity half-step update */
1662 update_coords(fplog,step,ir,mdatoms,state,fr->bMolPBC,f,
1663 bUpdateDoLR,fr->f_twin,fcd,
1664 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
1665 cr,nrnb,constr,&top->idef);
1668 /* Above, initialize just copies ekinh into ekin,
1669 * it doesn't copy position (for VV),
1670 * and entire integrator for MD.
1675 copy_rvecn(state->x,cbuf,0,state->natoms);
1677 bUpdateDoLR = (fr->bTwinRange && do_per_step(step,ir->nstcalclr));
1679 update_coords(fplog,step,ir,mdatoms,state,fr->bMolPBC,f,
1680 bUpdateDoLR,fr->f_twin,fcd,
1681 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1682 wallcycle_stop(wcycle,ewcUPDATE);
1684 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,
1685 fr->bMolPBC,graph,f,
1686 &top->idef,shake_vir,force_vir,
1687 cr,nrnb,wcycle,upd,constr,
1688 bInitStep,FALSE,bCalcVir,state->veta);
1692 /* erase F_EKIN and F_TEMP here? */
1693 /* just compute the kinetic energy at the half step to perform a trotter step */
1694 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1695 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1696 constr,NULL,FALSE,lastbox,
1697 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1698 cglo_flags | CGLO_TEMPERATURE
1700 wallcycle_start(wcycle,ewcUPDATE);
1701 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
1702 /* now we know the scaling, we can compute the positions again again */
1703 copy_rvecn(cbuf,state->x,0,state->natoms);
1705 bUpdateDoLR = (fr->bTwinRange && do_per_step(step,ir->nstcalclr));
1707 update_coords(fplog,step,ir,mdatoms,state,fr->bMolPBC,f,
1708 bUpdateDoLR,fr->f_twin,fcd,
1709 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
1710 wallcycle_stop(wcycle,ewcUPDATE);
1712 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
1713 /* are the small terms in the shake_vir here due
1714 * to numerical errors, or are they important
1715 * physically? I'm thinking they are just errors, but not completely sure.
1716 * For now, will call without actually constraining, constr=NULL*/
1717 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,
1718 state,fr->bMolPBC,graph,f,
1719 &top->idef,tmp_vir,force_vir,
1720 cr,nrnb,wcycle,upd,NULL,
1721 bInitStep,FALSE,bCalcVir,
1724 if (!bOK && !bFFscan)
1726 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
1729 if (fr->bSepDVDL && fplog && do_log)
1731 fprintf(fplog,sepdvdlformat,"Constraint dV/dl",0.0,dvdl);
1733 enerd->term[F_DVDL_BONDED] += dvdl;
1737 /* Need to unshift here */
1738 unshift_self(graph,state->box,state->x);
1741 GMX_BARRIER(cr->mpi_comm_mygroup);
1742 GMX_MPE_LOG(ev_update_finish);
1746 wallcycle_start(wcycle,ewcVSITECONSTR);
1749 shift_self(graph,state->box,state->x);
1751 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1752 top->idef.iparams,top->idef.il,
1753 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1757 unshift_self(graph,state->box,state->x);
1759 wallcycle_stop(wcycle,ewcVSITECONSTR);
1762 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
1763 /* With Leap-Frog we can skip compute_globals at
1764 * non-communication steps, but we need to calculate
1765 * the kinetic energy one step before communication.
1767 if (bGStat || do_per_step(step+1,nstglobalcomm) ||
1770 if (ir->nstlist == -1 && bFirstIterate)
1772 gs.sig[eglsNABNSB] = nlh.nabnsb;
1774 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1775 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1777 bFirstIterate ? &gs : NULL,
1778 (step_rel % gs.nstms == 0) &&
1779 (multisim_nsteps<0 || (step_rel<multisim_nsteps)),
1781 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1783 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
1784 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1785 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1786 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1787 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
1788 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
1791 if (ir->nstlist == -1 && bFirstIterate)
1793 nlh.nabnsb = gs.set[eglsNABNSB];
1794 gs.set[eglsNABNSB] = 0;
1797 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
1798 /* ############# END CALC EKIN AND PRESSURE ################# */
1800 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1801 the virial that should probably be addressed eventually. state->veta has better properies,
1802 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1803 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1806 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
1807 trace(shake_vir),&tracevir))
1811 bFirstIterate = FALSE;
1814 /* only add constraint dvdl after constraints */
1815 enerd->term[F_DVDL_BONDED] += dvdl;
1818 /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */
1819 sum_dhdl(enerd,state->lambda,ir->fepvals);
1821 update_box(fplog,step,ir,mdatoms,state,graph,f,
1822 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
1824 /* ################# END UPDATE STEP 2 ################# */
1825 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1827 /* The coordinates (x) were unshifted in update */
1828 if (bFFscan && (shellfc==NULL || bConverged))
1830 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
1832 &(top_global->mols),mdatoms->massT,pres))
1836 fprintf(stderr,"\n");
1842 /* We will not sum ekinh_old,
1843 * so signal that we still have to do it.
1845 bSumEkinhOld = TRUE;
1850 /* Only do GCT when the relaxation of shells (minimization) has converged,
1851 * otherwise we might be coupling to bogus energies.
1852 * In parallel we must always do this, because the other sims might
1856 /* Since this is called with the new coordinates state->x, I assume
1857 * we want the new box state->box too. / EL 20040121
1859 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
1861 mdatoms,&(top->idef),mu_aver,
1862 top_global->mols.nr,cr,
1863 state->box,total_vir,pres,
1864 mu_tot,state->x,f,bConverged);
1868 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1870 /* use the directly determined last velocity, not actually the averaged half steps */
1871 if (bTrotter && ir->eI==eiVV)
1873 enerd->term[F_EKIN] = last_ekin;
1875 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1879 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1883 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
1885 /* Check for excessively large energies */
1889 real etot_max = 1e200;
1891 real etot_max = 1e30;
1893 if (fabs(enerd->term[F_ETOT]) > etot_max)
1895 fprintf(stderr,"Energy too large (%g), giving up\n",
1896 enerd->term[F_ETOT]);
1899 /* ######### END PREPARING EDR OUTPUT ########### */
1901 /* Time for performance */
1902 if (((step % stepout) == 0) || bLastStep)
1904 runtime_upd_proc(runtime);
1910 gmx_bool do_dr,do_or;
1912 if (fplog && do_log && bDoExpanded)
1914 /* only needed if doing expanded ensemble */
1915 PrintFreeEnergyInfoToFile(fplog,ir->fepvals,ir->expandedvals,ir->bSimTemp?ir->simtempvals:NULL,
1916 &df_history,state->fep_state,ir->nstlog,step);
1918 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
1922 upd_mdebin(mdebin,bDoDHDL, TRUE,
1923 t,mdatoms->tmass,enerd,state,
1924 ir->fepvals,ir->expandedvals,lastbox,
1925 shake_vir,force_vir,total_vir,pres,
1926 ekind,mu_tot,constr);
1930 upd_mdebin_step(mdebin);
1933 do_dr = do_per_step(step,ir->nstdisreout);
1934 do_or = do_per_step(step,ir->nstorireout);
1936 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
1938 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
1940 if (ir->ePull != epullNO)
1942 pull_print_output(ir->pull,step,t);
1945 if (do_per_step(step,ir->nstlog))
1947 if(fflush(fplog) != 0)
1949 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
1955 /* Have to do this part after outputting the logfile and the edr file */
1956 state->fep_state = lamnew;
1957 for (i=0;i<efptNR;i++)
1959 state->lambda[i] = ir->fepvals->all_lambda[i][lamnew];
1962 /* Remaining runtime */
1963 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal()) && !bPMETuneRunning)
1967 fprintf(stderr,"\n");
1969 print_time(stderr,runtime,step,ir,cr);
1972 /* Replica exchange */
1974 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
1975 do_per_step(step,repl_ex_nst))
1977 bExchanged = replica_exchange(fplog,cr,repl_ex,
1981 if (bExchanged && DOMAINDECOMP(cr))
1983 dd_partition_system(fplog,step,cr,TRUE,1,
1984 state_global,top_global,ir,
1985 state,&f,mdatoms,top,fr,
1986 vsite,shellfc,constr,
1993 bStartingFromCpt = FALSE;
1995 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1996 /* With all integrators, except VV, we need to retain the pressure
1997 * at the current step for coupling at the next step.
1999 if ((state->flags & (1<<estPRES_PREV)) &&
2001 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
2003 /* Store the pressure in t_state for pressure coupling
2004 * at the next MD step.
2006 copy_mat(pres,state->pres_prev);
2009 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2011 if ( (membed!=NULL) && (!bLastStep) )
2013 rescale_membed(step_rel,membed,state_global->x);
2020 /* read next frame from input trajectory */
2021 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
2026 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
2030 if (!bRerunMD || !rerun_fr.bStep)
2032 /* increase the MD step number */
2037 cycles = wallcycle_stop(wcycle,ewcSTEP);
2038 if (DOMAINDECOMP(cr) && wcycle)
2040 dd_cycles_add(cr->dd,cycles,ddCyclStep);
2043 if (bPMETuneRunning || bPMETuneTry)
2045 /* PME grid + cut-off optimization with GPUs or PME nodes */
2047 /* Count the total cycles over the last steps */
2048 cycles_pmes += cycles;
2050 /* We can only switch cut-off at NS steps */
2051 if (step % ir->nstlist == 0)
2053 /* PME grid + cut-off optimization with GPUs or PME nodes */
2056 if (DDMASTER(cr->dd))
2058 /* PME node load is too high, start tuning */
2059 bPMETuneRunning = (dd_pme_f_ratio(cr->dd) >= 1.05);
2061 dd_bcast(cr->dd,sizeof(gmx_bool),&bPMETuneRunning);
2063 if (bPMETuneRunning || step_rel > ir->nstlist*50)
2065 bPMETuneTry = FALSE;
2068 if (bPMETuneRunning)
2070 /* init_step might not be a multiple of nstlist,
2071 * but the first cycle is always skipped anyhow.
2074 pme_load_balance(pme_loadbal,cr,
2075 (bVerbose && MASTER(cr)) ? stderr : NULL,
2077 ir,state,cycles_pmes,
2078 fr->ic,fr->nbv,&fr->pmedata,
2081 /* Update constants in forcerec/inputrec to keep them in sync with fr->ic */
2082 fr->ewaldcoeff = fr->ic->ewaldcoeff;
2083 fr->rlist = fr->ic->rlist;
2084 fr->rlistlong = fr->ic->rlistlong;
2085 fr->rcoulomb = fr->ic->rcoulomb;
2086 fr->rvdw = fr->ic->rvdw;
2092 if (step_rel == wcycle_get_reset_counters(wcycle) ||
2093 gs.set[eglsRESETCOUNTERS] != 0)
2095 /* Reset all the counters related to performance over the run */
2096 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime,
2097 fr->nbv != NULL && fr->nbv->bUseGPU ? fr->nbv->cu_nbv : NULL);
2098 wcycle_set_reset_counters(wcycle,-1);
2099 /* Correct max_hours for the elapsed time */
2100 max_hours -= run_time/(60.0*60.0);
2101 bResetCountersHalfMaxH = FALSE;
2102 gs.set[eglsRESETCOUNTERS] = 0;
2106 /* End of main MD loop */
2110 runtime_end(runtime);
2112 if (bRerunMD && MASTER(cr))
2117 if (!(cr->duty & DUTY_PME))
2119 /* Tell the PME only node to finish */
2120 gmx_pme_send_finish(cr);
2125 if (ir->nstcalcenergy > 0 && !bRerunMD)
2127 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
2128 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
2136 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
2138 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)));
2139 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
2142 if (pme_loadbal != NULL)
2144 pme_loadbal_done(pme_loadbal,fplog);
2147 if (shellfc && fplog)
2149 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
2150 (nconverged*100.0)/step_rel);
2151 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
2155 if (repl_ex_nst > 0 && MASTER(cr))
2157 print_replica_exchange_statistics(fplog,repl_ex);
2160 runtime->nsteps_done = step_rel;