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,2015,2016, 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 "md_support.h"
  41
  42 #include <climits>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/domdec/domdec.h"
  47 #include "gromacs/gmxlib/md_logging.h"
  48 #include "gromacs/gmxlib/network.h"
  49 #include "gromacs/gmxlib/nrnb.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdlib/mdrun.h"
  52 #include "gromacs/mdlib/sim_util.h"
  53 #include "gromacs/mdlib/simulationsignal.h"
  54 #include "gromacs/mdlib/tgroup.h"
  55 #include "gromacs/mdlib/update.h"
  56 #include "gromacs/mdlib/vcm.h"
  57 #include "gromacs/mdtypes/commrec.h"
  58 #include "gromacs/mdtypes/df_history.h"
  59 #include "gromacs/mdtypes/energyhistory.h"
  60 #include "gromacs/mdtypes/group.h"
  61 #include "gromacs/mdtypes/inputrec.h"
  62 #include "gromacs/mdtypes/md_enums.h"
  63 #include "gromacs/mdtypes/state.h"
  64 #include "gromacs/pbcutil/pbc.h"
  65 #include "gromacs/timing/wallcycle.h"
  66 #include "gromacs/topology/mtop_util.h"
  67 #include "gromacs/trajectory/trajectoryframe.h"
  68 #include "gromacs/utility/arrayref.h"
  69 #include "gromacs/utility/cstringutil.h"
  70 #include "gromacs/utility/fatalerror.h"
  71 #include "gromacs/utility/gmxassert.h"
  72 #include "gromacs/utility/smalloc.h"
  73 #include "gromacs/utility/snprintf.h"
  74
  75 // TODO move this to multi-sim module
  76 bool multisim_int_all_are_equal(const gmx_multisim_t *ms,
  77                                 gmx_int64_t           value)
  78 {
  79     bool         allValuesAreEqual = true;
  80     gmx_int64_t *buf;
  81
  82     GMX_RELEASE_ASSERT(ms, "Invalid use of multi-simulation pointer");
  83
  84     snew(buf, ms->nsim);
  85     /* send our value to all other master ranks, receive all of theirs */
  86     buf[ms->sim] = value;
  87     gmx_sumli_sim(ms->nsim, buf, ms);
  88
  89     for (int s = 0; s < ms->nsim; s++)
  90     {
  91         if (buf[s] != value)
  92         {
  93             allValuesAreEqual = false;
  94             break;
  95         }
  96     }
  97
  98     sfree(buf);
  99
 100     return allValuesAreEqual;
 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 void copy_coupling_state(t_state *statea, t_state *stateb,
 150                          gmx_ekindata_t *ekinda, gmx_ekindata_t *ekindb, t_grpopts* opts)
 151 {
 152
 153     /* MRS note -- might be able to get rid of some of the arguments.  Look over it when it's all debugged */
 154
 155     int i, j, nc;
 156
 157     /* Make sure we have enough space for x and v */
 158     if (statea->nalloc > stateb->nalloc)
 159     {
 160         stateb->nalloc = statea->nalloc;
 161         /* We need to allocate one element extra, since we might use
 162          * (unaligned) 4-wide SIMD loads to access rvec entries.
 163          */
 164         srenew(stateb->x, stateb->nalloc + 1);
 165         srenew(stateb->v, stateb->nalloc + 1);
 166     }
 167
 168     stateb->natoms     = statea->natoms;
 169     stateb->ngtc       = statea->ngtc;
 170     stateb->nnhpres    = statea->nnhpres;
 171     stateb->veta       = statea->veta;
 172     if (ekinda)
 173     {
 174         copy_mat(ekinda->ekin, ekindb->ekin);
 175         for (i = 0; i < stateb->ngtc; i++)
 176         {
 177             ekindb->tcstat[i].T  = ekinda->tcstat[i].T;
 178             ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
 179             copy_mat(ekinda->tcstat[i].ekinh, ekindb->tcstat[i].ekinh);
 180             copy_mat(ekinda->tcstat[i].ekinf, ekindb->tcstat[i].ekinf);
 181             ekindb->tcstat[i].ekinscalef_nhc =  ekinda->tcstat[i].ekinscalef_nhc;
 182             ekindb->tcstat[i].ekinscaleh_nhc =  ekinda->tcstat[i].ekinscaleh_nhc;
 183             ekindb->tcstat[i].vscale_nhc     =  ekinda->tcstat[i].vscale_nhc;
 184         }
 185     }
 186     copy_rvecn(statea->x, stateb->x, 0, stateb->natoms);
 187     copy_rvecn(statea->v, stateb->v, 0, stateb->natoms);
 188     copy_mat(statea->box, stateb->box);
 189     copy_mat(statea->box_rel, stateb->box_rel);
 190     copy_mat(statea->boxv, stateb->boxv);
 191
 192     for (i = 0; i < stateb->ngtc; i++)
 193     {
 194         nc = i*opts->nhchainlength;
 195         for (j = 0; j < opts->nhchainlength; j++)
 196         {
 197             stateb->nosehoover_xi[nc+j]  = statea->nosehoover_xi[nc+j];
 198             stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
 199         }
 200     }
 201     if (stateb->nhpres_xi != NULL)
 202     {
 203         for (i = 0; i < stateb->nnhpres; i++)
 204         {
 205             nc = i*opts->nhchainlength;
 206             for (j = 0; j < opts->nhchainlength; j++)
 207             {
 208                 stateb->nhpres_xi[nc+j]  = statea->nhpres_xi[nc+j];
 209                 stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
 210             }
 211         }
 212     }
 213 }
 214
 215 real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
 216 {
 217     real quantity = 0;
 218     switch (ir->etc)
 219     {
 220         case etcNO:
 221             break;
 222         case etcBERENDSEN:
 223             break;
 224         case etcNOSEHOOVER:
 225             quantity = NPT_energy(ir, state, MassQ);
 226             break;
 227         case etcVRESCALE:
 228             quantity = vrescale_energy(&(ir->opts), state->therm_integral);
 229             break;
 230         default:
 231             break;
 232     }
 233     return quantity;
 234 }
 235
 236 /* TODO Specialize this routine into init-time and loop-time versions?
 237    e.g. bReadEkin is only true when restoring from checkpoint */
 238 void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_inputrec *ir,
 239                      t_forcerec *fr, gmx_ekindata_t *ekind,
 240                      t_state *state, t_mdatoms *mdatoms,
 241                      t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
 242                      gmx_enerdata_t *enerd, tensor force_vir, tensor shake_vir, tensor total_vir,
 243                      tensor pres, rvec mu_tot, gmx_constr_t constr,
 244                      gmx::SimulationSignaller *signalCoordinator,
 245                      matrix box, int *totalNumberOfBondedInteractions,
 246                      gmx_bool *bSumEkinhOld, int flags)
 247 {
 248     tensor   corr_vir, corr_pres;
 249     gmx_bool bEner, bPres, bTemp;
 250     gmx_bool bStopCM, bGStat,
 251              bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
 252     real     prescorr, enercorr, dvdlcorr, dvdl_ekin;
 253
 254     /* translate CGLO flags to gmx_booleans */
 255     bStopCM       = flags & CGLO_STOPCM;
 256     bGStat        = flags & CGLO_GSTAT;
 257     bReadEkin     = (flags & CGLO_READEKIN);
 258     bScaleEkin    = (flags & CGLO_SCALEEKIN);
 259     bEner         = flags & CGLO_ENERGY;
 260     bTemp         = flags & CGLO_TEMPERATURE;
 261     bPres         = (flags & CGLO_PRESSURE);
 262     bConstrain    = (flags & CGLO_CONSTRAINT);
 263
 264     /* we calculate a full state kinetic energy either with full-step velocity verlet
 265        or half step where we need the pressure */
 266
 267     bEkinAveVel = (ir->eI == eiVV || (ir->eI == eiVVAK && bPres) || bReadEkin);
 268
 269     /* in initalization, it sums the shake virial in vv, and to
 270        sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
 271
 272     /* ########## Kinetic energy  ############## */
 273
 274     if (bTemp)
 275     {
 276         /* Non-equilibrium MD: this is parallellized, but only does communication
 277          * when there really is NEMD.
 278          */
 279
 280         if (PAR(cr) && (ekind->bNEMD))
 281         {
 282             accumulate_u(cr, &(ir->opts), ekind);
 283         }
 284         if (!bReadEkin)
 285         {
 286             calc_ke_part(state, &(ir->opts), mdatoms, ekind, nrnb, bEkinAveVel);
 287         }
 288     }
 289
 290     /* Calculate center of mass velocity if necessary, also parallellized */
 291     if (bStopCM)
 292     {
 293         calc_vcm_grp(0, mdatoms->homenr, mdatoms,
 294                      state->x, state->v, vcm);
 295     }
 296
 297     if (bTemp || bStopCM || bPres || bEner || bConstrain)
 298     {
 299         if (!bGStat)
 300         {
 301             /* We will not sum ekinh_old,
 302              * so signal that we still have to do it.
 303              */
 304             *bSumEkinhOld = TRUE;
 305
 306         }
 307         else
 308         {
 309             gmx::ArrayRef<real> signalBuffer = signalCoordinator->getCommunicationBuffer();
 310             if (PAR(cr))
 311             {
 312                 wallcycle_start(wcycle, ewcMoveE);
 313                 global_stat(gstat, cr, enerd, force_vir, shake_vir, mu_tot,
 314                             ir, ekind, constr, bStopCM ? vcm : NULL,
 315                             signalBuffer.size(), signalBuffer.data(),
 316                             totalNumberOfBondedInteractions,
 317                             *bSumEkinhOld, flags);
 318                 wallcycle_stop(wcycle, ewcMoveE);
 319             }
 320             signalCoordinator->finalizeSignals();
 321             *bSumEkinhOld = FALSE;
 322         }
 323     }
 324
 325     if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
 326     {
 327         correct_ekin(debug,
 328                      0, mdatoms->homenr,
 329                      state->v, vcm->group_p[0],
 330                      mdatoms->massT, mdatoms->tmass, ekind->ekin);
 331     }
 332
 333     /* Do center of mass motion removal */
 334     if (bStopCM)
 335     {
 336         check_cm_grp(fplog, vcm, ir, 1);
 337         do_stopcm_grp(0, mdatoms->homenr, mdatoms->cVCM,
 338                       state->x, state->v, vcm);
 339         inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 340     }
 341
 342     if (bEner)
 343     {
 344         /* Calculate the amplitude of the cosine velocity profile */
 345         ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
 346     }
 347
 348     if (bTemp)
 349     {
 350         /* Sum the kinetic energies of the groups & calc temp */
 351         /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with inputrecNptTrotter */
 352         /* three maincase:  VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
 353            Leap with AveVel is not supported; it's not clear that it will actually work.
 354            bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
 355            If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
 356          */
 357         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin,
 358                                        bEkinAveVel, bScaleEkin);
 359         enerd->dvdl_lin[efptMASS] = (double) dvdl_ekin;
 360
 361         enerd->term[F_EKIN] = trace(ekind->ekin);
 362     }
 363
 364     /* ##########  Long range energy information ###### */
 365
 366     if (bEner || bPres || bConstrain)
 367     {
 368         calc_dispcorr(ir, fr, box, state->lambda[efptVDW],
 369                       corr_pres, corr_vir, &prescorr, &enercorr, &dvdlcorr);
 370     }
 371
 372     if (bEner)
 373     {
 374         enerd->term[F_DISPCORR]  = enercorr;
 375         enerd->term[F_EPOT]     += enercorr;
 376         enerd->term[F_DVDL_VDW] += dvdlcorr;
 377     }
 378
 379     /* ########## Now pressure ############## */
 380     if (bPres || bConstrain)
 381     {
 382
 383         m_add(force_vir, shake_vir, total_vir);
 384
 385         /* Calculate pressure and apply LR correction if PPPM is used.
 386          * Use the box from last timestep since we already called update().
 387          */
 388
 389         enerd->term[F_PRES] = calc_pres(fr->ePBC, ir->nwall, box, ekind->ekin, total_vir, pres);
 390
 391         /* Calculate long range corrections to pressure and energy */
 392         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
 393            and computes enerd->term[F_DISPCORR].  Also modifies the
 394            total_vir and pres tesors */
 395
 396         m_add(total_vir, corr_vir, total_vir);
 397         m_add(pres, corr_pres, pres);
 398         enerd->term[F_PDISPCORR] = prescorr;
 399         enerd->term[F_PRES]     += prescorr;
 400     }
 401 }
 402
 403 void check_nst_param(FILE *fplog, t_commrec *cr,
 404                      const char *desc_nst, int nst,
 405                      const char *desc_p, int *p)
 406 {
 407     if (*p > 0 && *p % nst != 0)
 408     {
 409         /* Round up to the next multiple of nst */
 410         *p = ((*p)/nst + 1)*nst;
 411         md_print_warn(cr, fplog,
 412                       "NOTE: %s changes %s to %d\n", desc_nst, desc_p, *p);
 413     }
 414 }
 415
 416 void set_current_lambdas(gmx_int64_t step, t_lambda *fepvals, gmx_bool bRerunMD,
 417                          t_trxframe *rerun_fr, t_state *state_global, t_state *state, double lam0[])
 418 /* find the current lambdas.  If rerunning, we either read in a state, or a lambda value,
 419    requiring different logic. */
 420 {
 421     real frac;
 422     int  i, fep_state = 0;
 423     if (bRerunMD)
 424     {
 425         if (rerun_fr->bLambda)
 426         {
 427             if (fepvals->delta_lambda == 0)
 428             {
 429                 state_global->lambda[efptFEP] = rerun_fr->lambda;
 430                 for (i = 0; i < efptNR; i++)
 431                 {
 432                     if (i != efptFEP)
 433                     {
 434                         state->lambda[i] = state_global->lambda[i];
 435                     }
 436                 }
 437             }
 438             else
 439             {
 440                 /* find out between which two value of lambda we should be */
 441                 frac      = (step*fepvals->delta_lambda);
 442                 fep_state = static_cast<int>(floor(frac*fepvals->n_lambda));
 443                 /* interpolate between this state and the next */
 444                 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 445                 frac = (frac*fepvals->n_lambda)-fep_state;
 446                 for (i = 0; i < efptNR; i++)
 447                 {
 448                     state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
 449                         frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
 450                 }
 451             }
 452         }
 453         else if (rerun_fr->bFepState)
 454         {
 455             state_global->fep_state = rerun_fr->fep_state;
 456             for (i = 0; i < efptNR; i++)
 457             {
 458                 state_global->lambda[i] = fepvals->all_lambda[i][fep_state];
 459             }
 460         }
 461     }
 462     else
 463     {
 464         if (fepvals->delta_lambda != 0)
 465         {
 466             /* find out between which two value of lambda we should be */
 467             frac = (step*fepvals->delta_lambda);
 468             if (fepvals->n_lambda > 0)
 469             {
 470                 fep_state = static_cast<int>(floor(frac*fepvals->n_lambda));
 471                 /* interpolate between this state and the next */
 472                 /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 473                 frac = (frac*fepvals->n_lambda)-fep_state;
 474                 for (i = 0; i < efptNR; i++)
 475                 {
 476                     state_global->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
 477                         frac*(fepvals->all_lambda[i][fep_state+1]-fepvals->all_lambda[i][fep_state]);
 478                 }
 479             }
 480             else
 481             {
 482                 for (i = 0; i < efptNR; i++)
 483                 {
 484                     state_global->lambda[i] = lam0[i] + frac;
 485                 }
 486             }
 487         }
 488         else
 489         {
 490             /* if < 0, fep_state was never defined, and we should not set lambda from the state */
 491             if (state_global->fep_state > -1)
 492             {
 493                 state_global->fep_state = state->fep_state; /* state->fep_state is the one updated by bExpanded */
 494                 for (i = 0; i < efptNR; i++)
 495                 {
 496                     state_global->lambda[i] = fepvals->all_lambda[i][state_global->fep_state];
 497                 }
 498             }
 499         }
 500     }
 501     for (i = 0; i < efptNR; i++)
 502     {
 503         state->lambda[i] = state_global->lambda[i];
 504     }
 505 }
 506
 507 static void min_zero(int *n, int i)
 508 {
 509     if (i > 0 && (*n == 0 || i < *n))
 510     {
 511         *n = i;
 512     }
 513 }
 514
 515 static int lcd4(int i1, int i2, int i3, int i4)
 516 {
 517     int nst;
 518
 519     nst = 0;
 520     min_zero(&nst, i1);
 521     min_zero(&nst, i2);
 522     min_zero(&nst, i3);
 523     min_zero(&nst, i4);
 524     if (nst == 0)
 525     {
 526         gmx_incons("All 4 inputs for determining nstglobalcomm are <= 0");
 527     }
 528
 529     while (nst > 1 && ((i1 > 0 && i1 % nst != 0)  ||
 530                        (i2 > 0 && i2 % nst != 0)  ||
 531                        (i3 > 0 && i3 % nst != 0)  ||
 532                        (i4 > 0 && i4 % nst != 0)))
 533     {
 534         nst--;
 535     }
 536
 537     return nst;
 538 }
 539
 540 int check_nstglobalcomm(FILE *fplog, t_commrec *cr,
 541                         int nstglobalcomm, t_inputrec *ir)
 542 {
 543     if (!EI_DYNAMICS(ir->eI))
 544     {
 545         nstglobalcomm = 1;
 546     }
 547
 548     if (nstglobalcomm == -1)
 549     {
 550         // Set up the default behaviour
 551         if (!(ir->nstcalcenergy > 0 ||
 552               ir->nstlist > 0 ||
 553               ir->etc != etcNO ||
 554               ir->epc != epcNO))
 555         {
 556             /* The user didn't choose the period for anything
 557                important, so we just make sure we can send signals and
 558                write output suitably. */
 559             nstglobalcomm = 10;
 560             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
 561             {
 562                 nstglobalcomm = ir->nstenergy;
 563             }
 564         }
 565         else
 566         {
 567             /* The user has made a choice (perhaps implicitly), so we
 568              * ensure that we do timely intra-simulation communication
 569              * for (possibly) each of the four parts that care.
 570              *
 571              * TODO Does the Verlet scheme (+ DD) need any
 572              * communication at nstlist steps? Is the use of nstlist
 573              * here a leftover of the twin-range scheme? Can we remove
 574              * nstlist when we remove the group scheme?
 575              */
 576             nstglobalcomm = lcd4(ir->nstcalcenergy,
 577                                  ir->nstlist,
 578                                  ir->etc != etcNO ? ir->nsttcouple : 0,
 579                                  ir->epc != epcNO ? ir->nstpcouple : 0);
 580         }
 581     }
 582     else
 583     {
 584         // Check that the user's choice of mdrun -gcom will work
 585         if (ir->nstlist > 0 &&
 586             nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
 587         {
 588             nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
 589             md_print_warn(cr, fplog, "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n", nstglobalcomm);
 590         }
 591         if (ir->nstcalcenergy > 0)
 592         {
 593             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 594                             "nstcalcenergy", &ir->nstcalcenergy);
 595         }
 596         if (ir->etc != etcNO && ir->nsttcouple > 0)
 597         {
 598             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 599                             "nsttcouple", &ir->nsttcouple);
 600         }
 601         if (ir->epc != epcNO && ir->nstpcouple > 0)
 602         {
 603             check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 604                             "nstpcouple", &ir->nstpcouple);
 605         }
 606
 607         check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 608                         "nstenergy", &ir->nstenergy);
 609
 610         check_nst_param(fplog, cr, "-gcom", nstglobalcomm,
 611                         "nstlog", &ir->nstlog);
 612     }
 613
 614     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
 615     {
 616         md_print_warn(cr, fplog, "WARNING: Changing nstcomm from %d to %d\n",
 617                       ir->nstcomm, nstglobalcomm);
 618         ir->nstcomm = nstglobalcomm;
 619     }
 620
 621     if (fplog)
 622     {
 623         fprintf(fplog, "Intra-simulation communication will occur every %d steps.\n", nstglobalcomm);
 624     }
 625     return nstglobalcomm;
 626 }
 627
 628 void rerun_parallel_comm(t_commrec *cr, t_trxframe *fr,
 629                          gmx_bool *bNotLastFrame)
 630 {
 631     rvec    *xp, *vp;
 632
 633     if (MASTER(cr) && !*bNotLastFrame)
 634     {
 635         fr->natoms = -1;
 636     }
 637     xp = fr->x;
 638     vp = fr->v;
 639     gmx_bcast(sizeof(*fr), fr, cr);
 640     fr->x = xp;
 641     fr->v = vp;
 642
 643     *bNotLastFrame = (fr->natoms >= 0);
 644
 645 }
 646
 647 void set_state_entries(t_state *state, const t_inputrec *ir)
 648 {
 649     /* The entries in the state in the tpx file might not correspond
 650      * with what is needed, so we correct this here.
 651      */
 652     state->flags = 0;
 653     if (ir->efep != efepNO || ir->bExpanded)
 654     {
 655         state->flags |= (1<<estLAMBDA);
 656         state->flags |= (1<<estFEPSTATE);
 657     }
 658     state->flags |= (1<<estX);
 659     if (state->lambda == NULL)
 660     {
 661         snew(state->lambda, efptNR);
 662     }
 663     if (state->x == NULL)
 664     {
 665         /* We need to allocate one element extra, since we might use
 666          * (unaligned) 4-wide SIMD loads to access rvec entries.
 667          */
 668         snew(state->x, state->nalloc + 1);
 669     }
 670     if (EI_DYNAMICS(ir->eI))
 671     {
 672         state->flags |= (1<<estV);
 673         if (state->v == NULL)
 674         {
 675             snew(state->v, state->nalloc + 1);
 676         }
 677     }
 678     if (ir->eI == eiCG)
 679     {
 680         state->flags |= (1<<estCGP);
 681         if (state->cg_p == NULL)
 682         {
 683             /* cg_p is not stored in the tpx file, so we need to allocate it */
 684             snew(state->cg_p, state->nalloc + 1);
 685         }
 686     }
 687
 688     state->nnhpres = 0;
 689     if (ir->ePBC != epbcNONE)
 690     {
 691         state->flags |= (1<<estBOX);
 692         if (inputrecPreserveShape(ir))
 693         {
 694             state->flags |= (1<<estBOX_REL);
 695         }
 696         if ((ir->epc == epcPARRINELLORAHMAN) || (ir->epc == epcMTTK))
 697         {
 698             state->flags |= (1<<estBOXV);
 699         }
 700         if (ir->epc != epcNO)
 701         {
 702             if (inputrecNptTrotter(ir) || (inputrecNphTrotter(ir)))
 703             {
 704                 state->nnhpres = 1;
 705                 state->flags  |= (1<<estNHPRES_XI);
 706                 state->flags  |= (1<<estNHPRES_VXI);
 707                 state->flags  |= (1<<estSVIR_PREV);
 708                 state->flags  |= (1<<estFVIR_PREV);
 709                 state->flags  |= (1<<estVETA);
 710                 state->flags  |= (1<<estVOL0);
 711             }
 712             else
 713             {
 714                 state->flags |= (1<<estPRES_PREV);
 715             }
 716         }
 717     }
 718
 719     if (ir->etc == etcNOSEHOOVER)
 720     {
 721         state->flags |= (1<<estNH_XI);
 722         state->flags |= (1<<estNH_VXI);
 723     }
 724
 725     if (ir->etc == etcVRESCALE)
 726     {
 727         state->flags |= (1<<estTC_INT);
 728     }
 729
 730     init_gtc_state(state, state->ngtc, state->nnhpres, ir->opts.nhchainlength); /* allocate the space for nose-hoover chains */
 731     init_ekinstate(&state->ekinstate, ir);
 732     snew(state->enerhist, 1);
 733     init_energyhistory(state->enerhist);
 734     init_df_history(&state->dfhist, ir->fepvals->n_lambda);
 735     state->swapstate.eSwapCoords = ir->eSwapCoords;
 736 }