function run(;
    N::Int64,
    T::Int64,
    v::Float64,
    δt::Float64=1e-5,
    snapshot_at::Float64=0.1,
    framerate::Int64=0,
    save_data::Bool=true,
    n_steps_before_verlet_list_update::Int64=1000,
    packing_ratio::Float64=0.5,
    debug::Bool=false,
    seed::Int64=42,
)
    @assert N > 0
    @assert T > 0
    @assert v >= 0
    @assert δt > 0
    @assert snapshot_at ∉ 0.01:0.01:T
    @assert framerate >= 0
    @assert n_steps_before_verlet_list_update >= 0
    @assert packing_ratio > 0
    @assert seed > 0

    Random.seed!(seed)

    μ = 1.0
    D₀ = 1.0
    particle_diameter = 1.0
    Dᵣ = 3 * D₀ / (particle_diameter^2)
    σ = 1.0
    ϵ = 100.0
    interaction_r = 2^(1 / 6) * σ

    grid_n = round(Int64, ceil(sqrt(N)))

    N = grid_n^2

    l = sqrt(N * π / packing_ratio) * σ / 4
    grid_box_width = 2 * l / grid_n

    skin_r = 1.1 * (2 * v * n_steps_before_verlet_list_update * δt + 2 * interaction_r)

    integration_steps = floor(Int64, T / δt) + 1
    n_steps_before_save = round(Int64, snapshot_at / δt)

    n_frames = floor(Int64, integration_steps / n_steps_before_save) + 1

    args = (
        v=v,
        c₁=4 * ϵ * 6 * σ^6 * δt * μ,
        c₂=2 * σ^6,
        c₃=sqrt(2 * D₀ * δt),
        c₄=sqrt(2 * Dᵣ * δt),
        vδt=v * δt,
        μ=μ,
        interaction_r=interaction_r,
        interaction_r²=interaction_r^2,
        N=N,
        l=l,
        particle_diameter=particle_diameter,
        particles=generate_particles(grid_n, grid_box_width, l),
        skin_r=skin_r,
        skin_r²=skin_r^2,
        verlet_list=[PreVector(Int64, N - 1) for i in 1:(N - 1)],
        n_frames=n_frames,
        debug=debug,
    )

    sol, end_time = simulate(
        args, δt, T, n_steps_before_verlet_list_update, n_steps_before_save
    )

    args = nothing

    variables = (;
        # Input
        N,
        T,
        v,
        δt,
        snapshot_at,
        framerate,
        n_steps_before_verlet_list_update,
        packing_ratio,
        debug,
        # Calculated
        μ,
        D₀,
        particle_diameter,
        Dᵣ,
        σ,
        ϵ,
        interaction_r,
        grid_n,
        l,
        grid_box_width,
        skin_r,
        integration_steps,
        n_steps_before_save,
        n_frames,
        # Output
        end_time,
    )

    postprocess(sol, variables, save_data)

    return (; sol, variables)
end