From 5b6ce69788f28ec836a191ab70bfac8180c13dce Mon Sep 17 00:00:00 2001
From: "Michael F. Herbst" <info@michael-herbst.com>
Date: Wed, 22 Nov 2023 17:03:15 +0100
Subject: [PATCH] Fix a bug in spglib interface (#911)

---
 Project.toml               |  2 +-
 src/common/split_evenly.jl |  1 -
 src/external/spglib.jl     | 27 ++++++++++++++-------------
 src/symmetry.jl            |  4 ++--
 test/bzmesh_symmetry.jl    |  1 +
 test/symmetry_issues.jl    | 27 +++++++++++++++++++++++++++
 6 files changed, 45 insertions(+), 17 deletions(-)
 create mode 100644 test/symmetry_issues.jl

diff --git a/Project.toml b/Project.toml
index f027ee98af..61bb8d1384 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DFTK"
 uuid = "acf6eb54-70d9-11e9-0013-234b7a5f5337"
 authors = ["Michael F. Herbst <info@michael-herbst.com>", "Antoine Levitt <antoine.levitt@inria.fr>"]
-version = "0.6.13"
+version = "0.6.14"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
diff --git a/src/common/split_evenly.jl b/src/common/split_evenly.jl
index d6a56b5f2a..1b57cad1e6 100644
--- a/src/common/split_evenly.jl
+++ b/src/common/split_evenly.jl
@@ -1,4 +1,3 @@
-
 """
 Split an iterable evenly into N chunks, which will be returned.
 """
diff --git a/src/external/spglib.jl b/src/external/spglib.jl
index b4bf7847d7..5bbd23ddb5 100644
--- a/src/external/spglib.jl
+++ b/src/external/spglib.jl
@@ -5,24 +5,25 @@
 import Spglib
 
 function spglib_cell(lattice, atom_groups, positions, magnetic_moments)
-    offset  = 0
-    atoms   = zeros(Int, length(positions))
-    magmoms = zeros(length(positions))
+    spg_atoms     = Int[]
+    spg_magmoms   = Float64[]
+    spg_positions = Vector{Float64}[]
+
     arbitrary_spin = false
-    for (igroup, indices) in enumerate(atom_groups)
-        for iatom in indices
-            offset += 1
-            atoms[offset] = igroup
+    for (igroup, indices) in enumerate(atom_groups), iatom in indices
+        push!(spg_atoms, igroup)
+        push!(spg_positions, positions[iatom])
 
-            if !isempty(magnetic_moments)
-                magmom = magnetic_moments[iatom]
-                magmoms[offset] = magmom[3]
-                !iszero(magmom[1:2]) && (arbitrary_spin = true)
-            end
+        if isempty(magnetic_moments)
+            magmom = zeros(3)
+        else
+            magmom = magnetic_moments[iatom]
+            !iszero(magmom[1:2]) && (arbitrary_spin = true)
         end
+        push!(spg_magmoms, magmom[3])
     end
     @assert !arbitrary_spin
-    Spglib.SpglibCell(lattice, positions, atoms, magmoms)
+    Spglib.SpglibCell(lattice, spg_positions, spg_atoms, spg_magmoms)
 end
 function spglib_cell(system::AbstractSystem)
     parsed = parse_system(system)
diff --git a/src/symmetry.jl b/src/symmetry.jl
index 25a7fd3362..6cf4e8b47e 100644
--- a/src/symmetry.jl
+++ b/src/symmetry.jl
@@ -139,10 +139,10 @@ on each of the atoms. The symmetries are determined using spglib.
     end
     symmetries
 end
-function symmetry_operations(system::AbstractSystem; tol_symmetry=SYMMETRY_TOLERANCE)
+function symmetry_operations(system::AbstractSystem; kwargs...)
     parsed = parse_system(system)
     symmetry_operations(parsed.lattice, parsed.atoms, parsed.positions,
-                        parsed.magnetic_moments; tol_symmetry)
+                        parsed.magnetic_moments; kwargs...)
 end
 
 function _check_symmetries(symmetries::AbstractVector{<:SymOp}, lattice, atom_groups, positions;
diff --git a/test/bzmesh_symmetry.jl b/test/bzmesh_symmetry.jl
index 89b60fe236..ea38b4df0c 100644
--- a/test/bzmesh_symmetry.jl
+++ b/test/bzmesh_symmetry.jl
@@ -21,6 +21,7 @@
 
         model_sym = model_LDA(testcase.lattice, testcase.atoms, testcase.positions)
         basis = PlaneWaveBasis(model_sym; Ecut=5, case...)
+        DFTK.check_group(basis.symmetries)
         scfres = self_consistent_field(basis; is_converged=DFTK.ScfConvergenceDensity(1e-10))
         ρ2 = scfres.ρ
         E2 = scfres.energies.total
diff --git a/test/symmetry_issues.jl b/test/symmetry_issues.jl
new file mode 100644
index 0000000000..e39926a733
--- /dev/null
+++ b/test/symmetry_issues.jl
@@ -0,0 +1,27 @@
+# This file collects examples, where we had issues with symmetries (symmetry determination,
+# k-point reduction, etc.) which are now resolved. Should make sure we don't reintroduce bugs.
+
+@testitem "Symmetry issues" begin
+    using DFTK
+    using DFTK: spglib_dataset
+    using Unitful
+    using UnitfulAtomic
+    using AtomsBase
+
+    @testset "CuO2" begin
+        a = 4.474
+        lattice = [[0, a, a], [a, 0, a], [a, a, 0]]u"bohr"
+        x = 6.711
+        y = 2.237
+        atoms = [
+            Atom(:Cu, [0, 0, 0]u"bohr", magnetic_moment=0),
+            Atom(:O,  [x, y, x]u"bohr", magnetic_moment=0),
+            Atom(:O,  [x, y, y]u"bohr", magnetic_moment=0),
+        ]
+        system = periodic_system(atoms, lattice)
+
+        symmetries = DFTK.symmetry_operations(system; check_symmetry=true)
+        @test length(symmetries) == 48
+        @test spglib_dataset(system).spacegroup_number == 225
+    end
+end