on-the-fly rotation tracker following ParticlesMC architecture

Project.toml

@@ -15,6 +15,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Transducers = "28d57a85-8fef-5791-bfe6-a80928e7c999"
 
 [compat]
 Aqua = "0.8"
@@ -26,19 +27,20 @@ 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"
+Transducers = "0.4.85"
 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"]

src/ParticlesMC.jl

@@ -6,7 +6,7 @@ Exports commonly-used types (e.g., `Particles`, `Model`) and helper functions fo
 """
 module ParticlesMC
 
-using Arianna, StaticArrays
+using Arianna, StaticArrays, Transducers
 using Comonicon, TOML
 using Comonicon: @main
 
@@ -20,6 +20,7 @@ include("models.jl")
 include("molecules.jl")
 include("atoms.jl")
 include("moves.jl")
+include("rotation.jl")
 
 """Return the position of particle `i` in `system`.
 
@@ -125,6 +126,7 @@ export perform_action!, revert_action!
 include("IO/IO.jl")
 using .IO: XYZ, EXYZ, LAMMPS, load_configuration, load_chains
 export XYZ, EXYZ, LAMMPS, load_configuration, load_chains
+export ComputeRotation, StorePhiTrajectories, StoreLastPhiFrame
 
 
 """
@@ -250,6 +252,31 @@ ParticlesMC implemented in Comonicon.
     end
     push!(algorithm_list, (algorithm=Metropolis, pool=pool, seed=seed, parallel=parallel, sweepstep=length(chains[1])))
 
+    # Setup observables
+    for observable in get(sim, "observable", [])
+        alg = observable["algorithm"]
+        scheduler_params = observable["scheduler_params"]
+        interval = get(scheduler_params, "linear_interval", 1)
+        if "log_base" in keys(scheduler_params)
+            block = build_schedule(interval, 0, 2.0)
+            sched = build_schedule(steps, burn, block)
+        else
+            sched = build_schedule(steps, burn, interval)
+        end
+        if alg == "ComputeRotation"
+            parameters = get(observable, "parameters", Dict())
+            theta_T    = Float64.(get(parameters, "theta_T", [π/4]))
+            algorithm  = (
+                algorithm=ComputeRotation,
+                scheduler=sched,
+                theta_T=theta_T,
+            )
+        else
+            error("Unsupported observable algorithm: $alg")
+        end
+        push!(algorithm_list, algorithm)
+    end
+
     # Setup outputs
     for output in sim["output"]
         alg = output["algorithm"]
@@ -284,6 +311,12 @@ ParticlesMC implemented in Comonicon.
                 scheduler=sched,
                 fmt=eval(Meta.parse("$(fmt)()")),
             )
+        elseif alg == "StorePhiTrajectories" || alg == "StoreLastPhiFrame"
+            algorithm = (
+                algorithm=eval(Meta.parse(alg)),
+                scheduler=sched,
+                path=output_path,
+            )
         elseif alg == "PrintTimeSteps"
             algorithm = (
                 algorithm=eval(Meta.parse(alg)),

src/molecules.jl

@@ -9,6 +9,7 @@ models, and the neighbour list structure used for energy and neighbour queries.
 
 Fields
 - `position::Vector{SVector{D,T}}`: positions of particles.
+- `Φ::Vector{Vector{SVector{3,T}}}`: rotation vector of molecules per threshold index k.
 - `species::VS`: species/type of each particle.
 - `molecule::VS`: molecule identifier for each particle.
 - `molecule_species::VS`: species identifier for each molecule.
@@ -23,6 +24,7 @@ Fields
 """
 struct Molecules{D,  VS<:AbstractVector, C<:NeighbourList, T<:AbstractFloat, SM<:AbstractArray} <: Particles
     position::Vector{SVector{D,T}}
+    Φ::Vector{Vector{SVector{3,T}}}   # Φ[k][m] = rotation vector for molecule m, threshold k
     species::VS
     molecule::VS
     molecule_species::VS
@@ -76,6 +78,7 @@ Returns
 function System(position, species, molecule, density::T, temperature::T, model_matrix, bonds; molecule_species=nothing, list_type=EmptyList, list_parameters=nothing) where {T<:AbstractFloat}
     @assert length(position) == length(species)
     N = length(position)
+    Φ = Vector{Vector{SVector{3,T}}}()   # empty, filled by ComputeRotation
     Nmol = length(unique(molecule))
     start_mol, length_mol = get_first_and_counts(molecule)
     molecule_species = something(molecule_species, ones(Int, N))
@@ -84,7 +87,7 @@ function System(position, species, molecule, density::T, temperature::T, model_m
     energy = zeros(T, 1)
     maxcut = maximum([model.rcut for model in model_matrix])
     neighbour_list = list_type(box, maxcut, N; list_parameters=list_parameters)
-    system = Molecules(position, species, molecule, molecule_species,  start_mol, length_mol, density, temperature, energy, model_matrix, d, N, Nmol,box, neighbour_list, bonds)
+    system = Molecules(position, Φ, species, molecule, molecule_species,  start_mol, length_mol, density, temperature, energy, model_matrix, d, N, Nmol,box, neighbour_list, bonds)
     build_neighbour_list!(system)
     local_energy = [compute_energy_particle(system, i, neighbour_list) for i in eachindex(position)]
     energy = sum(local_energy) / 2