From 2cccfb84749ffd27f68c6558eaae72b3be3e5193 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 9 Jun 2026 10:10:57 +0200 Subject: [PATCH 1/4] try to improve efficiency without using eigen but direct calculation optimise allocation of bodyframes --- src/rotation.jl | 96 +++++++++++++++++++++-------------- test.jl | 131 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 188 insertions(+), 39 deletions(-) create mode 100644 test.jl diff --git a/src/rotation.jl b/src/rotation.jl index 64daf76..e6e51fe 100644 --- a/src/rotation.jl +++ b/src/rotation.jl @@ -27,40 +27,60 @@ function body_frame(r1::SVector{3,T}, r2::SVector{3,T}, r3::SVector{3,T}, L::T) 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 +function get_all_body_frames!(R_bodyframe, system::Molecules) + L = system.box[1] 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) + @inbounds for m in 1:system.Nmol + s = system.start_mol[m] + R_bodyframe[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} +# function rotation_vector(R::SMatrix{3,3,T}) where {T} - cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) +# cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) + +# 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]) - 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 - 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)) - # Rodrigues formula - sin_θ = clamp(-tr(n_skew * R) / 2, -one(T), one(T)) +# # theta +# θ = atan(sin_θ, cos_θ) - # theta - θ = atan(sin_θ, cos_θ) +# return θ * n +# end - return θ * n +function rotation_vector(R::SMatrix{3,3,T}) where {T} + cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) + θ = acos(cos_θ) + + ax = R[3,2] - R[2,3] + ay = R[1,3] - R[3,1] + az = R[2,1] - R[1,2] + sin_θ = sqrt(ax^2 + ay^2 + az^2) / 2 + + if θ < sqrt(eps(T)) # θ ≈ 0 + return zero(SVector{3,T}) + elseif π - θ < sqrt(eps(T)) # θ ≈ π + nx = sqrt(max((R[1,1] + 1) / 2, zero(T))) + ny = sqrt(max((R[2,2] + 1) / 2, zero(T))) + nz = sqrt(max((R[3,3] + 1) / 2, zero(T))) + n = SVector{3,T}(nx, ny, nz) + return θ * n / norm(n) + else # general case — zero allocs + inv_2s = θ / (2 * sin_θ) + return SVector{3,T}(ax * inv_2s, ay * inv_2s, az * inv_2s) + end end ##### Definition of the RotationState ie: the accumulating variables ##### @@ -68,13 +88,14 @@ end ########## ########## mutable struct RotationState{T} - R_ref::Vector{Vector{SMatrix{3,3,T,9}}} + R_ref::Vector{Vector{SMatrix{3,3,T,9}}} + R_bodyframe::Vector{SMatrix{3,3,T,9}} Φ_acc::Vector{Vector{SVector{3,T}}} initialized::Bool end # Start empty, filled during initialise when we first see the system -RotationState{T}() where {T} = RotationState{T}([], [], false) +RotationState{T}() where {T} = RotationState{T}([], [], [], false) ##### The Rotation computation algorithm ##### ########## ########## @@ -98,7 +119,6 @@ end ########## ########## function Arianna.initialise(algorithm::ComputeRotation, simulation::Simulation) - for c in eachindex(simulation.chains) system = simulation.chains[c] state = algorithm.states[c] @@ -106,15 +126,14 @@ function Arianna.initialise(algorithm::ComputeRotation, simulation::Simulation) N_mol = system.Nmol n_θ = length(algorithm.theta_T) - # compute initial body frames - R_all = get_all_body_frames(system) + state.R_bodyframe = Vector{SMatrix{3,3,T,9}}(undef, N_mol) # allocate once + get_all_body_frames!(state.R_bodyframe, system) # ← state.R_buf - # fill state - state.R_ref = [copy(R_all) for _ in 1:n_θ] - state.Φ_acc = [[zero(SVector{3,T}) for _ in 1:N_mol] for _ in 1:n_θ] - state.initialized = true - - resize!(system.Φ, n_θ) # a posteriori as we know n_θ + state.R_ref = [copy(state.R_bodyframe) for _ in 1:n_θ] + state.Φ_acc = [[zero(SVector{3,T}) for _ in 1:N_mol] for _ in 1:n_θ] + state.initialized = true + + resize!(system.Φ, n_θ) for k in 1:n_θ system.Φ[k] = [zero(SVector{3,T}) for _ in 1:N_mol] end @@ -125,20 +144,19 @@ end ########## ########## function Arianna.make_step!(simulation::Simulation, algorithm::ComputeRotation) - tcollect(eachindex(simulation.chains) |> Transducers.Map(c -> begin system = simulation.chains[c] state = algorithm.states[c] N_mol = system.Nmol - R_all = get_all_body_frames(system) + get_all_body_frames!(state.R_bodyframe, system) # ← in-place, no alloc # ← alias, no copy for (k, θ_T) in enumerate(algorithm.theta_T) for m in 1:N_mol - dR = state.R_ref[k][m]' * R_all[m] + dR = state.R_ref[k][m]' * state.R_bodyframe[m] Φ_current = rotation_vector(dR) Φ_total = state.Φ_acc[k][m] + Φ_current if norm(Φ_current) > θ_T state.Φ_acc[k][m] = Φ_total - state.R_ref[k][m] = R_all[m] + state.R_ref[k][m] = state.R_bodyframe[m] end system.Φ[k][m] = Φ_total end diff --git a/test.jl b/test.jl new file mode 100644 index 0000000..fd97b16 --- /dev/null +++ b/test.jl @@ -0,0 +1,131 @@ +using LinearAlgebra, StaticArrays, BenchmarkTools + +# ── original ────────────────────────────────────────────────────────────────── +function rotation_vector_eigen(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 + +# ── closed-form ─────────────────────────────────────────────────────────────── +function rotation_vector_closed(R::SMatrix{3,3,T}) where {T} + cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) + θ = acos(cos_θ) + + ax = R[3,2] - R[2,3] + ay = R[1,3] - R[3,1] + az = R[2,1] - R[1,2] + sin_θ = sqrt(ax^2 + ay^2 + az^2) / 2 # exact: norm of 2·sin(θ)·n / 2 + + if θ < sqrt(eps(T)) # θ ≈ 0 : no rotation + return zero(SVector{3,T}) + elseif π - θ < sqrt(eps(T)) # θ ≈ π : antisymmetric part vanishes + # use diagonal of symmetric part: (R+I)/2 = n⊗n + nx = sqrt(max((R[1,1] + 1) / 2, zero(T))) + ny = sqrt(max((R[2,2] + 1) / 2, zero(T))) + nz = sqrt(max((R[3,3] + 1) / 2, zero(T))) + n = SVector{3,T}(nx, ny, nz) + return θ * n / norm(n) + else # general case: zero allocs + inv_2s = θ / (2 * sin_θ) + return SVector{3,T}(ax * inv_2s, ay * inv_2s, az * inv_2s) + end +end + +# ── helper: build exact rotation matrix from axis + angle ───────────────────── +function make_rotation(axis::SVector{3,Float64}, θ::Float64) + n = axis / norm(axis) + # skew matrix stored column-major for SMatrix + N = SMatrix{3,3,Float64}( + 0, n[3], -n[2], + -n[3], 0, n[1], + n[2], -n[1], 0 + ) + return SMatrix{3,3,Float64}(I + sin(θ)*N + (1 - cos(θ))*(N*N)) +end + +# ── Test 1: 10_000 random rotations, θ ∈ (0.01, π-0.01) ───────────────────── +println("="^60) +println("Test 1 — 10_000 random rotations, θ ∈ (0.01, π-0.01)") +println("="^60) +max_err = 0.0 +for _ in 1:10_000 + axis = SVector{3,Float64}(randn(), randn(), randn()) + θ_true = 0.01 + (π - 0.02) * rand() + R = make_rotation(axis, θ_true) + v1 = rotation_vector_eigen(R) + v2 = rotation_vector_closed(R) + err = norm(v1 - v2) + global max_err = max(max_err, err) +end +println("Max |eigen - closed|: ", max_err, " (should be < 1e-12)") + +# ── Test 2: near-zero (θ → 0) ──────────────────────────────────────────────── +println() +println("="^60) +println("Test 2 — near-zero rotation θ ≈ 0") +println("="^60) +for θ_small in [1e-4, 1e-6, 1e-8, 1e-12] + axis = SVector{3,Float64}(1.0, 0.0, 0.0) + R = make_rotation(axis, θ_small) + v1 = rotation_vector_eigen(R) + v2 = rotation_vector_closed(R) + println("θ=%.0e |v_eigen|=%.6e |v_closed|=%.6e diff=%.2e\n", + θ_small, norm(v1), norm(v2), norm(v1-v2)) +end + +# ── Test 3: near-π — RELEVANT because theta_T=3.0 rad ──────────────────────── +println() +println("="^60) +println("Test 3 — near-π rotation (theta_T=3.0 means you DO reach this)") +println("="^60) +for δ in [1e-1, 1e-3, 1e-5, 1e-7] + axis = normalize(SVector{3,Float64}(1.0, 2.0, 3.0)) + R = make_rotation(axis, π - δ) + v1 = rotation_vector_eigen(R) + v2 = rotation_vector_closed(R) + println("θ=π-%.0e |v_eigen|=%.6f |v_closed|=%.6f diff=%.2e\n", + δ, norm(v1), norm(v2), norm(v1-v2)) +end + +# ── Test 4: benchmark ───────────────────────────────────────────────────────── +println() +println("="^60) +println("Test 4 — benchmark: 1000 calls each") +println("="^60) +Rs = [make_rotation(normalize(SVector{3,Float64}(randn(),randn(),randn())), + 0.1 + (π-0.2)*rand()) + for _ in 1:1000] + +t_eigen = @belapsed sum(norm, rotation_vector_eigen.($(Rs))) +t_closed = @belapsed sum(norm, rotation_vector_closed.($(Rs))) +println("eigen : %8.2f µs\n", t_eigen * 1e6) +println("closed : %8.2f µs\n", t_closed * 1e6) +println("speedup: %.1fx\n", t_eigen / t_closed) + +# ── Test 5: verify the identity R - Rᵀ = 2 sin(θ) N analytically ───────────── +println() +println("="^60) +println("Test 5 — verify R - Rᵀ = 2·sin(θ)·N identity directly") +println("="^60) +axis = normalize(SVector{3,Float64}(1.0, 2.0, 3.0)) +θ_ref = 1.23456 +R = make_rotation(axis, θ_ref) +antisym = (R - R') / 2 # = sin(θ)·N +n_from_antisym = SVector{3,Float64}( + antisym[3,2], # = sin(θ)·n[1] + antisym[1,3], # = sin(θ)·n[2] + antisym[2,1], # = sin(θ)·n[3] +) +n_recovered = n_from_antisym / norm(n_from_antisym) +println("True axis : ", round.(axis, sigdigits=8)) +println("Recovered axis: ", round.(n_recovered, sigdigits=8)) +println("Error on axis : ", norm(axis - n_recovered)) +println("True θ : ", θ_ref) +println("Recovered θ : ", atan(norm(n_from_antisym), (tr(R)-1)/2)) \ No newline at end of file From 7ca66ec02bb3caa0979f231e6c944e996a1b9d4a Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Mon, 15 Jun 2026 16:57:44 +0200 Subject: [PATCH 2/4] fix ill defined operation close to identity rotation --- src/rotation.jl | 54 ++++++++++++++++++------------------------------- 1 file changed, 20 insertions(+), 34 deletions(-) diff --git a/src/rotation.jl b/src/rotation.jl index e6e51fe..728cd8c 100644 --- a/src/rotation.jl +++ b/src/rotation.jl @@ -36,30 +36,6 @@ function get_all_body_frames!(R_bodyframe, system::Molecules) end 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 -# 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 θ * n -# end - function rotation_vector(R::SMatrix{3,3,T}) where {T} cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) θ = acos(cos_θ) @@ -69,17 +45,27 @@ function rotation_vector(R::SMatrix{3,3,T}) where {T} az = R[2,1] - R[1,2] sin_θ = sqrt(ax^2 + ay^2 + az^2) / 2 - if θ < sqrt(eps(T)) # θ ≈ 0 - return zero(SVector{3,T}) - elseif π - θ < sqrt(eps(T)) # θ ≈ π - nx = sqrt(max((R[1,1] + 1) / 2, zero(T))) - ny = sqrt(max((R[2,2] + 1) / 2, zero(T))) - nz = sqrt(max((R[3,3] + 1) / 2, zero(T))) - n = SVector{3,T}(nx, ny, nz) - return θ * n / norm(n) - else # general case — zero allocs + if sin_θ > sqrt(eps(T)) # guard sin_θ ≈ 0 + θ = atan(sin_θ, cos_θ) inv_2s = θ / (2 * sin_θ) return SVector{3,T}(ax * inv_2s, ay * inv_2s, az * inv_2s) + + elseif cos_θ > 0 # guard θ ≈ 0, near identity + return zero(SVector{3,T}) + + else # guard θ ≈ π + if R[1,1] >= R[2,2] && R[1,1] >= R[3,3] + nx = sqrt(max((R[1,1] + 1) / 2, zero(T))) + n = SVector{3,T}(nx, (R[1,2] + R[2,1]) / (4nx), (R[1,3] + R[3,1]) / (4nx)) + elseif R[2,2] >= R[3,3] + ny = sqrt(max((R[2,2] + 1) / 2, zero(T))) + n = SVector{3,T}((R[1,2] + R[2,1]) / (4ny), ny, (R[2,3] + R[3,2]) / (4ny)) + else + nz = sqrt(max((R[3,3] + 1) / 2, zero(T))) + n = SVector{3,T}((R[1,3] + R[3,1]) / (4nz), (R[2,3] + R[3,2]) / (4nz), nz) + end + θ = atan(sin_θ, cos_θ) + return θ * n / norm(n) end end @@ -175,7 +161,7 @@ function write_phi_frame(file::IOStream, t::Int, N_mol::Int, Φs::Vector{<:SVect for m in 1:N_mol println(file, "$m $(Φs[m][1]) $(Φs[m][2]) $(Φs[m][3])") end - flush(file) + #flush(file) not flushing at everyframe end struct StorePhiTrajectories <: AriannaAlgorithm From 3028f9a709060016d2b819ddaed97d5785662028 Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 16 Jun 2026 10:01:14 +0200 Subject: [PATCH 3/4] removed dead code --- src/rotation.jl | 1 - test.jl | 131 ------------------------------------------------ 2 files changed, 132 deletions(-) delete mode 100644 test.jl diff --git a/src/rotation.jl b/src/rotation.jl index 728cd8c..4ae2173 100644 --- a/src/rotation.jl +++ b/src/rotation.jl @@ -38,7 +38,6 @@ end function rotation_vector(R::SMatrix{3,3,T}) where {T} cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) - θ = acos(cos_θ) ax = R[3,2] - R[2,3] ay = R[1,3] - R[3,1] diff --git a/test.jl b/test.jl deleted file mode 100644 index fd97b16..0000000 --- a/test.jl +++ /dev/null @@ -1,131 +0,0 @@ -using LinearAlgebra, StaticArrays, BenchmarkTools - -# ── original ────────────────────────────────────────────────────────────────── -function rotation_vector_eigen(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 - -# ── closed-form ─────────────────────────────────────────────────────────────── -function rotation_vector_closed(R::SMatrix{3,3,T}) where {T} - cos_θ = clamp((tr(R) - 1) / 2, -one(T), one(T)) - θ = acos(cos_θ) - - ax = R[3,2] - R[2,3] - ay = R[1,3] - R[3,1] - az = R[2,1] - R[1,2] - sin_θ = sqrt(ax^2 + ay^2 + az^2) / 2 # exact: norm of 2·sin(θ)·n / 2 - - if θ < sqrt(eps(T)) # θ ≈ 0 : no rotation - return zero(SVector{3,T}) - elseif π - θ < sqrt(eps(T)) # θ ≈ π : antisymmetric part vanishes - # use diagonal of symmetric part: (R+I)/2 = n⊗n - nx = sqrt(max((R[1,1] + 1) / 2, zero(T))) - ny = sqrt(max((R[2,2] + 1) / 2, zero(T))) - nz = sqrt(max((R[3,3] + 1) / 2, zero(T))) - n = SVector{3,T}(nx, ny, nz) - return θ * n / norm(n) - else # general case: zero allocs - inv_2s = θ / (2 * sin_θ) - return SVector{3,T}(ax * inv_2s, ay * inv_2s, az * inv_2s) - end -end - -# ── helper: build exact rotation matrix from axis + angle ───────────────────── -function make_rotation(axis::SVector{3,Float64}, θ::Float64) - n = axis / norm(axis) - # skew matrix stored column-major for SMatrix - N = SMatrix{3,3,Float64}( - 0, n[3], -n[2], - -n[3], 0, n[1], - n[2], -n[1], 0 - ) - return SMatrix{3,3,Float64}(I + sin(θ)*N + (1 - cos(θ))*(N*N)) -end - -# ── Test 1: 10_000 random rotations, θ ∈ (0.01, π-0.01) ───────────────────── -println("="^60) -println("Test 1 — 10_000 random rotations, θ ∈ (0.01, π-0.01)") -println("="^60) -max_err = 0.0 -for _ in 1:10_000 - axis = SVector{3,Float64}(randn(), randn(), randn()) - θ_true = 0.01 + (π - 0.02) * rand() - R = make_rotation(axis, θ_true) - v1 = rotation_vector_eigen(R) - v2 = rotation_vector_closed(R) - err = norm(v1 - v2) - global max_err = max(max_err, err) -end -println("Max |eigen - closed|: ", max_err, " (should be < 1e-12)") - -# ── Test 2: near-zero (θ → 0) ──────────────────────────────────────────────── -println() -println("="^60) -println("Test 2 — near-zero rotation θ ≈ 0") -println("="^60) -for θ_small in [1e-4, 1e-6, 1e-8, 1e-12] - axis = SVector{3,Float64}(1.0, 0.0, 0.0) - R = make_rotation(axis, θ_small) - v1 = rotation_vector_eigen(R) - v2 = rotation_vector_closed(R) - println("θ=%.0e |v_eigen|=%.6e |v_closed|=%.6e diff=%.2e\n", - θ_small, norm(v1), norm(v2), norm(v1-v2)) -end - -# ── Test 3: near-π — RELEVANT because theta_T=3.0 rad ──────────────────────── -println() -println("="^60) -println("Test 3 — near-π rotation (theta_T=3.0 means you DO reach this)") -println("="^60) -for δ in [1e-1, 1e-3, 1e-5, 1e-7] - axis = normalize(SVector{3,Float64}(1.0, 2.0, 3.0)) - R = make_rotation(axis, π - δ) - v1 = rotation_vector_eigen(R) - v2 = rotation_vector_closed(R) - println("θ=π-%.0e |v_eigen|=%.6f |v_closed|=%.6f diff=%.2e\n", - δ, norm(v1), norm(v2), norm(v1-v2)) -end - -# ── Test 4: benchmark ───────────────────────────────────────────────────────── -println() -println("="^60) -println("Test 4 — benchmark: 1000 calls each") -println("="^60) -Rs = [make_rotation(normalize(SVector{3,Float64}(randn(),randn(),randn())), - 0.1 + (π-0.2)*rand()) - for _ in 1:1000] - -t_eigen = @belapsed sum(norm, rotation_vector_eigen.($(Rs))) -t_closed = @belapsed sum(norm, rotation_vector_closed.($(Rs))) -println("eigen : %8.2f µs\n", t_eigen * 1e6) -println("closed : %8.2f µs\n", t_closed * 1e6) -println("speedup: %.1fx\n", t_eigen / t_closed) - -# ── Test 5: verify the identity R - Rᵀ = 2 sin(θ) N analytically ───────────── -println() -println("="^60) -println("Test 5 — verify R - Rᵀ = 2·sin(θ)·N identity directly") -println("="^60) -axis = normalize(SVector{3,Float64}(1.0, 2.0, 3.0)) -θ_ref = 1.23456 -R = make_rotation(axis, θ_ref) -antisym = (R - R') / 2 # = sin(θ)·N -n_from_antisym = SVector{3,Float64}( - antisym[3,2], # = sin(θ)·n[1] - antisym[1,3], # = sin(θ)·n[2] - antisym[2,1], # = sin(θ)·n[3] -) -n_recovered = n_from_antisym / norm(n_from_antisym) -println("True axis : ", round.(axis, sigdigits=8)) -println("Recovered axis: ", round.(n_recovered, sigdigits=8)) -println("Error on axis : ", norm(axis - n_recovered)) -println("True θ : ", θ_ref) -println("Recovered θ : ", atan(norm(n_from_antisym), (tr(R)-1)/2)) \ No newline at end of file From 85fb2fa8ddf033c3430b01bd6c155b2de468431d Mon Sep 17 00:00:00 2001 From: cyfraysse Date: Tue, 16 Jun 2026 11:46:43 +0200 Subject: [PATCH 4/4] =?UTF-8?q?fixed=20comments=20as=20suggested=20by=20Fr?= =?UTF-8?q?an=C3=A7ois?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/rotation.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/rotation.jl b/src/rotation.jl index 4ae2173..845751e 100644 --- a/src/rotation.jl +++ b/src/rotation.jl @@ -31,7 +31,7 @@ function get_all_body_frames!(R_bodyframe, system::Molecules) L = system.box[1] pos = system.position @inbounds for m in 1:system.Nmol - s = system.start_mol[m] + s = system.start_mol[m] R_bodyframe[m] = body_frame(pos[s], pos[s+1], pos[s+2], L) end end @@ -112,7 +112,7 @@ function Arianna.initialise(algorithm::ComputeRotation, simulation::Simulation) n_θ = length(algorithm.theta_T) state.R_bodyframe = Vector{SMatrix{3,3,T,9}}(undef, N_mol) # allocate once - get_all_body_frames!(state.R_bodyframe, system) # ← state.R_buf + get_all_body_frames!(state.R_bodyframe, system) # fill state.R_bodyframe state.R_ref = [copy(state.R_bodyframe) for _ in 1:n_θ] state.Φ_acc = [[zero(SVector{3,T}) for _ in 1:N_mol] for _ in 1:n_θ] @@ -133,7 +133,7 @@ function Arianna.make_step!(simulation::Simulation, algorithm::ComputeRotation) system = simulation.chains[c] state = algorithm.states[c] N_mol = system.Nmol - get_all_body_frames!(state.R_bodyframe, system) # ← in-place, no alloc # ← alias, no copy + get_all_body_frames!(state.R_bodyframe, system) # no new vector created updates state.R_bodyframe for (k, θ_T) in enumerate(algorithm.theta_T) for m in 1:N_mol dR = state.R_ref[k][m]' * state.R_bodyframe[m]