From c291691754aa45442c5daa963f8df3e107264d6e Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Tue, 18 Nov 2025 14:07:46 +0100 Subject: [PATCH 1/3] Automatic stub file generation --- py-feos/Cargo.toml | 4 +++- py-feos/src/bin/stub_gen.rs | 18 ++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) create mode 100644 py-feos/src/bin/stub_gen.rs diff --git a/py-feos/Cargo.toml b/py-feos/Cargo.toml index e272ccdbd..c4418087e 100644 --- a/py-feos/Cargo.toml +++ b/py-feos/Cargo.toml @@ -18,7 +18,8 @@ pyo3 = { version = "0.27", features = [ "extension-module", "abi3-py310", "multiple-pymethods", - "indexmap" + "indexmap", + "experimental-inspect" ] } pythonize = "0.27" numpy = { version = "0.27" } @@ -43,6 +44,7 @@ feos = { workspace = true } feos-core = { workspace = true, features = ["ndarray"] } feos-derive = { workspace = true } feos-dft = { workspace = true, optional = true } +pyo3-introspection = "0.27.1" [features] default = [] diff --git a/py-feos/src/bin/stub_gen.rs b/py-feos/src/bin/stub_gen.rs new file mode 100644 index 000000000..5a5752b8c --- /dev/null +++ b/py-feos/src/bin/stub_gen.rs @@ -0,0 +1,18 @@ +use std::path::{Path, PathBuf}; + +fn main() { + println!("Generating stub files"); + let path = Path::new(".venv/lib/python3.12/site-packages/feos/feos.abi3.so"); + + assert!(path.is_file(), "Failed to locate cdylib"); + let main_module_name = "feos"; + let python_module = pyo3_introspection::introspect_cdylib(path, main_module_name) + .expect(format!("Failed introspection of {}", main_module_name).as_str()); + let type_stubs = pyo3_introspection::module_stub_files(&python_module); + + let stubst_string = type_stubs + .get(&PathBuf::from("__init__.pyi")) + .expect("Failed to get __init__.pyi"); + std::fs::write("feos.pyi", stubst_string).expect("Failed to write stubs file"); + println!("Generated stubs: {}", "feos.pyi") +} From 5384f7f0251f4abad35fa222ddbf4f9347b1f52a Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Wed, 19 Nov 2025 10:11:54 +0100 Subject: [PATCH 2/3] Change module creation --- py-feos/src/lib.rs | 113 ++++++++++++++++++--------------------------- 1 file changed, 45 insertions(+), 68 deletions(-) diff --git a/py-feos/src/lib.rs b/py-feos/src/lib.rs index b4f94441a..77b14604d 100644 --- a/py-feos/src/lib.rs +++ b/py-feos/src/lib.rs @@ -52,80 +52,57 @@ impl From for Verbosity { } #[pymodule] -fn feos(m: &Bound<'_, PyModule>) -> PyResult<()> { - m.add("__version__", env!("CARGO_PKG_VERSION"))?; +pub mod feos { + pub const __version__: &'static str = env!("CARGO_PKG_VERSION"); + + #[pymodule_export] + use super::PyVerbosity; + #[pymodule_export] + use crate::phase_equilibria::{PyPhaseDiagram, PyPhaseDiagramHetero, PyPhaseEquilibrium}; + #[pymodule_export] + use crate::state::{PyContributions, PyState, PyStateVec}; + + #[pymodule_export] + use crate::eos::PyEquationOfState; + #[pymodule_export] + use crate::parameter::{ + PyBinaryRecord, PyBinarySegmentRecord, PyChemicalRecord, PyGcParameters, PyIdentifier, + PyIdentifierOption, PyParameters, PyPureRecord, PySegmentRecord, PySmartsRecord, + }; - // Utility - m.add_class::()?; - - // State & phase equilibria. - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - - // Parameter - 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::()?; - - // Equation of state - m.add_class::()?; - - // // Estimator - // m.add_class::()?; - // m.add_class::()?; - // m.add_class::()?; - // m.add_class::()?; - - // AD #[cfg(feature = "ad")] - { - m.add_function(wrap_pyfunction!(ad::vapor_pressure_derivatives, m)?)?; - m.add_function(wrap_pyfunction!(ad::liquid_density_derivatives, m)?)?; - m.add_function(wrap_pyfunction!( - ad::equilibrium_liquid_density_derivatives, - m - )?)?; - m.add_function(wrap_pyfunction!(ad::bubble_point_pressure_derivatives, m)?)?; - m.add_function(wrap_pyfunction!(ad::dew_point_pressure_derivatives, m)?)?; - m.add_class::()?; - } + #[pymodule_export] + use crate::ad::{ + bubble_point_pressure_derivatives, dew_point_pressure_derivatives, + equilibrium_liquid_density_derivatives, liquid_density_derivatives, + vapor_pressure_derivatives, + }; #[cfg(feature = "dft")] - { - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - - // Solver - m.add_class::()?; - m.add_class::()?; - + #[pymodule_export] + use crate::dft::{ // Adsorption - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; - m.add_class::()?; + PyAdsorption1D, + PyAdsorption3D, + // Solver + PyDFTSolver, + PyDFTSolverLog, - // Interface - m.add_class::()?; - m.add_class::()?; + PyExternalPotential, + PyFMTVersion, + PyGeometry, + PyHelmholtzEnergyFunctional, // Solvation - m.add_class::()?; - m.add_class::()?; - } - Ok(()) + PyPairCorrelation, + PyPlanarInterface, + + PyPore1D, + PyPore2D, + PyPore3D, + + PySolvationProfile, + // Interface + PySurfaceTensionDiagram, + }; } From 38dd5944541bd53f775ae17ca2d92103a0db0b51 Mon Sep 17 00:00:00 2001 From: Gernot Bauer Date: Wed, 19 Nov 2025 14:32:30 +0100 Subject: [PATCH 3/3] Correct paths, stylesheets for catppuccin --- py-feos/docs/_static/custom_css.css | 427 ++++++++++++++++++ py-feos/docs/_static/latte.html | 282 ++++++++++++ py-feos/docs/_static/mocha.html | 282 ++++++++++++ py-feos/docs/_static/slate.scss | 150 ++++++ py-feos/docs/api/eos/eos.md | 1 + py-feos/docs/api/eos/epcsaft.md | 28 ++ py-feos/docs/api/eos/gc_pcsaft.md | 26 ++ py-feos/docs/api/eos/pcsaft.md | 17 + py-feos/docs/api/eos/pets.md | 25 + py-feos/docs/api/eos/saftvrmie.md | 17 + py-feos/docs/api/eos/saftvrqmie.md | 28 ++ py-feos/docs/api/eos/uvtheory.md | 28 ++ py-feos/docs/api/index.md | 35 ++ py-feos/docs/api/misc.md | 18 + py-feos/docs/api/parameters/identifier.md | 24 + py-feos/docs/api/parameters/parameters.md | 42 ++ py-feos/docs/api/parameters/records.md | 61 +++ py-feos/docs/api/state/state.md | 102 +++++ py-feos/docs/examples/basic-properties.md | 236 ++++++++++ py-feos/docs/examples/parameter-loading.md | 359 +++++++++++++++ py-feos/docs/getting-started/concepts.md | 265 +++++++++++ py-feos/docs/getting-started/installation.md | 151 +++++++ py-feos/docs/getting-started/quickstart.md | 185 ++++++++ py-feos/docs/hooks/typst_math.py | 131 ++++++ py-feos/docs/index.md | 102 +++++ py-feos/docs/javascript/mathjax.js | 20 + py-feos/docs/theory/dft/derivatives.md | 81 ++++ .../docs/theory/dft/enthalpy_of_adsorption.md | 113 +++++ .../theory/dft/euler_lagrange_equation.md | 121 +++++ .../docs/theory/dft/functional_derivatives.md | 44 ++ py-feos/docs/theory/dft/ideal_gas.md | 33 ++ py-feos/docs/theory/dft/pdgt.md | 52 +++ py-feos/docs/theory/dft/solver.md | 125 +++++ .../docs/user-guide/working-with-states.md | 328 ++++++++++++++ py-feos/mkdocs.yml | 192 ++++++++ py-feos/src/ad/mod.rs | 2 +- .../src/dft/adsorption/external_potential.rs | 4 +- py-feos/src/dft/adsorption/mod.rs | 8 +- py-feos/src/dft/adsorption/pore.rs | 12 +- py-feos/src/dft/interface/mod.rs | 2 +- .../dft/interface/surface_tension_diagram.rs | 2 +- py-feos/src/dft/mod.rs | 8 +- py-feos/src/dft/solvation.rs | 4 +- py-feos/src/dft/solver.rs | 4 +- py-feos/src/eos/mod.rs | 2 +- py-feos/src/lib.rs | 9 +- py-feos/src/parameter/chemical_record.rs | 4 +- py-feos/src/parameter/fragmentation.rs | 2 +- py-feos/src/parameter/identifier.rs | 4 +- py-feos/src/parameter/mod.rs | 4 +- py-feos/src/parameter/model_record.rs | 4 +- py-feos/src/parameter/segment.rs | 4 +- py-feos/src/phase_equilibria.rs | 10 +- py-feos/src/state.rs | 12 +- 54 files changed, 4184 insertions(+), 48 deletions(-) create mode 100644 py-feos/docs/_static/custom_css.css create mode 100644 py-feos/docs/_static/latte.html create mode 100644 py-feos/docs/_static/mocha.html create mode 100644 py-feos/docs/_static/slate.scss create mode 100644 py-feos/docs/api/eos/eos.md create mode 100644 py-feos/docs/api/eos/epcsaft.md create mode 100644 py-feos/docs/api/eos/gc_pcsaft.md create mode 100644 py-feos/docs/api/eos/pcsaft.md create mode 100644 py-feos/docs/api/eos/pets.md create mode 100644 py-feos/docs/api/eos/saftvrmie.md create mode 100644 py-feos/docs/api/eos/saftvrqmie.md create mode 100644 py-feos/docs/api/eos/uvtheory.md create mode 100644 py-feos/docs/api/index.md create mode 100644 py-feos/docs/api/misc.md create mode 100644 py-feos/docs/api/parameters/identifier.md create mode 100644 py-feos/docs/api/parameters/parameters.md create mode 100644 py-feos/docs/api/parameters/records.md create mode 100644 py-feos/docs/api/state/state.md create mode 100644 py-feos/docs/examples/basic-properties.md create mode 100644 py-feos/docs/examples/parameter-loading.md create mode 100644 py-feos/docs/getting-started/concepts.md create mode 100644 py-feos/docs/getting-started/installation.md create mode 100644 py-feos/docs/getting-started/quickstart.md create mode 100644 py-feos/docs/hooks/typst_math.py create mode 100644 py-feos/docs/index.md create mode 100644 py-feos/docs/javascript/mathjax.js create mode 100644 py-feos/docs/theory/dft/derivatives.md create mode 100644 py-feos/docs/theory/dft/enthalpy_of_adsorption.md create mode 100644 py-feos/docs/theory/dft/euler_lagrange_equation.md create mode 100644 py-feos/docs/theory/dft/functional_derivatives.md create mode 100644 py-feos/docs/theory/dft/ideal_gas.md create mode 100644 py-feos/docs/theory/dft/pdgt.md create mode 100644 py-feos/docs/theory/dft/solver.md create mode 100644 py-feos/docs/user-guide/working-with-states.md create mode 100644 py-feos/mkdocs.yml diff --git a/py-feos/docs/_static/custom_css.css b/py-feos/docs/_static/custom_css.css new file mode 100644 index 000000000..e05f89828 --- /dev/null +++ b/py-feos/docs/_static/custom_css.css @@ -0,0 +1,427 @@ +/* Fix /page#foo going to the top of the viewport and being hidden by the navbar */ +html { + scroll-padding-top: 50px; +} + +/* Fit the Twitter handle alongside the GitHub one in the top right. */ +div.md-header__source { + width: revert; + max-width: revert; +} + +a.md-source { + display: inline-block; +} + +.md-source__repository { + max-width: 100%; +} + +/* Emphasise sections of nav on left hand side */ + +nav.md-nav { + padding-left: 5px; +} + +nav.md-nav--secondary { + border-left: revert !important; +} + +.md-nav__title { + font-size: 0.9rem; +} + +.md-nav__item--section > .md-nav__link { + font-size: 0.9rem; +} + +/* Indent autogenerated documentation */ +div.doc-contents { + padding-left: 2.5em; + /* border-left: 0.5px solid; */ +} + +/* Increase visibility of splitters "---" */ + +/* [data-md-color-scheme="default"] .md-typeset hr { + border-bottom-color: rgb(0, 0, 0); + border-bottom-width: 1pt; +} */ + +/* [data-md-color-scheme="slate"] .md-typeset hr { + border-bottom-color: rgb(230, 230, 230); +} */ + +/* More space at the bottom of the page */ + +.md-main__inner { + margin-bottom: 1.5rem; +} + +/* Remove prev/next footer buttons */ + +.md-footer__inner { + display: none; +} + +/* Change font sizes */ + +html { + /* Decrease font size for overall webpage + Down from 137.5% which is the Material default */ + font-size: 120%; +} + +.md-typeset .admonition { + /* Increase font size in admonitions */ + font-size: 100% !important; +} + +.md-typeset details { + /* Increase font size in details */ + font-size: 100% !important; +} + +.md-typeset h1 { + font-size: 1.6rem; +} + +.md-typeset h2 { + font-size: 1.5rem; +} + +.md-typeset h3 { + font-size: 1.3rem; +} + +.md-typeset h4 { + font-size: 1.1rem; +} + +.md-typeset h5 { + font-size: 0.9rem; +} + +.md-typeset h6 { + font-size: 0.8rem; +} + +/* Bugfix: remove the superfluous parts generated when doing: + +??? Blah + + ::: library.something +*/ + +.md-typeset details .mkdocstrings > h4 { + display: none; +} + +.md-typeset details .mkdocstrings > h5 { + display: none; +} + +/* ============================================================================ + Catppuccin Mocha Theme - "feos" + Based on Catppuccin Mocha color palette + ============================================================================ */ + +@media screen { + [data-md-color-scheme="feos"] { + /* Indicate dark color scheme */ + color-scheme: dark; + + /* Catppuccin Mocha Color Palette */ + --ctp-rosewater: #f5e0dc; + --ctp-flamingo: #f2cdcd; + --ctp-pink: #f5c2e7; + --ctp-mauve: #cba6f7; + --ctp-red: #f38ba8; + --ctp-maroon: #eba0ac; + --ctp-peach: #fab387; + --ctp-yellow: #f9e2af; + --ctp-green: #a6e3a1; + --ctp-teal: #94e2d5; + --ctp-sky: #89dceb; + --ctp-sapphire: #74c7ec; + --ctp-blue: #89b4fa; + --ctp-lavender: #b4befe; + --ctp-text: #cdd6f4; + --ctp-subtext1: #bac2de; + --ctp-subtext0: #a6adc8; + --ctp-overlay2: #9399b2; + --ctp-overlay1: #7f849c; + --ctp-overlay0: #6c7086; + --ctp-surface2: #585b70; + --ctp-surface1: #45475a; + --ctp-surface0: #313244; + --ctp-base: #1e1e2e; + --ctp-mantle: #181825; + --ctp-crust: #11111b; + + /* Default color shades */ + --md-default-fg-color: var(--ctp-text); + --md-default-fg-color--light: var(--ctp-subtext1); + --md-default-fg-color--lighter: var(--ctp-subtext0); + --md-default-fg-color--lightest: var(--ctp-overlay2); + --md-default-bg-color: var(--ctp-base); + --md-default-bg-color--light: var(--ctp-surface0); + --md-default-bg-color--lighter: var(--ctp-surface1); + --md-default-bg-color--lightest: var(--ctp-surface2); + + /* Code color shades */ + --md-code-fg-color: var(--ctp-text); + --md-code-bg-color: var(--ctp-surface0); + --md-code-bg-color--light: var(--ctp-surface1); + --md-code-bg-color--lighter: var(--ctp-surface2); + + /* Code highlighting color shades */ + --md-code-hl-color: var(--ctp-blue); + --md-code-hl-color--light: rgba(137, 180, 250, 0.15); + + /* Code highlighting syntax colors */ + --md-code-hl-number-color: var(--ctp-peach); + --md-code-hl-special-color: var(--ctp-pink); + --md-code-hl-function-color: var(--ctp-blue); + --md-code-hl-constant-color: var(--ctp-peach); + --md-code-hl-keyword-color: var(--ctp-mauve); + --md-code-hl-string-color: var(--ctp-green); + --md-code-hl-name-color: var(--ctp-text); + --md-code-hl-operator-color: var(--ctp-sky); + --md-code-hl-punctuation-color: var(--ctp-overlay2); + --md-code-hl-comment-color: var(--ctp-overlay0); + --md-code-hl-generic-color: var(--ctp-subtext0); + --md-code-hl-variable-color: var(--ctp-text); + + /* Typeset color shades */ + --md-typeset-color: var(--ctp-text); + + /* Typeset link colors */ + --md-typeset-a-color: var(--ctp-sapphire); + + /* Typeset kbd color shades */ + --md-typeset-kbd-color: var(--ctp-surface1); + --md-typeset-kbd-accent-color: var(--ctp-surface2); + --md-typeset-kbd-border-color: var(--ctp-overlay0); + + /* Typeset mark color shades */ + --md-typeset-mark-color: rgba(249, 226, 175, 0.3); + + /* Typeset table color shades */ + --md-typeset-table-color: var(--ctp-surface0); + --md-typeset-table-color--light: var(--ctp-surface1); + + /* Admonition color shades */ + --md-admonition-fg-color: var(--ctp-text); + --md-admonition-bg-color: var(--ctp-base); + + /* Footer color shades */ + --md-footer-bg-color: var(--ctp-mantle); + --md-footer-bg-color--dark: var(--ctp-crust); + + /* Shadow depths */ + --md-shadow-z1: + 0 0.25rem 0.625rem rgba(0, 0, 0, 0.2), + 0 0 0.0625rem rgba(0, 0, 0, 0.15); + + --md-shadow-z2: + 0 0.25rem 0.625rem rgba(0, 0, 0, 0.3), + 0 0 0.0625rem rgba(0, 0, 0, 0.3); + + --md-shadow-z3: + 0 0.25rem 0.625rem rgba(0, 0, 0, 0.4), + 0 0 0.0625rem rgba(0, 0, 0, 0.4); + + /* Custom documentation heading colors */ + --doc-heading-color: var(--ctp-text); + --doc-heading-border-color: var(--ctp-surface1); + --doc-heading-color-alt: var(--ctp-surface0); + + /* Hide light mode images */ + img[src$="#only-light"], + img[src$="#gh-light-mode-only"] { + display: none; + } + + /* Primary color accents - adjust link colors */ + --md-primary-fg-color: var(--ctp-blue); + --md-primary-fg-color--light: var(--ctp-sapphire); + --md-primary-fg-color--dark: var(--ctp-lavender); + --md-primary-bg-color: var(--ctp-text); + --md-primary-bg-color--light: var(--ctp-subtext1); + + /* Accent colors */ + --md-accent-fg-color: var(--ctp-mauve); + --md-accent-fg-color--transparent: rgba(203, 166, 247, 0.1); + --md-accent-bg-color: var(--ctp-text); + --md-accent-bg-color--light: var(--ctp-subtext1); + } + + /* ======================================================================== + Catppuccin Latte Theme - "feos-light" + Based on Catppuccin Latte color palette (light mode) + ======================================================================== */ + + [data-md-color-scheme="feos-light"] { + /* Indicate light color scheme */ + color-scheme: light; + + /* Catppuccin Latte Color Palette */ + --ctp-rosewater: #dc8a78; + --ctp-flamingo: #dd7878; + --ctp-pink: #ea76cb; + --ctp-mauve: #8839ef; + --ctp-red: #d20f39; + --ctp-maroon: #e64553; + --ctp-peach: #fe640b; + --ctp-yellow: #df8e1d; + --ctp-green: #40a02b; + --ctp-teal: #179299; + --ctp-sky: #04a5e5; + --ctp-sapphire: #209fb5; + --ctp-blue: #1e66f5; + --ctp-lavender: #7287fd; + --ctp-text: #4c4f69; + --ctp-subtext1: #5c5f77; + --ctp-subtext0: #6c6f85; + --ctp-overlay2: #7c7f93; + --ctp-overlay1: #8c8fa1; + --ctp-overlay0: #9ca0b0; + --ctp-surface2: #acb0be; + --ctp-surface1: #bcc0cc; + --ctp-surface0: #ccd0da; + --ctp-base: #eff1f5; + --ctp-mantle: #e6e9ef; + --ctp-crust: #dce0e8; + + /* Default color shades */ + --md-default-fg-color: var(--ctp-text); + --md-default-fg-color--light: var(--ctp-subtext1); + --md-default-fg-color--lighter: var(--ctp-subtext0); + --md-default-fg-color--lightest: var(--ctp-overlay2); + --md-default-bg-color: var(--ctp-base); + --md-default-bg-color--light: var(--ctp-mantle); + --md-default-bg-color--lighter: var(--ctp-crust); + --md-default-bg-color--lightest: var(--ctp-surface0); + + /* Code color shades */ + --md-code-fg-color: var(--ctp-text); + --md-code-bg-color: var(--ctp-mantle); + --md-code-bg-color--light: var(--ctp-crust); + --md-code-bg-color--lighter: var(--ctp-surface0); + + /* Code highlighting color shades */ + --md-code-hl-color: var(--ctp-blue); + --md-code-hl-color--light: rgba(30, 102, 245, 0.15); + + /* Code highlighting syntax colors */ + --md-code-hl-number-color: var(--ctp-peach); + --md-code-hl-special-color: var(--ctp-pink); + --md-code-hl-function-color: var(--ctp-blue); + --md-code-hl-constant-color: var(--ctp-peach); + --md-code-hl-keyword-color: var(--ctp-mauve); + --md-code-hl-string-color: var(--ctp-green); + --md-code-hl-name-color: var(--ctp-text); + --md-code-hl-operator-color: var(--ctp-sky); + --md-code-hl-punctuation-color: var(--ctp-overlay1); + --md-code-hl-comment-color: var(--ctp-overlay0); + --md-code-hl-generic-color: var(--ctp-subtext0); + --md-code-hl-variable-color: var(--ctp-text); + + /* Typeset color shades */ + --md-typeset-color: var(--ctp-text); + + /* Typeset link colors */ + --md-typeset-a-color: var(--ctp-sapphire); + + /* Typeset kbd color shades */ + --md-typeset-kbd-color: var(--ctp-surface0); + --md-typeset-kbd-accent-color: var(--ctp-surface1); + --md-typeset-kbd-border-color: var(--ctp-overlay1); + + /* Typeset mark color shades */ + --md-typeset-mark-color: rgba(223, 142, 29, 0.3); + + /* Typeset table color shades */ + --md-typeset-table-color: var(--ctp-mantle); + --md-typeset-table-color--light: var(--ctp-crust); + + /* Admonition color shades */ + --md-admonition-fg-color: var(--ctp-text); + --md-admonition-bg-color: var(--ctp-base); + + /* Footer color shades */ + --md-footer-bg-color: var(--ctp-mantle); + --md-footer-bg-color--dark: var(--ctp-crust); + + /* Shadow depths */ + --md-shadow-z1: + 0 0.25rem 0.625rem rgba(76, 79, 105, 0.05), + 0 0 0.0625rem rgba(76, 79, 105, 0.1); + + --md-shadow-z2: + 0 0.25rem 0.625rem rgba(76, 79, 105, 0.1), + 0 0 0.0625rem rgba(76, 79, 105, 0.15); + + --md-shadow-z3: + 0 0.25rem 0.625rem rgba(76, 79, 105, 0.15), + 0 0 0.0625rem rgba(76, 79, 105, 0.2); + + /* Custom documentation heading colors */ + --doc-heading-color: var(--ctp-text); + --doc-heading-border-color: var(--ctp-surface1); + --doc-heading-color-alt: var(--ctp-surface0); + + /* Hide dark mode images */ + img[src$="#only-dark"], + img[src$="#gh-dark-mode-only"] { + display: none; + } + + /* Primary color accents */ + --md-primary-fg-color: var(--ctp-blue); + --md-primary-fg-color--light: var(--ctp-lavender); + --md-primary-fg-color--dark: var(--ctp-sapphire); + --md-primary-bg-color: var(--ctp-text); + --md-primary-bg-color--light: var(--ctp-subtext1); + + /* Accent colors */ + --md-accent-fg-color: var(--ctp-mauve); + --md-accent-fg-color--transparent: rgba(136, 57, 239, 0.1); + --md-accent-bg-color: var(--ctp-text); + --md-accent-bg-color--light: var(--ctp-subtext1); + } +} +/* Highlight functions, classes etc. type signatures. Really helps to make clear where + one item ends and another begins. */ + +/* [data-md-color-scheme="default"] { + --doc-heading-color: #DDD; + --doc-heading-border-color: #CCC; + --doc-heading-color-alt: #F0F0F0; +} */ +/* [data-md-color-scheme="slate"] { + --doc-heading-color: rgb(25,25,33); + --doc-heading-border-color: rgb(25,25,33); + --doc-heading-color-alt: rgb(33,33,44); + --md-code-bg-color: rgb(38,38,50); +} */ + +/* h4.doc-heading { + background-color: var(--doc-heading-color); + border: solid var(--doc-heading-border-color); + border-width: 1.5pt; + border-radius: 2pt; + padding: 0pt 5pt 2pt 5pt; +} + +h5.doc-heading, h6.heading { + background-color: var(--doc-heading-color-alt); + border-radius: 2pt; + padding: 0pt 5pt 2pt 5pt; +} */ + +/* Make errors in notebooks have scrolling */ +.output_error > pre { + overflow: auto; +} diff --git a/py-feos/docs/_static/latte.html b/py-feos/docs/_static/latte.html new file mode 100644 index 000000000..90628799b --- /dev/null +++ b/py-feos/docs/_static/latte.html @@ -0,0 +1,282 @@ +
+ 🌻 Latte + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
LabelsHexRGBHSL
+ + Rosewater#dc8a78rgb(220, 138, 120)hsl(11, 59%, 67%)
+ + Flamingo#dd7878rgb(221, 120, 120)hsl(0, 60%, 67%)
+ + Pink#ea76cbrgb(234, 118, 203)hsl(316, 73%, 69%)
+ + Mauve#8839efrgb(136, 57, 239)hsl(266, 85%, 58%)
+ + Red#d20f39rgb(210, 15, 57)hsl(347, 87%, 44%)
+ + Maroon#e64553rgb(230, 69, 83)hsl(355, 76%, 59%)
+ + Peach#fe640brgb(254, 100, 11)hsl(22, 99%, 52%)
+ + Yellow#df8e1drgb(223, 142, 29)hsl(35, 77%, 49%)
+ + Green#40a02brgb(64, 160, 43)hsl(109, 58%, 40%)
+ + Teal#179299rgb(23, 146, 153)hsl(183, 74%, 35%)
+ + Sky#04a5e5rgb(4, 165, 229)hsl(197, 97%, 46%)
+ + Sapphire#209fb5rgb(32, 159, 181)hsl(189, 70%, 42%)
+ + Blue#1e66f5rgb(30, 102, 245)hsl(220, 91%, 54%)
+ + Lavender#7287fdrgb(114, 135, 253)hsl(231, 97%, 72%)
+ + Text#4c4f69rgb(76, 79, 105)hsl(234, 16%, 35%)
+ + Subtext1#5c5f77rgb(92, 95, 119)hsl(233, 13%, 41%)
+ + Subtext0#6c6f85rgb(108, 111, 133)hsl(233, 10%, 47%)
+ + Overlay2#7c7f93rgb(124, 127, 147)hsl(232, 10%, 53%)
+ + Overlay1#8c8fa1rgb(140, 143, 161)hsl(231, 10%, 59%)
+ + Overlay0#9ca0b0rgb(156, 160, 176)hsl(228, 11%, 65%)
+ + Surface2#acb0bergb(172, 176, 190)hsl(227, 12%, 71%)
+ + Surface1#bcc0ccrgb(188, 192, 204)hsl(225, 14%, 77%)
+ + Surface0#ccd0dargb(204, 208, 218)hsl(223, 16%, 83%)
+ + Base#eff1f5rgb(239, 241, 245)hsl(220, 23%, 95%)
+ + Mantle#e6e9efrgb(230, 233, 239)hsl(220, 22%, 92%)
+ + Crust#dce0e8rgb(220, 224, 232)hsl(220, 21%, 89%)
+
diff --git a/py-feos/docs/_static/mocha.html b/py-feos/docs/_static/mocha.html new file mode 100644 index 000000000..ce8fe1b37 --- /dev/null +++ b/py-feos/docs/_static/mocha.html @@ -0,0 +1,282 @@ +
+ 🌿 Mocha + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
LabelsHexRGBHSL
+ + Rosewater#f5e0dcrgb(245, 224, 220)hsl(10, 56%, 91%)
+ + Flamingo#f2cdcdrgb(242, 205, 205)hsl(0, 59%, 88%)
+ + Pink#f5c2e7rgb(245, 194, 231)hsl(316, 72%, 86%)
+ + Mauve#cba6f7rgb(203, 166, 247)hsl(267, 84%, 81%)
+ + Red#f38ba8rgb(243, 139, 168)hsl(343, 81%, 75%)
+ + Maroon#eba0acrgb(235, 160, 172)hsl(350, 65%, 77%)
+ + Peach#fab387rgb(250, 179, 135)hsl(23, 92%, 75%)
+ + Yellow#f9e2afrgb(249, 226, 175)hsl(41, 86%, 83%)
+ + Green#a6e3a1rgb(166, 227, 161)hsl(115, 54%, 76%)
+ + Teal#94e2d5rgb(148, 226, 213)hsl(170, 57%, 73%)
+ + Sky#89dcebrgb(137, 220, 235)hsl(189, 71%, 73%)
+ + Sapphire#74c7ecrgb(116, 199, 236)hsl(199, 76%, 69%)
+ + Blue#89b4fargb(137, 180, 250)hsl(217, 92%, 76%)
+ + Lavender#b4befergb(180, 190, 254)hsl(232, 97%, 85%)
+ + Text#cdd6f4rgb(205, 214, 244)hsl(226, 64%, 88%)
+ + Subtext1#bac2dergb(186, 194, 222)hsl(227, 35%, 80%)
+ + Subtext0#a6adc8rgb(166, 173, 200)hsl(228, 24%, 72%)
+ + Overlay2#9399b2rgb(147, 153, 178)hsl(228, 17%, 64%)
+ + Overlay1#7f849crgb(127, 132, 156)hsl(230, 13%, 55%)
+ + Overlay0#6c7086rgb(108, 112, 134)hsl(231, 11%, 47%)
+ + Surface2#585b70rgb(88, 91, 112)hsl(233, 12%, 39%)
+ + Surface1#45475argb(69, 71, 90)hsl(234, 13%, 31%)
+ + Surface0#313244rgb(49, 50, 68)hsl(237, 16%, 23%)
+ + Base#1e1e2ergb(30, 30, 46)hsl(240, 21%, 15%)
+ + Mantle#181825rgb(24, 24, 37)hsl(240, 21%, 12%)
+ + Crust#11111brgb(17, 17, 27)hsl(240, 23%, 9%)
+
diff --git a/py-feos/docs/_static/slate.scss b/py-feos/docs/_static/slate.scss new file mode 100644 index 000000000..b170abf6f --- /dev/null +++ b/py-feos/docs/_static/slate.scss @@ -0,0 +1,150 @@ +//// +/// Copyright (c) 2016-2025 Martin Donath +/// +/// Permission is hereby granted, free of charge, to any person obtaining a +/// copy of this software and associated documentation files (the "Software"), +/// to deal in the Software without restriction, including without limitation +/// the rights to use, copy, modify, merge, publish, distribute, sublicense, +/// and/or sell copies of the Software, and to permit persons to whom the +/// Software is furnished to do so, subject to the following conditions: +/// +/// The above copyright notice and this permission notice shall be included in +/// all copies or substantial portions of the Software. +/// +/// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +/// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +/// FITNESS FOR A PARTICULAR PURPOSE AND NON-INFRINGEMENT. IN NO EVENT SHALL +/// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +/// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +/// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +/// DEALINGS +//// + +// ---------------------------------------------------------------------------- +// Rules +// ---------------------------------------------------------------------------- + +// Only use dark mode on screens +@media screen { + // Slate theme, i.e. dark mode + [data-md-color-scheme="slate"] { + // Indicate that the site is rendered with a dark color scheme + color-scheme: dark; + + // Default color shades + --md-default-fg-color: hsla(var(--md-hue), 15%, 90%, 0.82); + --md-default-fg-color--light: hsla(var(--md-hue), 15%, 90%, 0.56); + --md-default-fg-color--lighter: hsla(var(--md-hue), 15%, 90%, 0.32); + --md-default-fg-color--lightest: hsla(var(--md-hue), 15%, 90%, 0.12); + --md-default-bg-color: hsla(var(--md-hue), 15%, 14%, 1); + --md-default-bg-color--light: hsla(var(--md-hue), 15%, 14%, 0.54); + --md-default-bg-color--lighter: hsla(var(--md-hue), 15%, 14%, 0.26); + --md-default-bg-color--lightest: hsla(var(--md-hue), 15%, 14%, 0.07); + + // Code color shades + --md-code-fg-color: hsla(var(--md-hue), 18%, 86%, 0.82); + --md-code-bg-color: hsla(var(--md-hue), 15%, 18%, 1); + --md-code-bg-color--light: hsla(var(--md-hue), 15%, 18%, 0.9); + --md-code-bg-color--lighter: hsla(var(--md-hue), 15%, 18%, 0.54); + + // Code highlighting color shades + --md-code-hl-color: hsla(#{hex2hsl($clr-blue-a400)}, 1); + --md-code-hl-color--light: hsla(#{hex2hsl($clr-blue-a400)}, 0.1); + + // Code highlighting syntax color shades + --md-code-hl-number-color: hsla(6, 74%, 63%, 1); + --md-code-hl-special-color: hsla(340, 83%, 66%, 1); + --md-code-hl-function-color: hsla(291, 57%, 65%, 1); + --md-code-hl-constant-color: hsla(250, 62%, 70%, 1); + --md-code-hl-keyword-color: hsla(219, 66%, 64%, 1); + --md-code-hl-string-color: hsla(150, 58%, 44%, 1); + --md-code-hl-name-color: var(--md-code-fg-color); + --md-code-hl-operator-color: var(--md-default-fg-color--light); + --md-code-hl-punctuation-color: var(--md-default-fg-color--light); + --md-code-hl-comment-color: var(--md-default-fg-color--light); + --md-code-hl-generic-color: var(--md-default-fg-color--light); + --md-code-hl-variable-color: var(--md-default-fg-color--light); + + // Typeset color shades + --md-typeset-color: var(--md-default-fg-color); + + // Typeset `a` color shades + --md-typeset-a-color: var(--md-primary-fg-color); + + // Typeset `kbd` color shades + --md-typeset-kbd-color: hsla(var(--md-hue), 15%, 90%, 0.12); + --md-typeset-kbd-accent-color: hsla(var(--md-hue), 15%, 90%, 0.2); + --md-typeset-kbd-border-color: hsla(var(--md-hue), 15%, 14%, 1); + + // Typeset `mark` color shades + --md-typeset-mark-color: hsla(#{hex2hsl($clr-blue-a200)}, 0.3); + + // Typeset `table` color shades + --md-typeset-table-color: hsla(var(--md-hue), 15%, 95%, 0.12); + --md-typeset-table-color--light: hsla(var(--md-hue), 15%, 95%, 0.035); + + // Admonition color shades + --md-admonition-fg-color: var(--md-default-fg-color); + --md-admonition-bg-color: var(--md-default-bg-color); + + // Footer color shades + --md-footer-bg-color: hsla(var(--md-hue), 15%, 10%, 0.87); + --md-footer-bg-color--dark: hsla(var(--md-hue), 15%, 8%, 1); + + // Shadow depth 1 + --md-shadow-z1: + 0 #{px2rem(4px)} #{px2rem(10px)} hsla(0, 0%, 0%, 0.05), + 0 0 #{px2rem(1px)} hsla(0, 0%, 0%, 0.1); + + // Shadow depth 2 + --md-shadow-z2: + 0 #{px2rem(4px)} #{px2rem(10px)} hsla(0, 0%, 0%, 0.25), + 0 0 #{px2rem(1px)} hsla(0, 0%, 0%, 0.25); + + // Shadow depth 3 + --md-shadow-z3: + 0 #{px2rem(4px)} #{px2rem(10px)} hsla(0, 0%, 0%, 0.4), + 0 0 #{px2rem(1px)} hsla(0, 0%, 0%, 0.35); + + // Hide images for light mode + img[src$="#only-light"], + img[src$="#gh-light-mode-only"] { + display: none; + } + } + + // -------------------------------------------------------------------------- + + // Adjust link colors for dark primary colors + @each $name, + $color + in ( + "pink": hsl(340, 81%, 63%), + "purple": hsl(291, 53%, 63%), + "deep-purple": hsl(262, 73%, 70%), + "indigo": hsl(219, 76%, 62%), + "teal": hsl(174, 100%, 40%), + "green": hsl(122, 39%, 60%), + "deep-orange": hsl(14, 100%, 65%), + "brown": hsl(16, 45%, 56%), + // Set neutral colors to indigo + "grey": hsl(219, 66%, 62%), + "blue-grey": hsl(219, 66%, 62%), + "white": hsl(219, 66%, 62%), + "black": hsl(219, 66%, 62%) + ) + { + [data-md-color-scheme="slate"][data-md-color-primary="#{$name}"] { + --md-typeset-a-color: #{$color}; + } + } + + // -------------------------------------------------------------------------- + + // Switching in progress - disable all transitions temporarily + [data-md-color-switching] *, + [data-md-color-switching] *::before, + [data-md-color-switching] *::after { + transition-duration: 0ms !important; // stylelint-disable-line + } +} diff --git a/py-feos/docs/api/eos/eos.md b/py-feos/docs/api/eos/eos.md new file mode 100644 index 000000000..b19364b32 --- /dev/null +++ b/py-feos/docs/api/eos/eos.md @@ -0,0 +1 @@ +# Equation of State \ No newline at end of file diff --git a/py-feos/docs/api/eos/epcsaft.md b/py-feos/docs/api/eos/epcsaft.md new file mode 100644 index 000000000..79a81b640 --- /dev/null +++ b/py-feos/docs/api/eos/epcsaft.md @@ -0,0 +1,28 @@ +# ePC-SAFT + +Electrolyte PC-SAFT equation of state for ionic and non-ionic systems. + +```python +from feos.parameters import Parameters +from feos import EquationOfState + +# Load parameters for an ionic system +parameters = Parameters.from_json( + ["water", "sodium chloride"], + "parameters/epcsaft/example.json" +) + +# Create ePC-SAFT equation of state +epcsaft = EquationOfState.epcsaft( + parameters, + epcsaft_variant="advanced" +) +``` + +::: feos.EquationOfState.epcsaft + options: + show_root_full_path: true + docstring_section_style: spacy + summary: + attributes: false + functions: true \ No newline at end of file diff --git a/py-feos/docs/api/eos/gc_pcsaft.md b/py-feos/docs/api/eos/gc_pcsaft.md new file mode 100644 index 000000000..e23a1afff --- /dev/null +++ b/py-feos/docs/api/eos/gc_pcsaft.md @@ -0,0 +1,26 @@ +# GC-PC-SAFT + +Group contribution PC-SAFT equation of state for molecules described using segments. + +```python +from feos.parameters import GcParameters +from feos import EquationOfState + +# Load group contribution parameters +parameters = GcParameters.from_json_segments( + ["methanol", "ethanol"], + "parameters/gc_pcsaft/chemical_records.json", + "parameters/gc_pcsaft/segments.json" +) + +# Create GC-PC-SAFT equation of state +gc_pcsaft = EquationOfState.gc_pcsaft(parameters) +``` + +::: feos.EquationOfState.gc_pcsaft + options: + show_root_full_path: true + docstring_section_style: spacy + summary: + attributes: false + functions: true \ No newline at end of file diff --git a/py-feos/docs/api/eos/pcsaft.md b/py-feos/docs/api/eos/pcsaft.md new file mode 100644 index 000000000..3dd1bb251 --- /dev/null +++ b/py-feos/docs/api/eos/pcsaft.md @@ -0,0 +1,17 @@ +# PC-SAFT + +```python +from feos.parameters import Parameters +from feos import EquationOfState + +parameters = Parameters(...) +pcsaft = EquationOfState.pcsaft(parameters) +``` + +::: feos.EquationOfState.pcsaft + options: + show_root_full_path: true + docstring_section_style: spacy + summary: + attributes: false + functions: true diff --git a/py-feos/docs/api/eos/pets.md b/py-feos/docs/api/eos/pets.md new file mode 100644 index 000000000..39e405c9f --- /dev/null +++ b/py-feos/docs/api/eos/pets.md @@ -0,0 +1,25 @@ +# PeTS + +Perturbed truncated and shifted Lennard-Jones equation of state. + +```python +from feos.parameters import Parameters +from feos import EquationOfState + +# Load PeTS parameters +parameters = Parameters.from_json( + ["methane", "ethane"], + "parameters/pets/example.json" +) + +# Create PeTS equation of state +pets = EquationOfState.pets(parameters) +``` + +::: feos.EquationOfState.pets + options: + show_root_full_path: true + docstring_section_style: spacy + summary: + attributes: false + functions: true \ No newline at end of file diff --git a/py-feos/docs/api/eos/saftvrmie.md b/py-feos/docs/api/eos/saftvrmie.md new file mode 100644 index 000000000..17d3c2af7 --- /dev/null +++ b/py-feos/docs/api/eos/saftvrmie.md @@ -0,0 +1,17 @@ +# SAFT-VR Mie + +```python +from feos.parameters import Parameters +from feos import EquationOfState + +parameters = Parameters(...) +pcsaft = EquationOfState.saftvrmie(parameters) +``` + +::: feos.EquationOfState.saftvrmie + options: + show_root_full_path: true + docstring_section_style: spacy + summary: + attributes: false + functions: true diff --git a/py-feos/docs/api/eos/saftvrqmie.md b/py-feos/docs/api/eos/saftvrqmie.md new file mode 100644 index 000000000..d5cf6db3c --- /dev/null +++ b/py-feos/docs/api/eos/saftvrqmie.md @@ -0,0 +1,28 @@ +# SAFT-VRQ Mie + +SAFT-VRQ Mie equation of state for quantum fluids with Mie interactions. + +```python +from feos.parameters import Parameters +from feos import EquationOfState + +# Load SAFT-VRQ Mie parameters +parameters = Parameters.from_json( + ["hydrogen", "helium"], + "parameters/saftvrqmie/example.json" +) + +# Create SAFT-VRQ Mie equation of state +saftvrqmie = EquationOfState.saftvrqmie( + parameters, + inc_nonadd_term=True +) +``` + +::: feos.EquationOfState.saftvrqmie + options: + show_root_full_path: true + docstring_section_style: spacy + summary: + attributes: false + functions: true \ No newline at end of file diff --git a/py-feos/docs/api/eos/uvtheory.md b/py-feos/docs/api/eos/uvtheory.md new file mode 100644 index 000000000..74c44144c --- /dev/null +++ b/py-feos/docs/api/eos/uvtheory.md @@ -0,0 +1,28 @@ +# UV Theory + +UV theory equation of state for Mie fluids with different perturbation variants. + +```python +from feos.parameters import Parameters +from feos import EquationOfState + +# Load UV theory parameters +parameters = Parameters.from_json( + ["argon", "methane"], + "parameters/uvtheory/example.json" +) + +# Create UV theory equation of state +uvtheory = EquationOfState.uvtheory( + parameters, + perturbation="WCA" +) +``` + +::: feos.EquationOfState.uvtheory + options: + show_root_full_path: true + docstring_section_style: spacy + summary: + attributes: false + functions: true \ No newline at end of file diff --git a/py-feos/docs/api/index.md b/py-feos/docs/api/index.md new file mode 100644 index 000000000..0d51a6a83 --- /dev/null +++ b/py-feos/docs/api/index.md @@ -0,0 +1,35 @@ +# Introduction + +Working with FeO$_\text{s}$ broadly requires these steps: + +1. **Specifying parameters:** Create a `Parameters` or `GcParameters` object that contains parameters for the substance(s) of interest. +2. **Initializing a model:** That's either an `EquationOfState` or a `HelmholtzEnergyFunctional`. +3. **Specifying thermodynamic conditions:** That's either a `State` object, a `PhaseEquilibrium` object. There are several utility functions that build these objects in bulk, like a `PhaseDiagram`. +4. **Calculate properties:** Once the conditions are specified, you can calcuate properties. FeO$_\text{s}$ caches computations so repeated calls to the same method are very fast. + +FeO$_\text{s}$ uses dimensioned units as in- and ouputs. Make sure you have the `si_units` package installed. + +## Example + +```python +import si_units as si +from feos.parameters import Identifier, PureRecord, Parameters +from feos import Contributions, EquationOfState, State + +# Specify parameters for PC-SAFT +methane = PureRecord( + identifier=Identifier(name="methane"), + molarweight=16.043, + m=1.0, + sigma=3.7039, + epsilon_k=150.03 +) +parameters = Parameters.new_pure(methane) + +# Initialize model: at this point parameters are checked for consistency with the model. +model = EquationOfState.pcsaft(parameters) + +# Specify thermodynamic conditions and calculate property. +state = State(model, temperature=300 * si.KELVIN, pressure=5 * si.BAR) +entropy = state.molar_entropy(Contributions.Residual) +``` diff --git a/py-feos/docs/api/misc.md b/py-feos/docs/api/misc.md new file mode 100644 index 000000000..99b19ef4a --- /dev/null +++ b/py-feos/docs/api/misc.md @@ -0,0 +1,18 @@ +::: feos.Contributions + options: + members: + - IdealGas + - Residual + - Total + summary: + attributes: true + functions: false + +--- + +::: feos.Verbosity + options: + members: [] + summary: + attributes: true + functions: false \ No newline at end of file diff --git a/py-feos/docs/api/parameters/identifier.md b/py-feos/docs/api/parameters/identifier.md new file mode 100644 index 000000000..83a4669ee --- /dev/null +++ b/py-feos/docs/api/parameters/identifier.md @@ -0,0 +1,24 @@ +# Identifier + +::: feos.Identifier + options: + members: + - __init__ + summary: + attributes: false + functions: true + +--- + +::: feos.IdentifierOption + options: + members: + - Cas + - Name + - IupacName + - Smiles + - Inchi + - Formula + summary: + attributes: true + functions: false diff --git a/py-feos/docs/api/parameters/parameters.md b/py-feos/docs/api/parameters/parameters.md new file mode 100644 index 000000000..8dd388136 --- /dev/null +++ b/py-feos/docs/api/parameters/parameters.md @@ -0,0 +1,42 @@ +# Parameters + +```python +from feos.parameters import Parameters + +parameters = Parameters(...) +``` + +!!! info + + Parameters are note evaluated until provided to an equation of state. + +::: feos.Parameters + options: + members: + - from_records + - new_pure + - new_binary + - from_multiple_json + - from_json + - to_json_str + - pure_records + - binary_records + summary: + attributes: true + functions: true + +--- + +::: feos.GcParameters + options: + members: + - from_segments + - from_json_segments + - from_smiles + - from_json_smiles + - chemical_records + - segment_records + - binary_segment_records + summary: + attributes: true + functions: true diff --git a/py-feos/docs/api/parameters/records.md b/py-feos/docs/api/parameters/records.md new file mode 100644 index 000000000..732df5c5f --- /dev/null +++ b/py-feos/docs/api/parameters/records.md @@ -0,0 +1,61 @@ +# Records + +Record classes represent parameter sets for individual components, segments, and their interactions. + +::: feos.PureRecord + options: + members: + - from_json_str + - to_json_str + - to_dict + summary: + attributes: true + functions: true + +::: feos.SegmentRecord + options: + members: + - from_json_str + - to_json_str + - to_dict + - from_json + summary: + attributes: true + functions: true + +::: feos.BinaryRecord + options: + members: + - from_json_str + - to_json_str + - to_dict + summary: + attributes: true + functions: true + +::: feos.BinarySegmentRecord + options: + members: + - from_json_str + - to_json_str + summary: + attributes: true + functions: true + +::: feos.ChemicalRecord + options: + members: + - from_smiles + summary: + attributes: true + functions: true + +::: feos.SmartsRecord + options: + members: + - from_json_str + - to_json_str + - from_json + summary: + attributes: true + functions: true diff --git a/py-feos/docs/api/state/state.md b/py-feos/docs/api/state/state.md new file mode 100644 index 000000000..8e344e25b --- /dev/null +++ b/py-feos/docs/api/state/state.md @@ -0,0 +1,102 @@ +# State + +The `State` class represents a thermodynamic state at given conditions and provides access to all thermodynamic properties. + +```python +# Todo example +``` + +::: feos.State + options: + members: + - __init__ + - critical_point + - critical_point_pure + - critical_point_binary + - spinodal + - dp_dv + - dp_drho + - dp_dt + - dp_dni + - d2p_dv2 + - d2p_drho2 + - partial_molar_volume + - chemical_potential + - chemical_potential_contributions + - dmu_dt + - dmu_dni + - ln_phi + - ln_phi_pure_liquid + - ln_symmetric_activity_coefficient + - dln_phi_dt + - dln_phi_dp + - dln_phi_dnj + - thermodynamic_factor + - henrys_law_constant + - henrys_law_constant_binary + - molar_isochoric_heat_capacity + - dc_v_dt + - molar_isobaric_heat_capacity + - entropy + - ds_dt + - molar_entropy + - partial_molar_entropy + - enthalpy + - molar_enthalpy + - partial_molar_enthalpy + - helmholtz_energy + - molar_helmholtz_energy + - residual_helmholtz_energy_contributions + - gibbs_energy + - molar_gibbs_energy + - internal_energy + - molar_internal_energy + - joule_thomson + - isentropic_compressibility + - isothermal_compressibility + - isenthalpic_compressibility + - thermal_expansivity + - grueneisen_parameter + - structure_factor + - speed_of_sound + - mass + - total_mass + - mass_density + - massfracs + - specific_helmholtz_energy + - specific_entropy + - specific_internal_energy + - specific_gibbs_energy + - specific_enthalpy + - specific_isochoric_heat_capacity + - specific_isobaric_heat_capacity + - stability_analysis + - is_stable + summary: + attributes: true + functions: true +--- + +::: feos.StateVec + options: + members: + - __init__ + - __len__ + - __getitem__ + - molar_entropy + - specific_entropy + - molar_enthalpy + - specific_enthalpy + - to_dict + - temperature + - pressure + - compressibility + - density + - moles + - molefracs + - mass_density + - massfracs + summary: + attributes: true + functions: true + diff --git a/py-feos/docs/examples/basic-properties.md b/py-feos/docs/examples/basic-properties.md new file mode 100644 index 000000000..9e933ac07 --- /dev/null +++ b/py-feos/docs/examples/basic-properties.md @@ -0,0 +1,236 @@ +# Basic Properties + +This example shows how to calculate basic thermodynamic properties for pure components. + +## Pure Component Properties + +Let's calculate various properties for methane using PC-SAFT: + +```python +from feos import EquationOfState, State +from feos.parameters import PureRecord, Identifier, Parameters + +# PC-SAFT parameters for methane +methane = PureRecord( + identifier=Identifier(name="methane", cas="74-82-8"), + molarweight=16.04, + m=1.0, + sigma=3.7039, + epsilon_k=150.03, +) + +# Create equation of state +parameters = Parameters.new_pure(methane) +eos = EquationOfState.pcsaft(parameters) + +# Define conditions +T = 298.15 # K (room temperature) +p = 101325 # Pa (1 atm) + +# Create state +state = State(eos, temperature=T, pressure=p) + +print("=== Basic Properties ===") +print(f"Temperature: {state.temperature} K") +print(f"Pressure: {state.pressure()} Pa") +print(f"Density: {state.density}") +print(f"Molar volume: {state.molar_volume()}") +print(f"Compressibility factor: {state.compressibility()}") +``` + +## Derived Properties + +Calculate heat capacities, speed of sound, and other derived properties: + +```python +print("\n=== Derived Properties ===") +print(f"Isobaric heat capacity (Cp): {state.molar_isobaric_heat_capacity()}") +print(f"Isochoric heat capacity (Cv): {state.molar_isochoric_heat_capacity()}") +print(f"Heat capacity ratio (γ): {state.isentropic_compressibility()}") +print(f"Speed of sound: {state.speed_of_sound()}") +print(f"Joule-Thomson coefficient: {state.joule_thomson()}") +``` + +## Chemical Potential and Fugacity + +For phase equilibrium calculations, chemical potential and fugacity are important: + +```python +print("\n=== Chemical Properties ===") +print(f"Chemical potential: {state.chemical_potential()}") +print(f"Fugacity coefficient (ln φ): {state.ln_phi()}") +print(f"Fugacity: {state.fugacity()}") +``` + +## Properties at Different Conditions + +Calculate properties at various temperatures and pressures: + +```python +print("\n=== Properties at Different Conditions ===") + +# Temperature range at constant pressure +temperatures = [250, 300, 350, 400] # K +pressure = 101325 # Pa + +print("T [K]\tDensity [kg/m³]\tCp [J/mol/K]") +print("-" * 40) + +for T in temperatures: + state = State(eos, temperature=T, pressure=pressure) + density = state.density.value # Get numerical value + cp = state.molar_isobaric_heat_capacity().value + print(f"{T}\t{density:.2f}\t\t{cp:.2f}") +``` + +## Critical Properties + +Calculate critical point properties: + +```python +print("\n=== Critical Properties ===") + +# Calculate critical point +critical = State.critical_point_pure(eos) + +print(f"Critical temperature: {critical.temperature} K") +print(f"Critical pressure: {critical.pressure()} Pa") +print(f"Critical density: {critical.density}") +print(f"Critical compressibility: {critical.compressibility()}") +``` + +## Vapor Pressure + +Calculate vapor pressure at different temperatures: + +```python +print("\n=== Vapor Pressure ===") + +temperatures = [200, 250, 300, 350] # K + +print("T [K]\tP_sat [Pa]") +print("-" * 20) + +for T in temperatures: + try: + p_sat = State.vapor_pressure(eos, T) + print(f"{T}\t{p_sat.value:.0f}") + except RuntimeError: + print(f"{T}\tN/A (above critical)") +``` + +## Different Models + +Compare results from different equation of state models: + +```python +print("\n=== Model Comparison ===") + +# Create different EoS for the same component +eos_pcsaft = EquationOfState.pcsaft(parameters) + +# For comparison, we would need PETS parameters: +# pets_params = PureRecord( +# identifier=Identifier(name="methane"), +# molarweight=16.04, +# sigma=3.7, # Different parameter meaning in PETS +# epsilon_k=150, +# ) +# eos_pets = EquationOfState.pets(Parameters.new_pure(pets_params)) + +# Calculate density at same conditions +T, p = 298.15, 101325 +state_pcsaft = State(eos_pcsaft, temperature=T, pressure=p) + +print(f"PC-SAFT density: {state_pcsaft.density}") +# print(f"PETS density: {state_pets.density}") +``` + +## Working with Units + +Understanding how to work with units in FeOs: + +```python +print("\n=== Working with Units ===") + +state = State(eos, temperature=298.15, pressure=101325) +density = state.density + +# Access different aspects of the quantity +print(f"Full quantity: {density}") +print(f"Numerical value: {density.value}") +print(f"Units: {density.units}") + +# Convert to different units (if needed for calculations) +density_value = density.value # kg/m³ in SI units +print(f"Density in g/cm³: {density_value / 1000:.6f}") +``` + +## Error Handling + +Handle common calculation errors: + +```python +print("\n=== Error Handling ===") + +try: + # Try an extreme condition that might fail + extreme_state = State(eos, temperature=1000, pressure=1e8) # Very high T and p + print(f"Extreme condition density: {extreme_state.density}") +except RuntimeError as e: + print(f"Calculation failed: {e}") + +try: + # Try vapor pressure above critical temperature + high_temp_vp = State.vapor_pressure(eos, 500) # Above critical + print(f"Vapor pressure: {high_temp_vp}") +except RuntimeError as e: + print(f"Vapor pressure calculation failed: {e}") +``` + +## Complete Example Script + +Here's a complete script you can run: + +```python +#!/usr/bin/env python3 +""" +Basic thermodynamic property calculations with FeOs +""" + +from feos import EquationOfState, State +from feos.parameters import PureRecord, Identifier, Parameters + +def main(): + # Define component + methane = PureRecord( + identifier=Identifier(name="methane", cas="74-82-8"), + molarweight=16.04, + m=1.0, + sigma=3.7039, + epsilon_k=150.03, + ) + + # Create EoS + parameters = Parameters.new_pure(methane) + eos = EquationOfState.pcsaft(parameters) + + # Calculate properties at standard conditions + state = State(eos, temperature=298.15, pressure=101325) + + print("Methane properties at 298.15 K, 1 atm:") + print(f" Density: {state.density}") + print(f" Compressibility: {state.compressibility():.4f}") + print(f" Heat capacity (Cp): {state.molar_isobaric_heat_capacity()}") + + # Critical point + critical = State.critical_point_pure(eos) + print(f"\nCritical point:") + print(f" Temperature: {critical.temperature}") + print(f" Pressure: {critical.pressure()}") + +if __name__ == "__main__": + main() +``` + +This example demonstrates the fundamental property calculations available in FeOs. Next, explore [pure component VLE](pure-component-vle.md) for vapor-liquid equilibrium calculations. \ No newline at end of file diff --git a/py-feos/docs/examples/parameter-loading.md b/py-feos/docs/examples/parameter-loading.md new file mode 100644 index 000000000..a2a0c722b --- /dev/null +++ b/py-feos/docs/examples/parameter-loading.md @@ -0,0 +1,359 @@ +# Parameter Loading + +This example shows different ways to load and manage parameters in FeOs. + +## Loading Parameters from JSON + +FeOs can load parameters from JSON files, which is convenient for working with parameter databases: + +```python +from feos.parameters import Parameters, IdentifierOption + +# Load pure component parameters from JSON file +# (Note: You'll need actual JSON parameter files for this to work) +try: + parameters = Parameters.from_json( + substances=["methane", "ethane"], + pure_path="path/to/pure_parameters.json", + identifier_option=IdentifierOption.Name + ) + print("Parameters loaded successfully") +except FileNotFoundError: + print("Parameter file not found - creating parameters manually") + # Fall back to manual parameter creation + # ... (manual creation code) +``` + +## Creating Parameters Manually + +When you don't have JSON files, create parameters directly: + +```python +from feos.parameters import PureRecord, BinaryRecord, Identifier, Parameters + +# Create pure component records +methane = PureRecord( + identifier=Identifier(name="methane", cas="74-82-8"), + molarweight=16.04, + m=1.0, + sigma=3.7039, + epsilon_k=150.03, +) + +ethane = PureRecord( + identifier=Identifier(name="ethane", cas="74-84-0"), + molarweight=30.07, + m=1.6069, + sigma=3.5206, + epsilon_k=191.42, +) + +# Create binary interaction parameter +binary = BinaryRecord( + id1=Identifier(name="methane"), + id2=Identifier(name="ethane"), + k_ij=0.0, # No binary interaction +) + +# Combine into Parameters object +parameters = Parameters.from_records( + pure_records=[methane, ethane], + binary_records=[binary] +) + +print(f"Created parameters for {len(parameters.pure_records)} components") +``` + +## Working with Identifiers + +FeOs supports multiple identifier types for flexibility: + +```python +from feos.parameters import Identifier, IdentifierOption + +# Different ways to create identifiers +id_by_name = Identifier(name="water") +id_by_cas = Identifier(cas="7732-18-5") +id_by_smiles = Identifier(smiles="O") + +# Complete identifier with multiple fields +complete_id = Identifier( + name="water", + cas="7732-18-5", + iupac_name="oxidane", + smiles="O", + formula="H2O" +) + +print(f"Identifier: {complete_id}") + +# When loading from JSON, specify which identifier to match on +parameters_by_name = Parameters.from_json( + substances=["water"], + pure_path="parameters.json", + identifier_option=IdentifierOption.Name +) + +parameters_by_cas = Parameters.from_json( + substances=["7732-18-5"], + pure_path="parameters.json", + identifier_option=IdentifierOption.Cas +) +``` + +## Parameter Inspection + +Examine loaded parameters: + +```python +# Inspect pure records +for record in parameters.pure_records: + print(f"Component: {record.identifier.name}") + print(f" Molar weight: {record.molarweight} g/mol") + print(f" Model parameters: {record.model_record}") + print() + +# Inspect binary records +if hasattr(parameters, 'binary_records') and parameters.binary_records: + for binary in parameters.binary_records: + print(f"Binary interaction: {binary.id1.name} - {binary.id2.name}") + print(f" Parameters: {binary.model_record}") +``` + +## JSON Export + +Export parameters to JSON format: + +```python +# Convert parameters to JSON strings +pure_json, binary_json = parameters.to_json_str(pretty=True) + +print("Pure component parameters (JSON):") +print(pure_json[:200] + "..." if len(pure_json) > 200 else pure_json) + +if binary_json: + print("\nBinary parameters (JSON):") + print(binary_json[:200] + "..." if len(binary_json) > 200 else binary_json) + +# Save to files +with open("pure_params.json", "w") as f: + f.write(pure_json) + +if binary_json: + with open("binary_params.json", "w") as f: + f.write(binary_json) +``` + +## Different Parameter Sets + +Different equation of state models require different parameters: + +```python +# PC-SAFT parameters (associating component) +water_pcsaft = PureRecord( + identifier=Identifier(name="water"), + molarweight=18.015, + m=1.2047, # segment number + sigma=2.7927, # segment diameter [Å] + epsilon_k=353.944, # dispersion energy [K] + kappa_ab=0.0451, # association volume + epsilon_k_ab=2500.7, # association energy [K] + na=1, # association sites A + nb=1, # association sites B +) + +# PETS parameters (simpler model) +argon_pets = PureRecord( + identifier=Identifier(name="argon"), + molarweight=39.948, + sigma=3.4, # hard sphere diameter [Å] + epsilon_k=120, # dispersion energy [K] +) + +# Create different parameter sets +pcsaft_params = Parameters.new_pure(water_pcsaft) +pets_params = Parameters.new_pure(argon_pets) + +print("Created PC-SAFT parameters for water") +print("Created PETS parameters for argon") +``` + +## Group Contribution Parameters + +For group contribution methods like GC-PC-SAFT: + +```python +from feos.parameters import ( + GcParameters, ChemicalRecord, SegmentRecord, + BinarySegmentRecord, SmartsRecord +) + +# Define segments (functional groups) +ch3_segment = SegmentRecord( + identifier="CH3", + molarweight=15.035, + m=0.7745, + sigma=3.7729, + epsilon_k=229.0, +) + +ch2_segment = SegmentRecord( + identifier="CH2", + molarweight=14.027, + m=0.6744, + sigma=3.8395, + epsilon_k=247.0, +) + +# Define chemical structure (e.g., propane = CH3-CH2-CH3) +propane_structure = ChemicalRecord( + identifier=Identifier(name="propane"), + segments=["CH3", "CH2", "CH3"], + bonds=[[0, 1], [1, 2]] # Connect segments +) + +# Binary segment interaction +ch3_ch2_binary = BinarySegmentRecord( + id1="CH3", + id2="CH2", + model_record=0.0 # No interaction +) + +# Create GC parameters +gc_params = GcParameters.from_segments( + chemical_records=[propane_structure], + segment_records=[ch3_segment, ch2_segment], + binary_segment_records=[ch3_ch2_binary] +) + +print("Created group contribution parameters for propane") +``` + +## SMILES-based Parameter Generation + +Generate parameters from SMILES strings (requires rdkit): + +```python +try: + # This requires rdkit to be installed + from feos.parameters import SmartsRecord + + # Define SMARTS patterns for fragmentation + smarts = [ + SmartsRecord(group="CH3", smarts="[CH3]"), + SmartsRecord(group="CH2", smarts="[CH2]"), + # ... more patterns + ] + + # Generate parameters from SMILES + gc_params_smiles = GcParameters.from_smiles( + identifier="CCC", # SMILES for propane + smarts_records=smarts, + segment_records=[ch3_segment, ch2_segment], + binary_segment_records=[ch3_ch2_binary] + ) + + print("Generated parameters from SMILES") + +except ImportError: + print("rdkit not available - cannot generate from SMILES") +except Exception as e: + print(f"SMILES generation failed: {e}") +``` + +## Parameter Database Management + +Best practices for managing parameter databases: + +```python +import json +from pathlib import Path + +# Create a parameter database structure +def create_parameter_database(): + """Create a simple parameter database""" + + # Pure component database + pure_db = [ + { + "identifier": {"name": "methane", "cas": "74-82-8"}, + "molarweight": 16.04, + "model_record": { + "m": 1.0, + "sigma": 3.7039, + "epsilon_k": 150.03 + } + }, + { + "identifier": {"name": "ethane", "cas": "74-84-0"}, + "molarweight": 30.07, + "model_record": { + "m": 1.6069, + "sigma": 3.5206, + "epsilon_k": 191.42 + } + } + ] + + # Binary interaction database + binary_db = [ + { + "id1": {"name": "methane"}, + "id2": {"name": "ethane"}, + "model_record": {"k_ij": 0.0} + } + ] + + # Save to files + Path("pure_parameters.json").write_text( + json.dumps(pure_db, indent=2) + ) + Path("binary_parameters.json").write_text( + json.dumps(binary_db, indent=2) + ) + + return "pure_parameters.json", "binary_parameters.json" + +# Create and use database +pure_file, binary_file = create_parameter_database() + +# Load parameters from the database +params = Parameters.from_json( + substances=["methane", "ethane"], + pure_path=pure_file, + binary_path=binary_file, + identifier_option=IdentifierOption.Name +) + +print(f"Loaded parameters from database: {len(params.pure_records)} components") +``` + +## Error Handling + +Handle common parameter loading errors: + +```python +try: + # Try to load non-existent component + params = Parameters.from_json( + substances=["nonexistent_component"], + pure_path="pure_parameters.json", + identifier_option=IdentifierOption.Name + ) +except Exception as e: + print(f"Failed to load parameters: {e}") + +try: + # Invalid parameter values + invalid_record = PureRecord( + identifier=Identifier(name="invalid"), + molarweight=-10, # Invalid negative mass + m=1.0, + sigma=3.0, + epsilon_k=100, + ) +except ValueError as e: + print(f"Invalid parameter: {e}") +``` + +This example covers the various ways to work with parameters in FeOs. Next, explore [binary mixtures](binary-mixtures.md) to see how to use these parameters for mixture calculations. \ No newline at end of file diff --git a/py-feos/docs/getting-started/concepts.md b/py-feos/docs/getting-started/concepts.md new file mode 100644 index 000000000..e57edee83 --- /dev/null +++ b/py-feos/docs/getting-started/concepts.md @@ -0,0 +1,265 @@ +# Core Concepts + +Understanding the core concepts of FeOs will help you use the library effectively. This page explains the main components and how they work together. + +## Overview + +FeOs is built around a few key concepts: + +- **Parameters** define the molecular properties of chemical components +- **Equations of State (EoS)** describe how those components behave thermodynamically +- **States** represent specific thermodynamic conditions +- **Phase Equilibria** calculations find coexisting phases +- **DFT profiles** describe inhomogeneous fluid systems + +## Parameters + +Parameters contain the molecular information needed for thermodynamic calculations. + +### Pure Component Parameters + +Parameters for a single component are stored in a `PureRecord`: + +```python +from feos.parameters import PureRecord, Identifier + +# PC-SAFT parameters for water +water = PureRecord( + identifier=Identifier(name="water", cas="7732-18-5"), + molarweight=18.015, + m=1.2047, # segment number + sigma=2.7927, # segment diameter [Å] + epsilon_k=353.944, # dispersion energy [K] + kappa_ab=0.0451, # association volume + epsilon_k_ab=2500.7, # association energy [K] + na=1, # association sites of type A + nb=1, # association sites of type B +) +``` + +### Parameter Meanings + +Different equation of state models use different parameters: + +**PC-SAFT parameters:** +- `m`: Number of segments in the chain +- `sigma`: Segment diameter [Å] +- `epsilon_k`: Dispersion energy [K] +- `kappa_ab`: Association volume parameter +- `epsilon_k_ab`: Association energy [K] +- `na`, `nb`: Number of association sites + +**PETS parameters:** +- `sigma`: Hard sphere diameter [Å] +- `epsilon_k`: Dispersion energy [K] + +### Binary Interaction Parameters + +For mixtures, you may need binary interaction parameters: + +```python +from feos.parameters import BinaryRecord + +# Binary interaction for water-ethanol +binary = BinaryRecord( + id1=Identifier(name="water"), + id2=Identifier(name="ethanol"), + k_ij=0.02, # Binary interaction parameter +) +``` + +## Equations of State + +An equation of state describes the relationship between pressure, volume, temperature, and composition. + +### Creating an EoS + +```python +from feos import EquationOfState +from feos.parameters import Parameters + +# Create parameters object +parameters = Parameters.new_pure(water) + +# Create PC-SAFT equation of state +eos = EquationOfState.pcsaft(parameters) +``` + +### Available Models + +FeOs supports several equation of state models: + +- **PC-SAFT**: Perturbed-chain statistical associating fluid theory +- **EPC-SAFT**: Electrolyte PC-SAFT (for ionic systems) +- **GC-PC-SAFT**: Group contribution PC-SAFT +- **PETS**: Perturbed truncated and shifted Lennard-Jones +- **SAFT-VR Mie**: Variable range SAFT with Mie potential +- **UV Theory**: For Mie fluids + +Each model has its own parameter requirements and applicability range. + +## States + +A `State` represents a thermodynamic condition with known temperature, pressure, density, or other properties. + +### Creating States + +States can be created from different combinations of properties: + +```python +from feos import State + +# From temperature and pressure +state1 = State(eos, temperature=298.15, pressure=101325) + +# From temperature and density +state2 = State(eos, temperature=298.15, density=1000) + +# From temperature and molar volume +state3 = State(eos, temperature=298.15, volume=0.018) +``` + +### State Properties + +Once you have a state, you can calculate many properties: + +```python +# Basic properties +print(f"Pressure: {state.pressure()}") +print(f"Density: {state.density}") +print(f"Molar volume: {state.molar_volume()}") + +# Derived properties +print(f"Compressibility: {state.compressibility()}") +print(f"Speed of sound: {state.speed_of_sound()}") +print(f"Heat capacity (cp): {state.molar_isobaric_heat_capacity()}") +print(f"Heat capacity (cv): {state.molar_isochoric_heat_capacity()}") + +# Chemical potential and fugacity +print(f"Chemical potential: {state.chemical_potential()}") +print(f"Fugacity coefficient: {state.ln_phi()}") +``` + +## Phase Equilibria + +Phase equilibria calculations find conditions where multiple phases coexist. + +### Vapor-Liquid Equilibrium + +```python +from feos import PhaseEquilibrium + +# Calculate vapor pressure at given temperature +vle = PhaseEquilibrium.pure(eos, temperature=298.15, pressure=None) +print(f"Vapor pressure: {vle.vapor().pressure()}") +print(f"Liquid density: {vle.liquid().density}") +print(f"Vapor density: {vle.vapor().density}") +``` + +### Critical Points + +```python +# Critical point of pure component +critical = State.critical_point_pure(eos) +print(f"Critical temperature: {critical.temperature}") +print(f"Critical pressure: {critical.pressure()}") +print(f"Critical density: {critical.density}") +``` + +## Mixtures + +For multicomponent systems, you need parameters for all components plus any binary interactions: + +```python +# Multiple pure components +water_record = PureRecord(...) # water parameters +ethanol_record = PureRecord(...) # ethanol parameters + +# Binary interaction +binary_record = BinaryRecord(...) + +# Create mixture parameters +mixture_params = Parameters.from_records( + pure_records=[water_record, ethanol_record], + binary_records=[binary_record] +) + +# Create mixture EoS +mixture_eos = EquationOfState.pcsaft(mixture_params) + +# Create mixture state (need mole fractions) +mixture_state = State( + mixture_eos, + temperature=298.15, + pressure=101325, + molefracs=[0.5, 0.5] # 50% water, 50% ethanol +) +``` + +## Density Functional Theory (DFT) + +DFT calculations describe inhomogeneous fluid systems like interfaces and confined fluids. + +### Helmholtz Energy Functionals + +```python +from feos.dft import HelmholtzEnergyFunctional + +# Create functional (requires DFT-capable EoS) +functional = HelmholtzEnergyFunctional.pcsaft(parameters) +``` + +### Density Profiles + +DFT calculates density profiles in various geometries: + +```python +from feos.dft import PlanarInterface + +# Vapor-liquid interface +vle = PhaseEquilibrium.pure(eos, temperature=400, pressure=None) +interface = PlanarInterface.from_tanh(vle, n_grid=200, l_grid=10e-9) +interface.solve() + +print(f"Surface tension: {interface.surface_tension()}") +``` + +## Units and Quantities + +FeOs uses the `si-units` package for physical quantities: + +```python +# All inputs are in SI base units +T = 298.15 # Kelvin +p = 101325 # Pascal + +# Outputs are SIObject instances with units +density = state.density +print(f"Value: {density.value}") # Numerical value +print(f"Units: {density.units}") # Units string +print(f"Full: {density}") # Value with units +``` + +## Error Handling + +FeOs calculations can fail for various reasons: + +```python +try: + state = State(eos, temperature=298.15, pressure=1e10) # Very high pressure +except RuntimeError as e: + print(f"Calculation failed: {e}") +``` + +Common failure modes: +- **Convergence failure**: Iteration didn't converge +- **Invalid conditions**: Unphysical temperature/pressure +- **Phase detection**: Multiple solutions exist + +## Next Steps + +Now that you understand the core concepts: + +- **[Working with States](../user-guide/working-with-states.md)** - Advanced state calculations +- **[Parameters](../user-guide/parameters.md)** - Parameter management and databases +- **[Examples](../examples/basic-properties.md)** - Practical code examples \ No newline at end of file diff --git a/py-feos/docs/getting-started/installation.md b/py-feos/docs/getting-started/installation.md new file mode 100644 index 000000000..7e15c68ec --- /dev/null +++ b/py-feos/docs/getting-started/installation.md @@ -0,0 +1,151 @@ +# Installation + +This guide covers different ways to install the FeOs Python package. + +## Requirements + +- **Python**: 3.12 or higher +- **Operating System**: Windows, macOS, or Linux + +## Install from PyPI + +The easiest way to install FeOs is from PyPI using pip: + +```bash +pip install feos +``` + +This will install the latest stable release with all equation of state models and DFT functionality included. + +## Install with uv (Recommended for Development) + +If you're using [uv](https://docs.astral.sh/uv/) for Python package management: + +```bash +uv add feos +``` + +## Install from Source + +### Prerequisites for Building from Source + +To build from source, you need: + +- **Rust toolchain**: Install from [rustup.rs](https://rustup.rs/) +- **maturin**: Python package for building Rust extensions + +```bash +# Install Rust +curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh + +# Install maturin +pip install maturin +``` + +### Build and Install + +Clone the repository and build the package: + +```bash +git clone https://github.com/feos-org/feos.git +cd feos/py-feos + +# Install in development mode +maturin develop --release + +# Or build a wheel +maturin build --release +``` + +### Building with Specific Features + +You can build with only specific equation of state models: + +```bash +# Build with only PC-SAFT +maturin develop --release --features "python pcsaft" + +# Build with PC-SAFT and DFT +maturin develop --release --features "python pcsaft dft" +``` + +Available features: +- `pcsaft` - PC-SAFT equation of state +- `epcsaft` - Electrolyte PC-SAFT +- `pets` - PETS equation of state +- `saftvrmie` - SAFT-VR Mie +- `saftvrqmie` - SAFT-VR quantum Mie +- `uvtheory` - UV Theory +- `dft` - Density functional theory +- `all_models` - All available models (default) + +## Verify Installation + +Test your installation by running: + +```python +import feos +print(f"FeOs version: {feos.version}") + +# Test basic functionality +from feos import EquationOfState +from feos.parameters import PureRecord, Identifier, Parameters + +# This should run without errors +record = PureRecord( + Identifier(name="methane"), + molarweight=16.04, + m=1.0, + sigma=3.7039, + epsilon_k=150.03, +) +parameters = Parameters.new_pure(record) +eos = EquationOfState.pcsaft(parameters) +print("✓ FeOs installed successfully!") +``` + +## Optional Dependencies + +For enhanced functionality, you may want to install additional packages: + +```bash +# For Jupyter notebook support +pip install jupyter ipykernel + +# For plotting (used in some examples) +pip install matplotlib + +# For RDKit integration (SMILES support) +pip install rdkit +``` + +## Troubleshooting + +### Import Errors + +If you get import errors, ensure that: + +1. You're using the correct Python environment +2. The package was installed correctly (`pip list | grep feos`) +3. Your Python version is 3.12 or higher + +### Build Errors + +If building from source fails: + +1. **Check Rust installation**: `rustc --version` +2. **Update Rust**: `rustup update` +3. **Check maturin**: `maturin --version` +4. **Clean build**: Remove `target/` directory and rebuild + +### Performance Issues + +For optimal performance: + +1. **Use release builds**: Always include `--release` flag when building +2. **Enable LTO**: Use `maturin build --profile=release-lto` for maximum optimization +3. **Check CPU features**: The package is optimized for modern CPUs + +## Next Steps + +Once installed, continue with the [Quick Start](quickstart.md) guide to learn the basics of using FeOs. \ No newline at end of file diff --git a/py-feos/docs/getting-started/quickstart.md b/py-feos/docs/getting-started/quickstart.md new file mode 100644 index 000000000..6156904d5 --- /dev/null +++ b/py-feos/docs/getting-started/quickstart.md @@ -0,0 +1,185 @@ +# Quick Start + +This guide will get you up and running with FeOs in just a few minutes. We'll walk through the basic workflow of setting up parameters, creating an equation of state, and calculating thermodynamic properties. + +## Basic Workflow + +The typical FeOs workflow follows these steps: + +1. **Define parameters** for your chemical components +2. **Create an equation of state** using those parameters +3. **Define thermodynamic states** at specific conditions +4. **Calculate properties** for those states + +Let's walk through each step with a concrete example. + +## Step 1: Import FeOs + +Start by importing the main components: + +```python +from feos import EquationOfState, State +from feos.parameters import PureRecord, Identifier, Parameters +``` + +## Step 2: Define Component Parameters + +Create parameters for a pure component. Let's use methanol as an example: + +```python +# PC-SAFT parameters for methanol (Gross & Sadowski, 2002) +methanol_record = PureRecord( + identifier=Identifier(name="methanol", cas="67-56-1"), + molarweight=32.04, # g/mol + m=1.5255, # segment number + sigma=3.23, # segment diameter [Å] + epsilon_k=188.9, # dispersion energy [K] + kappa_ab=0.035176, # association volume + epsilon_k_ab=2899.5, # association energy [K] + na=1, # number of association sites A + nb=1, # number of association sites B +) + +# Create parameters object +parameters = Parameters.new_pure(methanol_record) +print(f"Parameters created for: {methanol_record.identifier.name}") +``` + +## Step 3: Create Equation of State + +Use the parameters to create a PC-SAFT equation of state: + +```python +# Create PC-SAFT equation of state +eos = EquationOfState.pcsaft(parameters) +print("PC-SAFT equation of state created") +``` + +## Step 4: Calculate Properties + +Now you can calculate thermodynamic properties at different conditions. + +### Critical Point + +```python +# Calculate critical point +critical_point = State.critical_point_pure(eos) + +print(f"Critical temperature: {critical_point.temperature}") +print(f"Critical pressure: {critical_point.pressure()}") +print(f"Critical density: {critical_point.density}") +``` + +### Properties at Specific Conditions + +```python +# Create state at room temperature and atmospheric pressure +T = 298.15 # K +p = 101325 # Pa +room_temp_state = State(eos, temperature=T, pressure=p) + +print(f"Density at {T} K, {p} Pa: {room_temp_state.density}") +print(f"Molar volume: {room_temp_state.molar_volume()}") +print(f"Compressibility factor: {room_temp_state.compressibility()}") +``` + +### Vapor Pressure + +```python +# Calculate vapor pressure at room temperature +vapor_pressure = State.vapor_pressure(eos, T) +print(f"Vapor pressure at {T} K: {vapor_pressure}") + +# Or calculate boiling temperature at atmospheric pressure +boiling_temp = State.boiling_temperature(eos, p) +print(f"Boiling temperature at {p} Pa: {boiling_temp}") +``` + +## Complete Example + +Here's the complete code from this quick start: + +```python +from feos import EquationOfState, State +from feos.parameters import PureRecord, Identifier, Parameters + +# Define methanol parameters +methanol_record = PureRecord( + identifier=Identifier(name="methanol", cas="67-56-1"), + molarweight=32.04, + m=1.5255, + sigma=3.23, + epsilon_k=188.9, + kappa_ab=0.035176, + epsilon_k_ab=2899.5, + na=1, + nb=1, +) + +# Create equation of state +parameters = Parameters.new_pure(methanol_record) +eos = EquationOfState.pcsaft(parameters) + +# Calculate critical point +critical_point = State.critical_point_pure(eos) +print(f"Critical point: T = {critical_point.temperature}, p = {critical_point.pressure()}") + +# Properties at room conditions +T = 298.15 # K +p = 101325 # Pa +state = State(eos, temperature=T, pressure=p) +print(f"Density at room conditions: {state.density}") + +# Vapor pressure +pv = State.vapor_pressure(eos, T) +print(f"Vapor pressure at {T} K: {pv}") +``` + +## Working with Units + +FeOs uses the `si-units` package for unit handling. All quantities are returned as `SIObject` instances: + +```python +# Temperature and pressure are unitless in SI base units (K and Pa) +T = 298.15 # K +p = 101325 # Pa + +state = State(eos, temperature=T, pressure=p) + +# Properties are returned with units +density = state.density +print(f"Density: {density}") # Will show value with units +print(f"Density value: {density.value}") # Just the numerical value +print(f"Density units: {density.units}") # Just the units +``` + +## Different Equation of State Models + +FeOs supports multiple equation of state models. Here's how to use different ones: + +```python +# PC-SAFT (as shown above) +eos_pcsaft = EquationOfState.pcsaft(parameters) + +# If you have PETS parameters +# eos_pets = EquationOfState.pets(pets_parameters) + +# If you have SAFT-VR Mie parameters +# eos_saftvrmie = EquationOfState.saftvrmie(saftvrmie_parameters) +``` + +## Next Steps + +Now that you understand the basics, explore: + +- **[Core Concepts](concepts.md)** - Deeper understanding of FeOs components +- **[Working with States](../user-guide/working-with-states.md)** - More advanced state calculations +- **[Parameters](../user-guide/parameters.md)** - Loading parameters from JSON files +- **[Examples](../examples/basic-properties.md)** - More practical examples + +## Common Issues + +- **Unit confusion**: Remember that temperature is in Kelvin and pressure in Pascal +- **Parameter units**: Molecular parameters have specific units (see parameter documentation) +- **Association**: Not all components have association sites (set `na=0, nb=0` for non-associating) +- **Convergence**: Some calculations may not converge for extreme conditions \ No newline at end of file diff --git a/py-feos/docs/hooks/typst_math.py b/py-feos/docs/hooks/typst_math.py new file mode 100644 index 000000000..0e317f776 --- /dev/null +++ b/py-feos/docs/hooks/typst_math.py @@ -0,0 +1,131 @@ +"""Render math with typst. + +# See + +https://forum.typst.app/t/guide-render-typst-math-in-mkdocs/4326 + +## Usage + +1. Install the markdown extensions pymdownx.arithmatex. +2. Add `math: typst` to pages' metadata. + +## Requirements + +- typst + +""" + +from __future__ import annotations + +import html +import re +from functools import cache +from subprocess import CalledProcessError, run +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from mkdocs.config.defaults import MkDocsConfig + from mkdocs.structure.files import Files + from mkdocs.structure.pages import Page + + +def should_render(page: Page) -> bool: + return page.meta.get("math") == "typst" + + +def on_page_markdown( + markdown: str, page: Page, config: MkDocsConfig, files: Files +) -> str | None: + if should_render(page): + assert "pymdownx.arithmatex" in config.markdown_extensions, ( + "Missing pymdownx.arithmatex in config.markdown_extensions. " + "Setting `math: typst` requires it to parse markdown." + ) + + +def on_post_page(output: str, page: Page, config: MkDocsConfig) -> str | None: + if should_render(page): + output = re.sub( + r'(.+?)', render_inline_math, output + ) + + output = re.sub( + r'
(.+?)
', + render_block_math, + output, + flags=re.MULTILINE | re.DOTALL, + ) + return output + + +def render_inline_math(match: re.Match[str]) -> str: + src = html.unescape(match.group(1)).removeprefix( + R"\(").removesuffix(R"\)").strip() + typ = f"${src}$" + return ( + '' + + fix_svg(typst_compile(typ)) + + for_screen_reader(typ) + + "" + ) + + +def render_block_math(match: re.Match[str]) -> str: + src = html.unescape(match.group(1)).removeprefix( + R"\[").removesuffix(R"\]").strip() + typ = f"$ {src} $" + return ( + '
' + + fix_svg(typst_compile(typ)) + + for_screen_reader(typ) + + "
" + ) + + +def for_screen_reader(typ: str) -> str: + return f'{html.escape(typ)}' + + +def fix_svg(svg: bytes) -> str: + """Fix the compiled SVG to be embedded in HTML + + - Strip trailing spaces + - Support dark theme + """ + return re.sub( + r' (fill|stroke)="#000000"', + r' \1="var(--md-typeset-color)"', + svg.decode().strip(), + ) + + +@cache +def typst_compile( + typ: str, + *, + prelude="#set page(width: auto, height: auto, margin: 0pt, fill: none)\n", + format="svg", +) -> bytes: + """Compile a Typst document + + https://github.com/marimo-team/marimo/discussions/2441 + """ + try: + return run( + ["typst", "compile", "-", "-", "--format", format], + input=(prelude + typ).encode(), + check=True, + capture_output=True, + ).stdout + except CalledProcessError as err: + raise RuntimeError( + f""" +Failed to render a typst math: + +```typst +{typ} +``` + +{err.stderr.decode()} +""".strip() + ) diff --git a/py-feos/docs/index.md b/py-feos/docs/index.md new file mode 100644 index 000000000..4426b27e3 --- /dev/null +++ b/py-feos/docs/index.md @@ -0,0 +1,102 @@ +# FeOs Python Documentation + +Welcome to the documentation for the **FeOs Python package** - Python bindings for the FeOs framework for equations of state and classical density functional theory. + +## What is FeOs? + +FeOs is a comprehensive framework that provides Rust implementations of different equation of state and Helmholtz energy functional models with corresponding Python bindings. The Python package offers a user-friendly interface to perform: + +- **Thermodynamic property calculations** for pure components and mixtures +- **Phase equilibrium calculations** including critical points and phase diagrams +- **Classical density functional theory** calculations for inhomogeneous systems +- **Parameter estimation** for equation of state models + +## Key Features + +- **Multiple EoS models**: PC-SAFT, EPC-SAFT, SAFT-VR Mie, PETS, UV Theory, and more +- **Group contribution methods**: GC-PC-SAFT for parameter prediction +- **DFT calculations**: Adsorption isotherms, surface tension, pair correlation functions +- **Parameter management**: JSON-based parameter databases and estimation tools +- **High performance**: Rust backend with automatic differentiation +- **Type safety**: Comprehensive type hints for excellent IDE support + +## Quick Example + +Get started with a simple thermodynamic calculation: + +```python +import feos + +# Define PC-SAFT parameters for methanol +record = feos.PureRecord( + feos.Identifier(name="methanol"), + molarweight=32.04, + m=1.5255, + sigma=3.23, + epsilon_k=188.9, + kappa_ab=0.035176, + epsilon_k_ab=2899.5, + na=1, + nb=1, +) + +# Create equation of state +parameters = feos.Parameters.new_pure(record) +eos = feos.EquationOfState.pcsaft(parameters) + +# Calculate critical point +critical_point = feos.State.critical_point_pure(eos) +print(f"Critical temperature: {critical_point.temperature}") +print(f"Critical pressure: {critical_point.pressure()}") +``` + +## Available Models + +The following models are currently implemented: + +| Model | Description | EoS | DFT | +|-------|-------------|-----|-----| +| `pcsaft` | Perturbed-chain statistical associating fluid theory | ✓ | ✓ | +| `epcsaft` | Electrolyte PC-SAFT | ✓ | | +| `gc-pcsaft` | Group contribution PC-SAFT | ✓ | ✓ | +| `pets` | Perturbed truncated and shifted Lennard-Jones | ✓ | ✓ | +| `uvtheory` | Equation of state for Mie fluids | ✓ | | +| `saftvrqmie` | SAFT-VR for quantum fluids | ✓ | ✓ | +| `saftvrmie` | SAFT-VR for variable range Mie interactions | ✓ | | + +## Getting Started + +Ready to dive in? Here's how to get started: + +1. **[Installation](getting-started/installation.md)** - Install the package via pip or build from source +2. **[Quick Start](getting-started/quickstart.md)** - Your first calculations with FeOs +3. **[Core Concepts](getting-started/concepts.md)** - Understand the main components + +## Documentation Structure + +- **[Getting Started](getting-started/installation.md)** - Installation and basic usage +- **[User Guide](user-guide/working-with-states.md)** - Detailed guides for different functionalities +- **[Examples](examples/basic-properties.md)** - Practical code examples for common tasks +- **[API Reference](api/feos.md)** - Complete API documentation + +## Need Help? + +- Browse the **[Examples](examples/basic-properties.md)** for code you can copy and adapt +- Check the **[API Reference](api/feos.md)** for detailed method documentation +- Visit the [main FeOs repository](https://github.com/feos-org/feos) for issues and discussions + +## Citation + +If you use FeOs in your research, please cite our publication: + +```bibtex +@article{rehner2023feos, + author = {Rehner, Philipp and Bauer, Gernot and Gross, Joachim}, + title = {FeOs: An Open-Source Framework for Equations of State and Classical Density Functional Theory}, + journal = {Industrial \& Engineering Chemistry Research}, + volume = {62}, + number = {12}, + pages = {5347-5357}, + year = {2023}, +} +``` diff --git a/py-feos/docs/javascript/mathjax.js b/py-feos/docs/javascript/mathjax.js new file mode 100644 index 000000000..819e7df5f --- /dev/null +++ b/py-feos/docs/javascript/mathjax.js @@ -0,0 +1,20 @@ +window.MathJax = { + tex: { + inlineMath: [["\\(", "\\)"]], + displayMath: [["\\[", "\\]"]], + processEscapes: true, + processEnvironments: true, + tags: 'ams' + }, + options: { + ignoreHtmlClass: ".*|", + processHtmlClass: "arithmatex" + } + }; + + document$.subscribe(() => { + MathJax.startup.output.clearCache() + MathJax.typesetClear() + MathJax.texReset() + MathJax.typesetPromise() + }) \ No newline at end of file diff --git a/py-feos/docs/theory/dft/derivatives.md b/py-feos/docs/theory/dft/derivatives.md new file mode 100644 index 000000000..c8e921625 --- /dev/null +++ b/py-feos/docs/theory/dft/derivatives.md @@ -0,0 +1,81 @@ +# Derivatives of density profiles + +For converged density profiles equilibrium properties can be calculated as partial derivatives of thermodynamic potentials analogous to classical (bulk) thermodynamics. The difference is that the derivatives have to be along a path of valid density profiles (solutions of the [Euler-Lagrange equation](euler_lagrange_equation.md)). + +The density profiles are calculated implicitly from the Euler-Lagrange equation, which can be written simplified as + +$$\Omega_{\rho_i}(T,\lbrace\mu_k\rbrace,[\lbrace\rho_k(\mathbf{r})\rbrace])=F_{\rho_i}(T,[\lbrace\rho_k(\mathbf{r})\rbrace])-\mu_i+V_i^\mathrm{ext}(\mathbf{r})=0$$ (eqn:euler_lagrange) + +Incorporating bond integrals can be done similar to the section on the [Newton solver](solver.md) but will not be discussed in this section. The derivatives of the density profiles can then be calculated from the total differential of eq. {eq}`eqn:euler_lagrange`, leading to + +$$\mathrm{d}\Omega_{\rho_i}(\mathbf{r})=\left(\frac{\partial\Omega_{\rho_i}(\mathbf{r})}{\partial T}\right)_{\mu_k,\rho_k}\mathrm{d}T+\sum_j\left(\frac{\partial\Omega_{\rho_i}(\mathbf{r})}{\partial\mu_j}\right)_{T,\mu_k,\rho_k}\mathrm{d}\mu_j+\int\sum_j\left(\frac{\delta\Omega_{\rho_i}(\mathbf{r})}{\delta\rho_j(\mathbf{r}')}\right)_{T,\mu_k,\rho_k}\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=0$$ + +Using eq. {eq}`eqn:euler_lagrange` and the shortened notation for derivatives of functionals in their natural variables, e.g., $F_T=\left(\frac{\partial F}{\partial T}\right)_{\rho_k}$, the expression can be simplified to + +$$F_{T\rho_i}(\mathbf{r})\mathrm{d}T-\mathrm{d}\mu_i+\int\sum_j F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=0$$ (eqn:gibbs_duhem) + +Similar to the Gibbs-Duhem relation for bulk phases, eq. {eq}`eqn:gibbs_duhem` shows how temperature, chemical potentials and the density profiles in an inhomogeneous system cannot be varied independently. The derivatives of the density profiles with respect to the intensive variables can be directly identified as + +$$\int\sum_j F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial T}\right)_{\mu_k}\mathrm{d}\mathbf{r}'=-F_{T\rho_i}(\mathbf{r})$$ + +and + +$$\int\sum_j F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial\mu_k}\right)_{T}\mathrm{d}\mathbf{r}'=\delta_{ik}$$ (eqn:drho_dmu) + +Both of these expressions are implicit (linear) equations for the derivatives. They can be solved rapidly analogously to the implicit expression appearing in the [Newton solver](solver.md). In practice, it is useful to explicitly cancel out the (often unknown) thermal de Broglie wavelength $\Lambda_i$ from the expression where it has no influence. This is done by splitting the intrinsic Helmholtz energy into an ideal gas and a residual part. + +$$F=k_\mathrm{B}T\int\sum_im_i\rho_i(\mathbf{r})\left(\ln\left(\rho_i(\mathbf{r})\Lambda_i^3\right)-1\right)\mathrm{d}\mathbf{r}+\mathcal{\hat F}^\mathrm{res}$$ + +Then $F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')=m_i\frac{k_\mathrm{B}T}{\rho_i(\mathbf{r})}\delta_{ij}\delta(\mathbf{r}-\mathbf{r}')+\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')$ and eq. {eq}`eqn:drho_dmu` can be rewritten as + +$$m_i\frac{k_\mathrm{B}T}{\rho_i(\mathbf{r})}\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\mu_k}\right)_T+\int\sum_j\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial\mu_k}\right)_{T}\mathrm{d}\mathbf{r}'=\delta_{ik}$$ + +In practice, the division by the density should be avoided for numerical reasons and the energetic properties are reduced with the factor $\beta=\frac{1}{k_\mathrm{B}T}$. The final expression is + +$$m_i\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\beta\mu_k}\right)_T+\rho_i(\mathbf{r})\int\sum_j\beta\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial\beta\mu_k}\right)_{T}\mathrm{d}\mathbf{r}'=\rho_i(\mathbf{r})\delta_{ik}$$ + +For the temperature derivative, it is more convenient to express eq. {eq}`eqn:gibbs_duhem` in terms of the pressure of a bulk phase that is in equilibrium with the inhomogeneous system. In the following, only paths along **constant bulk composition** are considered. With this constraint, the total differential of the chemical potential simplifies to + +$$\mathrm{d}\mu_i=-s_i\mathrm{d}T+v_i\mathrm{d}p$$ + +which can be used in eq. {eq}`eqn:gibbs_duhem` to give + +$$\left(F_{T\rho_i}(\mathbf{r})+s_i\right)\mathrm{d}T+\int\sum_j F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=v_i\mathrm{d}p$$ + +Even though $s_i$ is readily available in $\text{FeO}_\text{s}$ it is useful at this point to rewrite the partial molar entropy as + +$$s_i=v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}-F_{T\rho_i}^\mathrm{b}$$ + +Then, the intrinsic Helmholtz energy can be split into an ideal gas and a residual part again, and the de Broglie wavelength cancels. + +$$ +\begin{align*} +&\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+F_{T\rho_i}^\mathrm{res}(\mathbf{r})-F_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right)\mathrm{d}T\\ +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~+m_i\frac{k_\mathrm{B}T}{\rho_i(\mathbf{r})}\delta\rho_i(\mathbf{r})+\int\sum_j\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=v_i\mathrm{d}p +\end{align*} +$$ + +Finally the expressions for the derivatives with respect to pressure + +$$m_i\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\beta p}\right)_{T,x_k}+\rho_i(\mathbf{r})\int\sum_j\beta\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial\beta p}\right)_{T,x_k}\mathrm{d}\mathbf{r}'=v_i\rho_i(\mathbf{r})$$ + +and temperature + +$$ +\begin{align*} +&m_i\left(\frac{\partial\rho_i(\mathbf{r})}{\partial T}\right)_{p,x_k}+\rho_i(\mathbf{r})\int\sum_j\beta\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial T}\right)_{p,x_k}\mathrm{d}\mathbf{r}'\\ +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=-\frac{\rho_i(\mathbf{r})}{k_\mathrm{B}T}\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+F_{T\rho_i}^\mathrm{res}(\mathbf{r})-F_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right) +\end{align*} +$$ + +follow. All derivatives $x_i$ shown here can be calculated from the same linear equation + +$$m_ix_i+\rho_i(\mathbf{r})\int\sum_j\beta\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')x_i\mathrm{d}\mathbf{r}'=y_i$$ + +by just replacing the right hand side $y_i$. + +| derivative | right hand side | +| ------------------------------------------------------------------------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\beta\mu_k}\right)_T$ | $\rho_i(\mathbf{r})\delta_{ik}$ | +| $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\beta p}\right)_{T,x_k}$ | $\rho_i(\mathbf{r})v_i$ | +| $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial T}\right)_{p,x_k}$ | $-\frac{\rho_i(\mathbf{r})}{k_\mathrm{B}T}\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+F_{T\rho_i}^\mathrm{res}(\mathbf{r})-F_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right)$ | diff --git a/py-feos/docs/theory/dft/enthalpy_of_adsorption.md b/py-feos/docs/theory/dft/enthalpy_of_adsorption.md new file mode 100644 index 000000000..5ef18d69c --- /dev/null +++ b/py-feos/docs/theory/dft/enthalpy_of_adsorption.md @@ -0,0 +1,113 @@ +# Enthalpy of adsorption and the Clausius-Clapeyron relation + +## Enthalpy of adsorption + +The energy balance in differential form for a simple adsorption process can be written as + +$$\mathrm{d}U=h^\mathrm{in}\delta n^\mathrm{in}-h^\mathrm{b}\delta n^\mathrm{out}+\delta Q$$ (eqn:energy_balance) + +Here the balance is chosen to only include the fluid in the porous medium. The molar enthalpy $h^\mathrm{b}$ of the (bulk) fluid that leaves the adsorber is at a state that is in equilibrium with the porous medium. In contrast, the incoming stream can be at any condition. Analogously, the component balance is + +$$\mathrm{d}N_i=x_i^\mathrm{in}\delta n^\mathrm{in}-x_i\delta n^\mathrm{out}$$ (eqn:mass_balance) + +The differential of the internal energy can be replaced with the total differential in its variables temperature $T$ and number of particles $N_i$. The volume of the adsorber is fixed and thus not considered as a variable. + +$$\mathrm{d}U=\left(\frac{\partial U}{\partial T}\right)_{N_k}\mathrm{d}T+\sum_i\left(\frac{\partial U}{\partial N_i}\right)_{T,N_j}\mathrm{d}N_i$$ (eqn:U_differential) + +Eqs. {eq}`eqn:energy_balance`, {eq}`eqn:mass_balance` and {eq}`eqn:U_differential` can be combined into an expression for the heat of adsorption $\delta Q$ + +$$\delta Q=\left(\frac{\partial U}{\partial T}\right)_{N_k}\mathrm{d}T-\left(h^\mathrm{in}-\sum_ix_i^\mathrm{in}\left(\frac{\partial U}{\partial N_i}\right)_{T,N_j}\right)\delta n^\mathrm{in}+\left(h^\mathrm{b}-\sum_ix_i\left(\frac{\partial U}{\partial N_i}\right)_{T,N_j}\right)\delta n^\mathrm{out}$$ + +The heat of adsorption can thus be split into a sensible part that depends on the change in temperature, and a latent part that depends on the change in loading. The expression can be simplified by using the definitions of the isochoric heat capacity $C_v=\left(\frac{\partial U}{\partial T}\right)_{N_k}$ and the **partial molar enthalpy of adsorption** + +$$\Delta h_i^\mathrm{ads}=h_i^\mathrm{b}-\left(\frac{\partial U}{\partial N_i}\right)_{T,N_j}$$ + +yielding + +$$\delta Q=C_v\mathrm{d}T-\sum_ix_i^\mathrm{in}\left(h_i^\mathrm{in}-h_i^\mathrm{b}+\Delta h_i^\mathrm{ads}\right)\delta n^\mathrm{in}+\sum_ix_i\Delta h_i^\mathrm{ads}\delta n^\mathrm{out}$$ + +or + +$$\delta Q=C_v\mathrm{d}T-\sum_ix_i^\mathrm{in}\left(h_i^\mathrm{in}-h_i^\mathrm{b}+\Delta h_i^\mathrm{ads}\right)\delta n^\mathrm{in}+\Delta h^\mathrm{ads}\delta n^\mathrm{out}$$ + +with the **enthalpy of adsorption** + +$$\Delta h^\mathrm{ads}=\sum_ix_i\Delta h_i^\mathrm{ads}=h^\mathrm{b}-\sum_ix_i\left(\frac{\partial U}{\partial N_i}\right)_{T,N_j}$$ + +For **pure components** the balance equations simplify to + +$$\delta Q=C_v\mathrm{d}T-\left(h^\mathrm{in}-h^\mathrm{b}\right)\delta n^\mathrm{in}-\Delta h^\mathrm{ads}\mathrm{d}N$$ + +## Clausius-Clapeyron relation for porous media + +The Clausius-Clapeyron relation relates the $p-T$ slope of a pure component phase transition line to the corresponding enthalpy of phase change. For a vapor-liquid phase transition, the exact relation is + +$$\frac{\mathrm{d}p^\mathrm{sat}}{\mathrm{d}T}=\frac{s^\mathrm{V}-s^\mathrm{L}}{v^\mathrm{V}-v^\mathrm{L}}=\frac{h^\mathrm{V}-h^\mathrm{L}}{T\left(v^\mathrm{V}-v^\mathrm{L}\right)}$$ (eqn:temp_dep_press) + +In this expression, the enthalpy of vaporization $\Delta h^\mathrm{vap}=h^\mathrm{V}-h^\mathrm{L}$ can be identified. The molar volumes $v$ of the two phases can be replaced by the compressibility factor $Z=\frac{pv}{RT}$. Then, eq. {eq}`eqn:temp_dep_press` simplifies to + +$$\frac{\mathrm{d}p^\mathrm{sat}}{\mathrm{d}T}=\frac{p}{R T^2}\frac{\Delta h^\mathrm{vap}}{Z^\mathrm{V}-Z^\mathrm{L}}$$ + +which can be compactly written as + +$$\frac{\mathrm{d}\ln p^\mathrm{sat}}{\mathrm{d}\frac{1}{RT}}=-\frac{\Delta h^\mathrm{vap}}{Z^\mathrm{V}-Z^\mathrm{L}}$$ (eqn:Clausius_Clapeyron_exact) + +Eq. {eq}`eqn:Clausius_Clapeyron_exact` is still an exact expression. In practice, the volume (and hence the compressibility factor) of the liquid phase can often be neglected compared to the volume of the gas phase. Additionally assuming an ideal gas phase ($Z^\mathrm{V}\approx1$), leads to the expression commonly referred to as Clausius-Clapeyron relation: + +$$\frac{\mathrm{d}\ln p^\mathrm{sat}}{\mathrm{d}\frac{1}{RT}}=-\Delta h^\mathrm{vap}$$ (eqn:Clausius_Clapeyron) + +A similar relation can be derived for fluids adsorbed in a porous medium that is in equilibrium with a bulk phase. At this point it is important to clarify which variables describe the system + +- The adsorbed fluid and the bulk phase are in equilibrium. Therefore, the temperature $T$ and chemical potentials $\mu_i$ are the same for both phases. +- The density profiles and hence the number of particles $N_i$ in the porous medium is determined by $T$ and $\mu_i$. The volume of the porous medium is not considered as a thermodynamic variable but rather as a (constant) property of the adsorbent. +- All intensive properties of the bulk phase are fully determined by $T$ and $\mu_i$. In practice it can be useful to relate these properties to measurable properties like the pressure $p$ and the composition $x_i$. + +To find an expression of the slope of an isostere (constant $N_i$), the pressure, which is only defined for the bulk phase, has to be related to properties of the adsorbed fluid. + +$$\frac{\mathrm{d}\ln p}{\mathrm{d}\frac{1}{RT}}=-\frac{RT^2}{p}\frac{\mathrm{d}p}{\mathrm{d}T}$$ + +First, the pressure can be replaced using the Gibbs-Duhem relation for the bulk phase (index $\mathrm{b}$) + +$$\frac{\mathrm{d}\ln p}{\mathrm{d}\frac{1}{RT}}=-\frac{RT^2}{pv^\mathrm{b}}\left(s^\mathrm{b}+\sum_ix_i\left(\frac{\partial\mu_i}{\partial T}\right)_{N_k}\right)$$ (eqn:clausius_clapeyron_intermediate) + +Here the directional derivative $\frac{\mathrm{d}\mu_i}{\mathrm{d}T}$ could be replaced with a partial derivative amongst the variables describing the adsorbed fluid. The partial derivative can then be replaced using a Maxwell relation based on the Helmholtz energy $F$ as follows + +$$\left(\frac{\partial\mu_i}{\partial T}\right)_{N_k}=\left(\frac{\partial^2 F}{\partial T\partial N_i}\right)=-\left(\frac{\partial S}{\partial N_i}\right)_{T,N_j}$$ + +Using the Maxwell relation together with the compressibility factor of the bulk phase $Z^\mathrm{b}=\frac{pv^\mathrm{b}}{RT}$ in eq. {eq}`eqn:clausius_clapeyron_intermediate` results in + +$$\frac{\mathrm{d}\ln p}{\mathrm{d}\frac{1}{RT}}=-\frac{T}{Z^\mathrm{b}}\left(s^\mathrm{b}-\sum_ix_i\left(\frac{\partial S}{\partial N_i}\right)_{T,N_j}\right)$$ + +Finally, using $h^\mathrm{b}=Ts^\mathrm{b}+\sum_ix_i\mu_i$ and $\mathrm{d}U=T\mathrm{d}S+\sum_i\mu_i\mathrm{d}N_i$ leads to + +$$\frac{\mathrm{d}\ln p}{\mathrm{d}\frac{1}{RT}}=-\frac{1}{Z^\mathrm{b}}\left(h^\mathrm{b}-\sum_ix_i\left(\frac{\partial U}{\partial N_i}\right)_{T,N_j}\right)=-\frac{\Delta h^\mathrm{ads}}{Z^\mathrm{b}}$$ (eqn:deriv_relation_hads) + +The relation is exact and valid for an arbitrary number of components in the fluid phase. + +## Calculation of the enthalpy of adsorption from classical DFT + +In a DFT context, the introduction of entropies and internal energies are just unnecessary complications. The most useful definition of the (partial molar) enthalpy of adsorption is + +$$\Delta h_i^\mathrm{ads}=T\left(s_i^\mathrm{b}+\left(\frac{\partial\mu_i}{\partial T}\right)_{N_k}\right)$$ + +The derivative at constant number of particles is problematic and has to be replaced. This is done starting from the total differential of the number of particles + +$$\mathrm{d}N_i=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\mathrm{d}\mu_j+\left(\frac{\partial N_i}{\partial T}\right)_{\mu_k}\mathrm{d}T$$ (eqn:dn) + +Calculating the derivative with respect to $T$ at constant $N_i$ leads to + +$$0=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\left(\frac{\partial\mu_j}{\partial T}\right)_{N_k}+\left(\frac{\partial N_i}{\partial T}\right)_{\mu_k}$$ (eqn:dndt_1) + +from which the unknown derivative $\left(\frac{\partial\mu_i}{\partial T}\right)_{N_k}$ can be calculated. In practice the expression has the disadvantage that $\left(\frac{\partial N_i}{\partial T}\right)_{\mu_k}$ depends on the (sometimes unknown) thermal de Broglie wavelength which cancels later with $s_i^\mathrm{b}$. This can be remedied by first calculating the derivative of eq. {eq}`eqn:dn` with respect to $T$ at constant (bulk) pressure and composition. + +$$\left(\frac{\partial N_i}{\partial T}\right)_{p,x_k}=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\left(\frac{\partial\mu_j}{\partial T}\right)_{p,x_k}+\left(\frac{\partial N_i}{\partial T}\right)_{\mu_k}$$ (eqn:dndt_2) + +From classical bulk thermodynamics we know $\left(\frac{\partial\mu_j}{\partial T}\right)_{p,x_k}=-s_j^\mathrm{b}$ and therefore, eq. {eq}`eqn:dndt_2` can be used in eq. {eq}`eqn:dndt_1` to give + +$$0=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\left(s_j^\mathrm{b}+\left(\frac{\partial\mu_j}{\partial T}\right)_{N_k}\right)+\left(\frac{\partial N_i}{\partial T}\right)_{p,x_k}$$ + +After multiplying with $T$, the following elegant expression remains + +$$0=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\Delta h_j^\mathrm{ads}+T\left(\frac{\partial N_i}{\partial T}\right)_{p,x_k}$$ + +which is a symmetric linear system of equations due to $\left(\frac{\partial N_i}{\partial\mu_j}\right)_T=-\left(\frac{\partial^2\Omega}{\partial\mu_i\partial\mu_j}\right)_T$. The derivatives of the particle numbers are obtained by integrating over the respective derivatives of the density profiles which were discussed [previously](derivatives.md). diff --git a/py-feos/docs/theory/dft/euler_lagrange_equation.md b/py-feos/docs/theory/dft/euler_lagrange_equation.md new file mode 100644 index 000000000..b5b9b39fc --- /dev/null +++ b/py-feos/docs/theory/dft/euler_lagrange_equation.md @@ -0,0 +1,121 @@ +# Euler-Lagrange equation + +The fundamental expression in classical density functional theory is the relation between the grand potential functional $\Omega$ and the intrinsic Helmholtz energy functional $F$. + +$$\Omega(T,\mu,[\rho(r)])=F(T,[\rho(r)])-\sum_i\int\rho_i(r)\left(\mu_i-V_i^\mathrm{ext}(r)\right)\mathrm{d}r$$ + +What makes this expression so appealing is that the intrinsic Helmholtz energy functional only depends on the temperature $T$ and the density profiles $\rho_i(r)$ of the system and not on the external potentials $V_i^\mathrm{ext}(r)$. + +For a given temperature $T$, chemical potentials $\mu_i$ and external potentials $V_i^\mathrm{ext}(r)$ the grand potential reaches a minimum at equilibrium. Mathematically this condition can be written as + +\begin{equation} +\left.\frac{\delta\Omega}{\delta\rho*i(r)}\right|*{T,\mu}=F\_{\rho_i}(r)-\mu_i+V_i^{\mathrm{ext}}(r)=0 +\label{eqn:euler_lagrange_mu} +\end{equation} + +where $F_{\rho_i}(r)=\left.\frac{\delta F}{\delta\rho_i(r)}\right|_T$ is short for the functional derivative of the intrinsic Helmholtz energy. In this context, eq. $\eqref{eqn:euler_lagrange_mu}$ is commonly referred to as the Euler-Lagrange equation, an implicit nonlinear integral equation which needs to be solved for the equilibrium density profiles of the system. + +For a homogeneous (bulk) system, $V_i^\mathrm{ext}=0$ and we get + +$$F_{\rho_i}^\mathrm{b}-\mu_i=0$$ (eqn:euler_lagrange_bulk) + +The functional derivative of the Helmholtz energy of a bulk system $F_{\rho_i}^\mathrm{b}$ is a function of the temperature $T$ and bulk densities $\rho^\mathrm{b}$. At this point, it can be advantageous to relate the grand potential of an inhomogeneous system to the densities of a bulk system that is in equilibrium with the inhomogeneous system. This approach has several advantages: + +- The thermal de Broglie wavelength $\Lambda$ cancels out. +- If the chemical potential of the system is not known, all variables are the same quantity (densities). +- The bulk system is described by a Helmholtz energy which is explicit in the density, so there are no internal iterations required. + +Using eq. {eq}`eqn:euler_lagrange_bulk` in eq. {eq}`eqn:euler_lagrange_mu` leads to the Euler-Lagrange equation + +$$\left.\frac{\delta\Omega}{\delta\rho_i(r)}\right|_{T,\rho^\mathrm{b}}=F_{\rho_i}(r)-F_{\rho_i}^\mathrm{b}+V_i^{\mathrm{ext}}(r)=0$$ (eqn:euler_lagrange_rho) + +## Spherical molecules + +In the simplest case, the molecules under consideration can be described as spherical. Then the Helmholtz energy can be split into an ideal and a residual part: + +$$\beta F=\sum_i\int\rho_i(r)\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\right)\mathrm{d}r+\beta F^\mathrm{res}$$ + +with the thermal de Broglie wavelength $\Lambda_i$. The functional derivatives for an inhomogeneous and a bulk system follow as + +$$\beta F_{\rho_i}(r)=\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta F_{\rho_i}^\mathrm{res}$$ + +$$\beta F_{\rho_i}^\mathrm{b}=\ln\left(\rho_i^\mathrm{b}\Lambda_i^3\right)+\beta F_{\rho_i}^\mathrm{b,res}$$ + +Using these expressions in eq. {eq}`eqn:euler_lagrange_rho` results in + +$$\left.\frac{\delta\beta\Omega}{\delta\rho_i(r)}\right|_{T,\rho^\mathrm{b}}=\ln\left(\frac{\rho_i(r)}{\rho_i^\mathrm{b}}\right)+\beta\left(F_{\rho_i}^\mathrm{res}(r)-F_{\rho_i}^\mathrm{b,res}+V_i^{\mathrm{ext}}(r)\right)=0$$ + +The Euler-Lagrange equation can be recast as + +$$\rho_i(r)=\rho_i^\mathrm{b}e^{\beta\left(F_{\rho_i}^\mathrm{b,res}-F_{\rho_i}^\mathrm{res}(r)-V_i^\mathrm{ext}(r)\right)}$$ + +which is convenient because it leads directly to a recurrence relation known as Picard iteration. + +## Homosegmented chains + +For chain molecules that do not resolve individual segments (essentially the PC-SAFT Helmholtz energy functional) a chain contribution is introduced as + +$$\beta F^\mathrm{chain}=-\sum_i\int\rho_i(r)\left(m_i-1\right)\ln\left(\frac{y_{ii}\lambda_i(r)}{\rho_i(r)}\right)\mathrm{d}r$$ + +Where $m_i$ is the number of segments (i.e., the PC-SAFT chain length parameter), $y_{ii}$ is the cavity correlation function at contact in the reference fluid, and $\lambda_i$ is a weighted density. +The presence of $\rho_i(r)$ in the logarithm poses numerical problems. Therefore, it is convenient to rearrange the expression as + +$$ +\begin{align} +\beta F^\mathrm{chain}=&\sum_i\int\rho_i(r)\left(m_i-1\right)\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\right)\mathrm{d}r\\ +&\underbrace{-\sum_i\int\rho_i(r)\left(m_i-1\right)\left(\ln\left(y_{ii}\lambda_i(r)\Lambda_i^3\right)-1\right)\mathrm{d}r}_{\beta\hat{F}^\mathrm{chain}} +\end{align} +$$ + +Then the total Helmholtz energy + +$$\beta F=\sum_i\int\rho_i(r)\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\right)\mathrm{d}r+\beta F^\mathrm{chain}+\beta F^\mathrm{res}$$ + +can be rearranged to + +$$\beta F=\sum_i\int\rho_i(r)m_i\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\right)\mathrm{d}r+\underbrace{\beta\hat{F}^\mathrm{chain}+\beta F^\mathrm{res}}_{\beta\hat{F}^\mathrm{res}}$$ + +The functional derivatives are then similar to the spherical case + +$$\beta F_{\rho_i}(r)=m_i\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta\hat{F}_{\rho_i}^\mathrm{res}(r)$$ + +$$\beta F_{\rho_i}^\mathrm{b}=m_i\ln\left(\rho_i^\mathrm{b}\Lambda_i^3\right)+\beta\hat{F}_{\rho_i}^\mathrm{b,res}$$ + +and lead to a slightly modified Euler-Lagrange equation + +$$\rho_i(r)=\rho_i^\mathrm{b}e^{\frac{\beta}{m_i}\left(\hat F_{\rho_i}^\mathrm{b,res}-\hat F_{\rho_i}^\mathrm{res}(r)-V_i^\mathrm{ext}(r)\right)}$$ + +## Heterosegmented chains + +The expressions are more complex for models in which density profiles of individual segments are considered. A derivation is given in the appendix of [Rehner et al. (2022)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.105.034110). The resulting Euler-Lagrange equation is given as + +$$\rho_\alpha(r)=\Lambda_i^{-3}e^{\beta\left(\mu_i-\hat F_{\rho_\alpha}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)$$ + +with the bond integrals $I_{\alpha\alpha'}(r)$ that are calculated recursively from + +$$I_{\alpha\alpha'}(r)=\int e^{-\beta\left(\hat{F}_{\rho_{\alpha'}}(r')+V_{\alpha'}^\mathrm{ext}(r')\right)}\left(\prod_{\alpha''\neq\alpha}I_{\alpha'\alpha''}(r)\right)\omega_\mathrm{chain}^{\alpha\alpha'}(r-r')\mathrm{d}r$$ + +The index $\alpha$ is used for every segment on component $i$, $\alpha'$ refers to all segments bonded to segment $\alpha$ and $\alpha''$ to all segments bonded to $\alpha'$. +For bulk systems the expressions simplify to + +$$\rho_\alpha^\mathrm{b}=\Lambda_i^{-3}e^{\beta\left(\mu_i-\sum_\gamma\hat F_{\rho_\gamma}^\mathrm{b,res}\right)}$$ + +which shows that, by construction, the density of every segment in a molecule is identical in a bulk system. The index $\gamma$ refers to all segments on moecule $i$. The expressions can be combined in a similar way as for the molecular DFT: + +$$\rho_\alpha(r)=\rho_\alpha^\mathrm{b}e^{\beta\left(\sum_\gamma\hat F_{\rho_\gamma}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)$$ + +At this point it can be numerically useful to redistribute the bulk contributions back into the bond integrals + +$$\rho_\alpha(r)=\rho_\alpha^\mathrm{b}e^{\beta\left(\hat F_{\rho_\alpha}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)$$ + +$$I_{\alpha\alpha'}(r)=\int e^{\beta\left(\hat F_{\rho_{\alpha'}}^\mathrm{b,res}-\hat F_{\rho_{\alpha'}}^\mathrm{res}(r')-V_{\alpha'}^\mathrm{ext}(r')\right)}\left(\prod_{\alpha''\neq\alpha}I_{\alpha'\alpha''}(r)\right)\omega_\mathrm{chain}^{\alpha\alpha'}(r-r')\mathrm{d}r$$ + +## Combined expression + +To avoid having multiple implementations of the central part of the DFT code, the different descriptions of molecules can be combined in a single version of the Euler-Lagrange equation: + +$$\rho_\alpha(r)=\rho_\alpha^\mathrm{b}e^{\frac{\beta}{m_\alpha}\left(\hat F_{\rho_\alpha}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)$$ + +$$I_{\alpha\alpha'}(r)=\int e^{\frac{\beta}{m_{\alpha'}}\left(\hat F_{\rho_{\alpha'}}^\mathrm{b,res}-\hat F_{\rho_{\alpha'}}^\mathrm{res}(r')-V_{\alpha'}^\mathrm{ext}(r')\right)}\left(\prod_{\alpha''\neq\alpha}I_{\alpha'\alpha''}(r')\right)\omega_\mathrm{chain}^{\alpha\alpha'}(r-r')\mathrm{d}r'$$ + +If molecules consist of single (possibly non-spherical) segments, the Euler-Lagrange equation simplifies to that of the homosegmented chains shown above. For heterosegmented chains, the correct expression is obtained by setting $m_\alpha=1$. diff --git a/py-feos/docs/theory/dft/functional_derivatives.md b/py-feos/docs/theory/dft/functional_derivatives.md new file mode 100644 index 000000000..eb054b1f9 --- /dev/null +++ b/py-feos/docs/theory/dft/functional_derivatives.md @@ -0,0 +1,44 @@ +# Functional derivatives + +In the last section the functional derivative + +$$\hat F_{\rho_\alpha}^\mathrm{res}(r)=\left(\frac{\delta\hat F^\mathrm{res}}{\delta\rho_\alpha(r)}\right)_{T,\rho_{\alpha'\neq\alpha}}$$ + +was introduced as part of the Euler-Lagrange equation. The implementation of these functional derivatives can be a major difficulty during the development of a new Helmholtz energy model. In $\text{FeO}_\text{s}$ it is fully automated. The core assumption is that the residual Helmholtz energy functional $\hat F^\mathrm{res}$ can be written as a sum of contributions that each can be written in the following way: + +$$F=\int f[\rho(r)]\mathrm{d}r=\int f(\lbrace n_\gamma\rbrace)\mathrm{d}r$$ + +The Helmholtz energy density $f$ which would in general be a functional of the density itself can be expressed as a _function_ of weighted densities $n_\gamma$ which are obtained by convolving the density profiles with weight functions $\omega_\gamma^\alpha$ + +$$n_\gamma(r)=\sum_\alpha\int\rho_\alpha(r')\omega_\gamma^\alpha(r-r')\mathrm{d}r'\tag{1}$$ + +In practice the weight functions tend to have simple shapes like step functions (i.e. the weighted density is an average over a sphere) or Dirac distributions (i.e. the weighted density is an average over the surface of a sphere). + +For Helmholtz energy functionals that can be written in this form, the calculation of the functional derivative can be automated. In general the functional derivative can be written as + +$$F_{\rho_\alpha}(r)=\int\frac{\delta f(r')}{\delta\rho_\alpha(r)}\mathrm{d}r'=\int\sum_\gamma f_{n_\gamma}(r')\frac{\delta n_\gamma(r')}{\delta\rho_\alpha(r)}\mathrm{d}r'$$ + +with $f_{n_\gamma}$ as abbreviation for the _partial_ derivative $\frac{\partial f}{\partial n_\gamma}$. Using the definition of the weighted densities (1), the expression can be rewritten as + +$$ +\begin{align} +F_{\rho_\alpha}(r)&=\int\sum_\gamma f_{n_\gamma}(r')\frac{\delta n_\gamma(r')}{\delta\rho_\alpha(r)}\mathrm{d}r'=\int\sum_\gamma f_{n_\gamma}(r')\sum_{\alpha'}\int\underbrace{\frac{\delta\rho_{\alpha'}(r'')}{\delta\rho_\alpha(r)}}_{\delta(r-r'')\delta_{\alpha\alpha'}}\omega_\gamma^{\alpha'}(r'-r'')\mathrm{d}r''\mathrm{d}r'\\ +&=\sum_\gamma\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r'-r)\mathrm{d}r +\end{align} +$$ + +At this point the parity of the weight functions has to be taken into account. By construction scalar and spherically symmetric weight functions (the standard case) are even functions, i.e., $\omega(-r)=\omega(r)$. In contrast, vector valued weight functions, as they appear in fundamental measure theory, are odd functions, i.e., $\omega(-r)=-\omega(r)$. Therefore, the sum over the weight functions needs to be split into two contributions, as + +$$F_{\rho_\alpha}(r)=\sum_\gamma^\mathrm{scal}\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r-r')\mathrm{d}r-\sum_\gamma^\mathrm{vec}\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r-r')\mathrm{d}r\tag{2}$$ + +With this distinction, the calculation of the functional derivative is split into three steps + +1. Calculate the weighted densities $n_\gamma$ from eq. (1) +2. Evaluate the partial derivatives $f_{n_\gamma}$ +3. Use eq. (2) to obtain the functional derivative $F_{\rho_\alpha}$ + +A fast method to calculate the convolution integrals required in steps 1 and 3 is shown in the next section. + +The implementation of partial derivatives by hand can be cumbersome and error prone. $\text{FeO}_\text{s}$ uses automatic differentiation with dual numbers to facilitate this step. For details about dual numbers and their generalization, we refer to [our publication](https://www.frontiersin.org/articles/10.3389/fceng.2021.758090/full). The essential relation is that the Helmholtz energy density evaluated with a dual number as input for one of the weighted densities evaluates to a dual number with the function value as real part and the partial derivative as dual part. + +$$f(n_\gamma+\varepsilon,\lbrace n_{\gamma'\neq\gamma}\rbrace)=f+f_{n_\gamma}\varepsilon$$ diff --git a/py-feos/docs/theory/dft/ideal_gas.md b/py-feos/docs/theory/dft/ideal_gas.md new file mode 100644 index 000000000..a205c4f95 --- /dev/null +++ b/py-feos/docs/theory/dft/ideal_gas.md @@ -0,0 +1,33 @@ +# Ideal gas properties + +Classical DFT can be used to rapidly determine the ideal gas limit of fluids in porous media. In an ideal gas, there are no interactions between the fluid molecules and therefore the residual Helmholtz energy $F^\mathrm{res}$ and its derivatives vanish. Note that this is only the case for spherical or heterosegmented molecules ($m_\alpha=1$), as the chain contribution in the homosegmented model contains intramolecular interactions. The ideal gas density profile can then be obtained directly from the [Euler-Lagrange equation](euler_lagrange_equation.md): + +$$\rho_\alpha^\mathrm{ig}(\mathbf{r})=\rho_\alpha^\mathrm{b}e^{-\beta V_\alpha^\mathrm{ext}(\mathbf{r})}\prod_{\alpha'}I^\mathrm{ig}_{\alpha\alpha'}(\mathbf{r})$$ (eqn:rho_ideal_gas) + +$$I^\mathrm{ig}_{\alpha\alpha'}(\mathbf{r})=\int e^{-\beta V_{\alpha'}^\mathrm{ext}(\mathbf{r'})}\left(\prod_{\alpha''\neq\alpha}I^\mathrm{ig}_{\alpha'\alpha''}(\mathbf{r'})\right)\omega_\mathrm{chain}^{\alpha\alpha'}(\mathbf{r}-\mathbf{r'})\mathrm{d}\mathbf{r'}$$ + +Either from the derivatives presented [previously](derivatives.md), or from directly calculating derivatives of eq. {eq}`eqn:euler_lagrange_mu`, the **Henry coefficient** $H_\alpha$ can be calculated, as + +$$H_\alpha(T)=\left(\frac{\partial N_\alpha^\mathrm{ig}}{\partial p_\alpha}\right)_{T,x_k}=\int\left(\frac{\partial\rho_\alpha^\mathrm{ig}(\mathbf{r})}{\partial p_\alpha}\right)_{T,x_k}\mathrm{d}\mathbf{r}=\beta\int e^{-\beta V_\alpha^\mathrm{ext}(\mathbf{r})}\prod_{\alpha'}I^\mathrm{ig}_{\alpha\alpha'}(\mathbf{r})\mathrm{d}\mathbf{r}$$ + +By construction the Henry coefficients for all segments $\alpha$ of a molecule $i$ are identical. Therefore, the number of adsorbed molecules in an ideal gas state (the Henry regime) can be calculated from + +$$N_i^\mathrm{ig}=k_\mathrm{B}T\rho_i^\mathrm{b}H_i(T)$$ + +The expression can be used in the general equation for the **enthalpy of adsorption** (see [here](enthalpy_of_adsorption.md)) + +$$0=\sum_j\left(\frac{\partial N_i^\mathrm{ig}}{\partial\mu_j}\right)_T\Delta h_j^\mathrm{ads,ig}+T\left(\frac{\partial N_i^\mathrm{ig}}{\partial T}\right)_{p,x_k}$$ + +to simplify to + +$$0=\rho_i^\mathrm{b}H_i(T)\Delta h_i^\mathrm{ads,ig}+k_\mathrm{B}T^2\rho_i^\mathrm{b}H_i'(T)$$ + +and finally + +$$\Delta h_i^\mathrm{ads,ig}=-k_\mathrm{B}T^2\frac{H_i'(T)}{H_i(T)}$$ + +For a spherical molecule without bond integrals, the derivative can be evaluated straightforwardly to yield + +$$\Delta h_i^\mathrm{ads,ig}=\frac{\int\left(k_\mathrm{B}T-V_i^\mathrm{ext}(\mathbf{r})\right)e^{-\beta V_i^\mathrm{ext}(\mathbf{r})}\mathrm{d}\mathbf{r}}{\int e^{-\beta V_i^\mathrm{ext}(\mathbf{r})}\mathrm{d}\mathbf{r}}$$ + +Analytical derivatives of the bond integrals can be determined, however, in $\text{FeO}_\text{s}$ automatic differentiation with dual numbers is used to obtain correct derivatives with barely any implementation overhead. diff --git a/py-feos/docs/theory/dft/pdgt.md b/py-feos/docs/theory/dft/pdgt.md new file mode 100644 index 000000000..a6d09d172 --- /dev/null +++ b/py-feos/docs/theory/dft/pdgt.md @@ -0,0 +1,52 @@ +# Predictive density gradient theory + +Predictive density gradient theory (pDGT) is an efficient approach for the prediction of surface tensions, which is derived from non-local DFT, see [Rehner et al. (2018)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.98.063312). A gradient expansion is applied to the weighted densities of the Helmholtz energy functional to second order as well as to the Helmholtz energy density to first order. + +Weighted densities (in non-local DFT) are determined from + +$$ n*\alpha(\mathbf{r})=\sum_in*\alpha^i(\mathbf{r})=\sum*i\int\rho_i(\mathbf{r}- \mathbf{r}')\omega*\alpha^i(\mathbf{r}')\mathrm{d}\mathbf{r}'.$$ + +These convolutions are time-consuming calculations. Therefore, these equations are simplified by using a Taylor expansion around $\mathbf{r}$ for the density of each component $\rho_i$ as + +$$\rho_i(\mathbf{r}-\mathbf{r}')=\rho_i(\mathbf{r})-\nabla\rho_i(\mathbf{r})\cdot \mathbf{r}'+\frac{1}{2}\nabla\nabla\rho(\mathbf{r}):\mathbf{r}'\mathbf{r}'+\ldots$$ + +In the convolution integrals, the integration over angles can now be performed analytically for the spherically symmetric weight functions $\omega_\alpha^i(\mathbf{r})=\omega_\alpha^i(r)$ +which provides + +$$ +n_\alpha^i(\mathbf{r})=\rho_i(\mathbf{r})\underbrace{4\pi\int_0^\infty \omega_\alpha^i(r)r^2\mathrm{d} r}_{\omega_\alpha^{i0}} + +\nabla^2\rho_i(\mathbf{r})\underbrace{\frac{2}{3}\pi\int_0^\infty\omega_\alpha^i(r)r^4\mathrm{d} r}_{\omega_\alpha^{i2}}+\ldots +$$ + +with the weight constants $\omega_\alpha^{i0}$ and $\omega_\alpha^{i2}$. + +The resulting weighted densities can be split into a local part $n_\alpha^0(\mathbf{r})$ and an excess part $\Delta n_\alpha(\mathbf{r})$ as + +$$n_\alpha(\mathbf{r})=\underbrace{\sum_i\rho_i(\mathbf{r}) \omega_\alpha^{i0}}_{n_\alpha^0} +\underbrace{\sum_i\nabla^2\rho_i(\mathbf{r})\omega_\alpha^{i2}+\ldots}_{\Delta n_\alpha}.$$ + +The second simplification is the expansion of the reduced residual +Helmholtz energy density $\Phi(\{ n_\alpha\})$ around the local density approximation truncated after the second term + +$$ +\Phi(\lbrace n_\alpha\rbrace) + =\Phi(\lbrace n_\alpha^0\rbrace) + +\sum_i\sum_\alpha\frac{\partial\Phi}{\partial n_\alpha}\omega_\alpha^{i2}\nabla^2\rho_i + \ldots +$$ + +The Helmholtz energy functional (which was introduced in the section about the [Euler-Lagrange equation](euler_lagrange_equation.md)) then reads + +$$ F[\mathbf{\rho}(\mathbf{r})]=\int\left(f(\mathbf{\rho})+\sum*{ij}\frac{c*{ij}(\mathbf{\rho})}{2}\nabla\rho_i\cdot\nabla\rho_j\right)\mathrm{d}\mathbf{r}$$ + +with the density dependent influence parameter + +$$ \beta c*{ij}(\mathbf{\rho})=-\sum*{\alpha\beta}\frac{\partial^2\Phi}{\partial n*\alpha\partial n*\beta}\left(\omega*\alpha^{i2}\omega*\beta^{j0}+ \omega*\alpha^{i0}\omega*\beta^{j2}\right).$$ + +and the local Helmholtz energy density $f(\mathbf{\rho})$. + +For pure components, as derived in the original publication, the surface tension can be calculated from the surface excess grand potential per area according to + +$$ \gamma=\frac{F-\mu N+pV}{A}=\int\_{\rho^\mathrm{V}}^{\rho^\mathrm{L}} \sqrt{2c \left(f(\rho)-\rho\mu+p\right) } d\rho $$ + +Thus, no iterative solver is necessary to calculate the surface tension of pure components, which is a major advantage of pDGT. Finally, the density profile can be calculated from + +$$ z(\rho)=\int\_{\rho^\mathrm{V}}^{\rho} \sqrt{\frac{c/2}{ f(\rho)-\rho\mu+p} } d\rho $$ diff --git a/py-feos/docs/theory/dft/solver.md b/py-feos/docs/theory/dft/solver.md new file mode 100644 index 000000000..0dbc32997 --- /dev/null +++ b/py-feos/docs/theory/dft/solver.md @@ -0,0 +1,125 @@ +# DFT solvers + +Different solvers can be used to calculate the density profiles from the Euler-Lagrange equation introduced previously. The solvers differ in their stability, the rate of convergence, and the execution time. Unfortunately, the optimal solver and solver parameters depend on the studied system. + +## Picard iteration + +The form of the Euler-Lagrange equation + +$$\rho_\alpha(r)=\underbrace{\rho_\alpha^\mathrm{b}e^{\frac{\beta}{m_\alpha}\left(\hat F_{\rho_\alpha}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)}_{\equiv \mathcal{P}_\alpha(r;[\rho(r)])}$$ + +suggests using a simple fixed point iteration + +$$\rho_\alpha^{(k+1)}(r)=\mathcal{P}_\alpha\left(r;\left[\rho^{(k)}(r)\right]\right)$$ + +Except for some systems – typically at low densities – this iteration is unstable. Instead the new solution is obtained as combination of the old solution and the projected solution $\mathcal{P}$ + +$$\rho_\alpha^{(k+1)}(r)=(1-\nu)\rho_\alpha^{(k)}(r)+\nu\mathcal{P}_\alpha\left(r;\left[\rho^{(k)}(r)\right]\right)$$ + +The weighting between the old and projected solution is specified by the damping coefficient $\nu$. The expression can be rewritten as + +$$\rho_\alpha^{(k+1)}(r)=\rho_\alpha^{(k)}(r)+\nu\Delta\rho_\alpha^{(k)}(r)$$ + +with the search direction $\Delta\rho_\alpha(r)$ which is identical to the residual $\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)$ + +$$\Delta\rho_\alpha(r)=\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)\equiv\mathcal{P}_\alpha\left(r;\left[\rho(r)\right]\right)-\rho_\alpha(r)$$ + +The Euler-Lagrange equation can be reformulated as the "logarithmic" version + +$$\ln\rho_\alpha(r)=\ln\mathcal{P}_\alpha\left(r;\left[\rho(r)\right]\right)$$ + +Then repeating the same steps as above leads to the "logarithmic" Picard iteration + +$$\ln\rho_\alpha^{(k+1)}(r)=\ln\rho_\alpha^{(k)}(r)+\nu\Delta\ln\rho_\alpha^{(k)}(r)$$ + +or + +$$\rho_\alpha^{(k+1)}(r)=\rho_\alpha^{(k)}(r)e^{\nu\Delta\ln\rho_\alpha^{(k)}(r)}$$ + +with + +$$\Delta\ln\rho_\alpha(r)=\mathcal{\hat F}_\alpha\left(r;\left[\rho(r)\right]\right)\equiv\ln\mathcal{P}_\alpha\left(r;\left[\rho(r)\right]\right)-\ln\rho_\alpha(r)$$ + +## Newton algorithm + +A Newton iteration is a more refined approach to calculate the roots of the residual $\mathcal{F}$. From a Taylor expansion of the residual + +$$\mathcal{F}_\alpha\left(r;\left[\rho(r)+\Delta\rho(r)\right]\right)=\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)+\int\sum_\beta\frac{\delta\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)}{\delta\rho_\beta(r')}\Delta\rho_\beta(r')\mathrm{d}r'+\ldots$$ + +the Newton step is derived by setting the updated residual $\mathcal{F}_\alpha[\rho(r)+\Delta\rho(r)]$ to 0 and neglecting higher order terms. + +$$\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)=-\int\sum_\beta\frac{\delta\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)}{\delta\rho_\beta(r')}\Delta\rho_\beta(r')\mathrm{d}r'$$ (eqn:newton) + +The linear integral equation has to be solved for the step $\Delta\rho(r)$. Explicitly evaluating the functional derivatives of the residuals is not feasible due to their high dimensionality. Instead, a matrix-free linear solver like GMRES can be used. For GMRES only the action of the linear system on the variable is required (an evaluation of the right-hand side in the equation above for a given $\Delta\rho$). This action can be approximated numerically via + +$$\int\sum_\beta\frac{\delta\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)}{\delta\rho_\beta(r')}\Delta\rho_\beta(r')\mathrm{d}r'\approx\frac{\mathcal{F}_\alpha\left(r;\left[\rho(r)+s\Delta\rho(r)\right]\right)-\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)}{s}$$ + +However this approach requires the choice of an appropriate step size $s$ (something that we want to avoid in $\text{FeO}_\text{s}$) and also an evaluation of the full residual in every step of the linear solver. The solver can be sped up by doing parts of the functional derivative analytically beforehand. Using the definition of the residual in the rhs of eq. {eq}`eqn:newton` leads to + +$$ +\begin{align*} +q_\alpha(r)&\equiv-\int\sum_\beta\frac{\delta\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)}{\delta\rho_\beta(r')}\Delta\rho_\beta(r')\mathrm{d}r'\\ +&=\int\sum_\beta\frac{\delta}{\delta\rho_\beta(r')}\left(\rho_\alpha(r)-\rho_\alpha^\mathrm{b}e^{\frac{\beta}{m_\alpha}\left(\hat F_{\rho_\alpha}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)\right)\Delta\rho_\beta(r')\mathrm{d}r' +\end{align*} +$$ + +The functional derivative can be simplified using $\hat F_{\rho_\alpha\rho_\beta}^\mathrm{res}(r,r')=\frac{\delta \hat F_{\rho_\alpha}^\mathrm{b,res}(r)}{\delta\rho_\beta(r')}=\frac{\delta^2\hat F^\mathrm{b,res}}{\delta\rho_\alpha(r)\delta\rho_\beta(r')}$ + +$$ +\begin{align*} +q_\alpha(r)&=\int\sum_\beta\left(\delta_{\alpha\beta}\delta(r-r')+\left(\frac{\beta}{m_\alpha}\hat F_{\rho_\alpha\rho_\beta}^\mathrm{res}(r,r')-\sum_{\alpha'}\frac{1}{I_{\alpha\alpha'}(r)}\frac{\delta I_{\alpha\alpha'}(r)}{\delta\rho_\beta(r')}\right)\right.\\ +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\left.\times\rho_\alpha^\mathrm{b}e^{\frac{\beta}{m_\alpha}\left(\hat F_{\rho_\alpha}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)\right)\Delta\rho_\beta(r')\mathrm{d}r'\\ +&=\Delta\rho_\alpha(r)+\left(\frac{\beta}{m_\alpha}\underbrace{\int\sum_\beta\hat F_{\rho_\alpha\rho_\beta}^\mathrm{res}(r,r')\Delta\rho_\beta(r')\mathrm{d}r'}_{\Delta\hat F_{\rho_\alpha}^\mathrm{res}(r)}-\sum_{\alpha'}\frac{1}{I_{\alpha\alpha'}(r)}\underbrace{\int\sum_\beta\frac{\delta I_{\alpha\alpha'}(r)}{\delta\rho_\beta(r')}\Delta\rho_\beta(r')\mathrm{d}r'}_{\Delta I_{\alpha\alpha'}(r)}\right)\\ +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\times\rho_\alpha^\mathrm{b}e^{\frac{\beta}{m_\alpha}\left(\hat F_{\rho_\alpha}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r) +\end{align*} +$$ + +and finally + +$$q_\alpha(r)=\Delta\rho_\alpha(r)+\left(\frac{\beta}{m_\alpha}\Delta\hat F_{\rho_\alpha}^\mathrm{res}(r)-\sum_{\alpha'}\frac{\Delta I_{\alpha\alpha'}(r)}{I_{\alpha\alpha'}(r)}\right)\mathcal{P}_\alpha\left(r;\left[\rho(r)\right]\right)$$ (eqn:newton_rhs) + +Neglecting the second term in eq. {eq}`eqn:newton_rhs` leads to $\Delta_\alpha(r)=\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)$ which is the search direction of the Picard iteration. This observation implies that the Picard iteration is an approximation of the Newton solver that neglects the residual contribution to the Jacobian. Only using the ideal gas contribution to the Jacobian is a reasonable approximation for low densities and therefore, the Picard iteration converges quickly (with a large damping coefficient $\nu$) for low densities. + +The second functional derivative of the residual Helmholtz energy can be rewritten in terms of the weight functions. + +$$\hat F_{\rho_\alpha\rho_\beta}^\mathrm{res}(r,r')=\int\frac{\delta\hat f^\mathrm{res}(r'')}{\delta\rho_\alpha(r)\delta\rho_\beta(r')}\mathrm{d}r''=\int\sum_{\alpha\beta}\hat f^\mathrm{res}_{\alpha\beta}(r'')\frac{\delta n_\alpha(r'')}{\delta\rho_\alpha(r)}\frac{\delta n_\beta(r'')}{\delta\rho_\beta(r')}\mathrm{d}r''$$ + +Here $\hat f^\mathrm{res}_{\alpha\beta}=\frac{\partial^2\hat f^\mathrm{res}}{\partial n_\alpha\partial n_\beta}$ is the second partial derivative of the reduced Helmholtz energy density with respect to the weighted densities $n_\alpha$ and $n_\beta$. The definition of the weighted densities $n_\alpha(r)=\sum_\alpha\int\rho_\alpha(r')\omega_\alpha^i(r-r')\mathrm{d}r'$ is used to simplify the expression further. + +$$\hat F_{\rho_\alpha\rho_\beta}^\mathrm{res}(r,r')=\int\sum_{\alpha\beta}\hat f^\mathrm{res}_{\alpha\beta}(r'')\omega_\alpha^i(r''-r)\omega_\beta^j(r''-r')\mathrm{d}r''$$ + +With that, $\Delta\hat F_{\rho_\alpha}^\mathrm{res}(r)$ can be rewritten as + +$$ +\begin{align*} +\Delta\hat F_{\rho_\alpha}^\mathrm{res}(r)&=\int\sum_{\alpha\beta}\hat f^\mathrm{res}_{\alpha\beta}(r'')\omega_\alpha^i(r''-r)\underbrace{\sum_\beta\int\omega_\beta^j(r''-r')\Delta\rho_\beta(r')\mathrm{d}r'}_{\Delta n_\beta(r'')}\mathrm{d}r''\\ +&=\int\sum_\alpha\sum_\beta \hat f^\mathrm{res}_{\alpha\beta}(r'')\Delta n_\beta(r'')\omega_\alpha^i(r''-r)\mathrm{d}r'' +\end{align*}$$ (eqn:newton_F) + +To simplify the expression for $\Delta I_{\alpha\alpha'}(r)$, the recursive definition of the bond integrals is used. + +$$\begin{align*} +\Delta I_{\alpha\alpha'}(r)&=\iint\sum_\beta\frac{\delta}{\delta\rho_\beta(r'')}\left(e^{\frac{\beta}{m_{\alpha'}}\left(\hat F_{\rho_{\alpha'}}^\mathrm{b,res}-\hat F_{\rho_{\alpha'}}^\mathrm{res}(r')-V_{\alpha'}^\mathrm{ext}(r')\right)}\left(\prod_{\alpha''\neq\alpha}I_{\alpha'\alpha''}(r')\right)\omega_\mathrm{chain}^{\alpha\alpha'}(r-r')\right)\\ +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\times\Delta\rho_\beta(r'')\mathrm{d}r'\mathrm{d}r''\\ +&=\iint\sum_\beta\left(-\frac{\beta}{m_{\alpha'}}\hat F_{\rho_{\alpha'}\rho_\beta}^\mathrm{res}(r',r'')+\sum_{\alpha''\neq\alpha}\frac{1}{I_{\alpha'\alpha''}(r')}\frac{\delta I_{\alpha'\alpha''}(r')}{\delta\rho_\beta(r'')}\right)e^{\frac{\beta}{m_{\alpha'}}\left(\hat F_{\rho_{\alpha'}}^\mathrm{b,res}-\hat F_{\rho_{\alpha'}}^\mathrm{res}(r')-V_{\alpha'}^\mathrm{ext}(r')\right)}\\ +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\times\left(\prod_{\alpha''\neq\alpha}I_{\alpha'\alpha''}(r')\right)\omega_\mathrm{chain}^{\alpha\alpha'}(r-r')\Delta\rho_\beta(r'')\mathrm{d}r'\mathrm{d}r'' +\end{align*} +$$ + +Here, the definition of $\Delta\hat F_{\rho_\alpha}^\mathrm{res}(r)$ and $\Delta I_{\alpha\alpha'}(r)$ can be inserted leading to a recursive calculation of $\Delta I_{\alpha\alpha'}(r)$ similar to the original bond integrals. + +$$ +\begin{align*} +\Delta I_{\alpha\alpha'}(r)&=\int\left(-\frac{\beta}{m_{\alpha'}}\Delta\hat F_{\rho_{\alpha'}}^\mathrm{res}(r')+\sum_{\alpha''\neq\alpha}\frac{\Delta I_{\alpha'\alpha''}(r')}{I_{\alpha'\alpha''}(r')}\right)\\ +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\times e^{\frac{\beta}{m_{\alpha'}}\left(\hat F_{\rho_{\alpha'}}^\mathrm{b,res}-\hat F_{\rho_{\alpha'}}^\mathrm{res}(r')-V_{\alpha'}^\mathrm{ext}(r')\right)}\left(\prod_{\alpha''\neq\alpha}I_{\alpha'\alpha''}(r')\right)\omega_\mathrm{chain}^{\alpha\alpha'}(r-r')\mathrm{d}r' +\end{align*} +$$ + +(eqn:newton_I) + +In every iteration of GMRES, $q(r)$ needs to be evaluated from eqs. {eq}`eqn:newton_rhs`, {eq}`eqn:newton_F` and {eq}`eqn:newton_I`. The operations required for that are analogous to the calculation of weighted densities and the functional derivative in the Euler-Lagrange equation itself. Details of GMRES, including the pseudocode that the implementation in $\text{FeO}_\text{s}$ is based on, are given on [Wikipedia](https://de.wikipedia.org/wiki/GMRES-Verfahren) (German). + +The Newton solver converges exceptionally fast compared to a simple Picard iteration. The faster convergence comes at the cost of requiring multiple steps for solving the linear subsystem. With the algorithm outlined here, the evaluation of the second partial derivatives of the Helmholtz energy density is only required once for every Newton step. The GMRES algorithm only uses the very efficient convolution integrals and no additional evaluation of the model. + +## Anderson mixing + diff --git a/py-feos/docs/user-guide/working-with-states.md b/py-feos/docs/user-guide/working-with-states.md new file mode 100644 index 000000000..107d4b116 --- /dev/null +++ b/py-feos/docs/user-guide/working-with-states.md @@ -0,0 +1,328 @@ +# Working with States + +States are central to thermodynamic calculations in FeOs. This guide explains how to create states, calculate properties, and handle different scenarios. + +## State Creation + +States can be created from different combinations of intensive and extensive properties. + +### From Temperature and Pressure + +The most common way to create a state: + +```python +from feos import EquationOfState, State +from feos.parameters import PureRecord, Identifier, Parameters + +# Set up equation of state (see previous examples) +# ... eos creation code ... + +# Create state from T and p +state = State(eos, temperature=298.15, pressure=101325) +``` + +### From Temperature and Density + +When you know the density: + +```python +# Density in kg/m³ +state = State(eos, temperature=298.15, density=1000) +``` + +### From Temperature and Volume + +Specify molar volume: + +```python +# Molar volume in m³/mol +state = State(eos, temperature=298.15, volume=0.024) +``` + +### Mixed Specifications + +You can combine different property specifications: + +```python +# Temperature and enthalpy +state = State(eos, temperature=298.15, molar_enthalpy=1000) + +# Temperature and entropy +state = State(eos, temperature=298.15, molar_entropy=100) + +# Pressure and enthalpy (will iterate to find temperature) +state = State(eos, pressure=101325, molar_enthalpy=1000) +``` + +## State Properties + +Once you have a state, calculate various thermodynamic properties. + +### Basic Properties + +```python +# Primary state variables +print(f"Temperature: {state.temperature}") +print(f"Pressure: {state.pressure()}") +print(f"Density: {state.density}") +print(f"Molar volume: {state.molar_volume()}") + +# Derived properties +print(f"Compressibility factor: {state.compressibility()}") +print(f"Mass density: {state.mass_density()}") +``` + +### Energy Properties + +```python +# Internal energy and enthalpy +print(f"Molar internal energy: {state.molar_internal_energy()}") +print(f"Specific internal energy: {state.specific_internal_energy()}") +print(f"Molar enthalpy: {state.molar_enthalpy()}") +print(f"Specific enthalpy: {state.specific_enthalpy()}") + +# Entropy +print(f"Molar entropy: {state.molar_entropy()}") +print(f"Specific entropy: {state.specific_entropy()}") + +# Free energies +print(f"Molar Helmholtz energy: {state.molar_helmholtz_energy()}") +print(f"Molar Gibbs energy: {state.molar_gibbs_energy()}") +``` + +### Heat Capacities + +```python +# Heat capacities +print(f"Cp (isobaric): {state.molar_isobaric_heat_capacity()}") +print(f"Cv (isochoric): {state.molar_isochoric_heat_capacity()}") +print(f"Specific Cp: {state.specific_isobaric_heat_capacity()}") +print(f"Specific Cv: {state.specific_isochoric_heat_capacity()}") + +# Heat capacity ratio +gamma = state.molar_isobaric_heat_capacity() / state.molar_isochoric_heat_capacity() +print(f"Heat capacity ratio (γ): {gamma}") +``` + +### Transport Properties + +```python +# Speed of sound +print(f"Speed of sound: {state.speed_of_sound()}") + +# Compressibilities +print(f"Isothermal compressibility: {state.isothermal_compressibility()}") +print(f"Isentropic compressibility: {state.isentropic_compressibility()}") + +# Expansion coefficient +print(f"Thermal expansion: {state.thermal_expansion_coefficient()}") + +# Joule-Thomson coefficient +print(f"Joule-Thomson coefficient: {state.joule_thomson()}") +``` + +## Mixture States + +For multicomponent systems, you need to specify composition. + +### Binary Mixture + +```python +# Create binary mixture parameters (see parameter examples) +# ... mixture parameter setup ... + +# Specify composition with mole fractions +mixture_state = State( + mixture_eos, + temperature=298.15, + pressure=101325, + molefracs=[0.3, 0.7] # 30% component 1, 70% component 2 +) +``` + +### Partial Properties + +For mixtures, you can calculate partial molar properties: + +```python +# Partial molar volumes +partial_volumes = mixture_state.partial_molar_volume() +print(f"Partial molar volumes: {partial_volumes}") + +# Chemical potentials +chemical_potentials = mixture_state.chemical_potential() +print(f"Chemical potentials: {chemical_potentials}") + +# Fugacity coefficients +ln_phi = mixture_state.ln_phi() +print(f"ln(φ): {ln_phi}") +``` + +## Special States + +### Critical Points + +```python +# Pure component critical point +critical = State.critical_point_pure(eos) +print(f"Critical T: {critical.temperature}") +print(f"Critical p: {critical.pressure()}") +print(f"Critical ρ: {critical.density}") + +# For mixtures, critical points depend on composition +# mixture_critical = State.critical_point(mixture_eos, molefracs=[0.5, 0.5]) +``` + +### Saturation States + +```python +# Vapor pressure calculation +T = 298.15 +p_sat = State.vapor_pressure(eos, T) +print(f"Vapor pressure at {T} K: {p_sat}") + +# Boiling temperature +p = 101325 # Pa +T_boil = State.boiling_temperature(eos, p) +print(f"Boiling temperature at {p} Pa: {T_boil}") + +# Saturation density +liquid_density = State.saturated_liquid_density(eos, T) +vapor_density = State.saturated_vapor_density(eos, T) +print(f"Liquid density: {liquid_density}") +print(f"Vapor density: {vapor_density}") +``` + +## State Iterations + +Some state specifications require iterative solution. + +### Density Initialization + +When creating states, you can guide the density iteration: + +```python +# Specify initial guess for density iteration +state = State( + eos, + temperature=298.15, + pressure=101325, + density_initialization="liquid" # or "vapor" +) + +# Or provide specific initial density +state = State( + eos, + temperature=298.15, + pressure=101325, + density_initialization=1000 # kg/m³ +) +``` + +### Temperature Iterations + +For some property combinations, temperature is iterated: + +```python +# Provide initial temperature guess +state = State( + eos, + pressure=101325, + molar_enthalpy=1000, + initial_temperature=300 # K +) +``` + +## Error Handling and Troubleshooting + +### Common Failure Modes + +```python +try: + # This might fail if conditions are too extreme + extreme_state = State(eos, temperature=10000, pressure=1e10) +except RuntimeError as e: + print(f"State creation failed: {e}") + +try: + # Vapor pressure above critical temperature + high_temp_vp = State.vapor_pressure(eos, 1000) +except RuntimeError as e: + print(f"No vapor pressure above critical: {e}") +``` + +### Convergence Issues + +```python +# Adjust iteration parameters for difficult cases +state = State( + eos, + temperature=298.15, + pressure=101325, + max_iter=100, # Increase iterations + tol=1e-12, # Tighter tolerance + verbosity=feos.Verbosity.Iter # Show iteration progress +) +``` + +### Phase Detection + +```python +# When multiple solutions exist, specify which phase +liquid_state = State( + eos, + temperature=300, + pressure=p_sat, + density_initialization="liquid" +) + +vapor_state = State( + eos, + temperature=300, + pressure=p_sat, + density_initialization="vapor" +) +``` + +## Performance Considerations + +### Batch Calculations + +For many calculations at different conditions: + +```python +# More efficient than creating many individual states +temperatures = [250, 275, 300, 325, 350] +pressures = [1e5, 2e5, 3e5, 4e5, 5e5] + +results = [] +for T in temperatures: + for p in pressures: + try: + state = State(eos, temperature=T, pressure=p) + results.append({ + 'T': T, + 'p': p, + 'density': state.density.value, + 'cp': state.molar_isobaric_heat_capacity().value + }) + except RuntimeError: + # Skip failed calculations + continue + +print(f"Calculated {len(results)} state points") +``` + +### Property Caching + +Properties are calculated on demand and cached: + +```python +# First call calculates and caches +cp1 = state.molar_isobaric_heat_capacity() + +# Second call returns cached value (faster) +cp2 = state.molar_isobaric_heat_capacity() +``` + +This guide covers the fundamentals of working with states. Next, explore [phase equilibria](phase-equilibria.md) for vapor-liquid equilibrium calculations. \ No newline at end of file diff --git a/py-feos/mkdocs.yml b/py-feos/mkdocs.yml new file mode 100644 index 000000000..e96a08d5c --- /dev/null +++ b/py-feos/mkdocs.yml @@ -0,0 +1,192 @@ +site_name: FeOs +site_description: FeOs - framework for equations of state and classical density functional theory +site_url: https://feos-org.github.io/feos/ +repo_url: https://github.com/feos-org/feos +repo_name: feos-org/feos +edit_uri: edit/main/py-feos/docs/ + +extra_css: + - _static/custom_css.css + - _static/custom_css2.css + +theme: + name: material + font: + text: Fira Sans + code: Fira Sans Mono + features: + - navigation.tabs + - navigation.sections + - navigation.expand + - navigation.path + - navigation.top + - search.highlight + - search.suggest + - content.code.copy + - content.action.edit + - toc.integrate + palette: + # Palette toggle for automatic mode + - media: "(prefers-color-scheme)" + toggle: + icon: material/brightness-auto + name: Switch to light mode + + # Palette toggle for light mode + - media: "(prefers-color-scheme: light)" + scheme: default + primary: pink + accent: pink + toggle: + icon: material/brightness-7 + name: Switch to dark mode + + # Palette toggle for dark mode + - media: "(prefers-color-scheme: dark)" + scheme: feos + primary: black + toggle: + icon: material/brightness-4 + name: Switch to system preference + icon: + repo: fontawesome/brands/github + +plugins: + - search + - autorefs + - mkdocstrings: + handlers: + python: + # paths: ["feos"] + options: + # docstrings options + docstring_style: google + docstring_section_style: table # spacy or list + merge_init_into_class: true # put __init__ into class signature + show_if_no_docstring: false + docstring_options: + ignore_init_summary: false + trim_doctest_flags: true + + # general + find_stubs_package: true + show_source: false + + # headings + heading_level: 4 + show_root_heading: true + show_root_full_path: true + show_symbol_type_heading: true + show_symbol_type_toc: true + + # signatures + annotations_path: brief + show_signature_annotations: false + separate_signature: true + +markdown_extensions: + - pymdownx.arithmatex: # Render LaTeX via MathJax + generic: true + - pymdownx.superfences # Seems to enable syntax highlighting when used with the Material theme. + - pymdownx.details # Allowing hidden expandable regions denoted by ??? + - pymdownx.snippets: # Include one Markdown file into another + base_path: docs + - admonition + - tables + - toc: + toc_depth: 4 + - pymdownx.highlight: + pygments_lang_class: true + - pymdownx.extra + - pymdownx.tabbed: + alternate_style: true + +# Test: use typst +# hooks: +# - docs/hooks/typst_math.py + +# markdown_extensions: +# - pymdownx.arithmatex: # Render LaTeX via MathJax +# generic: true +# - pymdownx.highlight: +# anchor_linenums: true +# line_spans: __span +# pygments_lang_class: true +# - pymdownx.inlinehilite +# - pymdownx.snippets +# - pymdownx.superfences +# - pymdownx.tabbed: +# alternate_style: true +# - pymdownx.tasklist: +# custom_checkbox: true +# - admonition +# - pymdownx.details +# - attr_list +# - md_in_html +# - tables +# - toc: +# permalink: true + +extra_javascript: + - javascript/mathjax.js + - https://unpkg.com/mathjax@3/es5/tex-mml-chtml.js + +nav: + - Home: index.md + - Getting Started: + - Installation: getting-started/installation.md + - Quick Start: getting-started/quickstart.md + - Core Concepts: getting-started/concepts.md + - User Guide: + - Working with States: user-guide/working-with-states.md + - Equations of State: user-guide/equation-of-state.md + - Parameters: user-guide/parameters.md + - Phase Equilibria: user-guide/phase-equilibria.md + - Density Functional Theory: user-guide/density-functional-theory.md + - Parameter Estimation: user-guide/estimation.md + - Examples: + - Basic Properties: examples/basic-properties.md + - Pure Component VLE: examples/pure-component-vle.md + - Binary Mixtures: examples/binary-mixtures.md + - Parameter Loading: examples/parameter-loading.md + - DFT Profiles: examples/dft-profiles.md + - API Reference: + - api/index.md + - Misc: + - api/misc.md + - States: + - api/state/state.md + - api/state/phase_equilibrium.md + - Parameters: + - api/parameters/identifier.md + - api/parameters/records.md + - api/parameters/parameters.md + - Models: + - api/eos/eos.md + - SAFT-Type: + - api/eos/pcsaft.md + - api/eos/epcsaft.md + - api/eos/gc_pcsaft.md + - api/eos/saftvrmie.md + - api/eos/saftvrqmie.md + - Model Fluids: + - api/eos/uvtheory.md + - api/eos/pets.md + - Theory: + - Classical DFT: + - theory/dft/euler_lagrange_equation.md + - theory/dft/functional_derivatives.md + - theory/dft/solver.md + - theory/dft/derivatives.md + - theory/dft/enthalpy_of_adsorption.md + - theory/dft/ideal_gas.md + - theory/dft/pdgt.md + +extra: + social: + - icon: fontawesome/brands/github + link: https://github.com/feos-org/feos + # - icon: fontawesome/solid/paper-plane + # link: mailto:bauer@itt.uni-stuttgart.de + +copyright: Copyright © 2025 Gernot Bauer, Philipp Rehner diff --git a/py-feos/src/ad/mod.rs b/py-feos/src/ad/mod.rs index 1286f6c14..7b2e42a27 100644 --- a/py-feos/src/ad/mod.rs +++ b/py-feos/src/ad/mod.rs @@ -4,7 +4,7 @@ use numpy::{PyArray1, PyArray2, PyReadonlyArray2, ToPyArray}; use paste::paste; use pyo3::prelude::*; -#[pyclass(name = "EquationOfStateAD", eq, eq_int)] +#[pyclass(module = "feos.feos", name = "EquationOfStateAD", eq, eq_int)] #[derive(Clone, Copy, PartialEq)] pub enum PyEquationOfStateAD { PcSaftNonAssoc, diff --git a/py-feos/src/dft/adsorption/external_potential.rs b/py-feos/src/dft/adsorption/external_potential.rs index bf4d6e209..79ffa80c5 100644 --- a/py-feos/src/dft/adsorption/external_potential.rs +++ b/py-feos/src/dft/adsorption/external_potential.rs @@ -1,12 +1,12 @@ use feos_dft::adsorption::ExternalPotential; use ndarray::Array2; -use numpy::prelude::*; use numpy::PyArray1; +use numpy::prelude::*; use pyo3::prelude::*; use quantity::Length; /// A collection of external potentials. -#[pyclass(name = "ExternalPotential")] +#[pyclass(module = "feos.feos", name = "ExternalPotential")] #[derive(Clone)] pub struct PyExternalPotential(pub ExternalPotential); diff --git a/py-feos/src/dft/adsorption/mod.rs b/py-feos/src/dft/adsorption/mod.rs index 1f006be0a..18ca2045b 100644 --- a/py-feos/src/dft/adsorption/mod.rs +++ b/py-feos/src/dft/adsorption/mod.rs @@ -1,9 +1,9 @@ use super::PyDFTSolver; -use crate::eos::{parse_molefracs, PyEquationOfState}; +use crate::PyVerbosity; +use crate::eos::{PyEquationOfState, parse_molefracs}; use crate::error::PyFeosError; use crate::ideal_gas::IdealGasModel; use crate::residual::ResidualModel; -use crate::PyVerbosity; use feos_core::EquationOfState; use feos_dft::adsorption::{Adsorption, Adsorption1D, Adsorption3D}; use nalgebra::DMatrix; @@ -20,11 +20,11 @@ pub use external_potential::PyExternalPotential; pub use pore::{PyPore1D, PyPore2D, PyPore3D, PyPoreProfile1D, PyPoreProfile3D}; /// Container structure for adsorption isotherms in 1D pores. -#[pyclass(name = "Adsorption1D")] +#[pyclass(module = "feos.feos", name = "Adsorption1D")] pub struct PyAdsorption1D(Adsorption1D, ResidualModel>>>); /// Container structure for adsorption isotherms in 3D pores. -#[pyclass(name = "Adsorption3D")] +#[pyclass(module = "feos.feos", name = "Adsorption3D")] pub struct PyAdsorption3D(Adsorption3D, ResidualModel>>>); macro_rules! impl_adsorption_isotherm { diff --git a/py-feos/src/dft/adsorption/pore.rs b/py-feos/src/dft/adsorption/pore.rs index dc2cde9f9..3b2319a0c 100644 --- a/py-feos/src/dft/adsorption/pore.rs +++ b/py-feos/src/dft/adsorption/pore.rs @@ -76,10 +76,10 @@ macro_rules! impl_pore_profile { /// ------- /// Pore1D /// -#[pyclass(name = "Pore1D")] +#[pyclass(module = "feos.feos", name = "Pore1D")] pub struct PyPore1D(pub Pore1D); -#[pyclass(name = "PoreProfile1D")] +#[pyclass(module = "feos.feos", name = "PoreProfile1D")] pub struct PyPoreProfile1D( pub PoreProfile1D, ResidualModel>>>, ); @@ -175,10 +175,10 @@ impl PyPore1D { } } -#[pyclass(name = "Pore2D")] +#[pyclass(module = "feos.feos", name = "Pore2D")] pub struct PyPore2D(pub Pore2D); -#[pyclass(name = "PoreProfile2D")] +#[pyclass(module = "feos.feos", name = "PoreProfile2D")] pub struct PyPoreProfile2D( pub PoreProfile2D, ResidualModel>>>, ); @@ -266,10 +266,10 @@ impl PyPore2D { /// ------- /// Pore3D /// -#[pyclass(name = "Pore3D")] +#[pyclass(module = "feos.feos", name = "Pore3D")] pub struct PyPore3D(pub Pore3D); -#[pyclass(name = "PoreProfile3D")] +#[pyclass(module = "feos.feos", name = "PoreProfile3D")] pub struct PyPoreProfile3D( pub PoreProfile3D, ResidualModel>>>, ); diff --git a/py-feos/src/dft/interface/mod.rs b/py-feos/src/dft/interface/mod.rs index aae76c75c..f59b4f3b4 100644 --- a/py-feos/src/dft/interface/mod.rs +++ b/py-feos/src/dft/interface/mod.rs @@ -18,7 +18,7 @@ mod surface_tension_diagram; pub use surface_tension_diagram::PySurfaceTensionDiagram; /// A one-dimensional density profile of a vapor-liquid or liquid-liquid interface. -#[pyclass(name = "PlanarInterface")] +#[pyclass(module = "feos.feos", name = "PlanarInterface")] pub struct PyPlanarInterface( PlanarInterface, ResidualModel>>>, ); diff --git a/py-feos/src/dft/interface/surface_tension_diagram.rs b/py-feos/src/dft/interface/surface_tension_diagram.rs index e8a12461e..63e9a29a5 100644 --- a/py-feos/src/dft/interface/surface_tension_diagram.rs +++ b/py-feos/src/dft/interface/surface_tension_diagram.rs @@ -39,7 +39,7 @@ use std::sync::Arc; /// Returns /// ------- /// SurfaceTensionDiagram -#[pyclass(name = "SurfaceTensionDiagram")] +#[pyclass(module = "feos.feos", name = "SurfaceTensionDiagram")] pub struct PySurfaceTensionDiagram( SurfaceTensionDiagram, ResidualModel>>>, ); diff --git a/py-feos/src/dft/mod.rs b/py-feos/src/dft/mod.rs index 5986ee0d0..de13630c3 100644 --- a/py-feos/src/dft/mod.rs +++ b/py-feos/src/dft/mod.rs @@ -1,4 +1,4 @@ -use crate::eos::{parse_molefracs, PyEquationOfState}; +use crate::eos::{PyEquationOfState, parse_molefracs}; use crate::ideal_gas::IdealGasModel; use crate::residual::ResidualModel; use feos::hard_sphere::{FMTFunctional, FMTVersion}; @@ -24,7 +24,7 @@ pub(crate) use solver::{PyDFTSolver, PyDFTSolverLog}; /// Geometries of individual axes. #[derive(Clone, Copy, PartialEq)] -#[pyclass(name = "Geometry", eq, eq_int)] +#[pyclass(module = "feos.feos", name = "Geometry", eq, eq_int)] pub enum PyGeometry { Cartesian, Cylindrical, @@ -53,7 +53,7 @@ impl From for Geometry { /// Different versions of fundamental measure theory. #[derive(Clone, Copy, PartialEq)] -#[pyclass(name = "FMTVersion", eq, eq_int)] +#[pyclass(module = "feos.feos", name = "FMTVersion", eq, eq_int)] pub enum 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 WhiteBear, @@ -84,7 +84,7 @@ impl From for FMTVersion { } /// Collection of Helmholtz energy functionals. -#[pyclass(name = "HelmholtzEnergyFunctional")] +#[pyclass(module = "feos.feos", name = "HelmholtzEnergyFunctional")] #[derive(Clone)] pub struct PyHelmholtzEnergyFunctional; diff --git a/py-feos/src/dft/solvation.rs b/py-feos/src/dft/solvation.rs index 9e8706f3b..fbceb3c71 100644 --- a/py-feos/src/dft/solvation.rs +++ b/py-feos/src/dft/solvation.rs @@ -37,7 +37,7 @@ use std::sync::Arc; /// ------- /// SolvationProfile /// -#[pyclass(name = "SolvationProfile")] +#[pyclass(module = "feos.feos", name = "SolvationProfile")] pub struct PySolvationProfile( SolvationProfile, ResidualModel>>>, ); @@ -105,7 +105,7 @@ impl PySolvationProfile { /// ------- /// PairCorrelation /// -#[pyclass(name = "PairCorrelation")] +#[pyclass(module = "feos.feos", name = "PairCorrelation")] pub struct PyPairCorrelation( PairCorrelation, ResidualModel>>>, ); diff --git a/py-feos/src/dft/solver.rs b/py-feos/src/dft/solver.rs index 2c23d57c3..28b1bd417 100644 --- a/py-feos/src/dft/solver.rs +++ b/py-feos/src/dft/solver.rs @@ -16,7 +16,7 @@ use quantity::Time; /// Returns /// ------- /// DFTSolver -#[pyclass(name = "DFTSolver")] +#[pyclass(module = "feos.feos", name = "DFTSolver")] #[derive(Clone)] pub struct PyDFTSolver(pub DFTSolver); @@ -160,7 +160,7 @@ impl PyDFTSolver { } } -#[pyclass(name = "DFTSolverLog")] +#[pyclass(module = "feos.feos", name = "DFTSolverLog")] #[derive(Clone)] pub struct PyDFTSolverLog(pub DFTSolverLog); diff --git a/py-feos/src/eos/mod.rs b/py-feos/src/eos/mod.rs index 86b6a1447..44bbdc304 100644 --- a/py-feos/src/eos/mod.rs +++ b/py-feos/src/eos/mod.rs @@ -29,7 +29,7 @@ mod saftvrqmie; mod uvtheory; /// Collection of equations of state. -#[pyclass(name = "EquationOfState")] +#[pyclass(module = "feos.feos", name = "EquationOfState")] pub struct PyEquationOfState(pub Arc, ResidualModel>>); #[pymethods] diff --git a/py-feos/src/lib.rs b/py-feos/src/lib.rs index 77b14604d..9398adead 100644 --- a/py-feos/src/lib.rs +++ b/py-feos/src/lib.rs @@ -19,7 +19,7 @@ pub(crate) mod user_defined; /// Output level for phase equilibrium solvers. #[derive(Debug, Clone, Copy, PartialEq)] -#[pyclass(name = "Verbosity", eq, eq_int)] +#[pyclass(module = "feos.feos", name = "Verbosity", eq, eq_int)] pub(crate) enum PyVerbosity { /// Do not print output. None, @@ -53,7 +53,7 @@ impl From for Verbosity { #[pymodule] pub mod feos { - pub const __version__: &'static str = env!("CARGO_PKG_VERSION"); + use pyo3::prelude::*; #[pymodule_export] use super::PyVerbosity; @@ -105,4 +105,9 @@ pub mod feos { // Interface PySurfaceTensionDiagram, }; + + #[pymodule_init] + fn module_init(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add("__version__", env!("CARGO_PKG_VERSION")) + } } diff --git a/py-feos/src/parameter/chemical_record.rs b/py-feos/src/parameter/chemical_record.rs index b086a0f0f..879d18a28 100644 --- a/py-feos/src/parameter/chemical_record.rs +++ b/py-feos/src/parameter/chemical_record.rs @@ -1,4 +1,4 @@ -use super::fragmentation::{fragment_molecule, PySmartsRecord}; +use super::fragmentation::{PySmartsRecord, fragment_molecule}; use super::identifier::{PyIdentifier, PyIdentifierOption}; use crate::error::PyFeosError; use feos_core::parameter::{ChemicalRecord, Identifier}; @@ -7,7 +7,7 @@ use pyo3::prelude::*; use serde::{Deserialize, Serialize}; /// Information about segments and bonds of a molecule. -#[pyclass(name = "ChemicalRecord")] +#[pyclass(module = "feos.feos", name = "ChemicalRecord")] #[derive(Deserialize, Serialize, Debug, Clone)] pub(crate) struct PyChemicalRecord(ChemicalRecord); diff --git a/py-feos/src/parameter/fragmentation.rs b/py-feos/src/parameter/fragmentation.rs index b011f9fc5..d293d0fd0 100644 --- a/py-feos/src/parameter/fragmentation.rs +++ b/py-feos/src/parameter/fragmentation.rs @@ -9,7 +9,7 @@ use crate::error::PyFeosError; /// SMARTS code, required to fragmentize molecules into segments. #[derive(Clone, Serialize, Deserialize)] -#[pyclass(name = "SmartsRecord")] +#[pyclass(module = "feos.feos", name = "SmartsRecord")] pub(crate) struct PySmartsRecord { group: String, smarts: String, diff --git a/py-feos/src/parameter/identifier.rs b/py-feos/src/parameter/identifier.rs index 4baf195db..37277d81f 100644 --- a/py-feos/src/parameter/identifier.rs +++ b/py-feos/src/parameter/identifier.rs @@ -3,7 +3,7 @@ use pyo3::prelude::*; use serde::{Deserialize, Serialize}; /// Identifier to match on while reading parameters from files. -#[pyclass(name = "IdentifierOption", eq, eq_int)] +#[pyclass(module = "feos.feos", name = "IdentifierOption", eq, eq_int)] #[derive(Serialize, Deserialize, Debug, Clone, Copy, PartialEq)] pub enum PyIdentifierOption { Cas, @@ -43,7 +43,7 @@ impl From for IdentifierOption { } /// Different common identifiers for chemicals. -#[pyclass(name = "Identifier")] +#[pyclass(module = "feos.feos", name = "Identifier")] #[derive(Debug, Clone, Serialize, Deserialize)] pub struct PyIdentifier(pub Identifier); diff --git a/py-feos/src/parameter/mod.rs b/py-feos/src/parameter/mod.rs index 5f67e78b1..bb5364c0d 100644 --- a/py-feos/src/parameter/mod.rs +++ b/py-feos/src/parameter/mod.rs @@ -22,7 +22,7 @@ pub(crate) use model_record::{PyBinaryRecord, PyPureRecord}; pub(crate) use segment::{PyBinarySegmentRecord, PySegmentRecord}; /// Set of parameters that fully characterizes a mixture. -#[pyclass(name = "Parameters")] +#[pyclass(module = "feos.feos", name = "Parameters")] #[derive(Clone, Serialize, Deserialize)] pub struct PyParameters { pub pure_records: Vec>, @@ -497,7 +497,7 @@ impl PyParameters { /// Combination of chemical information and segment parameters that is used to /// parametrize a group-contribution model. -#[pyclass(name = "GcParameters")] +#[pyclass(module = "feos.feos", name = "GcParameters")] #[derive(Clone, Serialize, Deserialize)] pub struct PyGcParameters { chemical_records: Vec, diff --git a/py-feos/src/parameter/model_record.rs b/py-feos/src/parameter/model_record.rs index eeb06aa8e..a64cb732b 100644 --- a/py-feos/src/parameter/model_record.rs +++ b/py-feos/src/parameter/model_record.rs @@ -15,7 +15,7 @@ use serde_json::Value; #[derive(Serialize, Deserialize, Clone)] #[serde(from = "PureRecord")] #[serde(into = "PureRecord")] -#[pyclass(name = "PureRecord")] +#[pyclass(module = "feos.feos", name = "PureRecord")] pub struct PyPureRecord { #[pyo3(get)] pub identifier: PyIdentifier, @@ -128,7 +128,7 @@ impl PyPureRecord { #[derive(Serialize, Deserialize, Clone)] #[serde(from = "BinaryRecord")] #[serde(into = "BinaryRecord")] -#[pyclass(name = "BinaryRecord")] +#[pyclass(module = "feos.feos", name = "BinaryRecord")] pub struct PyBinaryRecord { #[pyo3(get)] pub id1: PyIdentifier, diff --git a/py-feos/src/parameter/segment.rs b/py-feos/src/parameter/segment.rs index 1816643ec..7a52fc449 100644 --- a/py-feos/src/parameter/segment.rs +++ b/py-feos/src/parameter/segment.rs @@ -10,7 +10,7 @@ use serde_json::Value; #[derive(Serialize, Deserialize, Clone)] #[serde(from = "SegmentRecord")] #[serde(into = "SegmentRecord")] -#[pyclass(name = "SegmentRecord")] +#[pyclass(module = "feos.feos", name = "SegmentRecord")] pub struct PySegmentRecord { #[pyo3(get)] identifier: String, @@ -113,7 +113,7 @@ impl PySegmentRecord { #[derive(Serialize, Deserialize, Clone)] #[serde(from = "BinaryRecord")] #[serde(into = "BinaryRecord")] -#[pyclass(name = "BinarySegmentRecord")] +#[pyclass(module = "feos.feos", name = "BinarySegmentRecord")] pub struct PyBinarySegmentRecord { #[pyo3(get)] pub id1: String, diff --git a/py-feos/src/phase_equilibria.rs b/py-feos/src/phase_equilibria.rs index 77f95df3d..b93622e5f 100644 --- a/py-feos/src/phase_equilibria.rs +++ b/py-feos/src/phase_equilibria.rs @@ -1,10 +1,10 @@ use crate::{ - eos::{parse_molefracs, PyEquationOfState}, + PyVerbosity, + eos::{PyEquationOfState, parse_molefracs}, error::PyFeosError, ideal_gas::IdealGasModel, residual::ResidualModel, state::{PyContributions, PyState, PyStateVec}, - PyVerbosity, }; use feos_core::{ Contributions, EquationOfState, PhaseDiagram, PhaseDiagramHetero, PhaseEquilibrium, ResidualDyn, @@ -20,7 +20,7 @@ use std::sync::Arc; use typenum::P3; /// A thermodynamic two phase equilibrium state. -#[pyclass(name = "PhaseEquilibrium")] +#[pyclass(module = "feos.feos", name = "PhaseEquilibrium")] #[derive(Clone)] pub struct PyPhaseEquilibrium( pub PhaseEquilibrium, ResidualModel>>, 2>, @@ -670,7 +670,7 @@ impl PyState { /// Returns /// ------- /// PhaseDiagram : the resulting phase diagram -#[pyclass(name = "PhaseDiagram")] +#[pyclass(module = "feos.feos", name = "PhaseDiagram")] pub struct PyPhaseDiagram(PhaseDiagram, ResidualModel>>, 2>); #[pymethods] @@ -1340,7 +1340,7 @@ impl PyPhaseDiagram { } /// Phase diagram for a binary mixture exhibiting a heteroazeotrope. -#[pyclass(name = "PhaseDiagramHetero")] +#[pyclass(module = "feos.feos", name = "PhaseDiagramHetero")] pub struct PyPhaseDiagramHetero( PhaseDiagramHetero, ResidualModel>>>, ); diff --git a/py-feos/src/state.rs b/py-feos/src/state.rs index 2b5492d85..b785648af 100644 --- a/py-feos/src/state.rs +++ b/py-feos/src/state.rs @@ -1,7 +1,7 @@ use crate::eos::parse_molefracs; use crate::{ - eos::PyEquationOfState, error::PyFeosError, ideal_gas::IdealGasModel, residual::ResidualModel, - PyVerbosity, + PyVerbosity, eos::PyEquationOfState, error::PyFeosError, ideal_gas::IdealGasModel, + residual::ResidualModel, }; use feos_core::{ Contributions, DensityInitialization, EquationOfState, FeosError, ResidualDyn, State, StateVec, @@ -15,7 +15,7 @@ use quantity::*; use std::collections::HashMap; use std::ops::{Deref, Neg, Sub}; use std::sync::Arc; -use typenum::{Quot, P3}; +use typenum::{P3, Quot}; type DpDn = Quantity>::Output>; type InvT = Quantity::Output>; @@ -24,7 +24,7 @@ type InvM = Quantity::Output>; /// Possible contributions that can be computed. #[derive(Clone, Copy, PartialEq)] -#[pyclass(name = "Contributions", eq, eq_int)] +#[pyclass(module = "feos.feos", name = "Contributions", eq, eq_int)] pub enum PyContributions { /// Only compute the ideal gas contribution IdealGas, @@ -101,7 +101,7 @@ impl From for Contributions { /// ------ /// Error /// When the state cannot be created using the combination of input. -#[pyclass(name = "State")] +#[pyclass(module = "feos.feos", name = "State")] #[derive(Clone)] pub struct PyState(pub State, ResidualModel>>>); @@ -1421,7 +1421,7 @@ impl PyState { /// Returns /// ------- /// StateVec -#[pyclass(name = "StateVec")] +#[pyclass(module = "feos.feos", name = "StateVec")] #[expect(clippy::type_complexity)] pub struct PyStateVec(Vec, ResidualModel>>>>);