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