From a4b346c0f23bc7fd546c8e7c69a78a7d586f48c1 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Sat, 2 Nov 2024 18:41:01 -0600
Subject: [PATCH 1/2] Add Wyckoff inference test

---
 src/Symmetry/Crystal.jl      |  6 ++----
 src/Symmetry/LatticeUtils.jl | 23 ++++++++++-------------
 src/Symmetry/WyckoffData.jl  | 17 +++++++++--------
 test/test_symmetry.jl        | 14 +++++++++++++-
 4 files changed, 34 insertions(+), 26 deletions(-)

diff --git a/src/Symmetry/Crystal.jl b/src/Symmetry/Crystal.jl
index 88a8e972f..e21cdf310 100644
--- a/src/Symmetry/Crystal.jl
+++ b/src/Symmetry/Crystal.jl
@@ -397,11 +397,9 @@ function crystal_from_spacegroup(latvecs::Mat3, positions::Vector{Vec3}, types::
     wyckoffs = find_wyckoff_for_position.(Ref(sg), positions; symprec)
     orbits = crystallographic_orbits_distinct(sg.symops, positions; symprec, wyckoffs)
 
-    # Check that inferred orbits match known multiplicities of the Wyckoffs
+    # Inferred orbits must match known multiplicities of the Wyckoffs
     foreach(orbits, wyckoffs) do orbit, wyckoff
-        if wyckoff.multiplicity ≉ length(orbit) / abs(det(sg.setting.R))
-            error("Internal error filling Wyckoff positions! Please report this.")
-        end
+        @assert wyckoff.multiplicity ≈ length(orbit) / abs(det(sg.setting.R))
     end
 
     all_positions = reduce(vcat, orbits)
diff --git a/src/Symmetry/LatticeUtils.jl b/src/Symmetry/LatticeUtils.jl
index 9b6477b5d..b7f5fa9d1 100644
--- a/src/Symmetry/LatticeUtils.jl
+++ b/src/Symmetry/LatticeUtils.jl
@@ -61,27 +61,25 @@ function is_standard_form(latvecs::Mat3)
     return latvecs ≈ conventional_latvecs
 end
 
-"""
-    CellType
-
-An enumeration over the different types of 3D Bravais unit cells.
-"""
+# Labels for the 7 lattice systems. Note that these are subtly distinct from the
+# 7 crystal systems. In particular, the rhombohedral lattice system (a=b=c,
+# α=β=γ) should not be confused with the trigonal crystal system. Trigonal
+# spacegroups are characterized by a 3-fold rotational symmetry. All trigonal
+# spacegroups (143-167) admit a hexagonal setting. Some of these (146, 148, 155,
+# 160, 161, 166, 167) additionally admit a rhombohedral setting.
 @enum CellType begin
     triclinic
     monoclinic
     orthorhombic
     tetragonal
-    # Rhombohedral is a special case. It is a lattice type (a=b=c, α=β=γ) but
-    # not a spacegroup type. Trigonal space groups are conventionally described
-    # using either hexagonal or rhombohedral lattices.
     rhombohedral
     hexagonal
     cubic
 end
 
-# Infer the `CellType` of a unit cell from its lattice vectors, i.e. the columns
-# of `latvecs`. Report an error if the unit cell is not in conventional form,
-# which would invalidate the table of symops for a given Hall number.
+# Infer the CellType (lattice system) from lattice vectors. Report an error if
+# the unit cell is not in conventional form, which would invalidate the table of
+# symops for a given Hall number.
 function cell_type(latvecs::Mat3)
     a, b, c, α, β, γ = lattice_params(latvecs)
 
@@ -121,8 +119,7 @@ function cell_type(latvecs::Mat3)
     return triclinic
 end
 
-# Return the standard cell convention for a given Hall number using the
-# convention of spglib, listed at
+# The lattice system that is expected for a given Hall number. Taken from:
 # http://pmsl.planet.sci.kobe-u.ac.jp/~seto/?page_id=37
 function cell_type(hall_number::Int)
     if 1 <= hall_number <= 2
diff --git a/src/Symmetry/WyckoffData.jl b/src/Symmetry/WyckoffData.jl
index 0d6a3bff2..9ca84dbcd 100644
--- a/src/Symmetry/WyckoffData.jl
+++ b/src/Symmetry/WyckoffData.jl
@@ -12,8 +12,9 @@ struct WyckoffExpr
 end
 
 # Extracts (F, c) from expressions of the form (F θ + c), where θ = [x, y, z].
-function WyckoffExpr(str)
+function WyckoffExpr(str::String)
     s = SymOp(strip(str, ['(', ')']))
+    @assert all(in((-2, -1, 0, 1, 2)), s.R)
     return WyckoffExpr(s.R, s.T)
 end
 
@@ -48,17 +49,17 @@ end
 # return nothing if this equation cannot be satisfied. Write F θ = b + Δb,
 # involving b = r - c, and arbitrary integers Δb. Given candidate Δb, one can
 # use the pseudo-inverse of F to find least squares solution θ, and then check
-# whether it satisfies the original linear system of equations. A search over
-# integer shifts Δb seems unavoidable.
+# whether it satisfies the original linear system of equations.
 function position_to_wyckoff_params(r::Vec3, w::WyckoffExpr; symprec=1e-8)
     (; F, c) = w
     # Wrapping components to [0, 1] simplifies the search in the next step.
     b = wrap_to_unit_cell(r - c; symprec)
-    # Loop over 8 vector shifts Δb ∈ [{-1,0}, {-1,0}, {-1,0}]. In the general
-    # case, one might need to search over all Δb ∈ ℕ³. This reduced search space
-    # has been empirically determined through consideration of all Wyckoffs of
-    # the 230 spacegroups. Similar logic appears in Crystalline.jl:
-    # https://github.com/thchr/Crystalline.jl/blob/31fa2f61a5db62f4dfd6b0e8fb994d8aacef3d5c/src/wyckoff.jl#L321-L335
+    # If F were arbitrary, one would need to search over all offsets Δb ∈ ℕ³.
+    # Fortunately, the actual matrices are constrained: matrix elements Fᵢⱼ are
+    # in {0, ±1, ±2} for hexagonal crystal systems, and in {0, ±1} for all other
+    # systems. Empirical tests indicate that just eight shifts are sufficient:
+    # Δb ∈ [{-1,0}, {-1,0}, {-1,0}]. Removing these shifts will cause failures
+    # in the "Wyckoff table" unit test.
     for Δb in Iterators.product((-1, 0), (-1, 0), (-1, 0))
         Δb = Vec3(Δb)
         θ = pinv(F) * (b + Δb)
diff --git a/test/test_symmetry.jl b/test/test_symmetry.jl
index 047f15a41..21d50e4aa 100644
--- a/test/test_symmetry.jl
+++ b/test/test_symmetry.jl
@@ -26,11 +26,23 @@
     end
     =#
 
+    # Test that the orbit of a Wyckoff matches the multiplicity data.
     for sgnum in 1:230
         sg = Sunny.Spacegroup(Sunny.standard_setting[sgnum])
         for (mult, letter, sitesym, pos) in Sunny.wyckoff_table[sgnum]
             orbit = Sunny.crystallographic_orbit(sg.symops, Sunny.WyckoffExpr(pos))
-            @assert length(orbit) == mult
+            @test length(orbit) == mult
+        end
+    end
+
+    # Test that Wyckoffs can be correctly inferred from a position
+    niters = 20
+    for sgnum in rand(1:230, niters)
+        for (_, letter, _, pos) in Sunny.wyckoff_table[sgnum]
+            w = Sunny.WyckoffExpr(pos)
+            θ = 10 * randn(3)
+            r = w.F * θ + w.c
+            @test letter == Sunny.find_wyckoff_for_position(sgnum, r; symprec=1e-8).letter
         end
     end
 end

From 9463d0d1d02adeb7ffaa8c3cadef2b1cc0e5bc32 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Sun, 3 Nov 2024 08:54:09 -0700
Subject: [PATCH 2/2] Tweak

---
 src/Symmetry/Crystal.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Symmetry/Crystal.jl b/src/Symmetry/Crystal.jl
index e21cdf310..d6e4d2a7e 100644
--- a/src/Symmetry/Crystal.jl
+++ b/src/Symmetry/Crystal.jl
@@ -397,7 +397,7 @@ function crystal_from_spacegroup(latvecs::Mat3, positions::Vector{Vec3}, types::
     wyckoffs = find_wyckoff_for_position.(Ref(sg), positions; symprec)
     orbits = crystallographic_orbits_distinct(sg.symops, positions; symprec, wyckoffs)
 
-    # Inferred orbits must match known multiplicities of the Wyckoffs
+    # Symmetry-propagated orbits must match known multiplicities of the Wyckoffs
     foreach(orbits, wyckoffs) do orbit, wyckoff
         @assert wyckoff.multiplicity ≈ length(orbit) / abs(det(sg.setting.R))
     end