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 gmx_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)
124 gmx_bool bPos,bEqual;
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,gmx_bool bInterSimGS,
310 matrix box, gmx_mtop_t *top_global, real *pcurr,
311 int natoms, gmx_bool *bSumEkinhOld, int flags)
315 tensor corr_vir,corr_pres,shakeall_vir;
316 gmx_bool bEner,bPres,bTemp, bVV;
317 gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
318 bFirstIterate,bReadEkin,bEkinAveVel,bScaleEkin, bConstrain;
319 real ekin,temp,prescorr,enercorr,dvdlcorr;
321 /* translate CGLO flags to gmx_booleans */
322 bRerunMD = flags & CGLO_RERUNMD;
323 bStopCM = flags & CGLO_STOPCM;
324 bGStat = flags & CGLO_GSTAT;
326 bReadEkin = (flags & CGLO_READEKIN);
327 bScaleEkin = (flags & CGLO_SCALEEKIN);
328 bEner = flags & CGLO_ENERGY;
329 bTemp = flags & CGLO_TEMPERATURE;
330 bPres = (flags & CGLO_PRESSURE);
331 bConstrain = (flags & CGLO_CONSTRAINT);
332 bIterate = (flags & CGLO_ITERATE);
333 bFirstIterate = (flags & CGLO_FIRSTITERATE);
335 /* we calculate a full state kinetic energy either with full-step velocity verlet
336 or half step where we need the pressure */
338 bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && bPres) || bReadEkin);
340 /* in initalization, it sums the shake virial in vv, and to
341 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
343 /* ########## Kinetic energy ############## */
347 /* Non-equilibrium MD: this is parallellized, but only does communication
348 * when there really is NEMD.
351 if (PAR(cr) && (ekind->bNEMD))
353 accumulate_u(cr,&(ir->opts),ekind);
358 restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
363 calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
368 /* Calculate center of mass velocity if necessary, also parallellized */
369 if (bStopCM && !bRerunMD && bEner)
371 calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
372 state->x,state->v,vcm);
376 if (bTemp || bPres || bEner || bConstrain)
380 /* We will not sum ekinh_old,
381 * so signal that we still have to do it.
383 *bSumEkinhOld = TRUE;
390 for(i=0; i<eglsNR; i++)
392 gs_buf[i] = gs->sig[i];
397 wallcycle_start(wcycle,ewcMoveE);
398 GMX_MPE_LOG(ev_global_stat_start);
399 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
401 gs != NULL ? eglsNR : 0,gs_buf,
403 *bSumEkinhOld,flags);
404 GMX_MPE_LOG(ev_global_stat_finish);
405 wallcycle_stop(wcycle,ewcMoveE);
409 if (MULTISIM(cr) && bInterSimGS)
413 /* Communicate the signals between the simulations */
414 gmx_sum_sim(eglsNR,gs_buf,cr->ms);
416 /* Communicate the signals form the master to the others */
417 gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
419 for(i=0; i<eglsNR; i++)
421 if (bInterSimGS || gs_simlocal[i])
423 /* Set the communicated signal only when it is non-zero,
424 * since signals might not be processed at each MD step.
426 gsi = (gs_buf[i] >= 0 ?
427 (int)(gs_buf[i] + 0.5) :
428 (int)(gs_buf[i] - 0.5));
433 /* Turn off the local signal */
438 *bSumEkinhOld = FALSE;
442 if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
445 mdatoms->start,mdatoms->start+mdatoms->homenr,
446 state->v,vcm->group_p[0],
447 mdatoms->massT,mdatoms->tmass,ekind->ekin);
451 /* Do center of mass motion removal */
452 if (bStopCM && !bRerunMD) /* is this correct? Does it get called too often with this logic? */
454 check_cm_grp(fplog,vcm,ir,1);
455 do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
456 state->x,state->v,vcm);
457 inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
463 /* Sum the kinetic energies of the groups & calc temp */
464 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
465 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
466 Leap with AveVel is not supported; it's not clear that it will actually work.
467 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
468 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
469 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
470 If FALSE, we go ahead and erase over it.
472 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,&(enerd->term[F_DKDL]),
473 bEkinAveVel,bIterate,bScaleEkin);
475 enerd->term[F_EKIN] = trace(ekind->ekin);
478 /* ########## Long range energy information ###### */
480 if (bEner || bPres || bConstrain)
482 calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda,
483 corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
486 if (bEner && bFirstIterate)
488 enerd->term[F_DISPCORR] = enercorr;
489 enerd->term[F_EPOT] += enercorr;
490 enerd->term[F_DVDL] += dvdlcorr;
491 if (fr->efep != efepNO) {
492 enerd->dvdl_lin += dvdlcorr;
496 /* ########## Now pressure ############## */
497 if (bPres || bConstrain)
500 m_add(force_vir,shake_vir,total_vir);
502 /* Calculate pressure and apply LR correction if PPPM is used.
503 * Use the box from last timestep since we already called update().
506 enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres,
507 (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
509 /* Calculate long range corrections to pressure and energy */
510 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
511 and computes enerd->term[F_DISPCORR]. Also modifies the
512 total_vir and pres tesors */
514 m_add(total_vir,corr_vir,total_vir);
515 m_add(pres,corr_pres,pres);
516 enerd->term[F_PDISPCORR] = prescorr;
517 enerd->term[F_PRES] += prescorr;
518 *pcurr = enerd->term[F_PRES];
519 /* calculate temperature using virial */
520 enerd->term[F_VTEMP] = calc_temp(trace(total_vir),ir->opts.nrdf[0]);
526 /* Definitions for convergence of iterated constraints */
528 /* iterate constraints up to 50 times */
529 #define MAXITERCONST 50
534 real f,fprev,x,xprev;
537 real allrelerr[MAXITERCONST+2];
538 int num_close; /* number of "close" violations, caused by limited precision. */
542 #define CONVERGEITER 0.000000001
543 #define CLOSE_ENOUGH 0.000001000
545 #define CONVERGEITER 0.0001
546 #define CLOSE_ENOUGH 0.0050
549 /* we want to keep track of the close calls. If there are too many, there might be some other issues.
550 so we make sure that it's either less than some predetermined number, or if more than that number,
551 only some small fraction of the total. */
552 #define MAX_NUMBER_CLOSE 50
553 #define FRACTION_CLOSE 0.001
555 /* maximum length of cyclic traps to check, emerging from limited numerical precision */
558 static void gmx_iterate_init(gmx_iterate_t *iterate,gmx_bool bIterate)
563 iterate->bIterate = bIterate;
564 iterate->num_close = 0;
565 for (i=0;i<MAXITERCONST+2;i++)
567 iterate->allrelerr[i] = 0;
571 static gmx_bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate_t *iterate, gmx_bool bFirstIterate, real fom, real *newf)
573 /* monitor convergence, and use a secant search to propose new
576 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
577 f(x_{i}) - f(x_{i-1})
579 The function we are trying to zero is fom-x, where fom is the
580 "figure of merit" which is the pressure (or the veta value) we
581 would get by putting in an old value of the pressure or veta into
582 the incrementor function for the step or half step. I have
583 verified that this gives the same answer as self consistent
584 iteration, usually in many fewer steps, especially for small tau_p.
586 We could possibly eliminate an iteration with proper use
587 of the value from the previous step, but that would take a bit
588 more bookkeeping, especially for veta, since tests indicate the
589 function of veta on the last step is not sufficiently close to
590 guarantee convergence this step. This is
591 good enough for now. On my tests, I could use tau_p down to
592 0.02, which is smaller that would ever be necessary in
593 practice. Generally, 3-5 iterations will be sufficient */
595 real relerr,err,xmin;
603 iterate->f = fom-iterate->x;
610 iterate->f = fom-iterate->x; /* we want to zero this difference */
611 if ((iterate->iter_i > 1) && (iterate->iter_i < MAXITERCONST))
613 if (iterate->f==iterate->fprev)
619 *newf = iterate->x - (iterate->x-iterate->xprev)*(iterate->f)/(iterate->f-iterate->fprev);
624 /* just use self-consistent iteration the first step to initialize, or
625 if it's not converging (which happens occasionally -- need to investigate why) */
629 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
630 difference between the closest of x and xprev to the new
631 value. To be 100% certain, we should check the difference between
632 the last result, and the previous result, or
634 relerr = (fabs((x-xprev)/fom));
636 but this is pretty much never necessary under typical conditions.
637 Checking numerically, it seems to lead to almost exactly the same
638 trajectories, but there are small differences out a few decimal
639 places in the pressure, and eventually in the v_eta, but it could
642 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
643 relerr = (fabs((*newf-xmin) / *newf));
646 err = fabs((iterate->f-iterate->fprev));
647 relerr = fabs(err/fom);
649 iterate->allrelerr[iterate->iter_i] = relerr;
651 if (iterate->iter_i > 0)
655 fprintf(debug,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
656 iterate->iter_i,fom,relerr,*newf);
659 if ((relerr < CONVERGEITER) || (err < CONVERGEITER) || (fom==0) || ((iterate->x == iterate->xprev) && iterate->iter_i > 1))
661 iterate->bIterate = FALSE;
664 fprintf(debug,"Iterating NPT constraints: CONVERGED\n");
668 if (iterate->iter_i > MAXITERCONST)
670 if (relerr < CLOSE_ENOUGH)
673 for (i=1;i<CYCLEMAX;i++) {
674 if ((iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-1]) &&
675 (iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-(1+2*i)])) {
679 fprintf(debug,"Exiting from an NPT iterating cycle of length %d\n",i);
686 /* step 1: trapped in a numerical attractor */
687 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
688 Better to give up convergence here than have the simulation die.
690 iterate->num_close++;
695 /* Step #2: test if we are reasonably close for other reasons, then monitor the number. If not, die */
697 /* how many close calls have we had? If less than a few, we're OK */
698 if (iterate->num_close < MAX_NUMBER_CLOSE)
700 sprintf(buf,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps,relerr);
701 md_print_warning(cr,fplog,buf);
702 iterate->num_close++;
704 /* if more than a few, check the total fraction. If too high, die. */
705 } else if (iterate->num_close/(double)nsteps > FRACTION_CLOSE) {
706 gmx_fatal(FARGS,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate->num_close/(double)nsteps);
712 gmx_fatal(FARGS,"Could not converge NPT constraints\n");
717 iterate->xprev = iterate->x;
719 iterate->fprev = iterate->f;
725 static void check_nst_param(FILE *fplog,t_commrec *cr,
726 const char *desc_nst,int nst,
727 const char *desc_p,int *p)
731 if (*p > 0 && *p % nst != 0)
733 /* Round up to the next multiple of nst */
734 *p = ((*p)/nst + 1)*nst;
735 sprintf(buf,"NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
736 md_print_warning(cr,fplog,buf);
740 static void reset_all_counters(FILE *fplog,t_commrec *cr,
741 gmx_large_int_t step,
742 gmx_large_int_t *step_rel,t_inputrec *ir,
743 gmx_wallcycle_t wcycle,t_nrnb *nrnb,
744 gmx_runtime_t *runtime)
746 char buf[STRLEN],sbuf[STEPSTRSIZE];
748 /* Reset all the counters related to performance over the run */
749 sprintf(buf,"Step %s: resetting all time and cycle counters\n",
750 gmx_step_str(step,sbuf));
751 md_print_warning(cr,fplog,buf);
753 wallcycle_stop(wcycle,ewcRUN);
754 wallcycle_reset_all(wcycle);
755 if (DOMAINDECOMP(cr))
757 reset_dd_statistics_counters(cr->dd);
760 ir->init_step += *step_rel;
761 ir->nsteps -= *step_rel;
763 wallcycle_start(wcycle,ewcRUN);
764 runtime_start(runtime);
765 print_date_and_time(fplog,cr->nodeid,"Restarted time",runtime);
768 static void min_zero(int *n,int i)
770 if (i > 0 && (*n == 0 || i < *n))
776 static int lcd4(int i1,int i2,int i3,int i4)
787 gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
790 while (nst > 1 && ((i1 > 0 && i1 % nst != 0) ||
791 (i2 > 0 && i2 % nst != 0) ||
792 (i3 > 0 && i3 % nst != 0) ||
793 (i4 > 0 && i4 % nst != 0)))
801 static int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
802 int nstglobalcomm,t_inputrec *ir)
806 if (!EI_DYNAMICS(ir->eI))
811 if (nstglobalcomm == -1)
813 if (!(ir->nstcalcenergy > 0 ||
819 if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
821 nstglobalcomm = ir->nstenergy;
826 /* Ensure that we do timely global communication for
827 * (possibly) each of the four following options.
829 nstglobalcomm = lcd4(ir->nstcalcenergy,
831 ir->etc != etcNO ? ir->nsttcouple : 0,
832 ir->epc != epcNO ? ir->nstpcouple : 0);
837 if (ir->nstlist > 0 &&
838 nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
840 nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
841 sprintf(buf,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
842 md_print_warning(cr,fplog,buf);
844 if (ir->nstcalcenergy > 0)
846 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
847 "nstcalcenergy",&ir->nstcalcenergy);
849 if (ir->etc != etcNO && ir->nsttcouple > 0)
851 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
852 "nsttcouple",&ir->nsttcouple);
854 if (ir->epc != epcNO && ir->nstpcouple > 0)
856 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
857 "nstpcouple",&ir->nstpcouple);
860 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
861 "nstenergy",&ir->nstenergy);
863 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
864 "nstlog",&ir->nstlog);
867 if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
869 sprintf(buf,"WARNING: Changing nstcomm from %d to %d\n",
870 ir->nstcomm,nstglobalcomm);
871 md_print_warning(cr,fplog,buf);
872 ir->nstcomm = nstglobalcomm;
875 return nstglobalcomm;
878 void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
879 t_inputrec *ir,gmx_mtop_t *mtop)
881 /* Check required for old tpx files */
882 if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
883 ir->nstcalcenergy % ir->nstlist != 0)
885 md_print_warning(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
887 if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
888 gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
889 ir->eConstrAlg == econtSHAKE)
891 md_print_warning(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
892 if (ir->epc != epcNO)
894 gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
897 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
898 "nstcalcenergy",&ir->nstcalcenergy);
899 if (ir->epc != epcNO)
901 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
902 "nstpcouple",&ir->nstpcouple);
904 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
905 "nstenergy",&ir->nstenergy);
906 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
907 "nstlog",&ir->nstlog);
908 if (ir->efep != efepNO)
910 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
911 "nstdhdl",&ir->nstdhdl);
917 gmx_bool bGStatEveryStep;
918 gmx_large_int_t step_ns;
919 gmx_large_int_t step_nscheck;
930 static void reset_nlistheuristics(gmx_nlheur_t *nlh,gmx_large_int_t step)
934 nlh->step_nscheck = step;
937 static void init_nlistheuristics(gmx_nlheur_t *nlh,
938 gmx_bool bGStatEveryStep,gmx_large_int_t step)
940 nlh->bGStatEveryStep = bGStatEveryStep;
947 reset_nlistheuristics(nlh,step);
950 static void update_nliststatistics(gmx_nlheur_t *nlh,gmx_large_int_t step)
952 gmx_large_int_t nl_lt;
953 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
955 /* Determine the neighbor list life time */
956 nl_lt = step - nlh->step_ns;
959 fprintf(debug,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh->nabnsb,gmx_step_str(nl_lt,sbuf));
963 nlh->s2 += nl_lt*nl_lt;
964 nlh->ab += nlh->nabnsb;
965 if (nlh->lt_runav == 0)
967 nlh->lt_runav = nl_lt;
968 /* Initialize the fluctuation average
969 * such that at startup we check after 0 steps.
971 nlh->lt_runav2 = sqr(nl_lt/2.0);
973 /* Running average with 0.9 gives an exp. history of 9.5 */
974 nlh->lt_runav2 = 0.9*nlh->lt_runav2 + 0.1*sqr(nlh->lt_runav - nl_lt);
975 nlh->lt_runav = 0.9*nlh->lt_runav + 0.1*nl_lt;
976 if (nlh->bGStatEveryStep)
978 /* Always check the nlist validity */
979 nlh->step_nscheck = step;
983 /* We check after: <life time> - 2*sigma
984 * The factor 2 is quite conservative,
985 * but we assume that with nstlist=-1 the user
986 * prefers exact integration over performance.
988 nlh->step_nscheck = step
989 + (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)) - 1;
993 fprintf(debug,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
994 gmx_step_str(nl_lt,sbuf),nlh->lt_runav,sqrt(nlh->lt_runav2),
995 gmx_step_str(nlh->step_nscheck-step+1,sbuf2),
996 (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)));
1000 static void set_nlistheuristics(gmx_nlheur_t *nlh,gmx_bool bReset,gmx_large_int_t step)
1006 reset_nlistheuristics(nlh,step);
1010 update_nliststatistics(nlh,step);
1013 nlh->step_ns = step;
1014 /* Initialize the cumulative coordinate scaling matrix */
1015 clear_mat(nlh->scale_tot);
1016 for(d=0; d<DIM; d++)
1018 nlh->scale_tot[d][d] = 1.0;
1022 static void rerun_parallel_comm(t_commrec *cr,t_trxframe *fr,
1023 gmx_bool *bNotLastFrame)
1028 bAlloc = (fr->natoms == 0);
1030 if (MASTER(cr) && !*bNotLastFrame)
1036 gmx_bcast(sizeof(*fr),fr,cr);
1040 *bNotLastFrame = (fr->natoms >= 0);
1042 if (*bNotLastFrame && PARTDECOMP(cr))
1044 /* x and v are the only variable size quantities stored in trr
1045 * that are required for rerun (f is not needed).
1049 snew(fr->x,fr->natoms);
1050 snew(fr->v,fr->natoms);
1054 gmx_bcast(fr->natoms*sizeof(fr->x[0]),fr->x[0],cr);
1058 gmx_bcast(fr->natoms*sizeof(fr->v[0]),fr->v[0],cr);
1063 double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
1064 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
1066 gmx_vsite_t *vsite,gmx_constr_t constr,
1067 int stepout,t_inputrec *ir,
1068 gmx_mtop_t *top_global,
1070 t_state *state_global,
1072 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
1073 gmx_edsam_t ed,t_forcerec *fr,
1074 int repl_ex_nst,int repl_ex_seed,
1075 real cpt_period,real max_hours,
1076 const char *deviceOptions,
1077 unsigned long Flags,
1078 gmx_runtime_t *runtime)
1081 gmx_large_int_t step,step_rel;
1084 gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
1085 gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
1086 bFirstStep,bStateFromTPX,bInitStep,bLastStep,
1087 bBornRadii,bStartingFromCpt;
1088 gmx_bool bDoDHDL=FALSE;
1089 gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
1090 bForceUpdate=FALSE,bCPT;
1092 gmx_bool bMasterState;
1093 int force_flags,cglo_flags;
1094 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
1096 t_trxstatus *status;
1099 t_state *bufstate=NULL;
1100 matrix *scale_tot,pcoupl_mu,M,ebox;
1102 t_trxframe rerun_fr;
1103 gmx_repl_ex_t repl_ex=NULL;
1106 gmx_localtop_t *top;
1107 t_mdebin *mdebin=NULL;
1108 t_state *state=NULL;
1109 rvec *f_global=NULL;
1112 gmx_enerdata_t *enerd;
1114 gmx_global_stat_t gstat;
1115 gmx_update_t upd=NULL;
1116 t_graph *graph=NULL;
1120 gmx_groups_t *groups;
1121 gmx_ekindata_t *ekind, *ekind_save;
1122 gmx_shellfc_t shellfc;
1123 int count,nconverged=0;
1126 gmx_bool bIonize=FALSE;
1127 gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
1129 gmx_bool bResetCountersHalfMaxH=FALSE;
1130 gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
1131 real temp0,mu_aver=0,dvdl;
1133 atom_id *grpindex=NULL;
1135 t_coupl_rec *tcr=NULL;
1136 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
1137 matrix boxcopy={{0}},lastbox;
1139 real fom,oldfom,veta_save,pcurr,scalevir,tracevir;
1142 real saved_conserved_quantity = 0;
1147 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
1148 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
1149 gmx_iterate_t iterate;
1151 /* Temporary addition for FAHCORE checkpointing */
1155 /* Check for special mdrun options */
1156 bRerunMD = (Flags & MD_RERUN);
1157 bIonize = (Flags & MD_IONIZE);
1158 bFFscan = (Flags & MD_FFSCAN);
1159 bAppend = (Flags & MD_APPENDFILES);
1160 if (Flags & MD_RESETCOUNTERSHALFWAY)
1164 /* Signal to reset the counters half the simulation steps. */
1165 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
1167 /* Signal to reset the counters halfway the simulation time. */
1168 bResetCountersHalfMaxH = (max_hours > 0);
1171 /* md-vv uses averaged full step velocities for T-control
1172 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1173 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1174 bVV = EI_VV(ir->eI);
1175 if (bVV) /* to store the initial velocities while computing virial */
1177 snew(cbuf,top_global->natoms);
1179 /* all the iteratative cases - only if there are constraints */
1180 bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
1181 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
1185 /* Since we don't know if the frames read are related in any way,
1186 * rebuild the neighborlist at every step.
1189 ir->nstcalcenergy = 1;
1193 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
1195 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
1196 bGStatEveryStep = (nstglobalcomm == 1);
1198 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
1201 "To reduce the energy communication with nstlist = -1\n"
1202 "the neighbor list validity should not be checked at every step,\n"
1203 "this means that exact integration is not guaranteed.\n"
1204 "The neighbor list validity is checked after:\n"
1205 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1206 "In most cases this will result in exact integration.\n"
1207 "This reduces the energy communication by a factor of 2 to 3.\n"
1208 "If you want less energy communication, set nstlist > 3.\n\n");
1211 if (bRerunMD || bFFscan)
1215 groups = &top_global->groups;
1217 /* Initial values */
1218 init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
1219 nrnb,top_global,&upd,
1220 nfile,fnm,&outf,&mdebin,
1221 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
1223 clear_mat(total_vir);
1225 /* Energy terms and groups */
1227 init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
1228 if (DOMAINDECOMP(cr))
1234 snew(f,top_global->natoms);
1237 /* Kinetic energy data */
1239 init_ekindata(fplog,top_global,&(ir->opts),ekind);
1240 /* needed for iteration of constraints */
1242 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
1243 /* Copy the cos acceleration to the groups struct */
1244 ekind->cosacc.cos_accel = ir->cos_accel;
1246 gstat = global_stat_init(ir);
1249 /* Check for polarizable models and flexible constraints */
1250 shellfc = init_shell_flexcon(fplog,
1251 top_global,n_flexible_constraints(constr),
1252 (ir->bContinuation ||
1253 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
1254 NULL : state_global->x);
1259 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
1261 set_deform_reference_box(upd,
1262 deform_init_init_step_tpx,
1263 deform_init_box_tpx);
1265 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
1270 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
1271 if ((io > 2000) && MASTER(cr))
1273 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1277 if (DOMAINDECOMP(cr)) {
1278 top = dd_init_local_top(top_global);
1281 dd_init_local_state(cr->dd,state_global,state);
1283 if (DDMASTER(cr->dd) && ir->nstfout) {
1284 snew(f_global,state_global->natoms);
1288 /* Initialize the particle decomposition and split the topology */
1289 top = split_system(fplog,top_global,ir,cr);
1291 pd_cg_range(cr,&fr->cg0,&fr->hcg);
1292 pd_at_range(cr,&a0,&a1);
1294 top = gmx_mtop_generate_local_top(top_global,ir);
1297 a1 = top_global->natoms;
1300 state = partdec_init_local_state(cr,state_global);
1303 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
1306 set_vsite_top(vsite,top,mdatoms,cr);
1309 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
1310 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
1314 make_local_shells(cr,mdatoms,shellfc);
1317 if (ir->pull && PAR(cr)) {
1318 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
1322 if (DOMAINDECOMP(cr))
1324 /* Distribute the charge groups over the nodes from the master node */
1325 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
1326 state_global,top_global,ir,
1327 state,&f,mdatoms,top,fr,
1328 vsite,shellfc,constr,
1332 update_mdatoms(mdatoms,state->lambda);
1336 /* Update mdebin with energy history if appending to output files */
1337 if ( Flags & MD_APPENDFILES )
1339 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
1341 /* Set the initial energy history in state to zero by updating once */
1342 update_energyhistory(&state_global->enerhist,mdebin);
1345 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
1346 /* Set the random state if we read a checkpoint file */
1347 set_stochd_state(upd,state);
1350 /* Initialize constraints */
1352 if (!DOMAINDECOMP(cr))
1353 set_constraints(constr,top,ir,mdatoms,cr);
1356 /* Check whether we have to GCT stuff */
1357 bTCR = ftp2bSet(efGCT,nfile,fnm);
1360 fprintf(stderr,"Will do General Coupling Theory!\n");
1362 gnx = top_global->mols.nr;
1364 for(i=0; (i<gnx); i++) {
1369 if (repl_ex_nst > 0 && MASTER(cr))
1370 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
1371 repl_ex_nst,repl_ex_seed);
1373 if (!ir->bContinuation && !bRerunMD)
1375 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
1377 /* Set the velocities of frozen particles to zero */
1378 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
1380 for(m=0; m<DIM; m++)
1382 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
1392 /* Constrain the initial coordinates and velocities */
1393 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
1394 graph,cr,nrnb,fr,top,shake_vir);
1398 /* Construct the virtual sites for the initial configuration */
1399 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
1400 top->idef.iparams,top->idef.il,
1401 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1407 /* I'm assuming we need global communication the first time! MRS */
1408 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
1409 | (bVV ? CGLO_PRESSURE:0)
1410 | (bVV ? CGLO_CONSTRAINT:0)
1411 | (bRerunMD ? CGLO_RERUNMD:0)
1412 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
1414 bSumEkinhOld = FALSE;
1415 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1416 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1417 constr,NULL,FALSE,state->box,
1418 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
1419 if (ir->eI == eiVVAK) {
1420 /* a second call to get the half step temperature initialized as well */
1421 /* we do the same call as above, but turn the pressure off -- internally to
1422 compute_globals, this is recognized as a velocity verlet half-step
1423 kinetic energy calculation. This minimized excess variables, but
1424 perhaps loses some logic?*/
1426 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1427 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1428 constr,NULL,FALSE,state->box,
1429 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1430 cglo_flags &~ CGLO_PRESSURE);
1433 /* Calculate the initial half step temperature, and save the ekinh_old */
1434 if (!(Flags & MD_STARTFROMCPT))
1436 for(i=0; (i<ir->opts.ngtc); i++)
1438 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
1443 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
1444 and there is no previous step */
1446 temp0 = enerd->term[F_TEMP];
1448 /* if using an iterative algorithm, we need to create a working directory for the state. */
1451 bufstate = init_bufstate(state);
1455 snew(xcopy,state->natoms);
1456 snew(vcopy,state->natoms);
1457 copy_rvecn(state->x,xcopy,0,state->natoms);
1458 copy_rvecn(state->v,vcopy,0,state->natoms);
1459 copy_mat(state->box,boxcopy);
1462 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
1463 temperature control */
1464 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
1468 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
1471 "RMS relative constraint deviation after constraining: %.2e\n",
1472 constr_rmsd(constr,FALSE));
1474 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
1477 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
1478 " input trajectory '%s'\n\n",
1479 *(top_global->name),opt2fn("-rerun",nfile,fnm));
1482 fprintf(stderr,"Calculated time to finish depends on nsteps from "
1483 "run input file,\nwhich may not correspond to the time "
1484 "needed to process input trajectory.\n\n");
1490 fprintf(stderr,"starting mdrun '%s'\n",
1491 *(top_global->name));
1492 if (ir->nsteps >= 0)
1494 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
1498 sprintf(tbuf,"%s","infinite");
1500 if (ir->init_step > 0)
1502 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
1503 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
1504 gmx_step_str(ir->init_step,sbuf2),
1505 ir->init_step*ir->delta_t);
1509 fprintf(stderr,"%s steps, %s ps.\n",
1510 gmx_step_str(ir->nsteps,sbuf),tbuf);
1513 fprintf(fplog,"\n");
1516 /* Set and write start time */
1517 runtime_start(runtime);
1518 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
1519 wallcycle_start(wcycle,ewcRUN);
1521 fprintf(fplog,"\n");
1523 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1525 chkpt_ret=fcCheckPointParallel( cr->nodeid,
1527 if ( chkpt_ret == 0 )
1528 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
1532 /***********************************************************
1534 * Loop over MD steps
1536 ************************************************************/
1538 /* if rerunMD then read coordinates and velocities from input trajectory */
1541 if (getenv("GMX_FORCE_UPDATE"))
1543 bForceUpdate = TRUE;
1546 rerun_fr.natoms = 0;
1549 bNotLastFrame = read_first_frame(oenv,&status,
1550 opt2fn("-rerun",nfile,fnm),
1551 &rerun_fr,TRX_NEED_X | TRX_READ_V);
1552 if (rerun_fr.natoms != top_global->natoms)
1555 "Number of atoms in trajectory (%d) does not match the "
1556 "run input file (%d)\n",
1557 rerun_fr.natoms,top_global->natoms);
1559 if (ir->ePBC != epbcNONE)
1563 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);
1565 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
1567 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
1574 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
1577 if (ir->ePBC != epbcNONE)
1579 /* Set the shift vectors.
1580 * Necessary here when have a static box different from the tpr box.
1582 calc_shifts(rerun_fr.box,fr->shift_vec);
1586 /* loop over MD steps or if rerunMD to end of input trajectory */
1588 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1589 bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
1590 bInitStep = bFirstStep && (bStateFromTPX || bVV);
1591 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
1593 bSumEkinhOld = FALSE;
1596 init_global_signals(&gs,cr,ir,repl_ex_nst);
1598 step = ir->init_step;
1601 if (ir->nstlist == -1)
1603 init_nlistheuristics(&nlh,bGStatEveryStep,step);
1606 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
1607 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
1609 wallcycle_start(wcycle,ewcSTEP);
1611 GMX_MPE_LOG(ev_timestep1);
1614 if (rerun_fr.bStep) {
1615 step = rerun_fr.step;
1616 step_rel = step - ir->init_step;
1618 if (rerun_fr.bTime) {
1628 bLastStep = (step_rel == ir->nsteps);
1629 t = t0 + step*ir->delta_t;
1632 if (ir->efep != efepNO)
1634 if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
1636 state_global->lambda = rerun_fr.lambda;
1640 state_global->lambda = lam0 + step*ir->delta_lambda;
1642 state->lambda = state_global->lambda;
1643 bDoDHDL = do_per_step(step,ir->nstdhdl);
1648 update_annealing_target_temp(&(ir->opts),t);
1653 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
1655 for(i=0; i<state_global->natoms; i++)
1657 copy_rvec(rerun_fr.x[i],state_global->x[i]);
1661 for(i=0; i<state_global->natoms; i++)
1663 copy_rvec(rerun_fr.v[i],state_global->v[i]);
1668 for(i=0; i<state_global->natoms; i++)
1670 clear_rvec(state_global->v[i]);
1674 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
1675 " Ekin, temperature and pressure are incorrect,\n"
1676 " the virial will be incorrect when constraints are present.\n"
1678 bRerunWarnNoV = FALSE;
1682 copy_mat(rerun_fr.box,state_global->box);
1683 copy_mat(state_global->box,state->box);
1685 if (vsite && (Flags & MD_RERUN_VSITE))
1687 if (DOMAINDECOMP(cr))
1689 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1693 /* Following is necessary because the graph may get out of sync
1694 * with the coordinates if we only have every N'th coordinate set
1696 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
1697 shift_self(graph,state->box,state->x);
1699 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1700 top->idef.iparams,top->idef.il,
1701 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1704 unshift_self(graph,state->box,state->x);
1709 /* Stop Center of Mass motion */
1710 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
1712 /* Copy back starting coordinates in case we're doing a forcefield scan */
1715 for(ii=0; (ii<state->natoms); ii++)
1717 copy_rvec(xcopy[ii],state->x[ii]);
1718 copy_rvec(vcopy[ii],state->v[ii]);
1720 copy_mat(boxcopy,state->box);
1725 /* for rerun MD always do Neighbour Searching */
1726 bNS = (bFirstStep || ir->nstlist != 0);
1731 /* Determine whether or not to do Neighbour Searching and LR */
1732 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
1734 bNS = (bFirstStep || bExchanged || bNStList ||
1735 (ir->nstlist == -1 && nlh.nabnsb > 0));
1737 if (bNS && ir->nstlist == -1)
1739 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
1743 /* < 0 means stop at next step, > 0 means stop at next NS step */
1744 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
1745 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
1750 /* Determine whether or not to update the Born radii if doing GB */
1751 bBornRadii=bFirstStep;
1752 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
1757 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
1758 do_verbose = bVerbose &&
1759 (step % stepout == 0 || bFirstStep || bLastStep);
1761 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
1765 bMasterState = TRUE;
1769 bMasterState = FALSE;
1770 /* Correct the new box if it is too skewed */
1771 if (DYNAMIC_BOX(*ir))
1773 if (correct_box(fplog,step,state->box,graph))
1775 bMasterState = TRUE;
1778 if (DOMAINDECOMP(cr) && bMasterState)
1780 dd_collect_state(cr->dd,state,state_global);
1784 if (DOMAINDECOMP(cr))
1786 /* Repartition the domain decomposition */
1787 wallcycle_start(wcycle,ewcDOMDEC);
1788 dd_partition_system(fplog,step,cr,
1789 bMasterState,nstglobalcomm,
1790 state_global,top_global,ir,
1791 state,&f,mdatoms,top,fr,
1792 vsite,shellfc,constr,
1793 nrnb,wcycle,do_verbose);
1794 wallcycle_stop(wcycle,ewcDOMDEC);
1795 /* If using an iterative integrator, reallocate space to match the decomposition */
1799 if (MASTER(cr) && do_log && !bFFscan)
1801 print_ebin_header(fplog,step,t,state->lambda);
1804 if (ir->efep != efepNO)
1806 update_mdatoms(mdatoms,state->lambda);
1809 if (bRerunMD && rerun_fr.bV)
1812 /* We need the kinetic energy at minus the half step for determining
1813 * the full step kinetic energy and possibly for T-coupling.*/
1814 /* This may not be quite working correctly yet . . . . */
1815 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1816 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
1817 constr,NULL,FALSE,state->box,
1818 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1819 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1821 clear_mat(force_vir);
1823 /* Ionize the atoms if necessary */
1826 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
1827 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
1830 /* Update force field in ffscan program */
1833 if (update_forcefield(fplog,
1835 mdatoms->nr,state->x,state->box)) {
1836 if (gmx_parallel_env_initialized())
1844 GMX_MPE_LOG(ev_timestep2);
1846 /* We write a checkpoint at this MD step when:
1847 * either at an NS step when we signalled through gs,
1848 * or at the last step (but not when we do not want confout),
1849 * but never at the first step or with rerun.
1851 bCPT = (((gs.set[eglsCHKPT] && bNS) ||
1852 (bLastStep && (Flags & MD_CONFOUT))) &&
1853 step > ir->init_step && !bRerunMD);
1856 gs.set[eglsCHKPT] = 0;
1859 /* Determine the energy and pressure:
1860 * at nstcalcenergy steps and at energy output steps (set below).
1862 bNstEner = do_per_step(step,ir->nstcalcenergy);
1865 (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)));
1867 /* Do we need global communication ? */
1868 bGStat = (bCalcEnerPres || bStopCM ||
1869 do_per_step(step,nstglobalcomm) ||
1870 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1872 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1874 if (do_ene || do_log)
1876 bCalcEnerPres = TRUE;
1880 /* these CGLO_ options remain the same throughout the iteration */
1881 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1882 (bStopCM ? CGLO_STOPCM : 0) |
1883 (bGStat ? CGLO_GSTAT : 0)
1886 force_flags = (GMX_FORCE_STATECHANGED |
1887 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1888 GMX_FORCE_ALLFORCES |
1889 (bNStList ? GMX_FORCE_DOLR : 0) |
1891 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
1892 (bDoDHDL ? GMX_FORCE_DHDL : 0)
1897 /* Now is the time to relax the shells */
1898 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
1900 bStopCM,top,top_global,
1902 state,f,force_vir,mdatoms,
1903 nrnb,wcycle,graph,groups,
1904 shellfc,fr,bBornRadii,t,mu_tot,
1905 state->natoms,&bConverged,vsite,
1916 /* The coordinates (x) are shifted (to get whole molecules)
1918 * This is parallellized as well, and does communication too.
1919 * Check comments in sim_util.c
1922 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1923 state->box,state->x,&state->hist,
1924 f,force_vir,mdatoms,enerd,fcd,
1925 state->lambda,graph,
1926 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1927 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1930 GMX_BARRIER(cr->mpi_comm_mygroup);
1934 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1935 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1938 if (bTCR && bFirstStep)
1940 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1941 fprintf(fplog,"Done init_coupling\n");
1945 if (bVV && !bStartingFromCpt && !bRerunMD)
1946 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1948 if (ir->eI==eiVV && bInitStep)
1950 /* if using velocity verlet with full time step Ekin,
1951 * take the first half step only to compute the
1952 * virial for the first step. From there,
1953 * revert back to the initial coordinates
1954 * so that the input is actually the initial step.
1956 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1958 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1959 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1962 update_coords(fplog,step,ir,mdatoms,state,
1963 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1964 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1965 cr,nrnb,constr,&top->idef);
1969 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1971 /* for iterations, we save these vectors, as we will be self-consistently iterating
1974 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1976 /* save the state */
1977 if (bIterations && iterate.bIterate) {
1978 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1981 bFirstIterate = TRUE;
1982 while (bFirstIterate || (bIterations && iterate.bIterate))
1984 if (bIterations && iterate.bIterate)
1986 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1987 if (bFirstIterate && bTrotter)
1989 /* The first time through, we need a decent first estimate
1990 of veta(t+dt) to compute the constraints. Do
1991 this by computing the box volume part of the
1992 trotter integration at this time. Nothing else
1993 should be changed by this routine here. If
1994 !(first time), we start with the previous value
1997 veta_save = state->veta;
1998 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
1999 vetanew = state->veta;
2000 state->veta = veta_save;
2005 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
2008 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2009 &top->idef,shake_vir,NULL,
2010 cr,nrnb,wcycle,upd,constr,
2011 bInitStep,TRUE,bCalcEnerPres,vetanew);
2013 if (!bOK && !bFFscan)
2015 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2020 { /* Need to unshift here if a do_force has been
2021 called in the previous step */
2022 unshift_self(graph,state->box,state->x);
2026 /* if VV, compute the pressure and constraints */
2027 /* For VV2, we strictly only need this if using pressure
2028 * control, but we really would like to have accurate pressures
2030 * Think about ways around this in the future?
2031 * For now, keep this choice in comments.
2033 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
2034 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
2036 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
2037 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2038 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2039 constr,NULL,FALSE,state->box,
2040 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2043 | (bTemp ? CGLO_TEMPERATURE:0)
2044 | (bPres ? CGLO_PRESSURE : 0)
2045 | (bPres ? CGLO_CONSTRAINT : 0)
2046 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
2047 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2050 /* explanation of above:
2051 a) We compute Ekin at the full time step
2052 if 1) we are using the AveVel Ekin, and it's not the
2053 initial step, or 2) if we are using AveEkin, but need the full
2054 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
2055 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2056 EkinAveVel because it's needed for the pressure */
2058 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2063 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
2067 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2072 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2073 state->veta,&vetanew))
2077 bFirstIterate = FALSE;
2080 if (bTrotter && !bInitStep) {
2081 copy_mat(shake_vir,state->svir_prev);
2082 copy_mat(force_vir,state->fvir_prev);
2083 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
2084 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2085 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
2086 enerd->term[F_EKIN] = trace(ekind->ekin);
2089 /* if it's the initial step, we performed this first step just to get the constraint virial */
2090 if (bInitStep && ir->eI==eiVV) {
2091 copy_rvecn(cbuf,state->v,0,state->natoms);
2094 if (fr->bSepDVDL && fplog && do_log)
2096 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2098 enerd->term[F_DHDL_CON] += dvdl;
2100 GMX_MPE_LOG(ev_timestep1);
2103 /* MRS -- now done iterating -- compute the conserved quantity */
2105 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
2108 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
2110 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
2112 saved_conserved_quantity -= enerd->term[F_DISPCORR];
2116 /* ######## END FIRST UPDATE STEP ############## */
2117 /* ######## If doing VV, we now have v(dt) ###### */
2119 /* ################## START TRAJECTORY OUTPUT ################# */
2121 /* Now we have the energies and forces corresponding to the
2122 * coordinates at time t. We must output all of this before
2124 * for RerunMD t is read from input trajectory
2126 GMX_MPE_LOG(ev_output_start);
2129 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
2130 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
2131 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
2132 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
2133 if (bCPT) { mdof_flags |= MDOF_CPT; };
2135 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
2138 /* Enforce writing positions and velocities at end of run */
2139 mdof_flags |= (MDOF_X | MDOF_V);
2144 fcReportProgress( ir->nsteps, step );
2146 /* sync bCPT and fc record-keeping */
2147 if (bCPT && MASTER(cr))
2148 fcRequestCheckPoint();
2151 if (mdof_flags != 0)
2153 wallcycle_start(wcycle,ewcTRAJ);
2156 if (state->flags & (1<<estLD_RNG))
2158 get_stochd_state(upd,state);
2164 state_global->ekinstate.bUpToDate = FALSE;
2168 update_ekinstate(&state_global->ekinstate,ekind);
2169 state_global->ekinstate.bUpToDate = TRUE;
2171 update_energyhistory(&state_global->enerhist,mdebin);
2174 write_traj(fplog,cr,outf,mdof_flags,top_global,
2175 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
2182 if (bLastStep && step_rel == ir->nsteps &&
2183 (Flags & MD_CONFOUT) && MASTER(cr) &&
2184 !bRerunMD && !bFFscan)
2186 /* x and v have been collected in write_traj,
2187 * because a checkpoint file will always be written
2190 fprintf(stderr,"\nWriting final coordinates.\n");
2191 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
2194 /* Make molecules whole only for confout writing */
2195 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
2197 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2198 *top_global->name,top_global,
2199 state_global->x,state_global->v,
2200 ir->ePBC,state->box);
2203 wallcycle_stop(wcycle,ewcTRAJ);
2205 GMX_MPE_LOG(ev_output_finish);
2207 /* kludge -- virial is lost with restart for NPT control. Must restart */
2208 if (bStartingFromCpt && bVV)
2210 copy_mat(state->svir_prev,shake_vir);
2211 copy_mat(state->fvir_prev,force_vir);
2213 /* ################## END TRAJECTORY OUTPUT ################ */
2215 /* Determine the pressure:
2216 * always when we want exact averages in the energy file,
2217 * at ns steps when we have pressure coupling,
2218 * otherwise only at energy output steps (set below).
2221 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2222 bCalcEnerPres = bNstEner;
2224 /* Do we need global communication ? */
2225 bGStat = (bGStatEveryStep || bStopCM || bNS ||
2226 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2228 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2230 if (do_ene || do_log)
2232 bCalcEnerPres = TRUE;
2236 /* Determine the wallclock run time up till now */
2237 run_time = gmx_gettime() - (double)runtime->real;
2239 /* Check whether everything is still allright */
2240 if (((int)gmx_get_stop_condition() > handled_stop_condition)
2246 /* this is just make gs.sig compatible with the hack
2247 of sending signals around by MPI_Reduce with together with
2249 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
2250 gs.sig[eglsSTOPCOND]=1;
2251 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
2252 gs.sig[eglsSTOPCOND]=-1;
2253 /* < 0 means stop at next step, > 0 means stop at next NS step */
2257 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2258 gmx_get_signal_name(),
2259 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2263 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2264 gmx_get_signal_name(),
2265 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2267 handled_stop_condition=(int)gmx_get_stop_condition();
2269 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
2270 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
2271 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
2273 /* Signal to terminate the run */
2274 gs.sig[eglsSTOPCOND] = 1;
2277 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2279 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2282 if (bResetCountersHalfMaxH && MASTER(cr) &&
2283 run_time > max_hours*60.0*60.0*0.495)
2285 gs.sig[eglsRESETCOUNTERS] = 1;
2288 if (ir->nstlist == -1 && !bRerunMD)
2290 /* When bGStatEveryStep=FALSE, global_stat is only called
2291 * when we check the atom displacements, not at NS steps.
2292 * This means that also the bonded interaction count check is not
2293 * performed immediately after NS. Therefore a few MD steps could
2294 * be performed with missing interactions.
2295 * But wrong energies are never written to file,
2296 * since energies are only written after global_stat
2299 if (step >= nlh.step_nscheck)
2301 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
2302 nlh.scale_tot,state->x);
2306 /* This is not necessarily true,
2307 * but step_nscheck is determined quite conservatively.
2313 /* In parallel we only have to check for checkpointing in steps
2314 * where we do global communication,
2315 * otherwise the other nodes don't know.
2317 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
2320 run_time >= nchkpt*cpt_period*60.0)) &&
2321 gs.set[eglsCHKPT] == 0)
2323 gs.sig[eglsCHKPT] = 1;
2328 gmx_iterate_init(&iterate,bIterations);
2331 /* for iterations, we save these vectors, as we will be redoing the calculations */
2332 if (bIterations && iterate.bIterate)
2334 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
2336 bFirstIterate = TRUE;
2337 while (bFirstIterate || (bIterations && iterate.bIterate))
2339 /* We now restore these vectors to redo the calculation with improved extended variables */
2342 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
2345 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2346 so scroll down for that logic */
2348 /* ######### START SECOND UPDATE STEP ################# */
2349 GMX_MPE_LOG(ev_update_start);
2350 /* Box is changed in update() when we do pressure coupling,
2351 * but we should still use the old box for energy corrections and when
2352 * writing it to the energy file, so it matches the trajectory files for
2353 * the same timestep above. Make a copy in a separate array.
2355 copy_mat(state->box,lastbox);
2358 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
2360 wallcycle_start(wcycle,ewcUPDATE);
2362 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
2365 if (bIterations && iterate.bIterate)
2373 /* we use a new value of scalevir to converge the iterations faster */
2374 scalevir = tracevir/trace(shake_vir);
2376 msmul(shake_vir,scalevir,shake_vir);
2377 m_add(force_vir,shake_vir,total_vir);
2378 clear_mat(shake_vir);
2380 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
2381 /* We can only do Berendsen coupling after we have summed
2382 * the kinetic energy or virial. Since the happens
2383 * in global_state after update, we should only do it at
2384 * step % nstlist = 1 with bGStatEveryStep=FALSE.
2389 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2390 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
2396 /* velocity half-step update */
2397 update_coords(fplog,step,ir,mdatoms,state,f,
2398 fr->bTwinRange && bNStList,fr->f_twin,fcd,
2399 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
2400 cr,nrnb,constr,&top->idef);
2403 /* Above, initialize just copies ekinh into ekin,
2404 * it doesn't copy position (for VV),
2405 * and entire integrator for MD.
2410 copy_rvecn(state->x,cbuf,0,state->natoms);
2413 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2414 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2415 wallcycle_stop(wcycle,ewcUPDATE);
2417 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2418 &top->idef,shake_vir,force_vir,
2419 cr,nrnb,wcycle,upd,constr,
2420 bInitStep,FALSE,bCalcEnerPres,state->veta);
2424 /* erase F_EKIN and F_TEMP here? */
2425 /* just compute the kinetic energy at the half step to perform a trotter step */
2426 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2427 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2428 constr,NULL,FALSE,lastbox,
2429 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2430 cglo_flags | CGLO_TEMPERATURE
2432 wallcycle_start(wcycle,ewcUPDATE);
2433 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
2434 /* now we know the scaling, we can compute the positions again again */
2435 copy_rvecn(cbuf,state->x,0,state->natoms);
2437 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2438 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2439 wallcycle_stop(wcycle,ewcUPDATE);
2441 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2442 /* are the small terms in the shake_vir here due
2443 * to numerical errors, or are they important
2444 * physically? I'm thinking they are just errors, but not completely sure.
2445 * For now, will call without actually constraining, constr=NULL*/
2446 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2447 &top->idef,tmp_vir,force_vir,
2448 cr,nrnb,wcycle,upd,NULL,
2449 bInitStep,FALSE,bCalcEnerPres,
2452 if (!bOK && !bFFscan)
2454 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2457 if (fr->bSepDVDL && fplog && do_log)
2459 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2461 enerd->term[F_DHDL_CON] += dvdl;
2465 /* Need to unshift here */
2466 unshift_self(graph,state->box,state->x);
2469 GMX_BARRIER(cr->mpi_comm_mygroup);
2470 GMX_MPE_LOG(ev_update_finish);
2474 wallcycle_start(wcycle,ewcVSITECONSTR);
2477 shift_self(graph,state->box,state->x);
2479 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
2480 top->idef.iparams,top->idef.il,
2481 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2485 unshift_self(graph,state->box,state->x);
2487 wallcycle_stop(wcycle,ewcVSITECONSTR);
2490 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
2491 if (ir->nstlist == -1 && bFirstIterate)
2493 gs.sig[eglsNABNSB] = nlh.nabnsb;
2495 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2496 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2498 bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
2500 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2502 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
2503 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
2504 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
2505 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
2506 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2509 if (ir->nstlist == -1 && bFirstIterate)
2511 nlh.nabnsb = gs.set[eglsNABNSB];
2512 gs.set[eglsNABNSB] = 0;
2514 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
2515 /* ############# END CALC EKIN AND PRESSURE ################# */
2517 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
2518 the virial that should probably be addressed eventually. state->veta has better properies,
2519 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
2520 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
2523 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2524 trace(shake_vir),&tracevir))
2528 bFirstIterate = FALSE;
2531 update_box(fplog,step,ir,mdatoms,state,graph,f,
2532 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
2534 /* ################# END UPDATE STEP 2 ################# */
2535 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2537 /* The coordinates (x) were unshifted in update */
2538 if (bFFscan && (shellfc==NULL || bConverged))
2540 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
2542 &(top_global->mols),mdatoms->massT,pres))
2544 if (gmx_parallel_env_initialized())
2548 fprintf(stderr,"\n");
2554 /* We will not sum ekinh_old,
2555 * so signal that we still have to do it.
2557 bSumEkinhOld = TRUE;
2562 /* Only do GCT when the relaxation of shells (minimization) has converged,
2563 * otherwise we might be coupling to bogus energies.
2564 * In parallel we must always do this, because the other sims might
2568 /* Since this is called with the new coordinates state->x, I assume
2569 * we want the new box state->box too. / EL 20040121
2571 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
2573 mdatoms,&(top->idef),mu_aver,
2574 top_global->mols.nr,cr,
2575 state->box,total_vir,pres,
2576 mu_tot,state->x,f,bConverged);
2580 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2582 /* sum up the foreign energy and dhdl terms */
2583 sum_dhdl(enerd,state->lambda,ir);
2585 /* use the directly determined last velocity, not actually the averaged half steps */
2586 if (bTrotter && ir->eI==eiVV)
2588 enerd->term[F_EKIN] = last_ekin;
2590 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
2594 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
2598 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
2600 /* Check for excessively large energies */
2604 real etot_max = 1e200;
2606 real etot_max = 1e30;
2608 if (fabs(enerd->term[F_ETOT]) > etot_max)
2610 fprintf(stderr,"Energy too large (%g), giving up\n",
2611 enerd->term[F_ETOT]);
2614 /* ######### END PREPARING EDR OUTPUT ########### */
2616 /* Time for performance */
2617 if (((step % stepout) == 0) || bLastStep)
2619 runtime_upd_proc(runtime);
2625 gmx_bool do_dr,do_or;
2627 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
2631 upd_mdebin(mdebin,bDoDHDL, TRUE,
2632 t,mdatoms->tmass,enerd,state,lastbox,
2633 shake_vir,force_vir,total_vir,pres,
2634 ekind,mu_tot,constr);
2638 upd_mdebin_step(mdebin);
2641 do_dr = do_per_step(step,ir->nstdisreout);
2642 do_or = do_per_step(step,ir->nstorireout);
2644 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
2646 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
2648 if (ir->ePull != epullNO)
2650 pull_print_output(ir->pull,step,t);
2653 if (do_per_step(step,ir->nstlog))
2655 if(fflush(fplog) != 0)
2657 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
2663 /* Remaining runtime */
2664 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
2668 fprintf(stderr,"\n");
2670 print_time(stderr,runtime,step,ir,cr);
2673 /* Replica exchange */
2675 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
2676 do_per_step(step,repl_ex_nst))
2678 bExchanged = replica_exchange(fplog,cr,repl_ex,
2679 state_global,enerd->term,
2682 if (bExchanged && DOMAINDECOMP(cr))
2684 dd_partition_system(fplog,step,cr,TRUE,1,
2685 state_global,top_global,ir,
2686 state,&f,mdatoms,top,fr,
2687 vsite,shellfc,constr,
2694 bStartingFromCpt = FALSE;
2696 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2697 /* Complicated conditional when bGStatEveryStep=FALSE.
2698 * We can not just use bGStat, since then the simulation results
2699 * would depend on nstenergy and nstlog or step_nscheck.
2701 if (((state->flags & (1<<estPRES_PREV)) ||
2702 (state->flags & (1<<estSVIR_PREV)) ||
2703 (state->flags & (1<<estFVIR_PREV))) &&
2705 (ir->nstlist > 0 && step % ir->nstlist == 0) ||
2706 (ir->nstlist < 0 && nlh.nabnsb > 0) ||
2707 (ir->nstlist == 0 && bGStat)))
2709 /* Store the pressure in t_state for pressure coupling
2710 * at the next MD step.
2712 if (state->flags & (1<<estPRES_PREV))
2714 copy_mat(pres,state->pres_prev);
2718 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2724 /* read next frame from input trajectory */
2725 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
2730 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
2734 if (!bRerunMD || !rerun_fr.bStep)
2736 /* increase the MD step number */
2741 cycles = wallcycle_stop(wcycle,ewcSTEP);
2742 if (DOMAINDECOMP(cr) && wcycle)
2744 dd_cycles_add(cr->dd,cycles,ddCyclStep);
2747 if (step_rel == wcycle_get_reset_counters(wcycle) ||
2748 gs.set[eglsRESETCOUNTERS] != 0)
2750 /* Reset all the counters related to performance over the run */
2751 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
2752 wcycle_set_reset_counters(wcycle,-1);
2753 /* Correct max_hours for the elapsed time */
2754 max_hours -= run_time/(60.0*60.0);
2755 bResetCountersHalfMaxH = FALSE;
2756 gs.set[eglsRESETCOUNTERS] = 0;
2759 /* End of main MD loop */
2763 runtime_end(runtime);
2765 if (bRerunMD && MASTER(cr))
2770 if (!(cr->duty & DUTY_PME))
2772 /* Tell the PME only node to finish */
2778 if (ir->nstcalcenergy > 0 && !bRerunMD)
2780 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
2781 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
2789 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
2791 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)));
2792 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
2795 if (shellfc && fplog)
2797 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
2798 (nconverged*100.0)/step_rel);
2799 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
2803 if (repl_ex_nst > 0 && MASTER(cr))
2805 print_replica_exchange_statistics(fplog,repl_ex);
2808 runtime->nsteps_done = step_rel;