From 409de331076584a62402344a54b2cd991ade0069 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Wed, 9 Mar 2022 16:06:59 +0100 Subject: [PATCH 01/14] Release v0.2.0 --- CHANGELOG.md | 8 ++++++++ Cargo.toml | 14 +++++++------- README.md | 2 +- build_wheel/Cargo.toml | 2 +- src/adsorption/mod.rs | 22 +++++++++++----------- src/functional.rs | 4 ++-- src/python/adsorption/mod.rs | 4 ++-- src/python/fundamental_measure_theory.rs | 1 - src/python/mod.rs | 2 +- src/python/profile.rs | 18 +++++++++--------- src/python/solver.rs | 3 --- 11 files changed, 42 insertions(+), 38 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f39271b..e2db505 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.2.0] - 2022-03-09 +### 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. + ## [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 index c2ac4dc..4d8ed43 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -16,22 +16,22 @@ 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 = { path = "../feos-core" } +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"] \ No newline at end of file +python = ["pyo3", "numpy", "feos-core/python"] diff --git a/README.md b/README.md index 1522431..3d581a1 100644 --- 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 index c012201..e05a1ba 100644 --- a/build_wheel/Cargo.toml +++ b/build_wheel/Cargo.toml @@ -9,5 +9,5 @@ crate-type = ["cdylib"] [dependencies] 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"] } diff --git a/src/adsorption/mod.rs b/src/adsorption/mod.rs index a629889..c5027fd 100644 --- a/src/adsorption/mod.rs +++ b/src/adsorption/mod.rs @@ -2,7 +2,7 @@ use super::functional::{HelmholtzEnergyFunctional, DFT}; use super::solver::DFTSolver; use feos_core::{ - Contributions, EosError, EosResult, EosUnit, EquationOfState, StateBuilder, VLEOptions, + Contributions, EosError, EosResult, EosUnit, EquationOfState, SolverOptions, StateBuilder, }; use ndarray::{arr1, Array1, Dimension, Ix1, Ix3}; use quantity::{QuantityArray1, QuantityArray2, QuantityScalar}; @@ -179,7 +179,7 @@ where pore, molefracs, solver, - VLEOptions::default(), + SolverOptions::default(), ); if let Ok(equilibrium) = equilibrium { let pressure = pressure.equilibrium(&equilibrium)?; @@ -259,9 +259,9 @@ where .pressure(pressure.get(0)) .moles(&moles) .build()?; - if functional.components() > 1 && !bulk.is_stable(VLEOptions::default())? { + if functional.components() > 1 && !bulk.is_stable(SolverOptions::default())? { bulk = bulk - .tp_flash(None, VLEOptions::default(), None)? + .tp_flash(None, SolverOptions::default(), None)? .vapor() .clone(); } @@ -273,9 +273,9 @@ where .pressure(pressure.get(i)) .moles(&moles) .build()?; - if functional.components() > 1 && !bulk.is_stable(VLEOptions::default())? { + if functional.components() > 1 && !bulk.is_stable(SolverOptions::default())? { bulk = bulk - .tp_flash(None, VLEOptions::default(), None)? + .tp_flash(None, SolverOptions::default(), None)? .vapor() .clone(); } @@ -299,7 +299,7 @@ where pore: &S, molefracs: Option<&Array1>, solver: Option<&DFTSolver>, - options: VLEOptions, + options: SolverOptions, ) -> EosResult> { let moles = functional.validate_moles(molefracs.map(|x| x * U::reference_moles()).as_ref())?; @@ -395,11 +395,11 @@ where QuantityArray1::from_shape_fn(self.profiles.len(), |i| match &self.profiles[i] { Ok(p) => { if p.profile.bulk.eos.components() > 1 - && !p.profile.bulk.is_stable(VLEOptions::default()).unwrap() + && !p.profile.bulk.is_stable(SolverOptions::default()).unwrap() { p.profile .bulk - .tp_flash(None, VLEOptions::default(), None) + .tp_flash(None, SolverOptions::default(), None) .unwrap() .vapor() .pressure(Contributions::Total) @@ -415,11 +415,11 @@ where QuantityArray1::from_shape_fn(self.profiles.len(), |i| match &self.profiles[i] { Ok(p) => { if p.profile.bulk.eos.components() > 1 - && !p.profile.bulk.is_stable(VLEOptions::default()).unwrap() + && !p.profile.bulk.is_stable(SolverOptions::default()).unwrap() { p.profile .bulk - .tp_flash(None, VLEOptions::default(), None) + .tp_flash(None, SolverOptions::default(), None) .unwrap() .vapor() .molar_gibbs_energy(Contributions::Total) diff --git a/src/functional.rs b/src/functional.rs index 4578f96..15a882c 100644 --- a/src/functional.rs +++ b/src/functional.rs @@ -186,8 +186,8 @@ impl DFT { D::Larger: Dimension, { // Calculate residual Helmholtz energy density and functional derivative - let t = temperature.to_reduced(U::reference_temperature()).unwrap(); - let rho = density.to_reduced(U::reference_density()).unwrap(); + let t = temperature.to_reduced(U::reference_temperature())?; + let rho = density.to_reduced(U::reference_density())?; let (mut f, dfdrho) = self.functional_derivative(t, &rho, convolver)?; // calculate the grand potential density diff --git a/src/python/adsorption/mod.rs b/src/python/adsorption/mod.rs index a09a8d2..e89e505 100644 --- a/src/python/adsorption/mod.rs +++ b/src/python/adsorption/mod.rs @@ -219,7 +219,7 @@ macro_rules! impl_adsorption_isotherm { solver: Option, max_iter: Option, tol: Option, - verbosity: Option, + verbosity: Option, ) -> PyResult { Ok(Self(Adsorption::phase_equilibrium( &functional.0, @@ -229,7 +229,7 @@ macro_rules! impl_adsorption_isotherm { &pore.0, molefracs.map(|x| x.to_owned_array()).as_ref(), solver.map(|s| s.0).as_ref(), - (max_iter, tol, verbosity.map(|v| v.0)).into(), + (max_iter, tol, verbosity).into(), )?)) } diff --git a/src/python/fundamental_measure_theory.rs b/src/python/fundamental_measure_theory.rs index 34625d9..b111a2a 100644 --- a/src/python/fundamental_measure_theory.rs +++ b/src/python/fundamental_measure_theory.rs @@ -4,7 +4,6 @@ use crate::functional::DFT; use crate::fundamental_measure_theory::{FMTFunctional, FMTVersion}; use crate::solvation::*; use crate::*; -use feos_core::python::{PyContributions, PyVerbosity}; use feos_core::*; use numpy::*; use pyo3::exceptions::PyValueError; diff --git a/src/python/mod.rs b/src/python/mod.rs index be762a5..2db7297 100644 --- a/src/python/mod.rs +++ b/src/python/mod.rs @@ -1,6 +1,6 @@ use pyo3::prelude::*; use pyo3::wrap_pymodule; -use quantity::python::PyInit_quantity; +use quantity::python::__PYO3_PYMODULE_DEF_QUANTITY; mod adsorption; mod fundamental_measure_theory; diff --git a/src/python/profile.rs b/src/python/profile.rs index 293679e..a73049b 100644 --- a/src/python/profile.rs +++ b/src/python/profile.rs @@ -117,14 +117,14 @@ macro_rules! impl_profile { /// Returns /// ------- /// SIArray - #[args(contributions = "PyContributions::Total()")] + #[args(contributions = "Contributions::Total")] #[pyo3(text_signature = "($self, contributions)")] fn entropy_density( &mut self, - contributions: PyContributions, + contributions: Contributions, ) -> PyResult<$si_arr> { Ok($si_arr::from( - self.0.profile.entropy_density(contributions.0)?, + self.0.profile.entropy_density(contributions)?, )) } @@ -139,14 +139,14 @@ macro_rules! impl_profile { /// Returns /// ------- /// SINumber - #[args(contributions = "PyContributions::Total()")] + #[args(contributions = "Contributions::Total")] #[pyo3(text_signature = "($self, contributions)")] fn entropy( &mut self, - contributions: PyContributions, + contributions: Contributions, ) -> PyResult { Ok(PySINumber::from( - self.0.profile.entropy(contributions.0)?, + self.0.profile.entropy(contributions)?, )) } @@ -161,14 +161,14 @@ macro_rules! impl_profile { /// Returns /// ------- /// SINumber - #[args(contributions = "PyContributions::Total()")] + #[args(contributions = "Contributions::Total")] #[pyo3(text_signature = "($self, contributions)")] fn internal_energy( &mut self, - contributions: PyContributions, + contributions: Contributions, ) -> PyResult { Ok(PySINumber::from( - self.0.profile.internal_energy(contributions.0)?, + self.0.profile.internal_energy(contributions)?, )) } } diff --git a/src/python/solver.rs b/src/python/solver.rs index d2a4be1..be2c03a 100644 --- a/src/python/solver.rs +++ b/src/python/solver.rs @@ -130,10 +130,7 @@ impl PyDFTSolver { fn _repr_markdown_(&self) -> String { self.0._repr_markdown_() } -} -#[pyproto] -impl pyo3::class::basic::PyObjectProtocol for PyDFTSolver { fn __repr__(&self) -> PyResult { Ok(self.0.to_string()) } From 86da06448913fa5fd5b4d6da8402a561d1c93f7e Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 10 Mar 2022 13:17:18 +0100 Subject: [PATCH 02/14] Add additional dual numbers for crit point calculation --- src/functional_contribution.rs | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/functional_contribution.rs b/src/functional_contribution.rs index d0bee33..b078820 100644 --- a/src/functional_contribution.rs +++ b/src/functional_contribution.rs @@ -34,10 +34,15 @@ macro_rules! impl_helmholtz_energy { impl_helmholtz_energy!(f64); impl_helmholtz_energy!(Dual64); +impl_helmholtz_energy!(Dual, f64>); impl_helmholtz_energy!(HyperDual64); impl_helmholtz_energy!(Dual3_64); impl_helmholtz_energy!(HyperDual); +impl_helmholtz_energy!(HyperDual, f64>); +impl_helmholtz_energy!(HyperDual, f64>); impl_helmholtz_energy!(Dual3); +impl_helmholtz_energy!(Dual3, f64>); +impl_helmholtz_energy!(Dual3, f64>); /// Individual functional contribution that can /// be evaluated using generalized (hyper) dual numbers. @@ -69,10 +74,15 @@ pub trait FunctionalContributionDual>: Display { pub trait FunctionalContribution: FunctionalContributionDual + FunctionalContributionDual + + FunctionalContributionDual, f64>> + FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + FunctionalContributionDual> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + Display { fn first_partial_derivatives( @@ -147,10 +157,15 @@ pub trait FunctionalContribution: impl FunctionalContribution for T where T: FunctionalContributionDual + FunctionalContributionDual + + FunctionalContributionDual, f64>> + FunctionalContributionDual + FunctionalContributionDual + FunctionalContributionDual> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + FunctionalContributionDual> + + FunctionalContributionDual, f64>> + + FunctionalContributionDual, f64>> + Display { } From 737c6036aaf37139dd02997964ade4afb2b48f6e Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 10 Mar 2022 14:06:46 +0100 Subject: [PATCH 03/14] change dependency to github --- Cargo.toml | 2 +- src/functional.rs | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 4d8ed43..adc3b64 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,7 +17,7 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ] [dependencies] quantity = { version = "0.5", features = ["linalg"] } -feos-core = { path = "../feos-core" } +feos-core = { git = "https://github.com/feos-org/feos-core" } num-dual = "0.5" ndarray = { version = "0.15", features = ["serde", "rayon"] } ndarray-stats = "0.5" diff --git a/src/functional.rs b/src/functional.rs index 15a882c..8bdad5a 100644 --- a/src/functional.rs +++ b/src/functional.rs @@ -296,8 +296,8 @@ impl DFT { Contributions::Total => { helmholtz_energy_density += &self.ideal_gas_contribution_dual::(temperature_dual, density); }, - Contributions::ResidualP|Contributions::IdealGas => panic!("Entropy density can only be calculated for Contributions::Residual or Contributions::Total"), - Contributions::Residual => (), + Contributions::ResidualNpt|Contributions::IdealGas => panic!("Entropy density can only be calculated for Contributions::Residual or Contributions::Total"), + Contributions::ResidualNvt => (), } Ok(helmholtz_energy_density.mapv(|f| -f.eps[0])) } @@ -367,8 +367,8 @@ impl DFT { Contributions::Total => { helmholtz_energy_density_dual += &self.ideal_gas_contribution_dual::(temperature_dual, density); }, - Contributions::ResidualP|Contributions::IdealGas => panic!("Internal energy density can only be calculated for Contributions::Residual or Contributions::Total"), - Contributions::Residual => (), + Contributions::ResidualNpt|Contributions::IdealGas => panic!("Internal energy density can only be calculated for Contributions::Residual or Contributions::Total"), + Contributions::ResidualNvt => (), } let helmholtz_energy_density = helmholtz_energy_density_dual .mapv(|f| f.re - f.eps[0] * temperature) From 10392b7ed0acc273fb100a3456878bb805f183c0 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 10 Mar 2022 16:29:26 +0100 Subject: [PATCH 04/14] moved creation of python module to build_wheel crate --- CHANGELOG.md | 8 +- Cargo.toml | 4 +- build_wheel/Cargo.toml | 7 +- build_wheel/src/lib.rs | 85 ++++++++++++++++++++- examples/FundamentalMeasureTheory.ipynb | 6 +- src/adsorption/external_potential.rs | 8 +- src/adsorption/fea_potential.rs | 18 ++--- src/adsorption/pore.rs | 24 +++--- src/convolver/mod.rs | 8 +- src/fundamental_measure_theory.rs | 1 + src/geometry.rs | 35 +++++---- src/lib.rs | 2 +- src/profile.rs | 13 +++- src/python/adsorption/external_potential.rs | 44 ----------- src/python/adsorption/mod.rs | 2 +- src/python/adsorption/pore.rs | 4 +- src/python/fundamental_measure_theory.rs | 77 ------------------- src/python/mod.rs | 44 +---------- src/python/solver.rs | 19 ++--- src/solver.rs | 72 +++++++---------- 20 files changed, 196 insertions(+), 285 deletions(-) delete mode 100644 src/python/fundamental_measure_theory.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index e2db505..8daf324 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,13 +6,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] -## [0.2.0] - 2022-03-09 +## [0.2.0] - 2022-03-?? +### Changed +- Renamed `AxisGeometry` to `Geometry`. +- Removed `PyGeometry` and `PyFMTVersion` in favor of a simpler implementation using `PyO3`'s new `#[pyclass]` for fieldless enums feature. +- `DFTSolver` now uses `Verbosity` instead of a `bool` to control its output. + ### 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 diff --git a/Cargo.toml b/Cargo.toml index adc3b64..848c35a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,13 +17,13 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ] [dependencies] quantity = { version = "0.5", features = ["linalg"] } -feos-core = { git = "https://github.com/feos-org/feos-core" } +#feos-core = { git = "https://github.com/feos-org/feos-core" } +feos-core = { path = "../feos-core" } 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.6" num-traits = "0.2" libc = "0.2" diff --git a/build_wheel/Cargo.toml b/build_wheel/Cargo.toml index e05a1ba..3d14965 100644 --- a/build_wheel/Cargo.toml +++ b/build_wheel/Cargo.toml @@ -1,13 +1,16 @@ [package] name = "feos_dft" version = "0.1.3" -authors = ["Philipp Rehner "] edition = "2018" [lib] crate-type = ["cdylib"] [dependencies] +quantity = "0.5" +#feos-core = { git = "https://github.com/feos-org/feos-core" } +feos-core = { path = "../../feos-core" } feos-dft = { path = "..", features = ["python"] } 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 index 69f7129..81bd936 100644 --- 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>); + +#[pymethods] +impl PyFMTFunctional { + #[new] + fn new(sigma: &PyArray1, version: FMTVersion) -> Self { + Self(Rc::new(FMTFunctional::new( + &sigma.to_owned_array(), + version, + ))) + } +} + +impl_equation_of_state!(PyFMTFunctional); + +impl_state!(DFT, 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::()?; + + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + + m.add_class::()?; + m.add_class::()?; + + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + + 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 index b7cba3d..17cd961 100644 --- 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.36 s, sys: 20.6 ms, total: 5.38 s\n", + "Wall time: 5.37 s\n" ] }, { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 3, diff --git a/src/adsorption/external_potential.rs b/src/adsorption/external_potential.rs index 5f84a72..2607d25 100644 --- a/src/adsorption/external_potential.rs +++ b/src/adsorption/external_potential.rs @@ -1,5 +1,5 @@ use crate::adsorption::fea_potential::calculate_fea_potential; -use crate::geometry::AxisGeometry; +use crate::geometry::Geometry; use feos_core::EosUnit; use libc::c_double; use ndarray::{Array1, Array2, Axis as Axis_nd}; @@ -184,7 +184,7 @@ impl ExternalPotential { system_size, n_grid, temperature, - AxisGeometry::Cartesian, + Geometry::Cartesian, ) } Self::Custom(_) => unreachable!(), @@ -330,7 +330,7 @@ impl ExternalPotential { system_size, n_grid, temperature, - AxisGeometry::Polar, + Geometry::Cylindrical, ) } Self::Custom(_) => unreachable!(), @@ -491,7 +491,7 @@ impl ExternalPotential { system_size, n_grid, temperature, - AxisGeometry::Spherical, + Geometry::Spherical, ) } Self::Custom(_) => unreachable!(), diff --git a/src/adsorption/fea_potential.rs b/src/adsorption/fea_potential.rs index 0a08a4a..da81444 100644 --- 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,7 +19,7 @@ pub fn calculate_fea_potential( system_size: &[QuantityScalar; 3], n_grid: &[usize; 2], temperature: f64, - geometry: AxisGeometry, + geometry: Geometry, ) -> Array1 { // allocate external potential let mut potential: Array1 = Array1::zeros(grid.len()); @@ -45,7 +45,7 @@ pub fn calculate_fea_potential( // 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 +54,7 @@ pub fn calculate_fea_potential( 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 +66,7 @@ pub fn calculate_fea_potential( // 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 +75,7 @@ pub fn calculate_fea_potential( 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 +97,13 @@ pub fn calculate_fea_potential( 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(), diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs index 3fba2c8..0db8e97 100644 --- a/src/adsorption/pore.rs +++ b/src/adsorption/pore.rs @@ -2,7 +2,7 @@ use crate::adsorption::{ExternalPotential, FluidParameters}; use crate::convolver::ConvolverFFT; use crate::functional::{HelmholtzEnergyFunctional, DFT}; use crate::functional_contribution::FunctionalContribution; -use crate::geometry::{Axis, AxisGeometry, Grid}; +use crate::geometry::{Axis, Geometry, Grid}; use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL}; use crate::solver::DFTSolver; use feos_core::{Contributions, EosResult, EosUnit, State, StateBuilder}; @@ -19,7 +19,7 @@ const DEFAULT_GRID_POINTS: usize = 2048; /// Parameters required to specify a 1D pore. pub struct Pore1D { // functional: Rc>, - geometry: AxisGeometry, + geometry: Geometry, pore_size: QuantityScalar, potential: ExternalPotential, n_grid: Option, @@ -28,7 +28,7 @@ pub struct Pore1D { impl Pore1D { pub fn new( - geometry: AxisGeometry, + geometry: Geometry, pore_size: QuantityScalar, potential: ExternalPotential, n_grid: Option, @@ -180,13 +180,13 @@ impl PoreSpecification for Pore1D { let n_grid = self.n_grid.unwrap_or(DEFAULT_GRID_POINTS); let axis = match self.geometry { - AxisGeometry::Cartesian => { + Geometry::Cartesian => { let potential_offset = POTENTIAL_OFFSET * bulk.eos.functional.sigma_ff().max().unwrap(); Axis::new_cartesian(n_grid, 0.5 * self.pore_size, Some(potential_offset))? } - AxisGeometry::Polar => Axis::new_polar(n_grid, self.pore_size)?, - AxisGeometry::Spherical => Axis::new_spherical(n_grid, self.pore_size)?, + Geometry::Cylindrical => Axis::new_polar(n_grid, self.pore_size)?, + Geometry::Spherical => Axis::new_spherical(n_grid, self.pore_size)?, }; // calculate external potential @@ -290,13 +290,13 @@ fn external_potential_1d( ) -> EosResult> { let potential_cutoff = potential_cutoff.unwrap_or(MAX_POTENTIAL); let effective_pore_size = match axis.geometry { - AxisGeometry::Spherical => pore_width.to_reduced(U::reference_length())?, - AxisGeometry::Polar => pore_width.to_reduced(U::reference_length())?, - AxisGeometry::Cartesian => 0.5 * pore_width.to_reduced(U::reference_length())?, + Geometry::Spherical => pore_width.to_reduced(U::reference_length())?, + Geometry::Cylindrical => pore_width.to_reduced(U::reference_length())?, + Geometry::Cartesian => 0.5 * pore_width.to_reduced(U::reference_length())?, }; let t = temperature.to_reduced(U::reference_temperature())?; let mut external_potential = match &axis.geometry { - AxisGeometry::Cartesian => { + Geometry::Cartesian => { potential.calculate_cartesian_potential( &(effective_pore_size + &axis.grid), fluid_parameters, @@ -307,13 +307,13 @@ fn external_potential_1d( t, ) } - AxisGeometry::Spherical => potential.calculate_spherical_potential( + Geometry::Spherical => potential.calculate_spherical_potential( &axis.grid, effective_pore_size, fluid_parameters, t, ), - AxisGeometry::Polar => potential.calculate_cylindrical_potential( + Geometry::Cylindrical => potential.calculate_cylindrical_potential( &axis.grid, effective_pore_size, fluid_parameters, diff --git a/src/convolver/mod.rs b/src/convolver/mod.rs index ac44ab9..9c49b09 100644 --- a/src/convolver/mod.rs +++ b/src/convolver/mod.rs @@ -1,4 +1,4 @@ -use crate::geometry::{Axis, AxisGeometry, Grid}; +use crate::geometry::{Axis, Geometry, Grid}; use crate::weight_functions::*; use ndarray::prelude::*; use ndarray::{Axis as Axis_nd, RemoveAxis, ScalarOperand, Slice}; @@ -132,9 +132,9 @@ where let mut lengths = Vec::with_capacity(cartesian_axes.len() + 1); let (transform, k_x) = match axis { Some(axis) => match axis.geometry { - AxisGeometry::Cartesian => CartesianTransform::new(axis), - AxisGeometry::Polar => PolarTransform::new(axis), - AxisGeometry::Spherical => SphericalTransform::new(axis), + Geometry::Cartesian => CartesianTransform::new(axis), + Geometry::Cylindrical => PolarTransform::new(axis), + Geometry::Spherical => SphericalTransform::new(axis), }, None => NoTransform::new(), }; diff --git a/src/fundamental_measure_theory.rs b/src/fundamental_measure_theory.rs index 675af62..7a8a1ab 100644 --- a/src/fundamental_measure_theory.rs +++ b/src/fundamental_measure_theory.rs @@ -23,6 +23,7 @@ pub trait FMTProperties { /// Different versions of fundamental measure theory #[derive(Clone, Copy)] +#[cfg_attr(feature = "python", pyo3::pyclass)] pub enum FMTVersion { /// White Bear ([Roth et al., 2002](https://doi.org/10.1088/0953-8984/14/46/313)) or modified ([Yu and Wu, 2002](https://doi.org/10.1063/1.1520530)) fundamental measure theory WhiteBear, diff --git a/src/geometry.rs b/src/geometry.rs index dfd9e7c..f4dae50 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -19,9 +19,9 @@ pub enum Grid { impl Grid { pub fn new_1d(axis: Axis) -> Self { match axis.geometry { - AxisGeometry::Cartesian => Self::Cartesian1(axis), - AxisGeometry::Polar => Self::Polar(axis), - AxisGeometry::Spherical => Self::Spherical(axis), + Geometry::Cartesian => Self::Cartesian1(axis), + Geometry::Cylindrical => Self::Polar(axis), + Geometry::Spherical => Self::Spherical(axis), } } @@ -66,18 +66,19 @@ impl Grid { /// Geometries of individual axes. #[derive(Copy, Clone)] -pub enum AxisGeometry { +#[cfg_attr(feature = "python", pyo3::pyclass)] +pub enum Geometry { Cartesian, - Polar, + Cylindrical, Spherical, } -impl AxisGeometry { +impl Geometry { /// Return the number of spatial dimensions for this geometry. pub fn dimension(&self) -> i32 { match self { Self::Cartesian => 1, - Self::Polar => 2, + Self::Cylindrical => 2, Self::Spherical => 3, } } @@ -86,7 +87,7 @@ impl AxisGeometry { /// An individual discretized axis. #[derive(Clone)] pub struct Axis { - pub geometry: AxisGeometry, + pub geometry: Geometry, pub grid: Array1, pub edges: Array1, integration_weights: Array1, @@ -110,7 +111,7 @@ impl Axis { let edges = Array1::linspace(0.0, l, points + 1); let integration_weights = Array1::from_elem(points, cell_size); Ok(Self { - geometry: AxisGeometry::Cartesian, + geometry: Geometry::Cartesian, grid, edges, integration_weights, @@ -128,7 +129,7 @@ impl Axis { 4.0 * FRAC_PI_3 * cell_size.powi(3) * (3 * k * k + 3 * k + 1) as f64 }); Ok(Self { - geometry: AxisGeometry::Spherical, + geometry: Geometry::Spherical, grid, edges, integration_weights, @@ -174,7 +175,7 @@ impl Axis { .collect(); Ok(Self { - geometry: AxisGeometry::Polar, + geometry: Geometry::Cylindrical, grid, edges, integration_weights, @@ -199,9 +200,9 @@ impl Axis { let length = (self.edges[self.grid.len()] - self.potential_offset - self.edges[0]) * U::reference_length(); (match self.geometry { - AxisGeometry::Cartesian => 1.0, - AxisGeometry::Polar => 4.0 * PI, - AxisGeometry::Spherical => 4.0 * FRAC_PI_3, + Geometry::Cartesian => 1.0, + Geometry::Cylindrical => 4.0 * PI, + Geometry::Spherical => 4.0 * FRAC_PI_3, }) * length.powi(self.geometry.dimension()) } @@ -219,10 +220,8 @@ impl Axis { n - 1 } else { match self.geometry { - AxisGeometry::Cartesian | AxisGeometry::Spherical => { - (x / self.edges[1]) as usize - } - AxisGeometry::Polar => { + Geometry::Cartesian | Geometry::Spherical => (x / self.edges[1]) as usize, + Geometry::Cylindrical => { if x < self.edges[1] { 0 } else { diff --git a/src/lib.rs b/src/lib.rs index 6ff73e9..3402b76 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -20,7 +20,7 @@ mod weight_functions; pub use convolver::{Convolver, ConvolverFFT}; pub use functional::{HelmholtzEnergyFunctional, DFT}; pub use functional_contribution::{FunctionalContribution, FunctionalContributionDual}; -pub use geometry::{Axis, AxisGeometry, Grid}; +pub use geometry::{Axis, Geometry, Grid}; pub use profile::{DFTProfile, DFTSpecification, DFTSpecifications}; pub use solver::DFTSolver; pub use weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape}; diff --git a/src/profile.rs b/src/profile.rs index a39e2f5..41432fa 100644 --- a/src/profile.rs +++ b/src/profile.rs @@ -3,8 +3,9 @@ use crate::functional::{HelmholtzEnergyFunctional, DFT}; use crate::geometry::Grid; use crate::solver::DFTSolver; use crate::weight_functions::WeightFunctionInfo; -use feos_core::{Contributions, EosError, EosResult, EosUnit, EquationOfState, State}; -use log::{info, warn}; +use feos_core::{ + log_result, Contributions, EosError, EosResult, EosUnit, EquationOfState, State, Verbosity, +}; use ndarray::{ s, Array, Array1, ArrayBase, ArrayViewMut, ArrayViewMut1, Axis as Axis_nd, Data, Dimension, Ix1, Ix2, Ix3, RemoveAxis, @@ -511,9 +512,13 @@ where // Call solver(s) let (converged, iterations) = solver.solve(&mut x, &mut residual)?; if converged { - info!("DFT solved in {} iterations", iterations); + log_result!(solver.verbosity, "DFT solved in {} iterations", iterations); } else if debug { - warn!("DFT not converged in {} iterations", iterations); + log_result!( + solver.verbosity, + "DFT not converged in {} iterations", + iterations + ); } else { return Err(EosError::NotConverged(String::from("DFT"))); } diff --git a/src/python/adsorption/external_potential.rs b/src/python/adsorption/external_potential.rs index 375952c..6d595b7 100644 --- a/src/python/adsorption/external_potential.rs +++ b/src/python/adsorption/external_potential.rs @@ -1,5 +1,4 @@ use crate::adsorption::ExternalPotential; -use crate::geometry::AxisGeometry; use numpy::PyArray1; use pyo3::prelude::*; use quantity::python::{PySIArray2, PySINumber}; @@ -216,46 +215,3 @@ impl PyExternalPotential { }) } } - -/// Geometry of a one-dimensional pore. -/// -/// Returns -/// ------- -/// Geometry -#[pyclass(name = "Geometry")] -#[derive(Clone)] -pub struct PyGeometry(pub AxisGeometry); - -#[pymethods] -#[allow(non_snake_case)] -impl PyGeometry { - /// Cartesian coordinates. - /// - /// Returns - /// ------- - /// AxisGeometry - #[classattr] - pub fn Cartesian() -> Self { - Self(AxisGeometry::Cartesian) - } - - /// Cylindrical coordinates. - /// - /// Returns - /// ------- - /// AxisGeometry - #[classattr] - pub fn Cylindrical() -> Self { - Self(AxisGeometry::Polar) - } - - /// Spherical coordinates. - /// - /// Returns - /// ------- - /// AxisGeometry - #[classattr] - pub fn Spherical() -> Self { - Self(AxisGeometry::Spherical) - } -} diff --git a/src/python/adsorption/mod.rs b/src/python/adsorption/mod.rs index e89e505..fbfe204 100644 --- a/src/python/adsorption/mod.rs +++ b/src/python/adsorption/mod.rs @@ -1,7 +1,7 @@ mod external_potential; mod pore; -pub use external_potential::{PyExternalPotential, PyGeometry}; +pub use external_potential::PyExternalPotential; #[macro_export] macro_rules! impl_adsorption { diff --git a/src/python/adsorption/pore.rs b/src/python/adsorption/pore.rs index b26611d..10a7c6f 100644 --- a/src/python/adsorption/pore.rs +++ b/src/python/adsorption/pore.rs @@ -34,14 +34,14 @@ macro_rules! impl_pore { impl PyPore1D { #[new] fn new( - geometry: PyGeometry, + geometry: Geometry, pore_size: PySINumber, potential: PyExternalPotential, n_grid: Option, potential_cutoff: Option, ) -> Self { Self(Pore1D::new( - geometry.0, + geometry, pore_size.into(), potential.0, n_grid, diff --git a/src/python/fundamental_measure_theory.rs b/src/python/fundamental_measure_theory.rs deleted file mode 100644 index b111a2a..0000000 --- a/src/python/fundamental_measure_theory.rs +++ /dev/null @@ -1,77 +0,0 @@ -use super::{PyDFTSolver, PyExternalPotential, PyGeometry}; -use crate::adsorption::*; -use crate::functional::DFT; -use crate::fundamental_measure_theory::{FMTFunctional, FMTVersion}; -use crate::solvation::*; -use crate::*; -use feos_core::*; -use numpy::*; -use pyo3::exceptions::PyValueError; -use pyo3::prelude::*; -use quantity::python::*; -use quantity::si::*; -use std::rc::Rc; - -/// Different versions of fundamental measure theory -#[pyclass(name = "FMTVersion")] -#[derive(Clone, Copy)] -pub struct PyFMTVersion(pub FMTVersion); - -#[pymethods] -#[allow(non_snake_case)] -impl PyFMTVersion { - /// White Bear ([Roth et al., 2002](https://doi.org/10.1088/0953-8984/14/46/313)) or modified ([Yu and Wu, 2002](https://doi.org/10.1063/1.1520530)) fundamental measure theory - #[classattr] - pub fn WhiteBear() -> Self { - Self(FMTVersion::WhiteBear) - } - - /// Scalar fundamental measure theory by [Kierlik and Rosinberg, 1990](https://doi.org/10.1103/PhysRevA.42.3382) - #[classattr] - pub fn KierlikRosinberg() -> Self { - Self(FMTVersion::KierlikRosinberg) - } - - /// Anti-symmetric White Bear fundamental measure theory ([Rosenfeld et al., 1997](https://doi.org/10.1103/PhysRevE.55.4245)) and SI of ([Kessler et al., 2021](https://doi.org/10.1016/j.micromeso.2021.111263)) - #[classattr] - pub fn AntiSymWhiteBear() -> Self { - Self(FMTVersion::AntiSymWhiteBear) - } -} - -/// 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>); - -#[pymethods] -impl PyFMTFunctional { - #[new] - fn new(sigma: &PyArray1, version: PyFMTVersion) -> Self { - Self(Rc::new(FMTFunctional::new( - &sigma.to_owned_array(), - version.0, - ))) - } -} - -impl_equation_of_state!(PyFMTFunctional); - -impl_state!(DFT, PyFMTFunctional); - -impl_pore!(FMTFunctional, PyFMTFunctional); -impl_adsorption!(FMTFunctional, PyFMTFunctional); - -impl_pair_correlation!(FMTFunctional); diff --git a/src/python/mod.rs b/src/python/mod.rs index 2db7297..9f98cda 100644 --- a/src/python/mod.rs +++ b/src/python/mod.rs @@ -1,50 +1,8 @@ -use pyo3::prelude::*; -use pyo3::wrap_pymodule; -use quantity::python::__PYO3_PYMODULE_DEF_QUANTITY; - mod adsorption; -mod fundamental_measure_theory; mod interface; mod profile; mod solvation; mod solver; -pub use adsorption::{PyExternalPotential, PyGeometry}; -pub use fundamental_measure_theory::PyFMTVersion; -use fundamental_measure_theory::*; +pub use adsorption::PyExternalPotential; pub use solver::PyDFTSolver; - -#[pymodule] -pub fn feos_dft(py: Python<'_>, m: &PyModule) -> PyResult<()> { - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - - m.add_class::()?; - m.add_class::()?; - - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - - 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/src/python/solver.rs b/src/python/solver.rs index be2c03a..49fbedb 100644 --- a/src/python/solver.rs +++ b/src/python/solver.rs @@ -1,32 +1,27 @@ use crate::DFTSolver; +use feos_core::Verbosity; use pyo3::prelude::*; /// Settings for the DFT solver. /// /// Parameters /// ---------- -/// output: bool, optional -/// Print the progress to the console. +/// verbosity: Verbosity, optional +/// The verbosity level of the solver. /// /// Returns /// ------- -/// empty solver: DFTSolver +/// DFTSolver #[pyclass(name = "DFTSolver")] #[derive(Clone)] -#[pyo3(text_signature = "(output=None)")] +#[pyo3(text_signature = "(verbosity=None)")] pub struct PyDFTSolver(pub DFTSolver); #[pymethods] impl PyDFTSolver { #[new] - fn new(output: Option) -> Self { - let mut solver = DFTSolver::new(); - if let Some(output) = output { - if output { - solver = solver.output(); - } - } - Self(solver) + fn new(verbosity: Option) -> Self { + Self(DFTSolver::new(verbosity.unwrap_or_default())) } /// The default solver. diff --git a/src/solver.rs b/src/solver.rs index 4ae52f1..68ffdd7 100644 --- a/src/solver.rs +++ b/src/solver.rs @@ -1,4 +1,4 @@ -use feos_core::{EosError, EosResult}; +use feos_core::{log_iter, EosError, EosResult, Verbosity}; use ndarray::prelude::*; use num_dual::linalg::{norm, LU}; use std::collections::VecDeque; @@ -45,24 +45,24 @@ enum DFTAlgorithm { #[derive(Clone)] pub struct DFTSolver { parameters: Vec, - output: bool, + pub verbosity: Verbosity, } impl Default for DFTSolver { fn default() -> Self { Self { parameters: vec![DEFAULT_PARAMS_ANDERSON_LOG, DEFAULT_PARAMS_ANDERSON], - output: false, + verbosity: Verbosity::None, } } } impl DFTSolver { /// Create a new empty `DFTSolver` object. - pub fn new() -> Self { + pub fn new(verbosity: Verbosity) -> Self { Self { parameters: Vec::new(), - output: false, + verbosity, } } @@ -110,23 +110,15 @@ impl DFTSolver { self } - /// Print the iteration to the console. - pub fn output(mut self) -> Self { - self.output = true; - self - } - pub(crate) fn solve(&self, x: &mut Array1, residual: &mut F) -> EosResult<(bool, usize)> where F: FnMut(&Array1, ArrayViewMut1, bool) -> EosResult<()>, { - if self.output { - println!("solver | iter | residual "); - } + log_iter!(self.verbosity, "solver | iter | residual "); let mut converged = false; let mut iterations = 0; for algorithm in &self.parameters { - let (c, i) = algorithm.solve(x, residual, self.output)?; + let (c, i) = algorithm.solve(x, residual, self.verbosity)?; converged = c; iterations += i; } @@ -139,16 +131,16 @@ impl SolverParameter { &self, x: &mut Array1, residual: &mut F, - output: bool, + verbosity: Verbosity, ) -> EosResult<(bool, usize)> where F: FnMut(&Array1, ArrayViewMut1, bool) -> EosResult<()>, { match self.solver { DFTAlgorithm::PicardIteration(max_rel) => { - self.solve_picard(max_rel, x, residual, output) + self.solve_picard(max_rel, x, residual, verbosity) } - DFTAlgorithm::AndersonMixing(mmax) => self.solve_anderson(mmax, x, residual, output), + DFTAlgorithm::AndersonMixing(mmax) => self.solve_anderson(mmax, x, residual, verbosity), } } @@ -157,14 +149,12 @@ impl SolverParameter { max_rel: f64, x: &mut Array1, residual: &mut F, - output: bool, + verbosity: Verbosity, ) -> EosResult<(bool, usize)> where F: FnMut(&Array1, ArrayViewMut1, bool) -> EosResult<()>, { - if output { - println!("{:-<43}", ""); - } + log_iter!(verbosity, "{:-<43}", ""); let mut resm = Array::zeros(x.raw_dim()); for k in 1..=self.max_iter { @@ -191,15 +181,14 @@ impl SolverParameter { // check for convergence let res = norm(&resm) / (resm.len() as f64).sqrt(); - if output { - println!( - "Picard iteration {:3} | {:>4} | {:.6e} | {}", - if self.log { "log" } else { "" }, - k, - res, - beta_min.unwrap_or(self.beta) - ); - } + log_iter!( + verbosity, + "Picard iteration {:3} | {:>4} | {:.6e} | {}", + if self.log { "log" } else { "" }, + k, + res, + beta_min.unwrap_or(self.beta) + ); if res.is_nan() { return Err(EosError::IterationFailed(String::from("Picard Iteration"))); @@ -216,14 +205,12 @@ impl SolverParameter { mmax: usize, x: &mut Array1, residual: &mut F, - output: bool, + verbosity: Verbosity, ) -> EosResult<(bool, usize)> where F: FnMut(&Array1, ArrayViewMut1, bool) -> EosResult<()>, { - if output { - println!("{:-<43}", ""); - } + log_iter!(verbosity, "{:-<43}", ""); let mut resm = VecDeque::with_capacity(mmax); let mut xm = VecDeque::with_capacity(mmax); let mut r; @@ -274,14 +261,13 @@ impl SolverParameter { // check for convergence let resv = &resm[m - 1]; let res = norm(resv) / (resv.len() as f64).sqrt(); - if output { - println!( - "Anderson mixing {:3} | {:>4} | {:.6e} ", - if self.log { "log" } else { "" }, - k, - res - ); - } + log_iter!( + verbosity, + "Anderson mixing {:3} | {:>4} | {:.6e} ", + if self.log { "log" } else { "" }, + k, + res + ); if res.is_nan() { return Err(EosError::IterationFailed(String::from("Anderson Mixing"))); From fe6407afc1a5e3162c0217ecdf6166f7f045bb9a Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 10 Mar 2022 16:47:02 +0100 Subject: [PATCH 05/14] back to github --- Cargo.toml | 3 +-- build_wheel/Cargo.toml | 3 +-- examples/FundamentalMeasureTheory.ipynb | 6 +++--- src/python/solver.rs | 12 ++++++++++++ 4 files changed, 17 insertions(+), 7 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 848c35a..cb9708d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,8 +17,7 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ] [dependencies] quantity = { version = "0.5", features = ["linalg"] } -#feos-core = { git = "https://github.com/feos-org/feos-core" } -feos-core = { path = "../feos-core" } +feos-core = { git = "https://github.com/feos-org/feos-core" } num-dual = "0.5" ndarray = { version = "0.15", features = ["serde", "rayon"] } ndarray-stats = "0.5" diff --git a/build_wheel/Cargo.toml b/build_wheel/Cargo.toml index 3d14965..64d6cf3 100644 --- a/build_wheel/Cargo.toml +++ b/build_wheel/Cargo.toml @@ -9,8 +9,7 @@ crate-type = ["cdylib"] [dependencies] quantity = "0.5" -#feos-core = { git = "https://github.com/feos-org/feos-core" } -feos-core = { path = "../../feos-core" } +feos-core = { git = "https://github.com/feos-org/feos-core" } feos-dft = { path = "..", features = ["python"] } pyo3 = { version = "0.16", features = ["extension-module", "abi3", "abi3-py37"] } numpy = "0.16" diff --git a/examples/FundamentalMeasureTheory.ipynb b/examples/FundamentalMeasureTheory.ipynb index 17cd961..1b79b9a 100644 --- a/examples/FundamentalMeasureTheory.ipynb +++ b/examples/FundamentalMeasureTheory.ipynb @@ -39,14 +39,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "CPU times: user 5.36 s, sys: 20.6 ms, total: 5.38 s\n", - "Wall time: 5.37 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": [ - "" + "" ] }, "execution_count": 3, diff --git a/src/python/solver.rs b/src/python/solver.rs index 49fbedb..0195918 100644 --- a/src/python/solver.rs +++ b/src/python/solver.rs @@ -77,6 +77,18 @@ impl PyDFTSolver { Self(solver) } + fn log_iter(&self) -> Self { + let mut solver = self.0.clone(); + solver.verbosity = Verbosity::Iter; + Self(solver) + } + + fn log_result(&self) -> Self { + let mut solver = self.0.clone(); + solver.verbosity = Verbosity::Result; + Self(solver) + } + /// Add Anderson mixing to the solver object. /// /// Parameters From c71605e42d5087de02e0ea76450e13514ab72f7d Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 25 Mar 2022 19:25:52 +0100 Subject: [PATCH 06/14] use StateVec in SurfaceTensionDiagram --- CHANGELOG.md | 7 +++-- Cargo.toml | 2 +- build_wheel/Cargo.toml | 2 +- src/interface/surface_tension_diagram.rs | 39 +++--------------------- 4 files changed, 11 insertions(+), 39 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8daf324..644d245 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,9 +8,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.2.0] - 2022-03-?? ### Changed -- Renamed `AxisGeometry` to `Geometry`. -- Removed `PyGeometry` and `PyFMTVersion` in favor of a simpler implementation using `PyO3`'s new `#[pyclass]` for fieldless enums feature. -- `DFTSolver` now uses `Verbosity` instead of a `bool` to control its output. +- 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) ### Packaging - Updated `pyo3` and `numpy` dependencies to 0.16. diff --git a/Cargo.toml b/Cargo.toml index cb9708d..e27e8a4 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,7 +17,7 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ] [dependencies] quantity = { version = "0.5", features = ["linalg"] } -feos-core = { git = "https://github.com/feos-org/feos-core" } +feos-core = { git = "https://github.com/feos-org/feos-core", branch = "state_vec" } num-dual = "0.5" ndarray = { version = "0.15", features = ["serde", "rayon"] } ndarray-stats = "0.5" diff --git a/build_wheel/Cargo.toml b/build_wheel/Cargo.toml index 64d6cf3..e765114 100644 --- a/build_wheel/Cargo.toml +++ b/build_wheel/Cargo.toml @@ -9,7 +9,7 @@ crate-type = ["cdylib"] [dependencies] quantity = "0.5" -feos-core = { git = "https://github.com/feos-org/feos-core" } +feos-core = { git = "https://github.com/feos-org/feos-core", branch = "state_vec" } feos-dft = { path = "..", features = ["python"] } pyo3 = { version = "0.16", features = ["extension-module", "abi3", "abi3-py37"] } numpy = "0.16" diff --git a/src/interface/surface_tension_diagram.rs b/src/interface/surface_tension_diagram.rs index e89ed35..28f4e6e 100644 --- a/src/interface/surface_tension_diagram.rs +++ b/src/interface/surface_tension_diagram.rs @@ -1,8 +1,7 @@ use super::PlanarInterface; use crate::functional::{HelmholtzEnergyFunctional, DFT}; use crate::solver::DFTSolver; -use feos_core::{Contributions, EosUnit, EquationOfState, PhaseEquilibrium}; -use ndarray::Array1; +use feos_core::{EosUnit, PhaseEquilibrium, StateVec}; use quantity::{QuantityArray1, QuantityScalar}; const DEFAULT_GRID_POINTS: usize = 2048; @@ -64,40 +63,12 @@ impl SurfaceTensionDiagram { Self { profiles } } - pub fn temperature(&self) -> QuantityArray1 { - QuantityArray1::from_shape_fn(self.profiles.len(), |i| { - self.profiles[i].profile.temperature - }) - } - - pub fn pressure(&self) -> QuantityArray1 { - QuantityArray1::from_shape_fn(self.profiles.len(), |i| { - self.profiles[i].vle.vapor().pressure(Contributions::Total) - }) + pub fn vapor(&self) -> StateVec<'_, U, DFT> { + self.profiles.iter().map(|p| p.vle.vapor()).collect() } - pub fn vapor_molefracs(&self) -> Array1 { - let mut x: Array1 = self - .profiles - .iter() - .map(|p| p.vle.vapor().molefracs[0]) - .collect(); - if self.profiles[0].vle.vapor().eos.components() == 1 { - x[0] = 0.0; - } - x - } - - pub fn liquid_molefracs(&self) -> Array1 { - let mut x: Array1 = self - .profiles - .iter() - .map(|p| p.vle.liquid().molefracs[0]) - .collect(); - if self.profiles[0].vle.liquid().eos.components() == 1 { - x[0] = 0.0; - } - x + pub fn liquid(&self) -> StateVec<'_, U, DFT> { + self.profiles.iter().map(|p| p.vle.liquid()).collect() } pub fn surface_tension(&mut self) -> QuantityArray1 { From 2bff8980c90425830a205bea35bbe7db7f744a58 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 25 Mar 2022 20:10:04 +0100 Subject: [PATCH 07/14] fix python macro --- .../interface/surface_tension_diagram.rs | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/src/python/interface/surface_tension_diagram.rs b/src/python/interface/surface_tension_diagram.rs index 2261b65..f45b9c5 100644 --- a/src/python/interface/surface_tension_diagram.rs +++ b/src/python/interface/surface_tension_diagram.rs @@ -62,23 +62,13 @@ macro_rules! impl_surface_tension_diagram { } #[getter] - pub fn get_temperature(&self) -> PySIArray1 { - self.0.temperature().into() + pub fn get_vapor(&self) -> PyStateVec { + self.0.vapor().into() } #[getter] - pub fn get_pressure(&self) -> PySIArray1 { - self.0.pressure().into() - } - - #[getter] - fn get_vapor_molefracs<'py>(&self, py: Python<'py>) -> &'py PyArray1 { - self.0.vapor_molefracs().view().to_pyarray(py) - } - - #[getter] - fn get_liquid_molefracs<'py>(&self, py: Python<'py>) -> &'py PyArray1 { - self.0.liquid_molefracs().view().to_pyarray(py) + pub fn get_liquid(&self) -> PyStateVec { + self.0.liquid().into() } #[getter] From 82a6a16805fa1f66ede6cc2a00d2699466761537 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sat, 26 Mar 2022 11:22:17 +0100 Subject: [PATCH 08/14] update Cargo.toml --- Cargo.toml | 2 +- build_wheel/Cargo.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index e27e8a4..cb9708d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,7 +17,7 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ] [dependencies] quantity = { version = "0.5", features = ["linalg"] } -feos-core = { git = "https://github.com/feos-org/feos-core", branch = "state_vec" } +feos-core = { git = "https://github.com/feos-org/feos-core" } num-dual = "0.5" ndarray = { version = "0.15", features = ["serde", "rayon"] } ndarray-stats = "0.5" diff --git a/build_wheel/Cargo.toml b/build_wheel/Cargo.toml index e765114..64d6cf3 100644 --- a/build_wheel/Cargo.toml +++ b/build_wheel/Cargo.toml @@ -9,7 +9,7 @@ crate-type = ["cdylib"] [dependencies] quantity = "0.5" -feos-core = { git = "https://github.com/feos-org/feos-core", branch = "state_vec" } +feos-core = { git = "https://github.com/feos-org/feos-core" } feos-dft = { path = "..", features = ["python"] } pyo3 = { version = "0.16", features = ["extension-module", "abi3", "abi3-py37"] } numpy = "0.16" From 30300daf1f1e9d6fb1845d743757725b680a78fb Mon Sep 17 00:00:00 2001 From: Philipp Rehner <69816385+prehner@users.noreply.github.com> Date: Thu, 31 Mar 2022 18:53:24 +0200 Subject: [PATCH 09/14] Add density as additional parameter to DFTProfile::new() (#24) * Add density as additional parameter to DFTProfile::new() * Fix calculation of pore volume * update changelog --- CHANGELOG.md | 3 ++- src/adsorption/mod.rs | 22 ++++++++++++-------- src/adsorption/pore.rs | 13 +++++++----- src/interface/mod.rs | 2 +- src/profile.rs | 34 ++++++++++++++++++------------- src/python/adsorption/pore.rs | 12 +++++++++-- src/solvation/mod.rs | 2 +- src/solvation/pair_correlation.rs | 2 +- 8 files changed, 57 insertions(+), 33 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 644d245..3e6ea38 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,12 +6,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] -## [0.2.0] - 2022-03-?? +## [0.2.0] - 2022-04-?? ### 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) ### Packaging - Updated `pyo3` and `numpy` dependencies to 0.16. diff --git a/src/adsorption/mod.rs b/src/adsorption/mod.rs index c5027fd..4b1a2b7 100644 --- a/src/adsorption/mod.rs +++ b/src/adsorption/mod.rs @@ -265,7 +265,10 @@ where .vapor() .clone(); } - let external_potential = pore.initialize(&bulk, None)?.profile.external_potential; + let external_potential = pore + .initialize(&bulk, None, None)? + .profile + .external_potential; for i in 0..pressure.len() { let mut bulk = StateBuilder::new(functional) @@ -279,11 +282,14 @@ where .vapor() .clone(); } - let mut p = pore.initialize(&bulk, Some(&external_potential))?; - let p2 = p.clone(); - if let Some(Ok(l)) = profiles.last() { - p.profile.density = l.profile.density.clone(); - } + let old_density = if let Some(Ok(l)) = profiles.last() { + Some(&l.profile.density) + } else { + None + }; + + let p = pore.initialize(&bulk, old_density, Some(&external_potential))?; + let p2 = pore.initialize(&bulk, None, Some(&external_potential))?; profiles.push(p.solve(solver).or_else(|_| p2.solve(solver))); } @@ -318,8 +324,8 @@ where .liquid() .build()?; - let mut vapor = pore.initialize(&vapor_bulk, None)?.solve(None)?; - let mut liquid = pore.initialize(&liquid_bulk, None)?.solve(solver)?; + let mut vapor = pore.initialize(&vapor_bulk, None, None)?.solve(None)?; + let mut liquid = pore.initialize(&liquid_bulk, None, None)?.solve(solver)?; // calculate initial value for the molar gibbs energy let nv = vapor.profile.bulk.density diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs index 0db8e97..eecd647 100644 --- a/src/adsorption/pore.rs +++ b/src/adsorption/pore.rs @@ -10,7 +10,7 @@ use ndarray::prelude::*; use ndarray::Axis as Axis_nd; use ndarray::Zip; use ndarray_stats::QuantileExt; -use quantity::{QuantityArray2, QuantityScalar}; +use quantity::{QuantityArray, QuantityArray2, QuantityArray4, QuantityScalar}; use std::rc::Rc; const POTENTIAL_OFFSET: f64 = 2.0; @@ -83,6 +83,7 @@ pub trait PoreSpecification { fn initialize( &self, bulk: &State>, + density: Option<&QuantityArray>, external_potential: Option<&Array>, ) -> EosResult>; @@ -96,9 +97,9 @@ pub trait PoreSpecification { { let bulk = StateBuilder::new(&Rc::new(Helium::new())) .temperature(298.0 * U::reference_temperature()) - .volume(U::reference_volume()) + .density(U::reference_density()) .build()?; - let pore = self.initialize(&bulk, None)?; + let pore = self.initialize(&bulk, None, None)?; let pot = pore .profile .external_potential @@ -174,6 +175,7 @@ impl PoreSpecification for Pore1D { fn initialize( &self, bulk: &State>, + density: Option<&QuantityArray2>, external_potential: Option<&Array2>, ) -> EosResult> { let dft = &bulk.eos; @@ -211,7 +213,7 @@ impl PoreSpecification for Pore1D { let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1)); Ok(PoreProfile { - profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?, + profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), density)?, grand_potential: None, interfacial_tension: None, }) @@ -226,6 +228,7 @@ impl PoreSpecification for Pore3D { fn initialize( &self, bulk: &State>, + density: Option<&QuantityArray4>, external_potential: Option<&Array4>, ) -> EosResult> { let dft = &bulk.eos; @@ -269,7 +272,7 @@ impl PoreSpecification for Pore3D { let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1)); Ok(PoreProfile { - profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?, + profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), density)?, grand_potential: None, interfacial_tension: None, }) diff --git a/src/interface/mod.rs b/src/interface/mod.rs index e2c7b66..d918bb2 100644 --- a/src/interface/mod.rs +++ b/src/interface/mod.rs @@ -82,7 +82,7 @@ impl PlanarInterface { let convolver = ConvolverFFT::plan(&grid, &weight_functions, None); Ok(Self { - profile: DFTProfile::new(grid, convolver, vle.vapor(), None)?, + profile: DFTProfile::new(grid, convolver, vle.vapor(), None, None)?, vle: vle.clone(), surface_tension: None, equimolar_radius: None, diff --git a/src/profile.rs b/src/profile.rs index 41432fa..8345407 100644 --- a/src/profile.rs +++ b/src/profile.rs @@ -189,6 +189,7 @@ where convolver: Rc>, bulk: &State>, external_potential: Option>, + density: Option<&QuantityArray>, ) -> EosResult { let dft = bulk.eos.clone(); @@ -201,26 +202,31 @@ where Array::zeros(n_grid).into_dimensionality().unwrap() }); - // intitialize density - let t = bulk.temperature.to_reduced(U::reference_temperature())?; - let bonds = dft - .bond_integrals(t, &external_potential, &convolver) - .mapv(f64::abs) - * (-&external_potential).mapv(f64::exp); - let mut density = Array::zeros(external_potential.raw_dim()); - let bulk_density = bulk.partial_density.to_reduced(U::reference_density())?; - for (s, &c) in dft.component_index.iter().enumerate() { - density - .index_axis_mut(Axis_nd(0), s) - .assign(&(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c])); - } + // initialize density + let density = if let Some(density) = density { + density.clone() + } else { + let t = bulk.temperature.to_reduced(U::reference_temperature())?; + let bonds = dft + .bond_integrals(t, &external_potential, &convolver) + .mapv(f64::abs) + * (-&external_potential).mapv(f64::exp); + let mut density = Array::zeros(external_potential.raw_dim()); + let bulk_density = bulk.partial_density.to_reduced(U::reference_density())?; + for (s, &c) in dft.component_index.iter().enumerate() { + density.index_axis_mut(Axis_nd(0), s).assign( + &(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]), + ); + } + density * U::reference_density() + }; Ok(Self { grid, convolver, dft: bulk.eos.clone(), temperature: bulk.temperature, - density: density * U::reference_density(), + density, specification: Rc::new(DFTSpecifications::ChemicalPotential), external_potential, bulk: bulk.clone(), diff --git a/src/python/adsorption/pore.rs b/src/python/adsorption/pore.rs index 10a7c6f..6f6f4f6 100644 --- a/src/python/adsorption/pore.rs +++ b/src/python/adsorption/pore.rs @@ -55,6 +55,8 @@ macro_rules! impl_pore { /// ---------- /// bulk : State /// The bulk state in equilibrium with the pore. + /// density : SIArray2, optional + /// Initial values for the density profile. /// external_potential : numpy.ndarray[float], optional /// The external potential in the pore. Used to /// save computation time in the case of costly @@ -63,14 +65,16 @@ macro_rules! impl_pore { /// Returns /// ------- /// PoreProfile1D - #[pyo3(text_signature = "($self, bulk, external_potential=None)")] + #[pyo3(text_signature = "($self, bulk, density=None, external_potential=None)")] fn initialize( &self, bulk: &PyState, + density: Option, external_potential: Option<&PyArray2>, ) -> PyResult { Ok(PyPoreProfile1D(self.0.initialize( &bulk.0, + density.as_deref(), external_potential.map(|e| e.to_owned_array()).as_ref(), )?)) } @@ -156,6 +160,8 @@ macro_rules! impl_pore { /// ---------- /// bulk : State /// The bulk state in equilibrium with the pore. + /// density : SIArray4, optional + /// Initial values for the density profile. /// external_potential : numpy.ndarray[float], optional /// The external potential in the pore. Used to /// save computation time in the case of costly @@ -164,14 +170,16 @@ macro_rules! impl_pore { /// Returns /// ------- /// PoreProfile3D - #[pyo3(text_signature = "($self, bulk, external_potential=None)")] + #[pyo3(text_signature = "($self, bulk, density=None, external_potential=None)")] fn initialize( &self, bulk: &PyState, + density: Option, external_potential: Option<&PyArray4>, ) -> PyResult { Ok(PyPoreProfile3D(self.0.initialize( &bulk.0, + density.as_deref(), external_potential.map(|e| e.to_owned_array()).as_ref(), )?)) } diff --git a/src/solvation/mod.rs b/src/solvation/mod.rs index 68a6688..3699edc 100644 --- a/src/solvation/mod.rs +++ b/src/solvation/mod.rs @@ -120,7 +120,7 @@ impl SolvationProfil let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1)); Ok(Self { - profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential))?, + profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), None)?, grand_potential: None, solvation_free_energy: None, }) diff --git a/src/solvation/pair_correlation.rs b/src/solvation/pair_correlation.rs index d9c6a5b..c6c39a6 100644 --- a/src/solvation/pair_correlation.rs +++ b/src/solvation/pair_correlation.rs @@ -59,7 +59,7 @@ impl PairCorrelation Date: Fri, 8 Apr 2022 16:29:28 +0200 Subject: [PATCH 10/14] Less redundant approach to the `DFT` wrapper struct (#27) --- CHANGELOG.md | 2 + src/adsorption/external_potential.rs | 4 +- src/adsorption/pore.rs | 27 +++-- src/functional.rs | 127 ++++++++++++++--------- src/fundamental_measure_theory.rs | 24 ++--- src/interface/mod.rs | 12 +-- src/interface/surface_tension_diagram.rs | 2 +- src/lib.rs | 2 +- src/pdgt.rs | 4 +- src/profile.rs | 20 ++-- src/solvation/mod.rs | 6 +- src/solvation/pair_correlation.rs | 4 +- 12 files changed, 129 insertions(+), 105 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3e6ea38..a795c53 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `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) ### Packaging - Updated `pyo3` and `numpy` dependencies to 0.16. diff --git a/src/adsorption/external_potential.rs b/src/adsorption/external_potential.rs index 2607d25..d23a9ec 100644 --- a/src/adsorption/external_potential.rs +++ b/src/adsorption/external_potential.rs @@ -1,4 +1,5 @@ use crate::adsorption::fea_potential::calculate_fea_potential; +use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::Geometry; use feos_core::EosUnit; use libc::c_double; @@ -55,10 +56,9 @@ pub enum ExternalPotential { } /// Parameters of the fluid required to evaluate the external potential. -pub trait FluidParameters { +pub trait FluidParameters: HelmholtzEnergyFunctional { fn epsilon_k_ff(&self) -> Array1; fn sigma_ff(&self) -> &Array1; - fn m(&self) -> Array1; } impl ExternalPotential { diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs index eecd647..cd1a849 100644 --- a/src/adsorption/pore.rs +++ b/src/adsorption/pore.rs @@ -1,6 +1,6 @@ use crate::adsorption::{ExternalPotential, FluidParameters}; use crate::convolver::ConvolverFFT; -use crate::functional::{HelmholtzEnergyFunctional, DFT}; +use crate::functional::{HelmholtzEnergyFunctional, MoleculeShape, DFT}; use crate::functional_contribution::FunctionalContribution; use crate::geometry::{Axis, Geometry, Grid}; use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL}; @@ -178,13 +178,12 @@ impl PoreSpecification for Pore1D { density: Option<&QuantityArray2>, external_potential: Option<&Array2>, ) -> EosResult> { - let dft = &bulk.eos; + let dft: &F = &bulk.eos; let n_grid = self.n_grid.unwrap_or(DEFAULT_GRID_POINTS); let axis = match self.geometry { Geometry::Cartesian => { - let potential_offset = - POTENTIAL_OFFSET * bulk.eos.functional.sigma_ff().max().unwrap(); + let potential_offset = POTENTIAL_OFFSET * bulk.eos.sigma_ff().max().unwrap(); Axis::new_cartesian(n_grid, 0.5 * self.pore_size, Some(potential_offset))? } Geometry::Cylindrical => Axis::new_polar(n_grid, self.pore_size)?, @@ -198,7 +197,7 @@ impl PoreSpecification for Pore1D { self.pore_size, bulk.temperature, &self.potential, - &bulk.eos.functional, + dft, &axis, self.potential_cutoff, ) @@ -209,7 +208,7 @@ impl PoreSpecification for Pore1D { // initialize convolver let grid = Grid::new_1d(axis); let t = bulk.temperature.to_reduced(U::reference_temperature())?; - let weight_functions = dft.functional.weight_functions(t); + let weight_functions = dft.weight_functions(t); let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1)); Ok(PoreProfile { @@ -231,7 +230,7 @@ impl PoreSpecification for Pore3D { density: Option<&QuantityArray4>, external_potential: Option<&Array4>, ) -> EosResult> { - let dft = &bulk.eos; + let dft: &F = &bulk.eos; // generate grid let x = Axis::new_cartesian(self.n_grid[0], self.system_size[0], None)?; @@ -252,7 +251,7 @@ impl PoreSpecification for Pore3D { let external_potential = external_potential.map_or_else( || { external_potential_3d( - &bulk.eos.functional, + dft, [&x, &y, &z], self.system_size, coordinates, @@ -268,7 +267,7 @@ impl PoreSpecification for Pore3D { // initialize convolver let grid = Grid::Periodical3(x, y, z); - let weight_functions = dft.functional.weight_functions(t); + let weight_functions = dft.weight_functions(t); let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1)); Ok(PoreProfile { @@ -453,7 +452,7 @@ impl Helium { fn new() -> DFT { let epsilon = arr1(&[EPSILON_HE]); let sigma = arr1(&[SIGMA_HE]); - DFT::new_homosegmented(Self { epsilon, sigma }, &Array1::ones(1)) + (Self { epsilon, sigma }).into() } } @@ -469,6 +468,10 @@ impl HelmholtzEnergyFunctional for Helium { fn compute_max_density(&self, _: &Array1) -> f64 { 1.0 } + + fn molecule_shape(&self) -> MoleculeShape { + MoleculeShape::Spherical(1) + } } impl FluidParameters for Helium { @@ -479,8 +482,4 @@ impl FluidParameters for Helium { fn sigma_ff(&self) -> &Array1 { &self.sigma } - - fn m(&self) -> Array1 { - arr1(&[1.0]) - } } diff --git a/src/functional.rs b/src/functional.rs index 8bdad5a..47114cb 100644 --- a/src/functional.rs +++ b/src/functional.rs @@ -12,50 +12,40 @@ use petgraph::graph::{Graph, UnGraph}; use petgraph::visit::EdgeRef; use petgraph::Directed; use quantity::{QuantityArray, QuantityArray1, QuantityScalar}; +use std::borrow::Cow; use std::fmt; -use std::ops::{AddAssign, MulAssign}; +use std::ops::{AddAssign, Deref, MulAssign}; use std::rc::Rc; /// Wrapper struct for the [HelmholtzEnergyFunctional] trait. +/// +/// Needed (for now) to generically implement the `EquationOfState` +/// trait for Helmholtz energy functionals. #[derive(Clone)] -pub struct DFT { - /// Helmholtz energy functional - pub functional: T, - /// map segment -> component - pub component_index: Array1, - /// chain lengths of individual components - pub m: Array1, - /// ideal chain contribution - pub ideal_chain_contribution: IdealChainContribution, -} +pub struct DFT(F); -impl DFT { - /// Create a new DFT struct for a homosegmented Helmholtz energy functional. - pub fn new_homosegmented(functional: T, m: &Array1) -> Self { - let component_index = Array1::from_shape_fn(m.len(), |i| i); - Self::new(functional, &component_index, m) +impl From for DFT { + fn from(functional: F) -> Self { + Self(functional) } +} - /// Create a new DFT struct for a heterosegmented Helmholtz energy functional. - pub fn new_heterosegmented(functional: T, component_index: &Array1) -> Self { - let m = Array1::ones(component_index.len()); - Self::new(functional, component_index, &m) +impl DFT { + pub fn into>(self) -> DFT { + DFT(self.0.into()) } +} - /// Create a new DFT struct for a general Helmholtz energy functional. - pub fn new(functional: T, component_index: &Array1, m: &Array1) -> Self { - Self { - functional, - component_index: component_index.clone(), - m: m.clone(), - ideal_chain_contribution: IdealChainContribution::new(component_index, m), - } +impl Deref for DFT { + type Target = F; + fn deref(&self) -> &F { + &self.0 } } impl, U: EosUnit> MolarWeight for DFT { fn molar_weight(&self) -> QuantityArray1 { - self.functional.molar_weight() + (self as &T).molar_weight() } } @@ -74,15 +64,15 @@ impl fmt::Display for DefaultIdealGasContribution { impl EquationOfState for DFT { fn components(&self) -> usize { - self.component_index[self.component_index.len() - 1] + 1 + self.component_index()[self.component_index().len() - 1] + 1 } fn subset(&self, component_list: &[usize]) -> Self { - self.functional.subset(component_list) + (self as &T).subset(component_list) } fn compute_max_density(&self, moles: &Array1) -> f64 { - self.functional.compute_max_density(moles) + (self as &T).compute_max_density(moles) } fn residual(&self) -> &[Box] { @@ -93,12 +83,11 @@ impl EquationOfState for DFT { where dyn HelmholtzEnergy: HelmholtzEnergyDual, { - self.functional - .contributions() + self.contributions() .iter() .map(|c| (c as &dyn HelmholtzEnergy).helmholtz_energy(state)) .sum::() - + self.ideal_chain_contribution.helmholtz_energy(state) + + self.ideal_chain_contribution().helmholtz_energy(state) } fn evaluate_residual_contributions>( @@ -109,7 +98,6 @@ impl EquationOfState for DFT { dyn HelmholtzEnergy: HelmholtzEnergyDual, { let mut res: Vec<(String, D)> = self - .functional .contributions() .iter() .map(|c| { @@ -120,23 +108,36 @@ impl EquationOfState for DFT { }) .collect(); res.push(( - self.ideal_chain_contribution.to_string(), - self.ideal_chain_contribution.helmholtz_energy(state), + self.ideal_chain_contribution().to_string(), + self.ideal_chain_contribution().helmholtz_energy(state), )); res } fn ideal_gas(&self) -> &dyn IdealGasContribution { - self.functional.ideal_gas() + (self as &T).ideal_gas() } } +/// Different representations for molecules within DFT. +pub enum MoleculeShape<'a> { + /// For spherical molecules, the number of components. + Spherical(usize), + /// For non-spherical molecules in a homosegmented approach, the chain length parameter $m$. + NonSpherical(&'a Array1), + /// For non-spherical molecules in a heterosegmented approach, the component index for every segment. + Heterosegmented(&'a Array1), +} + /// A general Helmholtz energy functional. pub trait HelmholtzEnergyFunctional: Sized { /// Return a slice of [FunctionalContribution]s. fn contributions(&self) -> &[Box]; - /// Return a [DFT] for the specified subset of components. + /// Return the shape of the molecules and the necessary specifications. + fn molecule_shape(&self) -> MoleculeShape; + + /// Return a functional for the specified subset of components. fn subset(&self, component_list: &[usize]) -> DFT; /// Return the maximum density in Angstrom^-3. @@ -170,6 +171,28 @@ pub trait HelmholtzEnergyFunctional: Sized { .map(|c| c.weight_functions(temperature)) .collect() } + + fn m(&self) -> Cow> { + match self.molecule_shape() { + MoleculeShape::Spherical(n) => Cow::Owned(Array1::ones(n)), + MoleculeShape::NonSpherical(m) => Cow::Borrowed(m), + MoleculeShape::Heterosegmented(component_index) => { + Cow::Owned(Array1::ones(component_index.len())) + } + } + } + + fn component_index(&self) -> Cow> { + match self.molecule_shape() { + MoleculeShape::Spherical(n) => Cow::Owned(Array1::from_shape_fn(n, |i| i)), + MoleculeShape::NonSpherical(m) => Cow::Owned(Array1::from_shape_fn(m.len(), |i| i)), + MoleculeShape::Heterosegmented(component_index) => Cow::Borrowed(component_index), + } + } + + fn ideal_chain_contribution(&self) -> IdealChainContribution { + IdealChainContribution::new(&self.component_index(), &self.m()) + } } impl DFT { @@ -191,11 +214,15 @@ impl DFT { let (mut f, dfdrho) = self.functional_derivative(t, &rho, convolver)?; // calculate the grand potential density - for ((rho, dfdrho), &m) in rho.outer_iter().zip(dfdrho.outer_iter()).zip(self.m.iter()) { + for ((rho, dfdrho), &m) in rho + .outer_iter() + .zip(dfdrho.outer_iter()) + .zip(self.m().iter()) + { f -= &((&dfdrho + m) * &rho); } - let bond_lengths = self.functional.bond_lengths(t); + let bond_lengths = self.bond_lengths(t); for segment in bond_lengths.node_indices() { let n = bond_lengths.neighbors(segment).count(); f += &(&rho.index_axis(Axis(0), segment.index()) * (0.5 * n as f64)); @@ -214,7 +241,7 @@ impl DFT { D::Larger: Dimension, { let n = self.components(); - let ig = self.functional.ideal_gas(); + let ig = self.ideal_gas(); let lambda = ig.de_broglie_wavelength(temperature, n); let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); for (i, rhoi) in density.outer_iter().enumerate() { @@ -233,7 +260,7 @@ impl DFT { D::Larger: Dimension, { let n = self.components(); - let ig = self.functional.ideal_gas(); + let ig = self.ideal_gas(); let lambda = ig.de_broglie_wavelength(temperature, n); let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); for (i, rhoi) in density.outer_iter().enumerate() { @@ -256,9 +283,9 @@ impl DFT { { let density_dual = density.mapv(N::from); let weighted_densities = convolver.weighted_densities(&density_dual); - let functional_contributions = self.functional.contributions(); + let functional_contributions = self.contributions(); let mut helmholtz_energy_density: Array = self - .ideal_chain_contribution + .ideal_chain_contribution() .calculate_helmholtz_energy_density(&density.mapv(N::from))?; for (c, wd) in functional_contributions.iter().zip(weighted_densities) { let nwd = wd.shape()[0]; @@ -319,11 +346,11 @@ impl DFT { let density_dual = density.mapv(Dual64::from); let temperature_dual = Dual64::from(temperature).derive(); let weighted_densities = convolver.weighted_densities(&density_dual); - let functional_contributions = self.functional.contributions(); + let functional_contributions = self.contributions(); let mut helmholtz_energy_density: Vec> = Vec::with_capacity(functional_contributions.len() + 1); helmholtz_energy_density.push( - self.ideal_chain_contribution + self.ideal_chain_contribution() .calculate_helmholtz_energy_density(&density.mapv(Dual64::from))?, ); @@ -389,7 +416,7 @@ impl DFT { D::Larger: Dimension, { let weighted_densities = convolver.weighted_densities(density); - let contributions = self.functional.contributions(); + let contributions = self.contributions(); let mut partial_derivatives = Vec::with_capacity(contributions.len()); let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0))); for (c, wd) in contributions.iter().zip(weighted_densities) { @@ -424,7 +451,7 @@ impl DFT { D::Larger: Dimension, { // calculate weight functions - let bond_lengths = self.functional.bond_lengths(temperature).into_edge_type(); + let bond_lengths = self.bond_lengths(temperature).into_edge_type(); let mut bond_weight_functions = bond_lengths.map( |_, _| (), |_, &l| WeightFunction::new_scaled(arr1(&[l]), WeightFunctionShape::Delta), diff --git a/src/fundamental_measure_theory.rs b/src/fundamental_measure_theory.rs index 7a8a1ab..a8279e1 100644 --- a/src/fundamental_measure_theory.rs +++ b/src/fundamental_measure_theory.rs @@ -1,6 +1,6 @@ //! Helmholtz energy functionals from fundamental measure theory. use crate::adsorption::FluidParameters; -use crate::functional::{HelmholtzEnergyFunctional, DFT}; +use crate::functional::{HelmholtzEnergyFunctional, MoleculeShape, DFT}; use crate::functional_contribution::*; use crate::solvation::PairPotential; use crate::weight_functions::{WeightFunction, WeightFunctionInfo, WeightFunctionShape}; @@ -293,14 +293,12 @@ impl FMTFunctional { }); let contributions: Vec> = vec![Box::new(FMTContribution::new(&properties, version))]; - DFT::new_homosegmented( - Self { - properties, - contributions, - version, - }, - &Array1::ones(sigma.len()), - ) + (Self { + properties, + contributions, + version, + }) + .into() } } @@ -320,6 +318,10 @@ impl HelmholtzEnergyFunctional for FMTFunctional { fn compute_max_density(&self, moles: &Array1) -> f64 { moles.sum() / (moles * &self.properties.sigma).sum() * 1.2 } + + fn molecule_shape(&self) -> MoleculeShape { + MoleculeShape::Spherical(self.properties.sigma.len()) + } } impl PairPotential for FMTFunctional { @@ -342,8 +344,4 @@ impl FluidParameters for FMTFunctional { fn sigma_ff(&self) -> &Array1 { &self.properties.sigma } - - fn m(&self) -> Array1 { - Array::ones(self.properties.sigma.len()) - } } diff --git a/src/interface/mod.rs b/src/interface/mod.rs index d918bb2..b535433 100644 --- a/src/interface/mod.rs +++ b/src/interface/mod.rs @@ -78,7 +78,7 @@ impl PlanarInterface { .vapor() .temperature .to_reduced(U::reference_temperature())?; - let weight_functions = dft.functional.weight_functions(t); + let weight_functions = dft.weight_functions(t); let convolver = ConvolverFFT::plan(&grid, &weight_functions, None); Ok(Self { @@ -98,7 +98,7 @@ impl PlanarInterface { let mut profile = Self::new(vle, n_grid, l_grid)?; // calculate segment indices - let indices = &profile.profile.dft.component_index; + let indices = &profile.profile.dft.component_index(); // calculate density profile let z0 = 0.5 * l_grid.to_reduced(U::reference_length())?; @@ -125,8 +125,8 @@ impl PlanarInterface { pub fn from_pdgt(vle: &PhaseEquilibrium, 2>, n_grid: usize) -> EosResult { let dft = &vle.vapor().eos; - if dft.component_index.len() != 1 { - panic!("Initialization from pDGT not possible for segment DFT"); + if dft.component_index().len() != 1 { + panic!("Initialization from pDGT not possible for segment DFT or mixtures"); } // calculate density profile from pDGT @@ -174,7 +174,7 @@ impl PlanarInterface { impl PlanarInterface { pub fn shift_equimolar_inplace(&mut self) { let s = self.profile.density.shape(); - let m = &self.profile.dft.m; + let m = &self.profile.dft.m(); let mut rho_l = 0.0 * U::reference_density(); let mut rho_v = 0.0 * U::reference_density(); let mut rho = Array::zeros(s[1]) * U::reference_density(); @@ -249,7 +249,7 @@ fn interp_symmetric( .unwrap() - 0.5 }); - let segments = vle_pdgt.vapor().eos.component_index.len(); + let segments = vle_pdgt.vapor().eos.component_index().len(); let mut reduced_density = interp( &z_pdgt.to_reduced(U::reference_length())?, &reduced_density, diff --git a/src/interface/surface_tension_diagram.rs b/src/interface/surface_tension_diagram.rs index 28f4e6e..0122432 100644 --- a/src/interface/surface_tension_diagram.rs +++ b/src/interface/surface_tension_diagram.rs @@ -34,7 +34,7 @@ impl SurfaceTensionDiagram { ) } else { // initialize with pDGT for single segments and tanh for mixtures and segment DFT - if vle.vapor().eos.component_index.len() == 1 { + if vle.vapor().eos.component_index().len() == 1 { PlanarInterface::from_pdgt(vle, n_grid) } else { PlanarInterface::from_tanh( diff --git a/src/lib.rs b/src/lib.rs index 3402b76..c6020cc 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -18,7 +18,7 @@ mod solver; mod weight_functions; pub use convolver::{Convolver, ConvolverFFT}; -pub use functional::{HelmholtzEnergyFunctional, DFT}; +pub use functional::{HelmholtzEnergyFunctional, MoleculeShape, DFT}; pub use functional_contribution::{FunctionalContribution, FunctionalContributionDual}; pub use geometry::{Axis, Geometry, Grid}; pub use profile::{DFTProfile, DFTSpecification, DFTSpecifications}; diff --git a/src/pdgt.rs b/src/pdgt.rs index 30bd262..2d92ae4 100644 --- a/src/pdgt.rs +++ b/src/pdgt.rs @@ -162,13 +162,13 @@ impl DFT { let mut delta_omega = Array::zeros(n_grid) * U::reference_pressure(); let mut influence_diagonal = Array::zeros(density.raw_dim()) * U::reference_influence_parameter(); - for contribution in self.functional.contributions() { + for contribution in self.contributions() { let (f, c) = contribution.influence_diagonal(vle.vapor().temperature, &density)?; delta_omega += &f; influence_diagonal += &c; } delta_omega += &self - .ideal_chain_contribution + .ideal_chain_contribution() .helmholtz_energy_density::<_, Ix1>(vle.vapor().temperature, &density)?; let t = vle diff --git a/src/profile.rs b/src/profile.rs index 8345407..e83ef86 100644 --- a/src/profile.rs +++ b/src/profile.rs @@ -98,7 +98,7 @@ impl DFTSpecification, _: &State>, ) -> EosResult> { - let m = &profile.dft.m; + let m: &Array1 = &profile.dft.m(); Ok(match self { Self::ChemicalPotential => chemical_potential.clone(), Self::Moles { moles } => (moles / z).mapv(f64::ln) * m, @@ -195,7 +195,7 @@ where // initialize external potential let external_potential = external_potential.unwrap_or_else(|| { - let mut n_grid = vec![dft.component_index.len()]; + let mut n_grid = vec![dft.component_index().len()]; grid.axes() .iter() .for_each(|&ax| n_grid.push(ax.grid.len())); @@ -213,7 +213,7 @@ where * (-&external_potential).mapv(f64::exp); let mut density = Array::zeros(external_potential.raw_dim()); let bulk_density = bulk.partial_density.to_reduced(U::reference_density())?; - for (s, &c) in dft.component_index.iter().enumerate() { + for (s, &c) in dft.component_index().iter().enumerate() { density.index_axis_mut(Axis_nd(0), s).assign( &(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]), ); @@ -287,7 +287,7 @@ where let mut d = rho.raw_dim(); d[0] = self.dft.components(); let mut density_comps = Array::zeros(d); - for (i, &j) in self.dft.component_index.iter().enumerate() { + for (i, &j) in self.dft.component_index().iter().enumerate() { density_comps .index_axis_mut(Axis_nd(0), j) .assign(&rho.index_axis(Axis_nd(0), i)); @@ -312,7 +312,6 @@ where .to_reduced(U::reference_temperature())?; let lambda_de_broglie = self .dft - .functional .ideal_gas() .de_broglie_wavelength(temperature, self.bulk.eos.components()); let mu_comp = self @@ -320,7 +319,7 @@ where .to_reduced(U::reference_molar_energy())? / temperature - lambda_de_broglie; - Ok(self.dft.component_index.mapv(|c| mu_comp[c])) + Ok(self.dft.component_index().mapv(|c| mu_comp[c])) } } @@ -399,11 +398,10 @@ where // Update bulk state let lambda_de_broglie = self .dft - .functional .ideal_gas() .de_broglie_wavelength(temperature, bulk.eos.components()); let mut mu_comp = Array::zeros(bulk.eos.components()); - for (s, &c) in self.dft.component_index.iter().enumerate() { + for (s, &c) in self.dft.component_index().iter().enumerate() { mu_comp[c] = chemical_potential[s]; } bulk.update_chemical_potential( @@ -424,7 +422,7 @@ where .bond_integrals(temperature, &dfdrho, &self.convolver); // Euler-Lagrange equation - let m = &self.dft.m; + let m = &self.dft.m(); res_rho .outer_iter_mut() .zip(dfdrho.outer_iter()) @@ -547,7 +545,7 @@ where pub fn entropy_density(&self, contributions: Contributions) -> EosResult> { // initialize convolver let t = self.temperature.to_reduced(U::reference_temperature())?; - let functional_contributions = self.dft.functional.contributions(); + let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions .iter() .map(|c| c.weight_functions(Dual64::from(t).derive())) @@ -569,7 +567,7 @@ where pub fn internal_energy(&self, contributions: Contributions) -> EosResult> { // initialize convolver let t = self.temperature.to_reduced(U::reference_temperature())?; - let functional_contributions = self.dft.functional.contributions(); + let functional_contributions = self.dft.contributions(); let weight_functions: Vec> = functional_contributions .iter() .map(|c| c.weight_functions(Dual64::from(t).derive())) diff --git a/src/solvation/mod.rs b/src/solvation/mod.rs index 3699edc..e739c8c 100644 --- a/src/solvation/mod.rs +++ b/src/solvation/mod.rs @@ -71,7 +71,7 @@ impl SolvationProfil cutoff_radius: Option>, potential_cutoff: Option, ) -> EosResult { - let dft = &bulk.eos; + let dft: &F = &bulk.eos; let system_size = system_size.unwrap_or([40.0 * U::reference_length(); 3]); @@ -104,7 +104,7 @@ impl SolvationProfil // calculate external potential let external_potential = external_potential_3d( - &dft.functional, + dft, [&x, &y, &z], coordinates, sigma_ss, @@ -116,7 +116,7 @@ impl SolvationProfil // initialize convolver let grid = Grid::Cartesian3(x, y, z); - let weight_functions = dft.functional.weight_functions(t); + let weight_functions = dft.weight_functions(t); let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1)); Ok(Self { diff --git a/src/solvation/pair_correlation.rs b/src/solvation/pair_correlation.rs index c6c39a6..2a895c7 100644 --- a/src/solvation/pair_correlation.rs +++ b/src/solvation/pair_correlation.rs @@ -46,7 +46,7 @@ impl PairCorrelation MAX_POTENTIAL { *x = MAX_POTENTIAL @@ -55,7 +55,7 @@ impl PairCorrelation Date: Mon, 11 Apr 2022 12:40:40 +0200 Subject: [PATCH 11/14] Add optional field cutoff_radius to FEA Potential (#25) --- CHANGELOG.md | 1 + src/adsorption/external_potential.rs | 7 +++++++ src/adsorption/fea_potential.rs | 3 ++- src/python/adsorption/external_potential.rs | 6 +++++- 4 files changed, 15 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a795c53..4d69d72 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `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. diff --git a/src/adsorption/external_potential.rs b/src/adsorption/external_potential.rs index d23a9ec..47a14d7 100644 --- a/src/adsorption/external_potential.rs +++ b/src/adsorption/external_potential.rs @@ -49,6 +49,7 @@ pub enum ExternalPotential { pore_center: [f64; 3], system_size: [QuantityScalar; 3], n_grid: [usize; 2], + cutoff_radius: Option, }, /// Custom potential @@ -168,6 +169,7 @@ impl ExternalPotential { pore_center, system_size, n_grid, + cutoff_radius, } => { // combining rules let epsilon_k_sf = @@ -185,6 +187,7 @@ impl ExternalPotential { n_grid, temperature, Geometry::Cartesian, + *cutoff_radius, ) } Self::Custom(_) => unreachable!(), @@ -314,6 +317,7 @@ impl ExternalPotential { pore_center, system_size, n_grid, + cutoff_radius, } => { // combining rules let epsilon_k_sf = @@ -331,6 +335,7 @@ impl ExternalPotential { n_grid, temperature, Geometry::Cylindrical, + *cutoff_radius, ) } Self::Custom(_) => unreachable!(), @@ -475,6 +480,7 @@ impl ExternalPotential { pore_center, system_size, n_grid, + cutoff_radius, } => { // combining rules let epsilon_k_sf = @@ -492,6 +498,7 @@ impl ExternalPotential { n_grid, temperature, Geometry::Spherical, + *cutoff_radius, ) } Self::Custom(_) => unreachable!(), diff --git a/src/adsorption/fea_potential.rs b/src/adsorption/fea_potential.rs index da81444..e1ef4a0 100644 --- a/src/adsorption/fea_potential.rs +++ b/src/adsorption/fea_potential.rs @@ -20,12 +20,13 @@ pub fn calculate_fea_potential( n_grid: &[usize; 2], temperature: f64, geometry: Geometry, + cutoff_radius: Option, ) -> Array1 { // allocate external potential let mut potential: Array1 = 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)| { diff --git a/src/python/adsorption/external_potential.rs b/src/python/adsorption/external_potential.rs index 6d595b7..9b23e91 100644 --- a/src/python/adsorption/external_potential.rs +++ b/src/python/adsorption/external_potential.rs @@ -185,13 +185,15 @@ impl PyExternalPotential { /// The size of the unit cell. /// n_grid : [int; 2] /// The number of grid points in each direction. + /// cutoff_radius : float, optional + /// The cutoff used in the calculation of fluid/wall interactions. /// Returns /// ------- /// ExternalPotential /// #[staticmethod] #[pyo3( - text_signature = "(coordinates, sigma_ss, epsilon_k_ss, pore_center, system_size, n_grid)" + text_signature = "(coordinates, sigma_ss, epsilon_k_ss, pore_center, system_size, n_grid, cutoff_radius=None)" )] pub fn FreeEnergyAveraged( coordinates: &PySIArray2, @@ -200,6 +202,7 @@ impl PyExternalPotential { pore_center: [f64; 3], system_size: [PySINumber; 3], n_grid: [usize; 2], + cutoff_radius: Option, ) -> Self { Self(ExternalPotential::FreeEnergyAveraged { coordinates: coordinates.clone().into(), @@ -212,6 +215,7 @@ impl PyExternalPotential { system_size[2].into(), ], n_grid, + cutoff_radius, }) } } From c68bc73b4a2ef386e47dbcc56f7b1d7c0d74c097 Mon Sep 17 00:00:00 2001 From: Rolf Stierle Date: Tue, 12 Apr 2022 15:36:18 +0200 Subject: [PATCH 12/14] Tangential pressure (#22) * Tangential pressure * fix python getters * remove partially duplicate grand_potential getter Co-authored-by: Philipp Rehner --- src/adsorption/mod.rs | 17 ++++++++++++++--- src/adsorption/pore.rs | 15 ++++++--------- src/functional.rs | 2 +- src/interface/mod.rs | 7 ++----- src/profile.rs | 9 +++++++++ src/python/profile.rs | 7 +++++++ src/solvation/mod.rs | 8 +------- src/solvation/pair_correlation.rs | 7 ++----- 8 files changed, 42 insertions(+), 30 deletions(-) diff --git a/src/adsorption/mod.rs b/src/adsorption/mod.rs index 4b1a2b7..b03f49a 100644 --- a/src/adsorption/mod.rs +++ b/src/adsorption/mod.rs @@ -4,7 +4,7 @@ use super::solver::DFTSolver; use feos_core::{ Contributions, EosError, EosResult, EosUnit, EquationOfState, SolverOptions, StateBuilder, }; -use ndarray::{arr1, Array1, Dimension, Ix1, Ix3}; +use ndarray::{arr1, Array1, Dimension, Ix1, Ix3, RemoveAxis}; use quantity::{QuantityArray1, QuantityArray2, QuantityScalar}; use std::rc::Rc; @@ -55,12 +55,17 @@ where }) } - fn equilibrium( + fn equilibrium< + D: Dimension + RemoveAxis + 'static, + F: HelmholtzEnergyFunctional + FluidParameters, + >( &self, equilibrium: &Adsorption, ) -> EosResult<(QuantityArray1, QuantityArray1)> where D::Larger: Dimension, + D::Smaller: Dimension, + ::Larger: Dimension, { let p_eq = equilibrium.pressure().get(0); match self { @@ -111,10 +116,16 @@ pub type Adsorption1D = Adsorption; /// Container structure for adsorption isotherms in 3D pores. pub type Adsorption3D = Adsorption; -impl Adsorption +impl< + U: EosUnit, + D: Dimension + RemoveAxis + 'static, + F: HelmholtzEnergyFunctional + FluidParameters, + > Adsorption where QuantityScalar: std::fmt::Display, D::Larger: Dimension, + D::Smaller: Dimension, + ::Larger: Dimension, { fn new>( functional: &Rc>, diff --git a/src/adsorption/pore.rs b/src/adsorption/pore.rs index cd1a849..e4b441c 100644 --- a/src/adsorption/pore.rs +++ b/src/adsorption/pore.rs @@ -8,7 +8,7 @@ use crate::solver::DFTSolver; use feos_core::{Contributions, EosResult, EosUnit, State, StateBuilder}; use ndarray::prelude::*; use ndarray::Axis as Axis_nd; -use ndarray::Zip; +use ndarray::{RemoveAxis, Zip}; use ndarray_stats::QuantileExt; use quantity::{QuantityArray, QuantityArray2, QuantityArray4, QuantityScalar}; use std::rc::Rc; @@ -133,22 +133,19 @@ impl Clone for PoreProfile { } } -impl PoreProfile +impl + PoreProfile where D::Larger: Dimension, + D::Smaller: Dimension, + ::Larger: Dimension, { pub fn solve_inplace(&mut self, solver: Option<&DFTSolver>, debug: bool) -> EosResult<()> { // Solve the profile self.profile.solve(solver, debug)?; // calculate grand potential density - let omega = self - .profile - .integrate(&self.profile.dft.grand_potential_density( - self.profile.temperature, - &self.profile.density, - &self.profile.convolver, - )?); + let omega = self.profile.grand_potential()?; self.grand_potential = Some(omega); // calculate interfacial tension diff --git a/src/functional.rs b/src/functional.rs index 47114cb..8fe1818 100644 --- a/src/functional.rs +++ b/src/functional.rs @@ -213,7 +213,7 @@ impl DFT { let rho = density.to_reduced(U::reference_density())?; let (mut f, dfdrho) = self.functional_derivative(t, &rho, convolver)?; - // calculate the grand potential density + // Calculate the grand potential density for ((rho, dfdrho), &m) in rho .outer_iter() .zip(dfdrho.outer_iter()) diff --git a/src/interface/mod.rs b/src/interface/mod.rs index b535433..0c98417 100644 --- a/src/interface/mod.rs +++ b/src/interface/mod.rs @@ -40,11 +40,8 @@ impl PlanarInterface { // postprocess self.surface_tension = Some(self.profile.integrate( - &(self.profile.dft.grand_potential_density( - self.profile.temperature, - &self.profile.density, - &self.profile.convolver, - )? + self.vle.vapor().pressure(Contributions::Total)), + &(self.profile.grand_potential_density()? + + self.vle.vapor().pressure(Contributions::Total)), )); let delta_rho = self.vle.liquid().density - self.vle.vapor().density; self.equimolar_radius = Some( diff --git a/src/profile.rs b/src/profile.rs index e83ef86..1d0b625 100644 --- a/src/profile.rs +++ b/src/profile.rs @@ -564,6 +564,15 @@ where Ok(self.integrate(&self.entropy_density(contributions)?)) } + pub fn grand_potential_density(&self) -> EosResult> { + self.dft + .grand_potential_density(self.temperature, &self.density, &self.convolver) + } + + pub fn grand_potential(&self) -> EosResult> { + Ok(self.integrate(&self.grand_potential_density()?)) + } + pub fn internal_energy(&self, contributions: Contributions) -> EosResult> { // initialize convolver let t = self.temperature.to_reduced(U::reference_temperature())?; diff --git a/src/python/profile.rs b/src/python/profile.rs index a73049b..5829ba1 100644 --- a/src/python/profile.rs +++ b/src/python/profile.rs @@ -171,6 +171,13 @@ macro_rules! impl_profile { self.0.profile.internal_energy(contributions)?, )) } + + #[getter] + fn get_grand_potential_density(&self) -> PyResult<$si_arr> { + Ok($si_arr::from( + self.0.profile.grand_potential_density()?, + )) + } } }; } diff --git a/src/solvation/mod.rs b/src/solvation/mod.rs index e739c8c..bb23fff 100644 --- a/src/solvation/mod.rs +++ b/src/solvation/mod.rs @@ -36,13 +36,7 @@ impl SolvationProfile { self.profile.solve(solver, debug)?; // calculate grand potential density - let omega = self - .profile - .integrate(&self.profile.dft.grand_potential_density( - self.profile.temperature, - &self.profile.density, - &self.profile.convolver, - )?); + let omega = self.profile.grand_potential()?; self.grand_potential = Some(omega); // calculate solvation free energy diff --git a/src/solvation/pair_correlation.rs b/src/solvation/pair_correlation.rs index 2a895c7..a9b5794 100644 --- a/src/solvation/pair_correlation.rs +++ b/src/solvation/pair_correlation.rs @@ -76,11 +76,8 @@ impl PairCorrelation Date: Tue, 12 Apr 2022 15:41:42 +0200 Subject: [PATCH 13/14] update changelog and version number --- CHANGELOG.md | 5 ++++- Cargo.toml | 2 +- build_wheel/Cargo.toml | 2 +- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4d69d72..9a601d1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,7 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] -## [0.2.0] - 2022-04-?? +## [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) diff --git a/Cargo.toml b/Cargo.toml index cb9708d..bf66355 100644 --- 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 "] edition = "2018" license = "MIT OR Apache-2.0" diff --git a/build_wheel/Cargo.toml b/build_wheel/Cargo.toml index 64d6cf3..a8e9eaf 100644 --- a/build_wheel/Cargo.toml +++ b/build_wheel/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "feos_dft" -version = "0.1.3" +version = "0.2.0" authors = ["Philipp Rehner "] edition = "2018" From 5e85d3c029790064a06db61eae8c92ce36e55de1 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Tue, 12 Apr 2022 15:51:59 +0200 Subject: [PATCH 14/14] udpate feos-core dependency --- Cargo.toml | 2 +- build_wheel/Cargo.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index bf66355..a974e19 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,7 +17,7 @@ rustdoc-args = [ "--html-in-header", "./docs-header.html" ] [dependencies] quantity = { version = "0.5", features = ["linalg"] } -feos-core = { git = "https://github.com/feos-org/feos-core" } +feos-core = "0.2" num-dual = "0.5" ndarray = { version = "0.15", features = ["serde", "rayon"] } ndarray-stats = "0.5" diff --git a/build_wheel/Cargo.toml b/build_wheel/Cargo.toml index a8e9eaf..762a796 100644 --- a/build_wheel/Cargo.toml +++ b/build_wheel/Cargo.toml @@ -9,7 +9,7 @@ crate-type = ["cdylib"] [dependencies] quantity = "0.5" -feos-core = { git = "https://github.com/feos-org/feos-core" } +feos-core = "0.2" feos-dft = { path = "..", features = ["python"] } pyo3 = { version = "0.16", features = ["extension-module", "abi3", "abi3-py37"] } numpy = "0.16"