prehner · Oct 21, 2026 · Oct 21, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 ### Fixed
 - Fixed bug in `_density_iteration` which used the residual Helmholtz energy as the pressure and the pressure as derivative. [#396](https://github.com/feos-org/feos/pull/396)
+- Fixed missing mole fraction factor in `State::thermodynamic_factor`. [#400](https://github.com/feos-org/feos/pull/400)
+- Fixed `State::henrys_law_constant` for mixed solvents. [#400](https://github.com/feos-org/feos/pull/400)
+- Fixed the acceleration in `PhaseEquilibrium::tp_flash` (Gibbs energy acceptance criterion not weighted by phase fractions and incomplete K-value history). [#400](https://github.com/feos-org/feos/pull/400)
 
 ## [0.10.3] - 2026-10-14
 ### Fixed

diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs
@@ -372,9 +372,9 @@ where
 
             self.rachford_rice_inplace(feed_state, &k)?;
             if let Some(k_vec) = k_vec
-                && i >= iterations - 3
+                && i + 4 >= iterations
             {
-                k_vec[i + 3 - iterations] = k.map(|ki| if ki > 0.0 { ki.ln() } else { 0.0 });
+                k_vec[i + 4 - iterations] = k.map(|ki| if ki > 0.0 { ki.ln() } else { 0.0 });
             }
         }
         Ok(false)
@@ -443,14 +443,16 @@ where
 
     // Total molar Gibbs energy excluding the constant contribution RT sum_i x_i ln(\Lambda_i^3)
     fn molar_gibbs_energy(&self) -> MolarEnergy {
-        self.states
-            .iter()
-            .fold(MolarEnergy::from_reduced(0.0), |acc, s| {
+        self.states.iter().zip(self.phase_fractions).fold(
+            MolarEnergy::from_reduced(0.0),
+            |acc, (s, phi)| {
                 let ln_rho_m1 = s.partial_density().to_reduced().map(|r| r.ln() - 1.0);
-                acc + s.residual_molar_helmholtz_energy()
+                acc + (s.residual_molar_helmholtz_energy()
                     + s.pressure(Contributions::Total) * s.molar_volume
-                    + RGAS * s.temperature * s.molefracs.dot(&ln_rho_m1)
-            })
+                    + RGAS * s.temperature * s.molefracs.dot(&ln_rho_m1))
+                    * phi
+            },
+        )
     }
 }
 

diff --git a/crates/feos-core/src/state/residual_properties.rs b/crates/feos-core/src/state/residual_properties.rs
@@ -449,7 +449,7 @@ impl<E: Residual + Subset> State<E> {
             .into_iter()
             .zip(&vle.vapor().molefracs)
             .for_each(|(i, &y)| molefracs_vapor[i] = y);
-        let vapor = State::new(eos, temperature, vle.vapor().density, molefracs.clone())?;
+        let vapor = State::new(eos, temperature, vle.vapor().density, molefracs_vapor)?;
 
         // Determine the Henry's law coefficients and return only those of the Henry components
         let p = vle.vapor().pressure(Contributions::Total).into_reduced();
@@ -473,9 +473,10 @@ impl<E: Residual> State<E> {
     /// Thermodynamic factor: $\Gamma_{ij}=\delta_{ij}+x_i\left(\frac{\partial\ln\varphi_i}{\partial x_j}\right)_{T,p,\Sigma}$
     pub fn thermodynamic_factor(&self) -> DMatrix<f64> {
         let dln_phi_dnj = self.n_dln_phi_dnj();
+        let x = &self.molefracs;
         let n = self.eos.components() - 1;
         DMatrix::from_fn(n, n, |i, j| {
-            dln_phi_dnj[(i, j)] - dln_phi_dnj[(i, n)] + if i == j { 1.0 } else { 0.0 }
+            x[i] * (dln_phi_dnj[(i, j)] - dln_phi_dnj[(i, n)]) + if i == j { 1.0 } else { 0.0 }
         })
     }
 }

diff --git a/crates/feos/tests/pcsaft/properties.rs b/crates/feos/tests/pcsaft/properties.rs
@@ -1,10 +1,12 @@
 use approx::assert_relative_eq;
 use feos::pcsaft::{PcSaft, PcSaftParameters};
 use feos_core::parameter::IdentifierOption;
-use feos_core::{DensityInitialization::Vapor, Residual, State};
+use feos_core::DensityInitialization::{Liquid, Vapor};
+use feos_core::{Contributions, PhaseEquilibrium, Residual, State};
 use nalgebra::dvector;
 use quantity::*;
 use std::error::Error;
+use std::sync::Arc;
 
 #[test]
 fn test_dln_phi_dp() -> Result<(), Box<dyn Error>> {
@@ -42,3 +44,42 @@ fn test_virial_is_not_nan() -> Result<(), Box<dyn Error>> {
     assert!(!virial_b.is_nan());
     Ok(())
 }
+
+fn esper2023(components: Vec<&str>) -> Result<PcSaft, Box<dyn Error>> {
+    Ok(PcSaft::new(PcSaftParameters::from_json(
+        components,
+        "../../parameters/pcsaft/esper2023.json",
+        Some("../../parameters/pcsaft/rehner2023_binary.json"),
+        IdentifierOption::Name,
+    )?))
+}
+
+#[test]
+fn test_thermodynamic_factor_component_order() -> Result<(), Box<dyn Error>> {
+    let ew = esper2023(vec!["ethanol", "water"])?;
+    let we = esper2023(vec!["water", "ethanol"])?;
+    let t = 300.0 * KELVIN;
+    for x in [0.2, 0.5, 1e-6] {
+        let s1 = State::new_npt(&&ew, t, BAR, dvector![x, 1.0 - x], Some(Liquid))?;
+        let s2 = State::new_npt(&&we, t, BAR, dvector![1.0 - x, x], Some(Liquid))?;
+        let g1 = s1.thermodynamic_factor()[(0, 0)];
+        let g2 = s2.thermodynamic_factor()[(0, 0)];
+        assert_relative_eq!(g1, g2, max_relative = 1e-6);
+    }
+    let s = State::new_npt(&&ew, t, BAR, dvector![1e-8, 1.0 - 1e-8], Some(Liquid))?;
+    assert_relative_eq!(s.thermodynamic_factor()[(0, 0)], 1.0, epsilon = 1e-6);
+    Ok(())
+}
+
+#[test]
+fn test_henrys_law_constant_mixed_solvent() -> Result<(), Box<dyn Error>> {
+    let eos = Arc::new(esper2023(vec!["methane", "propane", "hexane"])?);
+    let t = 300.0 * KELVIN;
+    let h = State::henrys_law_constant(&eos, t, &dvector![0.0, 0.5, 0.5])?[0];
+    let x1 = 1e-8;
+    let z = dvector![x1, 0.5 * (1.0 - x1), 0.5 * (1.0 - x1)];
+    let vle = PhaseEquilibrium::bubble_point(&eos, t, &z, None, None, Default::default())?;
+    let h_lim = vle.vapor().pressure(Contributions::Total) * vle.vapor().molefracs[0] / x1;
+    assert_relative_eq!(h, h_lim, max_relative = 1e-6);
+    Ok(())
+}
diff --git a/crates/feos/tests/pcsaft/tp_flash.rs b/crates/feos/tests/pcsaft/tp_flash.rs
@@ -102,3 +102,26 @@ fn test_tp_flash_zero_component() -> FeosResult<()> {
     );
     Ok(())
 }
+
+#[test]
+fn test_tp_flash_acceleration() -> FeosResult<()> {
+    // With working acceleration, this LLE converges in ~36 steps (previously ~180).
+    let eos = PcSaft::new(PcSaftParameters::from_json(
+        vec!["methanol", "hexane"],
+        "../../parameters/pcsaft/esper2023.json",
+        Some("../../parameters/pcsaft/rehner2023_binary.json"),
+        IdentifierOption::Name,
+    )?);
+    let options = SolverOptions::new().max_iter(15).tol(1e-12);
+    let vle = PhaseEquilibrium::tp_flash(
+        &&eos,
+        280.0 * KELVIN,
+        BAR,
+        &dvector![0.5, 0.5],
+        None,
+        options,
+        None,
+    )?;
+    assert!((vle.vapor().molefracs[0] - vle.liquid().molefracs[0]).abs() > 0.1);
+    Ok(())
+}