From d55e4495abd6dc5d28b0a2bf2209e19315f5304f Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 18 Sep 2026 15:03:11 +0200 Subject: [PATCH] Revive curved interfaces --- CHANGELOG.md | 4 + crates/feos-dft/src/adsorption/mod.rs | 9 +- crates/feos-dft/src/adsorption/pore.rs | 60 +++-- crates/feos-dft/src/geometry.rs | 18 ++ crates/feos-dft/src/interface/mod.rs | 197 ++++++++-------- .../src/interface/surface_tension_diagram.rs | 15 +- crates/feos-dft/src/lib.rs | 2 +- crates/feos-dft/src/profile/mod.rs | 214 +++++++++++------- crates/feos-dft/src/profile/properties.rs | 102 ++++----- .../src/solvation/pair_correlation.rs | 14 +- .../src/solvation/solvation_profile.rs | 8 +- crates/feos-dft/src/solver.rs | 37 +-- crates/feos/tests/gc_pcsaft/dft.rs | 8 +- crates/feos/tests/pcsaft/dft.rs | 24 +- docs/api/dft.md | 2 +- py-feos/src/dft/adsorption/pore.rs | 5 + py-feos/src/dft/interface/mod.rs | 138 ++++++----- .../dft/interface/surface_tension_diagram.rs | 16 +- py-feos/src/dft/mod.rs | 2 +- py-feos/src/dft/profile.rs | 27 +-- py-feos/src/dft/solvation.rs | 10 + py-feos/src/lib.rs | 2 +- 22 files changed, 517 insertions(+), 397 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 10aaae91b..e02d8970d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ 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) @@ -16,6 +17,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - 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 index fde4814fe..e16c6a4b6 100644 --- 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> { 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 index 7e58d4698..cdb75fcd4 100644 --- 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 { pub profile: DFTProfile, + pub bulk: State, pub grand_potential: Option, pub interfacial_tension: Option, } @@ -123,7 +124,18 @@ where density: Option<&Density>>, 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,6 +148,7 @@ where Self { profile, + bulk: bulk.clone(), grand_potential: None, interfacial_tension: None, } @@ -143,7 +156,7 @@ where 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()?; @@ -151,7 +164,7 @@ where // 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) -> 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>> { 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 { 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 + Copy + DctNum>( &self, temperature: N, ) -> DVector { - 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> = 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> { - 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> { - 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 { // 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 index cbdc88b36..bb48ff166 100644 --- 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. diff --git a/crates/feos-dft/src/interface/mod.rs b/crates/feos-dft/src/interface/mod.rs index c68349caa..9ff8b0bbc 100644 --- a/crates/feos-dft/src/interface/mod.rs +++ b/crates/feos-dft/src/interface/mod.rs @@ -1,12 +1,14 @@ //! Density profiles at planar interfaces and interfacial tensions. +use crate::DFTSpecification; use crate::functional::HelmholtzEnergyFunctional; -use crate::geometry::{Axis, Grid}; +use crate::geometry::{Axis, Geometry, Grid}; use crate::pdgt::PdgtFunctionalProperties; -use crate::profile::DFTProfile; +use crate::profile::{DFTProfile, InitialDensity}; use crate::solver::DFTSolver; use feos_core::{Contributions, FeosError, FeosResult, PhaseEquilibrium, ReferenceSystem}; use ndarray::{Array1, Array2, Axis as Axis_nd, Ix1, s}; use quantity::{Area, Density, Length, Moles, SurfaceTension, Temperature}; +use std::f64::consts::{FRAC_PI_3, PI}; mod surface_tension_diagram; pub use surface_tension_diagram::SurfaceTensionDiagram; @@ -16,32 +18,43 @@ const MIN_WIDTH: f64 = 100.0; /// Density profile and properties of a planar interface. #[derive(Clone)] -pub struct PlanarInterface { +pub struct Interface { + pub geometry: Geometry, pub profile: DFTProfile, pub vle: PhaseEquilibrium, pub surface_tension: Option, - pub equimolar_radius: Option, + pub surface_of_tension: Option, } -impl PlanarInterface { +impl Interface { pub fn solve_inplace(&mut self, solver: Option<&DFTSolver>, debug: bool) -> FeosResult<()> { // Solve the profile - self.profile.solve(solver, debug)?; + self.profile + .solve(self.vle.states.each_mut(), solver, debug)?; // postprocess - self.surface_tension = Some( - (self.profile.integrate( - &(self.profile.grand_potential_density()? - + self.vle.vapor().pressure(Contributions::Total)), - )) / Area::from_reduced(1.0), - ); - let delta_rho = self.vle.liquid().density - self.vle.vapor().density; - self.equimolar_radius = Some( - self.profile - .integrate(&(self.profile.density.sum_axis(Axis_nd(0)) - self.vle.vapor().density)) - / delta_rho - / Area::from_reduced(1.0), + let delta_p = self.vle.liquid().pressure(Contributions::Total) + - self.vle.vapor().pressure(Contributions::Total); + let delta_omega = self.profile.integrate( + &(self.profile.grand_potential_density()? + + self.vle.vapor().pressure(Contributions::Total)), ); + match self.geometry { + Geometry::Cartesian => { + self.surface_tension = Some(delta_omega / Area::from_reduced(1.0)); + } + Geometry::Cylindrical => { + self.surface_tension = + Some((delta_omega * delta_p / (PI * Length::from_reduced(1.0))).sqrt()); + self.surface_of_tension = + Some((delta_omega / (PI * Length::from_reduced(1.0) * delta_p)).sqrt()); + } + Geometry::Spherical => { + self.surface_tension = + Some((delta_omega * delta_p * delta_p / (16.0 * FRAC_PI_3)).cbrt()); + self.surface_of_tension = Some((delta_omega / (2.0 * FRAC_PI_3 * delta_p)).cbrt()); + } + }; Ok(()) } @@ -52,59 +65,71 @@ impl PlanarInterface { } } -impl PlanarInterface { - pub fn new(vle: &PhaseEquilibrium, n_grid: usize, l_grid: Length) -> Self { +impl Interface { + pub fn new(vle: &PhaseEquilibrium, axis: Axis, density: Density>) -> Self { // generate grid - let grid = Grid::Cartesian1(Axis::new_cartesian(n_grid, l_grid, None)); + let geometry = axis.geometry; + let grid = Grid::new_1d(axis); + + // initialize profile + let mut profile = DFTProfile::new( + grid, + &vle.vapor().eos, + vle.vapor().temperature, + None, + InitialDensity::Density(density), + DFTSpecification::from_state(vle.vapor()), + ); + profile.fix_total_moles(vle.vapor()); + + // shift the profile so that the equimolar surface is at x=0 + if matches!(geometry, Geometry::Cartesian) { + let indices = vle.vapor().eos.component_index(); + let rho_v = vle.vapor().partial_density(); + let excess = Density::from_shape_fn(profile.density.raw_dim(), |(i, z)| { + profile.density.get((i, z)) - rho_v.get(indices[i]) + }); + let delta_rho = vle.liquid().density - vle.vapor().density; + let ze = + (profile.integrate_segments(&excess).sum() / delta_rho) / Area::from_reduced(1.0); + profile.grid.axes_mut()[0].shift(ze); + } Self { - profile: DFTProfile::new(grid, vle.vapor(), None, None), + geometry, + profile, vle: vle.clone(), surface_tension: None, - equimolar_radius: None, + surface_of_tension: None, } } - pub fn from_tanh( + pub fn planar_from_tanh( vle: &PhaseEquilibrium, n_grid: usize, l_grid: Length, critical_temperature: Temperature, - fix_equimolar_surface: bool, ) -> Self { - let mut profile = Self::new(vle, n_grid, l_grid); - // calculate segment indices - let indices = &profile.profile.bulk.eos.component_index(); + let indices = &vle.vapor().eos.component_index(); + let ax = Axis::new_cartesian(n_grid, l_grid, None); // calculate density profile let z0 = 0.5 * l_grid.to_reduced(); let (z0, sign) = (z0.abs(), -z0.signum()); let reduced_temperature = (vle.vapor().temperature / critical_temperature).into_value(); - profile.profile.density = - Density::from_shape_fn(profile.profile.density.raw_dim(), |(i, z)| { - let rho_v = profile.vle.vapor().partial_density().get(indices[i]); - let rho_l = profile.vle.liquid().partial_density().get(indices[i]); - 0.5 * (rho_l - rho_v) - * (sign * (profile.profile.grid.grids()[0][z] - z0) / 3.0 - * (2.4728 - 2.3625 * reduced_temperature)) - .tanh() - + 0.5 * (rho_l + rho_v) - }); - - // specify specification - if fix_equimolar_surface { - profile.profile.fix_total_moles(); - } + let density = Density::from_shape_fn((indices.len(), ax.grid.len()), |(i, z)| { + let rho_v = vle.vapor().partial_density().get(indices[i]); + let rho_l = vle.liquid().partial_density().get(indices[i]); + 0.5 * (rho_l - rho_v) + * (sign * (ax.grid[z] - z0) / 3.0 * (2.4728 - 2.3625 * reduced_temperature)).tanh() + + 0.5 * (rho_l + rho_v) + }); - profile + Self::new(vle, ax, density) } - pub fn from_pdgt( - vle: &PhaseEquilibrium, - n_grid: usize, - fix_equimolar_surface: bool, - ) -> FeosResult { + pub fn planar_from_pdgt(vle: &PhaseEquilibrium, n_grid: usize) -> FeosResult { let dft = &vle.vapor().eos; if dft.component_index().len() != 1 { @@ -127,49 +152,38 @@ impl PlanarInterface { // create PlanarInterface let l_grid = Length::from_reduced(MIN_WIDTH).max(w_pdgt * RELATIVE_WIDTH); - let mut profile = Self::new(vle, n_grid, l_grid); // interpolate density profile from pDGT to DFT + let ax = Axis::new_cartesian(n_grid, l_grid, None); let r = l_grid * 0.5; - profile.profile.density = interp_symmetric( - vle, - z_pdgt, - rho_pdgt, - &profile.vle, - profile.profile.grid.grids()[0], - r, - )?; - - // specify specification - if fix_equimolar_surface { - profile.profile.fix_total_moles(); - } + let density = interp_symmetric(vle, z_pdgt, rho_pdgt, vle, &ax.grid, r)?; - Ok(profile) + Ok(Self::new(vle, ax, density)) } -} -impl PlanarInterface { - pub fn shift_equimolar_inplace(&mut self) { - let s = self.profile.density.shape(); - let m = &self.profile.bulk.eos.m(); - let mut rho_l = Density::from_reduced(0.0); - let mut rho_v = Density::from_reduced(0.0); - let mut rho = Density::zeros(s[1]); - for i in 0..s[0] { - rho_l += self.profile.density.get((i, 0)) * m[i]; - rho_v += self.profile.density.get((i, s[1] - 1)) * m[i]; - rho += &(&self.profile.density.index_axis(Axis_nd(0), i) * m[i]); + pub fn curved( + planar_interface: &Interface, + n_grid: usize, + l_grid: Length, + radius: Length, + geometry: Geometry, + ) -> FeosResult { + let ax = Axis::new(n_grid, l_grid, geometry); + let density = interp_symmetric( + &planar_interface.vle, + planar_interface.profile.axes()[0].clone(), + planar_interface.profile.density.clone(), + &planar_interface.vle, + &ax.grid, + radius, + )?; + let mut vle = planar_interface.vle.clone(); + if radius.is_sign_negative() { + let [vapor, liquid] = vle.states; + vle = PhaseEquilibrium::two_phase(liquid, vapor); } - let x = (rho - rho_v) / (rho_l - rho_v); - let ze = self.profile.grid.axes()[0].edges[0] + self.profile.integrate(&x).to_reduced(); - self.profile.grid.axes_mut()[0].grid -= ze; - } - - pub fn shift_equimolar(mut self) -> Self { - self.shift_equimolar_inplace(); - self + Ok(Self::new(&vle, ax, density)) } /// Relative adsorption of component `i' with respect to `j': \Gamma_i^(j) @@ -337,14 +351,16 @@ fn interp_symmetric( z: &Array1, radius: Length, ) -> FeosResult>> { + let indices = vle_pdgt.vapor().eos.component_index(); let reduced_density = Array2::from_shape_fn(rho_pdgt.raw_dim(), |(i, j)| { - ((rho_pdgt.get((i, j)) - vle_pdgt.vapor().partial_density().get(i)) - / (vle_pdgt.liquid().partial_density().get(i) - - vle_pdgt.vapor().partial_density().get(i))) + let c = indices[i]; + ((rho_pdgt.get((i, j)) - vle_pdgt.vapor().partial_density().get(c)) + / (vle_pdgt.liquid().partial_density().get(c) + - vle_pdgt.vapor().partial_density().get(c))) .into_value() - 0.5 }); - let segments = vle_pdgt.vapor().eos.component_index().len(); + let segments = indices.len(); let mut reduced_density = interp( &z_pdgt.to_reduced(), &reduced_density, @@ -366,9 +382,10 @@ fn interp_symmetric( Ok(Density::from_shape_fn( reduced_density.raw_dim(), |(i, j)| { + let c = indices[i]; reduced_density[(i, j)] - * (vle.liquid().partial_density().get(i) - vle.vapor().partial_density().get(i)) - + vle.vapor().partial_density().get(i) + * (vle.liquid().partial_density().get(c) - vle.vapor().partial_density().get(c)) + + vle.vapor().partial_density().get(c) }, )) } diff --git a/crates/feos-dft/src/interface/surface_tension_diagram.rs b/crates/feos-dft/src/interface/surface_tension_diagram.rs index cb32bf107..95d01dfe2 100644 --- a/crates/feos-dft/src/interface/surface_tension_diagram.rs +++ b/crates/feos-dft/src/interface/surface_tension_diagram.rs @@ -1,4 +1,4 @@ -use super::PlanarInterface; +use super::Interface; use crate::functional::HelmholtzEnergyFunctional; use crate::solver::DFTSolver; use feos_core::{PhaseEquilibrium, ReferenceSystem, StateVec}; @@ -9,7 +9,7 @@ const DEFAULT_GRID_POINTS: usize = 2048; /// Container structure for the efficient calculation of surface tension diagrams. pub struct SurfaceTensionDiagram { - pub profiles: Vec>, + pub profiles: Vec>, } // #[expect(clippy::ptr_arg)] @@ -20,32 +20,29 @@ impl SurfaceTensionDiagram { n_grid: Option, l_grid: Option, critical_temperature: Option, - fix_equimolar_surface: Option, solver: Option<&DFTSolver>, ) -> Self { let n_grid = n_grid.unwrap_or(DEFAULT_GRID_POINTS); - let mut profiles: Vec> = Vec::with_capacity(dia.len()); + let mut profiles: Vec> = Vec::with_capacity(dia.len()); for vle in dia.iter() { // check for a critical point let profile = if PhaseEquilibrium::is_trivial_solution(vle.vapor(), vle.liquid()) { - Ok(PlanarInterface::from_tanh( + Ok(Interface::planar_from_tanh( vle, 10, Length::from_reduced(100.0), Temperature::from_reduced(500.0), - fix_equimolar_surface.unwrap_or(false), )) } else { // initialize with pDGT for single segments and tanh for mixtures and segment DFT if vle.vapor().eos.component_index().len() == 1 { - PlanarInterface::from_pdgt(vle, n_grid, false) + Interface::planar_from_pdgt(vle, n_grid) } else { - Ok(PlanarInterface::from_tanh( + Ok(Interface::planar_from_tanh( vle, n_grid, l_grid.unwrap_or(Length::from_reduced(100.0)), critical_temperature.unwrap_or(Temperature::from_reduced(500.0)), - fix_equimolar_surface.unwrap_or(false), )) } .map(|mut profile| { diff --git a/crates/feos-dft/src/lib.rs b/crates/feos-dft/src/lib.rs index d9e93f99a..fc16b7105 100644 --- a/crates/feos-dft/src/lib.rs +++ b/crates/feos-dft/src/lib.rs @@ -19,6 +19,6 @@ pub use functional::{HelmholtzEnergyFunctional, HelmholtzEnergyFunctionalDyn, Mo pub use functional_contribution::FunctionalContribution; pub use geometry::{Axis, Geometry, Grid}; pub use pdgt::PdgtFunctionalProperties; -pub use profile::{DFTProfile, DFTSpecification}; +pub use profile::{DFTProfile, DFTSpecification, InitialDensity}; pub use solver::{DFTSolver, DFTSolverLog}; pub use weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape}; diff --git a/crates/feos-dft/src/profile/mod.rs b/crates/feos-dft/src/profile/mod.rs index 33726a661..70f6b31a9 100644 --- a/crates/feos-dft/src/profile/mod.rs +++ b/crates/feos-dft/src/profile/mod.rs @@ -6,7 +6,7 @@ use feos_core::{FeosError, FeosResult, ReferenceSystem, State}; use nalgebra::{DVector, Dyn, U1}; use ndarray::{Array, Array1, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix0, arr1}; use num_dual::DualNum; -use quantity::{_Volume, Density, Energy, Entropy, Length, Moles, Quantity, Volume}; +use quantity::{_Volume, Density, Energy, Entropy, Length, Moles, Quantity, Temperature, Volume}; use std::ops::{Add, MulAssign}; use std::sync::Arc; @@ -84,19 +84,44 @@ impl DFTSpecification { } } +/// Initial condition for the density profile of a [DFTProfile]. +pub enum InitialDensity { + /// Use the provided density profile. + Density(Density>), + /// Initialize the density profile from the given bulk partial densities, + /// accounting for the external potential and bond integrals. + Bulk(Density>), +} + +/// Intermediate results of the evaluation of the Euler-Lagrange equation. +pub(crate) struct Projection { + /// $e^{-\frac{\beta}{m_\alpha}\left(\hat F_{\rho_\alpha}^\mathrm{res}+V_\alpha^\mathrm{ext}\right)}$ + pub exp_dfdrho: Array, + /// Product of the bond integrals $\prod_{\alpha'}I_{\alpha\alpha'}$ + pub bonds: Array, + /// Configurational integrals $z_\alpha$ + pub z: Array1, + /// Fugacities $f_\alpha$ determined from the specification + pub fugacity: Array1, + /// Projected density profile $f_\alpha e^{\ldots}\prod_{\alpha'}I_{\alpha\alpha'}$ + pub density: Array, +} + #[derive(Clone)] /// A one-, two-, or three-dimensional density profile. pub struct DFTProfile { pub grid: Grid, pub convolver: Arc>, + pub functional: F, + pub temperature: Temperature, pub density: Density>, pub specification: DFTSpecification, pub external_potential: Array, - pub bulk: State, pub solver_log: DFTSolverLog, } impl DFTProfile { + /// Return the grid points of every axis. pub fn axes(&self) -> Vec>> { self.grid .grids() @@ -106,6 +131,7 @@ impl DFTProfile { .collect() } + /// Return the cell edges of every axis. pub fn edges(&self) -> Vec>> { self.grid .axes() @@ -122,25 +148,24 @@ where { /// Create a new density profile. /// - /// If no external potential is specified, it is set to 0. The density is - /// initialized based on the bulk state and the external potential. The - /// specification is set to `ChemicalPotential` and needs to be overriden - /// after this call if something else is required. + /// If no external potential is specified, it is set to 0. pub fn new( grid: Grid, - bulk: &State, + functional: &F, + temperature: Temperature, external_potential: Option<&Energy>>, - density: Option<&Density>>, + density: InitialDensity, + specification: DFTSpecification, ) -> Self { // initialize convolver - let t = bulk.temperature.to_reduced(); - let weight_functions = bulk.eos.weight_functions(t); + let t = temperature.to_reduced(); + let weight_functions = functional.weight_functions(t); let convolver = ConvolverFFT::plan(&grid, &weight_functions); // initialize external potential let external_potential = external_potential.map_or_else( || { - let mut n_grid = vec![bulk.eos.component_index().len()]; + let mut n_grid = vec![functional.component_index().len()]; grid.axes() .iter() .for_each(|&ax| n_grid.push(ax.grid.len())); @@ -158,34 +183,55 @@ where ); // initialize density - let density = density.map_or_else( - || { + let density = match density { + InitialDensity::Density(density) => density, + InitialDensity::Bulk(partial_density) => { let exp_dfdrho = (-&external_potential).mapv(f64::exp); - let mut bonds = bulk.eos.bond_integrals(t, &exp_dfdrho, convolver.as_ref()); + let mut bonds = functional.bond_integrals(t, &exp_dfdrho, convolver.as_ref()); bonds *= &exp_dfdrho; let mut density = Array::zeros(external_potential.raw_dim()); - let bulk_density = bulk.partial_density().into_reduced(); - for (s, &c) in bulk.eos.component_index().iter().enumerate() { + let bulk_density = partial_density.into_reduced(); + for (s, &c) in functional.component_index().iter().enumerate() { density.index_axis_mut(Axis_nd(0), s).assign( &(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]), ); } Density::from_reduced(density) - }, - Clone::clone, - ); + } + }; Self { grid, convolver, + functional: functional.clone(), + temperature, density, - specification: DFTSpecification::from_state(bulk), + specification, external_potential, - bulk: bulk.clone(), solver_log: DFTSolverLog::new(), } } + /// Create a new density profile in equilibrium with a bulk state. + /// + /// The density is initialized based on the bulk state and the external + /// potential. The specification is set to `ChemicalPotential` and needs + /// to be overriden after this call if something else is required. + pub fn from_bulk( + grid: Grid, + bulk: &State, + external_potential: Option<&Energy>>, + ) -> Self { + Self::new( + grid, + &bulk.eos, + bulk.temperature, + external_potential, + InitialDensity::Bulk(bulk.partial_density()), + DFTSpecification::from_state(bulk), + ) + } + /// Set a constraint to fix the number of particles of every component based on /// the current density profile. pub fn fix_moles(&mut self) { @@ -195,11 +241,12 @@ where /// Set a constraint to fix the total number of particles based on /// the current density profile. - pub fn fix_total_moles(&mut self) { + /// + /// The ratios of the fugacities are fixed to those of the bulk state. + pub fn fix_total_moles(&mut self, bulk: &State) { let rho = self.density.to_reduced(); let moles = self.grid.integrate_reduced_comp(&rho).sum(); - let DFTSpecification::ChemicalPotential(fugacity) = - DFTSpecification::from_state(&self.bulk) + let DFTSpecification::ChemicalPotential(fugacity) = DFTSpecification::from_state(bulk) else { unreachable!() }; @@ -208,7 +255,7 @@ where /// Return the external potential in SI units. pub fn external_potential(&self) -> Energy> { - Entropy::from_reduced(self.external_potential.clone()) * self.bulk.temperature + Entropy::from_reduced(self.external_potential.clone()) * self.temperature } } @@ -224,8 +271,8 @@ where profile: &ArrayBase, ) -> DVector { let integral = self.grid.integrate_reduced_comp(profile); - let mut integral_comp = DVector::zeros(self.bulk.eos.components()); - for (i, &j) in self.bulk.eos.component_index().iter().enumerate() { + let mut integral_comp = DVector::zeros(self.functional.components()); + for (i, &j) in self.functional.component_index().iter().enumerate() { integral_comp[j] = integral[i]; } integral_comp @@ -257,7 +304,7 @@ where Volume::from_reduced(functional_determinant) * value.sum() } - /// Integrate each component individually. + /// Integrate each segment individually. pub fn integrate_comp, U>( &self, profile: &Quantity, U>, @@ -270,7 +317,11 @@ where }) } - /// Integrate each segment individually and aggregate to components. + /// Integrate each segment individually and assign the result to the + /// corresponding component. + /// + /// Only valid for profiles for which all segments of a molecule have the + /// same integral (e.g., density profiles). pub fn integrate_segments, U>( &self, profile: &Quantity, U>, @@ -279,8 +330,8 @@ where _Volume: Add, { let integral = self.integrate_comp(profile); - let mut integral_comp = Quantity::new(DVector::zeros(self.bulk.eos.components())); - for (i, &j) in self.bulk.eos.component_index().iter().enumerate() { + let mut integral_comp = Quantity::new(DVector::zeros(self.functional.components())); + for (i, &j) in self.functional.component_index().iter().enumerate() { integral_comp.set(j, integral.get(i)); } integral_comp @@ -303,33 +354,26 @@ where ::Larger: Dimension, F: HelmholtzEnergyFunctional, { + /// Calculate the (reduced) weighted densities of every functional contribution. pub fn weighted_densities(&self) -> FeosResult>> { Ok(self .convolver .weighted_densities(&self.density.to_reduced())) } + /// Calculate the residual of the Euler-Lagrange equation and its norm. pub fn residual(&mut self, log: bool) -> FeosResult<(Array, f64)> { let density = self.density.to_reduced(); - let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&density, log)?; + let (res, res_norm, _) = self.euler_lagrange_equation(&density, log)?; Ok((res, res_norm)) } - #[expect(clippy::type_complexity)] - pub(crate) fn fugacity( - &mut self, - density: &Array, - ) -> FeosResult<( - Array, - Array1, - Array, - Array1, - )> { + pub(crate) fn project(&mut self, density: &Array) -> FeosResult> { // calculate reduced temperature - let temperature = self.bulk.temperature.to_reduced(); + let temperature = self.temperature.to_reduced(); // calculate intrinsic functional derivative - let (_, mut dfdrho) = self.bulk.eos.functional_derivative( + let (_, mut dfdrho) = self.functional.functional_derivative( temperature, density, self.convolver.as_ref(), @@ -341,16 +385,15 @@ where dfdrho .outer_iter_mut() - .zip(self.bulk.eos.m().iter()) + .zip(self.functional.m().iter()) .for_each(|(mut df, &m)| df /= m); // calculate bond integrals let exp_dfdrho = dfdrho.mapv(|x| (-x).exp()); - let bonds = self - .bulk - .eos - .bond_integrals(temperature, &exp_dfdrho, self.convolver.as_ref()); - let mut rho_projected = &exp_dfdrho * bonds; + let bonds = + self.functional + .bond_integrals(temperature, &exp_dfdrho, self.convolver.as_ref()); + let mut rho_projected = &exp_dfdrho * &bonds; let z = self.grid.integrate_reduced_comp(&rho_projected); // calculate fugacity based on the given specification @@ -364,7 +407,13 @@ where x *= f; }); - Ok((exp_dfdrho, z, rho_projected, fugacity)) + Ok(Projection { + exp_dfdrho, + bonds, + z, + fugacity, + density: rho_projected, + }) } #[expect(clippy::type_complexity)] @@ -372,21 +421,16 @@ where &mut self, density: &Array, log: bool, - ) -> FeosResult<( - Array, - f64, - Array, - Array1, - Array, - )> { + ) -> FeosResult<(Array, f64, Projection)> { // calculate functional derivatives and fugacity - let (exp_dfdrho, z, rho_projected, _) = self.fugacity(density)?; + let projection = self.project(density)?; + let rho_projected = &projection.density; // calculate residual let mut res = if log { rho_projected.mapv(f64::ln) - density.mapv(f64::ln) } else { - &rho_projected - density + rho_projected - density }; // set residual to 0 where external potentials are overwhelming @@ -397,16 +441,26 @@ where // calculate the norm of the residual let res_norm = - (density - &rho_projected).mapv(|x| x * x).sum().sqrt() / (res.len() as f64).sqrt(); + (density - rho_projected).mapv(|x| x * x).sum().sqrt() / (res.len() as f64).sqrt(); if res_norm.is_finite() { - Ok((res, res_norm, exp_dfdrho, z, rho_projected)) + Ok((res, res_norm, projection)) } else { Err(FeosError::IterationFailed("Euler-Lagrange equation".into())) } } - pub fn solve(&mut self, solver: Option<&DFTSolver>, debug: bool) -> FeosResult<()> { + /// Solve the Euler-Lagrange equation. + /// + /// For specifications other than `ChemicalPotential`, the provided bulk + /// states are updated to be in equilibrium with the converged profile. + /// If `debug` is true, no error is returned if the solver does not converge. + pub fn solve( + &mut self, + bulk_states: [&mut State; B], + solver: Option<&DFTSolver>, + debug: bool, + ) -> FeosResult<()> { // unwrap solver let solver = solver.cloned().unwrap_or_default(); @@ -418,22 +472,24 @@ where // Update bulk state if !matches!(self.specification, DFTSpecification::ChemicalPotential(_)) { - // solve a bulk profile with the Newton solver - let mut bulk_profile = DFTProfile::::new(Grid::Bulk, &self.bulk, None, None); - let (_, _, _, fugacity) = self.fugacity(&density)?; - bulk_profile.specification = DFTSpecification::ChemicalPotential(fugacity); - let solver = DFTSolver::new(None).newton(None, None, None, None); - bulk_profile.solve(Some(&solver), false)?; - - // create the state based on the results from the bulk profile - let component_index = self.bulk.eos.component_index(); - let mut partial_density = self.bulk.partial_density(); - bulk_profile - .density - .into_iter() - .enumerate() - .for_each(|(i, r)| partial_density.set(component_index[i], r)); - self.bulk = State::new_density(&self.bulk.eos, self.bulk.temperature, partial_density)?; + let fugacity = self.project(&density)?.fugacity; + for bulk in bulk_states { + // solve a bulk profile with the Newton solver + let mut bulk_profile = DFTProfile::::from_bulk(Grid::Bulk, bulk, None); + bulk_profile.specification = DFTSpecification::ChemicalPotential(fugacity.clone()); + let solver = DFTSolver::new(None).newton(None, None, None, None); + bulk_profile.solve([], Some(&solver), false)?; + + // create the state based on the results from the bulk profile + let component_index = self.functional.component_index(); + let mut partial_density = bulk.partial_density(); + bulk_profile + .density + .into_iter() + .enumerate() + .for_each(|(i, r)| partial_density.set(component_index[i], r)); + *bulk = State::new_density(&self.functional, self.temperature, partial_density)?; + } } // Update profile diff --git a/crates/feos-dft/src/profile/properties.rs b/crates/feos-dft/src/profile/properties.rs index fcdb18f9b..2162e62ec 100644 --- a/crates/feos-dft/src/profile/properties.rs +++ b/crates/feos-dft/src/profile/properties.rs @@ -4,7 +4,7 @@ use crate::convolver::{Convolver, PeriodicConvolver}; use crate::functional_contribution::FunctionalContribution; use crate::solver::DFTSolverLog; use crate::{ConvolverFFT, HelmholtzEnergyFunctional, WeightFunctionInfo}; -use feos_core::{Contributions, FeosResult, ReferenceSystem, Total}; +use feos_core::{Contributions, FeosResult, ReferenceSystem, State, Total}; use nalgebra::{DMatrix, DVector}; use ndarray::{Array, Array1, Axis, Dimension, Ix0, RemoveAxis}; use num_dual::{Dual64, DualNum}; @@ -29,10 +29,10 @@ where /// Calculate the grand potential density $\omega$. pub fn grand_potential_density(&self) -> FeosResult>> { // Calculate residual Helmholtz energy density and functional derivative - let t = self.bulk.temperature.to_reduced(); + let t = self.temperature.to_reduced(); let rho = self.density.to_reduced(); let mut solver_log = DFTSolverLog::new(); - let (mut f, dfdrho) = self.bulk.eos.functional_derivative( + let (mut f, dfdrho) = self.functional.functional_derivative( t, &rho, self.convolver.as_ref(), @@ -43,12 +43,12 @@ where for ((rho, dfdrho), &m) in rho .outer_iter() .zip(dfdrho.outer_iter()) - .zip(self.bulk.eos.m().iter()) + .zip(self.functional.m().iter()) { f -= &((&dfdrho + m) * &rho); } - let bond_lengths = self.bulk.eos.bond_lengths(t); + let bond_lengths = self.functional.bond_lengths(t); for segment in bond_lengths.node_indices() { let n = bond_lengths.neighbors(segment).count(); f += &(&rho.index_axis(Axis(0), segment.index()) * (0.5 * n as f64)); @@ -62,12 +62,13 @@ where Ok(self.integrate(&self.grand_potential_density()?)) } - /// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{F}}{\delta\rho_i(\mathbf{r})}$. + /// Calculate the reduced residual Helmholtz energy density $\beta\hat f^\mathrm{res}$ and + /// functional derivative $\beta\hat F_{\rho_\alpha}^\mathrm{res}(\mathbf{r})$. #[expect(clippy::type_complexity)] pub fn functional_derivative(&self) -> FeosResult<(Array, Array)> { let mut solver_log = DFTSolverLog::new(); - self.bulk.eos.functional_derivative( - self.bulk.temperature.to_reduced(), + self.functional.functional_derivative( + self.temperature.to_reduced(), &self.density.to_reduced(), self.convolver.as_ref(), &mut solver_log, @@ -92,10 +93,9 @@ where { let density_dual = density.mapv(N::from); let weighted_densities = convolver.weighted_densities(&density_dual); - let functional_contributions = self.bulk.eos.contributions(); + let functional_contributions = self.functional.contributions(); let mut helmholtz_energy_density: Array = self - .bulk - .eos + .functional .ideal_chain_contribution() .helmholtz_energy_density(&density.mapv(N::from))?; for (c, wd) in functional_contributions.into_iter().zip(weighted_densities) { @@ -118,9 +118,9 @@ where /// Untested with heterosegmented functionals. pub fn residual_entropy_density(&self) -> FeosResult>> { // initialize convolver - let temperature = self.bulk.temperature.to_reduced(); + let temperature = self.temperature.to_reduced(); let temperature_dual = Dual64::from(temperature).derivative(); - let functional_contributions = self.bulk.eos.contributions(); + let functional_contributions = self.functional.contributions(); let weight_functions: Vec> = functional_contributions .into_iter() .map(|c| c.weight_functions(temperature_dual)) @@ -151,11 +151,10 @@ where let density_dual = density.mapv(Dual64::from); let temperature_dual = Dual64::from(temperature).derivative(); let weighted_densities = convolver.weighted_densities(&density_dual); - let functional_contributions = self.bulk.eos.contributions(); + let functional_contributions = self.functional.contributions(); let mut helmholtz_energy_density: Vec> = Vec::new(); helmholtz_energy_density.push( - self.bulk - .eos + self.functional .ideal_chain_contribution() .helmholtz_energy_density(&density.mapv(Dual64::from))?, ); @@ -190,7 +189,7 @@ where temperature: Dual64, density: &Array, ) -> Array { - let lambda = self.bulk.eos.lift_total().ln_lambda3(temperature); + let lambda = self.functional.lift_total().ln_lambda3(temperature); let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); for (i, rhoi) in density.outer_iter().enumerate() { phi += &rhoi.mapv(|rhoi| (lambda[i] + rhoi.ln() - 1.0) * rhoi); @@ -206,9 +205,9 @@ where contributions: Contributions, ) -> FeosResult>> { // initialize convolver - let temperature = self.bulk.temperature.to_reduced(); + let temperature = self.temperature.to_reduced(); let temperature_dual = Dual64::from(temperature).derivative(); - let functional_contributions = self.bulk.eos.contributions(); + let functional_contributions = self.functional.contributions(); let weight_functions: Vec> = functional_contributions .into_iter() .map(|c| c.weight_functions(temperature_dual)) @@ -256,9 +255,9 @@ where D::Larger: Dimension, { // initialize convolver - let temperature = self.bulk.temperature.to_reduced(); + let temperature = self.temperature.to_reduced(); let temperature_dual = Dual64::from(temperature).derivative(); - let functional_contributions = self.bulk.eos.contributions(); + let functional_contributions = self.functional.contributions(); let weight_functions: Vec> = functional_contributions .into_iter() .map(|c| c.weight_functions(temperature_dual)) @@ -306,13 +305,12 @@ where &mut self, lhs: &Array, ) -> FeosResult> { - let t = self.bulk.temperature.into_reduced(); + let t = self.temperature.into_reduced(); let rho = self.density.to_reduced(); let second_partial_derivatives = - self.bulk - .eos + self.functional .second_partial_derivatives(t, &rho, self.convolver.as_ref())?; - let (_, _, exp_dfdrho, _, _) = self.euler_lagrange_equation(&rho, false)?; + let exp_dfdrho = self.project(&rho)?.exp_dfdrho; let rhs = |x: &_| { let delta_functional_derivative = Self::delta_functional_derivative( @@ -322,9 +320,9 @@ where ); let mut xm = x.clone(); xm.outer_iter_mut() - .zip(self.bulk.eos.m().iter()) + .zip(self.functional.m().iter()) .for_each(|(mut x, &m)| x *= m); - let delta_i = self.bulk.eos.delta_bond_integrals( + let delta_i = self.functional.delta_bond_integrals( t, &exp_dfdrho, &delta_functional_derivative, @@ -338,12 +336,12 @@ where /// Return the partial derivatives of the density profiles w.r.t. the chemical potentials $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\mu_k}\right)_T$ pub fn drho_dmu(&mut self) -> FeosResult> { - let shape: Vec<_> = std::iter::once(&self.bulk.eos.components()) + let shape: Vec<_> = std::iter::once(&self.functional.components()) .chain(self.density.shape()) .copied() .collect(); let mut drho_dmu = Array::zeros(shape).into_dimensionality().unwrap(); - let component_index = self.bulk.eos.component_index().into_owned(); + let component_index = self.functional.component_index().into_owned(); for (k, mut d) in drho_dmu.outer_iter_mut().enumerate() { let mut lhs = self.density.to_reduced(); for (i, mut l) in lhs.outer_iter_mut().enumerate() { @@ -354,7 +352,7 @@ where d.assign(&self.density_derivative(&lhs)?); } Ok(Quantity::from_reduced( - drho_dmu / self.bulk.temperature.to_reduced(), + drho_dmu / self.temperature.to_reduced(), )) } @@ -371,42 +369,42 @@ where } /// Return the partial derivatives of the density profiles w.r.t. the bulk pressure at constant temperature and bulk composition $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial p}\right)_{T,\mathbf{x}}$ - pub fn drho_dp(&mut self) -> FeosResult> { + pub fn drho_dp(&mut self, bulk: &State) -> FeosResult> { let mut lhs = self.density.to_reduced(); - let v = self.bulk.partial_molar_volume().to_reduced(); + let v = bulk.partial_molar_volume().to_reduced(); for (mut l, &c) in lhs .outer_iter_mut() - .zip(self.bulk.eos.component_index().iter()) + .zip(self.functional.component_index().iter()) { l *= v[c]; } Ok(Quantity::from_reduced( - self.density_derivative(&lhs)? / self.bulk.temperature.to_reduced(), + self.density_derivative(&lhs)? / self.temperature.to_reduced(), )) } /// Return the partial derivatives of the number of moles w.r.t. the bulk pressure at constant temperature and bulk composition $\left(\frac{\partial N_i}{\partial p}\right)_{T,\mathbf{x}}$ - pub fn dn_dp(&mut self) -> FeosResult { - let drho_dp = self.drho_dp()?; + pub fn dn_dp(&mut self, bulk: &State) -> FeosResult { + let drho_dp = self.drho_dp(bulk)?; Ok(self.integrate_segments(&drho_dp)) } /// Return the partial derivatives of the density profiles w.r.t. the temperature at constant bulk pressure and composition $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial T}\right)_{p,\mathbf{x}}$ - pub fn drho_dt(&mut self) -> FeosResult> { + pub fn drho_dt(&mut self, bulk: &State) -> FeosResult> { let rho = self.density.to_reduced(); - let t = self.bulk.temperature.to_reduced(); + let t = self.temperature.to_reduced(); let rho_dual = rho.mapv(Dual64::from); let t_dual = Dual64::from(t).derivative(); let mut solver_log = DFTSolverLog::new(); // calculate intrinsic functional derivative - let functional_contributions = self.bulk.eos.contributions(); + let functional_contributions = self.functional.contributions(); let weight_functions: Vec> = functional_contributions .into_iter() .map(|c| c.weight_functions(t_dual)) .collect(); let convolver: Arc> = ConvolverFFT::plan(&self.grid, &weight_functions); - let (_, mut dfdrho) = self.bulk.eos.functional_derivative( + let (_, mut dfdrho) = self.functional.functional_derivative( t_dual, &rho_dual, convolver.as_ref(), @@ -417,17 +415,16 @@ where dfdrho += &(&self.external_potential * t).mapv(|v| Dual64::from(v) / t_dual); // calculate bulk functional derivative - let partial_density = self.bulk.partial_density().into_reduced(); + let partial_density = bulk.partial_density().into_reduced(); let rho_bulk: Array1<_> = self - .bulk - .eos + .functional .component_index() .iter() .map(|&i| partial_density[i]) .collect(); let rho_bulk_dual = rho_bulk.mapv(Dual64::from); let bulk_convolver = PeriodicConvolver::<_, Ix0>::new_0d(&weight_functions); - let (_, dfdrho_bulk) = self.bulk.eos.functional_derivative( + let (_, dfdrho_bulk) = self.functional.functional_derivative( t_dual, &rho_bulk_dual, bulk_convolver.as_ref(), @@ -436,7 +433,7 @@ where dfdrho .outer_iter_mut() .zip(dfdrho_bulk) - .zip(self.bulk.eos.m().iter()) + .zip(self.functional.m().iter()) .for_each(|((mut df, df_b), &m)| { df.map_inplace(|df| *df = (*df - df_b) / Dual64::from(m)); }); @@ -444,17 +441,14 @@ where // calculate bond integrals let exp_dfdrho = dfdrho.mapv(|x| (-x).exp()); let bonds = self - .bulk - .eos + .functional .bond_integrals(t_dual, &exp_dfdrho, convolver.as_ref()); // solve for drho_dt let mut lhs = (exp_dfdrho * bonds).mapv(|x| (-x.ln() * t_dual).eps); - let x = - (self.bulk.partial_molar_volume() * self.bulk.dp_dt(Contributions::Total)).to_reduced(); + let x = (bulk.partial_molar_volume() * bulk.dp_dt(Contributions::Total)).to_reduced(); let x: Array1<_> = self - .bulk - .eos + .functional .component_index() .iter() .map(|&i| x[i]) @@ -462,7 +456,7 @@ where lhs.outer_iter_mut() .zip(rho.outer_iter()) .zip(rho_bulk) - .zip(self.bulk.eos.m().iter()) + .zip(self.functional.m().iter()) .zip(x) .for_each(|((((mut lhs, rho), rho_b), &m), x)| { lhs += &(&rho / rho_b).mapv(f64::ln); @@ -480,8 +474,8 @@ where } /// Return the partial derivatives of the number of moles w.r.t. the temperature at constant bulk pressure and composition $\left(\frac{\partial N_i}{\partial T}\right)_{p,\mathbf{x}}$ - pub fn dn_dt(&mut self) -> FeosResult { - let drho_dt = self.drho_dt()?; + pub fn dn_dt(&mut self, bulk: &State) -> FeosResult { + let drho_dt = self.drho_dt(bulk)?; Ok(self.integrate_segments(&drho_dt)) } } diff --git a/crates/feos-dft/src/solvation/pair_correlation.rs b/crates/feos-dft/src/solvation/pair_correlation.rs index 20eeab73b..6aba6309c 100644 --- a/crates/feos-dft/src/solvation/pair_correlation.rs +++ b/crates/feos-dft/src/solvation/pair_correlation.rs @@ -23,6 +23,7 @@ impl, T: PairPotential> PairPotential for C { /// Density profile and properties of a test particle system. pub struct PairCorrelation { pub profile: DFTProfile, + pub bulk: State, pub pair_correlation_function: Option>, pub self_solvation_free_energy: Option, pub structure_factor: Option, @@ -41,7 +42,8 @@ impl PairCorrelation { let grid = Grid::Spherical(axis); Self { - profile: DFTProfile::new(grid, bulk, Some(&external_potential), None), + profile: DFTProfile::from_bulk(grid, bulk, Some(&external_potential)), + bulk: bulk.clone(), pair_correlation_function: None, self_solvation_free_energy: None, structure_factor: None, @@ -50,10 +52,10 @@ impl PairCorrelation { 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 pair correlation function - let partial_density = self.profile.bulk.partial_density(); + let partial_density = self.bulk.partial_density(); self.pair_correlation_function = Some(Array::from_shape_fn( self.profile.density.raw_dim(), |(i, j)| (self.profile.density.get((i, j)) / partial_density.get(i)).into_value(), @@ -61,14 +63,12 @@ impl PairCorrelation { // calculate self solvation free energy self.self_solvation_free_energy = Some(self.profile.integrate( - &(self.profile.grand_potential_density()? - + self.profile.bulk.pressure(Contributions::Total)), + &(self.profile.grand_potential_density()? + self.bulk.pressure(Contributions::Total)), )); // calculate structure factor self.structure_factor = Some( - (self.profile.total_moles() - self.profile.bulk.density * self.profile.volume()) - .to_reduced() + (self.profile.total_moles() - self.bulk.density * self.profile.volume()).to_reduced() + 1.0, ); diff --git a/crates/feos-dft/src/solvation/solvation_profile.rs b/crates/feos-dft/src/solvation/solvation_profile.rs index 08893ad24..268751287 100644 --- a/crates/feos-dft/src/solvation/solvation_profile.rs +++ b/crates/feos-dft/src/solvation/solvation_profile.rs @@ -13,6 +13,7 @@ const CUTOFF_RADIUS: f64 = 14.0; /// Density profile and properties of a solute in a inhomogeneous bulk fluid. pub struct SolvationProfile { pub profile: DFTProfile, + pub bulk: State, pub grand_potential: Option, pub solvation_free_energy: Option, } @@ -20,7 +21,7 @@ pub struct SolvationProfile { impl SolvationProfile { 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()?; @@ -28,7 +29,7 @@ impl SolvationProfile { // calculate solvation free energy self.solvation_free_energy = Some( - (omega + self.profile.bulk.pressure(Contributions::Total) * self.profile.volume()) + (omega + self.bulk.pressure(Contributions::Total) * self.profile.volume()) / Moles::from_reduced(1.0), ); @@ -90,7 +91,8 @@ impl SolvationProfile { let grid = Grid::Cartesian3(x, y, z); Ok(Self { - profile: DFTProfile::new(grid, bulk, Some(&external_potential), None), + profile: DFTProfile::from_bulk(grid, bulk, Some(&external_potential)), + bulk: bulk.clone(), grand_potential: None, solvation_free_energy: None, }) diff --git a/crates/feos-dft/src/solver.rs b/crates/feos-dft/src/solver.rs index 9348b90e6..067d78f95 100644 --- a/crates/feos-dft/src/solver.rs +++ b/crates/feos-dft/src/solver.rs @@ -312,7 +312,7 @@ where for k in 0..picard.max_iter { // calculate residual - let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&*rho, picard.log)?; + let (res, res_norm, _) = self.euler_lagrange_equation(&*rho, picard.log)?; self.solver_log.add_residual(solver, k, res_norm, verbosity); // check for convergence @@ -354,7 +354,7 @@ where } else { rho + alpha * delta_rho }; - let Ok((_, res2, _, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { + let Ok((_, res2, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { continue; }; if res2 > res0 { @@ -367,7 +367,7 @@ where } else { rho + 0.5 * alpha * delta_rho }; - let Ok((_, res1, _, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { + let Ok((_, res1, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { continue; }; @@ -417,7 +417,7 @@ where let m = resm.len() + 1; // calculate residual - let (res, res_norm, _, _, _) = self.euler_lagrange_equation(&*rho, anderson.log)?; + let (res, res_norm, _) = self.euler_lagrange_equation(&*rho, anderson.log)?; self.solver_log.add_residual(solver, k, res_norm, verbosity); // check for convergence @@ -477,8 +477,7 @@ where let solver = if newton.log { "Newton (log)" } else { "Newton" }; for k in 0..newton.max_iter { // calculate initial residual - let (res, res_norm, exp_dfdrho, z, rho_p) = - self.euler_lagrange_equation(rho, newton.log)?; + let (res, res_norm, projection) = self.euler_lagrange_equation(rho, newton.log)?; self.solver_log.add_residual(solver, k, res_norm, verbosity); // check convergence @@ -490,8 +489,8 @@ where let second_partial_derivatives = self.solver_log .time_function("second partial derivatives", || { - self.bulk.eos.second_partial_derivatives( - self.bulk.temperature.into_reduced(), + self.functional.second_partial_derivatives( + self.temperature.into_reduced(), rho, self.convolver.as_ref(), ) @@ -506,20 +505,20 @@ where ); delta_functional_derivative .outer_iter_mut() - .zip(self.bulk.eos.m().iter()) + .zip(self.functional.m().iter()) .for_each(|(mut q, &m)| q /= m); - let delta_i = self.bulk.eos.delta_bond_integrals( - self.bulk.temperature.into_reduced(), - &exp_dfdrho, + let delta_i = self.functional.delta_bond_integrals( + self.temperature.into_reduced(), + &projection.exp_dfdrho, &delta_functional_derivative, self.convolver.as_ref(), ); let mut delta_exp_dfdrho = delta_functional_derivative - delta_i; - let delta_z = -self - .grid - .integrate_reduced_comp(&(&delta_exp_dfdrho * &exp_dfdrho)); + let delta_z = -self.grid.integrate_reduced_comp( + &(&delta_exp_dfdrho * &projection.exp_dfdrho * &projection.bonds), + ); - let delta_fugacity = self.specification.delta_fugacity(&z, &delta_z); + let delta_fugacity = self.specification.delta_fugacity(&projection.z, &delta_z); delta_exp_dfdrho .outer_iter_mut() .zip(delta_fugacity.iter()) @@ -527,7 +526,11 @@ where z -= f; }); - let rho = if newton.log { &*rho } else { &rho_p }; + let rho = if newton.log { + &*rho + } else { + &projection.density + }; delta_rho + delta_exp_dfdrho * rho }; diff --git a/crates/feos/tests/gc_pcsaft/dft.rs b/crates/feos/tests/gc_pcsaft/dft.rs index 6046ab4c3..82bca14fb 100644 --- a/crates/feos/tests/gc_pcsaft/dft.rs +++ b/crates/feos/tests/gc_pcsaft/dft.rs @@ -5,7 +5,7 @@ use feos::gc_pcsaft::{GcPcSaft, GcPcSaftFunctional, GcPcSaftParameters}; use feos_core::parameter::{ChemicalRecord, Identifier, IdentifierOption, SegmentRecord}; use feos_core::{PhaseEquilibrium, State, Verbosity}; use feos_dft::adsorption::{ExternalPotential, Pore1D, PoreSpecification}; -use feos_dft::interface::PlanarInterface; +use feos_dft::interface::Interface; use feos_dft::{DFTSolver, Geometry}; use nalgebra::dvector; use quantity::*; @@ -160,7 +160,7 @@ fn test_dft() -> Result<(), Box> { let points = 2048; let tc = State::critical_point(&&func, (), None, None, Default::default())?.temperature; let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?; - let profile = PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(None)?; + let profile = Interface::planar_from_tanh(&vle, points, w, tc).solve(None)?; println!( "hetero {} {} {}", profile.surface_tension.unwrap(), @@ -205,7 +205,7 @@ fn test_dft_assoc() -> Result<(), Box> { let w = 100.0 * ANGSTROM; let points = 4096; let vle = PhaseEquilibrium::pure(&func, t, None, Default::default())?; - let profile = PlanarInterface::from_tanh(&vle, points, w, 600.0 * KELVIN, false).solve(None)?; + let profile = Interface::planar_from_tanh(&vle, points, w, 600.0 * KELVIN).solve(None)?; println!( "hetero {} {} {}", profile.surface_tension.unwrap(), @@ -254,6 +254,6 @@ fn test_dft_newton() -> Result<(), Box> { let solver = DFTSolver::new(Some(Verbosity::Iter)) .picard_iteration(None, Some(10), None, None) .newton(None, None, None, None); - PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(Some(&solver))?; + Interface::planar_from_tanh(&vle, points, w, tc).solve(Some(&solver))?; Ok(()) } diff --git a/crates/feos/tests/pcsaft/dft.rs b/crates/feos/tests/pcsaft/dft.rs index b04a40aa1..59bea1705 100644 --- a/crates/feos/tests/pcsaft/dft.rs +++ b/crates/feos/tests/pcsaft/dft.rs @@ -6,7 +6,7 @@ use feos::ideal_gas::{Joback, JobackParameters}; use feos::pcsaft::{PcSaft, PcSaftFunctional, PcSaftParameters}; use feos_core::parameter::IdentifierOption; use feos_core::{Contributions, EquationOfState, FeosResult, PhaseEquilibrium, State, Verbosity}; -use feos_dft::interface::PlanarInterface; +use feos_dft::interface::Interface; use feos_dft::{DFTSolver, PdgtFunctionalProperties}; use nalgebra::dvector; use ndarray::Axis; @@ -107,10 +107,9 @@ fn test_dft_propane() -> Result<(), Box> { let vle_pure = PhaseEquilibrium::pure(&&func_pure, t, None, Default::default())?; let vle_full = PhaseEquilibrium::pure(&&func_full, t, None, Default::default())?; let vle_full_vec = PhaseEquilibrium::pure(&&func_full_vec, t, None, Default::default())?; - let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc, false).solve(None)?; - let profile_full = PlanarInterface::from_tanh(&vle_full, points, w, tc, false).solve(None)?; - let profile_full_vec = - PlanarInterface::from_tanh(&vle_full_vec, points, w, tc, false).solve(None)?; + let profile_pure = Interface::planar_from_tanh(&vle_pure, points, w, tc).solve(None)?; + let profile_full = Interface::planar_from_tanh(&vle_full, points, w, tc).solve(None)?; + let profile_full_vec = Interface::planar_from_tanh(&vle_full_vec, points, w, tc).solve(None)?; let _ = (&func_pure).solve_pdgt(&vle_pure, 198, 0, None)?; println!( "pure {} {} {} {}", @@ -216,9 +215,9 @@ fn test_dft_propane_newton() -> Result<(), Box> { let tc = State::critical_point(&&func, (), None, None, Default::default())?.temperature; let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?; let solver = DFTSolver::new(Some(Verbosity::Iter)).newton(None, None, None, None); - PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(Some(&solver))?; + Interface::planar_from_tanh(&vle, points, w, tc).solve(Some(&solver))?; let solver = DFTSolver::new(Some(Verbosity::Iter)).newton(Some(true), None, None, None); - PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(Some(&solver))?; + Interface::planar_from_tanh(&vle, points, w, tc).solve(Some(&solver))?; Ok(()) } @@ -237,9 +236,8 @@ fn test_dft_water() -> Result<(), Box> { let tc = State::critical_point(&&func_pure, (), None, None, Default::default())?.temperature; let vle_pure = PhaseEquilibrium::pure(&&func_pure, t, None, Default::default())?; let vle_full_vec = PhaseEquilibrium::pure(&&func_full_vec, t, None, Default::default())?; - let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc, false).solve(None)?; - let profile_full_vec = - PlanarInterface::from_tanh(&vle_full_vec, points, w, tc, false).solve(None)?; + let profile_pure = Interface::planar_from_tanh(&vle_pure, points, w, tc).solve(None)?; + let profile_full_vec = Interface::planar_from_tanh(&vle_full_vec, points, w, tc).solve(None)?; println!( "pure {} {} {}", profile_pure.surface_tension.unwrap(), @@ -319,7 +317,7 @@ fn test_entropy_bulk_values() -> Result<(), Box> { )?); let func = EquationOfState::new(joback, PcSaftFunctional::new(params)); let vle = PhaseEquilibrium::pure(&&func, 350.0 * KELVIN, None, Default::default())?; - let profile = PlanarInterface::from_pdgt(&vle, 2048, false)?.solve(None)?; + let profile = Interface::planar_from_pdgt(&vle, 2048)?.solve(None)?; let s_res = profile.profile.entropy_density(Contributions::Residual)?; let s_tot = profile.profile.entropy_density(Contributions::Total)?; println!( @@ -356,7 +354,7 @@ fn test_entropy_bulk_values() -> Result<(), Box> { s_res.get(2047), profile.vle.vapor().molar_entropy(Contributions::Residual) / profile.vle.vapor().molar_volume, - max_relative = 1e-8, + max_relative = 1e-6, ); assert_relative_eq!( s_tot.get(0), @@ -367,7 +365,7 @@ fn test_entropy_bulk_values() -> Result<(), Box> { assert_relative_eq!( s_tot.get(2047), profile.vle.vapor().molar_entropy(Contributions::Total) / profile.vle.vapor().molar_volume, - max_relative = 1e-8, + max_relative = 1e-6, ); Ok(()) } diff --git a/docs/api/dft.md b/docs/api/dft.md index 4505caa29..62faeec6c 100644 --- a/docs/api/dft.md +++ b/docs/api/dft.md @@ -31,7 +31,7 @@ Implementations of Helmholtz energy functionals for DFT. .. autosummary:: :toctree: generated/ - PlanarInterface + Interface SurfaceTensionDiagram ``` diff --git a/py-feos/src/dft/adsorption/pore.rs b/py-feos/src/dft/adsorption/pore.rs index 15eaa4120..2205f8ca4 100644 --- a/py-feos/src/dft/adsorption/pore.rs +++ b/py-feos/src/dft/adsorption/pore.rs @@ -164,6 +164,11 @@ impl PyPoreProfile { )) } + #[getter] + fn get_bulk(&self) -> PyState { + PyState(self.0.bulk.clone()) + } + #[getter] fn get_grand_potential(&self) -> Option { self.0.grand_potential diff --git a/py-feos/src/dft/interface/mod.rs b/py-feos/src/dft/interface/mod.rs index dd628605f..dc7ef8010 100644 --- a/py-feos/src/dft/interface/mod.rs +++ b/py-feos/src/dft/interface/mod.rs @@ -1,3 +1,5 @@ +use super::PyGeometry; +use super::adsorption::PyGrid; use super::profile::impl_profile; use super::{PyDFTSolver, PyDFTSolverLog}; use crate::error::PyFeosError; @@ -5,8 +7,9 @@ use crate::ideal_gas::IdealGasModel; use crate::phase_equilibria::PyPhaseEquilibrium; use crate::residual::ResidualModel; use crate::state::{PyContributions, PyState}; -use feos_core::EquationOfState; -use feos_dft::interface::PlanarInterface; +use feos_core::{EquationOfState, FeosError}; +use feos_dft::Grid; +use feos_dft::interface::Interface; use nalgebra::{DMatrix, DVector}; use ndarray::*; use numpy::*; @@ -17,16 +20,30 @@ use std::sync::Arc; mod surface_tension_diagram; pub use surface_tension_diagram::PySurfaceTensionDiagram; -/// A one-dimensional density profile of a vapor-liquid or liquid-liquid interface. -#[pyclass(name = "PlanarInterface")] -pub struct PyPlanarInterface( - PlanarInterface, ResidualModel>>>, -); +/// A one-dimensional density profile of a planar, cylindrical or spherical +/// vapor-liquid or liquid-liquid interface. +/// +/// Parameters +/// ---------- +/// grid : Grid +/// The grid on which the density is calculated +/// (`Grid.cartesian_1d`, `Grid.polar`, or `Grid.spherical`). +/// vle : PhaseEquilibrium +/// The bulk phase equilibrium. +/// density : SIArray2 +/// The initial density profile. +/// +/// Returns +/// ------- +/// Interface +/// +#[pyclass(name = "Interface")] +pub struct PyInterface(Interface, ResidualModel>>>); -impl_profile!(PyPlanarInterface); +impl_profile!(PyInterface); #[pymethods] -impl PyPlanarInterface { +impl PyInterface { /// Initialize a planar interface with a hyperbolic tangent. /// /// Parameters @@ -40,35 +57,22 @@ impl PyPlanarInterface { /// critical_temperature: SINumber /// An estimate for the critical temperature of the system. /// Used to guess the width of the interface. - /// fix_equimolar_surface: bool, optional - /// If True use additional constraints to fix the - /// equimolar surface of the system. - /// Defaults to False. /// /// Returns /// ------- - /// PlanarInterface + /// Interface /// #[staticmethod] - #[pyo3( - text_signature = "(vle, n_grid, l_grid, critical_temperature, fix_equimolar_surface=None)" - )] - #[pyo3(signature = (vle, n_grid, l_grid, critical_temperature, fix_equimolar_surface=None))] - fn from_tanh( + #[pyo3(text_signature = "(vle, n_grid, l_grid, critical_temperature)")] + #[pyo3(signature = (vle, n_grid, l_grid, critical_temperature))] + fn planar_from_tanh( vle: &PyPhaseEquilibrium, n_grid: usize, l_grid: Length, critical_temperature: Temperature, - fix_equimolar_surface: Option, ) -> Self { - let profile = PlanarInterface::from_tanh( - &vle.0, - n_grid, - l_grid, - critical_temperature, - fix_equimolar_surface.unwrap_or(false), - ); - PyPlanarInterface(profile) + let profile = Interface::planar_from_tanh(&vle.0, n_grid, l_grid, critical_temperature); + PyInterface(profile) } /// Initialize a planar interface with a pDGT calculation. @@ -79,69 +83,83 @@ impl PyPlanarInterface { /// The bulk phase equilibrium. /// n_grid : int /// The number of grid points. - /// fix_equimolar_surface: bool, optional - /// If True use additional constraints to fix the - /// equimolar surface of the system. - /// Defaults to False. /// /// Returns /// ------- - /// PlanarInterface + /// Interface /// #[staticmethod] - #[pyo3(text_signature = "(vle, n_grid, fix_equimolar_surface=None)")] - #[pyo3(signature = (vle, n_grid, fix_equimolar_surface=None))] - fn from_pdgt( - vle: &PyPhaseEquilibrium, - n_grid: usize, - fix_equimolar_surface: Option, - ) -> PyResult { - let profile = - PlanarInterface::from_pdgt(&vle.0, n_grid, fix_equimolar_surface.unwrap_or(false)) - .map_err(PyFeosError::from)?; - Ok(PyPlanarInterface(profile)) + #[pyo3(text_signature = "(vle, n_grid)")] + #[pyo3(signature = (vle, n_grid))] + fn planar_from_pdgt(vle: &PyPhaseEquilibrium, n_grid: usize) -> PyResult { + let profile = Interface::planar_from_pdgt(&vle.0, n_grid).map_err(PyFeosError::from)?; + Ok(Self(profile)) } - /// Initialize a planar interface with a provided density profile. + /// Initialize a curved interface from a converged planar interface. + /// + /// The planar density profile is interpolated onto the new grid so that + /// its equimolar surface is located at the specified radius. /// /// Parameters /// ---------- - /// vle : PhaseEquilibrium - /// The bulk phase equilibrium. + /// planar_interface : Interface + /// A (converged) planar interface. /// n_grid : int /// The number of grid points. /// l_grid: SINumber /// The width of the calculation domain. - /// density_profile: SIArray2 - /// Initial condition for the density profile iterations + /// radius: SINumber + /// The initial radius of the equimolar surface. A positive radius + /// corresponds to a droplet, a negative radius to a bubble. + /// geometry: Geometry + /// The geometry of the interface (cylindrical or spherical). /// /// Returns /// ------- - /// PlanarInterface + /// Interface /// #[staticmethod] - fn from_density_profile( - vle: &PyPhaseEquilibrium, + #[pyo3(text_signature = "(planar_interface, n_grid, l_grid, radius, geometry)")] + fn curved( + planar_interface: &PyInterface, n_grid: usize, l_grid: Length, - density_profile: Density>, - ) -> Self { - let mut profile = PlanarInterface::new(&vle.0, n_grid, l_grid); - profile.profile.density = density_profile; - PyPlanarInterface(profile) + radius: Length, + geometry: PyGeometry, + ) -> PyResult { + let profile = + Interface::curved(&planar_interface.0, n_grid, l_grid, radius, geometry.into()) + .map_err(PyFeosError::from)?; + Ok(Self(profile)) + } + + #[new] + fn new( + grid: PyGrid, + vle: &PyPhaseEquilibrium, + density: Density>, + ) -> PyResult { + let (Grid::Cartesian1(axis) | Grid::Polar(axis) | Grid::Spherical(axis)) = grid.0 else { + return Err(PyFeosError::from(FeosError::Error( + "An interface requires a 1D cartesian, polar, or spherical grid.".into(), + )) + .into()); + }; + Ok(Self(Interface::new(&vle.0, axis, density))) } } #[pymethods] -impl PyPlanarInterface { +impl PyInterface { #[getter] fn get_surface_tension(&mut self) -> Option { self.0.surface_tension } #[getter] - fn get_equimolar_radius(&mut self) -> Option { - self.0.equimolar_radius + fn get_surface_of_tension(&mut self) -> Option { + self.0.surface_of_tension } #[getter] diff --git a/py-feos/src/dft/interface/surface_tension_diagram.rs b/py-feos/src/dft/interface/surface_tension_diagram.rs index e8a12461e..64c68c4ce 100644 --- a/py-feos/src/dft/interface/surface_tension_diagram.rs +++ b/py-feos/src/dft/interface/surface_tension_diagram.rs @@ -1,4 +1,4 @@ -use super::PyPlanarInterface; +use super::PyInterface; use crate::dft::PyDFTSolver; use crate::phase_equilibria::PyPhaseEquilibrium; use crate::state::PyStateVec; @@ -29,10 +29,6 @@ use std::sync::Arc; /// critical_temperature: SINumber, optional /// An estimate for the critical temperature, used to initialize /// density profile (default: 500 K) -/// fix_equimolar_surface: bool, optional -/// If True use additional constraints to fix the -/// equimolar surface of the system. -/// Defaults to False. /// solver: DFTSolver, optional /// Custom solver options /// @@ -48,16 +44,15 @@ pub struct PySurfaceTensionDiagram( impl PySurfaceTensionDiagram { #[new] #[pyo3( - text_signature = "(dia, init_densities=None, n_grid=None, l_grid=None, critical_temperature=None, fix_equimolar_surface=None, solver=None)" + text_signature = "(dia, init_densities=None, n_grid=None, l_grid=None, critical_temperature=None, solver=None)" )] - #[pyo3(signature = (dia, init_densities=None, n_grid=None, l_grid=None, critical_temperature=None, fix_equimolar_surface=None, solver=None))] + #[pyo3(signature = (dia, init_densities=None, n_grid=None, l_grid=None, critical_temperature=None, solver=None))] pub fn isotherm( dia: Vec, init_densities: Option, n_grid: Option, l_grid: Option, critical_temperature: Option, - fix_equimolar_surface: Option, solver: Option, ) -> PyResult { let x: Vec<_> = dia.into_iter().map(|vle| vle.0).collect(); @@ -67,17 +62,16 @@ impl PySurfaceTensionDiagram { n_grid, l_grid, critical_temperature, - fix_equimolar_surface, solver.map(|s| s.0).as_ref(), ))) } #[getter] - fn get_profiles(&self) -> Vec { + fn get_profiles(&self) -> Vec { self.0 .profiles .iter() - .map(|p| PyPlanarInterface(p.clone())) + .map(|p| PyInterface(p.clone())) .collect() } diff --git a/py-feos/src/dft/mod.rs b/py-feos/src/dft/mod.rs index 07d1a6839..138f26b17 100644 --- a/py-feos/src/dft/mod.rs +++ b/py-feos/src/dft/mod.rs @@ -18,7 +18,7 @@ mod solver; pub(crate) use adsorption::{ PyAdsorption, PyExternalPotential, PyGrid, PyPore1D, PyPoreProfile, PyPoreSpecification, }; -pub(crate) use interface::{PyPlanarInterface, PySurfaceTensionDiagram}; +pub(crate) use interface::{PyInterface, PySurfaceTensionDiagram}; pub(crate) use solvation::PyPairCorrelation; #[cfg(feature = "rayon")] pub(crate) use solvation::PySolvationProfile; diff --git a/py-feos/src/dft/profile.rs b/py-feos/src/dft/profile.rs index 5de657f24..252d80dc6 100644 --- a/py-feos/src/dft/profile.rs +++ b/py-feos/src/dft/profile.rs @@ -67,7 +67,7 @@ macro_rules! impl_profile { #[getter] fn get_temperature(&self) -> Temperature { - self.0.profile.bulk.temperature + self.0.profile.temperature } #[getter] @@ -90,11 +90,6 @@ macro_rules! impl_profile { self.0.profile.external_potential().clone().into_dyn() } - #[getter] - fn get_bulk(&self) -> PyState { - PyState(self.0.profile.bulk.clone()) - } - #[getter] fn get_solver_log(&self) -> PyDFTSolverLog { PyDFTSolverLog(self.0.profile.solver_log.clone()) @@ -189,24 +184,20 @@ macro_rules! impl_profile { Ok(self.0.profile.dn_dmu().map_err(PyFeosError::from)?) } - #[getter] - fn get_drho_dp(&mut self) -> PyResult<> as std::ops::Div>::Output> { - Ok(self.0.profile.drho_dp().map_err(PyFeosError::from)?.into_dyn()) + fn drho_dp(&mut self, bulk: &PyState) -> PyResult<> as std::ops::Div>::Output> { + Ok(self.0.profile.drho_dp(&bulk.0).map_err(PyFeosError::from)?.into_dyn()) } - #[getter] - fn get_dn_dp(&mut self) -> PyResult<> as std::ops::Div>::Output> { - Ok(self.0.profile.dn_dp().map_err(PyFeosError::from)?) + fn dn_dp(&mut self, bulk: &PyState) -> PyResult<> as std::ops::Div>::Output> { + Ok(self.0.profile.dn_dp(&bulk.0).map_err(PyFeosError::from)?) } - #[getter] - fn get_drho_dt(&mut self) -> PyResult<> as std::ops::Div>::Output> { - Ok(self.0.profile.drho_dt().map_err(PyFeosError::from)?.into_dyn()) + fn drho_dt(&mut self, bulk: &PyState) -> PyResult<> as std::ops::Div>::Output> { + Ok(self.0.profile.drho_dt(&bulk.0).map_err(PyFeosError::from)?.into_dyn()) } - #[getter] - fn get_dn_dt(&mut self) -> PyResult<> as std::ops::Div>::Output> { - Ok(self.0.profile.dn_dt().map_err(PyFeosError::from)?) + fn dn_dt(&mut self, bulk: &PyState) -> PyResult<> as std::ops::Div>::Output> { + Ok(self.0.profile.dn_dt(&bulk.0).map_err(PyFeosError::from)?) } } }; diff --git a/py-feos/src/dft/solvation.rs b/py-feos/src/dft/solvation.rs index 5da6ac9f8..0a5d3bb5c 100644 --- a/py-feos/src/dft/solvation.rs +++ b/py-feos/src/dft/solvation.rs @@ -79,6 +79,11 @@ impl PySolvationProfile { )) } + #[getter] + fn get_bulk(&self) -> PyState { + PyState(self.0.bulk.clone()) + } + #[getter] fn get_grand_potential(&self) -> Option { self.0.grand_potential @@ -121,6 +126,11 @@ impl PyPairCorrelation { Self(PairCorrelation::new(&bulk.0, test_particle, n_grid, width)) } + #[getter] + fn get_bulk(&self) -> PyState { + PyState(self.0.bulk.clone()) + } + #[getter] fn get_pair_correlation_function<'py>( &self, diff --git a/py-feos/src/lib.rs b/py-feos/src/lib.rs index ee84f1fbb..dd8047f83 100644 --- a/py-feos/src/lib.rs +++ b/py-feos/src/lib.rs @@ -210,7 +210,7 @@ fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { // Interface m.add_class::()?; - m.add_class::()?; + m.add_class::()?; m.add_class::()?; #[cfg(feature = "rayon")]