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     snew(dh->lambda, nlambda);
  89     for (i = 0; i < nlambda; i++)
  90     {
  91         assert(lambda);
  92         dh->lambda[i] = lambda[i];
  93     }
  94
  95
  96     snew(dh->subblock_meta_d, dh->nlambda + 1);
  97
  98     dh->ndhmax = ndhmax + 2;
  99     for (i = 0; i < 2; i++)
 100     {
 101         dh->bin[i] = nullptr;
 102     }
 103
 104     snew(dh->dh, dh->ndhmax);
 105     snew(dh->dhf, dh->ndhmax);
 106
 107     if (nbins <= 0 || dx < GMX_REAL_EPS * 10)
 108     {
 109         dh->nhist = 0;
 110     }
 111     else
 112     {
 113         int i;
 114         /* pre-allocate the histogram */
 115         dh->nhist = 2; /* energies and derivatives histogram */
 116         dh->dx    = dx;
 117         dh->nbins = nbins;
 118         for (i = 0; i < dh->nhist; i++)
 119         {
 120             snew(dh->bin[i], dh->nbins);
 121         }
 122     }
 123     mde_delta_h_reset(dh);
 124 }
 125
 126 static void done_mde_delta_h(t_mde_delta_h* dh)
 127 {
 128     sfree(dh->lambda);
 129     sfree(dh->subblock_meta_d);
 130     sfree(dh->dh);
 131     sfree(dh->dhf);
 132     for (int i = 0; i < dh->nhist; i++)
 133     {
 134         sfree(dh->bin[i]);
 135     }
 136 }
 137
 138 /* Add a value to the delta_h list */
 139 static void mde_delta_h_add_dh(t_mde_delta_h* dh, double delta_h)
 140 {
 141     if (dh->ndh >= dh->ndhmax)
 142     {
 143         gmx_incons("delta_h array not big enough!");
 144     }
 145     dh->dh[dh->ndh] = delta_h;
 146     dh->ndh++;
 147 }
 148
 149 /* construct histogram with index hi */
 150 static void mde_delta_h_make_hist(t_mde_delta_h* dh, int hi, gmx_bool invert)
 151 {
 152     double       min_dh = FLT_MAX;
 153     double       max_dh = -FLT_MAX;
 154     unsigned int i;
 155     double       max_dh_hist; /* maximum binnable dh value */
 156     double       min_dh_hist; /* minimum binnable dh value */
 157     double       dx = dh->dx;
 158     double       f; /* energy mult. factor */
 159
 160     /* by applying a -1 scaling factor on the energies we get the same as
 161        having a negative dx, but we don't need to fix the min/max values
 162        beyond inverting x0 */
 163     f = invert ? -1 : 1;
 164
 165     /* first find min and max */
 166     for (i = 0; i < dh->ndh; i++)
 167     {
 168         if (f * dh->dh[i] < min_dh)
 169         {
 170             min_dh = f * dh->dh[i];
 171         }
 172         if (f * dh->dh[i] > max_dh)
 173         {
 174             max_dh = f * dh->dh[i];
 175         }
 176     }
 177
 178     /* reset the histogram */
 179     for (i = 0; i < dh->nbins; i++)
 180     {
 181         dh->bin[hi][i] = 0;
 182     }
 183     dh->maxbin[hi] = 0;
 184
 185     /* The starting point of the histogram is the lowest value found:
 186        that value has the highest contribution to the free energy.
 187
 188        Get this start value in number of histogram dxs from zero,
 189        as an integer.*/
 190
 191     dh->x0[hi] = static_cast<int64_t>(floor(min_dh / dx));
 192
 193     min_dh_hist = (dh->x0[hi]) * dx;
 194     max_dh_hist = (dh->x0[hi] + dh->nbins + 1) * dx;
 195
 196     /* and fill the histogram*/
 197     for (i = 0; i < dh->ndh; i++)
 198     {
 199         unsigned int bin;
 200
 201         /* Determine the bin number. If it doesn't fit into the histogram,
 202            add it to the last bin.
 203            We check the max_dh_int range because converting to integers
 204            might lead to overflow with unpredictable results.*/
 205         if ((f * dh->dh[i] >= min_dh_hist) && (f * dh->dh[i] <= max_dh_hist))
 206         {
 207             bin = static_cast<unsigned int>((f * dh->dh[i] - min_dh_hist) / dx);
 208         }
 209         else
 210         {
 211             bin = dh->nbins - 1;
 212         }
 213         /* double-check here because of possible round-off errors*/
 214         if (bin >= dh->nbins)
 215         {
 216             bin = dh->nbins - 1;
 217         }
 218         if (bin > dh->maxbin[hi])
 219         {
 220             dh->maxbin[hi] = bin;
 221         }
 222
 223         dh->bin[hi][bin]++;
 224     }
 225
 226     /* make sure we include a bin with 0 if we didn't use the full
 227        histogram width. This can then be used as an indication that
 228        all the data was binned. */
 229     if (dh->maxbin[hi] < dh->nbins - 1)
 230     {
 231         dh->maxbin[hi] += 1;
 232     }
 233 }
 234
 235
 236 static void mde_delta_h_handle_block(t_mde_delta_h* dh, t_enxblock* blk)
 237 {
 238     /* first check which type we should use: histogram or raw data */
 239     if (dh->nhist == 0)
 240     {
 241         int i;
 242
 243         /* We write raw data.
 244            Raw data consists of 3 subblocks: an int metadata block
 245            with type and derivative index, a foreign lambda block
 246            and and the data itself */
 247         add_subblocks_enxblock(blk, 3);
 248
 249         blk->id = enxDH;
 250
 251         /* subblock 1 */
 252         dh->subblock_meta_i[0] = dh->type;       /* block data type */
 253         dh->subblock_meta_i[1] = dh->derivative; /* derivative direction if
 254                                                     applicable (in indices
 255                                                     starting from first coord in
 256                                                     the main delta_h_coll) */
 257         blk->sub[0].nr   = 2;
 258         blk->sub[0].type = xdr_datatype_int;
 259         blk->sub[0].ival = dh->subblock_meta_i;
 260
 261         /* subblock 2 */
 262         for (i = 0; i < dh->nlambda; i++)
 263         {
 264             dh->subblock_meta_d[i] = dh->lambda[i];
 265         }
 266         blk->sub[1].nr   = dh->nlambda;
 267         blk->sub[1].type = xdr_datatype_double;
 268         blk->sub[1].dval = dh->subblock_meta_d;
 269
 270         /* subblock 3 */
 271         /* check if there's actual data to be written. */
 272         /*if (dh->ndh > 1)*/
 273         if (dh->ndh > 0)
 274         {
 275             unsigned int i;
 276
 277             blk->sub[2].nr = dh->ndh;
 278             /* Michael commented in 2012 that this use of explicit
 279                xdr_datatype_float was good for F@H for now.
 280                Apparently it's still good enough. */
 281             blk->sub[2].type = xdr_datatype_float;
 282             for (i = 0; i < dh->ndh; i++)
 283             {
 284                 dh->dhf[i] = static_cast<float>(dh->dh[i]);
 285             }
 286             blk->sub[2].fval = dh->dhf;
 287             dh->written      = TRUE;
 288         }
 289         else
 290         {
 291             blk->sub[2].nr   = 0;
 292             blk->sub[2].type = xdr_datatype_float;
 293             blk->sub[2].fval = nullptr;
 294         }
 295     }
 296     else
 297     {
 298         int nhist_written = 0;
 299         int i;
 300         int k;
 301
 302         /* TODO histogram metadata */
 303         /* check if there's actual data to be written. */
 304         if (dh->ndh > 1)
 305         {
 306             gmx_bool prev_complete = FALSE;
 307             /* Make the histogram(s) */
 308             for (i = 0; i < dh->nhist; i++)
 309             {
 310                 if (!prev_complete)
 311                 {
 312                     /* the first histogram is always normal, and the
 313                        second one is always reverse */
 314                     mde_delta_h_make_hist(dh, i, i == 1);
 315                     nhist_written++;
 316                     /* check whether this histogram contains all data: if the
 317                        last bin is 0, it does */
 318                     if (dh->bin[i][dh->nbins - 1] == 0)
 319                     {
 320                         prev_complete = TRUE;
 321                     }
 322                     if (!dh->derivative)
 323                     {
 324                         prev_complete = TRUE;
 325                     }
 326                 }
 327             }
 328             dh->written = TRUE;
 329         }
 330
 331         /* A histogram consists of 2, 3 or 4 subblocks:
 332            the foreign lambda value + histogram spacing, the starting point,
 333            and the histogram data (0, 1 or 2 blocks). */
 334         add_subblocks_enxblock(blk, nhist_written + 2);
 335         blk->id = enxDHHIST;
 336
 337         /* subblock 1: the lambda value + the histogram spacing */
 338         if (dh->nlambda == 1)
 339         {
 340             /* for backward compatibility */
 341             dh->subblock_meta_d[0] = dh->lambda[0];
 342         }
 343         else
 344         {
 345             dh->subblock_meta_d[0] = -1;
 346             for (i = 0; i < dh->nlambda; i++)
 347             {
 348                 dh->subblock_meta_d[2 + i] = dh->lambda[i];
 349             }
 350         }
 351         dh->subblock_meta_d[1] = dh->dx;
 352         blk->sub[0].nr         = 2 + ((dh->nlambda > 1) ? dh->nlambda : 0);
 353         blk->sub[0].type       = xdr_datatype_double;
 354         blk->sub[0].dval       = dh->subblock_meta_d;
 355
 356         /* subblock 2: the starting point(s) as a long integer */
 357         dh->subblock_meta_l[0] = nhist_written;
 358         dh->subblock_meta_l[1] = dh->type; /*dh->derivative ? 1 : 0;*/
 359         k                      = 2;
 360         for (i = 0; i < nhist_written; i++)
 361         {
 362             dh->subblock_meta_l[k++] = dh->x0[i];
 363         }
 364         /* append the derivative data */
 365         dh->subblock_meta_l[k++] = dh->derivative;
 366
 367         blk->sub[1].nr   = nhist_written + 3;
 368         blk->sub[1].type = xdr_datatype_int64;
 369         blk->sub[1].lval = dh->subblock_meta_l;
 370
 371         /* subblock 3 + 4 : the histogram data */
 372         for (i = 0; i < nhist_written; i++)
 373         {
 374             blk->sub[i + 2].nr = dh->maxbin[i] + 1; /* it's +1 because size=index+1
 375                                                        in C */
 376             blk->sub[i + 2].type = xdr_datatype_int;
 377             blk->sub[i + 2].ival = dh->bin[i];
 378         }
 379     }
 380 }
 381
 382 /* initialize the collection*/
 383 void mde_delta_h_coll_init(t_mde_delta_h_coll* dhc, const t_inputrec* ir)
 384 {
 385     int       i, j, n;
 386     double*   lambda_vec;
 387     int       ndhmax = ir->nstenergy / ir->nstcalcenergy;
 388     t_lambda* fep    = ir->fepvals.get();
 389
 390     dhc->temperature    = ir->opts.ref_t[0]; /* only store system temperature */
 391     dhc->start_time     = 0.;
 392     dhc->delta_time     = ir->delta_t * ir->fepvals->nstdhdl;
 393     dhc->start_time_set = FALSE;
 394
 395     /* this is the compatibility lambda value. If it is >=0, it is valid,
 396        and there is either an old-style lambda or a slow growth simulation. */
 397     dhc->start_lambda = ir->fepvals->init_lambda;
 398     /* for continuous change of lambda values */
 399     dhc->delta_lambda = ir->fepvals->delta_lambda * ir->fepvals->nstdhdl;
 400
 401     if (dhc->start_lambda < 0)
 402     {
 403         /* create the native lambda vectors */
 404         dhc->lambda_index = fep->init_fep_state;
 405         dhc->n_lambda_vec = 0;
 406         for (auto i : keysOf(fep->separate_dvdl))
 407         {
 408             if (fep->separate_dvdl[i])
 409             {
 410                 dhc->n_lambda_vec++;
 411             }
 412         }
 413         snew(dhc->native_lambda_vec, dhc->n_lambda_vec);
 414         snew(dhc->native_lambda_components, dhc->n_lambda_vec);
 415         j = 0;
 416         for (auto i : keysOf(fep->separate_dvdl))
 417         {
 418             if (fep->separate_dvdl[i])
 419             {
 420                 dhc->native_lambda_components[j] = static_cast<int>(i);
 421                 if (fep->init_fep_state >= 0 && fep->init_fep_state < fep->n_lambda)
 422                 {
 423                     dhc->native_lambda_vec[j] = fep->all_lambda[i][fep->init_fep_state];
 424                 }
 425                 else
 426                 {
 427                     dhc->native_lambda_vec[j] = -1;
 428                 }
 429                 j++;
 430             }
 431         }
 432     }
 433     else
 434     {
 435         /* don't allocate the meta-data subblocks for lambda vectors */
 436         dhc->native_lambda_vec        = nullptr;
 437         dhc->n_lambda_vec             = 0;
 438         dhc->native_lambda_components = nullptr;
 439         dhc->lambda_index             = -1;
 440     }
 441     /* allocate metadata subblocks */
 442     snew(dhc->subblock_d, c_subblockDNumPreEntries + dhc->n_lambda_vec);
 443     snew(dhc->subblock_i, c_subblockINumPreEntries + dhc->n_lambda_vec);
 444
 445     /* now decide which data to write out */
 446     dhc->nlambda     = 0;
 447     dhc->ndhdl       = 0;
 448     dhc->dh_expanded = nullptr;
 449     dhc->dh_energy   = nullptr;
 450     dhc->dh_pv       = nullptr;
 451
 452     /* total number of raw data point collections in the sample */
 453     dhc->ndh = 0;
 454
 455     {
 456         gmx_bool bExpanded           = FALSE;
 457         gmx_bool bEnergy             = FALSE;
 458         gmx_bool bPV                 = FALSE;
 459         int      n_lambda_components = 0;
 460
 461         /* first count the number of states */
 462
 463         /* add the dhdl's */
 464         if (fep->dhdl_derivatives == DhDlDerivativeCalculation::Yes)
 465         {
 466             for (auto i : keysOf(fep->separate_dvdl))
 467             {
 468                 if (ir->fepvals->separate_dvdl[i])
 469                 {
 470                     dhc->ndh += 1;
 471                     dhc->ndhdl += 1;
 472                 }
 473             }
 474         }
 475         /* add the lambdas */
 476         dhc->nlambda = ir->fepvals->lambda_stop_n - ir->fepvals->lambda_start_n;
 477         dhc->ndh += dhc->nlambda;
 478         /* another compatibility check */
 479         if (dhc->start_lambda < 0)
 480         {
 481             /* include one more for the specification of the state, by lambda or
 482                fep_state*/
 483             if (ir->expandedvals->elmcmove > LambdaMoveCalculation::No)
 484             {
 485                 dhc->ndh += 1;
 486                 bExpanded = TRUE;
 487             }
 488             /* whether to print energies */
 489             if (ir->fepvals->edHdLPrintEnergy != FreeEnergyPrintEnergy::No)
 490             {
 491                 dhc->ndh += 1;
 492                 bEnergy = TRUE;
 493             }
 494             if (ir->epc > PressureCoupling::No)
 495             {
 496                 dhc->ndh += 1; /* include pressure-volume work */
 497                 bPV = TRUE;
 498             }
 499         }
 500         /* allocate them */
 501         snew(dhc->dh, dhc->ndh);
 502
 503         /* now initialize them */
 504         /* the order, for now, must match that of the dhdl.xvg file because of
 505            how g_energy -odh is implemented */
 506         n = 0;
 507         if (bExpanded)
 508         {
 509             dhc->dh_expanded = dhc->dh + n;
 510             mde_delta_h_init(dhc->dh + n,
 511                              ir->fepvals->dh_hist_size,
 512                              ir->fepvals->dh_hist_spacing,
 513                              ndhmax,
 514                              dhbtEXPANDED,
 515                              0,
 516                              0,
 517                              nullptr);
 518             n++;
 519         }
 520         if (bEnergy)
 521         {
 522             dhc->dh_energy = dhc->dh + n;
 523             mde_delta_h_init(
 524                     dhc->dh + n, ir->fepvals->dh_hist_size, ir->fepvals->dh_hist_spacing, ndhmax, dhbtEN, 0, 0, nullptr);
 525             n++;
 526         }
 527         /* add the dhdl's */
 528         n_lambda_components = 0;
 529         if (fep->dhdl_derivatives == DhDlDerivativeCalculation::Yes)
 530         {
 531             dhc->dh_dhdl = dhc->dh + n;
 532             for (auto i : keysOf(fep->separate_dvdl))
 533             {
 534                 if (ir->fepvals->separate_dvdl[i])
 535                 {
 536                     /* we give it init_lambda for compatibility */
 537                     mde_delta_h_init(dhc->dh + n,
 538                                      ir->fepvals->dh_hist_size,
 539                                      ir->fepvals->dh_hist_spacing,
 540                                      ndhmax,
 541                                      dhbtDHDL,
 542                                      n_lambda_components,
 543                                      1,
 544                                      &(fep->init_lambda));
 545                     n++;
 546                     n_lambda_components++;
 547                 }
 548             }
 549         }
 550         else
 551         {
 552             for (auto i : keysOf(fep->separate_dvdl))
 553             {
 554                 if (fep->separate_dvdl[i])
 555                 {
 556                     n_lambda_components++; /* count the components */
 557                 }
 558             }
 559         }
 560         /* add the lambdas */
 561         dhc->dh_du = dhc->dh + n;
 562         snew(lambda_vec, n_lambda_components);
 563         for (i = ir->fepvals->lambda_start_n; i < ir->fepvals->lambda_stop_n; i++)
 564         {
 565             int k = 0;
 566
 567             for (auto j : keysOf(fep->separate_dvdl))
 568             {
 569                 if (fep->separate_dvdl[j])
 570                 {
 571                     lambda_vec[k++] = fep->all_lambda[j][i];
 572                 }
 573             }
 574
 575             mde_delta_h_init(dhc->dh + n,
 576                              ir->fepvals->dh_hist_size,
 577                              ir->fepvals->dh_hist_spacing,
 578                              ndhmax,
 579                              dhbtDH,
 580                              0,
 581                              n_lambda_components,
 582                              lambda_vec);
 583             n++;
 584         }
 585         sfree(lambda_vec);
 586         if (bPV)
 587         {
 588             dhc->dh_pv = dhc->dh + n;
 589             mde_delta_h_init(
 590                     dhc->dh + n, ir->fepvals->dh_hist_size, ir->fepvals->dh_hist_spacing, ndhmax, dhbtPV, 0, 0, nullptr);
 591             n++;
 592         }
 593     }
 594 }
 595
 596 void done_mde_delta_h_coll(t_mde_delta_h_coll* dhc)
 597 {
 598     if (dhc == nullptr)
 599     {
 600         return;
 601     }
 602     sfree(dhc->native_lambda_vec);
 603     sfree(dhc->native_lambda_components);
 604     sfree(dhc->subblock_d);
 605     sfree(dhc->subblock_i);
 606     for (int i = 0; i < dhc->ndh; ++i)
 607     {
 608         done_mde_delta_h(&dhc->dh[i]);
 609     }
 610     sfree(dhc->dh);
 611     sfree(dhc);
 612 }
 613
 614 /* add a bunch of samples - note fep_state is double to allow for better data storage */
 615 void mde_delta_h_coll_add_dh(t_mde_delta_h_coll*   dhc,
 616                              double                fep_state,
 617                              double                energy,
 618                              double                pV,
 619                              gmx::ArrayRef<double> dhdl,
 620                              double*               foreign_dU,
 621                              double                time)
 622 {
 623     int i;
 624
 625     if (!dhc->start_time_set)
 626     {
 627         dhc->start_time_set = TRUE;
 628         dhc->start_time     = time;
 629     }
 630
 631     for (i = 0; i < dhc->ndhdl; i++)
 632     {
 633         mde_delta_h_add_dh(dhc->dh_dhdl + i, dhdl[i]);
 634     }
 635     for (i = 0; i < dhc->nlambda; i++)
 636     {
 637         mde_delta_h_add_dh(dhc->dh_du + i, foreign_dU[i]);
 638     }
 639     if (dhc->dh_pv != nullptr)
 640     {
 641         mde_delta_h_add_dh(dhc->dh_pv, pV);
 642     }
 643     if (dhc->dh_energy != nullptr)
 644     {
 645         mde_delta_h_add_dh(dhc->dh_energy, energy);
 646     }
 647     if (dhc->dh_expanded != nullptr)
 648     {
 649         mde_delta_h_add_dh(dhc->dh_expanded, fep_state);
 650     }
 651 }
 652
 653 /* write the metadata associated with all the du blocks, and call
 654    handle_block to write out all the du blocks */
 655 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll* dhc, t_enxframe* fr, int nblock)
 656 {
 657     int         i;
 658     t_enxblock* blk;
 659
 660     /* add one block with one subblock as the collection's own data */
 661     nblock++;
 662     add_blocks_enxframe(fr, nblock);
 663     blk = fr->block + (nblock - 1);
 664
 665     /* only allocate lambda vector component blocks if they must be written out
 666        for backward compatibility */
 667     if (dhc->native_lambda_components != nullptr)
 668     {
 669         add_subblocks_enxblock(blk, 2);
 670     }
 671     else
 672     {
 673         add_subblocks_enxblock(blk, 1);
 674     }
 675
 676     dhc->subblock_d[0] = dhc->temperature;  /* temperature */
 677     dhc->subblock_d[1] = dhc->start_time;   /* time of first sample */
 678     dhc->subblock_d[2] = dhc->delta_time;   /* time difference between samples */
 679     dhc->subblock_d[3] = dhc->start_lambda; /* old-style lambda at starttime */
 680     dhc->subblock_d[4] = dhc->delta_lambda; /* lambda diff. between samples */
 681     /* set the lambda vector components if they exist */
 682     if (dhc->native_lambda_components != nullptr)
 683     {
 684         for (i = 0; i < dhc->n_lambda_vec; i++)
 685         {
 686             dhc->subblock_d[c_subblockDNumPreEntries + i] = dhc->native_lambda_vec[i];
 687         }
 688     }
 689     blk->id          = enxDHCOLL;
 690     blk->sub[0].nr   = c_subblockDNumPreEntries + dhc->n_lambda_vec;
 691     blk->sub[0].type = xdr_datatype_double;
 692     blk->sub[0].dval = dhc->subblock_d;
 693
 694     if (dhc->native_lambda_components != nullptr)
 695     {
 696         dhc->subblock_i[0] = dhc->lambda_index;
 697         /* set the lambda vector component IDs if they exist */
 698         dhc->subblock_i[1] = dhc->n_lambda_vec;
 699         for (i = 0; i < dhc->n_lambda_vec; i++)
 700         {
 701             dhc->subblock_i[c_subblockINumPreEntries + i] = dhc->native_lambda_components[i];
 702         }
 703         blk->sub[1].nr   = c_subblockINumPreEntries + dhc->n_lambda_vec;
 704         blk->sub[1].type = xdr_datatype_int;
 705         blk->sub[1].ival = dhc->subblock_i;
 706     }
 707
 708     for (i = 0; i < dhc->ndh; i++)
 709     {
 710         nblock++;
 711         add_blocks_enxframe(fr, nblock);
 712         blk = fr->block + (nblock - 1);
 713
 714         mde_delta_h_handle_block(dhc->dh + i, blk);
 715     }
 716 }
 717
 718 /* reset the data for a new round */
 719 void mde_delta_h_coll_reset(t_mde_delta_h_coll* dhc)
 720 {
 721     int i;
 722     for (i = 0; i < dhc->ndh; i++)
 723     {
 724         if (dhc->dh[i].written)
 725         {
 726             /* we can now throw away the data */
 727             mde_delta_h_reset(dhc->dh + i);
 728         }
 729     }
 730     dhc->start_time_set = FALSE;
 731 }
 732
 733 /* set the energyhistory variables to save state */
 734 void mde_delta_h_coll_update_energyhistory(const t_mde_delta_h_coll* dhc, energyhistory_t* enerhist)
 735 {
 736     if (enerhist->deltaHForeignLambdas == nullptr)
 737     {
 738         enerhist->deltaHForeignLambdas = std::make_unique<delta_h_history_t>();
 739         enerhist->deltaHForeignLambdas->dh.resize(dhc->ndh);
 740     }
 741
 742     delta_h_history_t* const deltaH = enerhist->deltaHForeignLambdas.get();
 743
 744     GMX_RELEASE_ASSERT(
 745             deltaH->dh.size() == static_cast<size_t>(dhc->ndh),
 746             "energy history number of delta_h histograms should match inputrec's number");
 747
 748     for (int i = 0; i < dhc->ndh; i++)
 749     {
 750         std::vector<real>& dh = deltaH->dh[i];
 751         dh.resize(dhc->dh[i].ndh);
 752         std::copy(dhc->dh[i].dh, dhc->dh[i].dh + dhc->dh[i].ndh, dh.begin());
 753     }
 754     deltaH->start_time   = dhc->start_time;
 755     deltaH->start_lambda = dhc->start_lambda;
 756 }
 757
 758
 759 /* restore the variables from an energyhistory */
 760 void mde_delta_h_coll_restore_energyhistory(t_mde_delta_h_coll* dhc, const delta_h_history_t* deltaH)
 761 {
 762     GMX_RELEASE_ASSERT(dhc, "Should have delta_h histograms");
 763     GMX_RELEASE_ASSERT(deltaH, "Should have delta_h histograms in energy history");
 764     GMX_RELEASE_ASSERT(
 765             deltaH->dh.size() == static_cast<size_t>(dhc->ndh),
 766             "energy history number of delta_h histograms should match inputrec's number");
 767
 768     for (gmx::index i = 0; i < gmx::ssize(deltaH->dh); i++)
 769     {
 770         dhc->dh[i].ndh = deltaH->dh[i].size();
 771         for (unsigned int j = 0; j < dhc->dh[i].ndh; j++)
 772         {
 773             dhc->dh[i].dh[j] = deltaH->dh[i][j];
 774         }
 775     }
 776     dhc->start_time = deltaH->start_time;
 777     if (deltaH->start_lambda_set)
 778     {
 779         dhc->start_lambda = deltaH->start_lambda;
 780     }
 781     dhc->start_time_set = (dhc->dh[0].ndh > 0);
 782 }