diff --git a/src/IO/IO.jl b/src/IO/IO.jl
index e4f6a43..99fa477 100644
--- a/src/IO/IO.jl
+++ b/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,6 +286,7 @@ 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])
@@ -293,9 +294,9 @@ function load_chains(init_path; args=Dict(), verbose=false)
         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
diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl
index 4199345..f4774c0 100644
--- a/src/ParticlesMC.jl
+++ b/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 = []
diff --git a/src/atoms.jl b/src/atoms.jl
index f528a4c..04dcf55 100644
--- a/src/atoms.jl
+++ b/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)
diff --git a/src/molecules.jl b/src/molecules.jl
index f174424..fe4a93c 100644
--- a/src/molecules.jl
+++ b/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)]
diff --git a/src/neighbours.jl b/src/neighbours.jl
index c737ba1..d77758a 100644
--- a/src/neighbours.jl
+++ b/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
diff --git a/test/runtests.jl b/test/runtests.jl
index 2b0f915..995c261 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -14,15 +14,17 @@ 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
@@ -30,8 +32,10 @@ end
     @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
-
-