Update on 03/14/26 at 19:08:21

Author: gnshb

polyprotic_pr_review.md

@@ -557,7 +557,10 @@ def formal_charge(particle_name):
                         terms = [10.0 ** (-cumsum_pK[i] + (i + 1) * pH) for i in range(n)]
                         denominator = 1.0 + sum(terms)
                         numerator = sum((i + 1) * terms[i] for i in range(n))
-                        charge = psi * numerator / denominator
+                        if acidity == "acidic":
+                            charge = -numerator / denominator
+                        else:
+                            charge = n - numerator / denominator
                     else:
                         pka = entry["pka_value"]
                         charge = psi / (1.0 + 10.0 ** (psi * (pH - pka)))

pyMBE/pyMBE.py

@@ -0,0 +1,137 @@
+"""
+Compare pyMBE's polyprotic Henderson-Hasselbalch implementation against:
+  1. The analytical calculate_Z formula from pbs_simulations
+  2. Ideal constant-pH simulation data from pbs_simulations sample_data
+
+Produces comparison plots for monoprotic, diprotic, and triprotic acids.
+"""
+import sys
+import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+sys.path.insert(0, os.path.join(os.path.dirname(__file__), ".."))
+import pyMBE
+
+# Add pbs_simulations analysis to path
+PBS_ROOT = os.path.expanduser("~/ntnu/pbs_simulations")
+sys.path.insert(0, os.path.join(PBS_ROOT, "scripts"))
+from analysis import analyze_time_series
+
+# ---------- Reference analytical formula from pbs_simulations ----------
+
+def calculate_Z_reference(pK, pH_array):
+    """Exact polyprotic HH formula from pbs_simulations/scripts/post_processing.py"""
+    pK = np.asarray(pK)
+    pH_array = np.atleast_1d(pH_array).astype(float)
+    y = pH_array[None, :]  # (1, N_pH)
+    i_range = np.arange(1, len(pK) + 1)[:, None]  # (n, 1)
+    cumsum_pK = np.cumsum(pK)[:, None]  # (n, 1)
+    fractions = 10.0 ** (-cumsum_pK + i_range * y)
+    numerator = np.sum(i_range * fractions, axis=0)
+    denominator = 1.0 + np.sum(fractions, axis=0)
+    return -numerator / denominator
+
+
+def compute_sim_Z(sample_dir, n):
+    """Compute average charge Z from simulation sample data using block binning."""
+    analyzed = analyze_time_series(path_to_datafolder=sample_dir,
+                                   ignore_files=["analyzed_data.csv"])
+    pH_vals = []
+    Z_vals = []
+    Z_err = []
+    for idx in range(len(analyzed)):
+        row = analyzed.iloc[idx]
+        pH_val = float(row[("pH", "value")])
+        n_ha = float(row[("mean", "n_ha")])
+        Z_num = -sum(float(row[("mean", f"n_a{i+1}")]) * (i + 1) for i in range(n))
+        Z_den = n_ha + sum(float(row[("mean", f"n_a{i+1}")]) for i in range(n))
+        pH_vals.append(pH_val)
+        Z_vals.append(Z_num / Z_den if Z_den > 0 else 0.0)
+    order = np.argsort(pH_vals)
+    return np.array(pH_vals)[order], np.array(Z_vals)[order]
+
+
+def compute_pyMBE_Z(pka_list, pH_array):
+    """Compute charge using pyMBE's calculate_HH."""
+    pmb = pyMBE.pymbe_library(seed=42)
+    n = len(pka_list)
+    if n == 1:
+        pmb.define_particle(name="acid",
+                            sigma=0.35 * pmb.units.nm,
+                            epsilon=1 * pmb.units("reduced_energy"),
+                            acidity="acidic",
+                            pka=pka_list[0])
+    else:
+        pmb.define_polyprotic_particle(name="acid",
+                                       sigma=0.35 * pmb.units.nm,
+                                       epsilon=1 * pmb.units("reduced_energy"),
+                                       n=n,
+                                       acidity="acidic",
+                                       pka_list=pka_list)
+    pmb.define_residue(name="res", central_bead="acid", side_chains=[])
+    pmb.define_molecule(name="mol", residue_list=["res"])
+    return pmb.calculate_HH(template_name="mol", pH_list=list(pH_array))
+
+
+# ---------- Configuration ----------
+
+pK_values = {
+    "monoprotic": [2.16],
+    "diprotic": [2.16, 7.21],
+    "triprotic": [2.16, 7.21, 12.32],
+}
+
+sample_dirs = {
+    "monoprotic": os.path.join(PBS_ROOT, "sample_data", "ideal-monoprotic"),
+    "diprotic": os.path.join(PBS_ROOT, "sample_data", "ideal-diprotic"),
+    "triprotic": os.path.join(PBS_ROOT, "sample_data", "ideal-triprotic"),
+}
+
+# ---------- Plot ----------
+
+fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharey=True)
+
+for ax, label in zip(axes, ["monoprotic", "diprotic", "triprotic"]):
+    pK = pK_values[label]
+    n = len(pK)
+    pH_fine = np.linspace(0.5, 14.5, 200)
+
+    # 1. Reference analytical formula
+    Z_ref = calculate_Z_reference(pK, pH_fine)
+
+    # 2. pyMBE calculate_HH
+    Z_pyMBE = compute_pyMBE_Z(pK, pH_fine)
+
+    # 3. Simulation data
+    pH_sim, Z_sim = compute_sim_Z(sample_dirs[label], n)
+
+    # Plot
+    ax.plot(pH_fine, Z_ref, "k-", lw=2, label="Analytical (pbs_simulations)")
+    ax.plot(pH_fine, Z_pyMBE, "r--", lw=2, label="pyMBE calculate_HH")
+    ax.scatter(pH_sim, Z_sim, c="blue", s=30, zorder=5, label="Simulation data (ideal)")
+
+    # pKa markers
+    for i, pk in enumerate(pK):
+        ax.axvline(x=pk, color="gray", ls=":", alpha=0.6)
+        ax.text(pk + 0.2, -n + 0.3, f"pK{i+1}={pk}", fontsize=8, color="gray")
+
+    ax.set_xlabel("pH")
+    ax.set_title(f"{n}-protic acid")
+    ax.grid(alpha=0.3)
+
+    # MSE between pyMBE and reference
+    mse = np.mean((np.array(Z_pyMBE) - Z_ref) ** 2)
+    ax.text(0.05, 0.05, f"MSE(pyMBE vs analytical) = {mse:.2e}",
+            transform=ax.transAxes, fontsize=8,
+            bbox=dict(boxstyle="round", facecolor="wheat", alpha=0.5))
+
+axes[0].set_ylabel("Average charge Z")
+axes[0].legend(loc="lower left", fontsize=8)
+
+plt.suptitle("Polyprotic HH: pyMBE vs Analytical vs Simulation", fontsize=14)
+plt.tight_layout()
+plt.savefig("tests/polyprotic_HH_comparison.png", dpi=150, bbox_inches="tight")
+print("Plot saved to tests/polyprotic_HH_comparison.png")
+plt.show()

tests/compare_polyprotic_HH.py

@@ -0,0 +1,93 @@
+"""
+Compare pyMBE's polyprotic HH for basic particles against the analytical formula.
+No simulation data available for bases — analytical comparison only.
+"""
+import sys
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+sys.path.insert(0, os.path.join(os.path.dirname(__file__), ".."))
+import pyMBE
+
+
+def calculate_Z_base_reference(pK, pH_array):
+    """
+    Analytical polyprotic HH for a basic particle.
+    Fully protonated has charge +n, each deprotonation removes +1.
+    Z = n - sum(i * 10^(-cumsum(pK)[i] + i*pH)) / (1 + sum(10^(-cumsum(pK)[i] + i*pH)))
+    """
+    pK = np.asarray(pK)
+    n = len(pK)
+    pH_array = np.atleast_1d(pH_array).astype(float)
+    y = pH_array[None, :]
+    i_range = np.arange(1, n + 1)[:, None]
+    cumsum_pK = np.cumsum(pK)[:, None]
+    fractions = 10.0 ** (-cumsum_pK + i_range * y)
+    numerator = np.sum(i_range * fractions, axis=0)
+    denominator = 1.0 + np.sum(fractions, axis=0)
+    return n - numerator / denominator
+
+
+def compute_pyMBE_Z_base(pka_list, pH_array):
+    """Compute charge using pyMBE's calculate_HH for a basic particle."""
+    pmb = pyMBE.pymbe_library(seed=42)
+    n = len(pka_list)
+    if n == 1:
+        pmb.define_particle(name="base",
+                            sigma=0.35 * pmb.units.nm,
+                            epsilon=1 * pmb.units("reduced_energy"),
+                            acidity="basic",
+                            pka=pka_list[0])
+    else:
+        pmb.define_polyprotic_particle(name="base",
+                                       sigma=0.35 * pmb.units.nm,
+                                       epsilon=1 * pmb.units("reduced_energy"),
+                                       n=n,
+                                       acidity="basic",
+                                       pka_list=pka_list)
+    pmb.define_residue(name="res", central_bead="base", side_chains=[])
+    pmb.define_molecule(name="mol", residue_list=["res"])
+    return pmb.calculate_HH(template_name="mol", pH_list=list(pH_array))
+
+
+pK_values = {
+    "monoprotic": [9.25],
+    "diprotic": [6.0, 10.0],
+    "triprotic": [4.0, 8.0, 12.0],
+}
+
+fig, axes = plt.subplots(1, 3, figsize=(15, 5), sharey=True)
+pH_fine = np.linspace(0.5, 14.5, 200)
+
+for ax, label in zip(axes, ["monoprotic", "diprotic", "triprotic"]):
+    pK = pK_values[label]
+    n = len(pK)
+
+    Z_ref = calculate_Z_base_reference(pK, pH_fine)
+    Z_pyMBE = compute_pyMBE_Z_base(pK, pH_fine)
+
+    ax.plot(pH_fine, Z_ref, "k-", lw=2, label="Analytical")
+    ax.plot(pH_fine, Z_pyMBE, "r--", lw=2, label="pyMBE calculate_HH")
+
+    for i, pk in enumerate(pK):
+        ax.axvline(x=pk, color="gray", ls=":", alpha=0.6)
+        ax.text(pk + 0.2, 0.3, f"pK{i+1}={pk}", fontsize=8, color="gray")
+
+    ax.set_xlabel("pH")
+    ax.set_title(f"{n}-protic base")
+    ax.grid(alpha=0.3)
+
+    mse = np.mean((np.array(Z_pyMBE) - Z_ref) ** 2)
+    ax.text(0.05, 0.05, f"MSE = {mse:.2e}",
+            transform=ax.transAxes, fontsize=8,
+            bbox=dict(boxstyle="round", facecolor="wheat", alpha=0.5))
+
+axes[0].set_ylabel("Average charge Z")
+axes[0].legend(loc="upper right", fontsize=8)
+
+plt.suptitle("Polyprotic HH (basic): pyMBE vs Analytical", fontsize=14)
+plt.tight_layout()
+plt.savefig("tests/polyprotic_HH_comparison_base.png", dpi=150, bbox_inches="tight")
+print("Plot saved to tests/polyprotic_HH_comparison_base.png")
+plt.show()

tests/compare_polyprotic_HH_base.py

@@ -300,5 +300,45 @@ def test_auto_naming_convention(self):
         self.assertEqual(names, ["H4Y", "H3Y", "H2Y", "HY", "Y"])
 
 
+    def test_calculate_HH_triprotic_acid(self):
+        """
+        Test Henderson-Hasselbalch charge for a triprotic acid at extreme pH values.
+        """
+        pmb = pyMBE.pymbe_library(seed=42)
+        pmb.define_polyprotic_particle(name="PO4",
+                                       sigma=0.35*pmb.units.nm,
+                                       epsilon=1*pmb.units("reduced_energy"),
+                                       n=3,
+                                       acidity="acidic",
+                                       pka_list=[2.15, 7.20, 12.35])
+        pmb.define_molecule(name="mol", residue_list=[])
+        pmb.define_residue(name="res", central_bead="PO4", side_chains=[])
+        pmb.define_molecule(name="test_mol", residue_list=["res"])
+        Z = pmb.calculate_HH(template_name="test_mol", pH_list=[0, 14])
+        # At pH 0: fully protonated, charge ≈ 0
+        self.assertAlmostEqual(Z[0], 0.0, places=1)
+        # At pH 14: fully deprotonated, charge ≈ -3
+        self.assertAlmostEqual(Z[1], -3.0, places=1)
+
+    def test_calculate_HH_diprotic_base(self):
+        """
+        Test Henderson-Hasselbalch charge for a diprotic base at extreme pH values.
+        """
+        pmb = pyMBE.pymbe_library(seed=42)
+        pmb.define_polyprotic_particle(name="B",
+                                       sigma=0.35*pmb.units.nm,
+                                       epsilon=1*pmb.units("reduced_energy"),
+                                       n=2,
+                                       acidity="basic",
+                                       pka_list=[6.0, 10.0])
+        pmb.define_residue(name="res", central_bead="B", side_chains=[])
+        pmb.define_molecule(name="test_mol", residue_list=["res"])
+        Z = pmb.calculate_HH(template_name="test_mol", pH_list=[0, 14])
+        # At pH 0: fully protonated, charge ≈ +2
+        self.assertAlmostEqual(Z[0], 2.0, places=2)
+        # At pH 14: fully deprotonated, charge ≈ 0
+        self.assertAlmostEqual(Z[1], 0.0, places=2)
+
+
 if __name__ == "__main__":
     ut.main()