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.231541177875 -0.48 0.0 2 -7.249024317152 -1.76 -0.89 0.0 637ms 3 -7.251138393957 -2.67 -1.31 0.0 84.4ms 4 -7.251287720328 -3.83 -1.62 0.0 84.8ms 5 -7.251325440039 -4.42 -1.92 0.0 87.9ms 6 -7.251335154306 -5.01 -2.22 0.0 209ms 7 -7.251337753935 -5.59 -2.51 0.0 84.1ms 8 -7.251338483911 -6.14 -2.79 0.0 85.3ms 9 -7.251338699518 -6.67 -3.06 0.0 85.1ms 10 -7.251338766276 -7.18 -3.33 0.0 84.9ms 11 -7.251338787794 -7.67 -3.59 0.0 85.1ms 12 -7.251338794955 -8.15 -3.84 0.0 90.4ms 13 -7.251338797396 -8.61 -4.09 0.0 219ms 14 -7.251338798243 -9.07 -4.33 0.0 81.3ms 15 -7.251338798541 -9.53 -4.56 0.0 84.0ms 16 -7.251338798646 -9.98 -4.80 0.0 83.6ms 17 -7.251338798684 -10.43 -5.03 0.0 85.6ms 18 -7.251338798697 -10.87 -5.26 0.0 84.8ms 19 -7.251338798702 -11.32 -5.48 0.0 88.2ms 20 -7.251338798704 -11.76 -5.71 0.0 90.9ms 21 -7.251338798704 -12.21 -5.93 0.0 211ms 22 -7.251338798704 -12.65 -6.16 0.0 85.3ms 23 -7.251338798704 -13.07 -6.38 0.0 85.2ms 24 -7.251338798705 -13.60 -6.60 0.0 85.6ms 25 -7.251338798705 -13.88 -6.82 0.0 85.2ms 26 -7.251338798705 + -14.75 -7.02 0.0 84.2ms 27 -7.251338798705 -14.27 -7.23 0.0 90.4ms 28 -7.251338798705 + -14.75 -7.34 0.0 201ms 29 -7.251338798705 -15.05 -7.54 0.0 84.4ms 30 -7.251338798705 + -14.75 -7.35 0.0 85.3ms 31 -7.251338798705 -14.45 -7.57 0.0 84.1ms 32 -7.251338798705 -14.75 -7.83 0.0 85.1ms 33 -7.251338798705 + -14.75 -7.64 0.0 84.3ms 34 -7.251338798705 + -14.75 -7.73 0.0 90.1ms 35 -7.251338798705 -14.45 -7.64 0.0 88.6ms 36 -7.251338798705 + -14.45 -7.48 0.0 190ms 37 -7.251338798705 -14.75 -7.67 0.0 81.5ms 38 -7.251338798705 + -Inf -7.93 0.0 85.1ms 39 -7.251338798705 + -Inf -7.78 0.0 82.9ms 40 -7.251338798705 + -14.75 -7.52 0.0 84.7ms 41 -7.251338798705 + -Inf -7.78 0.0 85.9ms 42 -7.251338798705 -14.45 -7.60 0.0 91.4ms 43 -7.251338798705 + -Inf -7.69 0.0 197ms 44 -7.251338798705 + -Inf -8.01 0.0 83.6ms
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
.