2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
47 #include "mtop_util.h"
48 #include "gmx_wallcycle.h"
51 #include "md_logging.h"
52 #include "md_support.h"
55 /* Is the signal in one simulation independent of other simulations? */
56 gmx_bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
58 /* check which of the multisim simulations has the shortest number of
59 steps and return that number of nsteps */
60 gmx_large_int_t get_multisim_nsteps(const t_commrec *cr,
61 gmx_large_int_t nsteps)
63 gmx_large_int_t steps_out;
70 snew(buf, cr->ms->nsim);
72 buf[cr->ms->sim] = nsteps;
73 gmx_sumli_sim(cr->ms->nsim, buf, cr->ms);
76 for (s = 0; s < cr->ms->nsim; s++)
78 /* find the smallest positive number */
79 if (buf[s] >= 0 && ((steps_out < 0) || (buf[s] < steps_out)) )
86 /* if we're the limiting simulation, don't do anything */
87 if (steps_out >= 0 && steps_out < nsteps)
90 snprintf(strbuf, 255, "Will stop simulation %%d after %s steps (another simulation will end then).\n", gmx_large_int_pfmt);
91 fprintf(stderr, strbuf, cr->ms->sim, steps_out);
94 /* broadcast to non-masters */
95 gmx_bcast(sizeof(gmx_large_int_t), &steps_out, cr);
99 int multisim_min(const gmx_multisim_t *ms, int nmin, int n)
102 gmx_bool bPos, bEqual;
107 gmx_sumi_sim(ms->nsim, buf, ms);
110 for (s = 0; s < ms->nsim; s++)
112 bPos = bPos && (buf[s] > 0);
113 bEqual = bEqual && (buf[s] == buf[0]);
119 nmin = min(nmin, buf[0]);
123 /* Find the least common multiple */
124 for (d = 2; d < nmin; d++)
127 while (s < ms->nsim && d % buf[s] == 0)
133 /* We found the LCM and it is less than nmin */
145 int multisim_nstsimsync(const t_commrec *cr,
146 const t_inputrec *ir, int repl_ex_nst)
153 nmin = multisim_min(cr->ms, nmin, ir->nstlist);
154 nmin = multisim_min(cr->ms, nmin, ir->nstcalcenergy);
155 nmin = multisim_min(cr->ms, nmin, repl_ex_nst);
158 gmx_fatal(FARGS, "Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
160 /* Avoid inter-simulation communication at every (second) step */
167 gmx_bcast(sizeof(int), &nmin, cr);
172 void init_global_signals(globsig_t *gs, const t_commrec *cr,
173 const t_inputrec *ir, int repl_ex_nst)
179 gs->nstms = multisim_nstsimsync(cr, ir, repl_ex_nst);
182 fprintf(debug, "Syncing simulations for checkpointing and termination every %d steps\n", gs->nstms);
190 for (i = 0; i < eglsNR; i++)
197 void copy_coupling_state(t_state *statea, t_state *stateb,
198 gmx_ekindata_t *ekinda, gmx_ekindata_t *ekindb, t_grpopts* opts)
201 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
205 /* Make sure we have enough space for x and v */
206 if (statea->nalloc > stateb->nalloc)
208 stateb->nalloc = statea->nalloc;
209 srenew(stateb->x, stateb->nalloc);
210 srenew(stateb->v, stateb->nalloc);
213 stateb->natoms = statea->natoms;
214 stateb->ngtc = statea->ngtc;
215 stateb->nnhpres = statea->nnhpres;
216 stateb->veta = statea->veta;
219 copy_mat(ekinda->ekin, ekindb->ekin);
220 for (i = 0; i < stateb->ngtc; i++)
222 ekindb->tcstat[i].T = ekinda->tcstat[i].T;
223 ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
224 copy_mat(ekinda->tcstat[i].ekinh, ekindb->tcstat[i].ekinh);
225 copy_mat(ekinda->tcstat[i].ekinf, ekindb->tcstat[i].ekinf);
226 ekindb->tcstat[i].ekinscalef_nhc = ekinda->tcstat[i].ekinscalef_nhc;
227 ekindb->tcstat[i].ekinscaleh_nhc = ekinda->tcstat[i].ekinscaleh_nhc;
228 ekindb->tcstat[i].vscale_nhc = ekinda->tcstat[i].vscale_nhc;
231 copy_rvecn(statea->x, stateb->x, 0, stateb->natoms);
232 copy_rvecn(statea->v, stateb->v, 0, stateb->natoms);
233 copy_mat(statea->box, stateb->box);
234 copy_mat(statea->box_rel, stateb->box_rel);
235 copy_mat(statea->boxv, stateb->boxv);
237 for (i = 0; i < stateb->ngtc; i++)
239 nc = i*opts->nhchainlength;
240 for (j = 0; j < opts->nhchainlength; j++)
242 stateb->nosehoover_xi[nc+j] = statea->nosehoover_xi[nc+j];
243 stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
246 if (stateb->nhpres_xi != NULL)
248 for (i = 0; i < stateb->nnhpres; i++)
250 nc = i*opts->nhchainlength;
251 for (j = 0; j < opts->nhchainlength; j++)
253 stateb->nhpres_xi[nc+j] = statea->nhpres_xi[nc+j];
254 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
260 real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
270 quantity = NPT_energy(ir, state, MassQ);
273 quantity = vrescale_energy(&(ir->opts), state->therm_integral);
281 void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
282 t_forcerec *fr, gmx_ekindata_t *ekind,
283 t_state *state, t_state *state_global, t_mdatoms *mdatoms,
284 t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
285 gmx_enerdata_t *enerd, tensor force_vir, tensor shake_vir, tensor total_vir,
286 tensor pres, rvec mu_tot, gmx_constr_t constr,
287 globsig_t *gs, gmx_bool bInterSimGS,
288 matrix box, gmx_mtop_t *top_global, real *pcurr,
289 int natoms, gmx_bool *bSumEkinhOld, int flags)
293 tensor corr_vir, corr_pres, shakeall_vir;
294 gmx_bool bEner, bPres, bTemp, bVV;
295 gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
296 bFirstIterate, bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
297 real ekin, temp, prescorr, enercorr, dvdlcorr, dvdl_ekin;
299 /* translate CGLO flags to gmx_booleans */
300 bRerunMD = flags & CGLO_RERUNMD;
301 bStopCM = flags & CGLO_STOPCM;
302 bGStat = flags & CGLO_GSTAT;
304 bReadEkin = (flags & CGLO_READEKIN);
305 bScaleEkin = (flags & CGLO_SCALEEKIN);
306 bEner = flags & CGLO_ENERGY;
307 bTemp = flags & CGLO_TEMPERATURE;
308 bPres = (flags & CGLO_PRESSURE);
309 bConstrain = (flags & CGLO_CONSTRAINT);
310 bIterate = (flags & CGLO_ITERATE);
311 bFirstIterate = (flags & CGLO_FIRSTITERATE);
313 /* we calculate a full state kinetic energy either with full-step velocity verlet
314 or half step where we need the pressure */
316 bEkinAveVel = (ir->eI == eiVV || (ir->eI == eiVVAK && bPres) || bReadEkin);
318 /* in initalization, it sums the shake virial in vv, and to
319 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
321 /* ########## Kinetic energy ############## */
325 /* Non-equilibrium MD: this is parallellized, but only does communication
326 * when there really is NEMD.
329 if (PAR(cr) && (ekind->bNEMD))
331 accumulate_u(cr, &(ir->opts), ekind);
336 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
341 calc_ke_part(state, &(ir->opts), mdatoms, ekind, nrnb, bEkinAveVel, bIterate);
347 /* Calculate center of mass velocity if necessary, also parallellized */
350 calc_vcm_grp(fplog, mdatoms->start, mdatoms->homenr, mdatoms,
351 state->x, state->v, vcm);
354 if (bTemp || bStopCM || bPres || bEner || bConstrain)
358 /* We will not sum ekinh_old,
359 * so signal that we still have to do it.
361 *bSumEkinhOld = TRUE;
368 for (i = 0; i < eglsNR; i++)
370 gs_buf[i] = gs->sig[i];
375 wallcycle_start(wcycle, ewcMoveE);
376 GMX_MPE_LOG(ev_global_stat_start);
377 global_stat(fplog, gstat, cr, enerd, force_vir, shake_vir, mu_tot,
378 ir, ekind, constr, bStopCM ? vcm : NULL,
379 gs != NULL ? eglsNR : 0, gs_buf,
381 *bSumEkinhOld, flags);
382 GMX_MPE_LOG(ev_global_stat_finish);
383 wallcycle_stop(wcycle, ewcMoveE);
387 if (MULTISIM(cr) && bInterSimGS)
391 /* Communicate the signals between the simulations */
392 gmx_sum_sim(eglsNR, gs_buf, cr->ms);
394 /* Communicate the signals form the master to the others */
395 gmx_bcast(eglsNR*sizeof(gs_buf[0]), gs_buf, cr);
397 for (i = 0; i < eglsNR; i++)
399 if (bInterSimGS || gs_simlocal[i])
401 /* Set the communicated signal only when it is non-zero,
402 * since signals might not be processed at each MD step.
404 gsi = (gs_buf[i] >= 0 ?
405 (int)(gs_buf[i] + 0.5) :
406 (int)(gs_buf[i] - 0.5));
411 /* Turn off the local signal */
416 *bSumEkinhOld = FALSE;
420 if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
423 mdatoms->start, mdatoms->start+mdatoms->homenr,
424 state->v, vcm->group_p[0],
425 mdatoms->massT, mdatoms->tmass, ekind->ekin);
428 /* Do center of mass motion removal */
431 check_cm_grp(fplog, vcm, ir, 1);
432 do_stopcm_grp(fplog, mdatoms->start, mdatoms->homenr, mdatoms->cVCM,
433 state->x, state->v, vcm);
434 inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
439 /* Calculate the amplitude of the cosine velocity profile */
440 ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
445 /* Sum the kinetic energies of the groups & calc temp */
446 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
447 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
448 Leap with AveVel is not supported; it's not clear that it will actually work.
449 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
450 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
451 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
452 If FALSE, we go ahead and erase over it.
454 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin,
455 bEkinAveVel, bIterate, bScaleEkin);
456 enerd->dvdl_lin[efptMASS] = (double) dvdl_ekin;
458 enerd->term[F_EKIN] = trace(ekind->ekin);
461 /* ########## Long range energy information ###### */
463 if (bEner || bPres || bConstrain)
465 calc_dispcorr(fplog, ir, fr, 0, top_global->natoms, box, state->lambda[efptVDW],
466 corr_pres, corr_vir, &prescorr, &enercorr, &dvdlcorr);
469 if (bEner && bFirstIterate)
471 enerd->term[F_DISPCORR] = enercorr;
472 enerd->term[F_EPOT] += enercorr;
473 enerd->term[F_DVDL_VDW] += dvdlcorr;
476 /* ########## Now pressure ############## */
477 if (bPres || bConstrain)
480 m_add(force_vir, shake_vir, total_vir);
482 /* Calculate pressure and apply LR correction if PPPM is used.
483 * Use the box from last timestep since we already called update().
486 enerd->term[F_PRES] = calc_pres(fr->ePBC, ir->nwall, box, ekind->ekin, total_vir, pres);
488 /* Calculate long range corrections to pressure and energy */
489 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
490 and computes enerd->term[F_DISPCORR]. Also modifies the
491 total_vir and pres tesors */
493 m_add(total_vir, corr_vir, total_vir);
494 m_add(pres, corr_pres, pres);
495 enerd->term[F_PDISPCORR] = prescorr;
496 enerd->term[F_PRES] += prescorr;
497 *pcurr = enerd->term[F_PRES];
501 void check_nst_param(FILE *fplog, t_commrec *cr,
502 const char *desc_nst, int nst,
503 const char *desc_p, int *p)
505 if (*p > 0 && *p % nst != 0)
507 /* Round up to the next multiple of nst */
508 *p = ((*p)/nst + 1)*nst;
509 md_print_warn(cr, fplog,
510 "NOTE: %s changes %s to %d\n", desc_nst, desc_p, *p);
514 void set_current_lambdas(gmx_large_int_t step, t_lambda *fepvals, gmx_bool bRerunMD,
515 t_trxframe *rerun_fr, t_state *state_global, t_state *state, double lam0[])
516 /* find the current lambdas. If rerunning, we either read in a state, or a lambda value,
517 requiring different logic. */
520 int i, fep_state = 0;
523 if (rerun_fr->bLambda)
525 if (fepvals->delta_lambda==0)
527 state_global->lambda[efptFEP] = rerun_fr->lambda;
528 for (i = 0; i < efptNR; i++)
532 state->lambda[i] = state_global->lambda[i];
538 /* find out between which two value of lambda we should be */
539 frac = (step*fepvals->delta_lambda);
540 fep_state = floor(frac*fepvals->n_lambda);
541 /* interpolate between this state and the next */
542 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
543 frac = (frac*fepvals->n_lambda)-fep_state;
544 for (i = 0; i < efptNR; i++)
546 state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
547 frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
551 else if (rerun_fr->bFepState)
553 state_global->fep_state = rerun_fr->fep_state;
554 for (i = 0; i < efptNR; i++)
556 state_global->lambda[i] = fepvals->all_lambda[i][fep_state];
562 if (fepvals->delta_lambda != 0)
564 /* find out between which two value of lambda we should be */
565 frac = (step*fepvals->delta_lambda);
566 if (fepvals->n_lambda > 0)
568 fep_state = floor(frac*fepvals->n_lambda);
569 /* interpolate between this state and the next */
570 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
571 frac = (frac*fepvals->n_lambda)-fep_state;
572 for (i = 0; i < efptNR; i++)
574 state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
575 frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
580 for (i = 0; i < efptNR; i++)
582 state_global->lambda[i] = lam0[i] + frac;
588 if (state->fep_state > 0) {
589 state_global->fep_state = state->fep_state; /* state->fep is the one updated by bExpanded */
590 for (i = 0; i < efptNR; i++)
592 state_global->lambda[i] = fepvals->all_lambda[i][state_global->fep_state];
597 for (i = 0; i < efptNR; i++)
599 state->lambda[i] = state_global->lambda[i];
603 static void min_zero(int *n, int i)
605 if (i > 0 && (*n == 0 || i < *n))
611 static int lcd4(int i1, int i2, int i3, int i4)
622 gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
625 while (nst > 1 && ((i1 > 0 && i1 % nst != 0) ||
626 (i2 > 0 && i2 % nst != 0) ||
627 (i3 > 0 && i3 % nst != 0) ||
628 (i4 > 0 && i4 % nst != 0)))
636 int check_nstglobalcomm(FILE *fplog, t_commrec *cr,
637 int nstglobalcomm, t_inputrec *ir)
639 if (!EI_DYNAMICS(ir->eI))
644 if (nstglobalcomm == -1)
646 if (!(ir->nstcalcenergy > 0 ||
652 if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
654 nstglobalcomm = ir->nstenergy;
659 /* Ensure that we do timely global communication for
660 * (possibly) each of the four following options.
662 nstglobalcomm = lcd4(ir->nstcalcenergy,
664 ir->etc != etcNO ? ir->nsttcouple : 0,
665 ir->epc != epcNO ? ir->nstpcouple : 0);
670 if (ir->nstlist > 0 &&
671 nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
673 nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
674 md_print_warn(cr, fplog, "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n", nstglobalcomm);
676 if (ir->nstcalcenergy > 0)
678 check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
679 "nstcalcenergy", &ir->nstcalcenergy);
681 if (ir->etc != etcNO && ir->nsttcouple > 0)
683 check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
684 "nsttcouple", &ir->nsttcouple);
686 if (ir->epc != epcNO && ir->nstpcouple > 0)
688 check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
689 "nstpcouple", &ir->nstpcouple);
692 check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
693 "nstenergy", &ir->nstenergy);
695 check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
696 "nstlog", &ir->nstlog);
699 if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
701 md_print_warn(cr, fplog, "WARNING: Changing nstcomm from %d to %d\n",
702 ir->nstcomm, nstglobalcomm);
703 ir->nstcomm = nstglobalcomm;
706 return nstglobalcomm;
709 void check_ir_old_tpx_versions(t_commrec *cr, FILE *fplog,
710 t_inputrec *ir, gmx_mtop_t *mtop)
712 /* Check required for old tpx files */
713 if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
714 ir->nstcalcenergy % ir->nstlist != 0)
716 md_print_warn(cr, fplog, "Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies\n");
718 if (gmx_mtop_ftype_count(mtop, F_CONSTR) +
719 gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0 &&
720 ir->eConstrAlg == econtSHAKE)
722 md_print_warn(cr, fplog, "With twin-range cut-off's and SHAKE the virial and pressure are incorrect\n");
723 if (ir->epc != epcNO)
725 gmx_fatal(FARGS, "Can not do pressure coupling with twin-range cut-off's and SHAKE");
728 check_nst_param(fplog, cr, "nstlist", ir->nstlist,
729 "nstcalcenergy", &ir->nstcalcenergy);
730 if (ir->epc != epcNO)
732 check_nst_param(fplog, cr, "nstlist", ir->nstlist,
733 "nstpcouple", &ir->nstpcouple);
735 check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
736 "nstenergy", &ir->nstenergy);
737 check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
738 "nstlog", &ir->nstlog);
739 if (ir->efep != efepNO)
741 check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
742 "nstdhdl", &ir->fepvals->nstdhdl);
746 if (EI_VV(ir->eI) && IR_TWINRANGE(*ir) && ir->nstlist > 1)
748 gmx_fatal(FARGS, "Twin-range multiple time stepping does not work with integrator %s.", ei_names[ir->eI]);
752 void rerun_parallel_comm(t_commrec *cr, t_trxframe *fr,
753 gmx_bool *bNotLastFrame)
758 bAlloc = (fr->natoms == 0);
760 if (MASTER(cr) && !*bNotLastFrame)
766 gmx_bcast(sizeof(*fr), fr, cr);
770 *bNotLastFrame = (fr->natoms >= 0);
772 if (*bNotLastFrame && PARTDECOMP(cr))
774 /* x and v are the only variable size quantities stored in trr
775 * that are required for rerun (f is not needed).
779 snew(fr->x, fr->natoms);
780 snew(fr->v, fr->natoms);
784 gmx_bcast(fr->natoms*sizeof(fr->x[0]), fr->x[0], cr);
788 gmx_bcast(fr->natoms*sizeof(fr->v[0]), fr->v[0], cr);