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
43 #if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
80 #include "mpelogging.h"
87 #include "compute_io.h"
89 #include "checkpoint.h"
90 #include "mtop_util.h"
91 #include "sighandler.h"
101 #include "corewrap.h"
105 /* simulation conditions to transmit. Keep in mind that they are
106 transmitted to other nodes through an MPI_Reduce after
107 casting them to a real (so the signals can be sent together with other
108 data). This means that the only meaningful values are positive,
110 enum { eglsNABNSB, eglsCHKPT, eglsSTOPCOND, eglsRESETCOUNTERS, eglsNR };
111 /* Is the signal in one simulation independent of other simulations? */
112 bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
115 int nstms; /* The frequency for intersimulation communication */
116 int sig[eglsNR]; /* The signal set by one process in do_md */
117 int set[eglsNR]; /* The communicated signal, equal for all processes */
121 static int multisim_min(const gmx_multisim_t *ms,int nmin,int n)
129 gmx_sumi_sim(ms->nsim,buf,ms);
132 for(s=0; s<ms->nsim; s++)
134 bPos = bPos && (buf[s] > 0);
135 bEqual = bEqual && (buf[s] == buf[0]);
141 nmin = min(nmin,buf[0]);
145 /* Find the least common multiple */
146 for(d=2; d<nmin; d++)
149 while (s < ms->nsim && d % buf[s] == 0)
155 /* We found the LCM and it is less than nmin */
167 static int multisim_nstsimsync(const t_commrec *cr,
168 const t_inputrec *ir,int repl_ex_nst)
175 nmin = multisim_min(cr->ms,nmin,ir->nstlist);
176 nmin = multisim_min(cr->ms,nmin,ir->nstcalcenergy);
177 nmin = multisim_min(cr->ms,nmin,repl_ex_nst);
180 gmx_fatal(FARGS,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
182 /* Avoid inter-simulation communication at every (second) step */
189 gmx_bcast(sizeof(int),&nmin,cr);
194 static void init_global_signals(globsig_t *gs,const t_commrec *cr,
195 const t_inputrec *ir,int repl_ex_nst)
201 gs->nstms = multisim_nstsimsync(cr,ir,repl_ex_nst);
204 fprintf(debug,"Syncing simulations for checkpointing and termination every %d steps\n",gs->nstms);
212 for(i=0; i<eglsNR; i++)
219 static void copy_coupling_state(t_state *statea,t_state *stateb,
220 gmx_ekindata_t *ekinda,gmx_ekindata_t *ekindb, t_grpopts* opts)
223 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
227 /* Make sure we have enough space for x and v */
228 if (statea->nalloc > stateb->nalloc)
230 stateb->nalloc = statea->nalloc;
231 srenew(stateb->x,stateb->nalloc);
232 srenew(stateb->v,stateb->nalloc);
235 stateb->natoms = statea->natoms;
236 stateb->ngtc = statea->ngtc;
237 stateb->nnhpres = statea->nnhpres;
238 stateb->veta = statea->veta;
241 copy_mat(ekinda->ekin,ekindb->ekin);
242 for (i=0; i<stateb->ngtc; i++)
244 ekindb->tcstat[i].T = ekinda->tcstat[i].T;
245 ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
246 copy_mat(ekinda->tcstat[i].ekinh,ekindb->tcstat[i].ekinh);
247 copy_mat(ekinda->tcstat[i].ekinf,ekindb->tcstat[i].ekinf);
248 ekindb->tcstat[i].ekinscalef_nhc = ekinda->tcstat[i].ekinscalef_nhc;
249 ekindb->tcstat[i].ekinscaleh_nhc = ekinda->tcstat[i].ekinscaleh_nhc;
250 ekindb->tcstat[i].vscale_nhc = ekinda->tcstat[i].vscale_nhc;
253 copy_rvecn(statea->x,stateb->x,0,stateb->natoms);
254 copy_rvecn(statea->v,stateb->v,0,stateb->natoms);
255 copy_mat(statea->box,stateb->box);
256 copy_mat(statea->box_rel,stateb->box_rel);
257 copy_mat(statea->boxv,stateb->boxv);
259 for (i = 0; i<stateb->ngtc; i++)
261 nc = i*opts->nhchainlength;
262 for (j=0; j<opts->nhchainlength; j++)
264 stateb->nosehoover_xi[nc+j] = statea->nosehoover_xi[nc+j];
265 stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
268 if (stateb->nhpres_xi != NULL)
270 for (i = 0; i<stateb->nnhpres; i++)
272 nc = i*opts->nhchainlength;
273 for (j=0; j<opts->nhchainlength; j++)
275 stateb->nhpres_xi[nc+j] = statea->nhpres_xi[nc+j];
276 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
282 static real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
292 quantity = NPT_energy(ir,state,MassQ);
295 quantity = vrescale_energy(&(ir->opts),state->therm_integral);
303 static void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
304 t_forcerec *fr, gmx_ekindata_t *ekind,
305 t_state *state, t_state *state_global, t_mdatoms *mdatoms,
306 t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
307 gmx_enerdata_t *enerd,tensor force_vir, tensor shake_vir, tensor total_vir,
308 tensor pres, rvec mu_tot, gmx_constr_t constr,
309 globsig_t *gs,bool bInterSimGS,
310 matrix box, gmx_mtop_t *top_global, real *pcurr,
311 int natoms, bool *bSumEkinhOld, int flags)
315 tensor corr_vir,corr_pres,shakeall_vir;
316 bool bEner,bPres,bTemp, bVV;
317 bool bRerunMD, bStopCM, bGStat, bIterate,
318 bFirstIterate,bReadEkin,bEkinAveVel,bScaleEkin, bConstrain;
319 real ekin,temp,prescorr,enercorr,dvdlcorr;
321 /* translate CGLO flags to booleans */
322 bRerunMD = flags & CGLO_RERUNMD;
323 bStopCM = flags & CGLO_STOPCM;
324 bGStat = flags & CGLO_GSTAT;
325 bReadEkin = flags & CGLO_READEKIN;
326 bScaleEkin = flags & CGLO_SCALEEKIN;
327 bEner = flags & CGLO_ENERGY;
328 bTemp = flags & CGLO_TEMPERATURE;
329 bPres = flags & CGLO_PRESSURE;
330 bConstrain = flags & CGLO_CONSTRAINT;
331 bIterate = flags & CGLO_ITERATE;
332 bFirstIterate = flags & CGLO_FIRSTITERATE;
334 /* we calculate a full state kinetic energy either with full-step velocity verlet
335 or half step where we need the pressure */
337 bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && bPres) || bReadEkin);
339 /* in initalization, it sums the shake virial in vv, and to
340 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
342 /* ########## Kinetic energy ############## */
346 /* Non-equilibrium MD: this is parallellized, but only does communication
347 * when there really is NEMD.
350 if (PAR(cr) && (ekind->bNEMD))
352 accumulate_u(cr,&(ir->opts),ekind);
357 restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
362 calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
367 /* Calculate center of mass velocity if necessary, also parallellized */
368 if (bStopCM && !bRerunMD && bEner)
370 calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
371 state->x,state->v,vcm);
375 if (bTemp || bPres || bEner || bConstrain)
379 /* We will not sum ekinh_old,
380 * so signal that we still have to do it.
382 *bSumEkinhOld = TRUE;
389 for(i=0; i<eglsNR; i++)
391 gs_buf[i] = gs->sig[i];
396 wallcycle_start(wcycle,ewcMoveE);
397 GMX_MPE_LOG(ev_global_stat_start);
398 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
400 gs != NULL ? eglsNR : 0,gs_buf,
402 *bSumEkinhOld,flags);
403 GMX_MPE_LOG(ev_global_stat_finish);
404 wallcycle_stop(wcycle,ewcMoveE);
408 if (MULTISIM(cr) && bInterSimGS)
412 /* Communicate the signals between the simulations */
413 gmx_sum_sim(eglsNR,gs_buf,cr->ms);
415 /* Communicate the signals form the master to the others */
416 gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
418 for(i=0; i<eglsNR; i++)
420 if (bInterSimGS || gs_simlocal[i])
422 /* Set the communicated signal only when it is non-zero,
423 * since signals might not be processed at each MD step.
425 gsi = (gs_buf[i] >= 0 ?
426 (int)(gs_buf[i] + 0.5) :
427 (int)(gs_buf[i] - 0.5));
432 /* Turn off the local signal */
437 *bSumEkinhOld = FALSE;
441 if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
444 mdatoms->start,mdatoms->start+mdatoms->homenr,
445 state->v,vcm->group_p[0],
446 mdatoms->massT,mdatoms->tmass,ekind->ekin);
450 /* Do center of mass motion removal */
451 if (bStopCM && !bRerunMD) /* is this correct? Does it get called too often with this logic? */
453 check_cm_grp(fplog,vcm,ir,1);
454 do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
455 state->x,state->v,vcm);
456 inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
462 /* Sum the kinetic energies of the groups & calc temp */
463 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
464 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
465 Leap with AveVel is not supported; it's not clear that it will actually work.
466 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
467 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
468 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
469 If FALSE, we go ahead and erase over it.
471 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,&(enerd->term[F_DKDL]),
472 bEkinAveVel,bIterate,bScaleEkin);
474 enerd->term[F_EKIN] = trace(ekind->ekin);
477 /* ########## Long range energy information ###### */
479 if (bEner || bPres || bConstrain)
481 calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda,
482 corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
485 if (bEner && bFirstIterate)
487 enerd->term[F_DISPCORR] = enercorr;
488 enerd->term[F_EPOT] += enercorr;
489 enerd->term[F_DVDL] += dvdlcorr;
490 if (fr->efep != efepNO) {
491 enerd->dvdl_lin += dvdlcorr;
495 /* ########## Now pressure ############## */
496 if (bPres || bConstrain)
499 m_add(force_vir,shake_vir,total_vir);
501 /* Calculate pressure and apply LR correction if PPPM is used.
502 * Use the box from last timestep since we already called update().
505 enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres,
506 (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
508 /* Calculate long range corrections to pressure and energy */
509 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
510 and computes enerd->term[F_DISPCORR]. Also modifies the
511 total_vir and pres tesors */
513 m_add(total_vir,corr_vir,total_vir);
514 m_add(pres,corr_pres,pres);
515 enerd->term[F_PDISPCORR] = prescorr;
516 enerd->term[F_PRES] += prescorr;
517 *pcurr = enerd->term[F_PRES];
518 /* calculate temperature using virial */
519 enerd->term[F_VTEMP] = calc_temp(trace(total_vir),ir->opts.nrdf[0]);
525 /* Definitions for convergence of iterated constraints */
527 /* iterate constraints up to 50 times */
528 #define MAXITERCONST 50
533 real f,fprev,x,xprev;
536 real allrelerr[MAXITERCONST+2];
537 int num_close; /* number of "close" violations, caused by limited precision. */
541 #define CONVERGEITER 0.000000001
542 #define CLOSE_ENOUGH 0.000001000
544 #define CONVERGEITER 0.0001
545 #define CLOSE_ENOUGH 0.0050
548 /* we want to keep track of the close calls. If there are too many, there might be some other issues.
549 so we make sure that it's either less than some predetermined number, or if more than that number,
550 only some small fraction of the total. */
551 #define MAX_NUMBER_CLOSE 50
552 #define FRACTION_CLOSE 0.001
554 /* maximum length of cyclic traps to check, emerging from limited numerical precision */
557 static void gmx_iterate_init(gmx_iterate_t *iterate,bool bIterate)
562 iterate->bIterate = bIterate;
563 iterate->num_close = 0;
564 for (i=0;i<MAXITERCONST+2;i++)
566 iterate->allrelerr[i] = 0;
570 static bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate_t *iterate, bool bFirstIterate, real fom, real *newf)
572 /* monitor convergence, and use a secant search to propose new
575 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
576 f(x_{i}) - f(x_{i-1})
578 The function we are trying to zero is fom-x, where fom is the
579 "figure of merit" which is the pressure (or the veta value) we
580 would get by putting in an old value of the pressure or veta into
581 the incrementor function for the step or half step. I have
582 verified that this gives the same answer as self consistent
583 iteration, usually in many fewer steps, especially for small tau_p.
585 We could possibly eliminate an iteration with proper use
586 of the value from the previous step, but that would take a bit
587 more bookkeeping, especially for veta, since tests indicate the
588 function of veta on the last step is not sufficiently close to
589 guarantee convergence this step. This is
590 good enough for now. On my tests, I could use tau_p down to
591 0.02, which is smaller that would ever be necessary in
592 practice. Generally, 3-5 iterations will be sufficient */
594 real relerr,err,xmin;
602 iterate->f = fom-iterate->x;
609 iterate->f = fom-iterate->x; /* we want to zero this difference */
610 if ((iterate->iter_i > 1) && (iterate->iter_i < MAXITERCONST))
612 if (iterate->f==iterate->fprev)
618 *newf = iterate->x - (iterate->x-iterate->xprev)*(iterate->f)/(iterate->f-iterate->fprev);
623 /* just use self-consistent iteration the first step to initialize, or
624 if it's not converging (which happens occasionally -- need to investigate why) */
628 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
629 difference between the closest of x and xprev to the new
630 value. To be 100% certain, we should check the difference between
631 the last result, and the previous result, or
633 relerr = (fabs((x-xprev)/fom));
635 but this is pretty much never necessary under typical conditions.
636 Checking numerically, it seems to lead to almost exactly the same
637 trajectories, but there are small differences out a few decimal
638 places in the pressure, and eventually in the v_eta, but it could
641 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
642 relerr = (fabs((*newf-xmin) / *newf));
645 err = fabs((iterate->f-iterate->fprev));
646 relerr = fabs(err/fom);
648 iterate->allrelerr[iterate->iter_i] = relerr;
650 if (iterate->iter_i > 0)
654 fprintf(debug,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
655 iterate->iter_i,fom,relerr,*newf);
658 if ((relerr < CONVERGEITER) || (err < CONVERGEITER) || (fom==0) || ((iterate->x == iterate->xprev) && iterate->iter_i > 1))
660 iterate->bIterate = FALSE;
663 fprintf(debug,"Iterating NPT constraints: CONVERGED\n");
667 if (iterate->iter_i > MAXITERCONST)
669 if (relerr < CLOSE_ENOUGH)
672 for (i=1;i<CYCLEMAX;i++) {
673 if ((iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-1]) &&
674 (iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-(1+2*i)])) {
678 fprintf(debug,"Exiting from an NPT iterating cycle of length %d\n",i);
685 /* step 1: trapped in a numerical attractor */
686 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
687 Better to give up convergence here than have the simulation die.
689 iterate->num_close++;
694 /* Step #2: test if we are reasonably close for other reasons, then monitor the number. If not, die */
696 /* how many close calls have we had? If less than a few, we're OK */
697 if (iterate->num_close < MAX_NUMBER_CLOSE)
699 sprintf(buf,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps,relerr);
700 md_print_warning(cr,fplog,buf);
701 iterate->num_close++;
703 /* if more than a few, check the total fraction. If too high, die. */
704 } else if (iterate->num_close/(double)nsteps > FRACTION_CLOSE) {
705 gmx_fatal(FARGS,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate->num_close/(double)nsteps);
711 gmx_fatal(FARGS,"Could not converge NPT constraints\n");
716 iterate->xprev = iterate->x;
718 iterate->fprev = iterate->f;
724 static void check_nst_param(FILE *fplog,t_commrec *cr,
725 const char *desc_nst,int nst,
726 const char *desc_p,int *p)
730 if (*p > 0 && *p % nst != 0)
732 /* Round up to the next multiple of nst */
733 *p = ((*p)/nst + 1)*nst;
734 sprintf(buf,"NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
735 md_print_warning(cr,fplog,buf);
739 static void reset_all_counters(FILE *fplog,t_commrec *cr,
740 gmx_large_int_t step,
741 gmx_large_int_t *step_rel,t_inputrec *ir,
742 gmx_wallcycle_t wcycle,t_nrnb *nrnb,
743 gmx_runtime_t *runtime)
745 char buf[STRLEN],sbuf[STEPSTRSIZE];
747 /* Reset all the counters related to performance over the run */
748 sprintf(buf,"Step %s: resetting all time and cycle counters\n",
749 gmx_step_str(step,sbuf));
750 md_print_warning(cr,fplog,buf);
752 wallcycle_stop(wcycle,ewcRUN);
753 wallcycle_reset_all(wcycle);
754 if (DOMAINDECOMP(cr))
756 reset_dd_statistics_counters(cr->dd);
759 ir->init_step += *step_rel;
760 ir->nsteps -= *step_rel;
762 wallcycle_start(wcycle,ewcRUN);
763 runtime_start(runtime);
764 print_date_and_time(fplog,cr->nodeid,"Restarted time",runtime);
767 static void min_zero(int *n,int i)
769 if (i > 0 && (*n == 0 || i < *n))
775 static int lcd4(int i1,int i2,int i3,int i4)
786 gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
789 while (nst > 1 && ((i1 > 0 && i1 % nst != 0) ||
790 (i2 > 0 && i2 % nst != 0) ||
791 (i3 > 0 && i3 % nst != 0) ||
792 (i4 > 0 && i4 % nst != 0)))
800 static int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
801 int nstglobalcomm,t_inputrec *ir)
805 if (!EI_DYNAMICS(ir->eI))
810 if (nstglobalcomm == -1)
812 if (!(ir->nstcalcenergy > 0 ||
818 if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
820 nstglobalcomm = ir->nstenergy;
825 /* Ensure that we do timely global communication for
826 * (possibly) each of the four following options.
828 nstglobalcomm = lcd4(ir->nstcalcenergy,
830 ir->etc != etcNO ? ir->nsttcouple : 0,
831 ir->epc != epcNO ? ir->nstpcouple : 0);
836 if (ir->nstlist > 0 &&
837 nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
839 nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
840 sprintf(buf,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
841 md_print_warning(cr,fplog,buf);
843 if (ir->nstcalcenergy > 0)
845 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
846 "nstcalcenergy",&ir->nstcalcenergy);
848 if (ir->etc != etcNO && ir->nsttcouple > 0)
850 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
851 "nsttcouple",&ir->nsttcouple);
853 if (ir->epc != epcNO && ir->nstpcouple > 0)
855 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
856 "nstpcouple",&ir->nstpcouple);
859 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
860 "nstenergy",&ir->nstenergy);
862 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
863 "nstlog",&ir->nstlog);
866 if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
868 sprintf(buf,"WARNING: Changing nstcomm from %d to %d\n",
869 ir->nstcomm,nstglobalcomm);
870 md_print_warning(cr,fplog,buf);
871 ir->nstcomm = nstglobalcomm;
874 return nstglobalcomm;
877 void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
878 t_inputrec *ir,gmx_mtop_t *mtop)
880 /* Check required for old tpx files */
881 if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
882 ir->nstcalcenergy % ir->nstlist != 0)
884 md_print_warning(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
886 if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
887 gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
888 ir->eConstrAlg == econtSHAKE)
890 md_print_warning(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
891 if (ir->epc != epcNO)
893 gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
896 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
897 "nstcalcenergy",&ir->nstcalcenergy);
898 if (ir->epc != epcNO)
900 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
901 "nstpcouple",&ir->nstpcouple);
903 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
904 "nstenergy",&ir->nstenergy);
905 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
906 "nstlog",&ir->nstlog);
907 if (ir->efep != efepNO)
909 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
910 "nstdhdl",&ir->nstdhdl);
916 bool bGStatEveryStep;
917 gmx_large_int_t step_ns;
918 gmx_large_int_t step_nscheck;
929 static void reset_nlistheuristics(gmx_nlheur_t *nlh,gmx_large_int_t step)
933 nlh->step_nscheck = step;
936 static void init_nlistheuristics(gmx_nlheur_t *nlh,
937 bool bGStatEveryStep,gmx_large_int_t step)
939 nlh->bGStatEveryStep = bGStatEveryStep;
946 reset_nlistheuristics(nlh,step);
949 static void update_nliststatistics(gmx_nlheur_t *nlh,gmx_large_int_t step)
951 gmx_large_int_t nl_lt;
952 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
954 /* Determine the neighbor list life time */
955 nl_lt = step - nlh->step_ns;
958 fprintf(debug,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh->nabnsb,gmx_step_str(nl_lt,sbuf));
962 nlh->s2 += nl_lt*nl_lt;
963 nlh->ab += nlh->nabnsb;
964 if (nlh->lt_runav == 0)
966 nlh->lt_runav = nl_lt;
967 /* Initialize the fluctuation average
968 * such that at startup we check after 0 steps.
970 nlh->lt_runav2 = sqr(nl_lt/2.0);
972 /* Running average with 0.9 gives an exp. history of 9.5 */
973 nlh->lt_runav2 = 0.9*nlh->lt_runav2 + 0.1*sqr(nlh->lt_runav - nl_lt);
974 nlh->lt_runav = 0.9*nlh->lt_runav + 0.1*nl_lt;
975 if (nlh->bGStatEveryStep)
977 /* Always check the nlist validity */
978 nlh->step_nscheck = step;
982 /* We check after: <life time> - 2*sigma
983 * The factor 2 is quite conservative,
984 * but we assume that with nstlist=-1 the user
985 * prefers exact integration over performance.
987 nlh->step_nscheck = step
988 + (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)) - 1;
992 fprintf(debug,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
993 gmx_step_str(nl_lt,sbuf),nlh->lt_runav,sqrt(nlh->lt_runav2),
994 gmx_step_str(nlh->step_nscheck-step+1,sbuf2),
995 (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)));
999 static void set_nlistheuristics(gmx_nlheur_t *nlh,bool bReset,gmx_large_int_t step)
1005 reset_nlistheuristics(nlh,step);
1009 update_nliststatistics(nlh,step);
1012 nlh->step_ns = step;
1013 /* Initialize the cumulative coordinate scaling matrix */
1014 clear_mat(nlh->scale_tot);
1015 for(d=0; d<DIM; d++)
1017 nlh->scale_tot[d][d] = 1.0;
1021 static void rerun_parallel_comm(t_commrec *cr,t_trxframe *fr,
1022 bool *bNotLastFrame)
1027 bAlloc = (fr->natoms == 0);
1029 if (MASTER(cr) && !*bNotLastFrame)
1035 gmx_bcast(sizeof(*fr),fr,cr);
1039 *bNotLastFrame = (fr->natoms >= 0);
1041 if (*bNotLastFrame && PARTDECOMP(cr))
1043 /* x and v are the only variable size quantities stored in trr
1044 * that are required for rerun (f is not needed).
1048 snew(fr->x,fr->natoms);
1049 snew(fr->v,fr->natoms);
1053 gmx_bcast(fr->natoms*sizeof(fr->x[0]),fr->x[0],cr);
1057 gmx_bcast(fr->natoms*sizeof(fr->v[0]),fr->v[0],cr);
1062 double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
1063 const output_env_t oenv, bool bVerbose,bool bCompact,
1065 gmx_vsite_t *vsite,gmx_constr_t constr,
1066 int stepout,t_inputrec *ir,
1067 gmx_mtop_t *top_global,
1069 t_state *state_global,
1071 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
1072 gmx_edsam_t ed,t_forcerec *fr,
1073 int repl_ex_nst,int repl_ex_seed,
1074 real cpt_period,real max_hours,
1075 const char *deviceOptions,
1076 unsigned long Flags,
1077 gmx_runtime_t *runtime)
1080 gmx_large_int_t step,step_rel;
1083 bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
1084 bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
1085 bFirstStep,bStateFromTPX,bInitStep,bLastStep,
1086 bBornRadii,bStartingFromCpt;
1088 bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
1089 bForceUpdate=FALSE,bCPT;
1092 int force_flags,cglo_flags;
1093 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
1095 t_trxstatus *status;
1098 t_state *bufstate=NULL;
1099 matrix *scale_tot,pcoupl_mu,M,ebox;
1101 t_trxframe rerun_fr;
1102 gmx_repl_ex_t repl_ex=NULL;
1105 gmx_localtop_t *top;
1106 t_mdebin *mdebin=NULL;
1107 t_state *state=NULL;
1108 rvec *f_global=NULL;
1111 gmx_enerdata_t *enerd;
1113 gmx_global_stat_t gstat;
1114 gmx_update_t upd=NULL;
1115 t_graph *graph=NULL;
1119 gmx_groups_t *groups;
1120 gmx_ekindata_t *ekind, *ekind_save;
1121 gmx_shellfc_t shellfc;
1122 int count,nconverged=0;
1126 bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
1128 bool bResetCountersHalfMaxH=FALSE;
1129 bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
1130 real temp0,mu_aver=0,dvdl;
1132 atom_id *grpindex=NULL;
1134 t_coupl_rec *tcr=NULL;
1135 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
1136 matrix boxcopy={{0}},lastbox;
1138 real fom,oldfom,veta_save,pcurr,scalevir,tracevir;
1141 real saved_conserved_quantity = 0;
1146 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
1147 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
1148 gmx_iterate_t iterate;
1150 /* Temporary addition for FAHCORE checkpointing */
1154 /* Check for special mdrun options */
1155 bRerunMD = (Flags & MD_RERUN);
1156 bIonize = (Flags & MD_IONIZE);
1157 bFFscan = (Flags & MD_FFSCAN);
1158 bAppend = (Flags & MD_APPENDFILES);
1159 if (Flags & MD_RESETCOUNTERSHALFWAY)
1163 /* Signal to reset the counters half the simulation steps. */
1164 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
1166 /* Signal to reset the counters halfway the simulation time. */
1167 bResetCountersHalfMaxH = (max_hours > 0);
1170 /* md-vv uses averaged full step velocities for T-control
1171 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1172 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1173 bVV = EI_VV(ir->eI);
1174 if (bVV) /* to store the initial velocities while computing virial */
1176 snew(cbuf,top_global->natoms);
1178 /* all the iteratative cases - only if there are constraints */
1179 bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
1180 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
1184 /* Since we don't know if the frames read are related in any way,
1185 * rebuild the neighborlist at every step.
1188 ir->nstcalcenergy = 1;
1192 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
1194 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
1195 bGStatEveryStep = (nstglobalcomm == 1);
1197 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
1200 "To reduce the energy communication with nstlist = -1\n"
1201 "the neighbor list validity should not be checked at every step,\n"
1202 "this means that exact integration is not guaranteed.\n"
1203 "The neighbor list validity is checked after:\n"
1204 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1205 "In most cases this will result in exact integration.\n"
1206 "This reduces the energy communication by a factor of 2 to 3.\n"
1207 "If you want less energy communication, set nstlist > 3.\n\n");
1210 if (bRerunMD || bFFscan)
1214 groups = &top_global->groups;
1216 /* Initial values */
1217 init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
1218 nrnb,top_global,&upd,
1219 nfile,fnm,&outf,&mdebin,
1220 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
1222 clear_mat(total_vir);
1224 /* Energy terms and groups */
1226 init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
1227 if (DOMAINDECOMP(cr))
1233 snew(f,top_global->natoms);
1236 /* Kinetic energy data */
1238 init_ekindata(fplog,top_global,&(ir->opts),ekind);
1239 /* needed for iteration of constraints */
1241 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
1242 /* Copy the cos acceleration to the groups struct */
1243 ekind->cosacc.cos_accel = ir->cos_accel;
1245 gstat = global_stat_init(ir);
1248 /* Check for polarizable models and flexible constraints */
1249 shellfc = init_shell_flexcon(fplog,
1250 top_global,n_flexible_constraints(constr),
1251 (ir->bContinuation ||
1252 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
1253 NULL : state_global->x);
1258 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
1260 set_deform_reference_box(upd,
1261 deform_init_init_step_tpx,
1262 deform_init_box_tpx);
1264 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
1269 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
1270 if ((io > 2000) && MASTER(cr))
1272 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1276 if (DOMAINDECOMP(cr)) {
1277 top = dd_init_local_top(top_global);
1280 dd_init_local_state(cr->dd,state_global,state);
1282 if (DDMASTER(cr->dd) && ir->nstfout) {
1283 snew(f_global,state_global->natoms);
1287 /* Initialize the particle decomposition and split the topology */
1288 top = split_system(fplog,top_global,ir,cr);
1290 pd_cg_range(cr,&fr->cg0,&fr->hcg);
1291 pd_at_range(cr,&a0,&a1);
1293 top = gmx_mtop_generate_local_top(top_global,ir);
1296 a1 = top_global->natoms;
1299 state = partdec_init_local_state(cr,state_global);
1302 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
1305 set_vsite_top(vsite,top,mdatoms,cr);
1308 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
1309 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
1313 make_local_shells(cr,mdatoms,shellfc);
1316 if (ir->pull && PAR(cr)) {
1317 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
1321 if (DOMAINDECOMP(cr))
1323 /* Distribute the charge groups over the nodes from the master node */
1324 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
1325 state_global,top_global,ir,
1326 state,&f,mdatoms,top,fr,
1327 vsite,shellfc,constr,
1331 update_mdatoms(mdatoms,state->lambda);
1335 /* Update mdebin with energy history if appending to output files */
1336 if ( Flags & MD_APPENDFILES )
1338 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
1340 /* Set the initial energy history in state to zero by updating once */
1341 update_energyhistory(&state_global->enerhist,mdebin);
1344 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
1345 /* Set the random state if we read a checkpoint file */
1346 set_stochd_state(upd,state);
1349 /* Initialize constraints */
1351 if (!DOMAINDECOMP(cr))
1352 set_constraints(constr,top,ir,mdatoms,cr);
1355 /* Check whether we have to GCT stuff */
1356 bTCR = ftp2bSet(efGCT,nfile,fnm);
1359 fprintf(stderr,"Will do General Coupling Theory!\n");
1361 gnx = top_global->mols.nr;
1363 for(i=0; (i<gnx); i++) {
1368 if (repl_ex_nst > 0 && MASTER(cr))
1369 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
1370 repl_ex_nst,repl_ex_seed);
1372 if (!ir->bContinuation && !bRerunMD)
1374 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
1376 /* Set the velocities of frozen particles to zero */
1377 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
1379 for(m=0; m<DIM; m++)
1381 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
1391 /* Constrain the initial coordinates and velocities */
1392 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
1393 graph,cr,nrnb,fr,top,shake_vir);
1397 /* Construct the virtual sites for the initial configuration */
1398 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
1399 top->idef.iparams,top->idef.il,
1400 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1406 /* I'm assuming we need global communication the first time! MRS */
1407 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
1408 | (bVV ? CGLO_PRESSURE:0)
1409 | (bVV ? CGLO_CONSTRAINT:0)
1410 | (bRerunMD ? CGLO_RERUNMD:0)
1411 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
1413 bSumEkinhOld = FALSE;
1414 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1415 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1416 constr,NULL,FALSE,state->box,
1417 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
1418 if (ir->eI == eiVVAK) {
1419 /* a second call to get the half step temperature initialized as well */
1420 /* we do the same call as above, but turn the pressure off -- internally to
1421 compute_globals, this is recognized as a velocity verlet half-step
1422 kinetic energy calculation. This minimized excess variables, but
1423 perhaps loses some logic?*/
1425 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1426 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1427 constr,NULL,FALSE,state->box,
1428 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1429 cglo_flags &~ CGLO_PRESSURE);
1432 /* Calculate the initial half step temperature, and save the ekinh_old */
1433 if (!(Flags & MD_STARTFROMCPT))
1435 for(i=0; (i<ir->opts.ngtc); i++)
1437 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
1442 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
1443 and there is no previous step */
1445 temp0 = enerd->term[F_TEMP];
1447 /* if using an iterative algorithm, we need to create a working directory for the state. */
1450 bufstate = init_bufstate(state);
1454 snew(xcopy,state->natoms);
1455 snew(vcopy,state->natoms);
1456 copy_rvecn(state->x,xcopy,0,state->natoms);
1457 copy_rvecn(state->v,vcopy,0,state->natoms);
1458 copy_mat(state->box,boxcopy);
1461 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
1462 temperature control */
1463 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
1467 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
1470 "RMS relative constraint deviation after constraining: %.2e\n",
1471 constr_rmsd(constr,FALSE));
1473 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
1476 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
1477 " input trajectory '%s'\n\n",
1478 *(top_global->name),opt2fn("-rerun",nfile,fnm));
1481 fprintf(stderr,"Calculated time to finish depends on nsteps from "
1482 "run input file,\nwhich may not correspond to the time "
1483 "needed to process input trajectory.\n\n");
1489 fprintf(stderr,"starting mdrun '%s'\n",
1490 *(top_global->name));
1491 if (ir->nsteps >= 0)
1493 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
1497 sprintf(tbuf,"%s","infinite");
1499 if (ir->init_step > 0)
1501 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
1502 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
1503 gmx_step_str(ir->init_step,sbuf2),
1504 ir->init_step*ir->delta_t);
1508 fprintf(stderr,"%s steps, %s ps.\n",
1509 gmx_step_str(ir->nsteps,sbuf),tbuf);
1512 fprintf(fplog,"\n");
1515 /* Set and write start time */
1516 runtime_start(runtime);
1517 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
1518 wallcycle_start(wcycle,ewcRUN);
1520 fprintf(fplog,"\n");
1522 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1524 chkpt_ret=fcCheckPointParallel( cr->nodeid,
1526 if ( chkpt_ret == 0 )
1527 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
1531 /***********************************************************
1533 * Loop over MD steps
1535 ************************************************************/
1537 /* if rerunMD then read coordinates and velocities from input trajectory */
1540 if (getenv("GMX_FORCE_UPDATE"))
1542 bForceUpdate = TRUE;
1545 rerun_fr.natoms = 0;
1548 bNotLastFrame = read_first_frame(oenv,&status,
1549 opt2fn("-rerun",nfile,fnm),
1550 &rerun_fr,TRX_NEED_X | TRX_READ_V);
1551 if (rerun_fr.natoms != top_global->natoms)
1554 "Number of atoms in trajectory (%d) does not match the "
1555 "run input file (%d)\n",
1556 rerun_fr.natoms,top_global->natoms);
1558 if (ir->ePBC != epbcNONE)
1562 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);
1564 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
1566 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
1573 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
1576 if (ir->ePBC != epbcNONE)
1578 /* Set the shift vectors.
1579 * Necessary here when have a static box different from the tpr box.
1581 calc_shifts(rerun_fr.box,fr->shift_vec);
1585 /* loop over MD steps or if rerunMD to end of input trajectory */
1587 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1588 bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
1589 bInitStep = bFirstStep && (bStateFromTPX || bVV);
1590 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
1592 bSumEkinhOld = FALSE;
1595 init_global_signals(&gs,cr,ir,repl_ex_nst);
1597 step = ir->init_step;
1600 if (ir->nstlist == -1)
1602 init_nlistheuristics(&nlh,bGStatEveryStep,step);
1605 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
1606 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
1608 wallcycle_start(wcycle,ewcSTEP);
1610 GMX_MPE_LOG(ev_timestep1);
1613 if (rerun_fr.bStep) {
1614 step = rerun_fr.step;
1615 step_rel = step - ir->init_step;
1617 if (rerun_fr.bTime) {
1627 bLastStep = (step_rel == ir->nsteps);
1628 t = t0 + step*ir->delta_t;
1631 if (ir->efep != efepNO)
1633 if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
1635 state_global->lambda = rerun_fr.lambda;
1639 state_global->lambda = lam0 + step*ir->delta_lambda;
1641 state->lambda = state_global->lambda;
1642 bDoDHDL = do_per_step(step,ir->nstdhdl);
1647 update_annealing_target_temp(&(ir->opts),t);
1652 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
1654 for(i=0; i<state_global->natoms; i++)
1656 copy_rvec(rerun_fr.x[i],state_global->x[i]);
1660 for(i=0; i<state_global->natoms; i++)
1662 copy_rvec(rerun_fr.v[i],state_global->v[i]);
1667 for(i=0; i<state_global->natoms; i++)
1669 clear_rvec(state_global->v[i]);
1673 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
1674 " Ekin, temperature and pressure are incorrect,\n"
1675 " the virial will be incorrect when constraints are present.\n"
1677 bRerunWarnNoV = FALSE;
1681 copy_mat(rerun_fr.box,state_global->box);
1682 copy_mat(state_global->box,state->box);
1684 if (vsite && (Flags & MD_RERUN_VSITE))
1686 if (DOMAINDECOMP(cr))
1688 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1692 /* Following is necessary because the graph may get out of sync
1693 * with the coordinates if we only have every N'th coordinate set
1695 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
1696 shift_self(graph,state->box,state->x);
1698 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1699 top->idef.iparams,top->idef.il,
1700 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1703 unshift_self(graph,state->box,state->x);
1708 /* Stop Center of Mass motion */
1709 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
1711 /* Copy back starting coordinates in case we're doing a forcefield scan */
1714 for(ii=0; (ii<state->natoms); ii++)
1716 copy_rvec(xcopy[ii],state->x[ii]);
1717 copy_rvec(vcopy[ii],state->v[ii]);
1719 copy_mat(boxcopy,state->box);
1724 /* for rerun MD always do Neighbour Searching */
1725 bNS = (bFirstStep || ir->nstlist != 0);
1730 /* Determine whether or not to do Neighbour Searching and LR */
1731 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
1733 bNS = (bFirstStep || bExchanged || bNStList ||
1734 (ir->nstlist == -1 && nlh.nabnsb > 0));
1736 if (bNS && ir->nstlist == -1)
1738 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
1742 /* < 0 means stop at next step, > 0 means stop at next NS step */
1743 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
1744 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
1749 /* Determine whether or not to update the Born radii if doing GB */
1750 bBornRadii=bFirstStep;
1751 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
1756 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
1757 do_verbose = bVerbose &&
1758 (step % stepout == 0 || bFirstStep || bLastStep);
1760 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
1764 bMasterState = TRUE;
1768 bMasterState = FALSE;
1769 /* Correct the new box if it is too skewed */
1770 if (DYNAMIC_BOX(*ir))
1772 if (correct_box(fplog,step,state->box,graph))
1774 bMasterState = TRUE;
1777 if (DOMAINDECOMP(cr) && bMasterState)
1779 dd_collect_state(cr->dd,state,state_global);
1783 if (DOMAINDECOMP(cr))
1785 /* Repartition the domain decomposition */
1786 wallcycle_start(wcycle,ewcDOMDEC);
1787 dd_partition_system(fplog,step,cr,
1788 bMasterState,nstglobalcomm,
1789 state_global,top_global,ir,
1790 state,&f,mdatoms,top,fr,
1791 vsite,shellfc,constr,
1792 nrnb,wcycle,do_verbose);
1793 wallcycle_stop(wcycle,ewcDOMDEC);
1794 /* If using an iterative integrator, reallocate space to match the decomposition */
1798 if (MASTER(cr) && do_log && !bFFscan)
1800 print_ebin_header(fplog,step,t,state->lambda);
1803 if (ir->efep != efepNO)
1805 update_mdatoms(mdatoms,state->lambda);
1808 if (bRerunMD && rerun_fr.bV)
1811 /* We need the kinetic energy at minus the half step for determining
1812 * the full step kinetic energy and possibly for T-coupling.*/
1813 /* This may not be quite working correctly yet . . . . */
1814 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1815 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
1816 constr,NULL,FALSE,state->box,
1817 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1818 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1820 clear_mat(force_vir);
1822 /* Ionize the atoms if necessary */
1825 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
1826 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
1829 /* Update force field in ffscan program */
1832 if (update_forcefield(fplog,
1834 mdatoms->nr,state->x,state->box)) {
1835 if (gmx_parallel_env_initialized())
1843 GMX_MPE_LOG(ev_timestep2);
1845 /* We write a checkpoint at this MD step when:
1846 * either at an NS step when we signalled through gs,
1847 * or at the last step (but not when we do not want confout),
1848 * but never at the first step or with rerun.
1850 bCPT = (((gs.set[eglsCHKPT] && bNS) ||
1851 (bLastStep && (Flags & MD_CONFOUT))) &&
1852 step > ir->init_step && !bRerunMD);
1855 gs.set[eglsCHKPT] = 0;
1858 /* Determine the energy and pressure:
1859 * at nstcalcenergy steps and at energy output steps (set below).
1861 bNstEner = do_per_step(step,ir->nstcalcenergy);
1864 (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)));
1866 /* Do we need global communication ? */
1867 bGStat = (bCalcEnerPres || bStopCM ||
1868 do_per_step(step,nstglobalcomm) ||
1869 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1871 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1873 if (do_ene || do_log)
1875 bCalcEnerPres = TRUE;
1879 /* these CGLO_ options remain the same throughout the iteration */
1880 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1881 (bStopCM ? CGLO_STOPCM : 0) |
1882 (bGStat ? CGLO_GSTAT : 0)
1885 force_flags = (GMX_FORCE_STATECHANGED |
1886 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1887 GMX_FORCE_ALLFORCES |
1888 (bNStList ? GMX_FORCE_DOLR : 0) |
1890 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
1891 (bDoDHDL ? GMX_FORCE_DHDL : 0)
1896 /* Now is the time to relax the shells */
1897 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
1899 bStopCM,top,top_global,
1901 state,f,force_vir,mdatoms,
1902 nrnb,wcycle,graph,groups,
1903 shellfc,fr,bBornRadii,t,mu_tot,
1904 state->natoms,&bConverged,vsite,
1915 /* The coordinates (x) are shifted (to get whole molecules)
1917 * This is parallellized as well, and does communication too.
1918 * Check comments in sim_util.c
1921 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1922 state->box,state->x,&state->hist,
1923 f,force_vir,mdatoms,enerd,fcd,
1924 state->lambda,graph,
1925 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1926 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1929 GMX_BARRIER(cr->mpi_comm_mygroup);
1933 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1934 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1937 if (bTCR && bFirstStep)
1939 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1940 fprintf(fplog,"Done init_coupling\n");
1944 if (bVV && !bStartingFromCpt && !bRerunMD)
1945 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1947 if (ir->eI==eiVV && bInitStep)
1949 /* if using velocity verlet with full time step Ekin,
1950 * take the first half step only to compute the
1951 * virial for the first step. From there,
1952 * revert back to the initial coordinates
1953 * so that the input is actually the initial step.
1955 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1957 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1958 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1961 update_coords(fplog,step,ir,mdatoms,state,
1962 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1963 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1964 cr,nrnb,constr,&top->idef);
1968 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1970 /* for iterations, we save these vectors, as we will be self-consistently iterating
1973 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1975 /* save the state */
1976 if (bIterations && iterate.bIterate) {
1977 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1980 bFirstIterate = TRUE;
1981 while (bFirstIterate || (bIterations && iterate.bIterate))
1983 if (bIterations && iterate.bIterate)
1985 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1986 if (bFirstIterate && bTrotter)
1988 /* The first time through, we need a decent first estimate
1989 of veta(t+dt) to compute the constraints. Do
1990 this by computing the box volume part of the
1991 trotter integration at this time. Nothing else
1992 should be changed by this routine here. If
1993 !(first time), we start with the previous value
1996 veta_save = state->veta;
1997 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1998 vetanew = state->veta;
1999 state->veta = veta_save;
2004 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
2007 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2008 &top->idef,shake_vir,NULL,
2009 cr,nrnb,wcycle,upd,constr,
2010 bInitStep,TRUE,bCalcEnerPres,vetanew);
2012 if (!bOK && !bFFscan)
2014 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2019 { /* Need to unshift here if a do_force has been
2020 called in the previous step */
2021 unshift_self(graph,state->box,state->x);
2025 /* if VV, compute the pressure and constraints */
2026 /* For VV2, we strictly only need this if using pressure
2027 * control, but we really would like to have accurate pressures
2029 * Think about ways around this in the future?
2030 * For now, keep this choice in comments.
2032 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
2033 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
2035 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
2036 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2037 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2038 constr,NULL,FALSE,state->box,
2039 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2042 | (bTemp ? CGLO_TEMPERATURE:0)
2043 | (bPres ? CGLO_PRESSURE : 0)
2044 | (bPres ? CGLO_CONSTRAINT : 0)
2045 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
2046 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2049 /* explanation of above:
2050 a) We compute Ekin at the full time step
2051 if 1) we are using the AveVel Ekin, and it's not the
2052 initial step, or 2) if we are using AveEkin, but need the full
2053 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
2054 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2055 EkinAveVel because it's needed for the pressure */
2057 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2062 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
2066 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2071 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2072 state->veta,&vetanew))
2076 bFirstIterate = FALSE;
2079 if (bTrotter && !bInitStep) {
2080 copy_mat(shake_vir,state->svir_prev);
2081 copy_mat(force_vir,state->fvir_prev);
2082 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
2083 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2084 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
2085 enerd->term[F_EKIN] = trace(ekind->ekin);
2088 /* if it's the initial step, we performed this first step just to get the constraint virial */
2089 if (bInitStep && ir->eI==eiVV) {
2090 copy_rvecn(cbuf,state->v,0,state->natoms);
2093 if (fr->bSepDVDL && fplog && do_log)
2095 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2097 enerd->term[F_DHDL_CON] += dvdl;
2099 GMX_MPE_LOG(ev_timestep1);
2102 /* MRS -- now done iterating -- compute the conserved quantity */
2104 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
2107 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
2109 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
2111 saved_conserved_quantity -= enerd->term[F_DISPCORR];
2115 /* ######## END FIRST UPDATE STEP ############## */
2116 /* ######## If doing VV, we now have v(dt) ###### */
2118 /* ################## START TRAJECTORY OUTPUT ################# */
2120 /* Now we have the energies and forces corresponding to the
2121 * coordinates at time t. We must output all of this before
2123 * for RerunMD t is read from input trajectory
2125 GMX_MPE_LOG(ev_output_start);
2128 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
2129 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
2130 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
2131 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
2132 if (bCPT) { mdof_flags |= MDOF_CPT; };
2134 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
2137 /* Enforce writing positions and velocities at end of run */
2138 mdof_flags |= (MDOF_X | MDOF_V);
2143 fcReportProgress( ir->nsteps, step );
2145 /* sync bCPT and fc record-keeping */
2146 if (bCPT && MASTER(cr))
2147 fcRequestCheckPoint();
2150 if (mdof_flags != 0)
2152 wallcycle_start(wcycle,ewcTRAJ);
2155 if (state->flags & (1<<estLD_RNG))
2157 get_stochd_state(upd,state);
2163 state_global->ekinstate.bUpToDate = FALSE;
2167 update_ekinstate(&state_global->ekinstate,ekind);
2168 state_global->ekinstate.bUpToDate = TRUE;
2170 update_energyhistory(&state_global->enerhist,mdebin);
2173 write_traj(fplog,cr,outf,mdof_flags,top_global,
2174 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
2181 if (bLastStep && step_rel == ir->nsteps &&
2182 (Flags & MD_CONFOUT) && MASTER(cr) &&
2183 !bRerunMD && !bFFscan)
2185 /* x and v have been collected in write_traj,
2186 * because a checkpoint file will always be written
2189 fprintf(stderr,"\nWriting final coordinates.\n");
2190 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
2193 /* Make molecules whole only for confout writing */
2194 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
2196 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2197 *top_global->name,top_global,
2198 state_global->x,state_global->v,
2199 ir->ePBC,state->box);
2202 wallcycle_stop(wcycle,ewcTRAJ);
2204 GMX_MPE_LOG(ev_output_finish);
2206 /* kludge -- virial is lost with restart for NPT control. Must restart */
2207 if (bStartingFromCpt && bVV)
2209 copy_mat(state->svir_prev,shake_vir);
2210 copy_mat(state->fvir_prev,force_vir);
2212 /* ################## END TRAJECTORY OUTPUT ################ */
2214 /* Determine the pressure:
2215 * always when we want exact averages in the energy file,
2216 * at ns steps when we have pressure coupling,
2217 * otherwise only at energy output steps (set below).
2220 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2221 bCalcEnerPres = bNstEner;
2223 /* Do we need global communication ? */
2224 bGStat = (bGStatEveryStep || bStopCM || bNS ||
2225 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2227 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2229 if (do_ene || do_log)
2231 bCalcEnerPres = TRUE;
2235 /* Determine the wallclock run time up till now */
2236 run_time = gmx_gettime() - (double)runtime->real;
2238 /* Check whether everything is still allright */
2239 if (((int)gmx_get_stop_condition() > handled_stop_condition)
2245 /* this is just make gs.sig compatible with the hack
2246 of sending signals around by MPI_Reduce with together with
2248 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
2249 gs.sig[eglsSTOPCOND]=1;
2250 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
2251 gs.sig[eglsSTOPCOND]=-1;
2252 /* < 0 means stop at next step, > 0 means stop at next NS step */
2256 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2257 gmx_get_signal_name(),
2258 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2262 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2263 gmx_get_signal_name(),
2264 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2266 handled_stop_condition=(int)gmx_get_stop_condition();
2268 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
2269 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
2270 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
2272 /* Signal to terminate the run */
2273 gs.sig[eglsSTOPCOND] = 1;
2276 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2278 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2281 if (bResetCountersHalfMaxH && MASTER(cr) &&
2282 run_time > max_hours*60.0*60.0*0.495)
2284 gs.sig[eglsRESETCOUNTERS] = 1;
2287 if (ir->nstlist == -1 && !bRerunMD)
2289 /* When bGStatEveryStep=FALSE, global_stat is only called
2290 * when we check the atom displacements, not at NS steps.
2291 * This means that also the bonded interaction count check is not
2292 * performed immediately after NS. Therefore a few MD steps could
2293 * be performed with missing interactions.
2294 * But wrong energies are never written to file,
2295 * since energies are only written after global_stat
2298 if (step >= nlh.step_nscheck)
2300 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
2301 nlh.scale_tot,state->x);
2305 /* This is not necessarily true,
2306 * but step_nscheck is determined quite conservatively.
2312 /* In parallel we only have to check for checkpointing in steps
2313 * where we do global communication,
2314 * otherwise the other nodes don't know.
2316 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
2319 run_time >= nchkpt*cpt_period*60.0)) &&
2320 gs.set[eglsCHKPT] == 0)
2322 gs.sig[eglsCHKPT] = 1;
2327 gmx_iterate_init(&iterate,bIterations);
2330 /* for iterations, we save these vectors, as we will be redoing the calculations */
2331 if (bIterations && iterate.bIterate)
2333 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
2335 bFirstIterate = TRUE;
2336 while (bFirstIterate || (bIterations && iterate.bIterate))
2338 /* We now restore these vectors to redo the calculation with improved extended variables */
2341 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
2344 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2345 so scroll down for that logic */
2347 /* ######### START SECOND UPDATE STEP ################# */
2348 GMX_MPE_LOG(ev_update_start);
2349 /* Box is changed in update() when we do pressure coupling,
2350 * but we should still use the old box for energy corrections and when
2351 * writing it to the energy file, so it matches the trajectory files for
2352 * the same timestep above. Make a copy in a separate array.
2354 copy_mat(state->box,lastbox);
2357 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
2359 wallcycle_start(wcycle,ewcUPDATE);
2361 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
2364 if (bIterations && iterate.bIterate)
2372 /* we use a new value of scalevir to converge the iterations faster */
2373 scalevir = tracevir/trace(shake_vir);
2375 msmul(shake_vir,scalevir,shake_vir);
2376 m_add(force_vir,shake_vir,total_vir);
2377 clear_mat(shake_vir);
2379 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
2380 /* We can only do Berendsen coupling after we have summed
2381 * the kinetic energy or virial. Since the happens
2382 * in global_state after update, we should only do it at
2383 * step % nstlist = 1 with bGStatEveryStep=FALSE.
2388 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2389 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
2395 /* velocity half-step update */
2396 update_coords(fplog,step,ir,mdatoms,state,f,
2397 fr->bTwinRange && bNStList,fr->f_twin,fcd,
2398 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
2399 cr,nrnb,constr,&top->idef);
2402 /* Above, initialize just copies ekinh into ekin,
2403 * it doesn't copy position (for VV),
2404 * and entire integrator for MD.
2409 copy_rvecn(state->x,cbuf,0,state->natoms);
2412 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2413 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2414 wallcycle_stop(wcycle,ewcUPDATE);
2416 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2417 &top->idef,shake_vir,force_vir,
2418 cr,nrnb,wcycle,upd,constr,
2419 bInitStep,FALSE,bCalcEnerPres,state->veta);
2423 /* erase F_EKIN and F_TEMP here? */
2424 /* just compute the kinetic energy at the half step to perform a trotter step */
2425 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2426 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2427 constr,NULL,FALSE,lastbox,
2428 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2429 cglo_flags | CGLO_TEMPERATURE
2431 wallcycle_start(wcycle,ewcUPDATE);
2432 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
2433 /* now we know the scaling, we can compute the positions again again */
2434 copy_rvecn(cbuf,state->x,0,state->natoms);
2436 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2437 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2438 wallcycle_stop(wcycle,ewcUPDATE);
2440 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2441 /* are the small terms in the shake_vir here due
2442 * to numerical errors, or are they important
2443 * physically? I'm thinking they are just errors, but not completely sure.
2444 * For now, will call without actually constraining, constr=NULL*/
2445 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2446 &top->idef,tmp_vir,force_vir,
2447 cr,nrnb,wcycle,upd,NULL,
2448 bInitStep,FALSE,bCalcEnerPres,
2451 if (!bOK && !bFFscan)
2453 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2456 if (fr->bSepDVDL && fplog && do_log)
2458 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2460 enerd->term[F_DHDL_CON] += dvdl;
2464 /* Need to unshift here */
2465 unshift_self(graph,state->box,state->x);
2468 GMX_BARRIER(cr->mpi_comm_mygroup);
2469 GMX_MPE_LOG(ev_update_finish);
2473 wallcycle_start(wcycle,ewcVSITECONSTR);
2476 shift_self(graph,state->box,state->x);
2478 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
2479 top->idef.iparams,top->idef.il,
2480 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2484 unshift_self(graph,state->box,state->x);
2486 wallcycle_stop(wcycle,ewcVSITECONSTR);
2489 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
2490 if (ir->nstlist == -1 && bFirstIterate)
2492 gs.sig[eglsNABNSB] = nlh.nabnsb;
2494 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2495 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2497 bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
2499 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2501 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
2502 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
2503 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
2504 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
2505 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2508 if (ir->nstlist == -1 && bFirstIterate)
2510 nlh.nabnsb = gs.set[eglsNABNSB];
2511 gs.set[eglsNABNSB] = 0;
2513 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
2514 /* ############# END CALC EKIN AND PRESSURE ################# */
2516 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
2517 the virial that should probably be addressed eventually. state->veta has better properies,
2518 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
2519 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
2522 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2523 trace(shake_vir),&tracevir))
2527 bFirstIterate = FALSE;
2530 update_box(fplog,step,ir,mdatoms,state,graph,f,
2531 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
2533 /* ################# END UPDATE STEP 2 ################# */
2534 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2536 /* The coordinates (x) were unshifted in update */
2537 if (bFFscan && (shellfc==NULL || bConverged))
2539 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
2541 &(top_global->mols),mdatoms->massT,pres))
2543 if (gmx_parallel_env_initialized())
2547 fprintf(stderr,"\n");
2553 /* We will not sum ekinh_old,
2554 * so signal that we still have to do it.
2556 bSumEkinhOld = TRUE;
2561 /* Only do GCT when the relaxation of shells (minimization) has converged,
2562 * otherwise we might be coupling to bogus energies.
2563 * In parallel we must always do this, because the other sims might
2567 /* Since this is called with the new coordinates state->x, I assume
2568 * we want the new box state->box too. / EL 20040121
2570 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
2572 mdatoms,&(top->idef),mu_aver,
2573 top_global->mols.nr,cr,
2574 state->box,total_vir,pres,
2575 mu_tot,state->x,f,bConverged);
2579 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2581 /* sum up the foreign energy and dhdl terms */
2582 sum_dhdl(enerd,state->lambda,ir);
2584 /* use the directly determined last velocity, not actually the averaged half steps */
2585 if (bTrotter && ir->eI==eiVV)
2587 enerd->term[F_EKIN] = last_ekin;
2589 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
2593 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
2597 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
2599 /* Check for excessively large energies */
2603 real etot_max = 1e200;
2605 real etot_max = 1e30;
2607 if (fabs(enerd->term[F_ETOT]) > etot_max)
2609 fprintf(stderr,"Energy too large (%g), giving up\n",
2610 enerd->term[F_ETOT]);
2613 /* ######### END PREPARING EDR OUTPUT ########### */
2615 /* Time for performance */
2616 if (((step % stepout) == 0) || bLastStep)
2618 runtime_upd_proc(runtime);
2626 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
2630 upd_mdebin(mdebin,bDoDHDL, TRUE,
2631 t,mdatoms->tmass,enerd,state,lastbox,
2632 shake_vir,force_vir,total_vir,pres,
2633 ekind,mu_tot,constr);
2637 upd_mdebin_step(mdebin);
2640 do_dr = do_per_step(step,ir->nstdisreout);
2641 do_or = do_per_step(step,ir->nstorireout);
2643 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
2645 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
2647 if (ir->ePull != epullNO)
2649 pull_print_output(ir->pull,step,t);
2652 if (do_per_step(step,ir->nstlog))
2654 if(fflush(fplog) != 0)
2656 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
2662 /* Remaining runtime */
2663 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
2667 fprintf(stderr,"\n");
2669 print_time(stderr,runtime,step,ir,cr);
2672 /* Replica exchange */
2674 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
2675 do_per_step(step,repl_ex_nst))
2677 bExchanged = replica_exchange(fplog,cr,repl_ex,
2678 state_global,enerd->term,
2681 if (bExchanged && DOMAINDECOMP(cr))
2683 dd_partition_system(fplog,step,cr,TRUE,1,
2684 state_global,top_global,ir,
2685 state,&f,mdatoms,top,fr,
2686 vsite,shellfc,constr,
2693 bStartingFromCpt = FALSE;
2695 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2696 /* Complicated conditional when bGStatEveryStep=FALSE.
2697 * We can not just use bGStat, since then the simulation results
2698 * would depend on nstenergy and nstlog or step_nscheck.
2700 if (((state->flags & (1<<estPRES_PREV)) ||
2701 (state->flags & (1<<estSVIR_PREV)) ||
2702 (state->flags & (1<<estFVIR_PREV))) &&
2704 (ir->nstlist > 0 && step % ir->nstlist == 0) ||
2705 (ir->nstlist < 0 && nlh.nabnsb > 0) ||
2706 (ir->nstlist == 0 && bGStat)))
2708 /* Store the pressure in t_state for pressure coupling
2709 * at the next MD step.
2711 if (state->flags & (1<<estPRES_PREV))
2713 copy_mat(pres,state->pres_prev);
2717 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2723 /* read next frame from input trajectory */
2724 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
2729 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
2733 if (!bRerunMD || !rerun_fr.bStep)
2735 /* increase the MD step number */
2740 cycles = wallcycle_stop(wcycle,ewcSTEP);
2741 if (DOMAINDECOMP(cr) && wcycle)
2743 dd_cycles_add(cr->dd,cycles,ddCyclStep);
2746 if (step_rel == wcycle_get_reset_counters(wcycle) ||
2747 gs.set[eglsRESETCOUNTERS] != 0)
2749 /* Reset all the counters related to performance over the run */
2750 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
2751 wcycle_set_reset_counters(wcycle,-1);
2752 /* Correct max_hours for the elapsed time */
2753 max_hours -= run_time/(60.0*60.0);
2754 bResetCountersHalfMaxH = FALSE;
2755 gs.set[eglsRESETCOUNTERS] = 0;
2758 /* End of main MD loop */
2762 runtime_end(runtime);
2764 if (bRerunMD && MASTER(cr))
2769 if (!(cr->duty & DUTY_PME))
2771 /* Tell the PME only node to finish */
2777 if (ir->nstcalcenergy > 0 && !bRerunMD)
2779 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
2780 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
2788 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
2790 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)));
2791 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
2794 if (shellfc && fplog)
2796 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
2797 (nconverged*100.0)/step_rel);
2798 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
2802 if (repl_ex_nst > 0 && MASTER(cr))
2804 print_replica_exchange_statistics(fplog,repl_ex);
2807 runtime->nsteps_done = step_rel;