unexpected illegal instruction, good case

illegal_instruction_DFTK.jl

using Brillouin
using Bravais
using DFTK
using Unitful
using UnitfulAtomic
using StaticArrays

a = 10.26  # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms     = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]

model = model_LDA(lattice, atoms, positions)
kline_density = austrip(40u"bohr")

if model.n_dim == 1  # Return fast for 1D model
    # TODO Is this special-casing of 1D is not needed for Brillouin.jl any more
    #
    # Just use irrfbz_path(1, DirectBasis{1}([1.0]))
    # (see https://github.com/JuliaMolSim/DFTK.jl/pull/496/files#r725205860)
    #
    # Length of the kpath is recip_lattice[1, 1] in 1D
    n_points = max(2, 1 + ceil(Int, kline_density * model.recip_lattice[1, 1]))
    kcoords  = [@SVector[coord, 0, 0] for coord in range(-1//2, 1//2, length=n_points)]
    klabels  = Dict("Γ" => zeros(3), "-½" => [-0.5, 0.0, 0.0], "½" => [0.5, 0, 0])
    return (kcoords=kcoords, klabels=klabels,
            kpath=[["-½", "½"]], kbranches=[1:length(kcoords)])
end

# - Brillouin.jl expects the input direct lattice to be in the conventional lattice
#   in the convention of the International Table of Crystallography Vol A (ITA).
# - spglib uses this convention for the returned conventional lattice,
#   so it can be directly used as input to Brillouin.jl
# - The output k-Points and reciprocal lattices will be in the CDML convention.
conv_latt = DFTK.spglib_standardize_cell(model; primitive=false,correct_symmetry=false).lattice
sgnum     = DFTK.spglib_spacegroup_number(model)  # Get ITA space-group number
direct_basis   = Bravais.DirectBasis(collect(eachcol(conv_latt)))
primitive_latt = Bravais.primitivize(direct_basis, Bravais.centering(sgnum, 3))

primitive_latt ≈ collect(eachcol(model.lattice)) || @warn(
    "DFTK's model.lattice and Brillouin's primitive lattice do not agree. " *
    "The kpath selected to plot the band structure might not be most appropriate."
)

kp     = Brillouin.irrfbz_path(sgnum, direct_basis)
kinter = Brillouin.interpolate(kp, density=kline_density)

# TODO Need to take care of time-reversal symmetry here!
#      See https://github.com/JuliaMolSim/DFTK.jl/pull/496/files#r725203554

# Need to double the points whenever a new path starts
# (for temporary compatibility with pymatgen)
# TODO Remove this later
kcoords = empty(first(kinter.kpaths))
for kbranch in kinter.kpaths
    idcs = findall(k -> any(sum(abs2, k - kcomp) < 1e-5
                            for kcomp in values(kp.points)), kbranch)
    @assert length(idcs) ≥ 2
    idcs = idcs[2:end-1]  # Don't duplicate first and last
    idcs = sort(append!(idcs, 1:length(kbranch)))
    append!(kcoords, kbranch[idcs])
end

T = eltype(kcoords[1])
klabels = Dict{String,Vector{T}}(string(key) => val for (key, val) in kp.points)
kpath   = [[string(el) for el in path] for path in kp.paths]
(; kcoords, klabels, kpath)