{ "cells": [ { "cell_type": "markdown", "source": [ "# Custom solvers\n", "In this example, we show how to define custom solvers. Our system\n", "will again be silicon, because we are not very imaginative" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using DFTK, LinearAlgebra\n", "\n", "a = 10.26\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", "# We take very (very) crude parameters\n", "model = model_LDA(lattice, atoms, positions)\n", "basis = PlaneWaveBasis(model; Ecut=5, kgrid=[1, 1, 1]);" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "We define our custom fix-point solver: simply a damped fixed-point" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function my_fp_solver(f, x0, max_iter; tol)\n", " mixing_factor = .7\n", " x = x0\n", " fx = f(x)\n", " for n = 1:max_iter\n", " inc = fx - x\n", " if norm(inc) < tol\n", " break\n", " end\n", " x = x + mixing_factor * inc\n", " fx = f(x)\n", " end\n", " (fixpoint=x, converged=norm(fx-x) < tol)\n", "end;" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Our eigenvalue solver just forms the dense matrix and diagonalizes\n", "it explicitly (this only works for very small systems)" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function my_eig_solver(A, X0; maxiter, tol, kwargs...)\n", " n = size(X0, 2)\n", " A = Array(A)\n", " E = eigen(A)\n", " λ = E.values[1:n]\n", " X = E.vectors[:, 1:n]\n", " (; λ, X, residual_norms=[], iterations=0, converged=true, n_matvec=0)\n", "end;" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "Finally we also define our custom mixing scheme. It will be a mixture\n", "of simple mixing (for the first 2 steps) and than default to Kerker mixing.\n", "In the mixing interface `δF` is $(ρ_\\text{out} - ρ_\\text{in})$, i.e.\n", "the difference in density between two subsequent SCF steps and the `mix`\n", "function returns $δρ$, which is added to $ρ_\\text{in}$ to yield $ρ_\\text{next}$,\n", "the density for the next SCF step." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "struct MyMixing\n", " n_simple # Number of iterations for simple mixing\n", "end\n", "MyMixing() = MyMixing(2)\n", "\n", "function DFTK.mix_density(mixing::MyMixing, basis, δF; n_iter, kwargs...)\n", " if n_iter <= mixing.n_simple\n", " return δF # Simple mixing -> Do not modify update at all\n", " else\n", " # Use the default KerkerMixing from DFTK\n", " DFTK.mix_density(KerkerMixing(), basis, δF; kwargs...)\n", " end\n", "end" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "That's it! Now we just run the SCF with these solvers" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n Energy log10(ΔE) log10(Δρ) Diag Δtime\n", "--- --------------- --------- --------- ---- ------\n", " 1 -7.235992833898 -0.50 0.0 \n", " 2 -7.249608665538 -1.87 -0.92 0.0 864ms\n", " 3 -7.251174940533 -2.81 -1.34 0.0 103ms\n", " 4 -7.251295969966 -3.92 -1.65 0.0 93.1ms\n", " 5 -7.251327319418 -4.50 -1.95 0.0 98.9ms\n", " 6 -7.251335595391 -5.08 -2.25 0.0 245ms\n", " 7 -7.251337861796 -5.64 -2.54 0.0 99.3ms\n", " 8 -7.251338511556 -6.19 -2.83 0.0 98.5ms\n", " 9 -7.251338706958 -6.71 -3.10 0.0 94.7ms\n", " 10 -7.251338768376 -7.21 -3.36 0.0 93.5ms\n", " 11 -7.251338788413 -7.70 -3.62 0.0 94.7ms\n", " 12 -7.251338795145 -8.17 -3.87 0.0 98.0ms\n", " 13 -7.251338797456 -8.64 -4.11 0.0 97.7ms\n", " 14 -7.251338798263 -9.09 -4.35 0.0 97.7ms\n", " 15 -7.251338798547 -9.55 -4.58 0.0 245ms\n", " 16 -7.251338798648 -10.00 -4.81 0.0 97.3ms\n", " 17 -7.251338798684 -10.44 -5.04 0.0 96.7ms\n", " 18 -7.251338798697 -10.89 -5.27 0.0 100ms\n", " 19 -7.251338798702 -11.33 -5.50 0.0 94.7ms\n", " 20 -7.251338798704 -11.78 -5.72 0.0 98.4ms\n", " 21 -7.251338798704 -12.22 -5.95 0.0 98.5ms\n", " 22 -7.251338798704 -12.66 -6.17 0.0 97.1ms\n", " 23 -7.251338798704 -13.11 -6.39 0.0 234ms\n", " 24 -7.251338798705 -13.62 -6.62 0.0 94.0ms\n", " 25 -7.251338798705 -13.97 -6.83 0.0 97.2ms\n", " 26 -7.251338798705 -14.21 -7.05 0.0 99.3ms\n", " 27 -7.251338798705 -15.05 -7.28 0.0 97.7ms\n", " 28 -7.251338798705 -14.75 -7.42 0.0 95.5ms\n", " 29 -7.251338798705 + -15.05 -7.66 0.0 103ms\n", " 30 -7.251338798705 + -13.65 -7.01 0.0 99.5ms\n", " 31 -7.251338798705 -13.75 -7.26 0.0 222ms\n", " 32 -7.251338798705 -14.75 -7.60 0.0 93.5ms\n", " 33 -7.251338798705 + -Inf -7.88 0.0 96.0ms\n", " 34 -7.251338798705 + -15.05 -7.42 0.0 96.5ms\n", " 35 -7.251338798705 -15.05 -7.67 0.0 104ms\n", " 36 -7.251338798705 -14.75 -7.96 0.0 95.5ms\n", " 37 -7.251338798705 + -14.75 -7.48 0.0 101ms\n", " 38 -7.251338798705 -14.75 -7.73 0.0 104ms\n", " 39 -7.251338798705 + -14.27 -7.32 0.0 105ms\n", " 40 -7.251338798705 -14.15 -7.53 0.0 216ms\n", " 41 -7.251338798705 + -14.27 -7.50 0.0 97.8ms\n", " 42 -7.251338798705 -14.75 -7.74 0.0 105ms\n", " 43 -7.251338798705 -14.45 -7.83 0.0 102ms\n", " 44 -7.251338798705 + -13.85 -7.15 0.0 98.8ms\n", " 45 -7.251338798705 -13.97 -7.39 0.0 100ms\n", " 46 -7.251338798705 + -Inf -7.71 0.0 101ms\n", " 47 -7.251338798705 -14.45 -7.82 0.0 111ms\n", " 48 -7.251338798704 + -13.36 -6.81 0.0 233ms\n", " 49 -7.251338798705 -13.44 -7.08 0.0 96.4ms\n", " 50 -7.251338798705 -14.45 -7.43 0.0 99.4ms\n", " 51 -7.251338798705 -15.05 -7.61 0.0 102ms\n", " 52 -7.251338798705 + -13.82 -7.08 0.0 94.7ms\n", " 53 -7.251338798705 -13.75 -7.32 0.0 94.3ms\n", " 54 -7.251338798705 + -14.75 -7.62 0.0 99.1ms\n", " 55 -7.251338798705 + -14.75 -7.46 0.0 100ms\n", " 56 -7.251338798705 + -14.75 -7.24 0.0 224ms\n", " 57 -7.251338798705 -14.35 -7.49 0.0 93.2ms\n", " 58 -7.251338798705 + -14.57 -7.51 0.0 98.1ms\n", " 59 -7.251338798705 -14.57 -7.73 0.0 97.5ms\n", " 60 -7.251338798705 + -15.05 -7.84 0.0 118ms\n", " 61 -7.251338798705 + -Inf -7.81 0.0 104ms\n", " 62 -7.251338798705 + -Inf -7.48 0.0 98.9ms\n", " 63 -7.251338798705 + -14.75 -7.57 0.0 101ms\n", " 64 -7.251338798705 -15.05 -7.47 0.0 98.9ms\n", " 65 -7.251338798705 + -15.05 -7.70 0.0 225ms\n", " 66 -7.251338798705 -14.75 -7.60 0.0 95.6ms\n", " 67 -7.251338798705 + -15.05 -7.73 0.0 93.9ms\n", " 68 -7.251338798705 -14.75 -7.84 0.0 98.1ms\n", " 69 -7.251338798705 + -14.35 -7.50 0.0 95.3ms\n", " 70 -7.251338798705 -14.45 -7.42 0.0 99.3ms\n", " 71 -7.251338798705 + -14.75 -7.65 0.0 103ms\n", " 72 -7.251338798705 -14.75 -7.86 0.0 101ms\n", " 73 -7.251338798705 + -14.10 -7.21 0.0 219ms\n", " 74 -7.251338798705 -14.10 -7.44 0.0 96.6ms\n", " 75 -7.251338798705 + -Inf -7.67 0.0 111ms\n", " 76 -7.251338798705 + -Inf -7.71 0.0 96.0ms\n", " 77 -7.251338798705 + -14.75 -7.46 0.0 99.3ms\n", " 78 -7.251338798705 -14.75 -7.74 0.0 96.6ms\n", " 79 -7.251338798705 + -14.45 -7.31 0.0 103ms\n", " 80 -7.251338798705 -14.57 -7.52 0.0 100ms\n", " 81 -7.251338798705 -15.05 -7.81 0.0 225ms\n", " 82 -7.251338798705 + -Inf -7.90 0.0 95.1ms\n", " 83 -7.251338798705 + -15.05 -7.31 0.0 97.9ms\n", " 84 -7.251338798705 + -14.57 -7.48 0.0 97.9ms\n", " 85 -7.251338798705 + -Inf -7.40 0.0 98.8ms\n", " 86 -7.251338798705 -14.45 -7.60 0.0 94.1ms\n", " 87 -7.251338798705 -14.57 -7.83 0.0 99.4ms\n", " 88 -7.251338798705 + -14.21 -7.33 0.0 101ms\n", " 89 -7.251338798705 -14.75 -7.51 0.0 99.3ms\n", " 90 -7.251338798705 -15.05 -7.65 0.0 225ms\n", " 91 -7.251338798705 + -14.75 -7.33 0.0 98.5ms\n", " 92 -7.251338798705 + -Inf -7.45 0.0 102ms\n", " 93 -7.251338798705 -14.35 -7.68 0.0 105ms\n", " 94 -7.251338798705 + -14.75 -7.87 0.0 103ms\n", " 95 -7.251338798705 + -13.80 -7.08 0.0 106ms\n", " 96 -7.251338798705 -13.75 -7.32 0.0 105ms\n", " 97 -7.251338798705 + -Inf -7.62 0.0 98.8ms\n", " 98 -7.251338798705 + -14.45 -7.45 0.0 231ms\n", " 99 -7.251338798705 -14.57 -7.68 0.0 93.8ms\n", " 100 -7.251338798705 + -13.97 -7.16 0.0 101ms\n", " 101 -7.251338798705 -14.01 -7.41 0.0 102ms\n", "┌ Warning: SCF not converged.\n", "└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/scf_callbacks.jl:38\n" ] } ], "cell_type": "code", "source": [ "scfres = self_consistent_field(basis;\n", " tol=1e-8,\n", " solver=my_fp_solver,\n", " eigensolver=my_eig_solver,\n", " mixing=MyMixing());" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "Note that the default convergence criterion is the difference in\n", "density. When this gets below `tol`, the\n", "\"driver\" `self_consistent_field` artificially makes the fixed-point\n", "solver think it's converged by forcing `f(x) = x`. You can customize\n", "this with the `is_converged` keyword argument to\n", "`self_consistent_field`." ], "metadata": {} } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.5" }, "kernelspec": { "name": "julia-1.8", "display_name": "Julia 1.8.5", "language": "julia" } }, "nbformat": 4 }