From 06999662a45af4cfbe65111c9833c8dcbd0d2881 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 8 Jul 2022 09:53:17 +0200 Subject: [PATCH 1/2] Add Steele potential with custom combining rules --- feos-core/src/phase_equilibria/mod.rs | 4 +- feos-dft/src/adsorption/external_potential.rs | 64 +++++++++++++++++++ .../python/adsorption/external_potential.rs | 35 ++++++++++ feos-dft/src/weight_functions.rs | 2 +- 4 files changed, 102 insertions(+), 3 deletions(-) diff --git a/feos-core/src/phase_equilibria/mod.rs b/feos-core/src/phase_equilibria/mod.rs index a98efe6b5..5bc5a5467 100644 --- a/feos-core/src/phase_equilibria/mod.rs +++ b/feos-core/src/phase_equilibria/mod.rs @@ -17,7 +17,7 @@ pub use phase_diagram_binary::PhaseDiagramHetero; pub use phase_diagram_pure::PhaseDiagram; /// Level of detail in the iteration output. -#[derive(Copy, Clone, PartialOrd, PartialEq)] +#[derive(Copy, Clone, PartialOrd, PartialEq, Eq)] #[cfg_attr(feature = "python", pyo3::pyclass)] pub enum Verbosity { /// Do not print output. @@ -37,7 +37,7 @@ impl Default for Verbosity { /// Options for the various phase equilibria solvers. /// /// If the values are [None], solver specific default -/// values are used. +/// values are used. #[derive(Copy, Clone, Default)] pub struct SolverOptions { /// Maximum number of iterations. diff --git a/feos-dft/src/adsorption/external_potential.rs b/feos-dft/src/adsorption/external_potential.rs index 47a14d7ff..9caee2d09 100644 --- a/feos-dft/src/adsorption/external_potential.rs +++ b/feos-dft/src/adsorption/external_potential.rs @@ -34,6 +34,13 @@ pub enum ExternalPotential { rho_s: f64, xi: Option, }, + /// Steele potential with custom combining rules: $V_i^\mathrm{ext}(z)=2\pi m_i\xi\varepsilon_{si}\sigma_{si}^2\Delta\rho_s\left(0.4\left(\frac{\sigma_{si}}{z}\right)^{10}-\left(\frac{\sigma_{si}}{z}\right)^4-\frac{\sigma_{si}^4}{3\Delta\left(z+0.61\Delta\right)^3}\right),~~~~\Delta=3.35$ + CustomSteele { + sigma_sf: Array1, + epsilon_k_sf: Array1, + rho_s: f64, + xi: Option, + }, /// Double well potential: $V_i^\mathrm{ext}(z)=\mathrm{min}\left(\frac{2\pi}{45} m_i\varepsilon_{2si}\sigma_{si}^3\rho_s\left(2\left(\frac{2\sigma_{si}}{z}\right)^9-15\left(\frac{2\sigma_{si}}{z}\right)^3\right),0\right)+\frac{2\pi}{45} m_i\varepsilon_{1si}\sigma_{si}^3\rho_s\left(2\left(\frac{\sigma_{si}}{z}\right)^9-15\left(\frac{\sigma_{si}}{z}\right)^3\right),~~~~\varepsilon_{1si}=\sqrt{\varepsilon_{1ss}\varepsilon_{ii}},~~~~\varepsilon_{2si}=\sqrt{\varepsilon_{2ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$ DoubleWell { sigma_ss: f64, @@ -138,6 +145,20 @@ impl ExternalPotential { / ((3.0 * DELTA_STEELE) * (z_grid + 0.61 * DELTA_STEELE).mapv(|x| x.powi(3)))) } + Self::CustomSteele { + sigma_sf, + epsilon_k_sf, + rho_s, + xi, + } => { + (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i]) + * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s) + * (0.4 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(10)) + - (sigma_sf[i] / z_grid).mapv(|x| x.powi(4)) + - sigma_sf[i].powi(4) + / ((3.0 * DELTA_STEELE) + * (z_grid + 0.61 * DELTA_STEELE).mapv(|x| x.powi(3)))) + } Self::DoubleWell { sigma_ss, epsilon1_k_ss, @@ -277,6 +298,23 @@ impl ExternalPotential { sigma_sf[i] / (pore_size + DELTA_STEELE * 0.61), )) } + Self::CustomSteele { + sigma_sf, + epsilon_k_sf, + rho_s, + xi, + } => { + (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i]) + * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s) + * (psi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size) + - psi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size) + - sigma_sf[i] / DELTA_STEELE + * phi( + 3, + &(r_grid / (pore_size + DELTA_STEELE * 0.61)), + sigma_sf[i] / (pore_size + DELTA_STEELE * 0.61), + )) + } Self::DoubleWell { sigma_ss, epsilon1_k_ss, @@ -440,6 +478,32 @@ impl ExternalPotential { pore_size + 0.61 * DELTA_STEELE, ))) } + Self::CustomSteele { + sigma_sf, + epsilon_k_sf, + rho_s, + xi, + } => { + (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i]) + * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s) + * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size) + - sum_n(4, r_grid, sigma_sf[i], pore_size) + - sigma_sf[i] / (3.0 * DELTA_STEELE) + * (sigma_sf[i].powi(3) + / r_grid + .mapv(|r| (pore_size + 0.61 * DELTA_STEELE - r).powi(3)) + + sigma_sf[i].powi(3) + / r_grid.mapv(|r| { + (pore_size + 0.61 * DELTA_STEELE + r).powi(3) + }) + + 1.5 + * sum_n( + 3, + r_grid, + sigma_sf[i], + pore_size + 0.61 * DELTA_STEELE, + ))) + } Self::DoubleWell { sigma_ss, epsilon1_k_ss, diff --git a/feos-dft/src/python/adsorption/external_potential.rs b/feos-dft/src/python/adsorption/external_potential.rs index 9b23e91d3..4226e7589 100644 --- a/feos-dft/src/python/adsorption/external_potential.rs +++ b/feos-dft/src/python/adsorption/external_potential.rs @@ -137,6 +137,41 @@ impl PyExternalPotential { }) } + /// Steele potential with custom combining rules + /// + /// .. math:: V_i^\mathrm{ext}(z)=2\pi m_i\xi\varepsilon_{si}\sigma_{si}^2\Delta\rho_s\left(0.4\left(\frac{\sigma_{si}}{z}\right)^{10}-\left(\frac{\sigma_{si}}{z}\right)^4-\frac{\sigma_{si}^4}{3\Delta\left(z+0.61\Delta\right)^3}\right),~~~~\Delta=3.35 + /// + /// Parameters + /// ---------- + /// sigma_sf : numpy.ndarray[float] + /// Solid-fluid interaction diameters. + /// epsilon_k_sf : numpy.ndarray[float] + /// Solid-fluid interaction energies. + /// rho_s : float + /// Density of the solid. + /// xi : float, optional + /// Binary wall-fluid interaction parameter. + /// + /// Returns + /// ------- + /// ExternalPotential + /// + #[staticmethod] + #[pyo3(text_signature = "(sigma_sf, epsilon_k_sf, rho_s, xi=None)")] + pub fn CustomSteele( + sigma_sf: &PyArray1, + epsilon_k_sf: &PyArray1, + rho_s: f64, + xi: Option, + ) -> Self { + Self(ExternalPotential::CustomSteele { + sigma_sf: sigma_sf.to_owned_array(), + epsilon_k_sf: epsilon_k_sf.to_owned_array(), + rho_s, + xi, + }) + } + /// Double well potential /// /// .. math:: V_i^\mathrm{ext}(z)=\mathrm{min}\left(\frac{2\pi}{45} m_i\varepsilon_{2si}\sigma_{si}^3\rho_s\left(2\left(\frac{2\sigma_{si}}{z}\right)^9-15\left(\frac{2\sigma_{si}}{z}\right)^3\right),0\right)+\frac{2\pi}{45} m_i\varepsilon_{1si}\sigma_{si}^3\rho_s\left(2\left(\frac{\sigma_{si}}{z}\right)^9-15\left(\frac{\sigma_{si}}{z}\right)^3\right),~~~~\varepsilon_{1si}=\sqrt{\varepsilon_{1ss}\varepsilon_{ii}},~~~~\varepsilon_{2si}=\sqrt{\varepsilon_{2ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right) diff --git a/feos-dft/src/weight_functions.rs b/feos-dft/src/weight_functions.rs index cb9c870dd..9184c8e01 100644 --- a/feos-dft/src/weight_functions.rs +++ b/feos-dft/src/weight_functions.rs @@ -150,7 +150,7 @@ impl> WeightFunction { } /// Possible weight function shapes. -#[derive(Clone, Copy, PartialEq, Debug)] +#[derive(Clone, Copy, PartialEq, Eq, Debug)] pub enum WeightFunctionShape { /// Heaviside step function Theta, From 7ca62d68b2f5bf53c15f96af3d61779f54d0c5dc Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 8 Jul 2022 10:01:03 +0200 Subject: [PATCH 2/2] update changelog --- feos-dft/CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md index 8c69a65d2..c3537a1ab 100644 --- a/feos-dft/CHANGELOG.md +++ b/feos-dft/CHANGELOG.md @@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added - Added getters for the fields of `Pore1D` in Python. [#30](https://github.com/feos-org/feos-dft/pull/30) +- Added Steele potential with custom combining rules. [#29](https://github.com/feos-org/feos/pull/29) ### Changed - Made FMT functional more flexible w.r.t. the shape of the weight functions. [#31](https://github.com/feos-org/feos-dft/pull/31)