feos-org · prehner · Jan 15, 2026 · Dec 17, 2025 · Jan 9, 2026 · Jan 9, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,6 +5,9 @@ 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 the calculation of critical points and tp-flashes when one or more of the mixture components have zero composition. [#331](https://github.com/feos-org/feos/pull/331)
+
 ### Packaging
 - Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323)
 

diff --git a/crates/feos-core/src/phase_equilibria/tp_flash.rs b/crates/feos-core/src/phase_equilibria/tp_flash.rs
@@ -265,17 +265,13 @@ impl<E: Residual> PhaseEquilibrium<E, 2> {
 
             // check for convergence
             *iter += 1;
-            let mut res_vec = ln_phi_l - ln_phi_v
-                + self
-                    .liquid()
-                    .molefracs
-                    .component_div(&self.vapor().molefracs)
-                    .map(|i| if i > 0.0 { i.ln() } else { 0.0 });
+            let mut res_vec = k.component_mul(&self.liquid().molefracs) - &self.vapor().molefracs;
 
             // Set residuum to 0 for non-volatile components
             if let Some(nvc) = non_volatile_components.as_ref() {
                 nvc.iter().for_each(|&c| res_vec[c] = 0.0);
             }
+
             let res = res_vec.norm();
             log_iter!(
                 verbosity,

diff --git a/crates/feos-core/src/state/critical_point.rs b/crates/feos-core/src/state/critical_point.rs
@@ -10,6 +10,7 @@ use num_dual::{
     implicit_derivative_binary, implicit_derivative_vec, jacobian, partial, partial2,
     third_derivative,
 };
+use num_traits::Zero;
 use quantity::{Density, Pressure, Temperature};
 
 const MAX_ITER_CRIT_POINT: usize = 50;
@@ -537,7 +538,13 @@ where
             let n = n + sqrt_z.component_mul(&u).map(Dual3::from_re) * s;
             let t = Dual3::from_re(temperature);
             let v = Dual3::from_re(molar_volume);
-            let ig = n.dot(&n.map(|n| (n / v).ln() - 1.0));
+            let ig = n.dot(&n.map(|n| {
+                if n.is_zero() {
+                    Dual3::zero()
+                } else {
+                    (n / v).ln() - 1.0
+                }
+            }));
             eos.lift().residual_helmholtz_energy(t, v, &n) / t + ig
         },
         D::from(0.0),

diff --git a/crates/feos/tests/pcsaft/critical_point.rs b/crates/feos/tests/pcsaft/critical_point.rs
@@ -1,10 +1,11 @@
 use approx::assert_relative_eq;
 use feos::pcsaft::{PcSaft, PcSaftParameters};
-use feos_core::State;
 use feos_core::parameter::IdentifierOption;
+use feos_core::{SolverOptions, State};
 use nalgebra::dvector;
 use quantity::*;
 use std::error::Error;
+use std::sync::Arc;
 
 #[test]
 fn test_critical_point_pure() -> Result<(), Box<dyn Error>> {
@@ -46,3 +47,31 @@ fn test_critical_point_mix() -> Result<(), Box<dyn Error>> {
     );
     Ok(())
 }
+
+#[test]
+fn test_critical_point_limits() -> Result<(), Box<dyn Error>> {
+    let params = PcSaftParameters::from_json(
+        vec!["propane", "butane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?;
+    let options = SolverOptions {
+        verbosity: feos_core::Verbosity::Iter,
+        ..Default::default()
+    };
+    let saft = Arc::new(PcSaft::new(params));
+    let cp_pure = State::critical_point_pure(&saft, None, None, options)?;
+    println!("{} {}", cp_pure[0], cp_pure[1]);
+    let molefracs = dvector![0.0, 1.0];
+    let cp_2 = State::critical_point(&saft, Some(&molefracs), None, None, options)?;
+    println!("{}", cp_2);
+    let molefracs = dvector![1.0, 0.0];
+    let cp_1 = State::critical_point(&saft, Some(&molefracs), None, None, options)?;
+    println!("{}", cp_1);
+    assert_eq!(cp_pure[0].temperature, cp_1.temperature);
+    assert_eq!(cp_pure[0].density, cp_1.density);
+    assert_eq!(cp_pure[1].temperature, cp_2.temperature);
+    assert_eq!(cp_pure[1].density, cp_2.density);
+    Ok(())
+}
diff --git a/crates/feos/tests/pcsaft/tp_flash.rs b/crates/feos/tests/pcsaft/tp_flash.rs
@@ -4,7 +4,6 @@ use feos_core::parameter::IdentifierOption;
 use feos_core::{Contributions, FeosResult, PhaseEquilibrium, SolverOptions};
 use nalgebra::dvector;
 use quantity::*;
-use std::error::Error;
 
 fn read_params(components: Vec<&str>) -> FeosResult<PcSaftParameters> {
     PcSaftParameters::from_json(
@@ -16,7 +15,7 @@ fn read_params(components: Vec<&str>) -> FeosResult<PcSaftParameters> {
 }
 
 #[test]
-fn test_tp_flash() -> Result<(), Box<dyn Error>> {
+fn test_tp_flash() -> FeosResult<()> {
     let propane = PcSaft::new(read_params(vec!["propane"])?);
     let butane = PcSaft::new(read_params(vec!["butane"])?);
     let t = 250.0 * KELVIN;
@@ -63,3 +62,51 @@ fn test_tp_flash() -> Result<(), Box<dyn Error>> {
     );
     Ok(())
 }
+
+#[test]
+fn test_tp_flash_zero_component() -> FeosResult<()> {
+    let eos_full = PcSaft::new(PcSaftParameters::from_json(
+        vec!["propane", "butane", "hexane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?);
+    let eos_binary = PcSaft::new(PcSaftParameters::from_json(
+        vec!["butane", "hexane"],
+        "tests/pcsaft/test_parameters.json",
+        None,
+        IdentifierOption::Name,
+    )?);
+    let options = SolverOptions {
+        verbosity: feos_core::Verbosity::Iter,
+        ..Default::default()
+    };
+    let vle_full = PhaseEquilibrium::tp_flash(
+        &&eos_full,
+        300.0 * KELVIN,
+        1.2 * BAR,
+        &(dvector![0.0, 0.5, 0.5] * MOL),
+        None,
+        options,
+        None,
+    )?;
+    let vle_binary = PhaseEquilibrium::tp_flash(
+        &&eos_binary,
+        300.0 * KELVIN,
+        1.2 * BAR,
+        &(dvector![0.5, 0.5] * MOL),
+        None,
+        options,
+        None,
+    )?;
+    println!("{vle_full}\n{vle_binary}");
+    assert_eq!(
+        vle_full.liquid().molefracs[1],
+        vle_binary.liquid().molefracs[0]
+    );
+    assert_eq!(
+        vle_full.vapor().molefracs[1],
+        vle_binary.vapor().molefracs[0]
+    );
+    Ok(())
+}