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