From eb72bf8be7a7a21f4c5567ef2299b8816151ede7 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Mon, 11 May 2026 17:44:53 +0200 Subject: [PATCH 01/24] feat: scaffold msad.jl with header comment --- src/ParticlesMC.jl | 1 + src/msad.jl | 10 ++++++++++ 2 files changed, 11 insertions(+) create mode 100644 src/msad.jl diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl index f652cb0..b33197e 100644 --- a/src/ParticlesMC.jl +++ b/src/ParticlesMC.jl @@ -20,6 +20,7 @@ include("models.jl") include("molecules.jl") include("atoms.jl") include("moves.jl") +include("msad.jl") """Return the position of particle `i` in `system`. diff --git a/src/msad.jl b/src/msad.jl new file mode 100644 index 0000000..b4039f1 --- /dev/null +++ b/src/msad.jl @@ -0,0 +1,10 @@ +# msad.jl +# +# On-the-fly Mean Squared Angular Displacement (MSAD) tracker. +# Implements three methods: +# - Integral : accumulates incremental rotation vectors every step +# - Threshold : relays the reference frame when rotation exceeds θ_T +# - Euler : compares directly to initial frame, computed on log schedule + + + From 19fd098de85d733446e0f9d37b67fa486215b709 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 12 May 2026 10:35:56 +0200 Subject: [PATCH 02/24] add MSADState and MSADTracker structs --- src/msad.jl | 77 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/src/msad.jl b/src/msad.jl index b4039f1..dbd46a3 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -6,5 +6,82 @@ # - Threshold : relays the reference frame when rotation exceeds θ_T # - Euler : compares directly to initial frame, computed on log schedule +using LinearAlgebra +using StaticArrays +##### Basic operations for rotation inspired fromo our python post processing ##### +########## ########## +function body_frame(r1::SVector{3,T}, r2::SVector{3,T}, r3::SVector{3,T}, L::T) where {T} + e1 = r2 - r1 + v = r3 - r1 + # minimum image convention + e1 = e1 - round.(e1 / L) .* L + v = v - round.( v / L) .* L + # Gram-Schmidt basis right handed + e1 = e1 / norm(e1) + e2 = cross(e1, v) + e2 = e2 / norm(e2) + e3 = cross(e1, e2) + # hcat stacks column vectors → 3×3 matrix + return hcat(e1, e2, e3) # SMatrix{3,3,T} +end + +function get_all_body_frames(system::Molecules) + L = system.box[1] # cubic box + T = typeof(L) + R = Vector{SMatrix{3,3,T,9}}(undef, system.Nmol) # pre-allocate + pos = system.position + @inbounds for m in 1:system.Nmol # iterate over the number of molecules + s = system.start_mol[m] # first bead index for the considered molecule + R[m] = body_frame(pos[s], pos[s+1], pos[s+2], L) + end + return R +end + +function rotation_vector(R::SMatrix{3,3,T}) where {T} + cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) + # skew-symmetric part encodes sin(θ)*n + skew = (R - R') / 2 + sin_n = SVector(skew[3,2], skew[1,3], skew[2,1]) # julia starts index at 1 + sin_θ = norm(sin_n) + θ = atan(sin_θ, cos_θ) + return sin_θ > 1e-10 ? (θ / sin_θ) * sin_n : zero(SVector{3,T}) ## might need to add the theta = pi case (unprobable for MC steps) +end + +##### Definition of the MSAD State ie: the accumulating variables ##### +# mutable because we need to update it +########## ########## + +mutable struct MSADState{T} + R0::Vector{SMatrix{3,3,T,9}} # (Nmol,3,3) initial body frames at t=0 + R_prev::Vector{SMatrix{3,3,T,9}} # (Nmol,3,3) body frames at previous time + phi_integral::Vector{SVector{3,T}} # (Nmol,3) Accumulated vector for integral method + R_ref_thresh::Vector{SMatrix{3,3,T,9}} # (Nmol,3,3) reference body frames for each molecule in the threshold method + phi_acc::Vector{SVector{3,T}} # (Nmol,3) Accumulated vector for threshold method + initialized::Bool +end + +# Start empty, filled during initialise when we first see the system +MSADState{T}() where {T} = MSADState{T}([], [], [], [], [], false) + +##### The MSAD computation algorithm ##### +########## ########## + +mutable struct MSADTracker{T} <: AriannaAlgorithm + states::Vector{MSADState{T}} # the states of the considered chain + theta_T::T # threshold angle in radians + output_schedule::Vector{Int} # schedule for writing MSAD to disk + paths::Vector{String} + files::Vector{IOStream} +end + +function MSADTracker(chains, theta_T::T, output_schedule::Vector{Int}, path::String) where {T} + n = length(chains) + states = [MSADState{T}() for _ in 1:n] # one empty state per chain filled during initialise + + dirs = [joinpath(path, "chains", "$c") for c in 1:n] + paths = [joinpath(d, "msad.dat") for d in dirs] + files = Vector{IOStream}(undef, n) + return MSADTracker{T}(states, theta_T, output_schedule, paths, files) +end \ No newline at end of file From 36838de2bcf345259e6136f42e3e5b052a4cf5d6 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 12 May 2026 10:44:45 +0200 Subject: [PATCH 03/24] add initialise and finalise for MSADTracker --- src/msad.jl | 49 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/src/msad.jl b/src/msad.jl index dbd46a3..1718e6d 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -79,9 +79,56 @@ end function MSADTracker(chains, theta_T::T, output_schedule::Vector{Int}, path::String) where {T} n = length(chains) states = [MSADState{T}() for _ in 1:n] # one empty state per chain filled during initialise - + dirs = [joinpath(path, "chains", "$c") for c in 1:n] paths = [joinpath(d, "msad.dat") for d in dirs] files = Vector{IOStream}(undef, n) return MSADTracker{T}(states, theta_T, output_schedule, paths, files) +end + +##### Initialisation of the simulation ##### +# called once at t = 0, set the system, the mutable struct and the storing paths +########## ########## + +function Arianna.initialise(tracker::MSADTracker, simulation::Simulation) + + for c in eachindex(simulation.chains) + system = simulation.chains[c] + state = tracker.states[c] + T = typeof(system.temperature) # get the float type from system so one can choose Float64 or 32 (generic) + + # create output directory if it doesn't exist + mkpath(dirname(tracker.paths[c])) + + # open output file — write header + tracker.files[c] = open(tracker.paths[c], "w") + println(tracker.files[c], "# t msad_euler msad_integral msad_thresh") + + # compute initial body frames + R_all = get_all_body_frames(system) + N_mol = system.Nmol + + # fill state, copy. so they are different pointers! + state.R0 = copy(R_all) + state.R_prev = copy(R_all) + state.phi_integral = [zero(SVector{3,T}) for _ in 1:N_mol] + state.R_ref_thresh = copy(R_all) + state.phi_acc = [zero(SVector{3,T}) for _ in 1:N_mol] + state.initialized = true + + # write t=0, all MSAD values are zero at start + println(tracker.files[c], "0 0.0 0.0 0.0") + flush(tracker.files[c]) + end + +end + +##### Finalise th simulation ##### +# close output to not keep unwritten data in memory +########## ########## + +function Arianna.finalise(tracker::MSADTracker, ::Simulation) + for f in tracker.files + close(f) + end end \ No newline at end of file From c19c9379fdae9a8453a943f4bfef5341f76a583e Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 12 May 2026 10:52:47 +0200 Subject: [PATCH 04/24] add compute_schedule to MSADTracker with subset assertion --- src/msad.jl | 85 +++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 83 insertions(+), 2 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index 1718e6d..3e1006a 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -71,19 +71,23 @@ MSADState{T}() where {T} = MSADState{T}([], [], [], [], [], false) mutable struct MSADTracker{T} <: AriannaAlgorithm states::Vector{MSADState{T}} # the states of the considered chain theta_T::T # threshold angle in radians + compute_schedule::Vector{Int} # schedule for computing rotation for integral and threshold methods output_schedule::Vector{Int} # schedule for writing MSAD to disk paths::Vector{String} files::Vector{IOStream} end -function MSADTracker(chains, theta_T::T, output_schedule::Vector{Int}, path::String) where {T} +function MSADTracker(chains, theta_T::T, compute_schedule::Vector{Int}, output_schedule::Vector{Int}, path::String) where {T} + + @assert all(t in compute_schedule for t in output_schedule) # output_schedule must be a subset of compute_schedule + n = length(chains) states = [MSADState{T}() for _ in 1:n] # one empty state per chain filled during initialise dirs = [joinpath(path, "chains", "$c") for c in 1:n] paths = [joinpath(d, "msad.dat") for d in dirs] files = Vector{IOStream}(undef, n) - return MSADTracker{T}(states, theta_T, output_schedule, paths, files) + return MSADTracker{T}(states, theta_T, compute_schedule, output_schedule, paths, files) end ##### Initialisation of the simulation ##### @@ -131,4 +135,81 @@ function Arianna.finalise(tracker::MSADTracker, ::Simulation) for f in tracker.files close(f) end +end + +##### Make a step in simulation ##### +########## ########## + +function Arianna.make_step!(simulation::Simulation, tracker::MSADTracker) + t = simulation.t + + for c in eachindex(simulation.chains) + system = simulation.chains[c] + state = tracker.states[c] + N_mol = system.Nmol + + # compute current body frames shared by all three methods + R_all = get_all_body_frames(system) + + # ── Integral update (every step) ────────────────────────── + # Python: dR = np.einsum('mij,mik->mjk', R_prev_all, R_all) + # phi_all += rotation_vector_batch(dR) + for m in 1:N_mol + dR = state.R_prev[m]' * R_all[m] + state.phi_integral[m] = state.phi_integral[m] + rotation_vector(dR) + end + + # ── Threshold update (every step) ───────────────────────── + # Python: dR_thresh = np.einsum('mij,mik->mjk', R_ref_thresh[k], R_all) + # phi_current = rotation_vector_batch(dR_thresh) + # phi_total = phi_acc_thresh[k] + phi_current + # relay = norms >= theta_T + phi_current = Vector{SVector{3,eltype(tracker.theta_T)}}(undef, N_mol) + phi_total = Vector{SVector{3,eltype(tracker.theta_T)}}(undef, N_mol) + for m in 1:N_mol + dR = state.R_ref_thresh[m]' * R_all[m] + phi_current[m] = rotation_vector(dR) + phi_total[m] = state.phi_acc[m] + phi_current[m] + # relay — Python: phi_acc_thresh[k][relay] = phi_total[relay] + # R_ref_thresh[k][relay] = R_all[relay] + if norm(phi_current[m]) >= tracker.theta_T + state.phi_acc[m] = phi_total[m] + state.R_ref_thresh[m] = R_all[m] + end + end + + # ── Advance R_prev (every step) ─────────────────────────── + # Python: R_prev_all = R_all.copy() + state.R_prev .= R_all + + # ── Output (only on output_schedule) ────────────────────── + # Python: times.append(frame_idx) + msad_euler/integral/thresh.append(...) + if t in tracker.output_schedule + # Euler — only computed here, not accumulated + # Python: R_rel = np.einsum('mij,mik->mjk', R0_all, R_all) + # thetas = rotation_vector_batch(R_rel) + # msad_e = np.sum(norm(thetas)**2) / N_mol + msad_euler = sum( + norm(rotation_vector(state.R0[m]' * R_all[m]))^2 + for m in 1:N_mol + ) / N_mol + + # Integral + # Python: msad_i = np.sum(norm(phi_all)**2) / N_mol + msad_integral = sum( + norm(state.phi_integral[m])^2 + for m in 1:N_mol + ) / N_mol + + # Threshold + # Python: np.sum(norm(phi_total)**2) / N_mol + msad_thresh = sum( + norm(phi_total[m])^2 + for m in 1:N_mol + ) / N_mol + + println(tracker.files[c], "$t $msad_euler $msad_integral $msad_thresh") + flush(tracker.files[c]) + end + end end \ No newline at end of file From 6bb1a4ff754aa3f44a138d6d55c43ffbcb64d18c Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 12 May 2026 10:58:31 +0200 Subject: [PATCH 05/24] add MSADTracker construction and assertion tests --- test/test_tracker.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 test/test_tracker.jl diff --git a/test/test_tracker.jl b/test/test_tracker.jl new file mode 100644 index 0000000..ca962c0 --- /dev/null +++ b/test/test_tracker.jl @@ -0,0 +1,21 @@ +using ParticlesMC, StaticArrays, LinearAlgebra + +# check MSADState constructs empty +s = ParticlesMC.MSADState{Float64}() +println("MSADState empty: ", s.initialized == false) + +# check MSADTracker constructor with valid schedules +compute = collect(0:10:100) +output = collect(0:20:100) +chains = [1, 2] +t = ParticlesMC.MSADTracker(chains, pi/4, compute, output, "/tmp/test_msad") +println("MSADTracker created: ", t.theta_T ≈ pi/4) + +# check assertion fires on invalid schedules +try + bad_output = [0, 5, 10] # 5 not in compute schedule + ParticlesMC.MSADTracker(chains, pi/4, compute, bad_output, "/tmp/test_msad") + println("ERROR: should have thrown") +catch e + println("Assertion caught correctly: ", e isa AssertionError) +end \ No newline at end of file From 42be0d7bf15815c4d8087b0c6ad04998cfd40fbd Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 12 May 2026 11:22:47 +0200 Subject: [PATCH 06/24] add the MSADTracker to the particlesmc loop --- src/ParticlesMC.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl index b33197e..4edd462 100644 --- a/src/ParticlesMC.jl +++ b/src/ParticlesMC.jl @@ -242,6 +242,31 @@ ParticlesMC implemented in Comonicon. else error("Unsupported policy: $policy for action: $action") end + elseif alg == "MSADTracker" + # read theta_T from parameters + parameters = get(output, "parameters", Dict()) + theta_T = Float64(get(parameters, "theta_T", π/4)) + + # build output schedule separately from compute schedule + output_scheduler_params = get(output, "output_scheduler_params", Dict()) + if "log_base" in keys(output_scheduler_params) + output_sched = build_schedule(steps, burn, 2.0) + elseif "linear_interval" in keys(output_scheduler_params) + output_sched = build_schedule(steps, burn, + output_scheduler_params["linear_interval"]) + else + output_sched = sched # default: same as compute + end + + # convert to Vector{Int} — build_schedule returns this already + algorithm = ( + algorithm=MSADTracker, + scheduler=sched, # compute schedule + theta_T=theta_T, + output_schedule=output_sched, + path=output_path, + ) + else error("Unsupported action: $action") end From c3e1e24711bbe60c500c9827af80e25d3aa5f894 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 12 May 2026 11:25:34 +0200 Subject: [PATCH 07/24] test operations in msad.jl --- test/test_msad.jl | 157 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 test/test_msad.jl diff --git a/test/test_msad.jl b/test/test_msad.jl new file mode 100644 index 0000000..0b12c2b --- /dev/null +++ b/test/test_msad.jl @@ -0,0 +1,157 @@ +using ParticlesMC +using StaticArrays +using LinearAlgebra +using Test + +@testset "body_frame" begin + + # known inputs — worked out by hand + r1 = SVector(0.0, 0.0, 0.0) + r2 = SVector(1.0, 0.0, 0.0) + r3 = SVector(0.0, 1.0, 0.0) + L = 10.0 + + R = ParticlesMC.body_frame(r1, r2, r3, L) + + # expected result from hand calculation: + # e1 = [1, 0, 0] + # e2 = [0, 0, 1] + # e3 = [0, -1, 0] + @test R[:, 1] ≈ SVector( 1.0, 0.0, 0.0) + @test R[:, 2] ≈ SVector( 0.0, 0.0, 1.0) + @test R[:, 3] ≈ SVector( 0.0, -1.0, 0.0) + + # R must be a valid rotation matrix: + # columns orthonormal → R'R = I + # right-handed → det(R) = +1 + @test R' * R ≈ I + @test det(R) ≈ 1.0 + +end + +@testset "body_frame minimum image" begin + + # same molecule but r2 is wrapped across the boundary + # r1 at 0.1, r2 at 9.9 — bond should be -0.2, not 9.8 + r1 = SVector(0.1, 0.0, 0.0) + r2 = SVector(9.9, 0.0, 0.0) + r3 = SVector(0.1, 1.0, 0.0) + L = 10.0 + + R_wrapped = ParticlesMC.body_frame(r1, r2, r3, L) + + # unwrapped version — same molecule, same frame + r2_unwrapped = SVector(-0.1, 0.0, 0.0) + R_unwrapped = ParticlesMC.body_frame(r1, r2_unwrapped, r3, L) + + @test R_wrapped ≈ R_unwrapped + +end + +@testset "rotation_vector" begin + + # identity matrix → zero rotation + R_id = SMatrix{3,3}(1.0I) + @test ParticlesMC.rotation_vector(R_id) ≈ SVector(0.0, 0.0, 0.0) + + # 90° rotation around z-axis + # R = [0 -1 0] + # [1 0 0] + # [0 0 1] + R_90z = SMatrix{3,3}(0.0, 1.0, 0.0, + -1.0, 0.0, 0.0, + 0.0, 0.0, 1.0) + Ω = ParticlesMC.rotation_vector(R_90z) + # should be θ=π/2 around z → Ω = [0, 0, π/2] + @test norm(Ω) ≈ π/2 + @test Ω / norm(Ω) ≈ SVector(0.0, 0.0, 1.0) + +end + +@testset "body_frame PBC - molecule split across each axis" begin + + L = 10.0 + + # --- x axis --- + # r1 near left wall, r2 wrapped to right wall + # true bond is [-0.2, 0, 0], not [9.8, 0, 0] + r1_x = SVector(0.1, 0.0, 0.0) + r2_x = SVector(9.9, 0.0, 0.0) # wrapped + r3_x = SVector(0.1, 1.0, 0.0) + R_wrapped_x = ParticlesMC.body_frame(r1_x, r2_x, r3_x, L) + R_unwrapped_x = ParticlesMC.body_frame(r1_x, SVector(-0.1, 0.0, 0.0), r3_x, L) + @test R_wrapped_x ≈ R_unwrapped_x + + # --- y axis --- + r1_y = SVector(0.0, 0.1, 0.0) + r2_y = SVector(1.0, 0.1, 0.0) + r3_y = SVector(0.0, 9.9, 0.0) # r3 wrapped across y + R_wrapped_y = ParticlesMC.body_frame(r1_y, r2_y, r3_y, L) + R_unwrapped_y = ParticlesMC.body_frame(r1_y, r2_y, SVector(0.0, -0.1, 0.0), L) + @test R_wrapped_y ≈ R_unwrapped_y + + # --- z axis --- + r1_z = SVector(0.0, 0.0, 0.1) + r2_z = SVector(1.0, 0.0, 0.1) + r3_z = SVector(0.0, 1.0, 0.1) + r2_z_wrapped = SVector(1.0, 0.0, 9.9 + 0.1) # wrapped across z... wait + # r1 near top wall, r2 wrapped to bottom + r1_z2 = SVector(0.0, 0.0, 0.1) + r2_z2 = SVector(0.0, 0.0, 9.9) # wrapped + r3_z2 = SVector(1.0, 0.0, 0.1) + R_wrapped_z = ParticlesMC.body_frame(r1_z2, r2_z2, r3_z2, L) + R_unwrapped_z = ParticlesMC.body_frame(r1_z2, SVector(0.0, 0.0, -0.1), r3_z2, L) + @test R_wrapped_z ≈ R_unwrapped_z + +end + +@testset "body_frame PBC - result is always a valid rotation matrix" begin + + L = 10.0 + # test many random molecules, some crossing boundaries + # a valid rotation matrix always satisfies R'R = I and det(R) = +1 + import Random + Random.seed!(42) + for _ in 1:100 + # random positions anywhere in box + r1 = SVector(rand()*L, rand()*L, rand()*L) + r2 = SVector(rand()*L, rand()*L, rand()*L) + r3 = SVector(rand()*L, rand()*L, rand()*L) + R = ParticlesMC.body_frame(r1, r2, r3, L) + @test R' * R ≈ I atol=1e-10 + @test det(R) ≈ 1.0 atol=1e-10 + end + +end + +@testset "rotation_vector PBC - roundtrip" begin + + # if we compute R between two frames, rotation_vector should + # give back an Ω whose norm equals the angle between them + # and applying it back should reconstruct the rotation + + # 180° rotation around x axis + R_180x = SMatrix{3,3}(1.0, 0.0, 0.0, + 0.0, -1.0, 0.0, + 0.0, 0.0, -1.0) + Ω = ParticlesMC.rotation_vector(R_180x) + @test norm(Ω) ≈ π atol=1e-10 + @test Ω / norm(Ω) ≈ SVector(1.0, 0.0, 0.0) atol=1e-10 + + # 180° rotation around y axis + R_180y = SMatrix{3,3}(-1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0,-1.0) + Ω = ParticlesMC.rotation_vector(R_180y) + @test norm(Ω) ≈ π atol=1e-10 + @test Ω / norm(Ω) ≈ SVector(0.0, 1.0, 0.0) atol=1e-10 + + # small rotation — near-identity, tests the sin≈0 branch + ε = 1e-8 + R_tiny = SMatrix{3,3}(1.0, ε, 0.0, + -ε, 1.0, 0.0, + 0.0, 0.0, 1.0) + Ω = ParticlesMC.rotation_vector(R_tiny) + @test norm(Ω) < 1e-6 # nearly zero rotation + +end \ No newline at end of file From 6163750c8c267390cbb8e2c745a976e9a49a7e0e Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 12 May 2026 14:27:16 +0200 Subject: [PATCH 08/24] cleaning msad.jl --- src/msad.jl | 24 ++++-------------------- 1 file changed, 4 insertions(+), 20 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index 3e1006a..92340dd 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -151,58 +151,42 @@ function Arianna.make_step!(simulation::Simulation, tracker::MSADTracker) # compute current body frames shared by all three methods R_all = get_all_body_frames(system) - # ── Integral update (every step) ────────────────────────── - # Python: dR = np.einsum('mij,mik->mjk', R_prev_all, R_all) - # phi_all += rotation_vector_batch(dR) + # Integral update for m in 1:N_mol dR = state.R_prev[m]' * R_all[m] state.phi_integral[m] = state.phi_integral[m] + rotation_vector(dR) end - # ── Threshold update (every step) ───────────────────────── - # Python: dR_thresh = np.einsum('mij,mik->mjk', R_ref_thresh[k], R_all) - # phi_current = rotation_vector_batch(dR_thresh) - # phi_total = phi_acc_thresh[k] + phi_current - # relay = norms >= theta_T + # Threshold update phi_current = Vector{SVector{3,eltype(tracker.theta_T)}}(undef, N_mol) phi_total = Vector{SVector{3,eltype(tracker.theta_T)}}(undef, N_mol) for m in 1:N_mol dR = state.R_ref_thresh[m]' * R_all[m] phi_current[m] = rotation_vector(dR) phi_total[m] = state.phi_acc[m] + phi_current[m] - # relay — Python: phi_acc_thresh[k][relay] = phi_total[relay] - # R_ref_thresh[k][relay] = R_all[relay] if norm(phi_current[m]) >= tracker.theta_T state.phi_acc[m] = phi_total[m] state.R_ref_thresh[m] = R_all[m] end end - # ── Advance R_prev (every step) ─────────────────────────── - # Python: R_prev_all = R_all.copy() state.R_prev .= R_all - # ── Output (only on output_schedule) ────────────────────── - # Python: times.append(frame_idx) + msad_euler/integral/thresh.append(...) + # output is time match output schedule if t in tracker.output_schedule - # Euler — only computed here, not accumulated - # Python: R_rel = np.einsum('mij,mik->mjk', R0_all, R_all) - # thetas = rotation_vector_batch(R_rel) - # msad_e = np.sum(norm(thetas)**2) / N_mol + # Euler only computed here not accumulated msad_euler = sum( norm(rotation_vector(state.R0[m]' * R_all[m]))^2 for m in 1:N_mol ) / N_mol # Integral - # Python: msad_i = np.sum(norm(phi_all)**2) / N_mol msad_integral = sum( norm(state.phi_integral[m])^2 for m in 1:N_mol ) / N_mol # Threshold - # Python: np.sum(norm(phi_total)**2) / N_mol msad_thresh = sum( norm(phi_total[m])^2 for m in 1:N_mol From f0e90fd89e787b95a9e5dfa1d96523697a9eff6f Mon Sep 17 00:00:00 2001 From: Cyprien Fraysse Date: Tue, 12 May 2026 11:29:08 +0200 Subject: [PATCH 09/24] feat: add MSADTracker branch to TOML parser --- src/ParticlesMC.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl index 4edd462..607204e 100644 --- a/src/ParticlesMC.jl +++ b/src/ParticlesMC.jl @@ -315,6 +315,31 @@ ParticlesMC implemented in Comonicon. algorithm=eval(Meta.parse(alg)), scheduler=sched, ) + elseif alg == "MSADTracker" + # read theta_T from parameters + parameters = get(output, "parameters", Dict()) + theta_T = Float64(get(parameters, "theta_T", π/4)) + + # build output schedule separately from compute schedule + output_scheduler_params = get(output, "output_scheduler_params", Dict()) + if "log_base" in keys(output_scheduler_params) + output_sched = build_schedule(steps, burn, 2.0) + elseif "linear_interval" in keys(output_scheduler_params) + output_sched = build_schedule(steps, burn, + output_scheduler_params["linear_interval"]) + else + output_sched = sched # default: same as compute + end + + # convert to Vector{Int} — build_schedule returns this already + algorithm = ( + algorithm=MSADTracker, + scheduler=sched, # compute schedule + theta_T=theta_T, + output_schedule=output_sched, + path=output_path, + ) + else error("Unsupported output algorithm: $alg") end From 3ecaa6101d0f1186a8445fd03783c93f62a64ee9 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 10:58:59 +0200 Subject: [PATCH 10/24] fix tracker --> algorithm, for consistency with Arianna.jl --- src/msad.jl | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index 92340dd..384ab3c 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -94,19 +94,19 @@ end # called once at t = 0, set the system, the mutable struct and the storing paths ########## ########## -function Arianna.initialise(tracker::MSADTracker, simulation::Simulation) +function Arianna.initialise(algorithm::MSADTracker, simulation::Simulation) for c in eachindex(simulation.chains) system = simulation.chains[c] - state = tracker.states[c] + state = algorithm.states[c] T = typeof(system.temperature) # get the float type from system so one can choose Float64 or 32 (generic) # create output directory if it doesn't exist - mkpath(dirname(tracker.paths[c])) + mkpath(dirname(algorithm.paths[c])) # open output file — write header - tracker.files[c] = open(tracker.paths[c], "w") - println(tracker.files[c], "# t msad_euler msad_integral msad_thresh") + algorithm.files[c] = open(algorithm.paths[c], "w") + println(algorithm.files[c], "# t msad_euler msad_integral msad_thresh") # compute initial body frames R_all = get_all_body_frames(system) @@ -121,8 +121,8 @@ function Arianna.initialise(tracker::MSADTracker, simulation::Simulation) state.initialized = true # write t=0, all MSAD values are zero at start - println(tracker.files[c], "0 0.0 0.0 0.0") - flush(tracker.files[c]) + println(algorithm.files[c], "0 0.0 0.0 0.0") + flush(algorithm.files[c]) end end @@ -131,8 +131,8 @@ end # close output to not keep unwritten data in memory ########## ########## -function Arianna.finalise(tracker::MSADTracker, ::Simulation) - for f in tracker.files +function Arianna.finalise(algorithm::MSADTracker, ::Simulation) + for f in algorithm.files close(f) end end @@ -140,12 +140,12 @@ end ##### Make a step in simulation ##### ########## ########## -function Arianna.make_step!(simulation::Simulation, tracker::MSADTracker) +function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) t = simulation.t for c in eachindex(simulation.chains) system = simulation.chains[c] - state = tracker.states[c] + state = algorithm.states[c] N_mol = system.Nmol # compute current body frames shared by all three methods @@ -158,13 +158,13 @@ function Arianna.make_step!(simulation::Simulation, tracker::MSADTracker) end # Threshold update - phi_current = Vector{SVector{3,eltype(tracker.theta_T)}}(undef, N_mol) - phi_total = Vector{SVector{3,eltype(tracker.theta_T)}}(undef, N_mol) + phi_current = Vector{SVector{3,eltype(algorithm.theta_T)}}(undef, N_mol) + phi_total = Vector{SVector{3,eltype(algorithm.theta_T)}}(undef, N_mol) for m in 1:N_mol dR = state.R_ref_thresh[m]' * R_all[m] phi_current[m] = rotation_vector(dR) phi_total[m] = state.phi_acc[m] + phi_current[m] - if norm(phi_current[m]) >= tracker.theta_T + if norm(phi_current[m]) >= algorithm.theta_T state.phi_acc[m] = phi_total[m] state.R_ref_thresh[m] = R_all[m] end @@ -173,7 +173,7 @@ function Arianna.make_step!(simulation::Simulation, tracker::MSADTracker) state.R_prev .= R_all # output is time match output schedule - if t in tracker.output_schedule + if t in algorithm.output_schedule # Euler only computed here not accumulated msad_euler = sum( norm(rotation_vector(state.R0[m]' * R_all[m]))^2 @@ -192,8 +192,8 @@ function Arianna.make_step!(simulation::Simulation, tracker::MSADTracker) for m in 1:N_mol ) / N_mol - println(tracker.files[c], "$t $msad_euler $msad_integral $msad_thresh") - flush(tracker.files[c]) + println(algorithm.files[c], "$t $msad_euler $msad_integral $msad_thresh") + flush(algorithm.files[c]) end end end \ No newline at end of file From 2d82dca7314fee69ae4d91b9e0fb54b8243dc6df Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 11:13:36 +0200 Subject: [PATCH 11/24] removed double space in output header --- src/msad.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/msad.jl b/src/msad.jl index 384ab3c..b4fbd2d 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -106,7 +106,7 @@ function Arianna.initialise(algorithm::MSADTracker, simulation::Simulation) # open output file — write header algorithm.files[c] = open(algorithm.paths[c], "w") - println(algorithm.files[c], "# t msad_euler msad_integral msad_thresh") + println(algorithm.files[c], "# t msad_euler msad_integral msad_thresh") # compute initial body frames R_all = get_all_body_frames(system) From 86108a418de3ae5b353905e29c753445d4f48e79 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 11:22:20 +0200 Subject: [PATCH 12/24] removed phi_total variable --- src/msad.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index b4fbd2d..5bed1ec 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -163,9 +163,9 @@ function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) for m in 1:N_mol dR = state.R_ref_thresh[m]' * R_all[m] phi_current[m] = rotation_vector(dR) - phi_total[m] = state.phi_acc[m] + phi_current[m] + #phi_total[m] = state.phi_acc[m] + phi_current[m] if norm(phi_current[m]) >= algorithm.theta_T - state.phi_acc[m] = phi_total[m] + state.phi_acc[m] = state.phi_acc[m] + phi_current[m] state.R_ref_thresh[m] = R_all[m] end end @@ -188,9 +188,14 @@ function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) # Threshold msad_thresh = sum( - norm(phi_total[m])^2 + norm(state.phi_acc[m] + phi_current[m])^2 for m in 1:N_mol ) / N_mol + + # msad_thresh = sum( + # norm(phi_total[m])^2 + # for m in 1:N_mol + # ) / N_mol println(algorithm.files[c], "$t $msad_euler $msad_integral $msad_thresh") flush(algorithm.files[c]) From c2f848d19c99037379a5efca23d8918973cce328 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 11:24:02 +0200 Subject: [PATCH 13/24] removed copy/paste mistake --- src/ParticlesMC.jl | 24 ------------------------ 1 file changed, 24 deletions(-) diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl index 607204e..a239e5d 100644 --- a/src/ParticlesMC.jl +++ b/src/ParticlesMC.jl @@ -242,30 +242,6 @@ ParticlesMC implemented in Comonicon. else error("Unsupported policy: $policy for action: $action") end - elseif alg == "MSADTracker" - # read theta_T from parameters - parameters = get(output, "parameters", Dict()) - theta_T = Float64(get(parameters, "theta_T", π/4)) - - # build output schedule separately from compute schedule - output_scheduler_params = get(output, "output_scheduler_params", Dict()) - if "log_base" in keys(output_scheduler_params) - output_sched = build_schedule(steps, burn, 2.0) - elseif "linear_interval" in keys(output_scheduler_params) - output_sched = build_schedule(steps, burn, - output_scheduler_params["linear_interval"]) - else - output_sched = sched # default: same as compute - end - - # convert to Vector{Int} — build_schedule returns this already - algorithm = ( - algorithm=MSADTracker, - scheduler=sched, # compute schedule - theta_T=theta_T, - output_schedule=output_sched, - path=output_path, - ) else error("Unsupported action: $action") From eab308448ed4a8659d76038bc20d7f938d1b8b9f Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 11:36:04 +0200 Subject: [PATCH 14/24] swap output and compute schedule for consistency with other amlgorithm ouput scheme --- src/ParticlesMC.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl index a239e5d..f876656 100644 --- a/src/ParticlesMC.jl +++ b/src/ParticlesMC.jl @@ -296,17 +296,19 @@ ParticlesMC implemented in Comonicon. parameters = get(output, "parameters", Dict()) theta_T = Float64(get(parameters, "theta_T", π/4)) - # build output schedule separately from compute schedule - output_scheduler_params = get(output, "output_scheduler_params", Dict()) - if "log_base" in keys(output_scheduler_params) - output_sched = build_schedule(steps, burn, 2.0) - elseif "linear_interval" in keys(output_scheduler_params) - output_sched = build_schedule(steps, burn, - output_scheduler_params["linear_interval"]) + # build output scheduler from scheduler_params for consistency + output_sched = sched + + # build compute schedule separately + compute_scheduler_params = get(output, "compute_scheduler_params", Dict()) + if "log_base" in keys(compute_scheduler_params) + compute_sched = build_schedule(steps, burn, Float64(compute_scheduler_params["log_base"])) + elseif "linear_interval" in keys(compute_scheduler_params) + compute_sched = build_schedule(steps, burn, compute_scheduler_params["linear_interval"]) else - output_sched = sched # default: same as compute + compute_sched = build_schedule(steps, burn, 1) # every step by default end - + # convert to Vector{Int} — build_schedule returns this already algorithm = ( algorithm=MSADTracker, From 39a863d1211e8d12580a962fc2aea89649f90d64 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 13:30:38 +0200 Subject: [PATCH 15/24] changed the rotation_vector function. eigenvector of R instead of skew matrix symmetric part. This version is more stable but might be slower --- Project.toml | 6 ++++-- src/msad.jl | 24 ++++++++++++++++----- test/bench_rotation.jl | 47 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 70 insertions(+), 7 deletions(-) create mode 100644 test/bench_rotation.jl diff --git a/Project.toml b/Project.toml index 4daa50b..2a011a6 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ version = "0.1.0" [deps] Arianna = "07692032-97b4-4f8d-80d7-e18df88d31a9" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Comonicon = "863f3e99-da2a-4334-8734-de3dacbe5542" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" @@ -19,6 +20,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" Arianna = "0.2" +BenchmarkTools = "1.8.0" Comonicon = "1.0" ComponentArrays = "0.15" ConcreteStructs = "0.2" @@ -26,19 +28,19 @@ DataStructures = "0.18" DelimitedFiles = "1.9" Distributions = "0.25" LinearAlgebra = "1.9" +Pkg = "1.9" Printf = "1.9" StaticArrays = "1.9" Statistics = "1.9" TOML = "1" Test = "1.9" julia = "1.9" -Pkg = "1.9" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test", "DelimitedFiles", "Aqua", "Pkg"] diff --git a/src/msad.jl b/src/msad.jl index 5bed1ec..9d83d14 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -40,13 +40,27 @@ function get_all_body_frames(system::Molecules) end function rotation_vector(R::SMatrix{3,3,T}) where {T} + cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) - # skew-symmetric part encodes sin(θ)*n - skew = (R - R') / 2 - sin_n = SVector(skew[3,2], skew[1,3], skew[2,1]) # julia starts index at 1 - sin_θ = norm(sin_n) + + vals, vecs = eigen(R) # eigenvalues and eigenvectors of R + idx = argmin(abs.(real.(vals) .- 1)) # find the idx corresponding to eigenvalue = 1 + n = real.(vecs[:, idx]) # extract the rotation axis vector + n = SVector{3,T}(n[1], n[2], n[3]) + + n_skew = SMatrix{3,3,T}( + 0, n[3], -n[2], + -n[3], 0, n[1], + n[2], -n[1], 0 + ) # skew matrix for n vector + + # Rodrigues formula + sin_θ = clamp(-tr(n_skew * R) / 2, -one(T), one(T)) + + # theta θ = atan(sin_θ, cos_θ) - return sin_θ > 1e-10 ? (θ / sin_θ) * sin_n : zero(SVector{3,T}) ## might need to add the theta = pi case (unprobable for MC steps) + + return θ * n end ##### Definition of the MSAD State ie: the accumulating variables ##### diff --git a/test/bench_rotation.jl b/test/bench_rotation.jl new file mode 100644 index 0000000..deb7aa0 --- /dev/null +++ b/test/bench_rotation.jl @@ -0,0 +1,47 @@ +using StaticArrays, LinearAlgebra, BenchmarkTools + +function rotation_vector_skew(R::SMatrix{3,3,T}) where {T} + cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) + skew = (R - R') / 2 + sin_n = SVector(skew[3,2], skew[1,3], skew[2,1]) + sin_θ = norm(sin_n) + θ = atan(sin_θ, cos_θ) + if sin_θ > 1e-6 + return (θ / sin_θ) * sin_n + elseif cos_θ > 0 + return sin_n + else + RpI = R + one(R) + c1 = SVector(RpI[1,1], RpI[2,1], RpI[3,1]) + c2 = SVector(RpI[1,2], RpI[2,2], RpI[3,2]) + c3 = SVector(RpI[1,3], RpI[2,3], RpI[3,3]) + n = norm(c1) >= norm(c2) ? c1 : c2 + n = norm(n) >= norm(c3) ? n : c3 + return π * (n / norm(n)) + end +end + +function rotation_vector_eig(R::SMatrix{3,3,T}) where {T} + cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) + vals, vecs = eigen(R) + idx = argmin(abs.(real.(vals) .- 1)) + n = real.(vecs[:, idx]) + n = SVector{3,T}(n[1], n[2], n[3]) + n_skew = SMatrix{3,3,T}( + 0, n[3], -n[2], + -n[3], 0, n[1], + n[2], -n[1], 0 + ) + sin_θ = clamp(-tr(n_skew * R) / 2, -one(T), one(T)) + θ = atan(sin_θ, cos_θ) + return θ * n +end + +# 90° rotation around z — normal case +R = SMatrix{3,3,Float64}(0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0) + +println("=== Skew-symmetric ===") +@btime rotation_vector_skew($R) + +println("=== Eigenvector ===") +@btime rotation_vector_eig($R) From e25b791459a54cdc2460fed533fede5fe6c7f5f4 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 14:04:08 +0200 Subject: [PATCH 16/24] changed the output from msad to phi vectors for integral and threshold method. I add to use again the phi_total variable... also I've dropped phi_euler as it can be computed directly from the configuration --- src/msad.jl | 77 ++++++++++++++++++++++++++++------------------------- 1 file changed, 40 insertions(+), 37 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index 9d83d14..bd9b0ee 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -40,7 +40,7 @@ function get_all_body_frames(system::Molecules) end function rotation_vector(R::SMatrix{3,3,T}) where {T} - + cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) vals, vecs = eigen(R) # eigenvalues and eigenvectors of R @@ -87,8 +87,9 @@ mutable struct MSADTracker{T} <: AriannaAlgorithm theta_T::T # threshold angle in radians compute_schedule::Vector{Int} # schedule for computing rotation for integral and threshold methods output_schedule::Vector{Int} # schedule for writing MSAD to disk - paths::Vector{String} - files::Vector{IOStream} + paths::Vector{String} # paths to store files + files_integral::Vector{IOStream} # file to write rotation vector in integral method + files_thresh::Vector{IOStream} # file to write rotation vector in threshold method end function MSADTracker(chains, theta_T::T, compute_schedule::Vector{Int}, output_schedule::Vector{Int}, path::String) where {T} @@ -99,11 +100,29 @@ function MSADTracker(chains, theta_T::T, compute_schedule::Vector{Int}, output_s states = [MSADState{T}() for _ in 1:n] # one empty state per chain filled during initialise dirs = [joinpath(path, "chains", "$c") for c in 1:n] - paths = [joinpath(d, "msad.dat") for d in dirs] - files = Vector{IOStream}(undef, n) - return MSADTracker{T}(states, theta_T, compute_schedule, output_schedule, paths, files) + + paths_integral = [joinpath(d, "phi_integral.dat") for d in dirs] + paths_thresh = [joinpath(d, "phi_thresh.dat") for d in dirs] + + files_integral = Vector{IOStream}(undef, n) + files_thresh = Vector{IOStream}(undef, n) + + return MSADTracker{T}(states, theta_T, compute_schedule, output_schedule, paths_integral, paths_thresh, files_integral, files_thresh) end +##### Function to write phi frame ##### +########## ########## + +function write_phi_frame(file::IOStream, t::Int, N_mol::Int, + phis::Vector{<:SVector}) + println(file, "t=$t") + for m in 1:N_mol + println(file, "$m $(phis[m][1]) $(phis[m][2]) $(phis[m][3])") + end + flush(file) +end + + ##### Initialisation of the simulation ##### # called once at t = 0, set the system, the mutable struct and the storing paths ########## ########## @@ -116,10 +135,12 @@ function Arianna.initialise(algorithm::MSADTracker, simulation::Simulation) T = typeof(system.temperature) # get the float type from system so one can choose Float64 or 32 (generic) # create output directory if it doesn't exist - mkpath(dirname(algorithm.paths[c])) + mkpath(dirname(algorithm.paths_integral[c])) # open output file — write header - algorithm.files[c] = open(algorithm.paths[c], "w") + algorithm.files_integral[c] = open(algorithm.paths_integral[c], "w") + algorithm.files_thresh[c] = open(algorithm.paths_thresh[c], "w") + println(algorithm.files[c], "# t msad_euler msad_integral msad_thresh") # compute initial body frames @@ -134,9 +155,9 @@ function Arianna.initialise(algorithm::MSADTracker, simulation::Simulation) state.phi_acc = [zero(SVector{3,T}) for _ in 1:N_mol] state.initialized = true - # write t=0, all MSAD values are zero at start - println(algorithm.files[c], "0 0.0 0.0 0.0") - flush(algorithm.files[c]) + # write t=0, all phi vectors + write_phi_frame(algorithm.files_integral[c],0,N_mol,state.phi_integral) + write_phi_frame(algorithm.files_thres[c],0,N_mol,state.phi_thresh) end end @@ -145,9 +166,10 @@ end # close output to not keep unwritten data in memory ########## ########## -function Arianna.finalise(algorithm::MSADTracker, ::Simulation) - for f in algorithm.files - close(f) +function Arianna.finalise(tracker::MSADTracker, ::Simulation) + for c in eachindex(tracker.files_integral) + close(tracker.files_integral[c]) + close(tracker.files_thresh[c]) end end @@ -177,9 +199,9 @@ function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) for m in 1:N_mol dR = state.R_ref_thresh[m]' * R_all[m] phi_current[m] = rotation_vector(dR) - #phi_total[m] = state.phi_acc[m] + phi_current[m] + phi_total[m] = state.phi_acc[m] + phi_current[m] if norm(phi_current[m]) >= algorithm.theta_T - state.phi_acc[m] = state.phi_acc[m] + phi_current[m] + state.phi_acc[m] = state.phi_total[m] state.R_ref_thresh[m] = R_all[m] end end @@ -188,31 +210,12 @@ function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) # output is time match output schedule if t in algorithm.output_schedule - # Euler only computed here not accumulated - msad_euler = sum( - norm(rotation_vector(state.R0[m]' * R_all[m]))^2 - for m in 1:N_mol - ) / N_mol # Integral - msad_integral = sum( - norm(state.phi_integral[m])^2 - for m in 1:N_mol - ) / N_mol + write_phi_frame(algorithm.files_integral[c],t,N_mol,state.phi_integral) # Threshold - msad_thresh = sum( - norm(state.phi_acc[m] + phi_current[m])^2 - for m in 1:N_mol - ) / N_mol - - # msad_thresh = sum( - # norm(phi_total[m])^2 - # for m in 1:N_mol - # ) / N_mol - - println(algorithm.files[c], "$t $msad_euler $msad_integral $msad_thresh") - flush(algorithm.files[c]) + write_phi_frame(algorithm.files_integral[c],t,N_mol,state.phi_total) end end end \ No newline at end of file From e98224d5d61042c019ef17f24de33335dd282d4f Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 14:25:29 +0200 Subject: [PATCH 17/24] fixed bugs in MSADTracker --- src/msad.jl | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index bd9b0ee..9f0a061 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -92,9 +92,16 @@ mutable struct MSADTracker{T} <: AriannaAlgorithm files_thresh::Vector{IOStream} # file to write rotation vector in threshold method end -function MSADTracker(chains, theta_T::T, compute_schedule::Vector{Int}, output_schedule::Vector{Int}, path::String) where {T} +function MSADTracker(chains; theta_T::T, output_schedule::Vector{Int}, path::String, scheduler::Vector{Int}=Int[], kwargs...) where {T} - @assert all(t in compute_schedule for t in output_schedule) # output_schedule must be a subset of compute_schedule + compute_schedule = scheduler + + if !isempty(compute_schedule) && !isempty(output_schedule) + @assert all(t in compute_schedule for t in output_schedule) """ + output_schedule contains steps not in compute_schedule. + You cannot write output at a step where make_step! was not called. + """ + end n = length(chains) states = [MSADState{T}() for _ in 1:n] # one empty state per chain filled during initialise @@ -208,7 +215,7 @@ function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) state.R_prev .= R_all - # output is time match output schedule + # output if time match output schedule if t in algorithm.output_schedule # Integral From ee390e9da2f82b7f01b8287bbfc40ae784036528 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 14:43:22 +0200 Subject: [PATCH 18/24] fix path mismatch --- src/msad.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index 9f0a061..43b21a7 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -87,7 +87,7 @@ mutable struct MSADTracker{T} <: AriannaAlgorithm theta_T::T # threshold angle in radians compute_schedule::Vector{Int} # schedule for computing rotation for integral and threshold methods output_schedule::Vector{Int} # schedule for writing MSAD to disk - paths::Vector{String} # paths to store files + path::Vector{String} # paths to store files files_integral::Vector{IOStream} # file to write rotation vector in integral method files_thresh::Vector{IOStream} # file to write rotation vector in threshold method end @@ -108,11 +108,8 @@ function MSADTracker(chains; theta_T::T, output_schedule::Vector{Int}, path::Str dirs = [joinpath(path, "chains", "$c") for c in 1:n] - paths_integral = [joinpath(d, "phi_integral.dat") for d in dirs] - paths_thresh = [joinpath(d, "phi_thresh.dat") for d in dirs] - - files_integral = Vector{IOStream}(undef, n) - files_thresh = Vector{IOStream}(undef, n) + algorithm.files_integral[c] = open(joinpath(dir, "phi_integral.dat"), "w") + algorithm.files_thresh[c] = open(joinpath(dir, "phi_thresh.dat"), "w") return MSADTracker{T}(states, theta_T, compute_schedule, output_schedule, paths_integral, paths_thresh, files_integral, files_thresh) end From 1f4527c25ae9c4f432d2c5731e2ad4e41f194c13 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 14:47:25 +0200 Subject: [PATCH 19/24] fix path mismatch bis --- src/msad.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index 43b21a7..66406db 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -108,8 +108,8 @@ function MSADTracker(chains; theta_T::T, output_schedule::Vector{Int}, path::Str dirs = [joinpath(path, "chains", "$c") for c in 1:n] - algorithm.files_integral[c] = open(joinpath(dir, "phi_integral.dat"), "w") - algorithm.files_thresh[c] = open(joinpath(dir, "phi_thresh.dat"), "w") + algorithm.files_integral[c] = open(joinpath(dirs, "phi_integral.dat"), "w") + algorithm.files_thresh[c] = open(joinpath(dirs, "phi_thresh.dat"), "w") return MSADTracker{T}(states, theta_T, compute_schedule, output_schedule, paths_integral, paths_thresh, files_integral, files_thresh) end From 2071047746aa26dccbac317968de3091b1056407 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 14:55:33 +0200 Subject: [PATCH 20/24] still fixing paths issues... --- src/msad.jl | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index 66406db..aac47b3 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -92,7 +92,12 @@ mutable struct MSADTracker{T} <: AriannaAlgorithm files_thresh::Vector{IOStream} # file to write rotation vector in threshold method end -function MSADTracker(chains; theta_T::T, output_schedule::Vector{Int}, path::String, scheduler::Vector{Int}=Int[], kwargs...) where {T} +function MSADTracker(chains; + theta_T::Float64, + output_schedule::Vector{Int}=Int[], + path::String=".", + scheduler::Vector{Int}=Int[], + kwargs...) compute_schedule = scheduler @@ -103,15 +108,15 @@ function MSADTracker(chains; theta_T::T, output_schedule::Vector{Int}, path::Str """ end - n = length(chains) - states = [MSADState{T}() for _ in 1:n] # one empty state per chain filled during initialise + n = length(chains) + states = [MSADState{Float64}() for _ in 1:n] - dirs = [joinpath(path, "chains", "$c") for c in 1:n] + files_integral = Vector{IOStream}(undef, n) + files_thresh = Vector{IOStream}(undef, n) - algorithm.files_integral[c] = open(joinpath(dirs, "phi_integral.dat"), "w") - algorithm.files_thresh[c] = open(joinpath(dirs, "phi_thresh.dat"), "w") - - return MSADTracker{T}(states, theta_T, compute_schedule, output_schedule, paths_integral, paths_thresh, files_integral, files_thresh) + return MSADTracker{Float64}(states, theta_T, compute_schedule, + output_schedule, path, + files_integral, files_thresh) end ##### Function to write phi frame ##### @@ -139,13 +144,12 @@ function Arianna.initialise(algorithm::MSADTracker, simulation::Simulation) T = typeof(system.temperature) # get the float type from system so one can choose Float64 or 32 (generic) # create output directory if it doesn't exist - mkpath(dirname(algorithm.paths_integral[c])) + dir = joinpath(algrithm.path, "chains", "$c") + mkpath(dir) # open output file — write header - algorithm.files_integral[c] = open(algorithm.paths_integral[c], "w") - algorithm.files_thresh[c] = open(algorithm.paths_thresh[c], "w") - - println(algorithm.files[c], "# t msad_euler msad_integral msad_thresh") + algorithm.files_integral[c] = open(joinpath(dir, "phi_integral.dat"), "w") + algorithm.files_thresh[c] = open(joinpath(dir, "phi_thresh.dat"), "w") # compute initial body frames R_all = get_all_body_frames(system) From dc6aa120f2d6f4c3c97807808bce0d0538f23ea4 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 15:12:19 +0200 Subject: [PATCH 21/24] still working on paths --- src/msad.jl | 50 +++++++++++++++++++++++++++----------------------- 1 file changed, 27 insertions(+), 23 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index aac47b3..06831fb 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -83,17 +83,18 @@ MSADState{T}() where {T} = MSADState{T}([], [], [], [], [], false) ########## ########## mutable struct MSADTracker{T} <: AriannaAlgorithm - states::Vector{MSADState{T}} # the states of the considered chain - theta_T::T # threshold angle in radians - compute_schedule::Vector{Int} # schedule for computing rotation for integral and threshold methods - output_schedule::Vector{Int} # schedule for writing MSAD to disk - path::Vector{String} # paths to store files - files_integral::Vector{IOStream} # file to write rotation vector in integral method - files_thresh::Vector{IOStream} # file to write rotation vector in threshold method + states::Vector{MSADState{T}} + theta_T::T + compute_schedule::Vector{Int} + output_schedule::Vector{Int} + paths_integral::Vector{String} + paths_thresh::Vector{String} + files_integral::Vector{IOStream} + files_thresh::Vector{IOStream} end function MSADTracker(chains; - theta_T::Float64, + theta_T::Float64=π/4, output_schedule::Vector{Int}=Int[], path::String=".", scheduler::Vector{Int}=Int[], @@ -110,12 +111,16 @@ function MSADTracker(chains; n = length(chains) states = [MSADState{Float64}() for _ in 1:n] + dirs = [joinpath(path, "chains", "$c") for c in 1:n] + + paths_integral = [joinpath(d, "phi_integral.dat") for d in dirs] + paths_thresh = [joinpath(d, "phi_thresh.dat") for d in dirs] files_integral = Vector{IOStream}(undef, n) files_thresh = Vector{IOStream}(undef, n) return MSADTracker{Float64}(states, theta_T, compute_schedule, - output_schedule, path, + output_schedule, paths_integral, paths_thresh, files_integral, files_thresh) end @@ -124,6 +129,7 @@ end function write_phi_frame(file::IOStream, t::Int, N_mol::Int, phis::Vector{<:SVector}) + println(file, N_mol) println(file, "t=$t") for m in 1:N_mol println(file, "$m $(phis[m][1]) $(phis[m][2]) $(phis[m][3])") @@ -131,7 +137,6 @@ function write_phi_frame(file::IOStream, t::Int, N_mol::Int, flush(file) end - ##### Initialisation of the simulation ##### # called once at t = 0, set the system, the mutable struct and the storing paths ########## ########## @@ -141,21 +146,20 @@ function Arianna.initialise(algorithm::MSADTracker, simulation::Simulation) for c in eachindex(simulation.chains) system = simulation.chains[c] state = algorithm.states[c] - T = typeof(system.temperature) # get the float type from system so one can choose Float64 or 32 (generic) + T = typeof(system.temperature) # create output directory if it doesn't exist - dir = joinpath(algrithm.path, "chains", "$c") - mkpath(dir) + mkpath(dirname(algorithm.paths_integral[c])) - # open output file — write header - algorithm.files_integral[c] = open(joinpath(dir, "phi_integral.dat"), "w") - algorithm.files_thresh[c] = open(joinpath(dir, "phi_thresh.dat"), "w") + # open output files + algorithm.files_integral[c] = open(algorithm.paths_integral[c], "w") + algorithm.files_thresh[c] = open(algorithm.paths_thresh[c], "w") # compute initial body frames R_all = get_all_body_frames(system) N_mol = system.Nmol - # fill state, copy. so they are different pointers! + # fill state state.R0 = copy(R_all) state.R_prev = copy(R_all) state.phi_integral = [zero(SVector{3,T}) for _ in 1:N_mol] @@ -163,14 +167,15 @@ function Arianna.initialise(algorithm::MSADTracker, simulation::Simulation) state.phi_acc = [zero(SVector{3,T}) for _ in 1:N_mol] state.initialized = true - # write t=0, all phi vectors - write_phi_frame(algorithm.files_integral[c],0,N_mol,state.phi_integral) - write_phi_frame(algorithm.files_thres[c],0,N_mol,state.phi_thresh) + # write t=0, all phi vectors are zero + write_phi_frame(algorithm.files_integral[c], 0, N_mol, + state.phi_integral) + write_phi_frame(algorithm.files_thresh[c], 0, N_mol, + [zero(SVector{3,T}) for _ in 1:N_mol]) end - end -##### Finalise th simulation ##### +##### Finalise the simulation ##### # close output to not keep unwritten data in memory ########## ########## @@ -180,7 +185,6 @@ function Arianna.finalise(tracker::MSADTracker, ::Simulation) close(tracker.files_thresh[c]) end end - ##### Make a step in simulation ##### ########## ########## From 8efb845b2fb452b72cc4766879c432fed898838f Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 15:20:31 +0200 Subject: [PATCH 22/24] phi_total is local not global memory --- src/msad.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/msad.jl b/src/msad.jl index 06831fb..4d52aa1 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -213,7 +213,7 @@ function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) phi_current[m] = rotation_vector(dR) phi_total[m] = state.phi_acc[m] + phi_current[m] if norm(phi_current[m]) >= algorithm.theta_T - state.phi_acc[m] = state.phi_total[m] + state.phi_acc[m] = phi_total[m] state.R_ref_thresh[m] = R_all[m] end end @@ -227,7 +227,7 @@ function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) write_phi_frame(algorithm.files_integral[c],t,N_mol,state.phi_integral) # Threshold - write_phi_frame(algorithm.files_integral[c],t,N_mol,state.phi_total) + write_phi_frame(algorithm.files_integral[c],t,N_mol,phi_total) end end end \ No newline at end of file From 8efddc84b5655fa3d00cd6d16df26f67b78fe4e3 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Fri, 15 May 2026 15:26:12 +0200 Subject: [PATCH 23/24] typo in files_thresh name --- src/msad.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/msad.jl b/src/msad.jl index 4d52aa1..7d4e94b 100644 --- a/src/msad.jl +++ b/src/msad.jl @@ -227,7 +227,7 @@ function Arianna.make_step!(simulation::Simulation, algorithm::MSADTracker) write_phi_frame(algorithm.files_integral[c],t,N_mol,state.phi_integral) # Threshold - write_phi_frame(algorithm.files_integral[c],t,N_mol,phi_total) + write_phi_frame(algorithm.files_thresh[c],t,N_mol,phi_total) end end end \ No newline at end of file From e2d510ff6ee29ba217645095582e4499838bec7f Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Wed, 20 May 2026 14:45:33 +0200 Subject: [PATCH 24/24] fix scheduler mismatch --- src/ParticlesMC.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl index f876656..6ff9a1b 100644 --- a/src/ParticlesMC.jl +++ b/src/ParticlesMC.jl @@ -312,7 +312,7 @@ ParticlesMC implemented in Comonicon. # convert to Vector{Int} — build_schedule returns this already algorithm = ( algorithm=MSADTracker, - scheduler=sched, # compute schedule + scheduler=compute_sched, # compute schedule theta_T=theta_T, output_schedule=output_sched, path=output_path,