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