Add new 3d_dft feature to avoid libc dependency

feos-org · g-bauer · Sep 13, 2022 · Sep 8, 2022 · Sep 8, 2022 · Sep 8, 2022
commit b014764a8b7a4b3f2ea3fdf24611d628875b49a7
diff --git a/Cargo.toml b/Cargo.toml
@@ -49,11 +49,12 @@ approx = "0.4"
 [features]
 default = []
 dft = ["feos-dft", "petgraph"]
+3d_dft = ["dft", "feos-dft/3d_dft"]
 estimator = []
 association = []
 pcsaft = ["association"]
 gc_pcsaft = ["association"]
 uvtheory = ["lazy_static"]
 pets = []
 python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python"]
-all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets"]
+all_models = ["estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets"]
diff --git a/feos-dft/Cargo.toml b/feos-dft/Cargo.toml
@@ -19,18 +19,18 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ]
 quantity = { version = "0.5", features = ["linalg"] }
 feos-core = { version = "0.3", path = "../feos-core" }
 num-dual = "0.5"
-ndarray = { version = "0.15", features = ["serde", "rayon"] }
-ndarray-stats = "0.5"
+ndarray = "0.15"
 rustdct = "0.7"
 rustfft = "6.0"
 ang = "0.6"
 num-traits = "0.2"
-libc = "0.2"
-gauss-quad = "0.1"
+libm = "0.2"
+gauss-quad = { version = "0.1", optional = true }
 petgraph = "0.6"
 numpy = { version = "0.16", optional = true }
 pyo3 = { version = "0.16", optional = true }
 
 [features]
 default = []
+3d_dft = ["gauss-quad", "ndarray/rayon"]
 python = ["pyo3", "numpy", "feos-core/python"]
diff --git a/feos-dft/src/adsorption/external_potential.rs b/feos-dft/src/adsorption/external_potential.rs
@@ -1,11 +1,14 @@
+#[cfg(feature = "3d_dft")]
 use crate::adsorption::fea_potential::calculate_fea_potential;
 use crate::functional::HelmholtzEnergyFunctional;
+#[cfg(feature = "3d_dft")]
 use crate::geometry::Geometry;
 use feos_core::EosUnit;
-use libc::c_double;
+use libm::tgamma;
 use ndarray::{Array1, Array2, Axis as Axis_nd};
+#[cfg(feature = "3d_dft")]
 use quantity::{QuantityArray2, QuantityScalar};
-use std::f64::consts::PI;
+use std::{f64::consts::PI, marker::PhantomData};
 
 const DELTA_STEELE: f64 = 3.35;
 
@@ -49,6 +52,7 @@ pub enum ExternalPotential<U> {
         rho_s: f64,
     },
     /// Free-energy averaged potential:
+    #[cfg(feature = "3d_dft")]
     FreeEnergyAveraged {
         coordinates: QuantityArray2<U>,
         sigma_ss: Array1<f64>,
@@ -61,6 +65,10 @@ pub enum ExternalPotential<U> {
 
     /// Custom potential
     Custom(Array2<f64>),
+
+    /// Needed to keep `FreeEnergyAveraged` optional
+    #[doc(hidden)]
+    Phantom(PhantomData<U>),
 }
 
 /// Parameters of the fluid required to evaluate the external potential.
@@ -69,6 +77,7 @@ pub trait FluidParameters: HelmholtzEnergyFunctional {
     fn sigma_ff(&self) -> &Array1<f64>;
 }
 
+#[allow(unused_variables)]
 impl<U: EosUnit> ExternalPotential<U> {
     // Evaluate the external potential in cartesian coordinates for a given grid and fluid parameters.
     pub fn calculate_cartesian_potential<P: FluidParameters>(
@@ -183,6 +192,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                             * (2.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
                                 - 15.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(3))))
                 }
+                #[cfg(feature = "3d_dft")]
                 Self::FreeEnergyAveraged {
                     coordinates,
                     sigma_ss,
@@ -211,7 +221,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         *cutoff_radius,
                     )
                 }
-                Self::Custom(_) => unreachable!(),
+                _ => unreachable!(),
             });
         }
         ext_pot
@@ -348,6 +358,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                             * sigma_sf[i].powi(3)
                             * *rho_s)
                 }
+                #[cfg(feature = "3d_dft")]
                 Self::FreeEnergyAveraged {
                     coordinates,
                     sigma_ss,
@@ -376,7 +387,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         *cutoff_radius,
                     )
                 }
-                Self::Custom(_) => unreachable!(),
+                _ => unreachable!(),
             });
         }
         ext_pot
@@ -537,6 +548,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                             * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
                                 - sum_n(4, r_grid, sigma_sf[i], pore_size)))
                 }
+                #[cfg(feature = "3d_dft")]
                 Self::FreeEnergyAveraged {
                     coordinates,
                     sigma_ss,
@@ -565,7 +577,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                         *cutoff_radius,
                     )
                 }
-                Self::Custom(_) => unreachable!(),
+                _ => unreachable!(),
             });
         }
         ext_pot
@@ -578,17 +590,17 @@ fn phi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
 
     (1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(m3n2)) * 4.0 * PI.sqrt() / n2m3
         * sigma_r.powf(n2m3)
-        * gamma(n as f64 - 0.5)
-        / gamma(n as f64)
+        * tgamma(n as f64 - 0.5)
+        / tgamma(n as f64)
         * taylor_2f1_phi(r_r, n)
 }
 
 fn psi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
     (1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(2.0 - 2.0 * n as f64))
         * 4.0
         * PI.sqrt()
-        * gamma(n as f64 - 0.5)
-        / gamma(n as f64)
+        * tgamma(n as f64 - 0.5)
+        / tgamma(n as f64)
         * sigma_r.powf(2.0 * n as f64 - 2.0)
         * taylor_2f1_psi(r_r, n)
 }
@@ -637,11 +649,3 @@ fn taylor_2f1_psi(x: &Array1<f64>, n: i32) -> Array1<f64> {
         _ => unreachable!(),
     }
 }
-
-extern "C" {
-    fn tgamma(x: c_double) -> c_double;
-}
-
-fn gamma(x: f64) -> f64 {
-    unsafe { tgamma(x) }
-}
diff --git a/feos-dft/src/adsorption/fea_potential.rs b/feos-dft/src/adsorption/fea_potential.rs
@@ -1,9 +1,9 @@
+use super::pore3d::{calculate_distance2, evaluate_lj_potential};
 use crate::profile::{CUTOFF_RADIUS, MAX_POTENTIAL};
 use crate::Geometry;
 use feos_core::EosUnit;
 use gauss_quad::GaussLegendre;
 use ndarray::{Array1, Array2, Zip};
-use ndarray_stats::SummaryStatisticsExt;
 use quantity::{QuantityArray2, QuantityScalar};
 use std::f64::consts::PI;
 use std::usize;
@@ -114,7 +114,7 @@ pub fn calculate_fea_potential<U: EosUnit>(
                 let distance2 = calculate_distance2(point, &coordinates, system_size);
                 let potential_sum: f64 = (0..sigma_sf.len())
                     .map(|alpha| {
-                        mi * evaluate(
+                        mi * evaluate_lj_potential(
                             distance2[alpha],
                             sigma_sf[alpha],
                             epsilon_k_sf[alpha],
@@ -125,44 +125,8 @@ pub fn calculate_fea_potential<U: EosUnit>(
                 potential_2d[[i1, i2]] = (-potential_sum.min(MAX_POTENTIAL)).exp();
             }
         }
-        *f = potential_2d.weighted_sum(&weights).unwrap();
+        *f = (potential_2d * &weights).sum();
     });
 
     -temperature * potential.map(|p| (p / weights_sum).ln())
 }
-
-// -temperature * potential.map(|p| (p / weights_sum).ln())
-
-/// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment
-fn evaluate(distance2: f64, sigma: f64, epsilon: f64, cutoff_radius2: f64) -> f64 {
-    let sigma_r = sigma.powi(2) / distance2;
-
-    let potential: f64 = if distance2 > cutoff_radius2 {
-        0.0
-    } else if distance2 == 0.0 {
-        f64::INFINITY
-    } else {
-        4.0 * epsilon * (sigma_r.powi(6) - sigma_r.powi(3))
-    };
-
-    potential
-}
-
-/// Evaluate the squared euclidian distance between a point and the coordinates of all solid atoms.
-fn calculate_distance2(
-    point: [f64; 3],
-    coordinates: &Array2<f64>,
-    system_size: [f64; 3],
-) -> Array1<f64> {
-    Array1::from_shape_fn(coordinates.ncols(), |i| {
-        let mut rx = coordinates[[0, i]] - point[0];
-        let mut ry = coordinates[[1, i]] - point[1];
-        let mut rz = coordinates[[2, i]] - point[2];
-
-        rx = rx - system_size[0] * (rx / system_size[0]).round();
-        ry = ry - system_size[1] * (ry / system_size[1]).round();
-        rz = rz - system_size[2] * (rz / system_size[2]).round();
-
-        rx.powi(2) + ry.powi(2) + rz.powi(2)
-    })
-}
diff --git a/feos-dft/src/adsorption/mod.rs b/feos-dft/src/adsorption/mod.rs
@@ -11,10 +11,16 @@ use std::iter;
 use std::rc::Rc;
 
 mod external_potential;
+#[cfg(feature = "3d_dft")]
 mod fea_potential;
 mod pore;
 pub use external_potential::{ExternalPotential, FluidParameters};
-pub use pore::{Pore1D, Pore3D, PoreProfile, PoreProfile1D, PoreProfile3D, PoreSpecification};
+pub use pore::{Pore1D, PoreProfile, PoreProfile1D, PoreSpecification};
+
+#[cfg(feature = "3d_dft")]
+mod pore3d;
+#[cfg(feature = "3d_dft")]
+pub use pore3d::{Pore3D, PoreProfile3D};
 
 const MAX_ITER_ADSORPTION_EQUILIBRIUM: usize = 50;
 const TOL_ADSORPTION_EQUILIBRIUM: f64 = 1e-8;