Choosing Grid Parameters (M and dr)

The Ornstein–Zernike solvers use radial Fourier/Hankel transforms to go between real and Fourier space. The grid is controlled by two parameters:

• M : number of grid points
• dr: approximate radial step size (may not be truly uniform, since Hankel transforms use non-uniform grids in arbitrary dimensions)

The maximum radius is

[ r_{\max} = M \cdot dr ]

which sets the size of the box in real space.

What Must Be Resolved?

• In real space, it is the direct correlation function (c(r)) that is Fourier transformed.

  • Therefore, dr and rmax must be chosen such that oscillations and decay of (c(r)) are well resolved.
  • Poor resolution of (c(r)) will produce inaccurate structure factors (S(k)).
  • In a three-dimensional system, if the potential has discontinuities, they must lie exactly between two gridpoints to ensure second order convergence.

• In Fourier space, the discretization also needs to be accurate.

  • If the (k)-grid is too coarse, long-wavelength properties (e.g. (S(k\to 0))) will be unreliable. Ensure that the direct correlation function (\gamma(k)) is well discretized.
  • If (r_{\max}) is too small, aliasing and finite-box effects may distort the Fourier transform.

Practical Guidelines

• Choose dr small enough to resolve short-range features of the potential and resulting (c(r)).
  • For hard spheres or Lennard–Jones: dr ≲ 0.05σ.
• Choose rmax large enough so that (c(r)) has decayed essentially to zero before the cutoff.
  • A good starting point: rmax ≈ 10σ.
• Increase M together with decreasing dr if necessary — this improves both real-space and Fourier-space resolution.
• Check convergence: repeat a calculation with finer grid and verify (g(r)), (S(k)), and thermodynamic routes are stable.

Example: Grid Convergence with dr

using OrnsteinZernike, Plots

system = SimpleFluid(3, 0.8, 1.0, LennardJones(1.0, 1.0))
closure = HypernettedChain()

Rmax = 10.0
p1 = plot(xlims=(0,2), xlabel="r", ylabel="c(r)")
p2 = plot(xlims=(0,20), xlabel="k", ylabel="c(k)")
for M in [1000, 100, 50]
    dr = Rmax/M
    sol  = @time solve(system, closure, NgIteration(M=M, dr=dr))
    kmax = round(maximum(sol.k), digits=2)
    plot!(p1, sol.r, sol.cr, label="M=$M, kmax=$(kmax)", ls=:dash, markers=:o, legend=false)
    plot!(p2, sol.k, sol.ck, label="M=$M, kmax=$(kmax)", ls=:dash, markers=:o, legend=:bottomright)
end
display(plot(p1,p2))

After iteration 0, the error is 8.13170801569.
After iteration 10, the error is 6.2030864521.
After iteration 20, the error is 0.38583419355.
After iteration 30, the error is 3.35326e-5.
After iteration 40, the error is 2.402e-8.
Converged after 46 iterations, the error is 8.5e-11.
  0.913790 seconds (737.22 k allocations: 37.367 MiB, 99.76% compilation time)
After iteration 0, the error is 2.70974354173.
After iteration 10, the error is 1.77538968751.
After iteration 20, the error is 0.10734410666.
After iteration 30, the error is 2.972445e-5.
After iteration 40, the error is 3.23e-9.
Converged after 45 iterations, the error is 6.5e-11.
  0.000425 seconds (562 allocations: 885.734 KiB)
After iteration 0, the error is 53.96901689238.
After iteration 10, the error is 0.93409779905.
After iteration 20, the error is 0.08035914126.
After iteration 30, the error is 0.00235930013.
After iteration 40, the error is 3.892e-8.
Converged after 50 iterations, the error is 6.0e-11.
  0.000305 seconds (592 allocations: 872.500 KiB)

Here the coarse grid under-resolves oscillations in (c(r)) and doesn't capture it's decay in Fourier space, which in turn leads to distorted structure factors (S(k)).