Ideal gas properties from DFT (#263)

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.

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
 ```
 

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<N: DualNum<f64> + Copy>(&self, temperature: N) -> UnGraph<(), N> {
                 match self {
                     #(#bond_lengths,)*
                     _ => Graph::with_capacity(0, 0),

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")]

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<T> = Quantity<T, _HenryCoefficient>;
+
 /// 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<N: DualNum<f64> + Copy + ScalarOperand + DctNum>(
+        &self,
+        temperature: N,
+    ) -> Array1<N> {
+        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<WeightFunctionInfo<N>> = 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<Array1<f64>> {
+        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<Array1<f64>> {
+        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<Ix1> for Pore1D {

feos-dft/src/functional.rs

@@ -32,7 +32,7 @@ impl<I: Components + Send + Sync, F: HelmholtzEnergyFunctional> HelmholtzEnergyF
         self.residual.compute_max_density(moles)
     }
 
-    fn bond_lengths(&self, temperature: f64) -> UnGraph<(), f64> {
+    fn bond_lengths<N: DualNum<f64> + 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>) -> f64;
 
     /// Overwrite this, if the functional consists of heterosegmented chains.
-    fn bond_lengths(&self, _temperature: f64) -> UnGraph<(), f64> {
+    fn bond_lengths<N: DualNum<f64> + 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<D>(
+    fn functional_derivative<D, N: DualNum<f64> + Copy + ScalarOperand>(
         &self,
-        temperature: f64,
-        density: &Array<f64, D::Larger>,
-        convolver: &Arc<dyn Convolver<f64, D>>,
-    ) -> EosResult<(Array<f64, D>, Array<f64, D::Larger>)>
+        temperature: N,
+        density: &Array<N, D::Larger>,
+        convolver: &Arc<dyn Convolver<N, D>>,
+    ) -> EosResult<(Array<N, D>, Array<N, D::Larger>)>
     where
         D: Dimension,
         D::Larger: Dimension<Smaller = D>,
@@ -226,50 +226,13 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync {
         ))
     }
 
-    #[expect(clippy::type_complexity)]
-    fn functional_derivative_dual<D>(
-        &self,
-        temperature: f64,
-        density: &Array<f64, D::Larger>,
-        convolver: &Arc<dyn Convolver<Dual64, D>>,
-    ) -> EosResult<(Array<Dual64, D>, Array<Dual64, D::Larger>)>
-    where
-        D: Dimension,
-        D::Larger: Dimension<Smaller = D>,
-    {
-        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 += &phi;
-        }
-        Ok((
-            helmholtz_energy_density,
-            convolver.functional_derivative(&partial_derivatives) * temperature_dual,
-        ))
-    }
-
     /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$
-    fn bond_integrals<D>(
+    fn bond_integrals<D, N: DualNum<f64> + Copy>(
         &self,
-        temperature: f64,
-        exponential: &Array<f64, D::Larger>,
-        convolver: &Arc<dyn Convolver<f64, D>>,
-    ) -> Array<f64, D::Larger>
+        temperature: N,
+        exponential: &Array<N, D::Larger>,
+        convolver: &Arc<dyn Convolver<N, D>>,
+    ) -> Array<N, D::Larger>
     where
         D: Dimension,
         D::Larger: Dimension<Smaller = D>,
@@ -290,7 +253,7 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync {
             }
         }
 
-        let mut i_graph: Graph<_, Option<Array<f64, D>>, Directed> =
+        let mut i_graph: Graph<_, Option<Array<N, D>>, 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<f64, D>, e| acc * e.weight().as_ref().unwrap(),
+                            |acc: Array<N, D>, e| acc * e.weight().as_ref().unwrap(),
                         );
                         i1 = Some(convolver.convolve(i0, &bond_weight_functions[edge.id()]));
                         break 'nodes;

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,49 +46,26 @@ pub trait FunctionalContribution: Display + Sync + Send {
             * state.volume
     }
 
-    fn first_partial_derivatives(
+    fn first_partial_derivatives<N: DualNum<f64> + Copy>(
         &self,
-        temperature: f64,
-        weighted_densities: Array2<f64>,
-        mut helmholtz_energy_density: ArrayViewMut1<f64>,
-        mut first_partial_derivative: ArrayViewMut2<f64>,
-    ) -> 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<Dual64>,
-        mut helmholtz_energy_density: ArrayViewMut1<Dual64>,
-        mut first_partial_derivative: ArrayViewMut2<Dual64>,
+        temperature: N,
+        weighted_densities: Array2<N>,
+        mut helmholtz_energy_density: ArrayViewMut1<N>,
+        mut first_partial_derivative: ArrayViewMut2<N>,
     ) -> EosResult<()> {
         let mut wd = weighted_densities.mapv(Dual::from_re);
         let t = Dual::from_re(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 = 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(())

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, D>) -> f64 {
+    fn integrate_reduced<N: DualNum<f64> + Copy>(&self, mut profile: Array<N, D>) -> 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<f64, D::Larger>) -> Array1<f64> {
+    fn integrate_reduced_comp<S: Data<Elem = N>, N: DualNum<f64> + Copy>(
+        &self,
+        profile: &ArrayBase<S, D::Larger>,
+    ) -> Array1<N> {
         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<S: Data<Elem = N>, N: DualNum<f64> + Copy>(
+        &self,
+        profile: &ArrayBase<S, D::Larger>,
+    ) -> Array1<N> {
+        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 Å.