Pseudospectral methods for density functional theory in bounded and unbounded domains

TL;DR

This paper introduces a novel pseudo-spectral collocation scheme for density functional theory that efficiently computes non-local terms with high accuracy, enabling better simulations of fluid phenomena in various domains.

Contribution

The paper presents an efficient pseudo-spectral method with specialized quadrature for DFT, reducing computational nodes needed and improving accuracy in complex fluid simulations.

Findings

Accurately computes non-local terms in DFT using fewer grid points.

Demonstrates effectiveness on two-dimensional fluid models with various interactions.

Results satisfy statistical mechanical sum rules.

Abstract

Classical Density Functional Theory (DFT) is a statistical-mechanical framework to analyze fluids, which accounts for nanoscale fluid inhomogeneities and non-local intermolecular interactions. DFT can be applied to a wide range of interfacial phenomena, as well as problems in adsorption, colloidal science and phase transitions in fluids. Typical DFT equations are highly non-linear, stiff and contain several convolution terms. We propose a novel, efficient pseudo-spectral collocation scheme for computing the non-local terms in real space with the help of a specialized Gauss quadrature. Due to the exponential accuracy of the quadrature and a convenient choice of collocation points near interfaces, we can use grids with a significantly lower number of nodes than most other reported methods. We demonstrate the capabilities of our numerical methodology by studying equilibrium and dynamic two-dimensional test cases with single- and multispecies hard-sphere and hard-disc particles modelled with fundamental measure theory, with and without van der Waals attractive forces, in bounded and unbounded physical domains. We show that our results satisfy statistical mechanical sum rules.

Tables

Table 1. Table 1 : Variable x 1 subscript 𝑥 1 x_{1} in the computational domain is discretized with Gauss-Lobatto-Chebyshev collocation points, x 1 ∈ [ − 1 , 1 ] subscript 𝑥 1 1 1 x_{1}\in[-1,1] , while variable x ^ 2 subscript ^ 𝑥 2 \hat{x}_{2} is discretized with an equispaced grid ( x ^ 2 ∈ [ 0 , 1 ) subscript ^ 𝑥 2 0 1 \hat{x}_{2}\in[0,1) ). Here, the physical domain represents polar coordinates ( r , φ ) 𝑟 𝜑 \left(r,\varphi\right) . If applicable, the red solid line denotes 𝒞 ​ { 𝒴 ​ ( 0 , x ^ 2 ) } 𝒞 𝒴 0 subscript ^ 𝑥 2 \mathcal{C}\left\{\mathcal{Y}\left(0,\hat{x}_{2}\right)\right\} , and each side of the red solid line contains 50 % percent 50 50\% of the collocation points.

Map from computational to physical domain	Discretisation in Cartesian coordinates
${𝒴}_{{𝒟}_1\left(R\right)}:\left(\begin{array}{c}{x}_1\\ {\hat{x}}_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}R{x}_1\\ 2\pi {\hat{x}}_2\end{array}\right)$
${𝒴}_{{𝒟}_2(R,L)}:\left(\begin{array}{c}{x}_1\\ {\hat{x}}_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}R+{𝒜}_2(x_1;L)\\ 2\pi {\hat{x}}_2\end{array}\right)$
${𝒴}_{{𝒟}_3(R_1,R_2)}:\left(\begin{array}{c}{x}_1\\ {\hat{x}}_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}{R}_1+\frac{x_1+1}{2}\left(R_2-R_1\right)\\ 2\pi {\hat{x}}_2\end{array}\right)$

Table 2. Table 2 : Variables x 1 , 2 subscript 𝑥 1 2 x_{1,2} in the computational domain are discretized with Gauss-Lobatto-Chebyshev collocation points, x 1 , 2 ∈ [ − 1 , 1 ] subscript 𝑥 1 2 1 1 x_{1,2}\in[-1,1] , while variable x ^ 2 subscript ^ 𝑥 2 \hat{x}_{2} is discretized with an equispaced grid ( x ^ 2 ∈ [ 0 , 1 ) subscript ^ 𝑥 2 0 1 \hat{x}_{2}\in[0,1) ). The physical domain represents angular spherical coordinates ( θ , ψ ) 𝜃 𝜓 \left(\theta,\psi\right) .

Map from computational to physical domain	Discretisation in Cartesian coordinates
${𝒴}_{{𝒮}_1(R,{\theta }_1,{\theta }_2)}:\left(\begin{array}{c}{x}_1\\ {\hat{x}}_2\end{array}\right)\to \left(\begin{array}{c}\theta \\ \psi \end{array}\right)=\left(\begin{array}{c}{\theta }_1+\left({\theta }_2-{\theta }_1\right)\frac{x_1+1}{2}\\ 2\pi {\hat{x}}_2\end{array}\right)$
${𝒴}_{{𝒮}_2(R,{\theta }_1,{\theta }_2)}:\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\to \left(\begin{array}{c}\theta \\ \psi \end{array}\right)=\left(\begin{array}{c}{\theta }_1+\left({\theta }_2-{\theta }_1\right)\frac{x_2+1}{2}\\ \pi \frac{1+x_2}{2}\end{array}\right)$

Table 3. Table 3 : Variables x 1 , 2 subscript 𝑥 1 2 x_{1,2} in the computational domain are discretized with Gauss-Lobatto-Chebyshev collocation points, x 1 , 2 ∈ [ − 1 , 1 ] subscript 𝑥 1 2 1 1 x_{1,2}\in[-1,1] . The physical domain represents polar coordinates ( r , φ ) 𝑟 𝜑 \left(r,\varphi\right) . We define L r = L ​ r d − R 3 ​ L + r d − R subscript 𝐿 𝑟 𝐿 subscript 𝑟 𝑑 𝑅 3 𝐿 subscript 𝑟 𝑑 𝑅 L_{r}=L\frac{r_{d}-R}{3L+r_{d}-R} and r d = | h sin ⁡ φ | subscript 𝑟 𝑑 ℎ 𝜑 r_{d}=\left|\frac{h}{\sin\varphi}\right| . Red line as in Table LABEL:tab:FourierShapes . 𝒜 2 subscript 𝒜 2 \mathcal{A}_{2} and 𝒜 2 , F subscript 𝒜 2 F \mathcal{A}_{2,\textnormal{F}} are defined in Eqs. ( 45 ) and ( 87 ), respectively.

Map from computational to physical domain	Discretisation in Cartesian coordinates
${𝒴}_{{𝒫}_1(R,L)}:\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}R+{𝒜}_2(x_1;L)\\ \pi \frac{x_2+1}{2}\end{array}\right)$
${𝒴}_{{𝒫}_2^{\pm }(R,L,h)}:\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}{𝒜}_{2,\text{F}}(x_1;R,r_d,L_r)\\ {\phi }_0^{\pm }+\left(\frac{\pi }{2}-\Delta \phi \right)\frac{x_2+1}{2}\end{array}\right)$
(${\phi }_0^{-}=\pi$, ${\phi }_0^{+}=\frac{3}{2}\pi +\Delta \phi$, $\Delta \phi =\arccos \left|\frac{h}{R}\right|$.)
${𝒴}_{{𝒫}_3(R,L,h)}:\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}{𝒜}_{2,\text{F}}(x_1;R,r_d,L_r)\\ \pi \frac{3+x_2}{2}\end{array}\right)$
${𝒴}_{{𝒫}_4^{\pm }(R_{\text{in}},R_{\text{out}},h)}:\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}{R}_{\text{in}}+\frac{1+x_1}{2}\Delta r\\ \frac{3}{2}\pi \pm {\phi }_0^{\pm }+\frac{1+x_2}{2}\left({\phi }_0^{-}-{\phi }_0^{+}\right)\end{array}\right)$
(${\phi }_0^{-}=\arccos \left|\frac{h}{R_{\text{out}}}\right|$, ${\phi }_0^{+}=\arccos \left|\frac{h}{R_{\text{in}}}\right|$, $\Delta r=\left|\frac{h}{\sin \phi }\right|-R_{\text{in}}$.)
${𝒴}_{{𝒫}_5(R_{\text{in}},R_{\text{out}},h)}:\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}{R}_{\text{in}}+\frac{1+x_1}{2}\Delta r\\ \frac{3}{2}\pi +\Delta \phi \end{array}\right)$
($\Delta \phi =\arccos \frac{h}{R_{\text{out}}}$, $\Delta r=\left|\frac{h}{\sin \phi }\right|-R_{\text{in}}$.)
${𝒴}_{{𝒫}_6(R_{\text{in}},R_{\text{out}},h)}:\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\to \left(\begin{array}{c}r\\ \phi \end{array}\right)=\left(\begin{array}{c}{R}_{\text{in}}+\frac{1+x_1}{2}\left(R_{\text{out}}-R_{\text{in}}\right)\\ \frac{\pi }{2}+x_2\left(\pi -\Delta \phi \right)\end{array}\right)$,  ($\Delta \phi =\arccos \left(\frac{h}{R_{\text{out}}}\right)$.)

Table 4. Table 4 : Assembly of geometries defined in Tabs.  LABEL:tab:FourierShapes - LABEL:tab:ShapeDiscretization to obtain the different intersections ℐ ​ ( y ) ℐ 𝑦 \mathcal{I}(y) of the support of weight χ 𝜒 \chi with the shifted half-space { ( x ′ , y ′ ) : y + y ′ > 0 } conditional-set superscript 𝑥 ′ superscript 𝑦 ′ 𝑦 superscript 𝑦 ′ 0 \{(x^{\prime},y^{\prime}):y+y^{\prime}>0\} . Here, Δ ​ φ = arcsin ⁡ y R Δ 𝜑 𝑦 𝑅 \Delta\varphi=\arcsin\frac{y}{R} , and L 𝐿 L is a mapping parameter.

Support of weight $\chi$	$ℐ\left(y\right)$
$\{𝐫:\left|𝐫\right|=R\}$
$\{𝐫:\left|𝐫\right|\le R\}$
$\{𝐫:\left|𝐫\right|\ge R\}$
$\{𝐫:R_{\text{in}}\le \left|𝐫\right|\le R_{\text{out}}\}$

Equations

$$\Omega [n]=\mathcal{F}[n]+\int n(\mathbf{r})\left\{V_{\text{ext}}(\mathbf{r})-\mu \right\}\text{d}\mathbf{r},$$

$$\frac{\delta \Omega [n]}{\delta n(\mathbf{r})}=0.$$

$$\mathcal{F}={\mathcal{F}}_{\text{id}}+{\mathcal{F}}_{\text{exc}},$$

$${\mathcal{F}}_{\text{id}}[n]={\beta }^{-1}\int n(\mathbf{r})\left(\log \left({\Lambda }^{3}n(\mathbf{r})\right)-1\right)\text{d}\mathbf{r},$$

$${\varphi }_{\text{LJ}}(r)=4\epsilon \left({\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right),$$

$${\mathcal{F}}_{\text{exc}}={\mathcal{F}}_{\text{HS}}+{\mathcal{F}}_{\text{attr}}.$$

$${\mathcal{F}}_{\text{HS}}[\{n^{(i)}\}]$$

$${{\rm Y}}_{\text{HS}}\left(\{n_{\alpha }\}\right)$$

$$n_{\alpha }\left(\mathbf{r}\right)=\sum _i\left(n^{(i)}\ast {\omega }_{\alpha ,i}\right)\left(\mathbf{r}\right)=\sum _i\int n^{(i)}\left({\mathbf{r}}^{\prime }\right){\omega }_{\alpha ,i}\left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)\text{d}{\mathbf{r}}^{\prime }.$$

$${\omega }_{2,i}\left(\mathbf{r}\right)=\delta \left(R_i-|\mathbf{r}|\right),{\omega }_{3,i}\left(\mathbf{r}\right)=\Theta \left(R_i-|\mathbf{r}|\right),{\mathbf{\omega }}_{2,i}\left(\mathbf{r}\right)=\frac{\mathbf{r}}{|\mathbf{r}|}\delta \left(R_i-|\mathbf{r}|\right),$$

$${\omega }_{1,i}\left(\mathbf{r}\right)=\frac{{\omega }_{2,i}\left(\mathbf{r}\right)}{4\pi R_i},{\omega }_{0,i}\left(\mathbf{r}\right)=\frac{{\omega }_{2,i}\left(\mathbf{r}\right)}{4\pi R_i^2},{\mathbf{\omega }}_{1,i}\left(\mathbf{r}\right)=\frac{{\mathbf{\omega }}_2\left(\mathbf{r}\right)}{4\pi R_i}.$$

$$\frac{\delta {\mathcal{F}}_{\text{HS}}[n]}{\delta n^{(i)}(\mathbf{r})}={\beta }^{-1}\sum _{\alpha }\left(\frac{\partial {{\rm Y}}_{\text{HS}}}{\partial n_{\alpha }}\ast {\hat{\omega }}_{\alpha ,i}\right)\left(\mathbf{r}\right),$$

$${{\rm Y}}_{\text{HD}}\left(\{n_{\alpha }\}\right)=$$

$$\frac{1}{4\pi \left(1-n_3^{\text{HD}}\right)}\left(\frac{19}{12}{\left(n_2^{\text{HD}}\right)}^2-\frac{5}{12}{\mathbf{n}}_2^{(1,\text{HD})}\cdot {\mathbf{n}}_2^{(1,\text{HD})}-\frac{7}{6}{\mathbf{n}}_2^{(2,\text{HD})}\cdot {\mathbf{n}}_2^{(2,\text{HD})}\right).$$

$${\omega }_{2,i}\left(\mathbf{r}\right)=\delta \left(R_i-|\mathbf{r}|\right)\text{and}{\omega }_{3,i}\left(\mathbf{r}\right)=\Theta \left(R_i-|\mathbf{r}|\right).$$

$${\omega }_{0,i}^{\text{HD}}\left(\mathbf{r}\right)=\frac{{\omega }_{2,i}\left(\mathbf{r}\right)}{2\pi R_i},{\mathbf{\omega }}_{2,i}^{(1,\text{HD})}\left(\mathbf{r}\right)=\mathbf{r}\delta \left(R_i-|\mathbf{r}|\right),{\mathbf{\omega }}_{2,i}^{(2,\text{HD})}\left(\mathbf{r}\right)=\left(\mathbf{r}\otimes \mathbf{r}\right)\delta \left(R_i-|\mathbf{r}|\right).$$

$${\mathcal{F}}_{\text{attr}}[n]$$

$$\displaystyle\phi_{\text{attr}}\left(r\right)=\varepsilon\left\{\begin{array}[]{ll}0&\text{for }r\leq\sigma\\ 4\left(\left(\displaystyle\frac{\sigma}{r}\right)^{12}-\left(\displaystyle\frac{\sigma}{r}\right)^{6}\right)&\text{for }r>\sigma\end{array}\right..$$

$$\frac{\delta {\mathcal{F}}_{\text{attr}}[n]}{\delta n(\mathbf{r})}=\left({\varphi }_{\text{attr}}\ast n\right)\left(\mathbf{r}\right)=\int {\varphi }_{\text{attr}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)n({\mathbf{r}}^{\prime })\text{d}{\mathbf{r}}^{\prime }.$$

$$\partial_{t}n({\bf r},t)={\boldsymbol{\nabla}}_{{\bf r}}\cdot\big{(}\mathbf{J}([n({\bf r},t)],{\bf r},t)\big{)},$$

$$\displaystyle\partial_{t}n({\bf r})=\frac{1}{\gamma m}{\boldsymbol{\nabla}}_{{\bf r}}\cdot\Big{[}n({\bf r},t){\boldsymbol{\nabla}}_{{\bf r}}\frac{\delta\mathcal{F}[n(\cdot,t)]}{\delta n({\bf r},t)}\Big{]},$$

$${\partial }_{t}n(\mathbf{r},t)$$

$${\partial }_t\mathbf{v}(\mathbf{r},t)$$

$$-\frac{32}{9}\pi {\epsilon }_{\text{D}}={\int }_{{\mathbb{R}}^3}{\varphi }_{\text{attr}}\left(|\mathbf{r}|\right)\text{d}\mathbf{r}.$$

$$\tilde{\mathbf{r}}$$

$$\tilde{\mu }$$

$$T_{\text{OD}}=\frac{\gamma m{\sigma }^2}{{\epsilon }_{\text{D}}}\text{and}{T}_{\text{IN}}=\sigma \sqrt{\frac{m}{{\epsilon }_{\text{D}}}},$$

$$\tilde{\gamma }=\gamma \sigma \sqrt{\frac{m}{{\epsilon }_{\text{D}}}}.$$

$$(n\ast \chi )(y)=\int \chi (y-\tilde{y})n(\tilde{y})\text{d}\tilde{y}$$

$$n\ast \chi ={\mathcal{F}}^{-1}\left\{\mathcal{F}\left\{n\right\}\cdot \mathcal{F}\left\{\chi \right\}\right\},$$

$$p_N(x)=\sum _{n=0}^{N}f(x_n)P_n(x),$$

$$P_n(x)=\prod _{m=0,m\ne n}^N\frac{x-x_m}{x_n-x_m},$$

Full text

1 Introduction

Nanoscopic effects associated with fluid interfaces play an important role in a wide range of natural phenomena and technological applications. These include the design of water-repellent surfaces, inkjet printing, oil recovery, as well as a growing number of applications in the rapidly developing fields of micro- and nanofluidics [48, 7, 21]. In the modern literature on soft-matter systems, particle-based approaches such as molecular dynamics (MD) or Monte-Carlo (MC) computations are popular tools which allow to obtain particle trajectories and phase-space configurations of statistical mechanical systems. However, the high numerical cost of direct simulations renders them intractable for large system sizes or long observation times.

In this context, classical density functional theory (DFT) has emerged in the past decades as a useful tool for investigations of nanoscale phenomena in soft-matter systems. DFT was first used to study fluid interfaces in the late 70’s [15, 17], and since then has been applied to cover a wide range of fields [72], from the modeling of atomic [60], molecular [10] and polyatomic [12, 77] systems to electrolyte solutions [31], and even water interfaces [34]. Within DFT, a fluid is described in terms of its one-body density alone. This is achieved by approximating the free energy of the system as a functional of the density and obtaining the density profile from a minimisation procedure, e.g. by employing the variational principle [17]. This is significantly less computationally demanding compared to MD-MC, whilst retaining nanoscale properties of the fluid [72].

The typically employed DFTs for atomic fluids, such as a Lennard–Jones (LJ) fluid, are based on the perturbation theory expansion of the free energy density in powers of the attractive intermolecular potential, with the reference free energy being that of a fluid with purely repulsive intermolecular interactions, such as a fluid of hard spheres [80, 65, 6]. A highly accurate free energy functional of a hard sphere fluid is provided by the Fundamental Measure Theory (FMT) introduced by Rosenfeld [56]. This employs a functional of weighted densities, defined as convolutions of the fluid density with weight functions of finite support. In the first order of the thermodynamic perturbation theory, the attractive intermolecular interactions are typically included via a mean-field approximation, as a convolution of the fluid density with the algebraically decaying attractive part of the total intermolecular potential.

The resulting integral equations are highly stiff and non-linear, and pose significant numerical challenges. Most commonly, the convolutions in these integral equations are solved by using fast Fourier transforms (FFT) [63, 46, 29, 35, 79, 40]. This ‘matrix-free’ approach has been employed in a variety of geometries in one- (1D), two- (2D) and three-dimensional (3D) domains (see Ref. [63] for a comparison of computations in different dimensions). A prerequisite for the use of FFT is a uniform cartesian grid [46], and accurate DFT computations therefore employ dense grids with 20 [40] to 50 [29] discretisation points per hard core diameter to achieve adequate resolutions of the fluid density near interfaces. An improvement can be obtained if the weight functions, which exhibit discontinuities, are Fourier-transformed analytically [35]. This was also employed in Ref. [79] in 3D with 10 discretisation points per hard sphere diameter. However, given the fact that far from the interfaces the fluid density is usually near-constant, such approaches are wasteful. Another drawback is that methods based on FFT are restricted to the use of periodic boundary conditions. Non-periodic scenarios therefore have to be performed in large periodic domains mimicking non-periodicity, which means that jumps in the density profile are included in the domain.

Alternatively, a real-space finite element approach, discretised on a cartesian grid with linear interpolation, was introduced in Refs. [20, 19, 21]. While this method is robust, it is also expensive [32]. Heroux et al. [32] reevaluated real-space approaches, realising that the high computational cost for solving DFT problems was also due to the use of preconditioners which are typically developed for partial differential equations (PDEs) but not nonlocal ones such as DFT equations. Therefore, an adapted solver was proposed in Ref. [32] and applied to hard-sphere, polymer and molecular fluids. The respective code was made publicly available with the TRAMONTO software package [1].

The aim of the present work is to introduce an alternative real-space quadrature based on the non-uniform pseudospectral discretisation, which can also be applied for unbounded physical domains. This allows the accurate discretisation of the fluid density profile with a small number of collocation points, by positioning collocation points close to fluid interfaces, therefore avoiding regions of near-constant fluid density. The convolutions are performed in real space, by separately discretising the intersection between the support of the density profile and the support of the respective weight functions, employing a quadrature scheme with spectral accuracy.

Our proposed scheme is highly accurate, efficient and fast. We have successfully applied it in a number of settings, for soft and hard-sphere fluids in 1D planar, spherically symmetric as well as 2D domains, for both equilibrium and non-equilibrium settings [23, 24, 22, 51, 52, 50]. A similar approach based on spectral methods was also used to compute the relaxation dynamics of planar films [73] and equilibrium phase transitions in confinement [74, 76, 75].

The paper is organised as follows. Secs. 2 and 3 contain succinct reviews of DFT for LJ fluids and pseudospectral methods in unbounded domains, respectively. The numerical aspects associated with the discretisation and quadrature in 2D domains are detailed in Sec. 3.4. We verify our numerical method in Sec. 4 and validate it with thermodynamic sum-rules in 1D and 2D settings, for single-fluid and multiple-species equilibrium and dynamic settings. Furthermore, we also include a comparison with stochastic sampling techniques. We conclude in Sec. 5 with the outlook and potential for further investigations.

2 Model Equations

Classical DFT is based on the theorem that the properties of an equilibrium many-body system can be uniquely described by a functional of the number density $n({\bf r})$ [17, 71, 47]. This functional has two main properties: (a) it is minimized by the equilibrium density distribution and (b) at equilibrium it corresponds to the grand potential of the system. For a single-component system, this functional is of the form

[TABLE]

where $\mathcal{F}$ is the Helmholtz free energy functional, $\mu$ is the chemical potential and $V_{\text{ext}}$ is the external potential of the system. In order to obtain the equilibrium density distribution, the minimising condition (a) may be formulated by the application of the variational principle, leading to the Euler-Lagrange equation

[TABLE]

The intrinsic Helmholtz free energy functional $\mathcal{F}[n]$ can be split into an ideal-gas contribution, $\mathcal{F}_{\text{id}}$, and an excess contribution due to the particle–particle interactions $\mathcal{F}_{\text{exc}}$

[TABLE]

with

[TABLE]

where $\beta^{-1}=k_{B}T$, $k_{B}$ is the Boltzmann constant, $T$ the temperature and $\Lambda$ the thermal wavelength. In contrast to the ideal-gas contribution, the excess free energy functional is not known exactly. When modelling a Lennard–Jones fluid with particles interacting via

[TABLE]

the excess free energy is typically split into contributions due to short-range repulsive particle–particle interactions, $\mathcal{F}_{\text{HS}}$, and long-range attractive ones, $\mathcal{F}_{\text{attr}}$, which are treated in a perturbative manner [80, 65, 6]:

[TABLE]

In (5), $\varepsilon$ and $\sigma$ denote the length and energy parameters of the Lennard–Jones potential, respectively.

As suggested by the notation, the repulsive contribution to the free energy $\mathcal{F}_{\text{HS}}$ is usually approximated by a hard-sphere potential. One of the most successful approaches to model hard-sphere fluids is FMT, first developed in 1989 by Rosenfeld [56, 59], which is based on the geometrical properties of hard spheres [18, 60]. Improvements to this original formulation have been presented in order to capture the freezing transition [57, 58, 66], or a more accurate equation of state [37, 60, 78]. These modified formulations, however, lead to equations with a functional form close to those of the original formulation, which is still widely used [40, 26, 4, 35, 18] and gives very accurate results if compared with MD or MC computations [59]. In the following, for the sake of simplicity we restrict ourselves to the implementation of the original FMT by Rosenfeld, and briefly discuss the formalism for 3D hard spheres and 2D hard disks. But, it should be noted, that our numerical framework can be easily adapted to other FMT formalisms.

2.1 Hard-sphere FMT

For a system of multiple species of hard spheres with radii $R_{i}$ and density $n^{(i)}$, Rosenfeld formulated the hard-sphere free energy as a functional of scalar- and vector-weighted densities $n_{\alpha}$, given by

[TABLE]

with

[TABLE]

The weighted densities $n_{\alpha}$ are defined through convolutions of weight functions $\omega_{\alpha,i}$ with the density field $n$:

[TABLE]

The vector-weighted densities ${\bf n}_{\alpha}$ are defined analogous to Eq. (9), with vector weights ${\boldsymbol{\omega}}_{\alpha,i}$. The weights, as obtained through the decomposition of the Mayer function [56], are

[TABLE]

where $\delta(x)$ is the Dirac delta and $\Theta(x)$ is the Heaviside step function. The remaining weights $\omega_{1},\omega_{0},{\boldsymbol{\omega}}_{1}$ depend linearly on the weights defined in (10) and are defined by

[TABLE]

We note that the system may be reduced by making use of the fact that the vector-weighted density ${\bf n}_{2}$ is the negative gradient of $n_{3}$: ${\bf n}_{2}=-{\boldsymbol{\nabla}}n_{3}$, such as applied by Merath [46]. Insertion of Eq. (7) into the Euler-Lagrange equation (2), leads to a contribution due to hard-sphere effects of

[TABLE]

where $\hat{\omega}_{\alpha,i}\left({\bf r}\right)=\omega_{\alpha,i}\left(-{\bf r}\right)$. The corresponding 2D convolution weight for a system which is invariant in one direction is given in Sec. A of the Appendix.

2.2 Hard-disk FMT

An analogous approach can be followed if hard disks are considered. Here, we follow the description by Roth et al. [61] and define the hard-disk free energy density as

[TABLE]

The weighted densities are defined analogously to (9) in 2D:

[TABLE]

Additional weights are

[TABLE]

We note that for the hard-disk case, the convolutions defining the weighted densities are defined in $\mathbb{R}^{2}$, as opposed to the weighted densities employed for the hard-sphere free energy, which are defined in $\mathbb{R}^{3}$.

2.3 Attractive free energy contribution

We follow the route of treating the attractive forces as perturbations to the repulsive forces [6], by including them in a mean-field manner as

[TABLE]

where $\phi_{\text{attr}}$ represents the attractive particle–particle interactions, modeled here with

[TABLE]

The contribution of the free energy functional $\mathcal{F}_{\text{attr}}$ to the Euler-Lagrange equation (2) is

[TABLE]

The corresponding 2D convolution weight for a system which is invariant in one direction is given in Sec. A of the Appendix.

2.4 Dynamic DFT

An analogous result to that underpinning DFT holds for the dynamics of many-body systems [11], leading to a dynamic DFT or DDFT. DDFTs typically take the form of a continuity equation

[TABLE]

where $t$ denotes the time and the challenge is to determine the flux $\mathbf{J}$, which is a functional of the density. Due to the success of DFT in describing nanoscale fluid properties, this functional is usually based on the free energy of a related equilibrium system. In particular, this ensures that the DDFT reduces to the correct DFT at equilibrium. In this way, DDFT can be regarded as a natural generalisation of DFT.

The first DDFTs for colloidal particles were obtained phenomenologically [17, 14], but there have since been a number of attempts to rigorously derive an overdamped/high friction DDFT from the Smoluchowski equation for the full $N$-body density. Marconi and Tarazona [42, 43] derived the first DDFT for pairwise interparticle potentials,

[TABLE]

where $m$ is the mass of one particle, $\gamma$ is the friction coefficient of the system and $\mathcal{F}$ is the free energy functional for a corresponding equilibrium system with the same one-body density, such as that given in (3), but including the external potential $V_{\text{ext}}$. This additional approximation is often called the adiabatic approximation. Later, this result was generalised to $N$-body interactions [3]. Hydrodynamic interactions, caused by the interplay between the flow in the bath and the colloidal movement, have also been included in more advanced DDFT models [55, 25, 23, 24].

In order to incorporate the effect of (microscopic) inertia in DDFT, two different approaches have been considered. The first uses a multiple timescale technique [41, 44, 67, 45], whilst the second takes momentum moments of the Kramers equation (the underdamped analogue of the Smoluchowski equation), leading to a continuity equation coupled to an evolution equation for the velocity [2],

[TABLE]

where $\mathbf{v}$ is the velocity distribution. This formalism has also recently been extended to include hydrodynamic interactions [25, 23, 24]. To close the infinite hierarchy of moment equations, one typically utilises the local equilibrium approximation in which the momentum distribution is assumed to be locally Maxwellian. DDFT has been successfully applied to a wide variety of problems including sedimentation, cluster formation and cell modelling, with results having been validated against the full underlying stochastic dynamics [55, 23, 24].

2.5 Nondimensionalization

We nondimensionalise the system with a length and an energy scale. For the sake of simplicity, we only consider hard-sphere or hard-disk mixtures where particles of all species have the same radius $R_{i}=R$. In particular, we take the length-scale of the Lennard–Jones potential $\sigma$ to be equal to the hard-sphere diameter $2R$, defining the length scale. As an energy scale for systems with interparticle attraction, we employ $\varepsilon_{\text{D}}$ such that

[TABLE]

For an attractive particle–particle interaction (19) without radial cutoff, we get $\varepsilon_{D}=\varepsilon$. If, however, a cutoff radius $r_{c}$ for the attractive interaction potential is employed for computational purposes, then (24) ensures that the fluids are represented by the same dimensionless bulk phase diagram, independent of the cutoff radius $r_{c}$, hence improving comparability. For pure hard-sphere systems, we set $\varepsilon_{\text{D}}=k_{B}T$. We thus introduce the following dimensionless quantities

[TABLE]

The second term in the nondimensionalisation of the chemical potential accounts for the contribution of the term $\Lambda$ in the ideal-gas free energy and leads to a constant shift of the dimensionless chemical potential.

For the dynamic computations, we introduce different time-scales for the overdamped and the inertial case, given by

[TABLE]

respectively. The dimensionless friction coefficient in the momentum equation (23) then becomes

[TABLE]

In what follows, we omit the tildes to simplify our notation.

3 Numerical Method

The accurate computation of convolutions

[TABLE]

is crucial for solving the Euler-Lagrange equation (2) and ultimately obtaining a reliable equilibrium density profile. Convolutions are needed to compute the weighted densities as given in Eq. (9), for the hard-sphere contribution to the free energy (see Eq. (12)) as well as the attractive contribution to the free energy (see Eq. (20)). The obvious choice is to employ the convolution theorem, and perform the computation of the convolutions in Fourier space [59] through

[TABLE]

where $\mathcal{F}$ is the Fourier transform. Employing FFT is of $O(N\log N)$ complexity, where $N$ is the number of Fourier modes.

Note, however, that a density profile of a liquid film adsorbed on a substrate exhibits a jump at the wall, and approaches continuously its bulk value with increasing distance to the wall. Meanwhile, any decomposition in Fourier modes requires a truncation of the physical domain and the assumption of periodicity. This entails two crucial difficulties: First, it requires a special treatment of the Gibbs phenomenon induced by the discontinuity of the density profile at the wall. Similarly, the representation of weight functions $\omega$ for FMT such as given in Eqs.(10) and (14) needs special treatment due to the short-range finite support properties of these functions [59]. Second, non-periodic scenarios have to be mimicked employing large periodic domains, therefore wasting precious computation power.

Several approaches have been presented to circumvent these difficulties. By computing the Fourier transform analytically, the accuracy may be improved substantially [35]. Alternatively, expansions of the equilibrium density profile in terms of the local curvatures at the wall have been proposed [36]. Here, however, we aim to avoid the complications introduced by computing the convolutions in Fourier space altogether by performing the convolutions in real-space. This approach, which we describe below, is based on pseudospectral methods.

3.1 Chebyshev pseudospectral method

Here we outline a real-space discretisation method for FMT-DFT employing a pseudospectral method [9]. In this framework, a function $f(x)$ defined on $x\in[-1,1]$ is represented by its functional values at the collocation points $x_{n}$. The function is then approximated by the interpolating polynomial, defined via

[TABLE]

where $P_{n}(x)$ denote the Lagrange polynomials of degree N, given by

[TABLE]

with the property

[TABLE]

To avoid the so-called Runge phenomenon that occurs in polynomial interpolation over equispaced grids, we choose the so-called Gauss-Lobatto-Chebyshev collocation points

[TABLE]

which are clustered at the endpoints of the interval. By doing so, we obtain a spectrally accurate representation of any smooth function $f(x)$ over the whole interval. Using Eq. (30), the value of $p_{N}(x)$ can, in principle, be computed anywhere in the domain. Nevertheless, it is generally known that the calculation of Eq. (31) is both costly and numerically unstable. These issues can be avoided by using the barycentric formula [5, 68]

[TABLE]

where ${\bar{w}}_{k}$ are the so-called barycentric weights, which, for the Gauss-Lobatto-Chebyshev points are given by [62]

[TABLE]

The integral of $f$ over $[-1,1]$ is computed through the use of the Clenshaw-Curtis quadrature [13]

[TABLE]

with weights

[TABLE]

3.2 Unbounded domains

When dealing with functions approaching a constant value in an unbounded domain (such as the density profile of a fluid in contact with a wall – the constant value being approached far from the wall) several numerical techniques can be used, such as [9, 64]:

1. (1)

domain truncation, 2. (2)

approximation by orthogonal systems in the unbounded domain, or 3. (3)

mapping of the unbounded domain onto a bounded domain and applying standard spectral methods.

Option (1) is appropriate when dealing with quickly-decaying functions. It has been widely used in DFT computations. As far as options (2) and (3) are concerned, the main difference is that in (2) the equations are studied in the physical (unbounded) domain, while in (3) the equations are considered in the mapped bounded domain, which is subsequently referred to as the computational domain. Here, we adopt (3) as its implementation provides a greater flexibility in studying different domain geometries. In other words, once the basic methods are implemented in the computational domain, we only have to adapt the mappings to the physical domain to study bounded, unbounded domains and combinations of both in 2D. It should also be emphasised that we have already successfully applied (3) in several physical settings, see Refs. [22, 23, 24, 51, 52, 73, 74, 75, 76].

We distinguish between algebraic, logarithmic and exponential maps from the computational domain $[-1,1]$ to the semi-infinite and infinite domains $[0,\infty]$ and $[-\infty,\infty]$, respectively. These have been compared extensively in [9], as well as in the earlier studies [28, 8]. In a minimalist approach, we choose the algebraic maps taking advantage of their robustness. In the direction normal and parallel to the wall, i.e. for a semi-infinite and infinite physical domains, the algebraic maps are

[TABLE]

and

[TABLE]

respectively. $L_{1,2}$ are mapping parameters linked to the length scale of the physical domain. In particular, $50\%$ of the discretisation points of the computational domain $[-1,1]$ are mapped onto the intervals $[-L_{1},L_{1}]$ and $[0,L_{2}]$ in the physical domain, respectively. One further advantage of using algebraic maps is that they retain their optimal properties as $N\to\infty$ for constant mapping parameters $L$ [9]. In Fig. 1, the low sensitivity of the algebraic map $\mathcal{A}_{1}$ with respect to $L_{1}$ is demonstrated by varying the number of collocation points, and comparing an analytical expression for the integral of a test function with the numerical result (see also Ref. [8]). Based on this, the parameters $L_{1},L_{2}$ are chosen empirically by taking into consideration the expected typical length scales of the density profile.

3.3 Computation of convolutions

We compute the convolution of the density profile $n$ and a weight function $\chi$ in real space. If the density profile is represented on a bounded domain by a relatively fine mesh, then the weight function $\chi$ may be represented on the same grid as the density distribution $n$ and the computation of Eq. (28) is straightforward.

Typically, however, we wish to compute the convolution of a relatively narrow function $\chi$, with support of the order of the particle diameter, with a density distribution $n$ which varies smoothly over a large domain. These density distributions can be accurately represented with a relatively small number of collocation points on the full domain. However, if one tries to represent the weight function $\chi$, shifted by the position of the collocation point $y_{k}$, on the same set of collocation points as the density profile, the result can be highly inaccurate. To illustrate this behavior, we plot in Fig. 2 the interpolation of three shifted Gaussian distributions on the same grid centred at $y=0$. It can clearly be seen how the quality of the representation rapidly deteriorates as the Gaussian is shifted to positions at which the grid becomes coarse.

In order to account for this behavior, we rewrite the convolution (28) as

[TABLE]

with the shifted density distribution

[TABLE]

and where $\mathcal{I}(y)$ is the intersection between the support of the weight function $\chi$ and the support of the shifted density distribution $n_{y}$, similar to the procedure employed in Ref. [16] in 1D. Assuming that the density profile changes slowly, we conclude that the behavior of the integrand in (46) is dominated by the weight $\chi$. We therefore discretise the domain $\mathcal{I}(y)$ for each collocation point $y=y_{k}$ based on the behavior of $\chi$, and interpolate the shifted density profile $n_{y}$ onto this grid.

In particular, we construct a convolution matrix ${\bf C}_{\chi}$ which acts on the vector $n_{k}$ of values of $n$ at the collocation points $y_{k}$ and returns $\left(\chi\ast n\right)_{k}$ at the collocation points $y_{k}$. The $k$-th row of ${\bf C}_{\chi}$, $\left({\bf C}_{\chi}\right)_{k}$, is then computed as follows:

Discretise $\mathcal{I}(y_{k})$ with $M$ collocation points $\hat{y}_{M}$ and compute the integration weights $\hat{w}_{M}$ for this domain. [row vector, length $M$]. 2. 2.

Evaluate the weight function $\chi_{M}=\chi(-\hat{y}_{M})$ on the discretisation points of $\mathcal{I}(y_{k})$ [column vector, length $M$]. 3. 3.

Using (36), compute the interpolation matrix ${\bf IP}_{k}$ of the original grid $y_{N}$ onto the grid of $\mathcal{I}(y_{k})$, given by $\hat{y}_{M}$, and computed in step 1. [$M\times N$ matrix]. 4. 4.

Set $\left({\bf C}_{\chi}\right)_{k}=\hat{w}_{M}\text{diag}\left(\chi_{M}\right){\bf IP}_{k}$.

Here, $\text{diag}(v)$ is the square matrix with the elements of the vector $v$ on the diagonal and zeros elsewhere.

This procedure apparently has a much higher complexity — determined by the computation of the interpolation matrices — of $O\left(MN^{2}\right)$, if compared to the computation of the convolution in Fourier space, which is of $O(N\log N)$ complexity. We note that this computation of the interpolation matrices is done point-by-point and is therefore easily parallelisable.

Also, due to its flexibility in describing the domain $\mathcal{I}(y_{k})$, the procedure described here naturally enables the accurate description of weight functions and density profiles with discontinuities and finite support. This allows us to use a smaller number of discretisation points for both the discretisation of the density profile and of $\mathcal{I}(y_{k})$.

Another way to bring down the complexity of the computation of the interpolation matrices, is to only consider points that are close to $y_{k}$. This use of a cutoff for the computation of the interpolation matrices does speed up computations, but the speed-up is not expected to be as dramatic because of the clustering of the points near high-gradient regions. Also, the introduction of such a cutoff will affect convergence properties and therefore has to be studied carefully. For this reason, we opted not to utilise such a cutoff in this work.

It is noteworthy that the procedure only has to be performed once for each geometry. In particular, in the process of solving the Euler-Lagrange equation (2), computing the convolution reduces to a matrix-vector multiplication of complexity $O(N^{2})$. We emphasise that a direct comparison of numerical procedures by means of their complexity is meaningful only if the number of points for both procedures is sufficiently large, of the same order of magnitude, and the results are of the same accuracy. We show in the following that with the real-space convolution procedure described here, a higher accuracy can be achieved while employing a considerably lower number of discretisation points, hence trading off its higher complexity.

3.4 Pseudospectral methods in 2D domains

The procedure is now generalised to 2D by appropriate application of tensor products [69], such as done e.g. in [38] to solve the Navier–Stokes (NS) equations or in [39] to solve 2D biharmonic boundary value problems. In particular, the vector space of 2D polynomials approximating a function $f\left(\bf x\right)$ defined on the unit cell $[-1,1]\times[-1,1]$, which here is the computational domain in 2D, is defined as the tensor product of the 1D polynomial vector spaces ${\bf P}_{N_{1}}\otimes{\bf P}_{N_{2}}$. Here, $N_{1},N_{2}$ are the number of collocation points along the first and second axes, respectively. The unit cell $[-1,1]\times[-1,1]$ is discretised with sets of 1D collocation points ${\bf x}_{1,2}\in[-1,1]^{N_{1,2}}$ in each direction, defined in Eq. (35) [69]. Details are given in Appendix B.

In Sec. 3.2, we have mapped the 1D computational domain $[-1,1]$ onto the whole line $[-\infty,\infty]$ and a semi-infinite space $[0,\infty]$, employing the algebraic maps $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$, respectively. Here, we proceed analogously for 2D and map the computational domain onto the physical domain through a general bijective map

[TABLE]

Here, $\mathcal{Y}_{1,2}$ are scalar functions which correspond, depending on the physical domain we wish to discretize, to $\mathcal{A}_{1},\mathcal{A}_{2}$ or to a simple scaling function. Depending on the geometry, the first and the second variables of the computational domain may be discretised employing either Gauss-Lobatto-Chebyshev collocation points in the interval $[-1,1]$ or an equispaced grid on $[0,1)$ for periodic functions.

The variables $(y_{1}^{\prime},y_{2}^{\prime})$ in the physical domain are then mapped to the Cartesian domain via

[TABLE]

In particular, depending on the geometry we wish to discretise, the physical domain may be represented in polar, spherical or skewed cartesian coordinates. The maps $\mathcal{C}$ for these cases are given in Sec. B.1 of the Appendix.

3.5 Domain discretisation for DFT-FMT equations in the half-space

Modeling a fluid with FMT requires the discretisation not only of the density $n$, but also of the weighted densities defined through Eq. (9). If the domain of the particles is confined by a hard wall, we obtain $n|_{y_{2}<0}=0$. Considering a hard-sphere fluid with spheres of diameter one, and employing the weights defined in Eq. (10), then it becomes clear that the support of the weighted densities is $\{(y_{1},y_{2}):y_{2}>-1/2\}$. Hence, a separate discretisation of the weighted densities is necessary.

Furthermore, if the contact density of a hard-sphere fluid at the wall is greater than zero, then the derivatives of the weighted densities in the direction normal to the wall, $\frac{\partial n_{\alpha}}{\partial y}$ exhibit a jump at $y_{2}=1/2$. This can best be seen for the behavior of $n_{2}$ and the second component of $\bf n_{2}$ in the computation of a hard-sphere fluid in contact with the hard wall, such as depicted in Fig. 6. There, the weighted densities $n_{2}$ and $\bf n_{2}$ are clearly non-differentiable at $y_{2}=1/2$.

In order to correctly describe this behavior, we split the half-space $\{(y_{1},y_{2}):y_{2}>-1/2\}$ used to describe the weighted densities, into the infinite slit $\{(y_{1},y_{2}):y_{2}\in[-1/2,1/2]\}$ and the half-space $\{(y_{1},y_{2}):y_{2}>1/2\}$, which are then discretised separately with pseudospectral methods. In particular, the slit is described by

[TABLE]

and the half-space is described by Eq. (84), but shifted by $1/2$ in the positive $y_{2}$-direction, such as shown in Fig. 3. In the computations which we present here, the number of collocation points in this strip, $N_{\text{strip}}$, is chosen empirically to be the next even integer to $N_{2}/3$. This value was found to provide an accurate representation of the density and the weighted densities in the vicinity of the wall.

In order to compute the convolutions as described in Sec. 3.3, one also has to discretise the intersections between the support of the weight function and the support of the density distribution. In Fig. 4, examples for such discretisations are given for intersections of the support of attractive particle–particle interactions and the half-space. A detailed description of the discretization of all geometries needed for the computation of the DFT-FMT model equations in the half-space is given in Sec. B.1.

The methods described here can also be used to describe a slit domain between two parallel walls. Similarly, one could generalize the method to a domain bounded by a nonplanar, but smooth surface, by adapting the maps from the computational to the physical domain. However, for more complex geometries, e.g. involving corners, the representation of the weighted densities and the representation of the overlap between the sphere and the wall pose considerable challenges. In this case, one way forward is to develop a spectral element-like approach, where the domain is split into subdomains. Local complexities of the geometry, where the weighted densities are known not to be smooth, are then accounted for by robust and simple discretization schemes, whilst the large part of the domain in which the density and the weighted densities are smooth, are still spectrally represented. Alternatively, the domain may be split into smaller subdomains, so that the solution is represented spectrally in each.

In general, the concepts described here can also be applied to the 1D case. In particular, the convolutions with weight functions with finite support will reduce to integrations over intervals, as opposed to the more complicated discretisations appearing in the 2D setting. Here, however, all computations for the cases exhibiting invariance in $y_{1}$-direction are performed by formally discretising the first computational variable $x_{1}$ with one collocation point only. This enables us to use the same numerical scheme as in 2D.

3.6 Iterative scheme

We have implemented both Newton [19, 20] and Picard [59] iterative schemes to solve the Euler-Lagrange equation (2). We note that if the density of the system is very small, then the ideal-gas contribution to (2) diverges logarithmically. For example, small densities can occur if the external potential $V_{\text{ext}}$ diverges. Therefore, in order to avoid a special treatment for the density profile at zero, we solve instead for

[TABLE]

Here, the external potential was split into a bounded time-dependent part $V_{\text{ext}}^{\prime}$ and a remainder ${\bar{V}}_{\text{ext}}$

[TABLE]

In other words, ${\bar{V}}_{\text{ext}}$ is the time-independent unbounded contribution to the external potential, such that $V_{\text{ext}}^{\prime}$ is bounded. If $V_{\text{ext}}$ is bounded, then ${\bar{V}}_{\text{ext}}$ is zero. We implemented the simple iterative scheme

[TABLE]

where $\lambda_{r}\in(0,1]$ is a relaxation parameter and the step size is defined by

[TABLE]

Here, ${\bf J}$ is the Jacobian of $\frac{\delta\Omega[n]}{\delta n({{\bf r}})}$ with respect to $z$. When solving for a DFT-FMT system with the Newton scheme, we include the weighted densities $n_{2},n_{3}$ and ${\bf n}_{2}$ as unknowns, similar to the explicit treatment of the nonlocal densities in [21]. We note that the Picard iterative scheme makes use of the functional form of $\frac{\delta\Omega[n]}{\delta n({{\bf r}})}$, which can be written as a fixed-point equation for $z$.

For the Newton iterations, we perform $O(10)$ steps with a relaxation parameter $\lambda_{r}=0.5$, and then set $\lambda_{r}=1$. For the Picard iterations, we initially set $\lambda_{r}=0.01$, and subsequently increase $\lambda_{r}$ to $0.2$. A Picard iterative step has a lower computational cost than one Newton iteration, given that in each Newton iteration, a linear system of equations has to be solved [see Eq. (54)]. This can be substantially improved by employing preconditioners [32]. The general idea is that the increased computational cost of each Newton iteration is compensated for by the smaller number of iterations needed compared to Picard iterations.

In 1D computations of a fluid in contact with a hard wall, the Newton scheme converges in $O(30)$ iterations, while the Picard scheme converges in $O(500-1000)$ iterative steps, leading to a comparable computational cost. In 2D computations of a contact line in contact with a wall, we find that the Picard iterative scheme is significantly more efficient than Newton iterations.

4 Verification and validation

The accuracy of the numerical scheme is validated in the following steps. First, we confirm that the procedure of Sec. 3.3 to compute convolutions does indeed converge. We then analyse the behavior of a hard-sphere fluid with and without long-range interactions in contact with a hard wall and confirm the accuracy of the results by checking whether the so-called thermodynamic sum rules apply. Similarly, in 2D, we present numerical results for a contact line and confirm the accuracy of the computations by checking the force balance normal to the interface. We then compare DFT results for multiple species of particles to the corresponding particle distributions computed by slice sampling (a Markov chain MC algorithm). In a final step, we perform DDFT computations of some model systems and check for the mass conservation of the systems.

4.1 Convergence test for convolutions

When solving the Euler-Langrange equation (2) with an FMT model for the hard-sphere free energy contribution and an attractive particle–particle free energy, convolutions have to be performed in various steps. In Fig. 5, the convergence of the convolution matrices ${\bf C}_{\chi}$ is tested. This is done by comparing the incremental change of ${\bf C}_{\chi}$ as the number of collocation points $M$ (of each of the computational variables used to discretise the domain $\mathcal{I}\left({\bf r}\right)$), is increased.

Let us first distinguish between the convolutions which need to be computed for the FMT part of the free energy and the convolutions for the attractive contributions, the main difference being that the support of the FMT weight functions $\omega_{\alpha}$ as defined in (10) is always finite and corresponds in 2D to a disk of diameter unity. In contrast, the support of attractive particle–particle interactions (19), which in 2D is given in Eqs. (68), (69), is a larger disk with the radius corresponding to the cutoff radius $r_{c}$, or, in the limiting case of $r_{c}=\infty$, a support corresponding to the full 2D space. In Fig. 5, we compare the convergence for a weight function $\phi_{\text{attr}}$ with cutoff radius $r_{c}=2.5$, such as used e.g. in [53, 40], and the full long-range particle–particle potential by setting $r_{c}=\infty$. As expected, the convergence for $\phi_{\text{attr}}(r_{c}=\infty)$ is not exponential. This is because the support of the weight function is unbounded, and no exponential convergence can be expected for unbounded domains other than for special cases such as exponentially decaying weight functions. In contrast, the convergence for a cutoff radius $r_{c}=2.5$ is exponential.

For the FMT part of the free energy, the weight functions $\omega_{\alpha}$ are first convolved with the density distribution $n$ to compute the weighted densities (9). Then, the weight functions are convolved with functions of the weighted densities given in Eq. (12) to compute the variation of the free energy. The main difference between the two cases is that the densities and the weighted densities are discretised using different domains, as shown in Fig. 3, and therefore the convolutions have to be computed using different intersections $\mathcal{I}$. As demonstrated in Fig. 5, the convergence is exponential for all cases.

We note that the convergence rates with respect to the number of collocation points $M$, depicted in Fig. 5, decrease with increasing number of collocation points $N$ of the main grid. This is because a higher number of modes of the density profile needs to be captured, therefore also requiring a higher value for $M$. In the following, we employ $M=20+N/4$ and $M=20+N/2$ collocation points for the computation of the convolution with the FMT weights $\omega_{\alpha}$ and the attractive potential with cutoff $r_{c}=2.5$, respectively. The prefactors $1/2$ and $1/4$ account for the different convergence rates shown in Fig. 5.

4.2 Planar films and DFT sum rules

We now proceed to compute planar wall-fluid interfaces for both pure hard-sphere fluids and hard-sphere fluids with attractive particle–particle interactions. Due to the convergence properties shown in Fig. 5, we choose for the latter case a Barker-Henderson potential $\phi_{\text{attr}}$ with radial cutoff $r_{c}=2.5$. We should note, however, that this effectively removes the long-range nature of the attractive particle–particle interactions.

In Fig. 6, results for a hard-sphere fluid in contact with a hard wall are compared with MC computations, giving a very good agreement. In the right subplots of Fig. 6, results for the weighted densities are depicted. It can be seen how the support of the weighted densities extends beyond the support of the density $n$, to $\{y_{2}\geq-0.5\}$. Also, the discontinuity of the weighted densities at $y_{2}=0.5$ can be observed.

For fluids with attractive particle–particle interactions, we employ a wall with an algebraically decaying 9-3 potential (e.g. [51, 52])

[TABLE]

which is obtained by integrating an LJ potential over the half-space $\{(y_{1},y_{2},y_{3})\in\mathbb{R}^{3}|y_{2}<1\}$ and where $\varepsilon_{\text{w}}$ is the strength of the wall-fluid interaction potential. This parameter is chosen such that the wall-liquid and the wall-vapor surface tensions are equivalent (with the surface tensions defined as the excess grand potential per unit area of the interface under consideration). All computations for simple fluids are performed at saturation chemical potential, such that bulk liquid and bulk vapor pressures are the same and hence both phases are equally stable. Hard-sphere fluids without attractive particle–particle interactions on the other hand are uniquely defined by their bulk densities and do not allow for two equally stable phases. Such fluids are always computed in contact with a hard, non-attractive wall. These restrictions are employed to reduce the number of pertinent physical parameters and to ensure comparability of the results.

The mapping parameter $L_{2}$ for the algebraic map (45) is chosen through a sensitivity study. In particular, for values $L_{2}\in[0.5,3.5]$ and for $N_{2}=100$ collocation points, the sum rule error is in the range of $5\pm 4\times 10^{-5}$ for the case of a fluid with attractions as well as for a hard-sphere fluid. This suggests a low sensitivity of the results with respect to the mapping parameter for all but very large or very small values of $L_{2}$, and we therefore choose $L_{2}=2$ for all further computations.

Let us now study the properties of a typical fluid density profile close to a hard wall. While wall-fluid density profiles do not exhibit any oscillations for very low bulk (or vapor) densities, the same is not true for higher bulk densities. In this case, we observe packing at the wall such as also depicted in Fig. 6. Figure 7 illustrates the difference between the quickly decaying density oscillations of a pure hard-sphere fluid in contact with a hard wall and the much slower decay of the oscillations if interparticle attractions and a long-range attractive wall potential are added to the free energy. In particular, the insets of Fig. 7 show the convergence of the density profiles with increasing number of collocation points.

Typically, the numerical accuracy of DFT computations is checked through sum rules, which represent statistical mechanical connections between thermodynamic properties of the system and the density profile [59]. The most prominent example is the so-called contact theorem, which relates the density at a planar hard wall with the bulk pressure of the system through

[TABLE]

As an aside, we note here that the sign of the second term on the right hand side differs from the expression given in the review article by Roth [59], but is consistent with Refs. [70, 33]. In Fig. 8, we show the convergence of the relative sum rule error with the number of collocation points.

For the fluid with interparticle attractions, the error for the wall-vapor interface is much lower than for the wall-liquid interface. This is due to the absence of oscillations in the density profile at the wall-vapor interface. Also, whilst the accuracy of the sum-rule for the wall-liquid interface is lower, the relative sum rule error converges at a comparable rate to that of a hard-sphere fluid of a similar bulk density. This result suggests that the main factor limiting the convergence at the wall are hard-sphere packing effects.

Whilst the sum rule provides an accurate check for consistency of the results with thermodynamic principles, it is limited to the behavior of the density close to the wall and therefore does not provide a convergence check for the full density profile. We therefore show in the right subplot of Fig. 8 the maximal incremental change of the density profile as the number of collocation points is increased. Evidently, the convergence is of a similar rate, but more smooth compared to the convergence of the sum rule error. In particular, the results for the Barker-Henderson wall-liquid interface converge slower than the results for a hard-sphere fluid. This can be linked to the slower decay of the oscillations for the Barker-Henderson fluid (see Fig. 7).

4.3 The liquid-vapor contact line

If a fluid with attractive particle–particle interactions is at saturation chemical potential, a liquid can be in stable contact with a vapor, forming an interface with a steep but smooth fluid density profile. When this liquid-vapor interface is brought in contact with a substrate, and the difference between the wall-liquid and the wall-vapor surface tensions is not greater than the liquid-vapor surface tension, a well-defined equilibrium contact angle is formed. This was also studied in [54] using the so-called local density approximation which approximates the repulsive contribution to the free energy by employing locally an equation of state (valid for the bulk fluid phase, but nevertheless used for the inhomogeneous case in the presence of a wall). Unlike FMT, this does not capture the oscillatory structure of the fluid in the vicinity of walls, which was later studied in Refs. [51, 52] for an equilibrium contact line using the original FMT framework introduced by Rosenfeld.

For the sake of simplicity, we focus on the case where the wall-vapor surface tension $\gamma_{\text{wv}}$ is equal to the wall-liquid surface tension $\gamma_{\text{wl}}$. In this case, the force balance parallel to the wall, known as the Young-Laplace equation, implies that the contact angle $\theta$ between the liquid-vapor interface and the substrate is $90^{\circ}$:

[TABLE]

where $\gamma_{\text{lv}}$ is the liquid-vapor surface tension and $\theta$ is the Young contact angle.

In our 2D computations, the contact angle is computed as follows. First, we define the mapping parameter $L_{1}$ of the algebraic map (44) used in (84). Knowing that the width of the liquid-vapor interface is about $4\sigma$ at temperature $k_{B}T=0.75\varepsilon$, we make use of the fact that $50\%$ of the collocation points are mapped onto the interval $[-L_{1},L_{1}]$ by the algebraic map (44). As we also aim to resolve properties such as the pressure of the fluid on the substrate, which vary at slightly larger distances from the contact line, we set $L_{1}=4$ for all further computations.

Furthermore, we enforce the density profile at distances $y_{2}>y_{2,\text{max}}$ from the substrate to be equivalent to an equilibrium liquid-vapor interface at an angle $\theta$ with respect to the substrate. This also assumes the validity of the Young-Laplace equation far from the substrate. Computationally, one can verify the validity of this assumption by gradually increasing $y_{2,\text{max}}$ and observing the behavior of the liquid-vapor interfaces for $y<y_{2,\text{max}}$. In all further computations, we set $y_{2,\text{max}}=25$.

In addition to the Young-Laplace equation, which represents the force balance across the contact line in the direction parallel to the wall, we can also verify the accuracy of our computations by checking the force balance normal to the wall. In short, the normal pressure acting from a hard attractive wall on a strip of fluid at position $y_{1}$ is also known as the disjoining pressure $\Pi$, which is defined as

[TABLE]

where $p$ is the bulk pressure. Insertion of Eq. (58) into Eq. (56) confirms that in the planar equilibrium case, the excess pressure acting on the substrate always vanishes. When the film height varies, however, the disjoining pressure does not vanish, and the integrated excess pressure has to correspond to the pulling force of the liquid-vapor interface in the direction normal to the substrate:

[TABLE]

In Fig. 9(a), we present results of 2D contact line computations for the disjoining pressure along the wall. It can be observed how the disjoining pressure approaches constant values

[TABLE]

away from the contact line, as the density profile normal to the wall converges to a planar equilibrium wall-liquid and wall-vapor interface. In agreement with the physical predictions, $\Pi_{\pm}$ vanishes as the number of collocation points is increased. In particular, the convergence of

[TABLE]

is shown with the dashed lines in Fig. 9(b).

We note, however, that the integration of the disjoining pressure in the normal force balance (59) is only well-defined if $\Pi$ vanishes for $y_{1}\to\pm\infty$. Consequently, the numerical error of $\Pi_{\pm}$ — regardless of its size — has to be accounted for when computing the normal force balance. Here, we do so by restricting the integration in the normal force balance to values of $y_{1}$ for which the absolute value of the disjoining pressure is greater than $2\Delta\Pi$: $\{y_{1}:|\Pi(y_{1})|>2\Delta\Pi\}$. The relative numerical error of the normal force balance (59) then becomes

[TABLE]

This is depicted in Fig. 9(b), where the Young contact angle is computed via the Young-Laplace equation.

4.4 Multiple species – comparison to slice sampling

Here we consider systems of a fixed number of particles of different species. The equilibrium states are computed similarly to the single-fluid case, by solving the Euler-Lagrange equation (2) for each species $k$ separately:

[TABLE]

The different species pairs $i,j$ may experience different interparticle potentials $\Phi_{\text{2D}}^{(i,j)}$. Also the external potentials $V_{\text{ext}}^{(k)}$ acting on species $k$ may differ. We demonstrate the convergence of the method with respect to the number of collocation points, and compare the results to those obtained by slice sampling. Details on the DDFT formulation for multi-species fluids as well as the physics of such fluids are given in [22].

Slice sampling [49] is a Markov chain MC method for sampling from a statistical distribution, in this case from the equilibrium configuration of the particles. A particular advantage of this method is that the results are invariant under scaling the distribution by a constant; hence one does not need to know the normalisation constant (or partition function). In particular, we consider two systems, soft particles (which repel each other with a Gaussian interparticle potential) in a hard wall box, and hard disks in an infinite planar geometry.

For the soft particles, the external potential is given by a quadratic

[TABLE]

along with a confining hard-walled box of side length 10. The system contains 20 ‘small’, 20 ‘medium’ and 20 ‘large’ particles, which represent different three different particle species. They repel each other via Gaussian potentials for each pair of species $i$, $j$:

[TABLE]

where $r$ is the interparticle distance. Here we define effective particle diameters $\alpha_{1}=0.5$, $\alpha_{2}=1$ and $\alpha_{3}=1.5$, which in turn define the length scale $\sigma_{i,j}$ of the interparticle potential via $\sigma_{i,j}=(\alpha_{i}+\alpha_{j})/2$. The density profiles for this system are given in Fig. 10, where it is also shown that the maximum relative error is about $2\%$ between DDFT and stochastic sampling for the parts of the density that have significant mass. Convergence with the number of collocation points is shown in Fig. 14 of Appendix C.

The external potential for the hard disk system consists of a weak confining potential $\bar{V}$ and three Gaussian wells, which are in different positions for the two species. Here, we set the confining background potential to

[TABLE]

and employ attractive Gaussian wells, leading to the external potential for the $k$th species

[TABLE]

Here, ${\bf y}=\left(y_{1},y_{2}\right)$ and ${\bf y}^{(k)}_{\ell}$ represents the position of the $\ell$-th external potential well for species $k$. For the first species we set ${\bf y}^{(1)}_{\{1,2,3\}}=\left\{\left(-1,0\right),\left(0.75,0\right),\left(0,0.75\right)\right\}$ and for the second species ${\bf y}^{(2)}_{\{1,2,3\}}=\left\{\left(-0.75,-0.25\right),\left(1,-0.25\right),\left(0,0.75\right)\right\}$. The system contains 10 particles from each of the two species. The FMT calculations are performed using 2D FMT for hard disks [61] (see Sec. 2.2). Computational results for the density profiles are depicted in Fig. 11, and convergence with the number of collocation points is shown in Fig. 15 of Appendix C.

In Fig. 12, we consider the convergence of the method by computing the relative $L^{2}$ error of the density $n$ for pairs of computations with $N$ and $N+2$ collocation points, with all other parameters fixed. This is done by interpolating each result onto the same equispaced grid of $100\times 100$ points and then computing the $L^{2}$ norm of the difference, $\|n_{N-2}-n_{N}\|_{2}$, and normalising by the norm of the most accurate computation $\|n(N_{\max})\|$ (which, as shown in Figs. 10 and 11, is virtually indistinguishable from the slice sampled ‘exact’ result). We plot the maximum error over all species. For the box calculation, the equispaced grid covers the whole box, whilst for the infinite domain it is restricted to $[-2,2]\times[-2,2]$, outside of which the densities are very small. It is clear that the method exhibits the expected spectral accuracy.

For the Gaussian case we use $M=20\times 20$ points to compute the convolution matrices accounting for the mean field contributions. For the FMT case we use $L=4$ in both directions for the infinite space, 10 collocation points for the circle FMT contributions and $10\times 10$ points for the disk contributions. The results show very good agreement with those obtained by doubling the number of points for the mean field and FMT computations, and the computations show very little dependence on $L$ in the region around the chosen values.

As mentioned, Figs. 10 and 11 compare the DFT computations with the slice sampled data. Slice sampling is performed using MATLAB’s ‘slicesample’ routine, with 5,000,000 samples; the resulting densities are then histogrammed into $50\times 50$ boxes. Evidently the agreement is very good. The slight deviations for higher densities in the FMT computations are to be expected as FMT is less accurate in the high-density regime.

4.5 Mass conservation of DDFT computations

In this section, we demonstrate the effectiveness of our numerical scheme in tackling the DDFT equations introduced in Sec. 2.4. In particular, we consider a simple scenario of the density changing between two different equilibrium configurations. The external potential is given by

[TABLE]

where $V_{\text{ext}}^{\text{BH}}$ is defined in (55).

Figure 13 shows the performance of our scheme as the number of collocation points is increased, for both the overdamped and inertial systems introduced in Sec. 2.4. We note that due to the unbounded nature of the geometry, the mass of the full system is infinite. In particular, local changes at the contact line lead to an influx of mass from the reservoir, thus increasing the net mass within a sub-domain near the contact line (see e.g. inset of left subplot of Fig. 13). For any finite domain, this change in mass has to be accounted for by a mass influx through its boundaries. In the left subplot of Fig. 13, this fundamental principle of mass conservation is used to validate the numerical scheme.

The equations are solved using the implicit variable-time-step solver for stiff problems ‘ode15s’ from MATLAB. It is noteworthy that the inclusion of inertia triples the number of unknowns in the problem as in addition to the density the two components of the velocity must be evaluated.

5 Conclusion

Historically, the main focus of computational fluid dynamics has been on analysing and finding numerical solutions to the NS equations and their simplifications. At the nanoscale, however, these equations do not apply as interparticle interactions, described by MD-MC, dominate. Unfortunately, despite drastic improvements in computational power, MD-MC simulations are only applicable for small fluid volumes. A compromise between the two is classical DFT which retains the microscopic details of macroscopic systems but at a cost significantly lower to that of macroscopic equations. Typically, DFT relies on the formulation of a free energy functional which is minimal at equilibrium, resulting in an integral equation for the (continuous) fluid density. The integral terms are due to the modelling of van der Waals attractive forces and hard-sphere effects via e.g. FMT [71]; the advantage of FMT being that it captures the oscillatory behavior of the fluid density in the vicinity of walls.

Typically, the corresponding DFT-FMT equations are solved using either Fourier methods for the convolutions [59], or real-space computations of the convolution grids [21, 1]. These approaches, however, formally restrict the applicability of the associated numerical schemes to bounded domains, where Fourier methods also assume periodicity of the solution. Here, we introduce an alternative numerical method based on a pseudospectral scheme with the collocation points mapped onto unbounded domains [9].

The resulting scheme is highly efficient, accurate and fast, and allows for the resolution of density profiles and nanoscale fluid behavior with a comparably lower number of collocation points. The convolutions required in DFT-FMT computations are performed in real space, which allows for a pre-computation of the convolution matrices. This process must be done only once for every geometry and can easily be parallelized.

We apply our scheme to four different DFT scenarios. First, we compute the density profile of a hard-sphere fluid with and without attractions close to a hard wall, and demonstrate convergence to the so-called thermodynamic sum rules [59] with increasing number of collocation points. We then apply the numerical scheme to the computation of an equilibrium contact angle and, in addition to the Young-Laplace equation, i.e. the tangential-stress balance at the contact line, we also confirm the normal-force balance. Subsequently, the scheme is applied to multi-species systems where it is successfully validated against Monte-Carlo computations. Finally, we present dynamic computations for a contact line, by using the DDFT framework developed in [24] and we confirm the accuracy of our scheme by testing for mass conservation.

We believe that the versatility of the outlined numerical methodology allows for the treatment of more involved geometries, such as capillaries and grooves, but also the inclusion of other nanoscale effects such as hydrodynamic interactions [23], which are neglected in the DDFT equations adopted here. Furthermore, of particular interest would be the development of a multiscale algorithm to bridge the nanoscale, described here by integral equations, with the macroscale described by NS. We shall address these and related questions in future studies.

Acknowledgements

We acknowledge financial support from Imperial College (IC) through a DTG International Studentship and from the Engineering and Physical Sciences Research Council (EPSRC) of the UK through Grants No. EP/L027186, EP/L025159 and EP/L020564 as well as an EPSRC-IC Pathways to Impact-Impact Acceleration Award, Grant No. EP/K503733.

Appendix A Convolution weights for 2D geometries

For a system of hard spheres, the DFT formulation includes at several instances convolutions in 3D: In the definition of the weighted densities $n_{\alpha}$ in (9), as well as in the hard-sphere and the attractive contributions to the free energy (12) and (20), respectively. If the system is invariant in one or more directions, these 3D convolutions may be rewritten as lower-dimensional convolutions. In particular, if the system is invariant in one direction in Cartesian coordinates, then a convolution $\left(n\ast\phi\right)_{\text{3D}}$ can be rewritten as a 2D convolution $\left(n\ast\Phi_{\text{2D}}\right)_{\text{2D}}$ in the $y_{1}$-$y_{2}$-plane with

[TABLE]

The corresponding 2D-weight function for the attractive contribution to the free energy as defined in (19) is for $r<1$

[TABLE]

and for $1\leq r<r_{c}$, with cutoff-length $r_{c}$

[TABLE]

Beyond the cutoff length, $r\geq r_{c}$, we set $\Phi_{\text{2D}}^{\text{attr}}(r)=0$. We note that for small values of $r$, we evaluate a Taylor expansion of $\Phi_{2D}^{\text{attr}}$ at zero for computational purposes. In dimensionless form, (24) becomes

[TABLE]

hence linking the dimensionless quantity $\frac{\varepsilon_{\text{D}}}{\varepsilon}$ to the dimensionless cutoff-radius $r_{c}$ through

[TABLE]

Finally, the convolution with FMT weight $\omega_{3,i}$ defined in (10) can be rewritten in 2D as a convolution with

[TABLE]

where $r=\sqrt{y_{1}^{2}+y_{2}^{2}}$.

Appendix B Implementation of pseudospectral methods in 2D

The collocation points ${\bf x}$ and the function $f({\bf x})$ in the 2D computational domain are represented in a block-structure as

[TABLE]

where $x_{j}^{i}$ denotes the $i$th collocation point of coordinate $j$. Any linear operator ${\bf A}=\left(a\right)_{ij}$ acting on $x_{1}$ can be written in 2D as a tensor product

[TABLE]

which can then be applied on the function vector (73). Similarly, an operator acting on $x_{2}$ can be written as

[TABLE]

This allows us to easily define linear interpolation, integration and differentiation operators in 2D.

B.1 Representation of the physical domain

The physical domain is mapped onto the Cartesian domain via $\mathcal{C}$ given in Eq. (49). In this case, the integration vector for the physical domain is computed through

[TABLE]

where $w_{N_{1,2}}$ are the 1D integration weights given in (43), and $\text{det}\left(\bf J\right)=\frac{\partial\mathcal{Y}_{1}}{\partial x_{1}}\frac{\partial\mathcal{Y}_{2}}{\partial x_{2}}-\frac{\partial\mathcal{Y}_{1}}{\partial x_{2}}\frac{\partial\mathcal{Y}_{2}}{\partial x_{1}}$ is the vector of the determinants of the Jacobian at the collocation points ${\bf x}$.

We distinguish between the following three cases:

Skewed Cartesian coordinates, by angle $\alpha$. In this case,

[TABLE]

and the integration weight vector has to be rescaled by $\sin\alpha$:

[TABLE]

This map has been applied when computing the density profile in the vicinity of a contact line [51, 52]. 2. 2.

Polar coordinates:

[TABLE]

$\varphi\in[0,2\pi)$ and $r\in[0,\infty]$ are the angular and the radial variables, respectively. The integration vector has to be multiplied with $r$:

[TABLE] 3. 3.

Projected spherical coordinates with radial parameter $R$:

[TABLE]

where $\theta\in[0,\pi]$ and $\psi\in[0,2\pi]$ are angular variables. We distinguish two cases: First, for an integration over the 2D disk defined by $\mathcal{C}_{\text{sph}}$, the integration weight vector has to be multiplied with the determinant of the Jacobian

[TABLE]

If, however, the map $\mathcal{C}_{\text{sph}}$ is complemented by a third cartesian component $y_{3}=R\sin\psi\sin\theta$, and the integration is done over the surface of the sphere with radius $R$ whilst assuming invariance of the density profile in $y_{3}$-direction, then the integration weight vector is multiplied with the classical weight for spherical coordinates

[TABLE]

B.2 Domain discretisation

For the computation of the convolutions introduced in Sec. 3.3, the intersections between the supports of the weight functions and of the density profile have to be discretised. Here, we present the maps from the computational domain to the physical domain needed for the discretisation of the different intersections.

First, we are interested in the density distribution of a fluid in contact with a planar hard wall. The domain which needs to be discretised is the half-space $H$. We therefore choose to map the domain $[-1,1]^{2}$, discretised by Gauss-Lobatto-Chebyshev collocation points in 2D, onto the half-space, using Eq. (48) and employing the algebraic maps defined in Eqs. (44) and (45) for the first and second variables $y_{1}^{\prime}$ and $y_{2}^{\prime}$, respectively:

[TABLE]

When computing the convolutions in real-space, the intersection $\mathcal{I}(y)$ between the support of the weight function $\chi$ and the shifted half-space needs to be discretised separately. For the weights $\Phi_{\text{2D}}^{\text{attr}}$, $\Phi_{\text{2D}}^{\omega_{3,i}}$ and $\omega_{2,i}$, given in Eqs. (68)-(69), (72), and (10), respectively, we distinguish four cases:

1. (1)

The support of the weight function corresponds to the surface of a sphere of radius $R$ in 3D: $A_{1}=\{{\bf r}\in\mathbb{R}^{3}:|{\bf r}|=R\}$, and the density under consideration is invariant in 1D. This is a special case of (2), and is applicable for weights $\omega_{2,i}$ and ${\boldsymbol{\omega}}_{2,i}$ given in Eqs. (10). 2. (2)

The support of the weight function corresponds to a disk of radius $R$ in 2D: $A_{2}=\{{\bf r}\in\mathbb{R}^{2}:|{\bf r}|<R\}$. This is used for weight $\omega_{3,i}$ in Eq. (10) if the density is invariant in one direction, and for the inner contribution to the 2D projection of the attractive particle–particle interaction (see Eqs. (68)). 3. (3)

The support of the weight function corresponds to an infinite annulus with inner radius $R$: $A_{3}=\{{\bf r}\in\mathbb{R}^{2}:|{\bf r}|>R\}$. This is applicable for the outer contribution to the the 2D projection of the attractive particle–particle interaction without radial cutoff, $\Phi_{\text{2D}}^{\text{attr}}$ (see Eqs. (69)). 4. (4)

The support of the weight function corresponds to a finite annulus with inner radius $R_{\text{in}}$ and outer radius $R_{\text{out}}$: $A_{4}=\{{\bf r}\in\mathbb{R}^{2}:R_{\text{in}}\leq|{\bf r}|\leq R_{\text{out}}\}$. This is applicable for the outer contribution to the the 2D projection of the attractive particle–particle interaction with a finite cutoff radius $r_{c}$ (see Eqs. (69)).

In Tabs. LABEL:tab:FourierShapes - LABEL:tab:ShapeDiscretization, we present discretisations of eleven different geometries. In Table 4, it is then described how these geometries are assembled to represent the intersections $\mathcal{I}\left({\bf r}\right)$ between $A_{1-4}$ with the shifted half-space.

For the discretisation of a disk, a finite annulus or an infinite annulus (see Table LABEL:tab:FourierShapes), we take advantage of the periodicity in the angular variable in the polar coordinates and employ an equispaced grid to describe the second computational variable $\hat{x}_{2}\in[0,1)$, which is then mapped to the angular variable $\varphi\in[0,2\pi)$. For the discretisation of a disk with radius $R$, $\mathcal{D}_{1}\left(R\right)$, the radial variable is described by an even number of Gauss-Lobatto-Chebyshev collocation points, and mapped from $[-1,1]\to[-R,R]$. Analogous to the procedure employed in Ref. [69], this avoids the treatment of the origin, but requires the explicit implementation of the symmetry $f(r,\varphi)=f(-r,\text{mod}(\varphi+\pi,2\pi))$, hence leading to an effective restriction of the computational variable $x_{1}$ to $[0,1]$, and restoring bijectivity of the map $\mathcal{Y}_{\mathcal{D}_{1}\left(R\right)}$. Alternatively, a Radau grid may be employed, which includes only one boundary point, therefore also avoiding the explicit treatment of the origin [9, 30].

The geometries described in Table LABEL:tab:ShapeDiscretization:Spherical map the computational domain onto spherical angular coordinates. The discretisation of $\mathcal{S}_{1}$ presented there uses the periodicity in the second variable, which is useful for the computation of the weighted density $n_{2}$ (see Eq. (9)). In particular, $\mathcal{S}_{1}\left(R,\theta_{1},\theta_{2}\right)$ is the surface of a sphere with radius $R$, capped at $\theta_{1}$ and $\theta_{2}$, respectively. If the full surface of the sphere is to be described, then $\theta_{1}=0$ and $\theta_{2}=\pi$. The convolution of the density distribution with the second FMT weight $\omega_{2}$ is then written as

[TABLE]

and similarly for the vector-weight ${\boldsymbol{\omega}}_{2}$. $n_{\bf r}$ is defined in Eq. (47). Here, the 2D density profile $n$ was extended into the third dimension by assuming invariance in $y_{3}$. It is now clear that the integrand is periodic in $\psi$.

In contrast to $\mathcal{S}_{1}$, $\mathcal{S}_{2}$ describes the 2D projection of $\mathcal{S}_{1}$ onto the $y_{1}$-$y_{2}$ plane. This is used to integrate the FMT weight $\Phi_{\text{2D}}^{\omega_{3}}$ for a system invariant in one direction (see Eq. (72)). In particular, the integration is written in angular spherical coordinates as:

[TABLE]

where $R^{2}\sin^{2}\left(\theta\right)\sin\left(\psi\right)$ is the determinant of the Jacobian of $\mathcal{C}_{\text{sph}}$.

The remaining discretisations $\mathcal{Y}_{\mathcal{P}_{1-6}}$ described in Table LABEL:tab:ShapeDiscretization employ a discretisation of the computational domain with Gauss-Lobatto-Chebyshev collocation points in both directions and map it onto polar coordinates. Crucially, the radial component for the maps $\mathcal{Y}_{\mathcal{P}_{2}}$ and $\mathcal{Y}_{\mathcal{P}_{3}}$ is mapped onto an interval, whose size depends on the angular variable $\varphi$ and which extends to infinity as $\varphi$ approaches the boundary. We therefore employ a generalisation of the algebraic map $\mathcal{A}_{2}$ by defining

[TABLE]

Finally, we note that for discretisations $\mathcal{Y}_{\mathcal{D}_{2}}$ and $\mathcal{Y}_{\mathcal{P}_{1-3}}$, the parameter $L$ has to be estimated. For unbounded domains, $L$ determines the spreading rate of the collocation points. In particular, it defines the position of the red solid lines in Tabs. LABEL:tab:FourierShapes-LABEL:tab:ShapeDiscretization. Here, the density profile is assumed to change relatively slowly when compared to the weight function, and consequently $L$ is obtained by minimizing the maximum interpolation error of the weight function $\chi$ on the respective discretization.

Appendix C Multi species - convergence