using DFTK
using ASEconvert

system = pyconvert(AbstractSystem, ase.build.bulk("Si"))
model  = model_LDA(attach_psp(system; Si="hgh/pbe/si-q4.hgh"))
basis  = PlaneWaveBasis(model; Ecut=5, kgrid=[3, 3, 3]);

using Plots
p = plot(yaxis=:log)
density_differences = Float64[];

using LinearAlgebra

function plot_callback(info)
    if info.stage == :finalize
        plot!(p, density_differences, label="|ρout - ρin|", markershape=:x)
    else
        push!(density_differences, norm(info.ρout - info.ρin))
    end
    info
end
callback = DFTK.ScfDefaultCallback() ∘ plot_callback;

scfres = self_consistent_field(basis; tol=1e-5, callback);

p