{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Comparison of DFT solvers"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We compare four different approaches for solving the DFT minimisation problem,\n",
"namely a density-based SCF, a potential-based SCF, direct minimisation and Newton."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"First we setup our problem"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "1.0e-6"
},
"metadata": {},
"execution_count": 1
}
],
"cell_type": "code",
"source": [
"using AtomsBuilder\n",
"using DFTK\n",
"using LinearAlgebra\n",
"using PseudoPotentialData\n",
"\n",
"pseudopotentials = PseudoFamily(\"dojo.nc.sr.pbesol.v0_4_1.standard.upf\")\n",
"model = model_DFT(bulk(:Si); functionals=PBEsol(), pseudopotentials)\n",
"basis = PlaneWaveBasis(model; Ecut=5, kgrid=[3, 3, 3])\n",
"\n",
"# Convergence we desire in the density\n",
"tol = 1e-6"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"## Density-based self-consistent field"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy log10(ΔE) log10(Δρ) Diag Δtime\n",
"--- --------------- --------- --------- ---- ------\n",
" 1 -8.397850398111 -0.90 5.2 27.0ms\n",
" 2 -8.400256071796 -2.62 -1.74 1.0 18.7ms\n",
" 3 -8.400406663600 -3.82 -2.96 1.5 19.8ms\n",
" 4 -8.400427833193 -4.67 -2.98 3.0 57.9ms\n",
" 5 -8.400427949655 -6.93 -3.05 1.0 18.9ms\n",
" 6 -8.400428149316 -6.70 -4.70 1.0 19.0ms\n",
" 7 -8.400428155862 -8.18 -4.45 3.2 26.0ms\n",
" 8 -8.400428156256 -9.40 -5.21 1.0 19.6ms\n",
" 9 -8.400428156275 -10.72 -6.13 1.0 19.7ms\n"
]
}
],
"cell_type": "code",
"source": [
"scfres_scf = self_consistent_field(basis; tol);"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"## Potential-based SCF"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -8.397726925764 -0.90 4.8 553ms\n",
" 2 -8.400374639656 -2.58 -1.77 0.80 2.0 495ms\n",
" 3 -8.400422170514 -4.32 -3.00 0.80 1.0 180ms\n",
" 4 -8.400428104597 -5.23 -3.46 0.80 2.2 19.9ms\n",
" 5 -8.400428151685 -7.33 -4.97 0.80 1.2 17.8ms\n",
" 6 -8.400428156269 -8.34 -5.35 0.80 3.0 23.0ms\n",
" 7 -8.400428156276 -11.12 -6.86 0.80 1.0 16.4ms\n"
]
}
],
"cell_type": "code",
"source": [
"scfres_scfv = DFTK.scf_potential_mixing(basis; tol);"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"## Direct minimization"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy log10(ΔE) log10(Δρ) Δtime\n",
"--- --------------- --------- --------- ------\n",
" 1 +0.926488639176 -1.09 2.77s\n",
" 2 -1.848359491024 0.44 -0.65 106ms\n",
" 3 -4.582656986494 0.44 -0.39 39.3ms\n",
" 4 -6.169239948250 0.20 -0.47 73.3ms\n",
" 5 -7.680944399394 0.18 -0.73 39.6ms\n",
" 6 -7.933004251426 -0.60 -1.45 29.2ms\n",
" 7 -8.218204309711 -0.54 -1.62 29.0ms\n",
" 8 -8.227933182068 -2.01 -1.94 28.8ms\n",
" 9 -8.275254161630 -1.32 -1.79 49.0ms\n",
" 10 -8.319166652806 -1.36 -1.67 39.4ms\n",
" 11 -8.353780606343 -1.46 -1.75 39.5ms\n",
" 12 -8.374376330935 -1.69 -2.43 29.0ms\n",
" 13 -8.389239941286 -1.83 -2.40 33.4ms\n",
" 14 -8.396952636531 -2.11 -2.50 29.1ms\n",
" 15 -8.398593985646 -2.78 -2.59 29.0ms\n",
" 16 -8.399631212256 -2.98 -3.11 29.0ms\n",
" 17 -8.400022884144 -3.41 -3.11 32.3ms\n",
" 18 -8.400262490188 -3.62 -3.37 29.0ms\n",
" 19 -8.400353215508 -4.04 -3.59 29.1ms\n",
" 20 -8.400407239636 -4.27 -4.05 32.2ms\n",
" 21 -8.400417995679 -4.97 -3.82 28.9ms\n",
" 22 -8.400425286419 -5.14 -4.36 29.0ms\n",
" 23 -8.400426555823 -5.90 -4.17 32.1ms\n",
" 24 -8.400427343278 -6.10 -4.71 29.1ms\n",
" 25 -8.400427745067 -6.40 -4.57 32.0ms\n",
" 26 -8.400427971511 -6.65 -5.01 29.0ms\n",
" 27 -8.400428076202 -6.98 -5.06 29.5ms\n",
" 28 -8.400428128024 -7.29 -5.21 33.7ms\n",
" 29 -8.400428141728 -7.86 -5.46 29.0ms\n",
" 30 -8.400428151034 -8.03 -5.36 29.0ms\n",
" 31 -8.400428154120 -8.51 -5.67 32.2ms\n",
" 32 -8.400428155458 -8.87 -5.85 29.0ms\n",
" 33 -8.400428155871 -9.38 -6.61 32.2ms\n"
]
}
],
"cell_type": "code",
"source": [
"scfres_dm = direct_minimization(basis; tol);"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"## Newton algorithm"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Start not too far from the solution to ensure convergence:\n",
"We run first a very crude SCF to get close and then switch to Newton."
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy log10(ΔE) log10(Δρ) Diag Δtime\n",
"--- --------------- --------- --------- ---- ------\n",
" 1 -8.397856585042 -0.90 5.0 26.7ms\n"
]
}
],
"cell_type": "code",
"source": [
"scfres_start = self_consistent_field(basis; tol=0.5);"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"Remove the virtual orbitals (which Newton cannot treat yet)"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy log10(ΔE) log10(Δρ) Δtime\n",
"--- --------------- --------- --------- ------\n",
" 1 -8.400427983085 -1.78 24.2s\n",
" 2 -8.400428156277 -6.76 -4.02 2.19s\n",
" 3 -8.400428156277 -14.75 -7.82 115ms\n"
]
}
],
"cell_type": "code",
"source": [
"ψ = DFTK.select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation).ψ\n",
"scfres_newton = newton(basis, ψ; tol);"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"## Comparison of results"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"|ρ_newton - ρ_scf| = 1.318651071213495e-6\n",
"|ρ_newton - ρ_scfv| = 1.2366889710809778e-7\n",
"|ρ_newton - ρ_dm| = 1.5628380220615619e-6\n"
]
}
],
"cell_type": "code",
"source": [
"println(\"|ρ_newton - ρ_scf| = \", norm(scfres_newton.ρ - scfres_scf.ρ))\n",
"println(\"|ρ_newton - ρ_scfv| = \", norm(scfres_newton.ρ - scfres_scfv.ρ))\n",
"println(\"|ρ_newton - ρ_dm| = \", norm(scfres_newton.ρ - scfres_dm.ρ))"
],
"metadata": {},
"execution_count": 7
}
],
"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
}