Theoretical background

Classical density functional theory is a powerful framework for the description of inhomogeneous fluids,^[1][2][3] based on a variational principle of the grand potential $\Omega([\rho], T)$ with respect to the density profile $\rho(\vec{r})$,

\[\frac{\delta \Omega([\rho], T)}{\delta \rho(\vec{r})} = 0.\]

Writing out external, ideal and excess contributions of $\Omega([\rho], T)$ for the case of a single species^[4] and evaluating the functional minimization yields the Euler-Lagrange equation

\[\rho(\vec{r}) = \exp\left(- \beta (V_\mathrm{ext}(\vec{r}) - \mu) + c_1(\vec{r}; [\rho], T) \right),\]

where $V_\mathrm{ext}(\vec{r})$ denotes an external potential, $\mu$ is the chemical potential and $\beta = 1 / (k_B T)$ with absolute temperature $T$ and Boltzmann's constant $k_B$.

The nontrivial effects of internal interactions are captured entirely by the one-body direct correlation functional

\[c_1(\vec{r}; [\rho], T) = - \beta \frac{\delta F_\mathrm{exc}([\rho], T)}{\delta \rho(\vec{r})},\]

which follows from functional differentiation of the excess Helmholtz free energy functional $F_\mathrm{exc}([\rho], T)$. Therefore, $c_1(\vec{r}; [\rho], T)$ is itself a functional of the density profile.

For practically all types of interacting fluids, only approximate forms of $c_1(\vec{r}; [\rho], T)$ (or, equivalently, of its generating functional $F_\mathrm{exc}([\rho], T)$) are known. A well-studied system is the hard-sphere fluid, for which highly accurate approximations of $F_\mathrm{exc}[\rho]$ can be constructed using fundamental measure theory (FMT).^[5] The hard-rod fluid in one dimension is one of the few exactly solvable models, for which the excess free energy functional has first been derived by Percus.^[6] Interparticle attraction is commonly incorporated via a mean-field contribution, which can account for arising fluid phase transitions.

With some (approximate) functional $c_1(\vec{r}; [\rho], T)$ at hand, the Euler-Lagrange equation can be solved self-consistently for the equilibrium density profile $\rho(\vec{r})$. A very simple numerical scheme is Picard iteration with mixing: The density profile is updated until convergence according to

\[\rho(\vec{r}) \leftarrow (1 - \alpha) \rho(\vec{r}) + \alpha \rho_\mathrm{EL}(\vec{r}),\]

where $\rho_\mathrm{EL}(\vec{r})$ is the right hand side of the Euler-Lagrange equation evaluated with the current density profile and $0 < \alpha < 1$ is a mixing parameter. Convergence is typically declared if $\Vert \Delta \rho_\mathrm{EL}(\vec{r}) \Vert_\infty$ becomes smaller than some predefined tolerance, whereby $\Delta \rho_\mathrm{EL}(\vec{r}) = \rho(\vec{r}) - \rho_\mathrm{EL}(\vec{r})$.