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