From 430139f85d22120cc8fefa3982702c3abba6ba9e Mon Sep 17 00:00:00 2001 From: Kun Chen Date: Wed, 9 Nov 2022 18:02:02 -0500 Subject: [PATCH] add irreducible_kcoord and its tests --- src/BaseMesh.jl | 1 + src/BrillouinZoneMeshes.jl | 1 + src/Cells.jl | 19 ++-- src/symmetry/PointSymmetry.jl | 12 +++ src/symmetry/SymOp.jl | 14 ++- src/symmetry/spglib.jl | 6 +- src/symmetry/symmetry.jl | 161 +++++++++++++--------------------- test/UniformMeshMap.jl | 18 ++++ 8 files changed, 115 insertions(+), 117 deletions(-) diff --git a/src/BaseMesh.jl b/src/BaseMesh.jl index a2a51e9..caf6a75 100644 --- a/src/BaseMesh.jl +++ b/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 diff --git a/src/BrillouinZoneMeshes.jl b/src/BrillouinZoneMeshes.jl index c8f6d13..24ebcc0 100644 --- a/src/BrillouinZoneMeshes.jl +++ b/src/BrillouinZoneMeshes.jl @@ -16,6 +16,7 @@ BaryCheb = CompositeGrids.BaryChebTools include("AbstractMeshes.jl") using .AbstractMeshes export AbstractMeshes, AbstractMesh +export fractional_coordinates include("printing.jl") diff --git a/src/Cells.jl b/src/Cells.jl index ec89879..14e5e58 100644 --- a/src/Cells.jl +++ b/src/Cells.jl @@ -72,11 +72,6 @@ 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[α] @@ -84,6 +79,7 @@ struct Cell{T,DIM} # 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 """ diff --git a/src/symmetry/PointSymmetry.jl b/src/symmetry/PointSymmetry.jl index 5189b70..6177639 100644 --- a/src/symmetry/PointSymmetry.jl +++ b/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) diff --git a/src/symmetry/SymOp.jl b/src/symmetry/SymOp.jl index d376830..871740a 100644 --- a/src/symmetry/SymOp.jl +++ b/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§ion=11 +# function axis(r::Mat3) + +# end diff --git a/src/symmetry/spglib.jl b/src/symmetry/spglib.jl index 17848e6..fd43675 100644 --- a/src/symmetry/spglib.jl +++ b/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) diff --git a/src/symmetry/symmetry.jl b/src/symmetry/symmetry.jl index c6be44d..b79f804 100644 --- a/src/symmetry/symmetry.jl +++ b/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 diff --git a/test/UniformMeshMap.jl b/test/UniformMeshMap.jl index d412294..b4c242d 100644 --- a/test/UniformMeshMap.jl +++ b/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 = [