Fix a bug in spglib interface (#911)

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"

src/common/split_evenly.jl

@@ -1,4 +1,3 @@
-
 """
 Split an iterable evenly into N chunks, which will be returned.
 """

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)

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;

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

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