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