src/gromacs/mdlib/md_support.cpp

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