Practical error bounds for the forces

DFTK includes an implementation of the strategy from ^[CDKL2022] to compute practical error bounds for forces and other quantities of interest.

This is an example showing how to compute error estimates for the forces on a ${\rm TiO}_2$ molecule, from which we can either compute asymptotically valid error bounds or increase the precision on the computation of the forces.

using DFTK
using Printf
using LinearAlgebra
using ForwardDiff

Setup

We setup manually the ${\rm TiO}_2$ configuration from Materials Project.

Ti = ElementPsp(:Ti, load_psp("hgh/lda/ti-q4.hgh"))
O  = ElementPsp(:O,  load_psp("hgh/lda/o-q6.hgh"))
atoms     = [Ti, Ti, O, O, O, O]
positions = [[0.5,     0.5,     0.5],  # Ti
             [0.0,     0.0,     0.0],  # Ti
             [0.19542, 0.80458, 0.5],  # O
             [0.80458, 0.19542, 0.5],  # O
             [0.30458, 0.30458, 0.0],  # O
             [0.69542, 0.69542, 0.0]]  # O
lattice   = [[8.79341  0.0      0.0];
             [0.0      8.79341  0.0];
             [0.0      0.0      5.61098]];

We apply a small displacement to one of the $\rm Ti$ atoms to get nonzero forces.

positions[1] .+= [-0.022, 0.028, 0.035]

3-element Vector{Float64}:
 0.478
 0.528
 0.535

We build a model with one $k$-point only, not too high Ecut_ref and small tolerance to limit computational time. These parameters can be increased for more precise results.

model = model_DFT(lattice, atoms, positions; functionals=LDA())
kgrid = [1, 1, 1]
Ecut_ref = 35
basis_ref = PlaneWaveBasis(model; Ecut=Ecut_ref, kgrid)
tol = 1e-5;

We also build a basis with smaller Ecut, to compute a variational approximation of the reference solution.

Ecut = 15
basis = PlaneWaveBasis(model; Ecut, kgrid);

Computations

Compute the solution on the smaller basis:

scfres = self_consistent_field(basis; tol, callback=identity);

Compute first order corrections refinement.δψ and refinement.δρ. Note that refinement.ψ and refinement.ρ are the quantities computed with Ecut and then extended to the reference grid. This step is roughly as expensive as the self_consistent_field call above.

refinement = refine_scfres(scfres, basis_ref; tol, callback=identity);

Error estimates

• Computation of the force from the variational solution without any post-processing:

f = compute_forces(scfres)

6-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [1.1355158087730388, -1.0152644365803742, 0.4001440405481873]
 [-0.5925123151954181, 0.8447323808567784, 0.46911993739625435]
 [-3.6627416238652746, 1.4344571509007669, -0.4095702055901477]
 [0.9250875188555172, 0.9711233982001515, -0.15000492263121368]
 [2.6531318898063825, 0.3993814404971774, -0.12657202389928535]
 [-0.45850443421553777, -2.634463072279953, -0.18270204499087583]

• Computation of the forces by a linearization argument when replacing the error $P-P_*$ by the modified residual $R_{\rm Schur}(P)$. The latter quantity is computable in practice.

force_refinement = refine_forces(refinement)
forces_refined = f + force_refinement.dF

6-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [1.2910450864872183, -1.1015921837641538, 0.6908021516912977]
 [-0.40995741882499254, 0.710968949769991, 0.5337610009082603]
 [-4.203159970920883, 1.625850324428251, -0.602708447644453]
 [0.8244264124763128, 1.2387832391065061, -0.34575956373124866]
 [3.01557295523826, 0.4905928527698006, -0.0943018689773653]
 [-0.5197570061664993, -2.962425215984985, -0.17946669823754965]

A practical estimate of the error on the forces is then the following:

dF_estimate = forces_refined - f

6-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [0.15552927771417946, -0.08632774718377956, 0.29065811114311035]
 [0.1825548963704256, -0.13376343108678734, 0.06464106351200594]
 [-0.5404183470556081, 0.19139317352748408, -0.1931382420543053]
 [-0.10066110637920445, 0.26765984090635464, -0.195754641100035]
 [0.3624410654318777, 0.0912114122726232, 0.03227015492192005]
 [-0.061252571950961565, -0.3279621437050322, 0.0032353467533261793]

Comparisons against non-practical estimates.

For practical computations one can stop at forces_refined and dF_estimate. We continue here with a comparison of different ways to obtain the refined forces, noting that the computational cost is much higher.

Computations

We compute the reference solution $P_*$ from which we will compute the references forces.

scfres_ref = self_consistent_field(basis_ref; tol, callback=identity)
ψ_ref = DFTK.select_occupied_orbitals(basis_ref, scfres_ref.ψ, scfres_ref.occupation).ψ;

• Compute the error $P-P_*$ on the associated orbitals $ϕ-ψ$ after aligning them: this is done by solving $\min |ϕ - ψU|$ for $U$ unitary matrix of size $N×N$ ($N$ being the number of electrons) whose solution is $U = S(S^*S)^{-1/2}$ where $S$ is the overlap matrix $ψ^*ϕ$.

function compute_error(ϕ, ψ)
    map(zip(ϕ, ψ)) do (ϕk, ψk)
        S = ψk'ϕk
        U = S*(S'S)^(-1/2)
        ϕk - ψk*U
    end
end
error = compute_error(refinement.ψ, ψ_ref);

Error estimates

We start with different estimations of the forces:

• Force from the reference solution

f_ref = compute_forces(scfres_ref)
forces   = Dict("F(P_*)" => f_ref)
relerror = Dict("F(P_*)" => 0.0)
compute_relerror(f) = norm(f - f_ref) / norm(f_ref);

• Force from the variational solution and relative error without any post-processing:

forces["F(P)"]   = f

relerror["F(P)"] = compute_relerror(f);

We then try to improve $F(P)$ using the first order linearization:

\[F(P) = F(P_*) + {\rm d}F(P)·(P-P_*).\]

To this end, we use the ForwardDiff.jl package to compute ${\rm d}F(P)$ using automatic differentiation.

function df(basis, occupation, ψ, δψ, ρ)
    δρ = DFTK.compute_δρ(basis, ψ, δψ, occupation)
    ForwardDiff.derivative(ε -> compute_forces(basis, ψ.+ε.*δψ, occupation; ρ=ρ+ε.*δρ), 0)
end;

• Computation of the forces by a linearization argument if we have access to the actual error $P-P_*$. Usually this is of course not the case, but this is the "best" improvement we can hope for with a linearisation, so we are aiming for this precision.

df_err = df(basis_ref, refinement.occupation, refinement.ψ,
            DFTK.proj_tangent(error, refinement.ψ), refinement.ρ)
forces["F(P) - df(P)⋅(P-P_*)"]   = f - df_err
relerror["F(P) - df(P)⋅(P-P_*)"] = compute_relerror(f - df_err);

• Computation of the forces by a linearization argument when replacing the error $P-P_*$ by the modified residual $R_{\rm Schur}(P)$. The latter quantity is computable in practice.

forces["F(P) - df(P)⋅Rschur(P)"]   = forces_refined
relerror["F(P) - df(P)⋅Rschur(P)"] = compute_relerror(forces_refined);

Summary of all forces on the first atom (Ti)

for (key, value) in pairs(forces)
    @printf("%30s = [%7.5f, %7.5f, %7.5f]   (rel. error: %7.5f)\n",
            key, (value[1])..., relerror[key])
end

                        F(P_*) = [1.47896, -1.25371, 0.81011]   (rel. error: 0.00000)
                          F(P) = [1.13552, -1.01526, 0.40014]   (rel. error: 0.20481)
        F(P) - df(P)⋅Rschur(P) = [1.29105, -1.10159, 0.69080]   (rel. error: 0.07836)
          F(P) - df(P)⋅(P-P_*) = [1.50903, -1.28630, 0.86147]   (rel. error: 0.08073)

Notice how close the computable expression $F(P) - {\rm d}F(P)⋅R_{\rm Schur}(P)$ is to the best linearization ansatz $F(P) - {\rm d}F(P)⋅(P-P_*)$.