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