feos-org · prehner · Nov 11, 2022 · Oct 6, 2022 · Oct 6, 2022 · Oct 6, 2022
diff --git a/docs/conf.py b/docs/conf.py
@@ -30,6 +30,7 @@
 myst_enable_extensions = [
     "colon_fence",
     "deflist",
+    "dollarmath"
 ]
 myst_heading_anchors = 3
 

diff --git a/docs/index.md b/docs/index.md
@@ -120,9 +120,9 @@ help_and_feedback
 :caption: Python
 :hidden:
 
-api/index
 tutorials/index
 recipes/index
+api/index
 ```
 
 ```{toctree}
@@ -131,4 +131,13 @@ recipes/index
 
 rustguide/index
 rust_api
+```
+
+```{toctree}
+:caption: Theory
+:hidden:
+
+theory/eos/index
+theory/dft/index
+theory/models/index
 ```
diff --git a/docs/theory/dft/euler_lagrange_equation.md b/docs/theory/dft/euler_lagrange_equation.md
@@ -0,0 +1,101 @@
+## Euler-Lagrange equation
+The fundamental expression in classical density functional theory is the relation between the grand potential $\Omega$ and the intrinsic Helmholtz energy $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 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
+
+$$\left.\frac{\delta\Omega}{\delta\rho_i(r)}\right|_{T,\mu}=F_{\rho_i}(r)-\mu_i+V_i^{\mathrm{ext}}(r)=0\tag{1}$$
+
+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
+
+$$F_{\rho_i}^\mathrm{b}-\mu_i=0$$
+
+which can be inserted into (1) to give
+
+$$F_{\rho_i}(r)=F_{\rho_i}^\mathrm{b}-V_i^\mathrm{ext}(r)\tag{2}$$
+
+### 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 de Broglie wavelength $\Lambda_i$. The functional derivatives for an inhomogeneous and a bulk system follow as
+
+$$\beta F_{\rho_i}=\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. (2) and solving for the density results in
+
+$$\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 the common form of the Euler-Lagrange equation for spherical molecules.
+
+### 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$$
+
+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}
+\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}=m_i\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta\hat{F}_{\rho_i}^\mathrm{res}$$
+
+$$\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
+
+$$I_{\alpha\alpha'}(r)=\int e^{-\beta\left(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)$$
+
+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/docs/theory/dft/functional_derivatives.md b/docs/theory/dft/functional_derivatives.md
@@ -0,0 +1,42 @@
+## 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)]dr=\int f(\lbrace n_\gamma\rbrace)dr$$
+
+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')dr'\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)}dr'=\int\sum_\gamma f_{n_\gamma}(r')\frac{\delta n_\gamma(r')}{\delta\rho_\alpha(r)}dr'$$
+
+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)}dr'=\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'')dr''dr'\\
+&=\sum_\gamma\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r'-r)dr
+\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, have odd parity, 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')dr-\sum_\gamma^\mathrm{vec}\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r-r')dr\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/docs/theory/dft/index.md b/docs/theory/dft/index.md
@@ -0,0 +1,12 @@
+# Classical density functional theory
+This section explains the implementation of the core expressions from classical density functional theory in $\text{FeO}_\text{s}$.
+
+```{eval-rst}
+.. toctree::
+   :maxdepth: 1
+
+   euler_lagrange_equation
+   functional_derivatives
+```
+
+It is currently still under construction. You can help by [contributing](https://github.com/feos-org/feos/issues/70).
diff --git a/docs/theory/eos/index.md b/docs/theory/eos/index.md
@@ -0,0 +1,9 @@
+# Equations of state
+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).
+
+```{eval-rst}
+.. toctree::
+   :maxdepth: 1
+```
diff --git a/docs/theory/models/index.md b/docs/theory/models/index.md
@@ -0,0 +1,9 @@
+# Models
+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).
+
+```{eval-rst}
+.. toctree::
+   :maxdepth: 1
+```
diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - Added `Send` and `Sync` as supertraits to `HelmholtzEnergyFunctional` and all related traits. [#57](https://github.com/feos-org/feos/pull/57)
 - Renamed the `3d_dft` feature to `rayon` in accordance to the other feos crates. [#62](https://github.com/feos-org/feos/pull/62)
+- 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)
 
 ## [0.3.2] - 2022-10-13
 ### Changed

diff --git a/feos-dft/src/convolver/mod.rs b/feos-dft/src/convolver/mod.rs
@@ -1,8 +1,10 @@
 use crate::geometry::{Axis, Geometry, Grid};
 use crate::weight_functions::*;
+use ndarray::linalg::Dot;
 use ndarray::prelude::*;
 use ndarray::{Axis as Axis_nd, RemoveAxis, ScalarOperand, Slice};
 use num_dual::*;
+use num_traits::Zero;
 use rustdct::DctNum;
 use std::ops::{AddAssign, MulAssign, SubAssign};
 use std::sync::Arc;
@@ -34,6 +36,45 @@ pub trait Convolver<T, D: Dimension>: Send + Sync {
     ) -> Array<T, D::Larger>;
 }
 
+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>
+where
+    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()
+    }
+}
+
 /// Base structure to hold either information about the weight function through
 /// `WeightFunctionInfo` or the weight functions themselves via
 /// `FFTWeightFunctions`.

diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs
@@ -445,7 +445,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
     pub fn bond_integrals<D>(
         &self,
         temperature: f64,
-        functional_derivative: &Array<f64, D::Larger>,
+        exponential: &Array<f64, D::Larger>,
         convolver: &Arc<dyn Convolver<f64, D>>,
     ) -> Array<f64, D::Larger>
     where
@@ -468,7 +468,6 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
             }
         }
 
-        let expdfdrho = functional_derivative.mapv(|x| (-x).exp());
         let mut i_graph: Graph<_, Option<Array<f64, D>>, Directed> =
             bond_weight_functions.map(|_, _| (), |_, _| None);
 
@@ -495,7 +494,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
                     if edges.clone().all(|e| e.weight().is_some()) {
                         edge_id = Some(edge.id());
                         let i0 = edges.fold(
-                            expdfdrho
+                            exponential
                                 .index_axis(Axis(0), edge.target().index())
                                 .to_owned(),
                             |acc: Array<f64, D>, e| acc * e.weight().as_ref().unwrap(),
@@ -514,7 +513,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
             }
         }
 
-        let mut i = Array::ones(functional_derivative.raw_dim());
+        let mut i = Array::ones(exponential.raw_dim());
         for node in i_graph.node_indices() {
             for edge in i_graph.edges(node) {
                 i.index_axis_mut(Axis(0), node.index())

diff --git a/feos-dft/src/interface/mod.rs b/feos-dft/src/interface/mod.rs
@@ -91,6 +91,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
         n_grid: usize,
         l_grid: QuantityScalar<U>,
         critical_temperature: QuantityScalar<U>,
+        fix_equimolar_surface: bool,
     ) -> EosResult<Self> {
         let mut profile = Self::new(vle, n_grid, l_grid)?;
 
@@ -113,13 +114,19 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
             });
 
         // specify specification
-        profile.profile.specification =
-            DFTSpecifications::total_moles_from_profile(&profile.profile)?;
+        if fix_equimolar_surface {
+            profile.profile.specification =
+                DFTSpecifications::total_moles_from_profile(&profile.profile)?;
+        }
 
         Ok(profile)
     }
 
-    pub fn from_pdgt(vle: &PhaseEquilibrium<U, DFT<F>, 2>, n_grid: usize) -> EosResult<Self> {
+    pub fn from_pdgt(
+        vle: &PhaseEquilibrium<U, DFT<F>, 2>,
+        n_grid: usize,
+        fix_equimolar_surface: bool,
+    ) -> EosResult<Self> {
         let dft = &vle.vapor().eos;
 
         if dft.component_index().len() != 1 {
@@ -161,8 +168,10 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
         )?;
 
         // specify specification
-        profile.profile.specification =
-            DFTSpecifications::total_moles_from_profile(&profile.profile)?;
+        if fix_equimolar_surface {
+            profile.profile.specification =
+                DFTSpecifications::total_moles_from_profile(&profile.profile)?;
+        }
 
         Ok(profile)
     }

diff --git a/feos-dft/src/interface/surface_tension_diagram.rs b/feos-dft/src/interface/surface_tension_diagram.rs
@@ -19,6 +19,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> SurfaceTensionDiagram<U, F> {
         n_grid: Option<usize>,
         l_grid: Option<QuantityScalar<U>>,
         critical_temperature: Option<QuantityScalar<U>>,
+        fix_equimolar_surface: Option<bool>,
         solver: Option<&DFTSolver>,
     ) -> Self {
         let n_grid = n_grid.unwrap_or(DEFAULT_GRID_POINTS);
@@ -31,17 +32,19 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> SurfaceTensionDiagram<U, F> {
                     10,
                     100.0 * U::reference_length(),
                     500.0 * U::reference_temperature(),
+                    fix_equimolar_surface.unwrap_or(false),
                 )
             } else {
                 // initialize with pDGT for single segments and tanh for mixtures and segment DFT
                 if vle.vapor().eos.component_index().len() == 1 {
-                    PlanarInterface::from_pdgt(vle, n_grid)
+                    PlanarInterface::from_pdgt(vle, n_grid, false)
                 } else {
                     PlanarInterface::from_tanh(
                         vle,
                         n_grid,
                         l_grid.unwrap_or(100.0 * U::reference_length()),
                         critical_temperature.unwrap_or(500.0 * U::reference_temperature()),
+                        fix_equimolar_surface.unwrap_or(false),
                     )
                 }
                 .map(|mut profile| {