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