Fix Redlich-Kister excess enthalpy: standard polynomial form, drop spurious 0.5 factor

Author: milanofthe

src/pathsim_chem/thermodynamics/enthalpy.py

@@ -301,22 +301,24 @@ class ExcessEnthalpyRedlichKister(Function):
     Computes the molar excess enthalpy using a Redlich-Kister polynomial
     expansion. This is a flexible, empirical model for representing binary
     excess properties. For a multi-component mixture, the excess enthalpy
-    is computed as a sum of binary pair contributions:
+    is summed over the unique binary pairs supplied:
 
     .. math::
 
-        h^E = \frac{1}{2} \sum_i \sum_j h^E_{ij}
+        h^E = \sum_{(i,j)} h^E_{ij}
 
-    Each binary pair contribution:
+    Each binary pair contribution follows the standard Redlich-Kister form
 
     .. math::
 
         h^E_{ij} = \frac{x_i x_j}{x_i + x_j}
-        \left(A(T) x_d + B(T) x_d^2 + C(T) x_d^3 + \ldots\right)
+        \left(A_0(T) + A_1(T) x_d + A_2(T) x_d^2 + \ldots\right)
 
-    where :math:`x_d = x_i - x_j` and the coefficients :math:`A(T)`,
-    :math:`B(T)`, ... are temperature-dependent polynomials:
-    :math:`A(T) = a_0 + a_1 T + a_2 T^2 + \ldots`
+    where :math:`x_d = x_i - x_j` and each :math:`A_k(T)` is a temperature
+    polynomial :math:`A_k(T) = a_{k,0} + a_{k,1} T + a_{k,2} T^2 + \ldots`.
+    For a true binary system :math:`x_i + x_j = 1` and the prefactor
+    reduces to :math:`x_i x_j`; in multicomponent mixtures the
+    :math:`/(x_i+x_j)` term provides the symmetric extension.
 
     **Input port:** ``T`` -- temperature [K].
 
@@ -328,12 +330,15 @@ class ExcessEnthalpyRedlichKister(Function):
         Mole fractions [N].
     coeffs : dict
         Redlich-Kister coefficients keyed by binary pair ``(i, j)`` as tuples.
-        Each value is a list of polynomial coefficient arrays, one per
-        Redlich-Kister term. Example for a single pair (0,1) with 3 terms::
+        Provide each unordered pair only once (e.g. ``(0, 1)``, not also
+        ``(1, 0)``). Each value is a list of polynomial coefficient arrays,
+        one per Redlich-Kister term. Example for a single pair (0,1) with
+        3 terms::
 
             {(0, 1): [[a0, a1], [b0, b1], [c0]]}
 
-        This gives A(T)=a0+a1*T, B(T)=b0+b1*T, C(T)=c0.
+        This gives A_0(T)=a0+a1*T, A_1(T)=b0+b1*T, A_2(T)=c0, where
+        :math:`A_k` multiplies :math:`x_d^k`.
     """
 
     input_port_labels = {"T": 0}
@@ -358,9 +363,9 @@ def _eval(self, T):
 
             xd = xi - xj
 
-            # evaluate each Redlich-Kister term
+            # evaluate Redlich-Kister polynomial: A_0 + A_1*xd + A_2*xd^2 + ...
             pair_hE = 0.0
-            xd_power = xd
+            xd_power = 1.0
             for poly_coeffs in terms:
                 # temperature polynomial: c0 + c1*T + c2*T^2 + ...
                 coeff_val = 0.0
@@ -373,7 +378,7 @@ def _eval(self, T):
 
             hE += xi * xj / x_sum * pair_hE
 
-        return 0.5 * hE
+        return hE
 
 
 # ENTHALPY DEPARTURE FROM EOS (9.7) =====================================================

tests/thermodynamics/test_enthalpy.py

@@ -201,28 +201,40 @@ def test_init(self):
         )
         self.assertEqual(hE.n, 2)
 
-    def test_symmetric_binary(self):
-        # Simplest case: one-term Redlich-Kister with constant A
-        # h^E_ij = x_i*x_j/(x_i+x_j) * A * (x_i-x_j)
-        # For x=[0.5, 0.5]: h^E = 0.5 * 0.5 / 1.0 * A * 0 = 0
+    def test_constant_coeff_only(self):
+        # One-term Redlich-Kister with constant A_0:
+        # h^E = x_i*x_j/(x_i+x_j) * A_0
+        # For x=[0.5, 0.5], A_0=1000: h^E = 0.25/1.0 * 1000 = 250
         hE = ExcessEnthalpyRedlichKister(
             x=[0.5, 0.5],
             coeffs={(0, 1): [[1000.0]]},
         )
         result = eval_block_T(hE, 300)
-        self.assertAlmostEqual(result, 0.0, places=10)
+        self.assertAlmostEqual(result, 250.0, places=10)
 
-    def test_asymmetric_composition(self):
-        # x=[0.3, 0.7], one-term A=2000
-        # h^E = 0.5 * 0.3*0.7/1.0 * 2000 * (0.3-0.7) = 0.5 * 0.21 * 2000 * (-0.4) = -84
+    def test_two_term_polynomial(self):
+        # x=[0.3, 0.7], coeffs A_0=2000, A_1=500
+        # h^E = 0.21 * (2000 + 500*(-0.4)) = 0.21 * 1800 = 378
         hE = ExcessEnthalpyRedlichKister(
             x=[0.3, 0.7],
-            coeffs={(0, 1): [[2000.0]]},
+            coeffs={(0, 1): [[2000.0], [500.0]]},
         )
         result = eval_block_T(hE, 300)
-        expected = 0.5 * 0.3 * 0.7 / 1.0 * 2000.0 * (0.3 - 0.7)
+        expected = 0.3 * 0.7 / 1.0 * (2000.0 + 500.0 * (0.3 - 0.7))
         self.assertAlmostEqual(result, expected, places=5)
 
+    def test_antisymmetric_under_pair_swap(self):
+        # Odd-order RK terms should flip sign when (i,j) <-> (j,i),
+        # since x_d = x_i - x_j flips. Verify by swapping x values.
+        hE_a = ExcessEnthalpyRedlichKister(
+            x=[0.3, 0.7], coeffs={(0, 1): [[0.0], [1000.0]]},  # A_1 only
+        )
+        hE_b = ExcessEnthalpyRedlichKister(
+            x=[0.7, 0.3], coeffs={(0, 1): [[0.0], [1000.0]]},
+        )
+        self.assertAlmostEqual(eval_block_T(hE_a, 300),
+                               -eval_block_T(hE_b, 300), places=8)
+
     def test_temperature_dependent_coeffs(self):
         # A(T) = 1000 + 2*T
         hE = ExcessEnthalpyRedlichKister(