The fundamental expression in classical density functional theory is the relation between the grand potential $\Omega$ and the intrinsic Helmholtz energy $F$.

What makes this expression so appealing is that the intrinsic Helmholtz energy does only depend on the temperature $T$ and the density profiles $\rho_i(r)$ of the system and not on the external potential $V_i^\mathrm{ext}$.

For a given temperature $T$, chemical potentials $\mu$ and external potentials $V^\mathrm{ext}(r)$ the grand potential reaches a minimum at equilibrium. Mathematically this condition can be written as

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. (1) is commonly referred to as the Euler-Lagrange equation, an implicit nonlinear integral equation which needs to be solved for the density profiles of the system.

For a homogeneous (bulk) system, $V^\mathrm{ext}=0$ and we get

which can be inserted into (1) to give

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:

with the de Broglie wavelength $\Lambda_i$. The functional derivatives for an inhomogeneous and a bulk system follow as

Using these expressions in eq. (2) and solving for the density results in

which is the common form of the Euler-Lagrange equation for spherical molecules.

For chain molecules that do not resolve individual segments (essentially the PC-SAFT Helmholtz energy functional) a chain contribution is introduced as

Here, $m_i$ is the number of segments (i.e., the PC-SAFT chain length parameter), $y_{ii}$ the cavity correlation function at contact in the reference fluid, and $\lambda_i$ a weighted density.

The presence of $\rho(r)$ in the logarithm poses numerical problems. Therefore, it is convenient to rearrange the expression as

$$\begin{align}

Then the total Helmholtz energy

can be rearranged to

The functional derivatives are then similar to the spherical case

and lead to a slightly modified Euler-Lagrange equation

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

with

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

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:

At this point it can be numerically useful to redistribute the bulk contributions back into the bond integrals

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:

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$.

In the last section the functional derivative

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:

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$

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

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}

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, have odd parity, i.e., $\omega(-r)=-\omega(r)$. Therefore, the sum over the weight functions needs to be split into two contributions, as

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.

This section explains the implementation of the core expressions from classical density functional theory in $\text{FeO}_\text{s}$.

It is currently still under construction. You can help by [contributing](https://github.com/feos-org/feos/issues/70).

This section explains the thermodynamic principles and algorithms used for equation of state modeling in $\text{FeO}_\text{s}$.

It is currently still under construction. You can help by [contributing](https://github.com/feos-org/feos/issues/70).

This section describes the thermodynamic models implemented in $\text{FeO}_\text{s}$.

It is currently still under construction. You can help by [contributing](https://github.com/feos-org/feos/issues/70).

- internally rewrote the implementation of the Euler-Lagrange equation to use a bulk density instead of the chemical potential as variable. [#60](https://github.com/feos-org/feos/pull/60)

use ndarray::linalg::Dot;

use num_traits::Zero;

pub(crate) struct BulkConvolver<T> {

weight_constants: Vec<Array2<T>>,

impl<T: DualNum<f64>> BulkConvolver<T> {

pub(crate) fn new(weight_functions: Vec<WeightFunctionInfo<T>>) -> Arc<dyn Convolver<T, Ix0>> {

let weight_constants = weight_functions

.into_iter()

.map(|w| w.weight_constants(Zero::zero(), 0))

.collect();

Arc::new(Self { weight_constants })

}

impl<T: DualNum<f64>> Convolver<T, Ix0> for BulkConvolver<T>

Array2<T>: Dot<Array1<T>, Output = Array1<T>>,

fn convolve(&self, _: Array0<T>, _: &WeightFunction<T>) -> Array0<T> {

unreachable!()

}

fn weighted_densities(&self, density: &Array1<T>) -> Vec<Array1<T>> {

self.weight_constants

.iter()

.map(|w| w.dot(density))

.collect()

}

fn functional_derivative(&self, partial_derivatives: &[Array1<T>]) -> Array1<T> {

self.weight_constants

.iter()

.zip(partial_derivatives.iter())

.map(|(w, pd)| pd.dot(w))

.reduce(|a, b| a + b)

.unwrap()

}

exponential: &Array<f64, D::Larger>,

exponential

let mut i = Array::ones(exponential.raw_dim());