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