\section{Introduction}\label{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 \cite{Mitchell2001,Bonn.20090527,FrinkSalinger:2002:NumericalChallenges}. 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~\cite{Ebner:1977sw,Evans}, and since then has been applied to cover a wide range of fields~\cite{Wu:2007fh}, from the modeling of atomic~\cite{RothEvansLang:2002},
 molecular~\cite{Cao_Wu:2004} and polyatomic~\cite{ChandlerMcCoy:1986,YuWu:2002} systems to electrolyte solutions~\cite{Haertel:2015:DFTDoublelayer}, and even water interfaces~\cite{Hughes:2013:Water}. 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 \cite{Evans}. This is significantly less computationally demanding compared to MD-MC, whilst retaining nanoscale properties of the fluid \cite{Wu:2007fh}.

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
\cite{Zwanzig:1954cq,TangWu:2003,Lasse:etal:2013}.  A highly accurate free
energy functional of a hard sphere fluid is provided by the
\revision{Fundamental Measure Theory (FMT) introduced by
Rosenfeld}~\cite{Rosenfeld:1989qc}. 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) \cite{Sears:2003uq,DissMerath:2008,Gross:2009:DFTSAFT,Knepley:2010,Zhou:2012three,Malijevsky:2013:CriticalPointWedgeFilling}.
This `matrix-free' approach has been employed in a variety of geometries in one- (1D), two- (2D) and three-dimensional (3D) domains (see Ref.~\cite{Sears:2003uq} for a comparison of computations in different dimensions).
A prerequisite for the use of FFT is a uniform cartesian grid~\cite{DissMerath:2008}, and accurate DFT computations therefore employ dense grids with 20~\cite{Malijevsky:2013:CriticalPointWedgeFilling} to 50~\cite{Gross:2009:DFTSAFT} 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~\cite{Knepley:2010}. This was also employed in Ref.~\cite{Zhou:2012three} 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.~\cite{DouglasFrink:2000425,FrinkSalinger:2000kx,FrinkSalinger:2002:NumericalChallenges}. While this method is robust, it is also expensive~\cite{heroux:2007vn}.
Heroux~\etal~\cite{heroux:2007vn} 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.~\cite{heroux:2007vn} and applied to hard-sphere, polymer and molecular fluids. The respective code was made publicly available with the TRAMONTO software
package~\cite{TRAMONTO:2007}.

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~\cite{Goddard:2012general,Goddard:2013Unification,Goddard:2013multi,Nold:FluidStructure:2014,nold2015nanoscale,Nold:2016:thesis}. A similar approach based on spectral methods was also used to compute the relaxation dynamics of planar films~\cite{YatsyshinSerafim:2012} and
equilibrium phase transitions in confinement~\cite{PeterPRE,PeterJChemPhys2,PeterJPhysCondMatt}.

The paper is organised as follows. Secs. \ref{sec:DFT} and \ref{sec:NumericalMethod} 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.~\ref{sec:SpectralTwoD}.
\revision{We verify our numerical method in Sec. \ref{sec:Validation} and validate} it with thermodynamic sum-rules in 1D and 2D settings, for single-fluid and multiple-species equilibrium and dynamic settings. \revision{Furthermore, we also include a} comparison with stochastic sampling techniques.
We conclude in Sec. \ref{sec:Conclusion} with the outlook and potential for further investigations.

\section{Model Equations \label{sec:DFT}}

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})$~\cite{Evans,Wu-DFT,Mermin:1965fk}. 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
\begin{align}
\Omega[n] = \mathcal{F}[n] + \int n(\pos) \klammCurl{V_{\text{ext}}(\pos) - \mu} \text{d}\pos,  \label{eq:GrandPotential}
\end{align}
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
\begin{align}
\frac{\delta \Omega[n]}{\delta n({\pos})} = 0.\label{eq:EulerLagrangeEquation}
\end{align}
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}}$
\begin{align}
\mathcal{F} = \mathcal{F}_{\text{id}} + \mathcal{F}_{\text{exc}}, \label{eq:FreeEnergyFunctional}
\end{align}
with
\begin{align}
\mathcal{F}_{\text{id}}[n] = \beta^{-1} \int  n(\pos) \klamm{ \log \klamm{\Lambda^3 n(\pos) } - 1} \text{d} \pos,
\end{align}
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
\begin{align}
\phi_{\text{LJ}}(r) = 4 \varepsilon \klamm{ \klamm{\frac{\sigma}{r}}^{12} - \klamm{\frac{\sigma}{r}}^{6} }, \label{eq:LennardJonesPotential}
\end{align}
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~\cite{Zwanzig:1954cq,TangWu:2003,Lasse:etal:2013}:
\begin{align}
\mathcal{F}_{\text{exc}} = \mathcal{F}_{\text{HS}} + \mathcal{F}_{\text{attr}}.
\end{align}
In \refe{eq:LennardJonesPotential}, $\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 \revision{FMT, first developed in
1989 by Rosenfeld}~\cite{Rosenfeld:1989qc,Roth:2010fk}, which is based on the geometrical properties of hard spheres~\cite{RosenfeldObituary:2002,RothEvansLang:2002}. Improvements to this original formulation have
been presented in order to capture the freezing
transition~\cite{rosenfeld1996dimensional,RosenfeldTarazona:1997,Tarazona:2000:Freezing}, or a more accurate equation of state~\cite{Lang:2001:thesis,RothEvansLang:2002,YuWu:2002:Modified}. These
modified formulations, however, lead to equations with a functional form
close to those of the original formulation, which is still widely used
\cite{Malijevsky:2013:CriticalPointWedgeFilling,Gor:2012:QuenchedSolid,ArcherEvans:2013,Knepley:2010,RosenfeldObituary:2002}
and gives very accurate results if compared with MD or MC computations
\cite{Roth:2010fk}. In the following, for the sake of simplicity we restrict ourselves to the implementation of the original \revision{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.

\subsection{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
\begin{align}
\mathcal{F}_{\text{HS}}[\{n^{(i)}\}] &= \beta^{-1} \int \Upsilon_{\text{HS}} \klamm{\klammCurl{n_\alpha\klamm{\pos}}} \text{d} \pos, \label{eq:DefRFMT_HardSphereContribution1}
\end{align}
with
\begin{align}
\Upsilon_{\text{HS}} \klamm{\{n_\alpha\}} &=
- n_0 \ln \klamm{1-n_3} + \frac{n_1 n_2}{1-n_3} + \frac{n_2^3}{24\pi \klamm{1-n_3}^2}
- \frac{{\bf n}_1\cdot{\bf n}_2}{1-n_3} - \frac{n_2\klamm{{\bf n}_2 \cdot {\bf n}_2}}{8\pi \klamm{1-n_3}^2}
. \label{eq:DefRFMT_HardSphereContribution2}
\end{align}
The weighted densities $n_\alpha$ are defined through convolutions of weight functions $\omega_{\alpha,i}$ with the density field $n$:
\begin{align}
n_\alpha \klamm{\bf r} =
\sum_i \klamm{n^{(i)} \ast \omega_{\alpha,i}}\klamm{\pos}
= \sum_i \int n^{(i)}\klamm{\pos '} \omega_{\alpha,i} \klamm{\pos - \pos'} \text{d} \pos'.\label{eq:nAlphaDef}
\end{align}
The vector-weighted densities ${\bf n}_\alpha$ are defined analogous to Eq.~\refe{eq:nAlphaDef}, with vector weights ${\boldsymbol \omega}_{\alpha,i}$. The weights, as obtained through the decomposition of the Mayer function~\cite{Rosenfeld:1989qc}, are
\begin{align}
\omega_{2,i}\klamm{\bf r} = \delta\klamm{R_i - |\pos|},\qquad
\omega_{3,i}\klamm{\bf r} = \Theta\klamm{R_i - |\pos|},\qquad
{\boldsymbol \omega}_{2,i}\klamm{\bf r} = \frac{\pos}{|\pos|} \delta\klamm{ R_i - |\pos| }, \label{eq:WeightedDensitiesRFMT2}
\end{align}
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 (\ref{eq:WeightedDensitiesRFMT2}) and are defined by
\begin{align}
\omega_{1,i}\klamm{\bf r} = \frac{\omega_{2,i}\klamm{\bf r}}{4\pi R_i} ,\qquad
\omega_{0,i}\klamm{\bf r} = \frac{\omega_{2,i}\klamm{\bf r}}{4 \pi R_i^2} ,\qquad
{\boldsymbol \omega}_{1,i}\klamm{\bf r} = \frac{{\boldsymbol \omega}_2\klamm{\bf r}}{4 \pi R_i}.
\end{align}
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~\cite{DissMerath:2008}.
Insertion of Eq.~(\ref{eq:DefRFMT_HardSphereContribution1}) into the
Euler-Lagrange equation \refe{eq:EulerLagrangeEquation}, leads to a
contribution due to hard-sphere effects of
\begin{align}
\frac{\delta \mathcal{F}_{\text{HS}}[n]}{\delta n^{(i)}(\pos)} = \beta^{-1} \sum_{\alpha} \klamm{\diff{\Upsilon_{\text{HS}}}{n_{\alpha}} \ast \hat \omega_{\alpha,i}}\klamm{\pos},
\label{eq:Rosenfeld_Variational}
\end{align}
where $\hat \omega_{\alpha,i}\klamm{\pos} = \omega_{\alpha,i}\klamm{-\pos}$.
The corresponding 2D convolution weight for a system which is invariant in one direction is given in Sec.~\ref{sec:ConvWeights2D} of the Appendix.

\subsection{Hard-disk FMT \label{sec:DiskFMT}}
An analogous approach can be followed if hard disks are considered. Here, we follow the description by Roth \etal~\cite{RothMeckeOettel2012} and define the hard-disk free energy density as
\begin{align}
\Upsilon_{\text{HD}} \klamm{\{n_\alpha\}} = &- n_0^{\text{HD}} \ln \klamm{1-n_3^{{\text{HD}}}}+\\
& \frac{1}{4\pi \klamm{1-n_3^{{\text{HD}}}}}\klamm{ \frac{19}{12} \klamm{n_2^{{\text{HD}}}}^2  - \frac{5}{12} {\bf n}_2^{(1,{\text{HD}})}\cdot{\bf n}_2^{(1,{\text{HD}})} - \frac{7}{6} {\bf n}_2^{(2,{\text{HD}})} \cdot {\bf n}_2^{(2,{\text{HD}})} }. \notag
\end{align}
The weighted densities are defined analogously to (\ref{eq:nAlphaDef}) in 2D:
\begin{align}
\omega_{2,i}\klamm{\bf r} = \delta\klamm{R_i - |\pos|}
\qquad\text{and}\qquad
\omega_{3,i}\klamm{\bf r} = \Theta\klamm{R_i - |\pos|}. \label{eq:w2w3:HD}
\end{align}
Additional weights are
\begin{align}
\omega_{0,i}^{\text{HD}}\klamm{\pos} = \frac{\omega_{2,i}\klamm{\pos}}{2\pi R_i},\qquad
{\boldsymbol \omega}_{2,i}^{(1,{\text{HD}})}\klamm{\pos} = {\bf \pos} \delta\klamm{R_i - |\pos|},\qquad
{\boldsymbol \omega}_{2,i}^{(2,{\text{HD}})}\klamm{\pos} = \klamm{{\pos} \otimes {\pos}} \delta\klamm{R_i - |\pos|}.
\end{align}
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$.

\subsection{Attractive free energy contribution}

We follow the route of treating the attractive forces as perturbations to the repulsive forces~\cite{Lasse:etal:2013},
by including them in a mean-field manner as
\begin{align}
\mathcal{F}_{\text{attr}}[n] &= \frac{1}{2 } \iint \phi_{\text{attr}}({|\pos - \pos'|}) n(\pos)n(\pos') \text{d}\pos' \text{d}\pos, \label{eq:FEattr}
\end{align}
where $\phi_{\text{attr}}$ represents the attractive particle--particle interactions, modeled here with
\begin{align}
\phi_{\text{attr}}\klamm{r} =
\varepsilon \left\{ \begin{array}{ll}
0 & \text{for } r \leq  \sigma\\
4 \klamm{
\klamm{\displaystyle\frac{\sigma}{r}}^{12} -
\klamm{\displaystyle\frac{\sigma}{r}}^{6}
} & \text{for } r > \sigma
\end{array}  \right. .\label{eq:pattr}
\end{align}
The contribution of the free energy functional $\mathcal{F}_{\text{attr}}$ to the Euler-Lagrange equation (\ref{eq:EulerLagrangeEquation}) is
\begin{align}
\frac{\delta \mathcal{F}_{\text{attr}}[n]}{\delta n({\pos})}
= \klamm{\phi_{\text{attr}} \ast n }\klamm{\pos} = \int \phi_{\text{attr}}({|\pos - \pos'|}) n(\pos') \text{d}\pos'. \label{eq:AttrContrEulerLagrange}
\end{align}
The corresponding 2D convolution weight for a system which is invariant in one direction is given in Sec.~\ref{sec:ConvWeights2D} of the Appendix.

\subsection{Dynamic DFT \label{DDFT}}

An analogous result to that underpinning DFT holds for the dynamics of many-body systems \cite{Chan:2005uq},
leading to a dynamic DFT or DDFT.  DDFTs typically take the form of a continuity equation
\[
	\partial_t n(\pos,t) = \grad{\pos} \cdot \big( \mathbf{J}([n(\pos,t)],\pos,t) \big),
\]
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~\cite{Evans,DietrichFrisch:1990}, 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 \cite{Marconi:1999ys,Marconi:2000ys} derived the first DDFT for pairwise interparticle potentials,
\begin{align}
	\partial_t n(\pos) =
	\frac{1}{\gamma m}
	\grad{\pos} \cdot \Big[ n(\pos,t) \grad{\pos}
	\frac{\delta \mathcal{F} [n(\cdot,t)]}{\delta n(\pos,t)} \Big],
\end{align}
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 (\ref{eq:FreeEnergyFunctional}), 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~\cite{Archer:2004mb}. 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~\cite{Rex:2009qf,BenDDFT,Goddard:2012general,Goddard:2013Unification}.

In order to incorporate the effect of (microscopic) inertia in DDFT, two
different approaches have been considered. The first uses a multiple
timescale technique
\cite{Marconi:2007zr,Marconi:2006zr,tarazona2008beyond,MaroniTarazona:2008:Beyond},
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~\cite{Archer.20090107},
\begin{align}
	\partial_t n(\pos,t) &= - \grad{\pos} \cdot \big( n(\pos,t) \mathbf{v}(\pos,t) \big) \\
	\partial_t  \mathbf{v}(\pos,t) &= -  \mathbf{v}(\pos,t) \cdot \grad{\pos} \mathbf{v}(\pos,t)
	- \gamma  \mathbf{v}(\pos,t)
	- \frac{1}{m} \grad{\pos} \frac{\delta \mathcal{F} [n(\pos,t)]}{\delta n(\pos,t)}, \label{eq:DDFT:MomentumEq}
\end{align}
where $\mathbf{v}$ is the velocity distribution.  This formalism has also recently
been extended to include hydrodynamic interactions~\cite{BenDDFT,Goddard:2012general,Goddard:2013Unification}.
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
\cite{Rex:2009qf,Goddard:2012general,Goddard:2013Unification}.

\subsection{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
\begin{align}
-\frac{32}{9} \pi \varepsilon_{\text{D}} = \int_{\mathbb{R}^3} \phi_{\text{attr}}\klamm{|\pos|} \text{d} \pos. \label{eq:defineEnergyScale}
\end{align}
For an attractive particle--particle interaction \refe{eq:pattr} 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
\refe{eq:defineEnergyScale} 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_BT$.
We thus introduce the following dimensionless quantities
\begin{align}
\tilde \pos &= \frac{\pos}{\sigma},\quad
\tilde T = \frac{k_B T}{\varepsilon_{\text{D}}},\quad
\tilde \phi_{\text{attr}} = \frac{\phi_{\text{attr}}}{\varepsilon_{\text{D}}},\quad
\tilde n = n \sigma^3,\notag\\
\tilde \mu &= \frac{\mu}{\varepsilon_{\text{D}}} - \frac{k_B T}{\varepsilon_{\text{D}}} \log \klamm{\frac{\Lambda^3}{\sigma^3}} \quad \text{and} \quad
\tilde V_{\text{ext}} = \frac{V_{\text{ext}}}{\varepsilon_{\text{D}}}.
\end{align}
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
\begin{align}
T_{\text{OD}} = \frac{\gamma m \sigma^2}{\varepsilon_{\text{D}}}
\qquad \text{and}\qquad
T_{\text{IN}} =\sigma \sqrt{\frac{m}{\varepsilon_{\text{D}}}},
\end{align}
respectively. The dimensionless friction coefficient in the momentum equation (\ref{eq:DDFT:MomentumEq}) then becomes
\begin{align}
\DL \gamma = \gamma \sigma \sqrt{\frac{m}{\varepsilon_{\text{D}}}}.
\end{align}
In what follows, we omit the tildes to simplify our notation.

\newcommand{\mathcal{F}}{\mathcal{F}}

\section{Numerical Method \label{sec:NumericalMethod}}

The accurate computation of convolutions
\begin{align}
(n \ast \chi)(y) = \int \chi(y- \tilde y) n(\tilde y) \text{d} \tilde y
	\label{convolutionDefn}
\end{align}
 is crucial for solving the
Euler-Lagrange equation \refe{eq:EulerLagrangeEquation} and ultimately
obtaining a reliable equilibrium density profile. Convolutions are needed to
compute the weighted densities as given in Eq.~\refe{eq:nAlphaDef}, for the
hard-sphere contribution to the free energy (see
Eq.~\refe{eq:Rosenfeld_Variational}) as well as the attractive contribution
to the free energy (see Eq.~\refe{eq:AttrContrEulerLagrange}). The obvious
choice is to employ the convolution theorem, and perform the computation of
the convolutions in Fourier space~\cite{Roth:2010fk} through
\begin{align}
n \ast \chi = \mathcal{F}^{-1}\klammCurl{ \mathcal{F}\klammCurl{n} \cdot \mathcal{F}\klammCurl{\chi}  }, \label{eq:ConvolutionFourier}
\end{align}
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.\refe{eq:WeightedDensitiesRFMT2} and \refe{eq:w2w3:HD} needs special treatment due to the short-range finite support properties of these functions~\cite{Roth:2010fk}. 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~\cite{Knepley:2010}. Alternatively, expansions of the
equilibrium density profile in terms of the local curvatures at the wall have
been proposed~\cite{konig:2005curvature}. 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.

\subsection{Chebyshev pseudospectral method}

Here we outline a real-space discretisation method for FMT-DFT employing a
pseudospectral method~\cite{Boyd:2001Cheb}. 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
\begin{align} 
p_N(x)= \sum_{n=0}^N f(x_n) P_n(x),
\label{eq:Decomposition} 
\end{align} 
where $P_n(x)$ denote the Lagrange polynomials of degree N,
given by
\begin{align}
P_n(x) =  \prod_{m=0,m\neq n}^{N}{ \frac{ x-x_m}{x_n-x_m}},
\label{eq:Legendre}
\end{align}
with the property
\begin{align}
\revision{
P_n\klamm{x_k} = \left\{ 
\begin{array}{ll}
1 & \text{for } k = n\\
0 & \text{for } k \neq n
\end{array}
\right. .}
\end{align} 
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
\begin{align}
x_{n} &= \cos t_n 
\qquad \text{with}\qquad
t_n = \frac{\pi n}{N} 
 \qquad \text{for}\qquad n = 0\ldots N, \label{eq:Spectral:ChebCollPts}
\end{align}
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.~\refe{eq:Decomposition}, 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.~\refe{eq:Legendre} is both costly and numerically unstable. These issues can be avoided by using the barycentric formula~\cite{Berrut:2004ly,Tee:2006Ad}
\begin{align}
p_N(x) = {\displaystyle\frac{{\mathlarger\sum_{k=0}^N} {\displaystyle\frac{\bar w_k}{x - x_k}f(x_k) }}{{\mathlarger\sum_{k=0}^N}  {\displaystyle\frac{\bar w_k}{x-x_k}}}}, \label{eq:NumericalMethods:BarycentricForm}
\end{align}
where ${\bar w}_k$ are the so-called barycentric weights, which, for the Gauss-Lobatto-Chebyshev points are given by~\cite{Salzer:1972}
\begin{align}
\bar w_j = \klamm{-1}^j d_j\qquad \text{with}\qquad
d_j = \left\{ 
\begin{array}{ll}
1/2 & \text{for } j\in \{0,N\} \\
1 &\text{otherwise}
\end{array}
\right. .\label{eq:WeightBaryIP}
\end{align}
The integral of $f$ over $[-1,1]$ is computed through the use of the Clenshaw-Curtis quadrature~\cite{clenshaw:1960}
\begin{align}
\int_{-1}^1 f(x) dx = \sum_{k=0}^N w_k f(x_k),\label{eq:Spectral:Clenshaw:Quadrature}
\end{align}
with weights
\begin{align}
w_j = \frac{2d_j}{N} \left\{
\begin{array}{ll}
 1 - {\mathlarger\sum_{k=1}^{(N-2)/2}} \fracd{ 2\cos\klamm{2kt_j}}{4k^2-1} - \fracd{\cos\klamm{\pi j}}{N^2-1}  &\qquad \text{for } N \text{ even}\\
 1 - {\mathlarger\sum_{k=1}^{(N-1)/2}} \fracd{ 2\cos\klamm{2kt_j}}{4k^2-1}&\qquad \text{for } N \text{ odd}
\end{array}\right..
\label{eq:Num:Cheb:ClenshawCurtisQuad1}
\end{align}

\subsection{Unbounded domains \label{sec:Spectral:UnboundedDomains}}
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~\cite{Boyd:2001Cheb,Shen2009}:
\begin{enumerate}[(1)]
\item domain truncation,
\item approximation by orthogonal systems in the unbounded domain, or
\item mapping of the unbounded domain onto a bounded domain and applying standard spectral methods.
\end{enumerate}
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.~\cite{Goddard:2013multi,Goddard:2012general,Goddard:2013Unification,Nold:FluidStructure:2014,nold2015nanoscale,YatsyshinSerafim:2012,PeterPRE,PeterJPhysCondMatt,PeterJChemPhys2}.

\begin{figure}[h]
\centering
	\includegraphics[width=8cm]{Figure1.pdf}
	\caption{\label{fig:LSensitity} Behavior of the $\log_{10}$-logarithm of the relative error for the integral of the test function
    $f(y) = (2/\pi)^{1/2} y^2 \exp\klamm{-\frac{y^2}{2}}$, given a discretisation
    of the whole line with mapped Gauss-Lobatto-Chebyshev collocation points using the algebraic map
    $\mathcal{A}_1$ in Eq.~\refe{eq:WholeLineAlgebraicMapping}, as a function of
    the number of collocation points and the mapping parameter $L_1$. The error is given in a logarithmic grayscale from white (corresponding to $10^0$) to black (corresponding to  $10^{-15})$.}
\end{figure}

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~\cite{Boyd:2001Cheb}, as well as in the earlier
studies~\cite{Grosch:1977273,Boyd:198243}. 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
\begin{align}
\mathcal{A}_1 &: [-1,1]\to [-\infty,\infty] \quad x \to L_1 \frac{x}{\sqrt{1-x^2}}, \label{eq:WholeLineAlgebraicMapping}
\end{align}
and
\begin{align}
\mathcal{A}_2 &: [-1,1]\to [0,\infty] \quad x \to L_2 \frac{1+x}{1-x}, \label{eq:SemiInfiniteAlgMapping}
\end{align}
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$~\cite{Boyd:2001Cheb}. In
Fig.~\ref{fig:LSensitity}, 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.~\cite{Boyd:198243}). 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.

\subsection{Computation of convolutions \label{S:convolution}}

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.~\refe{convolutionDefn} 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.~\ref{Fig:kernelAll} 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.

\begin{figure}[htbp]
\centering
\includegraphics{Figure2} 
\caption{Interpolation of three Gaussian distributions $\chi(y) = \frac{1}{2\pi}e^{-y^2/2}$, shifted by the position of
the $30^{\text{th}}$, $44^{\text{th}}$, and $47^{\text{th}}$ collocation
points for the left, middle and right subplot, respectively. The whole line is
discretized with $50$ Gauss-Lobatto-Chebyshev collocation points mapped from $[-1,1]$ to
$[-\infty,\infty]$ using the algebraic map $\mathcal{A}_1$ in Eq.~\refe{eq:WholeLineAlgebraicMapping} with mapping parameter $L_1 = 2$.}
\label{Fig:kernelAll}
\end{figure}

In order to account for this behavior, we rewrite the convolution \refe{convolutionDefn} as
\begin{align}
(n \ast \chi)(y) = \int_{\mathcal{I}(y)} \chi(- \tilde y) n_y(\tilde y) \text{d} \tilde y, \label{eq:ModifiedConvolutionExpression}
\end{align}
with the shifted density distribution 
\begin{align}
n_y(\tilde y) = n( y + \tilde y) \label{eq:Numerics:DefineShiftedDensity}
\end{align}
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.~\cite{elGendi:1969} in 1D. Assuming that the density profile changes slowly, we conclude that the behavior of the integrand in (\ref{eq:ModifiedConvolutionExpression}) 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 $\klamm{\chi \ast n}_k$ at the collocation points $y_k$. The $k$-th row of ${\bf C}_{\chi}$, $\klamm{{\bf C}_{\chi}}_k$, is then computed as follows:
\begin{enumerate}
\item 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$].
\item Evaluate the weight function $\chi_M = \chi(-\hat y_M)$ on the
    discretisation points of $\mathcal{I}(y_k)$ [column vector, length
    $M$].
\item Using \refe{eq:NumericalMethods:BarycentricForm}, 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].
\item Set $\klamm{{\bf C}_{\chi}}_k = \revision{\hat w_M} \text{diag}\klamm{\chi_M} {\bf IP}_k$.
\end{enumerate}
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\klamm{MN^2}$, if compared
to the computation of the convolution in Fourier space, which is of $O(N\log
N)$ complexity. \revision{We note that this computation of the interpolation matrices is
done point-by-point and is therefore easily parallelisable.}

\revision{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)$.

\revision{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.}

\revision{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
\refe{eq:EulerLagrangeEquation}, 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. 

\subsection{Pseudospectral methods in 2D domains \label{sec:SpectralTwoD}}

The procedure is now generalised to 2D by
appropriate application of tensor products~\cite{Trefethen_2000}, such as done e.g. in
\cite{Le:1992} to solve the Navier--Stokes (NS) equations or in
\cite{Mai:2007} to solve 2D biharmonic boundary value problems. In
particular, the vector space of 2D polynomials approximating a function 
$f\klamm{\bf x}$ 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.~\refe{eq:Spectral:ChebCollPts}~\cite{Trefethen_2000}. Details are given in Appendix \ref{sec:NumericalRepresentation}.

In Sec.~\ref{sec:Spectral:UnboundedDomains}, 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
\begin{align}
\mathcal{Y}: (x_1,x_2) \to (y_1',y_2') = \klamm{\mathcal{Y}_1(x_1,x_2),\mathcal{Y}_2(x_1,x_2)}. \label{eq:General2DMapFromCompDomain}
\end{align}
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',y_2')$ in the physical domain are then mapped to the Cartesian domain
via 
\begin{align}
\mathcal{C}: (y_1',y_2') \to (y_{1},y_{2}). \label{eq:DefCartMap}
\end{align}
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.~\ref{sec:App:RepPhys} of the Appendix.

\subsection{Domain discretisation for DFT-FMT equations in the half-space \label{sec:DomainDiscr}}

Modeling a fluid with FMT requires the discretisation not only of the density
$n$, but also of the weighted densities defined through
Eq.~(\ref{eq:nAlphaDef}). 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.~(\ref{eq:WeightedDensitiesRFMT2}), 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, $\diff{n_{\alpha}}{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.\ \ref{fig:Density_DFTFMT}. There, the
weighted densities $n_2$ and $\bf n_2$ are clearly
non-differentiable at $y_{2} = 1/2$.

\begin{figure}
	\centering
		\begin{subfigure}[t]{0.47\textwidth}
		\includegraphics{Figure3_a}
		\caption{Gridlines used to discretise the density distribution $n$ in the half space.}
		\end{subfigure}\quad
		\begin{subfigure}[t]{0.47\textwidth}
		\includegraphics{Figure3_b}
		\caption{Gridlines used to discretise the half-space composed by a strip of height unity and a half-space.}
		\end{subfigure}
	\caption{\label{fig:HalfSpaceGrid} Gridlines of the discretisation of density $n$ and weighted densities $n_\alpha$. In the subplots above, $20$ collocation points are employed in each direction with mapping parameters $L_1 = L_2 = 2$. $50\%$ of the collocation points in $y_1$ and $y_2$-direction are located between the two vertical red dashed lines, defined by $\mathcal{C}\klamm{\mathcal{Y}\klamm{\pm \frac{1}{\sqrt{2}},x_2}}$, and between the bottom of the half-space and the horizontal red dashed line, defined by $\mathcal{C}\klamm{\mathcal{Y}(x_1,-1}$ and $\mathcal{C}\klamm{\mathcal{Y}(x_1,0}$, respectively.}
\end{figure}

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
\begin{align}
\mathcal{Y}_{\text{slit}}(x_1,x_2) = \klamm{ \mathcal{A}_1(x_1),\frac{x_2}{2}}, \label{eq:SlitDiscretization}
\end{align}
and the half-space is described by Eq.~\refe{eq:HalfSpaceDiscretization}, but shifted by $1/2$ in the positive $y_2$-direction, such as shown in Fig.~\ref{fig:HalfSpaceGrid}.
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.~\ref{S:convolution}, one also has to discretise the intersections between the support of the weight function and the support of the density distribution. In Fig.~\ref{fig:IntersectionsExamples}, 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.~\ref{sec:App:RepPhys}.

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.

\begin{figure}
	\centering
	\begin{subfigure}[t]{0.47\textwidth}
	\includegraphics{Figure4_a}
	\caption{$\{{\bf r}: 1 \leq|\pos - \klamm{0,0.5}^T|\} \cap \{\pos:y_2 > 0\}$}
	\end{subfigure}\quad
	\begin{subfigure}[t]{0.47\textwidth}		
		\includegraphics{Figure4_b}
	\caption{$\{{\bf r}: 1 \leq|\pos - \klamm{0,1.1}^T| \} \cap \{\pos:y_2 > 0\}$}
	\end{subfigure}
	\begin{subfigure}[t]{0.47\textwidth}
		\includegraphics{Figure4_c}
		\caption{$\{{\bf r}: 1 \leq|\pos - \klamm{0,0.5}^T| \leq 2.5\} \cap \{\pos:y_2 > 0\}$}
	\end{subfigure}\quad
	\begin{subfigure}[t]{0.47\textwidth}		
		\includegraphics{Figure4_d}
		\caption{$\{{\bf r}: 1 \leq|\pos - \klamm{0,1.1}^T| \leq 2.5\} \cap \{\pos:y_2 > 0\}$}
	\end{subfigure}	
	\caption{\label{fig:IntersectionsExamples} Discretisation of intersections $\mathcal{I}$ between the support of $\Phi_{\text{2D}}^{\text{attr}}$ (see Eqs. (\ref{eq:phiattr:2D:1}), (\ref{eq:phiattr:2D:2})) with and without cutoff at $r_c = 2.5$ and the half-space, as given in Table \ref{tab:AssemblingShapesForIntersections} (see Sec.~\ref{Discretization} of the Appendix). Collocation points of different colors represent geometries which are discretised separately.}
\end{figure}

\subsection{Iterative scheme}

We have implemented both
Newton~\cite{FrinkSalinger:2000kx,DouglasFrink:2000425} and
Picard~\cite{Roth:2010fk} iterative schemes to solve the Euler-Lagrange
equation~(\ref{eq:EulerLagrangeEquation}). We note that if the density of the
system is very small, then the ideal-gas contribution to
(\ref{eq:EulerLagrangeEquation}) 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
\begin{align}
z = \beta^{-1} \log n + {\bar V}_{\text{ext}}.
\end{align}
Here, the external potential was split into a bounded time-dependent part $V_{\text{ext}}'$ and a remainder ${\bar V}_{\text{ext}}$
\begin{align}
V_{\text{ext}}\klamm{y_1,y_2,t} = V_{\text{ext}}'\klamm{y_1,y_2,t} + {\bar V}_{\text{ext}}\klamm{y_1,y_2} .
\end{align}
In other words, ${\bar V}_{\text{ext}}$ is the time-independent unbounded
contribution to the external potential, such that $V_{\text{ext}}'$ is bounded. If
$V_{\text{ext}}$ is bounded, then ${\bar V}_{\text{ext}}$ is zero. We implemented the
simple iterative scheme
\begin{align}
z_{n+1} = z_n + \lambda_r \Delta z,
\end{align}
where $\lambda_r \in (0,1]$ is a relaxation parameter and the step size is defined by
\begin{align}
{\bf J} \klamm{\Delta z}_{\text{Newton}} &= - \frac{\delta \Omega[n]}{\delta n({\pos})}
\qquad \text{and} \qquad
\klamm{\Delta z}_{\text{Picard}} = - \frac{\delta \Omega[n]}{\delta n({\pos})} \label{eq:IterativeSteps}.
\end{align}
Here, ${\bf J}$ is the Jacobian of $\frac{\delta \Omega[n]}{\delta
n({\pos})}$ 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~\cite{FrinkSalinger:2002:NumericalChallenges}. We note that the
Picard iterative scheme makes use of the functional form of $\frac{\delta
\Omega[n]}{\delta n({\pos})}$, 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.~(\ref{eq:IterativeSteps})]. This can be substantially improved by
employing preconditioners~\cite{heroux:2007vn}. 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.
\section{\revision{Verification and validation} \label{sec:Validation}}

The accuracy of the numerical scheme is validated in the following steps.
First, we confirm that the procedure of Sec.~\ref{S:convolution} 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.

\subsection{Convergence test for convolutions}
When solving the Euler-Langrange equation (\ref{eq:EulerLagrangeEquation})
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.~\ref{fig:ConvolutionMatrixError}, 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}\klamm{\pos}$), is increased.

\begin{figure}[ht]
	\centering
	\includegraphics{Figure5}	
	\caption[Validation: Convergence of the convolution operator ${\bf C}_{\chi}$]{ \label{fig:ConvolutionMatrixError}
	Convergence of convolution matrices ${\bf C}_{\chi}$, as the number of collocation points $M$ 
	used to discretise the intersections $\mathcal{I}$ is increased. For details on how 
	${\bf C}_{\chi}$ is computed, see Sec.~\ref{S:convolution}. Results are given for the particle--particle 
	interaction potential $\phi_{\text{attr}}$ without cutoff ($\ast$) and with cutoff radius $r_c = 2.5$ ($\bullet$) 
	(see Eq.~\refe{eq:pattr}), as well as for FMT weights $\omega_{2}$ ($\blacksquare$), $\omega_{3}$ ($\blacklozenge$), $\klamm{\boldsymbol \omega_{2}}_1$ ($\blacktriangledown$) and $\klamm{\boldsymbol \omega_{2}}_2$ ($\blacktriangleleft$) (see Eqs.~\refe{eq:WeightedDensitiesRFMT2}). Black and grey symbols denote the error of convolution matrices acting on functions defined on the half-spaces $\{(y_1,y_2):y_2 >0\}$
 and $\{(y_1,y_2):y_2 >-0.5\}$, respectively. The convergence is tested on the half-space introduced in Fig.~\ref{fig:HalfSpaceGrid} with $20 \times 20$ collocation points. \revision{As a guide to the eye, the blue solid line represents a steady convergence as $\sim e^{-M}$.}}
\end{figure}

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 (\ref{eq:WeightedDensitiesRFMT2}) is
always finite and corresponds in 2D to a disk of diameter unity. In contrast,
the support of attractive particle--particle interactions \refe{eq:pattr},
which in 2D is given in Eqs.~\refe{eq:phiattr:2D:1}, \refe{eq:phiattr:2D:2}, 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.~\ref{fig:ConvolutionMatrixError}, we compare the convergence for a
weight function $\phi_{\text{attr}}$ with cutoff radius $r_c = 2.5$, such as used e.g.
in~\cite{Parry:2014:WetGroove,Malijevsky:2013:CriticalPointWedgeFilling}, 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 \refe{eq:nAlphaDef}. Then, the weight functions are
convolved with functions of the weighted densities given in
Eq.~\refe{eq:Rosenfeld_Variational} 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.~\ref{fig:HalfSpaceGrid}, and therefore the convolutions have to be
computed using different intersections $\mathcal{I}$. As demonstrated in
Fig.~\ref{fig:ConvolutionMatrixError}, 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.~\ref{fig:ConvolutionMatrixError}, 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$. 
\revision{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.~\ref{fig:ConvolutionMatrixError}.}

\subsection{Planar films and DFT sum rules}

\begin{figure}
	\centering
	\includegraphics{Figure6}	
	\caption[Density profiles of a hard-sphere fluid in contact with a hard wall]{\label{fig:Density_DFTFMT} 
	The left subplot depicts the density profile of a hard-sphere fluid in contact with a hard wall, for bulk density $n_{\text{bulk}} = 0.7151$. The solid line represents results of a \revision{DFT-FMT} computation and the top and bottom dashed lines represent the values of the contact and bulk densities, respectively. Squares are MC computations by Groot \etal~\cite{Groot:1987uq}. The right subplots represent the weighted densities $n_2$, $n_3$ and the second component of the vector-weighted density ${\bf n}_2$ defined in Eq.~\refe{eq:nAlphaDef}. We employ $60$ collocation points in the direction normal to the wall.}
\end{figure}

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.~\ref{fig:ConvolutionMatrixError}, 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.~\ref{fig:Density_DFTFMT}, 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.~\ref{fig:Density_DFTFMT}, 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.

\begin{figure}
	\centering
	\begin{subfigure}[t]{0.47\textwidth}
	\includegraphics{Figure7_a}
	\caption{Hard-sphere fluid with bulk density $n_{\textnormal{b}} = 0.6220$ in contact with a hard wall.}
	\end{subfigure}
	\quad
	\begin{subfigure}[t]{0.47\textwidth}
	\includegraphics{Figure7_b}
	\caption{Barker-Henderson fluid with a cutoff radius $r_c = 2.5$ and an attractive wall (see Eq.~(\ref{eq:Vext:BH})) with parameter $\LJdepth_{\text{w}} = 0.865$ at $k_B T = 0.75 \varepsilon$.}
	\end{subfigure}
	\caption[Validation: Density profiles at the wall-fluid interface]{\label{fig:DensityProfiles_Convergence} 
	Density profiles of fluids in contact with a wall. The insets are density oscillations around the bulk density $n_{\textnormal{b}}$ as the number of discretisation points is increased from \revision{50 to 150 in steps of 25}, where darker shades of grey represent a higher number of discretisation points. \revision{The darkest lines representing the highest numbers of collocation points are indistinguishable.}}
\end{figure}

For fluids with attractive particle--particle interactions, we employ a wall
with an algebraically decaying 9-3 potential
(e.g.~\cite{Nold:FluidStructure:2014,nold2015nanoscale})
\begin{align}
V_{\text{ext}}^{\text{BH}}(y_2) &=
4 \pi \depthLJ_{\text{w}} \klamm{ \frac{1}{45}  \klamm{\frac{1}{y_2+1}}^9 - \frac{1}{6}  \klamm{\frac{1}{y_2+1}}^3}, \label{eq:Vext:BH}
\end{align}
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 $\depthLJ_{\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
(\ref{eq:SemiInfiniteAlgMapping}) 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.~\ref{fig:Density_DFTFMT}. Figure~\ref{fig:DensityProfiles_Convergence}
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.~\ref{fig:DensityProfiles_Convergence} show the convergence of
the density profiles with increasing number of collocation points.

\begin{figure}[ht]
	\centering
	\begin{subfigure}[t]{0.47\textwidth}
		\includegraphics{Figure8_a} 
		\caption{Relative sum rule error of the contact density $n_{\textnormal{C}}$.}
		\end{subfigure}		
		\quad
		\begin{subfigure}[t]{0.47\textwidth}
		\includegraphics{Figure8_b} 
		\caption{Maximum incremental change of the density profile. \label{fig:SumRuleError:b}}
		\end{subfigure}
	\caption[Validation: Convergence test employing thermodynamic sum rules]{\label{fig:SumRuleError}
	Convergence of contact density and maximum incremental change of density profile for a fluid in contact with a wall as function of the number of grid points in direction normal to the wall.
	Computations are done for:
	(1) a hard-sphere fluid with attractive interactions ($\bullet$: $\phi_{\text{attr}}\klamm{r_c = 2.5}$), in contact with an attractive wall (see Eq.~(\ref{eq:Vext:BH})), with $\LJdepth_{\text{w}} = 0.865$ and $k_B T = 0.75 \varepsilon$. Wall-liquid and wall-vapor interface computations are denoted with solid and dashed lines, respectively; (2) a hard sphere fluid in contact with a hard wall, with bulk densities corresponding to the vapor and liquid densities of the Barker Henderson fluid, $0.028$ ($\ast$) and $0.622$ ($\blacksquare$), respectively.}
\end{figure}

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~\cite{Roth:2010fk}. 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
\begin{align}
  \left.n\right|_{y_2 = 0}  = \beta p + \beta \int_0^\infty n\klamm{y_2} \frac{\text{d} V_{\text{ext}}}{\text{d} y_2} \text{d} y_2 . \label{eq:ContactTheorem}
\end{align}
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~\cite{Roth:2010fk}, but is consistent with Refs.~\cite{SwolHenderson:1989,Herring:2010vn}. In Fig.~\ref{fig:SumRuleError}, 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.~\ref{fig:SumRuleError} 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.~\ref{fig:DensityProfiles_Convergence}).

\subsection{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~\cite{Antonio2010} 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.~\cite{Nold:FluidStructure:2014,nold2015nanoscale} for an equilibrium contact line using the \revision{original FMT framework introduced by Rosenfeld}.

For the sake of simplicity, we focus on the case where the wall-vapor
surface tension $\surfaceTension_{\text{wv}}$ is equal to the wall-liquid surface tension
$\surfaceTension_{\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$:
\begin{align}
\surfaceTension_{\text{lv}} \cos \theta = \surfaceTension_{\text{wv}} - \surfaceTension_{\text{wl}},
\end{align}
where $\surfaceTension_{\text{lv}}$ is the liquid-vapor surface tension and $\theta$
is the Young contact angle.

\begin{figure}
	\centering
	 \begin{subfigure}[t]{0.47\textwidth}
	 	\includegraphics{Figure9_a} 
		\caption{\label{fig:SumRuleError2D:a}
		Disjoining pressure along the wall (see Eq.~(\ref{eq:DisjoiningPressureDefinition})). Darker shades of grey are results with a higher number of collocation points $N$.
		}
		\end{subfigure}
		\quad
	\begin{subfigure}[t]{0.47\textwidth}
	 	\includegraphics{Figure9_b} 
		\caption{\label{fig:SumRuleError2D:b}
		Relative error of the normal force balance (see Eq.~(\ref{eq:Numerics:defineDeltaDP}), solid lines), 
	and error in the planar limit, $\Delta \Pi$ (see Eq.~\refe{eq:Numerics:DeltaDP}, dashed lines) as a function of the number of collocation points $N$ in each direction.}		
		\end{subfigure}
	\caption[Validation: Normal force balance at the contact line.]{\label{fig:SumRuleError2D}
	Disjoining pressure and normal force balance at a $90^\circ$ contact angle.
	The top right inset of the right subplot shows a typical contour plot of the density profile, 
	with values $\nDensity_{\text{vap}} + \{0.05,0.5,0.95\}\klamm{\nDensity_{\text{liq}}-\nDensity_{\text{vap}}}$ represented by blue, green and red solid lines, respectively. The bottom left inset of the right subplot shows the 
	computation time of the Picard iterative scheme as a function of the number of collocation points.}
\end{figure}

In our 2D computations, the contact angle is computed as follows. First, we define the mapping parameter $L_1$ of the algebraic map (\ref{eq:WholeLineAlgebraicMapping}) used in (\ref{eq:HalfSpaceDiscretization}).
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 (\ref{eq:WholeLineAlgebraicMapping}). 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
\begin{align}
\Pi\klamm{y_1} = - \int_{0}^\infty n\klamm{y_1,y_2} \frac{\text{d} V_{\text{ext}}}{\text{d} y_2} \text{d} y_2 + \beta^{-1} n(y_1,0)- p, \label{eq:DisjoiningPressureDefinition}
\end{align}
where $p$ is the bulk pressure. Insertion of
Eq.~(\ref{eq:DisjoiningPressureDefinition}) into
Eq.~(\ref{eq:ContactTheorem}) 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:
\begin{align}
- \int_{- \infty}^\infty \Pi\klamm{y_1} \text{d} y_1 = \surfaceTension_{\text{lv}} \sin \theta.\label{eq:NormalForceBalance}
\end{align}

In Fig.~\ref{fig:SumRuleError2D: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 
\begin{align}
\Pi_\pm = \lim_{y_1 \to \pm\infty} \Pi(y_1)
\end{align}
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
\begin{align}
\Delta \Pi = \max |\Pi_\pm| \label{eq:Numerics:DeltaDP}
\end{align}
 is shown with the dashed lines in Fig.~\ref{fig:SumRuleError2D:b}. 

We note, however, that the integration of the disjoining pressure in the normal force balance 
(\ref{eq:NormalForceBalance}) 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 \refe{eq:NormalForceBalance} then becomes
\begin{align}
\tilde \Delta  = \left| 1 + \frac{ \int_{\{y_1:|\Pi(y_1)|>2\Delta \Pi\}} \Pi\klamm{y_1} \text{d} y_1 }{\surfaceTension_{\text{lv}} \sin \theta} \right|.
\label{eq:Numerics:defineDeltaDP}
\end{align}
This is depicted in Fig.~\ref{fig:SumRuleError2D:b}, where the Young contact angle is computed via the Young-Laplace equation. 

\subsection{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 (\ref{eq:EulerLagrangeEquation}) for
each species $k$ separately:
\begin{align}
\frac{\delta \Omega[\{n_k\}]}{\delta n_k({\pos})} = 0.
\end{align}
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~\cite{Goddard:2013multi}.

Slice sampling~\cite{Neal03} 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.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{Figure10.pdf}
	\caption{ \label{fig:multiSpeciesGaussian} Density profiles (top row), the absolute value of the difference between the DDFT result and that from stochastic sampling (second row) and contours (bottom row) for a system of three species interacting with Gaussian potentials (increasing size
	left--right) for 50 collocation points.}
\end{figure}

For the soft particles, the external potential is given by a quadratic
\[
	V_{\text{ext},\rm{S}}\klamm{y_1,y_2} = 2 \big( (y_1 - 8)^2 + (y_2 - 7)^2 \big)
\]
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$:
\[
	\Phi_{\text{2D},\rm{Gaussian}}^{(i,j)}(r) = 2 \exp \klamm{ -\frac{r^2}{\sigma_{i,j}^2}},
\]
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.~\ref{fig:multiSpeciesGaussian}, 
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.~\ref{fig:App:multiSpeciesGaussianN} of Appendix~\ref{sec:App:Multispecies}.

\begin{figure}
	\includegraphics[width=\textwidth]{Figure11.pdf}		
	\caption{ \label{fig:FMT2Species} 
	Comparison of density profiles (left subplot) and contours (right subplot) for a system of two hard disk species for 60 collocation points (solid lines: red, blue), virtually indistinguishable with results of slice sampling (dashed lines: green, magenta).}
\end{figure}

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
\[
	{\bar V}_{\rm{HD}}\klamm{y_1,y_2} = 0.01 (y_1^2 + y_2^2)
\]
and employ attractive Gaussian wells, leading to the external potential for the $k$th species
\[
	V_{\text{ext}}^{(k)}\klamm{y_1,y_2} = {\bar V}_{\rm{HD}}  - 3 \sum_{\ell=1}^3
	e^{-2|{\bf y} - {\bf y}^{(k)}_\ell|^2}.
\]
Here, ${\bf y} = \klamm{y_1,y_2}$ 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\{\klamm{-1,0},\klamm{0.75,0},\klamm{0,0.75}\right\}$ and for the second
species ${\bf y}^{(2)}_{\{1,2,3\}} =
\left\{\klamm{-0.75,-0.25},\klamm{1,-0.25},\klamm{0,0.75}\right\}$.
The system contains 10 particles from each of the two species. The FMT
calculations are performed using 2D FMT for hard
disks~\cite{RothMeckeOettel2012} (see Sec.~\ref{sec:DiskFMT}). Computational
results for the density profiles are depicted in Fig.~\ref{fig:FMT2Species}, and 
convergence with the number of collocation points is shown in Fig.~\ref{fig:App:FMT2SpeciesN} of Appendix~\ref{sec:App:Multispecies}.

\revision{In Fig.~\ref{fig:MultispeciesConvergence}, 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. \ref{fig:multiSpeciesGaussian} and \ref{fig:FMT2Species}, 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.~\ref{fig:multiSpeciesGaussian} and \ref{fig:FMT2Species}
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\times50$ 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.

\begin{figure}
\begin{center}
	\includegraphics[width=0.5\textwidth]{Figure12.pdf}
\end{center}
	\caption{\label{fig:MultispeciesConvergence} Maximum relative $L^2$ errors of successive pairs
	of computations for $N$ and $N+2$ collocation points for soft particles ($\blacksquare$) and FMT hard
	disks ($\bullet$).}
\end{figure}

\subsection{Mass conservation of DDFT computations}
\begin{figure}
	\centering
	 	 \begin{subfigure}[t]{0.47\textwidth}
 		\includegraphics{Figure13_a.pdf}
		\end{subfigure}\quad
	 	 \begin{subfigure}[t]{0.47\textwidth}
		\includegraphics{Figure13_b.pdf}
		\end{subfigure}\quad
	\caption{\label{fig:DynamicMaxMassError} Relative mass error per unit time in the rectangular subdomain $[-2,2]\times [0,2]$.
	Left: the maximum relative mass error per unit time $\frac{\text{d} }{\text{d} t}(m_{\text{err}}/m)$ as a function of the number of collocation points, where the inset depicts the mass in the
    subdomain as a function of time. Solid lines show results for the
    overdamped limit with $\Delta \depthLJ_{\text{w}} = 1$,  $\tau = 5$ and dashed lines are the inertial case with friction
    coefficient $\gamma = 2$. Right: contour plots for
    four different times. \revision{The blue, green and red solid curves show the locus of points in the domain for which $n = \nDensity_{\text{vap}} + i\klamm{\nDensity_{\text{liq}}-\nDensity_{\text{vap}}}$, where $i = 0.05,0.5$ and $0.95$, respectively.} All computations are performed for a Barker-Henderson-FMT fluid modelled with a cutoff potential $\phi_{\text{attr}}$ with dynamic external potential $V_1$ as defined in (\ref{eq:DynamicVext:1}).	
    The subdomain is discretized with $40 \times 40$ collocation points.}
	\end{figure}
	
\revision{In this section, we} demonstrate the effectiveness of our numerical scheme in tackling \revision{the DDFT equations introduced in Sec.~\ref{DDFT}. In particular,} we consider a simple scenario of the density changing between two different equilibrium configurations. The external potential is given by
\begin{align}
V_1(y,t;\depthLJ_{\text{w}},\Delta\depthLJ_{\text{w}},\tau) = V_{\text{ext}}^{\text{BH}}(y;\depthLJ_{\text{w}})
\left\{
\begin{array}{ll}
1 + \frac{\Delta\depthLJ_{\text{w}}}{\depthLJ_{\text{w}}} \sin\klamm{\frac{\pi t}{\tau}}& \text{for } t < \tau\\
1& \text{for } t \geq \tau
\end{array}
\right. , \label{eq:DynamicVext:1}
\end{align}
where $V_{\text{ext}}^{\text{BH}}$ is defined in (\ref{eq:Vext:BH}).

Figure~\ref{fig:DynamicMaxMassError} shows the performance of our scheme as the
number of collocation points is increased, for both the overdamped and
inertial systems introduced in Sec.~\ref{DDFT}. \revision{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.~\ref{fig:DynamicMaxMassError}). 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.~\ref{fig:DynamicMaxMassError}, 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.
\section{Conclusion \label{sec: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~\cite{Wu-DFT}; 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~\cite{Roth:2010fk}, or real-space
computations of the convolution grids~\cite{FrinkSalinger:2002:NumericalChallenges,TRAMONTO:2007}. 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~\cite{Boyd:2001Cheb}.

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~\cite{Roth:2010fk} 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~\cite{Goddard:2013Unification} 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~\cite{Goddard:2012general}, 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.
\section*{Acknowledgements}
We acknowledge financial support from Imperial College (IC) through a DTG
International Studentship \revision{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.}

\bibliographystyle{plain}
