feos-org · prehner · Jan 7, 2026 · Jan 27, 2026 · Jan 27, 2026 · Jan 27, 2026
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -2,9 +2,9 @@ name: Test
 
 on:
   push:
-    branches: [main]
+    branches: [main, development]
   pull_request:
-    branches: [main]
+    branches: [main, development]
 
 env:
   CARGO_TERM_COLOR: always

diff --git a/.github/workflows/wheels.yml b/.github/workflows/wheels.yml
@@ -1,9 +1,9 @@
 name: Build Wheels
 on:
   push:
-    branches: [main]
+    branches: [main, development]
   pull_request:
-    branches: [main]
+    branches: [main, development]
 jobs:
   linux:
     runs-on: ubuntu-latest

diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -4,6 +4,22 @@ All notable changes to this project will be documented in this file.
 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).
 
+## [Breaking]
+### Added
+- Extended tp-flash algorithm to static numbers of components and enabled automatic differentiation for binary systems. [#336](https://github.com/feos-org/feos/pull/336)
+- Rewrote `PhaseEquilibrium::pure_p` to mirror `pure_t` and enabled automatic differentiation. [#337](https://github.com/feos-org/feos/pull/337)
+- Added `boiling_temperature` to the list of properties for parallel evaluations of gradients. [#337](https://github.com/feos-org/feos/pull/337)
+- Added the `Composition` trait to allow more flexibility in the creation of states and phase equilibria. [#330](https://github.com/feos-org/feos/pull/330)
+
+### Changed
+- Removed any assumptions about the total number of moles in a `State` or `PhaseEquilibrium`. Evaluating extensive properties now returns a `Result`. [#330](https://github.com/feos-org/feos/pull/330)
+
+### Removed
+- Removed the `StateBuilder` struct, because it is mostly obsolete with the addition of the `Composition` trait. [#330](https://github.com/feos-org/feos/pull/330)
+
+### Packaging
+- Updated `quantity` dependency to 0.13 and removed the `typenum` dependency. [#323](https://github.com/feos-org/feos/pull/323)
+
 ## [Unreleased]
 ### Added
 - Added `py-feos` tests to GitHub Actions and moved `pyo3/extension-module` feature to `pyproject.toml`. [#334](https://github.com/feos-org/feos/pull/334)

diff --git a/Cargo.toml b/Cargo.toml
@@ -22,7 +22,7 @@ keywords = [
 categories = ["science"]
 
 [workspace.dependencies]
-quantity = "0.12"
+quantity = "0.13"
 num-dual = "0.13"
 ndarray = "0.17"
 nalgebra = "0.34"
@@ -33,7 +33,6 @@ serde = "1.0"
 serde_json = "1.0"
 indexmap = "2.0"
 itertools = "0.14"
-typenum = "1.16"
 rayon = "1.11"
 petgraph = "0.8"
 rustdct = "0.7"

diff --git a/crates/feos-core/Cargo.toml b/crates/feos-core/Cargo.toml
@@ -25,7 +25,6 @@ serde = { workspace = true, features = ["derive"] }
 serde_json = { workspace = true, features = ["preserve_order"] }
 indexmap = { workspace = true, features = ["serde"] }
 rayon = { workspace = true, optional = true }
-typenum = { workspace = true }
 itertools = { workspace = true }
 
 [dev-dependencies]

diff --git a/crates/feos-core/src/ad/mod.rs b/crates/feos-core/src/ad/mod.rs
@@ -1,10 +1,10 @@
 use crate::DensityInitialization::Liquid;
 use crate::density_iteration::density_iteration;
-use crate::{FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
+use crate::{Composition, FeosResult, PhaseEquilibrium, ReferenceSystem, Residual};
 use nalgebra::{Const, SVector, U1, U2};
 #[cfg(feature = "rayon")]
 use ndarray::{Array1, Array2, ArrayView2, Zip};
-use num_dual::{Derivative, DualSVec, DualStruct};
+use num_dual::{Derivative, DualNum, DualSVec, DualStruct, first_derivative, partial2};
 use quantity::{Density, Pressure, Temperature};
 #[cfg(feature = "rayon")]
 use quantity::{KELVIN, KILO, METER, MOL, PASCAL};
@@ -57,15 +57,59 @@ pub trait PropertiesAD {
             let v2 = Gradient::from(v2);
             let x = Self::pure_molefracs();
 
-            let a1 = self.residual_molar_helmholtz_energy(t, v1, &x);
-            let a2 = self.residual_molar_helmholtz_energy(t, v2, &x);
+            let a1 = self.residual_helmholtz_energy(t, v1, &x);
+            let a2 = self.residual_helmholtz_energy(t, v2, &x);
             (a1, a2)
         };
 
         let p = -(a1 - a2 + t * (v2 / v1).ln()) / (v1 - v2);
         Ok(Pressure::from_reduced(p))
     }
 
+    fn boiling_temperature<const P: usize>(
+        &self,
+        pressure: Pressure,
+    ) -> FeosResult<Temperature<Gradient<P>>>
+    where
+        Self: Residual<U1, Gradient<P>>,
+    {
+        let eos_f64 = self.re();
+        let (temperature, [vapor_density, liquid_density]) =
+            PhaseEquilibrium::pure_p(&eos_f64, pressure, None, Default::default())?;
+
+        // implicit differentiation is implemented here instead of just calling pure_t with dual
+        // numbers, because for the first derivative, we can avoid calculating density derivatives.
+        let t = temperature.into_reduced();
+        let v1 = 1.0 / liquid_density.to_reduced();
+        let v2 = 1.0 / vapor_density.to_reduced();
+        let p = pressure.into_reduced();
+        let t = Gradient::from(t);
+        let t = t + {
+            let v1 = Gradient::from(v1);
+            let v2 = Gradient::from(v2);
+            let p = Gradient::from(p);
+            let x = Self::pure_molefracs();
+
+            let residual_entropy = |v| {
+                let (a, s) = first_derivative(
+                    partial2(
+                        |t, &v, x| self.lift().residual_helmholtz_energy(t, v, x),
+                        &v,
+                        &x,
+                    ),
+                    t,
+                );
+                (a, -s)
+            };
+            let (a1, s1) = residual_entropy(v1);
+            let (a2, s2) = residual_entropy(v2);
+
+            let ln_rho = (v1 / v2).ln();
+            (p * (v2 - v1) + (a2 - a1 + t * ln_rho)) / (s2 - s1 - ln_rho)
+        };
+        Ok(Temperature::from_reduced(t))
+    }
+
     fn equilibrium_liquid_density<const P: usize>(
         &self,
         temperature: Temperature,
@@ -111,6 +155,26 @@ pub trait PropertiesAD {
         )
     }
 
+    #[cfg(feature = "rayon")]
+    fn boiling_temperature_parallel<const P: usize>(
+        parameter_names: [String; P],
+        parameters: ArrayView2<f64>,
+        input: ArrayView2<f64>,
+    ) -> (Array1<f64>, Array2<f64>, Array1<bool>)
+    where
+        Self: ParametersAD<1>,
+    {
+        parallelize::<_, Self, _, _>(
+            parameter_names,
+            parameters,
+            input,
+            |eos: &Self::Lifted<Gradient<P>>, inp| {
+                eos.boiling_temperature(inp[0] * PASCAL)
+                    .map(|p| p.convert_into(KELVIN))
+            },
+        )
+    }
+
     #[cfg(feature = "rayon")]
     fn liquid_density_parallel<const P: usize>(
         parameter_names: [String; P],
@@ -151,20 +215,21 @@ pub trait PropertiesAD {
         )
     }
 
-    fn bubble_point_pressure<const P: usize>(
+    fn bubble_point_pressure<const P: usize, X: Composition<f64, U2>>(
         &self,
         temperature: Temperature,
         pressure: Option<Pressure>,
-        liquid_molefracs: SVector<f64, 2>,
+        liquid_molefracs: X,
     ) -> FeosResult<Pressure<Gradient<P>>>
     where
         Self: Residual<U2, Gradient<P>>,
     {
         let eos_f64 = self.re();
+        let (liquid_molefracs, _) = liquid_molefracs.into_molefracs(&eos_f64)?;
         let vle = PhaseEquilibrium::bubble_point(
             &eos_f64,
             temperature,
-            &liquid_molefracs,
+            liquid_molefracs,
             pressure,
             None,
             Default::default(),
@@ -184,37 +249,38 @@ pub trait PropertiesAD {
             let y = y.map(Gradient::from);
             let x = liquid_molefracs.map(Gradient::from);
 
-            let a_v = self.residual_molar_helmholtz_energy(t, v_v, &y);
+            let a_v = self.residual_helmholtz_energy(t, v_v, &y);
             let (p_l, mu_res_l, dp_l, dmu_l) = self.dmu_dv(t, v_l, &x);
             let vi_l = dmu_l / dp_l;
             let v_l = vi_l.dot(&y);
             let a_l = (mu_res_l - vi_l * p_l).dot(&y);
             (a_l, a_v, v_l, v_v)
         };
-        let rho_l = vle.liquid().partial_density.to_reduced();
+        let rho_l = vle.liquid().partial_density().to_reduced();
         let rho_l = [rho_l[0], rho_l[1]];
-        let rho_v = vle.vapor().partial_density.to_reduced();
+        let rho_v = vle.vapor().partial_density().to_reduced();
         let rho_v = [rho_v[0], rho_v[1]];
         let p = -(a_v - a_l
             + t * (y[0] * (rho_v[0] / rho_l[0]).ln() + y[1] * (rho_v[1] / rho_l[1]).ln() - 1.0))
             / (v_v - v_l);
         Ok(Pressure::from_reduced(p))
     }
 
-    fn dew_point_pressure<const P: usize>(
+    fn dew_point_pressure<const P: usize, X: Composition<f64, U2>>(
         &self,
         temperature: Temperature,
         pressure: Option<Pressure>,
-        vapor_molefracs: SVector<f64, 2>,
+        vapor_molefracs: X,
     ) -> FeosResult<Pressure<Gradient<P>>>
     where
         Self: Residual<U2, Gradient<P>>,
     {
         let eos_f64 = self.re();
+        let (vapor_molefracs, _) = vapor_molefracs.into_molefracs(&eos_f64)?;
         let vle = PhaseEquilibrium::dew_point(
             &eos_f64,
             temperature,
-            &vapor_molefracs,
+            vapor_molefracs,
             pressure,
             None,
             Default::default(),
@@ -234,16 +300,16 @@ pub trait PropertiesAD {
             let x = x.map(Gradient::from);
             let y = vapor_molefracs.map(Gradient::from);
 
-            let a_l = self.residual_molar_helmholtz_energy(t, v_l, &x);
+            let a_l = self.residual_helmholtz_energy(t, v_l, &x);
             let (p_v, mu_res_v, dp_v, dmu_v) = self.dmu_dv(t, v_v, &y);
             let vi_v = dmu_v / dp_v;
             let v_v = vi_v.dot(&x);
             let a_v = (mu_res_v - vi_v * p_v).dot(&x);
             (a_l, a_v, v_l, v_v)
         };
-        let rho_l = vle.liquid().partial_density.to_reduced();
+        let rho_l = vle.liquid().partial_density().to_reduced();
         let rho_l = [rho_l[0], rho_l[1]];
-        let rho_v = vle.vapor().partial_density.to_reduced();
+        let rho_v = vle.vapor().partial_density().to_reduced();
         let rho_v = [rho_v[0], rho_v[1]];
         let p = -(a_l - a_v
             + t * (x[0] * (rho_l[0] / rho_v[0]).ln() + x[1] * (rho_l[1] / rho_v[1]).ln() - 1.0))
@@ -265,12 +331,8 @@ pub trait PropertiesAD {
             parameters,
             input,
             |eos: &Self::Lifted<Gradient<P>>, inp| {
-                eos.bubble_point_pressure(
-                    inp[0] * KELVIN,
-                    Some(inp[2] * PASCAL),
-                    SVector::from([inp[1], 1.0 - inp[1]]),
-                )
-                .map(|p| p.convert_into(PASCAL))
+                eos.bubble_point_pressure(inp[0] * KELVIN, Some(inp[2] * PASCAL), inp[1])
+                    .map(|p| p.convert_into(PASCAL))
             },
         )
     }
@@ -289,12 +351,8 @@ pub trait PropertiesAD {
             parameters,
             input,
             |eos: &Self::Lifted<Gradient<P>>, inp| {
-                eos.dew_point_pressure(
-                    inp[0] * KELVIN,
-                    Some(inp[2] * PASCAL),
-                    SVector::from([inp[1], 1.0 - inp[1]]),
-                )
-                .map(|p| p.convert_into(PASCAL))
+                eos.dew_point_pressure(inp[0] * KELVIN, Some(inp[2] * PASCAL), inp[1])
+                    .map(|p| p.convert_into(PASCAL))
             },
         )
     }

diff --git a/crates/feos-core/src/cubic.rs b/crates/feos-core/src/cubic.rs
@@ -221,7 +221,7 @@ mod tests {
         let parameters = PengRobinsonParameters::new_pure(propane)?;
         let pr = PengRobinson::new(parameters);
         let options = SolverOptions::new().verbosity(Verbosity::Iter);
-        let cp = State::critical_point(&&pr, None, None, None, options)?;
+        let cp = State::critical_point(&&pr, (), None, None, options)?;
         println!("{} {}", cp.temperature, cp.pressure(Contributions::Total));
         assert_relative_eq!(cp.temperature, tc * KELVIN, max_relative = 1e-4);
         assert_relative_eq!(

diff --git a/crates/feos-core/src/density_iteration.rs b/crates/feos-core/src/density_iteration.rs
@@ -87,7 +87,7 @@ where
     let (a_res, da_res) = first_derivative(
         |molar_volume| {
             eos.lift()
-                .residual_molar_helmholtz_energy(t, molar_volume, &x)
+                .residual_helmholtz_energy(t, molar_volume, &x)
         },
         molar_volume,
     );

diff --git a/crates/feos-core/src/equation_of_state/mod.rs b/crates/feos-core/src/equation_of_state/mod.rs
@@ -1,9 +1,10 @@
-use crate::{ReferenceSystem, StateHD};
+use crate::ReferenceSystem;
+use crate::state::StateHD;
 use nalgebra::{
     Const, DVector, DefaultAllocator, Dim, Dyn, OVector, SVector, U1, allocator::Allocator,
 };
 use num_dual::DualNum;
-use quantity::{Energy, MolarEnergy, Moles, Temperature, Volume};
+use quantity::{Dimensionless, MolarEnergy, MolarVolume, Temperature};
 use std::ops::Deref;
 
 mod residual;
@@ -164,17 +165,17 @@ where
     fn ideal_gas_helmholtz_energy<D2: DualNum<f64, Inner = D> + Copy>(
         &self,
         temperature: Temperature<D2>,
-        volume: Volume<D2>,
-        moles: &Moles<OVector<D2, N>>,
-    ) -> Energy<D2> {
+        volume: MolarVolume<D2>,
+        moles: &OVector<D2, N>,
+    ) -> MolarEnergy<D2> {
         let total_moles = moles.sum();
         let molefracs = moles / total_moles;
-        let molar_volume = volume / total_moles;
+        let molar_volume = volume.into_reduced() / total_moles;
         MolarEnergy::from_reduced(self.ideal_gas_molar_helmholtz_energy(
             temperature.into_reduced(),
-            molar_volume.into_reduced(),
+            molar_volume,
             &molefracs,
-        )) * total_moles
+        )) * Dimensionless::new(total_moles)
     }
 }