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,
62 gmx_bool derivative, double foreign_lambda)
66 dh->derivative=derivative;
67 dh->lambda=foreign_lambda;
75 snew(dh->dh, dh->ndhmax);
76 if ( nbins <= 0 || dx<GMX_REAL_EPS*10 )
83 /* pre-allocate the histogram */
90 for(i=0;i<dh->nhist;i++)
92 snew(dh->bin[i], dh->nbins);
95 mde_delta_h_reset(dh);
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)
101 if (dh->ndh >= dh->ndhmax)
103 gmx_incons("delta_h array not big enough!");
105 dh->dh[dh->ndh]=delta_h;
109 /* construct histogram with index hi */
110 static void mde_delta_h_make_hist(t_mde_delta_h *dh, int hi, gmx_bool invert)
112 double min_dh = FLT_MAX;
113 double max_dh = -FLT_MAX;
115 double max_dh_hist; /* maximum binnable dh value */
116 double min_dh_hist; /* minimum binnable dh value */
118 double f; /* energy mult. factor */
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 */
125 /* first find min and max */
126 for(i=0;i<dh->ndh;i++)
128 if (f*dh->dh[i] < min_dh)
130 if (f*dh->dh[i] > max_dh)
134 /* reset the histogram */
135 for(i=0;i<dh->nbins;i++)
141 /* The starting point of the histogram is the lowest value found:
142 that value has the highest contribution to the free energy.
144 Get this start value in number of histogram dxs from zero,
146 dh->x0[hi] = (gmx_large_int_t)floor(min_dh/dx);
148 min_dh_hist=(dh->x0[hi])*dx;
149 max_dh_hist=(dh->x0[hi] + dh->nbins + 1)*dx;
151 /* and fill the histogram*/
152 for(i=0;i<dh->ndh;i++)
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 ) )
162 bin = (unsigned int)( (f*dh->dh[i] - min_dh_hist)/dx );
168 /* double-check here because of possible round-off errors*/
169 if (bin >= dh->nbins)
173 if (bin > dh->maxbin[hi])
175 dh->maxbin[hi] = bin;
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)
189 void mde_delta_h_handle_block(t_mde_delta_h *dh, t_enxblock *blk)
191 /* first check which type we should use: histogram or raw data */
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);
202 dh->subblock_i[0]=dh->derivative ? 1 : 0; /* derivative type */
204 blk->sub[0].type=xdr_datatype_int;
205 blk->sub[0].ival=dh->subblock_i;
208 dh->subblock_d[0]=dh->lambda;
210 blk->sub[1].type=xdr_datatype_double;
211 blk->sub[1].dval=dh->subblock_d;
214 /* check if there's actual data to be written. */
217 blk->sub[2].nr=dh->ndh;
219 blk->sub[2].type=xdr_datatype_float;
220 blk->sub[2].fval=dh->dh;
222 blk->sub[2].type=xdr_datatype_double;
223 blk->sub[2].dval=dh->dh;
231 blk->sub[2].type=xdr_datatype_float;
233 blk->sub[2].type=xdr_datatype_double;
235 blk->sub[2].dval=NULL;
243 /* check if there's actual data to be written. */
246 gmx_bool prev_complete=FALSE;
247 /* Make the histogram(s) */
248 for(i=0;i<dh->nhist;i++)
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);
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)
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);
273 /* subblock 1: the foreign lambda value + the histogram spacing */
274 dh->subblock_d[0]=dh->lambda;
275 dh->subblock_d[1]=dh->dx;
277 blk->sub[0].type=xdr_datatype_double;
278 blk->sub[0].dval=dh->subblock_d;
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];
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;
290 /* subblock 3 + 4 : the histogram data */
291 for(i=0;i<nhist_written;i++)
293 blk->sub[i+2].nr=dh->maxbin[i]+1; /* it's +1 because size=index+1
295 blk->sub[i+2].type=xdr_datatype_int;
296 blk->sub[i+2].ival=dh->bin[i];
301 /* initialize the collection*/
302 void mde_delta_h_coll_init(t_mde_delta_h_coll *dhc, const t_inputrec *ir)
305 int ndhmax=ir->nstenergy/ir->nstcalcenergy;
307 dhc->temp=ir->opts.ref_t[0];
309 dhc->start_lambda=ir->init_lambda;
311 dhc->delta_time=ir->delta_t*ir->nstdhdl;
312 dhc->delta_lambda=ir->delta_lambda*ir->nstdhdl;
314 dhc->start_time_set=FALSE;
316 if (ir->dhdl_derivatives == dhdlderivativesYES)
325 dhc->ndh=ir->n_flambda+dhc->ndhdl;
326 snew(dhc->dh, dhc->ndh);
327 for(i=0;i<dhc->ndh;i++)
331 mde_delta_h_init(dhc->dh + i, ir->dh_hist_size,
332 ir->dh_hist_spacing, ndhmax,
333 TRUE, dhc->start_lambda);
337 mde_delta_h_init(dhc->dh + i, ir->dh_hist_size,
338 ir->dh_hist_spacing, ndhmax,
340 ir->flambda[i-dhc->ndhdl] );
345 /* add a bunch of samples */
346 void mde_delta_h_coll_add_dh(t_mde_delta_h_coll *dhc,
348 double *U, double time,
349 double native_lambda)
353 if (!dhc->start_time_set)
355 dhc->start_time_set=TRUE;
356 dhc->start_time=time;
357 dhc->start_lambda=native_lambda;
359 for(i=0;i<dhc->ndh;i++)
363 mde_delta_h_add_dh(dhc->dh + i, dhdl, time);
367 mde_delta_h_add_dh(dhc->dh + i, U[i+1-dhc->ndhdl] - U[0], time);
372 /* write the data associated with all the du blocks, but not the blocks
374 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll *dhc,
375 t_enxframe *fr, int nblock)
380 /* add one block with one subblock as the collection's own data */
382 add_blocks_enxframe(fr, nblock);
383 blk=fr->block + (nblock-1);
385 add_subblocks_enxblock(blk, 1);
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 */
395 blk->sub[0].type=xdr_datatype_double;
396 blk->sub[0].dval=dhc->subblock_d;
398 for(i=0;i<dhc->ndh;i++)
401 add_blocks_enxframe(fr, nblock);
402 blk=fr->block + (nblock-1);
404 mde_delta_h_handle_block(dhc->dh+i, blk);
408 /* reset the data for a new round */
409 void mde_delta_h_coll_reset(t_mde_delta_h_coll *dhc)
412 for(i=0;i<dhc->ndh;i++)
414 if (dhc->dh[i].written)
416 /* we can now throw away the data */
417 mde_delta_h_reset(dhc->dh + i);
420 dhc->start_time_set=FALSE;
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)
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;
437 if (enerhist->dht->nndh != dhc->ndh)
438 gmx_incons("energy history number of delta_h histograms != inputrec's number");
440 for(i=0;i<dhc->ndh;i++)
442 enerhist->dht->dh[i] = dhc->dh[i].dh;
443 enerhist->dht->ndh[i] = dhc->dh[i].ndh;
445 enerhist->dht->start_time=dhc->start_time;
446 enerhist->dht->start_lambda=dhc->start_lambda;
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)
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");
463 for(i=0;i<enerhist->dht->nndh;i++)
465 dhc->dh[i].ndh=enerhist->dht->ndh[i];
466 for(j=0;j<dhc->dh[i].ndh;j++)
468 dhc->dh[i].dh[j] = enerhist->dht->dh[i][j];
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;
477 dhc->start_time_set=FALSE;