src/gromacs/mdlib/mdebin_bar.c

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