choderalab / pymbar

Python implementation of the multistate Bennett acceptance ratio (MBAR)
http://pymbar.readthedocs.io
MIT License
240 stars 93 forks source link

`hybr` solver fails, then seeds `adaptive` solver with bad starting `f_i` #544

Open SimonBoothroyd opened 1 week ago

SimonBoothroyd commented 1 week ago

I've started to see some strange cases where MBAR (in pymbar 4.0) fails to converge (i.e. both the hybr and adaptive solvers in succesion), while UWHAM seems to converge ok and indicates reasonable overlap between states. In such cases, it seems that pymbar <4.0 also runs to convergence.

I have attached a set of u_kn and n_k from the gas phase leg of a hydration free energy calculation that exhibits this issue.

repro.zip

The issue can then be observed by running:

import numpy
import pymbar

print("PYMBAR", pymbar.__version__, flush=True)

u_kn = numpy.load("u_kn.npy")
n_k = numpy.load("n_k.npy")

mbar = pymbar.MBAR(u_kn, n_k)
print(mbar.W_nk)

it looks like the issue in 4.x is due to the default hybr solver basically failing to find a solution and returning non-sensical f_k_nonzero, which then seeds the fallback adaptive solver with with a very wrong start which then propagates through. If I instead just force MBAR to start with the adaptive solver things seem to converge ok.

In general, it would be good to throw optionally throw an error if the solvers fail to converge, or at least have an attribute on MBAR that can be queried for if the solver succeeded.