src/gromacs/mdlib/md_support.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "config.h"
  38
  39 #include "typedefs.h"
  40 #include "mdrun.h"
  41 #include "domdec.h"
  42 #include "gromacs/topology/mtop_util.h"
  43 #include "vcm.h"
  44 #include "nrnb.h"
  45 #include "macros.h"
  46 #include "md_logging.h"
  47 #include "md_support.h"
  48 #include "names.h"
  49
  50 #include "gromacs/legacyheaders/types/commrec.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/timing/wallcycle.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "gromacs/utility/smalloc.h"
  55
  56 /* Is the signal in one simulation independent of other simulations? */
  57 gmx_bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
  58
  59 /* check which of the multisim simulations has the shortest number of
  60    steps and return that number of nsteps */
  61 gmx_int64_t get_multisim_nsteps(const t_commrec *cr,
  62                                 gmx_int64_t      nsteps)
  63 {
  64     gmx_int64_t steps_out;
  65
  66     if (MASTER(cr))
  67     {
  68         gmx_int64_t     *buf;
  69         int              s;
  70
  71         snew(buf, cr->ms->nsim);
  72
  73         buf[cr->ms->sim] = nsteps;
  74         gmx_sumli_sim(cr->ms->nsim, buf, cr->ms);
  75
  76         steps_out = -1;
  77         for (s = 0; s < cr->ms->nsim; s++)
  78         {
  79             /* find the smallest positive number */
  80             if (buf[s] >= 0 && ((steps_out < 0) || (buf[s] < steps_out)) )
  81             {
  82                 steps_out = buf[s];
  83             }
  84         }
  85         sfree(buf);
  86
  87         /* if we're the limiting simulation, don't do anything */
  88         if (steps_out >= 0 && steps_out < nsteps)
  89         {
  90             char strbuf[255];
  91             snprintf(strbuf, 255, "Will stop simulation %%d after %s steps (another simulation will end then).\n", "%"GMX_PRId64);
  92             fprintf(stderr, strbuf, cr->ms->sim, steps_out);
  93         }
  94     }
  95     /* broadcast to non-masters */
  96     gmx_bcast(sizeof(gmx_int64_t), &steps_out, cr);
  97     return steps_out;
  98 }
  99
 100 int multisim_min(const gmx_multisim_t *ms, int nmin, int n)
 101 {
 102     int     *buf;
 103     gmx_bool bPos, bEqual;
 104     int      s, d;
 105
 106     snew(buf, ms->nsim);
 107     buf[ms->sim] = n;
 108     gmx_sumi_sim(ms->nsim, buf, ms);
 109     bPos   = TRUE;
 110     bEqual = TRUE;
 111     for (s = 0; s < ms->nsim; s++)
 112     {
 113         bPos   = bPos   && (buf[s] > 0);
 114         bEqual = bEqual && (buf[s] == buf[0]);
 115     }
 116     if (bPos)
 117     {
 118         if (bEqual)
 119         {
 120             nmin = min(nmin, buf[0]);
 121         }
 122         else
 123         {
 124             /* Find the least common multiple */
 125             for (d = 2; d < nmin; d++)
 126             {
 127                 s = 0;
 128                 while (s < ms->nsim && d % buf[s] == 0)
 129                 {
 130                     s++;
 131                 }
 132                 if (s == ms->nsim)
 133                 {
 134                     /* We found the LCM and it is less than nmin */
 135                     nmin = d;
 136                     break;
 137                 }
 138             }
 139         }
 140     }
 141     sfree(buf);
 142
 143     return nmin;
 144 }
 145
 146 int multisim_nstsimsync(const t_commrec *cr,
 147                         const t_inputrec *ir, int repl_ex_nst)
 148 {
 149     int nmin;
 150
 151     if (MASTER(cr))
 152     {
 153         nmin = INT_MAX;
 154         nmin = multisim_min(cr->ms, nmin, ir->nstlist);
 155         nmin = multisim_min(cr->ms, nmin, ir->nstcalcenergy);
 156         nmin = multisim_min(cr->ms, nmin, repl_ex_nst);
 157         if (nmin == INT_MAX)
 158         {
 159             gmx_fatal(FARGS, "Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
 160         }
 161         /* Avoid inter-simulation communication at every (second) step */
 162         if (nmin <= 2)
 163         {
 164             nmin = 10;
 165         }
 166     }
 167
 168     gmx_bcast(sizeof(int), &nmin, cr);
 169
 170     return nmin;
 171 }
 172
 173 void init_global_signals(globsig_t *gs, const t_commrec *cr,
 174                          const t_inputrec *ir, int repl_ex_nst)
 175 {
 176     int i;
 177
 178     if (MULTISIM(cr))
 179     {
 180         gs->nstms = multisim_nstsimsync(cr, ir, repl_ex_nst);
 181         if (debug)
 182         {
 183             fprintf(debug, "Syncing simulations for checkpointing and termination every %d steps\n", gs->nstms);
 184         }
 185     }
 186     else
 187     {
 188         gs->nstms = 1;
 189     }
 190
 191     for (i = 0; i < eglsNR; i++)
 192     {
 193         gs->sig[i] = 0;
 194         gs->set[i] = 0;
 195     }
 196 }
 197
 198 void copy_coupling_state(t_state *statea, t_state *stateb,
 199                          gmx_ekindata_t *ekinda, gmx_ekindata_t *ekindb, t_grpopts* opts)
 200 {
 201
 202     /* MRS note -- might be able to get rid of some of the arguments.  Look over it when it's all debugged */
 203
 204     int i, j, nc;
 205
 206     /* Make sure we have enough space for x and v */
 207     if (statea->nalloc > stateb->nalloc)
 208     {
 209         stateb->nalloc = statea->nalloc;
 210         srenew(stateb->x, stateb->nalloc);
 211         srenew(stateb->v, stateb->nalloc);
 212     }
 213
 214     stateb->natoms     = statea->natoms;
 215     stateb->ngtc       = statea->ngtc;
 216     stateb->nnhpres    = statea->nnhpres;
 217     stateb->veta       = statea->veta;
 218     if (ekinda)
 219     {
 220         copy_mat(ekinda->ekin, ekindb->ekin);
 221         for (i = 0; i < stateb->ngtc; i++)
 222         {
 223             ekindb->tcstat[i].T  = ekinda->tcstat[i].T;
 224             ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
 225             copy_mat(ekinda->tcstat[i].ekinh, ekindb->tcstat[i].ekinh);
 226             copy_mat(ekinda->tcstat[i].ekinf, ekindb->tcstat[i].ekinf);
 227             ekindb->tcstat[i].ekinscalef_nhc =  ekinda->tcstat[i].ekinscalef_nhc;
 228             ekindb->tcstat[i].ekinscaleh_nhc =  ekinda->tcstat[i].ekinscaleh_nhc;
 229             ekindb->tcstat[i].vscale_nhc     =  ekinda->tcstat[i].vscale_nhc;
 230         }
 231     }
 232     copy_rvecn(statea->x, stateb->x, 0, stateb->natoms);
 233     copy_rvecn(statea->v, stateb->v, 0, stateb->natoms);
 234     copy_mat(statea->box, stateb->box);
 235     copy_mat(statea->box_rel, stateb->box_rel);
 236     copy_mat(statea->boxv, stateb->boxv);
 237
 238     for (i = 0; i < stateb->ngtc; i++)
 239     {
 240         nc = i*opts->nhchainlength;
 241         for (j = 0; j < opts->nhchainlength; j++)
 242         {
 243             stateb->nosehoover_xi[nc+j]  = statea->nosehoover_xi[nc+j];
 244             stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
 245         }
 246     }
 247     if (stateb->nhpres_xi != NULL)
 248     {
 249         for (i = 0; i < stateb->nnhpres; i++)
 250         {
 251             nc = i*opts->nhchainlength;
 252             for (j = 0; j < opts->nhchainlength; j++)
 253             {
 254                 stateb->nhpres_xi[nc+j]  = statea->nhpres_xi[nc+j];
 255                 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
 256             }
 257         }
 258     }
 259 }
 260
 261 real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
 262 {
 263     real quantity = 0;
 264     switch (ir->etc)
 265     {
 266         case etcNO:
 267             break;
 268         case etcBERENDSEN:
 269             break;
 270         case etcNOSEHOOVER:
 271             quantity = NPT_energy(ir, state, MassQ);
 272             break;
 273         case etcVRESCALE:
 274             quantity = vrescale_energy(&(ir->opts), state->therm_integral);
 275             break;
 276         default:
 277             break;
 278     }
 279     return quantity;
 280 }
 281
 282 void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
 283                      t_forcerec *fr, gmx_ekindata_t *ekind,
 284                      t_state *state, t_state *state_global, t_mdatoms *mdatoms,
 285                      t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
 286                      gmx_enerdata_t *enerd, tensor force_vir, tensor shake_vir, tensor total_vir,
 287                      tensor pres, rvec mu_tot, gmx_constr_t constr,
 288                      globsig_t *gs, gmx_bool bInterSimGS,
 289                      matrix box, gmx_mtop_t *top_global,
 290                      gmx_bool *bSumEkinhOld, int flags)
 291 {
 292     int      i, gsi;
 293     real     gs_buf[eglsNR];
 294     tensor   corr_vir, corr_pres;
 295     gmx_bool bEner, bPres, bTemp, bVV;
 296     gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
 297              bFirstIterate, bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
 298     real     ekin, temp, prescorr, enercorr, dvdlcorr, dvdl_ekin;
 299
 300     /* translate CGLO flags to gmx_booleans */
 301     bRerunMD = flags & CGLO_RERUNMD;
 302     bStopCM  = flags & CGLO_STOPCM;
 303     bGStat   = flags & CGLO_GSTAT;
 304
 305     bReadEkin     = (flags & CGLO_READEKIN);
 306     bScaleEkin    = (flags & CGLO_SCALEEKIN);
 307     bEner         = flags & CGLO_ENERGY;
 308     bTemp         = flags & CGLO_TEMPERATURE;
 309     bPres         = (flags & CGLO_PRESSURE);
 310     bConstrain    = (flags & CGLO_CONSTRAINT);
 311     bIterate      = (flags & CGLO_ITERATE);
 312     bFirstIterate = (flags & CGLO_FIRSTITERATE);
 313
 314     /* we calculate a full state kinetic energy either with full-step velocity verlet
 315        or half step where we need the pressure */
 316
 317     bEkinAveVel = (ir->eI == eiVV || (ir->eI == eiVVAK && bPres) || bReadEkin);
 318
 319     /* in initalization, it sums the shake virial in vv, and to
 320        sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
 321
 322     /* ########## Kinetic energy  ############## */
 323
 324     if (bTemp)
 325     {
 326         /* Non-equilibrium MD: this is parallellized, but only does communication
 327          * when there really is NEMD.
 328          */
 329
 330         if (PAR(cr) && (ekind->bNEMD))
 331         {
 332             accumulate_u(cr, &(ir->opts), ekind);
 333         }
 334         debug_gmx();
 335         if (bReadEkin)
 336         {
 337             restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 338         }
 339         else
 340         {
 341
 342             calc_ke_part(state, &(ir->opts), mdatoms, ekind, nrnb, bEkinAveVel, bIterate);
 343         }
 344
 345         debug_gmx();
 346     }
 347
 348     /* Calculate center of mass velocity if necessary, also parallellized */
 349     if (bStopCM)
 350     {
 351         calc_vcm_grp(0, mdatoms->homenr, mdatoms,
 352                      state->x, state->v, vcm);
 353     }
 354
 355     if (bTemp || bStopCM || bPres || bEner || bConstrain)
 356     {
 357         if (!bGStat)
 358         {
 359             /* We will not sum ekinh_old,
 360              * so signal that we still have to do it.
 361              */
 362             *bSumEkinhOld = TRUE;
 363
 364         }
 365         else
 366         {
 367             if (gs != NULL)
 368             {
 369                 for (i = 0; i < eglsNR; i++)
 370                 {
 371                     gs_buf[i] = gs->sig[i];
 372                 }
 373             }
 374             if (PAR(cr))
 375             {
 376                 wallcycle_start(wcycle, ewcMoveE);
 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,
 380                             top_global, state,
 381                             *bSumEkinhOld, flags);
 382                 wallcycle_stop(wcycle, ewcMoveE);
 383             }
 384             if (gs != NULL)
 385             {
 386                 if (MULTISIM(cr) && bInterSimGS)
 387                 {
 388                     if (MASTER(cr))
 389                     {
 390                         /* Communicate the signals between the simulations */
 391                         gmx_sum_sim(eglsNR, gs_buf, cr->ms);
 392                     }
 393                     /* Communicate the signals form the master to the others */
 394                     gmx_bcast(eglsNR*sizeof(gs_buf[0]), gs_buf, cr);
 395                 }
 396                 for (i = 0; i < eglsNR; i++)
 397                 {
 398                     if (bInterSimGS || gs_simlocal[i])
 399                     {
 400                         /* Set the communicated signal only when it is non-zero,
 401                          * since signals might not be processed at each MD step.
 402                          */
 403                         gsi = (gs_buf[i] >= 0 ?
 404                                (int)(gs_buf[i] + 0.5) :
 405                                (int)(gs_buf[i] - 0.5));
 406                         if (gsi != 0)
 407                         {
 408                             gs->set[i] = gsi;
 409                         }
 410                         /* Turn off the local signal */
 411                         gs->sig[i] = 0;
 412                     }
 413                 }
 414             }
 415             *bSumEkinhOld = FALSE;
 416         }
 417     }
 418
 419     if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
 420     {
 421         correct_ekin(debug,
 422                      0, mdatoms->homenr,
 423                      state->v, vcm->group_p[0],
 424                      mdatoms->massT, mdatoms->tmass, ekind->ekin);
 425     }
 426
 427     /* Do center of mass motion removal */
 428     if (bStopCM)
 429     {
 430         check_cm_grp(fplog, vcm, ir, 1);
 431         do_stopcm_grp(0, mdatoms->homenr, mdatoms->cVCM,
 432                       state->x, state->v, vcm);
 433         inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 434     }
 435
 436     if (bEner)
 437     {
 438         /* Calculate the amplitude of the cosine velocity profile */
 439         ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
 440     }
 441
 442     if (bTemp)
 443     {
 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.
 452          */
 453         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin,
 454                                        bEkinAveVel, bScaleEkin);
 455         enerd->dvdl_lin[efptMASS] = (double) dvdl_ekin;
 456
 457         enerd->term[F_EKIN] = trace(ekind->ekin);
 458     }
 459
 460     /* ##########  Long range energy information ###### */
 461
 462     if (bEner || bPres || bConstrain)
 463     {
 464         calc_dispcorr(fplog, ir, fr, 0, top_global->natoms, box, state->lambda[efptVDW],
 465                       corr_pres, corr_vir, &prescorr, &enercorr, &dvdlcorr);
 466     }
 467
 468     if (bEner && bFirstIterate)
 469     {
 470         enerd->term[F_DISPCORR]  = enercorr;
 471         enerd->term[F_EPOT]     += enercorr;
 472         enerd->term[F_DVDL_VDW] += dvdlcorr;
 473     }
 474
 475     /* ########## Now pressure ############## */
 476     if (bPres || bConstrain)
 477     {
 478
 479         m_add(force_vir, shake_vir, total_vir);
 480
 481         /* Calculate pressure and apply LR correction if PPPM is used.
 482          * Use the box from last timestep since we already called update().
 483          */
 484
 485         enerd->term[F_PRES] = calc_pres(fr->ePBC, ir->nwall, box, ekind->ekin, total_vir, pres);
 486
 487         /* Calculate long range corrections to pressure and energy */
 488         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
 489            and computes enerd->term[F_DISPCORR].  Also modifies the
 490            total_vir and pres tesors */
 491
 492         m_add(total_vir, corr_vir, total_vir);
 493         m_add(pres, corr_pres, pres);
 494         enerd->term[F_PDISPCORR] = prescorr;
 495         enerd->term[F_PRES]     += prescorr;
 496     }
 497 }
 498
 499 void check_nst_param(FILE *fplog, t_commrec *cr,
 500                      const char *desc_nst, int nst,
 501                      const char *desc_p, int *p)
 502 {
 503     if (*p > 0 && *p % nst != 0)
 504     {
 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);
 509     }
 510 }
 511
 512 void set_current_lambdas(gmx_int64_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. */
 516 {
 517     real frac;
 518     int  i, fep_state = 0;
 519     if (bRerunMD)
 520     {
 521         if (rerun_fr->bLambda)
 522         {
 523             if (fepvals->delta_lambda == 0)
 524             {
 525                 state_global->lambda[efptFEP] = rerun_fr->lambda;
 526                 for (i = 0; i < efptNR; i++)
 527                 {
 528                     if (i != efptFEP)
 529                     {
 530                         state->lambda[i] = state_global->lambda[i];
 531                     }
 532                 }
 533             }
 534             else
 535             {
 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++)
 543                 {
 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]);
 546                 }
 547             }
 548         }
 549         else if (rerun_fr->bFepState)
 550         {
 551             state_global->fep_state = rerun_fr->fep_state;
 552             for (i = 0; i < efptNR; i++)
 553             {
 554                 state_global->lambda[i] = fepvals->all_lambda[i][fep_state];
 555             }
 556         }
 557     }
 558     else
 559     {
 560         if (fepvals->delta_lambda != 0)
 561         {
 562             /* find out between which two value of lambda we should be */
 563             frac = (step*fepvals->delta_lambda);
 564             if (fepvals->n_lambda > 0)
 565             {
 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++)
 571                 {
 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]);
 574                 }
 575             }
 576             else
 577             {
 578                 for (i = 0; i < efptNR; i++)
 579                 {
 580                     state_global->lambda[i] = lam0[i] + frac;
 581                 }
 582             }
 583         }
 584         else
 585         {
 586             if (state->fep_state > 0)
 587             {
 588                 state_global->fep_state = state->fep_state; /* state->fep is the one updated by bExpanded */
 589                 for (i = 0; i < efptNR; i++)
 590                 {
 591                     state_global->lambda[i] = fepvals->all_lambda[i][state_global->fep_state];
 592                 }
 593             }
 594         }
 595     }
 596     for (i = 0; i < efptNR; i++)
 597     {
 598         state->lambda[i] = state_global->lambda[i];
 599     }
 600 }
 601
 602 static void min_zero(int *n, int i)
 603 {
 604     if (i > 0 && (*n == 0 || i < *n))
 605     {
 606         *n = i;
 607     }
 608 }
 609
 610 static int lcd4(int i1, int i2, int i3, int i4)
 611 {
 612     int nst;
 613
 614     nst = 0;
 615     min_zero(&nst, i1);
 616     min_zero(&nst, i2);
 617     min_zero(&nst, i3);
 618     min_zero(&nst, i4);
 619     if (nst == 0)
 620     {
 621         gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
 622     }
 623
 624     while (nst > 1 && ((i1 > 0 && i1 % nst != 0)  ||
 625                        (i2 > 0 && i2 % nst != 0)  ||
 626                        (i3 > 0 && i3 % nst != 0)  ||
 627                        (i4 > 0 && i4 % nst != 0)))
 628     {
 629         nst--;
 630     }
 631
 632     return nst;
 633 }
 634
 635 int check_nstglobalcomm(FILE *fplog, t_commrec *cr,
 636                         int nstglobalcomm, t_inputrec *ir)
 637 {
 638     if (!EI_DYNAMICS(ir->eI))
 639     {
 640         nstglobalcomm = 1;
 641     }
 642
 643     if (nstglobalcomm == -1)
 644     {
 645         if (!(ir->nstcalcenergy > 0 ||
 646               ir->nstlist > 0 ||
 647               ir->etc != etcNO ||
 648               ir->epc != epcNO))
 649         {
 650             nstglobalcomm = 10;
 651             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
 652             {
 653                 nstglobalcomm = ir->nstenergy;
 654             }
 655         }
 656         else
 657         {
 658             /* Ensure that we do timely global communication for
 659              * (possibly) each of the four following options.
 660              */
 661             nstglobalcomm = lcd4(ir->nstcalcenergy,
 662                                  ir->nstlist,
 663                                  ir->etc != etcNO ? ir->nsttcouple : 0,
 664                                  ir->epc != epcNO ? ir->nstpcouple : 0);
 665         }
 666     }
 667     else
 668     {
 669         if (ir->nstlist > 0 &&
 670             nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
 671         {
 672             nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
 673             md_print_warn(cr, fplog, "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n", nstglobalcomm);
 674         }
 675         if (ir->nstcalcenergy > 0)
 676         {
 677             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 678                             "nstcalcenergy", &ir->nstcalcenergy);
 679         }
 680         if (ir->etc != etcNO && ir->nsttcouple > 0)
 681         {
 682             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 683                             "nsttcouple", &ir->nsttcouple);
 684         }
 685         if (ir->epc != epcNO && ir->nstpcouple > 0)
 686         {
 687             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 688                             "nstpcouple", &ir->nstpcouple);
 689         }
 690
 691         check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 692                         "nstenergy", &ir->nstenergy);
 693
 694         check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 695                         "nstlog", &ir->nstlog);
 696     }
 697
 698     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
 699     {
 700         md_print_warn(cr, fplog, "WARNING: Changing nstcomm from %d to %d\n",
 701                       ir->nstcomm, nstglobalcomm);
 702         ir->nstcomm = nstglobalcomm;
 703     }
 704
 705     return nstglobalcomm;
 706 }
 707
 708 void check_ir_old_tpx_versions(t_commrec *cr, FILE *fplog,
 709                                t_inputrec *ir, gmx_mtop_t *mtop)
 710 {
 711     /* Check required for old tpx files */
 712     if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
 713         ir->nstcalcenergy % ir->nstlist != 0)
 714     {
 715         md_print_warn(cr, fplog, "Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies\n");
 716
 717         if (gmx_mtop_ftype_count(mtop, F_CONSTR) +
 718             gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0 &&
 719             ir->eConstrAlg == econtSHAKE)
 720         {
 721             md_print_warn(cr, fplog, "With twin-range cut-off's and SHAKE the virial and pressure are incorrect\n");
 722             if (ir->epc != epcNO)
 723             {
 724                 gmx_fatal(FARGS, "Can not do pressure coupling with twin-range cut-off's and SHAKE");
 725             }
 726         }
 727         check_nst_param(fplog, cr, "nstlist", ir->nstlist,
 728                         "nstcalcenergy", &ir->nstcalcenergy);
 729         if (ir->epc != epcNO)
 730         {
 731             check_nst_param(fplog, cr, "nstlist", ir->nstlist,
 732                             "nstpcouple", &ir->nstpcouple);
 733         }
 734         check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
 735                         "nstenergy", &ir->nstenergy);
 736         check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
 737                         "nstlog", &ir->nstlog);
 738         if (ir->efep != efepNO)
 739         {
 740             check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy,
 741                             "nstdhdl", &ir->fepvals->nstdhdl);
 742         }
 743     }
 744
 745     if (EI_VV(ir->eI) && IR_TWINRANGE(*ir) && ir->nstlist > 1)
 746     {
 747         gmx_fatal(FARGS, "Twin-range multiple time stepping does not work with integrator %s.", ei_names[ir->eI]);
 748     }
 749 }
 750
 751 void rerun_parallel_comm(t_commrec *cr, t_trxframe *fr,
 752                          gmx_bool *bNotLastFrame)
 753 {
 754     gmx_bool bAlloc;
 755     rvec    *xp, *vp;
 756
 757     bAlloc = (fr->natoms == 0);
 758
 759     if (MASTER(cr) && !*bNotLastFrame)
 760     {
 761         fr->natoms = -1;
 762     }
 763     xp = fr->x;
 764     vp = fr->v;
 765     gmx_bcast(sizeof(*fr), fr, cr);
 766     fr->x = xp;
 767     fr->v = vp;
 768
 769     *bNotLastFrame = (fr->natoms >= 0);
 770
 771 }