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