function run_sim(
    dir::String;
    duration::Float64,
    snapshot_at::Float64=0.1,
    n_steps_before_verlet_list_update::Int64=1000,
    seed::Int64=42,
    n_bundle_snapshots::Int64=100,
)
    @assert length(dir) > 0
    @assert duration > 0
    @assert snapshot_at in 0.01:0.01:duration
    @assert n_steps_before_verlet_list_update >= 0
    @assert seed > 0
    @assert n_bundle_snapshots >= 0

    Random.seed!(seed)

    sim_consts = JSON3.read(read("$dir/sim_consts.json", String))

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

    integration_steps = floor(Int64, duration / sim_consts.δt) + 1
    n_steps_before_snapshot = round(Int64, snapshot_at / sim_consts.δt)

    n_snapshots = floor(Int64, integration_steps / n_steps_before_snapshot) + 1

    n_bundle_snapshots = min(n_snapshots, n_bundle_snapshots)

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

    status = JSON3.read(read("$dir/status.json", String))
    T0::Float64 = status.T

    T = T0 + duration

    start_datetime = now()

    run_vars = (;
        # Input
        duration,
        snapshot_at,
        n_steps_before_verlet_list_update,
        seed,
        n_bundle_snapshots,
        # Calculated
        skin_r,
        integration_steps,
        n_steps_before_snapshot,
        n_snapshots,
        T,
        # Read
        T0,
        start_datetime,
    )

    n_bundles::Int64 = status.n_bundles
    next_bundle = n_bundles + 1
    bundles_dir = "$dir/bundles"
    mkpath(bundles_dir)

    if n_bundle_snapshots > 0
        save_data = true

        open("$bundles_dir/run_vars_$next_bundle.json", "w") do f
            JSON3.pretty(f, run_vars)
        end
    else
        save_data = false
    end

    simulate(
        args,
        sim_consts.δt,
        T0,
        T,
        n_steps_before_verlet_list_update,
        n_steps_before_snapshot,
        n_bundles,
        dir,
        save_data,
    )

    return dir
end