add irreducible_kcoord and its tests

src/BaseMesh.jl

@@ -16,6 +16,7 @@ export UniformMesh, BaryChebMesh, CenteredMesh, EdgedMesh, UMesh, CompositeMesh
 abstract type AbstractUniformMesh{T,DIM} <: AbstractMesh{T,DIM} end
 export AbstractUniformMesh
 export inv_lattice_vector, lattice_vector, cell_volume
+export fractional_coordinates, cartesian_coordinates
 
 Base.length(mesh::AbstractUniformMesh) = prod(mesh.size)
 Base.size(mesh::AbstractUniformMesh) = mesh.size

src/BrillouinZoneMeshes.jl

@@ -16,6 +16,7 @@ BaryCheb = CompositeGrids.BaryChebTools
 include("AbstractMeshes.jl")
 using .AbstractMeshes
 export AbstractMeshes, AbstractMesh
+export fractional_coordinates
 
 include("printing.jl")
 

src/Cells.jl

@@ -72,18 +72,14 @@ struct Cell{T,DIM}
     cell_volume::T
     recip_cell_volume::T
 
-    # Particle types (elements) and particle positions and in fractional coordinates.
-    # Possibly empty. It's up to the `term_types` to make use of this (or not).
-    # `atom_groups` contains the groups of indices into atoms and positions, which
-    # point to identical atoms. It is computed automatically on Cell construction and may
-    # be used to optimise the term instantiation.
     atoms::Vector{Int}
     positions::Vector{SVector{DIM,T}}  # positions[i] is the location of atoms[i] in fract. coords
     atom_groups::Vector{Vector{Int}}  # atoms[i] == atoms[j] for all i, j in atom_group[α]
 
 
     # collection of all allowed G vectors
     G_vector::Vector{SVector{DIM,Int}}
+    symmetry::Vector{PointSymmetry.SymOp}
 end
 
 function Cell(;
@@ -117,7 +113,11 @@ function Cell(;
         G_vector = [SVector{DIM,Int}(zeros(DIM)),]
     end
 
-    return Cell{T,DIM}(lattice, recip_lattice, inv_lattice, inv_recip_lattice, cell_volume, recip_cell_volume, atoms, positions, atom_groups, G_vector)
+    cell = Cell{T,DIM}(lattice, recip_lattice, inv_lattice, inv_recip_lattice, cell_volume, recip_cell_volume, atoms, positions, atom_groups, G_vector, [])
+    for symop in default_symmetries(cell)
+        push!(cell.symmetry, symop)
+    end
+    return cell
 end
 
 function get_latvec(br::Cell, I::Int; isrecip=true)
@@ -128,10 +128,6 @@ function get_latvec(br::Cell, I::Int; isrecip=true)
     end
 end
 
-# normalize_magnetic_moment(::Nothing)::Vec3{Float64} = (0, 0, 0)
-# normalize_magnetic_moment(mm::Number)::Vec3{Float64} = (0, 0, mm)
-# normalize_magnetic_moment(mm::AbstractVector)::Vec3{Float64} = mm
-
 """
     function standard_cell(;
         dtype=Float64,
@@ -203,7 +199,8 @@ function standard_cell(;
         lattice=lattice,
         atoms=_atoms,
         positions=positions,
-        G_vector=G_vector)
+        G_vector=G_vector
+        )
 end
 
 """

src/symmetry/PointSymmetry.jl

@@ -17,6 +17,18 @@ const Mat3{T} = SMatrix{3,3,T,9} where {T}
 const Vec3{T} = SVector{3,T} where {T}
 const AbstractArray3{T} = AbstractArray{T,3}
 
+function _makeVec3(kcoord::AbstractVector{T}) where {T}
+    if length(kcoord) == 3
+        return Vec3{T}(kcoord)
+    elseif length(kcoord) == 2
+        return Vec3{T}(kcoord[1], kcoord[2], 0)
+    elseif length(kcoord) == 1
+        return Vec3{T}(kcoord[1], 0, 0)
+    else
+        error("illegal kcoord dimension!")
+    end
+end
+
 function _make3D(lattice::AbstractMatrix{T}, position::AbstractVector) where {T}
     @assert size(lattice, 1) == size(lattice, 2)
     DIM = size(lattice, 1)

src/symmetry/SymOp.jl

@@ -21,7 +21,7 @@ const SYMMETRY_TOLERANCE = 1e-5
 
 is_approx_integer(r; tol=SYMMETRY_TOLERANCE) = all(ri -> abs(ri - round(ri)) ≤ tol, r)
 
-struct SymOp{T <: Real}
+struct SymOp{T<:Real}
     # (Uu)(x) = u(W x + w) in real space
     W::Mat3{Int}
     w::Vec3{T}
@@ -37,7 +37,7 @@ function SymOp(W, w::AbstractVector{T}) where {T}
     SymOp{T}(W, w, S, τ)
 end
 
-function Base.convert(::Type{SymOp{T}}, S::SymOp{U}) where {U <: Real, T <: Real}
+function Base.convert(::Type{SymOp{T}}, S::SymOp{U}) where {U<:Real,T<:Real}
     SymOp{T}(S.W, T.(S.w), S.S, T.(S.τ))
 end
 
@@ -56,7 +56,7 @@ function Base.:*(op1::SymOp, op2::SymOp)
     w = op1.w + op1.W * op2.w
     SymOp(W, w)
 end
-Base.inv(op::SymOp) = SymOp(inv(op.W), -op.W\op.w)
+Base.inv(op::SymOp) = SymOp(inv(op.W), -op.W \ op.w)
 
 function check_group(symops::Vector; kwargs...)
     is_approx_in_symops(s1) = any(s -> isapprox(s, s1; kwargs...), symops)
@@ -66,10 +66,16 @@ function check_group(symops::Vector; kwargs...)
             error("check_group: symop $s with inverse $(inv(s)) is not in the group")
         end
         for s2 in symops
-            if !is_approx_in_symops(s*s2) || !is_approx_in_symops(s2*s)
+            if !is_approx_in_symops(s * s2) || !is_approx_in_symops(s2 * s)
                 error("check_group: product is not stable")
             end
         end
     end
     symops
 end
+
+# Get the rotation axis of a symmetry operation
+# https://en.wikipedia.org/w/index.php?title=Rotation_matrix&action=edit&section=11
+# function axis(r::Mat3)
+
+# end

src/symmetry/spglib.jl

@@ -127,15 +127,17 @@ function spglib_get_symmetry(lattice::AbstractMatrix{<:AbstractFloat}, atom_grou
         end
     end
 
-    Ws, ws
+    return Ws, ws
 end
 
-
+# The irreducible k-points are searched from unique k-point mesh grids from direct (real space) basis vectors 
+# and a set of rotation parts of symmetry operations in direct space with one or multiple stabilizers.
 function spglib_get_stabilized_reciprocal_mesh(kgrid_size, rotations::Vector;
     is_shift=Vec3(0, 0, 0),
     is_time_reversal=false,
     qpoints=[Vec3(0.0, 0.0, 0.0)])
     spg_rotations = cat([copy(Cint.(S')) for S in rotations]..., dims=3)
+
     nkpt = prod(kgrid_size)
     mapping = Vector{Cint}(undef, nkpt)
     grid_address = Matrix{Cint}(undef, 3, nkpt)

src/symmetry/symmetry.jl

@@ -24,19 +24,21 @@
 # (Uu)(G) = e^{-i G τ} u(S^-1 G)
 # In particular, we can choose the eigenvectors at Sk as u_{Sk} = U u_k
 
-# We represent then the BZ as a set of irreducible points `kpoints`,
-# and a set of weights `kweights` (summing to 1). The value of
-# observables is given by a weighted sum over the irreducible kpoints,
-# plus a symmetrization operation (which depends on the particular way
-# the observable transforms under the symmetry).
-
 # There is by decreasing cardinality
 # - The group of symmetry operations of the lattice
-# - The group of symmetry operations of the crystal (model.symmetries)
-# - The group of symmetry operations of the crystal that preserves the BZ mesh (basis.symmetries)
+# - The group of symmetry operations of the crystal (cell.symmetry)
+# - The group of symmetry operations of the crystal that preserves the BZ mesh (mesh.symmetry)
 
 # See https://juliamolsim.github.io/DFTK.jl/stable/advanced/symmetries for details.
 
+"""Bring ``k``-point coordinates into the range [-0.5, 0.5)"""
+function normalize_kpoint_coordinate(x::Real)
+    x = x - round(Int, x, RoundNearestTiesUp)
+    @assert -0.5 ≤ x < 0.5
+    x
+end
+normalize_kpoint_coordinate(k::AbstractVector) = normalize_kpoint_coordinate.(k)
+
 @doc raw"""
 Return the ``k``-point symmetry operations associated to a lattice and atoms.
 
@@ -50,41 +52,24 @@ function symmetry_operations(lattice, atoms, positions, magnetic_moments=[];
     @assert length(atoms) == length(positions)
     atom_groups = [findall(Ref(pot) .== atoms) for pot in Set(atoms)]
     Ws, ws = spglib_get_symmetry(lattice, atom_groups, positions, magnetic_moments; tol_symmetry)
-    [SymOp(W, w) for (W, w) in zip(Ws, ws)]
+    return [SymOp(W, w) for (W, w) in zip(Ws, ws)]
 end
 
 # Approximate in; can be performance-critical, so we optimize in case of rationals
 is_approx_in_(x::AbstractArray{<:Rational}, X) = any(isequal(x), X)
 is_approx_in_(x::AbstractArray{T}, X) where {T} = any(y -> isapprox(x, y; atol=sqrt(eps(T))), X)
 
-"""
-Filter out the symmetry operations that don't respect the symmetries of the discrete BZ grid
-"""
-function symmetries_preserving_kgrid(symmetries, kcoords)
-    kcoords_normalized = normalize_kpoint_coordinate.(kcoords)
-    function preserves_grid(symop)
-        all(is_approx_in_(normalize_kpoint_coordinate(symop.S * k), kcoords_normalized)
-            for k in kcoords_normalized)
-    end
-    filter(preserves_grid, symmetries)
-end
-
-"""
-Filter out the symmetry operations that don't respect the symmetries of the discrete real-space grid
-"""
-function symmetries_preserving_rgrid(symmetries, fft_size)
-    is_in_grid(r) =
-        all(zip(r, fft_size)) do (ri, size)
-            abs(ri * size - round(ri * size)) / size ≤ SYMMETRY_TOLERANCE
-        end
-
-    onehot3(i) = (x = zeros(Bool, 3); x[i] = true; Vec3(x))
-    function preserves_grid(symop)
-        all(is_in_grid(symop.W * onehot3(i) .// fft_size[i] + symop.w) for i = 1:3)
-    end
-
-    filter(preserves_grid, symmetries)
-end
+# """
+# Filter out the symmetry operations that don't respect the symmetries of the discrete BZ grid
+# """
+# function symmetries_preserving_kgrid(symmetries, kcoords)
+#     kcoords_normalized = normalize_kpoint_coordinate.(kcoords)
+#     function preserves_grid(symop)
+#         all(is_approx_in_(normalize_kpoint_coordinate(symop.S * k), kcoords_normalized)
+#             for k in kcoords_normalized)
+#     end
+#     filter(preserves_grid, symmetries)
+# end
 
 @doc raw"""
 Apply various standardisations to a lattice and a list of atoms. It uses spglib to detect
@@ -101,71 +86,47 @@ function standardize_atoms(lattice, atoms, positions, magnetic_moments=[]; kwarg
     (; ret.lattice, atoms, ret.positions, ret.magnetic_moments)
 end
 
-# find where in the irreducible basis `basis_irred` the k-point `kpt_unfolded` is handled
-function unfold_mapping(basis_irred, kpt_unfolded)
-    for ik_irred = 1:length(basis_irred.kpoints)
-        kpt_irred = basis_irred.kpoints[ik_irred]
-        for symop in basis_irred.symmetries
-            Sk_irred = normalize_kpoint_coordinate(symop.S * kpt_irred.coordinate)
-            k_unfolded = normalize_kpoint_coordinate(kpt_unfolded.coordinate)
-            if (Sk_irred ≈ k_unfolded) && (kpt_unfolded.spin == kpt_irred.spin)
-                return ik_irred, symop
-            end
-        end
-    end
-    error("Invalid unfolding of BZ")
-end
+# """"
+# Convert a `kcoords` into one that doesn't use BZ symmetry.
+# This is mainly useful for debug purposes (e.g. in cases we don't want to
+# bother thinking about symmetries).
+# """
+# function unfold_kcoords(kcoords, symmetries)
+#     all_kcoords = [normalize_kpoint_coordinate(symop.S * kcoord)
+#                    for kcoord in kcoords, symop in symmetries]
+
+#     # the above multiplications introduce an error
+#     unique(all_kcoords) do k
+#         digits = ceil(Int, -log10(SYMMETRY_TOLERANCE))
+#         normalize_kpoint_coordinate(round.(k; digits))
+#     end
+# end
 
-function unfold_array_(basis_irred, basis_unfolded, data, is_ψ)
-    if basis_irred == basis_unfolded
-        return data
-    end
-    if !(basis_irred.comm_kpts == basis_irred.comm_kpts == MPI.COMM_WORLD)
-        error("Brillouin zone symmetry unfolding not supported with MPI yet")
-    end
-    data_unfolded = similar(data, length(basis_unfolded.kpoints))
-    for ik_unfolded in 1:length(basis_unfolded.kpoints)
-        kpt_unfolded = basis_unfolded.kpoints[ik_unfolded]
-        ik_irred, symop = unfold_mapping(basis_irred, kpt_unfolded)
-        if is_ψ
-            # transform ψ_k from data into ψ_Sk in data_unfolded
-            kunfold_coord = kpt_unfolded.coordinate
-            @assert normalize_kpoint_coordinate(kunfold_coord) ≈ kunfold_coord
-            _, ψSk = apply_symop(symop, basis_irred,
-                basis_irred.kpoints[ik_irred], data[ik_irred])
-            data_unfolded[ik_unfolded] = ψSk
-        else
-            # simple copy
-            data_unfolded[ik_unfolded] = data[ik_irred]
+"""
+Filter out the kcoords that are irreducible under the given symmetries
+"""
+function irreducible_kcoord(symmetries::AbstractVector{SymOp}, kcoords::AbstractVector)
+    kcoords = _makeVec3.(kcoords)
+    kidxmap = [i for i in 1:length(kcoords)]
+    for (ki, k) in enumerate(kcoords)
+        if kidxmap[ki] != ki # already mapped to a symmetry equivalent
+            continue
+        end
+        symk = [symop.S * k for symop in symmetries]
+        for kj in ki+1:length(kcoords)
+            if kidxmap[kj] != kj # already mapped to a symmetry equivalent
+                continue
+            end
+            for s in 1:length(symmetries)
+                # If the difference between kred and W' * k == W^{-1} * k
+                # is only integer in fractional reciprocal-space coordinates, then
+                # kred and S' * k are equivalent k-points
+                if all(isinteger, kcoords[kj] - symk[s])
+                    kidxmap[kj] = ki
+                end
+            end
         end
     end
-    data_unfolded
-end
-
-function unfold_bz(scfres)
-    basis_unfolded = unfold_bz(scfres.basis)
-    ψ = unfold_array_(scfres.basis, basis_unfolded, scfres.ψ, true)
-    eigenvalues = unfold_array_(scfres.basis, basis_unfolded, scfres.eigenvalues, false)
-    occupation = unfold_array_(scfres.basis, basis_unfolded, scfres.occupation, false)
-    E, ham = energy_hamiltonian(basis_unfolded, ψ, occupation;
-        scfres.ρ, eigenvalues, scfres.εF)
-    @assert E.total ≈ scfres.energies.total
-    new_scfres = (; basis=basis_unfolded, ψ, ham, eigenvalues, occupation)
-    merge(scfres, new_scfres)
+    return kidxmap
 end
 
-""""
-Convert a `kcoords` into one that doesn't use BZ symmetry.
-This is mainly useful for debug purposes (e.g. in cases we don't want to
-bother thinking about symmetries).
-"""
-function unfold_kcoords(kcoords, symmetries)
-    all_kcoords = [normalize_kpoint_coordinate(symop.S * kcoord)
-                   for kcoord in kcoords, symop in symmetries]
-
-    # the above multiplications introduce an error
-    unique(all_kcoords) do k
-        digits = ceil(Int, -log10(SYMMETRY_TOLERANCE))
-        normalize_kpoint_coordinate(round.(k; digits))
-    end
-end

test/UniformMeshMap.jl

@@ -26,6 +26,23 @@
         end
         @test sort(klist) == sort(meshmap.irreducible_indices)
 
+        ##### test home-made irreducible_kcoord ##############
+        kidxmap = PointSymmetry.irreducible_kcoord(br.symmetry, [inv_lattice_vector(brmesh) * brmesh[i] for i in 1:length(brmesh)])
+        ir_klist = sort(unique(kidxmap))
+
+        # for kpoint in ir_klist
+        #     kvec = inv_lattice_vector(brmesh) * brmesh[kpoint]
+        #     println(kvec)
+        #     # _kvec = BrillouinZoneMeshes.Cells.recip_vector_red_to_cart(br, meshmap.kcoords_global[kpoint])
+        #     # @test kvec ≈ _kvec
+        # end
+
+        @test length(ir_klist) == length(meshmap.irreducible_indices)
+        for (idx, ir_idx) in enumerate(kidxmap)
+            # println(idx, ": ", fractional_coordinates(brmesh, idx), " -> ", ir_idx, " : ", fractional_coordinates(brmesh, ir_idx))
+            @test meshmap.map[idx] == meshmap.map[ir_idx]
+        end
+
         return meshmap, brmesh
     end
     size, shift = [4, 4, 4], false
@@ -40,6 +57,7 @@
     kweights = [0.015625, 0.125, 0.0625, 0.09375, 0.375, 0.1875, 0.046875, 0.09375]
     test(3, silicon.lattice, silicon.atoms, silicon.positions, size, shift;
         kcoords=kcoords, kweights=kweights)
+    exit(0)
 
     size, shift = [4, 4, 4], true
     kcoords = [