Added radial distribution function

analysis/mean_squared_displacement.jl

@@ -9,6 +9,8 @@ using ProgressMeter: ProgressMeter
 
 using ReCo: ReCo
 
+includet("../graphics/common_CairoMakie.jl")
+
 # IMPORTANT: Disable the periodic boundary conditions
 # The arguments types have to match for the function to be overwritten!
 ReCo.push_to_verlet_list!(::Any, ::Any, ::Any) = nothing
@@ -26,20 +28,6 @@ ReCo.minimum_image(v::SVector{2,Float64}, ::Float64) = v
 const δt = 1e-4
 const Dₜ = ReCo.DEFAULT_Dₜ
 
-function fill_with_bundle_property!(
-    v::Vector, property::Symbol, sim_dir::String, first_bundle::Int64=1
-)
-    bundle_paths = ReCo.sorted_bundle_paths(sim_dir)
-
-    for i in first_bundle:length(bundle_paths)
-        bundle::ReCo.Bundle = JLD2.load_object(bundle_paths[i])
-
-        append!(v, getproperty(bundle, property))
-    end
-
-    return nothing
-end
-
 function max_possible_displacement(T::Float64, v₀::Float64, δt::Float64=δt, Dₜ::Float64=Dₜ)
     return T * v₀ + T / δt * sqrt(2 * Dₜ * δt)
 end
@@ -52,6 +40,8 @@ function msd_simulation(
     parent_dir::String,
     comment::String="",
 )
+    Random.seed!(42)
+
     dir = ReCo.init_sim(;
         n_particles=1,
         v₀=v₀,
@@ -72,10 +62,8 @@ function msd_simulation(
 end
 
 function mean_squared_displacement(;
-    n_simulations::Int64, v₀s::AbstractVector{Float64}, T::Float64
-)
-    Random.seed!(42)
-
+    n_simulations::Int64, v₀s::NTuple{N,Float64}, T::Float64
+) where {N}
     n_v₀s = length(v₀s)
 
     main_parent_dir = "mean_squared_displacement_$(Dates.now())"
@@ -99,7 +87,9 @@ function mean_squared_displacement(;
     end
 
     ts = Float64[]
-    fill_with_bundle_property!(ts, :t, sim_dirs[1, 1], 2) # Skip the first bundle to avoid t = 0
+    ReCo.append_bundle_properties!(
+        (ts,), (:t,), sim_dirs[1, 1]; particle_slice=1, snapshot_slice=:, first_bundle=2
+    ) # Skip the first bundle to avoid t = 0
 
     mean_sq_displacements = zeros((length(ts), n_v₀s))
 
@@ -137,42 +127,9 @@ function expected_mean_squared_displacement(t::Float64, v₀::Float64)
     return (4 * Dₜ + 2 * v₀^2 / Dᵣ) * t + 2 * v₀^2 * (exp(-Dᵣ * t) - 1) / (Dᵣ^2)
 end
 
-function init_cairomakie!()
-    CairoMakie.activate!()
-    set_theme!()
-
-    return nothing
-end
-
-function gen_figure()
-    text_width_in_pt = 405
-
-    return Figure(;
-        resolution=(text_width_in_pt, 0.55 * text_width_in_pt),
-        fontsize=10,
-        figure_padding=1,
-    )
-end
-
-function set_gaps!(fig::Figure)
-    colgap!(fig.layout, 5)
-    rowgap!(fig.layout, 5)
-
-    return nothing
-end
-
-function save_fig(filename::String, fig::Figure)
-    parent_dir = "exports/graphics"
-    mkpath(parent_dir)
-
-    return save("$parent_dir/$filename", fig; pt_per_unit=1)
-end
-
 function plot_mean_sq_displacement_with_expectation(
-    ts::Vector{Float64},
-    mean_sq_displacements::Matrix{Float64},
-    v₀s::AbstractVector{Float64},
-)
+    ts::Vector{Float64}, mean_sq_displacements::Matrix{Float64}, v₀s::NTuple{N,Float64}
+) where {N}
     init_cairomakie!()
 
     fig = gen_figure()
@@ -205,7 +162,7 @@ function plot_mean_sq_displacement_with_expectation(
 end
 
 function run_msd_analysis()
-    v₀s = SVector(0.0, 20.0, 40.0, 60.0, 80.0)
+    v₀s = (0.0, 20.0, 40.0, 60.0, 80.0)
 
     ts, mean_sq_displacements = mean_squared_displacement(;
         n_simulations=200 * Threads.nthreads(), v₀s=v₀s, T=100.0
@@ -216,17 +173,17 @@ function run_msd_analysis()
     return nothing
 end
 
-function plot_random_walk(T::Float64, v₀::Float64, seed::Int64)
+function plot_random_walk(; T::Float64, v₀::Float64, seed::Int64)
     Random.seed!(seed)
 
     half_box_len = max_possible_displacement(T, v₀)
     dir = msd_simulation(v₀, half_box_len, T, 8.0, "random_walk_$(Dates.now())")
 
     ts = Float64[]
-    fill_with_bundle_property!(ts, :t, dir)
-
     cs = SVector{2,Float64}[]
-    fill_with_bundle_property!(cs, :c, dir)
+    ReCo.append_bundle_properties!(
+        (ts, cs), (:t, :c), dir; particle_slice=1, snapshot_slice=:
+    )
 
     init_cairomakie!()
 
@@ -281,7 +238,7 @@ function plot_random_walk(T::Float64, v₀::Float64, seed::Int64)
 end
 
 function run_random_walk()
-    plot_random_walk(100_000.0, 0.0, 12345)
+    plot_random_walk(; T=100_000.0, v₀=0.0, seed=12345)
 
     return nothing
 end

analysis/pair_correlation.jl

@@ -1,28 +0,0 @@
-if splitdir(pwd())[2] == "analysis"
-    cd("..")
-end
-
-if splitdir(pwd())[2] != "ReCo.jl"
-    error("You have to be in the main directeory ReCo.jl!")
-else
-    include("src/analysis/pair_correlation_function.jl")
-end
-
-##
-
-using JLD2: JLD2
-
-using CairoMakie
-CairoMakie.activate!()
-set_theme!(theme_black())
-
-##
-
-dirs = readdir("exports/2021_11_19"; join=true)
-
-##
-
-for dir in dirs
-    data = JLD2.load("$dir/data.jld2")
-    display(pair_correlation(data["sol"], data["variables"]).fig)
-end

analysis/radial_distribution_function/g_of_r_d_ratio_with_0_45_packing_ratio.csv

@@ -0,0 +1,149 @@
+# Source: https://www.researchgate.net/publication/223333736_Reconstruction_of_Ross'_linear-mixing_model_for_the_equation_of_state_of_deuterium
+r_d_ratio,g_of_r_d_ratio
+0.99392,0.08604
+0.99388,0.18421
+0.99384,0.28847
+0.99382,0.34463
+0.99378,0.42994
+0.99372,0.52947
+0.99365,0.63441
+0.99361,0.68891
+0.99354,0.77406
+0.99344,0.87468
+0.99332,0.98013
+0.99325,1.03317
+0.99314,1.11822
+0.99297,1.21964
+0.99278,1.32539
+0.99267,1.37722
+0.99249,1.4622
+0.99224,1.56412
+0.99194,1.67
+0.99179,1.72089
+0.99152,1.8058
+0.99116,1.90793
+0.99074,2.01372
+0.99053,2.06399
+0.99015,2.1488
+0.98965,2.25086
+0.98909,2.35636
+0.98881,2.40637
+0.98831,2.49108
+0.98767,2.5925
+0.98699,2.69664
+0.98665,2.74669
+0.98607,2.83091
+0.98538,2.92885
+0.98467,3.02834
+0.98395,3.12761
+0.98362,3.17414
+0.98304,3.25677
+0.98242,3.34772
+0.98183,3.43776
+0.98129,3.52581
+0.98081,3.61118
+0.98043,3.69358
+0.98015,3.77322
+0.98,3.85156
+0.98014,3.95169
+0.98089,4.03823
+0.98357,4.10894
+0.9944,4.032
+1.00014,3.93792
+1.00441,3.85868
+1.00862,3.77571
+1.01272,3.69218
+1.01669,3.61072
+1.0216,3.51176
+1.02638,3.42062
+1.03128,3.33517
+1.03664,3.25085
+1.04281,3.16285
+1.04804,3.09295
+1.05614,2.98982
+1.06293,2.90588
+1.07039,2.81547
+1.07549,2.75441
+1.08352,2.65945
+1.09202,2.56047
+1.09758,2.49642
+1.10625,2.39796
+1.11214,2.33189
+1.12071,2.23742
+1.12939,2.14379
+1.13793,2.05435
+1.14612,1.97171
+1.15388,1.89756
+1.16359,1.81176
+1.17456,1.72585
+1.18615,1.64833
+1.19996,1.56268
+1.21169,1.47967
+1.22345,1.39143
+1.23543,1.30583
+1.24773,1.23043
+1.26375,1.14911
+1.28235,1.06205
+1.30026,0.97628
+1.32001,0.89731
+1.34607,0.8185
+1.38029,0.74139
+1.42373,0.67706
+1.47928,0.62965
+1.54985,0.64005
+1.60887,0.68208
+1.65956,0.742
+1.70693,0.80174
+1.75549,0.86677
+1.7997,0.93281
+1.8408,0.99879
+1.88278,1.07151
+1.92418,1.14479
+1.96069,1.20961
+2.00089,1.27998
+2.06624,1.29427
+2.12188,1.24296
+2.16454,1.1732
+2.20365,1.10578
+2.2468,1.03273
+2.28968,0.97094
+2.34711,0.92061
+2.40659,0.87706
+2.47617,0.86598
+2.54506,0.87775
+2.6091,0.91495
+2.67198,0.9548
+2.73415,0.99707
+2.79524,1.03776
+2.86004,1.06782
+2.92949,1.08716
+2.99726,1.09694
+3.06572,1.07401
+3.13305,1.0512
+3.20406,1.01975
+3.26722,0.99225
+3.33163,0.96833
+3.40415,0.95312
+3.46976,0.95561
+3.53944,0.96948
+3.61034,0.98695
+3.68191,1.00826
+3.74431,1.02656
+3.81351,1.04126
+3.88832,1.04634
+3.95948,1.04227
+4.02721,1.03357
+4.09307,1.02294
+4.16757,1.01022
+4.23256,0.99784
+4.30367,0.98757
+4.3786,0.98803
+4.44527,0.99118
+4.51289,0.99726
+4.58223,1.00774
+4.65288,1.01261
+4.72314,1.01386
+4.79174,1.02311
+4.86131,1.03134
+4.93191,1.0253
+4.98585,1.01943

analysis/radial_distribution_function/radial_distribution_function.jl

@@ -0,0 +1,216 @@
+using CairoMakie
+using LaTeXStrings: @L_str
+using StaticArrays: SVector
+using JLD2: JLD2
+using Dates: Dates
+using CSV: CSV
+using DataFrames: DataFrames
+using Random: Random
+
+using ReCo: ReCo
+
+includet("../../graphics/common_CairoMakie.jl")
+
+function radial_distribution_simulation(;
+    n_particles::Int64, v₀s::NTuple{N,Float64}, T::Float64, packing_ratio::Float64
+) where {N}
+    Random.seed!(42)
+
+    n_v₀s = length(v₀s)
+
+    sim_dirs = Vector{String}(undef, n_v₀s)
+
+    parent_dir = "radial_distribution_$(Dates.now())"
+
+    Threads.@threads for v₀_ind in 1:n_v₀s
+        v₀ = v₀s[v₀_ind]
+
+        dir = ReCo.init_sim(;
+            n_particles=n_particles,
+            v₀=v₀,
+            packing_ratio=packing_ratio,
+            parent_dir=parent_dir,
+            comment="$v₀",
+        )
+
+        ReCo.run_sim(dir; duration=T, seed=v₀_ind)
+
+        sim_dirs[v₀_ind] = dir
+    end
+
+    return sim_dirs
+end
+
+function circular_shell_volume(lower_radius, Δradius)
+    return π * (2 * lower_radius * Δradius + Δradius^2)
+end
+
+function radial_distribution(;
+    sim_dirs::Vector{String}, n_radii::Int64, n_last_snapshots::Int64, n_particles::Int64
+)
+    sim_consts::ReCo.SimConsts = JLD2.load_object("$(sim_dirs[1])/sim_consts.jld2")
+
+    particle_radius = sim_consts.particle_radius
+
+    min_lower_radius = 0.0
+    max_lower_radius = 5 * (2 * particle_radius)
+    Δradius = (max_lower_radius - min_lower_radius) / n_radii
+
+    lower_radii = LinRange(min_lower_radius, max_lower_radius, n_radii)
+
+    box_volume = (2 * sim_consts.half_box_len)^2
+    volume_per_particle = box_volume / n_particles
+
+    n_sim_dirs = length(sim_dirs)
+
+    gs = Vector{Vector{Float64}}(undef, n_sim_dirs)
+
+    Threads.@threads for sim_dir_ind in 1:n_sim_dirs
+        sim_dir = sim_dirs[sim_dir_ind]
+
+        cs = Matrix{SVector{2,Float64}}(undef, (n_particles, n_last_snapshots))
+
+        bundle_paths = ReCo.sorted_bundle_paths(sim_dir; rev=true)
+
+        snapshot_conunter = 0
+        break_bundle_path_loop = false
+
+        for bundle_path in bundle_paths
+            bundle::ReCo.Bundle = JLD2.load_object(bundle_path)
+
+            for snapshot_ind in (bundle.n_snapshots):-1:1
+                snapshot_conunter += 1
+
+                @simd for particle_ind in 1:n_particles
+                    cs[particle_ind, snapshot_conunter] = bundle.c[
+                        particle_ind, snapshot_ind
+                    ]
+                end
+
+                if snapshot_conunter == n_last_snapshots
+                    break_bundle_path_loop = true
+                    break
+                end
+            end
+
+            if break_bundle_path_loop
+                break
+            end
+        end
+
+        if snapshot_conunter != n_last_snapshots
+            error("snapshot_conunter != n_last_snapshots")
+        end
+
+        g = zeros(n_radii)
+
+        for snapshot_ind in 1:n_last_snapshots
+            for p1_ind in 1:n_particles
+                for p2_ind in (p1_ind + 1):n_particles
+                    c1 = cs[p1_ind, snapshot_ind]
+                    c2 = cs[p2_ind, snapshot_ind]
+
+                    r⃗₁₂ = c2 - c1
+                    r⃗₁₂ = ReCo.minimum_image(r⃗₁₂, sim_consts.half_box_len)
+
+                    distance = ReCo.norm2d(r⃗₁₂)
+
+                    lower_radius_ind = ceil(Int64, distance / Δradius)
+
+                    if lower_radius_ind <= n_radii
+                        g[lower_radius_ind] += 2
+                    end
+                end
+            end
+        end
+
+        for (lower_radius_ind, lower_radius) in enumerate(lower_radii)
+            g[lower_radius_ind] *=
+                volume_per_particle / (
+                    n_last_snapshots *
+                    n_particles *
+                    circular_shell_volume(lower_radius, Δradius)
+                )
+        end
+
+        gs[sim_dir_ind] = g
+    end
+
+    return (lower_radii, gs, particle_radius)
+end
+
+function plot_radial_distributions(
+    v₀s::NTuple{N,Float64},
+    lower_radii::AbstractVector{Float64},
+    gs::Vector{Vector{Float64}},
+    particle_radius::Float64,
+) where {N}
+    println("Plotting the radial distributions")
+
+    init_cairomakie!()
+
+    fig = gen_figure()
+
+    max_lower_radius = maximum(lower_radii)
+    max_g = maximum(maximum.(gs))
+
+    ax = Axis(
+        fig[1:2, 1:2];
+        xticks=0:(2 * particle_radius):ceil(Int64, max_lower_radius),
+        yticks=0:ceil(Int64, max_g),
+        xlabel=L"r / d",
+        ylabel=L"g",
+    )
+
+    lines!(
+        ax,
+        SVector(0.0, max_lower_radius),
+        SVector(1.0, 1.0);
+        linestyle=:dash,
+        color=:red,
+        linewidth=1,
+    )
+
+    expected = CSV.read(
+        "analysis/g_of_r_d_ratio_with_0_45_packing_ratio.csv",
+        DataFrames.DataFrame;
+        header=2,
+    )
+    lines!(
+        ax,
+        expected.r_d_ratio,
+        expected.g_of_r_d_ratio;
+        label="Expectation for dense hard spheres",
+        color=:green,
+    )
+
+    for (g_ind, g) in enumerate(gs)
+        scatterlines!(
+            ax, lower_radii, g; markersize=3, linewidth=0.5, label=L"v_0 = %$(v₀s[g_ind])"
+        )
+    end
+
+    axislegend(ax; position=:rt, padding=3, rowgap=-3)
+
+    save_fig("radial_distribution.pdf", fig)
+
+    return nothing
+end
+
+function run_radial_distribution_analysis()
+    v₀s = (0.0, 80.0)
+
+    n_particles = 1000
+
+    sim_dirs = radial_distribution_simulation(;
+        n_particles=n_particles, v₀s=v₀s, T=100.0, packing_ratio=0.45
+    )
+
+    lower_radii, gs, particle_radius = radial_distribution(;
+        sim_dirs=sim_dirs, n_radii=75, n_last_snapshots=200, n_particles=n_particles
+    )
+
+    plot_radial_distributions(v₀s, lower_radii, gs, particle_radius)
+
+    return nothing
+end