And now the unit test works

Project.toml

@@ -59,7 +59,6 @@ FastLapackInterface = "1, 2"
 HDF5 = "0.15, 0.16, 0.17"
 Interpolations = "0.13, 0.14"
 MuladdMacro = "0.2"
-MKL_jll = "2023" # Pin to before v2024.X to ensure tests build successfully. See #324
 NQCBase = "0.2"
 NQCDistributions = "0.1"
 NQCModels = "0.8"
@@ -93,6 +92,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JuLIP = "945c410c-986d-556a-acb1-167a618e0462"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
@@ -104,4 +104,4 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SafeTestsets", "CSV", "DataFrames", "DiffEqNoiseProcess", "FiniteDiff", "DiffEqDevTools", "Logging", "MKL", "PyCall", "JuLIP", "Symbolics", "Statistics", "Plots"]
+test = ["Test", "SafeTestsets", "CSV", "DataFrames", "DiffEqNoiseProcess", "FiniteDiff", "DiffEqDevTools", "Logging", "MKL", "PyCall", "JuLIP", "Symbolics", "Statistics", "Plots", "JLD2"]

src/Analysis/diatomic.jl

@@ -6,17 +6,18 @@ using NQCDynamics: AbstractSimulation, Simulation, get_positions, Structure
 using NQCBase
 using Unitful, UnitfulAtomic
 using LinearAlgebra
+using Statistics
 
 # Surface normal projection
 surface_normal_height(x::AbstractVector, surface_normal::AbstractVector=[0,0,1])=norm(x.*normalize(surface_normal))
 
 """
-    surface_distance_condition(x::AbstractArray, indices::Vector{Int}, simulation::NQCDynamics.AbstractSimulation; surface_distance_threshold=5.0*u"Å")
+    surface_distance_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; surface_distance_threshold=5.0*u"Å")
 
     Checks that the diatomic molecule is at least `surface_distance_threshold` away from the highest substrate atom in the simulation.
 """
-function surface_distance_condition(x::AbstractArray, indices::Vector{Int}, simulation::NQCDynamics.AbstractSimulation; surface_distance_threshold=5.0*u"Å")
-        molecule_position=surface_normal_height(pbc_center_of_mass(x, indices..., simulation)) # molecule position wrt surface normal
+function surface_distance_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; surface_distance_threshold=5.0*u"Å")
+        molecule_position=surface_normal_height(Structure.pbc_center_of_mass(x, indices..., simulation)) # molecule position wrt surface normal
         substrate_heights=surface_normal_height.(eachcol(get_positions(x)[:, symdiff(1:length(simulation.atoms.masses), indices)])) # substrate positions along surface normal
         highest_z=max(substrate_heights[substrate_heights.≤molecule_position]...) # Ignore substrate above molecule. 
         #? surface_distance_threshold must be < vacuum above surface and > distance between substrate layers. 
@@ -29,12 +30,12 @@ function surface_distance_condition(x::AbstractArray, indices::Vector{Int}, simu
 end
 
 """
-    com_velocity_condition(x::AbstractArray, indices::Vector{Int}, simulation::NQCDynamics.AbstractSimulation; surface_normal::AbstractVector=[0,0,1])
+    com_velocity_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; surface_normal::AbstractVector=[0,0,1])
 
     Evaluates true if the centre of mass velocity vector of the diatomic molecule points to the surface.
 """
-function com_velocity_condition(x::AbstractArray, indices::Vector{Int}, simulation::NQCDynamics.AbstractSimulation; surface_normal::AbstractVector=[0,0,1])
-        if dot(velocity_center_of_mass(x, indices..., simulation), normalize(surface_normal))>0
+function com_velocity_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; surface_normal::AbstractVector=[0,0,1])
+        if dot(Structure.velocity_center_of_mass(x, indices..., simulation), normalize(surface_normal))>0
             return false
         else
             return true
@@ -53,14 +54,14 @@ This is evaluated using two conditions:
 2. Desorption begins at the turning point of the centre of mass velocity component along `surface_normal`, indicating overall movement away from the surface. 
 """
 function get_desorption_frame(trajectory::AbstractVector, diatomic_indices::Vector{Int}, simulation::AbstractSimulation;surface_normal::Vector = [0,0,1], surface_distance_threshold = 5.0*u"Å")
-    desorbed_frame=findfirst(surface_distance_condition.(trajectory, Ref(indices), Ref(simulation); surface_distance_threshold=surface_distance_threshold))
+    desorbed_frame=findfirst(surface_distance_condition.(trajectory, Ref(diatomic_indices), Ref(simulation); surface_distance_threshold=surface_distance_threshold))
 
     if isa(desorbed_frame, Nothing)
         @debug "No desorption event found."
         return nothing
     else
         @debug "Desorption observed in snapshot $(desorbed_frame)"
-        leaving_surface_frame=findlast(com_velocity_condition.(view(trajectory, 1:desorbed_frame), Ref(indices), Ref(simulation); surface_normal=surface_normal)) #ToDo testing views for memory efficiency, need to test time penalty. Also need to test if running on everything and findfirst-ing the Bool array is quicker. 
+        leaving_surface_frame=findlast(com_velocity_condition.(view(trajectory, 1:desorbed_frame), Ref(diatomic_indices), Ref(simulation); surface_normal=surface_normal)) #ToDo testing views for memory efficiency, need to test time penalty. Also need to test if running on everything and findfirst-ing the Bool array is quicker. 
         if leaving_surface_frame<0
             @warn "Desorption event detected, but no direction change detected."
             return nothing
@@ -81,7 +82,7 @@ function get_desorption_angle(trajectory::AbstractVector, indices::Vector{Int},
     # Determine the average centre of mass velocity to decrease error due to vibration and rotation orthogonal to true translational component. 
     com_velocities = zeros(eltype(trajectory[1]), length(surface_normal), length(trajectory)-desorption_frame)
     for i in 1:length(trajectory)-desorption_frame
-        com_velocities[:, i] .= velocity_center_of_mass(trajectory[i+desorption_frame], indices[1], indices[2], simulation)
+        com_velocities[:, i] .= Structure.velocity_center_of_mass(trajectory[i+desorption_frame], indices[1], indices[2], simulation)
     end
     average_velocity=mean(com_velocities;dims=2)
     # Now convert into an angle by arccos((a•b)/(|a|*|b|))

test/Analysis/diatomic.jl

@@ -2,20 +2,17 @@ using Test
 using NQCDynamics
 using Unitful, UnitfulAtomic
 using PyCall
-using NQCDynamics, NQCModels
+using NQCModels
 using JLD2
 
 @testset "Desorption Analysis" begin
     ase_io=pyimport("ase.io")
     desorption_trajectory_ase=ase_io.read("artifacts/desorption_test.xyz", index=":")
     desorption_dynamicsvariables=jldopen("artifacts/desorption_dynamics.jld2")["trajectory"]
     # Make a dummmy simulation out of the first structure (we don't need a working potential to test. )
-    atoms, initial_positions, cell=NQCDynamics.NQCBase.convert_from_ase_atoms(desorption_trajectory_ase[1])
+    atoms, initial_positions, cell=NQCDynamics.convert_from_ase_atoms(desorption_trajectory_ase[1])
     simulation=Simulation(atoms, AdiabaticASEModel(desorption_trajectory_ase[1]), cell=cell)
     diatomic_indices=[55,56]
-    @test @time "DesorptionFrame - No views" get_desorption_frame_noviews(desorption_dynamicsvariables, diatomic_indices, simulation; surface_distance_threshold=2.4u"Å") == 2675
-    @test @time "DesorptionFrame - With views" get_desorption_frame_views(desorption_dynamicsvariables, diatomic_indices, simulation; surface_distance_threshold=2.4u"Å") == 2675
-    @test @time "DesorptionAngle - No views" get_desorption_angle(desorption_dynamicsvariables, diatomic_indices, simulation; surface_distance_threshold=2.4u"Å", evaluation_function=get_desorption_frame_noviews) ≈ 28.05088202518
-    @test @time "DesorptionAngle - With views" get_desorption_angle(desorption_dynamicsvariables, diatomic_indices, simulation; surface_distance_threshold=2.4u"Å", evaluation_function=get_desorption_frame_views) ≈ 28.05088202518
-end
-
+    @test Analysis.Diatomic.get_desorption_frame(desorption_dynamicsvariables, diatomic_indices, simulation; surface_distance_threshold=2.4u"Å") == 2675
+    @test Analysis.Diatomic.get_desorption_angle(desorption_dynamicsvariables, diatomic_indices, simulation; surface_distance_threshold=2.4u"Å") ≈ 28.05088202518
+end

test/runtests.jl

@@ -17,6 +17,10 @@ if GROUP == "All" || GROUP == "Core"
     @safetestset "Ring Polymer Tests" begin
         include("Core/ring_polymers.jl")
     end
+
+end
+
+if GROUP == "All" || GROUP == "Analysis"
     @safetestset "Estimator tests" begin
         include("Core/estimators.jl")
     end