{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Practical error bounds for the forces\n",
"\n",
"DFTK includes an implementation of the strategy from [^CDKL2022] to compute\n",
"practical error bounds for forces and other quantities of interest.\n",
"\n",
"This is an example showing how to compute error estimates for the forces\n",
"on a ${\\rm TiO}_2$ molecule, from which we can either compute asymptotically\n",
"valid error bounds or increase the precision on the computation of the forces.\n",
"\n",
"[^CDKL2022]:\n",
" E. Cancès, G. Dusson, G. Kemlin, and A. Levitt\n",
" *Practical error bounds for properties in plane-wave electronic structure calculations*\n",
" [SIAM Journal on Scientific Computing 44 (5), B1312-B1340](https://doi.org/10.1137/21M1456224)"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using DFTK\n",
"using PseudoPotentialData\n",
"using Printf\n",
"using LinearAlgebra\n",
"using ForwardDiff"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"## Setup\n",
"We setup manually the ${\\rm TiO}_2$ configuration from\n",
"[Materials Project](https://materialsproject.org/materials/mp-2657/)."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"pseudopotentials = PseudoFamily(\"cp2k.nc.sr.lda.v0_1.largecore.gth\")\n",
"Ti = ElementPsp(:Ti, pseudopotentials)\n",
"O = ElementPsp(:O, pseudopotentials)\n",
"atoms = [Ti, Ti, O, O, O, O]\n",
"positions = [[0.5, 0.5, 0.5], # Ti\n",
" [0.0, 0.0, 0.0], # Ti\n",
" [0.19542, 0.80458, 0.5], # O\n",
" [0.80458, 0.19542, 0.5], # O\n",
" [0.30458, 0.30458, 0.0], # O\n",
" [0.69542, 0.69542, 0.0]] # O\n",
"lattice = [[8.79341 0.0 0.0];\n",
" [0.0 8.79341 0.0];\n",
" [0.0 0.0 5.61098]];"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"We apply a small displacement to one of the $\\rm Ti$ atoms to get nonzero\n",
"forces."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "3-element Vector{Float64}:\n 0.478\n 0.528\n 0.535"
},
"metadata": {},
"execution_count": 3
}
],
"cell_type": "code",
"source": [
"positions[1] .+= [-0.022, 0.028, 0.035]"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"We build a model with one $k$-point only, not too high `Ecut_ref` and small\n",
"tolerance to limit computational time. These parameters can be increased for\n",
"more precise results."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"model = model_DFT(lattice, atoms, positions; functionals=LDA())\n",
"kgrid = [1, 1, 1]\n",
"Ecut_ref = 35\n",
"basis_ref = PlaneWaveBasis(model; Ecut=Ecut_ref, kgrid)\n",
"tol = 1e-5;"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"We also build a basis with smaller `Ecut`, to compute a variational approximation of the reference solution.\n",
"\n",
"> **Choice of convergence parameters**\n",
">\n",
"> Be careful to choose `Ecut` not too close to `Ecut_ref`.\n",
"> Note also that the current choice `Ecut_ref = 35` is such that the\n",
"> reference solution is not converged and `Ecut = 15` is such that the\n",
"> asymptotic regime (crucial to validate the approach) is barely established."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"Ecut = 15\n",
"basis = PlaneWaveBasis(model; Ecut, kgrid);"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"## Computations\n",
"Compute the solution on the smaller basis:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"scfres = self_consistent_field(basis; tol, callback=identity);"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"Compute first order corrections `refinement.δψ` and `refinement.δρ`.\n",
"Note that `refinement.ψ` and `refinement.ρ` are the quantities computed with `Ecut`\n",
"and then extended to the reference grid.\n",
"This step is roughly as expensive as the `self_consistent_field` call above."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"refinement = refine_scfres(scfres, basis_ref; tol, callback=identity);"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"## Error estimates\n",
"- Computation of the force from the variational solution without any post-processing:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "6-element Vector{StaticArraysCore.SVector{3, Float64}}:\n [0.03418855950413535, -0.016276394305248854, 0.025767637977328173]\n [-0.5139216515045787, 0.6479201505659945, 0.2170681817193305]\n [-0.5238675354873454, 0.4801368623775526, -0.15971058202151486]\n [0.777613457067849, -0.8064496580824967, -0.12223539112573223]\n [0.6154633962454343, 0.3547562213960296, -0.0020629086600183477]\n [-0.3894877436373152, -0.6601071031231927, 0.04160786043271103]"
},
"metadata": {},
"execution_count": 8
}
],
"cell_type": "code",
"source": [
"f = compute_forces(scfres)"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"- Computation of the forces by a linearization argument when replacing the\n",
" error $P-P_*$ by the modified residual $R_{\\rm Schur}(P)$. The latter\n",
" quantity is computable in practice."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "6-element Vector{StaticArraysCore.SVector{3, Float64}}:\n [-0.21301722681824328, 0.2730113145648059, 0.29520434209420665]\n [-0.5681968786100761, 0.72046263411358, 0.2549407435641571]\n [-0.3622378313098461, 0.32855929296399394, -0.1598791007347069]\n [0.9488413396929462, -0.9818851370497104, -0.13380603475680183]\n [0.6787075799839407, 0.41328591303068946, -0.14153407483203972]\n [-0.48592754346491185, -0.7512568996296056, -0.11259874995984637]"
},
"metadata": {},
"execution_count": 9
}
],
"cell_type": "code",
"source": [
"force_refinement = refine_forces(refinement)\n",
"forces_refined = f + force_refinement.dF"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"A practical estimate of the error on the forces is then the following:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "6-element Vector{StaticArraysCore.SVector{3, Float64}}:\n [-0.24720578632237863, 0.2892877088700547, 0.2694367041168785]\n [-0.05427522710549737, 0.07254248354758552, 0.03787256184482657]\n [0.16162970417749933, -0.15157756941355865, -0.0001685187131920396]\n [0.17122788262509714, -0.17543547896721368, -0.011570643631069605]\n [0.06324418373850638, 0.058529691634659875, -0.13947116617202138]\n [-0.09643979982759665, -0.09114979650641297, -0.1542066103925574]"
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"dF_estimate = forces_refined - f"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"# Comparisons against non-practical estimates.\n",
"For practical computations one can stop at `forces_refined` and `dF_estimate`.\n",
"We continue here with a comparison of different ways to obtain the refined forces,\n",
"noting that the computational cost is much higher."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Computations\n",
"We compute the reference solution $P_*$ from which we will compute the\n",
"references forces."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"scfres_ref = self_consistent_field(basis_ref; tol, callback=identity)\n",
"ψ_ref = DFTK.select_occupied_orbitals(basis_ref, scfres_ref.ψ, scfres_ref.occupation).ψ;"
],
"metadata": {},
"execution_count": 11
},
{
"cell_type": "markdown",
"source": [
"- Compute the error $P-P_*$ on the associated orbitals $ϕ-ψ$ after aligning\n",
" them: this is done by solving $\\min |ϕ - ψU|$ for $U$ unitary matrix of\n",
" size $N×N$ ($N$ being the number of electrons) whose solution is\n",
" $U = S(S^*S)^{-1/2}$ where $S$ is the overlap matrix $ψ^*ϕ$."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function compute_error(ϕ, ψ)\n",
" map(zip(ϕ, ψ)) do (ϕk, ψk)\n",
" S = ψk'ϕk\n",
" U = S*(S'S)^(-1/2)\n",
" ϕk - ψk*U\n",
" end\n",
"end\n",
"error = compute_error(refinement.ψ, ψ_ref);"
],
"metadata": {},
"execution_count": 12
},
{
"cell_type": "markdown",
"source": [
"## Error estimates"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We start with different estimations of the forces:\n",
"- Force from the reference solution"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"f_ref = compute_forces(scfres_ref)\n",
"forces = Dict(\"F(P_*)\" => f_ref)\n",
"relerror = Dict(\"F(P_*)\" => 0.0)\n",
"compute_relerror(f) = norm(f - f_ref) / norm(f_ref);"
],
"metadata": {},
"execution_count": 13
},
{
"cell_type": "markdown",
"source": [
"- Force from the variational solution and relative error without\n",
" any post-processing:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"forces[\"F(P)\"] = f\n",
"\n",
"relerror[\"F(P)\"] = compute_relerror(f);"
],
"metadata": {},
"execution_count": 14
},
{
"cell_type": "markdown",
"source": [
"We then try to improve $F(P)$ using the first order linearization:\n",
"\n",
"$$\n",
"F(P) = F(P_*) + {\\rm d}F(P)·(P-P_*).\n",
"$$"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"To this end, we use the `ForwardDiff.jl` package to compute ${\\rm d}F(P)$\n",
"using automatic differentiation."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"function df(basis, occupation, ψ, δψ, ρ)\n",
" δρ = DFTK.compute_δρ(basis, ψ, δψ, occupation)\n",
" ForwardDiff.derivative(ε -> compute_forces(basis, ψ.+ε.*δψ, occupation; ρ=ρ+ε.*δρ), 0)\n",
"end;"
],
"metadata": {},
"execution_count": 15
},
{
"cell_type": "markdown",
"source": [
"- Computation of the forces by a linearization argument if we have access to\n",
" the actual error $P-P_*$. Usually this is of course not the case, but this\n",
" is the \"best\" improvement we can hope for with a linearisation, so we are\n",
" aiming for this precision."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"df_err = df(basis_ref, refinement.occupation, refinement.ψ,\n",
" DFTK.proj_tangent(error, refinement.ψ), refinement.ρ)\n",
"forces[\"F(P) - df(P)⋅(P-P_*)\"] = f - df_err\n",
"relerror[\"F(P) - df(P)⋅(P-P_*)\"] = compute_relerror(f - df_err);"
],
"metadata": {},
"execution_count": 16
},
{
"cell_type": "markdown",
"source": [
"- Computation of the forces by a linearization argument when replacing the\n",
" error $P-P_*$ by the modified residual $R_{\\rm Schur}(P)$. The latter\n",
" quantity is computable in practice."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"forces[\"F(P) - df(P)⋅Rschur(P)\"] = forces_refined\n",
"relerror[\"F(P) - df(P)⋅Rschur(P)\"] = compute_relerror(forces_refined);"
],
"metadata": {},
"execution_count": 17
},
{
"cell_type": "markdown",
"source": [
"Summary of all forces on the first atom (Ti)"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" F(P_*) = [-0.29996, 0.37533, 0.44083] (rel. error: 0.00000)\n",
" F(P) = [0.03419, -0.01628, 0.02577] (rel. error: 0.41550)\n",
" F(P) - df(P)⋅Rschur(P) = [-0.21302, 0.27301, 0.29520] (rel. error: 0.14681)\n",
" F(P) - df(P)⋅(P-P_*) = [-0.32017, 0.38917, 0.43183] (rel. error: 0.08583)\n"
]
}
],
"cell_type": "code",
"source": [
"for (key, value) in pairs(forces)\n",
" @printf(\"%30s = [%7.5f, %7.5f, %7.5f] (rel. error: %7.5f)\\n\",\n",
" key, (value[1])..., relerror[key])\n",
"end"
],
"metadata": {},
"execution_count": 18
},
{
"cell_type": "markdown",
"source": [
"Notice how close the computable expression $F(P) - {\\rm d}F(P)⋅R_{\\rm Schur}(P)$\n",
"is to the best linearization ansatz $F(P) - {\\rm d}F(P)⋅(P-P_*)$."
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.4"
},
"kernelspec": {
"name": "julia-1.11",
"display_name": "Julia 1.11.4",
"language": "julia"
}
},
"nbformat": 4
}