diff --git a/src/ReCo.jl b/src/ReCo.jl
index c0525a2..75b2e95 100644
--- a/src/ReCo.jl
+++ b/src/ReCo.jl
@@ -24,6 +24,5 @@ includet("Particle.jl")
 includet("setup.jl")
 includet("simulation.jl")
 includet("data.jl")
-includet("animation.jl")
-includet("postprocess.jl")
+# includet("animation.jl")
 includet("run.jl")
\ No newline at end of file
diff --git a/src/analysis/pair_correlation_function.jl b/src/analysis/pair_correlation_function.jl
index 95d6b8e..ff9ef3a 100644
--- a/src/analysis/pair_correlation_function.jl
+++ b/src/analysis/pair_correlation_function.jl
@@ -31,7 +31,7 @@ function pair_correlation(sol, variables)
             for j in 1:(variables.N)
                 if i != j
                     for k in 1:n_last_frames
-                        frame = variables.n_frames - k + 1
+                        frame = variables.n_snapshots - k + 1
                         c1 = sol.center[i, frame]
                         c2 = sol.center[j, frame]
 
diff --git a/src/animation.jl b/src/animation.jl
index d463cbf..5572321 100644
--- a/src/animation.jl
+++ b/src/animation.jl
@@ -41,7 +41,7 @@ function animate(directory::String, sol::Solution, variables; framerate=0)
         π2 = 2 * π
         particle_radius = variables.particle_diameter / 2
 
-        @showprogress 0.6 for frame in 1:(variables.n_frames)
+        @showprogress 0.6 for frame in 1:(variables.n_snapshots)
             @simd for i in 1:(variables.N)
                 circles[][i] = Circle(Point2(sol.center[i, frame]), particle_radius)
 
diff --git a/src/data.jl b/src/data.jl
index 61db8b9..d61b915 100644
--- a/src/data.jl
+++ b/src/data.jl
@@ -1,52 +1,60 @@
-struct Solution
+struct Bundle
     t::Vector{Float64}
-    center::Matrix{Vector{Float64}}
+    c::Matrix{Vector{Float64}}
     φ::Matrix{Float64}
-
-    function Solution(N::Int64, n_frames::Int64)
-        return new(
-            zeros(n_frames), [zeros(2) for i in 1:N, j in 1:n_frames], zeros((N, n_frames))
-        )
-    end
 end
 
-function save_frame!(sol, frame, t, particles)
-    frame += 1
+function Bundle(N::Int64, n_snapshots::Int64)
+    return Bundle(
+        zeros(n_snapshots),
+        [zeros(2) for i in 1:N, j in 1:n_snapshots],
+        zeros(N, n_snapshots),
+    )
+end
 
-    sol.t[frame] = t
+function first_n_snapshots(bundle::Bundle, n::Int64)
+    return Bundle(bundle.t[1:n], bundle.c[:, 1:n], bundle.φ[:, 1:n])
+end
+
+function save_snapshot!(
+    bundle::Bundle, n_snapshot::Int64, t::Float64, particles::Vector{Particle}
+)
+    bundle.t[n_snapshot] = t
+
+    bundle_c = bundle.c
+    bundle_φ = bundle.φ
 
     @simd for p in particles
         p_id = p.id
+        p_c = p.c
 
-        pre_pos = sol.center[p_id, frame]
+        bundle_c_id = bundle_c[p_id, n_snapshot]
 
-        @turbo for i in 1:2
-            pre_pos[i] = p.c[i]
+        @turbo for j in 1:2
+            bundle_c_id[j] = p_c[j]
         end
 
-        sol.φ[p_id, frame] = p.φ
+        bundle_φ[p_id, n_snapshot] = p.φ
     end
 
-    return frame
-end
-
-function save_data_jld(directory::String; kwargs...)
-    path = "$directory/data.jld2"
-    JLD2.jldsave(path; kwargs...)
-
-    println("Data saved to $path.")
-
     return nothing
 end
 
-function save_variables(directory::String; kwargs...)
-    path = "$directory/variables.json"
-
-    open(path, "w") do io
-        JSON3.pretty(io, kwargs...)
+function set_status(dir::String, n_bundles::Int64, T::Float64)
+    open("$dir/status.json", "w") do f
+        JSON3.write(f, Dict("n_bundles" => n_bundles, "T" => round(T; digits=3)))
     end
 
-    println("Variables saved to $path.")
+    return nothing
+end
+
+function save_bundle(dir::String, bundle::Bundle, n::Int64, T::Float64)
+    bundle_dir = "$dir/bundles/bundle$n"
+    mkpath(bundle_dir)
+
+    JLD2.save_object("$bundle_dir/bundle.jld2", bundle)
+
+    set_status(dir, n, T)
 
     return nothing
 end
\ No newline at end of file
diff --git a/src/postprocess.jl b/src/postprocess.jl
deleted file mode 100644
index 2f23db9..0000000
--- a/src/postprocess.jl
+++ /dev/null
@@ -1,23 +0,0 @@
-function postprocess(sol, variables, save_data)
-    if !isdir("exports")
-        mkdir("exports")
-    end
-
-    directory = "exports/$(variables.end_time)_T=$(variables.T)_N=$(variables.N)_v=$(variables.v)"
-
-    if save_data || (variables.framerate > 0)
-        mkdir(directory)
-
-        save_variables(directory; variables=variables)
-    end
-
-    if save_data
-        save_data_jld(directory; sol=sol, variables=variables)
-    end
-
-    if variables.framerate > 0
-        animate(directory, sol, variables)
-    end
-
-    return nothing
-end
\ No newline at end of file
diff --git a/src/run.jl b/src/run.jl
index f57777e..cd190c9 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -1,108 +1,110 @@
-function run(;
-    N::Int64,
-    T::Int64,
-    v::Float64,
-    δt::Float64=1e-5,
+function run_sim(
+    dir::String;
+    duration::Float64,
     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,
+    n_bundle_snapshots::Int64=100,
 )
-    @assert N > 0
-    @assert T > 0
-    @assert v >= 0
-    @assert δt > 0
-    @assert snapshot_at ∉ 0.01:0.01:T
-    @assert framerate >= 0
+    @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 packing_ratio > 0
     @assert seed > 0
+    @assert n_bundle_snapshots >= 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) * σ
+    sim_consts = JSON3.read(read("$dir/sim_consts.json", String))
 
-    grid_n = round(Int64, ceil(sqrt(N)))
+    skin_r =
+        1.1 * (
+            2 * sim_consts.v * n_steps_before_verlet_list_update * sim_consts.δt +
+            2 * sim_consts.interaction_r
+        )
 
-    N = grid_n^2
+    integration_steps = floor(Int64, duration / sim_consts.δt) + 1
+    n_steps_before_snapshot = round(Int64, snapshot_at / sim_consts.δt)
 
-    l = sqrt(N * π / packing_ratio) * σ / 4
-    grid_box_width = 2 * l / grid_n
+    n_snapshots = floor(Int64, integration_steps / n_steps_before_snapshot) + 1
 
-    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
+    n_bundle_snapshots = min(n_snapshots, n_bundle_snapshots)
 
     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),
+        v=sim_consts.v,
         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,
+        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),
     )
 
-    sol, end_time = simulate(
-        args, δt, T, n_steps_before_verlet_list_update, n_steps_before_save
-    )
+    status = JSON3.read(read("$dir/status.json", String))
+    T0::Float64 = status.T
 
-    args = nothing
+    T = T0 + duration
 
-    variables = (;
+    start_datetime = now()
+
+    run_vars = (;
         # Input
-        N,
-        T,
-        v,
-        δt,
+        duration,
         snapshot_at,
-        framerate,
         n_steps_before_verlet_list_update,
-        packing_ratio,
-        debug,
+        seed,
+        n_bundle_snapshots,
         # 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,
+        n_steps_before_snapshot,
+        n_snapshots,
+        T,
+        # Read
+        T0,
+        start_datetime,
     )
 
-    postprocess(sol, variables, save_data)
+    n_bundles::Int64 = status.n_bundles
+    next_bundle = n_bundles + 1
+    next_bundle_path = "$dir/bundles/bundle$next_bundle"
+    mkpath(next_bundle_path)
 
-    return (; sol, variables)
-end
+    if n_bundle_snapshots > 0
+        save_data = true
+
+        open("$next_bundle_path/run_vars.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
\ No newline at end of file
diff --git a/src/setup.jl b/src/setup.jl
index 68cfad4..4e26d5a 100644
--- a/src/setup.jl
+++ b/src/setup.jl
@@ -1,14 +1,11 @@
 function initial_particle_grid_pos(i, j; grid_box_width, l)
-    dim1_pos(x) = (x - 0.5) * grid_box_width - l
-
-    return SVector(dim1_pos(i), dim1_pos(j))
+    term = -0.5 * grid_box_width - l
+    return [k * grid_box_width + term for k in (i, j)]
 end
 
 function generate_particles(grid_n, grid_box_width, l)
     particles = Vector{Particle}(undef, grid_n^2)
 
-    particle_pos_in_grid_dim(i) = (i - 0.5) * grid_box_width - l
-
     id = 1
 
     for i in 1:grid_n
@@ -25,3 +22,61 @@ function generate_particles(grid_n, grid_box_width, l)
 
     return particles
 end
+
+function init_sim(;
+    N::Int64, v::Float64, δt::Float64=1e-5, packing_ratio::Float64=0.5, comment::String=""
+)
+    @assert N > 0
+    @assert v >= 0
+    @assert δt in 1e-7:1e-7:1e-4
+    @assert packing_ratio > 0
+
+    μ = 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
+
+    sim_consts = (;
+        # Input
+        N,
+        v,
+        δt,
+        packing_ratio,
+        # Calculated
+        μ,
+        D₀,
+        particle_diameter,
+        Dᵣ,
+        σ,
+        ϵ,
+        interaction_r,
+        grid_n,
+        l,
+        grid_box_width,
+    )
+
+    particles = generate_particles(grid_n, grid_box_width, l)
+    bundle = Bundle(N, 1)
+    save_snapshot!(bundle, 1, 0.0, particles)
+
+    dir = "exports/$(now())_N=$(N)_v=$(v)_#$(rand(1000:9999))$comment"
+    mkpath(dir)
+
+    open("$dir/sim_consts.json", "w") do f
+        JSON3.pretty(f, sim_consts)
+    end
+
+    save_bundle(dir, bundle, 1, 0.0)
+
+    return dir
+end
\ No newline at end of file
diff --git a/src/simulation.jl b/src/simulation.jl
index 64884d1..9e06308 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -66,22 +66,30 @@ end
 function simulate(
     args,
     δt::Float64,
-    T::Int64,
+    T0::Float64,
+    T::Float64,
     n_steps_before_verlet_list_update::Int64,
-    n_steps_before_save::Int64,
+    n_steps_before_snapshot::Int64,
+    n_bundles::Int64,
+    dir::String,
+    save_data::Bool,
 )
-    sol = Solution(args.N, args.n_frames)
-
-    frame = 0
-
-    frame = save_frame!(sol, frame, 0.0, args.particles)
+    bundle_snapshot_counter = 0
 
     start_time = now()
     println("Started simulation at $start_time.")
 
-    @showprogress 0.6 for (integration_step, t) in enumerate(0:δt:T)
-        if integration_step % n_steps_before_save == 0
-            frame = save_frame!(sol, frame, t, args.particles)
+    @showprogress 0.6 for (integration_step, t) in enumerate(T0:δt:T)
+        if (integration_step % n_steps_before_snapshot == 0) && save_data
+            bundle_snapshot_counter += 1
+            save_snapshot!(args.bundle, bundle_snapshot_counter, t, args.particles)
+
+            if bundle_snapshot_counter == args.n_bundle_snapshots
+                bundle_snapshot_counter = 0
+                n_bundles += 1
+
+                save_bundle(dir, args.bundle, n_bundles, t)
+            end
         end
 
         if integration_step % n_steps_before_verlet_list_update == 0
@@ -91,9 +99,17 @@ function simulate(
         euler!(args)
     end
 
+    if bundle_snapshot_counter > 0
+        bundle = first_n_snapshots(args.bundle, bundle_snapshot_counter)
+
+        n_bundles += 1
+
+        save_bundle(dir, bundle, n_bundles, T)
+    end
+
     end_time = now()
     elapsed_time = canonicalize(CompoundPeriod(end_time - start_time))
     println("Simulation done at $end_time and took $elapsed_time.")
 
-    return (; sol, end_time)
+    return nothing
 end