Skip to content

Commit 2415513

Browse files
authored
Ideal gas properties from DFT (#263)
1 parent 283fde2 commit 2415513

File tree

11 files changed

+205
-170
lines changed

11 files changed

+205
-170
lines changed

docs/theory/dft/ideal_gas.md

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
# Ideal gas properties
2+
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):
3+
4+
$$\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)
5+
6+
$$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'}$$
7+
8+
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
9+
10+
$$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}$$
11+
12+
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
13+
14+
$$N_i^\mathrm{ig}=k_\mathrm{B}T\rho_i^\mathrm{b}H_i(T)$$
15+
16+
The expression can be used in the general equation for the **enthalpy of adsorption** (see [here](enthalpy_of_adsorption.md))
17+
18+
$$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}$$
19+
20+
to simplify to
21+
22+
$$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)$$
23+
24+
and finally
25+
26+
$$\Delta h_i^\mathrm{ads,ig}=-k_\mathrm{B}T^2\frac{H_i'(T)}{H_i(T)}$$
27+
28+
For a spherical molecule without bond integrals, the derivative can be evaluated straightforwardly to yield
29+
30+
$$\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}}$$
31+
32+
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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ This section explains the implementation of the core expressions from classical
1010
solver
1111
derivatives
1212
enthalpy_of_adsorption
13+
ideal_gas
1314
pdgt
1415
```
1516

feos-derive/src/dft.rs

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -141,7 +141,7 @@ fn impl_helmholtz_energy_functional(
141141
#(#contributions,)*
142142
}
143143
}
144-
fn bond_lengths(&self, temperature: f64) -> UnGraph<(), f64> {
144+
fn bond_lengths<N: DualNum<f64> + Copy>(&self, temperature: N) -> UnGraph<(), N> {
145145
match self {
146146
#(#bond_lengths,)*
147147
_ => Graph::with_capacity(0, 0),

feos-dft/src/adsorption/mod.rs

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ mod fea_potential;
1616
mod pore;
1717
mod pore2d;
1818
pub use external_potential::{ExternalPotential, FluidParameters};
19-
pub use pore::{Pore1D, PoreProfile, PoreProfile1D, PoreSpecification};
19+
pub use pore::{HenryCoefficient, Pore1D, PoreProfile, PoreProfile1D, PoreSpecification};
2020
pub use pore2d::{Pore2D, PoreProfile2D};
2121

2222
#[cfg(feature = "rayon")]

feos-dft/src/adsorption/pore.rs

Lines changed: 49 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -7,18 +7,25 @@ use crate::profile::{DFTProfile, MAX_POTENTIAL};
77
use crate::solver::DFTSolver;
88
use crate::WeightFunctionInfo;
99
use feos_core::{Components, Contributions, EosResult, ReferenceSystem, State, StateBuilder};
10-
use ndarray::prelude::*;
11-
use ndarray::Axis as Axis_nd;
12-
use ndarray::RemoveAxis;
10+
use ndarray::{prelude::*, ScalarOperand};
11+
use ndarray::{Axis as Axis_nd, RemoveAxis};
1312
use num_dual::linalg::LU;
14-
use num_dual::DualNum;
15-
use quantity::{Density, Dimensionless, Energy, Length, MolarEnergy, Temperature, Volume, KELVIN};
13+
use num_dual::{Dual64, DualNum};
14+
use quantity::{
15+
Density, Dimensionless, Energy, Length, MolarEnergy, Quantity, Temperature, Volume, _Moles,
16+
_Pressure, KELVIN, RGAS,
17+
};
18+
use rustdct::DctNum;
1619
use std::fmt::Display;
1720
use std::sync::Arc;
21+
use typenum::Diff;
1822

1923
const POTENTIAL_OFFSET: f64 = 2.0;
2024
const DEFAULT_GRID_POINTS: usize = 2048;
2125

26+
pub type _HenryCoefficient = Diff<_Moles, _Pressure>;
27+
pub type HenryCoefficient<T> = Quantity<T, _HenryCoefficient>;
28+
2229
/// Parameters required to specify a 1D pore.
2330
pub struct Pore1D {
2431
pub geometry: Geometry,
@@ -144,6 +151,43 @@ where
144151
* Dimensionless::new(&self.profile.bulk.molefracs))
145152
.sum())
146153
}
154+
155+
fn _henry_coefficients<N: DualNum<f64> + Copy + ScalarOperand + DctNum>(
156+
&self,
157+
temperature: N,
158+
) -> Array1<N> {
159+
if self.profile.dft.m().iter().any(|&m| m != 1.0) {
160+
panic!("Henry coefficients can only be calculated for spherical and heterosegmented molecules!")
161+
};
162+
let pot = self.profile.external_potential.mapv(N::from)
163+
* self.profile.temperature.to_reduced()
164+
/ temperature;
165+
let exp_pot = pot.mapv(|v| (-v).exp());
166+
let functional_contributions = self.profile.dft.contributions();
167+
let weight_functions: Vec<WeightFunctionInfo<N>> = functional_contributions
168+
.map(|c| c.weight_functions(temperature))
169+
.collect();
170+
let convolver = ConvolverFFT::<_, D>::plan(&self.profile.grid, &weight_functions, None);
171+
let bonds = self
172+
.profile
173+
.dft
174+
.bond_integrals(temperature, &exp_pot, &convolver);
175+
self.profile.integrate_reduced_segments(&(exp_pot * bonds))
176+
}
177+
178+
pub fn henry_coefficients(&self) -> HenryCoefficient<Array1<f64>> {
179+
let t = self.profile.temperature.to_reduced();
180+
Volume::from_reduced(self._henry_coefficients(t)) / (RGAS * self.profile.temperature)
181+
}
182+
183+
pub fn ideal_gas_enthalpy_of_adsorption(&self) -> MolarEnergy<Array1<f64>> {
184+
let t = Dual64::from(self.profile.temperature.to_reduced()).derivative();
185+
let h_dual = self._henry_coefficients(t);
186+
let h = h_dual.mapv(|h| h.re);
187+
let dh = h_dual.mapv(|h| h.eps);
188+
let t = self.profile.temperature.to_reduced();
189+
RGAS * self.profile.temperature * Dimensionless::from_reduced((&h - t * dh) / h)
190+
}
147191
}
148192

149193
impl PoreSpecification<Ix1> for Pore1D {

feos-dft/src/functional.rs

Lines changed: 15 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ impl<I: Components + Send + Sync, F: HelmholtzEnergyFunctional> HelmholtzEnergyF
3232
self.residual.compute_max_density(moles)
3333
}
3434

35-
fn bond_lengths(&self, temperature: f64) -> UnGraph<(), f64> {
35+
fn bond_lengths<N: DualNum<f64> + Copy>(&self, temperature: N) -> UnGraph<(), N> {
3636
self.residual.bond_lengths(temperature)
3737
}
3838
}
@@ -158,7 +158,7 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync {
158158
fn compute_max_density(&self, moles: &Array1<f64>) -> f64;
159159

160160
/// Overwrite this, if the functional consists of heterosegmented chains.
161-
fn bond_lengths(&self, _temperature: f64) -> UnGraph<(), f64> {
161+
fn bond_lengths<N: DualNum<f64> + Copy>(&self, _temperature: N) -> UnGraph<(), N> {
162162
Graph::with_capacity(0, 0)
163163
}
164164

@@ -190,14 +190,14 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync {
190190
IdealChainContribution::new(&self.component_index(), &self.m())
191191
}
192192

193-
/// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{F}}{\delta\rho_i(\mathbf{r})}$.
193+
/// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{\beta F}}{\delta\rho_i(\mathbf{r})}$.
194194
#[expect(clippy::type_complexity)]
195-
fn functional_derivative<D>(
195+
fn functional_derivative<D, N: DualNum<f64> + Copy + ScalarOperand>(
196196
&self,
197-
temperature: f64,
198-
density: &Array<f64, D::Larger>,
199-
convolver: &Arc<dyn Convolver<f64, D>>,
200-
) -> EosResult<(Array<f64, D>, Array<f64, D::Larger>)>
197+
temperature: N,
198+
density: &Array<N, D::Larger>,
199+
convolver: &Arc<dyn Convolver<N, D>>,
200+
) -> EosResult<(Array<N, D>, Array<N, D::Larger>)>
201201
where
202202
D: Dimension,
203203
D::Larger: Dimension<Smaller = D>,
@@ -226,50 +226,13 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync {
226226
))
227227
}
228228

229-
#[expect(clippy::type_complexity)]
230-
fn functional_derivative_dual<D>(
231-
&self,
232-
temperature: f64,
233-
density: &Array<f64, D::Larger>,
234-
convolver: &Arc<dyn Convolver<Dual64, D>>,
235-
) -> EosResult<(Array<Dual64, D>, Array<Dual64, D::Larger>)>
236-
where
237-
D: Dimension,
238-
D::Larger: Dimension<Smaller = D>,
239-
{
240-
let temperature_dual = Dual64::from(temperature).derivative();
241-
let density_dual = density.mapv(Dual64::from);
242-
let weighted_densities = convolver.weighted_densities(&density_dual);
243-
let contributions = self.contributions();
244-
let mut partial_derivatives = Vec::new();
245-
let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
246-
for (c, wd) in contributions.zip(weighted_densities) {
247-
let nwd = wd.shape()[0];
248-
let ngrid = wd.len() / nwd;
249-
let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
250-
let mut pd = Array::zeros(wd.raw_dim());
251-
c.first_partial_derivatives_dual(
252-
temperature_dual,
253-
wd.into_shape_with_order((nwd, ngrid)).unwrap(),
254-
phi.view_mut().into_shape_with_order(ngrid).unwrap(),
255-
pd.view_mut().into_shape_with_order((nwd, ngrid)).unwrap(),
256-
)?;
257-
partial_derivatives.push(pd);
258-
helmholtz_energy_density += &phi;
259-
}
260-
Ok((
261-
helmholtz_energy_density,
262-
convolver.functional_derivative(&partial_derivatives) * temperature_dual,
263-
))
264-
}
265-
266229
/// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$
267-
fn bond_integrals<D>(
230+
fn bond_integrals<D, N: DualNum<f64> + Copy>(
268231
&self,
269-
temperature: f64,
270-
exponential: &Array<f64, D::Larger>,
271-
convolver: &Arc<dyn Convolver<f64, D>>,
272-
) -> Array<f64, D::Larger>
232+
temperature: N,
233+
exponential: &Array<N, D::Larger>,
234+
convolver: &Arc<dyn Convolver<N, D>>,
235+
) -> Array<N, D::Larger>
273236
where
274237
D: Dimension,
275238
D::Larger: Dimension<Smaller = D>,
@@ -290,7 +253,7 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync {
290253
}
291254
}
292255

293-
let mut i_graph: Graph<_, Option<Array<f64, D>>, Directed> =
256+
let mut i_graph: Graph<_, Option<Array<N, D>>, Directed> =
294257
bond_weight_functions.map(|_, _| (), |_, _| None);
295258

296259
let bonds = i_graph.edge_count();
@@ -319,7 +282,7 @@ pub trait HelmholtzEnergyFunctional: Components + Sized + Send + Sync {
319282
exponential
320283
.index_axis(Axis(0), edge.target().index())
321284
.to_owned(),
322-
|acc: Array<f64, D>, e| acc * e.weight().as_ref().unwrap(),
285+
|acc: Array<N, D>, e| acc * e.weight().as_ref().unwrap(),
323286
);
324287
i1 = Some(convolver.convolve(i0, &bond_weight_functions[edge.id()]));
325288
break 'nodes;

feos-dft/src/functional_contribution.rs

Lines changed: 8 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ use feos_core::{EosResult, StateHD};
33
use ndarray::prelude::*;
44
use ndarray::{RemoveAxis, ScalarOperand};
55
use num_dual::*;
6-
use num_traits::{One, Zero};
6+
use num_traits::Zero;
77
use std::fmt::Display;
88

99
/// Individual functional contribution that can be evaluated using generalized (hyper) dual numbers.
@@ -46,49 +46,26 @@ pub trait FunctionalContribution: Display + Sync + Send {
4646
* state.volume
4747
}
4848

49-
fn first_partial_derivatives(
49+
fn first_partial_derivatives<N: DualNum<f64> + Copy>(
5050
&self,
51-
temperature: f64,
52-
weighted_densities: Array2<f64>,
53-
mut helmholtz_energy_density: ArrayViewMut1<f64>,
54-
mut first_partial_derivative: ArrayViewMut2<f64>,
55-
) -> EosResult<()> {
56-
let mut wd = weighted_densities.mapv(Dual64::from);
57-
let t = Dual64::from(temperature);
58-
let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0)));
59-
60-
for i in 0..wd.shape()[0] {
61-
wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps = 1.0);
62-
phi = self.helmholtz_energy_density(t, wd.view())?;
63-
first_partial_derivative
64-
.index_axis_mut(Axis(0), i)
65-
.assign(&phi.mapv(|p| p.eps));
66-
wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps = 0.0);
67-
}
68-
helmholtz_energy_density.assign(&phi.mapv(|p| p.re));
69-
Ok(())
70-
}
71-
72-
fn first_partial_derivatives_dual(
73-
&self,
74-
temperature: Dual64,
75-
weighted_densities: Array2<Dual64>,
76-
mut helmholtz_energy_density: ArrayViewMut1<Dual64>,
77-
mut first_partial_derivative: ArrayViewMut2<Dual64>,
51+
temperature: N,
52+
weighted_densities: Array2<N>,
53+
mut helmholtz_energy_density: ArrayViewMut1<N>,
54+
mut first_partial_derivative: ArrayViewMut2<N>,
7855
) -> EosResult<()> {
7956
let mut wd = weighted_densities.mapv(Dual::from_re);
8057
let t = Dual::from_re(temperature);
8158
let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0)));
8259

8360
for i in 0..wd.shape()[0] {
8461
wd.index_axis_mut(Axis(0), i)
85-
.map_inplace(|x| x.eps = Dual::one());
62+
.map_inplace(|x| x.eps = N::one());
8663
phi = self.helmholtz_energy_density(t, wd.view())?;
8764
first_partial_derivative
8865
.index_axis_mut(Axis(0), i)
8966
.assign(&phi.mapv(|p| p.eps));
9067
wd.index_axis_mut(Axis(0), i)
91-
.map_inplace(|x| x.eps = Dual::zero());
68+
.map_inplace(|x| x.eps = N::zero());
9269
}
9370
helmholtz_energy_density.assign(&phi.mapv(|p| p.re));
9471
Ok(())

feos-dft/src/profile/mod.rs

Lines changed: 19 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ use feos_core::{Components, EosError, EosResult, ReferenceSystem, State};
66
use ndarray::{
77
Array, Array1, Array2, Array3, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix1, Ix2, Ix3,
88
};
9+
use num_dual::DualNum;
910
use quantity::{Density, Length, Moles, Quantity, Temperature, Volume, _Volume, DEGREES};
1011
use std::ops::{Add, MulAssign};
1112
use std::sync::Arc;
@@ -247,23 +248,38 @@ where
247248
}
248249
}
249250

250-
fn integrate_reduced(&self, mut profile: Array<f64, D>) -> f64 {
251+
fn integrate_reduced<N: DualNum<f64> + Copy>(&self, mut profile: Array<N, D>) -> N {
251252
let (integration_weights, functional_determinant) = self.grid.integration_weights();
252253

253254
for (i, w) in integration_weights.into_iter().enumerate() {
254255
for mut l in profile.lanes_mut(Axis_nd(i)) {
255-
l.mul_assign(w);
256+
l.mul_assign(&w.mapv(N::from));
256257
}
257258
}
258259
profile.sum() * functional_determinant
259260
}
260261

261-
fn integrate_reduced_comp(&self, profile: &Array<f64, D::Larger>) -> Array1<f64> {
262+
fn integrate_reduced_comp<S: Data<Elem = N>, N: DualNum<f64> + Copy>(
263+
&self,
264+
profile: &ArrayBase<S, D::Larger>,
265+
) -> Array1<N> {
262266
Array1::from_shape_fn(profile.shape()[0], |i| {
263267
self.integrate_reduced(profile.index_axis(Axis_nd(0), i).to_owned())
264268
})
265269
}
266270

271+
pub(crate) fn integrate_reduced_segments<S: Data<Elem = N>, N: DualNum<f64> + Copy>(
272+
&self,
273+
profile: &ArrayBase<S, D::Larger>,
274+
) -> Array1<N> {
275+
let integral = self.integrate_reduced_comp(profile);
276+
let mut integral_comp = Array1::zeros(self.dft.components());
277+
for (i, &j) in self.dft.component_index().iter().enumerate() {
278+
integral_comp[j] = integral[i];
279+
}
280+
integral_comp
281+
}
282+
267283
/// Return the volume of the profile.
268284
///
269285
/// In periodic directions, the length is assumed to be 1 Å.

0 commit comments

Comments
 (0)