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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
40 #include "gromacs/utility/arrayref.h"
41 #include "mdebin_bar.h"
50 #include "gromacs/fileio/enxio.h"
51 #include "gromacs/mdlib/energyoutput.h"
52 #include "gromacs/mdtypes/energyhistory.h"
53 #include "gromacs/mdtypes/inputrec.h"
54 #include "gromacs/mdtypes/md_enums.h"
55 #include "gromacs/trajectory/energyframe.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/gmxassert.h"
58 #include "gromacs/utility/smalloc.h"
60 /* Number of entries in subblock_d preceding the lambda components */
61 constexpr int c_subblockDNumPreEntries = 5;
62 /* Number of entries in subblock_i preceding the lambda components */
63 constexpr int c_subblockINumPreEntries = 2;
65 /* reset the delta_h list to prepare it for new values */
66 static void mde_delta_h_reset(t_mde_delta_h* dh)
72 /* initialize the delta_h list */
73 static void mde_delta_h_init(t_mde_delta_h* dh,
85 dh->derivative = derivative;
86 dh->nlambda = nlambda;
88 dh->lambda.resize(nlambda);
89 for (i = 0; i < nlambda; i++)
92 dh->lambda[i] = lambda[i];
96 dh->subblock_meta_d.resize(dh->nlambda + 1);
98 dh->ndhmax = ndhmax + 2;
100 dh->dh.resize(dh->ndhmax);
101 dh->dhf.resize(dh->ndhmax);
103 if (nbins <= 0 || dx < GMX_REAL_EPS * 10)
110 /* pre-allocate the histogram */
111 dh->nhist = 2; /* energies and derivatives histogram */
114 for (i = 0; i < dh->nhist; i++)
116 dh->bin[i].resize(dh->nbins);
119 mde_delta_h_reset(dh);
122 /* Add a value to the delta_h list */
123 static void mde_delta_h_add_dh(t_mde_delta_h* dh, double delta_h)
125 if (dh->ndh >= dh->ndhmax)
127 gmx_incons("delta_h array not big enough!");
129 dh->dh[dh->ndh] = delta_h;
133 /* construct histogram with index hi */
134 static void mde_delta_h_make_hist(t_mde_delta_h* dh, int hi, gmx_bool invert)
136 double min_dh = FLT_MAX;
137 double max_dh = -FLT_MAX;
139 double max_dh_hist; /* maximum binnable dh value */
140 double min_dh_hist; /* minimum binnable dh value */
142 double f; /* energy mult. factor */
144 /* by applying a -1 scaling factor on the energies we get the same as
145 having a negative dx, but we don't need to fix the min/max values
146 beyond inverting x0 */
149 /* first find min and max */
150 for (i = 0; i < dh->ndh; i++)
152 if (f * dh->dh[i] < min_dh)
154 min_dh = f * dh->dh[i];
156 if (f * dh->dh[i] > max_dh)
158 max_dh = f * dh->dh[i];
162 /* reset the histogram */
163 for (i = 0; i < dh->nbins; i++)
169 /* The starting point of the histogram is the lowest value found:
170 that value has the highest contribution to the free energy.
172 Get this start value in number of histogram dxs from zero,
175 dh->x0[hi] = static_cast<int64_t>(floor(min_dh / dx));
177 min_dh_hist = (dh->x0[hi]) * dx;
178 max_dh_hist = (dh->x0[hi] + dh->nbins + 1) * dx;
180 /* and fill the histogram*/
181 for (i = 0; i < dh->ndh; i++)
185 /* Determine the bin number. If it doesn't fit into the histogram,
186 add it to the last bin.
187 We check the max_dh_int range because converting to integers
188 might lead to overflow with unpredictable results.*/
189 if ((f * dh->dh[i] >= min_dh_hist) && (f * dh->dh[i] <= max_dh_hist))
191 bin = static_cast<unsigned int>((f * dh->dh[i] - min_dh_hist) / dx);
197 /* double-check here because of possible round-off errors*/
198 if (bin >= dh->nbins)
202 if (bin > dh->maxbin[hi])
204 dh->maxbin[hi] = bin;
210 /* make sure we include a bin with 0 if we didn't use the full
211 histogram width. This can then be used as an indication that
212 all the data was binned. */
213 if (dh->maxbin[hi] < dh->nbins - 1)
220 static void mde_delta_h_handle_block(t_mde_delta_h* dh, t_enxblock* blk)
222 /* first check which type we should use: histogram or raw data */
227 /* We write raw data.
228 Raw data consists of 3 subblocks: an int metadata block
229 with type and derivative index, a foreign lambda block
230 and and the data itself */
231 add_subblocks_enxblock(blk, 3);
236 dh->subblock_meta_i[0] = dh->type; /* block data type */
237 dh->subblock_meta_i[1] = dh->derivative; /* derivative direction if
238 applicable (in indices
239 starting from first coord in
240 the main delta_h_coll) */
242 blk->sub[0].type = XdrDataType::Int;
243 blk->sub[0].ival = dh->subblock_meta_i.data();
246 for (i = 0; i < dh->nlambda; i++)
248 dh->subblock_meta_d[i] = dh->lambda[i];
250 blk->sub[1].nr = dh->nlambda;
251 blk->sub[1].type = XdrDataType::Double;
252 blk->sub[1].dval = dh->subblock_meta_d.data();
255 /* check if there's actual data to be written. */
261 blk->sub[2].nr = dh->ndh;
262 /* Michael commented in 2012 that this use of explicit
263 XdrDataType::Float was good for F@H for now.
264 Apparently it's still good enough. */
265 blk->sub[2].type = XdrDataType::Float;
266 for (i = 0; i < dh->ndh; i++)
268 dh->dhf[i] = static_cast<float>(dh->dh[i]);
270 blk->sub[2].fval = dh->dhf.data();
276 blk->sub[2].type = XdrDataType::Float;
277 blk->sub[2].fval = nullptr;
282 int nhist_written = 0;
286 /* TODO histogram metadata */
287 /* check if there's actual data to be written. */
290 gmx_bool prev_complete = FALSE;
291 /* Make the histogram(s) */
292 for (i = 0; i < dh->nhist; i++)
296 /* the first histogram is always normal, and the
297 second one is always reverse */
298 mde_delta_h_make_hist(dh, i, i == 1);
300 /* check whether this histogram contains all data: if the
301 last bin is 0, it does */
302 if (dh->bin[i][dh->nbins - 1] == 0)
304 prev_complete = TRUE;
308 prev_complete = TRUE;
315 /* A histogram consists of 2, 3 or 4 subblocks:
316 the foreign lambda value + histogram spacing, the starting point,
317 and the histogram data (0, 1 or 2 blocks). */
318 add_subblocks_enxblock(blk, nhist_written + 2);
321 /* subblock 1: the lambda value + the histogram spacing */
322 if (dh->nlambda == 1)
324 /* for backward compatibility */
325 dh->subblock_meta_d[0] = dh->lambda[0];
329 dh->subblock_meta_d[0] = -1;
330 for (i = 0; i < dh->nlambda; i++)
332 dh->subblock_meta_d[2 + i] = dh->lambda[i];
335 dh->subblock_meta_d[1] = dh->dx;
336 blk->sub[0].nr = 2 + ((dh->nlambda > 1) ? dh->nlambda : 0);
337 blk->sub[0].type = XdrDataType::Double;
338 blk->sub[0].dval = dh->subblock_meta_d.data();
340 /* subblock 2: the starting point(s) as a long integer */
341 dh->subblock_meta_l[0] = nhist_written;
342 dh->subblock_meta_l[1] = dh->type; /*dh->derivative ? 1 : 0;*/
344 for (i = 0; i < nhist_written; i++)
346 dh->subblock_meta_l[k++] = dh->x0[i];
348 /* append the derivative data */
349 dh->subblock_meta_l[k++] = dh->derivative;
351 blk->sub[1].nr = nhist_written + 3;
352 blk->sub[1].type = XdrDataType::Int64;
353 blk->sub[1].lval = dh->subblock_meta_l.data();
355 /* subblock 3 + 4 : the histogram data */
356 for (i = 0; i < nhist_written; i++)
358 blk->sub[i + 2].nr = dh->maxbin[i] + 1; /* it's +1 because size=index+1
360 blk->sub[i + 2].type = XdrDataType::Int;
361 blk->sub[i + 2].ival = dh->bin[i].data();
366 /* initialize the collection*/
367 t_mde_delta_h_coll::t_mde_delta_h_coll(const t_inputrec& inputrec)
371 int ndhmax = inputrec.nstenergy / inputrec.nstcalcenergy;
372 t_lambda* fep = inputrec.fepvals.get();
374 temperature = inputrec.opts.ref_t[0]; /* only store system temperature */
376 delta_time = inputrec.delta_t * inputrec.fepvals->nstdhdl;
377 start_time_set = FALSE;
379 /* this is the compatibility lambda value. If it is >=0, it is valid,
380 and there is either an old-style lambda or a slow growth simulation. */
381 start_lambda = inputrec.fepvals->init_lambda;
382 /* for continuous change of lambda values */
383 delta_lambda = inputrec.fepvals->delta_lambda * inputrec.fepvals->nstdhdl;
385 if (start_lambda < 0)
387 /* create the native lambda vectors */
388 lambda_index = fep->init_fep_state;
390 for (auto i : keysOf(fep->separate_dvdl))
392 if (fep->separate_dvdl[i])
397 native_lambda_vec.resize(n_lambda_vec);
398 native_lambda_components.resize(n_lambda_vec);
400 for (auto i : keysOf(fep->separate_dvdl))
402 if (fep->separate_dvdl[i])
404 native_lambda_components[j] = static_cast<int>(i);
405 if (fep->init_fep_state >= 0 && fep->init_fep_state < fep->n_lambda)
407 native_lambda_vec[j] = fep->all_lambda[i][fep->init_fep_state];
411 native_lambda_vec[j] = -1;
419 /* don't allocate the meta-data subblocks for lambda vectors */
423 /* allocate metadata subblocks */
424 subblock_d.resize(c_subblockDNumPreEntries + n_lambda_vec);
425 subblock_i.resize(c_subblockINumPreEntries + n_lambda_vec);
427 /* now decide which data to write out */
430 dh_expanded_index = -1;
431 dh_energy_index = -1;
434 /* total number of raw data point collections in the sample */
438 gmx_bool bExpanded = FALSE;
439 gmx_bool bEnergy = FALSE;
440 gmx_bool bPV = FALSE;
441 int n_lambda_components = 0;
443 /* first count the number of states */
446 if (fep->dhdl_derivatives == DhDlDerivativeCalculation::Yes)
448 for (auto i : keysOf(fep->separate_dvdl))
450 if (inputrec.fepvals->separate_dvdl[i])
457 /* add the lambdas */
458 nlambda = inputrec.fepvals->lambda_stop_n - inputrec.fepvals->lambda_start_n;
460 /* another compatibility check */
461 if (start_lambda < 0)
463 /* include one more for the specification of the state, by lambda or
465 if (inputrec.expandedvals->elmcmove > LambdaMoveCalculation::No)
470 /* whether to print energies */
471 if (inputrec.fepvals->edHdLPrintEnergy != FreeEnergyPrintEnergy::No)
476 if (inputrec.epc > PressureCoupling::No)
478 ndh += 1; /* include pressure-volume work */
485 /* now initialize them */
486 /* the order, for now, must match that of the dhdl.xvg file because of
487 how g_energy -odh is implemented */
491 dh_expanded_index = n;
492 mde_delta_h_init(&dh[n],
493 inputrec.fepvals->dh_hist_size,
494 inputrec.fepvals->dh_hist_spacing,
505 mde_delta_h_init(&dh[n],
506 inputrec.fepvals->dh_hist_size,
507 inputrec.fepvals->dh_hist_spacing,
516 n_lambda_components = 0;
517 if (fep->dhdl_derivatives == DhDlDerivativeCalculation::Yes)
520 for (auto i : keysOf(fep->separate_dvdl))
522 if (inputrec.fepvals->separate_dvdl[i])
524 /* we give it init_lambda for compatibility */
525 mde_delta_h_init(&dh[n],
526 inputrec.fepvals->dh_hist_size,
527 inputrec.fepvals->dh_hist_spacing,
532 &(fep->init_lambda));
534 n_lambda_components++;
540 for (auto i : keysOf(fep->separate_dvdl))
542 if (fep->separate_dvdl[i])
544 n_lambda_components++; /* count the components */
548 /* add the lambdas */
550 snew(lambda_vec, n_lambda_components);
551 for (i = inputrec.fepvals->lambda_start_n; i < inputrec.fepvals->lambda_stop_n; i++)
555 for (auto j : keysOf(fep->separate_dvdl))
557 if (fep->separate_dvdl[j])
559 lambda_vec[k++] = fep->all_lambda[j][i];
563 mde_delta_h_init(&dh[n],
564 inputrec.fepvals->dh_hist_size,
565 inputrec.fepvals->dh_hist_spacing,
577 mde_delta_h_init(&dh[n],
578 inputrec.fepvals->dh_hist_size,
579 inputrec.fepvals->dh_hist_spacing,
590 /* add a bunch of samples - note fep_state is double to allow for better data storage */
591 void mde_delta_h_coll_add_dh(t_mde_delta_h_coll* dhc,
595 gmx::ArrayRef<double> dhdl,
601 if (!dhc->start_time_set)
603 dhc->start_time_set = TRUE;
604 dhc->start_time = time;
607 for (i = 0; i < dhc->ndhdl; i++)
609 mde_delta_h_add_dh(&dhc->dh[dhc->dh_dhdl_index + i], dhdl[i]);
611 for (i = 0; i < dhc->nlambda; i++)
613 mde_delta_h_add_dh(&dhc->dh[dhc->dh_du_index + i], foreign_dU[i]);
615 if (dhc->dh_pv_index >= 0)
617 mde_delta_h_add_dh(&dhc->dh[dhc->dh_pv_index], pV);
619 if (dhc->dh_energy_index >= 0)
621 mde_delta_h_add_dh(&dhc->dh[dhc->dh_energy_index], energy);
623 if (dhc->dh_expanded_index >= 0)
625 mde_delta_h_add_dh(&dhc->dh[dhc->dh_expanded_index], fep_state);
629 /* write the metadata associated with all the du blocks, and call
630 handle_block to write out all the du blocks */
631 void mde_delta_h_coll_handle_block(t_mde_delta_h_coll* dhc, t_enxframe* fr, int nblock)
636 /* add one block with one subblock as the collection's own data */
638 add_blocks_enxframe(fr, nblock);
639 blk = fr->block + (nblock - 1);
641 /* only allocate lambda vector component blocks if they must be written out
642 for backward compatibility */
643 if (!dhc->native_lambda_components.empty())
645 add_subblocks_enxblock(blk, 2);
649 add_subblocks_enxblock(blk, 1);
652 dhc->subblock_d[0] = dhc->temperature; /* temperature */
653 dhc->subblock_d[1] = dhc->start_time; /* time of first sample */
654 dhc->subblock_d[2] = dhc->delta_time; /* time difference between samples */
655 dhc->subblock_d[3] = dhc->start_lambda; /* old-style lambda at starttime */
656 dhc->subblock_d[4] = dhc->delta_lambda; /* lambda diff. between samples */
657 /* set the lambda vector components if they exist */
658 if (!dhc->native_lambda_components.empty())
660 for (i = 0; i < dhc->n_lambda_vec; i++)
662 dhc->subblock_d[c_subblockDNumPreEntries + i] = dhc->native_lambda_vec[i];
666 blk->sub[0].nr = c_subblockDNumPreEntries + dhc->n_lambda_vec;
667 blk->sub[0].type = XdrDataType::Double;
668 blk->sub[0].dval = dhc->subblock_d.data();
670 if (!dhc->native_lambda_components.empty())
672 dhc->subblock_i[0] = dhc->lambda_index;
673 /* set the lambda vector component IDs if they exist */
674 dhc->subblock_i[1] = dhc->n_lambda_vec;
675 for (i = 0; i < dhc->n_lambda_vec; i++)
677 dhc->subblock_i[c_subblockINumPreEntries + i] = dhc->native_lambda_components[i];
679 blk->sub[1].nr = c_subblockINumPreEntries + dhc->n_lambda_vec;
680 blk->sub[1].type = XdrDataType::Int;
681 blk->sub[1].ival = dhc->subblock_i.data();
684 for (i = 0; i < dhc->ndh; i++)
687 add_blocks_enxframe(fr, nblock);
688 blk = fr->block + (nblock - 1);
690 mde_delta_h_handle_block(&dhc->dh[i], blk);
694 /* reset the data for a new round */
695 void mde_delta_h_coll_reset(t_mde_delta_h_coll* dhc)
698 for (i = 0; i < dhc->ndh; i++)
700 if (dhc->dh[i].written)
702 /* we can now throw away the data */
703 mde_delta_h_reset(&dhc->dh[i]);
706 dhc->start_time_set = FALSE;
709 /* set the energyhistory variables to save state */
710 void mde_delta_h_coll_update_energyhistory(const t_mde_delta_h_coll* dhc, energyhistory_t* enerhist)
712 if (enerhist->deltaHForeignLambdas == nullptr)
714 enerhist->deltaHForeignLambdas = std::make_unique<delta_h_history_t>();
715 enerhist->deltaHForeignLambdas->dh.resize(dhc->ndh);
718 delta_h_history_t* const deltaH = enerhist->deltaHForeignLambdas.get();
721 deltaH->dh.size() == static_cast<size_t>(dhc->ndh),
722 "energy history number of delta_h histograms should match inputrec's number");
724 for (int i = 0; i < dhc->ndh; i++)
726 std::vector<real>& dh = deltaH->dh[i];
727 for (unsigned int j = 0; j < dhc->dh[i].ndh; j++)
729 dh.emplace_back(dhc->dh[i].dh[j]);
732 deltaH->start_time = dhc->start_time;
733 deltaH->start_lambda = dhc->start_lambda;
737 /* restore the variables from an energyhistory */
738 void mde_delta_h_coll_restore_energyhistory(t_mde_delta_h_coll* dhc, const delta_h_history_t* deltaH)
740 GMX_RELEASE_ASSERT(dhc, "Should have delta_h histograms");
741 GMX_RELEASE_ASSERT(deltaH, "Should have delta_h histograms in energy history");
743 deltaH->dh.size() == static_cast<size_t>(dhc->ndh),
744 "energy history number of delta_h histograms should match inputrec's number");
746 for (gmx::index i = 0; i < gmx::ssize(deltaH->dh); i++)
748 dhc->dh[i].ndh = deltaH->dh[i].size();
749 for (unsigned int j = 0; j < dhc->dh[i].ndh; j++)
751 dhc->dh[i].dh[j] = deltaH->dh[i][j];
754 dhc->start_time = deltaH->start_time;
755 if (deltaH->start_lambda_set)
757 dhc->start_lambda = deltaH->start_lambda;
759 dhc->start_time_set = (dhc->dh[0].ndh > 0);