From 37610cf2b08943917121b2bdd2f5cf95fa963f0a Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 5 Dec 2024 17:52:05 +0100 Subject: [PATCH 1/5] Update to PyO3 v0.23 --- Cargo.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Cargo.toml b/Cargo.toml index afba6891d..081c4fc0b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -51,7 +51,11 @@ typenum = "1.16" [dependencies.pyo3] version = "0.23" +<<<<<<< HEAD features = ["extension-module", "abi3", "abi3-py39"] +======= +features = ["extension-module", "abi3", "abi3-py37"] +>>>>>>> 13057cfa (Update to PyO3 v0.23) optional = true [dev-dependencies] From c74f2d3daf64547cf4edcb51e14c1b468fc7d898 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 25 Jan 2024 14:39:00 +0100 Subject: [PATCH 2/5] Ideal gas properties from DFT --- Cargo.toml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 081c4fc0b..afba6891d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -51,11 +51,7 @@ typenum = "1.16" [dependencies.pyo3] version = "0.23" -<<<<<<< HEAD features = ["extension-module", "abi3", "abi3-py39"] -======= -features = ["extension-module", "abi3", "abi3-py37"] ->>>>>>> 13057cfa (Update to PyO3 v0.23) optional = true [dev-dependencies] From a6dfc8d7018e43a28c591e93970fcb18565ce33a Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 25 Jan 2024 14:39:00 +0100 Subject: [PATCH 3/5] Ideal gas properties from DFT --- Cargo.toml | 2 +- feos-dft/src/adsorption/pore.rs | 41 +++++++++++++- feos-dft/src/python/adsorption/pore.rs | 76 +++++++++----------------- 3 files changed, 65 insertions(+), 54 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index afba6891d..430da3561 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -70,7 +70,7 @@ codegen-units = 1 [features] -default = [] +default = ["python", "pcsaft", "dft"] dft = ["feos-dft", "petgraph"] estimator = [] association = [] diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs index 57fbc017a..1afbce042 100644 --- a/feos-dft/src/adsorption/pore.rs +++ b/feos-dft/src/adsorption/pore.rs @@ -6,19 +6,27 @@ use crate::geometry::{Axis, Geometry, Grid}; use crate::profile::{DFTProfile, MAX_POTENTIAL}; use crate::solver::DFTSolver; use crate::WeightFunctionInfo; +use feos_core::si::{ + Density, Dimensionless, Energy, Length, MolarEnergy, MolarWeight, Quantity, Temperature, + Volume, _Moles, _Pressure, KELVIN, RGAS, +}; use feos_core::{Components, Contributions, EosResult, ReferenceSystem, State, StateBuilder}; +use feos_core::{Components, Contributions, EosResult, State, StateBuilder}; use ndarray::prelude::*; -use ndarray::Axis as Axis_nd; -use ndarray::RemoveAxis; +use ndarray::{Axis as Axis_nd, RemoveAxis}; use num_dual::linalg::LU; use num_dual::DualNum; use quantity::{Density, Dimensionless, Energy, Length, MolarEnergy, Temperature, Volume, KELVIN}; use std::fmt::Display; use std::sync::Arc; +use typenum::Diff; const POTENTIAL_OFFSET: f64 = 2.0; const DEFAULT_GRID_POINTS: usize = 2048; +pub type _HenryCoefficient = Diff<_Moles, _Pressure>; +pub type HenryCoefficient = Quantity; + /// Parameters required to specify a 1D pore. pub struct Pore1D { pub geometry: Geometry, @@ -144,6 +152,35 @@ where * Dimensionless::new(&self.profile.bulk.molefracs)) .sum()) } + + pub fn henry_coefficient(&self) -> EosResult>> { + let mut pot = self.profile.external_potential.clone(); + // pot.outer_iter_mut() + // .zip(self.profile.dft.m().iter()) + // .for_each(|(mut v, &m)| v /= m); + let exp_pot = Dimensionless::from_reduced(pot.mapv(|v| (-v).exp())); + let m = self.profile.dft.m().into_owned(); + Ok(self.profile.integrate_comp(&exp_pot) + / (RGAS * self.profile.temperature * Dimensionless::from_reduced(m))) + } + + pub fn ideal_gas_enthalpy_of_adsorption(&self) -> EosResult>> + where + D::Larger: Dimension, + { + let mut pot = self.profile.external_potential.clone(); + pot.outer_iter_mut() + .zip(self.profile.dft.m().iter()) + .for_each(|(mut v, &m)| v /= m); + let potm1 = Dimensionless::from_reduced(self.profile.external_potential.mapv(|v| 1.0 - v)); + let exp_pot = Dimensionless::from_reduced(pot.mapv(|v| (-v).exp())); + let h = self + .profile + .integrate_comp(&(potm1 * &exp_pot)) + .convert_into(self.profile.integrate_comp(&exp_pot)); + let m = self.profile.dft.m().into_owned(); + Ok(h * RGAS * self.profile.temperature / Dimensionless::from_reduced(m)) + } } impl PoreSpecification for Pore1D { diff --git a/feos-dft/src/python/adsorption/pore.rs b/feos-dft/src/python/adsorption/pore.rs index 1dffe4b3c..976aea922 100644 --- a/feos-dft/src/python/adsorption/pore.rs +++ b/feos-dft/src/python/adsorption/pore.rs @@ -28,6 +28,7 @@ macro_rules! impl_pore { pub struct PyPoreProfile1D(PoreProfile1D<$func>); impl_1d_profile!(PyPoreProfile1D, [get_r, get_z]); + impl_pore_profile!(PyPoreProfile1D); #[pymethods] impl PyPore1D { @@ -113,29 +114,6 @@ macro_rules! impl_pore { } } - #[pymethods] - impl PyPoreProfile1D { - #[getter] - fn get_grand_potential(&self) -> Option { - self.0.grand_potential - } - - #[getter] - fn get_interfacial_tension(&self) -> Option { - self.0.interfacial_tension - } - - #[getter] - fn get_partial_molar_enthalpy_of_adsorption(&self) -> PyResult>> { - Ok(self.0.partial_molar_enthalpy_of_adsorption()?.into()) - } - - #[getter] - fn get_enthalpy_of_adsorption(&self) -> PyResult { - Ok(self.0.enthalpy_of_adsorption()?.into()) - } - } - #[pyclass(name = "Pore2D")] pub struct PyPore2D(Pore2D); @@ -143,6 +121,7 @@ macro_rules! impl_pore { pub struct PyPoreProfile2D(PoreProfile2D<$func>); impl_2d_profile!(PyPoreProfile2D, get_x, get_y); + impl_pore_profile!(PyPoreProfile2D); #[pymethods] impl PyPore2D { @@ -198,38 +177,12 @@ macro_rules! impl_pore { } } - #[pymethods] - impl PyPoreProfile2D { - #[getter] - fn get_grand_potential(&self) -> Option { - self.0.grand_potential - } - - #[getter] - fn get_interfacial_tension(&self) -> Option { - self.0.interfacial_tension - } - - #[getter] - fn get_partial_molar_enthalpy_of_adsorption(&self) -> PyResult>> { - Ok(self.0.partial_molar_enthalpy_of_adsorption()?) - } - - #[getter] - fn get_enthalpy_of_adsorption(&self) -> PyResult { - Ok(self.0.enthalpy_of_adsorption()?) - } - } - /// Parameters required to specify a 3D pore. /// /// Parameters /// ---------- /// system_size : [SINumber; 3] /// The size of the unit cell. - /// angles : [Angle; 3] - /// The angles of the unit cell or `None` if the unit cell - /// is orthorombic /// n_grid : [int; 3] /// The number of grid points in each direction. /// coordinates : numpy.ndarray[float] @@ -238,6 +191,9 @@ macro_rules! impl_pore { /// The size parameters of all interaction sites. /// epsilon_k_ss : numpy.ndarray[float] /// The energy parameter of all interaction sites. + /// angles : [Angle; 3], optional + /// The angles of the unit cell or `None` if the unit cell + /// is orthorombic /// potential_cutoff: float, optional /// Maximum value for the external potential. /// cutoff_radius: SINumber, optional @@ -254,6 +210,7 @@ macro_rules! impl_pore { pub struct PyPoreProfile3D(PoreProfile3D<$func>); impl_3d_profile!(PyPoreProfile3D, get_x, get_y, get_z); + impl_pore_profile!(PyPoreProfile3D); #[pymethods] impl PyPore3D { @@ -320,9 +277,14 @@ macro_rules! impl_pore { Ok(self.0.pore_volume()?.into()) } } + }; +} +#[macro_export] +macro_rules! impl_pore_profile { + ($py_profile:ty) => { #[pymethods] - impl PyPoreProfile3D { + impl $py_profile { #[getter] fn get_grand_potential(&self) -> Option { self.0.grand_potential @@ -334,7 +296,9 @@ macro_rules! impl_pore { } #[getter] - fn get_partial_molar_enthalpy_of_adsorption(&self) -> PyResult>> { + fn get_partial_molar_enthalpy_of_adsorption( + &self, + ) -> PyResult>> { Ok(self.0.partial_molar_enthalpy_of_adsorption()?) } @@ -342,6 +306,16 @@ macro_rules! impl_pore { fn get_enthalpy_of_adsorption(&self) -> PyResult { Ok(self.0.enthalpy_of_adsorption()?) } + + #[getter] + fn get_henry_coefficient(&self) -> PyResult { + Ok(self.0.henry_coefficient()?.into()) + } + + #[getter] + fn get_ideal_gas_enthalpy_of_adsorption(&self) -> PyResult { + Ok(self.0.ideal_gas_enthalpy_of_adsorption()?.into()) + } } }; } From 12bc6730ae3c3b475bb67cbe9ac4c71d1ff42da4 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sun, 8 Dec 2024 19:47:04 +0100 Subject: [PATCH 4/5] Implementation improvements --- Cargo.toml | 2 +- docs/theory/dft/ideal_gas.md | 32 +++++++++++ docs/theory/dft/index.md | 1 + feos-derive/src/dft.rs | 2 +- feos-dft/src/adsorption/mod.rs | 2 +- feos-dft/src/adsorption/pore.rs | 70 ++++++++++++----------- feos-dft/src/functional.rs | 67 +++++----------------- feos-dft/src/functional_contribution.rs | 39 +++---------- feos-dft/src/profile/mod.rs | 22 ++++++- feos-dft/src/profile/properties.rs | 76 +++++++++++++++++-------- feos-dft/src/python/adsorption/pore.rs | 8 +-- src/gc_pcsaft/dft/mod.rs | 4 +- 12 files changed, 173 insertions(+), 152 deletions(-) create mode 100644 docs/theory/dft/ideal_gas.md diff --git a/Cargo.toml b/Cargo.toml index 430da3561..afba6891d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -70,7 +70,7 @@ codegen-units = 1 [features] -default = ["python", "pcsaft", "dft"] +default = [] dft = ["feos-dft", "petgraph"] estimator = [] association = [] diff --git a/docs/theory/dft/ideal_gas.md b/docs/theory/dft/ideal_gas.md new file mode 100644 index 000000000..ebed3634d --- /dev/null +++ b/docs/theory/dft/ideal_gas.md @@ -0,0 +1,32 @@ +# Ideal gas properties +Classical DFT can be used to rapidly determine the ideal gas limit of fluids in porous media. In an ideal gas, there are no interactions between the fluid molecules and therefore the residual Helmholtz energy $F^\mathrm{res}$ and its derivatives vanish. Note that this is only the case for spherical or heterosegmented molecules ($m_\alpha=1$), as the chain contribution in the homosegmented model contains intramolecular interactions. The ideal gas density profile can then be obtained directly from the [Euler-Lagrange equation](euler_lagrange_equation.md): + +$$\rho_\alpha^\mathrm{ig}(\mathbf{r})=\rho_\alpha^\mathrm{b}e^{-\beta V_\alpha^\mathrm{ext}(\mathbf{r})}\prod_{\alpha'}I^\mathrm{ig}_{\alpha\alpha'}(\mathbf{r})$$ (eqn:rho_ideal_gas) + +$$I^\mathrm{ig}_{\alpha\alpha'}(\mathbf{r})=\int e^{-\beta V_{\alpha'}^\mathrm{ext}(\mathbf{r'})}\left(\prod_{\alpha''\neq\alpha}I^\mathrm{ig}_{\alpha'\alpha''}(\mathbf{r'})\right)\omega_\mathrm{chain}^{\alpha\alpha'}(\mathbf{r}-\mathbf{r'})\mathrm{d}\mathbf{r'}$$ + + Either from the derivatives presented [previously](derivatives.md), or from directly calculating derivatives of eq. {eq}`eqn:euler_lagrange_mu`, the **Henry coefficient** $H_\alpha$ can be calculated, as + + $$H_\alpha(T)=\left(\frac{\partial N_\alpha^\mathrm{ig}}{\partial p_\alpha}\right)_{T,x_k}=\int\left(\frac{\partial\rho_\alpha^\mathrm{ig}(\mathbf{r})}{\partial p_\alpha}\right)_{T,x_k}\mathrm{d}\mathbf{r}=\beta\int e^{-\beta V_\alpha^\mathrm{ext}(\mathbf{r})}\prod_{\alpha'}I^\mathrm{ig}_{\alpha\alpha'}(\mathbf{r})\mathrm{d}\mathbf{r}$$ + +By construction the Henry coefficients for all segments $\alpha$ of a molecule $i$ are identical. Therefore, the number of adsorbed molecules in an ideal gas state (the Henry regime) can be calculated from + +$$N_i^\mathrm{ig}=k_\mathrm{B}T\rho_i^\mathrm{b}H_i(T)$$ + +The expression can be used in the general equation for the **enthalpy of adsorption** (see [here](enthalpy_of_adsorption.md)) + +$$0=\sum_j\left(\frac{\partial N_i^\mathrm{ig}}{\partial\mu_j}\right)_T\Delta h_j^\mathrm{ads,ig}+T\left(\frac{\partial N_i^\mathrm{ig}}{\partial T}\right)_{p,x_k}$$ + +to simplify to + +$$0=\rho_i^\mathrm{b}H_i(T)\Delta h_i^\mathrm{ads,ig}+k_\mathrm{B}T^2\rho_i^\mathrm{b}H_i'(T)$$ + +and finally + +$$\Delta h_i^\mathrm{ads,ig}=-k_\mathrm{B}T^2\frac{H_i'(T)}{H_i(T)}$$ + +For a spherical molecule without bond integrals, the derivative can be evaluated straightforwardly to yield + +$$\Delta h_i^\mathrm{ads,ig}=\frac{\int\left(k_\mathrm{B}T-V_i^\mathrm{ext}(\mathbf{r})\right)e^{-\beta V_i^\mathrm{ext}(\mathbf{r})}\mathrm{d}\mathbf{r}}{\int e^{-\beta V_i^\mathrm{ext}(\mathbf{r})}\mathrm{d}\mathbf{r}}$$ + +Analytical derivatives of the bond integrals can be determined, however, in $\text{FeO}_\text{s}$ automatic differentiation with dual numbers is used to obtain correct derivatives with barely any implementation overhead. \ No newline at end of file diff --git a/docs/theory/dft/index.md b/docs/theory/dft/index.md index 3fb26446d..efd0331a0 100644 --- a/docs/theory/dft/index.md +++ b/docs/theory/dft/index.md @@ -10,6 +10,7 @@ This section explains the implementation of the core expressions from classical solver derivatives enthalpy_of_adsorption + ideal_gas pdgt ``` diff --git a/feos-derive/src/dft.rs b/feos-derive/src/dft.rs index 91225a8b9..d71ec58ba 100644 --- a/feos-derive/src/dft.rs +++ b/feos-derive/src/dft.rs @@ -141,7 +141,7 @@ fn impl_helmholtz_energy_functional( #(#contributions,)* } } - fn bond_lengths(&self, temperature: f64) -> UnGraph<(), f64> { + fn bond_lengths + Copy>(&self, temperature: N) -> UnGraph<(), N> { match self { #(#bond_lengths,)* _ => Graph::with_capacity(0, 0), diff --git a/feos-dft/src/adsorption/mod.rs b/feos-dft/src/adsorption/mod.rs index c96f3628b..447227fbc 100644 --- a/feos-dft/src/adsorption/mod.rs +++ b/feos-dft/src/adsorption/mod.rs @@ -16,7 +16,7 @@ mod fea_potential; mod pore; mod pore2d; pub use external_potential::{ExternalPotential, FluidParameters}; -pub use pore::{Pore1D, PoreProfile, PoreProfile1D, PoreSpecification}; +pub use pore::{HenryCoefficient, Pore1D, PoreProfile, PoreProfile1D, PoreSpecification}; pub use pore2d::{Pore2D, PoreProfile2D}; #[cfg(feature = "rayon")] diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs index 1afbce042..126949380 100644 --- a/feos-dft/src/adsorption/pore.rs +++ b/feos-dft/src/adsorption/pore.rs @@ -6,17 +6,16 @@ use crate::geometry::{Axis, Geometry, Grid}; use crate::profile::{DFTProfile, MAX_POTENTIAL}; use crate::solver::DFTSolver; use crate::WeightFunctionInfo; -use feos_core::si::{ - Density, Dimensionless, Energy, Length, MolarEnergy, MolarWeight, Quantity, Temperature, - Volume, _Moles, _Pressure, KELVIN, RGAS, -}; use feos_core::{Components, Contributions, EosResult, ReferenceSystem, State, StateBuilder}; -use feos_core::{Components, Contributions, EosResult, State, StateBuilder}; -use ndarray::prelude::*; +use ndarray::{prelude::*, ScalarOperand}; use ndarray::{Axis as Axis_nd, RemoveAxis}; use num_dual::linalg::LU; -use num_dual::DualNum; -use quantity::{Density, Dimensionless, Energy, Length, MolarEnergy, Temperature, Volume, KELVIN}; +use num_dual::{Dual64, DualNum}; +use quantity::{ + Density, Dimensionless, Energy, Length, MolarEnergy, Quantity, Temperature, Volume, _Moles, + _Pressure, KELVIN, RGAS, +}; +use rustdct::DctNum; use std::fmt::Display; use std::sync::Arc; use typenum::Diff; @@ -153,33 +152,38 @@ where .sum()) } - pub fn henry_coefficient(&self) -> EosResult>> { - let mut pot = self.profile.external_potential.clone(); - // pot.outer_iter_mut() - // .zip(self.profile.dft.m().iter()) - // .for_each(|(mut v, &m)| v /= m); - let exp_pot = Dimensionless::from_reduced(pot.mapv(|v| (-v).exp())); - let m = self.profile.dft.m().into_owned(); - Ok(self.profile.integrate_comp(&exp_pot) - / (RGAS * self.profile.temperature * Dimensionless::from_reduced(m))) + fn _henry_coefficients + Copy + ScalarOperand + DctNum>( + &self, + temperature: N, + ) -> Array1 { + let pot = self.profile.external_potential.mapv(N::from) + * self.profile.temperature.to_reduced() + / temperature; + let exp_pot = pot.mapv(|v| (-v).exp()); + let functional_contributions = self.profile.dft.contributions(); + let weight_functions: Vec> = functional_contributions + .map(|c| c.weight_functions(temperature)) + .collect(); + let convolver = ConvolverFFT::<_, D>::plan(&self.profile.grid, &weight_functions, None); + let bonds = self + .profile + .dft + .bond_integrals(temperature, &exp_pot, &convolver); + self.profile.integrate_reduced_segments(&(exp_pot * bonds)) } - pub fn ideal_gas_enthalpy_of_adsorption(&self) -> EosResult>> - where - D::Larger: Dimension, - { - let mut pot = self.profile.external_potential.clone(); - pot.outer_iter_mut() - .zip(self.profile.dft.m().iter()) - .for_each(|(mut v, &m)| v /= m); - let potm1 = Dimensionless::from_reduced(self.profile.external_potential.mapv(|v| 1.0 - v)); - let exp_pot = Dimensionless::from_reduced(pot.mapv(|v| (-v).exp())); - let h = self - .profile - .integrate_comp(&(potm1 * &exp_pot)) - .convert_into(self.profile.integrate_comp(&exp_pot)); - let m = self.profile.dft.m().into_owned(); - Ok(h * RGAS * self.profile.temperature / Dimensionless::from_reduced(m)) + pub fn henry_coefficients(&self) -> HenryCoefficient> { + 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.temperature.to_reduced()).derivative(); + let h_dual = self._henry_coefficients(t); + let h = h_dual.mapv(|h| h.re); + let dh = h_dual.mapv(|h| h.eps); + let t = self.profile.temperature.to_reduced(); + RGAS * self.profile.temperature * Dimensionless::from_reduced((&h - t * dh) / h) } } diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs index 8492daf7c..4f8460303 100644 --- a/feos-dft/src/functional.rs +++ b/feos-dft/src/functional.rs @@ -32,7 +32,7 @@ impl HelmholtzEnergyF self.residual.compute_max_density(moles) } - fn bond_lengths(&self, temperature: f64) -> UnGraph<(), f64> { + fn bond_lengths + Copy>(&self, temperature: N) -> UnGraph<(), N> { self.residual.bond_lengths(temperature) } } @@ -158,7 +158,7 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { fn compute_max_density(&self, moles: &Array1) -> f64; /// Overwrite this, if the functional consists of heterosegmented chains. - fn bond_lengths(&self, _temperature: f64) -> UnGraph<(), f64> { + fn bond_lengths + Copy>(&self, _temperature: N) -> UnGraph<(), N> { Graph::with_capacity(0, 0) } @@ -190,14 +190,14 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { IdealChainContribution::new(&self.component_index(), &self.m()) } - /// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{F}}{\delta\rho_i(\mathbf{r})}$. + /// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{\beta F}}{\delta\rho_i(\mathbf{r})}$. #[expect(clippy::type_complexity)] - fn functional_derivative( + fn functional_derivative + Copy + ScalarOperand>( &self, - temperature: f64, - density: &Array, - convolver: &Arc>, - ) -> EosResult<(Array, Array)> + temperature: N, + density: &Array, + convolver: &Arc>, + ) -> EosResult<(Array, Array)> where D: Dimension, D::Larger: Dimension, @@ -226,50 +226,13 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { )) } - #[expect(clippy::type_complexity)] - fn functional_derivative_dual( - &self, - temperature: f64, - density: &Array, - convolver: &Arc>, - ) -> EosResult<(Array, Array)> - where - D: Dimension, - D::Larger: Dimension, - { - let temperature_dual = Dual64::from(temperature).derivative(); - let density_dual = density.mapv(Dual64::from); - let weighted_densities = convolver.weighted_densities(&density_dual); - let contributions = self.contributions(); - let mut partial_derivatives = Vec::new(); - let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0))); - for (c, wd) in contributions.zip(weighted_densities) { - let nwd = wd.shape()[0]; - let ngrid = wd.len() / nwd; - let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); - let mut pd = Array::zeros(wd.raw_dim()); - c.first_partial_derivatives_dual( - temperature_dual, - wd.into_shape_with_order((nwd, ngrid)).unwrap(), - phi.view_mut().into_shape_with_order(ngrid).unwrap(), - pd.view_mut().into_shape_with_order((nwd, ngrid)).unwrap(), - )?; - partial_derivatives.push(pd); - helmholtz_energy_density += φ - } - Ok(( - helmholtz_energy_density, - convolver.functional_derivative(&partial_derivatives) * temperature_dual, - )) - } - /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$ - fn bond_integrals( + fn bond_integrals + Copy>( &self, - temperature: f64, - exponential: &Array, - convolver: &Arc>, - ) -> Array + temperature: N, + exponential: &Array, + convolver: &Arc>, + ) -> Array where D: Dimension, D::Larger: Dimension, @@ -290,7 +253,7 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { } } - let mut i_graph: Graph<_, Option>, Directed> = + let mut i_graph: Graph<_, Option>, Directed> = bond_weight_functions.map(|_, _| (), |_, _| None); let bonds = i_graph.edge_count(); @@ -319,7 +282,7 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync { exponential .index_axis(Axis(0), edge.target().index()) .to_owned(), - |acc: Array, e| acc * e.weight().as_ref().unwrap(), + |acc: Array, e| acc * e.weight().as_ref().unwrap(), ); i1 = Some(convolver.convolve(i0, &bond_weight_functions[edge.id()])); break 'nodes; diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs index b714fa48f..f88f89ffe 100644 --- a/feos-dft/src/functional_contribution.rs +++ b/feos-dft/src/functional_contribution.rs @@ -3,7 +3,7 @@ use feos_core::{EosResult, StateHD}; use ndarray::prelude::*; use ndarray::{RemoveAxis, ScalarOperand}; use num_dual::*; -use num_traits::{One, Zero}; +use num_traits::Zero; use std::fmt::Display; /// Individual functional contribution that can be evaluated using generalized (hyper) dual numbers. @@ -46,35 +46,12 @@ pub trait FunctionalContribution: Display + Sync + Send { * state.volume } - fn first_partial_derivatives( + fn first_partial_derivatives + Copy>( &self, - temperature: f64, - weighted_densities: Array2, - mut helmholtz_energy_density: ArrayViewMut1, - mut first_partial_derivative: ArrayViewMut2, - ) -> EosResult<()> { - let mut wd = weighted_densities.mapv(Dual64::from); - let t = Dual64::from(temperature); - let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0))); - - for i in 0..wd.shape()[0] { - wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps = 1.0); - phi = self.helmholtz_energy_density(t, wd.view())?; - first_partial_derivative - .index_axis_mut(Axis(0), i) - .assign(&phi.mapv(|p| p.eps)); - wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps = 0.0); - } - helmholtz_energy_density.assign(&phi.mapv(|p| p.re)); - Ok(()) - } - - fn first_partial_derivatives_dual( - &self, - temperature: Dual64, - weighted_densities: Array2, - mut helmholtz_energy_density: ArrayViewMut1, - mut first_partial_derivative: ArrayViewMut2, + temperature: N, + weighted_densities: Array2, + mut helmholtz_energy_density: ArrayViewMut1, + mut first_partial_derivative: ArrayViewMut2, ) -> EosResult<()> { let mut wd = weighted_densities.mapv(Dual::from_re); let t = Dual::from_re(temperature); @@ -82,13 +59,13 @@ pub trait FunctionalContribution: Display + Sync + Send { for i in 0..wd.shape()[0] { wd.index_axis_mut(Axis(0), i) - .map_inplace(|x| x.eps = Dual::one()); + .map_inplace(|x| x.eps = N::one()); phi = self.helmholtz_energy_density(t, wd.view())?; first_partial_derivative .index_axis_mut(Axis(0), i) .assign(&phi.mapv(|p| p.eps)); wd.index_axis_mut(Axis(0), i) - .map_inplace(|x| x.eps = Dual::zero()); + .map_inplace(|x| x.eps = N::zero()); } helmholtz_energy_density.assign(&phi.mapv(|p| p.re)); Ok(()) diff --git a/feos-dft/src/profile/mod.rs b/feos-dft/src/profile/mod.rs index fc6c675ee..9a33991f3 100644 --- a/feos-dft/src/profile/mod.rs +++ b/feos-dft/src/profile/mod.rs @@ -6,6 +6,7 @@ use feos_core::{Components, EosError, EosResult, ReferenceSystem, State}; use ndarray::{ Array, Array1, Array2, Array3, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix1, Ix2, Ix3, }; +use num_dual::DualNum; use quantity::{Density, Length, Moles, Quantity, Temperature, Volume, _Volume, DEGREES}; use std::ops::{Add, MulAssign}; use std::sync::Arc; @@ -247,23 +248,38 @@ where } } - fn integrate_reduced(&self, mut profile: Array) -> f64 { + fn integrate_reduced + Copy>(&self, mut profile: Array) -> N { let (integration_weights, functional_determinant) = self.grid.integration_weights(); for (i, w) in integration_weights.into_iter().enumerate() { for mut l in profile.lanes_mut(Axis_nd(i)) { - l.mul_assign(w); + l.mul_assign(&w.mapv(N::from)); } } profile.sum() * functional_determinant } - fn integrate_reduced_comp(&self, profile: &Array) -> Array1 { + fn integrate_reduced_comp, N: DualNum + Copy>( + &self, + profile: &ArrayBase, + ) -> Array1 { Array1::from_shape_fn(profile.shape()[0], |i| { self.integrate_reduced(profile.index_axis(Axis_nd(0), i).to_owned()) }) } + pub(crate) fn integrate_reduced_segments, N: DualNum + Copy>( + &self, + profile: &ArrayBase, + ) -> Array1 { + let integral = self.integrate_reduced_comp(profile); + let mut integral_comp = Array1::zeros(self.dft.components()); + for (i, &j) in self.dft.component_index().iter().enumerate() { + integral_comp[j] = integral[i]; + } + integral_comp + } + /// Return the volume of the profile. /// /// In periodic directions, the length is assumed to be 1 Å. diff --git a/feos-dft/src/profile/properties.rs b/feos-dft/src/profile/properties.rs index 2fd169024..d5d9cfc55 100644 --- a/feos-dft/src/profile/properties.rs +++ b/feos-dft/src/profile/properties.rs @@ -296,13 +296,16 @@ 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(&self) -> EosResult> { - let shape = self.density.shape(); - let shape: Vec<_> = std::iter::once(&shape[0]).chain(shape).copied().collect(); + let shape: Vec<_> = std::iter::once(&self.dft.components()) + .chain(self.density.shape()) + .copied() + .collect(); let mut drho_dmu = Array::zeros(shape).into_dimensionality().unwrap(); + let component_index = self.dft.component_index(); 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() { - if i != k { + if component_index[i] != k { l.fill(0.0); } } @@ -315,12 +318,14 @@ where /// Return the partial derivatives of the number of moles w.r.t. the chemical potentials $\left(\frac{\partial N_i}{\partial\mu_k}\right)_T$ pub fn dn_dmu(&self) -> EosResult { - let drho_dmu = self.drho_dmu()?; + let drho_dmu = self.drho_dmu()?.into_reduced(); let n = drho_dmu.shape()[0]; - let dn_dmu = DnDmu::from_shape_fn([n; 2], |(i, j)| { - self.integrate(&drho_dmu.index_axis(Axis(0), i).index_axis(Axis(0), j)) - }); - Ok(dn_dmu) + let mut dn_dmu = Array2::zeros([n; 2]); + dn_dmu + .outer_iter_mut() + .zip(drho_dmu.outer_iter()) + .for_each(|(mut dn, drho)| dn.assign(&self.integrate_reduced_segments(&drho))); + Ok(DnDmu::from_reduced(dn_dmu)) } /// 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}}$ @@ -341,50 +346,73 @@ where } /// 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}}$ - /// - /// Not compatible with heterosegmented DFT. pub fn drho_dt(&self) -> EosResult> { let rho = self.density.to_reduced(); let t = self.temperature.to_reduced(); + let rho_dual = rho.mapv(Dual64::from); + let t_dual = Dual64::from(t).derivative(); - // calculate temperature derivative of functional derivative + // calculate intrinsic functional derivative let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions - .map(|c| c.weight_functions(Dual64::from(t).derivative())) + .map(|c| c.weight_functions(t_dual)) .collect(); let convolver: Arc> = ConvolverFFT::plan(&self.grid, &weight_functions, None); - let (_, dfdrhodt) = self.dft.functional_derivative_dual(t, &rho, &convolver)?; + let (_, mut dfdrho) = self + .dft + .functional_derivative(t_dual, &rho_dual, &convolver)?; - // calculate temperature derivative of bulk functional derivative + // calculate total functional derivative + dfdrho += &((&self.external_potential * t).mapv(Dual64::from) / t_dual); + + // calculate bulk functional derivative let partial_density = self.bulk.partial_density.to_reduced(); let rho_bulk = self.dft.component_index().mapv(|i| partial_density[i]); + let rho_bulk_dual = rho_bulk.mapv(Dual64::from); let bulk_convolver = BulkConvolver::new(weight_functions); - let (_, dfdrhodt_bulk) = + let (_, dfdrho_bulk) = self.dft - .functional_derivative_dual(t, &rho_bulk, &bulk_convolver)?; + .functional_derivative(t_dual, &rho_bulk_dual, &bulk_convolver)?; + dfdrho + .outer_iter_mut() + .zip(dfdrho_bulk) + .zip(self.dft.m().iter()) + .for_each(|((mut df, df_b), &m)| { + df -= df_b; + df /= Dual64::from(m) + }); + + // calculate bond integrals + let exp_dfdrho = dfdrho.mapv(|x| (-x).exp()); + let bonds = self.dft.bond_integrals(t_dual, &exp_dfdrho, &convolver); // solve for drho_dt + let mut lhs = ((exp_dfdrho * bonds).mapv(|x| -x.ln()) * t_dual).mapv(|d| d.eps); let x = (self.bulk.partial_molar_volume() * self.bulk.dp_dt(Contributions::Total)).to_reduced(); - let mut lhs = dfdrhodt.mapv(|d| d.eps); - lhs.outer_iter_mut() - .zip(dfdrhodt_bulk) - .zip(x) - .for_each(|((mut lhs, d), x)| lhs -= d.eps - x); + let x = self.dft.component_index().mapv(|i| x[i]); lhs.outer_iter_mut() .zip(rho.outer_iter()) .zip(rho_bulk) .zip(self.dft.m().iter()) - .for_each(|(((mut lhs, rho), rho_b), &m)| lhs += &((&rho / rho_b).mapv(f64::ln) * m)); + .zip(x) + .for_each(|((((mut lhs, rho), rho_b), &m), x)| { + lhs += &(&rho / rho_b).mapv(f64::ln); + lhs *= m; + lhs += x; + }); lhs *= &(-&rho / t); + lhs.iter_mut().for_each(|l| { + if !l.is_finite() { + *l = 0.0 + } + }); Ok(Quantity::from_reduced(self.density_derivative(&lhs)?)) } /// 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}}$ - /// - /// Not compatible with heterosegmented DFT. pub fn dn_dt(&self) -> EosResult { Ok(self.integrate_segments(&self.drho_dt()?)) } diff --git a/feos-dft/src/python/adsorption/pore.rs b/feos-dft/src/python/adsorption/pore.rs index 976aea922..583d062de 100644 --- a/feos-dft/src/python/adsorption/pore.rs +++ b/feos-dft/src/python/adsorption/pore.rs @@ -308,13 +308,13 @@ macro_rules! impl_pore_profile { } #[getter] - fn get_henry_coefficient(&self) -> PyResult { - Ok(self.0.henry_coefficient()?.into()) + fn get_henry_coefficients(&self) -> HenryCoefficient> { + self.0.henry_coefficients() } #[getter] - fn get_ideal_gas_enthalpy_of_adsorption(&self) -> PyResult { - Ok(self.0.ideal_gas_enthalpy_of_adsorption()?.into()) + fn get_ideal_gas_enthalpy_of_adsorption(&self) -> MolarEnergy> { + self.0.ideal_gas_enthalpy_of_adsorption() } } }; diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 3f753dcc4..50191c9cd 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -110,7 +110,7 @@ impl HelmholtzEnergyFunctional for GcPcSaftFunctional { Box::new(contributions.into_iter()) } - fn bond_lengths(&self, temperature: f64) -> UnGraph<(), f64> { + fn bond_lengths + Copy>(&self, temperature: N) -> UnGraph<(), N> { // temperature dependent segment diameter let d = self.parameters.hs_diameter(temperature); @@ -120,7 +120,7 @@ impl HelmholtzEnergyFunctional for GcPcSaftFunctional { let (i, j) = self.parameters.bonds.edge_endpoints(e).unwrap(); let di = d[i.index()]; let dj = d[j.index()]; - 0.5 * (di + dj) + (di + dj) * 0.5 }, ) } From 9857203b8030fc237eec8d454eec33d882b261e7 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Tue, 10 Dec 2024 14:56:35 +0100 Subject: [PATCH 5/5] add error message, if ideal gas properties are calculated for homosegmented chains --- feos-dft/src/adsorption/pore.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs index 126949380..90eb1c3bd 100644 --- a/feos-dft/src/adsorption/pore.rs +++ b/feos-dft/src/adsorption/pore.rs @@ -156,6 +156,9 @@ where &self, temperature: N, ) -> Array1 { + if self.profile.dft.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.temperature.to_reduced() / temperature;