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