1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
45 #include "gmx_fatal.h"
50 #include "mdebin_bar.h"
52 /* reset the delta_h list to prepare it for new values */
53 static void mde_delta_h_reset(t_mde_delta_h *dh)
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)
71 snew(dh->dh, dh->ndhmax);
72 snew(dh->dhf,dh->ndhmax);
74 if ( nbins <= 0 || dx<GMX_REAL_EPS*10 )
81 /* pre-allocate the histogram */
82 dh->nhist=2; /* energies and derivatives histogram */
85 for(i=0;i<dh->nhist;i++)
87 snew(dh->bin[i], dh->nbins);
90 mde_delta_h_reset(dh);
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)
96 if (dh->ndh >= dh->ndhmax)
98 gmx_incons("delta_h array not big enough!");
100 dh->dh[dh->ndh]=delta_h;
104 /* construct histogram with index hi */
105 static void mde_delta_h_make_hist(t_mde_delta_h *dh, int hi, gmx_bool invert)
107 double min_dh = FLT_MAX;
108 double max_dh = -FLT_MAX;
110 double max_dh_hist; /* maximum binnable dh value */
111 double min_dh_hist; /* minimum binnable dh value */
113 double f; /* energy mult. factor */
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 */
120 /* first find min and max */
121 for(i=0;i<dh->ndh;i++)
123 if (f*dh->dh[i] < min_dh)
125 if (f*dh->dh[i] > max_dh)
129 /* reset the histogram */
130 for(i=0;i<dh->nbins;i++)
136 /* The starting point of the histogram is the lowest value found:
137 that value has the highest contribution to the free energy.
139 Get this start value in number of histogram dxs from zero,
142 dh->x0[hi] = (gmx_large_int_t)floor(min_dh/dx);
144 min_dh_hist=(dh->x0[hi])*dx;
145 max_dh_hist=(dh->x0[hi] + dh->nbins + 1)*dx;
147 /* and fill the histogram*/
148 for(i=0;i<dh->ndh;i++)
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 ) )
158 bin = (unsigned int)( (f*dh->dh[i] - min_dh_hist)/dx );
164 /* double-check here because of possible round-off errors*/
165 if (bin >= dh->nbins)
169 if (bin > dh->maxbin[hi])
171 dh->maxbin[hi] = bin;
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)
185 void mde_delta_h_handle_block(t_mde_delta_h *dh, t_enxblock *blk)
187 /* first check which type we should use: histogram or raw data */
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);
200 dh->subblock_i[0]=dh->derivative ? 1 : 0; /* derivative type */
202 blk->sub[0].type=xdr_datatype_int;
203 blk->sub[0].ival=dh->subblock_i;
206 dh->subblock_d[0]=dh->lambda;
208 blk->sub[1].type=xdr_datatype_double;
209 blk->sub[1].dval=dh->subblock_d;
212 /* check if there's actual data to be written. */
216 blk->sub[2].nr=dh->ndh;
217 /* For F@H for now. */
220 blk->sub[2].type=xdr_datatype_float;
221 for(i=0;i<dh->ndh;i++)
223 dh->dhf[i] = (float)dh->dh[i];
225 blk->sub[2].fval=dh->dhf;
227 blk->sub[2].type=xdr_datatype_double;
228 blk->sub[2].dval=dh->dh;
236 blk->sub[2].type=xdr_datatype_float;
237 blk->sub[2].fval=NULL;
239 blk->sub[2].type=xdr_datatype_double;
241 blk->sub[2].dval=NULL;
249 /* check if there's actual data to be written. */
252 gmx_bool prev_complete=FALSE;
253 /* Make the histogram(s) */
254 for(i=0;i<dh->nhist;i++)
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);
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)
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);
279 /* subblock 1: the foreign lambda value + the histogram spacing */
280 dh->subblock_d[0]=dh->lambda;
281 dh->subblock_d[1]=dh->dx;
283 blk->sub[0].type=xdr_datatype_double;
284 blk->sub[0].dval=dh->subblock_d;
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];
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;
296 /* subblock 3 + 4 : the histogram data */
297 for(i=0;i<nhist_written;i++)
299 blk->sub[i+2].nr=dh->maxbin[i]+1; /* it's +1 because size=index+1
301 blk->sub[i+2].type=xdr_datatype_int;
302 blk->sub[i+2].ival=dh->bin[i];
307 /* initialize the collection*/
308 void mde_delta_h_coll_init(t_mde_delta_h_coll *dhc, const t_inputrec *ir)
312 int ndhmax=ir->nstenergy/ir->nstcalcenergy;
314 dhc->temperature=ir->opts.ref_t[0]; /* only store system temperature */
316 dhc->delta_time=ir->delta_t*ir->fepvals->nstdhdl;
317 dhc->start_time_set=FALSE;
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;
323 /* total number of raw data points in the sample */
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) {
331 /* whether to print energies */
332 if (ir->fepvals->bPrintEnergy) {
337 for (i=0;i<efptNR;i++)
339 if (ir->fepvals->separate_dvdl[i])
345 /* add the lambdas */
346 dhc->ndh += ir->fepvals->n_lambda;
348 if (ir->epc > epcNO) {
349 dhc->ndh += 1; /* include pressure-volume work */
352 snew(dhc->dh, dhc->ndh);
353 for(i=0;i<dhc->ndh;i++)
355 mde_delta_h_init(dhc->dh+i, ir->fepvals->dh_hist_size,
356 ir->fepvals->dh_hist_spacing, ndhmax);
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,
376 if (!dhc->start_time_set)
378 dhc->start_time_set=TRUE;
379 dhc->start_time=time;
385 mde_delta_h_add_dh(dhc->dh+n,fep_state,time);
390 mde_delta_h_add_dh(dhc->dh+n,energy,time);
393 for (i=0;i<ndhdl;i++)
395 mde_delta_h_add_dh(dhc->dh+n, dhdl[i], time);
398 for (i=0;i<nlambda;i++)
400 mde_delta_h_add_dh(dhc->dh+n, foreign_dU[i], time);
405 mde_delta_h_add_dh(dhc->dh+n, pV, time);
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?*/
413 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll *dhc,
414 t_enxframe *fr, int nblock)
419 /* add one block with one subblock as the collection's own data */
421 add_blocks_enxframe(fr, nblock);
422 blk=fr->block + (nblock-1);
424 add_subblocks_enxblock(blk, 1);
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 */
433 blk->sub[0].type=xdr_datatype_double;
434 blk->sub[0].dval=dhc->subblock_d;
436 for(i=0;i<dhc->ndh;i++)
439 add_blocks_enxframe(fr, nblock);
440 blk=fr->block + (nblock-1);
442 mde_delta_h_handle_block(dhc->dh+i, blk);
446 /* reset the data for a new round */
447 void mde_delta_h_coll_reset(t_mde_delta_h_coll *dhc)
450 for(i=0;i<dhc->ndh;i++)
452 if (dhc->dh[i].written)
454 /* we can now throw away the data */
455 mde_delta_h_reset(dhc->dh + i);
458 dhc->start_time_set=FALSE;
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)
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;
475 if (enerhist->dht->nndh != dhc->ndh)
476 gmx_incons("energy history number of delta_h histograms != inputrec's number");
478 for(i=0;i<dhc->ndh;i++)
480 enerhist->dht->dh[i] = dhc->dh[i].dh;
481 enerhist->dht->ndh[i] = dhc->dh[i].ndh;
483 enerhist->dht->start_time=dhc->start_time;
484 enerhist->dht->start_lambda=dhc->start_lambda;
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)
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");
501 for(i=0;i<enerhist->dht->nndh;i++)
503 dhc->dh[i].ndh=enerhist->dht->ndh[i];
504 for(j=0;j<dhc->dh[i].ndh;j++)
506 dhc->dh[i].dh[j] = enerhist->dht->dh[i][j];
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;
515 dhc->start_time_set=FALSE;