using DFTK using LinearAlgebra 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) basis = PlaneWaveBasis(model; Ecut=5, kgrid=[3, 3, 3]) # Convergence we desire in the density tol = 1e-6 scfres_scf = self_consistent_field(basis; tol); scfres_scfv = DFTK.scf_potential_mixing(basis; tol); scfres_dm = direct_minimization(basis; tol=tol^2); scfres_start = self_consistent_field(basis; tol=0.5); ψ = DFTK.select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation).ψ scfres_newton = newton(basis, ψ; tol); println("|ρ_newton - ρ_scf| = ", norm(scfres_newton.ρ - scfres_scf.ρ)) println("|ρ_newton - ρ_scfv| = ", norm(scfres_newton.ρ - scfres_scfv.ρ)) println("|ρ_newton - ρ_dm| = ", norm(scfres_newton.ρ - scfres_dm.ρ))