Theory

In this section, we describe our conventions and notation.

Single component systems

The Ornstein zernike equation $h(r) = c(r) + \rho \int d\textbf{r}' c(\textbf{r}')h(|\textbf{r} - \textbf{r}'|) $ links the pair correlation function $g(r) = h(r)+1$ to the number density $\rho$ and the direct correlation function $c(r)$. In order to solve it, a second relation between $h(r)$ and $c(r)$ must be specified. This is called the closure relation.

In practise, the closure relation typically takes the form $c(r) = f(\gamma(r), r, u(r))$, where $u(r)$ is the interaciton potential and $\gamma(r) = h(r) - c(r)$ is the indirect correlation function. Sometimes, instead the closure is defined for the bridge function $b(r)$, and the closure relation is then given by $h(r) - 1 = \exp\left(-\beta u(r) + h(r) - c(r) + b(r) \right)$, where $\beta = 1/k_BT$.

Mixtures

Everything above generalizes to the mixture case:

The Ornstein zernike equation

\[h_{ij}(r) = c_{ij}(r) + \sum_l \rho_l \int d\textbf{r}' c_{il}(\textbf{r}')h_{lj}(|\textbf{r} - \textbf{r}'|)\]

links the pair correlation function $g_{ij}(r) = h_{ij}(r)+1$ to the species specific number density $\rho_{i}$ and the direct correlation function $c_{ij}(r)$. In order to solve it, a second relation between $h_{ij}(r)$ and $c_{ij}(r)$ must be specified. This is called the closure relation.

In practise, the closure relation typically takes the form $c_{ij}(r) = f(\gamma_{ij}(r), r, u_{ij}(r))$, where $u_{ij}(r)$ is the interaciton potential and $\gamma_{ij}(r) = h_{ij}(r) - c_{ij}(r)$ is the indirect correlation function. Sometimes, instead the closure is defined for the bridge function $b(r)$, and the closure relation is then given by $h_{ij}(r) - 1 = \exp\left(-\beta u_{ij}(r) + h_{ij}(r) - c_{ij}(r) + b_{ij}(r) \right)$, where $\beta = 1/k_BT$.

Fourier Transforms

The ability to numerically solve the Ornstein-Zernike equation relies heavily on doing repeated Fourier Transforms. In arbitrary dimensions, these Fourier transforms can be written as Hankel transforms in the case that the argument is a radially symmetric function. In particular, in $d$ dimensions, the radial Fourier transform and its inverse are given by

\[\hat{F}(k) = (2\pi)^{d/2} k ^{1-d/2}\int_0^\infty dr r^{d/2}J_{d/2-1}(kr)F(r)\]

\[F(r) = (2\pi)^{-d/2} r ^{1-d/2}\int_0^\infty dk k^{d/2}J_{d/2-1}(kr)\hat{F}(k),\]

in which $J_{d/2-1}(x)$ is the bessel function of order $d/2-1$. In the special cases of 1 and 3 dimensions, the transform simplifies into:

\[\hat{F}(k) = 2\int_0^\infty dr \cos(kr)F(r)\]

\[F(r) = \frac{1}{\pi}\int_0^\infty dk \cos(kr)\hat{F}(k),\]

in 1$d$, and

\[\hat{F}(k) = \frac{4\pi}{k}\int_0^\infty dr r \sin(kr)F(r)\]

\[F(r) = \frac{1}{2\pi^2r}\int_0^\infty dk k\sin(kr)\hat{F}(k),\]

in 3$d$. This package uses discrete versions of all of the above.

Thermodynamic properties

Using the structure as determined by this package, several thermodynamic properies can be computed. In particular, this package contains methods to compute the (virial) pressure $p$, the isothermal compressibility $\chi$, and the excess internal energy per particle $E_x$.

For mixtures, they are computed respectively from the following definitions

\[p = k_BT \rho_0\sum_i x_i - 1/6 \rho_0^2 \sum_{ij} x_i x_j \int d\textbf{r} r g_{ij}(r) u'_{ij}(r)\]

\[1/(\rho_0 k_BT \chi) = 1 - ρ_0 \sum_{ij} x_i x_j \hat{c}_{ij}(k\to0)\]

,

\[E_x = 1/2 \rho_0 \sum_{ij} x_i x_j \int d\textbf{r} g_{ij}(r) u_{ij}(r)\]

in which $\rho_0=N/V$. The functions to use are compute_virial_pressure, compute_compressibility, and, compute_excess_energy.

Thermodynamic Routes

The Ornstein–Zernike framework gives several thermodynamic quantities:

• Virial route (pressure): from (g(r)) and (u'(r)).
• Compressibility route:
  • compute_compressibility returns isothermal compressibility (\chi).
  • Pressure via compressibility route requires integration over density using [ \frac{\partial (\beta p)}{\partial \rho} \;=\; \frac{1}{S(0)}, \qquad S(0) = \rho\,kBT\,\kappaT, ] so [ \beta p(\rho) \;=\; \int_0^\rho \frac{d\rho'}{S(0;\rho')}. ]
• Energy route: excess internal energy from (g(r)).

Because closures are approximate, routes generally disagree; the gap is a measure of closure error.

Example: Virial Pressure vs. Compressibility-Route Pressure

Below we:

1. compute virial pressure directly,
2. compute (\chi) at a sequence of densities,
3. integrate (1/S(0)) to obtain the compressibility-route pressure.

using OrnsteinZernike

# --- System definition at target temperature and potential ---
kBT     = 1.0
σ, ϵ    = 1.0, 1.0
closure = Verlet()   # any closure works; V shown here
potential = PowerLaw(ϵ, σ, 8)

# Choose a target density and also a ramp to integrate from 0 → ρ_target
ρ_target   = 0.8
delta_ρ = 0.01
ρ_grid     = collect((delta_ρ/2):delta_ρ:(ρ_target-delta_ρ/2))  # avoid 0 to keep numerics stable
system_at = ρ -> SimpleFluid(3, ρ, kBT, potential)

# (1) Virial-route pressure at ρ_target
sol_target = solve(system_at(ρ_target), closure, NgIteration(M=5000, dr=0.01))
p_virial   = compute_virial_pressure(sol_target, system_at(ρ_target))   # this is "virial pressure"

# (2) Compressibility κ_T along a density grid
chi_values = similar(ρ_grid)
sols = solve(system_at(ρ_grid[end]), closure, DensityRamp(NgIteration(M=5000, dr=0.01), ρ_grid))

for (i, ρ) in enumerate(ρ_grid)
    chi_values[i] = compute_compressibility(sols[i], system_at(ρ))  # Isothermal compressibility κ_T
end

# (3) Integrate to get compressibility-route pressure:
#     d(βp)/dρ = 1 / S(0) with S(0) = ρ * kBT * κ_T
S0 = ρ_grid .* kBT .* chi_values

# Midpoint integration of 1/S0 over ρ to get β p(ρ)
βp  = delta_ρ * sum(1 ./ S0)
p_comp = kBT * βp    # convert βp → p

println("Virial-route pressure at          ρ = $(ρ_target): p_virial = $(p_virial)")
println("Compressibility-route pressure at ρ = $(ρ_target): p_comp   = $(p_comp)")

After iteration 0, the error is 2.1507597035.
After iteration 10, the error is 0.00251783142.
After iteration 20, the error is 4.0e-11.
Converged after 21 iterations, the error is 4.1e-11.

Solving the system at ρ = 0.005.

After iteration 0, the error is 1.62401e-6.
Converged after 2 iterations, the error is 3.0e-13.

Solving the system at ρ = 0.015.

After iteration 0, the error is 7.674102e-5.
Converged after 3 iterations, the error is 4.2e-13.

Solving the system at ρ = 0.025.

After iteration 0, the error is 0.00031568423.
Converged after 3 iterations, the error is 8.9e-13.

Solving the system at ρ = 0.035.

After iteration 0, the error is 0.00077128395.
Converged after 4 iterations, the error is 1.0e-12.

Solving the system at ρ = 0.045.

After iteration 0, the error is 0.00146249159.
Converged after 4 iterations, the error is 1.3e-11.

Solving the system at ρ = 0.055.

After iteration 0, the error is 0.00238687308.
Converged after 4 iterations, the error is 8.9e-11.

Solving the system at ρ = 0.065.

After iteration 0, the error is 0.00352841026.
Converged after 5 iterations, the error is 6.3e-14.

Solving the system at ρ = 0.075.

After iteration 0, the error is 0.00486302617.
Converged after 5 iterations, the error is 1.8e-13.

Solving the system at ρ = 0.085.

After iteration 0, the error is 0.00636244919.
Converged after 5 iterations, the error is 4.2e-13.

Solving the system at ρ = 0.095.

After iteration 0, the error is 0.00799686893.
Converged after 5 iterations, the error is 8.9e-13.

Solving the system at ρ = 0.105.

After iteration 0, the error is 0.00973671746.
Converged after 5 iterations, the error is 1.7e-12.

Solving the system at ρ = 0.115.

After iteration 0, the error is 0.01155382136.
Converged after 5 iterations, the error is 3.0e-12.

Solving the system at ρ = 0.125.

After iteration 0, the error is 0.01342210474.
Converged after 5 iterations, the error is 3.7e-12.

Solving the system at ρ = 0.135.

After iteration 0, the error is 0.01531797528.
Converged after 5 iterations, the error is 4.5e-12.

Solving the system at ρ = 0.145.

After iteration 0, the error is 0.01722049002.
Converged after 5 iterations, the error is 7.1e-12.

Solving the system at ρ = 0.155.

After iteration 0, the error is 0.01911137079.
Converged after 5 iterations, the error is 7.7e-12.

Solving the system at ρ = 0.165.

After iteration 0, the error is 0.02097492048.
Converged after 5 iterations, the error is 7.4e-11.

Solving the system at ρ = 0.175.

After iteration 0, the error is 0.02279787633.
Converged after 6 iterations, the error is 3.7e-12.

Solving the system at ρ = 0.185.

After iteration 0, the error is 0.02456922673.
Converged after 6 iterations, the error is 2.6e-11.

Solving the system at ρ = 0.195.

After iteration 0, the error is 0.02628000935.
Converged after 6 iterations, the error is 5.1e-11.

Solving the system at ρ = 0.205.

After iteration 0, the error is 0.02792310381.
Converged after 6 iterations, the error is 7.9e-11.

Solving the system at ρ = 0.215.

After iteration 0, the error is 0.029493027.
Converged after 6 iterations, the error is 8.4e-11.

Solving the system at ρ = 0.225.

After iteration 0, the error is 0.03098573703.
Converged after 7 iterations, the error is 7.7e-12.

Solving the system at ρ = 0.235.

After iteration 0, the error is 0.03239844903.
Converged after 7 iterations, the error is 1.4e-11.

Solving the system at ρ = 0.245.

After iteration 0, the error is 0.03372946491.
Converged after 7 iterations, the error is 3.3e-11.

Solving the system at ρ = 0.255.

After iteration 0, the error is 0.03497801778.
Converged after 7 iterations, the error is 8.2e-12.

Solving the system at ρ = 0.265.

After iteration 0, the error is 0.0361441314.
Converged after 8 iterations, the error is 4.8e-12.

Solving the system at ρ = 0.275.

After iteration 0, the error is 0.03722849404.
Converged after 8 iterations, the error is 6.2e-12.

Solving the system at ρ = 0.285.

After iteration 0, the error is 0.03823234631.
Converged after 8 iterations, the error is 8.5e-12.

Solving the system at ρ = 0.295.

After iteration 0, the error is 0.03915738191.
Converged after 8 iterations, the error is 1.3e-11.

Solving the system at ρ = 0.305.

After iteration 0, the error is 0.04000566041.
Converged after 8 iterations, the error is 1.9e-11.

Solving the system at ρ = 0.315.

After iteration 0, the error is 0.04077953108.
Converged after 8 iterations, the error is 2.9e-11.

Solving the system at ρ = 0.325.

After iteration 0, the error is 0.04148156658.
Converged after 8 iterations, the error is 4.5e-11.

Solving the system at ρ = 0.335.

After iteration 0, the error is 0.04211450575.
Converged after 8 iterations, the error is 6.7e-11.

Solving the system at ρ = 0.345.

After iteration 0, the error is 0.0426812044.
Converged after 8 iterations, the error is 1.0e-10.

Solving the system at ρ = 0.355.

After iteration 0, the error is 0.04318459328.
Converged after 9 iterations, the error is 6.4e-12.

Solving the system at ρ = 0.365.

After iteration 0, the error is 0.04362764241.
Converged after 9 iterations, the error is 1.1e-11.

Solving the system at ρ = 0.375.

After iteration 0, the error is 0.04401333103.
Converged after 9 iterations, the error is 1.8e-11.

Solving the system at ρ = 0.385.

After iteration 0, the error is 0.04434462252.
Converged after 9 iterations, the error is 3.0e-11.

Solving the system at ρ = 0.395.

After iteration 0, the error is 0.04462444362.
Converged after 9 iterations, the error is 5.0e-11.

Solving the system at ρ = 0.405.

After iteration 0, the error is 0.04485566746.
Converged after 9 iterations, the error is 8.1e-11.

Solving the system at ρ = 0.415.

After iteration 0, the error is 0.04504109986.
Converged after 10 iterations, the error is 2.3e-13.

Solving the system at ρ = 0.425.

After iteration 0, the error is 0.04518346853.
Converged after 10 iterations, the error is 4.9e-13.

Solving the system at ρ = 0.435.

After iteration 0, the error is 0.04528541474.
Converged after 10 iterations, the error is 8.4e-13.

Solving the system at ρ = 0.445.

After iteration 0, the error is 0.04534948707.
Converged after 10 iterations, the error is 1.2e-12.

Solving the system at ρ = 0.455.

After iteration 0, the error is 0.04537813711.
Converged after 10 iterations, the error is 1.7e-12.

Solving the system at ρ = 0.465.

After iteration 0, the error is 0.04537371661.
Converged after 10 iterations, the error is 2.5e-12.

Solving the system at ρ = 0.475.

After iteration 0, the error is 0.04533847608.
Converged after 10 iterations, the error is 3.8e-12.

Solving the system at ρ = 0.485.

After iteration 0, the error is 0.04527456451.
Converged after 10 iterations, the error is 6.1e-12.

Solving the system at ρ = 0.495.

After iteration 0, the error is 0.04518403005.
Converged after 10 iterations, the error is 1.0e-11.

Solving the system at ρ = 0.505.

After iteration 0, the error is 0.04506882161.
Converged after 10 iterations, the error is 1.7e-11.

Solving the system at ρ = 0.515.

After iteration 0, the error is 0.04493079108.
Converged after 10 iterations, the error is 3.1e-11.

Solving the system at ρ = 0.525.

After iteration 0, the error is 0.04477169631.
Converged after 10 iterations, the error is 5.6e-11.

Solving the system at ρ = 0.535.

After iteration 0, the error is 0.04459320452.
After iteration 10, the error is 2.0e-11.
Converged after 11 iterations, the error is 1.7e-11.

Solving the system at ρ = 0.545.

After iteration 0, the error is 0.0443968963.
After iteration 10, the error is 2.0e-11.
Converged after 11 iterations, the error is 2.2e-11.

Solving the system at ρ = 0.555.

After iteration 0, the error is 0.04418426999.
After iteration 10, the error is 2.0e-11.
Converged after 11 iterations, the error is 1.9e-11.

Solving the system at ρ = 0.565.

After iteration 0, the error is 0.04395674654.
After iteration 10, the error is 1.0e-11.
Converged after 11 iterations, the error is 6.3e-12.

Solving the system at ρ = 0.575.

After iteration 0, the error is 0.04371567475.
After iteration 10, the error is 0.0.
Converged after 11 iterations, the error is 4.8e-12.

Solving the system at ρ = 0.585.

After iteration 0, the error is 0.04346233693.
After iteration 10, the error is 1.0e-11.
Converged after 11 iterations, the error is 7.8e-12.

Solving the system at ρ = 0.595.

After iteration 0, the error is 0.04319795502.
After iteration 10, the error is 0.0.
Converged after 11 iterations, the error is 4.6e-12.

Solving the system at ρ = 0.605.

After iteration 0, the error is 0.04292369721.
Converged after 10 iterations, the error is 6.0e-11.

Solving the system at ρ = 0.615.

After iteration 0, the error is 0.04264068498.
After iteration 10, the error is 1.0e-11.
Converged after 11 iterations, the error is 6.6e-12.

Solving the system at ρ = 0.625.

After iteration 0, the error is 0.04235000088.
After iteration 10, the error is 2.0e-11.
Converged after 11 iterations, the error is 1.6e-11.

Solving the system at ρ = 0.635.

After iteration 0, the error is 0.04205269689.
After iteration 10, the error is 2.0e-11.
Converged after 11 iterations, the error is 1.6e-11.

Solving the system at ρ = 0.645.

After iteration 0, the error is 0.04174980354.
After iteration 10, the error is 1.0e-11.
Converged after 11 iterations, the error is 1.5e-11.

Solving the system at ρ = 0.655.

After iteration 0, the error is 0.04144234002.
After iteration 10, the error is 1.0e-11.
Converged after 11 iterations, the error is 1.5e-11.

Solving the system at ρ = 0.665.

After iteration 0, the error is 0.04113132524.
After iteration 10, the error is 1.0e-11.
Converged after 11 iterations, the error is 1.4e-11.

Solving the system at ρ = 0.675.

After iteration 0, the error is 0.04081779007.
After iteration 10, the error is 1.0e-11.
Converged after 11 iterations, the error is 1.1e-11.

Solving the system at ρ = 0.685.

After iteration 0, the error is 0.04050279099.
After iteration 10, the error is 0.0.
Converged after 11 iterations, the error is 4.9e-12.

Solving the system at ρ = 0.695.

After iteration 0, the error is 0.04018742519.
After iteration 10, the error is 1.0e-11.
Converged after 11 iterations, the error is 9.1e-12.

Solving the system at ρ = 0.705.

After iteration 0, the error is 0.03987284753.
After iteration 10, the error is 2.0e-11.
Converged after 11 iterations, the error is 2.4e-11.

Solving the system at ρ = 0.715.

After iteration 0, the error is 0.03956028934.
After iteration 10, the error is 4.0e-11.
Converged after 11 iterations, the error is 4.2e-11.

Solving the system at ρ = 0.725.

After iteration 0, the error is 0.0392510795.
After iteration 10, the error is 6.0e-11.
Converged after 11 iterations, the error is 6.1e-11.

Solving the system at ρ = 0.735.

After iteration 0, the error is 0.03894666784.
After iteration 10, the error is 8.0e-11.
Converged after 11 iterations, the error is 7.7e-11.

Solving the system at ρ = 0.745.

After iteration 0, the error is 0.03864865123.
After iteration 10, the error is 9.0e-11.
Converged after 11 iterations, the error is 8.7e-11.

Solving the system at ρ = 0.755.

After iteration 0, the error is 0.0383588024.
After iteration 10, the error is 9.0e-11.
Converged after 11 iterations, the error is 9.2e-11.

Solving the system at ρ = 0.765.

After iteration 0, the error is 0.03807910178.
After iteration 10, the error is 1.0e-10.
Converged after 12 iterations, the error is 5.7e-12.

Solving the system at ρ = 0.775.

After iteration 0, the error is 0.03781177219.
After iteration 10, the error is 9.0e-11.
Converged after 11 iterations, the error is 8.7e-11.

Solving the system at ρ = 0.785.

After iteration 0, the error is 0.03755931664.
After iteration 10, the error is 5.9e-10.
Converged after 12 iterations, the error is 2.1e-11.

Solving the system at ρ = 0.795.

After iteration 0, the error is 0.0373245587.
After iteration 10, the error is 1.11e-9.
Converged after 12 iterations, the error is 6.7e-11.

Solving the system at ρ = 0.795.

After iteration 0, the error is 0.0.
Converged after 1 iterations, the error is 2.3e-12.
Virial-route pressure at          ρ = 0.8: p_virial = 6.170761036688983
Compressibility-route pressure at ρ = 0.8: p_comp   = 7.291206648140269

Notes

• compute_compressibility(sol, system) returns (\chi) (units of inverse pressure).
• We used a density ramp inside the loop to help convergence and reuse previous solutions.

Why Routes Differ

• Virial route is sensitive to short-range structure (contact region).
• Compressibility route probes long-wavelength fluctuations via (S(0)).
• Discrepancies reflect closure approximation error.

Making Routes Agree

Closures like Rogers–Young include a tunable parameter (e.g., (\alpha)) that can be chosen so that ( p{\mathrm{vir}} \approx p{\mathrm{comp}} ). See the Thermodynamic Consistency tutorial for a bisection example.