feos-org · prehner · Feb 12, 2026 · Jan 27, 2026
diff --git a/crates/feos-core/src/ad/mod.rs b/crates/feos-core/src/ad/mod.rs
@@ -4,7 +4,7 @@ use crate::{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};
@@ -66,6 +66,50 @@ pub trait PropertiesAD {
         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_molar_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],

diff --git a/crates/feos-core/src/equation_of_state/residual.rs b/crates/feos-core/src/equation_of_state/residual.rs
@@ -1,6 +1,9 @@
 use crate::{FeosError, FeosResult, ReferenceSystem, state::StateHD};
+use nalgebra::SVector;
 use nalgebra::{DVector, DefaultAllocator, Dim, Dyn, OMatrix, OVector, U1, allocator::Allocator};
-use num_dual::{DualNum, Gradients, partial, partial2, second_derivative, third_derivative};
+use num_dual::{
+    DualNum, Gradients, hessian, partial, partial2, second_derivative, third_derivative,
+};
 use quantity::ad::first_derivative;
 use quantity::*;
 use std::ops::{Deref, Div};
@@ -313,12 +316,41 @@ where
             molar_volume,
         );
         (
-            a * density,
+            a,
             -da + temperature * density,
             molar_volume * molar_volume * d2a + temperature,
         )
     }
 
+    /// calculates a_res, p, s_res, dp_drho, dp_dt
+    fn p_dpdrho_dpdt(
+        &self,
+        temperature: D,
+        density: D,
+        molefracs: &OVector<D, N>,
+    ) -> (D, D, D, D, D) {
+        let molar_volume = density.recip();
+        let (a, da, d2a) = hessian::<_, _, _, nalgebra::U2, _>(
+            partial(
+                |vt: SVector<_, 2>, x: &OVector<_, N>| {
+                    let [[v, t]] = vt.data.0;
+                    self.lift().residual_molar_helmholtz_energy(t, v, x)
+                },
+                molefracs,
+            ),
+            &SVector::from([molar_volume, temperature]),
+        );
+        let [[da_dv, da_dt]] = da.data.0;
+        let [[d2a_dv2, d2a_dvdt], _] = d2a.data.0;
+        (
+            a,
+            -da_dv + temperature * density,
+            -da_dt,
+            molar_volume * molar_volume * d2a_dv2 + temperature,
+            -d2a_dvdt + density,
+        )
+    }
+
     /// calculates p, dp_drho, d2p_drho2
     fn p_dpdrho_d2pdrho2(
         &self,

diff --git a/crates/feos-core/src/phase_equilibria/mod.rs b/crates/feos-core/src/phase_equilibria/mod.rs
@@ -1,5 +1,5 @@
 use crate::equation_of_state::Residual;
-use crate::errors::{FeosError, FeosResult};
+use crate::errors::FeosResult;
 use crate::state::{DensityInitialization, State};
 use crate::{Contributions, ReferenceSystem};
 use nalgebra::allocator::Allocator;
@@ -44,12 +44,6 @@ pub struct PhaseEquilibrium<E, const P: usize, N: Dim = Dyn, D: DualNum<f64> + C
 where
     DefaultAllocator: Allocator<N>;
 
-// impl<E, const P: usize> Clone for PhaseEquilibrium<E, P> {
-//     fn clone(&self) -> Self {
-//         Self(self.0.clone())
-//     }
-// }
-
 impl<E: Residual, const P: usize> fmt::Display for PhaseEquilibrium<E, P> {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         for (i, s) in self.0.iter().enumerate() {
@@ -224,14 +218,6 @@ impl<E: Residual<N>, N: Dim> PhaseEquilibrium<E, 2, N>
 where
     DefaultAllocator: Allocator<N>,
 {
-    pub(super) fn check_trivial_solution(self) -> FeosResult<Self> {
-        if Self::is_trivial_solution(self.vapor(), self.liquid()) {
-            Err(FeosError::TrivialSolution)
-        } else {
-            Ok(self)
-        }
-    }
-
     /// Check if the two states form a trivial solution
     pub fn is_trivial_solution(state1: &State<E, N>, state2: &State<E, N>) -> bool {
         let rho1 = state1.molefracs.clone() * state1.density.into_reduced();