Restyling for reweighting implementation

[invemichele]

Hi everybody,
I think there is room for improvement in the current implementation of the reweighting, and I would be happy to contribute.
In particular the way the integral for c(t) is now calculated is, I believe, unnecessarily cumbersome. It comes directly from the Grid::integrate() function, which is not used anywhere else in the code (or at least not in the master branch).
I think that a simple direct summation over the grid points is faster, works in any dimension (solving #328 ), and does not have any downside effect since the reweighting is only available when using a grid bias.
This also avoids calling getBiasAndDerivatives which unnecessarily returns the derivatives.

I implemented this changes in a pull request, together with the support for the non well-tempered case (I don't see why not having it) and a fix for the overflow of c(t) in case of very high bias values.
Please let me know if I am missing a piece of the picture, and I am misunderstanding the code.

More in general I believe that also the user interface should be changed. In particular it does not make much sense to calculate the reweighting on a grid different from the bias one, and thus ask the user for REWEIGHTING_NGRID. especially because (if I understand the code correctly) the method called is Grid::getValue(const vector<double> & x), which calls Grid::getIndices(const vector<double> & x), which returns floor((x[i]-min_[i])/dx_[i]), thus as a matter of fact, if we double the number of reweighting bins, there is no difference in the value of c(t), except for some small border effect.
Thus REWEIGHTING_NGRID could be dropped in favor of a flag like CALC_RCT or CALC_REWEIGHTING. It is fine to leave the REWEIGHTING_NHILLS option, but I would change the default frequency to 1, and let the user choose if increasing it.
Of course these are just suggestions, driven by my typical usage, but might not be best suited for all the possible scenarios.

I hope this can be useful.

[valsson]

Hi @invemichele,

I would agree with simplifying the integral calculation as you suggest.

I just have one suggestion on the changes. Currently the integration is just done by summing over the grid points (and as c(t) is a ratio of two integrals there is no need to consider the volume of the grid bins). However, strictly speaking this is not correct numerical integration, it would be better to use something like trapezoidal rule for the integration so that the points at the edges have different weights. In many cases the integrated functions might be zero at the edges of the grid so this would not matter. But anyway I think it would be better to do it exactly.

In the VES code I have done all integrals by using trapezoidal rule. For the multidimensional case is just a product of one-dimensional trapezoidal rules. There is a function in src/ves/GridIntegrationWeights.cpp that for a given grid just gives a std::vector of the size of the grid with all the integration weights. It also handles the periodic case correctly. I could add that function to the Grid class. Then one could do something like:

std::vector weights = BiasGrid_->getIntegrationWeights();
for (Grid::index_t t=rank; t<BiasGrid_->getSize(); t+=stride) {
const double val=BiasGrid_->getValue(t);
Z_0+= weights[t] * std::exp(minusBetaFval-big_number);
Z_V+= weights[t] * std::exp(minusBetaFplusVval-big_number);
}

As for the changing the keywords, it would not make sense to keep the REWEIGHTING_NGRID keywords. But the developers should comment what is the exact policy of "breaking inputs".

As for changing the default for REWEIGHTING_NHILLS from 50 to 1. I am not sure about that, there is of course some overhead in doing the integration (especially for more than 2 dimension) and perhaps doing the calculation of c(t) every time a Gaussian is added it is not needed as the changes in c(t) are not that great. Or do you have experience where it is really important to calculate c(t) at every Gaussian addition?

Regards,
Omar

[invemichele]

Ciao Omar,
thanks for the quick reply. I also thought about using trapezoidal rule for the integral, but then as you wrote I thought that from a practical point of view the integrand will always vanish at the border. Also in the case of a periodic CV the values at the border are the same thus trapezoidal weights sum up to 1, and since in plumed grid the max value is not set in case of periodic CV, so it is correct to give weight 1 to the min value.
I don't know if I have been clear, but the idea is that in this specific case the simple summation is equivalent to the trapezoidal rule.
I also thought about Simpson's rule, but then I read this on wikipedia, where they claim that for our purposes (periodic or rapidly vanishing functions) trapezoidal rule is actually more accurate.

Anyway, it is probably a good idea to implement a proper integration method in the grid class, as you suggested, and in that case it could be used here as well.

Regarding the choice for default REWEIGHTING_NHILLS: no, I do not have a strong reason to suggest to put the default to 1. It is mainly because that is the correct thing to do in theory and in my simulations it does not give big overhead. But as I said this might not be a good choice for the average user. I leave that decision to more experienced MetaD users.

[carlocamilloni]

Ciao,

Sounds very good.
I would go for a CALC_RCT flag, while I would not change the default for REWEIGHTING_NHILLS at the moment, but its documentation can be expanded a bit to suggest to try the effect of lower values on the performances.
About the pull request it should be astyled, and the two relevant regtests should be reset.

[carlocamilloni]

I agree it would be good to have the integrate method in the grid class.

[GiovanniBussi]

Hi all, please also check #399 which looks related for what I understand (or maybe not).

I am not sure I am following everything about the integration stuff. Anyway, the way grids are implemented in plumed you can write the potential as an analytical function within each grid "cube". Maybe you might want to try to integrate that function? Or maybe the methods you proposed are equally good or better.

[invemichele]

Hi @GiovanniBussi
yes, it is related to #399, because mine would be an alternative solution to the same issue #328, and the part of the code referred to in #399 would be cancelled.

You are right, my integration method does not consider the spline interpolation. Then probably the small difference I saw in 1D when doubling REWEIGHTING_NGRID was due to the spline and not just to some border effect.
Still I don't know if the best thing is to define a new grid for the calculation of c(t), as it is now, since it is much slower and most of the users keep it equal to the bias one anyway. Maybe it is not too hard to include in an exact way the spline in the integration (i'm not familiar with splines, but I can try to look into this). Or we could keep this small (at least in my tests) numerical difference, as we do with the default REWEIGHTING_NHILLS=50 for the increase in performance.

A last possibility would be to keep both implementation, and switch to the faster one only in case the user sets REWEIGHTING_NGRID to be equal to GRID_BIN, since in that case they are equivalent. In this scenario the correction from #399 should be included.

[invemichele]

Hi everybody,
I did some quick testing, just to get a better idea of the differences. I ran 20 ns alanine dipeptide with bias on phi and psi, as in Belfast tutorial, with a 63x63 bias grid (GRID_SPACING=0.1,0.1).
I compared the current "old" implementation, with and without bug #399, and the proposed "new" implementation. As a reference I used old implementation without bug and with REWEIGHTING_NGRID=630,630, which should give a very good estimate, accounting for spline interpolation.
On y axis is the relative error in percentage, as a function of time: abs(rct-rct_ref)/rct_ref*100

Old implementation without bug, and new implementation are numerically indistinguishable when REWEIGHTING_NGRID=GRID_BIN .
Here some the wall times (I ran only once, so this numbers are not too reliable):

• 09:37 without rct calculation
• 09:58 new implementation, NHILLS=1 (cyano points)
• 10:47 old implementation, no bug, NHILLS=50, NGRID=63,63 (violet points)
• 12:29 old implementation, with bug, NHILLS=1, NGRID=63,63 (green points)
• 1:57:17 old implementation, no bug, NHILLS=1, NGRID=630,630 (reference)

I don't know if there are cases in which ignoring the spline interpolation can give rise to big differences, for alanine dipeptide this is not the case, there is virtually no difference at all.

I also ran some simulations with a much bigger bias grid GRID_BIN=300,300 (again, I ran only once each simulation):

• 12:27 without rct calculation
• 12:43 new implementation, NHILLS=1
• 54:40 old implementation, NHILLS=1

From these tests looks like with the new implementation there is no overhead in keeping NHILLS=1.

I made some changes to my previous pull request:

• changed user input form REWEIGHTING_NGRID and REWEIGHTING_NHILLS to CALC_RCT and RCT_USTRIDE
• made c(t) update stride=1 by default instead of 50
• added some comments and doc to the code
• astyle and regtest update

To be honest I don't know, probably I exaggerated since I broke backward compatibility, and I arbitrary chose to take away the possibility to increase the precision of the rct calculation (to not confuse the user, since it is not a perceivable difference, as far as I understand). Of course feel free to change these things!

Probably the bug fix #399 should be included in old plumed versions.
Regards,

Michele

[carlocamilloni]

@GiovanniBussi @invemichele I think this is good. While about #399 I think it still need to be double check. Anyway it would be good to have the feature improved and fully working for plumed 2.5

[GiovanniBussi]

My only doubt is that I discussed exactly this possibility (using a single grid) with @gtribello before he merged the first implementation and apparently there was some reason to use a double grid. So I would feel better if @gtribello confirms it's ok before we remove the feature.