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