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