|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "markdown", |
| 5 | + "id": "90e862bb-783a-4f97-a227-8f6e0b4012b8", |
| 6 | + "metadata": {}, |
| 7 | + "source": [ |
| 8 | + "# The Dirac equation for the Hydrogen atom\n", |
| 9 | + "\n", |
| 10 | + "In this notebook we will illustrate how one can solve the Dirac equation for a Hydrogen atom using the Multiwavelet framework provided by VAMPyR\n", |
| 11 | + "\n", |
| 12 | + "The Dirac equation can be coincisely written as follows:\n", |
| 13 | + "\n", |
| 14 | + "\\begin{equation}\n", |
| 15 | + "(\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p} + V) \\phi = \\epsilon \\phi \n", |
| 16 | + "\\end{equation}\n", |
| 17 | + "\n", |
| 18 | + "where $\\phi = (\\phi^{L\\alpha}, \\phi^{L\\beta}, \\phi^{S\\alpha}, \\phi^{S\\beta})^t$ represents a 4-component spinor, $\\boldsymbol{\\alpha} = \n", |
| 19 | + "\\begin{pmatrix}\n", |
| 20 | + "0_2 & \\boldsymbol{\\sigma} \\\\\n", |
| 21 | + "\\boldsymbol{\\sigma} & 0_2 & \n", |
| 22 | + "\\end{pmatrix}\n", |
| 23 | + "$ and\n", |
| 24 | + "$\\beta = \n", |
| 25 | + "\\begin{pmatrix}\n", |
| 26 | + "1_2 & 0_2 \\\\\n", |
| 27 | + "0_2 & -1_2\n", |
| 28 | + "\\end{pmatrix}\n", |
| 29 | + "$ are the 4x4 Dirac matrices, $\\boldsymbol{\\sigma}$ is a cartesian vector collecting the three 2x2 Pauli matrices, $\\mathbf{p}$ is the momentum operator, $c$ is the speed of light, $m$ is the electron mass and $V$ is the nuclear potential.\n", |
| 30 | + "\n", |
| 31 | + "As for the non-relativistic case, the equation is first rewritten in its integral formulation:\n", |
| 32 | + "$$\\phi = \\frac{1}{2mc^2}(\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p} + \\epsilon) \\left[ G_\\mu \\star (V \\psi) \\right]$$\n", |
| 33 | + "\n", |
| 34 | + "where $G_\\mu(x) = \\frac{e^{-\\mu |x|}}{4 \\pi |x|}$ is the Helmholtz Green's kernel and $\\mu = \\sqrt{\\frac{m^2c^4-\\epsilon}{mc^2}}$. An initial guess can be obtained by taking a Slater orbital or a Gaussian function for the $\\psi^{L\\alpha}$ component and then applying the restricted kinetic balance:\n", |
| 35 | + "\n", |
| 36 | + "$$\n", |
| 37 | + "\\begin{pmatrix}\n", |
| 38 | + "\\phi^{S\\alpha} \\\\\n", |
| 39 | + "\\phi^{S\\beta}\n", |
| 40 | + "\\end{pmatrix}\n", |
| 41 | + "= \\frac{1}{2c}\\boldsymbol{\\sigma} \\cdot \\mathbf{p} \n", |
| 42 | + "\\begin{pmatrix}\n", |
| 43 | + "\\phi^{L\\alpha} \\\\\n", |
| 44 | + "0\n", |
| 45 | + "\\end{pmatrix}\n", |
| 46 | + "$$\n", |
| 47 | + "The guess obtained is then plugged into the integral form of the Dirac equation, which can then be iterated until convergence" |
| 48 | + ] |
| 49 | + }, |
| 50 | + { |
| 51 | + "cell_type": "markdown", |
| 52 | + "id": "54c75006-0358-42c5-9075-37a5b8a909ba", |
| 53 | + "metadata": {}, |
| 54 | + "source": [ |
| 55 | + "We start by loading the relevant packages: the 3d version of `vampyr`, `numpy`, the `complex_function` class which deals with complex funtions and the `orbital` class which deals with 4-component spinors. Each complex function is handled as a pair of `function_tree`s and each spinor is handled as a 4c vector of complex functions. The `nuclear_potential` package is self-explanatory" |
| 56 | + ] |
| 57 | + }, |
| 58 | + { |
| 59 | + "cell_type": "code", |
| 60 | + "execution_count": 13, |
| 61 | + "id": "37e0e6b2-886e-4415-8612-d2a652f24c4f", |
| 62 | + "metadata": {}, |
| 63 | + "outputs": [], |
| 64 | + "source": [ |
| 65 | + "from vampyr import vampyr3d as vp\n", |
| 66 | + "from orbital4c import orbital as orb\n", |
| 67 | + "from orbital4c import nuclear_potential as nucpot\n", |
| 68 | + "from orbital4c import complex_fcn as cf\n", |
| 69 | + "import numpy as np" |
| 70 | + ] |
| 71 | + }, |
| 72 | + { |
| 73 | + "cell_type": "markdown", |
| 74 | + "id": "7fbe4360-9936-47d4-bb6f-d51bb3b34b69", |
| 75 | + "metadata": {}, |
| 76 | + "source": [ |
| 77 | + "As a reference, the exact Dirac energy for the ground state Hydrogen atom is computed in the function below." |
| 78 | + ] |
| 79 | + }, |
| 80 | + { |
| 81 | + "cell_type": "code", |
| 82 | + "execution_count": 14, |
| 83 | + "id": "0ef3d827-810b-4411-8ba3-bd63271033e4", |
| 84 | + "metadata": {}, |
| 85 | + "outputs": [ |
| 86 | + { |
| 87 | + "name": "stdout", |
| 88 | + "output_type": "stream", |
| 89 | + "text": [ |
| 90 | + "Exact Energy -0.5000066565989982\n" |
| 91 | + ] |
| 92 | + } |
| 93 | + ], |
| 94 | + "source": [ |
| 95 | + "def analytic_1s(light_speed, n, k, Z):\n", |
| 96 | + " alpha = 1/light_speed\n", |
| 97 | + " gamma = orb.compute_gamma(k,Z,alpha)\n", |
| 98 | + " tmp1 = n - np.abs(k) + gamma\n", |
| 99 | + " tmp2 = Z * alpha / tmp1\n", |
| 100 | + " tmp3 = 1 + tmp2**2\n", |
| 101 | + " return light_speed**2 / np.sqrt(tmp3)\n", |
| 102 | + "\n", |
| 103 | + "light_speed = 137.03599913900001\n", |
| 104 | + "alpha = 1/light_speed\n", |
| 105 | + "k = -1\n", |
| 106 | + "l = 0\n", |
| 107 | + "n = 1\n", |
| 108 | + "m = 0.5\n", |
| 109 | + "Z = 1\n", |
| 110 | + "atom = \"H\"\n", |
| 111 | + "\n", |
| 112 | + "energy_1s = analytic_1s(light_speed, n, k, Z)\n", |
| 113 | + "print('Exact Energy',energy_1s - light_speed**2, flush = True)" |
| 114 | + ] |
| 115 | + }, |
| 116 | + { |
| 117 | + "cell_type": "markdown", |
| 118 | + "id": "6277d776-2236-496a-bdd5-41929d11ac3d", |
| 119 | + "metadata": {}, |
| 120 | + "source": [ |
| 121 | + "The `MultiResolutionAnalysis` object defining the simulation box is constructed and passed to the classes for complex functions and spinors" |
| 122 | + ] |
| 123 | + }, |
| 124 | + { |
| 125 | + "cell_type": "code", |
| 126 | + "execution_count": 15, |
| 127 | + "id": "5ac41574-6f43-411d-841f-f8fdf0749b66", |
| 128 | + "metadata": {}, |
| 129 | + "outputs": [], |
| 130 | + "source": [ |
| 131 | + "mra = vp.MultiResolutionAnalysis(box=[-30,30], order=6)\n", |
| 132 | + "prec = 1.0e-4\n", |
| 133 | + "origin = [0.1, 0.2, 0.3] # origin moved to avoid placing the nuclar charge on a node\n", |
| 134 | + "\n", |
| 135 | + "orb.orbital4c.light_speed = light_speed\n", |
| 136 | + "orb.orbital4c.mra = mra\n", |
| 137 | + "cf.complex_fcn.mra = mra" |
| 138 | + ] |
| 139 | + }, |
| 140 | + { |
| 141 | + "cell_type": "markdown", |
| 142 | + "id": "d5abdbe4-ba5b-49d3-ae7e-febab5ced4fe", |
| 143 | + "metadata": {}, |
| 144 | + "source": [ |
| 145 | + "We construct a starting guess by taking a simple Gaussian function and initialize the real part of the $\\phi^{L\\alpha}$ component of the spinor with it. Thereafter the restricted kinetic balance is employed. This is implemented in the `init_small_components` method of the `orbital` class" |
| 146 | + ] |
| 147 | + }, |
| 148 | + { |
| 149 | + "cell_type": "code", |
| 150 | + "execution_count": 16, |
| 151 | + "id": "95826623-b759-4318-afc5-e0c225cb5f1d", |
| 152 | + "metadata": {}, |
| 153 | + "outputs": [], |
| 154 | + "source": [ |
| 155 | + "a_coeff = 3.0\n", |
| 156 | + "b_coeff = np.sqrt(a_coeff/np.pi)**3\n", |
| 157 | + "gauss = vp.GaussFunc(b_coeff, a_coeff, origin)\n", |
| 158 | + "gauss_tree = vp.FunctionTree(mra)\n", |
| 159 | + "vp.advanced.build_grid(out=gauss_tree, inp=gauss)\n", |
| 160 | + "vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss)\n", |
| 161 | + "gauss_tree.normalize()\n", |
| 162 | + "\n", |
| 163 | + "spinor_H = orb.orbital4c()\n", |
| 164 | + "La_comp = cf.complex_fcn()\n", |
| 165 | + "La_comp.copy_fcns(real = gauss_tree)\n", |
| 166 | + "spinor_H.copy_components(La = La_comp)\n", |
| 167 | + "spinor_H.init_small_components(prec/10)\n", |
| 168 | + "spinor_H.normalize()" |
| 169 | + ] |
| 170 | + }, |
| 171 | + { |
| 172 | + "cell_type": "markdown", |
| 173 | + "id": "f6a40ddd-7933-46ee-8128-96caf89f3682", |
| 174 | + "metadata": {}, |
| 175 | + "source": [ |
| 176 | + "The nuclear potential is defined and projected onto the `V_tree`" |
| 177 | + ] |
| 178 | + }, |
| 179 | + { |
| 180 | + "cell_type": "code", |
| 181 | + "execution_count": 17, |
| 182 | + "id": "4f0cda6f-217c-4941-8ce5-fa52576c2d98", |
| 183 | + "metadata": {}, |
| 184 | + "outputs": [], |
| 185 | + "source": [ |
| 186 | + "Peps = vp.ScalingProjector(mra, prec)\n", |
| 187 | + "f = lambda x: nucpot.coulomb_HFYGB(x, origin, Z, prec)\n", |
| 188 | + "V_tree = Peps(f)\n", |
| 189 | + "\n", |
| 190 | + "default_der = 'BS'" |
| 191 | + ] |
| 192 | + }, |
| 193 | + { |
| 194 | + "cell_type": "markdown", |
| 195 | + "id": "1678bc97-d85d-4a48-9aea-13901631534d", |
| 196 | + "metadata": {}, |
| 197 | + "source": [ |
| 198 | + "The orbital is optimized by iterating the integral version of the Dirac equation as follows:\n", |
| 199 | + "1. Application of the Dirac Hamiltonian $f^n = \\hat{h}_D \\phi^n = (\\beta m c^2+ c \\boldsymbol{\\alpha} \\cdot \\mathbf{p}) \\phi^n$\n", |
| 200 | + "2. Application of the potnetial operator $g^n = \\hat{V} \\phi^n$\n", |
| 201 | + "3. Sum $h^n = f^n + g^n$\n", |
| 202 | + "4. Expectation value of the energy $\\left\\langle \\phi^n | h^n \\right\\rangle$\n", |
| 203 | + "5. Calculation of the Helmholtz parameter $\\mu$\n", |
| 204 | + "6. Convolution with the Helmholtz kernel $t^n = G_\\mu \\star g^n$\n", |
| 205 | + "7. application of the shifted Dirac Hamiltonian $\\tilde{\\phi}^{n+1} = (\\hat{h}_D + \\epsilon) t^n$\n", |
| 206 | + "8. normalization of the new iterate\n", |
| 207 | + "9. calculation of the change in the orbital $\\delta \\phi^n = \\phi^{n+1}-\\phi^n$\n", |
| 208 | + "\n", |
| 209 | + "Once the orbital error is below the requested threshold the iteration is interrupted and the final energy is computed." |
| 210 | + ] |
| 211 | + }, |
| 212 | + { |
| 213 | + "cell_type": "code", |
| 214 | + "execution_count": 18, |
| 215 | + "id": "40a443b6-33f4-4771-8b06-ccd12c9da35d", |
| 216 | + "metadata": {}, |
| 217 | + "outputs": [ |
| 218 | + { |
| 219 | + "name": "stdout", |
| 220 | + "output_type": "stream", |
| 221 | + "text": [ |
| 222 | + "Energy -0.14180720307558659\n", |
| 223 | + "Error 0.3648353655184581\n", |
| 224 | + "Energy -0.49611934473068686\n", |
| 225 | + "Error 0.0024839139737050653\n", |
| 226 | + "Energy -0.4997319277244969\n", |
| 227 | + "Error 0.00023304055836475225\n", |
| 228 | + "Energy -0.49998214000152075\n", |
| 229 | + "Error 2.3987883028320866e-05\n", |
| 230 | + "Final Energy -0.5000042447172746\n", |
| 231 | + "Exact Energy -0.5000066565989982\n", |
| 232 | + "Difference -2.411881723674014e-06\n" |
| 233 | + ] |
| 234 | + } |
| 235 | + ], |
| 236 | + "source": [ |
| 237 | + "orbital_error = 1\n", |
| 238 | + "while orbital_error > prec:\n", |
| 239 | + " # 1. Application of the Dirac Hamiltonian\n", |
| 240 | + " hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)\n", |
| 241 | + " # 2. Application of the potnetial operator\n", |
| 242 | + " v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)\n", |
| 243 | + " # 3. Sum\n", |
| 244 | + " add_psi = hd_psi + v_psi\n", |
| 245 | + " # 4. Expectation value of the energy\n", |
| 246 | + " energy = (spinor_H.dot(add_psi)).real\n", |
| 247 | + " print('Energy',energy-light_speed**2)\n", |
| 248 | + " # 5. Calculation of the Helmholtz parameter\n", |
| 249 | + " mu = orb.calc_dirac_mu(energy, light_speed)\n", |
| 250 | + " # 6. Convolution with the Helmholtz kernel\n", |
| 251 | + " tmp = orb.apply_helmholtz(v_psi, mu, prec)\n", |
| 252 | + " tmp.crop(prec/10)\n", |
| 253 | + " # 7. application of the shifted Dirac Hamiltonian\n", |
| 254 | + " new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = default_der) \n", |
| 255 | + " new_orbital.crop(prec/10)\n", |
| 256 | + " # 8. normalization of the new iterate\n", |
| 257 | + " new_orbital.normalize()\n", |
| 258 | + " delta_psi = new_orbital - spinor_H\n", |
| 259 | + " # 9. calculation of the change in the orbital\n", |
| 260 | + " orbital_error = (delta_psi.dot(delta_psi)).real\n", |
| 261 | + " print('Error',orbital_error, flush = True)\n", |
| 262 | + " spinor_H = new_orbital\n", |
| 263 | + "\n", |
| 264 | + "# Computing the final energy\n", |
| 265 | + "hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)\n", |
| 266 | + "v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)\n", |
| 267 | + "add_psi = hd_psi + v_psi\n", |
| 268 | + "energy = (spinor_H.dot(add_psi)).real\n", |
| 269 | + "print('Final Energy',energy - light_speed**2)\n", |
| 270 | + "\n", |
| 271 | + "energy_1s = analytic_1s(light_speed, n, k, Z)\n", |
| 272 | + "\n", |
| 273 | + "print('Exact Energy',energy_1s - light_speed**2)\n", |
| 274 | + "print('Difference',energy_1s - energy)" |
| 275 | + ] |
| 276 | + } |
| 277 | + ], |
| 278 | + "metadata": { |
| 279 | + "kernelspec": { |
| 280 | + "display_name": "Python 3 (ipykernel)", |
| 281 | + "language": "python", |
| 282 | + "name": "python3" |
| 283 | + }, |
| 284 | + "language_info": { |
| 285 | + "codemirror_mode": { |
| 286 | + "name": "ipython", |
| 287 | + "version": 3 |
| 288 | + }, |
| 289 | + "file_extension": ".py", |
| 290 | + "mimetype": "text/x-python", |
| 291 | + "name": "python", |
| 292 | + "nbconvert_exporter": "python", |
| 293 | + "pygments_lexer": "ipython3", |
| 294 | + "version": "3.10.13" |
| 295 | + } |
| 296 | + }, |
| 297 | + "nbformat": 4, |
| 298 | + "nbformat_minor": 5 |
| 299 | +} |
0 commit comments