feos_dft/adsorption/
pore.rs

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

feos_dft/adsorption/
pore.rs
