From faa080479c53f380f199ff623f3f655b2b829071 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Thu, 10 Oct 2024 14:31:17 -0600
Subject: [PATCH] Tests passing

---
 src/Symmetry/Crystal.jl     | 14 +++++++-------
 src/Symmetry/WyckoffData.jl | 18 ++++++++++++++++++
 test/test_symmetry.jl       |  9 +++++----
 3 files changed, 30 insertions(+), 11 deletions(-)

diff --git a/src/Symmetry/Crystal.jl b/src/Symmetry/Crystal.jl
index 723f6febf..76b682281 100644
--- a/src/Symmetry/Crystal.jl
+++ b/src/Symmetry/Crystal.jl
@@ -412,7 +412,8 @@ function crystal_from_symops(latvecs::Mat3, positions::Vector{Vec3}, types::Vect
     all_types = String[]
     classes = Int[]
 
-    # Fill atom positions by applying all symops
+    # Fill atom positions by applying all symops. See crystallographic_orbit
+    # for similar logic.
     for i = eachindex(positions)
         for s = symops
             x = wrap_to_unit_cell(transform(s, positions[i]); symprec)
@@ -455,11 +456,11 @@ end
 
 function get_wyckoff(cryst::Crystal, i::Int)
     (; classes, positions, sg_number, sg_setting, symprec) = cryst
-    c = classes[i]
-    rs = positions[findall(==(c), classes)]
-    rs_std = [transform(sg_setting, r) for r in rs]
-    wyckoff = find_wyckoff_for_orbit(sg_number, rs_std; symprec)
-    @assert wyckoff.multiplicity ≈ length(rs) / det(sg_setting.R)
+    isnothing(sg_number) || isnothing(sg_setting) && return nothing
+    r = transform(sg_setting, positions[i])
+    wyckoff = find_wyckoff_for_position(sg_number, r; symprec)
+    cell_multiplicity = count(==(classes[i]), classes)
+    @assert wyckoff.multiplicity ≈ cell_multiplicity / abs(det(sg_setting.R))
     return wyckoff
 end
 
@@ -777,7 +778,6 @@ function hyperkagome_crystal(; a=1.0)
     # https://materials.springer.com/isp/crystallographic/docs/sd_1723194
     x = 0.141
     cryst = Crystal(latvecs, [[1/8, x, x+1/4]], 213)
-    @assert get_wyckoff(cryst, 1).letter = 'Z'
     return cryst
 end
 
diff --git a/src/Symmetry/WyckoffData.jl b/src/Symmetry/WyckoffData.jl
index f61379c72..a0a6f077f 100644
--- a/src/Symmetry/WyckoffData.jl
+++ b/src/Symmetry/WyckoffData.jl
@@ -4,6 +4,24 @@ struct Wyckoff
     sitesym      :: String
 end
 
+function crystallographic_orbit(symops, r; symprec)
+    orbit = Vec3[]
+    for s in symops
+        r′ = transform(s, r)
+        if !any(is_periodic_copy.(orbit, Ref(r′); symprec))
+            push!(orbit, wrap_to_unit_cell(r′; symprec))
+        end
+    end
+    return orbit
+end
+
+function find_wyckoff_for_position(sg_number, r; symprec)
+    Rs, Ts = Spglib.get_symmetry_from_database(standard_setting[sg_number])
+    symops = SymOp.(Rs, Ts)
+    rs = crystallographic_orbit(symops, r; symprec)
+    return find_wyckoff_for_orbit(sg_number, rs; symprec)
+end
+
 function find_wyckoff_for_orbit(sg_number, rs; symprec)
     for w in reverse(wyckoff_table[sg_number])
         (mult, letter, sitesym, pos) = w
diff --git a/test/test_symmetry.jl b/test/test_symmetry.jl
index 59ca013dc..2efb7ea86 100644
--- a/test/test_symmetry.jl
+++ b/test/test_symmetry.jl
@@ -178,7 +178,7 @@ end
     # Inference of Wyckoff symbols
     lat_vecs = lattice_vectors(1, 1, 1.2, 90, 90, 120)
     cryst = Crystal(lat_vecs, [[0.2, 0.2, 1/2]], 164)
-    @test Sunny.get_wyckoff(cryst, 1) == Wyckoff(6, 'h', ".2.")
+    @test Sunny.get_wyckoff(cryst, 1) == Sunny.Wyckoff(6, 'h', ".2.")
 
     ### Check settings for monoclinic spacegroup
 
@@ -186,13 +186,13 @@ end
     latvecs = lattice_vectors(1, 1.1, 1.2, 90, 100, 90)
     cryst = Crystal(latvecs, [[0, 0.2, 1/2]], "C 1 2 1")
     @test cryst.sg_label == "'C 2 = C 1 2 1' (5)"
-    @test Sunny.get_wyckoff(cryst, 1).letter == 'Z'
+    @test Sunny.get_wyckoff(cryst, 1) == Sunny.Wyckoff(2, 'b', "2")
 
     # Alternative setting
     latvecs2 = reduce(hcat, eachcol(latvecs)[[3, 1, 2]])
     cryst2 = Crystal(latvecs2, [[1/2, 0, 0.2]], "A 1 1 2")
     @test cryst2.sg_label == "'C 2 = A 1 1 2' (5)"
-    @test Sunny.get_wyckoff(cryst, 1).letter == 'Z'
+    @test Sunny.get_wyckoff(cryst, 1) == Sunny.Wyckoff(2, 'b', "2")
 
     # Verify `cryst` is already in standard setting
     @test cryst.sg_setting.R ≈ I
@@ -390,8 +390,9 @@ end
         Modified reference frame! Transform using `rotate_operator(op, R)`.
         """
 
+    cryst = Sunny.hyperkagome_crystal()
+    @assert Sunny.get_wyckoff(cryst, 1) == Sunny.Wyckoff(12, 'd', "..2")
     capt = IOCapture.capture() do
-        cryst = Sunny.hyperkagome_crystal()
         print_suggested_frame(cryst, 2)
     end
     @test capt.output == """