First steps

In order to solve an Ornstein-Zernike equation, there are a number of things that need to be specified.

1. The interaction potential,
2. The system parameters: density, temperature and dimensionality
3. The closure relation
4. The solver

Let's start with a very simple example: a three-dimensional 1-component system, where the particles interact according to an inverse power law potential $u(r)=\epsilon (\sigma/r)^n$. In this case, we can make use of the built-in potential

using OrnsteinZernike
ϵ = 1.0
σ = 1.0
n = 8
potential = PowerLaw(ϵ, σ, n)

PowerLaw{Float64, Float64, Int64}(1.0, 1.0, 8)

Now that we have the potential, we define the system

dims = 3 # we consider a 3D system
ρ = 0.6 # number density
kBT = 1.0 # thermal energy
system = SimpleLiquid(dims, ρ, kBT, potential)

3 dimensional SimpleLiquid:
 ρ = 0.6
 kBT = 1.0
 potential = PowerLaw{Float64, Float64, Int64}(1.0, 1.0, 8)

The SimpleLiquid object is meant to be used when dealing with systems that have spherically symmetric interaction potentials and no external fields.

The third step is to define a closure relation. For now, let's stick to the simple Hypernetted Chain closure

closure = HypernettedChain()

HypernettedChain()

We can now solve the system.

sol = solve(system, closure)

OZSolution{Vector{Float64}, Vector{Float64}}:
 r = [0.005 0.015 … 9.985 9.995000000000001]'
 k = [0.15707963267948966 0.47123889803846897 … 313.6880264609408 314.00218572629984]'
 cr = [-7.699093104727448 -7.695706108519444 … -1.012081396467579e-8 -1.0040090092644505e-8]'
 gr = [-1.0243272896559574e-10 -1.0231193670051653e-10 … 0.9999999938362887 0.9999999938605444]'
 ck = [-16.014517909323175 -15.680903323809734 … 1.552960283565252e-16 -5.634497317451505e-17]'
 Sk = [0.09426216097138174 0.09607493542793866 … 1.0 1.0]

Note that we have not specified the method by which we do so. In this case, a simple default method is chosen that works well in most cases.

The sol object contains fields for the radial distribution function gr, direct correlation function (in real and Fourier space) cr and ck, the static structure factor Sk, and arrays for r and k. For example, we can plot the radial distribution function as follows:

using Plots
plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")

Full code:

using OrnsteinZernike
ϵ = 1.0
σ = 1.0
n = 8
potential = PowerLaw(ϵ, σ, n)
dims = 3 # we consider a 3D system
ρ = 0.6 # number density
kBT = 1.0 # thermal energy
system = SimpleLiquid(dims, ρ, kBT, potential)
closure = HypernettedChain()
sol = solve(system, closure)
using Plots
plot(sol.r, sol.gr, xlims=(0,5), xlabel="r", ylabel="g(r)")