Add entropy and enthalpy

CHANGES

@@ -12,11 +12,12 @@ The rules for this file:
   * accompany each entry with github issue/PR number (Issue #xyz)
   * release numbers follow "Semantic Versioning" https://semver.org
 
-**/**/**** xiki-tempula
+**/**/**** xiki-tempula, jaclark5
 
   * 2.5.0
 
 Enhancements
+  - Added matrices for entropy and enthalpy for MBAR (PR #406)
   - Parallelise read and preprocess for ABFE workflow. (PR #371)
 
 
@@ -27,7 +28,7 @@ Enhancements
 Fixes
   - [doc] tutorial: use alchemlyb.concat (PR #399)
   - Resolve pandas FutureWarnings in bar_.py and mbar_.py (issue #408 PR #407)
- 
+
 
 09/17/2024 jaclark5, orbeckst
 

src/alchemlyb/estimators/base.py

@@ -5,6 +5,10 @@ class _EstimatorMixOut:
     _d_delta_f_ = None
     _delta_f_ = None
     _states_ = None
+    _d_delta_h_ = None
+    _delta_h_ = None
+    _d_delta_s_ = None
+    _delta_s_ = None
 
     @property
     def d_delta_f_(self):
@@ -14,6 +18,22 @@ def d_delta_f_(self):
     def delta_f_(self):
         return self._delta_f_
 
+    @property
+    def d_delta_h_(self):
+        return self._d_delta_h_
+
+    @property
+    def delta_h_(self):
+        return self._delta_h_
+
+    @property
+    def d_delta_s_(self):
+        return self._d_delta_s_
+
+    @property
+    def delta_s_(self):
+        return self._delta_s_
+
     @property
     def states_(self):
         return self._states_

src/alchemlyb/estimators/mbar_.py

@@ -57,6 +57,20 @@ class MBAR(BaseEstimator, _EstimatorMixOut):
         The estimated statistical uncertainty (one standard deviation) in
         dimensionless free energy differences.
 
+    delta_h_ : DataFrame
+        The estimated dimensionless enthalpy difference between each state.
+
+    d_delta_h_ : DataFrame
+        The estimated statistical uncertainty (one standard deviation) in
+        dimensionless enthalpy differences.
+
+    delta_s_ : DataFrame, optional
+        The estimated dimensionless entropy difference between each state.
+
+    d_delta_s_ : DataFrame
+        The estimated statistical uncertainty (one standard deviation) in
+        dimensionless entropy differences.
+
     theta_ : DataFrame
         The theta matrix.
 
@@ -80,9 +94,11 @@ class MBAR(BaseEstimator, _EstimatorMixOut):
     .. versionchanged:: 2.1.0
         `n_bootstraps` option added.
     .. versionchanged:: 2.4.0
-       Handle initial estimate, initial_f_k, from bar in the instance
-       that not all lambda states represented as column headers are
-       represented in the indices of u_nk.
+        Handle initial estimate, initial_f_k, from bar in the instance
+        that not all lambda states represented as column headers are
+        represented in the indices of u_nk.
+    .. versionchanged:: 2.5.0
+        Added computation of enthalpy and entropy
 
     """
 
@@ -110,7 +126,7 @@ def __init__(
         # handle for pymbar.MBAR object
         self._mbar = None
 
-    def fit(self, u_nk):
+    def fit(self, u_nk, compute_entropy_enthalpy=False):
         """
         Compute overlap matrix of reduced potentials using multi-state
         Bennett acceptance ratio.
@@ -121,6 +137,9 @@ def fit(self, u_nk):
             ``u_nk[n, k]`` is the reduced potential energy of uncorrelated
             configuration ``n`` evaluated at state ``k``.
 
+        compute_entropy_enthalpy : bool, optional, default=False
+            Compute entropy and enthalpy.
+
         """
         # sort by state so that rows from same state are in contiguous blocks
         u_nk = u_nk.sort_index(level=u_nk.index.names[1:])
@@ -171,13 +190,18 @@ def fit(self, u_nk):
             solver_protocol=self.method,
             n_bootstraps=self.n_bootstraps,
         )
-        if self.n_bootstraps == 0:
-            uncertainty_method = None
-        else:
-            uncertainty_method = "bootstrap"
+
+        uncertainty_method = None if self.n_bootstraps == 0 else "bootstrap"
         out = self._mbar.compute_free_energy_differences(
             return_theta=True, uncertainty_method=uncertainty_method
         )
+        if compute_entropy_enthalpy:
+            out.update(
+                self._mbar.compute_entropy_and_enthalpy(
+                    uncertainty_method=uncertainty_method
+                )
+            )
+
         self._delta_f_ = pd.DataFrame(
             out["Delta_f"], columns=self._states_, index=self._states_
         )
@@ -187,9 +211,29 @@ def fit(self, u_nk):
         self.theta_ = pd.DataFrame(
             out["Theta"], columns=self._states_, index=self._states_
         )
+        if compute_entropy_enthalpy:
+            self._delta_h_ = pd.DataFrame(
+                out["Delta_u"], columns=self._states_, index=self._states_
+            )
+            self._d_delta_h_ = pd.DataFrame(
+                out["dDelta_u"], columns=self._states_, index=self._states_
+            )
+            self._delta_s_ = pd.DataFrame(
+                out["Delta_s"], columns=self._states_, index=self._states_
+            )
+            self._d_delta_s_ = pd.DataFrame(
+                out["dDelta_s"], columns=self._states_, index=self._states_
+            )
 
         self._delta_f_.attrs = u_nk.attrs
         self._d_delta_f_.attrs = u_nk.attrs
+        if compute_entropy_enthalpy:
+            self._delta_h_.attrs = u_nk.attrs
+            self._d_delta_h_.attrs = u_nk.attrs
+            self._delta_s_.attrs = u_nk.attrs
+            self._d_delta_s_.attrs = u_nk.attrs
+            self._delta_s_.attrs["energy_unit"] = "k"
+            self._d_delta_s_.attrs["energy_unit"] = "k"
 
         return self
 

src/alchemlyb/tests/test_fep_estimators.py

@@ -198,6 +198,31 @@ def test_bootstrap(gmx_benzene_Coulomb_u_nk):
     assert mbar_bootstrap_err != mbar_err
 
 
+def test_enthalpy_entropy_mbar(gmx_benzene_Coulomb_u_nk):
+    mbar = MBAR()
+    u_nk = alchemlyb.concat(gmx_benzene_Coulomb_u_nk)
+    mbar.fit(u_nk, compute_entropy_enthalpy=True)
+
+    assert mbar.delta_f_.iloc[0, :].to_numpy() == pytest.approx(
+        np.array([0.0, 1.619069, 2.557990, 2.986302, 3.041156]), abs=1e-6
+    )
+    assert mbar.delta_h_.iloc[0, :].to_numpy() == pytest.approx(
+        np.array([0.0, 1.241970, 1.695000, 1.706555, 1.388906]), abs=1e-6
+    )
+    assert mbar.delta_s_.iloc[0, :].to_numpy() == pytest.approx(
+        np.array([0.0, -0.377099, -0.862990, -1.279746, -1.652250]), abs=1e-6
+    )
+    assert mbar.d_delta_f_.iloc[0, :].to_numpy() == pytest.approx(
+        np.array([0.0, 0.008802, 0.014432, 0.018097, 0.020879]), abs=1e-6
+    )
+    assert mbar.d_delta_h_.iloc[0, :].to_numpy() == pytest.approx(
+        np.array([0.0, 0.011598, 0.016538, 0.018077, 0.018940]), abs=1e-6
+    )
+    assert mbar.d_delta_s_.iloc[0, :].to_numpy() == pytest.approx(
+        np.array([0.0, 0.012242, 0.019519, 0.023606, 0.026107]), abs=1e-6
+    )
+
+
 def test_wrong_initial_f_k():
     with pytest.raises(
         ValueError, match="Only `BAR` is supported as string input to `initial_f_k`"