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,2016,2017,2018, 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"
47 #include "gromacs/fileio/enxio.h"
48 #include "gromacs/mdlib/mdebin.h"
49 #include "gromacs/mdtypes/energyhistory.h"
50 #include "gromacs/mdtypes/inputrec.h"
51 #include "gromacs/mdtypes/md_enums.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/gmxassert.h"
54 #include "gromacs/utility/smalloc.h"
56 /* reset the delta_h list to prepare it for new values */
57 static void mde_delta_h_reset(t_mde_delta_h
*dh
)
63 /* initialize the delta_h list */
64 static void mde_delta_h_init(t_mde_delta_h
*dh
, int nbins
,
65 double dx
, unsigned int ndhmax
,
66 int type
, int derivative
, int nlambda
,
72 dh
->derivative
= derivative
;
73 dh
->nlambda
= nlambda
;
75 snew(dh
->lambda
, nlambda
);
76 for (i
= 0; i
< nlambda
; i
++)
79 dh
->lambda
[i
] = lambda
[i
];
83 snew(dh
->subblock_meta_d
, dh
->nlambda
+1);
85 dh
->ndhmax
= ndhmax
+2;
86 for (i
= 0; i
< 2; i
++)
91 snew(dh
->dh
, dh
->ndhmax
);
92 snew(dh
->dhf
, dh
->ndhmax
);
94 if (nbins
<= 0 || dx
< GMX_REAL_EPS
*10)
101 /* pre-allocate the histogram */
102 dh
->nhist
= 2; /* energies and derivatives histogram */
105 for (i
= 0; i
< dh
->nhist
; i
++)
107 snew(dh
->bin
[i
], dh
->nbins
);
110 mde_delta_h_reset(dh
);
113 static void done_mde_delta_h(t_mde_delta_h
*dh
)
116 sfree(dh
->subblock_meta_d
);
119 for (int i
= 0; i
< dh
->nhist
; i
++)
125 /* Add a value to the delta_h list */
126 static void mde_delta_h_add_dh(t_mde_delta_h
*dh
, double delta_h
)
128 if (dh
->ndh
>= dh
->ndhmax
)
130 gmx_incons("delta_h array not big enough!");
132 dh
->dh
[dh
->ndh
] = delta_h
;
136 /* construct histogram with index hi */
137 static void mde_delta_h_make_hist(t_mde_delta_h
*dh
, int hi
, gmx_bool invert
)
139 double min_dh
= FLT_MAX
;
140 double max_dh
= -FLT_MAX
;
142 double max_dh_hist
; /* maximum binnable dh value */
143 double min_dh_hist
; /* minimum binnable dh value */
145 double f
; /* energy mult. factor */
147 /* by applying a -1 scaling factor on the energies we get the same as
148 having a negative dx, but we don't need to fix the min/max values
149 beyond inverting x0 */
152 /* first find min and max */
153 for (i
= 0; i
< dh
->ndh
; i
++)
155 if (f
*dh
->dh
[i
] < min_dh
)
157 min_dh
= f
*dh
->dh
[i
];
159 if (f
*dh
->dh
[i
] > max_dh
)
161 max_dh
= f
*dh
->dh
[i
];
165 /* reset the histogram */
166 for (i
= 0; i
< dh
->nbins
; i
++)
172 /* The starting point of the histogram is the lowest value found:
173 that value has the highest contribution to the free energy.
175 Get this start value in number of histogram dxs from zero,
178 dh
->x0
[hi
] = (gmx_int64_t
)floor(min_dh
/dx
);
180 min_dh_hist
= (dh
->x0
[hi
])*dx
;
181 max_dh_hist
= (dh
->x0
[hi
] + dh
->nbins
+ 1)*dx
;
183 /* and fill the histogram*/
184 for (i
= 0; i
< dh
->ndh
; i
++)
188 /* Determine the bin number. If it doesn't fit into the histogram,
189 add it to the last bin.
190 We check the max_dh_int range because converting to integers
191 might lead to overflow with unpredictable results.*/
192 if ( (f
*dh
->dh
[i
] >= min_dh_hist
) && (f
*dh
->dh
[i
] <= max_dh_hist
) )
194 bin
= (unsigned int)( (f
*dh
->dh
[i
] - min_dh_hist
)/dx
);
200 /* double-check here because of possible round-off errors*/
201 if (bin
>= dh
->nbins
)
205 if (bin
> dh
->maxbin
[hi
])
207 dh
->maxbin
[hi
] = bin
;
213 /* make sure we include a bin with 0 if we didn't use the full
214 histogram width. This can then be used as an indication that
215 all the data was binned. */
216 if (dh
->maxbin
[hi
] < dh
->nbins
-1)
223 static void mde_delta_h_handle_block(t_mde_delta_h
*dh
, t_enxblock
*blk
)
225 /* first check which type we should use: histogram or raw data */
230 /* We write raw data.
231 Raw data consists of 3 subblocks: an int metadata block
232 with type and derivative index, a foreign lambda block
233 and and the data itself */
234 add_subblocks_enxblock(blk
, 3);
239 dh
->subblock_meta_i
[0] = dh
->type
; /* block data type */
240 dh
->subblock_meta_i
[1] = dh
->derivative
; /* derivative direction if
241 applicable (in indices
242 starting from first coord in
243 the main delta_h_coll) */
245 blk
->sub
[0].type
= xdr_datatype_int
;
246 blk
->sub
[0].ival
= dh
->subblock_meta_i
;
249 for (i
= 0; i
< dh
->nlambda
; i
++)
251 dh
->subblock_meta_d
[i
] = dh
->lambda
[i
];
253 blk
->sub
[1].nr
= dh
->nlambda
;
254 blk
->sub
[1].type
= xdr_datatype_double
;
255 blk
->sub
[1].dval
= dh
->subblock_meta_d
;
258 /* check if there's actual data to be written. */
264 blk
->sub
[2].nr
= dh
->ndh
;
265 /* Michael commented in 2012 that this use of explicit
266 xdr_datatype_float was good for F@H for now.
267 Apparently it's still good enough. */
268 blk
->sub
[2].type
= xdr_datatype_float
;
269 for (i
= 0; i
< dh
->ndh
; i
++)
271 dh
->dhf
[i
] = (float)dh
->dh
[i
];
273 blk
->sub
[2].fval
= dh
->dhf
;
279 blk
->sub
[2].type
= xdr_datatype_float
;
280 blk
->sub
[2].fval
= nullptr;
285 int nhist_written
= 0;
289 /* TODO histogram metadata */
290 /* check if there's actual data to be written. */
293 gmx_bool prev_complete
= FALSE
;
294 /* Make the histogram(s) */
295 for (i
= 0; i
< dh
->nhist
; i
++)
299 /* the first histogram is always normal, and the
300 second one is always reverse */
301 mde_delta_h_make_hist(dh
, i
, i
== 1);
303 /* check whether this histogram contains all data: if the
304 last bin is 0, it does */
305 if (dh
->bin
[i
][dh
->nbins
-1] == 0)
307 prev_complete
= TRUE
;
311 prev_complete
= TRUE
;
318 /* A histogram consists of 2, 3 or 4 subblocks:
319 the foreign lambda value + histogram spacing, the starting point,
320 and the histogram data (0, 1 or 2 blocks). */
321 add_subblocks_enxblock(blk
, nhist_written
+2);
324 /* subblock 1: the lambda value + the histogram spacing */
325 if (dh
->nlambda
== 1)
327 /* for backward compatibility */
328 dh
->subblock_meta_d
[0] = dh
->lambda
[0];
332 dh
->subblock_meta_d
[0] = -1;
333 for (i
= 0; i
< dh
->nlambda
; i
++)
335 dh
->subblock_meta_d
[2+i
] = dh
->lambda
[i
];
338 dh
->subblock_meta_d
[1] = dh
->dx
;
339 blk
->sub
[0].nr
= 2+ ((dh
->nlambda
> 1) ? dh
->nlambda
: 0);
340 blk
->sub
[0].type
= xdr_datatype_double
;
341 blk
->sub
[0].dval
= dh
->subblock_meta_d
;
343 /* subblock 2: the starting point(s) as a long integer */
344 dh
->subblock_meta_l
[0] = nhist_written
;
345 dh
->subblock_meta_l
[1] = dh
->type
; /*dh->derivative ? 1 : 0;*/
347 for (i
= 0; i
< nhist_written
; i
++)
349 dh
->subblock_meta_l
[k
++] = dh
->x0
[i
];
351 /* append the derivative data */
352 dh
->subblock_meta_l
[k
++] = dh
->derivative
;
354 blk
->sub
[1].nr
= nhist_written
+3;
355 blk
->sub
[1].type
= xdr_datatype_int64
;
356 blk
->sub
[1].lval
= dh
->subblock_meta_l
;
358 /* subblock 3 + 4 : the histogram data */
359 for (i
= 0; i
< nhist_written
; i
++)
361 blk
->sub
[i
+2].nr
= dh
->maxbin
[i
]+1; /* it's +1 because size=index+1
363 blk
->sub
[i
+2].type
= xdr_datatype_int
;
364 blk
->sub
[i
+2].ival
= dh
->bin
[i
];
369 /* initialize the collection*/
370 void mde_delta_h_coll_init(t_mde_delta_h_coll
*dhc
, const t_inputrec
*ir
)
374 int ndhmax
= ir
->nstenergy
/ir
->nstcalcenergy
;
375 t_lambda
*fep
= ir
->fepvals
;
377 dhc
->temperature
= ir
->opts
.ref_t
[0]; /* only store system temperature */
378 dhc
->start_time
= 0.;
379 dhc
->delta_time
= ir
->delta_t
*ir
->fepvals
->nstdhdl
;
380 dhc
->start_time_set
= FALSE
;
382 /* this is the compatibility lambda value. If it is >=0, it is valid,
383 and there is either an old-style lambda or a slow growth simulation. */
384 dhc
->start_lambda
= ir
->fepvals
->init_lambda
;
385 /* for continuous change of lambda values */
386 dhc
->delta_lambda
= ir
->fepvals
->delta_lambda
*ir
->fepvals
->nstdhdl
;
388 if (dhc
->start_lambda
< 0)
390 /* create the native lambda vectors */
391 dhc
->lambda_index
= fep
->init_fep_state
;
392 dhc
->n_lambda_vec
= 0;
393 for (i
= 0; i
< efptNR
; i
++)
395 if (fep
->separate_dvdl
[i
])
400 snew(dhc
->native_lambda_vec
, dhc
->n_lambda_vec
);
401 snew(dhc
->native_lambda_components
, dhc
->n_lambda_vec
);
403 for (i
= 0; i
< efptNR
; i
++)
405 if (fep
->separate_dvdl
[i
])
407 dhc
->native_lambda_components
[j
] = i
;
408 if (fep
->init_fep_state
>= 0 &&
409 fep
->init_fep_state
< fep
->n_lambda
)
411 dhc
->native_lambda_vec
[j
] =
412 fep
->all_lambda
[i
][fep
->init_fep_state
];
416 dhc
->native_lambda_vec
[j
] = -1;
424 /* don't allocate the meta-data subblocks for lambda vectors */
425 dhc
->native_lambda_vec
= nullptr;
426 dhc
->n_lambda_vec
= 0;
427 dhc
->native_lambda_components
= nullptr;
428 dhc
->lambda_index
= -1;
430 /* allocate metadata subblocks */
431 snew(dhc
->subblock_d
, 5 + dhc
->n_lambda_vec
);
432 snew(dhc
->subblock_i
, 1 + dhc
->n_lambda_vec
);
434 /* now decide which data to write out */
437 dhc
->dh_expanded
= nullptr;
438 dhc
->dh_energy
= nullptr;
439 dhc
->dh_pv
= nullptr;
441 /* total number of raw data point collections in the sample */
445 gmx_bool bExpanded
= FALSE
;
446 gmx_bool bEnergy
= FALSE
;
447 gmx_bool bPV
= FALSE
;
448 int n_lambda_components
= 0;
450 /* first count the number of states */
453 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
455 for (i
= 0; i
< efptNR
; i
++)
457 if (ir
->fepvals
->separate_dvdl
[i
])
464 /* add the lambdas */
465 dhc
->nlambda
= ir
->fepvals
->lambda_stop_n
- ir
->fepvals
->lambda_start_n
;
466 dhc
->ndh
+= dhc
->nlambda
;
467 /* another compatibility check */
468 if (dhc
->start_lambda
< 0)
470 /* include one more for the specification of the state, by lambda or
472 if (ir
->expandedvals
->elmcmove
> elmcmoveNO
)
477 /* whether to print energies */
478 if (ir
->fepvals
->edHdLPrintEnergy
!= edHdLPrintEnergyNO
)
485 dhc
->ndh
+= 1; /* include pressure-volume work */
490 snew(dhc
->dh
, dhc
->ndh
);
492 /* now initialize them */
493 /* the order, for now, must match that of the dhdl.xvg file because of
494 how g_energy -odh is implemented */
498 dhc
->dh_expanded
= dhc
->dh
+n
;
499 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
500 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
501 dhbtEXPANDED
, 0, 0, nullptr);
506 dhc
->dh_energy
= dhc
->dh
+n
;
507 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
508 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
509 dhbtEN
, 0, 0, nullptr);
513 n_lambda_components
= 0;
514 if (fep
->dhdl_derivatives
== edhdlderivativesYES
)
516 dhc
->dh_dhdl
= dhc
->dh
+ n
;
517 for (i
= 0; i
< efptNR
; i
++)
519 if (ir
->fepvals
->separate_dvdl
[i
])
521 /* we give it init_lambda for compatibility */
522 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
523 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
524 dhbtDHDL
, n_lambda_components
, 1,
525 &(fep
->init_lambda
));
527 n_lambda_components
++;
533 for (i
= 0; i
< efptNR
; i
++)
535 if (ir
->fepvals
->separate_dvdl
[i
])
537 n_lambda_components
++; /* count the components */
542 /* add the lambdas */
543 dhc
->dh_du
= dhc
->dh
+ n
;
544 snew(lambda_vec
, n_lambda_components
);
545 for (i
= ir
->fepvals
->lambda_start_n
; i
< ir
->fepvals
->lambda_stop_n
; i
++)
549 for (j
= 0; j
< efptNR
; j
++)
551 if (ir
->fepvals
->separate_dvdl
[j
])
553 lambda_vec
[k
++] = fep
->all_lambda
[j
][i
];
557 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
558 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
559 dhbtDH
, 0, n_lambda_components
, lambda_vec
);
565 dhc
->dh_pv
= dhc
->dh
+n
;
566 mde_delta_h_init(dhc
->dh
+n
, ir
->fepvals
->dh_hist_size
,
567 ir
->fepvals
->dh_hist_spacing
, ndhmax
,
568 dhbtPV
, 0, 0, nullptr);
574 void done_mde_delta_h_coll(t_mde_delta_h_coll
*dhc
)
580 sfree(dhc
->native_lambda_vec
);
581 sfree(dhc
->native_lambda_components
);
582 sfree(dhc
->subblock_d
);
583 sfree(dhc
->subblock_i
);
584 for (int i
= 0; i
< dhc
->ndh
; ++i
)
586 done_mde_delta_h(&dhc
->dh
[i
]);
592 /* add a bunch of samples - note fep_state is double to allow for better data storage */
593 void mde_delta_h_coll_add_dh(t_mde_delta_h_coll
*dhc
,
603 if (!dhc
->start_time_set
)
605 dhc
->start_time_set
= TRUE
;
606 dhc
->start_time
= time
;
609 for (i
= 0; i
< dhc
->ndhdl
; i
++)
611 mde_delta_h_add_dh(dhc
->dh_dhdl
+i
, dhdl
[i
]);
613 for (i
= 0; i
< dhc
->nlambda
; i
++)
615 mde_delta_h_add_dh(dhc
->dh_du
+i
, foreign_dU
[i
]);
617 if (dhc
->dh_pv
!= nullptr)
619 mde_delta_h_add_dh(dhc
->dh_pv
, pV
);
621 if (dhc
->dh_energy
!= nullptr)
623 mde_delta_h_add_dh(dhc
->dh_energy
, energy
);
625 if (dhc
->dh_expanded
!= nullptr)
627 mde_delta_h_add_dh(dhc
->dh_expanded
, fep_state
);
632 /* write the metadata associated with all the du blocks, and call
633 handle_block to write out all the du blocks */
634 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll
*dhc
,
635 t_enxframe
*fr
, int nblock
)
640 /* add one block with one subblock as the collection's own data */
642 add_blocks_enxframe(fr
, nblock
);
643 blk
= fr
->block
+ (nblock
-1);
645 /* only allocate lambda vector component blocks if they must be written out
646 for backward compatibility */
647 if (dhc
->native_lambda_components
!= nullptr)
649 add_subblocks_enxblock(blk
, 2);
653 add_subblocks_enxblock(blk
, 1);
656 dhc
->subblock_d
[0] = dhc
->temperature
; /* temperature */
657 dhc
->subblock_d
[1] = dhc
->start_time
; /* time of first sample */
658 dhc
->subblock_d
[2] = dhc
->delta_time
; /* time difference between samples */
659 dhc
->subblock_d
[3] = dhc
->start_lambda
; /* old-style lambda at starttime */
660 dhc
->subblock_d
[4] = dhc
->delta_lambda
; /* lambda diff. between samples */
661 /* set the lambda vector components if they exist */
662 if (dhc
->native_lambda_components
!= nullptr)
664 for (i
= 0; i
< dhc
->n_lambda_vec
; i
++)
666 dhc
->subblock_d
[5+i
] = dhc
->native_lambda_vec
[i
];
670 blk
->sub
[0].nr
= 5 + dhc
->n_lambda_vec
;
671 blk
->sub
[0].type
= xdr_datatype_double
;
672 blk
->sub
[0].dval
= dhc
->subblock_d
;
674 if (dhc
->native_lambda_components
!= nullptr)
676 dhc
->subblock_i
[0] = dhc
->lambda_index
;
677 /* set the lambda vector component IDs if they exist */
678 dhc
->subblock_i
[1] = dhc
->n_lambda_vec
;
679 for (i
= 0; i
< dhc
->n_lambda_vec
; i
++)
681 dhc
->subblock_i
[i
+2] = dhc
->native_lambda_components
[i
];
683 blk
->sub
[1].nr
= 2 + dhc
->n_lambda_vec
;
684 blk
->sub
[1].type
= xdr_datatype_int
;
685 blk
->sub
[1].ival
= dhc
->subblock_i
;
688 for (i
= 0; i
< dhc
->ndh
; i
++)
691 add_blocks_enxframe(fr
, nblock
);
692 blk
= fr
->block
+ (nblock
-1);
694 mde_delta_h_handle_block(dhc
->dh
+i
, blk
);
698 /* reset the data for a new round */
699 void mde_delta_h_coll_reset(t_mde_delta_h_coll
*dhc
)
702 for (i
= 0; i
< dhc
->ndh
; i
++)
704 if (dhc
->dh
[i
].written
)
706 /* we can now throw away the data */
707 mde_delta_h_reset(dhc
->dh
+ i
);
710 dhc
->start_time_set
= FALSE
;
713 /* set the energyhistory variables to save state */
714 void mde_delta_h_coll_update_energyhistory(const t_mde_delta_h_coll
*dhc
,
715 energyhistory_t
*enerhist
)
717 if (enerhist
->deltaHForeignLambdas
== nullptr)
719 enerhist
->deltaHForeignLambdas
.reset(new delta_h_history_t
);
720 enerhist
->deltaHForeignLambdas
->dh
.resize(dhc
->ndh
);
723 delta_h_history_t
* const deltaH
= enerhist
->deltaHForeignLambdas
.get();
725 GMX_RELEASE_ASSERT(deltaH
->dh
.size() == static_cast<size_t>(dhc
->ndh
), "energy history number of delta_h histograms should match inputrec's number");
727 for (int i
= 0; i
< dhc
->ndh
; i
++)
729 std::vector
<real
> &dh
= deltaH
->dh
[i
];
730 dh
.resize(dhc
->dh
[i
].ndh
);
731 std::copy(dh
.begin(), dh
.end(), dhc
->dh
[i
].dh
);
733 deltaH
->start_time
= dhc
->start_time
;
734 deltaH
->start_lambda
= dhc
->start_lambda
;
739 /* restore the variables from an energyhistory */
740 void mde_delta_h_coll_restore_energyhistory(t_mde_delta_h_coll
*dhc
,
741 const delta_h_history_t
*deltaH
)
743 GMX_RELEASE_ASSERT(dhc
, "Should have delta_h histograms");
744 GMX_RELEASE_ASSERT(deltaH
, "Should have delta_h histograms in energy history");
745 GMX_RELEASE_ASSERT(deltaH
->dh
.size() == static_cast<size_t>(dhc
->ndh
), "energy history number of delta_h histograms should match inputrec's number");
747 for (unsigned int i
= 0; i
< deltaH
->dh
.size(); i
++)
749 dhc
->dh
[i
].ndh
= deltaH
->dh
[i
].size();
750 for (unsigned int j
= 0; j
< dhc
->dh
[i
].ndh
; j
++)
752 dhc
->dh
[i
].dh
[j
] = deltaH
->dh
[i
][j
];
755 dhc
->start_time
= deltaH
->start_time
;
756 if (deltaH
->start_lambda_set
)
758 dhc
->start_lambda
= deltaH
->start_lambda
;
760 if (dhc
->dh
[0].ndh
> 0)
762 dhc
->start_time_set
= TRUE
;
766 dhc
->start_time_set
= FALSE
;