src/gromacs/mdlib/mdebin_bar.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.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "gromacs/utility/arrayref.h"
  41 #include "mdebin_bar.h"
  42
  43 #include <cassert>
  44 #include <cfloat>
  45 #include <cmath>
  46 #include <cstring>
  47
  48 #include <memory>
  49
  50 #include "gromacs/fileio/enxio.h"
  51 #include "gromacs/mdlib/energyoutput.h"
  52 #include "gromacs/mdtypes/energyhistory.h"
  53 #include "gromacs/mdtypes/inputrec.h"
  54 #include "gromacs/mdtypes/md_enums.h"
  55 #include "gromacs/trajectory/energyframe.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/gmxassert.h"
  58 #include "gromacs/utility/smalloc.h"
  59
  60 /* Number of entries in subblock_d preceding the lambda components */
  61 constexpr int c_subblockDNumPreEntries = 5;
  62 /* Number of entries in subblock_i preceding the lambda components */
  63 constexpr int c_subblockINumPreEntries = 2;
  64
  65 /* reset the delta_h list to prepare it for new values */
  66 static void mde_delta_h_reset(t_mde_delta_h* dh)
  67 {
  68     dh->ndh     = 0;
  69     dh->written = FALSE;
  70 }
  71
  72 /* initialize the delta_h list */
  73 static void mde_delta_h_init(t_mde_delta_h* dh,
  74                              int            nbins,
  75                              double         dx,
  76                              unsigned int   ndhmax,
  77                              int            type,
  78                              int            derivative,
  79                              int            nlambda,
  80                              const double*  lambda)
  81 {
  82     int i;
  83
  84     dh->type       = type;
  85     dh->derivative = derivative;
  86     dh->nlambda    = nlambda;
  87
  88     dh->lambda.resize(nlambda);
  89     for (i = 0; i < nlambda; i++)
  90     {
  91         assert(lambda);
  92         dh->lambda[i] = lambda[i];
  93     }
  94
  95
  96     dh->subblock_meta_d.resize(dh->nlambda + 1);
  97
  98     dh->ndhmax = ndhmax + 2;
  99
 100     dh->dh.resize(dh->ndhmax);
 101     dh->dhf.resize(dh->ndhmax);
 102
 103     if (nbins <= 0 || dx < GMX_REAL_EPS * 10)
 104     {
 105         dh->nhist = 0;
 106     }
 107     else
 108     {
 109         int i;
 110         /* pre-allocate the histogram */
 111         dh->nhist = 2; /* energies and derivatives histogram */
 112         dh->dx    = dx;
 113         dh->nbins = nbins;
 114         for (i = 0; i < dh->nhist; i++)
 115         {
 116             dh->bin[i].resize(dh->nbins);
 117         }
 118     }
 119     mde_delta_h_reset(dh);
 120 }
 121
 122 /* Add a value to the delta_h list */
 123 static void mde_delta_h_add_dh(t_mde_delta_h* dh, double delta_h)
 124 {
 125     if (dh->ndh >= dh->ndhmax)
 126     {
 127         gmx_incons("delta_h array not big enough!");
 128     }
 129     dh->dh[dh->ndh] = delta_h;
 130     dh->ndh++;
 131 }
 132
 133 /* construct histogram with index hi */
 134 static void mde_delta_h_make_hist(t_mde_delta_h* dh, int hi, gmx_bool invert)
 135 {
 136     double       min_dh = FLT_MAX;
 137     double       max_dh = -FLT_MAX;
 138     unsigned int i;
 139     double       max_dh_hist; /* maximum binnable dh value */
 140     double       min_dh_hist; /* minimum binnable dh value */
 141     double       dx = dh->dx;
 142     double       f; /* energy mult. factor */
 143
 144     /* by applying a -1 scaling factor on the energies we get the same as
 145        having a negative dx, but we don't need to fix the min/max values
 146        beyond inverting x0 */
 147     f = invert ? -1 : 1;
 148
 149     /* first find min and max */
 150     for (i = 0; i < dh->ndh; i++)
 151     {
 152         if (f * dh->dh[i] < min_dh)
 153         {
 154             min_dh = f * dh->dh[i];
 155         }
 156         if (f * dh->dh[i] > max_dh)
 157         {
 158             max_dh = f * dh->dh[i];
 159         }
 160     }
 161
 162     /* reset the histogram */
 163     for (i = 0; i < dh->nbins; i++)
 164     {
 165         dh->bin[hi][i] = 0;
 166     }
 167     dh->maxbin[hi] = 0;
 168
 169     /* The starting point of the histogram is the lowest value found:
 170        that value has the highest contribution to the free energy.
 171
 172        Get this start value in number of histogram dxs from zero,
 173        as an integer.*/
 174
 175     dh->x0[hi] = static_cast<int64_t>(floor(min_dh / dx));
 176
 177     min_dh_hist = (dh->x0[hi]) * dx;
 178     max_dh_hist = (dh->x0[hi] + dh->nbins + 1) * dx;
 179
 180     /* and fill the histogram*/
 181     for (i = 0; i < dh->ndh; i++)
 182     {
 183         unsigned int bin;
 184
 185         /* Determine the bin number. If it doesn't fit into the histogram,
 186            add it to the last bin.
 187            We check the max_dh_int range because converting to integers
 188            might lead to overflow with unpredictable results.*/
 189         if ((f * dh->dh[i] >= min_dh_hist) && (f * dh->dh[i] <= max_dh_hist))
 190         {
 191             bin = static_cast<unsigned int>((f * dh->dh[i] - min_dh_hist) / dx);
 192         }
 193         else
 194         {
 195             bin = dh->nbins - 1;
 196         }
 197         /* double-check here because of possible round-off errors*/
 198         if (bin >= dh->nbins)
 199         {
 200             bin = dh->nbins - 1;
 201         }
 202         if (bin > dh->maxbin[hi])
 203         {
 204             dh->maxbin[hi] = bin;
 205         }
 206
 207         dh->bin[hi][bin]++;
 208     }
 209
 210     /* make sure we include a bin with 0 if we didn't use the full
 211        histogram width. This can then be used as an indication that
 212        all the data was binned. */
 213     if (dh->maxbin[hi] < dh->nbins - 1)
 214     {
 215         dh->maxbin[hi] += 1;
 216     }
 217 }
 218
 219
 220 static void mde_delta_h_handle_block(t_mde_delta_h* dh, t_enxblock* blk)
 221 {
 222     /* first check which type we should use: histogram or raw data */
 223     if (dh->nhist == 0)
 224     {
 225         int i;
 226
 227         /* We write raw data.
 228            Raw data consists of 3 subblocks: an int metadata block
 229            with type and derivative index, a foreign lambda block
 230            and and the data itself */
 231         add_subblocks_enxblock(blk, 3);
 232
 233         blk->id = enxDH;
 234
 235         /* subblock 1 */
 236         dh->subblock_meta_i[0] = dh->type;       /* block data type */
 237         dh->subblock_meta_i[1] = dh->derivative; /* derivative direction if
 238                                                     applicable (in indices
 239                                                     starting from first coord in
 240                                                     the main delta_h_coll) */
 241         blk->sub[0].nr   = 2;
 242         blk->sub[0].type = XdrDataType::Int;
 243         blk->sub[0].ival = dh->subblock_meta_i.data();
 244
 245         /* subblock 2 */
 246         for (i = 0; i < dh->nlambda; i++)
 247         {
 248             dh->subblock_meta_d[i] = dh->lambda[i];
 249         }
 250         blk->sub[1].nr   = dh->nlambda;
 251         blk->sub[1].type = XdrDataType::Double;
 252         blk->sub[1].dval = dh->subblock_meta_d.data();
 253
 254         /* subblock 3 */
 255         /* check if there's actual data to be written. */
 256         /*if (dh->ndh > 1)*/
 257         if (dh->ndh > 0)
 258         {
 259             unsigned int i;
 260
 261             blk->sub[2].nr = dh->ndh;
 262             /* Michael commented in 2012 that this use of explicit
 263                XdrDataType::Float was good for F@H for now.
 264                Apparently it's still good enough. */
 265             blk->sub[2].type = XdrDataType::Float;
 266             for (i = 0; i < dh->ndh; i++)
 267             {
 268                 dh->dhf[i] = static_cast<float>(dh->dh[i]);
 269             }
 270             blk->sub[2].fval = dh->dhf.data();
 271             dh->written      = TRUE;
 272         }
 273         else
 274         {
 275             blk->sub[2].nr   = 0;
 276             blk->sub[2].type = XdrDataType::Float;
 277             blk->sub[2].fval = nullptr;
 278         }
 279     }
 280     else
 281     {
 282         int nhist_written = 0;
 283         int i;
 284         int k;
 285
 286         /* TODO histogram metadata */
 287         /* check if there's actual data to be written. */
 288         if (dh->ndh > 1)
 289         {
 290             gmx_bool prev_complete = FALSE;
 291             /* Make the histogram(s) */
 292             for (i = 0; i < dh->nhist; i++)
 293             {
 294                 if (!prev_complete)
 295                 {
 296                     /* the first histogram is always normal, and the
 297                        second one is always reverse */
 298                     mde_delta_h_make_hist(dh, i, i == 1);
 299                     nhist_written++;
 300                     /* check whether this histogram contains all data: if the
 301                        last bin is 0, it does */
 302                     if (dh->bin[i][dh->nbins - 1] == 0)
 303                     {
 304                         prev_complete = TRUE;
 305                     }
 306                     if (!dh->derivative)
 307                     {
 308                         prev_complete = TRUE;
 309                     }
 310                 }
 311             }
 312             dh->written = TRUE;
 313         }
 314
 315         /* A histogram consists of 2, 3 or 4 subblocks:
 316            the foreign lambda value + histogram spacing, the starting point,
 317            and the histogram data (0, 1 or 2 blocks). */
 318         add_subblocks_enxblock(blk, nhist_written + 2);
 319         blk->id = enxDHHIST;
 320
 321         /* subblock 1: the lambda value + the histogram spacing */
 322         if (dh->nlambda == 1)
 323         {
 324             /* for backward compatibility */
 325             dh->subblock_meta_d[0] = dh->lambda[0];
 326         }
 327         else
 328         {
 329             dh->subblock_meta_d[0] = -1;
 330             for (i = 0; i < dh->nlambda; i++)
 331             {
 332                 dh->subblock_meta_d[2 + i] = dh->lambda[i];
 333             }
 334         }
 335         dh->subblock_meta_d[1] = dh->dx;
 336         blk->sub[0].nr         = 2 + ((dh->nlambda > 1) ? dh->nlambda : 0);
 337         blk->sub[0].type       = XdrDataType::Double;
 338         blk->sub[0].dval       = dh->subblock_meta_d.data();
 339
 340         /* subblock 2: the starting point(s) as a long integer */
 341         dh->subblock_meta_l[0] = nhist_written;
 342         dh->subblock_meta_l[1] = dh->type; /*dh->derivative ? 1 : 0;*/
 343         k                      = 2;
 344         for (i = 0; i < nhist_written; i++)
 345         {
 346             dh->subblock_meta_l[k++] = dh->x0[i];
 347         }
 348         /* append the derivative data */
 349         dh->subblock_meta_l[k++] = dh->derivative;
 350
 351         blk->sub[1].nr   = nhist_written + 3;
 352         blk->sub[1].type = XdrDataType::Int64;
 353         blk->sub[1].lval = dh->subblock_meta_l.data();
 354
 355         /* subblock 3 + 4 : the histogram data */
 356         for (i = 0; i < nhist_written; i++)
 357         {
 358             blk->sub[i + 2].nr = dh->maxbin[i] + 1; /* it's +1 because size=index+1
 359                                                        in C */
 360             blk->sub[i + 2].type = XdrDataType::Int;
 361             blk->sub[i + 2].ival = dh->bin[i].data();
 362         }
 363     }
 364 }
 365
 366 /* initialize the collection*/
 367 t_mde_delta_h_coll::t_mde_delta_h_coll(const t_inputrec& inputrec)
 368 {
 369     int       i, j, n;
 370     double*   lambda_vec;
 371     int       ndhmax = inputrec.nstenergy / inputrec.nstcalcenergy;
 372     t_lambda* fep    = inputrec.fepvals.get();
 373
 374     temperature    = inputrec.opts.ref_t[0]; /* only store system temperature */
 375     start_time     = 0.;
 376     delta_time     = inputrec.delta_t * inputrec.fepvals->nstdhdl;
 377     start_time_set = FALSE;
 378
 379     /* this is the compatibility lambda value. If it is >=0, it is valid,
 380        and there is either an old-style lambda or a slow growth simulation. */
 381     start_lambda = inputrec.fepvals->init_lambda;
 382     /* for continuous change of lambda values */
 383     delta_lambda = inputrec.fepvals->delta_lambda * inputrec.fepvals->nstdhdl;
 384
 385     if (start_lambda < 0)
 386     {
 387         /* create the native lambda vectors */
 388         lambda_index = fep->init_fep_state;
 389         n_lambda_vec = 0;
 390         for (auto i : keysOf(fep->separate_dvdl))
 391         {
 392             if (fep->separate_dvdl[i])
 393             {
 394                 n_lambda_vec++;
 395             }
 396         }
 397         native_lambda_vec.resize(n_lambda_vec);
 398         native_lambda_components.resize(n_lambda_vec);
 399         j = 0;
 400         for (auto i : keysOf(fep->separate_dvdl))
 401         {
 402             if (fep->separate_dvdl[i])
 403             {
 404                 native_lambda_components[j] = static_cast<int>(i);
 405                 if (fep->init_fep_state >= 0 && fep->init_fep_state < fep->n_lambda)
 406                 {
 407                     native_lambda_vec[j] = fep->all_lambda[i][fep->init_fep_state];
 408                 }
 409                 else
 410                 {
 411                     native_lambda_vec[j] = -1;
 412                 }
 413                 j++;
 414             }
 415         }
 416     }
 417     else
 418     {
 419         /* don't allocate the meta-data subblocks for lambda vectors */
 420         n_lambda_vec = 0;
 421         lambda_index = -1;
 422     }
 423     /* allocate metadata subblocks */
 424     subblock_d.resize(c_subblockDNumPreEntries + n_lambda_vec);
 425     subblock_i.resize(c_subblockINumPreEntries + n_lambda_vec);
 426
 427     /* now decide which data to write out */
 428     nlambda           = 0;
 429     ndhdl             = 0;
 430     dh_expanded_index = -1;
 431     dh_energy_index   = -1;
 432     dh_pv_index       = -1;
 433
 434     /* total number of raw data point collections in the sample */
 435     ndh = 0;
 436
 437     {
 438         gmx_bool bExpanded           = FALSE;
 439         gmx_bool bEnergy             = FALSE;
 440         gmx_bool bPV                 = FALSE;
 441         int      n_lambda_components = 0;
 442
 443         /* first count the number of states */
 444
 445         /* add the dhdl's */
 446         if (fep->dhdl_derivatives == DhDlDerivativeCalculation::Yes)
 447         {
 448             for (auto i : keysOf(fep->separate_dvdl))
 449             {
 450                 if (inputrec.fepvals->separate_dvdl[i])
 451                 {
 452                     ndh += 1;
 453                     ndhdl += 1;
 454                 }
 455             }
 456         }
 457         /* add the lambdas */
 458         nlambda = inputrec.fepvals->lambda_stop_n - inputrec.fepvals->lambda_start_n;
 459         ndh += nlambda;
 460         /* another compatibility check */
 461         if (start_lambda < 0)
 462         {
 463             /* include one more for the specification of the state, by lambda or
 464                fep_state*/
 465             if (inputrec.expandedvals->elmcmove > LambdaMoveCalculation::No)
 466             {
 467                 ndh += 1;
 468                 bExpanded = TRUE;
 469             }
 470             /* whether to print energies */
 471             if (inputrec.fepvals->edHdLPrintEnergy != FreeEnergyPrintEnergy::No)
 472             {
 473                 ndh += 1;
 474                 bEnergy = TRUE;
 475             }
 476             if (inputrec.epc > PressureCoupling::No)
 477             {
 478                 ndh += 1; /* include pressure-volume work */
 479                 bPV = TRUE;
 480             }
 481         }
 482         /* allocate them */
 483         dh.resize(ndh);
 484
 485         /* now initialize them */
 486         /* the order, for now, must match that of the dhdl.xvg file because of
 487            how g_energy -odh is implemented */
 488         n = 0;
 489         if (bExpanded)
 490         {
 491             dh_expanded_index = n;
 492             mde_delta_h_init(&dh[n],
 493                              inputrec.fepvals->dh_hist_size,
 494                              inputrec.fepvals->dh_hist_spacing,
 495                              ndhmax,
 496                              dhbtEXPANDED,
 497                              0,
 498                              0,
 499                              nullptr);
 500             n++;
 501         }
 502         if (bEnergy)
 503         {
 504             dh_energy_index = n;
 505             mde_delta_h_init(&dh[n],
 506                              inputrec.fepvals->dh_hist_size,
 507                              inputrec.fepvals->dh_hist_spacing,
 508                              ndhmax,
 509                              dhbtEN,
 510                              0,
 511                              0,
 512                              nullptr);
 513             n++;
 514         }
 515         /* add the dhdl's */
 516         n_lambda_components = 0;
 517         if (fep->dhdl_derivatives == DhDlDerivativeCalculation::Yes)
 518         {
 519             dh_dhdl_index = n;
 520             for (auto i : keysOf(fep->separate_dvdl))
 521             {
 522                 if (inputrec.fepvals->separate_dvdl[i])
 523                 {
 524                     /* we give it init_lambda for compatibility */
 525                     mde_delta_h_init(&dh[n],
 526                                      inputrec.fepvals->dh_hist_size,
 527                                      inputrec.fepvals->dh_hist_spacing,
 528                                      ndhmax,
 529                                      dhbtDHDL,
 530                                      n_lambda_components,
 531                                      1,
 532                                      &(fep->init_lambda));
 533                     n++;
 534                     n_lambda_components++;
 535                 }
 536             }
 537         }
 538         else
 539         {
 540             for (auto i : keysOf(fep->separate_dvdl))
 541             {
 542                 if (fep->separate_dvdl[i])
 543                 {
 544                     n_lambda_components++; /* count the components */
 545                 }
 546             }
 547         }
 548         /* add the lambdas */
 549         dh_du_index = n;
 550         snew(lambda_vec, n_lambda_components);
 551         for (i = inputrec.fepvals->lambda_start_n; i < inputrec.fepvals->lambda_stop_n; i++)
 552         {
 553             int k = 0;
 554
 555             for (auto j : keysOf(fep->separate_dvdl))
 556             {
 557                 if (fep->separate_dvdl[j])
 558                 {
 559                     lambda_vec[k++] = fep->all_lambda[j][i];
 560                 }
 561             }
 562
 563             mde_delta_h_init(&dh[n],
 564                              inputrec.fepvals->dh_hist_size,
 565                              inputrec.fepvals->dh_hist_spacing,
 566                              ndhmax,
 567                              dhbtDH,
 568                              0,
 569                              n_lambda_components,
 570                              lambda_vec);
 571             n++;
 572         }
 573         sfree(lambda_vec);
 574         if (bPV)
 575         {
 576             dh_pv_index = n;
 577             mde_delta_h_init(&dh[n],
 578                              inputrec.fepvals->dh_hist_size,
 579                              inputrec.fepvals->dh_hist_spacing,
 580                              ndhmax,
 581                              dhbtPV,
 582                              0,
 583                              0,
 584                              nullptr);
 585             n++;
 586         }
 587     }
 588 }
 589
 590 /* add a bunch of samples - note fep_state is double to allow for better data storage */
 591 void mde_delta_h_coll_add_dh(t_mde_delta_h_coll*   dhc,
 592                              double                fep_state,
 593                              double                energy,
 594                              double                pV,
 595                              gmx::ArrayRef<double> dhdl,
 596                              double*               foreign_dU,
 597                              double                time)
 598 {
 599     int i;
 600
 601     if (!dhc->start_time_set)
 602     {
 603         dhc->start_time_set = TRUE;
 604         dhc->start_time     = time;
 605     }
 606
 607     for (i = 0; i < dhc->ndhdl; i++)
 608     {
 609         mde_delta_h_add_dh(&dhc->dh[dhc->dh_dhdl_index + i], dhdl[i]);
 610     }
 611     for (i = 0; i < dhc->nlambda; i++)
 612     {
 613         mde_delta_h_add_dh(&dhc->dh[dhc->dh_du_index + i], foreign_dU[i]);
 614     }
 615     if (dhc->dh_pv_index >= 0)
 616     {
 617         mde_delta_h_add_dh(&dhc->dh[dhc->dh_pv_index], pV);
 618     }
 619     if (dhc->dh_energy_index >= 0)
 620     {
 621         mde_delta_h_add_dh(&dhc->dh[dhc->dh_energy_index], energy);
 622     }
 623     if (dhc->dh_expanded_index >= 0)
 624     {
 625         mde_delta_h_add_dh(&dhc->dh[dhc->dh_expanded_index], fep_state);
 626     }
 627 }
 628
 629 /* write the metadata associated with all the du blocks, and call
 630    handle_block to write out all the du blocks */
 631 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll* dhc, t_enxframe* fr, int nblock)
 632 {
 633     int         i;
 634     t_enxblock* blk;
 635
 636     /* add one block with one subblock as the collection's own data */
 637     nblock++;
 638     add_blocks_enxframe(fr, nblock);
 639     blk = fr->block + (nblock - 1);
 640
 641     /* only allocate lambda vector component blocks if they must be written out
 642        for backward compatibility */
 643     if (!dhc->native_lambda_components.empty())
 644     {
 645         add_subblocks_enxblock(blk, 2);
 646     }
 647     else
 648     {
 649         add_subblocks_enxblock(blk, 1);
 650     }
 651
 652     dhc->subblock_d[0] = dhc->temperature;  /* temperature */
 653     dhc->subblock_d[1] = dhc->start_time;   /* time of first sample */
 654     dhc->subblock_d[2] = dhc->delta_time;   /* time difference between samples */
 655     dhc->subblock_d[3] = dhc->start_lambda; /* old-style lambda at starttime */
 656     dhc->subblock_d[4] = dhc->delta_lambda; /* lambda diff. between samples */
 657     /* set the lambda vector components if they exist */
 658     if (!dhc->native_lambda_components.empty())
 659     {
 660         for (i = 0; i < dhc->n_lambda_vec; i++)
 661         {
 662             dhc->subblock_d[c_subblockDNumPreEntries + i] = dhc->native_lambda_vec[i];
 663         }
 664     }
 665     blk->id          = enxDHCOLL;
 666     blk->sub[0].nr   = c_subblockDNumPreEntries + dhc->n_lambda_vec;
 667     blk->sub[0].type = XdrDataType::Double;
 668     blk->sub[0].dval = dhc->subblock_d.data();
 669
 670     if (!dhc->native_lambda_components.empty())
 671     {
 672         dhc->subblock_i[0] = dhc->lambda_index;
 673         /* set the lambda vector component IDs if they exist */
 674         dhc->subblock_i[1] = dhc->n_lambda_vec;
 675         for (i = 0; i < dhc->n_lambda_vec; i++)
 676         {
 677             dhc->subblock_i[c_subblockINumPreEntries + i] = dhc->native_lambda_components[i];
 678         }
 679         blk->sub[1].nr   = c_subblockINumPreEntries + dhc->n_lambda_vec;
 680         blk->sub[1].type = XdrDataType::Int;
 681         blk->sub[1].ival = dhc->subblock_i.data();
 682     }
 683
 684     for (i = 0; i < dhc->ndh; i++)
 685     {
 686         nblock++;
 687         add_blocks_enxframe(fr, nblock);
 688         blk = fr->block + (nblock - 1);
 689
 690         mde_delta_h_handle_block(&dhc->dh[i], blk);
 691     }
 692 }
 693
 694 /* reset the data for a new round */
 695 void mde_delta_h_coll_reset(t_mde_delta_h_coll* dhc)
 696 {
 697     int i;
 698     for (i = 0; i < dhc->ndh; i++)
 699     {
 700         if (dhc->dh[i].written)
 701         {
 702             /* we can now throw away the data */
 703             mde_delta_h_reset(&dhc->dh[i]);
 704         }
 705     }
 706     dhc->start_time_set = FALSE;
 707 }
 708
 709 /* set the energyhistory variables to save state */
 710 void mde_delta_h_coll_update_energyhistory(const t_mde_delta_h_coll* dhc, energyhistory_t* enerhist)
 711 {
 712     if (enerhist->deltaHForeignLambdas == nullptr)
 713     {
 714         enerhist->deltaHForeignLambdas = std::make_unique<delta_h_history_t>();
 715         enerhist->deltaHForeignLambdas->dh.resize(dhc->ndh);
 716     }
 717
 718     delta_h_history_t* const deltaH = enerhist->deltaHForeignLambdas.get();
 719
 720     GMX_RELEASE_ASSERT(
 721             deltaH->dh.size() == static_cast<size_t>(dhc->ndh),
 722             "energy history number of delta_h histograms should match inputrec's number");
 723
 724     for (int i = 0; i < dhc->ndh; i++)
 725     {
 726         std::vector<real>& dh = deltaH->dh[i];
 727         for (unsigned int j = 0; j < dhc->dh[i].ndh; j++)
 728         {
 729             dh.emplace_back(dhc->dh[i].dh[j]);
 730         }
 731     }
 732     deltaH->start_time   = dhc->start_time;
 733     deltaH->start_lambda = dhc->start_lambda;
 734 }
 735
 736
 737 /* restore the variables from an energyhistory */
 738 void mde_delta_h_coll_restore_energyhistory(t_mde_delta_h_coll* dhc, const delta_h_history_t* deltaH)
 739 {
 740     GMX_RELEASE_ASSERT(dhc, "Should have delta_h histograms");
 741     GMX_RELEASE_ASSERT(deltaH, "Should have delta_h histograms in energy history");
 742     GMX_RELEASE_ASSERT(
 743             deltaH->dh.size() == static_cast<size_t>(dhc->ndh),
 744             "energy history number of delta_h histograms should match inputrec's number");
 745
 746     for (gmx::index i = 0; i < gmx::ssize(deltaH->dh); i++)
 747     {
 748         dhc->dh[i].ndh = deltaH->dh[i].size();
 749         for (unsigned int j = 0; j < dhc->dh[i].ndh; j++)
 750         {
 751             dhc->dh[i].dh[j] = deltaH->dh[i][j];
 752         }
 753     }
 754     dhc->start_time = deltaH->start_time;
 755     if (deltaH->start_lambda_set)
 756     {
 757         dhc->start_lambda = deltaH->start_lambda;
 758     }
 759     dhc->start_time_set = (dhc->dh[0].ndh > 0);
 760 }