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 57fbc017a..90eb1c3bd 100644 --- a/feos-dft/src/adsorption/pore.rs +++ b/feos-dft/src/adsorption/pore.rs @@ -7,18 +7,25 @@ use crate::profile::{DFTProfile, MAX_POTENTIAL}; use crate::solver::DFTSolver; use crate::WeightFunctionInfo; use feos_core::{Components, Contributions, EosResult, ReferenceSystem, State, StateBuilder}; -use ndarray::prelude::*; -use ndarray::Axis as Axis_nd; -use ndarray::RemoveAxis; +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; 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 +151,43 @@ where * Dimensionless::new(&self.profile.bulk.molefracs)) .sum()) } + + fn _henry_coefficients + Copy + ScalarOperand + DctNum>( + &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; + 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 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) + } } impl PoreSpecification for Pore1D { 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 1dffe4b3c..583d062de 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_coefficients(&self) -> HenryCoefficient> { + self.0.henry_coefficients() + } + + #[getter] + 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 }, ) }