Examples

Concrete, runnable examples that exercise the major code paths.

Open-BC Coulomb, 5 charges in 3D

using MultilevelSummation
using StaticArrays
using Random
Random.seed!(42)

positions = [SVector{3,Float64}((rand(3) .* 2)...) for _ in 1:5]
charges   = randn(5) ./ sqrt(5)

cell     = zero(SMatrix{3,3,Float64})            # unused for open BC
periodic = (false, false, false)

# Baseline: direct O(N²) Coulomb
using MultilevelSummation.Reference: naive_energy
K = Coulomb()
U_naive = naive_energy(positions, charges, cell, periodic, K)

# MSM with a fairly small a — splitting error is the dominant approximation here.
calc = MSMCalculator(HardyC2Cubic(2.0, 3), CubicC1(), 0.2)
U_msm = msm_energy(positions, charges, cell, periodic, calc)

(U_naive, U_msm, abs(U_msm - U_naive) / abs(U_naive))

(-0.08508664602342171, -0.07634694482575935, 0.10271530970037993)

Convergence in a

Holding h fixed and growing a should reduce the splitting error monotonically.

using MultilevelSummation, StaticArrays, Random
using MultilevelSummation.Reference: naive_energy
Random.seed!(0xC0FFEE)

positions = [SVector{3,Float64}((rand(3) .* 2)...) for _ in 1:8]
charges   = randn(8); charges .-= sum(charges)/length(charges)
cell, periodic = zero(SMatrix{3,3,Float64}), (false, false, false)

U_naive = naive_energy(positions, charges, cell, periodic, Coulomb())

for a in (0.8, 1.6, 3.2, 6.4)
    calc = MSMCalculator(HardyC2Cubic(a, 3), CubicC1(), 0.2)
    U = msm_energy(positions, charges, cell, periodic, calc)
    @show a, abs(U - U_naive) / abs(U_naive)
end

(a, abs(U - U_naive) / abs(U_naive)) = (0.8, 0.07483594153438516)
(a, abs(U - U_naive) / abs(U_naive)) = (1.6, 0.013746410333744958)
(a, abs(U - U_naive) / abs(U_naive)) = (3.2, 0.002401370669336198)
(a, abs(U - U_naive) / abs(U_naive)) = (6.4, 0.0005538276368636638)

Periodic 3D Coulomb

using MultilevelSummation, StaticArrays, Random
Random.seed!(7)

h, L = 0.5, 4                        # 8 fine points per axis ⇒ top is 1×1×1
cell_L = h * 8
cell = SMatrix{3,3,Float64}(cell_L * one(SMatrix{3,3,Float64}))
periodic = (true, true, true)

N = 8
positions = [SVector{3,Float64}((rand(3) .* cell_L)...) for _ in 1:N]
charges   = randn(N); charges .-= sum(charges)/N

calc = MSMCalculator(HardyC2Cubic(2.0, L), CubicC1(), h)
U, F = msm_energy_forces(positions, charges, cell, periodic, calc)
(U, sum(F))                          # ΣF ≈ 0 expected

(-1.772514777870421, [0.029026986702012814, -0.03124331469721546, 0.01615475477658279])

Through the AtomsCalculators interface

using MultilevelSummation, AtomsBase, AtomsCalculators, Unitful, StaticArrays

L = 4.0u"Å"
cell_vec = (SVector(L, 0u"Å", 0u"Å"),
            SVector(0u"Å", L, 0u"Å"),
            SVector(0u"Å", 0u"Å", L))
sys = periodic_system([
    Atom(:Na, SVector(1.0u"Å", 2.0u"Å", 0.5u"Å"); charge = +1.0),
    Atom(:Cl, SVector(2.5u"Å", 1.5u"Å", 3.0u"Å"); charge = -1.0),
], cell_vec)

calc = MSMCalculator(HardyC2Cubic(2.0, 4), CubicC1(), 0.5)
ef = AtomsCalculators.energy_forces(sys, calc)
(ef.energy, ef.forces[1])

(-0.5859876127899326, [0.0707274850556901, -0.03910912379936127, -0.0698562289975313])

Mixed BC: periodic in xy, open in z

using MultilevelSummation, StaticArrays, Random
Random.seed!(11)

# 4×4 periodic in xy, open in z
L_xy   = 4.0
cell   = SMatrix{3,3,Float64}([L_xy 0 0; 0 L_xy 0; 0 0 1.0])
periodic = (true, true, false)

positions = [SVector(2.0, 1.5, 0.3), SVector(0.5, 3.2, -0.4)]
charges   = [1.0, -1.0]

calc = MSMCalculator(HardyC2Cubic(2.0, 3), CubicC1(), 0.5)
msm_energy(positions, charges, cell, periodic, calc)

-0.5355393303057003

Naive Ewald reference

A naive 3D Ewald reference lives in the public MultilevelSummation.Reference submodule. Used in tests for periodic-Coulomb accuracy checks and as the reference inside MultilevelSummation.Tune.sweep:

using MultilevelSummation.Reference: ewald_energy, ewald_energy_forces

U      = ewald_energy(positions, charges, cell; α=0.7, R_cut=10.0, k_cut=12.0)
U, F   = ewald_energy_forces(positions, charges, cell; α=0.7, R_cut=10.0, k_cut=12.0)

For auto-tuned parameters use the Tune wrapper:

using MultilevelSummation.Tune
U = Tune.ewald_reference(positions, charges, cell; tol = 1e-9)

Both perform the standard erfc / reciprocal-Gaussian split and self- validate via α-invariance in test/test_ewald.jl.

Programmatic hyperparameter sweep

using MultilevelSummation
using MultilevelSummation.Tune
using StaticArrays, Random
Random.seed!(0xABCD)

L = 8.0; cell = SMatrix{3,3,Float64}(L * one(SMatrix{3,3,Float64}))
periodic = (true, true, true)
positions = [SVector{3,Float64}((rand(3) .* L)...) for _ in 1:32]
charges = randn(32); charges .-= sum(charges) / length(charges)

U_ref = Tune.ewald_reference(positions, charges, cell)
results = Tune.sweep(positions, charges, cell, periodic;
                    reference  = U_ref,
                    h_values   = (0.5, 1.0),
                    a_values   = (1.0, 2.0),
                    L_strategy = :all)
length(results), Tune.recommend(results; max_rel_err = 0.1)

(14, SweepResult(h=1, a=2, L=3, n_grid=8, energy=-5.105100e+00, rel_err=4.811e-02, t_msm=0.0012s))

The shipped tuning/tune_NaCl.jl is a thin caller of this same API on realistic NaCl supercells — see tuning/README.md for how to run it.