feos-org · prehner · Jul 13, 2022 · Jul 3, 2022 · Jul 3, 2022 · Jul 3, 2022
diff --git a/.github/workflows/release_saft.yml b/.github/workflows/release_saft.yml
@@ -0,0 +1,24 @@
+name: Release feos-saft
+
+on:
+  push:
+    tags: ["feos-saft-v*"]
+
+jobs:
+  release-crates-io:
+    name: Release crates.io
+    runs-on: ubuntu-latest
+    defaults:
+      run:
+        working-directory: ./feos-saft
+    steps:
+      - uses: actions/checkout@v2
+      - uses: actions-rs/toolchain@v1
+        with:
+          profile: minimal
+          toolchain: stable
+          override: true
+      - uses: katyo/publish-crates@v1
+        with:
+          registry-token: ${{ secrets.CRATES_IO_TOKEN }}
+          path: './feos-saft'
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -15,7 +15,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        crate: [feos-core, feos-dft]
+        crate: [feos-core, feos-dft, feos-saft]
     steps:
       - uses: actions/checkout@v3
       - name: Build

diff --git a/Cargo.toml b/Cargo.toml
@@ -16,7 +16,7 @@ categories = ["science"]
 features = ["all_models"]
 
 [workspace]
-members = ["feos-core", "feos-dft"]
+members = ["feos-core", "feos-dft", "feos-saft"]
 
 [lib]
 crate-type = ["rlib", "cdylib"]
@@ -26,7 +26,7 @@ quantity = "0.5"
 num-dual = "0.5"
 feos-core = { version = "0.3", path = "feos-core" }
 feos-dft = { version = "0.2", path = "feos-dft", optional = true }
-#feos-pets = { version = "0.1", features = ["python"] }
+feos-saft = { version = "0.1", path = "feos-saft", optional = true }
 numpy = { version = "0.16", optional = true }
 ndarray = { version = "0.15", features = ["approx"] }
 petgraph = { version = "0.6", optional = true }
@@ -48,11 +48,11 @@ approx = "0.4"
 
 [features]
 default = []
-dft = ["feos-dft", "petgraph"]
+dft = ["feos-dft", "petgraph", "feos-saft?/dft"]
 estimator = []
-pcsaft = []
-gc_pcsaft = []
+pcsaft = ["feos-saft/association"]
+gc_pcsaft = ["feos-saft/association"]
 uvtheory = ["lazy_static"]
-pets = []
-python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python"]
+pets = ["feos-saft"]
+python = ["pyo3", "numpy", "feos-core/python", "feos-dft?/python", "feos-saft?/python"]
 all_models = ["dft", "estimator", "pcsaft", "gc_pcsaft", "uvtheory", "pets"]
diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md
@@ -13,6 +13,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Made FMT functional more flexible w.r.t. the shape of the weight functions. [#31](https://github.com/feos-org/feos-dft/pull/31)
 - Changed interface of `PairCorrelationFunction` to facilitate the calculation of pair correlation functions in mixtures. [#29](https://github.com/feos-org/feos-dft/pull/29)
 
+### Removed
+- Moved the implementation of fundamental measure theory to the new `feos-saft` crate. [#27](https://github.com/feos-org/feos/pull/27)
+
 ## [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)

diff --git a/feos-dft/src/lib.rs b/feos-dft/src/lib.rs
@@ -7,7 +7,6 @@ pub mod adsorption;
 mod convolver;
 mod functional;
 mod functional_contribution;
-pub mod fundamental_measure_theory;
 mod geometry;
 mod ideal_chain_contribution;
 pub mod interface;

diff --git a/feos-saft/Cargo.toml b/feos-saft/Cargo.toml
@@ -0,0 +1,32 @@
+[package]
+name = "feos-saft"
+version = "0.1.0"
+authors = ["Gernot Bauer <bauer@itt.uni-stuttgart.de>",
+           "Philipp Rehner <prehner@ethz.ch"]
+edition = "2021"
+license = "MIT OR Apache-2.0"
+description = "Core traits and functionalities for the `feos` project."
+homepage = "https://github.com/feos-org"
+readme = "README.md"
+repository = "https://github.com/feos-org/feos"
+keywords = ["physics", "thermodynamics", "equations_of_state", "phase_equilibria"]
+categories = ["science"]
+workspace = ".."
+
+[package.metadata.docs.rs]
+rustdoc-args = [ "--html-in-header", "./docs-header.html" ]
+
+[dependencies]
+num-dual = { version = "0.5", features = ["linalg"] }
+feos-core = { version = "0.3", path = "../feos-core" }
+feos-dft = { version = "0.2", path = "../feos-dft", optional = true }
+ndarray = { version = "0.15", features = ["serde"] }
+pyo3 = { version = "0.16", optional = true }
+serde = "1.0"
+serde_json = "1.0"
+
+[features]
+default = []
+association = []
+dft = ["feos-dft"]
+python = ["pyo3"]
diff --git a/feos-saft/docs-header.html b/feos-saft/docs-header.html
@@ -0,0 +1,15 @@
+<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/katex@0.10.0/dist/katex.min.css" integrity="sha384-9eLZqc9ds8eNjO3TmqPeYcDj8n+Qfa4nuSiGYa6DjLNcv9BtN69ZIulL9+8CqC9Y" crossorigin="anonymous">
+<script src="https://cdn.jsdelivr.net/npm/katex@0.10.0/dist/katex.min.js"                  integrity="sha384-K3vbOmF2BtaVai+Qk37uypf7VrgBubhQreNQe9aGsz9lB63dIFiQVlJbr92dw2Lx" crossorigin="anonymous"></script>
+<script src="https://cdn.jsdelivr.net/npm/katex@0.10.0/dist/contrib/auto-render.min.js"    integrity="sha384-kmZOZB5ObwgQnS/DuDg6TScgOiWWBiVt0plIRkZCmE6rDZGrEOQeHM5PcHi+nyqe" crossorigin="anonymous"></script>
+<script>
+    document.addEventListener("DOMContentLoaded", function() {
+        renderMathInElement(document.body, {
+            delimiters: [
+                {left: "$$", right: "$$", display: true},
+                {left: "\\(", right: "\\)", display: false},
+                {left: "$", right: "$", display: false},
+                {left: "\\[", right: "\\]", display: true}
+            ]
+        });
+    });
+</script>
diff --git a/feos-saft/license-apache b/feos-saft/license-apache
@@ -0,0 +1 @@
+../license-apache
diff --git a/feos-saft/license-mit b/feos-saft/license-mit
@@ -0,0 +1 @@
+../license-mit
diff --git a/src/gc_pcsaft/dft/association.rs → feos-saft/src/association/dft.rs b/src/gc_pcsaft/dft/association.rs → feos-saft/src/association/dft.rs
@@ -1,60 +1,40 @@
-use super::parameter::GcPcSaftFunctionalParameters;
-use crate::gc_pcsaft::eos::association::{
-    assoc_site_frac_a, assoc_site_frac_ab, helmholtz_energy_density_cross_association,
-};
-use feos_core::EosError;
+use super::*;
+use crate::HardSphereProperties;
+use feos_core::EosResult;
 use feos_dft::{
     FunctionalContributionDual, WeightFunction, WeightFunctionInfo, WeightFunctionShape,
 };
-use ndarray::*;
 use num_dual::DualNum;
 use std::f64::consts::PI;
-use std::fmt;
 use std::ops::MulAssign;
-use std::rc::Rc;
 
 pub const N0_CUTOFF: f64 = 1e-9;
 
-#[derive(Clone)]
-pub struct AssociationFunctional {
-    parameters: Rc<GcPcSaftFunctionalParameters>,
-    max_iter: usize,
-    tol: f64,
-}
-
-impl AssociationFunctional {
-    pub fn new(parameters: &Rc<GcPcSaftFunctionalParameters>, max_iter: usize, tol: f64) -> Self {
-        Self {
-            parameters: parameters.clone(),
-            max_iter,
-            tol,
-        }
-    }
-}
-
-impl<N> FunctionalContributionDual<N> for AssociationFunctional
+impl<N, P> FunctionalContributionDual<N> for Association<P>
 where
     N: DualNum<f64> + ScalarOperand,
+    P: HardSphereProperties,
 {
     fn weight_functions(&self, temperature: N) -> WeightFunctionInfo<N> {
         let p = &self.parameters;
         let r = p.hs_diameter(temperature) * 0.5;
-        WeightFunctionInfo::new(p.component_index.clone(), false)
+        let [_, _, _, c3] = p.geometry_coefficients(temperature);
+        WeightFunctionInfo::new(p.component_index().into_owned(), false)
             .add(
                 WeightFunction::new_scaled(r.clone(), WeightFunctionShape::Delta),
                 false,
             )
             .add(
                 WeightFunction {
-                    prefactor: p.m.mapv(N::from),
+                    prefactor: c3.clone(),
                     kernel_radius: r.clone(),
                     shape: WeightFunctionShape::DeltaVec,
                 },
                 false,
             )
             .add(
                 WeightFunction {
-                    prefactor: p.m.mapv(N::from),
+                    prefactor: c3,
                     kernel_radius: r,
                     shape: WeightFunctionShape::Theta,
                 },
@@ -66,43 +46,40 @@ where
         &self,
         temperature: N,
         weighted_densities: ArrayView2<N>,
-    ) -> Result<Array1<N>, EosError> {
+    ) -> EosResult<Array1<N>> {
         let p = &self.parameters;
 
         // number of segments
-        let segments = p.m.len();
+        let n = self.association_parameters.component_index.len();
 
         // number of associating segments
-        let nassoc = p.assoc_segment.len();
+        let nassoc = self.association_parameters.assoc_comp.len();
 
         // number of dimensions
-        let dim = (weighted_densities.shape()[0] - 1) / segments - 1;
+        let dim = (weighted_densities.shape()[0] - 1) / n - 1;
 
         // weighted densities
-        let n0i = weighted_densities.slice_axis(Axis(0), Slice::new(0, Some(segments as isize), 1));
+        let n0i = weighted_densities.slice_axis(Axis(0), Slice::new(0, Some(n as isize), 1));
         let n2vi: Vec<_> = (0..dim)
             .map(|i| {
                 weighted_densities.slice_axis(
                     Axis(0),
-                    Slice::new(
-                        (segments * (i + 1)) as isize,
-                        Some((segments * (i + 2)) as isize),
-                        1,
-                    ),
+                    Slice::new((n * (i + 1)) as isize, Some((n * (i + 2)) as isize), 1),
                 )
             })
             .collect();
-        let n3 = weighted_densities.index_axis(Axis(0), segments * (dim + 1));
+        let n3 = weighted_densities.index_axis(Axis(0), n * (dim + 1));
 
         // calculate rho0 (only associating segments)
+        let [_, _, c2, _] = p.geometry_coefficients(temperature);
         let diameter = p.hs_diameter(temperature);
         let mut n2i = n0i.to_owned();
         for (i, mut n2i) in n2i.outer_iter_mut().enumerate() {
-            n2i.mul_assign(diameter[i].powi(2) * p.m[i] * PI);
+            n2i.mul_assign(diameter[i].powi(2) * c2[i] * PI);
         }
         let mut rho0: Array2<N> = (n2vi
             .iter()
-            .fold(Array2::zeros(n2i.raw_dim()), |acc, n2vi| acc + n2vi * n2vi)
+            .fold(Array::zeros(n0i.raw_dim()), |acc, n2vi| acc + n2vi * n2vi)
             / -(&n2i * &n2i)
             + 1.0)
             * n0i;
@@ -111,8 +88,9 @@ where
                 *rho0 = n0i;
             }
         });
-        let rho0 =
-            Array2::from_shape_fn((nassoc, n3.len()), |(i, j)| rho0[(p.assoc_segment[i], j)]);
+        let rho0 = Array2::from_shape_fn((nassoc, n3.len()), |(i, j)| {
+            rho0[(self.association_parameters.assoc_comp[i], j)]
+        });
 
         // calculate xi
         let n2v: Vec<_> = n2vi.iter().map(|n2vi| n2vi.sum_axis(Axis(0))).collect();
@@ -136,23 +114,24 @@ where
         // only one associating component
         if nassoc == 1 {
             // association strength
-            let k = &n2 * &n3i * diameter[p.assoc_segment[0]] * 0.5;
+            let k = &n2 * &n3i * diameter[self.association_parameters.assoc_comp[0]] * 0.5;
             let deltarho = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i)
-                * ((temperature.recip() * p.epsilon_k_aibj[(0, 0)]).exp_m1()
-                    * p.sigma3_kappa_aibj[(0, 0)])
+                * ((temperature.recip() * self.association_parameters.epsilon_k_aibj[(0, 0)])
+                    .exp_m1()
+                    * self.association_parameters.sigma3_kappa_aibj[(0, 0)])
                 * rho0.index_axis(Axis(0), 0);
 
-            let na = p.na[0];
-            let nb = p.nb[0];
+            let na = self.association_parameters.na[0];
+            let nb = self.association_parameters.nb[0];
             let f = |x: N| x.ln() - x * 0.5 + 0.5;
             if nb > 0.0 {
                 // no cross association, two association sites
-                let xa = deltarho.mapv(|d| assoc_site_frac_ab(d, na, nb));
+                let xa = deltarho.mapv(|d| Self::assoc_site_frac_ab(d, na, nb));
                 let xb = (&xa - 1.0) * (na / nb) + 1.0;
                 Ok((xa.mapv(f) * na + xb.mapv(f) * nb) * rho0.index_axis(Axis(0), 0))
             } else {
                 // no cross association, one association site
-                let xa = deltarho.mapv(|d| assoc_site_frac_a(d, na));
+                let xa = deltarho.mapv(|d| Self::assoc_site_frac_a(d, na));
 
                 Ok(xa.mapv(f) * na * rho0.index_axis(Axis(0), 0))
             }
@@ -167,12 +146,7 @@ where
                 .zip(n3i.iter())
                 .zip(xi.iter())
                 .map(|(((rho0, &n2), &n3i), &xi)| {
-                    helmholtz_energy_density_cross_association(
-                        &p.assoc_segment,
-                        &p.sigma3_kappa_aibj,
-                        &p.epsilon_k_aibj,
-                        &p.na,
-                        &p.nb,
+                    self.helmholtz_energy_density_cross_association(
                         temperature,
                         &rho0,
                         &diameter,
@@ -190,9 +164,3 @@ where
         }
     }
 }
-
-impl fmt::Display for AssociationFunctional {
-    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        write!(f, "Association functional")
-    }
-}