src/mdlib/md_support.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  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.
  15
  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.
  20  *
  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.
  27  *
  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.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include "typedefs.h"
  41 #include "string2.h"
  42 #include "smalloc.h"
  43 #include "mdrun.h"
  44 #include "domdec.h"
  45 #include "mtop_util.h"
  46 #include "gmx_wallcycle.h"
  47 #include "vcm.h"
  48 #include "nrnb.h"
  49 #include "md_logging.h"
  50 #include "md_support.h"
  51
  52 /* Is the signal in one simulation independent of other simulations? */
  53 gmx_bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
  54
  55 /* check which of the multisim simulations has the shortest number of
  56    steps and return that number of nsteps */
  57 gmx_large_int_t get_multisim_nsteps(const t_commrec *cr,
  58                                     gmx_large_int_t nsteps)
  59 {
  60     gmx_large_int_t steps_out;
  61
  62     if MASTER(cr)
  63     {
  64         gmx_large_int_t *buf;
  65         int s;
  66
  67         snew(buf,cr->ms->nsim);
  68
  69         buf[cr->ms->sim] = nsteps;
  70         gmx_sumli_sim(cr->ms->nsim, buf, cr->ms);
  71
  72         steps_out=-1;
  73         for(s=0; s<cr->ms->nsim; s++)
  74         {
  75             /* find the smallest positive number */
  76             if (buf[s]>= 0 && ((steps_out < 0) || (buf[s]<steps_out)) )
  77             {
  78                 steps_out=buf[s];
  79             }
  80         }
  81         sfree(buf);
  82
  83         /* if we're the limiting simulation, don't do anything */
  84         if (steps_out>=0 && steps_out<nsteps)
  85         {
  86             char strbuf[255];
  87             snprintf(strbuf, 255, "Will stop simulation %%d after %s steps (another simulation will end then).\n", gmx_large_int_pfmt);
  88             fprintf(stderr, strbuf, cr->ms->sim, steps_out);
  89         }
  90     }
  91     /* broadcast to non-masters */
  92     gmx_bcast(sizeof(gmx_large_int_t), &steps_out, cr);
  93     return steps_out;
  94 }
  95
  96 int multisim_min(const gmx_multisim_t *ms,int nmin,int n)
  97 {
  98     int  *buf;
  99     gmx_bool bPos,bEqual;
 100     int  s,d;
 101
 102     snew(buf,ms->nsim);
 103     buf[ms->sim] = n;
 104     gmx_sumi_sim(ms->nsim,buf,ms);
 105     bPos   = TRUE;
 106     bEqual = TRUE;
 107     for(s=0; s<ms->nsim; s++)
 108     {
 109         bPos   = bPos   && (buf[s] > 0);
 110         bEqual = bEqual && (buf[s] == buf[0]);
 111     }
 112     if (bPos)
 113     {
 114         if (bEqual)
 115         {
 116             nmin = min(nmin,buf[0]);
 117         }
 118         else
 119         {
 120             /* Find the least common multiple */
 121             for(d=2; d<nmin; d++)
 122             {
 123                 s = 0;
 124                 while (s < ms->nsim && d % buf[s] == 0)
 125                 {
 126                     s++;
 127                 }
 128                 if (s == ms->nsim)
 129                 {
 130                     /* We found the LCM and it is less than nmin */
 131                     nmin = d;
 132                     break;
 133                 }
 134             }
 135         }
 136     }
 137     sfree(buf);
 138
 139     return nmin;
 140 }
 141
 142 int multisim_nstsimsync(const t_commrec *cr,
 143                         const t_inputrec *ir,int repl_ex_nst)
 144 {
 145     int nmin;
 146
 147     if (MASTER(cr))
 148     {
 149         nmin = INT_MAX;
 150         nmin = multisim_min(cr->ms,nmin,ir->nstlist);
 151         nmin = multisim_min(cr->ms,nmin,ir->nstcalcenergy);
 152         nmin = multisim_min(cr->ms,nmin,repl_ex_nst);
 153         if (nmin == INT_MAX)
 154         {
 155             gmx_fatal(FARGS,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
 156         }
 157         /* Avoid inter-simulation communication at every (second) step */
 158         if (nmin <= 2)
 159         {
 160             nmin = 10;
 161         }
 162     }
 163
 164     gmx_bcast(sizeof(int),&nmin,cr);
 165
 166     return nmin;
 167 }
 168
 169 void init_global_signals(globsig_t *gs,const t_commrec *cr,
 170                          const t_inputrec *ir,int repl_ex_nst)
 171 {
 172     int i;
 173
 174     if (MULTISIM(cr))
 175     {
 176         gs->nstms = multisim_nstsimsync(cr,ir,repl_ex_nst);
 177         if (debug)
 178         {
 179             fprintf(debug,"Syncing simulations for checkpointing and termination every %d steps\n",gs->nstms);
 180         }
 181     }
 182     else
 183     {
 184         gs->nstms = 1;
 185     }
 186
 187     for(i=0; i<eglsNR; i++)
 188     {
 189         gs->sig[i] = 0;
 190         gs->set[i] = 0;
 191     }
 192 }
 193
 194 void copy_coupling_state(t_state *statea,t_state *stateb,
 195                          gmx_ekindata_t *ekinda,gmx_ekindata_t *ekindb, t_grpopts* opts)
 196 {
 197
 198     /* MRS note -- might be able to get rid of some of the arguments.  Look over it when it's all debugged */
 199
 200     int i,j,nc;
 201
 202     /* Make sure we have enough space for x and v */
 203     if (statea->nalloc > stateb->nalloc)
 204     {
 205         stateb->nalloc = statea->nalloc;
 206         srenew(stateb->x,stateb->nalloc);
 207         srenew(stateb->v,stateb->nalloc);
 208     }
 209
 210     stateb->natoms     = statea->natoms;
 211     stateb->ngtc       = statea->ngtc;
 212     stateb->nnhpres    = statea->nnhpres;
 213     stateb->veta       = statea->veta;
 214     if (ekinda)
 215     {
 216         copy_mat(ekinda->ekin,ekindb->ekin);
 217         for (i=0; i<stateb->ngtc; i++)
 218         {
 219             ekindb->tcstat[i].T = ekinda->tcstat[i].T;
 220             ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
 221             copy_mat(ekinda->tcstat[i].ekinh,ekindb->tcstat[i].ekinh);
 222             copy_mat(ekinda->tcstat[i].ekinf,ekindb->tcstat[i].ekinf);
 223             ekindb->tcstat[i].ekinscalef_nhc =  ekinda->tcstat[i].ekinscalef_nhc;
 224             ekindb->tcstat[i].ekinscaleh_nhc =  ekinda->tcstat[i].ekinscaleh_nhc;
 225             ekindb->tcstat[i].vscale_nhc =  ekinda->tcstat[i].vscale_nhc;
 226         }
 227     }
 228     copy_rvecn(statea->x,stateb->x,0,stateb->natoms);
 229     copy_rvecn(statea->v,stateb->v,0,stateb->natoms);
 230     copy_mat(statea->box,stateb->box);
 231     copy_mat(statea->box_rel,stateb->box_rel);
 232     copy_mat(statea->boxv,stateb->boxv);
 233
 234     for (i = 0; i<stateb->ngtc; i++)
 235     {
 236         nc = i*opts->nhchainlength;
 237         for (j=0; j<opts->nhchainlength; j++)
 238         {
 239             stateb->nosehoover_xi[nc+j]  = statea->nosehoover_xi[nc+j];
 240             stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
 241         }
 242     }
 243     if (stateb->nhpres_xi != NULL)
 244     {
 245         for (i = 0; i<stateb->nnhpres; i++)
 246         {
 247             nc = i*opts->nhchainlength;
 248             for (j=0; j<opts->nhchainlength; j++)
 249             {
 250                 stateb->nhpres_xi[nc+j]  = statea->nhpres_xi[nc+j];
 251                 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
 252             }
 253         }
 254     }
 255 }
 256
 257 real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
 258 {
 259     real quantity = 0;
 260     switch (ir->etc)
 261     {
 262     case etcNO:
 263         break;
 264     case etcBERENDSEN:
 265         break;
 266     case etcNOSEHOOVER:
 267         quantity = NPT_energy(ir,state,MassQ);
 268         break;
 269     case etcVRESCALE:
 270         quantity = vrescale_energy(&(ir->opts),state->therm_integral);
 271         break;
 272     default:
 273         break;
 274     }
 275     return quantity;
 276 }
 277
 278 void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
 279                      t_forcerec *fr, gmx_ekindata_t *ekind,
 280                      t_state *state, t_state *state_global, t_mdatoms *mdatoms,
 281                      t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
 282                      gmx_enerdata_t *enerd,tensor force_vir, tensor shake_vir, tensor total_vir,
 283                      tensor pres, rvec mu_tot, gmx_constr_t constr,
 284                      globsig_t *gs,gmx_bool bInterSimGS,
 285                      matrix box, gmx_mtop_t *top_global, real *pcurr,
 286                      int natoms, gmx_bool *bSumEkinhOld, int flags)
 287 {
 288     int  i,gsi;
 289     real gs_buf[eglsNR];
 290     tensor corr_vir,corr_pres,shakeall_vir;
 291     gmx_bool bEner,bPres,bTemp, bVV;
 292     gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
 293         bFirstIterate,bReadEkin,bEkinAveVel,bScaleEkin, bConstrain;
 294     real ekin,temp,prescorr,enercorr,dvdlcorr;
 295
 296     /* translate CGLO flags to gmx_booleans */
 297     bRerunMD = flags & CGLO_RERUNMD;
 298     bStopCM = flags & CGLO_STOPCM;
 299     bGStat = flags & CGLO_GSTAT;
 300
 301     bReadEkin = (flags & CGLO_READEKIN);
 302     bScaleEkin = (flags & CGLO_SCALEEKIN);
 303     bEner = flags & CGLO_ENERGY;
 304     bTemp = flags & CGLO_TEMPERATURE;
 305     bPres  = (flags & CGLO_PRESSURE);
 306     bConstrain = (flags & CGLO_CONSTRAINT);
 307     bIterate = (flags & CGLO_ITERATE);
 308     bFirstIterate = (flags & CGLO_FIRSTITERATE);
 309
 310     /* we calculate a full state kinetic energy either with full-step velocity verlet
 311        or half step where we need the pressure */
 312
 313     bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && bPres) || bReadEkin);
 314
 315     /* in initalization, it sums the shake virial in vv, and to
 316        sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
 317
 318     /* ########## Kinetic energy  ############## */
 319
 320     if (bTemp)
 321     {
 322         /* Non-equilibrium MD: this is parallellized, but only does communication
 323          * when there really is NEMD.
 324          */
 325
 326         if (PAR(cr) && (ekind->bNEMD))
 327         {
 328             accumulate_u(cr,&(ir->opts),ekind);
 329         }
 330         debug_gmx();
 331         if (bReadEkin)
 332         {
 333             restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
 334         }
 335         else
 336         {
 337
 338             calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
 339         }
 340
 341         debug_gmx();
 342     }
 343
 344     /* Calculate center of mass velocity if necessary, also parallellized */
 345     if (bStopCM)
 346     {
 347         calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
 348                      state->x,state->v,vcm);
 349     }
 350
 351     if (bTemp || bStopCM || bPres || bEner || bConstrain)
 352     {
 353         if (!bGStat)
 354         {
 355             /* We will not sum ekinh_old,
 356              * so signal that we still have to do it.
 357              */
 358             *bSumEkinhOld = TRUE;
 359
 360         }
 361         else
 362         {
 363             if (gs != NULL)
 364             {
 365                 for(i=0; i<eglsNR; i++)
 366                 {
 367                     gs_buf[i] = gs->sig[i];
 368                 }
 369             }
 370             if (PAR(cr))
 371             {
 372                 wallcycle_start(wcycle,ewcMoveE);
 373                 GMX_MPE_LOG(ev_global_stat_start);
 374                 global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
 375                             ir,ekind,constr,bStopCM ? vcm : NULL,
 376                             gs != NULL ? eglsNR : 0,gs_buf,
 377                             top_global,state,
 378                             *bSumEkinhOld,flags);
 379                 GMX_MPE_LOG(ev_global_stat_finish);
 380                 wallcycle_stop(wcycle,ewcMoveE);
 381             }
 382             if (gs != NULL)
 383             {
 384                 if (MULTISIM(cr) && bInterSimGS)
 385                 {
 386                     if (MASTER(cr))
 387                     {
 388                         /* Communicate the signals between the simulations */
 389                         gmx_sum_sim(eglsNR,gs_buf,cr->ms);
 390                     }
 391                     /* Communicate the signals form the master to the others */
 392                     gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
 393                 }
 394                 for(i=0; i<eglsNR; i++)
 395                 {
 396                     if (bInterSimGS || gs_simlocal[i])
 397                     {
 398                         /* Set the communicated signal only when it is non-zero,
 399                          * since signals might not be processed at each MD step.
 400                          */
 401                         gsi = (gs_buf[i] >= 0 ?
 402                                (int)(gs_buf[i] + 0.5) :
 403                                (int)(gs_buf[i] - 0.5));
 404                         if (gsi != 0)
 405                         {
 406                             gs->set[i] = gsi;
 407                         }
 408                         /* Turn off the local signal */
 409                         gs->sig[i] = 0;
 410                     }
 411                 }
 412             }
 413             *bSumEkinhOld = FALSE;
 414         }
 415     }
 416
 417     if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
 418     {
 419         correct_ekin(debug,
 420                      mdatoms->start,mdatoms->start+mdatoms->homenr,
 421                      state->v,vcm->group_p[0],
 422                      mdatoms->massT,mdatoms->tmass,ekind->ekin);
 423     }
 424
 425     /* Do center of mass motion removal */
 426     if (bStopCM)
 427     {
 428         check_cm_grp(fplog,vcm,ir,1);
 429         do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
 430                       state->x,state->v,vcm);
 431         inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
 432     }
 433
 434     if (bEner)
 435     {
 436         /* Calculate the amplitude of the cosine velocity profile */
 437         ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
 438     }
 439
 440     if (bTemp)
 441     {
 442         /* Sum the kinetic energies of the groups & calc temp */
 443         /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
 444         /* three maincase:  VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
 445            Leap with AveVel is not supported; it's not clear that it will actually work.
 446            bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
 447            If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
 448            bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
 449            If FALSE, we go ahead and erase over it.
 450         */
 451         enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,&(enerd->term[F_DKDL]),
 452                                        bEkinAveVel,bIterate,bScaleEkin);
 453
 454         enerd->term[F_EKIN] = trace(ekind->ekin);
 455     }
 456
 457     /* ##########  Long range energy information ###### */
 458
 459     if (bEner || bPres || bConstrain)
 460     {
 461         calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda[efptVDW],
 462                       corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
 463     }
 464
 465     if (bEner && bFirstIterate)
 466     {
 467         enerd->term[F_DISPCORR] = enercorr;
 468         enerd->term[F_EPOT] += enercorr;
 469         enerd->term[F_DVDL_VDW] += dvdlcorr;
 470     }
 471
 472     /* ########## Now pressure ############## */
 473     if (bPres || bConstrain)
 474     {
 475
 476         m_add(force_vir,shake_vir,total_vir);
 477
 478         /* Calculate pressure and apply LR correction if PPPM is used.
 479          * Use the box from last timestep since we already called update().
 480          */
 481
 482         enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres);
 483
 484         /* Calculate long range corrections to pressure and energy */
 485         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
 486            and computes enerd->term[F_DISPCORR].  Also modifies the
 487            total_vir and pres tesors */
 488
 489         m_add(total_vir,corr_vir,total_vir);
 490         m_add(pres,corr_pres,pres);
 491         enerd->term[F_PDISPCORR] = prescorr;
 492         enerd->term[F_PRES] += prescorr;
 493         *pcurr = enerd->term[F_PRES];
 494         /* calculate temperature using virial */
 495         enerd->term[F_VTEMP] = calc_temp(trace(total_vir),ir->opts.nrdf[0]);
 496
 497     }
 498 }
 499
 500 void check_nst_param(FILE *fplog,t_commrec *cr,
 501                      const char *desc_nst,int nst,
 502                      const char *desc_p,int *p)
 503 {
 504     if (*p > 0 && *p % nst != 0)
 505     {
 506         /* Round up to the next multiple of nst */
 507         *p = ((*p)/nst + 1)*nst;
 508         md_print_warn(cr,fplog,
 509                       "NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
 510     }
 511 }
 512
 513 void set_current_lambdas(gmx_large_int_t step, t_lambda *fepvals, gmx_bool bRerunMD,
 514                          t_trxframe *rerun_fr,t_state *state_global, t_state *state, double lam0[])
 515 /* find the current lambdas.  If rerunning, we either read in a state, or a lambda value,
 516    requiring different logic. */
 517 {
 518     real frac;
 519     int i,fep_state=0;
 520     if (bRerunMD)
 521     {
 522         if (rerun_fr->bLambda)
 523         {
 524             if (fepvals->delta_lambda!=0)
 525             {
 526                 state_global->lambda[efptFEP] = rerun_fr->lambda;
 527                 for (i=0;i<efptNR;i++)
 528                 {
 529                     if (i!= efptFEP)
 530                     {
 531                         state->lambda[i] = state_global->lambda[i];
 532                     }
 533                 }
 534             }
 535             else
 536             {
 537                 /* find out between which two value of lambda we should be */
 538                 frac = (step*fepvals->delta_lambda);
 539                 fep_state = floor(frac*fepvals->n_lambda);
 540                 /* interpolate between this state and the next */
 541                 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 542                 frac = (frac*fepvals->n_lambda)-fep_state;
 543                 for (i=0;i<efptNR;i++)
 544                 {
 545                     state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
 546                         frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
 547                 }
 548             }
 549         }
 550         else if (rerun_fr->bFepState)
 551         {
 552             state_global->fep_state = rerun_fr->fep_state;
 553             for (i=0;i<efptNR;i++)
 554             {
 555                 state_global->lambda[i] = fepvals->all_lambda[i][fep_state];
 556             }
 557         }
 558     }
 559     else
 560     {
 561         if (fepvals->delta_lambda!=0)
 562         {
 563             /* find out between which two value of lambda we should be */
 564             frac = (step*fepvals->delta_lambda);
 565             if (fepvals->n_lambda > 0)
 566             {
 567                 fep_state = floor(frac*fepvals->n_lambda);
 568                 /* interpolate between this state and the next */
 569                 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 570                 frac = (frac*fepvals->n_lambda)-fep_state;
 571                 for (i=0;i<efptNR;i++)
 572                 {
 573                     state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
 574                         frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
 575                 }
 576             }
 577             else
 578             {
 579                 for (i=0;i<efptNR;i++)
 580                 {
 581                     state_global->lambda[i] = lam0[i] + frac;
 582                 }
 583             }
 584         }
 585     }
 586     for (i=0;i<efptNR;i++)
 587     {
 588         state->lambda[i] = state_global->lambda[i];
 589     }
 590 }
 591
 592 static void min_zero(int *n,int i)
 593 {
 594     if (i > 0 && (*n == 0 || i < *n))
 595     {
 596         *n = i;
 597     }
 598 }
 599
 600 static int lcd4(int i1,int i2,int i3,int i4)
 601 {
 602     int nst;
 603
 604     nst = 0;
 605     min_zero(&nst,i1);
 606     min_zero(&nst,i2);
 607     min_zero(&nst,i3);
 608     min_zero(&nst,i4);
 609     if (nst == 0)
 610     {
 611         gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
 612     }
 613
 614     while (nst > 1 && ((i1 > 0 && i1 % nst != 0)  ||
 615                        (i2 > 0 && i2 % nst != 0)  ||
 616                        (i3 > 0 && i3 % nst != 0)  ||
 617                        (i4 > 0 && i4 % nst != 0)))
 618     {
 619         nst--;
 620     }
 621
 622     return nst;
 623 }
 624
 625 int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
 626                         int nstglobalcomm,t_inputrec *ir)
 627 {
 628     if (!EI_DYNAMICS(ir->eI))
 629     {
 630         nstglobalcomm = 1;
 631     }
 632
 633     if (nstglobalcomm == -1)
 634     {
 635         if (!(ir->nstcalcenergy > 0 ||
 636               ir->nstlist > 0 ||
 637               ir->etc != etcNO ||
 638               ir->epc != epcNO))
 639         {
 640             nstglobalcomm = 10;
 641             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
 642             {
 643                 nstglobalcomm = ir->nstenergy;
 644             }
 645         }
 646         else
 647         {
 648             /* Ensure that we do timely global communication for
 649              * (possibly) each of the four following options.
 650              */
 651             nstglobalcomm = lcd4(ir->nstcalcenergy,
 652                                  ir->nstlist,
 653                                  ir->etc != etcNO ? ir->nsttcouple : 0,
 654                                  ir->epc != epcNO ? ir->nstpcouple : 0);
 655         }
 656     }
 657     else
 658     {
 659         if (ir->nstlist > 0 &&
 660             nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
 661         {
 662             nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
 663             md_print_warn(cr,fplog,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
 664         }
 665         if (ir->nstcalcenergy > 0)
 666         {
 667             check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 668                             "nstcalcenergy",&ir->nstcalcenergy);
 669         }
 670         if (ir->etc != etcNO && ir->nsttcouple > 0)
 671         {
 672             check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 673                             "nsttcouple",&ir->nsttcouple);
 674         }
 675         if (ir->epc != epcNO && ir->nstpcouple > 0)
 676         {
 677             check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 678                             "nstpcouple",&ir->nstpcouple);
 679         }
 680
 681         check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 682                         "nstenergy",&ir->nstenergy);
 683
 684         check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
 685                         "nstlog",&ir->nstlog);
 686     }
 687
 688     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
 689     {
 690         md_print_warn(cr,fplog,"WARNING: Changing nstcomm from %d to %d\n",
 691                       ir->nstcomm,nstglobalcomm);
 692         ir->nstcomm = nstglobalcomm;
 693     }
 694
 695     return nstglobalcomm;
 696 }
 697
 698 void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
 699                                t_inputrec *ir,gmx_mtop_t *mtop)
 700 {
 701     /* Check required for old tpx files */
 702     if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
 703         ir->nstcalcenergy % ir->nstlist != 0)
 704     {
 705         md_print_warn(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies\n");
 706
 707         if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
 708             gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
 709             ir->eConstrAlg == econtSHAKE)
 710         {
 711             md_print_warn(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect\n");
 712             if (ir->epc != epcNO)
 713             {
 714                 gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
 715             }
 716         }
 717         check_nst_param(fplog,cr,"nstlist",ir->nstlist,
 718                         "nstcalcenergy",&ir->nstcalcenergy);
 719         if (ir->epc != epcNO)
 720         {
 721             check_nst_param(fplog,cr,"nstlist",ir->nstlist,
 722                             "nstpcouple",&ir->nstpcouple);
 723         }
 724         check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
 725                         "nstenergy",&ir->nstenergy);
 726         check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
 727                         "nstlog",&ir->nstlog);
 728         if (ir->efep != efepNO)
 729         {
 730             check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
 731                             "nstdhdl",&ir->fepvals->nstdhdl);
 732         }
 733     }
 734 }
 735
 736 void rerun_parallel_comm(t_commrec *cr,t_trxframe *fr,
 737                          gmx_bool *bNotLastFrame)
 738 {
 739     gmx_bool bAlloc;
 740     rvec *xp,*vp;
 741
 742     bAlloc = (fr->natoms == 0);
 743
 744     if (MASTER(cr) && !*bNotLastFrame)
 745     {
 746         fr->natoms = -1;
 747     }
 748     xp = fr->x;
 749     vp = fr->v;
 750     gmx_bcast(sizeof(*fr),fr,cr);
 751     fr->x = xp;
 752     fr->v = vp;
 753
 754     *bNotLastFrame = (fr->natoms >= 0);
 755
 756     if (*bNotLastFrame && PARTDECOMP(cr))
 757     {
 758         /* x and v are the only variable size quantities stored in trr
 759          * that are required for rerun (f is not needed).
 760          */
 761         if (bAlloc)
 762         {
 763             snew(fr->x,fr->natoms);
 764             snew(fr->v,fr->natoms);
 765         }
 766         if (fr->bX)
 767         {
 768             gmx_bcast(fr->natoms*sizeof(fr->x[0]),fr->x[0],cr);
 769         }
 770         if (fr->bV)
 771         {
 772             gmx_bcast(fr->natoms*sizeof(fr->v[0]),fr->v[0],cr);
 773         }
 774     }
 775 }