function run(;
    N::Int64,
    T::Float64,
    v::Float64=20.0,
    δt::Float64=2e-5,
    save_at::Float64=0.1,
    framerate::Int64=0,
    save_data::Bool=false,
    n_steps_before_verlet_list_update::Int64=100,
    )

    Random.seed!(42)

    μ = 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 = float(sqrt(N))
    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, save_at / δt)

    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],
        n_frames = floor(Int64, integration_steps / n_steps_before_save) + 1,
    )

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

    name_part = "$(end_time)_T=$(T)_N=$(N)_v=$(v)_dt=$(δt)"

    if save_data
        save_data_jld(sol, args, name_part)
    end

    if framerate > 0
        animate(sol, args, name_part; framerate=framerate)
    end

    return (sol = sol, args = args, name_part = name_part)
end