We compare four different approaches for solving the DFT minimisation problem, namely a density-based SCF, a potential-based SCF, direct minimisation and Newton.
First we setup our problem
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
1.0e-6
scfres_scf = self_consistent_field(basis; tol);
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -7.846818949646 -0.70 4.5 2 -7.852317784126 -2.26 -1.53 1.0 36.6ms 3 -7.852617645060 -3.52 -2.56 1.5 40.2ms 4 -7.852645858657 -4.55 -2.89 2.5 47.8ms 5 -7.852646504983 -6.19 -3.21 1.2 170ms 6 -7.852646677529 -6.76 -4.23 1.0 35.2ms 7 -7.852646686666 -8.04 -4.95 2.2 49.7ms 8 -7.852646686723 -10.24 -5.44 1.2 38.2ms 9 -7.852646686729 -11.28 -5.71 1.2 37.8ms 10 -7.852646686730 -11.96 -6.60 1.0 37.8ms
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -7.846925223248 -0.70 4.8 2 -7.852529038765 -2.25 -1.63 0.80 2.0 346ms 3 -7.852636784512 -3.97 -2.72 0.80 1.0 32.2ms 4 -7.852646502333 -5.01 -3.29 0.80 2.2 43.6ms 5 -7.852646682143 -6.75 -4.13 0.80 1.2 34.4ms 6 -7.852646686366 -8.37 -4.83 0.80 1.2 33.7ms 7 -7.852646686722 -9.45 -5.89 0.80 2.0 42.1ms 8 -7.852646686730 -11.13 -6.67 0.80 2.0 41.5ms
Note: Unlike the other algorithms, tolerance for this one is in the energy, thus we square the density tolerance value to be roughly equivalent.
scfres_dm = direct_minimization(basis; tol=tol^2);
Iter Function value Gradient norm 0 1.381312e+01 3.058315e+00 * time: 0.5095341205596924 1 9.712065e-01 1.717580e+00 * time: 0.7737569808959961 2 -1.935801e+00 1.916885e+00 * time: 0.8076651096343994 3 -3.728783e+00 1.842075e+00 * time: 0.8577439785003662 4 -5.139869e+00 1.920588e+00 * time: 0.9129831790924072 5 -6.710954e+00 1.367151e+00 * time: 0.9677610397338867 6 -7.452315e+00 6.212700e-01 * time: 1.0199639797210693 7 -7.676187e+00 4.592504e-01 * time: 1.0547881126403809 8 -7.753233e+00 2.034475e-01 * time: 1.0879709720611572 9 -7.788749e+00 7.622030e-02 * time: 1.1221320629119873 10 -7.815776e+00 6.348562e-02 * time: 1.1550359725952148 11 -7.833394e+00 5.858491e-02 * time: 1.1883571147918701 12 -7.842064e+00 4.015647e-02 * time: 1.2197380065917969 13 -7.850047e+00 1.962371e-02 * time: 1.253767967224121 14 -7.851678e+00 1.101054e-02 * time: 1.2899491786956787 15 -7.852405e+00 6.464870e-03 * time: 1.3239161968231201 16 -7.852580e+00 5.886294e-03 * time: 1.3566551208496094 17 -7.852629e+00 1.622977e-03 * time: 1.390160083770752 18 -7.852642e+00 1.000632e-03 * time: 1.4240131378173828 19 -7.852645e+00 4.842165e-04 * time: 1.4560701847076416 20 -7.852646e+00 3.415662e-04 * time: 1.4911260604858398 21 -7.852646e+00 2.065074e-04 * time: 1.5242431163787842 22 -7.852647e+00 1.141188e-04 * time: 1.5592751502990723 23 -7.852647e+00 6.159923e-05 * time: 1.593641996383667 24 -7.852647e+00 5.006607e-05 * time: 1.6258301734924316 25 -7.852647e+00 2.300981e-05 * time: 1.657628059387207 26 -7.852647e+00 1.126302e-05 * time: 1.6911160945892334 27 -7.852647e+00 5.934221e-06 * time: 1.724276065826416 28 -7.852647e+00 3.654801e-06 * time: 1.7565810680389404 29 -7.852647e+00 2.198017e-06 * time: 1.7894861698150635 30 -7.852647e+00 1.361253e-06 * time: 1.8241140842437744 31 -7.852647e+00 6.561026e-07 * time: 1.8576619625091553 32 -7.852647e+00 3.630724e-07 * time: 1.891535997390747 33 -7.852647e+00 1.658644e-07 * time: 1.927828073501587 34 -7.852647e+00 1.058262e-07 * time: 1.9621570110321045 35 -7.852647e+00 5.879169e-08 * time: 1.9952220916748047 36 -7.852647e+00 4.118221e-08 * time: 2.033506155014038 37 -7.852647e+00 2.182735e-08 * time: 2.0687620639801025 38 -7.852647e+00 1.155588e-08 * time: 2.1025850772857666 39 -7.852647e+00 9.669505e-09 * time: 2.1523780822753906 40 -7.852647e+00 6.365953e-09 * time: 2.185908079147339 41 -7.852647e+00 6.365953e-09 * time: 2.368324041366577
Start not too far from the solution to ensure convergence: We run first a very crude SCF to get close and then switch to Newton.
scfres_start = self_consistent_field(basis; tol=0.5);
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -7.846757130857 -0.70 4.5
Remove the virtual orbitals (which Newton cannot treat yet)
ψ = DFTK.select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation).ψ
scfres_newton = newton(basis, ψ; tol);
n Energy log10(ΔE) log10(Δρ) Δtime --- --------------- --------- --------- ------ 1 -7.852645641253 -1.64 2 -7.852646686730 -5.98 -3.68 2.26s 3 -7.852646686730 -12.84 -7.06 328ms
println("|ρ_newton - ρ_scf| = ", norm(scfres_newton.ρ - scfres_scf.ρ))
println("|ρ_newton - ρ_scfv| = ", norm(scfres_newton.ρ - scfres_scfv.ρ))
println("|ρ_newton - ρ_dm| = ", norm(scfres_newton.ρ - scfres_dm.ρ))
|ρ_newton - ρ_scf| = 2.6589651775334955e-7 |ρ_newton - ρ_scfv| = 8.547977389425358e-8 |ρ_newton - ρ_dm| = 1.0808370690073526e-9