{ "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 DFTK\n", "using LinearAlgebra\n", "\n", "a = 10.26 # Silicon lattice constant in Bohr\n", "lattice = a / 2 * [[0 1 1.];\n", " [1 0 1.];\n", " [1 1 0.]]\n", "Si = ElementPsp(:Si; psp=load_psp(\"hgh/lda/Si-q4\"))\n", "atoms = [Si, Si]\n", "positions = [ones(3)/8, -ones(3)/8]\n", "\n", "model = model_LDA(lattice, atoms, positions)\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 -7.847787614592 -0.70 4.8 91.5ms\n", " 2 -7.852409462876 -2.34 -1.53 1.0 18.9ms\n", " 3 -7.852610187022 -3.70 -2.54 1.0 18.8ms\n", " 4 -7.852645422740 -4.45 -2.78 2.2 25.3ms\n", " 5 -7.852646089369 -6.18 -2.87 1.2 20.3ms\n", " 6 -7.852646665562 -6.24 -3.81 1.0 19.3ms\n", " 7 -7.852646684815 -7.72 -4.66 1.5 21.7ms\n", " 8 -7.852646686683 -8.73 -5.02 2.0 24.4ms\n", " 9 -7.852646686723 -10.40 -5.56 1.0 19.4ms\n", " 10 -7.852646686728 -11.34 -5.58 1.8 23.5ms\n", " 11 -7.852646686730 -11.74 -6.20 1.0 19.6ms\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 -7.847756480124 -0.70 4.5 469ms\n", " 2 -7.852560493068 -2.32 -1.62 0.80 2.2 2.40s\n", " 3 -7.852641694852 -4.09 -2.69 0.80 1.0 234ms\n", " 4 -7.852646498055 -5.32 -3.42 0.80 1.8 21.8ms\n", " 5 -7.852646681896 -6.74 -4.42 0.80 2.0 22.6ms\n", " 6 -7.852646686624 -8.33 -4.89 0.80 2.2 24.6ms\n", " 7 -7.852646686725 -10.00 -5.50 0.80 1.2 19.3ms\n", " 8 -7.852646686730 -11.36 -6.50 0.80 1.5 20.6ms\n" ] } ], "cell_type": "code", "source": [ "scfres_scfv = DFTK.scf_potential_mixing(basis; tol);" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "## Direct minimization\n", "Note: Unlike the other algorithms, tolerance for this one is in the energy,\n", "thus we square the density tolerance value to be roughly equivalent." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n Energy log10(ΔE) log10(Δρ) Δtime\n", "--- --------------- --------- --------- ------\n", " 1 +1.124538990136 -1.01 4.52s\n", " 2 -1.345591143722 0.39 -0.69 144ms\n", " 3 -3.451642713493 0.32 -0.47 34.0ms\n", " 4 -4.562430912188 0.05 -0.58 34.1ms\n", " 5 -6.432544925586 0.27 -0.68 33.8ms\n", " 6 -7.223242472318 -0.10 -1.12 34.7ms\n", " 7 -7.532616065722 -0.51 -1.25 26.0ms\n", " 8 -7.680630806733 -0.83 -1.69 25.8ms\n", " 9 -7.750715360337 -1.15 -1.96 26.3ms\n", " 10 -7.790830884447 -1.40 -2.12 26.1ms\n", " 11 -7.819231479371 -1.55 -2.04 26.0ms\n", " 12 -7.837961942154 -1.73 -2.34 26.2ms\n", " 13 -7.848980106297 -1.96 -2.32 26.0ms\n", " 14 -7.851241415069 -2.65 -3.50 25.7ms\n", " 15 -7.852174911668 -3.03 -2.95 25.9ms\n", " 16 -7.852459472415 -3.55 -3.52 25.7ms\n", " 17 -7.852584481891 -3.90 -3.75 25.4ms\n", " 18 -7.852628127098 -4.36 -3.89 25.4ms\n", " 19 -7.852640254273 -4.92 -4.39 25.5ms\n", " 20 -7.852644677252 -5.35 -4.47 25.9ms\n", " 21 -7.852646109306 -5.84 -4.47 26.6ms\n", " 22 -7.852646478109 -6.43 -5.00 26.3ms\n", " 23 -7.852646618850 -6.85 -4.98 26.0ms\n", " 24 -7.852646660364 -7.38 -5.43 26.1ms\n", " 25 -7.852646675945 -7.81 -5.46 26.0ms\n", " 26 -7.852646682952 -8.15 -5.68 25.5ms\n", " 27 -7.852646685294 -8.63 -6.03 25.3ms\n", " 28 -7.852646686276 -9.01 -6.47 25.5ms\n", " 29 -7.852646686576 -9.52 -6.57 25.3ms\n", " 30 -7.852646686670 -10.03 -6.85 25.6ms\n", " 31 -7.852646686709 -10.40 -7.40 25.8ms\n", " 32 -7.852646686723 -10.86 -7.52 44.0ms\n", " 33 -7.852646686728 -11.33 -7.54 25.8ms\n", " 34 -7.852646686729 -11.85 -7.71 25.5ms\n", " 35 -7.852646686730 -12.35 -8.22 25.6ms\n", " 36 -7.852646686730 -12.78 -8.41 25.8ms\n", " 37 -7.852646686730 -13.29 -8.51 26.2ms\n", " 38 -7.852646686730 -13.75 -8.52 26.4ms\n", " 39 -7.852646686730 -14.10 -9.14 25.9ms\n", " 40 -7.852646686730 -14.10 -9.37 25.6ms\n", " 41 -7.852646686730 + -Inf -9.60 74.6ms\n", "┌ Warning: DM not converged.\n", "└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/scf_callbacks.jl:60\n" ] } ], "cell_type": "code", "source": [ "scfres_dm = direct_minimization(basis; tol=tol^2);" ], "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 -7.847715983527 -0.70 4.8 36.4ms\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 -7.852645825570 -1.63 14.2s\n", " 2 -7.852646686730 -6.06 -3.69 3.63s\n", " 3 -7.852646686730 -13.16 -7.18 136ms\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| = 6.427403060808416e-7\n", "|ρ_newton - ρ_scfv| = 7.50386640588271e-8\n", "|ρ_newton - ρ_dm| = 1.8391242184086972e-9\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.10.4" }, "kernelspec": { "name": "julia-1.10", "display_name": "Julia 1.10.4", "language": "julia" } }, "nbformat": 4 }