using DFTK

Si = ElementCohenBergstresser(:Si)
atoms = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]
lattice = Si.lattice_constant / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]];

model = Model(lattice, atoms, positions; terms=[Kinetic(), AtomicLocal()])
basis = PlaneWaveBasis(model, Ecut=10.0, kgrid=(2, 2, 2));

ham = Hamiltonian(basis)
eigres = diagonalize_all_kblocks(DFTK.lobpcg_hyper, ham, 6)
εF = DFTK.compute_occupation(basis, eigres.λ).εF

using Plots
using Unitful

p = plot_bandstructure(basis; n_bands=8, εF, kline_density=10, unit=u"eV")
ylims!(p, (-5, 6))