feos-org / feos

FeOs - A Framework for Equations of State and Classical Density Functional Theory

FeOs - A Framework for Equations of State and Classical Density Functional Theory

The FeOs package provides Rust implementations of different equation of state and Helmholtz energy functional models and corresponding Python bindings.

from feos.eos import EquationOfState, State
from feos.pcsaft import PcSaftParameters, PcSaftRecord

# PC-SAFT parameters for methanol (Gross and Sadowski 2002)
record = PcSaftRecord(1.5255, 3.23, 188.9, kappa_ab=0.035176, epsilon_k_ab=2899.5, na=1, nb=1)

# Build an equation of state
parameters = PcSaftParameters.from_model_records([record])
eos = EquationOfState.pcsaft(parameters)

# Define thermodynamic conditions
critical_point = State.critical_point(eos)

# Compute properties
p = critical_point.pressure()
t = critical_point.temperature
print(f'Critical point for methanol: T={t}, p={p}.')

Critical point for methanol: T=531.5 K, p=10.7 MPa.

Models

The following models are currently published as part of the FeOs framework

name	description	eos	dft
pcsaft	perturbed-chain (polar) statistical associating fluid theory	✓	✓
epcsaft	electrolyte PC-SAFT	✓
gc-pcsaft	(heterosegmented) group contribution PC-SAFT	✓	✓
pets	perturbed truncated and shifted Lennard-Jones mixtures	✓	✓
uvtheory	equation of state for Mie fluids and mixtures	✓
saftvrqmie	equation of state for quantum fluids and mixtures	✓	✓
saftvrmie	statistical associating fluid theory for variable range interactions of Mie form	✓

The list is being expanded continuously. Currently under development are implementations of ePC-SAFT and a Helmholtz energy functional for the UV theory.

Other public repositories that implement models within the FeOs framework, but are currently not part of the feos Python package, are

name	description	eos	dft
feos-fused-chains	heterosegmented fused-sphere chain functional		✓

Parameters

In addition to the source code for the Rust and Python packages, this repository contains JSON files with previously published parameters for the different models including group contribution methods. The parameter files can be read directly from Rust or Python.

Properties and phase equilibria

The crate makes use of generalized (hyper-) dual numbers to generically calculate exact partial derivatives from Helmholtz energy equations of state. The derivatives are used to calculate

• equilibrium properties (pressure, heat capacity, fugacity, and many more),
• transport properties (viscosity, thermal conductivity, diffusion coefficients) using the entropy scaling approach
• critical points and phase equilibria for pure components and mixtures.

In addition to that, utilities are provided to assist in the handling of parameters for both molecular equations of state and (homosegmented) group contribution methods and for the generation of phase diagrams for pure components and binary mixtures.

Classical density functional theory

FeOs uses efficient numerical methods to calculate density profiles in inhomogeneous systems. Highlights include:

• Fast calculation of convolution integrals in cartesian (1D, 2D and 3D), polar, cylindrical, and spherical coordinate systems using FFT and related algorithms.
• Automatic calculation of partial derivatives of Helmholtz energy densities (including temperature derivatives) using automatic differentiation with generalized (hyper-) dual numbers.
• Modeling of heterosegmented molecules, including branched molecules.
• Functionalities for calculating surface tensions, adsorption isotherms, pair correlation functions, and solvation free energies.

Cargo features

Without additional features activated, the command

cargo test --release

will only build and test the core functionalities of the crate. To run unit and integration tests for specific models, run

cargo test --release --features pcsaft

to test, e.g., the implementation of PC-SAFT or

cargo test --release --features all_models

to run tests on all implemented models.

Python package

FeOs uses the PyO3 framework to provide Python bindings. The Python package can be installed via pip and runs on Windows, Linux and macOS:

pip install feos

If there is no compiled package for your system available from PyPI and you have a Rust compiler installed, you can instead build the python package from source using

pip install git+https://github.com/feos-org/feos

This command builds the package without link-time optimization (LTO) that can be used to increase the performance further. See the Building from source section for information about building the wheel including LTO.

Building from source

To compile the code you need the Rust compiler and maturin (>=0.13,<0.14) installed. To install the package directly into the active environment (virtualenv or conda), use

maturin develop --release

which uses the python and all_models feature as specified in the pyproject.toml file.

Alternatively, you can specify the models or features that you want to include in the python package explicitly, e.g.

maturin develop --release --features "python pcsaft dft"

for the PC-SAFT equation of state and Helmholtz energy functional.

To build wheels including link-time optimization (LTO), use

maturin build --profile="release-lto"

which will use the python and all_models features specified in the pyproject.toml file. Use the following command to build a wheel with specific features:

maturin build --profile="release-lto" --features "python ..."

LTO increases compile times measurably but the resulting wheel is more performant and has a smaller size. For development however, we recommend using the --release flag.

Documentation

For a documentation of the Python API, Python examples, and a guide to the underlying Rust framework check out the documentation.

Benchmarks

Check out the benches directory for information about provided Rust benchmarks and how to run them.

Developers

This software is currently maintained by members of the groups of

• Prof. Joachim Gross, Institute of Thermodynamics and Thermal Process Engineering (ITT), University of Stuttgart
• Prof. André Bardow, Energy and Process Systems Engineering (EPSE), ETH Zurich

Contributing

FeOs grew from the need to maintain a common codebase used within the scientific work done in our groups. We share the code publicly as a platform to publish our own research but also encourage other researchers and developers to contribute their own models or implementations of existing equations of state.

If you want to contribute to FeOs, there are several ways to go: improving the documentation and helping with language issues, testing the code on your systems to find bugs, adding new models or algorithms, or providing feature requests. Feel free to message us if you have questions or open an issue to discuss improvements.

Cite us

If you find FeOs useful for your own scientific studies, consider citing our publication accompanying this library.

@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},
}

feos's Issues

Access to I&EC Paper

Hi,
I was reading of your work, and wanted to read your in-print paper "FeOs: An Open-Source Framework for Equations of State and Classical Density Functional Theory".
But it's behind a paywall I hve no access to.

Could you perhaps post an author-copy?
But maybe I'm jumping the gun?

No offence, but it seems a contradiction in terms to have a paper on Open Source that is not itself Open source....:}
I tried SciHub and LibraryGenesis but no luck.

Add version to package

We should add __version__ to the package. I.e. add this to lib.rs

m.add("__version__", env!("CARGO_PKG_VERSION"))?;

Add second non-mixed derivative to `PartialDerivative` and `Dual2_64` to cache?

At the moment, we calculate all properties of a State that use second derivatives via HyperDual64. Second derivatives w.r.t a single property, i.e. non-mixed derivatives, can be calculated using Dual2_64 which are a bit more efficient.

Should we add Dual2_64 to the PartialDerivative enum and to the cache? E.g. we could do

// in state/mod.rs

#[derive(Clone, Copy, Eq, Hash, PartialEq, Debug)]
pub(crate) enum PartialDerivative {
    Zeroth,
    First(Derivative),
    Second(Derivative),
    SecondMixed(Derivative, Derivative),
    Third(Derivative),
}

/// Creates a [StateHD] taking the first and second (non-mixed) derivatives.
pub fn derive2(&self, derivative: Derivative) -> StateHD<Dual2_64> {
    let mut t = Dual2_64::from(self.reduced_temperature);
    let mut v = Dual2_64::from(self.reduced_volume);
    let mut n = self.reduced_moles.mapv(Dual2_64::from);
    match derivative {
        Derivative::DT => t = t.derive(),
        Derivative::DV => v = v.derive(),
        Derivative::DN(i) => n[i] = n[i].derive(),
    }
    StateHD::new(t, v, n)
}

/// Creates a [StateHD] taking the first and second mixed partial derivatives.
pub fn derive2_mixed(
    &self,
    derivative1: Derivative,
    derivative2: Derivative,
) -> StateHD<HyperDual64> {
    let mut t = HyperDual64::from(self.reduced_temperature);
    let mut v = HyperDual64::from(self.reduced_volume);
    let mut n = self.reduced_moles.mapv(HyperDual64::from);
    match derivative1 {
        Derivative::DT => t.eps1[0] = 1.0,
        Derivative::DV => v.eps1[0] = 1.0,
        Derivative::DN(i) => n[i].eps1[0] = 1.0,
    }
    match derivative2 {
        Derivative::DT => t.eps2[0] = 1.0,
        Derivative::DV => v.eps2[0] = 1.0,
        Derivative::DN(i) => n[i].eps2[0] = 1.0,
    }
    StateHD::new(t, v, n)
}

// in state/cache.rs

pub fn get_or_insert_with_d2_64<F: FnOnce() -> Dual2_64>(
    &mut self,
    derivative: Derivative,
    f: F,
) -> f64 {
    if let Some(&value) = self.map.get(&PartialDerivative::Second(derivative)) {
        self.hit += 1;
        value
    } else {
        self.miss += 1;
        let value = f();
        self.map.insert(PartialDerivative::Zeroth, value.re);
        self.map
            .insert(PartialDerivative::First(derivative), value.v1[0]);
        self.map
            .insert(PartialDerivative::Second(derivative), value.v2[0]);
        // also fill PartialDerivtave::SecondMixed(derivative, derivative)?
        value.v2[0]
    }
}

// in state/properties.rs

fn dp_dv_(&self, evaluate: Evaluate) -> QuantityScalar<U> {
        -self.get_or_compute_derivative(PartialDerivative::Second(DV), evaluate)
    }

fn dp_dt_(&self, evaluate: Evaluate) -> QuantityScalar<U> {
    -self.get_or_compute_derivative(PartialDerivative::SecondMixed(DV, DT), evaluate)
}

Background: This will speed up algorithms that use those non-mixed, second derivatives, e.g. the density iteration. First implementation shows speedup of ~9 % for density iteration and 6 % for tp_flash.

Initial temperature for bubble/dew point calculations at given pressure

Currently an initial value for the temperature is required when calculating bubble or dew points at given pressure. This is in contrast to the calculations at given temperature for which robust estimations of starting values for the pressure are available.

If some application requires bubble or dew points at a given pressure without any knowledge about the temperature of the system, heuristics need to be implemented.

Add way to specify ideal gas contribution for EoS implemented in Python

Issues with release workflows

The workflow used to publish workspace crates seems to be not working currently. For the main feos crate https://github.com/actions-rs/cargo should work but apparently the repository is not being actively maintained any more.

A possible (much simpler) option would be https://github.com/dtolnay/rust-toolchain

Enhancements for Python EoS tutorial

These points came up in #106

• we should add information about internally stored units of the State object to Python tutorial
• we should mention how the ideal gas contribution works with an EoS implemented in Python

Add benchmarks

Following #65, we should add benchmarks for possible future performance investigations.
Feel free to add to the list of functions to benchmark.

• direct calls to the residual Helmholtz energy function given a StateHD (circumventing closure construction and cache) [#89 ].
  • Helmholtz energy
  • first derivative using Dual64
  • second derivative using HyperDual64
  • third derivative using Dual3
• Properties of State (creating a State inside benchmark or manually clearing cache) [#89]
• functions including algorithms / solver
  • PhaseEquilibrium: constructors
  • critical point
• DFT
  • surface tension
  • pores
• Python
  • identifying overhead

Notes

• Compare impact of compiler options
  • lto = "thin" versus lto = true
  • target-cpu=native
  • and more?

`DataSet` with no datapoints crash

Currently, the unit of the predicted property is "inferred" by taking the unit of the 0'th element of the target vector. This panics if the vector is empty. See here.

In predict, instead of allocating an array and iterating over the index, we can use iter().map().collect().

Using the estimator tool

Is there a tutorial on how to use the estimator tool from python? I couldn't seem to find one, but I might be missing something.

Missing Contributions objects from the feos.pcsaft module

Cannot run this notebook because of above error pcsaft_phase_diagram.ipynb

Modeling copolymer systems using PC-SAFT eos

Dear doctor:
I am trying to write a python program for copolymer PC-SAFT equation of state.
In the article "Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules", Professor Gross gave all expression of residual Helmholtz energy, density, compressibility factor, fugacity coefficient, residual enthalpy and entropy of homopolymer system.
In the other article "Modeling Copolymer Systems Using the Perturbed-Chain SAFT Equation of State", he gave expressions of residual Helmholtz energy of comopolymer system, but did not provide expressions of other physical quantities.
Can you provide expressions or codes of density, compressibility factor, fugacity coefficient, residual enthalpy and entropy of copolymer system without derivative term? Thank you very much.

Theory guide: PC-SAFT

Documentation bugs and missing docstrings

API

• Add short example for each model in the api/model.md

UV Theory

• EquationOfState.UVTheory:
  • wrong data type for parameters and wrong docstring
  • missing docstring for Perturbation
  • feos.uvtheory.Perturbation: missing docstring
  • feos.uvtheory.UVParameters: "binary saft parameter records"
  • feos.uvtheory.UVRecord: describe attributes and/or default constructor.

Switch to autodiff?

In implementing teqp (https://github.com/usnistgov/teqp), I found that automatic differentiation was MUCH faster than complex step derivatives and multicomplex numbers. As you wrote the library in Rust, have you benchmarked your hyper dual derivatives with autodiff? I used this package in C++: https://github.com/autodiff/autodiff. I think you could also use this in Rust as Rust is C++-derived?

Transparency of ideal gas model selection

When using PC-SAFT it is somewhat unclear which ideal gas model is used (see, e.g., #68). Should the selection be more explicit? Some thoughts:

1. The default ideal gas model is nice to have because it leads to correct pressures and derivatives of pressures. However, calculating caloric properties with the default model should be discouraged. At the same time, residual caloric properties are fine.
2. Caloric properties could be hidden behind a trait similar to mass specific properties. That would not solve the problem for PC-SAFT, though, which clearly should implement that trait.
3. For PC-SAFT specifically there should probably be either one single method, or an explicit choice.

Then there is also the issue with models implemented in Python:

• #38

New design for `DataSet` and `Estimator`

The prediction of a target property for a given model and a set of thermodynamic states is a problem that can easily be parallelized.
Our current DataSet trait's predict method is very flexible in that a user can freely decide how to iterate over data and possibly calculate properties that are needed during that iteration. The downside of the current design is that we cannot implement the parallelization in a generic way.

This issue discusses ways to change the DataSet trait so that a generic parallel evaluation of data points can be provided.

In- and ouput of the prediction of a single data point

We could add a method into the trait that returns the prediction for a single data point.

// in trait DataSet
fn predict_datapoint(&self, eos: &Arc<E>, input: SINumber) -> Result<SINumber, EstimatorError>;

The prediction method for all data can then loop over this method. The same method could be reused for a parallel and a sequential prediction method.

However, different properties require different inputs both in terms of the number of inputs as well as the data type(s). E.g.

• vapor pressure
  • Input: temperature as SINumber
  • Output: pressure as SINumber
• binary vle using chemical potential
  • Input: x, y as f64, t, p as SINumber
  • Output: delta in chemical potentials for both species (2 SINumber)

The above predict_datapoint method would work for the vapor pressure but not for the chemical potential.

Option 1: Unify in- and output of the prediction method for a data point

We could get rid of the unit in the in- and output and be generic in the length of in- and output data.

Implementation (click)

pub trait DataSet<E: EquationOfState>: Send + Sync {
    /// This function is called prior to calculation of predictions.
    ///
    /// Can be used to calculate an internal state needed for the predictions,
    /// such as critical temperature or Antoine Parameters for vapor pressure
    /// extrapolation.
    fn prepare(&mut self, eos: &Arc<E>) -> Result<(), EstimatorError>;

    /// Prediction of the target given the model and a single data point.
    fn predict_datapoint(&self, eos: &Arc<E>, input: &[f64]) -> Vec<Result<f64, EstimatorError>>;

    /// Return an iterator over input elements.
    fn input(&self) -> &[Vec<f64>];

    /// Prediction of the target given the model for all stored data points.
    fn predict(&mut self, eos: &Arc<E>) -> Result<Vec<f64>, EstimatorError> {
        self.prepare(eos)?;
        self.input()
            .par_iter()
            .flat_map(|input| self.predict_datapoint(eos, input))
            .collect()
    }
}

• Advantage:
  • relative_difference (and thus the cost function) can be generically implemented since the result of predict is just a Vec<f64>.
  • DataSets can be stored in a Vec<Arc<dyn DataSet<E>>> in the Estimator which then can be used as is.
  • Return values of predict are always the same which is nice to check/plot how a model's parameters perform.
• Disadvantage: We have non-intuitive data types in the interfaces and (possibly) lower performance since lots of Vecs have to be allocated during the iteration (most of the time with a single value)

Option 2: Add associated types to the DataSet trait

The interface can be written in a nicer way using associated types and constants. However, the trait is no longer object safe. We'd have to build an enum with all possible implementors of DataSet so that those can be used within an Estimator inside a Vec.

Implementation (click)

pub trait DataSet<E: EquationOfState, const NTARGET: usize>: Send + Sync {
    type Input: Send + Sync;   // data type of input
    type Target: Send + Sync; // data type of output

    fn target(&self) -> &[Self::Target];

    /// This function is called prior to calculation of predictions.
    ///
    /// Can be used to calculate an internal state needed for the predictions,
    /// such as critical temperature or Antoine Parameters for vapor pressure
    /// extrapolation.
    fn prepare(&mut self, eos: &Arc<E>) -> Result<(), EstimatorError>;

    /// Prediction of the target given the model and a single data point.
    fn predict_datapoint(
        &self,
        eos: &Arc<E>,
        input: &Self::Input,
    ) -> Result<Self::Target, EstimatorError>;

    /// Return an iterator over input elements.
    fn input(&self) -> &[Self::Input];

    /// Calculate relative difference
    fn relative_difference_datapoint(
        &self,
        eos: &Arc<E>,
        input: &Self::Input,
        target: &Self::Target,
    ) -> Result<[f64; NTARGET], EstimatorError>;

    /// Prediction of the target given the model for all stored data points.
    fn predict(&mut self, eos: &Arc<E>) -> Result<Vec<Self::Target>, EstimatorError> {
        self.prepare(eos)?;
        self.input()
            .par_iter()
            .map(|input| self.predict_datapoint(eos, input))
            .collect()
    }

    fn relative_difference(&mut self, eos: &Arc<E>) -> Result<Vec<f64>, EstimatorError> {
        self.prepare(eos)?;

        // can be done without collecting twice with par_extend?
        let rel_dif = self
            .input()
            .par_iter()
            .zip(self.target())
            .map(|(input, target)| self.relative_difference_datapoint(eos, input, target))
            .collect::<Result<Vec<[f64; NTARGET]>, EstimatorError>>()?;
        Ok(rel_dif.par_iter().cloned().flatten().collect())
    }
}

• Advantage:
  • Implementation is convenient because we can use dimensioned types.
• Disadvantage:
  • The user has to implement how the relative difference between a prediction and a target can be calculated.
  • No longer object safe
  • Requires new design for Estimator.

An implementation for the vapor pressure could look like this:

Implementation (click)

pub struct VaporPressure {
    pub target: Vec<SINumber>,
    temperature: Vec<SINumber>,
    max_temperature: SINumber,
    antoine: Option<Antoine>,
    solver_options: SolverOptions,
}

impl VaporPressure {
    /// Create a new data set for vapor pressure.
    ///
    /// If the equation of state fails to compute the vapor pressure
    /// (e.g. when it underestimates the critical point) the vapor
    /// pressure can be estimated.
    /// If `extrapolate` is `true`, the vapor pressure is estimated by
    /// calculating the slope of ln(p) over 1/T.
    /// If `extrapolate` is `false`, it is set to `NAN`.
    pub fn new(
        target: SIArray1,
        temperature: SIArray1,
        extrapolate: bool,
        solver_options: Option<SolverOptions>,
    ) -> Result<Self, EstimatorError> {
        let max_temperature = temperature
            .to_reduced(KELVIN)?
            .into_iter()
            .reduce(|a, b| a.max(b))
            .unwrap()
            * KELVIN;

        let antoine = if extrapolate {
            Some(Antoine::default())
        } else {
            None
        };

        Ok(Self {
            target: target.into_iter().collect(),
            temperature: temperature.into_iter().collect(),
            max_temperature,
            antoine,
            solver_options: solver_options.unwrap_or_default(),
        })
    }
}

impl<E: EquationOfState> DataSet<E, 1> for VaporPressure {
    type Input = SINumber;
    type Target = SINumber;

    fn target(&self) -> &[Self::Target] {
        &self.target
    }

    fn prepare(&mut self, eos: &Arc<E>) -> Result<(), EstimatorError> {
        if self.antoine.is_some() {
            let critical_point = State::critical_point(
                eos,
                None,
                Some(self.max_temperature * KELVIN),
                self.solver_options,
            )?;
            let tc = critical_point.temperature;
            let pc = critical_point.pressure(Contributions::Total);

            let t0 = 0.9 * tc;
            let p0 = PhaseEquilibrium::pure(eos, t0, None, self.solver_options)?
                .vapor()
                .pressure(Contributions::Total);
            self.antoine = Some(Antoine::new(pc, p0, tc, t0)?);
        }
        Ok(())
    }

    fn input(&self) -> &[Self::Input] {
        &self.temperature
    }

    fn predict_datapoint(
        &self,
        eos: &Arc<E>,
        input: &Self::Input,
    ) -> Result<Self::Target, EstimatorError> {
        let vle = PhaseEquilibrium::pure(eos, *input, None, self.solver_options);
        if let Ok(vle) = vle {
            Ok(vle.vapor().pressure(Contributions::Total))
        } else if let Some(antoine) = &self.antoine {
            antoine.pressure(*input)
        } else {
            Ok(f64::NAN * PASCAL)
        }
    }

    fn relative_difference_datapoint(
        &self,
        eos: &Arc<E>,
        input: &Self::Input,
        target: &Self::Target,
    ) -> Result<[f64; 1], EstimatorError> {
        self.predict_datapoint(eos, input)
            .and_then(|prediction| Ok([((target - prediction) / target).into_value()?]))
    }
}

For more complex predictions, e.g. BinaryVleChemicalPotential this could look like so:

click

pub struct BinaryVleChemicalPotential {
    input: Vec<(f64, f64, SINumber, SINumber)>,
    target: Vec<(SINumber, SINumber)>,
}

impl BinaryVleChemicalPotential {
    pub fn new(
        temperature: SIArray1,
        pressure: SIArray1,
        liquid_molefracs: Array1<f64>,
        vapor_molefracs: Array1<f64>,
    ) -> Self {
        let target: Vec<(SINumber, SINumber)> = (0..temperature.len())
            .map(|_| (500.0 * JOULE / MOL, 500.0 * JOULE / MOL))
            .collect();
        let mut input = Vec::new();
        for (((&xi, &yi), t), p) in liquid_molefracs
            .iter()
            .zip(vapor_molefracs.iter())
            .zip(temperature.into_iter())
            .zip(pressure.into_iter())
        {
            input.push((xi, yi, t, p));
        }
        assert_eq!(input.len(), target.len());
        Self { input, target }
    }
}

impl<E: EquationOfState> DataSet<E, 2> for BinaryVleChemicalPotential {
    type Input = (f64, f64, SINumber, SINumber);
    type Target = (SINumber, SINumber);

    fn target(&self) -> &[Self::Target] {
        &self.target
    }

    fn input(&self) -> &[Self::Input] {
        &self.input
    }

    #[inline]
    fn prepare(&mut self, eos: &Arc<E>) -> Result<(), EstimatorError> {
        Ok(())
    }

    fn predict_datapoint(
        &self,
        eos: &Arc<E>,
        input: &Self::Input,
    ) -> Result<Self::Target, EstimatorError> {
        let xi = input.0;
        let yi = input.1;
        let t = input.2;
        let p = input.3;
        let liquid_moles = arr1(&[xi, 1.0 - xi]) * MOL;
        let liquid = State::new_npt(eos, t, p, &liquid_moles, DensityInitialization::Liquid)?;
        let mu_liquid = liquid.chemical_potential(Contributions::Total);
        let vapor_moles = arr1(&[yi, 1.0 - yi]) * MOL;
        let vapor = State::new_npt(eos, t, p, &vapor_moles, DensityInitialization::Vapor)?;
        let mu_vapor = vapor.chemical_potential(Contributions::Total);

        Ok((
            mu_liquid.get(0) - mu_vapor.get(0) + 500.0 * JOULE / MOL,
            mu_liquid.get(1) - mu_vapor.get(1) + 500.0 * JOULE / MOL,
        ))
    }

    fn relative_difference_datapoint(
        &self,
        eos: &Arc<E>,
        input: &Self::Input,
        target: &Self::Target,
    ) -> Result<[f64; 2], EstimatorError> {
        self.predict_datapoint(eos, input).and_then(|prediction| {
            Ok([
                ((target.0 - prediction.0) / target.0).into_value()?,
                ((target.1 - prediction.1) / target.1).into_value()?,
            ])
        })
    }
}

Option 3: Add parallel version of predict leaving the implementation to the user

• Advantage:
  • requires the minimal amount of change to the current DataSet trait
  • flexible
  • Estimator stays almost the same
• Disadvantage:
  • the user has to implement the parallel iterators
  • code has to be written twice

Readme and documentation for feos-saft

Theory guide: Critical points

Convergence problems with SAFT-VRQ-Mie in pores

The SAFT-VRQ-Mie functional converges less robustly than PC-SAFT due to the calculation of local mole fractions. Due to the Gibbs phenomenon, it is almost unavoidable that the density within the pore wall becomes 0 at some point. In that case the calculation of local mole fractions generates NaN.

One remedy for pure components is to write a dedicated implementation of the functional specifically for pure components.

Other remedies would be to "clean" the values of the Helmholtz energy density in places where the (total) density is zero or rewrite the functional without calling the equation of state.

Difference in Molar Enthalpies between EoS and DFT

Contributions.IdealGas uses different methods for EoS and DFT (see below).

Joback gives very small values for the IdealGas molar enthalpy (as fat as I tested all in the range ~1e-13 * JOULE).

import numpy as np

from feos.si import *
from feos.pcsaft import *

from feos.eos import EquationOfState, Contributions
from feos.eos import State as EosState
from feos.eos import PhaseEquilibrium as EosPhaseEquilibrium

from feos.dft import HelmholtzEnergyFunctional, PhaseEquilibrium
from feos.dft import State as DftState
from feos.dft import PhaseEquilibrium as DftPhaseEquilibrium

substances = ['nitrogen', 'methane']
params = PcSaftParameters.from_json(substances=substances, pure_path='pure_parameters.json', binary_path='binary_parameters.json', search_option=IdentifierOption.Name)

eos = EquationOfState.pcsaft(params)
dft = HelmholtzEnergyFunctional.pcsaft(params)

vle_eos = EosPhaseEquilibrium.bubble_point(eos, temperature_or_pressure=110*KELVIN, liquid_molefracs=np.array([0.1, 0.9]))
vle_dft = DftPhaseEquilibrium.bubble_point(dft, temperature_or_pressure=110*KELVIN, liquid_molefracs=np.array([0.1, 0.9]))

print(vle_eos.liquid.helmholtz_energy_contributions()) # Gives QSPR as IdealGas contribution
print(vle_dft.liquid.helmholtz_energy_contributions()) # Gives Joback as IdealGas contribution

print(vle_eos.liquid.molar_enthalpy(contributions=Contributions.ResidualNvt) - vle_dft.liquid.molar_enthalpy(contributions=Contributions.ResidualNvt)) #Gives very good agreement in the ResidualNvt contribution

print(vle_eos.liquid.molar_enthalpy(contributions=Contributions.IdealGas) - vle_dft.liquid.molar_enthalpy(contributions=Contributions.IdealGas)) # Gives large discrepancy in IdealGas contribution

Add example notebooks to documentation

Feel free to add more. I'd favor more separate notebooks with focused scope - makes it easier for user to find what they look for.

• For each equation of state: simple setup / parameter handling
• Working with units
• Dual numbers
• Phase equilibria, for pure substances and for binary mixtures
  • different constructor methods
  • phase diagram with plots
• Critical point calculation for pure substances and mixtures
• DFT / PDGT
  • Setup, density profile for vle, surface tension
  • Setup pdgt, surface tension
  • external potential, radial distribution function
  • external potential, adsorption, adsorption isotherm
  • external potential, solvation

Theory guide: Convolution integrals

The matrix appears to be singular

When trying to use the package to simulate hydrogen adsorption with saftvrqmie, the codes gives the error: "The matrix appears to be singular"

from feos.si import *  
from feos.dft import *
from feos.saftvrqmie import *
import matplotlib.pyplot as plt

parameters= SaftVRQMieParameters.from_json(
    substances=['hydrogen'], 
    pure_path='.../parameters/saftvrqmie/hammer2023.json')

func = HelmholtzEnergyFunctional.saftvrqmie(parameters, FMTVersion.WhiteBear)

temp = 77*KELVIN
psize= 200*ANGSTROM

solver = DFTSolver().picard_iteration(tol=1e-5, beta=0.05).anderson_mixing()

potential = ExternalPotential.LJ93(3.0, 100.0, 0.08)

bulk = State(func, temp, pressure=30*BAR)
porex = Pore1D(Geometry.Cartesian, psize, potential).initialize(bulk).solve(solver)

plt.plot(porex.z/(ANGSTROM), (porex.density*(METER**3/MOL)/1000).T)

The error is:

Cell In [20], line 9
      6 potential = ExternalPotential.LJ93(3.0, 100.0, 0.08)
      8 bulk = State(func, temp, pressure=30*BAR)
----> 9 porex = Pore1D(Geometry.Cartesian, psize, potential).initialize(bulk).solve(solver)
     11 plt.plot(porex.z/(ANGSTROM), (porex.density*(METER**3/MOL)/1000).T)
     13 plt.xlabel(r"$z$ / A")

RuntimeError: The matrix appears to be singular.

When I try with a bigger pore (300 A) it runs with no problems, and if I try an even smaller pore the error changes to:

Cell In [21], line 9
      6 potential = ExternalPotential.LJ93(3.0, 100.0, 0.08)
      8 bulk = State(func, temp, pressure=30*BAR)
----> 9 porex = Pore1D(Geometry.Cartesian, psize, potential).initialize(bulk).solve(solver)
     11 plt.plot(porex.z/(ANGSTROM), (porex.density*(METER**3/MOL)/1000).T)
     13 plt.xlabel(r"$z$ / A")

RuntimeError: `Picard Iteration` encountered illegal values during the iteration.

I also tried the same configuration using ´pcsaft´ and had no issues for any pore sizes. Could you please look into this errors?

Fix `pcsaft_phase_diagram` notebook

The Important objects from FeOs section mentions that Contributions can be imported from the pcsaft module but they are defined in the feos.eos module.

More general association schemes

Some recent and some not so recent publications show the advantage of incorporating association contributions for binary mixtures in which both pure components are not self-associating. Conceptually this is straightforward to implement in any SAFT model. The question is, how the additional parameters can be specified, especially because they have no effect on pure component properties.

`PcSaftParameters.from_records` can not be used without binary parameters

c1 = PureRecord(Identifier(), 1.0, PcSaftRecord(1, 3.0, 150.0))
c2 = PureRecord(Identifier(), 1.0, PcSaftRecord(1, 3.0, 150.0))
func = HelmholtzEnergyFunctional.pcsaft(PcSaftParameters.from_records([c1, c2]))

can not be run due to a missing argument. Passing an empty list results in an unwrap error.

Fix Python compatible `StateHD<D>`s and respective dual numbers

For the python equation of state, users can define their own helmholtz_energy function that takes a StateHD object as input. In the function, properties of the state can be accessed, e.g. one can write state.moles.
To make this work, there has to be a non-generic state object for each (hyper) dual number.

Currently, missing dual numbers are simply defined in feos-core. Using these definitions, for each a non-generic StateHD version is built.
There are currently two issues:

1. simply defining the dual numbers is not enough - we need arithmetic operations, etc. from num-dual (which currently does not export all variants we need to Python)
2. we do not call impl_state_hd! (which implement the getters) for all states that are built. Hence, there will be errors in Python when one tries to call the getters.

To fix this, we

1. have to properly export all dual numbers that are needed from num-dual. (Check class names in macros!)
2. properly build all state definitions
3. build state definition and implementation for HelmholtzEnergyDual in single macro?

Move `EosVariant` and `FunctionalVariant` out of `python` module

Having these enums outside of the python feature makes them usable in other contexts as they can be imported in other Rust libraries without having to compile with python.

FeOs: Defining parameters including kij for general fitting of cubic EoS

Hello,

I would like to fit other cubic EoS to experimental data using FeOs. For this I have used as a starting point core_user_define_eos.ipynb and binary_parameter_optimization.ipynb. This works fine for PC-SAFT, but reading multiple *.json files including binary parameters for a general EoS is causing an error which I do not understand. Do you have any insights? I am working with the Peng-Robinson example to start. Please refer to the attachment and the following section of code.

Code

parameters = PengRobinsonParameters.from_multiple_json(
[
(['argon'], file_pure),
(['ammonia'], file_pure)
],
binary_path=file_binary,
)
parameters

Code-generated Error

RuntimeError Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_3844\3446048090.py in
----> 1 parameters = PengRobinsonParameters.from_multiple_json(
2 [
3 (['argon'], file_pure),
4 (['ammonia'], file_pure)
5 ],

RuntimeError: invalid type: map, expected f64 at line 19 column 25

T=374K.zip

Second order DFT solver

This paper describes a matrix-free Newton algorithm for DFT. The action of the Jacobian on the solution is calculated efficiently using FFT. This can be combined with a linear solver like GMRES that requires only the product of the Jacobian with the current x value.

Convergence could be improved, but as of right now it looks like we would need to implement the GMRES algorithm ourselves.

Export `DQVariants` to Python

DQVariants appears in the pc_saft constructor but it is currently not exported in the pcsaft module

Add example for `Estimator` and `DataSet` objects

We should add an example notebook on how to use the DataSet and Estimator objects to optimize EoS and entropy scaling parameters. Probably when possible changes to the API (due to parallelization?) are done.

Overhaul variants of `EosError`

EosError grew organically while writing core. For 0.2. we should tidy up this enum a bit and maybe introduce a variant that can be used to report more generic errors.

Theory guide

I want to include a theory guide with the documentation of FeOs. The scope should be clearly separated from the user and developer guides that are already part of the documentation. That means, it is supposed to focus on the theory underlying the code and not on implementation details. Basically, it is supposed to outline which equations the code aims to solve and which models are included.

The following topics could at some point be part of the guide:

• EOS

  • #99
  • Entropy scaling
  • #101
  • Tp-flash
  • bubble and dew points

• DFT

  • Euler-Lagrange equation
  • Functional derivatives
  • #121
  • #74
  • grand potential
  • 3D DFT (in particular the calculation of external potentials)
  • planar interfaces
  • #116

• models

  • #148
  • #149
  • #72
  • UV?

If you would like to contribute to any of these or other topics, leave a comment below or open a pull request :)

Theory guide: pDGT

core_dual_numbers example - confusion with units and references

hi,
I was reading and trying the first example 'core_dual_numbers'
and got confused with some aspects:
(some of which apply to the example "core_user_defined_eos.ipynb as well")

see attachment notebook examples.zip

in the rev1 attachment i highlight some points of confusion
in the rev2 attachment i try to refactor the Helmholtz function of the EOS to S.I. and Reid/Poling/Prausnitz notation, and try to get the same results as before for the dimensionless results (a_kt)

main issues:

1. results for a_kt of the pyfunc and a_kt of the "state.helmholtz_energy()/(KB300KELVIN)" method does not seem to match - what does that mean?

see values in the attachment

2. the reference of the chemical potential, and its relation with fugacity coefficient is not clear for me,
   also lnphi and lnphi_pure do not seem to match, what does that mean?
   (in my opinion the chemical potential / lnphi is a very motivating application for autodiffing)

see values in the attachment

3. i observed that the pyobject STATE, in the notebook body, contains quantities with units,
   but the type of the object state that arrives at the Helmholtz function of the class is of the type "<class 'builtins.StateF'>"
   and contains floats with implicit units - that also caused me some confusion when trying to take advantage of the units system inside the Helmholtz function, perhaps rename the argument of the "def helmholtz_energy(self, state):" to
   "def helmholtz_energy(self, state_floats):" / warn in the docstring of the class in the example that the received object contains only floats and array of floats (mol fracs)

suggestions:

1. the example represents a bulk system (no confinement), but the values used for the example are of few angstrom**3 and few molecules size

please see suggestion of macroscopic values

2. the example uses values inside the helmholtz energy not in S.I., and with nonstandard notation for the a and a/(bt) or a/(brt) term

please see suggestion of my attempt at an eos reprogrammed in S.I. and with notation from reid/poling/prausnitz
where a has R**2 instead of R and helmholtz uses a term like a/(brt) instead of a/(bt)

thanks for your attention

PhaseDiagram.binary_vle() Unable to Calculation Binary Mixtures without Critical Point

When creating a binary phase diagram using PhaseDiagram.binary_vle() the current algorithm is not able to handle systems where a critical point does not exist. This would be the case, for example, for the mixture Ar-NH3 at 298.15 K.

T=298K.zip

Missing argument in

Signature:
State.critical_point_binary(
    eos,
    temperature_or_pressure,
    #initial_temperature=None,  # this line is missing
    initial_molefracs=None,
    max_iter=None,
    tol=None,
    verbosity=None,
)
Docstring:
Create a thermodynamic state at critical conditions for a binary system.

Parameters
----------
eos: EquationOfState
    The equation of state to use.
temperature_or_pressure: SINumber
    temperature_or_pressure.
initial_temperature: SINumber, optional
    An initial guess for the temperature.
initial_molefracs: [float], optional
    An initial guess for the composition.
max_iter : int, optional
    The maximum number of iterations.
tol: float, optional
    The solution tolerance.
verbosity : Verbosity, optional
    The verbosity.

Returns
-------
State : State at critical conditions.
Type:      builtin_function_or_method

Third Virial Coefficient of uv-BH implementation returns NaN

Example in Python (as in lj_models.ipynb):

sigma = 3.7039
eps_k = 150.03
parameters = UVParameters.new_simple(12.0, 6.0, sigma, eps_k)
uvtheory_bh = EquationOfState.uvtheory(parameters, perturbation=Perturbation.BarkerHenderson)
uvtheory_bh.third_virial_coefficient(6 * KELVIN)

returns NaN m³/mol²
Works for uv-B3 & uv-WCA implementation.

C/C++ interface

I would like to experiment with the very exciting uv-theory stuff in my code, and I can call that from Python, but to integrate with the tools I am developing, I need a C/C++ interface. Is it possible to develop such an interface easily with rust?

Add more methods to SINumber

Thanks for this package! The SINumber wrapper in python does not seem to have the same methods as the equivalent object in rust (see here). For example, I cannot call .value on a SINumber. Would it be possible to add that functionality? Below is an example:

from feos.si import *
T = 300*KELVIN
T.value

AttributeError                            Traceback (most recent call last)
[<ipython-input-55-7f7315c6633b>](https://localhost:8080/#) in <module>
      1 from feos.si import *
      2 T = 300*KELVIN
----> 3 T.value

AttributeError: 'si_units.SINumber' object has no attribute 'value'

Changes needed for maturin >= 0.13

With version 0.13 of maturin, there are some changes that will crash our current CI workflow for building wheels and recommendations we give in our README.

• --no-sdist is no longer needed
• cargo-extra-args is no longer needed. Args to cargo build work out of the box for maturin build.

Theory guide: DFT solvers

`PairPotential` in PC-SAFT and PeTS

After the changes in feos-dft that should simplify the calculation of pair correlation functions in mixtures, the new interface for pair potentials has not yet been implemented for PC-SAFT and PeTS.

Solving `DFTProfile` for binary mixture with `x1 = 1` results in singular matrix

parameters = PcSaftParameters.from_json(
    ['toluene', 'acetone'], 
    '../parameters/pcsaft/loetgeringlin2018.json'
)
pcsaft = HelmholtzEnergyFunctional.pcsaft(parameters)
vle = PhaseDiagram.binary_vle(pcsaft, 300 * KELVIN, 251)
gamma = [PlanarInterface.from_tanh(pe, 1024, 100 * ANGSTROM, 500 * KELVIN).solve().surface_tension() for pe in vle.states]

results in RuntimeError: The matrix appears to be singular. for both end-points (pure substances) of the phase diagram.

The above functionality can be done with a SurfaceTensionDiagram but would it be helpful to handle PhaseEquilibrium inputs for mixtures with x=1?

"Hello world" example doesn't run

The main example in the README doesn't work:

from feos.eos import EquationOfState, State
from feos.pcsaft import PcSaftParameters

# Build an equation of state
parameters = PcSaftParameters.from_json(['methanol'], 'parameters.json')
eos = EquationOfState.pcsaft(parameters)

# Define thermodynamic conditions
critical_point = State.critical_point(eos)

# Compute properties
p = critical_point.pressure()
t = critical_point.temperature
print(f'Critical point for methanol: T={t}, p={p}.')

I think because the path to "parameters.json" is invalid. Looking into the wheel, it looks like the parameter files are not provided with the wheel. Could you please include them with the wheel?

Citations

• Add CITATION file.
• "How to cite" section in README?
• Add link to Hammer et al. for SAFT-VRQ-Mie parameters (here)

Reuse `HardSphere` implementation within PeTS and UV models

The implementation of the PeTS and UV equations of state might be simplified by reusing the implementation of the BMCSL equation of state that is used in the SAFT models.

Theory guide: Properties