From 7b5e5bd5a96a1ee9cd8311e5a92e19300dcc086b Mon Sep 17 00:00:00 2001 From: jac16 Date: Wed, 18 Sep 2024 11:15:55 -0400 Subject: [PATCH 1/3] Add computation of entropic/enthalpy contributions from BAR --- pymbar/__init__.py | 2 +- pymbar/other_estimators.py | 146 ++++++++++++++++++++++++++++++++++++- pymbar/tests/test_bar.py | 47 ++++++++++++ 3 files changed, 192 insertions(+), 3 deletions(-) diff --git a/pymbar/__init__.py b/pymbar/__init__.py index d692f01b..4bf1ad86 100644 --- a/pymbar/__init__.py +++ b/pymbar/__init__.py @@ -30,7 +30,7 @@ from . import timeseries, testsystems, confidenceintervals from .mbar import MBAR -from .other_estimators import bar, bar_overlap, bar_zero, exp, exp_gauss +from .other_estimators import bar, bar_overlap, bar_enthalpy_entropy, bar_zero, exp, exp_gauss from .fes import FES diff --git a/pymbar/other_estimators.py b/pymbar/other_estimators.py index d9539df0..5fd54820 100644 --- a/pymbar/other_estimators.py +++ b/pymbar/other_estimators.py @@ -153,6 +153,148 @@ def bar_zero(w_F, w_R, DeltaF): return fzero +def bar_enthalpy_entropy(u11, u12, u22, u21, DeltaF, dDeltaF=None): + """Compute the enthalpy and entropy components with BAR + + This function calculates the enthalpy difference between states using a + derivation of the entropy from the Bennett Acceptance Ratio (BAR) method + at constant temperature. + + from DOI: 10.1021/jp103050u + + M = log( N_F / N_R ) + + g_1(x) = 1 / [1 + exp(M + w_F(x) - Delta F)] + g_2(x) = 1 / [1 + exp(-M - w_R(x) + Delta F)] + + a_F = _F - _F * _F + _0 + a_R = _R - _R * _R - _1 + + we want: + + H = (N_F * a_F - N_R * a_R) / (N_F * _F + N_R * _R) + + The entropy is then taken as the difference between the enthalpy and free energy. + + Note that: + w_F = u_kln[0, 1] - u_kln[0, 0] + w_R = u_kln[1, 0] - u_kln[1, 1] + + Parameters + ---------- + u11 : np.ndarray + Reduced potential energy of state 1 in a configuration + sampled from state 1 in snapshot t. + t = 0...(T_1-1) Length T_1 is deduced from vector. + u12 : np.ndarray + Reduced potential energy of state 2 in a configuration + sampled from state 1 in snapshot t. + t = 0...(T_1-1) Length must be equal to T_1. + u22 : np.ndarray + Reduced potential energy of state 2 in a configuration + sampled from state 2 in snapshot t. + t = 0...(T_2-1) Length T_2 is deduced from vector. + u21 : np.ndarray + Reduced potential energy of state 1 in a configuration + sampled from state 2 in snapshot t. + t = 0...(T_2-1) Length must be equal to T_2. + Delta_f : float + Free energy difference + dDelta_f : float, default=None + Estimated standard deviation of free energy difference + + Returns + ------- + dict + 'Delta_f' : float + Free energy difference + 'dDelta_f' : float + Estimated standard error of free energy difference + 'Delta_h' : float + Enthalpy difference + 'dDelta_h' : float + Estimated standard error of enthalpy difference + 'Delta_s' : float + Enthalpy difference + 'dDelta_s' : float + Estimated standard error of enthalpy difference if + the standard deviation of the free energy, ``dDelta_f``, + is provided. + + Examples + -------- + Compute free energy difference between two specified samples of work values. + + >>> from pymbar import testsystems + >>> [u11, u12, u22, u21] = ??? # testsystems.gaussian_work_example(mu_F=None, DeltaF=1.0, seed=0) + >>> DeltaH = bar_enthalpy(u11, u12, u22, u21, DeltaF) + + """ + + np.seterr(over="raise") # raise exceptions to overflows + results = {'Delta_f': DeltaF} + if dDeltaF is not None: + results["dDelta_f"] = dDeltaF + + u11, u12 = np.array(u11, np.float64), np.array(u12, np.float64) + u22, u21 = np.array(u22, np.float64), np.array(u21, np.float64) + DeltaF = float(DeltaF) + + # Recommended stable implementation of bar. + + # Determine number of forward and reverse work values provided. + T_1 = float(u11.size) # number of forward work values + T_2 = float(u22.size) # number of reverse work values + if len(u12) != T_1: + raise ValueError("The length of u12 must be equal to the length of u11.") + if len(u21) != T_2: + raise ValueError("The length of u21 must be equal to the length of u22.") + + # Compute log ratio of forward and reverse counts. + M = np.log(T_1 / T_2) + + g_A1 = 1 / (1 + np.exp(u12 - u11 + M - DeltaF)) + g_A2 = 1 / (1 + np.exp(u22 - u21 + M - DeltaF)) + g_B1 = 1 / (1 + np.exp(-u12 + u11 - M + DeltaF)) + g_B2 = 1 / (1 + np.exp(-u22 + u21 - M + DeltaF)) + + a_1 = np.mean(g_A1 * u11) - np.mean(g_A1) * np.mean(u11) + np.mean(g_A1 * g_B1 * (u12 - u11)) + a_2 = np.mean(g_B2 * u22) - np.mean(g_B2) * np.mean(u22) - np.mean(g_A2 * g_B2 * (u22 - u21)) + + tmp1, tmp2 = np.mean(g_A1 * g_B1), np.mean(g_A2 * g_B2) + results["Delta_h"] = (T_1 * a_1 - T_2 * a_2) / (T_1 * tmp1 + T_2 * tmp2) + + # Calculate the uncertainty as standard errors + da_1 = np.sqrt( + np.std(g_A1 * u11)**2 + + np.sqrt((np.std(g_A1) * np.mean(u11))**2 + (np.mean(g_A1) * np.std(u11))**2) + + np.std(g_A1 * g_B1 * (u12 - u11))**2 + ) / np.sqrt(T_1) + da_2 = np.sqrt( + np.std(g_B2 * u22)**2 + + np.sqrt((np.std(g_B2) * np.mean(u22))**2 + (np.mean(g_B2) * np.std(u22))**2) + + np.std(g_A2 * g_B2 * (u22 - u21))**2 + ) / np.sqrt(T_2) + dtmp1, dtmp2 = np.std(g_A1 * g_B1) / np.sqrt(T_1), np.std(g_A2 * g_B2) / np.sqrt(T_2) + + dHda1 = T_1 / (T_1 * tmp1 + T_2 * tmp2) + dHda2 = -T_2 / (T_1 * tmp1 + T_2 * tmp2) + dHdtmp1 = T_1 * (T_2 * a_2 - T_1 * a_1) / (T_1 * tmp1 + T_2 * tmp2)**2 + dHdtmp2 = T_2 * (T_2 * a_2 - T_1 * a_1) / (T_1 * tmp1 + T_2 * tmp2)**2 + + results["dDelta_h"] = np.sqrt( + (dHda1 * da_1)**2 + + (dHda2 * da_2)**2 + + (dHdtmp1 * dtmp1)**2 + + (dHdtmp2 * dtmp2)**2 + ) + + results["Delta_s"] = results["Delta_h"] - results["Delta_f"] + if 'dDelta_f' in results: + results["dDelta_s"] = np.sqrt(results["dDelta_h"]**2 + results["dDelta_f"]**2) + + return results + def bar( w_F, w_R, @@ -189,7 +331,7 @@ def bar( relative_tolerance : float, optional, default=1E-12 Can be set to determine the relative tolerance convergence criteria (default 1.0e-12) verbose : bool - Should be set to True if verbse debug output is desired (default False) + Should be set to True if verbose debug output is desired (default False) method: str, optional, default='false-position' Choice of method to solve bar nonlinear equations: one of 'bisection', 'self-consistent-iteration' or 'false-position' (default : 'false-position'). iterated_solution: bool, optional, default=True @@ -218,7 +360,7 @@ def bar( -------- Compute free energy difference between two specified samples of work values. - >>> from pymbar import testsystems + >>> from pymbar import testsystems, bar >>> [w_F, w_R] = testsystems.gaussian_work_example(mu_F=None, DeltaF=1.0, seed=0) >>> results = bar(w_F, w_R) >>> print('Free energy difference is {:.3f} +- {:.3f} kT'.format(results['Delta_f'], results['dDelta_f'])) diff --git a/pymbar/tests/test_bar.py b/pymbar/tests/test_bar.py index 29ffda5c..c8090ba8 100644 --- a/pymbar/tests/test_bar.py +++ b/pymbar/tests/test_bar.py @@ -98,6 +98,53 @@ def test_bar_free_energies(bar_and_test): assert_almost_equal(dfe_bar, dfe_mbar, decimal=3) +@pytest.mark.parametrize("system_generator", system_generators) +def test_bar_enthalpy_entropy(system_generator): + """Test Bar estimate of enthalpy and entropy in harmonic oscillator""" + + name, test = system_generator() + + Nk = int(1e+5) + _, u_kln, _ = test.sample(mode="u_kln", N_k=(Nk, Nk)) + w_F = u_kln[0, 1] - u_kln[0, 0] + w_R = u_kln[1, 0] - u_kln[1, 1] + + _result = estimators.bar(w_F, w_R) + results_bar = estimators.bar_enthalpy_entropy( + u_kln[0,0], + u_kln[0,1], + u_kln[1,1], + u_kln[1,0], + _result['Delta_f'], + dDeltaF=_result['dDelta_f'], + ) + + G = test.analytical_free_energies() + H = test.analytical_observable(observable="potential energy") + S = test.analytical_entropies() + results_analytical = { + "Delta_f": G[1]-G[0], + "Delta_h": H[1]-H[0], + "Delta_s": S[1]-S[0], + } + + assert_almost_equal( + results_bar["Delta_f"], results_analytical["Delta_f"], decimal=2 + ) + assert_almost_equal( + results_bar["Delta_h"], results_analytical["Delta_h"], decimal=2 + ) + assert_almost_equal( + results_bar["Delta_s"], results_analytical["Delta_s"], decimal=2 + ) + assert_almost_equal( + results_bar["dDelta_f"], 0.0, decimal=2 + ) + assert_almost_equal( + results_bar["dDelta_h"], 0.0, decimal=2 + ) + + def test_bar_overlap(): for system_generator in system_generators: name, test = system_generator() From b878c0cc2bd6e1e2bf54fb32e658a0d40f2424fe Mon Sep 17 00:00:00 2001 From: jac16 Date: Wed, 18 Sep 2024 11:17:57 -0400 Subject: [PATCH 2/3] Remove tracking of project settings --- .gitignore | 3 +++ .vscode/settings.json | 5 ----- 2 files changed, 3 insertions(+), 5 deletions(-) delete mode 100644 .vscode/settings.json diff --git a/.gitignore b/.gitignore index 6b353eff..cb5a29e5 100644 --- a/.gitignore +++ b/.gitignore @@ -118,6 +118,9 @@ venv.bak/ .spyderproject .spyproject +# VSCode project settings +.vscode/ + # Rope project settings .ropeproject diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index 36fe635c..00000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,5 +0,0 @@ -{ - "python.linting.enabled": true, - "python.linting.pylintEnabled": true, - "python.pythonPath": "/home/jaime/.conda/envs/pymbar/bin/python" -} \ No newline at end of file From 8ffe4c991b1448b569341e23957a4ab8745f3a09 Mon Sep 17 00:00:00 2001 From: jac16 Date: Wed, 18 Sep 2024 11:18:35 -0400 Subject: [PATCH 3/3] Fix docstring for HarmonicOscillator --- pymbar/testsystems/harmonic_oscillators.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pymbar/testsystems/harmonic_oscillators.py b/pymbar/testsystems/harmonic_oscillators.py index 32ac4cd6..cf19d1ad 100644 --- a/pymbar/testsystems/harmonic_oscillators.py +++ b/pymbar/testsystems/harmonic_oscillators.py @@ -10,7 +10,7 @@ class HarmonicOscillatorsTestCase(object): Generate energy samples with default parameters. >>> testcase = HarmonicOscillatorsTestCase() - >>> x_kn, u_kln, N_k, s_n = testcase.sample() + >>> x_n, u_kn, N_k, s_n = testcase.sample() Retrieve analytical properties. @@ -22,12 +22,12 @@ class HarmonicOscillatorsTestCase(object): Generate energy samples with default parameters in one line. - >>> x_kn, u_kln, N_k, s_n = HarmonicOscillatorsTestCase().sample() + >>> x_n, u_kn, N_k, s_n = HarmonicOscillatorsTestCase().sample() Generate energy samples with specified parameters. >>> testcase = HarmonicOscillatorsTestCase(O_k=[0, 1, 2, 3, 4], K_k=[1, 2, 4, 8, 16]) - >>> x_kn, u_kln, N_k, s_n = testcase.sample(N_k=[10, 20, 30, 40, 50]) + >>> x_n, u_kn, N_k, s_n = testcase.sample(N_k=[10, 20, 30, 40, 50]) Test sampling in different output modes. @@ -131,9 +131,11 @@ def sample(self, N_k=(10, 20, 30, 40, 50), mode="u_kn", seed=None): s_n : np.ndarray, shape=(n_samples), dtype='int' s_n is the state of origin of x_n + if mode == 'u_kln': + x_kn : np.ndarray, shape=(n_states, n_samples), dtype=float 1D harmonic oscillator positions - u_kln : np.ndarray, shape=(n_states, n_states, n_samples), dytpe=float, only if mode='u_kln' + u_kln : np.ndarray, shape=(n_states, n_states, n_samples), dytpe=float u_kln[k,l,n] is reduced potential of sample n from state k evaluated at state l. N_k : np.ndarray, shape=(n_states), dtype=int32 N_k[k] is the number of samples generated from state k