feos-org · prehner · Jul 8, 2022 · Jul 8, 2022 · Jul 8, 2022
diff --git 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/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)

diff --git a/feos-dft/src/adsorption/external_potential.rs b/feos-dft/src/adsorption/external_potential.rs
@@ -34,6 +34,13 @@ pub enum ExternalPotential<U> {
         rho_s: f64,
         xi: Option<f64>,
     },
+    /// 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<f64>,
+        epsilon_k_sf: Array1<f64>,
+        rho_s: f64,
+        xi: Option<f64>,
+    },
     /// 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<U: EosUnit> ExternalPotential<U> {
                                 / ((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<U: EosUnit> ExternalPotential<U> {
                                     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<U: EosUnit> ExternalPotential<U> {
                                             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
@@ -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<f64>,
+        epsilon_k_sf: &PyArray1<f64>,
+        rho_s: f64,
+        xi: Option<f64>,
+    ) -> 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
@@ -150,7 +150,7 @@ impl<T: DualNum<f64>> WeightFunction<T> {
 }
 
 /// Possible weight function shapes.
-#[derive(Clone, Copy, PartialEq, Debug)]
+#[derive(Clone, Copy, PartialEq, Eq, Debug)]
 pub enum WeightFunctionShape {
     /// Heaviside step function
     Theta,