The Dirac equation for the Hydrogen atom#
+In this notebook we will illustrate how one can solve the Dirac equation for a Hydrogen atom using the Multiwavelet framework provided by VAMPyR
+The Dirac equation can be coincisely written as follows:
+where \(\phi = (\phi^{L\alpha}, \phi^{L\beta}, \phi^{S\alpha}, \phi^{S\beta})^t\) represents a 4-component spinor, $\boldsymbol{\alpha} =
+\( and +\)\beta =
+\( 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.
+As for the non-relativistic case, the equation is first rewritten in its integral formulation: +$\(\phi = \frac{1}{2mc^2}(\beta m c^2+ c \boldsymbol{\alpha} \cdot \mathbf{p} + \epsilon) \left[ G_\mu \star (V \psi) \right]\)$
+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:
+The guess obtained is then plugged into the integral form of the Dirac equation, which can then be iterated until convergence
+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
from vampyr import vampyr3d as vp
+from orbital4c import orbital as orb
+from orbital4c import nuclear_potential as nucpot
+from orbital4c import complex_fcn as cf
+import numpy as np
+
As a reference, the exact Dirac energy for the ground state Hydrogen atom is computed in the function below.
+def analytic_1s(light_speed, n, k, Z):
+ alpha = 1/light_speed
+ gamma = orb.compute_gamma(k,Z,alpha)
+ tmp1 = n - np.abs(k) + gamma
+ tmp2 = Z * alpha / tmp1
+ tmp3 = 1 + tmp2**2
+ return light_speed**2 / np.sqrt(tmp3)
+
+light_speed = 137.03599913900001
+alpha = 1/light_speed
+k = -1
+l = 0
+n = 1
+m = 0.5
+Z = 1
+atom = "H"
+
+energy_1s = analytic_1s(light_speed, n, k, Z)
+print('Exact Energy',energy_1s - light_speed**2, flush = True)
+
Exact Energy -0.5000066565989982
+
The MultiResolutionAnalysis
object defining the simulation box is constructed and passed to the classes for complex functions and spinors
mra = vp.MultiResolutionAnalysis(box=[-30,30], order=6)
+prec = 1.0e-4
+origin = [0.1, 0.2, 0.3] # origin moved to avoid placing the nuclar charge on a node
+
+orb.orbital4c.light_speed = light_speed
+orb.orbital4c.mra = mra
+cf.complex_fcn.mra = mra
+
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
a_coeff = 3.0
+b_coeff = np.sqrt(a_coeff/np.pi)**3
+gauss = vp.GaussFunc(b_coeff, a_coeff, origin)
+gauss_tree = vp.FunctionTree(mra)
+vp.advanced.build_grid(out=gauss_tree, inp=gauss)
+vp.advanced.project(prec=prec, out=gauss_tree, inp=gauss)
+gauss_tree.normalize()
+
+spinor_H = orb.orbital4c()
+La_comp = cf.complex_fcn()
+La_comp.copy_fcns(real = gauss_tree)
+spinor_H.copy_components(La = La_comp)
+spinor_H.init_small_components(prec/10)
+spinor_H.normalize()
+
The nuclear potential is defined and projected onto the V_tree
Peps = vp.ScalingProjector(mra, prec)
+f = lambda x: nucpot.coulomb_HFYGB(x, origin, Z, prec)
+V_tree = Peps(f)
+
+default_der = 'BS'
+
The orbital is optimized by iterating the integral version of the Dirac equation as follows:
+-
+
Application of the Dirac Hamiltonian \(f^n = \hat{h}_D \phi^n = (\beta m c^2+ c \boldsymbol{\alpha} \cdot \mathbf{p}) \phi^n\)
+Application of the potnetial operator \(g^n = \hat{V} \phi^n\)
+Sum \(h^n = f^n + g^n\)
+Expectation value of the energy \(\left\langle \phi^n | h^n \right\rangle\)
+Calculation of the Helmholtz parameter \(\mu\)
+Convolution with the Helmholtz kernel \(t^n = G_\mu \star g^n\)
+application of the shifted Dirac Hamiltonian \(\tilde{\phi}^{n+1} = (\hat{h}_D + \epsilon) t^n\)
+normalization of the new iterate
+calculation of the change in the orbital \(\delta \phi^n = \phi^{n+1}-\phi^n\)
+
Once the orbital error is below the requested threshold the iteration is interrupted and the final energy is computed.
+orbital_error = 1
+while orbital_error > prec:
+ # 1. Application of the Dirac Hamiltonian
+ hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)
+ # 2. Application of the potnetial operator
+ v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
+ # 3. Sum
+ add_psi = hd_psi + v_psi
+ # 4. Expectation value of the energy
+ energy = (spinor_H.dot(add_psi)).real
+ print('Energy',energy-light_speed**2)
+ # 5. Calculation of the Helmholtz parameter
+ mu = orb.calc_dirac_mu(energy, light_speed)
+ # 6. Convolution with the Helmholtz kernel
+ tmp = orb.apply_helmholtz(v_psi, mu, prec)
+ tmp.crop(prec/10)
+ # 7. application of the shifted Dirac Hamiltonian
+ new_orbital = orb.apply_dirac_hamiltonian(tmp, prec, energy, der = default_der)
+ new_orbital.crop(prec/10)
+ # 8. normalization of the new iterate
+ new_orbital.normalize()
+ delta_psi = new_orbital - spinor_H
+ # 9. calculation of the change in the orbital
+ orbital_error = (delta_psi.dot(delta_psi)).real
+ print('Error',orbital_error, flush = True)
+ spinor_H = new_orbital
+
+# Computing the final energy
+hd_psi = orb.apply_dirac_hamiltonian(spinor_H, prec, der = default_der)
+v_psi = orb.apply_potential(-1.0, V_tree, spinor_H, prec)
+add_psi = hd_psi + v_psi
+energy = (spinor_H.dot(add_psi)).real
+print('Final Energy',energy - light_speed**2)
+
+energy_1s = analytic_1s(light_speed, n, k, Z)
+
+print('Exact Energy',energy_1s - light_speed**2)
+print('Difference',energy_1s - energy)
+
Energy -0.14180720307558659
+Error 0.3648353655184581
+Energy -0.49611934473068686
+Error 0.0024839139737050653
+Energy -0.4997319277244969
+Error 0.00023304055836475225
+Energy -0.49998214000152075
+Error 2.3987883028320866e-05
+Final Energy -0.5000042447172746
+Exact Energy -0.5000066565989982
+Difference -2.411881723674014e-06
+