Development guide

General notes

ClassicalDFT.jl makes heavy use of Julia's type system and its multiple dispatch capabilities. In particular, various aspects of a classical DFT calculation, such as

• the geometry of the system,
• the specific density functional one wants to use,
• the numerical solution method for the self-consistent Euler-Lagrange equation,

are implemented as concrete types (or "classes" if you are used to the usual object-oriented terminology). Common functions which are required in the calculation, e.g.

• how to compute integrals and convolutions in the considered geometry,
• how to evaluate $F_\mathrm{exc}([\rho], T)$ and $c_1(\vec{r}; [\rho], T)$ for the given functional,
• how to perform an iteration step,

are specialized on these types. Julia's multiple dispatch mechanism then selects the correct implementation by choosing the method with a signature that matches as specifically as possible the types of the given arguments.

Example

Subtypes of AbstractGeometry represent the geometrical properties of a system. Concrete geometries must provide implementations for certain functions. For instance, coordinate_labels(geom) should return a tuple with labels of the different coordinate axes, which depends on the concrete type of the given geometry geom:

using ClassicalDFT

coordinate_labels(geom::CartesianGeometry{3}) = ("x", "y", "z")
coordinate_labels(geom::PlanarGeometry) = ("z",)
coordinate_labels(geom::SphericalGeometry) = ("r",)

Implementing custom functionals

Fundamental measure theory (FMT)

If your new functional is derived from fundamental measure theory (FMT), you merely need to provide a specialization of the functions Φ! and dΦ!, which calculate the FMT integrand and its derivatives with respect to all weighted densities. Refer to src/functionals/fmt.jl to see how this is done for the available FMT functionals.

Mean-field

If your functional is of mean-field type, you may use the constructor MeanFieldFunctional and provide the mean-field part $\phi_\mathrm{MF}(r)$ of your interaction potential as a function in the keyword argument ϕ_MF_func. However, be aware that, depending on the considered geometry, $\phi_\mathrm{MF}(r)$ may need to be integrated out first. Refer to src/functionals/meanfield.jl to see how this is done for the available mean-field functionals.

General functionals

Your custom functional MyFunctional{T,Geom,Nspecies} must be a concrete subtype of AbstractFunctional{T,Geom,Nspecies}, and it must hence also have the type parameters T<:Real, Geom<:AbstractGeometry{T}, and Nspecies (you are free to define additional parameters, of course). The fields geom::Geom and σ::NTuple{Nspecies,T} are required and should be set respectively to the considered geometry and the utilized length scales (usually the particle sizes). The constructor should have the signature MyFunctional(geom::AbstractGeometry; σ=1, <further keyword arguments>) and it may dispatch on concrete geometries, if necessary. The number of species Nspecies that the functional applies to should be inferred from the provided keyword arguments (e.g. from σ).

To make the functional work in a classical DFT calculation, the functions Fexc_func and c1!_func must be specialized for your concrete functional type. Make sure to specialize only the two-argument methods Fexc_func(::AbstractFunctional, ::Val{Nspecies}) and c1!_func(::AbstractFunctional, ::Val{Nspecies}) for your custom functional.

The following is a rough sketch of a functional implementing the local density approximation (LDA) for a single species, where $F_\mathrm{exc}([ρ], T) = \int \mathrm{d}\vec{r} ρ(\vec{r}) f_\mathrm{exc}(ρ(\vec{r}), T)$ with some function $f_\mathrm{exc}(\rho_b, T)$ that is associated to the bulk equation of state of the considered fluid.

using ClassicalDFT
import ClassicalDFT: AbstractGeometry, AbstractFunctional, Fexc_func, c1!_func, integrate_func

struct LDAFunctional{T<:Real,Geom<:AbstractGeometry{T},Nspecies} <: AbstractFunctional{T,Geom,Nspecies}
    geom::Geom
    σ::NTuple{Nspecies,T}
    function LDAFunctional(geom::Geom; σ=1) where {T<:Real,Geom<:AbstractGeometry{T}}
        σ isa Real || throw(ArgumentError("Nspecies > 1 left as an exercise for the reader"))
        new{T,Geom,1}(geom, (T(σ),))
    end
end

fexc(functional::LDAFunctional, ρb, T) = ρb # your fexc function

dfexc(functional::LDAFunctional, ρb, T) = 1.0 # your fexc derivative

function Fexc_func(functional::LDAFunctional, ::Val{Nspecies}) where {Nspecies}
    Nspecies == 1 || throw(ArgumentError("Nspecies > 1 left as an exercise for the reader"))
    geom = geometry(functional)
    integrate = integrate_func(geom)
    integrand = zeros(coordinate_size(geom))
    function Fexc(ρ::NTuple{Nspecies,<:AbstractArray}, T::Real=1.0)
        integrand .= ρ[1] .* fexc.(functional, ρ[1], T)
        integrate(integrand)
    end
end

function c1!_func(functional::LDAFunctional, ::Val{Nspecies}) where {Nspecies}
    Nspecies == 1 || throw(ArgumentError("Nspecies > 1 left as an exercise for the reader"))
    function c1!(result::NTuple{Nspecies,<:AbstractArray}, ρ::NTuple{Nspecies,<:AbstractArray}, T::Real=1.0)
        result[1] .= .-(fexc.(functional, ρ[1], T) .+ ρ[1] .* dfexc.(functional, ρ[1], T)) ./ T
    end
end