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