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