Skip to content

uwham #

Estimate free energies using the UWHAM method [1]

References

[1]: Zhiqiang Tan, Emilio Gallicchio, Mauro Lapelosa, and Ronald M. Levy, "Theory of binless multi-state free energy estimation with applications to protein-ligand binding", J. Chem. Phys. 136, 144102 (2012)

Functions:

  • estimate_f_i –

    Estimates the free energies of a set of sampled states.

estimate_f_i #

estimate_f_i(
    u_kn: ndarray, n_k: ndarray
) -> tuple[ndarray, ndarray, ndarray]

Estimates the free energies of a set of sampled states.

Parameters:

  • u_kn (ndarray) –

    The uncorrelated reduced potentials sampled at k states with shape=(n_states, n_samples).

  • n_k (ndarray) –

    The number of samples at state k.

Returns:

  • tuple[ndarray, ndarray, ndarray] –

    The estimated reduced free energies and their estimated variance.

Source code in femto/fe/uwham.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def estimate_f_i(
    u_kn: numpy.ndarray, n_k: numpy.ndarray
) -> tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]:
    """Estimates the free energies of a set of *sampled* states.

    Args:
        u_kn: The uncorrelated reduced potentials sampled at ``k`` states with
            ``shape=(n_states, n_samples)``.
        n_k: The number of samples at state ``k``.

    Returns:
        The estimated reduced free energies and their estimated variance.
    """

    u_kn = numpy.array(u_kn)
    n_k = numpy.array(n_k)

    ln_q = -u_kn.T

    n_samples, n_states = ln_q.shape

    if n_states != len(n_k):
        raise RuntimeError("The number of states do not match")
    if n_samples != n_k.sum():
        raise RuntimeError("The number of samples do not match")

    ln_z = numpy.zeros(len(n_k) - 1)  # ln_z_0 is always fixed at 0.0
    ln_q -= ln_q[:, :1]

    n = n_k.sum()
    factor = n_k / n

    result = scipy.optimize.minimize(
        functools.partial(_compute_kappa, ln_q=ln_q, n=n, factor=factor),
        ln_z,
        method="trust-ncg",
        jac=True,
        hess=functools.partial(_compute_kappa_hessian, ln_q=ln_q, n=n, factor=factor),
    )

    if not result.success:
        raise RuntimeError("The UWHAM minimization failed to converge.")

    f_i = numpy.insert(-result.x, 0, 0.0)
    ln_z = numpy.insert(result.x, 0, 0.0)

    weights = _compute_weights(ln_z, ln_q, factor)

    if not numpy.allclose(weights.sum(axis=0) / n, 1.0, atol=1e-3):
        raise RuntimeError("The UWHAM weights do not sum to 1.0")

    df_i = _compute_variance(ln_z, weights, factor, n)

    return f_i, df_i, weights / n