2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "mdebin_bar.h"
45 #include "gromacs/fileio/enxio.h"
46 #include "gromacs/mdlib/mdebin.h"
47 #include "gromacs/mdtypes/energyhistory.h"
48 #include "gromacs/utility/fatalerror.h"
49 #include "gromacs/utility/smalloc.h"
51 /* reset the delta_h list to prepare it for new values */
52 static void mde_delta_h_reset(t_mde_delta_h
*dh
)
58 /* initialize the delta_h list */
59 static void mde_delta_h_init(t_mde_delta_h
*dh
, int nbins
,
60 double dx
, unsigned int ndhmax
,
61 int type
, int derivative
, int nlambda
,
67 dh
->derivative
= derivative
;
69 dh
->nlambda
= nlambda
;
71 snew(dh
->lambda
, nlambda
);
72 for (i
= 0; i
< nlambda
; i
++)
74 dh
->lambda
[i
] = lambda
[i
];
78 snew(dh
->subblock_meta_d
, dh
->nlambda
+1);
80 dh
->ndhmax
= ndhmax
+2;
81 for (i
= 0; i
< 2; i
++)
86 snew(dh
->dh
, dh
->ndhmax
);
87 snew(dh
->dhf
, dh
->ndhmax
);
89 if (nbins
<= 0 || dx
< GMX_REAL_EPS
*10)
96 /* pre-allocate the histogram */
97 dh
->nhist
= 2; /* energies and derivatives histogram */
100 for (i
= 0; i
< dh
->nhist
; i
++)
102 snew(dh
->bin
[i
], dh
->nbins
);
105 mde_delta_h_reset(dh
);
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
)
111 if (dh
->ndh
>= dh
->ndhmax
)
113 gmx_incons("delta_h array not big enough!");
115 dh
->dh
[dh
->ndh
] = delta_h
;
119 /* construct histogram with index hi */
120 static void mde_delta_h_make_hist(t_mde_delta_h
*dh
, int hi
, gmx_bool invert
)
122 double min_dh
= FLT_MAX
;
123 double max_dh
= -FLT_MAX
;
125 double max_dh_hist
; /* maximum binnable dh value */
126 double min_dh_hist
; /* minimum binnable dh value */
128 double f
; /* energy mult. factor */
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 */
135 /* first find min and max */
136 for (i
= 0; i
< dh
->ndh
; i
++)
138 if (f
*dh
->dh
[i
] < min_dh
)
140 min_dh
= f
*dh
->dh
[i
];
142 if (f
*dh
->dh
[i
] > max_dh
)
144 max_dh
= f
*dh
->dh
[i
];
148 /* reset the histogram */
149 for (i
= 0; i
< dh
->nbins
; i
++)
155 /* The starting point of the histogram is the lowest value found:
156 that value has the highest contribution to the free energy.
158 Get this start value in number of histogram dxs from zero,
161 dh
->x0
[hi
] = (gmx_int64_t
)floor(min_dh
/dx
);
163 min_dh_hist
= (dh
->x0
[hi
])*dx
;
164 max_dh_hist
= (dh
->x0
[hi
] + dh
->nbins
+ 1)*dx
;
166 /* and fill the histogram*/
167 for (i
= 0; i
< dh
->ndh
; i
++)
171 /* Determine the bin number. If it doesn't fit into the histogram,
172 add it to the last bin.
173 We check the max_dh_int range because converting to integers
174 might lead to overflow with unpredictable results.*/
175 if ( (f
*dh
->dh
[i
] >= min_dh_hist
) && (f
*dh
->dh
[i
] <= max_dh_hist
) )
177 bin
= (unsigned int)( (f
*dh
->dh
[i
] - min_dh_hist
)/dx
);
183 /* double-check here because of possible round-off errors*/
184 if (bin
>= dh
->nbins
)
188 if (bin
> dh
->maxbin
[hi
])
190 dh
->maxbin
[hi
] = bin
;
196 /* make sure we include a bin with 0 if we didn't use the full
197 histogram width. This can then be used as an indication that
198 all the data was binned. */
199 if (dh
->maxbin
[hi
] < dh
->nbins
-1)
206 void mde_delta_h_handle_block(t_mde_delta_h
*dh
, t_enxblock
*blk
)
208 /* first check which type we should use: histogram or raw data */
213 /* We write raw data.
214 Raw data consists of 3 subblocks: an int metadata block
215 with type and derivative index, a foreign lambda block
216 and and the data itself */
217 add_subblocks_enxblock(blk
, 3);
222 dh
->subblock_meta_i
[0] = dh
->type
; /* block data type */
223 dh
->subblock_meta_i
[1] = dh
->derivative
; /* derivative direction if
224 applicable (in indices
225 starting from first coord in
226 the main delta_h_coll) */
228 blk
->sub
[0].type
= xdr_datatype_int
;
229 blk
->sub
[0].ival
= dh
->subblock_meta_i
;
232 for (i
= 0; i
< dh
->nlambda
; i
++)
234 dh
->subblock_meta_d
[i
] = dh
->lambda
[i
];
236 blk
->sub
[1].nr
= dh
->nlambda
;
237 blk
->sub
[1].type
= xdr_datatype_double
;
238 blk
->sub
[1].dval
= dh
->subblock_meta_d
;
241 /* check if there's actual data to be written. */
247 blk
->sub
[2].nr
= dh
->ndh
;
248 /* For F@H for now. */
251 blk
->sub
[2].type
= xdr_datatype_float
;
252 for (i
= 0; i
< dh
->ndh
; i
++)
254 dh
->dhf
[i
] = (float)dh
->dh
[i
];
256 blk
->sub
[2].fval
= dh
->dhf
;
258 blk
->sub
[2].type
= xdr_datatype_double
;
259 blk
->sub
[2].dval
= dh
->dh
;
267 blk
->sub
[2].type
= xdr_datatype_float
;
268 blk
->sub
[2].fval
= NULL
;
270 blk
->sub
[2].type
= xdr_datatype_double
;
271 blk
->sub
[2].dval
= NULL
;
277 int nhist_written
= 0;
281 /* TODO histogram metadata */
282 /* check if there's actual data to be written. */
285 gmx_bool prev_complete
= FALSE
;
286 /* Make the histogram(s) */
287 for (i
= 0; i
< dh
->nhist
; i
++)
291 /* the first histogram is always normal, and the
292 second one is always reverse */
293 mde_delta_h_make_hist(dh
, i
, i
== 1);
295 /* check whether this histogram contains all data: if the
296 last bin is 0, it does */
297 if (dh
->bin
[i
][dh
->nbins
-1] == 0)
299 prev_complete
= TRUE
;
303 prev_complete
= TRUE
;
310 /* A histogram consists of 2, 3 or 4 subblocks:
311 the foreign lambda value + histogram spacing, the starting point,
312 and the histogram data (0, 1 or 2 blocks). */
313 add_subblocks_enxblock(blk
, nhist_written
+2);
316 /* subblock 1: the lambda value + the histogram spacing */
317 if (dh
->nlambda
== 1)
319 /* for backward compatibility */
320 dh
->subblock_meta_d
[0] = dh
->lambda
[0];
324 dh
->subblock_meta_d
[0] = -1;
325 for (i
= 0; i
< dh
->nlambda
; i
++)
327 dh
->subblock_meta_d
[2+i
] = dh
->lambda
[i
];
330 dh
->subblock_meta_d
[1] = dh
->dx
;
331 blk
->sub
[0].nr
= 2+ ((dh
->nlambda
> 1) ? dh
->nlambda
: 0);
332 blk
->sub
[0].type
= xdr_datatype_double
;
333 blk
->sub
[0].dval
= dh
->subblock_meta_d
;
335 /* subblock 2: the starting point(s) as a long integer */
336 dh
->subblock_meta_l
[0] = nhist_written
;
337 dh
->subblock_meta_l
[1] = dh
->type
; /*dh->derivative ? 1 : 0;*/
339 for (i
= 0; i
< nhist_written
; i
++)
341 dh
->subblock_meta_l
[k
++] = dh
->x0
[i
];
343 /* append the derivative data */
344 dh
->subblock_meta_l
[k
++] = dh
->derivative
;
346 blk
->sub
[1].nr
= nhist_written
+3;
347 blk
->sub
[1].type
= xdr_datatype_int64
;
348 blk
->sub
[1].lval
= dh
->subblock_meta_l
;
350 /* subblock 3 + 4 : the histogram data */
351 for (i
= 0; i
< nhist_written
; i
++)
353 blk
->sub
[i
+2].nr
= dh
->maxbin
[i
]+1; /* it's +1 because size=index+1
355 blk
->sub
[i
+2].type
= xdr_datatype_int
;
356 blk
->sub
[i
+2].ival
= dh
->bin
[i
];
361 /* initialize the collection*/
362 void mde_delta_h_coll_init(t_mde_delta_h_coll
*dhc
, const t_inputrec
*ir
)
366 int ndhmax
= ir
->nstenergy
/ir
->nstcalcenergy
;
367 t_lambda
*fep
= ir
->fepvals
;
369 dhc
->temperature
= ir
->opts
.ref_t
[0]; /* only store system temperature */
370 dhc
->start_time
= 0.;
371 dhc
->delta_time
= ir
->delta_t
*ir
->fepvals
->nstdhdl
;
372 dhc
->start_time_set
= FALSE
;
374 /* this is the compatibility lambda value. If it is >=0, it is valid,
375 and there is either an old-style lambda or a slow growth simulation. */
376 dhc
->start_lambda
= ir
->fepvals
->init_lambda
;
377 /* for continuous change of lambda values */
378 dhc
->delta_lambda
= ir
->fepvals
->delta_lambda
*ir
->fepvals
->nstdhdl
;
380 if (dhc
->start_lambda
< 0)
382 /* create the native lambda vectors */
383 dhc
->lambda_index
= fep
->init_fep_state
;
384 dhc
->n_lambda_vec
= 0;
385 for (i
= 0; i
< efptNR
; i
++)
387 if (fep
->separate_dvdl
[i
])
392 snew(dhc
->native_lambda_vec
, dhc
->n_lambda_vec
);
393 snew(dhc
->native_lambda_components
, dhc
->n_lambda_vec
);
395 for (i
= 0; i
< efptNR
; i
++)
397 if (fep
->separate_dvdl
[i
])
399 dhc
->native_lambda_components
[j
] = i
;
400 if (fep
->init_fep_state
>= 0 &&
401 fep
->init_fep_state
< fep
->n_lambda
)
403 dhc
->native_lambda_vec
[j
] =
404 fep
->all_lambda
[i
][fep
->init_fep_state
];
408 dhc
->native_lambda_vec
[j
] = -1;
416 /* don't allocate the meta-data subblocks for lambda vectors */
417 dhc
->native_lambda_vec
= NULL
;
418 dhc
->n_lambda_vec
= 0;
419 dhc
->native_lambda_components
= 0;
420 dhc
->lambda_index
= -1;
422 /* allocate metadata subblocks */
423 snew(dhc
->subblock_d
, 5 + dhc
->n_lambda_vec
);
424 snew(dhc
->subblock_i
, 1 + dhc
->n_lambda_vec
);
426 /* now decide which data to write out */
429 dhc
->dh_expanded
= NULL
;
430 dhc
->dh_energy
= NULL
;
433 /* total number of raw data point collections in the sample */
437 gmx_bool bExpanded
= FALSE
;
438 gmx_bool bEnergy
= FALSE
;
439 gmx_bool bPV
= FALSE
;
440 int n_lambda_components
= 0;
442 /* first count the number of states */
445 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
447 for (i
= 0; i
< efptNR
; i
++)
449 if (ir
->fepvals
->separate_dvdl
[i
])
456 /* add the lambdas */
457 dhc
->nlambda
= ir
->fepvals
->lambda_stop_n
- ir
->fepvals
->lambda_start_n
;
458 dhc
->ndh
+= dhc
->nlambda
;
459 /* another compatibility check */
460 if (dhc
->start_lambda
< 0)
462 /* include one more for the specification of the state, by lambda or
464 if (ir
->expandedvals
->elmcmove
> elmcmoveNO
)
469 /* whether to print energies */
470 if (ir
->fepvals
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
477 dhc
->ndh
+= 1; /* include pressure-volume work */
482 snew(dhc
->dh
, dhc
->ndh
);
484 /* now initialize them */
485 /* the order, for now, must match that of the dhdl.xvg file because of
486 how g_energy -odh is implemented */
490 dhc
->dh_expanded
= dhc
->dh
+n
;
491 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
492 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
493 dhbtEXPANDED
, 0, 0, NULL
);
498 dhc
->dh_energy
= dhc
->dh
+n
;
499 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
500 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
505 n_lambda_components
= 0;
506 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
508 dhc
->dh_dhdl
= dhc
->dh
+ n
;
509 for (i
= 0; i
< efptNR
; i
++)
511 if (ir
->fepvals
->separate_dvdl
[i
])
513 /* we give it init_lambda for compatibility */
514 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
515 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
516 dhbtDHDL
, n_lambda_components
, 1,
517 &(fep
->init_lambda
));
519 n_lambda_components
++;
525 for (i
= 0; i
< efptNR
; i
++)
527 if (ir
->fepvals
->separate_dvdl
[i
])
529 n_lambda_components
++; /* count the components */
534 /* add the lambdas */
535 dhc
->dh_du
= dhc
->dh
+ n
;
536 snew(lambda_vec
, n_lambda_components
);
537 for (i
= ir
->fepvals
->lambda_start_n
; i
< ir
->fepvals
->lambda_stop_n
; i
++)
541 for (j
= 0; j
< efptNR
; j
++)
543 if (ir
->fepvals
->separate_dvdl
[j
])
545 lambda_vec
[k
++] = fep
->all_lambda
[j
][i
];
549 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
550 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
551 dhbtDH
, 0, n_lambda_components
, lambda_vec
);
557 dhc
->dh_pv
= dhc
->dh
+n
;
558 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
559 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
566 /* add a bunch of samples - note fep_state is double to allow for better data storage */
567 void mde_delta_h_coll_add_dh(t_mde_delta_h_coll
*dhc
,
577 if (!dhc
->start_time_set
)
579 dhc
->start_time_set
= TRUE
;
580 dhc
->start_time
= time
;
583 for (i
= 0; i
< dhc
->ndhdl
; i
++)
585 mde_delta_h_add_dh(dhc
->dh_dhdl
+i
, dhdl
[i
]);
587 for (i
= 0; i
< dhc
->nlambda
; i
++)
589 mde_delta_h_add_dh(dhc
->dh_du
+i
, foreign_dU
[i
]);
591 if (dhc
->dh_pv
!= NULL
)
593 mde_delta_h_add_dh(dhc
->dh_pv
, pV
);
595 if (dhc
->dh_energy
!= NULL
)
597 mde_delta_h_add_dh(dhc
->dh_energy
, energy
);
599 if (dhc
->dh_expanded
!= NULL
)
601 mde_delta_h_add_dh(dhc
->dh_expanded
, fep_state
);
606 /* write the metadata associated with all the du blocks, and call
607 handle_block to write out all the du blocks */
608 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll
*dhc
,
609 t_enxframe
*fr
, int nblock
)
614 /* add one block with one subblock as the collection's own data */
616 add_blocks_enxframe(fr
, nblock
);
617 blk
= fr
->block
+ (nblock
-1);
619 /* only allocate lambda vector component blocks if they must be written out
620 for backward compatibility */
621 if (dhc
->native_lambda_components
!= NULL
)
623 add_subblocks_enxblock(blk
, 2);
627 add_subblocks_enxblock(blk
, 1);
630 dhc
->subblock_d
[0] = dhc
->temperature
; /* temperature */
631 dhc
->subblock_d
[1] = dhc
->start_time
; /* time of first sample */
632 dhc
->subblock_d
[2] = dhc
->delta_time
; /* time difference between samples */
633 dhc
->subblock_d
[3] = dhc
->start_lambda
; /* old-style lambda at starttime */
634 dhc
->subblock_d
[4] = dhc
->delta_lambda
; /* lambda diff. between samples */
635 /* set the lambda vector components if they exist */
636 if (dhc
->native_lambda_components
!= NULL
)
638 for (i
= 0; i
< dhc
->n_lambda_vec
; i
++)
640 dhc
->subblock_d
[5+i
] = dhc
->native_lambda_vec
[i
];
644 blk
->sub
[0].nr
= 5 + dhc
->n_lambda_vec
;
645 blk
->sub
[0].type
= xdr_datatype_double
;
646 blk
->sub
[0].dval
= dhc
->subblock_d
;
648 if (dhc
->native_lambda_components
!= NULL
)
650 dhc
->subblock_i
[0] = dhc
->lambda_index
;
651 /* set the lambda vector component IDs if they exist */
652 dhc
->subblock_i
[1] = dhc
->n_lambda_vec
;
653 for (i
= 0; i
< dhc
->n_lambda_vec
; i
++)
655 dhc
->subblock_i
[i
+2] = dhc
->native_lambda_components
[i
];
657 blk
->sub
[1].nr
= 2 + dhc
->n_lambda_vec
;
658 blk
->sub
[1].type
= xdr_datatype_int
;
659 blk
->sub
[1].ival
= dhc
->subblock_i
;
662 for (i
= 0; i
< dhc
->ndh
; i
++)
665 add_blocks_enxframe(fr
, nblock
);
666 blk
= fr
->block
+ (nblock
-1);
668 mde_delta_h_handle_block(dhc
->dh
+i
, blk
);
672 /* reset the data for a new round */
673 void mde_delta_h_coll_reset(t_mde_delta_h_coll
*dhc
)
676 for (i
= 0; i
< dhc
->ndh
; i
++)
678 if (dhc
->dh
[i
].written
)
680 /* we can now throw away the data */
681 mde_delta_h_reset(dhc
->dh
+ i
);
684 dhc
->start_time_set
= FALSE
;
687 /* set the energyhistory variables to save state */
688 void mde_delta_h_coll_update_energyhistory(t_mde_delta_h_coll
*dhc
,
689 energyhistory_t
*enerhist
)
694 snew(enerhist
->dht
, 1);
695 snew(enerhist
->dht
->ndh
, dhc
->ndh
);
696 snew(enerhist
->dht
->dh
, dhc
->ndh
);
697 enerhist
->dht
->nndh
= dhc
->ndh
;
701 if (enerhist
->dht
->nndh
!= dhc
->ndh
)
703 gmx_incons("energy history number of delta_h histograms != inputrec's number");
706 for (i
= 0; i
< dhc
->ndh
; i
++)
708 enerhist
->dht
->dh
[i
] = dhc
->dh
[i
].dh
;
709 enerhist
->dht
->ndh
[i
] = dhc
->dh
[i
].ndh
;
711 enerhist
->dht
->start_time
= dhc
->start_time
;
712 enerhist
->dht
->start_lambda
= dhc
->start_lambda
;
717 /* restore the variables from an energyhistory */
718 void mde_delta_h_coll_restore_energyhistory(t_mde_delta_h_coll
*dhc
,
719 energyhistory_t
*enerhist
)
726 gmx_incons("No delta_h histograms found");
730 gmx_incons("No delta_h histograms found in energy history");
732 if (enerhist
->dht
->nndh
!= dhc
->ndh
)
734 gmx_incons("energy history number of delta_h histograms != inputrec's number");
737 for (i
= 0; i
< enerhist
->dht
->nndh
; i
++)
739 dhc
->dh
[i
].ndh
= enerhist
->dht
->ndh
[i
];
740 for (j
= 0; j
< dhc
->dh
[i
].ndh
; j
++)
742 dhc
->dh
[i
].dh
[j
] = enerhist
->dht
->dh
[i
][j
];
745 dhc
->start_time
= enerhist
->dht
->start_time
;
746 if (enerhist
->dht
->start_lambda_set
)
748 dhc
->start_lambda
= enerhist
->dht
->start_lambda
;
750 if (dhc
->dh
[0].ndh
> 0)
752 dhc
->start_time_set
= TRUE
;
756 dhc
->start_time_set
= FALSE
;