feos-org · prehner · Apr 12, 2022 · Mar 9, 2022 · Mar 10, 2022 · Mar 10, 2022
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -6,6 +6,28 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.2.0] - 2022-04-12
+### Added
+- Added `grand_potential_density` getter for DFT profiles in Python. [#22](https://github.com/feos-org/feos-dft/pull/22)
+
+### Changed
+- Renamed `AxisGeometry` to `Geometry`. [#19](https://github.com/feos-org/feos-dft/pull/19)
+- Removed `PyGeometry` and `PyFMTVersion` in favor of a simpler implementation using `PyO3`'s new `#[pyclass]` for fieldless enums feature. [#19](https://github.com/feos-org/feos-dft/pull/19)
+- `DFTSolver` now uses `Verbosity` instead of a `bool` to control its output. [#19](https://github.com/feos-org/feos-dft/pull/19)
+- `SurfaceTensionDiagram` now uses the new `StateVec` struct to access properties of the bulk phases. [#19](https://github.com/feos-org/feos-dft/pull/19)
+- `Pore1D::initialize` and `Pore3D::initialize` now accept initial values for the density profiles as optional arguments. [#24](https://github.com/feos-org/feos-dft/pull/24)
+- Internally restructured the `DFT` structure to avoid redundant data. [#24](https://github.com/feos-org/feos-dft/pull/24)
+- Removed the `m` function in `FluidParameters`, it is instead inferred from `HelmholtzEnergyFunctional` which is now a supertrait of `FluidParameters`. [#24](https://github.com/feos-org/feos-dft/pull/24)
+- Added optional field `cutoff_radius` to `ExternalPotential::FreeEnergyAveraged`. [#25](https://github.com/feos-org/feos-dft/pull/25)
+
+### Packaging
+- Updated `pyo3` and `numpy` dependencies to 0.16.
+- Updated `quantity` dependency to 0.5.
+- Updated `num-dual` dependency to 0.5.
+- Updated `feos-core` dependency to 0.2.
+- Updated `ang` dependency to 0.6.
+- Removed `log` dependency.
+
 ## [0.1.3] - 2022-02-17
 ### Fixed
 - The pore volume for `Pore3D` is now also accesible from Python. [#16](https://github.com/feos-org/feos-dft/pull/16)

diff --git a/Cargo.toml b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "feos-dft"
-version = "0.1.3"
+version = "0.2.0"
 authors = ["Philipp Rehner <prehner@ethz.ch>"]
 edition = "2018"
 license = "MIT OR Apache-2.0"
@@ -16,22 +16,21 @@ exclude = ["/.github/*", "*.ipynb"]
 rustdoc-args = [ "--html-in-header", "./docs-header.html" ]
 
 [dependencies]
-quantity = { version = "0.4", features = ["linalg"] }
-feos-core = "0.1"
-num-dual = "0.4"
+quantity = { version = "0.5", features = ["linalg"] }
+feos-core = "0.2"
+num-dual = "0.5"
 ndarray = { version = "0.15", features = ["serde", "rayon"] }
 ndarray-stats = "0.5"
 rustdct = "0.7"
 rustfft = "6.0"
-log = "0.4"
-ang = "0.5"
+ang = "0.6"
 num-traits = "0.2"
 libc = "0.2"
 gauss-quad = "0.1"
 petgraph = "0.6"
-numpy = { version = "0.15", optional = true }
-pyo3 = { version = "0.15", optional = true }
+numpy = { version = "0.16", optional = true }
+pyo3 = { version = "0.16", optional = true }
 
 [features]
 default = []
-python = ["pyo3", "numpy", "feos-core/python"]
+python = ["pyo3", "numpy", "feos-core/python"]
diff --git a/README.md b/README.md
@@ -12,7 +12,7 @@ Add this to your `Cargo.toml`
 
 ```toml
 [dependencies]
-feos-dft = "0.1"
+feos-dft = "0.2"
 ```
 
 ## Test building python wheel

diff --git a/build_wheel/Cargo.toml b/build_wheel/Cargo.toml
@@ -1,13 +1,15 @@
 [package]
 name = "feos_dft"
-version = "0.1.3"
-authors = ["Philipp Rehner <prehner@ethz.ch"]
+version = "0.2.0"
+authors = ["Philipp Rehner <prehner@ethz.ch>"]
 edition = "2018"
 
 [lib]
 crate-type = ["cdylib"]
 
 [dependencies]
+quantity = "0.5"
+feos-core = "0.2"
 feos-dft = { path = "..", features = ["python"] }
-pyo3 = { version = "0.15", features = ["extension-module", "abi3", "abi3-py36"] }
-
+pyo3 = { version = "0.16", features = ["extension-module", "abi3", "abi3-py37"] }
+numpy = "0.16"
diff --git a/build_wheel/src/lib.rs b/build_wheel/src/lib.rs
@@ -1,7 +1,86 @@
-use feos_dft::python::feos_dft;
+use feos_core::*;
+use feos_dft::adsorption::*;
+use feos_dft::fundamental_measure_theory::{FMTFunctional, FMTVersion};
+use feos_dft::python::{PyDFTSolver, PyExternalPotential};
+use feos_dft::solvation::PairCorrelation;
+use feos_dft::*;
+use numpy::*;
+use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
+use pyo3::wrap_pymodule;
+use quantity::python::*;
+use quantity::si::SIUnit;
+use std::rc::Rc;
+
+/// Helmholtz energy functional for hard sphere systems.
+///
+/// Parameters
+/// ----------
+/// sigma : numpy.ndarray[float]
+///     The diameters of the hard spheres in Angstrom.
+/// version : FMTVersion
+///     The specific version of FMT to be used.
+///
+/// Returns
+/// -------
+/// FMTFunctional
+#[pyclass(name = "FMTFunctional", unsendable)]
+#[pyo3(text_signature = "(sigma, version)")]
+#[derive(Clone)]
+pub struct PyFMTFunctional(Rc<DFT<FMTFunctional>>);
+
+#[pymethods]
+impl PyFMTFunctional {
+    #[new]
+    fn new(sigma: &PyArray1<f64>, version: FMTVersion) -> Self {
+        Self(Rc::new(FMTFunctional::new(
+            &sigma.to_owned_array(),
+            version,
+        )))
+    }
+}
+
+impl_equation_of_state!(PyFMTFunctional);
+
+impl_state!(DFT<FMTFunctional>, PyFMTFunctional);
+
+impl_pore!(FMTFunctional, PyFMTFunctional);
+impl_adsorption!(FMTFunctional, PyFMTFunctional);
+
+impl_pair_correlation!(FMTFunctional);
 
 #[pymodule]
-pub fn build_wheel(py: Python<'_>, m: &PyModule) -> PyResult<()> {
-    feos_dft(py, m)
+pub fn feos_dft(py: Python<'_>, m: &PyModule) -> PyResult<()> {
+    m.add_class::<Verbosity>()?;
+
+    m.add_class::<PyExternalPotential>()?;
+    m.add_class::<Geometry>()?;
+    m.add_class::<PyDFTSolver>()?;
+
+    m.add_class::<FMTVersion>()?;
+    m.add_class::<PyFMTFunctional>()?;
+
+    m.add_class::<PyState>()?;
+    m.add_class::<PyPore1D>()?;
+    m.add_class::<PyPore3D>()?;
+    m.add_class::<PyPairCorrelation>()?;
+    m.add_class::<PyAdsorption1D>()?;
+    m.add_class::<PyAdsorption3D>()?;
+
+    m.add_wrapped(wrap_pymodule!(quantity))?;
+
+    py.run(
+        "\
+import sys
+quantity.SINumber.__module__ = 'feos_dft.si'
+quantity.SIArray1.__module__ = 'feos_dft.si'
+quantity.SIArray2.__module__ = 'feos_dft.si'
+quantity.SIArray3.__module__ = 'feos_dft.si'
+quantity.SIArray4.__module__ = 'feos_dft.si'
+sys.modules['feos_dft.si'] = quantity
+    ",
+        None,
+        Some(m.dict()),
+    )?;
+    Ok(())
 }
diff --git a/examples/FundamentalMeasureTheory.ipynb b/examples/FundamentalMeasureTheory.ipynb
@@ -39,14 +39,14 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "CPU times: user 5.33 s, sys: 7.47 ms, total: 5.34 s\n",
-      "Wall time: 5.33 s\n"
+      "CPU times: user 5.27 s, sys: 13.5 ms, total: 5.28 s\n",
+      "Wall time: 5.26 s\n"
      ]
     },
     {
      "data": {
       "text/plain": [
-       "<matplotlib.legend.Legend at 0x7fe9ce507400>"
+       "<matplotlib.legend.Legend at 0x7f1e31ad7490>"
       ]
      },
      "execution_count": 3,

diff --git a/src/adsorption/external_potential.rs b/src/adsorption/external_potential.rs
@@ -1,5 +1,6 @@
 use crate::adsorption::fea_potential::calculate_fea_potential;
-use crate::geometry::AxisGeometry;
+use crate::functional::HelmholtzEnergyFunctional;
+use crate::geometry::Geometry;
 use feos_core::EosUnit;
 use libc::c_double;
 use ndarray::{Array1, Array2, Axis as Axis_nd};
@@ -48,17 +49,17 @@ pub enum ExternalPotential<U> {
         pore_center: [f64; 3],
         system_size: [QuantityScalar<U>; 3],
         n_grid: [usize; 2],
+        cutoff_radius: Option<f64>,
     },
 
     /// Custom potential
     Custom(Array2<f64>),
 }
 
 /// Parameters of the fluid required to evaluate the external potential.
-pub trait FluidParameters {
+pub trait FluidParameters: HelmholtzEnergyFunctional {
     fn epsilon_k_ff(&self) -> Array1<f64>;
     fn sigma_ff(&self) -> &Array1<f64>;
-    fn m(&self) -> Array1<f64>;
 }
 
 impl<U: EosUnit> ExternalPotential<U> {
@@ -168,6 +169,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                     pore_center,
                     system_size,
                     n_grid,
+                    cutoff_radius,
                 } => {
                     // combining rules
                     let epsilon_k_sf =
@@ -184,7 +186,8 @@ impl<U: EosUnit> ExternalPotential<U> {
                         system_size,
                         n_grid,
                         temperature,
-                        AxisGeometry::Cartesian,
+                        Geometry::Cartesian,
+                        *cutoff_radius,
                     )
                 }
                 Self::Custom(_) => unreachable!(),
@@ -314,6 +317,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                     pore_center,
                     system_size,
                     n_grid,
+                    cutoff_radius,
                 } => {
                     // combining rules
                     let epsilon_k_sf =
@@ -330,7 +334,8 @@ impl<U: EosUnit> ExternalPotential<U> {
                         system_size,
                         n_grid,
                         temperature,
-                        AxisGeometry::Polar,
+                        Geometry::Cylindrical,
+                        *cutoff_radius,
                     )
                 }
                 Self::Custom(_) => unreachable!(),
@@ -475,6 +480,7 @@ impl<U: EosUnit> ExternalPotential<U> {
                     pore_center,
                     system_size,
                     n_grid,
+                    cutoff_radius,
                 } => {
                     // combining rules
                     let epsilon_k_sf =
@@ -491,7 +497,8 @@ impl<U: EosUnit> ExternalPotential<U> {
                         system_size,
                         n_grid,
                         temperature,
-                        AxisGeometry::Spherical,
+                        Geometry::Spherical,
+                        *cutoff_radius,
                     )
                 }
                 Self::Custom(_) => unreachable!(),

diff --git a/src/adsorption/fea_potential.rs b/src/adsorption/fea_potential.rs
@@ -1,5 +1,5 @@
 use crate::profile::{CUTOFF_RADIUS, MAX_POTENTIAL};
-use crate::AxisGeometry;
+use crate::Geometry;
 use feos_core::EosUnit;
 use gauss_quad::GaussLegendre;
 use ndarray::{Array1, Array2, Zip};
@@ -19,13 +19,14 @@ pub fn calculate_fea_potential<U: EosUnit>(
     system_size: &[QuantityScalar<U>; 3],
     n_grid: &[usize; 2],
     temperature: f64,
-    geometry: AxisGeometry,
+    geometry: Geometry,
+    cutoff_radius: Option<f64>,
 ) -> Array1<f64> {
     // allocate external potential
     let mut potential: Array1<f64> = Array1::zeros(grid.len());
 
     // calculate squared cutoff radius
-    let cutoff_radius2 = CUTOFF_RADIUS.powi(2);
+    let cutoff_radius2 = cutoff_radius.unwrap_or(CUTOFF_RADIUS).powi(2);
 
     // dimensionless solid coordinates
     let coordinates = Array2::from_shape_fn(coordinates.raw_dim(), |(i, j)| {
@@ -45,7 +46,7 @@ pub fn calculate_fea_potential<U: EosUnit>(
     // Cylindrical coordinates => phi
     // Spherical coordinates => phi
     let (nodes1, weights1) = match geometry {
-        AxisGeometry::Cartesian => {
+        Geometry::Cartesian => {
             let nodes = Array1::linspace(
                 0.5 * system_size[1] / n_grid[0] as f64,
                 system_size[1] - 0.5 * system_size[1] / n_grid[0] as f64,
@@ -54,7 +55,7 @@ pub fn calculate_fea_potential<U: EosUnit>(
             let weights = Array1::from_elem(n_grid[0], system_size[1] / n_grid[0] as f64);
             (nodes, weights)
         }
-        AxisGeometry::Spherical | AxisGeometry::Polar => {
+        Geometry::Spherical | Geometry::Cylindrical => {
             let nodes = PI + Array1::from_vec(GaussLegendre::nodes_and_weights(n_grid[0]).0) * PI;
             let weights = Array1::from_vec(GaussLegendre::nodes_and_weights(n_grid[0]).1) * PI;
             (nodes, weights)
@@ -66,7 +67,7 @@ pub fn calculate_fea_potential<U: EosUnit>(
     // Cylindrical coordinates => z
     // Spherical coordinates => theta
     let (nodes2, weights2) = match geometry {
-        AxisGeometry::Polar | AxisGeometry::Cartesian => {
+        Geometry::Cylindrical | Geometry::Cartesian => {
             let nodes = Array1::linspace(
                 0.5 * system_size[2] / n_grid[1] as f64,
                 system_size[2] - 0.5 * system_size[2] / n_grid[1] as f64,
@@ -75,7 +76,7 @@ pub fn calculate_fea_potential<U: EosUnit>(
             let weights = Array1::from_elem(n_grid[1], system_size[2] / n_grid[1] as f64);
             (nodes, weights)
         }
-        AxisGeometry::Spherical => {
+        Geometry::Spherical => {
             let nodes = PI / 2.0
                 + Array1::from_vec(GaussLegendre::nodes_and_weights(n_grid[1]).0) * PI / 2.0;
             let weights = Array1::from_vec(GaussLegendre::nodes_and_weights(n_grid[1]).1) * PI
@@ -97,13 +98,13 @@ pub fn calculate_fea_potential<U: EosUnit>(
         for (i1, &n1) in nodes1.iter().enumerate() {
             for (i2, &n2) in nodes2.iter().enumerate() {
                 let point = match geometry {
-                    AxisGeometry::Cartesian => [grid[i0], n1, n2],
-                    AxisGeometry::Polar => [
+                    Geometry::Cartesian => [grid[i0], n1, n2],
+                    Geometry::Cylindrical => [
                         pore_center[0] + grid[i0] * n1.cos(),
                         pore_center[1] + grid[i0] * n1.sin(),
                         n2,
                     ],
-                    AxisGeometry::Spherical => [
+                    Geometry::Spherical => [
                         pore_center[0] + grid[i0] * n2.sin() * n1.cos(),
                         pore_center[1] + grid[i0] * n2.sin() * n1.sin(),
                         pore_center[2] + grid[i0] * n2.cos(),