diff --git a/src/DFTK.jl b/src/DFTK.jl index 05f8aea724..83782204a7 100644 --- a/src/DFTK.jl +++ b/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 diff --git a/src/Model.jl b/src/Model.jl index b8218b3059..5211293c33 100644 --- a/src/Model.jl +++ b/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 diff --git a/src/PlaneWaveBasis.jl b/src/PlaneWaveBasis.jl index 2ee83d69bc..827a187680 100644 --- a/src/PlaneWaveBasis.jl +++ b/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, diff --git a/src/SymOp.jl b/src/SymOp.jl index 55fb9b3296..2c0f875295 100644 --- a/src/SymOp.jl +++ b/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) diff --git a/src/bzmesh.jl b/src/bzmesh.jl index a386b40817..81c1c558cd 100644 --- a/src/bzmesh.jl +++ b/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.") diff --git a/src/postprocess/current.jl b/src/postprocess/current.jl index 499d6213ef..c941ac6550 100644 --- a/src/postprocess/current.jl +++ b/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])) diff --git a/src/postprocess/dos.jl b/src/postprocess/dos.jl index 638ff8e8e0..99fb573d5e 100644 --- a/src/postprocess/dos.jl +++ b/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 ε """ diff --git a/src/symmetry.jl b/src/symmetry.jl index a1573f087e..a135f2fd31 100644 --- a/src/symmetry.jl +++ b/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 diff --git a/src/terms/terms.jl b/src/terms/terms.jl index 79c7aea0a4..04e5a3cf60 100644 --- a/src/terms/terms.jl +++ b/src/terms/terms.jl @@ -32,9 +32,11 @@ end include("Hamiltonian.jl") # breaks_symmetries on a term builder answers true if this term breaks -# the symmetries of the lattice/atoms (in which case k-point reduction -# is invalid) +# the geometrical symmetries of the lattice/atoms (in which case k-point +# reduction is invalid) breaks_symmetries(::Any) = false +# answers true if the term breaks the time-reversal symmetry +breaks_time_reversal(::Any) = false include("kinetic.jl") @@ -53,9 +55,11 @@ include("pairwise.jl") include("magnetic.jl") breaks_symmetries(::Magnetic) = true +breaks_time_reversal(::Magnetic) = true include("anyonic.jl") breaks_symmetries(::Anyonic) = true +breaks_time_reversal(::Anyonic) = true # forces computes either nothing or an array forces[el][at][α] compute_forces(::Term, ::AbstractBasis, ψ, occ; kwargs...) = nothing # by default, no force diff --git a/src/workarounds/intervals.jl b/src/workarounds/intervals.jl index 435034fcb3..000c20b1b4 100644 --- a/src/workarounds/intervals.jl +++ b/src/workarounds/intervals.jl @@ -27,10 +27,10 @@ function _is_well_conditioned(A::AbstractArray{<:Interval}; kwargs...) end function symmetry_operations(lattice::AbstractMatrix{<:Interval}, atoms, magnetic_moments=[]; - tol_symmetry=max(SYMMETRY_TOLERANCE, maximum(radius, lattice))) + is_time_reversal=false, tol_symmetry=max(SYMMETRY_TOLERANCE, maximum(radius, lattice))) @assert tol_symmetry < 1e-2 symmetry_operations(IntervalArithmetic.mid.(lattice), atoms, magnetic_moments; - tol_symmetry) + is_time_reversal, tol_symmetry) end function local_potential_fourier(el::ElementCohenBergstresser, q::T) where {T <: Interval} diff --git a/test/chi0.jl b/test/chi0.jl index 45b7b400ee..bbaed2a756 100644 --- a/test/chi0.jl +++ b/test/chi0.jl @@ -73,7 +73,6 @@ function test_chi0(testcase; symmetry=false, temperature=0, # just to cover it here if temperature > 0 - N = compute_nos(εF, basis, res.eigenvalues) D = compute_dos(εF, basis, res.eigenvalues) LDOS = compute_ldos(εF, basis, res.eigenvalues, res.ψ) end diff --git a/test/compute_density.jl b/test/compute_density.jl index 8396134f75..9dea53ebe0 100644 --- a/test/compute_density.jl +++ b/test/compute_density.jl @@ -107,7 +107,7 @@ if mpi_nprocs() == 1 # not easy to distribute for (ik, k) in enumerate(kcoords) Hk_ir = ham_ir.blocks[ik] for symop in ksymops[ik] - Skpoint, ψSk = DFTK.apply_ksymop(symop, ham_ir.basis, Hk_ir.kpoint, ψ_ir[ik]) + Skpoint, ψSk = DFTK.apply_symop(symop, ham_ir.basis, Hk_ir.kpoint, ψ_ir[ik]) ikfull = findfirst(1:length(kfull)) do idx all(isinteger, round.(kfull[idx] - Skpoint.coordinate, digits=10))