Add conformer-ensemble thermochemistry helpers

ThermoScreening/thermo/__init__.py

@@ -9,3 +9,4 @@
 from .screening import screen
 from .conformers import generate as generate_conformers, write_conformers
 from .reactions import reaction_free_energy, reduction_potential
+from .ensemble import boltzmann_weights, ensemble_free_energy, lowest_gibbs

ThermoScreening/thermo/_units.py

@@ -0,0 +1,35 @@
+"""Shared Hartree-based unit conversions for the post-processing helpers."""
+
+from ..utils.physicalConstants import PhysicalConstants
+
+HARTREE_TO_EV = PhysicalConstants["H"] / PhysicalConstants["eV"]
+HARTREE_TO_KCAL_PER_MOL = (
+    PhysicalConstants["H"] * PhysicalConstants["N_A"] / (PhysicalConstants["cal"] * 1000.0)
+)
+HARTREE_TO_KJ_PER_MOL = PhysicalConstants["H"] * PhysicalConstants["N_A"] / 1000.0
+
+# Output units for energies given internally in Hartree (per particle).
+# "eV" and "H" are per particle; "kcal" and "kJ" are per mole.
+UNIT_FACTORS = {
+    "H": 1.0,
+    "eV": HARTREE_TO_EV,
+    "kcal": HARTREE_TO_KCAL_PER_MOL,
+    "kJ": HARTREE_TO_KJ_PER_MOL,
+}
+
+
+def convert_from_hartree(value, unit):
+    """
+    Convert an energy in Hartree to ``unit`` (one of ``H``/``eV``/``kcal``/``kJ``).
+
+    Raises
+    ------
+    ValueError
+        If ``unit`` is not supported.
+    """
+    try:
+        factor = UNIT_FACTORS[unit]
+    except KeyError:
+        known = ", ".join(UNIT_FACTORS)
+        raise ValueError(f"Unknown unit {unit!r}; choose one of: {known}.")
+    return value * factor

ThermoScreening/thermo/ensemble.py

@@ -0,0 +1,119 @@
+"""Conformer-ensemble thermochemistry from computed ``Thermo`` objects.
+
+Pure post-processing helpers that combine the absolute Gibbs free energies
+(``Thermo.total_EeGtot()`` = electronic + thermal, in Hartree) of several
+conformers of the *same* molecule into ensemble properties: Boltzmann
+populations, the Boltzmann-averaged (ensemble) free energy, and the
+lowest-free-energy conformer. Compute each conformer with the same engine and
+conditions (see :func:`ThermoScreening.thermo.conformers.generate`), then
+combine them here.
+"""
+
+import math
+
+from ..utils.physicalConstants import PhysicalConstants
+from ._units import convert_from_hartree
+
+# Boltzmann constant in Hartree per kelvin (k_B [J/K] converted to Hartree).
+_BOLTZMANN_HARTREE_PER_K = PhysicalConstants["kB"] / PhysicalConstants["H"]
+
+
+def _energies_and_kt(thermos, temperature):
+    """Return the list of conformer free energies (Hartree) and k_B*T (Hartree)."""
+    if not thermos:
+        raise ValueError("thermos must contain at least one conformer.")
+    if temperature <= 0:
+        raise ValueError("temperature must be positive.")
+    energies = [thermo.total_EeGtot() for thermo in thermos]
+    return energies, _BOLTZMANN_HARTREE_PER_K * temperature
+
+
+def boltzmann_weights(thermos, temperature=298.15):
+    """
+    Boltzmann populations of a conformer ensemble at ``temperature`` (K).
+
+    Parameters
+    ----------
+    thermos : iterable of Thermo
+        Conformers of the same molecule, computed with matching settings.
+    temperature : float
+        Temperature in kelvin. Default 298.15.
+
+    Returns
+    -------
+    list of float
+        Normalised populations (summing to 1) in the order of ``thermos``.
+
+    Raises
+    ------
+    ValueError
+        If ``thermos`` is empty or ``temperature`` is not positive.
+    """
+    thermos = list(thermos)
+    energies, kt = _energies_and_kt(thermos, temperature)
+    # Shift by the minimum for numerical stability (the ratios are unchanged).
+    reference = min(energies)
+    boltzmann = [math.exp(-(energy - reference) / kt) for energy in energies]
+    partition = sum(boltzmann)
+    return [weight / partition for weight in boltzmann]
+
+
+def ensemble_free_energy(thermos, temperature=298.15, unit="H"):
+    """
+    Ensemble (Boltzmann) free energy of a conformer set.
+
+    ``G_ensemble = -k_B T ln( sum_i exp(-G_i / k_B T) )``, which lies at or below
+    the lowest conformer free energy by the mixing (conformational) entropy.
+
+    Parameters
+    ----------
+    thermos : iterable of Thermo
+        Conformers of the same molecule, computed with matching settings.
+    temperature : float
+        Temperature in kelvin. Default 298.15.
+    unit : str
+        Output unit: ``"H"`` (Hartree per particle, default), ``"eV"`` (per
+        particle), ``"kcal"`` (kcal/mol), or ``"kJ"`` (kJ/mol).
+
+    Returns
+    -------
+    float
+        The ensemble free energy in the requested unit.
+
+    Raises
+    ------
+    ValueError
+        If ``thermos`` is empty, ``temperature`` is not positive, or ``unit`` is
+        not supported.
+    """
+    thermos = list(thermos)
+    energies, kt = _energies_and_kt(thermos, temperature)
+    reference = min(energies)
+    partition = sum(math.exp(-(energy - reference) / kt) for energy in energies)
+    free_energy = reference - kt * math.log(partition)
+    return convert_from_hartree(free_energy, unit)
+
+
+def lowest_gibbs(thermos):
+    """
+    Return the conformer with the lowest absolute Gibbs free energy.
+
+    Parameters
+    ----------
+    thermos : iterable of Thermo
+        Conformers of the same molecule, computed with matching settings.
+
+    Returns
+    -------
+    Thermo
+        The lowest-free-energy conformer.
+
+    Raises
+    ------
+    ValueError
+        If ``thermos`` is empty.
+    """
+    thermos = list(thermos)
+    if not thermos:
+        raise ValueError("thermos must contain at least one conformer.")
+    return min(thermos, key=lambda thermo: thermo.total_EeGtot())

ThermoScreening/thermo/reactions.py

@@ -8,26 +8,13 @@
 appropriate -- then combine them here.
 """
 
-from ..utils.physicalConstants import PhysicalConstants
+from ._units import HARTREE_TO_EV as _HARTREE_TO_EV, convert_from_hartree
 
 # Absolute potential of the standard hydrogen electrode (V). Convention-dependent
 # (values from ~4.28 to ~4.44 V are used); this is the IUPAC-recommended value.
 # Override ``reference_potential`` for a calibrated reference.
 SHE_ABSOLUTE_POTENTIAL = 4.44
 
-_HARTREE_TO_EV = PhysicalConstants["H"] / PhysicalConstants["eV"]
-_HARTREE_TO_KCAL_PER_MOL = (
-    PhysicalConstants["H"] * PhysicalConstants["N_A"] / (PhysicalConstants["cal"] * 1000.0)
-)
-_HARTREE_TO_KJ_PER_MOL = PhysicalConstants["H"] * PhysicalConstants["N_A"] / 1000.0
-
-_UNIT_FACTORS = {
-    "H": 1.0,
-    "eV": _HARTREE_TO_EV,
-    "kcal": _HARTREE_TO_KCAL_PER_MOL,
-    "kJ": _HARTREE_TO_KJ_PER_MOL,
-}
-
 
 def _total_gibbs(species):
     """
@@ -67,17 +54,11 @@ def reaction_free_energy(reactants, products, unit="kcal"):
     ValueError
         If ``unit`` is not supported.
     """
-    try:
-        factor = _UNIT_FACTORS[unit]
-    except KeyError:
-        known = ", ".join(_UNIT_FACTORS)
-        raise ValueError(f"Unknown unit {unit!r}; choose one of: {known}.")
-
     delta = (
         sum(_total_gibbs(product) for product in products)
         - sum(_total_gibbs(reactant) for reactant in reactants)
     )
-    return delta * factor
+    return convert_from_hartree(delta, unit)
 
 
 def reduction_potential(

docs/api.rst

@@ -29,6 +29,13 @@ Reactions and redox
 .. autofunction:: ThermoScreening.thermo.reactions.reaction_free_energy
 .. autofunction:: ThermoScreening.thermo.reactions.reduction_potential
 
+Conformer ensembles
+-------------------
+
+.. autofunction:: ThermoScreening.thermo.ensemble.boltzmann_weights
+.. autofunction:: ThermoScreening.thermo.ensemble.ensemble_free_energy
+.. autofunction:: ThermoScreening.thermo.ensemble.lowest_gibbs
+
 Thermochemistry core
 --------------------
 

docs/usage.rst

@@ -72,6 +72,32 @@ Conformers, then screen the ensemble:
    write_conformers(conformers, "butanol_confs")
    results = screen("butanol_confs", engine="xtb-cli")
 
+Conformer ensembles
+-------------------
+
+A flexible molecule is a Boltzmann ensemble of conformers, not a single
+structure. Generate and compute each conformer, then combine their free energies
+into ensemble properties:
+
+.. code-block:: python
+
+   from ThermoScreening.thermo.api import xtb_cli_thermo
+   from ThermoScreening.thermo.conformers import generate
+   from ThermoScreening.thermo import (
+       boltzmann_weights, ensemble_free_energy, lowest_gibbs,
+   )
+
+   conformers = generate("CCCCO", max_conformers=10)          # n-butanol
+   thermos = [xtb_cli_thermo(c, solvent="water") for c in conformers]
+
+   weights = boltzmann_weights(thermos)          # populations at 298.15 K
+   G = ensemble_free_energy(thermos, unit="kcal")  # Boltzmann-averaged free energy
+   best = lowest_gibbs(thermos)                   # the single lowest-G conformer
+
+The ensemble free energy lies at or below the lowest conformer's by the mixing
+(conformational) entropy; use ``lowest_gibbs`` when you instead want the single
+dominant structure to carry forward (e.g. into :func:`reaction_free_energy`).
+
 Reactions and redox
 -------------------
 

tests/thermo/test_ensemble.py

@@ -0,0 +1,107 @@
+import math
+
+import pytest
+
+from ThermoScreening.thermo.ensemble import (
+    boltzmann_weights,
+    ensemble_free_energy,
+    lowest_gibbs,
+)
+
+_KB_HARTREE_PER_K = 1.380649e-23 / 4.3597447222071e-18  # k_B in Hartree/K
+_H_TO_KCAL = 627.5094740631
+
+
+class _FakeThermo:
+    """A stand-in exposing only the free energy the ensemble helpers use."""
+
+    def __init__(self, eegtot):
+        self._eegtot = eegtot
+
+    def total_EeGtot(self):
+        return self._eegtot
+
+
+def test_boltzmann_weights_degenerate():
+    weights = boltzmann_weights([_FakeThermo(-1.0), _FakeThermo(-1.0)])
+    assert weights == pytest.approx([0.5, 0.5])
+
+
+def test_boltzmann_weights_single_conformer():
+    assert boltzmann_weights([_FakeThermo(-3.14)]) == pytest.approx([1.0])
+
+
+def test_boltzmann_weights_normalised_and_ordered():
+    weights = boltzmann_weights([_FakeThermo(-1.0), _FakeThermo(-1.001), _FakeThermo(-0.999)])
+    assert sum(weights) == pytest.approx(1.0)
+    # lower free energy -> larger population
+    assert weights[1] > weights[0] > weights[2]
+
+
+def test_boltzmann_weights_known_ratio():
+    # gap of exactly k_B*T -> population ratio of e:1
+    temperature = 300.0
+    gap = _KB_HARTREE_PER_K * temperature
+    weights = boltzmann_weights([_FakeThermo(0.0), _FakeThermo(gap)], temperature=temperature)
+    assert weights[0] / weights[1] == pytest.approx(math.e)
+
+
+def test_ensemble_free_energy_single_is_identity():
+    assert ensemble_free_energy([_FakeThermo(-2.0)], unit="H") == pytest.approx(-2.0)
+
+
+def test_ensemble_free_energy_degenerate_mixing_entropy():
+    # two degenerate conformers -> G_ens = e - k_B*T*ln(2), below either one
+    temperature = 298.15
+    energy = -5.0
+    kt = _KB_HARTREE_PER_K * temperature
+    expected = energy - kt * math.log(2)
+    result = ensemble_free_energy([_FakeThermo(energy), _FakeThermo(energy)], temperature=temperature)
+    assert result == pytest.approx(expected)
+    assert result < energy
+
+
+def test_ensemble_free_energy_at_or_below_minimum():
+    thermos = [_FakeThermo(-1.0), _FakeThermo(-1.2), _FakeThermo(-0.8)]
+    g_ens = ensemble_free_energy(thermos, unit="H")
+    assert g_ens <= min(t.total_EeGtot() for t in thermos)
+
+
+def test_ensemble_free_energy_unit_conversion():
+    h = ensemble_free_energy([_FakeThermo(-1.0), _FakeThermo(-1.2)], unit="H")
+    kcal = ensemble_free_energy([_FakeThermo(-1.0), _FakeThermo(-1.2)], unit="kcal")
+    assert kcal == pytest.approx(h * _H_TO_KCAL, rel=1e-4)
+
+
+def test_ensemble_free_energy_rejects_unknown_unit():
+    with pytest.raises(ValueError, match="Unknown unit"):
+        ensemble_free_energy([_FakeThermo(0.0)], unit="furlong")
+
+
+def test_lowest_gibbs_returns_min_conformer():
+    a, b, c = _FakeThermo(-1.0), _FakeThermo(-1.5), _FakeThermo(-0.5)
+    assert lowest_gibbs([a, b, c]) is b
+
+
+@pytest.mark.parametrize("func", [boltzmann_weights, ensemble_free_energy, lowest_gibbs])
+def test_empty_ensemble_raises(func):
+    with pytest.raises(ValueError, match="at least one conformer"):
+        func([])
+
+
+@pytest.mark.parametrize("func", [boltzmann_weights, ensemble_free_energy])
+@pytest.mark.parametrize("temperature", [0.0, -10.0])
+def test_non_positive_temperature_raises(func, temperature):
+    with pytest.raises(ValueError, match="temperature must be positive"):
+        func([_FakeThermo(0.0)], temperature=temperature)
+
+
+def test_public_api_exported():
+    from ThermoScreening.thermo import boltzmann_weights as bw
+    from ThermoScreening.thermo import ensemble_free_energy as efe
+    from ThermoScreening.thermo import lowest_gibbs as lg
+    from ThermoScreening.thermo import ensemble
+
+    assert bw is ensemble.boltzmann_weights
+    assert efe is ensemble.ensemble_free_energy
+    assert lg is ensemble.lowest_gibbs