[WIP] Time-reversal symmetry for kpoint reduction

src/DFTK.jl

@@ -175,7 +175,6 @@ export compute_stresses_cart
 include("postprocess/stresses.jl")
 export compute_dos
 export compute_ldos
-export compute_nos
 export plot_dos
 include("postprocess/dos.jl")
 export compute_χ0

src/Model.jl

@@ -94,7 +94,7 @@ function Model(lattice::AbstractMatrix{T};
                temperature=T(0.0),
                smearing=nothing,
                spin_polarization=default_spin_polarization(magnetic_moments),
-               symmetries=default_symmetries(lattice, atoms, magnetic_moments, terms, spin_polarization),
+               symmetries=true,
                ) where {T <: Real}
     lattice = Mat3{T}(lattice)
     temperature = T(austrip(temperature))
@@ -181,7 +181,7 @@ end
 Default logic to determine the symmetry operations to be used in the model.
 """
 function default_symmetries(lattice, atoms, magnetic_moments, terms, spin_polarization;
-                        tol_symmetry=1e-5)
+                        tol_symmetry=SYMMETRY_TOLERANCE)
     dimension = count(!iszero, eachcol(lattice))
     if spin_polarization == :full || dimension != 3
         return [one(SymOp)]  # Symmetry not supported in spglib
@@ -193,8 +193,9 @@ function default_symmetries(lattice, atoms, magnetic_moments, terms, spin_polari
     else
         magnetic_moments = [el => normalize_magnetic_moment.(magmoms)
                             for (el, magmoms) in magnetic_moments]
-        return symmetry_operations(lattice, atoms, magnetic_moments,
-                                   tol_symmetry=tol_symmetry)
+        is_time_reversal = !any(breaks_time_reversal, terms)
+        return symmetry_operations(lattice, atoms, magnetic_moments;
+                                   is_time_reversal, tol_symmetry)
     end
 end
 

src/PlaneWaveBasis.jl

@@ -93,7 +93,7 @@ struct PlaneWaveBasis{T} <: AbstractBasis{T}
     #                                       respective rank in comm_kpts
 
     # Symmetry operations that leave the reducible Brillouin zone invariant.
-    # Subset of model.symmetries, and superset of all the ksymops.
+    # Subset of model.symmetries.
     # Nearly all computations will be done inside this symmetry group;
     # the exception is inexact operations on the FFT grid (ie xc),
     # which doesn't respect the symmetry
@@ -240,7 +240,7 @@ end
 @timing function PlaneWaveBasis(model::Model{T}, Ecut::Number,
                                 kcoords::AbstractVector, ksymops,
                                 symmetries=symmetries_preserving_kgrid(model.symmetries,
-                                                                       kcoords, ksymops);
+                                                                       unfold_kcoords(kcoords, vcat(ksymops...)));
                                 fft_size=nothing, variational=true,
                                 fft_size_algorithm=:fast, supersampling=2,
                                 kgrid=nothing, kshift=nothing,
@@ -276,7 +276,8 @@ Creates a new basis identical to `basis`, but with a custom set of kpoints
 @timing function PlaneWaveBasis(basis::PlaneWaveBasis, kcoords::AbstractVector,
                                 ksymops::AbstractVector)
     kgrid = kshift = nothing
-    symmetries = symmetries_preserving_kgrid(basis.model.symmetries, kcoords, ksymops)
+    symmetries = symmetries_preserving_kgrid(basis.model.symmetries,
+                                             unfold_kcoords(kcoords, basis.model.symmetries))
     PlaneWaveBasis(basis.model, basis.Ecut,
                    basis.fft_size, basis.variational,
                    kcoords, ksymops, kgrid, kshift,

src/SymOp.jl

@@ -12,34 +12,40 @@
 # (all these formulas are valid both in reduced and cartesian coordinates)
 
 # Time-reversal symmetries are the anti-unitaries
-# (Uu)(x) = conj(u(Wx+w))
+# (Uu)(x) = (θu)(Wx+w)
 # or in Fourier space
-# (Uu)(G) = e^{i G τ} conj(u(-S^-1 G))
+# (Uu)(G) = e^{-i G τ} (θu)(S^-1 G), with the modified S = -W', τ = W^-1 w
+# where θ is the conjugation operator
 
 # Tolerance to consider two atomic positions as equal (in relative coordinates)
 const SYMMETRY_TOLERANCE = 1e-5
 
 struct SymOp{T <: Real}
-    # (Uu)(x) = u(W x + w) in real space
+    # (Uu)(x) = (θu)(W x + w) in real space
     W::Mat3{Int}
     w::Vec3{T}
+    θ::Bool # true for TRS, false for no TRS
 
-    # (Uu)(G) = e^{-i G τ} u(S^-1 G) in reciprocal space
+    # (Uu)(G) = e^{-i G τ} (θu)(S^-1 G) in reciprocal space
     S::Mat3{Int}
     τ::Vec3{T}
 
-    function SymOp(W, w)
+    function SymOp(W, w, θ=false)
         w = mod.(w, 1)
         S = W'
         τ = -W \w
-        new{eltype(τ)}(W, w, S, τ)
+        if θ
+            S = -S
+            τ = -τ
+        end
+        new{eltype(τ)}(W, w, θ, S, τ)
     end
 end
 
-Base.:(==)(op1::SymOp, op2::SymOp) = op1.W == op2.W && op1.w == op2.w
+Base.:(==)(op1::SymOp, op2::SymOp) = op1.W == op2.W && op1.w == op2.w && op1.θ == op2.θ
 function Base.isapprox(op1::SymOp, op2::SymOp; atol=SYMMETRY_TOLERANCE)
     is_approx_integer(r) = all(ri -> abs(ri - round(ri)) ≤ atol, r)
-    op1.W == op2.W && is_approx_integer(op1.w - op2.w)
+    op1.W == op2.W && is_approx_integer(op1.w - op2.w) && op1.θ == op2.θ
 end
 Base.one(::Type{SymOp}) = SymOp(Mat3{Int}(I), Vec3(zeros(Bool, 3)))
 Base.one(::SymOp) = one(SymOp)
@@ -48,9 +54,9 @@ Base.one(::SymOp) = one(SymOp)
 function Base.:*(op1, op2)
     W = op1.W * op2.W
     w = op1.w + op1.W * op2.w
-    SymOp(W, w)
+    SymOp(W, w, xor(op1.θ, op2.θ))
 end
-Base.inv(op) = SymOp(inv(op.W), -op.W\op.w)
+Base.inv(op) = SymOp(inv(op.W), -op.W\op.w, op.θ)
 
 function check_group(symops::Vector; kwargs...)
     is_approx_in_symops(s1) = any(s -> isapprox(s, s1; kwargs...), symops)

src/bzmesh.jl

@@ -62,13 +62,13 @@ function bzmesh_ir_wedge(kgrid_size, symmetries; kshift=[0, 0, 0])
     # Filter those symmetry operations (S,τ) that preserve the MP grid
     symmetries = symmetries_preserving_kgrid(symmetries, kpoints_mp)
 
+    is_time_reversal = any(symop -> symop.θ, symmetries)
+
     # Give the remaining symmetries to spglib to compute an irreducible k-point mesh
-    # TODO implement time-reversal symmetry and turn the flag to true
     is_shift = Int.(2 * kshift)
-    Ws = [symop.S' for symop in symmetries]
+    Ws = unique([symop.W for symop in symmetries])
     _, mapping, grid = spglib_get_stabilized_reciprocal_mesh(
-        kgrid_size, Ws, is_shift=is_shift, is_time_reversal=false
-    )
+        kgrid_size, Ws; is_shift, is_time_reversal)
     # Convert irreducible k-points to DFTK conventions
     kgrid_size = Vec3{Int}(kgrid_size)
     kirreds = [(kshift .+ grid[ik + 1]) .// kgrid_size
@@ -113,9 +113,7 @@ function bzmesh_ir_wedge(kgrid_size, symmetries; kshift=[0, 0, 0])
 
     if !isempty(kreds_notmapped)
         # add them as reducible anyway
-        Ws = [ symop.S'           for symop in symmetries]
-        ws = [-symop.S' * symop.τ for symop in symmetries]
-        eirreds, esymops = find_irreducible_kpoints(kreds_notmapped, Ws, ws)
+        eirreds, esymops = find_irreducible_kpoints(kreds_notmapped, symmetries)
         @info("$(length(kreds_notmapped)) reducible kpoints could not be generated from " *
               "the irreducible kpoints returned by spglib. $(length(eirreds)) of " *
               "these are added as extra irreducible kpoints.")

src/postprocess/current.jl

@@ -3,7 +3,7 @@
 Computes the *probability* (not charge) current, ∑ fn Im(ψn* ∇ψn)
 """
 function compute_current(basis::PlaneWaveBasis, ψ, occupation)
-    @assert all(symops -> length(symops) == 1, basis.ksymops) == 1  # TODO lift this
+    @assert length(basis.symmetries) == 1 # TODO lift this
     current = [zeros(eltype(basis), basis.fft_size) for α = 1:3]
     for (ik, kpt) in enumerate(basis.kpoints)
         for (n, ψnk) in enumerate(eachcol(ψ[ik]))

src/postprocess/dos.jl

@@ -7,45 +7,6 @@
 # LDOS (local density of states)
 # LD = sum_n f_n ψn
 
-@doc raw"""
-    compute_nos(ε, basis, eigenvalues; smearing=basis.model.smearing,
-                temperature=basis.model.temperature)
-
-The number of Kohn-Sham states in a temperature window of width `temperature` around the
-energy `ε` contributing to the DOS at temperature `T`.
-
-This quantity is not a physical quantity, but rather a dimensionless approximate measure
-for how well properties near the Fermi surface are sampled with the passed `smearing`
-and temperature `T`. It increases with both `T` and better sampling of the BZ with
-``k``-points. A value ``\gg 1`` indicates a good sampling of properties near the
-Fermi surface.
-"""
-function compute_nos(ε, basis, eigenvalues; smearing=basis.model.smearing,
-                     temperature=basis.model.temperature)
-    N = zeros(typeof(ε), basis.model.n_spin_components)
-    if (temperature == 0) || smearing isa Smearing.None
-        error("compute_nos only supports finite temperature")
-    end
-
-    # Note the differences to the DOS and LDOS functions: We are not counting states
-    # per BZ volume (like in DOS), but absolute number of states. Therefore n_symeqk
-    # is used instead of kweigths. We count the states inside a temperature window
-    # `temperature` centred about εik. For this the states are weighted by the distribution
-    # -f'((εik - ε)/temperature).
-    #
-    # To explicitly show the similarity with DOS and the temperature dependence we employ
-    # -f'((εik - ε)/temperature) = temperature * ( d/dε f_τ(εik - ε') )|_{ε' = ε}
-    for σ in 1:basis.model.n_spin_components, ik = krange_spin(basis, σ)
-        n_symeqk = length(basis.ksymops[ik])  # Number of symmetry-equivalent k-points
-        for (iband, εnk) in enumerate(eigenvalues[ik])
-            enred = (εnk - ε) / temperature
-            N[σ] -= n_symeqk * Smearing.occupation_derivative(smearing, enred)
-        end
-    end
-    N = mpi_sum(N, basis.comm_kpts)
-end
-
-
 """
 Total density of states at energy ε
 """

src/symmetry.jl

@@ -24,29 +24,31 @@
 # (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`, a
-# list of symmetry operations `ksymops` allowing the reconstruction of
-# the full (reducible) BZ, and a set of weights `kweights` (summing to
-# 1). The value of an observable (eg energy) per unit cell is given as
-# the value of that observable at each irreducible k-point, weighted by
-# kweight
+# 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 set of symmetry operations that we use to reduce the
-#   reducible Brillouin zone (RBZ) to the irreducible (IBZ) (basis.ksymops)
 
 # See https://juliamolsim.github.io/DFTK.jl/stable/advanced/symmetries for details.
 
 @doc raw"""
-Return the ``k``-point symmetry operations associated to a lattice and atoms.
+Return the symmetry operations associated to a lattice and atoms.
 """
-function symmetry_operations(lattice, atoms, magnetic_moments=[]; tol_symmetry=SYMMETRY_TOLERANCE)
+function symmetry_operations(lattice, atoms, magnetic_moments=[];
+                             is_time_reversal=false, tol_symmetry=SYMMETRY_TOLERANCE)
     Ws, ws = spglib_get_symmetry(lattice, atoms, magnetic_moments; tol_symmetry)
-    symmetries = [SymOp(W, w) for (W, w) in zip(Ws, ws)]
-    unique(symmetries)
+    symmetries = ([SymOp(W, w) for (W, w) in zip(Ws, ws)])
+    if is_time_reversal
+        symmetries = vcat(symmetries,
+                          [SymOp(symop.W, symop.w, true) for symop in symmetries])
+    end
+    symmetries
 end
 
 """
@@ -64,25 +66,11 @@ function symmetries_preserving_kgrid(symmetries, kcoords)
     filter(preserves_grid, symmetries)
 end
 
-
-"""
-Find the subset of symmetries compatible with the grid induced by the given kcoords and ksymops
-"""
-function symmetries_preserving_kgrid(symmetries, kcoords, ksymops)
-    new_symmetries = symmetries_preserving_kgrid(symmetries, unfold_kcoords(kcoords, ksymops))
-    # check for inconsistent ksymops/symmetries
-    if !all(s1 -> any(s2 -> isapprox(s1, s2), new_symmetries), Iterators.flatten(ksymops))
-        error("symmetries_preserving_kgrid: ksymops must be a subset of symmetries")
-    end
-    new_symmetries
-end
-
-
 """
 Implements a primitive search to find an irreducible subset of kpoints
 amongst the provided kpoints.
 """
-function find_irreducible_kpoints(kcoords, Ws, ws)
+function find_irreducible_kpoints(kcoords, symmetries)
     # This function is required because spglib sometimes flags kpoints
     # as reducible, where we cannot find a symmetry operation to
     # generate them from the provided irreducible kpoints. This
@@ -103,16 +91,16 @@ function find_irreducible_kpoints(kcoords, Ws, ws)
         kcoords_mapped[ik] = true
 
         for jk in findall(.!kcoords_mapped)
-            isym = findfirst(1:length(Ws)) do isym
-                # If the difference between kred and W' * k == W^{-1} * k
+            isym = findfirst(1:length(symmetries)) do isym
+                # If the difference between kred and S*k
                 # is only integer in fractional reciprocal-space coordinates, then
-                # kred and S' * k are equivalent k-points
-                all(isinteger, kcoords[jk] - (Ws[isym]' * kcoords[ik]))
+                # kred and S * k are equivalent k-points
+                all(isinteger, kcoords[jk] - (symmetries[isym].S * kcoords[ik]))
             end
 
             if !isnothing(isym)  # Found a reducible k-point
                 kcoords_mapped[jk] = true
-                push!(thisk_symops, SymOp(Ws[isym], ws[isym]))
+                push!(thisk_symops, symmetries[isym])
             end
         end  # jk
 
@@ -126,9 +114,9 @@ Apply a symmetry operation to eigenvectors `ψk` at a given `kpoint` to obtain a
 equivalent point in [-0.5, 0.5)^3 and associated eigenvectors (expressed in the
 basis of the new ``k``-point).
 """
-function apply_ksymop(ksymop::SymOp, basis, kpoint, ψk::AbstractVecOrMat)
-    S, τ = ksymop.S, ksymop.τ
-    ksymop == one(SymOp) && return kpoint, ψk
+function apply_symop(symop::SymOp, basis, kpoint, ψk::AbstractVecOrMat)
+    S, τ, θ= symop.S, symop.τ, symop.θ
+    symop == one(SymOp) && return kpoint, ψk
 
     # Apply S and reduce coordinates to interval [-0.5, 0.5)
     # Doing this reduction is important because
@@ -162,7 +150,8 @@ function apply_ksymop(ksymop::SymOp, basis, kpoint, ψk::AbstractVecOrMat)
         for (ig, G_full) in enumerate(Gs_full)
             igired = index_G_vectors(basis, kpoint, invS * G_full)
             @assert igired !== nothing
-            ψSk[ig, iband] = cis(-2π * dot(G_full, τ)) * ψk[igired, iband]
+            ψSk[ig, iband] = cis(-2π * dot(G_full, τ)) * 
+                             (θ ? conj(ψk[igired, iband]) : ψk[igired, iband])
         end
     end
 
@@ -172,7 +161,7 @@ end
 """
 Apply a symmetry operation to a density.
 """
-function apply_ksymop(symop::SymOp, basis, ρin)
+function apply_symop(symop::SymOp, basis, ρin)
     symop == one(SymOp) && return ρin
     symmetrize_ρ(basis, ρin; symmetries=[symop])
 end
@@ -201,7 +190,8 @@ function accumulate_over_symmetries!(ρaccu, ρin, basis, symmetries)
         for (ig, G) in enumerate(G_vectors_generator(basis.fft_size))
             igired = index_G_vectors(basis, invS * G)
             if igired !== nothing
-                @inbounds ρaccu[ig] += cis(-2T(π) * T(dot(G, symop.τ))) * ρin[igired]
+                @inbounds ρaccu[ig] += cis(-2T(π) * T(dot(G, symop.τ))) *
+                                       (symop.θ ? conj(ρin[igired]) : ρin[igired])
             end
         end
     end  # symop
@@ -285,10 +275,10 @@ This is mainly useful for debug purposes (e.g. in cases we don't want to
 bother thinking about symmetries).
 """
 function unfold_bz(basis::PlaneWaveBasis)
-    if all(length.(basis.ksymops_global) .== 1)
+    if length(basis.symmetries) == 1
         return basis
     else
-        kcoords = unfold_kcoords(basis.kcoords_global, basis.ksymops_global)
+        kcoords = unfold_kcoords(basis.kcoords_global, basis.symmetries)
         new_basis = PlaneWaveBasis(basis.model,
                                    basis.Ecut, basis.fft_size, basis.variational,
                                    kcoords, [[one(SymOp)] for _ in 1:length(kcoords)],
@@ -300,7 +290,7 @@ end
 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.ksymops[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)
@@ -326,7 +316,7 @@ function unfold_array_(basis_irred, basis_unfolded, data, 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_ksymop(symop, basis_irred,
+            _, ψSk = apply_symop(symop, basis_irred,
                                   basis_irred.kpoints[ik_irred], data[ik_irred])
             data_unfolded[ik_unfolded] = ψSk
         else
@@ -349,12 +339,9 @@ function unfold_bz(scfres)
     merge(scfres, new_scfres)
 end
 
-function unfold_kcoords(kcoords, ksymops)
-    all_kcoords = eltype(kcoords)[]
-    for ik = 1:length(kcoords)
-        for symop in ksymops[ik]
-            push!(all_kcoords, symop.S * kcoords[ik])
-        end
-    end
-    normalize_kpoint_coordinate.(all_kcoords)
+function unfold_kcoords(kcoords, symmetries)
+    all_kcoords = [symop.S * kcoord for kcoord in kcoords, symop in symmetries]
+    # the above multiplications introduce an error
+    unique(k -> normalize_kpoint_coordinate(round.(k; digits=ceil(Int, -log10(SYMMETRY_TOLERANCE)))),
+           normalize_kpoint_coordinate.(all_kcoords))
 end