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