Custom solvers

In this example, we show how to define custom solvers. Our system will again be silicon, because we are not very imaginative

using DFTK, LinearAlgebra

a = 10.26
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]

# We take very (very) crude parameters
model = model_LDA(lattice, atoms, positions)
basis = PlaneWaveBasis(model; Ecut=5, kgrid=[1, 1, 1]);

We define our custom fix-point solver: simply a damped fixed-point

function my_fp_solver(f, x0, max_iter; tol)
    mixing_factor = .7
    x = x0
    fx = f(x)
    for n = 1:max_iter
        inc = fx - x
        if norm(inc) < tol
            break
        end
        x = x + mixing_factor * inc
        fx = f(x)
    end
    (fixpoint=x, converged=norm(fx-x) < tol)
end;

Our eigenvalue solver just forms the dense matrix and diagonalizes it explicitly (this only works for very small systems)

function my_eig_solver(A, X0; maxiter, tol, kwargs...)
    n = size(X0, 2)
    A = Array(A)
    E = eigen(A)
    λ = E.values[1:n]
    X = E.vectors[:, 1:n]
    (; λ, X, residual_norms=[], iterations=0, converged=true, n_matvec=0)
end;

Finally we also define our custom mixing scheme. It will be a mixture of simple mixing (for the first 2 steps) and than default to Kerker mixing. In the mixing interface δF is $(ρ_\text{out} - ρ_\text{in})$, i.e. the difference in density between two subsequent SCF steps and the mix function returns $δρ$, which is added to $ρ_\text{in}$ to yield $ρ_\text{next}$, the density for the next SCF step.

struct MyMixing
    n_simple  # Number of iterations for simple mixing
end
MyMixing() = MyMixing(2)

function DFTK.mix_density(mixing::MyMixing, basis, δF; n_iter, kwargs...)
    if n_iter <= mixing.n_simple
        return δF  # Simple mixing -> Do not modify update at all
    else
        # Use the default KerkerMixing from DFTK
        DFTK.mix_density(KerkerMixing(), basis, δF; kwargs...)
    end
end

That's it! Now we just run the SCF with these solvers

scfres = self_consistent_field(basis;
                               tol=1e-8,
                               solver=my_fp_solver,
                               eigensolver=my_eig_solver,
                               mixing=MyMixing());
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -7.235992833898                   -0.50    0.0
  2   -7.249608665538       -1.87       -0.92    0.0    628ms
  3   -7.251174940533       -2.81       -1.34    0.0   94.6ms
  4   -7.251295969966       -3.92       -1.65    0.0   92.4ms
  5   -7.251327319418       -4.50       -1.95    0.0   94.6ms
  6   -7.251335595391       -5.08       -2.25    0.0   94.8ms
  7   -7.251337861796       -5.64       -2.54    0.0   97.1ms
  8   -7.251338511556       -6.19       -2.83    0.0    103ms
  9   -7.251338706958       -6.71       -3.10    0.0    261ms
 10   -7.251338768376       -7.21       -3.36    0.0   94.3ms
 11   -7.251338788413       -7.70       -3.62    0.0   93.0ms
 12   -7.251338795145       -8.17       -3.87    0.0   96.8ms
 13   -7.251338797456       -8.64       -4.11    0.0   91.2ms
 14   -7.251338798263       -9.09       -4.35    0.0   93.0ms
 15   -7.251338798547       -9.55       -4.58    0.0   97.4ms
 16   -7.251338798648      -10.00       -4.81    0.0    106ms
 17   -7.251338798684      -10.44       -5.04    0.0    251ms
 18   -7.251338798697      -10.89       -5.27    0.0    105ms
 19   -7.251338798702      -11.33       -5.50    0.0   97.6ms
 20   -7.251338798704      -11.78       -5.72    0.0   94.3ms
 21   -7.251338798704      -12.22       -5.95    0.0   94.4ms
 22   -7.251338798704      -12.66       -6.17    0.0   97.8ms
 23   -7.251338798704      -13.11       -6.39    0.0    106ms
 24   -7.251338798705      -13.62       -6.62    0.0    247ms
 25   -7.251338798705      -13.97       -6.83    0.0    104ms
 26   -7.251338798705      -14.21       -7.05    0.0   95.0ms
 27   -7.251338798705      -15.05       -7.28    0.0   97.3ms
 28   -7.251338798705      -14.75       -7.42    0.0   95.4ms
 29   -7.251338798705   +  -15.05       -7.66    0.0   96.2ms
 30   -7.251338798705   +  -13.65       -7.01    0.0   97.7ms
 31   -7.251338798705      -13.75       -7.26    0.0    106ms
 32   -7.251338798705      -14.75       -7.60    0.0    279ms
 33   -7.251338798705   +    -Inf       -7.88    0.0    100ms
 34   -7.251338798705   +  -15.05       -7.42    0.0    105ms
 35   -7.251338798705      -15.05       -7.67    0.0    101ms
 36   -7.251338798705      -14.75       -7.96    0.0   98.1ms
 37   -7.251338798705   +  -14.75       -7.48    0.0   97.7ms
 38   -7.251338798705      -14.75       -7.73    0.0   97.4ms
 39   -7.251338798705   +  -14.27       -7.32    0.0    101ms
 40   -7.251338798705      -14.15       -7.53    0.0    233ms
 41   -7.251338798705   +  -14.27       -7.50    0.0    103ms
 42   -7.251338798705      -14.75       -7.74    0.0   98.8ms
 43   -7.251338798705      -14.45       -7.83    0.0    104ms
 44   -7.251338798705   +  -13.85       -7.15    0.0    103ms
 45   -7.251338798705      -13.97       -7.39    0.0   96.7ms
 46   -7.251338798705   +    -Inf       -7.71    0.0   94.9ms
 47   -7.251338798705      -14.45       -7.82    0.0    102ms
 48   -7.251338798704   +  -13.36       -6.81    0.0    226ms
 49   -7.251338798705      -13.44       -7.08    0.0   94.8ms
 50   -7.251338798705      -14.45       -7.43    0.0    105ms
 51   -7.251338798705      -15.05       -7.61    0.0   97.8ms
 52   -7.251338798705   +  -13.82       -7.08    0.0    105ms
 53   -7.251338798705      -13.75       -7.32    0.0   94.2ms
 54   -7.251338798705   +  -14.75       -7.62    0.0   99.0ms
 55   -7.251338798705   +  -14.75       -7.46    0.0    104ms
 56   -7.251338798705   +  -14.75       -7.24    0.0    227ms
 57   -7.251338798705      -14.35       -7.49    0.0    102ms
 58   -7.251338798705   +  -14.57       -7.51    0.0   99.2ms
 59   -7.251338798705      -14.57       -7.73    0.0   96.8ms
 60   -7.251338798705   +  -15.05       -7.84    0.0    101ms
 61   -7.251338798705   +    -Inf       -7.81    0.0   95.3ms
 62   -7.251338798705   +    -Inf       -7.48    0.0    109ms
 63   -7.251338798705   +  -14.75       -7.57    0.0    102ms
 64   -7.251338798705      -15.05       -7.47    0.0    221ms
 65   -7.251338798705   +  -15.05       -7.70    0.0    100ms
 66   -7.251338798705      -14.75       -7.60    0.0    100ms
 67   -7.251338798705   +  -15.05       -7.73    0.0   97.1ms
 68   -7.251338798705      -14.75       -7.84    0.0   99.1ms
 69   -7.251338798705   +  -14.35       -7.50    0.0   97.5ms
 70   -7.251338798705      -14.45       -7.42    0.0   99.1ms
 71   -7.251338798705   +  -14.75       -7.65    0.0    103ms
 72   -7.251338798705      -14.75       -7.86    0.0    231ms
 73   -7.251338798705   +  -14.10       -7.21    0.0    102ms
 74   -7.251338798705      -14.10       -7.44    0.0   99.3ms
 75   -7.251338798705   +    -Inf       -7.67    0.0   98.7ms
 76   -7.251338798705   +    -Inf       -7.71    0.0    102ms
 77   -7.251338798705   +  -14.75       -7.46    0.0    102ms
 78   -7.251338798705      -14.75       -7.74    0.0   98.7ms
 79   -7.251338798705   +  -14.45       -7.31    0.0    102ms
 80   -7.251338798705      -14.57       -7.52    0.0    228ms
 81   -7.251338798705      -15.05       -7.81    0.0   97.6ms
 82   -7.251338798705   +    -Inf       -7.90    0.0   98.7ms
 83   -7.251338798705   +  -15.05       -7.31    0.0   95.2ms
 84   -7.251338798705   +  -14.57       -7.48    0.0    106ms
 85   -7.251338798705   +    -Inf       -7.40    0.0   97.4ms
 86   -7.251338798705      -14.45       -7.60    0.0    105ms
 87   -7.251338798705      -14.57       -7.83    0.0    257ms
 88   -7.251338798705   +  -14.21       -7.33    0.0   96.5ms
 89   -7.251338798705      -14.75       -7.51    0.0   98.8ms
 90   -7.251338798705      -15.05       -7.65    0.0    108ms
 91   -7.251338798705   +  -14.75       -7.33    0.0   97.9ms
 92   -7.251338798705   +    -Inf       -7.45    0.0   97.6ms
 93   -7.251338798705      -14.35       -7.68    0.0    101ms
 94   -7.251338798705   +  -14.75       -7.87    0.0    106ms
 95   -7.251338798705   +  -13.80       -7.08    0.0    229ms
 96   -7.251338798705      -13.75       -7.32    0.0    104ms
 97   -7.251338798705   +    -Inf       -7.62    0.0    101ms
 98   -7.251338798705   +  -14.45       -7.45    0.0   97.9ms
 99   -7.251338798705      -14.57       -7.68    0.0   97.7ms
 100   -7.251338798705   +  -13.97       -7.16    0.0    100ms
 101   -7.251338798705      -14.01       -7.41    0.0   97.3ms
┌ Warning: SCF not converged.
└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/scf_callbacks.jl:38

Note that the default convergence criterion is the difference in density. When this gets below tol, the "driver" self_consistent_field artificially makes the fixed-point solver think it's converged by forcing f(x) = x. You can customize this with the is_converged keyword argument to self_consistent_field.