using DFTK

a = 10.263141334305942  # 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];

pystruct = pymatgen_structure(lattice, atoms, positions)
pystruct.make_supercell([2, 2, 2])
lattice   = load_lattice(pystruct)
positions = load_positions(pystruct)
atoms     = fill(Si, length(positions));

model = model_LDA(lattice, atoms, positions)
basis = PlaneWaveBasis(model; Ecut=5, kgrid=(1, 1, 1))

scfres = direct_minimization(basis; tol=1e-3);
ψ, _ = DFTK.select_occupied_orbitals(basis, scfres.ψ, scfres.occupation)
scfres_newton = newton(basis, ψ; tol=1e-12)

scfres_newton.energies