{ "cells": [ { "cell_type": "markdown", "source": [ "# Analysing SCF convergence\n", "\n", "The goal of this example is to explain the differing convergence behaviour\n", "of SCF algorithms depending on the choice of the mixing.\n", "For this we look at the eigenpairs of the Jacobian governing the SCF convergence,\n", "that is\n", "$$\n", "1 - α P^{-1} \\varepsilon^\\dagger \\qquad \\text{with} \\qquad \\varepsilon^\\dagger = (1-\\chi_0 K).\n", "$$\n", "where $α$ is the damping $P^{-1}$ is the mixing preconditioner\n", "(e.g. `KerkerMixing`, `LdosMixing`)\n", "and $\\varepsilon^\\dagger$ is the dielectric operator.\n", "\n", "We thus investigate the largest and smallest eigenvalues of\n", "$(P^{-1} \\varepsilon^\\dagger)$ and $\\varepsilon^\\dagger$.\n", "The ratio of largest to smallest eigenvalue of this operator is the condition number\n", "$$\n", "\\kappa = \\frac{\\lambda_\\text{max}}{\\lambda_\\text{min}},\n", "$$\n", "which can be related to the rate of convergence of the SCF.\n", "The smaller the condition number, the faster the convergence.\n", "For more details on SCF methods, see Self-consistent field methods.\n", "\n", "For our investigation we consider a crude aluminium setup:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "FlexibleSystem(Al₁₆, periodicity = TTT):\n cell_vectors : [ 16.2 0 0;\n 0 4.05 0;\n 0 0 4.05]u\"Å\"\n" }, "metadata": {}, "execution_count": 1 } ], "cell_type": "code", "source": [ "using AtomsBuilder\n", "using DFTK\n", "\n", "system_Al = bulk(:Al; cubic=true) * (4, 1, 1)" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "and we discretise:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using PseudoPotentialData\n", "\n", "pseudopotentials = PseudoFamily(\"dojo.nc.sr.lda.v0_4_1.standard.upf\")\n", "model_Al = model_DFT(system_Al; functionals=LDA(), temperature=1e-3,\n", " symmetries=false, pseudopotentials)\n", "basis_Al = PlaneWaveBasis(model_Al; Ecut=7, kgrid=[1, 1, 1]);" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "On aluminium (a metal) already for moderate system sizes (like the 8 layers\n", "we consider here) the convergence without mixing / preconditioner is slow:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n Energy log10(ΔE) log10(Δρ) Diag Δtime\n", "--- --------------- --------- --------- ---- ------\n", " 1 -36.73372409678 -0.88 11.0 480ms\n", " 2 -36.54927769127 + -0.73 -1.33 1.0 157ms\n", " 3 +61.16604959193 + 1.99 -0.06 22.0 379ms\n", " 4 -35.82323484490 1.99 -0.96 10.0 333ms\n", " 5 -25.75016816170 + 1.00 -0.54 5.0 185ms\n", " 6 -36.36466582549 1.03 -1.19 4.0 171ms\n", " 7 -36.69438399901 -0.48 -1.57 3.0 170ms\n", " 8 -36.71542333898 -1.68 -1.73 2.0 113ms\n", " 9 -36.73973626349 -1.61 -2.06 2.0 108ms\n", " 10 -36.73920968223 + -3.28 -2.06 3.0 131ms\n", " 11 -36.74133415779 -2.67 -2.28 2.0 110ms\n", " 12 -36.74137938893 -4.34 -2.30 2.0 118ms\n", " 13 -36.74152077867 -3.85 -2.42 1.0 99.0ms\n", " 14 -36.74221927551 -3.16 -2.62 1.0 98.8ms\n", " 15 -36.74132123693 + -3.05 -2.47 2.0 124ms\n", " 16 -36.73596244634 + -2.27 -2.11 3.0 174ms\n", " 17 -36.73736756507 -2.85 -2.21 3.0 157ms\n", " 18 -36.74144784101 -2.39 -2.53 3.0 146ms\n", " 19 -36.74150148859 -4.27 -2.50 3.0 147ms\n", " 20 -36.74244708404 -3.02 -2.94 2.0 118ms\n", " 21 -36.74247259177 -4.59 -3.09 2.0 136ms\n", " 22 -36.74251298891 -4.39 -3.68 2.0 117ms\n", " 23 -36.74251206447 + -6.03 -3.76 2.0 159ms\n", " 24 -36.74251434697 -5.64 -4.10 1.0 99.0ms\n", " 25 -36.74251450683 -6.80 -4.01 2.0 127ms\n", " 26 -36.74251455560 -7.31 -4.33 2.0 115ms\n", " 27 -36.74251474605 -6.72 -4.43 2.0 124ms\n", " 28 -36.74251453687 + -6.68 -4.36 2.0 127ms\n", " 29 -36.74251476405 -6.64 -4.97 2.0 117ms\n", " 30 -36.74251476752 -8.46 -5.13 3.0 147ms\n", " 31 -36.74251477028 -8.56 -5.20 3.0 178ms\n", " 32 -36.74251477111 -9.08 -5.33 1.0 95.7ms\n", " 33 -36.74251477251 -8.85 -5.66 2.0 117ms\n", " 34 -36.74251477235 + -9.79 -5.55 3.0 148ms\n", " 35 -36.74251477230 + -10.25 -5.56 3.0 139ms\n", " 36 -36.74251477302 -9.14 -6.29 2.0 117ms\n", " 37 -36.74251477290 + -9.91 -5.98 3.0 148ms\n", " 38 -36.74251477303 -9.89 -6.59 3.0 146ms\n", " 39 -36.74251477300 + -10.46 -6.15 3.0 176ms\n", " 40 -36.74251477303 -10.45 -6.60 3.0 135ms\n", " 41 -36.74251477304 -11.42 -6.98 2.0 117ms\n", " 42 -36.74251477304 -12.58 -7.06 3.0 147ms\n", " 43 -36.74251477304 -12.56 -7.12 1.0 98.3ms\n", " 44 -36.74251477304 -12.89 -7.36 2.0 117ms\n", " 45 -36.74251477304 -12.77 -7.77 1.0 98.4ms\n", " 46 -36.74251477304 + -12.87 -7.39 3.0 155ms\n", " 47 -36.74251477304 -13.30 -7.40 3.0 173ms\n", " 48 -36.74251477304 + -Inf -7.55 3.0 177ms\n", " 49 -36.74251477304 -13.15 -7.81 3.0 143ms\n", " 50 -36.74251477304 -13.67 -8.21 2.0 114ms\n", " 51 -36.74251477304 -14.15 -8.18 3.0 140ms\n", " 52 -36.74251477304 + -14.15 -8.16 2.0 115ms\n", " 53 -36.74251477304 + -Inf -8.22 3.0 123ms\n", " 54 -36.74251477304 -14.15 -8.28 2.0 149ms\n", " 55 -36.74251477304 + -Inf -8.90 1.0 99.3ms\n", " 56 -36.74251477304 + -Inf -8.68 3.0 156ms\n", " 57 -36.74251477304 + -Inf -8.67 2.0 119ms\n", " 58 -36.74251477304 + -Inf -9.29 2.0 117ms\n", " 59 -36.74251477304 + -13.85 -8.60 3.0 166ms\n", " 60 -36.74251477304 -14.15 -9.27 4.0 167ms\n", " 61 -36.74251477304 + -Inf -8.59 4.0 169ms\n", " 62 -36.74251477304 -14.15 -9.10 3.0 180ms\n", " 63 -36.74251477304 + -13.85 -9.24 2.0 129ms\n", " 64 -36.74251477304 -14.15 -10.09 2.0 114ms\n", " 65 -36.74251477304 -14.15 -9.70 4.0 176ms\n", " 66 -36.74251477304 + -13.85 -10.20 3.0 146ms\n", " 67 -36.74251477304 + -Inf -10.24 1.0 101ms\n", " 68 -36.74251477304 + -Inf -10.57 1.0 163ms\n", " 69 -36.74251477304 -14.15 -10.49 3.0 692ms\n", " 70 -36.74251477304 + -14.15 -10.76 2.0 115ms\n", " 71 -36.74251477304 + -Inf -10.80 2.0 115ms\n", " 72 -36.74251477304 -13.85 -10.98 2.0 108ms\n", " 73 -36.74251477304 + -13.85 -11.19 2.0 138ms\n", " 74 -36.74251477304 -13.85 -10.63 3.0 168ms\n", " 75 -36.74251477304 + -Inf -11.35 3.0 165ms\n", " 76 -36.74251477304 + -13.85 -11.02 3.0 161ms\n", " 77 -36.74251477304 -13.85 -11.74 3.0 159ms\n", " 78 -36.74251477304 + -14.15 -11.93 2.0 138ms\n", " 79 -36.74251477304 + -Inf -12.23 2.0 121ms\n" ] } ], "cell_type": "code", "source": [ "# Note: DFTK uses the self-adapting LdosMixing() by default, so to truly disable\n", "# any preconditioning, we need to supply `mixing=SimpleMixing()` explicitly.\n", "scfres_Al = self_consistent_field(basis_Al; tol=1e-12, mixing=SimpleMixing());" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "while when using the Kerker preconditioner it is much faster:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "n Energy log10(ΔE) log10(Δρ) Diag Δtime\n", "--- --------------- --------- --------- ---- ------\n", " 1 -36.73113129843 -0.88 11.0 412ms\n", " 2 -36.73922462098 -2.09 -1.36 1.0 657ms\n", " 3 -36.73897761884 + -3.61 -1.61 3.0 128ms\n", " 4 -36.74219400837 -2.49 -2.16 2.0 112ms\n", " 5 -36.74231044516 -3.93 -2.40 5.0 114ms\n", " 6 -36.74244084393 -3.88 -2.45 2.0 103ms\n", " 7 -36.74248499908 -4.36 -3.12 1.0 90.2ms\n", " 8 -36.74251177551 -4.57 -3.22 3.0 151ms\n", " 9 -36.74251295375 -5.93 -3.45 1.0 92.0ms\n", " 10 -36.74251418318 -5.91 -3.76 2.0 101ms\n", " 11 -36.74251471578 -6.27 -4.21 6.0 124ms\n", " 12 -36.74251476845 -7.28 -4.66 2.0 129ms\n", " 13 -36.74251476124 + -8.14 -4.66 2.0 107ms\n", " 14 -36.74251477267 -7.94 -5.35 1.0 96.0ms\n", " 15 -36.74251477297 -9.52 -5.53 3.0 140ms\n", " 16 -36.74251477301 -10.40 -5.71 2.0 98.3ms\n", " 17 -36.74251477292 + -10.04 -5.84 2.0 104ms\n", " 18 -36.74251477304 -9.95 -6.58 1.0 114ms\n", " 19 -36.74251477304 -12.10 -6.66 8.0 172ms\n", " 20 -36.74251477304 + -11.98 -6.82 1.0 95.9ms\n", " 21 -36.74251477304 -12.13 -7.08 2.0 131ms\n", " 22 -36.74251477304 -12.28 -7.57 2.0 106ms\n", " 23 -36.74251477304 + -14.15 -7.74 3.0 139ms\n", " 24 -36.74251477304 -14.15 -7.94 1.0 96.1ms\n", " 25 -36.74251477304 + -Inf -8.52 2.0 103ms\n", " 26 -36.74251477304 + -13.85 -8.41 3.0 143ms\n", " 27 -36.74251477304 -13.85 -8.65 1.0 112ms\n", " 28 -36.74251477304 + -13.85 -9.20 2.0 105ms\n", " 29 -36.74251477304 + -Inf -9.49 3.0 144ms\n", " 30 -36.74251477304 -13.85 -9.91 1.0 96.2ms\n", " 31 -36.74251477304 + -Inf -10.16 3.0 142ms\n", " 32 -36.74251477304 + -Inf -10.56 1.0 96.0ms\n", " 33 -36.74251477304 + -14.15 -10.86 2.0 111ms\n", " 34 -36.74251477304 + -Inf -11.31 3.0 128ms\n", " 35 -36.74251477304 + -Inf -11.36 2.0 131ms\n", " 36 -36.74251477304 -14.15 -11.80 2.0 104ms\n", " 37 -36.74251477304 + -14.15 -12.08 4.0 132ms\n" ] } ], "cell_type": "code", "source": [ "scfres_Al = self_consistent_field(basis_Al; tol=1e-12, mixing=KerkerMixing());" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "Given this `scfres_Al` we construct functions representing\n", "$\\varepsilon^\\dagger$ and $P^{-1}$:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "epsilon (generic function with 1 method)" }, "metadata": {}, "execution_count": 5 } ], "cell_type": "code", "source": [ "# Function, which applies P^{-1} for the case of KerkerMixing\n", "Pinv_Kerker(δρ) = DFTK.mix_density(KerkerMixing(), basis_Al, δρ)\n", "\n", "# Function which applies ε† = 1 - χ0 K\n", "function epsilon(δρ)\n", " δV = apply_kernel(basis_Al, δρ; ρ=scfres_Al.ρ)\n", " χ0δV = apply_χ0(scfres_Al, δV)\n", " δρ - χ0δV\n", "end" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "With these functions available we can now compute the desired eigenvalues.\n", "For simplicity we only consider the first few largest ones." ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "┌ Info: Arnoldi eigsolve finished after 1 iterations:\n", "│ * 7 eigenvalues converged\n", "│ * norm of residuals = (8.18e-06, 2.11e-08, 6.48e-09, 1.67e-05, 9.36e-06, 2.68e-05, 2.73e-05)\n", "└ * number of operations = 15\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "44.02448779307413" }, "metadata": {}, "execution_count": 6 } ], "cell_type": "code", "source": [ "using KrylovKit\n", "λ_Simple, X_Simple = eigsolve(epsilon, randn(size(scfres_Al.ρ)), 3, :LM;\n", " tol=1e-3, eager=true, verbosity=2)\n", "λ_Simple_max = maximum(real.(λ_Simple))" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "The smallest eigenvalue is a bit more tricky to obtain, so we will just assume" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "0.952" }, "metadata": {}, "execution_count": 7 } ], "cell_type": "code", "source": [ "λ_Simple_min = 0.952" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "This makes the condition number around 30:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "46.24420986667451" }, "metadata": {}, "execution_count": 8 } ], "cell_type": "code", "source": [ "cond_Simple = λ_Simple_max / λ_Simple_min" ], "metadata": {}, "execution_count": 8 }, { "cell_type": "markdown", "source": [ "This does not sound large compared to the condition numbers you might know\n", "from linear systems.\n", "\n", "However, this is sufficient to cause a notable slowdown, which would be even more\n", "pronounced if we did not use Anderson, since we also would need to drastically\n", "reduce the damping (try it!)." ], "metadata": {} }, { "cell_type": "markdown", "source": [ " Having computed the eigenvalues of the dielectric matrix\n", "we can now also look at the eigenmodes, which are responsible for\n", "the bad convergence behaviour. The largest eigenmode for example:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=1}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de5xU9X3/8e9cd/bCLuyyF3aRyGUjFfyJVVhUHisGiBcCUazEREyLMbV5WNQW46N9pE1JrZaHNmoSNU1SiTHapKH6QJMHRDRQUCs3EY3EW1EIt+W+7G3uc35/TLIlzvvYmd0ZZtbv6/nXznfPnvnOnDPz2dn9vs/H4ziOAQDAVt5iTwAAgGKiEAIArEYhBABYjUIIALAahRAAYDUKIQDAahRCAIDVKIQAAKtRCAEAVqMQAgCsVhKF8N577+3o6Cj2LAAANiqJQvj000/v2bOn2LMAANioJAohAADFQiEEAFiNQggAsBqFEABgNQohAMBqFEIAgNUohAAAq1EIAQBWoxACAKxGIQQAWI1CCACwGoUQAGA1CiEAwGoUQgCA1SiEAACrUQgBAFajEAIArEYhBABYjUIIALAahRAAYDUKIQDAahRCAIDVKIQAAKtRCAEAVqMQAgCsRiEEAFiNQggAsBqFEABgNQohAMBqFEIAgNUohAAAq1EIAQBWoxACAKxGIQQAWI1CCACwGoUQAGA1CiEAwGoUQgCA1SiEAACrUQgBAFajEAIArEYhBABYjUIIALAahRAAYDUKIQDAahRCAIDVKIQAAKv5B7+LAwcO/PznP9+1a1dFRcWcOXMuvvji/m+tXbv2hRdeaGhouOmmm4YPHz74+wIAIL/y8Ilw06ZNr732WnNzczKZnD9//ve///30+KOPPrp48eKWlpatW7e2t7cnEonB3xcAAPnlcRwnj7v79re//dRTT23YsCGVSn3yk5+8//7758+fn0qlJk+efNddd11zzTXyp6ZPn/6tb32rra0tjzMBACAb+fwfYTwe37Zt2+TJk40x+/bte//99+fMmWOM8Xq9s2fP3rBhQx7vCwCAvMjD/wiNMW+//fbcuXMPHTo0ZcqUtWvXGmM6OjqqqqrKy8vTGzQ0NOzYscPtxw8fPrx8+fKGhob0zVGjRi1btiwvEwMA4KPl5xNha2vr1q1bX3zxxfLy8iVLlhhj/H5/Mpns3yCRSASDQbcfD4VCEydOPP/3zj333LzMCgCA/1N+PhH6fL7a2tra2tp77rnnkksuefTRR5ubm/v6+jo7O9OLRffv39/c3Oz249XV1VdddRX/IwQAnH55+EQYDof7v96+ffsZZ5xhjGlqarrgggtWrlxpjOnt7V2zZs28efMGf18AAORXHj4R/tmf/dmRI0c+8YlP7N+//9VXX/33f//39Pg999zz+c9//pVXXtmxY8cFF1zQ3t4++PsCACC/8hCfCIfDmzZt2r9/f11d3UUXXVRTU9P/rd/+9rcbN25sbm6eOXOm1+v66ZP4BACgWPLwibC8vPzSSy+V3xozZsyiRYsGfxcAABQI1xoFAFiNQggAsBqFEABgNQohAMBqFEIAgNUohAAAq1EIAQBWoxACAKxGIQQAWI1CCACwGoUQAGA1CiEAwGoUQgCA1SiEAACrUQgBAFajEAIArEYhBABYjUIIALAahRAAYDUKIQDAahRCAIDVKIQAAKtRCAEAVqMQAgCsRiEEAFiNQggAsBqFEABgNQohAMBqFEIAgNUohAAAq1EIAQBWoxACAKxGIQQAWM1f7AkYY8xZ514QrZ+w8f1jHxr/n+O9mRtvev945uCv3/vwz6Yd6+jOHAx3dmYOJmPhzEGvPyB36yurEKOOI+7rxKHMwc49b8rdFp0vWJ45WDZshNxYPmPR7hN5npMxxpjyEU2Zg41nT8scbJlQK/cweczwzMEzR1ZmDtZXBTMH68rFoDGmrkKcISNC4mVVXeYTW5bp30R9XR1isFsMxvftyhwMv7NT7rZj61ticIc4RU/sEy+cznhS7lad+KbC58kcDHjFYG8yJXf7Pz3xzMHD0YTcuOjq1fH1e8TjPRYTT2MspZ7EXIh7Mub84SG58dn/ryFzcNT5Z2QOVo9rzhwsbxgpd+sdUZ856BuuBtWWqQqXt5oqsXGvRzy0Y2Hx3LY2DJO7PRWfCAEAVqMQAgCsRiEEAFiNQggAsBqFEABgNQohAMBqFEIAgNVKIkdojJNSWaSkytaoJJIJBEWIxxgTUjGvVKI6czAWFntwki7BqaQIM6USscxBGbYrWXK2iUiZ3DiVFDGvwfMFxN0FKsQh86hTIeWSx4olRFgtoTaWqTiVBzPGGJ/6hlcPyh24pcfUuHxoKXGKOvIxGOMNiNd72TARkaysVrnJrqjcbTipnka1pUzLyR9327hk9SXEbAPqg0aBHladCjKGfPqjjjxDknH1tibDo46Ofupz3Kvm4FXv2HLQ6JeffBYH/NTyiRAAYDUKIQDAahRCAIDVKIQAAKtRCAEAVqMQAgCsVhLxCcfRC9nlsmq/WoxbVqEfSMUwvfQ/k1yLH+vrkxsnoqI/VLxPdK6J953McgIlK9YrulYZYzwecSA8MjngF2vxZSLCGBOsEv2SQjWqk4tq9JNSMQljTJ/qfRNVS8Pl6na3+IRfnTYyKSEzFR6XZegyw+GkxOp2J5FDWyJ/uXg5lNeJ3lvyafS4rMX3dUYyB3vUHuRgpzo0xpg+l/ZMpUk2k/KoRyZPJZnDMS7drGoCImZQpyJk1TLAYXTOIBlRAbCoGHQ762R6waNCEY4c9LjFJ9TG6uxwCw79n/hECACwGoUQAGA1CiEAwGoUQgCA1SiEAACrUQgBAFYrifhEyrUPgBgsUwuCR1Soi+Ubk1QJDLnEX68nTuiF3bGeE2JQxQxkpuLjwVHrl/1lFZmDZTUjMwdlIsIYE6wU8YlghditX60Xd+s+EY6JBd8RHZ9QPU9cLqzvU1EJPSh3oHpHGGM8jhqXa9ZVIxR9hhsTqAhlDlbUVcqNs9+tzFr0ZJ2p6IzrmMTQ6j4hZf8Aal1yDvVl4i06+6REqEbnx2RaLBkV51IiLFqOpFSfCmOMTj/41LBMSrh1n1DjKfXePuBThk+EAACrUQgBAFajEAIArEYhBABYjUIIALBa3grhiRMnYjFxedYjR45Eo2LdEQAApSAP8YmHH374n/7pn+LxeDgcvuKKKx577LGqqipjzO7du6+++urDhw/39vbeddddS5Yscd+HI9esS2V+UbyrywNyY53KUINxtbze16uX86bUqvdkNCwG44P9JUC2bvCoFhzJmFiwfpr5gmKBftmw2szB0PAGuYeyctWqQq0X96tr8LuJJdRBVyuw5ZkoG0q4jWfffUJfQt8YkxLjjsxaqD145IJ1l+4TwWoRTUmpXgpJlzYRkS5xkntVfCKqXnqD7zIhn/CgyzGTK+xPZ1RDngnDXOMT4lAOrxRvd8FK9cIp1+/wXvUuKg+6TErILY1LYw35ZmXUoGxJYYxxVKOblFGzLWL3iVGjRm3cuPHo0aMHDhzYt2/f3XffnR5funRpe3v7/v37N2/e/LWvfW3Xrl2Dvy8AAPIrD4VwwYIFra2txpiamprLLrvsvffeM8acPHny2Wefvf32240xZ5111mWXXfbkk08O/r4AAMivfC6WicViq1atmjVrljFm7969Xq937Nix6W+dddZZH3zwgdsPJpP6Ty4AABRa3gqh4zhLliypq6v78z//c2NMV1dXxSmXxaqsrOzs1I3OjTFHjx7N1zQAAMhJ3q41unTp0h07djz//PM+n88YU19f393dnUqlvF6vMaazs7OhQS+OMMY0NjbmaxoAAOQkP58Iv/a1r61fv37NmjXV1dXpkTFjxlRVVW3fvj19c+vWreecc05e7gsAgDzKwyfC5cuXP/jgg9/85jd/9atfGWPq6uo+9alPlZWV3XjjjXfeeecjjzzy0ksv7dix46mnnhr8fQEAkF95KISO48ydO3fdunXpm5/85Cc/9alPGWPuueeev//7v//c5z43atSoNWvWjBgxwn0XOrwlUymyx03Apz/aBlVcxqc2loNufWdkexUd88qF1ycOh7+sPHPQ4xVbus02oQKOBeJRSSCvX8TX/AF97sl0oE8dR69buE+RLb1k6sjJoXmOS3Aq+593jT3JMyzribk8Mx6/OjrqQPiCatAluClDaTo2qR7BYFOExoRUKC2kG1/pHGFcvXgLFC2U0wq4vHjL1KH0h9QbRUgd3KB+Y5RtmFze1uTTkMtzox6ak8urxOVkynIsK3kohH/7t38rx0Oh0H333XffffcN/i4AACgQrjUKALAahRAAYDUKIQDAahRCAIDVKIQAAKvl7coyg+LJfn2sScomSi5tQWIJ1U1GtphRg3LNvTF6BbRMDuQklRQdT2T4oWTbMMkMSSoh2vQkVHsX45KUkKeH1yUzI8lsiWyI48kp/pD1oOaWz9FnWNYTc2kq5CTU0VEHIqm6kiXjOiCUUq8yHYiSx1HuNBcR1bUq5XIc5XNz+powudxX3OWtRjauKo+IoyNfDm6RKqc827yaDlrkFH5QD82T01Oedboul2n9AT4RAgCsRiEEAFiNQggAsBqFEABgNQohAMBqFEIAgNVKIz5hPHIhuxRVa7W7wnG5cWdvLHMw0is2joXVenG13NwY41VJCZ9qE+GLiMYLybiIE7hJJcRDKFkywhHtPp456JY2ScWHi91WVKg95LBYOuiXTUtUfELt1SWPoMddBuUScJffRFU8Rj9jag9OUp+3ibA48WJdfZmDkU4R2ol26vM2rl5QMs8keylUuMRgolm3c5FPeCR5OjMROZBnQndcB8COeFXiRT2y6qjYMlQj3n+MMYHKQOagDGDI5iRuySWdJlLhFo8e1Efc44iN5QyyryPZ7A0AAFtQCAEAVqMQAgCsRiEEAFiNQggAsBqFEABgtZKIT3iN8at11fLS6VG1zvhEn84Y9HWJBd993WIwrIIWyZheL+4NBDMHg5Vi3X9KJSWSJ3OIT5Qsj1q4n4yJZfeRE4fElqqrhjEmXtUlNq5pyBwMBOszB70umYryoDjVQwERSJArsFMuF8uXvVD0oNyBT2dIHI8Y9/jVq9UnBt26psT7RLil71hv5mDvITHYd1QfsmhPtgmfKtVaJO6STelVDWFibkGWkpR944YTLvGJaEpEU7pVhKwuKM6ZFpeJVdaLPJKvTJxL/nIRwJCZClcyzOOoQbfAjBr3elSuY6DtJ/hECACwGoUQAGA1CiEAwGoUQgCA1SiEAACrUQgBAFajEAIArFYSOUKPR+cIVZMck1D9O6J9oomScYkMynBhLCyCU27tbHyBUOagp1Lm6kRyK3LyqNxtaZL5SGNMKikSTomIeBpl26nkySNyt/E+kSOU4avKupGZg16VVDPGVKiUVZnOEYofdwnmmYSKtWXfhslxa8Oksower3i16nChC9mGKXxMpANljrBXvXCMMWGVkZSBrmHq6LgFA2W0LpZ1b6bTrFJ1Jgqok6lHpQATLmdYd0L1bFJ7iKg3RvmEG2Pq1OHxhUQ82lcmBt3OOnnQHZkCVIMeGS40OnToUTlCGT3PBp8IAQBWoxACAKxGIQQAWI1CCACwGoUQAGA1CiEAwGolEZ8wxiN73/jU4mO5Nj0e0+tuI31iiX+0RzX6Uf2DvP6A3K03oJqYeMVafF+wXO6hNMnZ+kPiwRqXZ0zGJ3IisxYyU+GoU8GtDVNQrSN3af4lftwtPpFU35BJCZf2QW6rvdW4fGjqrHNbRJ6Ki5RRtFs0UertEoOdLq2C5JNTocJPMk5QLmNSxgQH3FOnGCr8YrZ+dSD65HvVoLtLHYuK/UbK9SGTZ4hPNVfyqoiRcYn96Aehch2645JbNkadYfIcH/AZwydCAIDVKIQAAKtRCAEAVqMQAgCsRiEEAFiNQggAsFpJxCc8LstePXJYXtrfZXF6Sl2mPZUQC/RlLwWPz+Ui63p5uuwYMJR+1fDItfhqMP2Nws7mFKmEODqOXlStV1DLfE72clrcnsPGuU1KLhjPetAlcJKKizXrMbXkPZ5bBEQMy0CEz2UHLqmKEuVTz/npDIDIY+PaqkNOTL5Z5RZTkLMY5OBHjOfNUHqbBgAg7yiEAACrUQgBAFajEAIArEYhBABYjUIIALBaScQnHNe1tGpYXoLfZZ2yV/Uc8PrLxH2p9eJua/Hlwn15iXS525LlqKu/y8H0Nwo7m1PINiDy6OhD49IRIns5LYPPYePcJpXtWefWLMOjXiayvUBQraQPePVus39qZf4i6fIsJAu+Zj6fZB8SHQArDHlPLsknlxNPvlnpEyynWQxy8CPG84ZPhAAAq1EIAQBWoxACAKxGIQQAWI1CCACwGoUQAGA1CiEAwGolkSM0xpExr6SKHcnEYCCoAzOhCpE/SyWqMwdjYbEHJ6kjdE4yoXYbyxxMxsJyD6VJzjYREbFL49K4avB8AXF3gQpxyGQqzq0hV0w15EqojWVuyq2Jk0vznew78uTSd0Y+NBn9dEn2eQPi9V42LJg5WFktBk2X6F9mjAmrxJ+cQUw9BPnjbhuXrL6EmG1AfdAo0MOqKxPvYCGf/qgjz5BkXL2tqS5dbhlilzS3moNs7uba8U0lhtWGA35q+UQIALAahRAAYDUKIQDAahRCAIDVKIQAAKvlYdVoNBrdtm3ba6+9FovF/vqv//rUb/3kJz95/vnnm5qabr311qampsHfFwAA+ZWHQvjSSy/ddttto0ePfvHFF08thA8++OBDDz20bNmyl19+ecaMGTt37iwr0wvxHUcvZJfLqv1qMW5ZhX4gFcP0PWaSa/FjfX1y40S0N3Mw3tetBk9mOYGSFevtlOMejzgQsjWS1y/W4stEhDEmWDU8czBUU5856POp+ISKSRhj+mJiFXhULQ2Xq9vd4hN+2dhIdgpTu/C4tbLSLb3E6nYnIQbd+MvFy6G8rjxzUD6NHpe1+L7OSOZgj9qDHOxUh8YY05ccSi3MetVsPeqR6X5JLmdYhTrJa1TnrDoVIauWAQ6jcwbJiAqARcWg21kn0wseFYpw5KDHLT6hNtY9owYYoMjDn0ZnzZr15ptv3n333acOJpPJ+++//5FHHlm0aNF3v/vdqqqqp59+evD3BQBAfhXqf4R79+7dt2/fzJkz0zcvvfTSl19+uUD3BQDAgBWqEHZ0dFRXVweDv/ubWH19/cGDB902Pn78eIGmAQDARytUIQyFQtHo/16QKRqNlpeLf0WkVVRUFGgaAAB8tEIVwpaWlkgkcvTo0fTNvXv3trS0uG0cCoUKNA0AAD5aoQphfX39jBkznnzySWPMiRMnVq9evWDBggLdFwAAA5aH+MShQ4cuuuiiaDQaDofHjx8/evToDRs2GGOWL19+9dVXr1u3bufOnVdccUVbW5vbHlKufQDEYJlaEDyiQl0s35ikSmDIJf56PXFCL+yO9ZwQgypmIDMVHw+OWr/sLxN/5S6rGZk5KBMRxphgpYhPBNUfz/1qvbhb94lwTCz4juj4hOp54nJhfZ+KSuhBuQPVO8IY43HUuFyzrhqh6DPcmECF+NNLRV2l3Dj73cqsRU/WmYrOuI5JDK3uE1L2D6DWJedQXybeorNPSoRqdH5MpsWSUXEuJcKi5UhK9akwxuj0g08Ny6SEW/cJNZ5S7+0DPmXyUAjr6uqef/75/92j/3f7vPjii99+++0tW7Y0NTVNmTJl8HcEAEDe5aEQ+v3+cePGyW/V1tZefvnlg78LAAAKhGuNAgCsRiEEAFiNQggAsBqFEABgtTwslskHR65Zl8r8onhXlwfkxjqVoQbjanm9r1cv502pVe/JaFgMxsXi45zI1g0e1YIjGRML1k8zX1As0C8bVps5GBreIPdQVq5aVaj14n51DX43sYQ66GoFtjwTZUMJt/Hsu0/oS+gbY1Ji3JFZC7UHj1yw7tJ9Ilgtoikp1Ush6dImItIlTnKvik9E1Utv8F0m5BMedDlmcoX96YxqyDNhmGt8QhzK4ZXi7S5YqV445fod3qveReVBl0kJuaVxaawh36yMGpQtKYwxjmp0kzJqtkXsPgEAwNBFIQQAWI1CCACwGoUQAGA1CiEAwGoUQgCA1SiEAACrlUaO0NHhLZlKkT1uAj5d0YMqLuNTG8tBt74zsr2KjnnlwusTh8NfVp456PGKLd1mm1ABxwLxqCSQ1y/ia/6APvdkOtCnjqPXLdynyJZeMnXk5NA8xyU4lf3Pu8ae5BmW9cRcnhmPXx0ddSB8QTXoEtyUoTQdm1SPYLApQmNCKpQW0o2vdI4wrl68BYoWymkFXF68ZepQ+kPqjSKkDm5QvzHKNkwub2vyacjluVEPzcnlVeJyMmU5lhU+EQIArEYhBABYjUIIALAahRAAYDUKIQDAahRCAIDVSiM+4cl+faxJyiZKLm1BYgnVTUa2mFGDcs29MXoFtEwO5CSVFB1PZPihZNswyQxJKiHa9CRUexfjkpSQp4fXJTMjyWyJbIjjySn+kPWg5pbP0WdY1hNzaSrkJNTRUQciqbqSJeM6IJRSrzIdiJLHUe40FxHVtSrlchzlc3P6mjC53Ffc5a1GNq4qj4ijI18ObpEqpzzbvJoOWuQUflAPzZPTU551ui6Xaf0BPhECAKxGIQQAWI1CCACwGoUQAGA1CiEAwGoUQgCA1UojPmE8ciG7FFVrtbvCcblxZ28sczDSKzaOhdV6cbXc3BjjVUkJn2oT4YuIxgvJuIgTuEklxEMoWTLCEe0+njnoljZJxYeL3VZUqD3ksFg66JdNS1R8Qu3VJY+gx10G5RJwl99EVTxGP2NqD05Sn7eJsDjxYl19mYORThHaiXbq8zauXlAyzyR7KVS4xGCiWbdzkU94JHk6MxE5kGdCd1wHwI54VeJFPbLqqNgyVCPef4wxgcpA5qAMYMjmJG7JJZ0mUuEWjx7UR9zjiI3lDLKvI9nsDQAAW1AIAQBWoxACAKxGIQQAWI1CCACwGoUQAGC1kohPeI3xq3XV8tLpUbXO+ESfzhj0dYkF333dYjCsghbJmF4v7g0EMweDlWLdf0olJZInc4hPlCyPWrifjIll95ETh8SWqquGMSZe1SU2rmnIHAwE6zMHvS6ZivKgONVDARFIkCuwUy4Xy5e9UPSg3IFPZ0gcjxj3+NWr1ScG3bqmxPtEuKXvWG/mYO8hMdh3VB+yaE+2CZ8q1Vok7pJN6VUNYWJuQZaSlH3jhhMu8YloSkRTulWErC4ozpkWl4lV1os8kq9MnEv+chHAkJkKVzLM46hBt8CMGvd6VK5joO0n+EQIALAahRAAYDUKIQDAahRCAIDVKIQAAKtRCAEAVqMQAgCsVhI5Qo9H5whVkxyTUP07on2iiZJxiQzKcGEsLIJTbu1sfIFQ5qCnUubqRHIrcvKo3G1pkvlIY0wqKRJOiYh4GmXbqeTJI3K38T6RI5Thq8q6kZmDXpVUM8ZUqJRVmc4Rih93CeaZhIq1Zd+GyXFrw6SyjB6veLXqcKEL2YYpfEykA2WOsFe9cIwxYZWRlIGuYerouAUDZbQulnVvptOsUnUmCqiTqUelABMuZ1h3QvVsUnuIqDdG+YQbY+rU4fGFRDzaVyYG3c46edAdmQJUgx4ZLjQ6dOhROUIZPc8GnwgBAFajEAIArEYhBABYjUIIALAahRAAYDUKIQDAaiURnzDGI3vf+NTiY7k2PR7T624jfWKJf7RHNfpR/YO8/oDcrTegmph4xVp8X7Bc7qE0ydn6Q+LBGpdnTMYnciKzFjJT4ahTwa0NU1CtI3dp/iV+3C0+kVTfkEkJl/ZBbqu91bh8aOqsc1tEnoqLlFG0WzRR6u0Sg50urYLkk1Ohwk8yTlAuY1LGBAfcU6cYKvxitn51IPrke9Wgu0sdi4r9Rsr1IZNniE81V/KqiJFxif3oB6FyHbrjkls2Rp1h8hwf8BnDJ0IAgNUohAAAq1EIAQBWoxACAKxGIQQAWI1CCACwWknEJzwuy149clhe2t9lcXpKXaY9lRAL9GUvBY/P5SLrenm67BgwlH7V8Mi1+Gow/Y3CzuYUqYQ4Oo5eVK1XUMt8TvZyWtyew8a5TUouGM960CVwkoqLNesxteQ9nlsERAzLQITPZQcuqYoS5VPP+ekMgMhj49qqQ05MvlnlFlOQsxjk4EeM581QepsGACDvKIQAAKtRCAEAVqMQAgCsRiEEAFitsIUwmUy+9dZbx44dK+i9AAAwYB65Bj0v3nrrrblz51ZVVe3fv//WW2/9h3/4B7ctp0+f/q1vfautra1AMwEAwE0BPxHecccdX/jCF954443XXnvtgQce+M1vflO4+wIAYGAKVQhPnDjxy1/+8i/+4i+MMWPGjLnyyit/8pOfFOi+AAAYsEJdWWbv3r2BQGD06NHpm+PHj//tb3/rtnEsFnv77bf9/t9Nprq6urW1tUATAwDgVIUqhD09PeXl/9vuvKKioqtLNBlP6+joWL58eWVlZfpma2srHx8BAKdHoQphY2NjV1dXMpn0+XzGmOPHjzc1NbltPGbMGBbLAACKolD/IzzjjDNqa2s3b96cvvnKK6+cd955BbovAAAGrFCfCIPB4M0333z77bc/8MADL7300rvvvvv5z3++QPcFAMCAFbAN07JlyyorK7/+9a83NjauX79+2LBhhbsvAAAGpoCB+uwRqAcAFAvXGgUAWI1CCACwGoUQAGA1CiEAwGoUQgCA1SiEAACrUQgBAFajEAIArEYhBABYjUIIALAahRAAYDUKIQDAahRCAIDVKIQAAKtRCAEAVqMQAgCsRiEEAFiNQggAsBqFEABgNQohAMBqFEIAgNUohAAAq1EIAQBWoxACAKxGIQQAWI1CCACwGoUQAGA1CiEAwGoUQgCA1SiEAACrUQgBAFajEAIArEYhBABYjUIIALAahRAAYDUKIQDAahRCAIDVKIQAAKtRCAEAVqMQAgCsRiEEAFiNQggAsBqFEABgNQohAMBqFEIAgNUohAAAq1EIAQBWoxACAKxGIQQAWI1CCACwGoUQAGA1CiEAwGoUQgCA1SiEAACrUQgBAFajEAIArEYhBABYjcJtMa0AAA0kSURBVEIIALAahRAAYDUKIQDAahRCAIDVKIQAAKv5B7+Lnp6edevW7dixIx6P33XXXf3jyWTyu9/97tq1a5ubm7/61a+OHz9+8PcFAEB+5eET4fbt27/5zW/u3Lnz/vvvP3V8+fLl3//+97/yla/U1dVdcsklfX19g78vAADyKw+FsL29fcOGDXfeeeepg/F4/KGHHvrOd75zxRVX3H333S0tLf/xH/8x+PsCACC/CvU/wr179x4+fPjiiy9O35wxY8bWrVsLdF8AAAxYVv8jjEQiBw4cyBw/44wzAoGA/JFDhw4NGzbM7//d/uvq6rZt2+a2/3379i1ZsqSmpiZ9c/To0T/84Q+zmRgAAIOUVSF8/fXXb7jhhszx1atXT5gwQf5IZWVlJBLpvxkOh6uqqtz2P3LkyC9+8YsTJ05M32xqaspmVgAADF5WhbCtre3dd9/Nab8tLS2xWOzgwYOjRo0yxuzevXvMmDFuG4dCoalTp7a1teV0FwAADF6h/kdYV1c3a9asFStWGGM6OjpWr1597bXXFui+AAAYsDzkCA8ePDhp0qRkMhkOh2traz/xiU+89tprxpj77rtv7ty5a9as2bVr1+LFi6dMmTL4+wIAIL/yUAgbGxt37drVf9Pr/d2nzClTpuzatevXv/51Y2PjR/xdFACAIspDIfR6vSNGjJDfSv/zb/B3AQBAgXCtUQCA1SiEAACrUQgBAFajEAIArEYhBABYjUKYf//8z/+cTCaLPYv8e/nll9evX1/sWeSf4zh33313sWdREGvXrt2yZUuxZ5F/sVjs3nvvLfYsCuLpp5/+zW9+U+xZ5F9XV9e3v/3tYs/CFYUw/5YvX97b21vsWeTff/3Xf/3qV78q9izyL5lMLlu2rNizKIjnnnvupZdeKvYs8u/EiRMf6n76sfHss89+LBv17Nu373vf+16xZ+GKQggAsBqFEABgNQohAMBqHsdxij0H09zc7PP5gsFgsSeSH3v27BkzZozH4yn2RPLs5MmTqVTK7XJ6Q9ru3bvPPPPMYs8i/44fP+73+6urq4s9kTxLpVL79u37WF7B+OjRo6FQ6CO6tw5R8Xj88OHDLS0tp/+uv/CFL9x1110fvU1JFMJDhw59nFaXRKPRsrKyYs8i/9JLYX0+X7Enkn8f10OWSCQ8Hg+HbAiJx+M+n6+/dcHHSbEO2ahRo8rLyz96m5IohAAAFMvH8PcOAACyRyEEAFiNQggAsBqFEABgtTx0qLfcgQMHXn311YMHD15++eWnrufetWvXj370o1gsdt11102ZMqWIMxyYY8eO/eIXv3jzzTerq6sXLFgwadKk/m8999xzzz33XFNT00033VRbW1vESQ7Ms88+++qrr3Z3d7e2ti5atGjYsGHp8Y6OjhUrVhw7duyzn/1se3t7cSc5GDt27NiyZcvChQuHDx+eHtm4ceMzzzxTV1d34403NjU1FXd6A/D4449HIpH012PHjp0zZ0766z179jz22GO9vb3XXnvt1KlTizfBgUsmkz/96U+3bds2YsSIBQsWTJ482RjjOM7PfvazzZs3n3nmmV/+8pf/z0WPpeaDDz54/vnnTx2ZP39++sR76623nnjiCWPM9ddff/bZZxdnfhn4RDhY06ZNu/fee5cuXfrGG2/0D+7Zs2fq1KmxWGz48OGXXHLJULx44B133PHss882NjaePHly2rRp/af1D3/4wxtvvHHs2LFvvPHGjBkzYrFYcec5AE888UQgEBg7duyqVaumT5+efoft7u5ua2vbvXv3mDFjrrnmmlWrVhV7mgPU3d19ww033HzzzYcOHUqPrFq1asGCBWecccaePXva2tq6urqKO8MBuOOOO7Zt2/b++++///77hw8fTg92dHRMnTq1s7OzoaFhzpw5GzduLO4kByCZTM6dO/ehhx5qaWlJpVKvvPJKenzZsmXf+MY3Wltb16xZM2/evOJOcgB6enre/73169f/5V/+ZToQ8tZbb02fPj0QCJSVlV144YU7d+4s9kx/z8HgJJNJx3HGjx//85//vH/wzjvvXLRoUfrrr3/96wsXLizO5AYhHA73f7106dJrr73WcZxUKtXa2vr000+nvz7nnHN++tOfFm2KgxaPx6uqql555RXHcR555JEZM2akx1esWDF16tSiTm3gvvKVrzz88MPGmLfffjs9Mm3atH/7t39Lf93e3v7www8Xb3YDVF9f3/9w+n3jG9+46qqr0l/fe++9V1555Wmf12D94Ac/mDRpUiwWO3Wwp6enpqZm+/btjuNEIpG6urpNmzYVaYJ58NWvfnXBggXpr2+++eZbbrkl/fVtt9120003FW9ef4BPhIMlo68bN2789Kc/nf56zpw5GzZsOL2TyoNQKNT/dTQaTV/q4uDBg++99176D1Mej2f27NlD8aH12759u9frHTdunMk4ZFu3bu3r6yvq7AZiw4YNr7/++pe//OX+kUgksmXLlv6/JQ7Rs9EY8+STTz744IOnNtP4GLzKVq9evWjRol/+8pcPPPDAf//3f6cHd+zYEQgEzjvvPGNMWVlZe3v7UHxoaYlE4oknnrjxxhvTNzdu3Nh/Kn76058uncdFISyIgwcP1tfXp79uaGg4cuRIIpEo7pQG7PXXX3/88cf/6q/+yhhz8ODBioqK/us/NTY2HjhwoKizG6DFixc3NzfPnDnzxz/+cUNDg8k4ZOmRYk4xd319fbfccssPfvCDU385Sz+K9CMyQ/aQXXjhhdFodPfu3QsWLLj99tvTgx86ZL29vUPur74ffPDBo48+unLlyp6enmuvvfbBBx80xnR0dPQfLzNkD1naL37xC6/Xe/nll6dvfuiQlc5LjMUyBeH3+/srXyKRGLrXTNqzZ89VV111//33n3POOcaYQCBwakWPx+ND9DJX3/nOd7q6ulavXr148eItW7aMHz/+Q4fMGDPkLn77d3/3d9ddd93ZZ599al/oQCBgfv+IzJA9ZM8880z6i1tuueWP/uiPbr311nHjxmUesvSDHUK8Xm9ra+vjjz9ujGlra1u4cOFtt9126uMyxsTj8aF7tdgVK1b86Z/+af9F/j50yErnJTYk351LX0tLS/8vcfv37x81atRQLIR79uyZOXPm3/zN33zpS19KjzQ3N8disSNHjqRvph9a8SY4cFVVVc3NzTfddNO55567evVq84eHbN++fT6fr7GxsahzzNmTTz65cuXKCy64YNq0acaYa665Jr3cyefz7d+/P73N0D1kaa2trfX19R988IHJeJXV1tYOudWVLS0t/SsnJ02adPLkyZMnTzY3N3d0dPT/NjN0D1lHR8dzzz23ePHi/pEPHbLm5uYiTe3Dht6785Awb968lStXOo5jjFm5cuVQXPe1b9++WbNmLV269Oabb+4fHDly5IUXXrhy5UpjTF9f3+rVq+fPn1+8OQ5ENBpNpVLpr7u6ut5999106GXevHmrVq1KL4L9z//8z8suu6x0fl3N0tq1a1esWPG9733vX//1X40x//iP/3jRRRcFAoHLL788fcji8fgzzzwz5A5ZJBJxfn9J5C1bthw9enTixInGmHnz5j311FPpozlEX2VXX3315s2b049u06ZNzc3Nw4cPP++886qrq9euXWuMOXTo0IsvvviZz3ym2DMdiMcee+zCCy+cMGFC/0j6jTH99c9+9rMSOmRFXqwz9C1evPj8888vKyubMGHC+eefv3PnTsdxTpw4MWnSpFmzZl1zzTUtLS27d+8u9jRz9id/8iehUOj837vhhhvS4+vWraurq/viF784ZcqU+fPnp1Kp4s4zV5s3b04HJK677rqmpqaFCxem1/3GYrGZM2dOmzbt+uuvHzly5NatW4s904FL//Wpf5nltm3b6urqrr/++ra2tpkzZ35ojWLpW7Nmzfjx4z/3uc999rOfHTZs2L/8y7+kx3t6ev74j/+4vb194cKFjY2N77zzTnHnOQDhcPjiiy9ub2//0pe+VF9fn16P7TjOE088UV9fv3jx4gkTJixZsqS4kxywiRMnPv7446eOdHR0jB8//sorr/zMZz4zbty4gwcPFmtuH0L3icF65513enp6+m9OnDixsrLSGBOJRF544YVoNDp79uyampriTXCAdu3a1dnZ2X+zsrIy/Wu4MWb//v0vvvhiY2PjJZdcMhT/5Pvee++9+eabyWTy7LPPPjXSm0gk1q9ff/z48UsvvfTU1QpDjuM427dvnzRpUv/S38OHD69fv762tvbSSy/1+4fYyoBkMvnGG2+8++676d/MRo8e3f+taDS6bt26np6e2bNnD9FOmemzrqura/r06ae263vvvfe2bds2duzY6dOnF3F6A5ZIJF5//fXJkyd/6H/Svb29L7zwguM4s2fPLp22ixRCAIDVht6v8wAA5BGFEABgNQohAMBqFEIAgNUohAAAq1EIAQBWoxACAKxGIQQAWI1CCACwGoUQAGA1CiEAwGr/H6umurnGGsDgAAAAAElFTkSuQmCC", "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "execution_count": 9 } ], "cell_type": "code", "source": [ "using Statistics\n", "using Plots\n", "mode_xy = mean(real.(X_Simple[1]), dims=3)[:, :, 1, 1] # Average along z axis\n", "heatmap(mode_xy', c=:RdBu_11, aspect_ratio=1, grid=false,\n", " legend=false, clim=(-0.006, 0.006))" ], "metadata": {}, "execution_count": 9 }, { "cell_type": "markdown", "source": [ "This mode can be physically interpreted as the reason why this SCF converges\n", "slowly. For example in this case it displays a displacement of electron\n", "density from the centre to the extremal parts of the unit cell. This\n", "phenomenon is called charge-sloshing." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "We repeat the exercise for the Kerker-preconditioned dielectric operator:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "┌ Info: Arnoldi eigsolve finished after 1 iterations:\n", "│ * 3 eigenvalues converged\n", "│ * norm of residuals = (2.12e-04, 4.21e-04, 2.33e-04)\n", "└ * number of operations = 12\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Plot{Plots.GRBackend() n=1}", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de3zU1Z3/8ZOZ3O/3K4lAEkRuohDui1jAa7GK1dqq3aJ23T5cL7tYH7uP7nbturo8dKu2VbttV+u62nZLddH2B4oWlrsQQK6iwUAg5AohkHsmmZnfH9Om1O97JCQThnhez7+SD5PzPfP9nu/5ZMj5fE+E3+83AADYyhXuDgAAEE4kQgCA1UiEAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAaiRCAIDVLohE+OSTT9bX14e7FwAAG10QifCNN944cuRIuHsBALDRBZEIAQAIFxIhAMBqJEIAgNVIhAAAq5EIAQBWIxECAKxGIgQAWI1ECACwGokQAGA1EiEAwGokQgCA1UiEAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAaiRCAIDVSIQAAKuRCAEAViMRAgCsRiIEAFiNRAgAsBqJEABgNRIhAMBqJEIAgNVIhAAAq5EIAQBWIxECAKxGIgQAWI1ECACwGokQAGA1EiEAwGokQgCA1UiEAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAaiRCAIDVSIQAAKuRCAEAViMRAgCsRiIEAFgtcvBN1NbW/va3v62srIyPj1+4cOHs2bP7/mn16tXvvfdednb2Pffck5qaOvhjAQAQWiH4RPj+++9/8MEH+fn5Xq/3hhtu+OlPfxqIv/jii0uWLCkoKCgvL587d25vb+/gjwUAQGhF+P3+EDb3wx/+8PXXX1+3bp3P5xszZszTTz99ww03+Hy+CRMmPPbYYzfffLP8qRkzZvzgBz+YPn16CHsCAEB/hPJvhD09Pdu3b58wYYIx5tixY4cOHVq4cKExxuVyLViwYN26dSE8FgAAIRGCvxEaYz766KPrr7++oaFh8uTJq1evNsbU19cnJibGxcUFXpCdnb1r165gP97Y2Lhs2bLs7OzAt3l5eY8++mhIOgYAwGcLzSfC0tLS8vLyDRs2xMXF3X///caYyMhIr9fb94Le3t7o6OhgPx4bGzt27Ngpf3TppZeGpFcAAJxVaD4Rut3u9PT09PT0J5544oorrnjxxRfz8/M7OjpOnToVWCxaU1OTn58f7MeTk5NvvPFG/kYIADj/QvCJsLOzs+/rnTt3FhYWGmNyc3OnTp26fPlyY0x7e/uqVasWLVo0+GMBABBaIfhE+I1vfOP48eMXXXRRTU3Njh07fvGLXwTiTzzxxFe/+tUtW7bs2rVr6tSpc+fOHfyxAAAIrRCUT3R2dr7//vs1NTUZGRmzZs1KSUnp+6ejR4+uX78+Pz9/3rx5LlfQT5+UTwAAwiUEnwjj4uKuvPJK+U9FRUV33HHH4A8BAMAQ4VmjAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAaiRCAIDVSIQAAKuRCAEAViMRAgCsRiIEAFiNRAgAsBqJEABgNRIhAMBqJEIAgNVIhAAAq5EIAQBWIxECAKxGIgQAWI1ECACwGokQAGA1EiEAwGokQgCA1UiEAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAaiRCAIDVSIQAAKtFhrsDxhgzbXTepPr3O998/1PxiLhE54sjs0c4gxEZ+bJlb2K2M9gdk+wMNnd6ncGmzl7Z7LGWLtVCjzPockU4g9kJ0bLZgqQYZzAtVlyjtDi3MxjZdkI262ptdAZ9x6udQe+JOmfQ3+vRzSaliT7kFIkWUnPFsZLEpTHGtBtxEpq7xNU53iFOeE1Lt2y21SMuZWyk+EUwL1F0IDdRX7K0WHEhUqL8zqC8Cm4VNMb01h91Br3N+sWi2TR9biNzxdWRF8Kngqd7xGA2Qa5OfZsYNnVt4up09fpks0nRYuQXJIurkxUf5QzKS2OMSTCiD/JCRJyqdwZ7G8SlMcb4WptFC5Fi2Lgz85xBV1ahblZdiN7ETGdQzmDNXXoGq2kVJ6GxXVwyn08M5rQ4ccKNMSOSY53BjLj+zmAx3S2yWXebuDr+plpnsLfxmDMYe9Xdstkz8YkQAGA1EiEAwGokQgCA1UiEAACrkQgBAFYjEQIArEYiBABY7YKoI+yoPVHx6tvOeHyWqCPMmlzqDCZcOk22HDEywRnsjUpyBk+qgpvttadls5s+aXIGqxpanUGXW/yqMaEwVTY7e3S6M3hpruhtqqqRiujUve09vM8ZPL1zhzPYtP+wM9ijyvWMMUmFopgpe+pYZzB24kxnMCIuRTbbHSFKr2pbRYXTtmOnnMEtn+hiyuPNnc5gfLw41tRR4irMHCmKJo0x47LEAEtxi4oud7sYM90HymWzTdv3OIOnPhGFU1JqiS6rzZg6yRmMGT/DGfQlZDiDHb26Mq9SndstVaKubvvhk6LZDl2ompUW5wzOLBGjbtoIcUPFR4maNmNMol+U0EU01ziDXXu3OION2z+SzbZWi4EXpQocM8aPcgZTLp8im3WXirg/QZyEFo8Ydfsa22Szmw6JC7GvWtxQPq+o8hyZIyYlY8zsEjFsygrEnZ4UIybGWE+7bNZXW+kMtu/e5gwe33XQGSyljhAAgM9GIgQAWI1ECACwGokQAGA1EiEAwGokQgCA1S6I8onqhrYXN4jtMy5Si49nfSi25LgkQm8Qk5RzkTPoic1xBj86IVbuLi8XvTLGHNj6iTPYdFAshXerNdzVk2bLZk93iF1yMtUS/8JktStQs15e37RZrAL/8JeitxsPifXTp3rEsmxjzES15UqZWkReGhPvDEZmjpTNdrhEQcKuerE/y+sbq5zBynJReGCMaan52BmMTRP7Q9VeLkpxvL6RstmCJHES8qJEKU5PdYUzWLN6k2x233JR8bK5qUO+2GlWxhEZn3BSVPgUJotyEZM52hk7LTYfM8aYjYdEZcj/21DlDB7dKZa8dzaL3Y6MMSkFFzuDtWWiAiR6nuhtfpLeOSvTJ06jp2K3M1j5vxudwfJ3q2Sze9XWbKlRouBkzihxn47r0jUk2WrPJn9msTMoNyBbuU+f281bxGZS9XvEaPT2iPdVWVomm62dXuIMJsWILJOXKOZ2V5ACsNY9YrI68AvR2807xEZyS1+Qrf75oc/+EgAAPr9IhAAAq5EIAQBWIxECAKxGIgQAWI1ECACw2gVRPuHzG4/P74wfbBNLilveF8+JTxstlpsbYxKnXeEMepLF/hV7a8UC/YqdYpGxMaZh7zoZd+rtEmu1j2z6rXxxXMLNzuBflGY5g9MLxL4cvXV60fzRNQecwd98qHdp6L9yteeAUYvLM8aJB/bnTl0gm+1wiwfYb6wQva3YLBZVt9SIKoVg2hvF9a1YI+pzktNvkS3MLRa9vSReDOauw6J+45OVImiMebtBbxrQT8F+PFYdLnvKGPHKSfOdwbo2URZijNmwv8EZrFjzO2fQ6wlSgaE0V+0VzfaIIoGN+cnO4LyRqizEGJdXFEqd3COGqKyU0MM+iHqvOGO/OSAG8zfUTWqMSS8TN7W/ZI4zeFh17IM9unwi2BTUT8EmQHeUKFnZqzbVmV0otqTwnz4um63fIqb3t1UiaOjWQ/Ss+EQIALAaiRAAYDUSIQDAaiRCAIDVSIQAAKuFLBE2Nzd7PGKR5/Hjx7u7xUIvAAAuBCEon3j++ef/9V//taenp7Oz89prr3355ZcTExONMVVVVTfddFNjY2N7e/tjjz12//33D/5YJsgC2cb9et3t6JZmZ7BXlWpU1IsH8zfu15sDDJHGykpnsLpZPGlevAFjvCf1UunqICdnKMjF5TM+EQud89Q6eGNMq1dsdnFUFbecU6VE/8kl/o3VogPGmKaOHhH1ib1QWo+KGoMddYMqkzhX8nATVcdifeIqNHXIcadPzjlVSvSfvOhHa8VmCK0evWtKhE8MvFNqiJ5TpcQgBbtJJ6qbWl6G6mZRqSWnlKEjJ8yKOWIGk5OwT03XJsj0PuBKCSkEnwjz8vLWr19/4sSJ2traY8eOPf7444H40qVL586dW1NTs3Xr1u985zuV5/eSAADQHyFIhIsXLy4tLTXGpKSkXH311QcPHjTGnD59+q233nrooYeMMRdffPHVV1/92muvDf5YAACEVigXy3g8nhUrVsyfP98YU11d7XK5Ro0aFfiniy+++PDhw8F+0Kv+NwwAgPMgZInQ7/fff//9GRkZf/VXf2WMaWlpiY//077kCQkJp06J3c8DTpwY7OO+AAAYmJAlwqVLl+7atWvFihVut9sYk5WV1dra6vP5Av966tSp7OzsYD+bk5MTqm4AAHBOQpMIv/Od76xdu3bVqlXJyX94+m1RUVFiYuLOnTsD35aXl0+cODEkxwIAIIRCUD6xbNmyZ5999vvf//7vf/97Y0xGRsYXvvCFmJiYu+6665FHHnnhhRc2bty4a9eu119/ffDHAgAgtEKQCP1+//XXX79mzZrAt2PGjPnCF75gjHniiSf+6Z/+6Stf+UpeXt6qVavS0tIGf6xguk7pojS/KlZTFSympVMUhPV2i9KcodPd0uQMtve7XMbbqU/CySAFVedNd6t40oJRlWrGmB6/T7Qgy/XOo+52PRK6e0VvI/ziP1p620Vd3QlPKGuhzkoeTnYswi9uEvlmTfCTc97I4dHj1b01fjHw9BA9j4LdpPKmFmWqQSYKOaUMHTlhyqlVTsJyujbBp/cQCkEi/Id/+AcZj42Nfeqpp5566qnBHwIAgCHCs0YBAFYjEQIArEYiBABYjUQIALAaiRAAYLUQrBq9EMSmxsh4RJSIu9Tq4+S4KGcwMibeGTRDVlYRk5zhDCbE9PcaueP0SUiPdg+8T6EQkxQtoi7dqyiX+OUsJl5cnfMpJkGPhJhI0Vt/hBhhkQmxzmBmtL641WrF+eDJw8mOybcg36wJfnLOGzk8otxBfsv3iYGnh+h5FOwmlTe1rAuRE4WcUowx7cer+9+3/pMTppxa5SQsp2sTfHoPIT4RAgCsRiIEAFiNRAgAsBqJEABgNRIhAMBqJEIAgNWGX/lEjlolnD0+S77YlSy2vIhUS3fH5CY5g+XjZ8tma3e++1ldHKjs4mJnsDBNrEiWj593p+fKZgvlyWlsP5eu9VdZWpwzmFpS4Az6gyyVTnKLdeRF+cnO4MGCMc5gS03FWbp4Nu5oUU6QXSg6YIzJkHUdLvFo/aQisfv0lLxE2Wz1oebP6OGAycPJjsnilox4Oe70yTmqTqPXI3a6OCfJ6qLL4ZEUpCDB7xUDTw7RsrRKZ7C8ufMsXRwQfZMGuanVzg16opBTijHmZOWu/vet/7LVhCmnVjkJy+naBJnec3Y3OoMN/d6r59OHHtiPAQDw+UAiBABYjUQIALAaiRAAYDUSIQDAaiRCAIDVLojyCVeEiVaraS9Sa9NnTclzBnNnTpAtR6SIdbfRbnGsiWoF9pjLi2Sz3p4rnMGmg+XOoDtKLCLPnaSrMi6bJJZKj1IFCRFyf4O8i2SzRV+4xBn8ckObM7jx0Cln8FSPVzY7MVm8tbKFI53B9EljnUFfdIJsNl7tPjFnTKYzWD2rzBmsLNdVGadrPnYG49LECS+6fJoz+BfjVY2BMXmJausGt1jDHTPqYmew5Dpd7HHN8n3O4Oam/u55MitDbwdRcp3og+yY3y1uvbxEPV3Ik9Ny8ovO4NGd25zBzuZ62WxKgehYcdkkZ1AOj/go/Vu+zy0GnhyiZQuPip9/t0o2u7dFVIakRokSjjmjUp1BeZOaIDe1V93+cqKQU4oxprN9kTNYv2eTOFaPeF8ZpeLWM0EmTDm1yklYTtcmyPR+jZqsNu+oky2cFZ8IAQBWIxECAKxGIgQAWI1ECACwGokQAGA1EiEAwGokQgCA1S6IOsLCnMS7J4h6l/gssWtM1uRSZzBhkq5r8cWlOIOyhGVspqguuqVshGx2kyrZqbpspDPocotfNSYUikIiY8zs0enOYEGyKIzT2+Gk5ctwxqyZzuCk2GhxrP2HncGejh7ZbFKhKN7KnirqsaLHXOoMeqN1rVt8hDhjk3NFKZJnzkhncEuu3tjoePNl4ljx4iRMHSWuwsyReoOYlFi114/69TKqUOwfVHCVriiNTRc715R8Uitf7JRaEmQkTBVFeLJjPS4xM+g3a8yc0RnOoFvVg24vEiO/o8Mjm81Sd9nMEjHq5PAIVkfo94uBJ4do8U2icFMOe2NMWfUJZzBKVUJnjB/lDKZcPkU2K29qefvLieI6Na8aY1LUyN+nyjF9Xp8zODJHjE9jzOwSMRLk1ConYTldmyDT+yV+sRtV5riDsoWz4hMhAMBqJEIAgNVIhAAAq5EIAQBWIxECAKxGIgQAWO2CKJ+Iz88c85V5znhEnFgKH5ktShoiMvR6ca/a6ydSbfmUHitOxdR8vZw3N1GsVG7uFGUGLnWs7ASxdtkYU5Akmk2OFmvW1TYsxh9k8XHkKLGJSVqieHHyOLGJib9Xr253JYmKgsgcsQ+LP1Ws4fZH6v2SYtTi8PwkccZmXyQ6MDJNV2W0esTWSLGR4hfBPHVxcxP1JYuPVFfCJS6ZN0GsLI+5RJf9ZKeJjY0ypjbKFzu507JlPDJXXB3ZMfkW9Js1pljVOSSoHYjG54ll9129YoG+MSYpWtySskggS1UpxKgF+sYYvxEt+NMKnMHYSXOcwXw1wo0xvtZmZzAiUgwbd6bYSM6VVaibVTe1vP3lRDEhW1cTpcWJMzZN1Qj5fKJKQf64MWaE2ppNTq1yEvYH2ZrNnV/sDCbGiFEXO1K8sj/4RAgAsBqJEABgNRIhAMBqJEIAgNVIhAAAq5EIAQBWuyDKJ1yxCZEjSkQ8RizGNYlqjW+QrQz8bvEG5apq+TT0JLUi2RhToFYJp6hVwi610jklRp92+bz8qCCrwJ38br2m2R8jllDLNdwRakWy6dW7T0TEqyfQJ4jtBfxR6jqqXSaMMW51xmLVwv3Ufp9wY0xnj9cZjFZ7g8jrGKcKLUywVeDqlbJcxBcrtk0wxkTmiBohV4J+5L94paqNCXY42TG/ukvkmzVBTk6WqhGSJ9yj9jcwxsSpAozkGBGUw8MdpLfGL/ogh2iEmmoigwwwf7LYtMREilvSlSTuEZ+6SU3wm9pJThTBtuCQNULx6oT71CYPiaqyxQSZMOXUKk+inK5NkOndlZrlDEZFDjCj8YkQAGA1EiEAwGokQgCA1UiEAACrkQgBAFYjEQIArHZBlE/4utp7j33ijAfZfUJsI+AKVjkQLVrwqzft8YpVwq0esebeGFPT0uUMDtHuEzFqxbkU4Q1S59Dd5gx6T4iNJmTwHHefEKcxQm40EasXzXvVcu2uXhE81SVGQk1Lt2z2HHaf8IreuoPsPiEX7ke41UnoFR1zdbXIZnsbjjmD3ubB7z7R78oBVQMSZJcI06n+4Xi7GDZ1beIkDH73iSiXuP3jImUZizFGHC6iR9zRpk1sKNHbcFQ3OtjdJ/RUHBGkMMypR81gHT363MoL0agu2TntPhGh6spk+YS8NhFecZMaY1yeDtHCqePOYG+juHHcYvcdxyHO/hIAAD6/SIQAAKuRCAEAViMRAgCsRiIEAFiNRAgAsBqJEABgtQuijrCj9kTFq2874/FZogowa3KpM5hw6TTZcsTIBGewN0psZ3NSFaVtrz0tm930SZMzWNXQ6gy6VBXghEKxD4sxZvZosZPLpbmit6mxYruTiE7d297D+5zB0zt3OINN+w87gz0dujwxqTDTGcyeOtYZjJ040xmMiNNbBXVHiNKr2lZR4bTt2ClncMsnJ2Szx5s7ncH4eHGsqaPEVZg5UhRNGmPGZYkBluIW5afudjFmug+Uy2abtu9xBk99Uitf7JRaki/jGVMnOYMx42c4g76EDGewo1fvSlapzu2WKlFXt/3wSdFshy5UzUoTm4LNLBGjbtoIcUPFy82/jEn0ixK6iOYaZ7Br7xZnsHH7R7LZ1mox8KLiRb1dxvhRzmDK5VNks+5SEfcniJPQooqe9zWKGmJjzKZD4kLsqxY3lE9tkjUyR+8INrtEDJuyAnGnJ8WIiTHW0y6b9dVWOoPtu7c5g8d3HXQGS6+6WzZ7Jj4RAgCsRiIEAFiNRAgAsBqJEABgNRIhAMBqIVg12t3dvX379g8++MDj8fzd3/3dmf/0y1/+8t13383NzX3ggQdyc3MHfywAAEIrBIlw48aNDz744IgRIzZs2HBmInz22Wefe+65Rx99dNOmTXPmzNm/f39MjNqOx5jqhrYXN4jtMy5Si49nfSg2o7kkQuz0YYxJyrnIGfTE5jiDH50QK3eXl4teGWMObBWbRjUdFEvh3WoNd/Wk2bLZ0x1FzmCmWuJfmKx2BWrWy+ubNotV4B/+UvR24yGxfvpUj96LaqLacqVMLSIvjRFbyURmjpTNdrhEQcKuerFd0esbq5zBynJReGCMOV3zsTMYlyZ+Oau9XJTieH0jZbMFSeIk5EWJUpye6gpnsGb1JtnsvuWi4mVzk9iMRpqVcUTGJ5wUFT6FyaJcxGSOdsZOq62KjDEbD4nKkP+3ocoZPLpTLHnvbK6XzaYUXOwM1paJCpDoeaK3+Ul656xMnziNnordzmDl/250BsvfrZLN7lVbs6VGiYKTOaPEfTquS9eQZKs9m/yZxc6g3IBs5T59bjdvEZtJ1e8Ro9GrNqiqLC2TzdZOL3EGk2JElslLFHO7K0gBWOseMVkd+IXo7eYdYiO5pS/IVv/80Gd/ydnMnz9/3759jz/++JlBr9f79NNPv/DCC3fcccePf/zjxMTEN954Y/DHAgAgtIbqb4TV1dXHjh2bN29e4Nsrr7xy0yb9yy8AAGE0VImwvr4+OTk5OvoP/zuRlZVVVyc+tAacPCkecwAAwHkwVIkwNja2u/tP/2fd3d0dFycemBQQHy/+hgQAwHkwVImwoKCgq6vrxIk/LJ2orq4uKCgI9uLYWP1UQAAAhtpQJcKsrKw5c+a89tprxpjm5uaVK1cuXrx4iI4FAMCAhaB8oqGhYdasWd3d3Z2dncXFxSNGjFi3bp0xZtmyZTfddNOaNWv2799/7bXXTp8+PVgLPr/x+PzO+ME2saS45X3xnPi00WK5uTEmcdoVzqAnWexfsbdWLNCv2CkWGRtjGvauk3Gn3i6xVvvIpt/KF8cl3OwM/kVpljM4vUDsy9FbpxfNH11zwBn8zYd6l4b+K1d7Dhi1uDxjnHhgf+7UBbLZDrd4gP3GCtHbis1iUXVLjahSCKa9UVzfijWiPic5/RbZwtxi0dtL4sVg7jos6jc+WSmCxpi3G/SmAf0U7Mdj1eGyp4wRr5w03xmsaxNlIcaYDfsbnMGKNb9zBr2eIBUYSnPVXtFsjygS2Jif7AzOG6nKQoxxeUWh1Mk9YojKSgk97IOo94oz9psDYjB/Q92kxpj0MnFT+0vmOIOHVcc+2KPLJ4JNQf0UbAJ0R4mSlb1qU53ZhWJLCv/p47LZ+i1ien9bJYKGbj1EzyoEiTAjI+Pdd9/9U4uRf2hz9uzZH3300bZt23JzcydPnjz4AwEAEHIhSISRkZGjR4uCVmNMenr6NddcM/hDAAAwRHjWKADAaiRCAIDVSIQAAKuRCAEAVgvBYpnzTC6Qbdyv192Obml2BntVqUZFvXgwf+P+8/p81MbKSmewulk8aV68AWO8J/VS6eogJ2coyMXlMz4RC53z1Dp4Y0yrV2x2cVQVt5xTpUT/ySX+jdWiA8aYpo4eEfWJvVBaj4oagx11gyqTOFfycBNVx2J94io0dchxp0/OOVVK9J+86EdrxWYIrR69a0qETwy8U2qInlOlxCAFu0knqptaXobqZlGpJaeUoSMnzIo5YgaTk7BPTdcmyPQ+4EoJiU+EAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAasOvjlDqOqWL0vyqWE1VsJiWTlEQ1tstSnOGTndLkzPY3u9yGW+nPgkngxRUnTfdrWI7LaMq1YwxPX6faEGW651H3e16JHT3it5G+MXvl73toq7uhCeUtVBnJQ8nOxbhFzeJfLMm+Mk5b+Tw6PHq3hq/GHh6iJ5HwW5SeVOLMtUgE4WcUoaOnDDl1ConYTldm+DTewjxiRAAYDUSIQDAaiRCAIDVSIQAAKuRCAEAViMRAgCs9jkpn4hNjZHxiCgRd6nVx8lxUc5gZEy8bHaIyipikjOcwYSY/l4jd5w+CenR7oH3KRRikqJF1KV7FeUSv5zFxIurcz7FJOiREBMpeuuPECMsMiHWGcyM1he3Wq04Hzx5ONkx+RbkmzXBT855I4dHlDvIb/k+MfD0ED2Pgt2k8qaWdSFyopBTijGm/Xh1//vWf3LClFOrnITldG2CT+8hxCdCAIDVSIQAAKuRCAEAViMRAgCsRiIEAFiNRAgAsNrwK5/IUauEs8dnyRe7ktOcwUi1dHdMbpIzWD5+tmy2due7n9XFgcouLnYGC9PEimT5+Hl3eq5stlCenMb2c+laf5WlxTmDqSUFzqA/yFLpJLdYR16Un+wMHiwY4wy21FScpYtn444W5QTZhaIDxpgMWdfhEo/WTyrKcQan5CXKZqsPNX9GDwdMHk52TBa3ZMTLcadPzlF1Gr0esdPFOUlWF10Oj6QgBQl+rxh4coiWpVU6g+XNnWfp4oDomzTITa12btAThZxSjDEnK3f1v2/9l60mTDm1yklYTtcmyPSes7vRGWzo9149nz70wH4MAIDPBxIhAMBqJEIAgNVIhAAAq5EIAQBWIxECAKx2QZRPuCJMtFpNe5Famz5rSp4zmDtzgmw5IkWsu412i2NNVCuwx1xeJJv19lzhDDYdLHcG3VFiEXnuJF2VcdkksVR6lCpIiJD7G+RdJJst+sIlzuCXG9qcwY2HTjmDp3q8stmJyeKtlS0c6QymTxrrDPqiE2Sz8Wr3iTljMp3B6lllzmBlua7KOF3zsTMYlyZOeHMzwHgAABRtSURBVNHl05zBvxivagyMyUtUWze4xRrumFEXO4Ml1+lij2uW73MGNzf1d8+TWRl6O4iS60QfZMf8bnHr5SXq6UKenJaTX3QGj+7c5gx2NtfLZlMKRMeKyyY5g3J4xEfp3/J9bjHw5BAtW3hU/Py7VbLZvS2iMiQ1SpRwzBmV6gzKm9QEuam96vaXE4WcUowxne2LnMH6PZvEsXrE+8ooFbeeCTJhyqlVTsJyujZBpvdr1GS1eUedbOGs+EQIALAaiRAAYDUSIQDAaiRCAIDVSIQAAKuRCAEAViMRAgCsdkHUERbmJN49QdS7xGeJXWOyJpc6gwmTdF2LLy7FGZQlLGMzRXXRLWUjZLObVMlO1WUjnUGXW/yqMaFQFBIZY2aPTncGC5JFYZzeDictX4YzZs10BifFRotj7T/sDPZ09MhmkwpF8Vb2VFGPFT3mUmfQG61r3eIjxBmbnCtKkTxzRjqDW3L1xkbHmy8Tx4oXJ2HqKHEVZo7UG8SkxKq9ftSvl1GFYv+ggqt0RWlsuti5puSTWvlip9SSICNhqijCkx3rcYmZQb9ZY+aMznAG3aoedHuRGPkdHR7ZbJa6y2aWiFEnh0ewOkK/Xww8OUSLbxKFm3LYG2PKqk84g1GqEjpj/ChnMOXyKbJZeVPL219OFNepedUYk6JG/n5Vjun1+pzBkTlifBpjZpeIkSCnVjkJy+naBJneL/GL3agyxx2ULZwVnwgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAahdE+UR8fuaYr8xzxiPixFL4yGxR0hCRodeLe9VeP5Fqy6f0WHEqpubr5by5iWKlcnOnKDNwqWNlJ4i1y8aYgiTRbHK0WLOutmEx/iCLjyNHiU1M0hLFi5PHiU1M/L16dbsrSVQUROaIfVj8qWINtz9S75cUoxaH5yeJMzb7ItGBkWm6KqPVI7ZGio0UvwjmqYubm6gvWXykuhIuccm8CWJlecwluuwnO01sbJQxtVG+2Mmdli3jkbni6siOybeg36wxxarOIUHtQDQ+Tyy77+oVC/SNMUnR4paURQJZqkohRi3QN8b4jWjBn1bgDMZOmuMM5qsRbozxtTY7gxGRYti4M8VGcq6sQt2suqnl7S8nignZupooLU6csWmqRsjnE1UK8seNMSPU1mxyapWTsD/I1mzu/GJnMDFGjLrYkeKV/cEnQgCA1UiEAACrkQgBAFYjEQIArEYiBABYjUQIALDaBVE+4YpNiBxRIuIxYjGuSVRrfINsZeB3izcoV1XLp6EnqRXJxpgCtUo4Ra0SdqmVzikx+rTL5+VHBVkF7uR36zXN/hixhFqu4Y5QK5JNr959IiJePYE+QWwv4I9S11HtMmGMcaszFqsW7qf2+4QbYzp7vM5gtNobRF7HOFVoYYKtAlevlOUivlixbYIxJjJH1Ai5EvQj/8UrVW1MsMPJjvnVXSLfrAlycrJUjZA84R61v4ExJk4VYCTHiKAcHu4gvTV+0Qc5RCPUVBMZZID5k8WmJSZS3JKuJHGP+NRNaoLf1E5yogi2BYesEYpXJ9ynNnlIVJUtJsiEKadWeRLldG2CTO+u1CxnMCpygBmNT4QAAKuRCAEAViMRAgCsRiIEAFiNRAgAsNrQJkKv13vgwIGmpqYhPQoAAAMW4VerY0PiwIED119/fWJiYk1NzQMPPPDP//zPwV45Y8aMH/zgB9OnTx+ingAAEMwQfiJ8+OGHv/a1r+3Zs+eDDz545plnPvzww6E7FgAAAzNUibC5ufntt9/+67/+a2NMUVHRdddd98tf/nKIjgUAwIAN1ZNlqquro6KiRoz4wwMyiouLjx49GuzFHo/no48+ivzjQwGSk5NLS0uHqGMAAJxpqBJhW1tbXNyfntcVHx/f0tIS7MX19fXLli1LSPjD9sSlpaV8fAQAnB9DlQhzcnJaWlq8Xq/b7TbGnDx5Mjc3N9iLi4qKWCwDAAiLofobYWFhYXp6+tatWwPfbtmy5bLLLhuiYwEAMGBD9YkwOjr63nvvfeihh5555pmNGzdWVFR89atfHaJjAQAwYEO4DdOjjz6akJDw3e9+NycnZ+3atUlJ/d1EBgCA82YIC+r7j4J6AEC48KxRAIDVSIQAAKuRCAEAViMRAgCsRiIEAFiNRAgAsBqJEABgNRIhAMBqJEIAgNVIhAAAq5EIAQBWIxECAKxGIgQAWI1ECACwGokQAGA1EiEAwGokQgCA1UiEAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAaiRCAIDVSIQAAKuRCAEAViMRAgCsRiIEAFiNRAgAsBqJEABgNRIhAMBqJEIAgNVIhAAAq5EIAQBWIxECAKxGIgQAWI1ECACwGokQAGA1EiEAwGokQgCA1UiEAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAaiRCAIDVSIQAAKuRCAEAViMRAgCsRiIEAFiNRAgAsBqJEABgNRIhAMBqJEIAgNVIhAAAq5EIAQBWIxECAKxGIgQAWI1ECACwWuTgm2hra1uzZs2uXbt6enoee+yxvrjX6/3xj3+8evXq/Pz8b3/728XFxYM/FgAAoRWCT4Q7d+78/ve/v3///qeffvrM+LJly376059+61vfysjIuOKKKzo6OgZ/LAAAQisEiXDu3Lnr1q175JFHzgz29PQ899xzP/rRj6699trHH3+8oKDgf/7nfwZ/LAAAQmuo/kZYXV3d2Ng4e/bswLdz5swpLy8fomMBADBg/fobYVdXV21trTNeWFgYFRUlf6ShoSEpKSky8g/tZ2RkbN++PVj7x44du//++1NSUgLfjhgx4uc//3l/OgYAwCD1KxHu3r37zjvvdMZXrlxZUlIifyQhIaGrq6vv287OzsTExGDtZ2Zmfv3rXx87dmzg29zc3P70CgCAwetXIpw+fXpFRcU5tVtQUODxeOrq6vLy8owxVVVVRUVFwV4cGxtbVlY2ffr0czoEAACDN1R/I8zIyJg/f/5LL71kjKmvr1+5cuUtt9wyRMcCAGDAQlBHWFdXN378eK/X29nZmZ6eftFFF33wwQfGmKeeeur6669ftWpVZWXlkiVLJk+ePPhjAQAQWiFIhDk5OZWVlX3fulx/+JQ5efLkysrKvXv35uTkfMb/iwIAEEYhSIQulystLU3+U+CPf4M/BAAAQ4RnjQIArEYiBABYjUQIALAaiRAAYDUSIQDAaiTC0Pu3f/s3r9cb7l6E3qZNm9auXRvuXoSe3+9//PHHw92LIbF69ept27aFuxeh5/F4nnzyyXD3Yki88cYbH374Ybh7EXotLS0//OEPw92LoEiEobds2bL29vZw9yL0/u///u/3v/99uHsRel6v99FHHw13L4bEO++8s3HjxnD3IvSam5s/tfvp58Zbb731udyo59ixYz/5yU/C3YugSIQAAKuRCAEAViMRAgCsFuH3+8PdB5Ofn+92u6Ojo8PdkdA4cuRIUVFRREREuDsSYqdPn/b5fMEepzesVVVVjRw5Mty9CL2TJ09GRkYmJyeHuyMh5vP5jh079rl8gvGJEydiY2M/Y/fWYaqnp6exsbGgoOD8H/prX/vaY4899tmvuSASYUNDw+dpdUl3d3dMTEy4exF6gaWwbrc73B0Jvc/rJevt7Y2IiOCSDSM9PT1ut7tv64LPk3Bdsry8vLi4uM9+zQWRCAEACJfP4e8dAAD0H4kQAGA1EiEAwGokQgCA1UKwQ73lamtrd+zYUVdXd80115y5nruysvK//uu/PB7PbbfdNnny5DD2cGCampp+97vf7du3Lzk5efHixePHj+/7p3feeeedd97Jzc2955570tPTw9jJgXnrrbd27NjR2tpaWlp6xx13JCUlBeL19fUvvfRSU1PTl770pblz54a3k4Oxa9eubdu23XrrrampqYHI+vXr33zzzYyMjLvuuis3Nze83RuAV155paurK/D1qFGjFi5cGPj6yJEjL7/8cnt7+y233FJWVha+Dg6c1+v91a9+tX379rS0tMWLF0+YMMEY4/f7f/3rX2/dunXkyJHf/OY3z7ro8UJz+PDhd99998zIDTfcEBh4Bw4cePXVV40xt99++7hx48LTPwc+EQ7WtGnTnnzyyaVLl+7Zs6cveOTIkbKyMo/Hk5qaesUVVwzHhwc+/PDDb731Vk5OzunTp6dNm9Y3rH/+85/fddddo0aN2rNnz5w5czweT3j7OQCvvvpqVFTUqFGjVqxYMWPGjMAM29raOn369KqqqqKioptvvnnFihXh7uYAtba23nnnnffee29DQ0MgsmLFisWLFxcWFh45cmT69OktLS3h7eEAPPzww9u3bz906NChQ4caGxsDwfr6+rKyslOnTmVnZy9cuHD9+vXh7eQAeL3e66+//rnnnisoKPD5fFu2bAnEH3300e9973ulpaWrVq1atGhReDs5AG1tbYf+aO3atX/zN38TKAg5cODAjBkzoqKiYmJiZs6cuX///nD39I/8GByv1+v3+4uLi3/729/2BR955JE77rgj8PV3v/vdW2+9NTydG4TOzs6+r5cuXXrLLbf4/X6fz1daWvrGG28Evp44ceKvfvWrsHVx0Hp6ehITE7ds2eL3+1944YU5c+YE4i+99FJZWVlYuzZw3/rWt55//nljzEcffRSITJs27T//8z8DX8+dO/f5558PX+8GKCsrq+/t9Pne97534403Br5+8sknr7vuuvPer8H62c9+Nn78eI/Hc2awra0tJSVl586dfr+/q6srIyPj/fffD1MHQ+Db3/724sWLA1/fe++99913X+DrBx988J577glfv/4MnwgHS5a+rl+//qqrrgp8vXDhwnXr1p3fToVAbGxs39fd3d2BR13U1dUdPHgw8B9TERERCxYsGI5vrc/OnTtdLtfo0aON45KVl5d3dHSEtXcDsW7dut27d3/zm9/si3R1dW3btq3v/xKH6Wg0xrz22mvPPvvsmZtpfA7uspUrV95xxx1vv/32M888s3nz5kBw165dUVFRl112mTEmJiZm7ty5w/GtBfT29r766qt33XVX4Nv169f3DcWrrrrqwnlfJMIhUVdXl5WVFfg6Ozv7+PHjvb294e3SgO3evfuVV17527/9W2NMXV1dfHx83/OfcnJyamtrw9q7AVqyZEl+fv68efP++7//Ozs72zguWSASzi6eu46Ojvvuu+9nP/vZmb+cBd5F4B2ZYXvJZs6c2d3dXVVVtXjx4oceeigQ/NQla29vH3b/63v48OEXX3xx+fLlbW1tt9xyy7PPPmuMqa+v77teZthesoDf/e53LpfrmmuuCXz7qUt24dxiLJYZEpGRkX2Zr7e3d/g+M+nIkSM33njj008/PXHiRGNMVFTUmRm9p6dnmD7m6kc/+lFLS8vKlSuXLFmybdu24uLiT10yY8ywe/jtP/7jP952223jxo07c1/oqKgo88d3ZIbtJXvzzTcDX9x3332XXHLJAw88MHr0aOclC7zZYcTlcpWWlr7yyivGmOnTp996660PPvjgme/LGNPT0zN8nxb70ksv/eVf/mXfQ/4+dckunFtsWM7OF76CgoK+X+Jqamry8vKGYyI8cuTIvHnz/v7v//7uu+8ORPLz8z0ez/HjxwPfBt5a+Do4cImJifn5+ffcc8+ll166cuVK8+eX7NixY263OycnJ6x9PGevvfba8uXLp06dOm3aNGPMzTffHFju5Ha7a2pqAq8ZvpcsoLS0NCsr6/Dhw8Zxl6Wnpw+71ZUFBQV9KyfHjx9/+vTp06dP5+fn19fX9/02M3wvWX19/TvvvLNkyZK+yKcuWX5+fpi69mnDb3YeFhYtWrR8+XK/32+MWb58+XBc93Xs2LH58+cvXbr03nvv7QtmZmbOnDlz+fLlxpiOjo6VK1fecMMN4evjQHR3d/t8vsDXLS0tFRUVgaKXRYsWrVixIrAI9je/+c3VV1994fy62k+rV69+6aWXfvKTn/zHf/yHMeZf/uVfZs2aFRUVdc011wQuWU9Pz5tvvjnsLllXV5f/j49E3rZt24kTJ8aOHWuMWbRo0euvvx64msP0Lrvpppu2bt0aeHfvv/9+fn5+amrqZZddlpycvHr1amNMQ0PDhg0bvvjFL4a7pwPx8ssvz5w5s6SkpC8SmBgDX//617++gC5ZmBfrDH9LliyZMmVKTExMSUnJlClT9u/f7/f7m5ubx48fP3/+/JtvvrmgoKCqqirc3TxnX/7yl2NjY6f80Z133hmIr1mzJiMj4+tf//rkyZNvuOEGn88X3n6eq61btwYKJG677bbc3Nxbb701sO7X4/HMmzdv2rRpt99+e2ZmZnl5ebh7OnCB/33qW2a5ffv2jIyM22+/ffr06fPmzfvUGsUL36pVq4qLi7/yla986UtfSkpK+vd///dAvK2t7fLLL587d+6tt96ak5Pz8ccfh7efA9DZ2Tl79uy5c+fefffdWVlZgfXYfr//1VdfzcrKWrJkSUlJyf333x/eTg7Y2LFjX3nllTMj9fX1xcXF11133Re/+MXRo0fX1dWFq2+fwu4Tg/Xxxx+3tbX1fTt27NiEhARjTFdX13vvvdfd3b1gwYKUlJTwdXCAKisrT5061fdtQkJC4NdwY0xNTc2GDRtycnKuuOKK4fhfvgcPHty3b5/X6x03btyZJb29vb1r1649efLklVdeeeZqhWHH7/fv3Llz/PjxfUt/Gxsb165dm56efuWVV0ZGDrOVAV6vd8+ePRUVFYHfzEaMGNH3T93d3WvWrGlra1uwYMEw3SkzMOpaWlpmzJhx5nZ9Bw8e3L59+6hRo2bMmBHG7g1Yb2/v7t27J0yY8Km/Sbe3t7/33nt+v3/BggUXzraLJEIAgNWG36/zAACEEIkQAGA1EiEAwGokQgCA1UiEAACrkQgBAFYjEQIArEYiBABYjUQIALAaiRAAYDUSIQDAav8fbtMQfAQrAGMAAAAASUVORK5CYII=", "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "execution_count": 10 } ], "cell_type": "code", "source": [ "λ_Kerker, X_Kerker = eigsolve(Pinv_Kerker ∘ epsilon,\n", " randn(size(scfres_Al.ρ)), 3, :LM;\n", " tol=1e-3, eager=true, verbosity=2)\n", "\n", "mode_xy = mean(real.(X_Kerker[1]), dims=3)[:, :, 1, 1] # Average along z axis\n", "heatmap(mode_xy', c=:RdBu_11, aspect_ratio=1, grid=false,\n", " legend=false, clim=(-0.006, 0.006))" ], "metadata": {}, "execution_count": 10 }, { "cell_type": "markdown", "source": [ "Clearly the charge-sloshing mode is no longer dominating.\n", "\n", "The largest eigenvalue is now" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "4.7236110498662836" }, "metadata": {}, "execution_count": 11 } ], "cell_type": "code", "source": [ "maximum(real.(λ_Kerker))" ], "metadata": {}, "execution_count": 11 }, { "cell_type": "markdown", "source": [ "Since the smallest eigenvalue in this case remains of similar size (it is now\n", "around 0.8), this implies that the conditioning improves noticeably when\n", "`KerkerMixing` is used.\n", "\n", "**Note:** Since LdosMixing requires solving a linear system at each\n", "application of $P^{-1}$, determining the eigenvalues of\n", "$P^{-1} \\varepsilon^\\dagger$ is slightly more expensive and thus not shown. The\n", "results are similar to `KerkerMixing`, however." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "We could repeat the exercise for an insulating system (e.g. a Helium chain).\n", "In this case you would notice that the condition number without mixing\n", "is actually smaller than the condition number with Kerker mixing. In other\n", "words employing Kerker mixing makes the convergence *worse*. A closer\n", "investigation of the eigenvalues shows that Kerker mixing reduces the\n", "smallest eigenvalue of the dielectric operator this time, while keeping\n", "the largest value unchanged. Overall the conditioning thus workens." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "**Takeaways:**\n", "- For metals the conditioning of the dielectric matrix increases steeply with system size.\n", "- The Kerker preconditioner tames this and makes SCFs on large metallic\n", " systems feasible by keeping the condition number of order 1.\n", "- For insulating systems the best approach is to not use any mixing.\n", "- **The ideal mixing** strongly depends on the dielectric properties of\n", " system which is studied (metal versus insulator versus semiconductor)." ], "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 }