src/gromacs/mdlib/mdebin_bar.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <string.h>
  41 #include <float.h>
  42 #include <math.h>
  43 #include "typedefs.h"
  44 #include "string2.h"
  45 #include "gmx_fatal.h"
  46 #include "mdebin.h"
  47 #include "smalloc.h"
  48 #include "enxio.h"
  49 #include "gmxfio.h"
  50 #include "mdebin_bar.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, double time)
 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             min_dh=f*dh->dh[i];
 141         if (f*dh->dh[i] > max_dh)
 142             max_dh=f*dh->dh[i];
 143     }
 144
 145     /* reset the histogram */
 146     for(i=0;i<dh->nbins;i++)
 147     {
 148         dh->bin[hi][i]=0;
 149     }
 150     dh->maxbin[hi]=0;
 151
 152     /* The starting point of the histogram is the lowest value found:
 153        that value has the highest contribution to the free energy.
 154
 155        Get this start value in number of histogram dxs from zero,
 156        as an integer.*/
 157
 158     dh->x0[hi] = (gmx_large_int_t)floor(min_dh/dx);
 159
 160     min_dh_hist=(dh->x0[hi])*dx;
 161     max_dh_hist=(dh->x0[hi] + dh->nbins + 1)*dx;
 162
 163     /* and fill the histogram*/
 164     for(i=0;i<dh->ndh;i++)
 165     {
 166         unsigned int bin;
 167
 168         /* Determine the bin number. If it doesn't fit into the histogram,
 169            add it to the last bin.
 170            We check the max_dh_int range because converting to integers
 171            might lead to overflow with unpredictable results.*/
 172         if ( (f*dh->dh[i] >= min_dh_hist) && (f*dh->dh[i] <= max_dh_hist ) )
 173         {
 174             bin = (unsigned int)( (f*dh->dh[i] - min_dh_hist)/dx );
 175         }
 176         else
 177         {
 178             bin = dh->nbins-1;
 179         }
 180         /* double-check here because of possible round-off errors*/
 181         if (bin >= dh->nbins)
 182         {
 183             bin = dh->nbins-1;
 184         }
 185         if (bin > dh->maxbin[hi])
 186         {
 187             dh->maxbin[hi] = bin;
 188         }
 189
 190         dh->bin[hi][bin]++;
 191     }
 192
 193     /* make sure we include a bin with 0 if we didn't use the full
 194        histogram width. This can then be used as an indication that
 195        all the data was binned. */
 196     if (dh->maxbin[hi] < dh->nbins-1)
 197         dh->maxbin[hi] += 1;
 198 }
 199
 200
 201 void mde_delta_h_handle_block(t_mde_delta_h *dh, t_enxblock *blk)
 202 {
 203     /* first check which type we should use: histogram or raw data */
 204     if (dh->nhist == 0)
 205     {
 206         int i;
 207
 208         /* We write raw data.
 209            Raw data consists of 3 subblocks: an int metadata block
 210            with type and derivative index, a foreign lambda block
 211            and and the data itself */
 212         add_subblocks_enxblock(blk, 3);
 213
 214         blk->id=enxDH;
 215
 216         /* subblock 1 */
 217         dh->subblock_meta_i[0]=dh->type; /* block data type */
 218         dh->subblock_meta_i[1]=dh->derivative; /* derivative direction if
 219                                                   applicable (in indices
 220                                                   starting from first coord in
 221                                                   the main delta_h_coll) */
 222         blk->sub[0].nr=2;
 223         blk->sub[0].type=xdr_datatype_int;
 224         blk->sub[0].ival=dh->subblock_meta_i;
 225
 226         /* subblock 2 */
 227         for(i=0;i<dh->nlambda;i++)
 228         {
 229             dh->subblock_meta_d[i]=dh->lambda[i];
 230         }
 231         blk->sub[1].nr=dh->nlambda;
 232         blk->sub[1].type=xdr_datatype_double;
 233         blk->sub[1].dval=dh->subblock_meta_d;
 234
 235         /* subblock 3 */
 236         /* check if there's actual data to be written. */
 237         /*if (dh->ndh > 1)*/
 238         if (dh->ndh > 0)
 239         {
 240             unsigned int i;
 241
 242             blk->sub[2].nr=dh->ndh;
 243 /* For F@H for now. */
 244 #undef GMX_DOUBLE
 245 #ifndef GMX_DOUBLE
 246             blk->sub[2].type=xdr_datatype_float;
 247             for(i=0;i<dh->ndh;i++)
 248             {
 249                 dh->dhf[i] = (float)dh->dh[i];
 250             }
 251             blk->sub[2].fval=dh->dhf;
 252 #else
 253             blk->sub[2].type=xdr_datatype_double;
 254             blk->sub[2].dval=dh->dh;
 255 #endif
 256             dh->written=TRUE;
 257         }
 258         else
 259         {
 260             blk->sub[2].nr=0;
 261 #ifndef GMX_DOUBLE
 262             blk->sub[2].type=xdr_datatype_float;
 263             blk->sub[2].fval=NULL;
 264 #else
 265             blk->sub[2].type=xdr_datatype_double;
 266             blk->sub[2].dval=NULL;
 267 #endif
 268         }
 269     }
 270     else
 271     {
 272         int nhist_written=0;
 273         int i;
 274         int k;
 275
 276         /* TODO histogram metadata */
 277         /* check if there's actual data to be written. */
 278         if (dh->ndh > 1)
 279         {
 280             gmx_bool prev_complete=FALSE;
 281             /* Make the histogram(s) */
 282             for(i=0;i<dh->nhist;i++)
 283             {
 284                 if (!prev_complete)
 285                 {
 286                     /* the first histogram is always normal, and the
 287                        second one is always reverse */
 288                     mde_delta_h_make_hist(dh, i, i==1);
 289                     nhist_written++;
 290                     /* check whether this histogram contains all data: if the
 291                        last bin is 0, it does */
 292                     if (dh->bin[i][dh->nbins-1] == 0)
 293                         prev_complete=TRUE;
 294                     if (!dh->derivative)
 295                         prev_complete=TRUE;
 296                 }
 297             }
 298             dh->written=TRUE;
 299         }
 300
 301         /* A histogram consists of 2, 3 or 4 subblocks:
 302            the foreign lambda value + histogram spacing, the starting point,
 303            and the histogram data (0, 1 or 2 blocks). */
 304         add_subblocks_enxblock(blk, nhist_written+2);
 305         blk->id=enxDHHIST;
 306
 307         /* subblock 1: the lambda value + the histogram spacing */
 308         if (dh->nlambda == 1)
 309         {
 310             /* for backward compatibility */
 311             dh->subblock_meta_d[0]=dh->lambda[0];
 312         }
 313         else
 314         {
 315             dh->subblock_meta_d[0]=-1;
 316             for(i=0;i<dh->nlambda;i++)
 317             {
 318                 dh->subblock_meta_d[2+i]=dh->lambda[i];
 319             }
 320         }
 321         dh->subblock_meta_d[1]=dh->dx;
 322         blk->sub[0].nr = 2+ ((dh->nlambda>1) ? dh->nlambda : 0);
 323         blk->sub[0].type=xdr_datatype_double;
 324         blk->sub[0].dval=dh->subblock_meta_d;
 325
 326         /* subblock 2: the starting point(s) as a long integer */
 327         dh->subblock_meta_l[0]=nhist_written;
 328         dh->subblock_meta_l[1]=dh->type; /*dh->derivative ? 1 : 0;*/
 329         k=2;
 330         for(i=0;i<nhist_written;i++)
 331             dh->subblock_meta_l[k++]=dh->x0[i];
 332         /* append the derivative data */
 333         dh->subblock_meta_l[k++]=dh->derivative;
 334
 335         blk->sub[1].nr=nhist_written+3;
 336         blk->sub[1].type=xdr_datatype_large_int;
 337         blk->sub[1].lval=dh->subblock_meta_l;
 338
 339         /* subblock 3 + 4 : the histogram data */
 340         for(i=0;i<nhist_written;i++)
 341         {
 342             blk->sub[i+2].nr=dh->maxbin[i]+1; /* it's +1 because size=index+1
 343                                                  in C */
 344             blk->sub[i+2].type=xdr_datatype_int;
 345             blk->sub[i+2].ival=dh->bin[i];
 346         }
 347     }
 348 }
 349
 350 /* initialize the collection*/
 351 void mde_delta_h_coll_init(t_mde_delta_h_coll *dhc, const t_inputrec *ir)
 352 {
 353     int i,j,n;
 354     double lambda;
 355     double *lambda_vec;
 356     int ndhmax=ir->nstenergy/ir->nstcalcenergy;
 357     t_lambda *fep=ir->fepvals;
 358
 359     dhc->temperature=ir->opts.ref_t[0];  /* only store system temperature */
 360     dhc->start_time=0.;
 361     dhc->delta_time=ir->delta_t*ir->fepvals->nstdhdl;
 362     dhc->start_time_set=FALSE;
 363
 364     /* this is the compatibility lambda value. If it is >=0, it is valid,
 365        and there is either an old-style lambda or a slow growth simulation. */
 366     dhc->start_lambda=ir->fepvals->init_lambda;
 367     /* for continuous change of lambda values */
 368     dhc->delta_lambda=ir->fepvals->delta_lambda*ir->fepvals->nstdhdl;
 369
 370     if (dhc->start_lambda < 0)
 371     {
 372         /* create the native lambda vectors */
 373         dhc->lambda_index=fep->init_fep_state;
 374         dhc->n_lambda_vec=0;
 375         for(i=0;i<efptNR;i++)
 376         {
 377             if (fep->separate_dvdl[i])
 378             {
 379                 dhc->n_lambda_vec++;
 380             }
 381         }
 382         snew(dhc->native_lambda_vec, dhc->n_lambda_vec);
 383         snew(dhc->native_lambda_components, dhc->n_lambda_vec);
 384         j=0;
 385         for(i=0;i<efptNR;i++)
 386         {
 387             if (fep->separate_dvdl[i])
 388             {
 389                 dhc->native_lambda_components[j]=i;
 390                 if (fep->init_fep_state >=0 &&
 391                     fep->init_fep_state < fep->n_lambda)
 392                 {
 393                     dhc->native_lambda_vec[j]=
 394                                 fep->all_lambda[i][fep->init_fep_state];
 395                 }
 396                 else
 397                 {
 398                     dhc->native_lambda_vec[j]=-1;
 399                 }
 400                 j++;
 401             }
 402         }
 403     }
 404     else
 405     {
 406         /* don't allocate the meta-data subblocks for lambda vectors */
 407         dhc->native_lambda_vec=NULL;
 408         dhc->n_lambda_vec=0;
 409         dhc->native_lambda_components=0;
 410         dhc->lambda_index=-1;
 411     }
 412     /* allocate metadata subblocks */
 413     snew(dhc->subblock_d, 5 + dhc->n_lambda_vec);
 414     snew(dhc->subblock_i, 1 + dhc->n_lambda_vec);
 415
 416     /* now decide which data to write out */
 417     dhc->nlambda=0;
 418     dhc->ndhdl=0;
 419     dhc->dh_expanded=NULL;
 420     dhc->dh_energy=NULL;
 421     dhc->dh_pv=NULL;
 422
 423     /* total number of raw data point collections in the sample */
 424     dhc->ndh = 0;
 425
 426     {
 427         gmx_bool bExpanded=FALSE;
 428         gmx_bool bEnergy=FALSE;
 429         gmx_bool bPV=FALSE;
 430         int n_lambda_components=0;
 431
 432         /* first count the number of states */
 433
 434         /* add the dhdl's */
 435         if (fep->dhdl_derivatives == edhdlderivativesYES)
 436         {
 437             for (i=0;i<efptNR;i++)
 438             {
 439                 if (ir->fepvals->separate_dvdl[i])
 440                 {
 441                     dhc->ndh+=1;
 442                     dhc->ndhdl+=1;
 443                 }
 444             }
 445         }
 446         /* add the lambdas */
 447         dhc->nlambda = ir->fepvals->lambda_stop_n - ir->fepvals->lambda_start_n;
 448         dhc->ndh += dhc->nlambda;
 449         /* another compatibility check */
 450         if (dhc->start_lambda < 0)
 451         {
 452             /* include one more for the specification of the state, by lambda or
 453                fep_state*/
 454             if (ir->expandedvals->elmcmove > elmcmoveNO) {
 455                 dhc->ndh +=1;
 456                 bExpanded=TRUE;
 457             }
 458             /* whether to print energies */
 459             if (ir->fepvals->bPrintEnergy) {
 460                 dhc->ndh += 1;
 461                 bEnergy=TRUE;
 462             }
 463             if (ir->epc > epcNO) {
 464                 dhc->ndh += 1;  /* include pressure-volume work */
 465                 bPV=TRUE;
 466             }
 467         }
 468         /* allocate them */
 469         snew(dhc->dh, dhc->ndh);
 470
 471         /* now initialize them */
 472         /* the order, for now, must match that of the dhdl.xvg file because of
 473            how g_energy -odh is implemented */
 474         n=0;
 475         if (bExpanded)
 476         {
 477             dhc->dh_expanded=dhc->dh+n;
 478             mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
 479                              ir->fepvals->dh_hist_spacing, ndhmax,
 480                              dhbtEXPANDED, 0, 0, NULL);
 481             n++;
 482         }
 483         if (bEnergy)
 484         {
 485             dhc->dh_energy=dhc->dh+n;
 486             mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
 487                              ir->fepvals->dh_hist_spacing, ndhmax,
 488                              dhbtEN, 0, 0, NULL);
 489             n++;
 490         }
 491         /* add the dhdl's */
 492         n_lambda_components=0;
 493         if (fep->dhdl_derivatives == edhdlderivativesYES)
 494         {
 495             dhc->dh_dhdl = dhc->dh + n;
 496             for (i=0;i<efptNR;i++)
 497             {
 498                 if (ir->fepvals->separate_dvdl[i])
 499                 {
 500                     /* we give it init_lambda for compatibility */
 501                     mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
 502                                      ir->fepvals->dh_hist_spacing, ndhmax,
 503                                      dhbtDHDL, n_lambda_components, 1,
 504                                      &(fep->init_lambda));
 505                     n++;
 506                     n_lambda_components++;
 507                 }
 508             }
 509         }
 510         else
 511         {
 512             for (i=0;i<efptNR;i++)
 513             {
 514                 if (ir->fepvals->separate_dvdl[i])
 515                 {
 516                     n_lambda_components++; /* count the components */
 517                 }
 518             }
 519
 520         }
 521         /* add the lambdas */
 522         dhc->dh_du = dhc->dh + n;
 523         snew(lambda_vec, n_lambda_components);
 524         for(i=ir->fepvals->lambda_start_n;i<ir->fepvals->lambda_stop_n;i++)
 525         {
 526             int k=0;
 527
 528             for(j=0;j<efptNR;j++)
 529             {
 530                 if (ir->fepvals->separate_dvdl[j])
 531                 {
 532                     lambda_vec[k++] = fep->all_lambda[j][i];
 533                 }
 534             }
 535
 536             mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
 537                              ir->fepvals->dh_hist_spacing, ndhmax,
 538                              dhbtDH, 0, n_lambda_components, lambda_vec);
 539             n++;
 540         }
 541         sfree(lambda_vec);
 542         if (bPV)
 543         {
 544             dhc->dh_pv=dhc->dh+n;
 545             mde_delta_h_init(dhc->dh+n, ir->fepvals->dh_hist_size,
 546                              ir->fepvals->dh_hist_spacing, ndhmax,
 547                              dhbtPV, 0, 0, NULL);
 548             n++;
 549         }
 550     }
 551 }
 552
 553 /* add a bunch of samples - note fep_state is double to allow for better data storage */
 554 void mde_delta_h_coll_add_dh(t_mde_delta_h_coll *dhc,
 555                              double fep_state,
 556                              double energy,
 557                              double pV,
 558                              double *dhdl,
 559                              double *foreign_dU,
 560                              double time)
 561 {
 562     int i;
 563
 564     if (!dhc->start_time_set)
 565     {
 566         dhc->start_time_set=TRUE;
 567         dhc->start_time=time;
 568     }
 569
 570     for (i=0;i<dhc->ndhdl;i++)
 571     {
 572         mde_delta_h_add_dh(dhc->dh_dhdl+i, dhdl[i], time);
 573     }
 574     for (i=0;i<dhc->nlambda;i++)
 575     {
 576         mde_delta_h_add_dh(dhc->dh_du+i, foreign_dU[i], time);
 577     }
 578     if (dhc->dh_pv != NULL)
 579     {
 580         mde_delta_h_add_dh(dhc->dh_pv, pV, time);
 581     }
 582     if (dhc->dh_energy != NULL)
 583     {
 584         mde_delta_h_add_dh(dhc->dh_energy,energy,time);
 585     }
 586     if (dhc->dh_expanded != NULL)
 587     {
 588         mde_delta_h_add_dh(dhc->dh_expanded,fep_state,time);
 589     }
 590
 591 }
 592
 593 /* write the metadata associated with all the du blocks, and call
 594    handle_block to write out all the du blocks */
 595 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll *dhc,
 596                                    t_enxframe *fr, int nblock)
 597 {
 598     int i;
 599     t_enxblock *blk;
 600
 601     /* add one block with one subblock as the collection's own data */
 602     nblock++;
 603     add_blocks_enxframe(fr, nblock);
 604     blk=fr->block + (nblock-1);
 605
 606     /* only allocate lambda vector component blocks if they must be written out
 607        for backward compatibility */
 608     if (dhc->native_lambda_components!=NULL)
 609     {
 610         add_subblocks_enxblock(blk, 2);
 611     }
 612     else
 613     {
 614         add_subblocks_enxblock(blk, 1);
 615     }
 616
 617     dhc->subblock_d[0] = dhc->temperature; /* temperature */
 618     dhc->subblock_d[1] = dhc->start_time; /* time of first sample */
 619     dhc->subblock_d[2] = dhc->delta_time; /* time difference between samples */
 620     dhc->subblock_d[3] = dhc->start_lambda; /* old-style lambda at starttime */
 621     dhc->subblock_d[4] = dhc->delta_lambda; /* lambda diff. between samples */
 622     /* set the lambda vector components if they exist */
 623     if (dhc->native_lambda_components!=NULL)
 624     {
 625         for(i=0;i<dhc->n_lambda_vec;i++)
 626         {
 627             dhc->subblock_d[5+i] = dhc->native_lambda_vec[i];
 628         }
 629     }
 630     blk->id=enxDHCOLL;
 631     blk->sub[0].nr=5 + dhc->n_lambda_vec;
 632     blk->sub[0].type=xdr_datatype_double;
 633     blk->sub[0].dval=dhc->subblock_d;
 634
 635     if (dhc->native_lambda_components != NULL)
 636     {
 637         dhc->subblock_i[0] = dhc->lambda_index;
 638         /* set the lambda vector component IDs if they exist */
 639         dhc->subblock_i[1] = dhc->n_lambda_vec;
 640         for(i=0;i<dhc->n_lambda_vec;i++)
 641         {
 642             dhc->subblock_i[i+2] = dhc->native_lambda_components[i];
 643         }
 644         blk->sub[1].nr=2 + dhc->n_lambda_vec;
 645         blk->sub[1].type=xdr_datatype_int;
 646         blk->sub[1].ival=dhc->subblock_i;
 647     }
 648
 649     for(i=0;i<dhc->ndh;i++)
 650     {
 651         nblock++;
 652         add_blocks_enxframe(fr, nblock);
 653         blk=fr->block + (nblock-1);
 654
 655         mde_delta_h_handle_block(dhc->dh+i, blk);
 656     }
 657 }
 658
 659 /* reset the data for a new round */
 660 void mde_delta_h_coll_reset(t_mde_delta_h_coll *dhc)
 661 {
 662     int i;
 663     for(i=0;i<dhc->ndh;i++)
 664     {
 665         if (dhc->dh[i].written)
 666         {
 667             /* we can now throw away the data */
 668             mde_delta_h_reset(dhc->dh + i);
 669         }
 670     }
 671     dhc->start_time_set=FALSE;
 672 }
 673
 674 /* set the energyhistory variables to save state */
 675 void mde_delta_h_coll_update_energyhistory(t_mde_delta_h_coll *dhc,
 676                                            energyhistory_t *enerhist)
 677 {
 678     int i;
 679     if (!enerhist->dht)
 680     {
 681         snew(enerhist->dht, 1);
 682         snew(enerhist->dht->ndh, dhc->ndh);
 683         snew(enerhist->dht->dh, dhc->ndh);
 684         enerhist->dht->nndh=dhc->ndh;
 685     }
 686     else
 687     {
 688         if (enerhist->dht->nndh != dhc->ndh)
 689             gmx_incons("energy history number of delta_h histograms != inputrec's number");
 690     }
 691     for(i=0;i<dhc->ndh;i++)
 692     {
 693         enerhist->dht->dh[i] = dhc->dh[i].dh;
 694         enerhist->dht->ndh[i] = dhc->dh[i].ndh;
 695     }
 696     enerhist->dht->start_time=dhc->start_time;
 697     enerhist->dht->start_lambda=dhc->start_lambda;
 698 }
 699
 700
 701
 702 /* restore the variables from an energyhistory */
 703 void mde_delta_h_coll_restore_energyhistory(t_mde_delta_h_coll *dhc,
 704                                             energyhistory_t *enerhist)
 705 {
 706     int i;
 707     unsigned int j;
 708
 709     if (dhc && !enerhist->dht)
 710         gmx_incons("No delta_h histograms in energy history");
 711     if (enerhist->dht->nndh != dhc->ndh)
 712         gmx_incons("energy history number of delta_h histograms != inputrec's number");
 713
 714     for(i=0;i<enerhist->dht->nndh;i++)
 715     {
 716         dhc->dh[i].ndh=enerhist->dht->ndh[i];
 717         for(j=0;j<dhc->dh[i].ndh;j++)
 718         {
 719             dhc->dh[i].dh[j] = enerhist->dht->dh[i][j];
 720         }
 721     }
 722     dhc->start_time=enerhist->dht->start_time;
 723     if (enerhist->dht->start_lambda_set)
 724         dhc->start_lambda=enerhist->dht->start_lambda;
 725     if (dhc->dh[0].ndh > 0)
 726         dhc->start_time_set=TRUE;
 727     else
 728         dhc->start_time_set=FALSE;
 729 }
 730
 731
 732
 733