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"
54 /* Is the signal in one simulation independent of other simulations? */
55 gmx_bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
57 /* check which of the multisim simulations has the shortest number of
58 steps and return that number of nsteps */
59 gmx_large_int_t get_multisim_nsteps(const t_commrec *cr,
60 gmx_large_int_t nsteps)
62 gmx_large_int_t steps_out;
69 snew(buf,cr->ms->nsim);
71 buf[cr->ms->sim] = nsteps;
72 gmx_sumli_sim(cr->ms->nsim, buf, cr->ms);
75 for(s=0; s<cr->ms->nsim; s++)
77 /* find the smallest positive number */
78 if (buf[s]>= 0 && ((steps_out < 0) || (buf[s]<steps_out)) )
85 /* if we're the limiting simulation, don't do anything */
86 if (steps_out>=0 && steps_out<nsteps)
89 snprintf(strbuf, 255, "Will stop simulation %%d after %s steps (another simulation will end then).\n", gmx_large_int_pfmt);
90 fprintf(stderr, strbuf, cr->ms->sim, steps_out);
93 /* broadcast to non-masters */
94 gmx_bcast(sizeof(gmx_large_int_t), &steps_out, cr);
98 int multisim_min(const gmx_multisim_t *ms,int nmin,int n)
101 gmx_bool bPos,bEqual;
106 gmx_sumi_sim(ms->nsim,buf,ms);
109 for(s=0; s<ms->nsim; s++)
111 bPos = bPos && (buf[s] > 0);
112 bEqual = bEqual && (buf[s] == buf[0]);
118 nmin = min(nmin,buf[0]);
122 /* Find the least common multiple */
123 for(d=2; d<nmin; d++)
126 while (s < ms->nsim && d % buf[s] == 0)
132 /* We found the LCM and it is less than nmin */
144 int multisim_nstsimsync(const t_commrec *cr,
145 const t_inputrec *ir,int repl_ex_nst)
152 nmin = multisim_min(cr->ms,nmin,ir->nstlist);
153 nmin = multisim_min(cr->ms,nmin,ir->nstcalcenergy);
154 nmin = multisim_min(cr->ms,nmin,repl_ex_nst);
157 gmx_fatal(FARGS,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
159 /* Avoid inter-simulation communication at every (second) step */
166 gmx_bcast(sizeof(int),&nmin,cr);
171 void init_global_signals(globsig_t *gs,const t_commrec *cr,
172 const t_inputrec *ir,int repl_ex_nst)
178 gs->nstms = multisim_nstsimsync(cr,ir,repl_ex_nst);
181 fprintf(debug,"Syncing simulations for checkpointing and termination every %d steps\n",gs->nstms);
189 for(i=0; i<eglsNR; i++)
196 void copy_coupling_state(t_state *statea,t_state *stateb,
197 gmx_ekindata_t *ekinda,gmx_ekindata_t *ekindb, t_grpopts* opts)
200 /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
204 /* Make sure we have enough space for x and v */
205 if (statea->nalloc > stateb->nalloc)
207 stateb->nalloc = statea->nalloc;
208 srenew(stateb->x,stateb->nalloc);
209 srenew(stateb->v,stateb->nalloc);
212 stateb->natoms = statea->natoms;
213 stateb->ngtc = statea->ngtc;
214 stateb->nnhpres = statea->nnhpres;
215 stateb->veta = statea->veta;
218 copy_mat(ekinda->ekin,ekindb->ekin);
219 for (i=0; i<stateb->ngtc; i++)
221 ekindb->tcstat[i].T = ekinda->tcstat[i].T;
222 ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
223 copy_mat(ekinda->tcstat[i].ekinh,ekindb->tcstat[i].ekinh);
224 copy_mat(ekinda->tcstat[i].ekinf,ekindb->tcstat[i].ekinf);
225 ekindb->tcstat[i].ekinscalef_nhc = ekinda->tcstat[i].ekinscalef_nhc;
226 ekindb->tcstat[i].ekinscaleh_nhc = ekinda->tcstat[i].ekinscaleh_nhc;
227 ekindb->tcstat[i].vscale_nhc = ekinda->tcstat[i].vscale_nhc;
230 copy_rvecn(statea->x,stateb->x,0,stateb->natoms);
231 copy_rvecn(statea->v,stateb->v,0,stateb->natoms);
232 copy_mat(statea->box,stateb->box);
233 copy_mat(statea->box_rel,stateb->box_rel);
234 copy_mat(statea->boxv,stateb->boxv);
236 for (i = 0; i<stateb->ngtc; i++)
238 nc = i*opts->nhchainlength;
239 for (j=0; j<opts->nhchainlength; j++)
241 stateb->nosehoover_xi[nc+j] = statea->nosehoover_xi[nc+j];
242 stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
245 if (stateb->nhpres_xi != NULL)
247 for (i = 0; i<stateb->nnhpres; i++)
249 nc = i*opts->nhchainlength;
250 for (j=0; j<opts->nhchainlength; j++)
252 stateb->nhpres_xi[nc+j] = statea->nhpres_xi[nc+j];
253 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
259 real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
269 quantity = NPT_energy(ir,state,MassQ);
272 quantity = vrescale_energy(&(ir->opts),state->therm_integral);
280 void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
281 t_forcerec *fr, gmx_ekindata_t *ekind,
282 t_state *state, t_state *state_global, t_mdatoms *mdatoms,
283 t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
284 gmx_enerdata_t *enerd,tensor force_vir, tensor shake_vir, tensor total_vir,
285 tensor pres, rvec mu_tot, gmx_constr_t constr,
286 globsig_t *gs,gmx_bool bInterSimGS,
287 matrix box, gmx_mtop_t *top_global, real *pcurr,
288 int natoms, gmx_bool *bSumEkinhOld, int flags)
292 tensor corr_vir,corr_pres,shakeall_vir;
293 gmx_bool bEner,bPres,bTemp, bVV;
294 gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
295 bFirstIterate,bReadEkin,bEkinAveVel,bScaleEkin, bConstrain;
296 real ekin,temp,prescorr,enercorr,dvdlcorr;
298 /* translate CGLO flags to gmx_booleans */
299 bRerunMD = flags & CGLO_RERUNMD;
300 bStopCM = flags & CGLO_STOPCM;
301 bGStat = flags & CGLO_GSTAT;
303 bReadEkin = (flags & CGLO_READEKIN);
304 bScaleEkin = (flags & CGLO_SCALEEKIN);
305 bEner = flags & CGLO_ENERGY;
306 bTemp = flags & CGLO_TEMPERATURE;
307 bPres = (flags & CGLO_PRESSURE);
308 bConstrain = (flags & CGLO_CONSTRAINT);
309 bIterate = (flags & CGLO_ITERATE);
310 bFirstIterate = (flags & CGLO_FIRSTITERATE);
312 /* we calculate a full state kinetic energy either with full-step velocity verlet
313 or half step where we need the pressure */
315 bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && bPres) || bReadEkin);
317 /* in initalization, it sums the shake virial in vv, and to
318 sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
320 /* ########## Kinetic energy ############## */
324 /* Non-equilibrium MD: this is parallellized, but only does communication
325 * when there really is NEMD.
328 if (PAR(cr) && (ekind->bNEMD))
330 accumulate_u(cr,&(ir->opts),ekind);
335 restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
340 calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
346 /* Calculate center of mass velocity if necessary, also parallellized */
349 calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
350 state->x,state->v,vcm);
353 if (bTemp || bStopCM || bPres || bEner || bConstrain)
357 /* We will not sum ekinh_old,
358 * so signal that we still have to do it.
360 *bSumEkinhOld = TRUE;
367 for(i=0; i<eglsNR; i++)
369 gs_buf[i] = gs->sig[i];
374 wallcycle_start(wcycle,ewcMoveE);
375 GMX_MPE_LOG(ev_global_stat_start);
376 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
377 ir,ekind,constr,bStopCM ? vcm : NULL,
378 gs != NULL ? eglsNR : 0,gs_buf,
380 *bSumEkinhOld,flags);
381 GMX_MPE_LOG(ev_global_stat_finish);
382 wallcycle_stop(wcycle,ewcMoveE);
386 if (MULTISIM(cr) && bInterSimGS)
390 /* Communicate the signals between the simulations */
391 gmx_sum_sim(eglsNR,gs_buf,cr->ms);
393 /* Communicate the signals form the master to the others */
394 gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
396 for(i=0; i<eglsNR; i++)
398 if (bInterSimGS || gs_simlocal[i])
400 /* Set the communicated signal only when it is non-zero,
401 * since signals might not be processed at each MD step.
403 gsi = (gs_buf[i] >= 0 ?
404 (int)(gs_buf[i] + 0.5) :
405 (int)(gs_buf[i] - 0.5));
410 /* Turn off the local signal */
415 *bSumEkinhOld = FALSE;
419 if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
422 mdatoms->start,mdatoms->start+mdatoms->homenr,
423 state->v,vcm->group_p[0],
424 mdatoms->massT,mdatoms->tmass,ekind->ekin);
427 /* Do center of mass motion removal */
430 check_cm_grp(fplog,vcm,ir,1);
431 do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
432 state->x,state->v,vcm);
433 inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
438 /* Calculate the amplitude of the cosine velocity profile */
439 ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
444 /* Sum the kinetic energies of the groups & calc temp */
445 /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
446 /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
447 Leap with AveVel is not supported; it's not clear that it will actually work.
448 bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
449 If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
450 bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
451 If FALSE, we go ahead and erase over it.
453 enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,&(enerd->term[F_DKDL]),
454 bEkinAveVel,bIterate,bScaleEkin);
456 enerd->term[F_EKIN] = trace(ekind->ekin);
459 /* ########## Long range energy information ###### */
461 if (bEner || bPres || bConstrain)
463 calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda[efptVDW],
464 corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
467 if (bEner && bFirstIterate)
469 enerd->term[F_DISPCORR] = enercorr;
470 enerd->term[F_EPOT] += enercorr;
471 enerd->term[F_DVDL_VDW] += dvdlcorr;
474 /* ########## Now pressure ############## */
475 if (bPres || bConstrain)
478 m_add(force_vir,shake_vir,total_vir);
480 /* Calculate pressure and apply LR correction if PPPM is used.
481 * Use the box from last timestep since we already called update().
484 enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres);
486 /* Calculate long range corrections to pressure and energy */
487 /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
488 and computes enerd->term[F_DISPCORR]. Also modifies the
489 total_vir and pres tesors */
491 m_add(total_vir,corr_vir,total_vir);
492 m_add(pres,corr_pres,pres);
493 enerd->term[F_PDISPCORR] = prescorr;
494 enerd->term[F_PRES] += prescorr;
495 *pcurr = enerd->term[F_PRES];
499 void check_nst_param(FILE *fplog,t_commrec *cr,
500 const char *desc_nst,int nst,
501 const char *desc_p,int *p)
503 if (*p > 0 && *p % nst != 0)
505 /* Round up to the next multiple of nst */
506 *p = ((*p)/nst + 1)*nst;
507 md_print_warn(cr,fplog,
508 "NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
512 void set_current_lambdas(gmx_large_int_t step, t_lambda *fepvals, gmx_bool bRerunMD,
513 t_trxframe *rerun_fr,t_state *state_global, t_state *state, double lam0[])
514 /* find the current lambdas. If rerunning, we either read in a state, or a lambda value,
515 requiring different logic. */
521 if (rerun_fr->bLambda)
523 if (fepvals->delta_lambda!=0)
525 state_global->lambda[efptFEP] = rerun_fr->lambda;
526 for (i=0;i<efptNR;i++)
530 state->lambda[i] = state_global->lambda[i];
536 /* find out between which two value of lambda we should be */
537 frac = (step*fepvals->delta_lambda);
538 fep_state = floor(frac*fepvals->n_lambda);
539 /* interpolate between this state and the next */
540 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
541 frac = (frac*fepvals->n_lambda)-fep_state;
542 for (i=0;i<efptNR;i++)
544 state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
545 frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
549 else if (rerun_fr->bFepState)
551 state_global->fep_state = rerun_fr->fep_state;
552 for (i=0;i<efptNR;i++)
554 state_global->lambda[i] = fepvals->all_lambda[i][fep_state];
560 if (fepvals->delta_lambda!=0)
562 /* find out between which two value of lambda we should be */
563 frac = (step*fepvals->delta_lambda);
564 if (fepvals->n_lambda > 0)
566 fep_state = floor(frac*fepvals->n_lambda);
567 /* interpolate between this state and the next */
568 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
569 frac = (frac*fepvals->n_lambda)-fep_state;
570 for (i=0;i<efptNR;i++)
572 state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
573 frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
578 for (i=0;i<efptNR;i++)
580 state_global->lambda[i] = lam0[i] + frac;
585 for (i=0;i<efptNR;i++)
587 state->lambda[i] = state_global->lambda[i];
591 static void min_zero(int *n,int i)
593 if (i > 0 && (*n == 0 || i < *n))
599 static int lcd4(int i1,int i2,int i3,int i4)
610 gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
613 while (nst > 1 && ((i1 > 0 && i1 % nst != 0) ||
614 (i2 > 0 && i2 % nst != 0) ||
615 (i3 > 0 && i3 % nst != 0) ||
616 (i4 > 0 && i4 % nst != 0)))
624 int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
625 int nstglobalcomm,t_inputrec *ir)
627 if (!EI_DYNAMICS(ir->eI))
632 if (nstglobalcomm == -1)
634 if (!(ir->nstcalcenergy > 0 ||
640 if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
642 nstglobalcomm = ir->nstenergy;
647 /* Ensure that we do timely global communication for
648 * (possibly) each of the four following options.
650 nstglobalcomm = lcd4(ir->nstcalcenergy,
652 ir->etc != etcNO ? ir->nsttcouple : 0,
653 ir->epc != epcNO ? ir->nstpcouple : 0);
658 if (ir->nstlist > 0 &&
659 nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
661 nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
662 md_print_warn(cr,fplog,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
664 if (ir->nstcalcenergy > 0)
666 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
667 "nstcalcenergy",&ir->nstcalcenergy);
669 if (ir->etc != etcNO && ir->nsttcouple > 0)
671 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
672 "nsttcouple",&ir->nsttcouple);
674 if (ir->epc != epcNO && ir->nstpcouple > 0)
676 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
677 "nstpcouple",&ir->nstpcouple);
680 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
681 "nstenergy",&ir->nstenergy);
683 check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
684 "nstlog",&ir->nstlog);
687 if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
689 md_print_warn(cr,fplog,"WARNING: Changing nstcomm from %d to %d\n",
690 ir->nstcomm,nstglobalcomm);
691 ir->nstcomm = nstglobalcomm;
694 return nstglobalcomm;
697 void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
698 t_inputrec *ir,gmx_mtop_t *mtop)
700 /* Check required for old tpx files */
701 if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
702 ir->nstcalcenergy % ir->nstlist != 0)
704 md_print_warn(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies\n");
706 if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
707 gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
708 ir->eConstrAlg == econtSHAKE)
710 md_print_warn(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect\n");
711 if (ir->epc != epcNO)
713 gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
716 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
717 "nstcalcenergy",&ir->nstcalcenergy);
718 if (ir->epc != epcNO)
720 check_nst_param(fplog,cr,"nstlist",ir->nstlist,
721 "nstpcouple",&ir->nstpcouple);
723 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
724 "nstenergy",&ir->nstenergy);
725 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
726 "nstlog",&ir->nstlog);
727 if (ir->efep != efepNO)
729 check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
730 "nstdhdl",&ir->fepvals->nstdhdl);
735 void rerun_parallel_comm(t_commrec *cr,t_trxframe *fr,
736 gmx_bool *bNotLastFrame)
741 bAlloc = (fr->natoms == 0);
743 if (MASTER(cr) && !*bNotLastFrame)
749 gmx_bcast(sizeof(*fr),fr,cr);
753 *bNotLastFrame = (fr->natoms >= 0);
755 if (*bNotLastFrame && PARTDECOMP(cr))
757 /* x and v are the only variable size quantities stored in trr
758 * that are required for rerun (f is not needed).
762 snew(fr->x,fr->natoms);
763 snew(fr->v,fr->natoms);
767 gmx_bcast(fr->natoms*sizeof(fr->x[0]),fr->x[0],cr);
771 gmx_bcast(fr->natoms*sizeof(fr->v[0]),fr->v[0],cr);