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,2017,2018,2019, 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 #include <cmath>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/domdec/domdec.h"
  48 #include "gromacs/gmxlib/network.h"
  49 #include "gromacs/gmxlib/nrnb.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdlib/dispersioncorrection.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/enerdata.h"
  60 #include "gromacs/mdtypes/energyhistory.h"
  61 #include "gromacs/mdtypes/forcerec.h"
  62 #include "gromacs/mdtypes/group.h"
  63 #include "gromacs/mdtypes/inputrec.h"
  64 #include "gromacs/mdtypes/md_enums.h"
  65 #include "gromacs/mdtypes/state.h"
  66 #include "gromacs/pbcutil/pbc.h"
  67 #include "gromacs/pulling/pull.h"
  68 #include "gromacs/timing/wallcycle.h"
  69 #include "gromacs/topology/mtop_util.h"
  70 #include "gromacs/trajectory/trajectoryframe.h"
  71 #include "gromacs/utility/arrayref.h"
  72 #include "gromacs/utility/cstringutil.h"
  73 #include "gromacs/utility/fatalerror.h"
  74 #include "gromacs/utility/gmxassert.h"
  75 #include "gromacs/utility/logger.h"
  76 #include "gromacs/utility/smalloc.h"
  77 #include "gromacs/utility/snprintf.h"
  78
  79 // TODO move this to multi-sim module
  80 bool multisim_int_all_are_equal(const gmx_multisim_t *ms,
  81                                 int64_t               value)
  82 {
  83     bool         allValuesAreEqual = true;
  84     int64_t     *buf;
  85
  86     GMX_RELEASE_ASSERT(ms, "Invalid use of multi-simulation pointer");
  87
  88     snew(buf, ms->nsim);
  89     /* send our value to all other master ranks, receive all of theirs */
  90     buf[ms->sim] = value;
  91     gmx_sumli_sim(ms->nsim, buf, ms);
  92
  93     for (int s = 0; s < ms->nsim; s++)
  94     {
  95         if (buf[s] != value)
  96         {
  97             allValuesAreEqual = false;
  98             break;
  99         }
 100     }
 101
 102     sfree(buf);
 103
 104     return allValuesAreEqual;
 105 }
 106
 107 int multisim_min(const gmx_multisim_t *ms, int nmin, int n)
 108 {
 109     int     *buf;
 110     gmx_bool bPos, bEqual;
 111     int      s, d;
 112
 113     snew(buf, ms->nsim);
 114     buf[ms->sim] = n;
 115     gmx_sumi_sim(ms->nsim, buf, ms);
 116     bPos   = TRUE;
 117     bEqual = TRUE;
 118     for (s = 0; s < ms->nsim; s++)
 119     {
 120         bPos   = bPos   && (buf[s] > 0);
 121         bEqual = bEqual && (buf[s] == buf[0]);
 122     }
 123     if (bPos)
 124     {
 125         if (bEqual)
 126         {
 127             nmin = std::min(nmin, buf[0]);
 128         }
 129         else
 130         {
 131             /* Find the least common multiple */
 132             for (d = 2; d < nmin; d++)
 133             {
 134                 s = 0;
 135                 while (s < ms->nsim && d % buf[s] == 0)
 136                 {
 137                     s++;
 138                 }
 139                 if (s == ms->nsim)
 140                 {
 141                     /* We found the LCM and it is less than nmin */
 142                     nmin = d;
 143                     break;
 144                 }
 145             }
 146         }
 147     }
 148     sfree(buf);
 149
 150     return nmin;
 151 }
 152
 153 /* TODO Specialize this routine into init-time and loop-time versions?
 154    e.g. bReadEkin is only true when restoring from checkpoint */
 155 void compute_globals(FILE *fplog, gmx_global_stat *gstat, t_commrec *cr, t_inputrec *ir,
 156                      t_forcerec *fr, gmx_ekindata_t *ekind,
 157                      t_state *state, t_mdatoms *mdatoms,
 158                      t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
 159                      gmx_enerdata_t *enerd, tensor force_vir, tensor shake_vir, tensor total_vir,
 160                      tensor pres, rvec mu_tot, gmx::Constraints *constr,
 161                      gmx::SimulationSignaller *signalCoordinator,
 162                      matrix box, int *totalNumberOfBondedInteractions,
 163                      gmx_bool *bSumEkinhOld, int flags)
 164 {
 165     tensor   corr_vir, corr_pres;
 166     gmx_bool bEner, bPres, bTemp;
 167     gmx_bool bStopCM, bGStat,
 168              bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
 169     gmx_bool bCheckNumberOfBondedInteractions;
 170     real     prescorr, enercorr, dvdlcorr, dvdl_ekin;
 171
 172     /* translate CGLO flags to gmx_booleans */
 173     bStopCM                          = ((flags & CGLO_STOPCM) != 0);
 174     bGStat                           = ((flags & CGLO_GSTAT) != 0);
 175     bReadEkin                        = ((flags & CGLO_READEKIN) != 0);
 176     bScaleEkin                       = ((flags & CGLO_SCALEEKIN) != 0);
 177     bEner                            = ((flags & CGLO_ENERGY) != 0);
 178     bTemp                            = ((flags & CGLO_TEMPERATURE) != 0);
 179     bPres                            = ((flags & CGLO_PRESSURE) != 0);
 180     bConstrain                       = ((flags & CGLO_CONSTRAINT) != 0);
 181     bCheckNumberOfBondedInteractions = ((flags & CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS) != 0);
 182
 183     /* we calculate a full state kinetic energy either with full-step velocity verlet
 184        or half step where we need the pressure */
 185
 186     bEkinAveVel = (ir->eI == eiVV || (ir->eI == eiVVAK && bPres) || bReadEkin);
 187
 188     /* in initalization, it sums the shake virial in vv, and to
 189        sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
 190
 191     /* ########## Kinetic energy  ############## */
 192
 193     if (bTemp)
 194     {
 195         /* Non-equilibrium MD: this is parallellized, but only does communication
 196          * when there really is NEMD.
 197          */
 198
 199         if (PAR(cr) && (ekind->bNEMD))
 200         {
 201             accumulate_u(cr, &(ir->opts), ekind);
 202         }
 203         if (!bReadEkin)
 204         {
 205             calc_ke_part(state, &(ir->opts), mdatoms, ekind, nrnb, bEkinAveVel);
 206         }
 207     }
 208
 209     /* Calculate center of mass velocity if necessary, also parallellized */
 210     if (bStopCM)
 211     {
 212         calc_vcm_grp(0, mdatoms->homenr, mdatoms,
 213                      state->x.rvec_array(), state->v.rvec_array(), vcm);
 214     }
 215
 216     if (bTemp || bStopCM || bPres || bEner || bConstrain || bCheckNumberOfBondedInteractions)
 217     {
 218         if (!bGStat)
 219         {
 220             /* We will not sum ekinh_old,
 221              * so signal that we still have to do it.
 222              */
 223             *bSumEkinhOld = TRUE;
 224
 225         }
 226         else
 227         {
 228             gmx::ArrayRef<real> signalBuffer = signalCoordinator->getCommunicationBuffer();
 229             if (PAR(cr))
 230             {
 231                 wallcycle_start(wcycle, ewcMoveE);
 232                 global_stat(gstat, cr, enerd, force_vir, shake_vir, mu_tot,
 233                             ir, ekind, constr, bStopCM ? vcm : nullptr,
 234                             signalBuffer.size(), signalBuffer.data(),
 235                             totalNumberOfBondedInteractions,
 236                             *bSumEkinhOld, flags);
 237                 wallcycle_stop(wcycle, ewcMoveE);
 238             }
 239             signalCoordinator->finalizeSignals();
 240             *bSumEkinhOld = FALSE;
 241         }
 242     }
 243
 244     /* Do center of mass motion removal */
 245     if (bStopCM)
 246     {
 247         check_cm_grp(fplog, vcm, ir, 1);
 248         /* At initialization, do not pass x with acceleration-correction mode
 249          * to avoid (incorrect) correction of the initial coordinates.
 250          */
 251         rvec *xPtr = nullptr;
 252         if (vcm->mode == ecmANGULAR || (vcm->mode == ecmLINEAR_ACCELERATION_CORRECTION && !(flags & CGLO_INITIALIZATION)))
 253         {
 254             xPtr = state->x.rvec_array();
 255         }
 256         do_stopcm_grp(*mdatoms,
 257                       xPtr, state->v.rvec_array(), *vcm);
 258         inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 259     }
 260
 261     if (bEner)
 262     {
 263         /* Calculate the amplitude of the cosine velocity profile */
 264         ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
 265     }
 266
 267     if (bTemp)
 268     {
 269         /* Sum the kinetic energies of the groups & calc temp */
 270         /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with inputrecNptTrotter */
 271         /* three maincase:  VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
 272            Leap with AveVel is not supported; it's not clear that it will actually work.
 273            bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
 274            If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
 275          */
 276         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin,
 277                                        bEkinAveVel, bScaleEkin);
 278         enerd->dvdl_lin[efptMASS] = static_cast<double>(dvdl_ekin);
 279
 280         enerd->term[F_EKIN] = trace(ekind->ekin);
 281     }
 282
 283     /* ##########  Long range energy information ###### */
 284
 285     if (bEner || bPres || bConstrain)
 286     {
 287         calc_dispcorr(ir, fr, box, state->lambda[efptVDW],
 288                       corr_pres, corr_vir, &prescorr, &enercorr, &dvdlcorr);
 289     }
 290
 291     if (bEner)
 292     {
 293         enerd->term[F_DISPCORR]  = enercorr;
 294         enerd->term[F_EPOT]     += enercorr;
 295         enerd->term[F_DVDL_VDW] += dvdlcorr;
 296     }
 297
 298     /* ########## Now pressure ############## */
 299     if (bPres || bConstrain)
 300     {
 301
 302         m_add(force_vir, shake_vir, total_vir);
 303
 304         /* Calculate pressure and apply LR correction if PPPM is used.
 305          * Use the box from last timestep since we already called update().
 306          */
 307
 308         enerd->term[F_PRES] = calc_pres(fr->ePBC, ir->nwall, box, ekind->ekin, total_vir, pres);
 309
 310         /* Calculate long range corrections to pressure and energy */
 311         /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
 312            and computes enerd->term[F_DISPCORR].  Also modifies the
 313            total_vir and pres tesors */
 314
 315         m_add(total_vir, corr_vir, total_vir);
 316         m_add(pres, corr_pres, pres);
 317         enerd->term[F_PDISPCORR] = prescorr;
 318         enerd->term[F_PRES]     += prescorr;
 319     }
 320 }
 321
 322 /* check whether an 'nst'-style parameter p is a multiple of nst, and
 323    set it to be one if not, with a warning. */
 324 static void check_nst_param(const gmx::MDLogger &mdlog,
 325                             const char *desc_nst, int nst,
 326                             const char *desc_p, int *p)
 327 {
 328     if (*p > 0 && *p % nst != 0)
 329     {
 330         /* Round up to the next multiple of nst */
 331         *p = ((*p)/nst + 1)*nst;
 332         GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
 333                 "NOTE: %s changes %s to %d", desc_nst, desc_p, *p);
 334     }
 335 }
 336
 337 void setCurrentLambdasRerun(int64_t step, const t_lambda *fepvals,
 338                             const t_trxframe *rerun_fr, const double *lam0,
 339                             t_state *globalState)
 340 {
 341     GMX_RELEASE_ASSERT(globalState != nullptr, "setCurrentLambdasGlobalRerun should be called with a valid state object");
 342
 343     if (rerun_fr->bLambda)
 344     {
 345         if (fepvals->delta_lambda == 0)
 346         {
 347             globalState->lambda[efptFEP] = rerun_fr->lambda;
 348         }
 349         else
 350         {
 351             /* find out between which two value of lambda we should be */
 352             real frac      = step*fepvals->delta_lambda;
 353             int  fep_state = static_cast<int>(std::floor(frac*fepvals->n_lambda));
 354             /* interpolate between this state and the next */
 355             /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 356             frac           = frac*fepvals->n_lambda - fep_state;
 357             for (int i = 0; i < efptNR; i++)
 358             {
 359                 globalState->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
 360                     frac*(fepvals->all_lambda[i][fep_state+1] - fepvals->all_lambda[i][fep_state]);
 361             }
 362         }
 363     }
 364     else if (rerun_fr->bFepState)
 365     {
 366         globalState->fep_state = rerun_fr->fep_state;
 367         for (int i = 0; i < efptNR; i++)
 368         {
 369             globalState->lambda[i] = fepvals->all_lambda[i][globalState->fep_state];
 370         }
 371     }
 372 }
 373
 374 void setCurrentLambdasLocal(int64_t step, const t_lambda *fepvals,
 375                             const double *lam0, t_state *state)
 376 /* find the current lambdas.  If rerunning, we either read in a state, or a lambda value,
 377    requiring different logic. */
 378 {
 379     if (fepvals->delta_lambda != 0)
 380     {
 381         /* find out between which two value of lambda we should be */
 382         real frac = step*fepvals->delta_lambda;
 383         if (fepvals->n_lambda > 0)
 384         {
 385             int fep_state = static_cast<int>(std::floor(frac*fepvals->n_lambda));
 386             /* interpolate between this state and the next */
 387             /* this assumes that the initial lambda corresponds to lambda==0, which is verified in grompp */
 388             frac          = frac*fepvals->n_lambda - fep_state;
 389             for (int i = 0; i < efptNR; i++)
 390             {
 391                 state->lambda[i] = lam0[i] + (fepvals->all_lambda[i][fep_state]) +
 392                     frac*(fepvals->all_lambda[i][fep_state + 1] - fepvals->all_lambda[i][fep_state]);
 393             }
 394         }
 395         else
 396         {
 397             for (int i = 0; i < efptNR; i++)
 398             {
 399                 state->lambda[i] = lam0[i] + frac;
 400             }
 401         }
 402     }
 403     else
 404     {
 405         /* if < 0, fep_state was never defined, and we should not set lambda from the state */
 406         if (state->fep_state > -1)
 407         {
 408             for (int i = 0; i < efptNR; i++)
 409             {
 410                 state->lambda[i] = fepvals->all_lambda[i][state->fep_state];
 411             }
 412         }
 413     }
 414 }
 415
 416 static void min_zero(int *n, int i)
 417 {
 418     if (i > 0 && (*n == 0 || i < *n))
 419     {
 420         *n = i;
 421     }
 422 }
 423
 424 static int lcd4(int i1, int i2, int i3, int i4)
 425 {
 426     int nst;
 427
 428     nst = 0;
 429     min_zero(&nst, i1);
 430     min_zero(&nst, i2);
 431     min_zero(&nst, i3);
 432     min_zero(&nst, i4);
 433     if (nst == 0)
 434     {
 435         gmx_incons("All 4 inputs for determining nstglobalcomm are <= 0");
 436     }
 437
 438     while (nst > 1 && ((i1 > 0 && i1 % nst != 0)  ||
 439                        (i2 > 0 && i2 % nst != 0)  ||
 440                        (i3 > 0 && i3 % nst != 0)  ||
 441                        (i4 > 0 && i4 % nst != 0)))
 442     {
 443         nst--;
 444     }
 445
 446     return nst;
 447 }
 448
 449 int check_nstglobalcomm(const gmx::MDLogger &mdlog, int nstglobalcomm, t_inputrec *ir, const t_commrec * cr)
 450 {
 451     if (!EI_DYNAMICS(ir->eI))
 452     {
 453         nstglobalcomm = 1;
 454     }
 455
 456     if (nstglobalcomm == -1)
 457     {
 458         // Set up the default behaviour
 459         if (!(ir->nstcalcenergy > 0 ||
 460               ir->nstlist > 0 ||
 461               ir->etc != etcNO ||
 462               ir->epc != epcNO))
 463         {
 464             /* The user didn't choose the period for anything
 465                important, so we just make sure we can send signals and
 466                write output suitably. */
 467             nstglobalcomm = 10;
 468             if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
 469             {
 470                 nstglobalcomm = ir->nstenergy;
 471             }
 472         }
 473         else
 474         {
 475             /* The user has made a choice (perhaps implicitly), so we
 476              * ensure that we do timely intra-simulation communication
 477              * for (possibly) each of the four parts that care.
 478              *
 479              * TODO Does the Verlet scheme (+ DD) need any
 480              * communication at nstlist steps? Is the use of nstlist
 481              * here a leftover of the twin-range scheme? Can we remove
 482              * nstlist when we remove the group scheme?
 483              */
 484             nstglobalcomm = lcd4(ir->nstcalcenergy,
 485                                  ir->nstlist,
 486                                  ir->etc != etcNO ? ir->nsttcouple : 0,
 487                                  ir->epc != epcNO ? ir->nstpcouple : 0);
 488         }
 489     }
 490     else
 491     {
 492         // Check that the user's choice of mdrun -gcom will work
 493         if (ir->nstlist > 0 &&
 494             nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
 495         {
 496             nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
 497             GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
 498                     "WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d",
 499                     nstglobalcomm);
 500         }
 501         if (ir->nstcalcenergy > 0)
 502         {
 503             check_nst_param(mdlog, "-gcom", nstglobalcomm,
 504                             "nstcalcenergy", &ir->nstcalcenergy);
 505         }
 506         if (ir->etc != etcNO && ir->nsttcouple > 0)
 507         {
 508             check_nst_param(mdlog, "-gcom", nstglobalcomm,
 509                             "nsttcouple", &ir->nsttcouple);
 510         }
 511         if (ir->epc != epcNO && ir->nstpcouple > 0)
 512         {
 513             check_nst_param(mdlog, "-gcom", nstglobalcomm,
 514                             "nstpcouple", &ir->nstpcouple);
 515         }
 516
 517         check_nst_param(mdlog, "-gcom", nstglobalcomm,
 518                         "nstenergy", &ir->nstenergy);
 519
 520         check_nst_param(mdlog, "-gcom", nstglobalcomm,
 521                         "nstlog", &ir->nstlog);
 522     }
 523
 524     if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
 525     {
 526         GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
 527                 "WARNING: Changing nstcomm from %d to %d",
 528                 ir->nstcomm, nstglobalcomm);
 529         ir->nstcomm = nstglobalcomm;
 530     }
 531
 532     if (cr->nnodes > 1)
 533     {
 534         GMX_LOG(mdlog.info).appendTextFormatted(
 535                 "Intra-simulation communication will occur every %d steps.\n", nstglobalcomm);
 536     }
 537     return nstglobalcomm;
 538
 539 }
 540
 541 void rerun_parallel_comm(t_commrec *cr, t_trxframe *fr,
 542                          gmx_bool *bLastStep)
 543 {
 544     rvec    *xp, *vp;
 545
 546     if (MASTER(cr) && *bLastStep)
 547     {
 548         fr->natoms = -1;
 549     }
 550     xp = fr->x;
 551     vp = fr->v;
 552     gmx_bcast(sizeof(*fr), fr, cr);
 553     fr->x = xp;
 554     fr->v = vp;
 555
 556     *bLastStep = (fr->natoms < 0);
 557
 558 }
 559
 560 // TODO Most of this logic seems to belong in the respective modules
 561 void set_state_entries(t_state *state, const t_inputrec *ir)
 562 {
 563     /* The entries in the state in the tpx file might not correspond
 564      * with what is needed, so we correct this here.
 565      */
 566     state->flags = 0;
 567     if (ir->efep != efepNO || ir->bExpanded)
 568     {
 569         state->flags |= (1<<estLAMBDA);
 570         state->flags |= (1<<estFEPSTATE);
 571     }
 572     state->flags |= (1<<estX);
 573     GMX_RELEASE_ASSERT(state->x.size() == state->natoms, "We should start a run with an initialized state->x");
 574     if (EI_DYNAMICS(ir->eI))
 575     {
 576         state->flags |= (1<<estV);
 577     }
 578
 579     state->nnhpres = 0;
 580     if (ir->ePBC != epbcNONE)
 581     {
 582         state->flags |= (1<<estBOX);
 583         if (inputrecPreserveShape(ir))
 584         {
 585             state->flags |= (1<<estBOX_REL);
 586         }
 587         if ((ir->epc == epcPARRINELLORAHMAN) || (ir->epc == epcMTTK))
 588         {
 589             state->flags |= (1<<estBOXV);
 590             state->flags |= (1<<estPRES_PREV);
 591         }
 592         if (inputrecNptTrotter(ir) || (inputrecNphTrotter(ir)))
 593         {
 594             state->nnhpres = 1;
 595             state->flags  |= (1<<estNHPRES_XI);
 596             state->flags  |= (1<<estNHPRES_VXI);
 597             state->flags  |= (1<<estSVIR_PREV);
 598             state->flags  |= (1<<estFVIR_PREV);
 599             state->flags  |= (1<<estVETA);
 600             state->flags  |= (1<<estVOL0);
 601         }
 602         if (ir->epc == epcBERENDSEN)
 603         {
 604             state->flags  |= (1<<estBAROS_INT);
 605         }
 606     }
 607
 608     if (ir->etc == etcNOSEHOOVER)
 609     {
 610         state->flags |= (1<<estNH_XI);
 611         state->flags |= (1<<estNH_VXI);
 612     }
 613
 614     if (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)
 615     {
 616         state->flags |= (1<<estTHERM_INT);
 617     }
 618
 619     init_gtc_state(state, state->ngtc, state->nnhpres, ir->opts.nhchainlength); /* allocate the space for nose-hoover chains */
 620     init_ekinstate(&state->ekinstate, ir);
 621
 622     if (ir->bExpanded)
 623     {
 624         snew(state->dfhist, 1);
 625         init_df_history(state->dfhist, ir->fepvals->n_lambda);
 626     }
 627
 628     if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
 629     {
 630         state->flags |= (1<<estPULLCOMPREVSTEP);
 631     }
 632 }