feat: implement Verlet neighbour list (#65)

**Summary**
This PR implements a flavor of Verlet list. It addresses #63. 

The VerletList object contains for each particle in the system, an array
of all its neighbours (defined as all particles within the interaction
cutoff + some arbitrary parameter dr).
These arrays of neighbours are updated only when a given particle moves
a distance greater than dr/2, ensuring that the list of neighbours
always contains all particles within the interaction cutoff radius.

**Test plan**
The CI was updated to check over a trajectory that the energy is
identical to that obtained using different neighbour list methods.

**Performances**
The VerletList implementation was compared to the LinkedList one, for a
system of ortho-terphenyl molecules, at density 1.2 and temperature 2.
In the attached plot, the simulation time is reported for different
values of the parameter dr. The dashed line at the top is the timing
with the LinkedList implementation. For reference, the interaction
cutoff is 1.23.
The timing is done once for each Verlet list, and 4 times for the
LinkedList, with values ranging from 75s to 105s.
There's quite a lot of variability between repeats, but we still have a
~2x speedup for this system, for 0.3 < dr < 0.6.
<img width="640" height="480" alt="verlet"
src="https://github.com/user-attachments/assets/a2175cb8-1746-4092-b414-528dab85a5fe"
/>

Author: V-Francois

src/IO/IO.jl

@@ -2,7 +2,7 @@ module IO
 
 using ..ParticlesMC: Particles, Atoms, Molecules, System
 using ..ParticlesMC: fold_back, volume_sphere
-using ..ParticlesMC: EmptyList, LinkedList, CellList
+using ..ParticlesMC: EmptyList, LinkedList, CellList, VerletList
 using ..ParticlesMC: Model, GeneralKG, JBB, BHHP, SoftSpheres, KobAndersen, Trimer, LennardJones
 using Arianna
 using Distributions, LinearAlgebra, StaticArrays, Printf
@@ -286,16 +286,17 @@ function load_chains(init_path; args=Dict(), verbose=false)
     if haskey(args, "list_type") && !isnothing(args["list_type"])
         list_type = eval(Meta.parse(args["list_type"]))
     end
+    list_parameters = get(args, "list_parameters", nothing)
     verbose && println("Using $list_type as cell list type")
     # Create independent chains
     bool_molecule = :molecule in keys(config_dict[1])
     if bool_molecule
         initial_molecule_array = broadcast_dict(config_dict, :molecule)
         initial_bond_array = broadcast_dict(config_dict, :bond)
         #initial_btype_array = broadcast_dict(config_dict, :btype)
-        chains = [System(initial_position_array[k], initial_species_array[k], initial_molecule_array[k], initial_density_array[k], initial_temperature_array[k], model_matrix, initial_bond_array[k], list_type=list_type) for k in eachindex(initial_position_array)]
+        chains = [System(initial_position_array[k], initial_species_array[k], initial_molecule_array[k], initial_density_array[k], initial_temperature_array[k], model_matrix, initial_bond_array[k], list_type=list_type, list_parameters=list_parameters) for k in eachindex(initial_position_array)]
     else
-        chains = [System(initial_position_array[k], initial_species_array[k], initial_density_array[k], initial_temperature_array[k], model_matrix, list_type=list_type) for k in eachindex(initial_position_array)]
+        chains = [System(initial_position_array[k], initial_species_array[k], initial_density_array[k], initial_temperature_array[k], model_matrix, list_type=list_type, list_parameters=list_parameters) for k in eachindex(initial_position_array)]
     end
     verbose && println("$(length(chains)) chains created")
     return chains

src/ParticlesMC.jl

@@ -115,7 +115,7 @@ end
 export callback_energy
 #export nearest_image_distance
 export Model, GeneralKG, JBB, BHHP, SoftSpheres, KobAndersen, Trimer
-export NeighbourList, LinkedList, CellList, EmptyList
+export NeighbourList, LinkedList, CellList, EmptyList, VerletList
 export Atoms, Molecules
 export Displacement, DiscreteSwap, MoleculeFlip
 export fold_back, System
@@ -153,6 +153,7 @@ ParticlesMC implemented in Comonicon.
         error("Configuration file '$config' does not exist in the current path.")
     end
     list_type = get(system, "list_type", "LinkedList")  # optional field
+    list_parameters = get(system, "list_parameters", nothing)  # optional field
     bonds = get(system, "bonds", nothing)
 
     # Extract simulation parameters
@@ -173,6 +174,7 @@ ParticlesMC implemented in Comonicon.
             "density" => [density],
             "model" => [model],
             "list_type" => list_type,
+            "list_parameters" => list_parameters,
             "bonds" => bonds,
         ))
     else
@@ -181,6 +183,7 @@ ParticlesMC implemented in Comonicon.
             "density" => [density],
             "model" => [model],
             "list_type" => list_type,
+            "list_parameters" => list_parameters,
         ))
     end
     algorithm_list = []

src/atoms.jl

@@ -37,14 +37,14 @@ Create an `Atoms` system, initialize its neighbour list, and compute initial ene
 constructs an `Atoms` instance, builds the neighbour list of type `list_type`,
 and computes the initial total energy (stored in `system.energy[1]`).
 """
-function System(position, species, density::T, temperature::T, model_matrix; list_type=EmptyList) where{T<:AbstractFloat}
+function System(position, species, density::T, temperature::T, model_matrix; list_type=EmptyList, list_parameters=nothing) where{T<:AbstractFloat}
     @assert length(position) == length(species)
     N = length(position)
     d = length(Array(position)[1])
     energy = Vector{T}(undef, 1)
     box = @SVector fill(T((N / density)^(1 / d)), d)
     maxcut = maximum([model.rcut for model in model_matrix])
-    neighbour_list = list_type(box, maxcut, N)
+    neighbour_list = list_type(box, maxcut, N; list_parameters=list_parameters)
     species_list = isa(species[1], Integer) ? SpeciesList(species) : nothing
     system = Atoms(position, species, density, energy, temperature, model_matrix, N, d, box, neighbour_list, species_list)
     build_neighbour_list!(system)

src/molecules.jl

@@ -73,7 +73,7 @@ Arguments
 Returns
 - `Molecules` instance with neighbour list built and `energy[1]` set.
 """
-function System(position, species, molecule, density::T, temperature::T, model_matrix, bonds; molecule_species=nothing, list_type=EmptyList) where {T<:AbstractFloat}
+function System(position, species, molecule, density::T, temperature::T, model_matrix, bonds; molecule_species=nothing, list_type=EmptyList, list_parameters=nothing) where {T<:AbstractFloat}
     @assert length(position) == length(species)
     N = length(position)
     Nmol = length(unique(molecule))
@@ -83,7 +83,7 @@ function System(position, species, molecule, density::T, temperature::T, model_m
     box = @SVector fill(T((N / density)^(1 / d)), d)
     energy = zeros(T, 1)
     maxcut = maximum([model.rcut for model in model_matrix])
-    neighbour_list = list_type(box, maxcut, N)
+    neighbour_list = list_type(box, maxcut, N; list_parameters=list_parameters)
     system = Molecules(position, species, molecule, molecule_species,  start_mol, length_mol, density, temperature, energy, model_matrix, d, N, Nmol,box, neighbour_list, bonds)
     build_neighbour_list!(system)
     local_energy = [compute_energy_particle(system, i, neighbour_list) for i in eachindex(position)]

src/neighbours.jl

@@ -18,7 +18,7 @@ struct EmptyList <: NeighbourList end
 
 """Construct an `EmptyList` (ignores box, rcut, N).
 """
-function EmptyList(box, rcut, N)
+function EmptyList(box, rcut, N; list_parameters=nothing)
     return EmptyList()
 end
 
@@ -114,7 +114,7 @@ end
 
 Cells are chosen so that their dimensions are at least `rcut`, and neighbour cells are precomputed.
 """
-function CellList(box::SVector{d,T}, rcut::T, N::Int) where {d,T<:AbstractFloat}
+function CellList(box::SVector{d,T}, rcut::T, N::Int; list_parameters=nothing) where {d,T<:AbstractFloat}
 
     # Calculate cell dimensions ensuring they're >= rcut
     cell_box = @. box / floor(Int, box / rcut)
@@ -233,7 +233,7 @@ end
 
 """Construct a `LinkedList` neighbour list given box, cutoff `rcut`, and number of particles `N`.
 """
-function LinkedList(box, rcut, N)
+function LinkedList(box, rcut, N; list_parameters=nothing)
     cell = box ./ fld.(box, rcut)
     ncells = ntuple(a -> Int(box[a] / cell[a]), length(box))
     head = -ones(Int, prod(ncells))
@@ -375,3 +375,169 @@ function (linked_list::LinkedList)(system::Particles, i::Int)
 
     return LinkedIterator(neighbour_cells, linked_list.head, linked_list.list)
 end
+
+
+# Verlet list
+
+"""Verlet-list neighbour list implementation.
+
+Uses a cell list to find all particles within rcut + dr.
+Keeps an array of neighbour ids for each particle
+"""
+struct VerletList{T<:AbstractFloat, d} <: NeighbourList
+    cs::Vector{Int}
+    box::SVector{d,T}
+    ncells::NTuple{d,Int}
+    cells::Vector{Vector{Int}}  # List of particles in each cell
+    rcut::T
+    dr::T
+    dr2o4::T
+    neighbours::Vector{Vector{Int}}
+    neighbour_cells::Vector{Vector{Int}}  # List of neighbouring cells
+    positions_at_last_update::Vector{SVector{d,T}}
+end
+
+"""Construct a `VerletList` neighbour list given box, cutoff `rcut`, cutoff padding `dr`, and number of particles `N`.
+"""
+function VerletList(box, rcut, N; list_parameters=nothing)
+    if list_parameters == nothing || !haskey(list_parameters, "dr")
+        error("Verlet list must have dr as a parameter")
+    end
+    dr = list_parameters["dr"]
+    cell = box ./ fld.(box, rcut + dr)
+    ncells = ntuple(a -> Int(box[a] / cell[a]), length(box))
+    head = -ones(Int, prod(ncells))
+    list = zeros(Int, N)
+    cs = zeros(Int, N)
+    neighbour_cells = build_neighbour_cells(ncells)
+    neighbours = [Int[] for _ in 1:N]
+    positions_at_last_update = [zeros(SVector{length(box), typeof(box[1])}) for _ in 1:N]
+    cells = [Int[] for _ in 1:prod(ncells)]
+    return VerletList(cs, cell, ncells, cells, rcut, dr, (dr^2)/4, neighbours, neighbour_cells, positions_at_last_update)
+end
+
+"""Calling a VerletList objects return an object which can be iterated upon.
+
+This iteration will return the indices of the neighbours of particle i.
+"""
+function (verlet_list::VerletList)(::Particles, i::Int)
+    return (j for j in verlet_list.neighbours[i])
+end
+
+"""Populate `neighbour_list` (a `VerletList`) by constructing the list of neighbours for each particle.
+
+These neighbours are all particles within a distance `rcut` + `dr`
+"""
+function build_neighbour_list!(system::Particles, neighbour_list::VerletList)
+    # Populate cell list
+    for (i, position_i) in enumerate(system)
+        c = get_cell_index(position_i, neighbour_list)
+        neighbour_list.cs[i] = c
+        push!(neighbour_list.cells[c], i)  # Directly append particle index
+    end
+
+    # Now iterate over all pairs i,j, and see if they are within the cutoff
+    r_cutoff2 = (neighbour_list.rcut + neighbour_list.dr)^2
+    for (i, position_i) in enumerate(system)
+        c = get_cell_index(position_i, neighbour_list)
+        neighbour_cells = neighbour_list.neighbour_cells[c]
+        # Scan the neighbourhood of cell mc (including itself)
+        # and from there scan atoms in cell c2
+        for c2 in neighbour_cells
+            @inbounds for j in neighbour_list.cells[c2]
+                # Don't double count, or self interact
+                if j <= i
+                    continue
+                end
+                position_j = get_position(system, j)
+                r2 = nearest_image_distance_squared(position_i, position_j, get_box(system))
+                if r2 < r_cutoff2
+                    push!(neighbour_list.neighbours[i], j)
+                    push!(neighbour_list.neighbours[j], i)
+                end
+            end
+        end
+
+        neighbour_list.positions_at_last_update[i] = copy(position_i)
+    end
+
+    return nothing
+end
+
+
+"""Update particle `i` moving from cell `c` to cell `c2` in a `VerletList`.
+
+Performs an in-place removal from the old cell and appends to the new one.
+"""
+function update_cell_list!(i, c, c2, neighbour_list::VerletList)
+    # Remove from old cell using in-place filter
+    filter!(x -> x != i, neighbour_list.cells[c])
+    # Add to new cell
+    push!(neighbour_list.cells[c2], i)
+    # Update particle's cell index
+    neighbour_list.cs[i] = c2
+    return nothing
+end
+
+"""Return the old and new cell indices for particle `i` in `neighbour_list` (VerletList).
+
+Used to detect whether a particle has moved between cells.
+"""
+function old_new_cell(system::Particles, i, neighbour_list::VerletList)
+    c = neighbour_list.cs[i]
+    # New cell index
+    mc2 = get_cell(get_position(system, i), neighbour_list)
+    c2 = cell_index(neighbour_list, mc2)
+    return c, c2
+end
+
+"""Update Verlet neighbour list
+
+The cell and list for any given particle is only update if it has moved a distance > dr/2 since the last update"""
+function update_neighbour_list!(system::Particles, i::Int, neighbour_list::VerletList)
+    # Check if particle moved for more than half of dr since last update
+    position_i = get_position(system, i)
+    # WARNING: assumes that positions are never wrapped
+    dx = position_i - neighbour_list.positions_at_last_update[i]
+    if sum(abs2, dx) < neighbour_list.dr2o4
+        return nothing
+    end
+
+    # Need to recompute neighbour list for the particle
+    # Alternatively, we could redo everything and reset all positions, unclear which one is the fastest
+
+    # First, remove i for all other lists
+    @inbounds for j in neighbour_list.neighbours[i]
+        # TODO: benchmark against filter
+        deleteat!(neighbour_list.neighbours[j], neighbour_list.neighbours[j] .== i)
+    end
+    neighbour_list.neighbours[i] = Int[]
+
+    # Update i's cell if needed
+    c, c2 = old_new_cell(system, i, neighbour_list)
+    if c != c2
+        update_cell_list!(i, c, c2, neighbour_list)
+    end
+
+    # Then, recreate list for i
+    r_cutoff2 = (neighbour_list.rcut + neighbour_list.dr)^2
+    c = get_cell_index(position_i, neighbour_list)
+    neighbour_cells = neighbour_list.neighbour_cells[c]
+    for c2 in neighbour_cells
+        @inbounds for j in neighbour_list.cells[c2]
+            if i == j
+                continue
+            end
+            position_j = get_position(system, j)
+            r2 = nearest_image_distance_squared(position_i, position_j, get_box(system))
+            if r2 < r_cutoff2
+                push!(neighbour_list.neighbours[i], j)
+                push!(neighbour_list.neighbours[j], i)
+            end
+        end
+    end
+
+    neighbour_list.positions_at_last_update[i] = copy(position_i)
+
+    return nothing
+end

test/runtests.jl

@@ -14,24 +14,28 @@ using Pkg
     ENV["PATH"] = ENV["PATH"] * path_sep * julia_bin
     @test success(`bash -c "command -v particlesmc"`)
     @test success(`particlesmc params.toml`)
-end 
+end
 
 
 @testset "Potential energy test" begin
     # Test inital configuration
     chains_el = load_chains("config_0.exyz", args=Dict("temperature" => [0.231], "model" => ["JBB"], "list_type" => "EmptyList"))
     chains_ll = load_chains("config_0.lmp", args=Dict("temperature" => [0.231], "model" => ["JBB"], "list_type" => "LinkedList"))
+    chains_vl = load_chains("config_0.lmp", args=Dict("temperature" => [0.231], "model" => ["JBB"], "list_type" => "VerletList", "list_parameters" => Dict("dr" => 0.2)))
     system_el = chains_el[1]
     system_ll = chains_ll[1]
+    system_vl = chains_vl[1]
     @test system_el.N == system_ll.N
     @test system_el.d == system_ll.d
     @test system_el.temperature == system_ll.temperature
     @test all.(system_el.position == system_ll.position)
     @test all.(system_el.species == system_ll.species)
     energy_el = system_el.energy[1] / length(system_el)
     energy_ll = system_ll.energy[1] / length(system_ll)
+    energy_vl = system_vl.energy[1] / length(system_vl)
     @test isapprox(energy_el, -2.676832, atol=1e-6)
     @test isapprox(energy_ll, -2.676832, atol=1e-6)
+    @test isapprox(energy_vl, -2.676832, atol=1e-6)
 
     # Test simulation energy
     M = 1
@@ -63,19 +67,28 @@ end
     path_el = "data/noswap/empty_list/"
     simulation = Simulation(chains_el, algorithm_list, steps; path=path_el, verbose=true)
     run!(simulation)
-    
+
     ## Linked List simulation
     chains_ll = [deepcopy(system_ll)]
     path_ll = "data/noswap/linked_list/"
     simulation = Simulation(chains_ll, algorithm_list, steps; path=path_ll, verbose=true)
     run!(simulation)
 
+    ## Verlet List simulation
+    chains_vl = [deepcopy(system_vl)]
+    path_vl = "data/noswap/verlet_list/"
+    simulation = Simulation(chains_vl, algorithm_list, steps; path=path_vl, verbose=true)
+    run!(simulation)
+
     ## Read energy data and compare
     path_energy_el = joinpath(path_el, "energy.dat")
     path_energy_ll = joinpath(path_ll, "energy.dat")
+    path_energy_vl = joinpath(path_vl, "energy.dat")
     energy_el= readdlm(path_energy_el)[:, 2]
     energy_ll = readdlm(path_energy_ll)[:, 2]
+    energy_vl = readdlm(path_energy_vl)[:, 2]
     @test isapprox(energy_el, energy_ll, atol=1e-6)
+    @test isapprox(energy_ll, energy_vl, atol=1e-6)
 
     # SWAPS
     pswap = 0.8
@@ -100,7 +113,7 @@ end
     path_el = "data/swap/empty_list/"
     simulation = Simulation(chains_el, algorithm_list, steps; path=path_el, verbose=true)
     run!(simulation)
-    
+
     ## Linked List simulation
     chains_ll = [deepcopy(system_ll)]
     path_ll = "data/swap/linked_list/"
@@ -162,7 +175,7 @@ end
     path_el = "data/noswap/empty_list/"
     simulation = Simulation(chains_el, algorithm_list, steps; path=path_el, verbose=true)
     run!(simulation)
-    
+
     ## Linked List simulation
     chains_ll = [deepcopy(system_ll)]
     path_ll = "data/noswap/linked_list/"
@@ -177,5 +190,3 @@ end
     @test isapprox(energy_el, energy_ll, atol=1e-6)
 
 end
-
-