prehner · Oct 24, 2026 · Oct 23, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,6 +5,12 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+### Fixed
+- Fixed `State::critical_point_binary` at fixed pressure, which evaluated the pressure at the wrong volume and never converged. [#402](https://github.com/feos-org/feos/pull/402)
+- Fixed `PhaseEquilibrium::tp_flash` with an initial state returning the previous phase split if only the feed composition changes. [#402](https://github.com/feos-org/feos/pull/402)
+- Fixed implicit derivatives in `State::new_nph`, `new_nps`, `new_nth`, `new_nts`, and `new_nvu` being lost if the initial guess is already converged.
+- Fixed `State::residual_chemical_potential_contributions` passing the temperature as molar volume. [#402](https://github.com/feos-org/feos/pull/402)
+- Fixed duplicate point in `PhaseDiagram::dew_point_line`. [#402](https://github.com/feos-org/feos/pull/402)
 
 ## [0.10.4] - 2026-10-21
 ### Fixed

diff --git a/crates/feos-core/src/phase_equilibria/phase_envelope.rs b/crates/feos-core/src/phase_equilibria/phase_envelope.rs
@@ -107,9 +107,9 @@ impl<E: Residual> PhaseDiagram<E, 2> {
 
         let min_pressure = vle.as_ref().unwrap().vapor().pressure(Contributions::Total);
         let p_c = sc.pressure(Contributions::Total);
-        let max_pressure =
-            min_pressure + (p_c - min_pressure) * ((n_p - 2) as f64 / (n_p - 1) as f64);
-        let pressures = Pressure::linspace(min_pressure, max_pressure, n_p);
+        // the first point (min_pressure) is already part of the temperature sweep
+        let dp = (p_c - min_pressure) / ((n_p + 1) as f64);
+        let pressures = Pressure::linspace(min_pressure + dp, p_c - dp, n_p);
 
         for pi in &pressures {
             let t_init = vle.as_ref().map(|vle| vle.vapor().temperature);

diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs
@@ -127,16 +127,17 @@ where
     ) -> FeosResult<PhaseEquilibrium<E, 2, N>> {
         // initialization
         if let Some(initial_state) = initial_state {
+            // the feed may have changed: split it with the old K values
             let mut init = initial_state.clone();
-            init.update_states(
-                self,
-                initial_state.vapor().molefracs.clone(),
-                initial_state.liquid().molefracs.clone(),
-                initial_state.vapor_phase_fraction(),
-            )?;
-            let vle = self.tp_flash_(init, options, non_volatile_components.clone());
-            if vle.is_ok() {
-                return vle;
+            let k = initial_state
+                .vapor()
+                .molefracs
+                .component_div(&initial_state.liquid().molefracs);
+            if init.rachford_rice_inplace(self, &k).is_ok() {
+                let vle = self.tp_flash_(init, options, non_volatile_components.clone());
+                if vle.is_ok() {
+                    return vle;
+                }
             }
         }
 

diff --git a/crates/feos-core/src/state/critical_point.rs b/crates/feos-core/src/state/critical_point.rs
@@ -570,7 +570,7 @@ where
         &temperature,
         &molefracs,
     );
-    let (_, da) = first_derivative(a, D::one());
+    let (_, da) = first_derivative(a, density.recip());
     let p_calc = -da + density * temperature;
     SVector::from([c1, c2, -p_calc + pressure])
 }

diff --git a/crates/feos-core/src/state/mod.rs b/crates/feos-core/src/state/mod.rs
@@ -620,16 +620,15 @@ where
     let maxiter = 50;
 
     for _ in 0..maxiter {
-        let (fx, dfx, mut state) = f(x0)?;
+        let (fx, dfx, _) = f(x0)?;
         let x = x0 - fx / dfx;
         if is_close(x.re(), x0.re(), atol, rtol) {
             // Ensure that at least NDERIV iterations are performed (for implicit AD)
             for _ in 0..D::NDERIV {
-                let (fx, dfx, s) = f(x0)?;
+                let (fx, dfx, _) = f(x0)?;
                 x0 -= fx / dfx;
-                state = s;
             }
-            return Ok(state);
+            return Ok(f(x0)?.2);
         }
         x0 = x;
     }

diff --git a/crates/feos-core/src/state/residual_properties.rs b/crates/feos-core/src/state/residual_properties.rs
@@ -560,7 +560,7 @@ where
         component: usize,
     ) -> Vec<(&'static str, MolarEnergy<D>)> {
         let t = Dual::from_re(self.temperature.into_reduced());
-        let v = Dual::from_re(self.temperature.into_reduced());
+        let v = Dual::from_re(self.density.into_reduced().recip());
         let mut x = self.molefracs.map(Dual::from_re);
         x[component].eps = D::one();
         let contributions = self.eos.lift().helmholtz_energy_contributions(t, v, &x);

diff --git a/crates/feos/tests/pcsaft/critical_point.rs b/crates/feos/tests/pcsaft/critical_point.rs
@@ -1,7 +1,7 @@
 use approx::assert_relative_eq;
 use feos::pcsaft::{PcSaft, PcSaftParameters};
 use feos_core::parameter::IdentifierOption;
-use feos_core::{SolverOptions, State};
+use feos_core::{Contributions, SolverOptions, State};
 use nalgebra::dvector;
 use quantity::*;
 use std::error::Error;
@@ -75,3 +75,29 @@ fn test_critical_point_limits() -> Result<(), Box<dyn Error>> {
     assert_eq!(cp_pure[1].density, cp_2.density);
     Ok(())
 }
+
+#[test]
+fn test_critical_point_binary_pressure() -> Result<(), Box<dyn Error>> {
+    let eos = PcSaft::new(PcSaftParameters::from_json(
+        vec!["methane", "propane"],
+        "../../parameters/pcsaft/esper2023.json",
+        Some("../../parameters/pcsaft/rehner2023_binary.json"),
+        IdentifierOption::Name,
+    )?);
+    let options = SolverOptions::default();
+    let cp_t = State::critical_point_binary(&&eos, 250.0 * KELVIN, None, None, None, options)?;
+    let p = cp_t.pressure(Contributions::Total);
+    let x = cp_t.molefracs[0];
+    let cp_p = State::critical_point_binary(
+        &&eos,
+        p,
+        Some(cp_t.temperature),
+        Some([x, 1.0 - x]),
+        Some(cp_t.density),
+        options,
+    )?;
+    assert_relative_eq!(cp_p.temperature, cp_t.temperature, max_relative = 1e-8);
+    assert_relative_eq!(cp_p.molefracs[0], x, max_relative = 1e-6);
+    assert_relative_eq!(cp_p.density, cp_t.density, max_relative = 1e-6);
+    Ok(())
+}
diff --git a/crates/feos/tests/pcsaft/properties.rs b/crates/feos/tests/pcsaft/properties.rs
@@ -1,7 +1,7 @@
 use approx::assert_relative_eq;
 use feos::pcsaft::{PcSaft, PcSaftParameters};
-use feos_core::parameter::IdentifierOption;
 use feos_core::DensityInitialization::{Liquid, Vapor};
+use feos_core::parameter::IdentifierOption;
 use feos_core::{Contributions, PhaseEquilibrium, Residual, State};
 use nalgebra::dvector;
 use quantity::*;
@@ -83,3 +83,18 @@ fn test_henrys_law_constant_mixed_solvent() -> Result<(), Box<dyn Error>> {
     assert_relative_eq!(h, h_lim, max_relative = 1e-6);
     Ok(())
 }
+
+#[test]
+fn test_residual_chemical_potential_contributions() -> Result<(), Box<dyn Error>> {
+    let eos = esper2023(vec!["ethanol", "water"])?;
+    let s = State::new_npt(&&eos, 300.0 * KELVIN, BAR, dvector![0.3, 0.7], Some(Liquid))?;
+    let mu_res = s.residual_chemical_potential();
+    for i in 0..2 {
+        let mu_i = s
+            .residual_chemical_potential_contributions(i)
+            .into_iter()
+            .fold(0.0 * JOULE / MOL, |acc, (_, mu)| acc + mu);
+        assert_relative_eq!(mu_i, mu_res.get(i), max_relative = 1e-10);
+    }
+    Ok(())
+}
diff --git a/crates/feos/tests/pcsaft/px_flashes.rs b/crates/feos/tests/pcsaft/px_flashes.rs
@@ -2,8 +2,8 @@ use approx::assert_relative_eq;
 use feos::ideal_gas::Joback;
 use feos::pcsaft::PcSaftBinary;
 use feos_core::{
-    Contributions, EquationOfState, FeosResult, IdealGasAD, PhaseEquilibrium, ReferenceSystem,
-    SolverOptions, Verbosity, ad::ParametersAD,
+    Contributions, DensityInitialization, EquationOfState, FeosResult, IdealGasAD,
+    PhaseEquilibrium, ReferenceSystem, SolverOptions, State, Verbosity, ad::ParametersAD,
 };
 use nalgebra::U1;
 use num_dual::{DualStruct, DualVec};
@@ -142,3 +142,53 @@ fn test_ps_flash() -> FeosResult<()> {
 
     Ok(())
 }
+
+#[test]
+fn test_new_nph_warm_start_derivative() -> FeosResult<()> {
+    let params = [
+        [1.5, 3.4, 180.0, 2.2, 0.03, 2500., 2.0, 1.0],
+        [4.5, 3.6, 250.0, 1.2, 0.015, 1500., 1.0, 2.0],
+    ];
+    let kij = 0.15;
+    let joback = [
+        Joback([380., 0.0, 0.0, 0.0, 0.0]),
+        Joback([210., 0.0, 0.0, 0.0, 0.0]),
+    ];
+    let eos = EquationOfState::new(joback.clone(), PcSaftBinary::new(params, kij));
+    let (p, t, x) = (5.0 * BAR, 600.0 * KELVIN, 0.3);
+    let h = State::new_npt(&eos, t, p, x, Some(DensityInitialization::Vapor))?
+        .molar_enthalpy(Contributions::Total);
+
+    let mut flat_params: Vec<f64> = params[0].to_vec();
+    flat_params.extend_from_slice(&params[1]);
+    flat_params.push(kij);
+    let pcsaft_ad = PcSaftBinary::<f64, 8>::seed_derivatives(&flat_params, ["k_ij"]);
+    let joback_ad = joback.each_ref().map(|j| j.lift());
+    let eos_ad = EquationOfState::new(joback_ad, pcsaft_ad);
+
+    let dkij = 1e-7;
+    let eos_h = EquationOfState::new(joback.clone(), PcSaftBinary::new(params, kij + dkij));
+    let t_h = State::new_nph(&eos_h, p, h, x, Some(DensityInitialization::Vapor), Some(t))?;
+    let dt_h = (t_h.temperature - t).into_reduced() / dkij;
+
+    // starting from the converged temperature, the derivative must not be lost
+    for t0 in [500.0 * KELVIN, t] {
+        let s = State::new_nph(
+            &eos_ad,
+            Pressure::from_inner(&p),
+            MolarEnergy::from_inner(&h),
+            DualVec::from_inner(&x),
+            Some(DensityInitialization::Vapor),
+            Some(Temperature::from_inner(&t0)),
+        )?;
+        let [[dt]] = s
+            .temperature
+            .into_reduced()
+            .eps
+            .unwrap_generic(U1, U1)
+            .data
+            .0;
+        assert_relative_eq!(dt, dt_h, max_relative = 1e-4);
+    }
+    Ok(())
+}
diff --git a/crates/feos/tests/pcsaft/tp_flash.rs b/crates/feos/tests/pcsaft/tp_flash.rs
@@ -125,3 +125,33 @@ fn test_tp_flash_acceleration() -> FeosResult<()> {
     assert!((vle.vapor().molefracs[0] - vle.liquid().molefracs[0]).abs() > 0.1);
     Ok(())
 }
+
+#[test]
+fn test_tp_flash_initial_state_new_feed() -> FeosResult<()> {
+    let eos = PcSaft::new(PcSaftParameters::from_json(
+        vec!["propane", "n-butane", "pentane"],
+        "../../parameters/pcsaft/esper2023.json",
+        Some("../../parameters/pcsaft/rehner2023_binary.json"),
+        IdentifierOption::Name,
+    )?);
+    let (t, p) = (300.0 * KELVIN, 4.0 * BAR);
+    let options = SolverOptions::new().tol(1e-12);
+    let mut vle: Option<PhaseEquilibrium<_, 2>> = None;
+    for z1 in [0.4, 0.5, 0.6, 0.7] {
+        let z = dvector![z1, 0.5 * (1.0 - z1), 0.5 * (1.0 - z1)];
+        let cold = PhaseEquilibrium::tp_flash(&&eos, t, p, &z, None, options, None)?;
+        let warm = PhaseEquilibrium::tp_flash(&&eos, t, p, &z, vle.as_ref(), options, None)?;
+        assert_relative_eq!(
+            warm.vapor_phase_fraction(),
+            cold.vapor_phase_fraction(),
+            max_relative = 1e-8
+        );
+        assert_relative_eq!(
+            warm.liquid().molefracs[0],
+            cold.liquid().molefracs[0],
+            max_relative = 1e-8
+        );
+        vle = Some(warm);
+    }
+    Ok(())
+}