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