feos_dft/adsorption/
pore.rs

1use crate::WeightFunctionInfo;
2use crate::adsorption::{ExternalPotential, FluidParameters};
3use crate::convolver::ConvolverFFT;
4use crate::functional::{HelmholtzEnergyFunctional, HelmholtzEnergyFunctionalDyn, MoleculeShape};
5use crate::functional_contribution::FunctionalContribution;
6use crate::geometry::{Axis, Geometry, Grid};
7use crate::profile::{DFTProfile, MAX_POTENTIAL};
8use crate::solver::DFTSolver;
9use feos_core::{
10    Contributions, FeosResult, ReferenceSystem, ResidualDyn, State, StateBuilder, StateHD,
11};
12use nalgebra::{DVector, dvector};
13use ndarray::prelude::*;
14use ndarray::{Axis as Axis_nd, RemoveAxis};
15use num_dual::linalg::LU;
16use num_dual::{Dual64, DualNum};
17use quantity::{
18    _Moles, _Pressure, Density, Dimensionless, Energy, KELVIN, Length, MolarEnergy, Quantity, RGAS,
19    Temperature, Volume,
20};
21use rustdct::DctNum;
22use typenum::Diff;
23
24const POTENTIAL_OFFSET: f64 = 2.0;
25const DEFAULT_GRID_POINTS: usize = 2048;
26
27pub type _HenryCoefficient = Diff<_Moles, _Pressure>;
28pub type HenryCoefficient<T> = Quantity<T, _HenryCoefficient>;
29
30/// Parameters required to specify a 1D pore.
31pub struct Pore1D {
32    pub geometry: Geometry,
33    pub pore_size: Length,
34    pub potential: ExternalPotential,
35    pub n_grid: Option<usize>,
36    pub potential_cutoff: Option<f64>,
37}
38
39impl Pore1D {
40    pub fn new(
41        geometry: Geometry,
42        pore_size: Length,
43        potential: ExternalPotential,
44        n_grid: Option<usize>,
45        potential_cutoff: Option<f64>,
46    ) -> Self {
47        Self {
48            geometry,
49            pore_size,
50            potential,
51            n_grid,
52            potential_cutoff,
53        }
54    }
55}
56
57/// Trait for the generic implementation of adsorption applications.
58pub trait PoreSpecification<D: Dimension> {
59    /// Initialize a new single pore.
60    fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
61        &self,
62        bulk: &State<F>,
63        density: Option<&Density<Array<f64, D::Larger>>>,
64        external_potential: Option<&Array<f64, D::Larger>>,
65    ) -> FeosResult<PoreProfile<D, F>>;
66
67    /// Return the pore volume using Helium at 298 K as reference.
68    fn pore_volume(&self) -> FeosResult<Volume>
69    where
70        D::Larger: Dimension<Smaller = D>,
71    {
72        let bulk = StateBuilder::new(&&Helium)
73            .temperature(298.0 * KELVIN)
74            .density(Density::from_reduced(1.0))
75            .build()?;
76        let pore = self.initialize(&bulk, None, None)?;
77        let pot = Dimensionless::from_reduced(
78            pore.profile
79                .external_potential
80                .index_axis(Axis(0), 0)
81                .mapv(|v| (-v).exp()),
82        );
83        Ok(pore.profile.integrate(&pot))
84    }
85}
86
87/// Density profile and properties of a confined system in arbitrary dimensions.
88#[derive(Clone)]
89pub struct PoreProfile<D: Dimension, F> {
90    pub profile: DFTProfile<D, F>,
91    pub grand_potential: Option<Energy>,
92    pub interfacial_tension: Option<Energy>,
93}
94
95/// Density profile and properties of a 1D confined system.
96pub type PoreProfile1D<F> = PoreProfile<Ix1, F>;
97
98impl<D: Dimension + RemoveAxis + 'static, F: HelmholtzEnergyFunctional> PoreProfile<D, F>
99where
100    D::Larger: Dimension<Smaller = D>,
101    D::Smaller: Dimension<Larger = D>,
102    <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
103{
104    pub fn solve_inplace(&mut self, solver: Option<&DFTSolver>, debug: bool) -> FeosResult<()> {
105        // Solve the profile
106        self.profile.solve(solver, debug)?;
107
108        // calculate grand potential density
109        let omega = self.profile.grand_potential()?;
110        self.grand_potential = Some(omega);
111
112        // calculate interfacial tension
113        self.interfacial_tension =
114            Some(omega + self.profile.bulk.pressure(Contributions::Total) * self.profile.volume());
115
116        Ok(())
117    }
118
119    pub fn solve(mut self, solver: Option<&DFTSolver>) -> FeosResult<Self> {
120        self.solve_inplace(solver, false)?;
121        Ok(self)
122    }
123
124    pub fn update_bulk(mut self, bulk: &State<F>) -> Self {
125        self.profile.bulk = bulk.clone();
126        self.grand_potential = None;
127        self.interfacial_tension = None;
128        self
129    }
130
131    pub fn partial_molar_enthalpy_of_adsorption(&self) -> FeosResult<MolarEnergy<DVector<f64>>> {
132        let a = self.profile.dn_dmu()?;
133        let a_unit = a.get2(0, 0);
134        let b = -self.profile.temperature * self.profile.dn_dt()?;
135        let b_unit = b.get(0);
136
137        let h_ads = LU::new((a / a_unit).into_value())?.solve(&(b / b_unit).into_value());
138        Ok(&h_ads * b_unit / a_unit)
139    }
140
141    pub fn enthalpy_of_adsorption(&self) -> FeosResult<MolarEnergy> {
142        Ok((self.partial_molar_enthalpy_of_adsorption()?
143            * Dimensionless::new(&self.profile.bulk.molefracs))
144        .sum())
145    }
146
147    fn _henry_coefficients<N: DualNum<f64> + Copy + DctNum>(&self, temperature: N) -> DVector<N> {
148        if self.profile.bulk.eos.m().iter().any(|&m| m != 1.0) {
149            panic!(
150                "Henry coefficients can only be calculated for spherical and heterosegmented molecules!"
151            )
152        };
153        let pot = (self.profile.external_potential.mapv(N::from)
154            * self.profile.temperature.to_reduced())
155        .mapv(|v| v / temperature);
156        let exp_pot = pot.mapv(|v| (-v).exp());
157        let functional_contributions = self.profile.bulk.eos.contributions();
158        let weight_functions: Vec<WeightFunctionInfo<N>> = functional_contributions
159            .into_iter()
160            .map(|c| c.weight_functions(temperature))
161            .collect();
162        let convolver =
163            ConvolverFFT::<_, D>::plan(&self.profile.grid, &weight_functions, self.profile.lanczos);
164        let bonds = self
165            .profile
166            .bulk
167            .eos
168            .bond_integrals(temperature, &exp_pot, convolver.as_ref());
169        self.profile.integrate_reduced_segments(&(exp_pot * bonds))
170    }
171
172    pub fn henry_coefficients(&self) -> HenryCoefficient<DVector<f64>> {
173        let t = self.profile.temperature.to_reduced();
174        Volume::from_reduced(self._henry_coefficients(t)) / (RGAS * self.profile.temperature)
175    }
176
177    pub fn ideal_gas_enthalpy_of_adsorption(&self) -> MolarEnergy<DVector<f64>> {
178        let t = Dual64::from(self.profile.temperature.to_reduced()).derivative();
179        let h_dual = self._henry_coefficients(t);
180        let h = h_dual.map(|h| h.re);
181        let dh = h_dual.map(|h| h.eps);
182        let t = self.profile.temperature.to_reduced();
183        RGAS * self.profile.temperature
184            * Dimensionless::from_reduced((&h - t * dh).component_div(&h))
185    }
186}
187
188impl PoreSpecification<Ix1> for Pore1D {
189    fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
190        &self,
191        bulk: &State<F>,
192        density: Option<&Density<Array2<f64>>>,
193        external_potential: Option<&Array2<f64>>,
194    ) -> FeosResult<PoreProfile1D<F>> {
195        let dft: &F = &bulk.eos;
196        let n_grid = self.n_grid.unwrap_or(DEFAULT_GRID_POINTS);
197
198        let axis = match self.geometry {
199            Geometry::Cartesian => {
200                let potential_offset = POTENTIAL_OFFSET
201                    * bulk
202                        .eos
203                        .sigma_ff()
204                        .iter()
205                        .max_by(|a, b| a.total_cmp(b))
206                        .unwrap();
207                Axis::new_cartesian(n_grid, 0.5 * self.pore_size, Some(potential_offset))
208            }
209            Geometry::Cylindrical => Axis::new_polar(n_grid, self.pore_size),
210            Geometry::Spherical => Axis::new_spherical(n_grid, self.pore_size),
211        };
212
213        // calculate external potential
214        let external_potential = external_potential.map_or_else(
215            || {
216                external_potential_1d(
217                    self.pore_size,
218                    bulk.temperature,
219                    &self.potential,
220                    dft,
221                    &axis,
222                    self.potential_cutoff,
223                )
224            },
225            |e| e.clone(),
226        );
227
228        // initialize grid
229        let grid = Grid::new_1d(axis);
230
231        Ok(PoreProfile {
232            profile: DFTProfile::new(grid, bulk, Some(external_potential), density, Some(1)),
233            grand_potential: None,
234            interfacial_tension: None,
235        })
236    }
237}
238
239fn external_potential_1d<P: HelmholtzEnergyFunctional + FluidParameters>(
240    pore_width: Length,
241    temperature: Temperature,
242    potential: &ExternalPotential,
243    fluid_parameters: &P,
244    axis: &Axis,
245    potential_cutoff: Option<f64>,
246) -> Array2<f64> {
247    let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL);
248    let effective_pore_size = match axis.geometry {
249        Geometry::Spherical => pore_width.to_reduced(),
250        Geometry::Cylindrical => pore_width.to_reduced(),
251        Geometry::Cartesian => 0.5 * pore_width.to_reduced(),
252    };
253    let t = temperature.to_reduced();
254    let mut external_potential = match &axis.geometry {
255        Geometry::Cartesian => {
256            potential.calculate_cartesian_potential(
257                &(effective_pore_size + &axis.grid),
258                fluid_parameters,
259                t,
260            ) + &potential.calculate_cartesian_potential(
261                &(effective_pore_size - &axis.grid),
262                fluid_parameters,
263                t,
264            )
265        }
266        Geometry::Spherical => potential.calculate_spherical_potential(
267            &axis.grid,
268            effective_pore_size,
269            fluid_parameters,
270            t,
271        ),
272        Geometry::Cylindrical => potential.calculate_cylindrical_potential(
273            &axis.grid,
274            effective_pore_size,
275            fluid_parameters,
276            t,
277        ),
278    } / t;
279
280    for (i, &z) in axis.grid.iter().enumerate() {
281        if z > effective_pore_size {
282            external_potential
283                .index_axis_mut(Axis_nd(1), i)
284                .fill(potential_cutoff);
285        }
286    }
287    external_potential.map_inplace(|x| {
288        if *x > potential_cutoff {
289            *x = potential_cutoff
290        }
291    });
292    external_potential
293}
294
295const EPSILON_HE: f64 = 10.9;
296const SIGMA_HE: f64 = 2.64;
297
298#[derive(Clone, Copy)]
299struct Helium;
300
301impl ResidualDyn for Helium {
302    fn components(&self) -> usize {
303        1
304    }
305    fn compute_max_density<D: DualNum<f64> + Copy>(&self, _: &DVector<D>) -> D {
306        D::from(1.0)
307    }
308
309    fn reduced_helmholtz_energy_density_contributions<D: DualNum<f64> + Copy>(
310        &self,
311        state: &StateHD<D>,
312    ) -> Vec<(&'static str, D)> {
313        self.evaluate_bulk(state)
314    }
315}
316
317impl HelmholtzEnergyFunctionalDyn for Helium {
318    type Contribution<'a>
319        = HeliumContribution
320    where
321        Self: 'a;
322
323    fn contributions<'a>(&'a self) -> impl Iterator<Item = Self::Contribution<'a>> {
324        std::iter::empty()
325    }
326
327    fn molecule_shape(&self) -> MoleculeShape<'_> {
328        MoleculeShape::Spherical(1)
329    }
330}
331
332impl FluidParameters for &Helium {
333    fn epsilon_k_ff(&self) -> DVector<f64> {
334        dvector![EPSILON_HE]
335    }
336
337    fn sigma_ff(&self) -> DVector<f64> {
338        dvector![SIGMA_HE]
339    }
340}
341
342struct HeliumContribution;
343
344impl FunctionalContribution for HeliumContribution {
345    fn weight_functions<N: DualNum<f64> + Copy>(&self, _: N) -> WeightFunctionInfo<N> {
346        unreachable!()
347    }
348
349    fn helmholtz_energy_density<N: DualNum<f64> + Copy>(
350        &self,
351        _: N,
352        _: ArrayView2<N>,
353    ) -> FeosResult<Array1<N>> {
354        unreachable!()
355    }
356
357    fn name(&self) -> &'static str {
358        unreachable!()
359    }
360}
feos_dft/adsorption/pore.rs

feos_dft/adsorption/
pore.rs
