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