feos-org · prehner · Dec 4, 2022 · Nov 11, 2022 · Nov 11, 2022 · Nov 11, 2022
diff --git a/docs/index.md b/docs/index.md
@@ -27,7 +27,7 @@ p = critical_point.pressure()
 t = critical_point.temperature
 print(f'Critical point for methanol: T={t}, p={p}.')
 ```
-```terminal
+```output
 Critical point for methanol: T=531.5 K, p=10.7 MPa.
 ```
 ````
@@ -56,7 +56,7 @@ let p = critical_point.pressure(Contributions::Total);
 let t = critical_point.temperature;
 println!("Critical point for methanol: T={}, p={}.", t, p);
 ```
-```terminal
+```output
 Critical point for methanol: T=531.5 K, p=10.7 MPa.
 ```
 ````

diff --git a/docs/theory/dft/euler_lagrange_equation.md b/docs/theory/dft/euler_lagrange_equation.md
@@ -1,44 +1,51 @@
-## Euler-Lagrange equation
+# 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
-$$
+$$\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}$.
+What makes this expression so appealing is that the intrinsic Helmholtz energy only depends on the temperature $T$ and the density profiles $\rho_i(r)$ of the system and not on the external potential $V_i^\mathrm{ext}(r)$.
 
 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}$$
+$$\left.\frac{\delta\Omega}{\delta\rho_i(r)}\right|_{T,\mu}=F_{\rho_i}(r)-\mu_i+V_i^{\mathrm{ext}}(r)=0$$ (eqn:euler_lagrange_mu)
 
 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$$
+$$F_{\rho_i}^\mathrm{b}-\mu_i=0$$ (eqn:euler_lagrange_bulk)
 
-which can be inserted into (1) to give
+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.
 
-$$F_{\rho_i}(r)=F_{\rho_i}^\mathrm{b}-V_i^\mathrm{ext}(r)\tag{2}$$
+Using eq. {eq}`eqn:euler_lagrange_bulk` in eq. {eq}`eqn:euler_lagrange_mu` leads to the Euler-Lagrange equation
 
-### Spherical molecules
+$$\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 de Broglie wavelength $\Lambda_i$. The functional derivatives for an inhomogeneous and a bulk system follow as
+with the thermal 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}(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. (2) and solving for the density results in
+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 the common form of the Euler-Lagrange equation for spherical molecules.
+which is convenient because it leads directly to a recurrence relation known as Picard iteration.
 
-### Homosegmented chains
+## 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$$
@@ -61,22 +68,22 @@ $$\beta F=\sum_i\int\rho_i(r)m_i\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\rig
 
 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}(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
+## 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
+with the bond integrals $I_{\alpha\alpha'}(r)$ that are calculated recursively from
 
-$$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$$
+$$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
@@ -93,9 +100,11 @@ $$\rho_\alpha(r)=\rho_\alpha^\mathrm{b}e^{\beta\left(\hat F_{\rho_\alpha}^\mathr
 
 $$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
+## 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/docs/theory/dft/functional_derivatives.md b/docs/theory/dft/functional_derivatives.md
@@ -1,33 +1,33 @@
-## Functional derivatives
+# 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$$
+$$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')dr'\tag{1}$$
+$$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)}dr'=\int\sum_\gamma f_{n_\gamma}(r')\frac{\delta n_\gamma(r')}{\delta\rho_\alpha(r)}dr'$$
+$$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)}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
+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, 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}$$
+$$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
 

diff --git a/docs/theory/dft/index.md b/docs/theory/dft/index.md
@@ -7,6 +7,7 @@ This section explains the implementation of the core expressions from classical
 
    euler_lagrange_equation
    functional_derivatives
+   solver
 ```
 
 It is currently still under construction. You can help by [contributing](https://github.com/feos-org/feos/issues/70).