using DFTK
using LinearAlgebra
using ASEconvert

function aluminium_setup(repeat=1; Ecut=7.0, kgrid=[2, 2, 2])
    a = 7.65339
    lattice = a * Matrix(I, 3, 3)
    Al = ElementPsp(:Al, psp=load_psp("hgh/lda/al-q3"))
    atoms     = [Al, Al, Al, Al]
    positions = [[0.0, 0.0, 0.0], [0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]
    unit_cell = periodic_system(lattice, atoms, positions)

    # Make supercell in ASE:
    # We convert our lattice to the conventions used in ASE, make the supercell
    # and then convert back ...
    supercell_ase = convert_ase(unit_cell) * pytuple((repeat, 1, 1))
    supercell     = pyconvert(AbstractSystem, supercell_ase)

    # Unfortunately right now the conversion to ASE drops the pseudopotential information,
    # so we need to reattach it:
    supercell = attach_psp(supercell, Al="hgh/lda/al-q3")

    # Construct an LDA model and discretise
    # Note: We disable symmetries explicitly here. Otherwise the problem sizes
    #       we are able to run on the CI are too simple to observe the numerical
    #       instabilities we want to trigger here.
    model = model_LDA(supercell; temperature=1e-3, symmetries=false)
    PlaneWaveBasis(model; Ecut, kgrid)
end;

setup = aluminium_setup(5)
convert_ase(periodic_system(setup.model)).write("al_supercell.png")

self_consistent_field(aluminium_setup(1); tol=1e-4);

self_consistent_field(aluminium_setup(2); tol=1e-4);

self_consistent_field(aluminium_setup(4); tol=1e-4);

self_consistent_field(aluminium_setup(1); tol=1e-4, mixing=SimpleMixing());

self_consistent_field(aluminium_setup(4); tol=1e-4, mixing=SimpleMixing());