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 /* FIX ME after 4.5 */
327 /* temporary hack because we are using gmx_bool (unsigned char) */
329 bReadEkin = (flags & CGLO_READEKIN) != 0;
330 bScaleEkin = (flags & CGLO_SCALEEKIN) != 0;
331 bEner = flags & CGLO_ENERGY;
332 bTemp = flags & CGLO_TEMPERATURE;
333 bPres = (flags & CGLO_PRESSURE) != 0;
334 bConstrain = (flags & CGLO_CONSTRAINT) != 0;
335 bIterate = (flags & CGLO_ITERATE) != 0;
336 bFirstIterate = (flags & CGLO_FIRSTITERATE) != 0;
338 /* we calculate a full state kinetic energy either with full-step velocity verlet
339 or half step where we need the pressure */
341 bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && bPres) || bReadEkin);
343 /* in initalization, it sums the shake virial in vv, and to
344 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
346 /* ########## Kinetic energy ############## */
350 /* Non-equilibrium MD: this is parallellized, but only does communication
351 * when there really is NEMD.
354 if (PAR(cr) && (ekind->bNEMD))
356 accumulate_u(cr,&(ir->opts),ekind);
361 restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
366 calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
371 /* Calculate center of mass velocity if necessary, also parallellized */
372 if (bStopCM && !bRerunMD && bEner)
374 calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
375 state->x,state->v,vcm);
379 if (bTemp || bPres || bEner || bConstrain)
383 /* We will not sum ekinh_old,
384 * so signal that we still have to do it.
386 *bSumEkinhOld = TRUE;
393 for(i=0; i<eglsNR; i++)
395 gs_buf[i] = gs->sig[i];
400 wallcycle_start(wcycle,ewcMoveE);
401 GMX_MPE_LOG(ev_global_stat_start);
402 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
404 gs != NULL ? eglsNR : 0,gs_buf,
406 *bSumEkinhOld,flags);
407 GMX_MPE_LOG(ev_global_stat_finish);
408 wallcycle_stop(wcycle,ewcMoveE);
412 if (MULTISIM(cr) && bInterSimGS)
416 /* Communicate the signals between the simulations */
417 gmx_sum_sim(eglsNR,gs_buf,cr->ms);
419 /* Communicate the signals form the master to the others */
420 gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
422 for(i=0; i<eglsNR; i++)
424 if (bInterSimGS || gs_simlocal[i])
426 /* Set the communicated signal only when it is non-zero,
427 * since signals might not be processed at each MD step.
429 gsi = (gs_buf[i] >= 0 ?
430 (int)(gs_buf[i] + 0.5) :
431 (int)(gs_buf[i] - 0.5));
436 /* Turn off the local signal */
441 *bSumEkinhOld = FALSE;
445 if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
448 mdatoms->start,mdatoms->start+mdatoms->homenr,
449 state->v,vcm->group_p[0],
450 mdatoms->massT,mdatoms->tmass,ekind->ekin);
454 /* Do center of mass motion removal */
455 if (bStopCM && !bRerunMD) /* is this correct? Does it get called too often with this logic? */
457 check_cm_grp(fplog,vcm,ir,1);
458 do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
459 state->x,state->v,vcm);
460 inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
466 /* Sum the kinetic energies of the groups & calc temp */
467 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
468 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
469 Leap with AveVel is not supported; it's not clear that it will actually work.
470 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
471 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
472 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
473 If FALSE, we go ahead and erase over it.
475 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,&(enerd->term[F_DKDL]),
476 bEkinAveVel,bIterate,bScaleEkin);
478 enerd->term[F_EKIN] = trace(ekind->ekin);
481 /* ########## Long range energy information ###### */
483 if (bEner || bPres || bConstrain)
485 calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda,
486 corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
489 if (bEner && bFirstIterate)
491 enerd->term[F_DISPCORR] = enercorr;
492 enerd->term[F_EPOT] += enercorr;
493 enerd->term[F_DVDL] += dvdlcorr;
494 if (fr->efep != efepNO) {
495 enerd->dvdl_lin += dvdlcorr;
499 /* ########## Now pressure ############## */
500 if (bPres || bConstrain)
503 m_add(force_vir,shake_vir,total_vir);
505 /* Calculate pressure and apply LR correction if PPPM is used.
506 * Use the box from last timestep since we already called update().
509 enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres,
510 (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
512 /* Calculate long range corrections to pressure and energy */
513 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
514 and computes enerd->term[F_DISPCORR]. Also modifies the
515 total_vir and pres tesors */
517 m_add(total_vir,corr_vir,total_vir);
518 m_add(pres,corr_pres,pres);
519 enerd->term[F_PDISPCORR] = prescorr;
520 enerd->term[F_PRES] += prescorr;
521 *pcurr = enerd->term[F_PRES];
522 /* calculate temperature using virial */
523 enerd->term[F_VTEMP] = calc_temp(trace(total_vir),ir->opts.nrdf[0]);
529 /* Definitions for convergence of iterated constraints */
531 /* iterate constraints up to 50 times */
532 #define MAXITERCONST 50
537 real f,fprev,x,xprev;
540 real allrelerr[MAXITERCONST+2];
541 int num_close; /* number of "close" violations, caused by limited precision. */
545 #define CONVERGEITER 0.000000001
546 #define CLOSE_ENOUGH 0.000001000
548 #define CONVERGEITER 0.0001
549 #define CLOSE_ENOUGH 0.0050
552 /* we want to keep track of the close calls. If there are too many, there might be some other issues.
553 so we make sure that it's either less than some predetermined number, or if more than that number,
554 only some small fraction of the total. */
555 #define MAX_NUMBER_CLOSE 50
556 #define FRACTION_CLOSE 0.001
558 /* maximum length of cyclic traps to check, emerging from limited numerical precision */
561 static void gmx_iterate_init(gmx_iterate_t *iterate,gmx_bool bIterate)
566 iterate->bIterate = bIterate;
567 iterate->num_close = 0;
568 for (i=0;i<MAXITERCONST+2;i++)
570 iterate->allrelerr[i] = 0;
574 static gmx_bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate_t *iterate, gmx_bool bFirstIterate, real fom, real *newf)
576 /* monitor convergence, and use a secant search to propose new
579 The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
580 f(x_{i}) - f(x_{i-1})
582 The function we are trying to zero is fom-x, where fom is the
583 "figure of merit" which is the pressure (or the veta value) we
584 would get by putting in an old value of the pressure or veta into
585 the incrementor function for the step or half step. I have
586 verified that this gives the same answer as self consistent
587 iteration, usually in many fewer steps, especially for small tau_p.
589 We could possibly eliminate an iteration with proper use
590 of the value from the previous step, but that would take a bit
591 more bookkeeping, especially for veta, since tests indicate the
592 function of veta on the last step is not sufficiently close to
593 guarantee convergence this step. This is
594 good enough for now. On my tests, I could use tau_p down to
595 0.02, which is smaller that would ever be necessary in
596 practice. Generally, 3-5 iterations will be sufficient */
598 real relerr,err,xmin;
606 iterate->f = fom-iterate->x;
613 iterate->f = fom-iterate->x; /* we want to zero this difference */
614 if ((iterate->iter_i > 1) && (iterate->iter_i < MAXITERCONST))
616 if (iterate->f==iterate->fprev)
622 *newf = iterate->x - (iterate->x-iterate->xprev)*(iterate->f)/(iterate->f-iterate->fprev);
627 /* just use self-consistent iteration the first step to initialize, or
628 if it's not converging (which happens occasionally -- need to investigate why) */
632 /* Consider a slight shortcut allowing us to exit one sooner -- we check the
633 difference between the closest of x and xprev to the new
634 value. To be 100% certain, we should check the difference between
635 the last result, and the previous result, or
637 relerr = (fabs((x-xprev)/fom));
639 but this is pretty much never necessary under typical conditions.
640 Checking numerically, it seems to lead to almost exactly the same
641 trajectories, but there are small differences out a few decimal
642 places in the pressure, and eventually in the v_eta, but it could
645 if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
646 relerr = (fabs((*newf-xmin) / *newf));
649 err = fabs((iterate->f-iterate->fprev));
650 relerr = fabs(err/fom);
652 iterate->allrelerr[iterate->iter_i] = relerr;
654 if (iterate->iter_i > 0)
658 fprintf(debug,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
659 iterate->iter_i,fom,relerr,*newf);
662 if ((relerr < CONVERGEITER) || (err < CONVERGEITER) || (fom==0) || ((iterate->x == iterate->xprev) && iterate->iter_i > 1))
664 iterate->bIterate = FALSE;
667 fprintf(debug,"Iterating NPT constraints: CONVERGED\n");
671 if (iterate->iter_i > MAXITERCONST)
673 if (relerr < CLOSE_ENOUGH)
676 for (i=1;i<CYCLEMAX;i++) {
677 if ((iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-1]) &&
678 (iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-(1+2*i)])) {
682 fprintf(debug,"Exiting from an NPT iterating cycle of length %d\n",i);
689 /* step 1: trapped in a numerical attractor */
690 /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
691 Better to give up convergence here than have the simulation die.
693 iterate->num_close++;
698 /* Step #2: test if we are reasonably close for other reasons, then monitor the number. If not, die */
700 /* how many close calls have we had? If less than a few, we're OK */
701 if (iterate->num_close < MAX_NUMBER_CLOSE)
703 sprintf(buf,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps,relerr);
704 md_print_warning(cr,fplog,buf);
705 iterate->num_close++;
707 /* if more than a few, check the total fraction. If too high, die. */
708 } else if (iterate->num_close/(double)nsteps > FRACTION_CLOSE) {
709 gmx_fatal(FARGS,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate->num_close/(double)nsteps);
715 gmx_fatal(FARGS,"Could not converge NPT constraints\n");
720 iterate->xprev = iterate->x;
722 iterate->fprev = iterate->f;
728 static void check_nst_param(FILE *fplog,t_commrec *cr,
729 const char *desc_nst,int nst,
730 const char *desc_p,int *p)
734 if (*p > 0 && *p % nst != 0)
736 /* Round up to the next multiple of nst */
737 *p = ((*p)/nst + 1)*nst;
738 sprintf(buf,"NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
739 md_print_warning(cr,fplog,buf);
743 static void reset_all_counters(FILE *fplog,t_commrec *cr,
744 gmx_large_int_t step,
745 gmx_large_int_t *step_rel,t_inputrec *ir,
746 gmx_wallcycle_t wcycle,t_nrnb *nrnb,
747 gmx_runtime_t *runtime)
749 char buf[STRLEN],sbuf[STEPSTRSIZE];
751 /* Reset all the counters related to performance over the run */
752 sprintf(buf,"Step %s: resetting all time and cycle counters\n",
753 gmx_step_str(step,sbuf));
754 md_print_warning(cr,fplog,buf);
756 wallcycle_stop(wcycle,ewcRUN);
757 wallcycle_reset_all(wcycle);
758 if (DOMAINDECOMP(cr))
760 reset_dd_statistics_counters(cr->dd);
763 ir->init_step += *step_rel;
764 ir->nsteps -= *step_rel;
766 wallcycle_start(wcycle,ewcRUN);
767 runtime_start(runtime);
768 print_date_and_time(fplog,cr->nodeid,"Restarted time",runtime);
771 static void min_zero(int *n,int i)
773 if (i > 0 && (*n == 0 || i < *n))
779 static int lcd4(int i1,int i2,int i3,int i4)
790 gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
793 while (nst > 1 && ((i1 > 0 && i1 % nst != 0) ||
794 (i2 > 0 && i2 % nst != 0) ||
795 (i3 > 0 && i3 % nst != 0) ||
796 (i4 > 0 && i4 % nst != 0)))
804 static int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
805 int nstglobalcomm,t_inputrec *ir)
809 if (!EI_DYNAMICS(ir->eI))
814 if (nstglobalcomm == -1)
816 if (!(ir->nstcalcenergy > 0 ||
822 if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
824 nstglobalcomm = ir->nstenergy;
829 /* Ensure that we do timely global communication for
830 * (possibly) each of the four following options.
832 nstglobalcomm = lcd4(ir->nstcalcenergy,
834 ir->etc != etcNO ? ir->nsttcouple : 0,
835 ir->epc != epcNO ? ir->nstpcouple : 0);
840 if (ir->nstlist > 0 &&
841 nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
843 nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
844 sprintf(buf,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
845 md_print_warning(cr,fplog,buf);
847 if (ir->nstcalcenergy > 0)
849 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
850 "nstcalcenergy",&ir->nstcalcenergy);
852 if (ir->etc != etcNO && ir->nsttcouple > 0)
854 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
855 "nsttcouple",&ir->nsttcouple);
857 if (ir->epc != epcNO && ir->nstpcouple > 0)
859 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
860 "nstpcouple",&ir->nstpcouple);
863 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
864 "nstenergy",&ir->nstenergy);
866 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
867 "nstlog",&ir->nstlog);
870 if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
872 sprintf(buf,"WARNING: Changing nstcomm from %d to %d\n",
873 ir->nstcomm,nstglobalcomm);
874 md_print_warning(cr,fplog,buf);
875 ir->nstcomm = nstglobalcomm;
878 return nstglobalcomm;
881 void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
882 t_inputrec *ir,gmx_mtop_t *mtop)
884 /* Check required for old tpx files */
885 if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
886 ir->nstcalcenergy % ir->nstlist != 0)
888 md_print_warning(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
890 if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
891 gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
892 ir->eConstrAlg == econtSHAKE)
894 md_print_warning(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
895 if (ir->epc != epcNO)
897 gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
900 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
901 "nstcalcenergy",&ir->nstcalcenergy);
902 if (ir->epc != epcNO)
904 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
905 "nstpcouple",&ir->nstpcouple);
907 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
908 "nstenergy",&ir->nstenergy);
909 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
910 "nstlog",&ir->nstlog);
911 if (ir->efep != efepNO)
913 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
914 "nstdhdl",&ir->nstdhdl);
920 gmx_bool bGStatEveryStep;
921 gmx_large_int_t step_ns;
922 gmx_large_int_t step_nscheck;
933 static void reset_nlistheuristics(gmx_nlheur_t *nlh,gmx_large_int_t step)
937 nlh->step_nscheck = step;
940 static void init_nlistheuristics(gmx_nlheur_t *nlh,
941 gmx_bool bGStatEveryStep,gmx_large_int_t step)
943 nlh->bGStatEveryStep = bGStatEveryStep;
950 reset_nlistheuristics(nlh,step);
953 static void update_nliststatistics(gmx_nlheur_t *nlh,gmx_large_int_t step)
955 gmx_large_int_t nl_lt;
956 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
958 /* Determine the neighbor list life time */
959 nl_lt = step - nlh->step_ns;
962 fprintf(debug,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh->nabnsb,gmx_step_str(nl_lt,sbuf));
966 nlh->s2 += nl_lt*nl_lt;
967 nlh->ab += nlh->nabnsb;
968 if (nlh->lt_runav == 0)
970 nlh->lt_runav = nl_lt;
971 /* Initialize the fluctuation average
972 * such that at startup we check after 0 steps.
974 nlh->lt_runav2 = sqr(nl_lt/2.0);
976 /* Running average with 0.9 gives an exp. history of 9.5 */
977 nlh->lt_runav2 = 0.9*nlh->lt_runav2 + 0.1*sqr(nlh->lt_runav - nl_lt);
978 nlh->lt_runav = 0.9*nlh->lt_runav + 0.1*nl_lt;
979 if (nlh->bGStatEveryStep)
981 /* Always check the nlist validity */
982 nlh->step_nscheck = step;
986 /* We check after: <life time> - 2*sigma
987 * The factor 2 is quite conservative,
988 * but we assume that with nstlist=-1 the user
989 * prefers exact integration over performance.
991 nlh->step_nscheck = step
992 + (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)) - 1;
996 fprintf(debug,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
997 gmx_step_str(nl_lt,sbuf),nlh->lt_runav,sqrt(nlh->lt_runav2),
998 gmx_step_str(nlh->step_nscheck-step+1,sbuf2),
999 (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)));
1003 static void set_nlistheuristics(gmx_nlheur_t *nlh,gmx_bool bReset,gmx_large_int_t step)
1009 reset_nlistheuristics(nlh,step);
1013 update_nliststatistics(nlh,step);
1016 nlh->step_ns = step;
1017 /* Initialize the cumulative coordinate scaling matrix */
1018 clear_mat(nlh->scale_tot);
1019 for(d=0; d<DIM; d++)
1021 nlh->scale_tot[d][d] = 1.0;
1025 static void rerun_parallel_comm(t_commrec *cr,t_trxframe *fr,
1026 gmx_bool *bNotLastFrame)
1031 bAlloc = (fr->natoms == 0);
1033 if (MASTER(cr) && !*bNotLastFrame)
1039 gmx_bcast(sizeof(*fr),fr,cr);
1043 *bNotLastFrame = (fr->natoms >= 0);
1045 if (*bNotLastFrame && PARTDECOMP(cr))
1047 /* x and v are the only variable size quantities stored in trr
1048 * that are required for rerun (f is not needed).
1052 snew(fr->x,fr->natoms);
1053 snew(fr->v,fr->natoms);
1057 gmx_bcast(fr->natoms*sizeof(fr->x[0]),fr->x[0],cr);
1061 gmx_bcast(fr->natoms*sizeof(fr->v[0]),fr->v[0],cr);
1066 double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
1067 const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
1069 gmx_vsite_t *vsite,gmx_constr_t constr,
1070 int stepout,t_inputrec *ir,
1071 gmx_mtop_t *top_global,
1073 t_state *state_global,
1075 t_nrnb *nrnb,gmx_wallcycle_t wcycle,
1076 gmx_edsam_t ed,t_forcerec *fr,
1077 int repl_ex_nst,int repl_ex_seed,
1078 real cpt_period,real max_hours,
1079 const char *deviceOptions,
1080 unsigned long Flags,
1081 gmx_runtime_t *runtime)
1084 gmx_large_int_t step,step_rel;
1087 gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
1088 gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
1089 bFirstStep,bStateFromTPX,bInitStep,bLastStep,
1090 bBornRadii,bStartingFromCpt;
1091 gmx_bool bDoDHDL=FALSE;
1092 gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
1093 bForceUpdate=FALSE,bCPT;
1095 gmx_bool bMasterState;
1096 int force_flags,cglo_flags;
1097 tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
1099 t_trxstatus *status;
1102 t_state *bufstate=NULL;
1103 matrix *scale_tot,pcoupl_mu,M,ebox;
1105 t_trxframe rerun_fr;
1106 gmx_repl_ex_t repl_ex=NULL;
1109 gmx_localtop_t *top;
1110 t_mdebin *mdebin=NULL;
1111 t_state *state=NULL;
1112 rvec *f_global=NULL;
1115 gmx_enerdata_t *enerd;
1117 gmx_global_stat_t gstat;
1118 gmx_update_t upd=NULL;
1119 t_graph *graph=NULL;
1123 gmx_groups_t *groups;
1124 gmx_ekindata_t *ekind, *ekind_save;
1125 gmx_shellfc_t shellfc;
1126 int count,nconverged=0;
1129 gmx_bool bIonize=FALSE;
1130 gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
1132 gmx_bool bResetCountersHalfMaxH=FALSE;
1133 gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
1134 real temp0,mu_aver=0,dvdl;
1136 atom_id *grpindex=NULL;
1138 t_coupl_rec *tcr=NULL;
1139 rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
1140 matrix boxcopy={{0}},lastbox;
1142 real fom,oldfom,veta_save,pcurr,scalevir,tracevir;
1145 real saved_conserved_quantity = 0;
1150 char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
1151 int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
1152 gmx_iterate_t iterate;
1154 /* Temporary addition for FAHCORE checkpointing */
1158 /* Check for special mdrun options */
1159 bRerunMD = (Flags & MD_RERUN);
1160 bIonize = (Flags & MD_IONIZE);
1161 bFFscan = (Flags & MD_FFSCAN);
1162 bAppend = (Flags & MD_APPENDFILES);
1163 if (Flags & MD_RESETCOUNTERSHALFWAY)
1167 /* Signal to reset the counters half the simulation steps. */
1168 wcycle_set_reset_counters(wcycle,ir->nsteps/2);
1170 /* Signal to reset the counters halfway the simulation time. */
1171 bResetCountersHalfMaxH = (max_hours > 0);
1174 /* md-vv uses averaged full step velocities for T-control
1175 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
1176 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
1177 bVV = EI_VV(ir->eI);
1178 if (bVV) /* to store the initial velocities while computing virial */
1180 snew(cbuf,top_global->natoms);
1182 /* all the iteratative cases - only if there are constraints */
1183 bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
1184 bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
1188 /* Since we don't know if the frames read are related in any way,
1189 * rebuild the neighborlist at every step.
1192 ir->nstcalcenergy = 1;
1196 check_ir_old_tpx_versions(cr,fplog,ir,top_global);
1198 nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
1199 bGStatEveryStep = (nstglobalcomm == 1);
1201 if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
1204 "To reduce the energy communication with nstlist = -1\n"
1205 "the neighbor list validity should not be checked at every step,\n"
1206 "this means that exact integration is not guaranteed.\n"
1207 "The neighbor list validity is checked after:\n"
1208 " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
1209 "In most cases this will result in exact integration.\n"
1210 "This reduces the energy communication by a factor of 2 to 3.\n"
1211 "If you want less energy communication, set nstlist > 3.\n\n");
1214 if (bRerunMD || bFFscan)
1218 groups = &top_global->groups;
1220 /* Initial values */
1221 init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
1222 nrnb,top_global,&upd,
1223 nfile,fnm,&outf,&mdebin,
1224 force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
1226 clear_mat(total_vir);
1228 /* Energy terms and groups */
1230 init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
1231 if (DOMAINDECOMP(cr))
1237 snew(f,top_global->natoms);
1240 /* Kinetic energy data */
1242 init_ekindata(fplog,top_global,&(ir->opts),ekind);
1243 /* needed for iteration of constraints */
1245 init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
1246 /* Copy the cos acceleration to the groups struct */
1247 ekind->cosacc.cos_accel = ir->cos_accel;
1249 gstat = global_stat_init(ir);
1252 /* Check for polarizable models and flexible constraints */
1253 shellfc = init_shell_flexcon(fplog,
1254 top_global,n_flexible_constraints(constr),
1255 (ir->bContinuation ||
1256 (DOMAINDECOMP(cr) && !MASTER(cr))) ?
1257 NULL : state_global->x);
1262 tMPI_Thread_mutex_lock(&deform_init_box_mutex);
1264 set_deform_reference_box(upd,
1265 deform_init_init_step_tpx,
1266 deform_init_box_tpx);
1268 tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
1273 double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
1274 if ((io > 2000) && MASTER(cr))
1276 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
1280 if (DOMAINDECOMP(cr)) {
1281 top = dd_init_local_top(top_global);
1284 dd_init_local_state(cr->dd,state_global,state);
1286 if (DDMASTER(cr->dd) && ir->nstfout) {
1287 snew(f_global,state_global->natoms);
1291 /* Initialize the particle decomposition and split the topology */
1292 top = split_system(fplog,top_global,ir,cr);
1294 pd_cg_range(cr,&fr->cg0,&fr->hcg);
1295 pd_at_range(cr,&a0,&a1);
1297 top = gmx_mtop_generate_local_top(top_global,ir);
1300 a1 = top_global->natoms;
1303 state = partdec_init_local_state(cr,state_global);
1306 atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
1309 set_vsite_top(vsite,top,mdatoms,cr);
1312 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
1313 graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
1317 make_local_shells(cr,mdatoms,shellfc);
1320 if (ir->pull && PAR(cr)) {
1321 dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
1325 if (DOMAINDECOMP(cr))
1327 /* Distribute the charge groups over the nodes from the master node */
1328 dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
1329 state_global,top_global,ir,
1330 state,&f,mdatoms,top,fr,
1331 vsite,shellfc,constr,
1335 update_mdatoms(mdatoms,state->lambda);
1339 /* Update mdebin with energy history if appending to output files */
1340 if ( Flags & MD_APPENDFILES )
1342 restore_energyhistory_from_state(mdebin,&state_global->enerhist);
1344 /* Set the initial energy history in state to zero by updating once */
1345 update_energyhistory(&state_global->enerhist,mdebin);
1348 if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
1349 /* Set the random state if we read a checkpoint file */
1350 set_stochd_state(upd,state);
1353 /* Initialize constraints */
1355 if (!DOMAINDECOMP(cr))
1356 set_constraints(constr,top,ir,mdatoms,cr);
1359 /* Check whether we have to GCT stuff */
1360 bTCR = ftp2bSet(efGCT,nfile,fnm);
1363 fprintf(stderr,"Will do General Coupling Theory!\n");
1365 gnx = top_global->mols.nr;
1367 for(i=0; (i<gnx); i++) {
1372 if (repl_ex_nst > 0 && MASTER(cr))
1373 repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
1374 repl_ex_nst,repl_ex_seed);
1376 if (!ir->bContinuation && !bRerunMD)
1378 if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
1380 /* Set the velocities of frozen particles to zero */
1381 for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
1383 for(m=0; m<DIM; m++)
1385 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
1395 /* Constrain the initial coordinates and velocities */
1396 do_constrain_first(fplog,constr,ir,mdatoms,state,f,
1397 graph,cr,nrnb,fr,top,shake_vir);
1401 /* Construct the virtual sites for the initial configuration */
1402 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
1403 top->idef.iparams,top->idef.il,
1404 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1410 /* I'm assuming we need global communication the first time! MRS */
1411 cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
1412 | (bVV ? CGLO_PRESSURE:0)
1413 | (bVV ? CGLO_CONSTRAINT:0)
1414 | (bRerunMD ? CGLO_RERUNMD:0)
1415 | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
1417 bSumEkinhOld = FALSE;
1418 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1419 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1420 constr,NULL,FALSE,state->box,
1421 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
1422 if (ir->eI == eiVVAK) {
1423 /* a second call to get the half step temperature initialized as well */
1424 /* we do the same call as above, but turn the pressure off -- internally to
1425 compute_globals, this is recognized as a velocity verlet half-step
1426 kinetic energy calculation. This minimized excess variables, but
1427 perhaps loses some logic?*/
1429 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1430 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
1431 constr,NULL,FALSE,state->box,
1432 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1433 cglo_flags &~ CGLO_PRESSURE);
1436 /* Calculate the initial half step temperature, and save the ekinh_old */
1437 if (!(Flags & MD_STARTFROMCPT))
1439 for(i=0; (i<ir->opts.ngtc); i++)
1441 copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
1446 enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
1447 and there is no previous step */
1449 temp0 = enerd->term[F_TEMP];
1451 /* if using an iterative algorithm, we need to create a working directory for the state. */
1454 bufstate = init_bufstate(state);
1458 snew(xcopy,state->natoms);
1459 snew(vcopy,state->natoms);
1460 copy_rvecn(state->x,xcopy,0,state->natoms);
1461 copy_rvecn(state->v,vcopy,0,state->natoms);
1462 copy_mat(state->box,boxcopy);
1465 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
1466 temperature control */
1467 trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
1471 if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
1474 "RMS relative constraint deviation after constraining: %.2e\n",
1475 constr_rmsd(constr,FALSE));
1477 fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
1480 fprintf(stderr,"starting md rerun '%s', reading coordinates from"
1481 " input trajectory '%s'\n\n",
1482 *(top_global->name),opt2fn("-rerun",nfile,fnm));
1485 fprintf(stderr,"Calculated time to finish depends on nsteps from "
1486 "run input file,\nwhich may not correspond to the time "
1487 "needed to process input trajectory.\n\n");
1493 fprintf(stderr,"starting mdrun '%s'\n",
1494 *(top_global->name));
1495 if (ir->nsteps >= 0)
1497 sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
1501 sprintf(tbuf,"%s","infinite");
1503 if (ir->init_step > 0)
1505 fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
1506 gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
1507 gmx_step_str(ir->init_step,sbuf2),
1508 ir->init_step*ir->delta_t);
1512 fprintf(stderr,"%s steps, %s ps.\n",
1513 gmx_step_str(ir->nsteps,sbuf),tbuf);
1516 fprintf(fplog,"\n");
1519 /* Set and write start time */
1520 runtime_start(runtime);
1521 print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
1522 wallcycle_start(wcycle,ewcRUN);
1524 fprintf(fplog,"\n");
1526 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
1528 chkpt_ret=fcCheckPointParallel( cr->nodeid,
1530 if ( chkpt_ret == 0 )
1531 gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
1535 /***********************************************************
1537 * Loop over MD steps
1539 ************************************************************/
1541 /* if rerunMD then read coordinates and velocities from input trajectory */
1544 if (getenv("GMX_FORCE_UPDATE"))
1546 bForceUpdate = TRUE;
1549 rerun_fr.natoms = 0;
1552 bNotLastFrame = read_first_frame(oenv,&status,
1553 opt2fn("-rerun",nfile,fnm),
1554 &rerun_fr,TRX_NEED_X | TRX_READ_V);
1555 if (rerun_fr.natoms != top_global->natoms)
1558 "Number of atoms in trajectory (%d) does not match the "
1559 "run input file (%d)\n",
1560 rerun_fr.natoms,top_global->natoms);
1562 if (ir->ePBC != epbcNONE)
1566 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);
1568 if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
1570 gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
1577 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
1580 if (ir->ePBC != epbcNONE)
1582 /* Set the shift vectors.
1583 * Necessary here when have a static box different from the tpr box.
1585 calc_shifts(rerun_fr.box,fr->shift_vec);
1589 /* loop over MD steps or if rerunMD to end of input trajectory */
1591 /* Skip the first Nose-Hoover integration when we get the state from tpx */
1592 bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
1593 bInitStep = bFirstStep && (bStateFromTPX || bVV);
1594 bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
1596 bSumEkinhOld = FALSE;
1599 init_global_signals(&gs,cr,ir,repl_ex_nst);
1601 step = ir->init_step;
1604 if (ir->nstlist == -1)
1606 init_nlistheuristics(&nlh,bGStatEveryStep,step);
1609 bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
1610 while (!bLastStep || (bRerunMD && bNotLastFrame)) {
1612 wallcycle_start(wcycle,ewcSTEP);
1614 GMX_MPE_LOG(ev_timestep1);
1617 if (rerun_fr.bStep) {
1618 step = rerun_fr.step;
1619 step_rel = step - ir->init_step;
1621 if (rerun_fr.bTime) {
1631 bLastStep = (step_rel == ir->nsteps);
1632 t = t0 + step*ir->delta_t;
1635 if (ir->efep != efepNO)
1637 if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
1639 state_global->lambda = rerun_fr.lambda;
1643 state_global->lambda = lam0 + step*ir->delta_lambda;
1645 state->lambda = state_global->lambda;
1646 bDoDHDL = do_per_step(step,ir->nstdhdl);
1651 update_annealing_target_temp(&(ir->opts),t);
1656 if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
1658 for(i=0; i<state_global->natoms; i++)
1660 copy_rvec(rerun_fr.x[i],state_global->x[i]);
1664 for(i=0; i<state_global->natoms; i++)
1666 copy_rvec(rerun_fr.v[i],state_global->v[i]);
1671 for(i=0; i<state_global->natoms; i++)
1673 clear_rvec(state_global->v[i]);
1677 fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
1678 " Ekin, temperature and pressure are incorrect,\n"
1679 " the virial will be incorrect when constraints are present.\n"
1681 bRerunWarnNoV = FALSE;
1685 copy_mat(rerun_fr.box,state_global->box);
1686 copy_mat(state_global->box,state->box);
1688 if (vsite && (Flags & MD_RERUN_VSITE))
1690 if (DOMAINDECOMP(cr))
1692 gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
1696 /* Following is necessary because the graph may get out of sync
1697 * with the coordinates if we only have every N'th coordinate set
1699 mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
1700 shift_self(graph,state->box,state->x);
1702 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
1703 top->idef.iparams,top->idef.il,
1704 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
1707 unshift_self(graph,state->box,state->x);
1712 /* Stop Center of Mass motion */
1713 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
1715 /* Copy back starting coordinates in case we're doing a forcefield scan */
1718 for(ii=0; (ii<state->natoms); ii++)
1720 copy_rvec(xcopy[ii],state->x[ii]);
1721 copy_rvec(vcopy[ii],state->v[ii]);
1723 copy_mat(boxcopy,state->box);
1728 /* for rerun MD always do Neighbour Searching */
1729 bNS = (bFirstStep || ir->nstlist != 0);
1734 /* Determine whether or not to do Neighbour Searching and LR */
1735 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
1737 bNS = (bFirstStep || bExchanged || bNStList ||
1738 (ir->nstlist == -1 && nlh.nabnsb > 0));
1740 if (bNS && ir->nstlist == -1)
1742 set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
1746 /* < 0 means stop at next step, > 0 means stop at next NS step */
1747 if ( (gs.set[eglsSTOPCOND] < 0 ) ||
1748 ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
1753 /* Determine whether or not to update the Born radii if doing GB */
1754 bBornRadii=bFirstStep;
1755 if (ir->implicit_solvent && (step % ir->nstgbradii==0))
1760 do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
1761 do_verbose = bVerbose &&
1762 (step % stepout == 0 || bFirstStep || bLastStep);
1764 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
1768 bMasterState = TRUE;
1772 bMasterState = FALSE;
1773 /* Correct the new box if it is too skewed */
1774 if (DYNAMIC_BOX(*ir))
1776 if (correct_box(fplog,step,state->box,graph))
1778 bMasterState = TRUE;
1781 if (DOMAINDECOMP(cr) && bMasterState)
1783 dd_collect_state(cr->dd,state,state_global);
1787 if (DOMAINDECOMP(cr))
1789 /* Repartition the domain decomposition */
1790 wallcycle_start(wcycle,ewcDOMDEC);
1791 dd_partition_system(fplog,step,cr,
1792 bMasterState,nstglobalcomm,
1793 state_global,top_global,ir,
1794 state,&f,mdatoms,top,fr,
1795 vsite,shellfc,constr,
1796 nrnb,wcycle,do_verbose);
1797 wallcycle_stop(wcycle,ewcDOMDEC);
1798 /* If using an iterative integrator, reallocate space to match the decomposition */
1802 if (MASTER(cr) && do_log && !bFFscan)
1804 print_ebin_header(fplog,step,t,state->lambda);
1807 if (ir->efep != efepNO)
1809 update_mdatoms(mdatoms,state->lambda);
1812 if (bRerunMD && rerun_fr.bV)
1815 /* We need the kinetic energy at minus the half step for determining
1816 * the full step kinetic energy and possibly for T-coupling.*/
1817 /* This may not be quite working correctly yet . . . . */
1818 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
1819 wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
1820 constr,NULL,FALSE,state->box,
1821 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
1822 CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
1824 clear_mat(force_vir);
1826 /* Ionize the atoms if necessary */
1829 ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
1830 mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
1833 /* Update force field in ffscan program */
1836 if (update_forcefield(fplog,
1838 mdatoms->nr,state->x,state->box)) {
1839 if (gmx_parallel_env_initialized())
1847 GMX_MPE_LOG(ev_timestep2);
1849 /* We write a checkpoint at this MD step when:
1850 * either at an NS step when we signalled through gs,
1851 * or at the last step (but not when we do not want confout),
1852 * but never at the first step or with rerun.
1854 bCPT = (((gs.set[eglsCHKPT] && bNS) ||
1855 (bLastStep && (Flags & MD_CONFOUT))) &&
1856 step > ir->init_step && !bRerunMD);
1859 gs.set[eglsCHKPT] = 0;
1862 /* Determine the energy and pressure:
1863 * at nstcalcenergy steps and at energy output steps (set below).
1865 bNstEner = do_per_step(step,ir->nstcalcenergy);
1868 (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)));
1870 /* Do we need global communication ? */
1871 bGStat = (bCalcEnerPres || bStopCM ||
1872 do_per_step(step,nstglobalcomm) ||
1873 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
1875 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
1877 if (do_ene || do_log)
1879 bCalcEnerPres = TRUE;
1883 /* these CGLO_ options remain the same throughout the iteration */
1884 cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
1885 (bStopCM ? CGLO_STOPCM : 0) |
1886 (bGStat ? CGLO_GSTAT : 0)
1889 force_flags = (GMX_FORCE_STATECHANGED |
1890 ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1891 GMX_FORCE_ALLFORCES |
1892 (bNStList ? GMX_FORCE_DOLR : 0) |
1894 (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
1895 (bDoDHDL ? GMX_FORCE_DHDL : 0)
1900 /* Now is the time to relax the shells */
1901 count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
1903 bStopCM,top,top_global,
1905 state,f,force_vir,mdatoms,
1906 nrnb,wcycle,graph,groups,
1907 shellfc,fr,bBornRadii,t,mu_tot,
1908 state->natoms,&bConverged,vsite,
1919 /* The coordinates (x) are shifted (to get whole molecules)
1921 * This is parallellized as well, and does communication too.
1922 * Check comments in sim_util.c
1925 do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
1926 state->box,state->x,&state->hist,
1927 f,force_vir,mdatoms,enerd,fcd,
1928 state->lambda,graph,
1929 fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
1930 (bNS ? GMX_FORCE_NS : 0) | force_flags);
1933 GMX_BARRIER(cr->mpi_comm_mygroup);
1937 mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
1938 mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
1941 if (bTCR && bFirstStep)
1943 tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
1944 fprintf(fplog,"Done init_coupling\n");
1948 if (bVV && !bStartingFromCpt && !bRerunMD)
1949 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1951 if (ir->eI==eiVV && bInitStep)
1953 /* if using velocity verlet with full time step Ekin,
1954 * take the first half step only to compute the
1955 * virial for the first step. From there,
1956 * revert back to the initial coordinates
1957 * so that the input is actually the initial step.
1959 copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
1961 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1962 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
1965 update_coords(fplog,step,ir,mdatoms,state,
1966 f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
1967 ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
1968 cr,nrnb,constr,&top->idef);
1972 gmx_iterate_init(&iterate,bIterations && !bInitStep);
1974 /* for iterations, we save these vectors, as we will be self-consistently iterating
1977 /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
1979 /* save the state */
1980 if (bIterations && iterate.bIterate) {
1981 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
1984 bFirstIterate = TRUE;
1985 while (bFirstIterate || (bIterations && iterate.bIterate))
1987 if (bIterations && iterate.bIterate)
1989 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
1990 if (bFirstIterate && bTrotter)
1992 /* The first time through, we need a decent first estimate
1993 of veta(t+dt) to compute the constraints. Do
1994 this by computing the box volume part of the
1995 trotter integration at this time. Nothing else
1996 should be changed by this routine here. If
1997 !(first time), we start with the previous value
2000 veta_save = state->veta;
2001 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
2002 vetanew = state->veta;
2003 state->veta = veta_save;
2008 if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
2011 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2012 &top->idef,shake_vir,NULL,
2013 cr,nrnb,wcycle,upd,constr,
2014 bInitStep,TRUE,bCalcEnerPres,vetanew);
2016 if (!bOK && !bFFscan)
2018 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2023 { /* Need to unshift here if a do_force has been
2024 called in the previous step */
2025 unshift_self(graph,state->box,state->x);
2029 /* if VV, compute the pressure and constraints */
2030 /* For VV2, we strictly only need this if using pressure
2031 * control, but we really would like to have accurate pressures
2033 * Think about ways around this in the future?
2034 * For now, keep this choice in comments.
2036 /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
2037 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
2039 bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
2040 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2041 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2042 constr,NULL,FALSE,state->box,
2043 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2046 | (bTemp ? CGLO_TEMPERATURE:0)
2047 | (bPres ? CGLO_PRESSURE : 0)
2048 | (bPres ? CGLO_CONSTRAINT : 0)
2049 | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
2050 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2053 /* explanation of above:
2054 a) We compute Ekin at the full time step
2055 if 1) we are using the AveVel Ekin, and it's not the
2056 initial step, or 2) if we are using AveEkin, but need the full
2057 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
2058 b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
2059 EkinAveVel because it's needed for the pressure */
2061 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
2066 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
2070 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2075 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2076 state->veta,&vetanew))
2080 bFirstIterate = FALSE;
2083 if (bTrotter && !bInitStep) {
2084 copy_mat(shake_vir,state->svir_prev);
2085 copy_mat(force_vir,state->fvir_prev);
2086 if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
2087 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
2088 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
2089 enerd->term[F_EKIN] = trace(ekind->ekin);
2092 /* if it's the initial step, we performed this first step just to get the constraint virial */
2093 if (bInitStep && ir->eI==eiVV) {
2094 copy_rvecn(cbuf,state->v,0,state->natoms);
2097 if (fr->bSepDVDL && fplog && do_log)
2099 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2101 enerd->term[F_DHDL_CON] += dvdl;
2103 GMX_MPE_LOG(ev_timestep1);
2106 /* MRS -- now done iterating -- compute the conserved quantity */
2108 saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
2111 last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
2113 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
2115 saved_conserved_quantity -= enerd->term[F_DISPCORR];
2119 /* ######## END FIRST UPDATE STEP ############## */
2120 /* ######## If doing VV, we now have v(dt) ###### */
2122 /* ################## START TRAJECTORY OUTPUT ################# */
2124 /* Now we have the energies and forces corresponding to the
2125 * coordinates at time t. We must output all of this before
2127 * for RerunMD t is read from input trajectory
2129 GMX_MPE_LOG(ev_output_start);
2132 if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
2133 if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
2134 if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
2135 if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
2136 if (bCPT) { mdof_flags |= MDOF_CPT; };
2138 #if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
2141 /* Enforce writing positions and velocities at end of run */
2142 mdof_flags |= (MDOF_X | MDOF_V);
2147 fcReportProgress( ir->nsteps, step );
2149 /* sync bCPT and fc record-keeping */
2150 if (bCPT && MASTER(cr))
2151 fcRequestCheckPoint();
2154 if (mdof_flags != 0)
2156 wallcycle_start(wcycle,ewcTRAJ);
2159 if (state->flags & (1<<estLD_RNG))
2161 get_stochd_state(upd,state);
2167 state_global->ekinstate.bUpToDate = FALSE;
2171 update_ekinstate(&state_global->ekinstate,ekind);
2172 state_global->ekinstate.bUpToDate = TRUE;
2174 update_energyhistory(&state_global->enerhist,mdebin);
2177 write_traj(fplog,cr,outf,mdof_flags,top_global,
2178 step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
2185 if (bLastStep && step_rel == ir->nsteps &&
2186 (Flags & MD_CONFOUT) && MASTER(cr) &&
2187 !bRerunMD && !bFFscan)
2189 /* x and v have been collected in write_traj,
2190 * because a checkpoint file will always be written
2193 fprintf(stderr,"\nWriting final coordinates.\n");
2194 if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
2197 /* Make molecules whole only for confout writing */
2198 do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
2200 write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
2201 *top_global->name,top_global,
2202 state_global->x,state_global->v,
2203 ir->ePBC,state->box);
2206 wallcycle_stop(wcycle,ewcTRAJ);
2208 GMX_MPE_LOG(ev_output_finish);
2210 /* kludge -- virial is lost with restart for NPT control. Must restart */
2211 if (bStartingFromCpt && bVV)
2213 copy_mat(state->svir_prev,shake_vir);
2214 copy_mat(state->fvir_prev,force_vir);
2216 /* ################## END TRAJECTORY OUTPUT ################ */
2218 /* Determine the pressure:
2219 * always when we want exact averages in the energy file,
2220 * at ns steps when we have pressure coupling,
2221 * otherwise only at energy output steps (set below).
2224 bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
2225 bCalcEnerPres = bNstEner;
2227 /* Do we need global communication ? */
2228 bGStat = (bGStatEveryStep || bStopCM || bNS ||
2229 (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
2231 do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
2233 if (do_ene || do_log)
2235 bCalcEnerPres = TRUE;
2239 /* Determine the wallclock run time up till now */
2240 run_time = gmx_gettime() - (double)runtime->real;
2242 /* Check whether everything is still allright */
2243 if (((int)gmx_get_stop_condition() > handled_stop_condition)
2249 /* this is just make gs.sig compatible with the hack
2250 of sending signals around by MPI_Reduce with together with
2252 if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
2253 gs.sig[eglsSTOPCOND]=1;
2254 if ( gmx_get_stop_condition() == gmx_stop_cond_next )
2255 gs.sig[eglsSTOPCOND]=-1;
2256 /* < 0 means stop at next step, > 0 means stop at next NS step */
2260 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2261 gmx_get_signal_name(),
2262 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2266 "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
2267 gmx_get_signal_name(),
2268 gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
2270 handled_stop_condition=(int)gmx_get_stop_condition();
2272 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
2273 (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
2274 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
2276 /* Signal to terminate the run */
2277 gs.sig[eglsSTOPCOND] = 1;
2280 fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2282 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
2285 if (bResetCountersHalfMaxH && MASTER(cr) &&
2286 run_time > max_hours*60.0*60.0*0.495)
2288 gs.sig[eglsRESETCOUNTERS] = 1;
2291 if (ir->nstlist == -1 && !bRerunMD)
2293 /* When bGStatEveryStep=FALSE, global_stat is only called
2294 * when we check the atom displacements, not at NS steps.
2295 * This means that also the bonded interaction count check is not
2296 * performed immediately after NS. Therefore a few MD steps could
2297 * be performed with missing interactions.
2298 * But wrong energies are never written to file,
2299 * since energies are only written after global_stat
2302 if (step >= nlh.step_nscheck)
2304 nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
2305 nlh.scale_tot,state->x);
2309 /* This is not necessarily true,
2310 * but step_nscheck is determined quite conservatively.
2316 /* In parallel we only have to check for checkpointing in steps
2317 * where we do global communication,
2318 * otherwise the other nodes don't know.
2320 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
2323 run_time >= nchkpt*cpt_period*60.0)) &&
2324 gs.set[eglsCHKPT] == 0)
2326 gs.sig[eglsCHKPT] = 1;
2331 gmx_iterate_init(&iterate,bIterations);
2334 /* for iterations, we save these vectors, as we will be redoing the calculations */
2335 if (bIterations && iterate.bIterate)
2337 copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
2339 bFirstIterate = TRUE;
2340 while (bFirstIterate || (bIterations && iterate.bIterate))
2342 /* We now restore these vectors to redo the calculation with improved extended variables */
2345 copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
2348 /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
2349 so scroll down for that logic */
2351 /* ######### START SECOND UPDATE STEP ################# */
2352 GMX_MPE_LOG(ev_update_start);
2353 /* Box is changed in update() when we do pressure coupling,
2354 * but we should still use the old box for energy corrections and when
2355 * writing it to the energy file, so it matches the trajectory files for
2356 * the same timestep above. Make a copy in a separate array.
2358 copy_mat(state->box,lastbox);
2361 if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
2363 wallcycle_start(wcycle,ewcUPDATE);
2365 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
2368 if (bIterations && iterate.bIterate)
2376 /* we use a new value of scalevir to converge the iterations faster */
2377 scalevir = tracevir/trace(shake_vir);
2379 msmul(shake_vir,scalevir,shake_vir);
2380 m_add(force_vir,shake_vir,total_vir);
2381 clear_mat(shake_vir);
2383 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
2384 /* We can only do Berendsen coupling after we have summed
2385 * the kinetic energy or virial. Since the happens
2386 * in global_state after update, we should only do it at
2387 * step % nstlist = 1 with bGStatEveryStep=FALSE.
2392 update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
2393 update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
2399 /* velocity half-step update */
2400 update_coords(fplog,step,ir,mdatoms,state,f,
2401 fr->bTwinRange && bNStList,fr->f_twin,fcd,
2402 ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
2403 cr,nrnb,constr,&top->idef);
2406 /* Above, initialize just copies ekinh into ekin,
2407 * it doesn't copy position (for VV),
2408 * and entire integrator for MD.
2413 copy_rvecn(state->x,cbuf,0,state->natoms);
2416 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2417 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2418 wallcycle_stop(wcycle,ewcUPDATE);
2420 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2421 &top->idef,shake_vir,force_vir,
2422 cr,nrnb,wcycle,upd,constr,
2423 bInitStep,FALSE,bCalcEnerPres,state->veta);
2427 /* erase F_EKIN and F_TEMP here? */
2428 /* just compute the kinetic energy at the half step to perform a trotter step */
2429 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2430 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2431 constr,NULL,FALSE,lastbox,
2432 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2433 cglo_flags | CGLO_TEMPERATURE
2435 wallcycle_start(wcycle,ewcUPDATE);
2436 trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
2437 /* now we know the scaling, we can compute the positions again again */
2438 copy_rvecn(cbuf,state->x,0,state->natoms);
2440 update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
2441 ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
2442 wallcycle_stop(wcycle,ewcUPDATE);
2444 /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
2445 /* are the small terms in the shake_vir here due
2446 * to numerical errors, or are they important
2447 * physically? I'm thinking they are just errors, but not completely sure.
2448 * For now, will call without actually constraining, constr=NULL*/
2449 update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
2450 &top->idef,tmp_vir,force_vir,
2451 cr,nrnb,wcycle,upd,NULL,
2452 bInitStep,FALSE,bCalcEnerPres,
2455 if (!bOK && !bFFscan)
2457 gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
2460 if (fr->bSepDVDL && fplog && do_log)
2462 fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
2464 enerd->term[F_DHDL_CON] += dvdl;
2468 /* Need to unshift here */
2469 unshift_self(graph,state->box,state->x);
2472 GMX_BARRIER(cr->mpi_comm_mygroup);
2473 GMX_MPE_LOG(ev_update_finish);
2477 wallcycle_start(wcycle,ewcVSITECONSTR);
2480 shift_self(graph,state->box,state->x);
2482 construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
2483 top->idef.iparams,top->idef.il,
2484 fr->ePBC,fr->bMolPBC,graph,cr,state->box);
2488 unshift_self(graph,state->box,state->x);
2490 wallcycle_stop(wcycle,ewcVSITECONSTR);
2493 /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
2494 if (ir->nstlist == -1 && bFirstIterate)
2496 gs.sig[eglsNABNSB] = nlh.nabnsb;
2498 compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
2499 wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
2501 bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
2503 top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
2505 | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
2506 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
2507 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
2508 | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
2509 | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
2512 if (ir->nstlist == -1 && bFirstIterate)
2514 nlh.nabnsb = gs.set[eglsNABNSB];
2515 gs.set[eglsNABNSB] = 0;
2517 /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
2518 /* ############# END CALC EKIN AND PRESSURE ################# */
2520 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
2521 the virial that should probably be addressed eventually. state->veta has better properies,
2522 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
2523 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
2526 done_iterating(cr,fplog,step,&iterate,bFirstIterate,
2527 trace(shake_vir),&tracevir))
2531 bFirstIterate = FALSE;
2534 update_box(fplog,step,ir,mdatoms,state,graph,f,
2535 ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
2537 /* ################# END UPDATE STEP 2 ################# */
2538 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
2540 /* The coordinates (x) were unshifted in update */
2541 if (bFFscan && (shellfc==NULL || bConverged))
2543 if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
2545 &(top_global->mols),mdatoms->massT,pres))
2547 if (gmx_parallel_env_initialized())
2551 fprintf(stderr,"\n");
2557 /* We will not sum ekinh_old,
2558 * so signal that we still have to do it.
2560 bSumEkinhOld = TRUE;
2565 /* Only do GCT when the relaxation of shells (minimization) has converged,
2566 * otherwise we might be coupling to bogus energies.
2567 * In parallel we must always do this, because the other sims might
2571 /* Since this is called with the new coordinates state->x, I assume
2572 * we want the new box state->box too. / EL 20040121
2574 do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
2576 mdatoms,&(top->idef),mu_aver,
2577 top_global->mols.nr,cr,
2578 state->box,total_vir,pres,
2579 mu_tot,state->x,f,bConverged);
2583 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
2585 /* sum up the foreign energy and dhdl terms */
2586 sum_dhdl(enerd,state->lambda,ir);
2588 /* use the directly determined last velocity, not actually the averaged half steps */
2589 if (bTrotter && ir->eI==eiVV)
2591 enerd->term[F_EKIN] = last_ekin;
2593 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
2597 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
2601 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
2603 /* Check for excessively large energies */
2607 real etot_max = 1e200;
2609 real etot_max = 1e30;
2611 if (fabs(enerd->term[F_ETOT]) > etot_max)
2613 fprintf(stderr,"Energy too large (%g), giving up\n",
2614 enerd->term[F_ETOT]);
2617 /* ######### END PREPARING EDR OUTPUT ########### */
2619 /* Time for performance */
2620 if (((step % stepout) == 0) || bLastStep)
2622 runtime_upd_proc(runtime);
2628 gmx_bool do_dr,do_or;
2630 if (!(bStartingFromCpt && (EI_VV(ir->eI))))
2634 upd_mdebin(mdebin,bDoDHDL, TRUE,
2635 t,mdatoms->tmass,enerd,state,lastbox,
2636 shake_vir,force_vir,total_vir,pres,
2637 ekind,mu_tot,constr);
2641 upd_mdebin_step(mdebin);
2644 do_dr = do_per_step(step,ir->nstdisreout);
2645 do_or = do_per_step(step,ir->nstorireout);
2647 print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
2649 eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
2651 if (ir->ePull != epullNO)
2653 pull_print_output(ir->pull,step,t);
2656 if (do_per_step(step,ir->nstlog))
2658 if(fflush(fplog) != 0)
2660 gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
2666 /* Remaining runtime */
2667 if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
2671 fprintf(stderr,"\n");
2673 print_time(stderr,runtime,step,ir,cr);
2676 /* Replica exchange */
2678 if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
2679 do_per_step(step,repl_ex_nst))
2681 bExchanged = replica_exchange(fplog,cr,repl_ex,
2682 state_global,enerd->term,
2685 if (bExchanged && DOMAINDECOMP(cr))
2687 dd_partition_system(fplog,step,cr,TRUE,1,
2688 state_global,top_global,ir,
2689 state,&f,mdatoms,top,fr,
2690 vsite,shellfc,constr,
2697 bStartingFromCpt = FALSE;
2699 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
2700 /* Complicated conditional when bGStatEveryStep=FALSE.
2701 * We can not just use bGStat, since then the simulation results
2702 * would depend on nstenergy and nstlog or step_nscheck.
2704 if (((state->flags & (1<<estPRES_PREV)) ||
2705 (state->flags & (1<<estSVIR_PREV)) ||
2706 (state->flags & (1<<estFVIR_PREV))) &&
2708 (ir->nstlist > 0 && step % ir->nstlist == 0) ||
2709 (ir->nstlist < 0 && nlh.nabnsb > 0) ||
2710 (ir->nstlist == 0 && bGStat)))
2712 /* Store the pressure in t_state for pressure coupling
2713 * at the next MD step.
2715 if (state->flags & (1<<estPRES_PREV))
2717 copy_mat(pres,state->pres_prev);
2721 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
2727 /* read next frame from input trajectory */
2728 bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
2733 rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
2737 if (!bRerunMD || !rerun_fr.bStep)
2739 /* increase the MD step number */
2744 cycles = wallcycle_stop(wcycle,ewcSTEP);
2745 if (DOMAINDECOMP(cr) && wcycle)
2747 dd_cycles_add(cr->dd,cycles,ddCyclStep);
2750 if (step_rel == wcycle_get_reset_counters(wcycle) ||
2751 gs.set[eglsRESETCOUNTERS] != 0)
2753 /* Reset all the counters related to performance over the run */
2754 reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
2755 wcycle_set_reset_counters(wcycle,-1);
2756 /* Correct max_hours for the elapsed time */
2757 max_hours -= run_time/(60.0*60.0);
2758 bResetCountersHalfMaxH = FALSE;
2759 gs.set[eglsRESETCOUNTERS] = 0;
2762 /* End of main MD loop */
2766 runtime_end(runtime);
2768 if (bRerunMD && MASTER(cr))
2773 if (!(cr->duty & DUTY_PME))
2775 /* Tell the PME only node to finish */
2781 if (ir->nstcalcenergy > 0 && !bRerunMD)
2783 print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
2784 eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
2792 if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
2794 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)));
2795 fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
2798 if (shellfc && fplog)
2800 fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
2801 (nconverged*100.0)/step_rel);
2802 fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
2806 if (repl_ex_nst > 0 && MASTER(cr))
2808 print_replica_exchange_statistics(fplog,repl_ex);
2811 runtime->nsteps_done = step_rel;