prehner · Sep 28, 2026 · Sep 18, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -9,13 +9,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added `PoreSpecification` enum to specify the state of the fluid in the pore (currently chemical potential or moles). [#371](https://github.com/feos-org/feos/pull/371)
 - Added the `Grid` and `PoreProfile` classes as an entry points to DFT in various coordinate systems in Python. [#376](https://github.com/feos-org/feos/pull/376)
 - `PoreProfile.solver_log` now stores a timing profile of the different solver steps. [#385](https://github.com/feos-org/feos/pull/385)
+- Added `Interface::curved` to calculate cylindrical and spherical interfaces from a planar interface. [#387](https://github.com/feos-org/feos/pull/387)
 
 ### Changed
 - Changed data type of initial temperatures or pressure for phase equilibrium calculations (`TemperatureOrPressure::Other`) from `D` to `f64`. [#369](https://github.com/feos-org/feos/pull/369)
 - Reworked DFT solution algorithms slightly for the cases in which additional specifications are given. [#371](https://github.com/feos-org/feos/pull/371)
 - External potentials are passed and returned as quantities (energies) instead of reduced units. [#372](https://github.com/feos-org/feos/pull/372)
 - Merged the `Adsorption1D` and `Adsorption3D` classes in Python into `Adsorption` using dynamically dimensioned arrays. [#376](https://github.com/feos-org/feos/pull/376)
 - If the `rayon` feature is activated (default for the released Python wheels), the evaluation of the helmholtz energy density and its derivatives will now be run in parallel. The number of CPUs can be set via the global rayon threadpool (see changelog for `0.9.5`) [#385](https://github.com/feos-org/feos/pull/385)
+- Renamed `PlanarInterface` to `Interface` and its constructors `from_tanh` and `from_pdgt` to `planar_from_tanh` and `planar_from_pdgt`. [#387](https://github.com/feos-org/feos/pull/387)
+- Replaced `PlanarInterface.from_density_profile` in Python with the `Interface(grid, vle, density)` constructor in analogy to `PoreProfile`. [#387](https://github.com/feos-org/feos/pull/387)
+- `drho_dp`, `dn_dp`, `drho_dt`, and `dn_dt` of DFT profiles now take the bulk state as an argument (methods instead of properties in Python). [#387](https://github.com/feos-org/feos/pull/387)
 
 ### Removed
 - Removed the `DFTSpecification` trait in favor of only using the `DFTSpecification` enum (renamed from `DFTSpecifications`). [#371](https://github.com/feos-org/feos/pull/371)

diff --git a/crates/feos-dft/src/adsorption/mod.rs b/crates/feos-dft/src/adsorption/mod.rs
@@ -353,17 +353,16 @@ where
     pub fn pressure(&self) -> Pressure<Array1<f64>> {
         Pressure::from_shape_fn(self.profiles.len(), |i| match &self.profiles[i] {
             Ok(p) => {
-                if p.profile.bulk.eos.components() > 1
-                    && !p.profile.bulk.is_stable(SolverOptions::default()).unwrap()
+                if p.profile.functional.components() > 1
+                    && !p.bulk.is_stable(SolverOptions::default()).unwrap()
                 {
-                    p.profile
-                        .bulk
+                    p.bulk
                         .tp_flash(None, SolverOptions::default(), None)
                         .unwrap()
                         .vapor()
                         .pressure(Contributions::Total)
                 } else {
-                    p.profile.bulk.pressure(Contributions::Total)
+                    p.bulk.pressure(Contributions::Total)
                 }
             }
             Err(_) => Pressure::from_reduced(f64::NAN),

diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs
@@ -3,7 +3,7 @@ use crate::convolver::ConvolverFFT;
 use crate::functional::{HelmholtzEnergyFunctional, HelmholtzEnergyFunctionalDyn, MoleculeShape};
 use crate::functional_contribution::FunctionalContribution;
 use crate::geometry::{Axis, Geometry, Grid};
-use crate::profile::DFTProfile;
+use crate::profile::{DFTProfile, InitialDensity};
 use crate::solver::DFTSolver;
 use crate::{DFTSpecification, WeightFunctionInfo};
 use feos_core::{Contributions, FeosResult, ReferenceSystem, ResidualDyn, State, StateHD};
@@ -106,6 +106,7 @@ pub enum PoreSpecification {
 #[derive(Clone)]
 pub struct PoreProfile<D: Dimension, F> {
     pub profile: DFTProfile<D, F>,
+    pub bulk: State<F>,
     pub grand_potential: Option<Energy>,
     pub interfacial_tension: Option<Energy>,
 }
@@ -123,7 +124,18 @@ where
         density: Option<&Density<Array<f64, D::Larger>>>,
         specification: PoreSpecification,
     ) -> Self {
-        let mut profile = DFTProfile::new(grid, bulk, Some(external_potential), density);
+        let density = match density {
+            Some(density) => InitialDensity::Density(density.clone()),
+            None => InitialDensity::Bulk(bulk.partial_density()),
+        };
+        let mut profile = DFTProfile::new(
+            grid,
+            &bulk.eos,
+            bulk.temperature,
+            Some(external_potential),
+            density,
+            DFTSpecification::from_state(bulk),
+        );
 
         // fix the number of particles
         match specification {
@@ -136,22 +148,23 @@ where
 
         Self {
             profile,
+            bulk: bulk.clone(),
             grand_potential: None,
             interfacial_tension: None,
         }
     }
 
     pub fn solve_inplace(&mut self, solver: Option<&DFTSolver>, debug: bool) -> FeosResult<()> {
         // Solve the profile
-        self.profile.solve(solver, debug)?;
+        self.profile.solve([&mut self.bulk], solver, debug)?;
 
         // calculate grand potential density
         let omega = self.profile.grand_potential()?;
         self.grand_potential = Some(omega);
 
         // calculate interfacial tension
         self.interfacial_tension =
-            Some(omega + self.profile.bulk.pressure(Contributions::Total) * self.profile.volume());
+            Some(omega + self.bulk.pressure(Contributions::Total) * self.profile.volume());
 
         Ok(())
     }
@@ -163,7 +176,7 @@ where
 
     pub fn update_bulk(mut self, bulk: &State<F>) -> Self {
         self.profile.specification = DFTSpecification::from_state(bulk);
-        self.profile.bulk = bulk.clone();
+        self.bulk = bulk.clone();
         self.grand_potential = None;
         self.interfacial_tension = None;
         self
@@ -174,7 +187,7 @@ where
     ) -> FeosResult<MolarEnergy<DVector<f64>>> {
         let a = self.profile.dn_dmu()?;
         let a_unit = a.get2(0, 0);
-        let b = -self.profile.bulk.temperature * self.profile.dn_dt()?;
+        let b = -self.profile.temperature * self.profile.dn_dt(&self.bulk)?;
         let b_unit = b.get(0);
 
         let h_ads = LU::new((a / a_unit).into_value())?.solve(&(b / b_unit).into_value());
@@ -184,67 +197,68 @@ where
     pub fn enthalpy_of_adsorption(&mut self) -> FeosResult<MolarEnergy> {
         Ok(self
             .partial_molar_enthalpy_of_adsorption()?
-            .dot(&Dimensionless::new(self.profile.bulk.molefracs.clone())))
+            .dot(&Dimensionless::new(self.bulk.molefracs.clone())))
     }
 
     fn _henry_coefficients<N: DualNum<Primitive = f64> + Copy + DctNum>(
         &self,
         temperature: N,
     ) -> DVector<N> {
-        if self.profile.bulk.eos.m().iter().any(|&m| m != 1.0) {
+        if self.profile.functional.m().iter().any(|&m| m != 1.0) {
             panic!(
                 "Henry coefficients can only be calculated for spherical and heterosegmented molecules!"
             )
         };
         let pot = (self.profile.external_potential.mapv(N::from)
-            * self.profile.bulk.temperature.to_reduced())
+            * self.profile.temperature.to_reduced())
         .mapv(|v| v / temperature);
         let exp_pot = pot.mapv(|v| (-v).exp());
-        let functional_contributions = self.profile.bulk.eos.contributions();
+        let functional_contributions = self.profile.functional.contributions();
         let weight_functions: Vec<WeightFunctionInfo<N>> = functional_contributions
             .into_iter()
             .map(|c| c.weight_functions(temperature))
             .collect();
         let convolver = ConvolverFFT::<_, D>::plan(&self.profile.grid, &weight_functions);
-        let bonds = self
-            .profile
-            .bulk
-            .eos
-            .bond_integrals(temperature, &exp_pot, convolver.as_ref());
+        let bonds =
+            self.profile
+                .functional
+                .bond_integrals(temperature, &exp_pot, convolver.as_ref());
         self.profile.integrate_reduced_segments(&(exp_pot * bonds))
     }
 
     pub fn henry_coefficients(&self) -> HenryCoefficient<DVector<f64>> {
-        let t = self.profile.bulk.temperature.to_reduced();
-        Volume::from_reduced(self._henry_coefficients(t)) / (RGAS * self.profile.bulk.temperature)
+        let t = self.profile.temperature.to_reduced();
+        Volume::from_reduced(self._henry_coefficients(t)) / (RGAS * self.profile.temperature)
     }
 
     pub fn ideal_gas_enthalpy_of_adsorption(&self) -> MolarEnergy<DVector<f64>> {
-        let t = Dual64::from(self.profile.bulk.temperature.to_reduced()).derivative();
+        let t = Dual64::from(self.profile.temperature.to_reduced()).derivative();
         let h_dual = self._henry_coefficients(t);
         let h = h_dual.map(|h| h.re);
         let dh = h_dual.map(|h| h.eps);
-        let t = self.profile.bulk.temperature.to_reduced();
-        RGAS * self.profile.bulk.temperature
+        let t = self.profile.temperature.to_reduced();
+        RGAS * self.profile.temperature
             * Dimensionless::from_reduced((&h - t * dh).component_div(&h))
     }
 
     pub fn into_dyn(self) -> PoreProfile<IxDyn, F> {
         // initialize convolver
-        let t = self.profile.bulk.temperature.to_reduced();
-        let weight_functions = self.profile.bulk.eos.weight_functions(t);
+        let t = self.profile.temperature.to_reduced();
+        let weight_functions = self.profile.functional.weight_functions(t);
         let convolver = ConvolverFFT::plan(&self.profile.grid, &weight_functions);
 
         PoreProfile {
             profile: DFTProfile {
                 grid: self.profile.grid,
                 convolver,
+                functional: self.profile.functional,
+                temperature: self.profile.temperature,
                 density: self.profile.density.into_dyn(),
                 specification: self.profile.specification,
                 external_potential: self.profile.external_potential.into_dyn(),
-                bulk: self.profile.bulk,
                 solver_log: self.profile.solver_log,
             },
+            bulk: self.bulk,
             grand_potential: self.grand_potential,
             interfacial_tension: self.interfacial_tension,
         }

diff --git a/crates/feos-dft/src/geometry.rs b/crates/feos-dft/src/geometry.rs
@@ -179,6 +179,15 @@ pub struct Axis {
 }
 
 impl Axis {
+    /// Create a new axis for the given geometry.
+    pub fn new(points: usize, length: Length, geometry: Geometry) -> Self {
+        match geometry {
+            Geometry::Cartesian => Self::new_cartesian(points, length, None),
+            Geometry::Cylindrical => Self::new_polar(points, length),
+            Geometry::Spherical => Self::new_spherical(points, length),
+        }
+    }
+
     /// Create a new (equidistant) cartesian axis.
     ///
     /// The potential_offset is required to make sure that particles
@@ -260,6 +269,15 @@ impl Axis {
         }
     }
 
+    pub fn shift(&mut self, x0: Length) {
+        if !matches!(self.geometry, Geometry::Cartesian) {
+            panic!("Cannot shift cylindrical or spherical axes!");
+        }
+        let x0 = x0.into_reduced();
+        self.grid -= x0;
+        self.edges -= x0;
+    }
+
     /// Returns the total length of the axis.
     ///
     /// This includes the `potential_offset` and used e.g.