Public interface

solve(functional::AbstractFunctional; <keyword arguments>)

Find the density profile $\rho_i(\vec{r})$ of each species $i$ from self-consistent solution of the Euler-Lagrange equations

\[\rho_i(\vec{r}) = A_i \exp\left(-\beta V_{\mathrm{ext},i}(\vec{r}) + c_{1,i}(\vec{r}; [\{\rho_i\}], T)\right)\]

of classical DFT. The first argument functional determines how to evaluate the species-resolved direct correlation functionals $c_{1,i}(\vec{r}; [\{\rho_i\}], T)$. The prefactors $A_i$ are determined per species either from the chemical potential $\mu_i$ according to $A_i = \exp(\beta \mu_i)$ or from the target mean density $\bar{\rho}_i$, in which case $A_i$ sets the normalization. Specify either μ or ρ̄ for each species.

The function returns a NamedTuple with the following keys:

• ρ: Latest state of the density profiles of all species.
• converged: Convergence indicator.
• iteration: Iteration count.
• ΔρELabsmax: Latest species-wise residues.
• Further settings and quantities for reproducibility: functional, solver, T, Vext, μ, ρ̄.

Arguments

• T: Temperature. Default: 1.0
• μ: Chemical potential of each species. Can be a number or missing. Default: missing for all species
• ρ̄: Mean density of each species. Can be a number or missing. Default: missing for all species
• Vext: External potential of each species. Must be compatible with the underlying geometry and can be given either as a function or an array. Default: 0 for all species
• ρ_init: Initial density profile of each species. Must be compatible with the underlying geometry and can be given either as a function or an array. Default: 0 for all species
• solver: Solve method. Default: PicardSolver with default keyword arguments
• tol: Value of $\Vert \Delta \rho_\mathrm{EL}(\vec{r}) \Vert_\infty$ for which convergence is declared, whereby $\Delta \rho_\mathrm{EL}(\vec{r})$ is the difference of right and left hand side of the Euler-Lagrange equation. Default: a reasonable value for the utilized floating-point precision, but no less than 1e-8
• maxiter: Maximum number of iterations. Default: 10000
• iterate_logρ: Whether to take the logarithm of the Euler-Lagrange equation and use $\log\rho(\vec{r})$ for iteration. Default: false
• boundary_to_NaN: Whether to set the values of all output arrays to NaN within regions that are controlled by the boundary conditions of the considered geometry. Default: false
• callbacks: Callbacks that are invoked in each iteration. Default: none

See also Geometries, Functionals, Solvers, Callbacks.

AndersonSolver(functional::AbstractFunctional; m::Int=10, α::Real=0.15, cond_threshold::Real=1000, rel_absmax_threshold::Real=0.1)

Construct a solver for the functional functional that uses Anderson acceleration with a history of at most m previous iterates. α controls the mixing ratio of new guess and old states similar to the corresponding setting in PicardSolver. To ensure numerical stability, the matrix of previous residue differences must have a condition number of at least cond_threshold, otherwise old iterates are successively dropped until this threshold is reached. Additionally, Anderson acceleration is only invoked if the relative deviation between $\Vert \rho(\vec{r}) \Vert_\infty$ and $\Vert \rho_\mathrm{EL}(\vec{r}) \Vert_\infty$ is less than rel_absmax_threshold (i.e. for states sufficiently close to the solution), otherwise simple mixed Picard iterates are performed.

PicardSolver(functional::AbstractFunctional; α::Real=0.05, adaptive::Bool=false, αmin::Real=0.001)

Construct a solver for the functional functional that uses Picard iteration with mixing parameter α.

One iteration step has the form

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

where $\rho_\mathrm{EL}(\vec{r})$ is determined by evaluating the right hand side of the Euler-Lagrange equation with $\rho(\vec{r})$.

If adaptive is true, the solver will reduce $α$ based on the value of $\Vert \Delta \rho_\mathrm{EL}(\vec{r}) \Vert_\infty$ (but no lower than αmin) to ensure stability. This might be useful especially during the first few iterations, but note that the adaptation mechanism is rather heuristic.

FixedPointAccelerationJLSolver(functional::AbstractFunctional; Algorithm::Symbol=:Anderson, m::Int=10, default_kwargs=(Dampening=0.15, Dampening_With_Input=true), kwargs...)

Construct a solver for the functional functional which uses the package FixedPointAcceleration to perform iterations. Algorithm chooses the algorithm to use and m sets the maximum history size of previous iterations to keep track of. The remaining keyword arguments override the sane defaults in default_kwargs and are forwarded to FixedPointAcceleration.fixed_point_new_input.

CartesianGeometry{D}(L⃗, dr⃗; boundary=:periodic)
CartesianGeometry{D}(L⃗::SVector{D,<:Real}, dr⃗::SVector{D,<:Real}; boundary::NTuple{D,Symbol}=:periodic)

Construct a D-dimensional Cartesian coordinate system with box size L⃗, discretization dr⃗ and boundary conditions boundary (if D is 1, one of :periodic, :bulk, :walls, otherwise one of :periodic, :walls). Scalar arguments are interpreted to apply equally for each Cartesian axis.

PlanarGeometry(L::Real, dz::Real; boundary=:periodic)

Construct a planar three-dimensional coordinate system with length L, discretization dz and boundary conditions boundary (one of :periodic, :bulk, :walls). All quantities are assumed to depend only on the coordinate $z$ and to be translationally invariant along the lateral directions $x$ and $y$.

SphericalGeometry(R::Real, dr::Real; boundary=:bulk)

Construct a spherical three-dimensional coordinate system from the origin to the outer radius R with discretization dr and boundary conditions boundary (one of :bulk, :wall). All quantities are assumed to depend only on the radial coordinate $r$.

SphericalShellGeometry(R_inner::Real, R_outer::Real, dr::Real; boundary=:bulk)

Construct a spherical three-dimensional coordinate system around a solute in a shell between inner radius R_inner and outer radius R_outer with discretization dr and boundary conditions boundary (one of :bulk, :walls). All quantities are assumed to depend only on the radial coordinate $r$.

coordinate_axes(geom::AbstractGeometry)

Return the coordinate axes of geom.

coordinate_labels(::Type{<:AbstractGeometry})

Return a tuple with the coordinate labels of the given geometry type.

coordinate_mesh(geom::AbstractGeometry)

Return an iterator over the coordinate mesh of geom. The iterator will visit every coordinate r⃗ in the system (including those within the boundary regions).

If only the size of the mesh is required, use coordinate_size. To obtain the individual coordinate axes, which make up the components of the mesh, see coordinate_axes.

coordinate_size(geom::AbstractGeometry)

Return the size of the coordinate mesh of geom.

total_dims(::Type{<:AbstractGeometry})
total_dims(geom::AbstractGeometry)

Return the number of dimensions of the underlying geometry.

CompositeFunctional(geom::AbstractGeometry, functionals...)

Construct a functional that is a sum of other functionals in the geometry geom.

A popular example is the FMT + mean-field functional $F_\mathrm{exc} = F_\mathrm{exc}^\mathrm{FMT} + F_\mathrm{exc}^\mathrm{MF}$.

The provided functionals should be constructors of <:AbstractFunctionals. To override the default keyword arguments, a NamedTuple can alternatively be given, which contains the constructor of the functional as the first element and suitable keyword arguments in the remaining entries.

Examples

# Pass only the type of the functional to use default keyword arguments
CompositeFunctional(geom, RosenfeldFunctional, LennardJonesMeanFieldFunctional)

# Pass a NamedTuple to provide custom keyword arguments
CompositeFunctional(geom, RosenfeldFunctional, (LennardJonesMeanFieldFunctional, rcut=3.0))

NullFunctional(geom::AbstractGeometry)

Construct a dummy functional which is identically zero. Using this functional corresponds to the system being an ideal gas (this may be useful for testing purposes).

Fexc_func(functional::AbstractFunctional)

Return a function to evaluate the excess free energy $F_\mathrm{exc}([\{\rho_i\}], T)$ for a given functional.

The returned function defines the method

Fexc(ρ::NTuple{Nspecies,<:AbstractArray}, T::Real=1.0)::Real

where Nspecies is the number of species. For the case of a single species, the following method is additionally defined:

Fexc(ρ::AbstractArray, T::Real=1.0)::Real

c1!_func(functional::AbstractFunctional)

Return a function to evaluate the species-resolved one-body direct correlation functional $c_{1,i}(\vec{r}; [\{\rho_i\}], T)$ according to the given functional.

The returned function defines the method

c1!(result::NTuple{Nspecies,<:AbstractArray}, ρ::NTuple{Nspecies,<:AbstractArray}, T::Real=1.0)

where Nspecies is the number of species. For the case of a single species, the following method is additionally defined:

c1!(result::AbstractArray, ρ::AbstractArray, T::Real=1.0)

c1_func(functional::AbstractFunctional)

Return a function to evaluate the species-resolved one-body direct correlation functional $c_{1,i}(\vec{r}; [\{\rho_i\}], T)$ according to the given functional. Contrary to c1!_func, the result is returned instead of mutated in-place.

The returned function defines the method

c1(ρ::NTuple{Nspecies,<:AbstractArray}, T::Real=1.0)::NTuple{Nspecies,<:AbstractArray}

where Nspecies is the number of species. For the case of a single species, the following method is additionally defined:

c1(ρ::AbstractArray, T::Real=1.0)::AbstractArray

geometry(functional::AbstractFunctional)

Return the underlying geometry for which the functional was constructed.

Ω_func(functional::AbstractFunctional)

Return a function to evaluate the grand potential

\[\Omega = F_\mathrm{id}([\{\rho_i\}], T) + F_\mathrm{exc}([\{\rho_i\}], T) + \sum_i \int \mathrm{d} \vec{r} \rho_i(\vec{r}) (V_{\mathrm{ext},i}(\vec{r}) - \mu_i)\]

where $F_\mathrm{id}([\{\rho_i\}], T) = k_B T \sum_i \int \mathrm{d} \vec{r} \rho_i(\vec{r}) (\log \rho_i(\vec{r}) - 1)$ and $F_\mathrm{exc}([\{\rho_i\}], T)$ is evaluated according to the given functional (see also Fexc_func).

The returned function defines the method

Ω(ρ::NTuple{Nspecies,<:AbstractArray}, Vext::NTuple{Nspecies,<:AbstractArray}, μ::NTuple{Nspecies,<:Real}, T::Real=1.0)::Real

where Nspecies is the number of species. For the case of a single species, the following method is additionally defined:

Ω(ρ::AbstractArray, Vext::AbstractArray, μ::Real, T::Real=1.0)::Real

Class of functionals that describe hard-core systems with fundamental measure theory (FMT).

The excess free energy possesses the general form

\[\beta F_\mathrm{exc}[\{\rho_i\}] = \int \mathrm{d}\vec{r} \Phi(\{n_\alpha(\vec{r})\})\]

with (scalar and vectorial) weighted densities $n_\alpha(\vec{r}) = \sum_i (\omega_{\alpha,i} * \rho_i)(\vec{r})$ that are obtained from convolutions of the density profiles $\rho_i(\vec{r})$ with appropriate weight functions $\{\omega_{\alpha,i}(\vec{r})\}$.

For constructing specific kinds of well-known FMT functionals, see PercusFunctional, RosenfeldFunctional, WhiteBearFunctional, WhiteBearMkIIFunctional.

PercusFunctional(geom::AbstractGeometry{<:Real,1}; σ=1, dΦ_mode=:analytic)

Construct the Percus functional [J. K. Percus, J. Stat. Phys. 15, 505 (1976)] for hard rods with size σ in the one-dimensional geometry geom.

RosenfeldFunctional(geom::AbstractGeometry{<:Real,3}; σ=1, dΦ_mode=:analytic)

Construct the Rosenfeld FMT functional [Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989)] for hard spheres with diameter σ in the geometry geom. Note that the current implementation fixes a sign error [S. Phan et al., Phys. Rev. E 48, 618 (1993)] in one of the terms of $\Phi(\{n_\alpha(\vec{r})\})$ in Rosenfeld's original publication.

WhiteBearFunctional(geom::AbstractGeometry{<:Real,3}; σ=1, dΦ_mode=:analytic)

Construct the White Bear FMT functional [R. Roth et. al., J. Phys. Condens. Matter 14, 12063 (2002)] for hard spheres with diameter σ in the geometry geom.

WhiteBearMkIIFunctional(geom::AbstractGeometry{<:Real,3}; σ=1, dΦ_mode=:analytic)

Construct the White Bear MkII FMT functional [H. Hansen-Goos and R. Roth, J. Phys. Condens. Matter 18, 8413 (2006)] for hard spheres with diameter σ in the geometry geom.

Class of mean-field functionals with excess free energy

\[F_\mathrm{exc}[\{\rho_i\}] = \frac{1}{2} \sum_i \sum_j \int \mathrm{d}\vec{r} \int \mathrm{d}\vec{r}' \rho_i(\vec{r}) \rho_j(\vec{r}') \phi_{\mathrm{MF},ij}(\vert\vec{r} - \vec{r}'\vert)\]

where $\phi_{\mathrm{MF},ij}(\vert\vec{r} - \vec{r}'\vert)$ is the contribution of the pair potential between species $i$ and $j$ that is treated via the mean-field/random-phase approximation. The argument σ sets the length scales (particle sizes), ϕ_MF_func are the functions in the above integrand, and pad should be chosen manually to provide sufficient padding.

See also CompositeFunctional for using a mean-field functional in combination with another functional (e.g. Fundamental measure theory (FMT) functionals for describing hard-core repulsion) and LennardJonesMeanFieldFunctional for a specialization to truncated Lennard-Jones interactions.

MeanFieldFunctional{T,Geom,Nspecies}(geom::Geom; σ::NTuple{Nspecies,<:Real}, ϕ_MF_func::NTuple{Nspecies,NTuple{Nspecies,<:Function}}, pad::Real) where {T<:Real,Geom<:AbstractGeometry{T},Nspecies}

Construct a mean-field functional for Nspecies species where the argument σ sets the length scales (particle sizes), ϕ_MF_func are functions describing the mean-field contribution of interparticle interactions, and pad should be chosen manually to provide sufficient padding.

While this constructor is suitable for generic mean-field functionals, consider using specializations for specific interaction types: LennardJonesMeanFieldFunctional.

LennardJonesMeanFieldFunctional(geom::AbstractGeometry; σ::Real=1, rcut::Real=2.5 .* σ, ϵ::Real=1)

Construct the mean-field functional for the truncated Lennard-Jones fluid in the geometry geom.

The truncated Lennard-Jones pair potential between particles of species $i$ and $j$ has the form $\phi_{ij}(r) = 4 \epsilon_{ij} ((\sigma_{ij} / r)^{-12} - (\sigma_{ij} / r)^{-6})$ for interparticle distance $r = \vert\vec{r} - \vec{r}'\vert < r_{\mathrm{cut},ij}$ and is $0$ for $r \geq r_{\mathrm{cut},ij}$. The attractive part relevant for the mean-field functional is defined as $\phi_{\mathrm{MF},ij}(r) = \phi_{ij}(r)$ for $r \geq r_{\mathrm{min},ij} = 2^{1 / 6} \sigma_{ij}$ and $\phi_{\mathrm{MF},ij}(r) = -\epsilon_{ij}$ for $r < r_{\mathrm{min},ij}$.

The arguments σ, ϵ and rcut set the respective parameters.

PrintCallback(quantities::Vector{<:String}=["iteration", "‖ΔρEL‖"]; every::Int=1, colwidth::Int=15, io::IO=stdout)

Print various quantities during iteration.

Possible elements of quantities are "iteration", "‖ΔρEL‖" (or "ΔρELabsmax"), "Ω", "μ", "ρ̄". every sets the interval in terms of the iteration count to print output. colwidth sets the column width in terms of the number of characters. io specifies the output stream to use (this is stdout by default, but could also be a writeable file stream for instance).

SleepCallback(seconds)

Sleep for a time of seconds after each iteration step.

PlotCallback(; every::Int=1, plot_kwargs...)

Show an interactive plot of the density profile throughout the iteration. every sets the interval in terms of the iteration count to update the plot. The remaining keyword arguments are forwarded to plot.