Skip to content

Commit d63bb09

Browse files
Alexsp32Alexsp32
committed
Add surface-site-classification to branch (blocking PR for this PR: #437)
1 parent d46e6b2 commit d63bb09

File tree

5 files changed

+299
-0
lines changed

5 files changed

+299
-0
lines changed

‎src/Analysis/Analysis.jl‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,10 @@ using Reexport: @reexport
88
include("diatomic.jl")
99
export Diatomic
1010

11+
# Surface facet symmetry site classification
12+
include("high-symmetry_sites.jl")
13+
export HighSymmetrySites
14+
1115
# Allow re-processing of simulation data.
1216
include("postprocess.jl")
1317
export Postprocess

‎src/Analysis/high-symmetry_sites.jl‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
"""
2+
These analysis functions contain some postprocessing code to associate positions of adsorbate molecules on a surface slab with high-symmetry locations on the surface slab.
3+
4+
Base definitions for fcc metal surface facets are included, but there are options to expand this further.
5+
6+
**This code only works for 2D positions in X and Y, so the surface needs to lie in the XY plane.**
7+
8+
"""
9+
module HighSymmetrySites
10+
11+
12+
using NQCDynamics: AbstractSimulation, PeriodicCell, get_positions
13+
14+
export SlabStructure, FCC100Sites, FCC110Sites, FCC111Sites, FCC211Sites, positions_to_category, classify_every_frame
15+
16+
struct SlabStructure
17+
adsorbate_indices::Vector{Int}
18+
symmetry_sites::Dict{Symbol, Vector{Vector{Float64}}}
19+
supercell_size::Vector{Float64}
20+
end
21+
22+
const FCC100Sites = Dict(
23+
:top => vec([vcat(i...) for i in Iterators.product([0.0,1.0], [0.0,1.0])]),
24+
:bridge => [
25+
[0.5, 0.0],
26+
[0.5, 1.0],
27+
[0.0, 0.5],
28+
[1.0, 0.5],
29+
],
30+
:fcc => [
31+
[0.5, 0.5],
32+
],
33+
)
34+
const FCC110Sites = Dict(
35+
:top => vec([vcat(i...) for i in Iterators.product([0.0,1.0], [0.0,1.0])]),
36+
:long_bridge => [
37+
[0.5, 0.0],
38+
[0.5, 1.0],
39+
],
40+
:short_bridge => [
41+
[0.0, 0.5],
42+
[1.0, 0.5],
43+
],
44+
:center_hollow => [
45+
[0.5, 0.5],
46+
],
47+
:step_hollow => [
48+
[0.25, 0.5],
49+
[0.75, 0.5],
50+
],
51+
)
52+
const FCC111Sites = Dict(
53+
:top => vec([vcat(i...) for i in Iterators.product([0.0,1.0], [0.0,1.0])]),
54+
:fcc => [[0.33,0.33]],
55+
:hcp => [[0.66,0.66]],
56+
:bridge => [
57+
[0.5,0.0],
58+
[0.5,0.5],
59+
[0.0,0.5],
60+
[1.0,0.5],
61+
[0.5,1.0],
62+
],
63+
)
64+
const FCC211Sites = Dict( # Based on Cao2018 site definitions
65+
:step_edge => vec([vcat(i...) for i in Iterators.product([0.0,1.0], [0.0,1.0])]), # Top sites, marking desorption bridged over the step edge
66+
:short_step => [ # This is the higher hollow site close to the step edge and would count for desorption via the step edge. Site 1 in Cao2018
67+
[0.125, 0.5, ] #
68+
],
69+
:fcc_high => [ # Site 2 in Cao2018, equivalent to desorption not over the step edge
70+
[0.5, 0, ],
71+
[0.5, 1.0, ],
72+
],
73+
:fcc_low => [ # This is hollow site 3 in Cao2018, equivalent to desorption not over the step edge
74+
[0.75,0.0, ],
75+
[0.75,1.0, ],
76+
],
77+
:long_step => [ # This is hollow site 4, the lower one with the higher step edge, counting for desorption over the higher step.
78+
[0.875, 0.5, ],
79+
], #Not really sure where b2 goes.
80+
)
81+
82+
"""
83+
positions_to_category(position, categories, cell::PeriodicCell; fractional::Bool = false, snap_to_site::Float64 = 0.03)
84+
85+
Classifies a given 2D position according to its proximity to high-symmetry surface sites.
86+
87+
# Arguments
88+
- `position`: A 2-element (or longer, only first 2D used) vector representing the coordinates to classify.
89+
- `categories`: A dictionary mapping site names (as Symbols) to lists of site coordinates.
90+
- `cell::PeriodicCell`: The surface unit cell used to calculate distances and perform coordinate transformations.
91+
- `fractional::Bool = false`: Whether to use fractional coordinates to determine distances to sites. If `true`, `position` is mapped directly to the sites; if `false`, converts to fractional using `cell`.
92+
- `snap_to_site::Float64 = 0.03`: The distance tolerance (in the same units as position/cell) for snapping to a given site category.
93+
94+
# Returns
95+
- The key of the closest site category (from `categories`) if within `snap_to_site` of any defined site. If not within this range or outside the cell, returns `:other`.
96+
97+
# Notes
98+
- Positions are always truncated to 2D before classification.
99+
- If the position is not inside the periodic cell, returns `:other`.
100+
- Site assignment is based on the minimum Euclidean distance to any site in each category.
101+
"""
102+
function positions_to_category(
103+
position,
104+
categories,
105+
cell::PeriodicCell;
106+
fractional::Bool=false,
107+
snap_to_site::Float64=0.03,
108+
)
109+
position = first(position,2) # Ensure we only ever carry X and Y
110+
category_names = collect(keys(categories))
111+
if !check_atoms_in_cell(cell, hcat(position))
112+
return :other
113+
else
114+
position_copy = position
115+
# Fractionalise positions if requested.
116+
fractional ? position_copy = cell.inverse * position_copy : nothing
117+
site_distances = Float64[]
118+
for site in category_names
119+
if fractional
120+
push!(site_distances, minimum([norm(position_copy - point) for point in categories[site]])) # Gets closest distance for each category
121+
else
122+
push!(site_distances, minimum([norm(position_copy - cell.vectors * point) for point in categories[site]])) # Gets closest distance for each category
123+
end
124+
end
125+
end
126+
if any(site_distances .≤ snap_to_site) # Check if any category is sufficiently close
127+
return category_names[argmin(site_distances)] # Return the classified category
128+
else
129+
return :other
130+
end
131+
end
132+
133+
"""
134+
classify_every_frame(
135+
trajectory,
136+
cell::PeriodicCell,
137+
slab::SlabStructure,
138+
fractional::Bool = false,
139+
snap_to_site::Float64 = 0.03,
140+
)
141+
142+
Classifies the adsorption site of adsorbate atoms for every frame in a trajectory.
143+
144+
# Arguments
145+
- `trajectory`: Vector{DynamicsVariables type} of structures to classify.
146+
- `cell::PeriodicCell`: The periodic cell representation of the simulation.
147+
- `slab::SlabStructure`: Slab structure object containing adsorbate atom indices, supercell size, and symmetry site information.
148+
- `fractional::Bool = false`: If `true`, uses fractional coordinates for site matching; otherwise, uses Cartesian coordinates.
149+
- `snap_to_site::Float64 = 0.03`: Distance tolerance for assigning a position to a high-symmetry site.
150+
151+
# Returns
152+
- A vector of vectors, where each inner vector contains the classified site (as a `Symbol`) for an adsorbate atom at one frame.
153+
154+
# Notes
155+
- Each adsorbate atom is assigned to a symmetry site or `:other` for each trajectory frame.
156+
- Positions are wrapped back into the primitive cell before site classification.
157+
- Uses `positions_to_category` to perform the categorization for each atom in each frame.
158+
"""
159+
function classify_every_frame(
160+
trajectory,
161+
cell::PeriodicCell,
162+
slab::SlabStructure;
163+
fractional::Bool = false,
164+
snap_to_site::Float64 = 0.03,
165+
)
166+
# Determine a 2D representation of the unit cell first
167+
fullsize_cell = cell.vectors
168+
primitive_cell = PeriodicCell(fullsize_cell ./ slab.supercell_size)
169+
primitive_cell_2d = PeriodicCell(getindex(fullsize_cell ./ slab.supercell_size, 1:2, 1:2))
170+
171+
site_associations = [Symbol[] for _ in eachindex(slab.adsorbate_indices)] # Generate site association vector containing info for each adsorbate atom.
172+
@inbounds for trajectory_frame in eachindex(trajectory)
173+
@inbounds for (output_list_index, atom_index) in enumerate(slab.adsorbate_indices)
174+
position_H = hcat(
175+
get_positions(
176+
trajectory[trajectory_frame]
177+
)[:, atom_index]
178+
)
179+
# Move into primitive unit cell
180+
apply_cell_boundaries!(primitive_cell, position_H)
181+
position_H = vec(position_H)
182+
# Then categorise
183+
category = positions_to_category(
184+
position_H,
185+
slab.symmetry_sites,
186+
primitive_cell_2d,
187+
fractional=fractional,
188+
snap_to_site=snap_to_site
189+
)
190+
push!(site_associations[output_list_index], category)
191+
end
192+
end
193+
return site_associations
194+
end
195+
196+
197+
end

‎src/DynamicsOutputs.jl‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -547,6 +547,30 @@ end
547547

548548
export OutputNoise
549549

550+
struct OutputSurfaceSiteClassification
551+
slab_info::Analysis.HighSymmetrySites.SlabStructure
552+
fractional::Bool
553+
snap_to_site::Float64
554+
end
555+
function OutputSurfaceSiteClassification(slab_info::Analysis.HighSymmetrySites.SlabStructure; fractional=false, snap_to_sites=0.03)
556+
return OutputSurfaceSiteClassification(
557+
slab_info,
558+
fractional,
559+
snap_to_sites
560+
)
561+
end
562+
563+
function (output::OutputSurfaceSiteClassification)(sol, i)
564+
cell = sol.prob.p.cell # Simulation needs to have an associated PeriodicCell
565+
return Analysis.HighSymmetrySites.classify_every_frame(
566+
sol.u,
567+
cell,
568+
output.slab_info,
569+
output.fractional,
570+
output.snap_to_site,
571+
)
572+
end
573+
550574
"""
551575
OutputEverything(sol,i)
552576

‎test/Analysis/high_symmetry.jl‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
using Test
2+
using NQCDynamics
3+
using LinearAlgebra
4+
5+
@testset "Surface site classification" begin
6+
test_cell = PeriodicCell(I(3) .* 3) # Unit length cell
7+
dummy_sites = Dict(
8+
:type_1 => [
9+
[0.0,0.0],
10+
],
11+
:type_2 => [
12+
[1.0,1.0],
13+
],
14+
)
15+
testing_positions = [rand(3,10) for _ in 1:4]
16+
# Indices 9 and 10 are the atom to test
17+
testing_positions[1][:, 9] .= [0,0,0]
18+
testing_positions[1][:, 10] .= [0,0,0]
19+
testing_positions[2][:, 9] .= [1,1,0]
20+
testing_positions[2][:, 10] .= [1,1,0]
21+
testing_positions[3][:, 9] .= [0.5,0,0]
22+
testing_positions[3][:, 10] .= [0.5,0,0]
23+
testing_positions[4][:, 9] .= [0.25,0,0]
24+
testing_positions[4][:, 10] .= [0.25,0,0]
25+
# Case 1: 0,0 should snap to site type 1
26+
result_case1 = Analysis.HighSymmetrySites.positions_to_category(
27+
testing_positions[1][:, 9],
28+
dummy_sites,
29+
PeriodicCell(I(2))
30+
)
31+
@test result_case1 == :type_1
32+
# Case 2: 1,1 should snap to site type 2
33+
result_case2 = Analysis.HighSymmetrySites.positions_to_category(
34+
testing_positions[2][:, 9],
35+
dummy_sites,
36+
PeriodicCell(I(2))
37+
)
38+
@test result_case2 == :type_2
39+
# Case 3: 0.5,0 should map to :other if snap_to_site isn't particularly large.
40+
result_case3 = Analysis.HighSymmetrySites.positions_to_category(
41+
testing_positions[3][:, 9],
42+
dummy_sites,
43+
PeriodicCell(I(2))
44+
)
45+
@test result_case3== :other
46+
# Case 4: 0.25,0 should map to site type 1 if snap_to_site is set to 0.5.
47+
result_case4 = Analysis.HighSymmetrySites.positions_to_category(
48+
testing_positions[2][:, 9],
49+
dummy_sites,
50+
PeriodicCell(I(2)),
51+
snap_to_site = 0.5,
52+
)
53+
@test result_case4 == :type_1
54+
# Case 5: Analysing the full trajectory should give the correct results.
55+
result_case5 = Analysis.HighSymmetrySites.classify_every_frame(
56+
[DynamicsVariables(v = rand(3,10), r = i) for i in testing_positions],
57+
test_cell,
58+
Analysis.HighSymmetrySites.SlabStructure(
59+
[9,10],
60+
dummy_sites,
61+
[3.0,3.0,1.0],
62+
),
63+
snap_to_site = 0.45,
64+
)
65+
for i in 1:2
66+
@test result_case5[1][i] == :site_1
67+
@test result_case5[2][i] == :site_2
68+
@test result_case5[3][i] == :other
69+
@test result_case5[4][i] == :site_1
70+
end
71+
end

‎test/runtests.jl‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,9 @@ if GROUP == "All" || GROUP == "Analysis"
3030
@safetestset "Diatomic Analysis" begin
3131
include("Analysis/diatomic.jl")
3232
end
33+
@safetestset "Surface site selectivity" begin
34+
include("Analysis/high_symmetry.jl")
35+
end
3336
end
3437

3538
if GROUP == "All" || GROUP == "InitialConditions"

0 commit comments

Comments
 (0)