diff --git a/Manifest.toml b/Manifest.toml
index d319c1e..9f9cef4 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -294,9 +294,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[deps.Distributions]]
 deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
-git-tree-sha1 = "24d26ca2197c158304ab2329af074fbe14c988e4"
+git-tree-sha1 = "2e97190dfd4382499a4ac349e8d316491c9db341"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.45"
+version = "0.25.46"
 
 [[deps.DocStringExtensions]]
 deps = ["LibGit2"]
@@ -328,9 +328,9 @@ version = "1.3.0"
 
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
+git-tree-sha1 = "ae13fcbc7ab8f16b0856729b050ef0c446aa3492"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.2.10+0"
+version = "2.4.4+0"
 
 [[deps.ExprTools]]
 git-tree-sha1 = "56559bbef6ca5ea0c0818fa5c90320398a6fbf8d"
@@ -363,9 +363,9 @@ version = "3.3.10+0"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "67551df041955cc6ee2ed098718c8fcd7fc7aebe"
+git-tree-sha1 = "80ced645013a5dbdc52cf70329399c35ce007fae"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.12.0"
+version = "1.13.0"
 
 [[deps.FilePathsBase]]
 deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
@@ -579,10 +579,10 @@ uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 version = "0.9.3"
 
 [[deps.ImageIO]]
-deps = ["FileIO", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"]
-git-tree-sha1 = "816fc866edd8307a6e79a575e6585bfab8cef27f"
+deps = ["FileIO", "JpegTurbo", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"]
+git-tree-sha1 = "464bdef044df52e6436f8c018bea2d48c40bb27b"
 uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
-version = "0.6.0"
+version = "0.6.1"
 
 [[deps.Imath_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -668,9 +668,9 @@ version = "1.0.0"
 
 [[deps.JLD2]]
 deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
-git-tree-sha1 = "bcb31db46795eeb64480c89d854615bc78a13289"
+git-tree-sha1 = "b528d68220e2aba1d2d0c0461b6f7eda8c5c1e33"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.4.19"
+version = "0.4.20"
 
 [[deps.JLLWrappers]]
 deps = ["Preferences"]
@@ -690,11 +690,17 @@ git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
 uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
 version = "1.9.2"
 
+[[deps.JpegTurbo]]
+deps = ["CEnum", "FileIO", "ImageCore", "JpegTurbo_jll", "TOML"]
+git-tree-sha1 = "a77b273f1ddec645d1b7c4fd5fb98c8f90ad10a5"
+uuid = "b835a17e-a41a-41e7-81f0-2f016b05efe0"
+version = "0.1.1"
+
 [[deps.JpegTurbo_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
+git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
-version = "2.1.0+0"
+version = "2.1.2+0"
 
 [[deps.Juno]]
 deps = ["Base64", "Logging", "Media", "Profile"]
@@ -882,9 +888,9 @@ uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 version = "0.2.1"
 
 [[deps.ManualMemory]]
-git-tree-sha1 = "9cb207b18148b2199db259adfa923b45593fe08e"
+git-tree-sha1 = "bcaef4fc7a0cfe2cba636d84cda54b5e4e4ca3cd"
 uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667"
-version = "0.1.6"
+version = "0.1.8"
 
 [[deps.MappedArrays]]
 git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142"
@@ -924,9 +930,9 @@ version = "0.4.10"
 
 [[deps.MethodAnalysis]]
 deps = ["AbstractTrees"]
-git-tree-sha1 = "40c1181bf7943b176c4a11edd67e72ab81fa3b1d"
+git-tree-sha1 = "81e123ea81d6081fe8b733dbe79e1291d55cfb0f"
 uuid = "85b6ec6f-f7df-4429-9514-a64bcd9ee824"
-version = "0.4.4"
+version = "0.4.6"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
@@ -1058,9 +1064,9 @@ version = "0.11.5"
 
 [[deps.PNGFiles]]
 deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
-git-tree-sha1 = "6d105d40e30b635cfed9d52ec29cf456e27d38f8"
+git-tree-sha1 = "eb4dbb8139f6125471aa3da98fb70f02dc58e49c"
 uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
-version = "0.3.12"
+version = "0.3.14"
 
 [[deps.Packing]]
 deps = ["GeometryBasics"]
@@ -1076,9 +1082,9 @@ version = "0.5.11"
 
 [[deps.Pango_jll]]
 deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "9bc1871464b12ed19297fbc56c4fb4ba84988b0d"
+git-tree-sha1 = "3a121dfbba67c94a5bec9dde613c3d0cbcf3a12b"
 uuid = "36c8627f-9965-5494-a995-c6b170f724f3"
-version = "1.47.0+0"
+version = "1.50.3+0"
 
 [[deps.Parameters]]
 deps = ["OrderedCollections", "UnPack"]
@@ -1153,9 +1159,9 @@ uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 
 [[deps.ProfileView]]
 deps = ["Cairo", "Colors", "FileIO", "FlameGraphs", "Graphics", "Gtk", "GtkObservables", "InteractiveUtils", "IntervalSets", "MethodAnalysis", "Preferences", "Profile", "UUIDs"]
-git-tree-sha1 = "93b795da19d3668ff3ba5463481d137ba7fb3c2a"
+git-tree-sha1 = "1638391c64d5ae2853b54823df3f98150dc202f3"
 uuid = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
-version = "1.4.0"
+version = "1.5.0"
 
 [[deps.ProgressMeter]]
 deps = ["Distributed", "Printf"]
@@ -1377,9 +1383,9 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [[deps.SpecialFunctions]]
 deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "e6bf188613555c78062842777b116905a9f9dd49"
+git-tree-sha1 = "a4116accb1c84f0a8e1b9932d873654942b2364b"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.1.0"
+version = "2.1.1"
 
 [[deps.StackViews]]
 deps = ["OffsetArrays"]
@@ -1395,9 +1401,9 @@ version = "0.4.1"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "2884859916598f974858ff01df7dfc6c708dd895"
+git-tree-sha1 = "a635a9333989a094bddc9f940c04c549cd66afcf"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.3.3"
+version = "1.3.4"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
diff --git a/README.adoc b/README.adoc
index b337e71..f93ca01 100644
--- a/README.adoc
+++ b/README.adoc
@@ -31,6 +31,7 @@ Pkg.activate(".")
 ----
 
 === Install dependencies
+
 After activating the package environment, run the follwing to install the package dependencies:
 
 [source, julia]
@@ -39,7 +40,31 @@ Pkg.instantiate()
 ----
 
 == Run simulation
-// TODO
+
+Import the package:
+
+[source, julia]
+----
+using ReCo
+----
+
+It might take a moment to precompile.
+
+Initialize the simulation with 100 particles having a self-propulsion velocity of 40.0 and return the relative path to the simulation directory:
+
+[source, julia]
+----
+sim_dir = init_sim(100, 40.0)
+----
+
+Run the simulation:
+
+[source, julia]
+----
+run_sim(sim_dir, duration=20.0)
+----
+
+The values for the number of particles, self-propulsion velocity and simulation duration are used here as an example. For more information about possible values and other optional arguments, press `?` in the REPL after running `using ReCo`. Then type `init_sim` or `run_sim` followed by pressing enter. This will show the method's documention.
 
 == Run reinforcement learning
 // TODO
diff --git a/analysis/mean_squared_displacement.jl b/analysis/mean_squared_displacement.jl
index 4a62ed3..4a9fe66 100644
--- a/analysis/mean_squared_displacement.jl
+++ b/analysis/mean_squared_displacement.jl
@@ -41,7 +41,7 @@ function msd_simulation(
 )
     Random.seed!(seed)
 
-    dir = ReCo.init_sim(;
+    sim_dir = ReCo.init_sim(;
         n_particles=1,
         v₀=v₀,
         parent_dir=parent_dir,
@@ -51,14 +51,14 @@ function msd_simulation(
     )
 
     ReCo.run_sim(
-        dir;
+        sim_dir;
         duration=T,
         seed=rand(1:typemax(Int64)),
         snapshot_at=snapshot_at,
         n_bundle_snapshots=1000,
     )
 
-    return dir
+    return sim_dir
 end
 
 function mean_squared_displacement(;
@@ -78,9 +78,11 @@ function mean_squared_displacement(;
         parent_dir = main_parent_dir * "/$v₀"
 
         Threads.@threads for sim_ind in 1:n_simulations
-            dir = msd_simulation(v₀, half_box_len, T, 0.5, parent_dir; comment="$sim_ind")
+            sim_dir = msd_simulation(
+                v₀, half_box_len, T, 0.5, parent_dir; comment="$sim_ind"
+            )
 
-            sim_dirs[sim_ind, v₀_ind] = dir
+            sim_dirs[sim_ind, v₀_ind] = sim_dir
 
             ProgressMeter.next!(progress; showvalues=[(:v₀, v₀)])
         end
@@ -175,12 +177,14 @@ end
 
 function plot_random_walk(; T::Float64, v₀::Float64, seed::Int64)
     half_box_len = max_possible_displacement(T, v₀)
-    dir = msd_simulation(v₀, half_box_len, T, 8.0, "random_walk_$(Dates.now())"; seed=seed)
+    sim_dir = msd_simulation(
+        v₀, half_box_len, T, 8.0, "random_walk_$(Dates.now())"; seed=seed
+    )
 
     ts = Float64[]
     cs = SVector{2,Float64}[]
     ReCo.append_bundle_properties!(
-        (ts, cs), (:t, :c), dir; particle_slice=1, snapshot_slice=:
+        (ts, cs), (:t, :c), sim_dir; particle_slice=1, snapshot_slice=:
     )
 
     init_cairomakie!()
diff --git a/analysis/radial_distribution_function/radial_distribution_function.jl b/analysis/radial_distribution_function/radial_distribution_function.jl
index 629861e..beadcbd 100644
--- a/analysis/radial_distribution_function/radial_distribution_function.jl
+++ b/analysis/radial_distribution_function/radial_distribution_function.jl
@@ -12,7 +12,7 @@ using ReCo: ReCo
 includet("../../src/Visualization/common_CairoMakie.jl")
 
 function radial_distribution_simulation(;
-    n_particles::Int64, v₀s::NTuple{N,Float64}, T::Float64, packing_ratio::Float64
+    n_particles::Int64, v₀s::NTuple{N,Float64}, T::Float64, packing_fraction::Float64
 ) where {N}
     Random.seed!(42)
 
@@ -25,17 +25,17 @@ function radial_distribution_simulation(;
     Threads.@threads for v₀_ind in 1:n_v₀s
         v₀ = v₀s[v₀_ind]
 
-        dir = ReCo.init_sim(;
+        sim_dir = ReCo.init_sim(;
             n_particles=n_particles,
             v₀=v₀,
-            packing_ratio=packing_ratio,
+            packing_fraction=packing_fraction,
             parent_dir=parent_dir,
             comment="$v₀",
         )
 
-        ReCo.run_sim(dir; duration=T, seed=v₀_ind)
+        ReCo.run_sim(sim_dir; duration=T, seed=v₀_ind)
 
-        sim_dirs[v₀_ind] = dir
+        sim_dirs[v₀_ind] = sim_dir
     end
 
     return sim_dirs
@@ -191,7 +191,7 @@ function plot_radial_distributions(;
 
     if include_comparison
         comparison_curve = CSV.read(
-            "analysis/radial_distribution_function/g_of_r_d_ratio_with_0_37_packing_ratio.csv",
+            "analysis/radial_distribution_function/g_of_r_d_ratio_with_0_37_packing_fraction.csv",
             DataFrames.DataFrame;
             header=3,
         )
@@ -217,7 +217,7 @@ function run_radial_distribution_analysis()
     n_particles = 1000
 
     sim_dirs = radial_distribution_simulation(;
-        n_particles=n_particles, v₀s=v₀s, T=100.0, packing_ratio=0.37
+        n_particles=n_particles, v₀s=v₀s, T=100.0, packing_fraction=0.37
     )
 
     lower_radii, gs, particle_radius = radial_distribution(;
diff --git a/src/RL/RL.jl b/src/RL/RL.jl
index c24a447..34ff165 100644
--- a/src/RL/RL.jl
+++ b/src/RL/RL.jl
@@ -74,7 +74,7 @@ function run_rl(;
     seed::Int64=42,
     ϵ_stable::Float64=0.00001,
     skin_to_interaction_radius_ratio::Float64=ReCo.DEFAULT_SKIN_TO_INTERACTION_RADIUS_RATIO,
-    packing_ratio::Float64=0.15,
+    packing_fraction::Float64=0.15,
     show_progress::Bool=true,
     reward_discount::Float64=0.1,
     show_simulation_progress::Bool=true,
@@ -96,7 +96,7 @@ function run_rl(;
         n_particles,
         0.0;
         skin_to_interaction_radius_ratio=skin_to_interaction_radius_ratio,
-        packing_ratio=packing_ratio,
+        packing_fraction=packing_fraction,
     )
     n_particles = sim_consts.n_particles # Not always equal to the input!
 
@@ -137,7 +137,7 @@ function run_rl(;
     progress = ProgressMeter.Progress(n_episodes; dt=2, enabled=show_progress, desc="RL: ")
 
     for episode in 1:n_episodes
-        dir = ReCo.init_sim_with_sim_consts(sim_consts; parent_dir=parent_dir)
+        sim_dir = ReCo.init_sim_with_sim_consts(sim_consts; parent_dir=parent_dir)
 
         # Reset
         reset!(env)
@@ -148,7 +148,7 @@ function run_rl(;
 
         # Episode
         ReCo.run_sim(
-            dir;
+            sim_dir;
             duration=episode_duration,
             seed=episode,
             env_helper=env_helper,
diff --git a/src/data.jl b/src/data.jl
index b425f57..c04cc53 100644
--- a/src/data.jl
+++ b/src/data.jl
@@ -2,7 +2,7 @@ struct SimConsts
     n_particles::Int64
     v₀::Float64
     δt::Float64
-    packing_ratio::Float64
+    packing_fraction::Float64
     μₜ::Float64
     Dₜ::Float64
     particle_radius::Float64
@@ -72,19 +72,19 @@ function save_snapshot!(
     return nothing
 end
 
-function save_bundle(dir::String, bundle::Bundle, n_bundle::Int64, T::Float64)
-    bundles_dir = "$dir/bundles"
+function save_bundle(sim_dir::String, bundle::Bundle, n_bundle::Int64, T::Float64)
+    bundles_dir = "$sim_dir/bundles"
     mkpath(bundles_dir)
 
     JLD2.save_object("$bundles_dir/$n_bundle.jld2", bundle)
 
-    JLD2.save_object("$dir/sim_state.jld2", SimState(n_bundle, round(T; digits=3)))
+    JLD2.save_object("$sim_dir/sim_state.jld2", SimState(n_bundle, round(T; digits=3)))
 
     return nothing
 end
 
-function sorted_bundle_paths(dir::String; rev::Bool=false)
-    bundles_dir = "$dir/bundles"
+function sorted_bundle_paths(sim_dir::String; rev::Bool=false)
+    bundles_dir = "$sim_dir/bundles"
     bundle_paths = readdir(bundles_dir; join=true, sort=false)
 
     n_bundles = length(bundle_paths)
@@ -160,8 +160,8 @@ function BundlesInfo(sim_dir::String)
     return BundlesInfo(n_bundles, total_n_snapshots, bundle_n_snapshots, bundle_paths)
 end
 
-function load_sim_consts(dir::String)
-    sim_consts::ReCo.SimConsts = JLD2.load_object("$dir/sim_consts.jld2")
+function load_sim_consts(sim_dir::String)
+    sim_consts::ReCo.SimConsts = JLD2.load_object("$sim_dir/sim_consts.jld2")
 
     return sim_consts
 end
diff --git a/src/run.jl b/src/run.jl
index efa49d7..8456582 100644
--- a/src/run.jl
+++ b/src/run.jl
@@ -1,3 +1,8 @@
+const DEFAULT_SNAPSHOT_AT = 0.1
+const DEFAULT_SEED = 42
+const DEFAULT_N_BUNDLE_SNAPSHOTS = 100
+const DEFAULT_SHOW_PROGRESS = true
+
 function empty_hook(args...)
     return nothing
 end
@@ -6,16 +11,34 @@ function gen_cell_list_box(half_box_len::Float64, skin_radius::Float64)
     return CellListMap.Box(SVector(2 * half_box_len, 2 * half_box_len), skin_radius)
 end
 
+"""
+    run_sim(sim_dir::String; <keyword arguments>)
+
+Run the initialized simulation in its directory `sim_dir`.
+
+This method starts or resumes a simulation. For long simulations, the simulation can be stopped by pressing `Ctrl` + `c`. This stopped simulation can be resumed later by running this method again with the same simulation directory.
+
+Some of the last snapshots might be lost if the simulations is stopped (see the argument `n_bundle_snapshots`).
+
+# Arguments
+- `sim_dir::String`: Relative path of the initialized simulation directory.
+- `duration::Float64`: Duration of the simulation.
+- `snapshot_at::Float64=$DEFAULT_SNAPSHOT_AT`: Snapshot time interval.
+- `seed::Int64=$DEFAULT_SEED`: Random number generator seed.
+- `n_bundle_snapshots::Int64=$DEFAULT_N_BUNDLE_SNAPSHOTS`: Number of snapshots in a bundle. This number is relevant for long simulations that can be stopped while running. A simulation can be continued from the last bundle of snapshots. If the number of snapshots in a bundle is too high and the simulation is stopped, many of the last snapshots can be lost. A low number results in high IO since snapshots are then bundled and stored more often. For example, setting this number to 1 results in saving every snapshot immidiately without bundeling it with other snapshots which would be more efficient. Setting the number to 1000 could mean loosing 999 snapshots in the worst case if the simulation is stopped before having 1000 snapshots to bundle and save.
+- `env_helper::Union{RL.EnvHelper,Nothing}=nothing`: Environment helper. It should be left as the default `nothing` unless this method is used internally for reinforcement learning.
+- `show_progress::Bool=$DEFAULT_SHOW_PROGRESS`: Show simulation progress bar.
+"""
 function run_sim(
-    dir::String;
+    sim_dir::String;
     duration::Float64,
-    snapshot_at::Float64=0.1,
-    seed::Int64=42,
-    n_bundle_snapshots::Int64=100,
+    snapshot_at::Float64=DEFAULT_SNAPSHOT_AT,
+    seed::Int64=DEFAULT_SEED,
+    n_bundle_snapshots::Int64=DEFAULT_N_BUNDLE_SNAPSHOTS,
     env_helper::Union{RL.EnvHelper,Nothing}=nothing,
-    show_progress::Bool=true,
+    show_progress::Bool=DEFAULT_SHOW_PROGRESS,
 )
-    @assert length(dir) > 0
+    @assert length(sim_dir) > 0
     @assert duration > 0
     @assert snapshot_at in 0.001:0.001:duration
     @assert seed > 0
@@ -23,7 +46,7 @@ function run_sim(
 
     Random.seed!(seed)
 
-    sim_consts = load_sim_consts(dir)
+    sim_consts = load_sim_consts(sim_dir)
 
     integration_steps = floor(Int64, duration / sim_consts.δt) + 1
 
@@ -34,7 +57,7 @@ function run_sim(
 
     n_bundle_snapshots = min(n_snapshots, n_bundle_snapshots)
 
-    sim_state::SimState = JLD2.load_object("$dir/sim_state.jld2")
+    sim_state::SimState = JLD2.load_object("$sim_dir/sim_state.jld2")
     n_bundles = sim_state.n_bundles
 
     T0::Float64 = sim_state.T
@@ -66,14 +89,14 @@ function run_sim(
 
         next_bundle = n_bundles + 1
 
-        runs_dir = "$dir/runs"
+        runs_dir = "$sim_dir/runs"
 
         if save_data
             JLD2.save_object("$runs_dir/run_params_$next_bundle.jld2", run_params)
         end
     end
 
-    bundles_dir = "$dir/bundles"
+    bundles_dir = "$sim_dir/bundles"
     bundle::Bundle = JLD2.load_object("$bundles_dir/$n_bundles.jld2")
 
     particles = gen_particles(bundle, sim_consts.n_particles)
@@ -113,7 +136,7 @@ function run_sim(
         sim_consts.n_steps_before_verlet_list_update,
         n_steps_before_snapshot,
         n_bundles,
-        dir,
+        sim_dir,
         save_data,
         env_helper,
     )
diff --git a/src/setup.jl b/src/setup.jl
index f7f1d4a..3c0e010 100644
--- a/src/setup.jl
+++ b/src/setup.jl
@@ -1,4 +1,4 @@
-const DEFAULT_PACKING_RATIO = 0.5
+const DEFAULT_PACKING_FRACTION = 0.5
 const DEFAULT_δt = 1e-5
 const DEFAULT_SKIN_TO_INTERACTION_RADIUS_RATIO = 2.5
 const DEFAULT_EXPORTS_DIR = "exports"
@@ -8,6 +8,7 @@ const DEFAULT_PARTICLE_RADIUS = 0.5
 const DEFAULT_Dₜ = 1.0
 const DEFAULT_μₜ = 1.0
 const DEFAULT_ϵ = 100.0
+const DEFAULT_HALF_BOX_LEN = 0.0
 
 function initial_particle_grid_pos(
     i::Int64, j::Int64, grid_box_width::Float64, half_box_len::Float64
@@ -50,14 +51,14 @@ function gen_sim_consts(
     n_particles::Int64,
     v₀::Float64;
     δt::Float64=DEFAULT_δt,
-    packing_ratio::Float64=DEFAULT_PACKING_RATIO,
+    packing_fraction::Float64=DEFAULT_PACKING_FRACTION,
     skin_to_interaction_radius_ratio::Float64=DEFAULT_SKIN_TO_INTERACTION_RADIUS_RATIO,
     half_box_len::Float64=0.0,
 )
     @assert n_particles > 0
     @assert v₀ >= 0
     @assert δt in 1e-7:1e-7:1e-4
-    @assert packing_ratio > 0
+    @assert packing_fraction > 0
 
     μₜ = DEFAULT_μₜ
     Dₜ = DEFAULT_Dₜ
@@ -87,9 +88,11 @@ function gen_sim_consts(
     n_particles = grid_n^2
 
     if half_box_len == 0.0
-        half_box_len = sqrt(n_particles * π / packing_ratio) * σ / 4
-    elseif packing_ratio != DEFAULT_PACKING_RATIO
-        error("You can not specify half_box_len and packing_ratio at the same time!")
+        half_box_len = sqrt(n_particles * π / packing_fraction) * σ / 4
+    elseif packing_fraction != DEFAULT_PACKING_FRACTION
+        error("You can not specify half_box_len and packing_fraction at the same time!")
+    else
+        @assert half_box_len > 0.0
     end
 
     grid_box_width = 2 * half_box_len / grid_n
@@ -99,7 +102,7 @@ function gen_sim_consts(
         n_particles,
         v₀,
         δt,
-        packing_ratio,
+        packing_fraction,
         # Internal
         μₜ,
         Dₜ,
@@ -128,49 +131,65 @@ function init_sim_with_sim_consts(
     bundle = Bundle(sim_consts.n_particles, 1)
     save_snapshot!(bundle, 1, 0.0, particles)
 
-    dir = exports_dir
+    sim_dir = exports_dir
 
     if length(parent_dir) > 0
-        dir *= "/$parent_dir"
+        sim_dir *= "/$parent_dir"
     end
 
     start_datetime = Dates.now()
-    dir *= "/$(start_datetime)_N=$(sim_consts.n_particles)_v=$(sim_consts.v₀)_#$(rand(1000:9999))"
+    sim_dir *= "/$(start_datetime)_N=$(sim_consts.n_particles)_v=$(sim_consts.v₀)_#$(rand(1000:9999))"
 
     if length(comment) > 0
-        dir *= "_$comment"
+        sim_dir *= "_$comment"
     end
 
-    mkpath(dir)
+    mkpath(sim_dir)
 
-    task = @async JLD2.save_object("$dir/sim_consts.jld2", sim_consts)
+    task = @async JLD2.save_object("$sim_dir/sim_consts.jld2", sim_consts)
 
-    save_bundle(dir, bundle, 1, 0.0)
+    save_bundle(sim_dir, bundle, 1, 0.0)
 
-    runs_dir = "$dir/runs"
+    runs_dir = "$sim_dir/runs"
     mkpath(runs_dir)
 
     wait(task)
 
-    return dir
+    return sim_dir
 end
 
+"""
+    init_sim(n_particles::Int64, v₀::Float64; <keyword arguments>)
+
+Initialize simulation and return the relative path of the simulation directory.
+
+# Arguments
+- `n_particles::Int64`: Number of particles.
+- `v₀::Float64`: Self-propulsion velocity. Only values in the interval [0.0, 80.0] are tested.
+- `δt::Float64=$DEFAULT_δt`: Integration time step.
+- `packing_fraction::Float64=$DEFAULT_PACKING_FRACTION`: Packing fraction.
+- `skin_to_interaction_radius_ratio::Float64=$DEFAULT_SKIN_TO_INTERACTION_RADIUS_RATIO`: Ratio of skin radius to interaction radius.
+- `exports_dir::String="$DEFAULT_EXPORTS_DIR"`: Relative path to exports directory.
+- `parent_dir::String="$DEFAULT_PARENT_DIR"`: Name of the directory within the simulation directory is placed.
+- `comment::String="$DEFAULT_COMMENT"`: Comment to append to the directory name.
+- `half_box_len::Float64=$DEFAULT_HALF_BOX_LEN` Half box length. The default of 0.0 means that the half box length will be calculated from the packing fraction. Otherwise, the provided half box length will be used. It is not possible to provide a half box length and a packing fraction at the same time.
+"""
 function init_sim(;
     n_particles::Int64,
     v₀::Float64,
     δt::Float64=DEFAULT_δt,
-    packing_ratio::Float64=DEFAULT_PACKING_RATIO,
+    packing_fraction::Float64=DEFAULT_PACKING_FRACTION,
     skin_to_interaction_radius_ratio::Float64=DEFAULT_SKIN_TO_INTERACTION_RADIUS_RATIO,
     exports_dir::String=DEFAULT_EXPORTS_DIR,
     parent_dir::String=DEFAULT_PARENT_DIR,
     comment::String=DEFAULT_COMMENT,
-    half_box_len::Float64=0.0,
+    half_box_len::Float64=DEFAULT_HALF_BOX_LEN,
 )
     sim_consts = gen_sim_consts(
         n_particles,
         v₀;
         δt=δt,
-        packing_ratio=packing_ratio,
+        packing_fraction=packing_fraction,
         skin_to_interaction_radius_ratio=skin_to_interaction_radius_ratio,
         half_box_len,
     )
diff --git a/src/simulation.jl b/src/simulation.jl
index 4fd6213..1c93479 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -114,7 +114,7 @@ function simulate!(
     n_steps_before_verlet_list_update::Int64,
     n_steps_before_snapshot::Int64,
     n_bundles::Int64,
-    dir::String,
+    sim_dir::String,
     save_data::Bool,
     env_helper::Union{RL.EnvHelper,Nothing},
 )
@@ -148,7 +148,7 @@ function simulate!(
                     bundle_snapshot_counter = 0
                     n_bundles += 1
 
-                    save_bundle(dir, args.bundle, n_bundles, t)
+                    save_bundle(sim_dir, args.bundle, n_bundles, t)
                 end
             end
         end
@@ -193,7 +193,7 @@ function simulate!(
 
         n_bundles += 1
 
-        save_bundle(dir, bundle, n_bundles, T)
+        save_bundle(sim_dir, bundle, n_bundles, T)
     end
 
     return nothing