Notebook

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.

[^CDKL2022]: E. Cancès, G. Dusson, G. Kemlin, and A. Levitt Practical error bounds for properties in plane-wave electronic structure calculations SIAM Journal on Scientific Computing 44 (5), B1312-B1340

In [1]:

using DFTK
using PseudoPotentialData
using Printf
using LinearAlgebra
using ForwardDiff

Setup¶

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

In [2]:

pseudopotentials = PseudoFamily("cp2k.nc.sr.lda.v0_1.largecore.gth")
Ti = ElementPsp(:Ti, pseudopotentials)
O  = ElementPsp(:O,  pseudopotentials)
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.

In [3]:

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

Out[3]:

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.

In [4]:

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.

Choice of convergence parameters

Be careful to choose Ecut not too close to Ecut_ref. Note also that the current choice Ecut_ref = 35 is such that the reference solution is not converged and Ecut = 15 is such that the asymptotic regime (crucial to validate the approach) is barely established.

In [5]:

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

Computations¶

Compute the solution on the smaller basis:

In [6]:

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.

In [7]:

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

Error estimates¶

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

In [8]:

f = compute_forces(scfres)

Out[8]:

6-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [0.03418855950413535, -0.016276394305248854, 0.025767637977328173]
 [-0.5139216515045787, 0.6479201505659945, 0.2170681817193305]
 [-0.5238675354873454, 0.4801368623775526, -0.15971058202151486]
 [0.777613457067849, -0.8064496580824967, -0.12223539112573223]
 [0.6154633962454343, 0.3547562213960296, -0.0020629086600183477]
 [-0.3894877436373152, -0.6601071031231927, 0.04160786043271103]

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.

In [9]:

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

Out[9]:

6-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [-0.21301722681824328, 0.2730113145648059, 0.29520434209420665]
 [-0.5681968786100761, 0.72046263411358, 0.2549407435641571]
 [-0.3622378313098461, 0.32855929296399394, -0.1598791007347069]
 [0.9488413396929462, -0.9818851370497104, -0.13380603475680183]
 [0.6787075799839407, 0.41328591303068946, -0.14153407483203972]
 [-0.48592754346491185, -0.7512568996296056, -0.11259874995984637]

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

In [10]:

dF_estimate = forces_refined - f

Out[10]:

6-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [-0.24720578632237863, 0.2892877088700547, 0.2694367041168785]
 [-0.05427522710549737, 0.07254248354758552, 0.03787256184482657]
 [0.16162970417749933, -0.15157756941355865, -0.0001685187131920396]
 [0.17122788262509714, -0.17543547896721368, -0.011570643631069605]
 [0.06324418373850638, 0.058529691634659875, -0.13947116617202138]
 [-0.09643979982759665, -0.09114979650641297, -0.1542066103925574]

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.

In [11]:

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 $ψ^*ϕ$ .

In [12]:

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

In [13]:

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:

In [14]:

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.

In [15]:

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.

In [16]:

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.

In [17]:

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)

In [18]:

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_*) = [-0.29996, 0.37533, 0.44083]   (rel. error: 0.00000)
                          F(P) = [0.03419, -0.01628, 0.02577]   (rel. error: 0.41550)
        F(P) - df(P)⋅Rschur(P) = [-0.21302, 0.27301, 0.29520]   (rel. error: 0.14681)
          F(P) - df(P)⋅(P-P_*) = [-0.32017, 0.38917, 0.43183]   (rel. error: 0.08583)

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_*)$ .