src/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                              gmx_bool derivative, double foreign_lambda)
  63 {
  64     int i;
  65
  66     dh->derivative=derivative;
  67     dh->lambda=foreign_lambda;
  68
  69     dh->ndhmax=ndhmax+2;
  70     for(i=0;i<2;i++)
  71     {
  72         dh->bin[i]=NULL;
  73     }
  74
  75     snew(dh->dh, ndhmax);
  76     if ( nbins <= 0 || dx<GMX_REAL_EPS*10 )
  77     {
  78         dh->nhist=0;
  79     }
  80     else
  81     {
  82         int i;
  83         /* pre-allocate the histogram */
  84         if (derivative)
  85             dh->nhist=2;
  86         else
  87             dh->nhist=1;
  88         dh->dx=dx;
  89         dh->nbins=nbins;
  90         for(i=0;i<dh->nhist;i++)
  91         {
  92             snew(dh->bin[i], nbins);
  93         }
  94     }
  95     mde_delta_h_reset(dh);
  96 }
  97
  98 /* free the contents of the delta_h list */
  99 static void mde_delta_h_free(t_mde_delta_h *dh)
 100 {
 101     int i;
 102     for(i=0;i<dh->nhist;i++)
 103     {
 104         sfree(dh->bin[i]);
 105     }
 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, double time)
 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             min_dh=f*dh->dh[i];
 140         if (f*dh->dh[i] > max_dh)
 141             max_dh=f*dh->dh[i];
 142     }
 143
 144     /* reset the histogram */
 145     for(i=0;i<dh->nbins;i++)
 146     {
 147         dh->bin[hi][i]=0;
 148     }
 149     dh->maxbin[hi]=0;
 150
 151     /* The starting point of the histogram is the lowest value found:
 152        that value has the highest contribution to the free energy.
 153
 154        Get this start value in number of histogram dxs from zero,
 155        as an integer.*/
 156     dh->x0[hi] = (gmx_large_int_t)floor(min_dh/dx);
 157
 158     min_dh_hist=(dh->x0[hi])*dx;
 159     max_dh_hist=(dh->x0[hi] + dh->nbins + 1)*dx;
 160
 161     /* and fill the histogram*/
 162     for(i=0;i<dh->ndh;i++)
 163     {
 164         unsigned int bin;
 165
 166         /* Determine the bin number. If it doesn't fit into the histogram,
 167            add it to the last bin.
 168            We check the max_dh_int range because converting to integers
 169            might lead to overflow with unpredictable results.*/
 170         if ( (f*dh->dh[i] >= min_dh_hist) && (f*dh->dh[i] <= max_dh_hist ) )
 171         {
 172             bin = (unsigned int)( (f*dh->dh[i] - min_dh_hist)/dx );
 173         }
 174         else
 175         {
 176             bin = dh->nbins-1;
 177         }
 178         /* double-check here because of possible round-off errors*/
 179         if (bin >= dh->nbins)
 180         {
 181             bin = dh->nbins-1;
 182         }
 183         if (bin > dh->maxbin[hi])
 184         {
 185             dh->maxbin[hi] = bin;
 186         }
 187
 188         dh->bin[hi][bin]++;
 189     }
 190
 191     /* make sure we include a bin with 0 if we didn't use the full
 192        histogram width. This can then be used as an indication that
 193        all the data was binned. */
 194     if (dh->maxbin[hi] < dh->nbins-1)
 195         dh->maxbin[hi] += 1;
 196 }
 197
 198
 199 void mde_delta_h_handle_block(t_mde_delta_h *dh, t_enxblock *blk)
 200 {
 201     /* first check which type we should use: histogram or raw data */
 202     if (dh->nhist == 0)
 203     {
 204         /* We write raw data.
 205            Raw data consists of 3 subblocks: a block with the
 206            the foreign lambda, and the data itself */
 207         add_subblocks_enxblock(blk, 3);
 208
 209         blk->id=enxDH;
 210
 211         /* subblock 1 */
 212         dh->subblock_i[0]=dh->derivative ? 1 : 0; /* derivative type */
 213         blk->sub[0].nr=1;
 214         blk->sub[0].type=xdr_datatype_int;
 215         blk->sub[0].ival=dh->subblock_i;
 216
 217         /* subblock 2 */
 218         dh->subblock_d[0]=dh->lambda;
 219         blk->sub[1].nr=1;
 220         blk->sub[1].type=xdr_datatype_double;
 221         blk->sub[1].dval=dh->subblock_d;
 222
 223         /* subblock 3 */
 224         /* check if there's actual data to be written. */
 225         if (dh->ndh > 1)
 226         {
 227             blk->sub[2].nr=dh->ndh;
 228 #ifndef GMX_DOUBLE
 229             blk->sub[2].type=xdr_datatype_float;
 230             blk->sub[2].fval=dh->dh;
 231 #else
 232             blk->sub[2].type=xdr_datatype_double;
 233             blk->sub[2].dval=dh->dh;
 234 #endif
 235             dh->written=TRUE;
 236         }
 237         else
 238         {
 239             blk->sub[2].nr=0;
 240 #ifndef GMX_DOUBLE
 241             blk->sub[2].type=xdr_datatype_float;
 242 #else
 243             blk->sub[2].type=xdr_datatype_double;
 244 #endif
 245             blk->sub[2].dval=NULL;
 246         }
 247     }
 248     else
 249     {
 250         int nhist_written=0;
 251         int i;
 252
 253         /* check if there's actual data to be written. */
 254         if (dh->ndh > 1)
 255         {
 256             gmx_bool prev_complete=FALSE;
 257             /* Make the histogram(s) */
 258             for(i=0;i<dh->nhist;i++)
 259             {
 260                 if (!prev_complete)
 261                 {
 262                     /* the first histogram is always normal, and the
 263                        second one is always reverse */
 264                     mde_delta_h_make_hist(dh, i, i==1);
 265                     nhist_written++;
 266                     /* check whether this histogram contains all data: if the
 267                        last bin is 0, it does */
 268                     if (dh->bin[i][dh->nbins-1] == 0)
 269                         prev_complete=TRUE;
 270                     if (!dh->derivative)
 271                         prev_complete=TRUE;
 272                 }
 273             }
 274             dh->written=TRUE;
 275         }
 276
 277         /* A histogram consists of 2, 3 or 4 subblocks:
 278            the foreign lambda value + histogram spacing, the starting point,
 279            and the histogram data (0, 1 or 2 blocks). */
 280         add_subblocks_enxblock(blk, nhist_written+2);
 281         blk->id=enxDHHIST;
 282
 283         /* subblock 1: the foreign lambda value + the histogram spacing */
 284         dh->subblock_d[0]=dh->lambda;
 285         dh->subblock_d[1]=dh->dx;
 286         blk->sub[0].nr=2;
 287         blk->sub[0].type=xdr_datatype_double;
 288         blk->sub[0].dval=dh->subblock_d;
 289
 290         /* subblock 2: the starting point(s) as a long integer */
 291         dh->subblock_l[0]=nhist_written;
 292         dh->subblock_l[1]=dh->derivative ? 1 : 0;
 293         for(i=0;i<nhist_written;i++)
 294             dh->subblock_l[2+i]=dh->x0[i];
 295
 296         blk->sub[1].nr=nhist_written+2;
 297         blk->sub[1].type=xdr_datatype_large_int;
 298         blk->sub[1].lval=dh->subblock_l;
 299
 300         /* subblock 3 + 4 : the histogram data */
 301         for(i=0;i<nhist_written;i++)
 302         {
 303             blk->sub[i+2].nr=dh->maxbin[i]+1; /* it's +1 because size=index+1
 304                                                  in C */
 305             blk->sub[i+2].type=xdr_datatype_int;
 306             blk->sub[i+2].ival=dh->bin[i];
 307         }
 308     }
 309 }
 310
 311 /* initialize the collection*/
 312 void mde_delta_h_coll_init(t_mde_delta_h_coll *dhc, const t_inputrec *ir)
 313 {
 314     int i;
 315     int ndhmax=ir->nstenergy/ir->nstcalcenergy;
 316
 317     dhc->temp=ir->opts.ref_t[0];
 318     dhc->start_time=0.;
 319     dhc->start_lambda=ir->init_lambda;
 320
 321     dhc->delta_time=ir->delta_t*ir->nstdhdl;
 322     dhc->delta_lambda=ir->delta_lambda*ir->nstdhdl;
 323
 324     dhc->start_time_set=FALSE;
 325
 326     if (ir->dhdl_derivatives == dhdlderivativesYES)
 327     {
 328         dhc->ndhdl=1;
 329     }
 330     else
 331     {
 332         dhc->ndhdl=0;
 333     }
 334
 335     dhc->ndh=ir->n_flambda+dhc->ndhdl;
 336     snew(dhc->dh, dhc->ndh);
 337     for(i=0;i<dhc->ndh;i++)
 338     {
 339         if (i<dhc->ndhdl)
 340         {
 341             mde_delta_h_init(dhc->dh + i, ir->dh_hist_size,
 342                              ir->dh_hist_spacing, ndhmax,
 343                              TRUE, dhc->start_lambda);
 344         }
 345         else
 346         {
 347             mde_delta_h_init(dhc->dh + i, ir->dh_hist_size,
 348                              ir->dh_hist_spacing, ndhmax,
 349                              FALSE,
 350                              ir->flambda[i-dhc->ndhdl] );
 351         }
 352     }
 353 }
 354
 355 /* add a bunch of samples */
 356 void mde_delta_h_coll_add_dh(t_mde_delta_h_coll *dhc,
 357                              double dhdl,
 358                              double *U, double time,
 359                              double native_lambda)
 360 {
 361     int i;
 362
 363     if (!dhc->start_time_set)
 364     {
 365         dhc->start_time_set=TRUE;
 366         dhc->start_time=time;
 367         dhc->start_lambda=native_lambda;
 368     }
 369     for(i=0;i<dhc->ndh;i++)
 370     {
 371         if (i<dhc->ndhdl)
 372         {
 373             mde_delta_h_add_dh(dhc->dh + i, dhdl, time);
 374         }
 375         else
 376         {
 377             mde_delta_h_add_dh(dhc->dh + i, U[i+1-dhc->ndhdl] - U[0], time);
 378         }
 379     }
 380 }
 381
 382 /* write the data associated with all the du blocks, but not the blocks
 383    themselves. */
 384 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll *dhc,
 385                                    t_enxframe *fr, int nblock)
 386 {
 387     int i;
 388     t_enxblock *blk;
 389
 390     /* add one block with one subblock as the collection's own data */
 391     nblock++;
 392     add_blocks_enxframe(fr, nblock);
 393     blk=fr->block + (nblock-1);
 394
 395     add_subblocks_enxblock(blk, 1);
 396
 397     dhc->subblock_d[0] = dhc->temp; /* temperature */
 398     dhc->subblock_d[1] = dhc->start_time; /* time of first sample */
 399     dhc->subblock_d[2] = dhc->delta_time; /* time difference between samples */
 400     dhc->subblock_d[3] = dhc->start_lambda; /* lambda at starttime */
 401     dhc->subblock_d[4] = dhc->delta_lambda; /* lambda diff. between samples */
 402
 403     blk->id=enxDHCOLL;
 404     blk->sub[0].nr=5;
 405     blk->sub[0].type=xdr_datatype_double;
 406     blk->sub[0].dval=dhc->subblock_d;
 407
 408     for(i=0;i<dhc->ndh;i++)
 409     {
 410         nblock++;
 411         add_blocks_enxframe(fr, nblock);
 412         blk=fr->block + (nblock-1);
 413
 414         mde_delta_h_handle_block(dhc->dh+i, blk);
 415     }
 416 }
 417
 418 /* reset the data for a new round */
 419 void mde_delta_h_coll_reset(t_mde_delta_h_coll *dhc)
 420 {
 421     int i;
 422     for(i=0;i<dhc->ndh;i++)
 423     {
 424         if (dhc->dh[i].written)
 425         {
 426             /* we can now throw away the data */
 427             mde_delta_h_reset(dhc->dh + i);
 428         }
 429     }
 430     dhc->start_time_set=FALSE;
 431 }
 432
 433 /* set the energyhistory variables to save state */
 434 void mde_delta_h_coll_update_energyhistory(t_mde_delta_h_coll *dhc,
 435                                            energyhistory_t *enerhist)
 436 {
 437     int i;
 438     if (!enerhist->dht)
 439     {
 440         snew(enerhist->dht, 1);
 441         snew(enerhist->dht->ndh, dhc->ndh);
 442         snew(enerhist->dht->dh, dhc->ndh);
 443         enerhist->dht->nndh=dhc->ndh;
 444     }
 445     else
 446     {
 447         if (enerhist->dht->nndh != dhc->ndh)
 448             gmx_incons("energy history number of delta_h histograms != inputrec's number");
 449     }
 450     for(i=0;i<dhc->ndh;i++)
 451     {
 452         enerhist->dht->dh[i] = dhc->dh[i].dh;
 453         enerhist->dht->ndh[i] = dhc->dh[i].ndh;
 454     }
 455     enerhist->dht->start_time=dhc->start_time;
 456     enerhist->dht->start_lambda=dhc->start_lambda;
 457 }
 458
 459
 460
 461 /* restore the variables from an energyhistory */
 462 void mde_delta_h_coll_restore_energyhistory(t_mde_delta_h_coll *dhc,
 463                                             energyhistory_t *enerhist)
 464 {
 465     int i;
 466     unsigned int j;
 467
 468     if (dhc && !enerhist->dht)
 469         gmx_incons("No delta_h histograms in energy history");
 470     if (enerhist->dht->nndh != dhc->ndh)
 471         gmx_incons("energy history number of delta_h histograms != inputrec's number");
 472
 473     for(i=0;i<enerhist->dht->nndh;i++)
 474     {
 475         dhc->dh[i].ndh=enerhist->dht->ndh[i];
 476         for(j=0;j<dhc->dh[i].ndh;j++)
 477         {
 478             dhc->dh[i].dh[j] = enerhist->dht->dh[i][j];
 479         }
 480     }
 481     dhc->start_time=enerhist->dht->start_time;
 482     if (enerhist->dht->start_lambda_set)
 483         dhc->start_lambda=enerhist->dht->start_lambda;
 484     if (dhc->dh[0].ndh > 0)
 485         dhc->start_time_set=TRUE;
 486     else
 487         dhc->start_time_set=FALSE;
 488 }
 489
 490
 491
 492