Scalar Dirac equation, and numerical considerations

Table of Contents

• Preliminaries
  • Scalar Dirac equations
  • Boundary Conditions near the origin
  • Boundary Conditions at the augmentation radius
• Solving the Dirac equation Radial Mesh
  • Legacy and modern implementations
  • Input file switches controlling how the integrations are performed: the RSEDEF token

Preliminaries

Questaal codes normally solve the Dirac equation in a scalar relativistic approximation for partial waves inside the augmentation spheres, though some parts of the code (mostly ASA codes) have extensions to the full Dirac equation.

These integrations are performed numerically on a discrete radial mesh. For the most part you don’t need to worry about how integrations are carried out; for the most part Questaal takes care of the numerics pretty well without you needing to delve into the details. If you are haveing an issue and looking for a quick fix, see documentation for the RSEDEF token.

There are details that at times can be important, (occasionally the integrators can trip up or be a little inaccurate), which this page documents.

Scalar Dirac equations

The scalar relativistic approximation reduces the four component Dirac equation to two components (or two pairs of components in a spin polarized case). The spin orbit couping that connects the two pairs can be included perturbatively and has a structure of the form $L{\cdot}S$. The two components consist of a large component $g$ and a small component $f$. In a spherical potential the two satisfy these coupled equations $$\frac{dg}{dr} = \frac{g}{r} + u f \ \ \text{ and }\ \ \frac{df}{dr} = -\frac{f}{r} - y g$$

where $$u(r) = c + \phi(r), \ \ y(r) = -\frac{l(l+1)}{r^2 u(r)} + \phi(r), \ \ \text{ and }\ \ \phi(r) = \frac{1}{c} \left( \epsilon + \frac{2Z}{r} - v(r) \right)$$

$c$ is the speed of light, $\epsilon$ is the energy, $2Z/r$ is the nuclear potential (in Ry, the units Questaal uses) and $v(r)$ the external potential. Note that the partial waves have an $l$ quantum number, but the index is suppressed.

Boundary Conditions near the origin

A second order differential equation (or two coupled first order ones) has two solutions, one regular at the origin and one irregular. On physical grounds the solution must be regular at the origin, and for the scalar Dirac case the behavior at the origin is defined by $$g = C r^{\gamma}, \quad \text{where} \quad \gamma = \left[1 + l(l+1) - \left(2Z/c\right)^2\right]^{1/2}$$

Boundary Conditions at the augmentation radius

We are left with one free boundary condition. For the free atom, there is a physical requirement that the wave function is normalizable, which requires the solution decay as $r{\rightarrow}\infty$. Solutions that are regular both at the origin and at $\infty$ occur only at discrete energies $\epsilon$; this is how the quantization condition manifests itself.

In the solid with an all-electron method, the partial waves are evaluated only up to a “muffin tin” radius (aka augmentation radius) $r_m$, and joined to envelope functions there, which for Questaal are smoothed Hankel functions, jigsaw puzzle orbitals, or possibly plane waves. The partial wave must join smoothly and differentiably to the envelope function: matching value can always accomplished by normalization of the partial wave, but to match the slope (or logarithmic derivative, $D=r_m{\times}\mathrm{slope/value}$) is valid only for a particular $\epsilon$. Thus we have freedom to choose either $\epsilon$ or $D$, and there is a 1-1 correspondence between the two. (More precisely there is a multiplicity of energies for a given $D$, each with a different number of nodes $n$: there is a 1-1 mapping between $D$ and the ($\epsilon, n$) pair.)

Solving the Dirac equation Radial Mesh

$g$ and $f$ must be integrated numerically on a radial mesh. These functions will oscillate very rapidly for small $r$ where they must be orthogonalized to core levels; thus meshes are typically logarithmic. Questaal uses a modified form: a shifted logarithmic mesh. For points $i=1{\cdots}m$, the mesh is $$r_i = b \cdot \exp(a(i-1)-1)$$

Since $r_m$ must be the augmentation radius, there is only one degree of freedom in the choice of $a$ and $b$.

When $a{\cdot}i$ is small compared to unity the $r_i$ are linearly spaced, with spacing $ab$; as $i$ increases it smoothly tums into a logarithmic mesh where $r_{i+1}/r_i$ is constant.

To control stability there are several considerations:

1. When $a{\cdot}i$ is much larger than one, the mesh is logarithmic with fixed $r_{i+1}/r_i=a$. If $a$ is too large, the integration near the augmentation radius is inaccurate.

2. When $a{\cdot}i$ is small compared to unity the $r_i$ are linearly spaced, with spacing $ab$ If $ab$ is too small, the integrator becomes inaccurte as numerical problems with points too closely spaced appears.

3. In a linear method not only the partial wave, but its energy derivative, must be integrated.

4. Core states with $\epsilon$ very negative decay exponentially fast for large $r$.

To address the last issue, the partial wave is integrating outward from the origin, and inward from $r_m$ to some central meeting point $r_c$. The logarithmic derivative of the two solutions must match at $r_c$, which is done by fixing $D$ and varying $\epsilon$, or the reverse. It is a nonlinear problem either way, and either $D$ or $\epsilon$ must be determined by successive iterations. (Most of the time in the Questaal codes, $D$ is fixed and $\epsilon$ is varied.)

Legacy and modern implementations

As for points 1, and 2, it has been found empirically that $a=0.015$ is quite reliable; this is the typical default value. (You can explicitly specify $a$ for a particular species in the input file.) $b$ is fixed by an internal formula (see subs/rmesh.f in the distribution). In the legacy code the formula was a simple function of the nuclear atomic number $Z$. However it can happen that $b$ is too small, especially when $a$ is small, and the mesh generator can adjust for this.

Regarding point 3, $\dot{\phi}$ (actually the relativistic analog $\dot{g},\,\dot{f}$) can be computed by finite difference at several energies, or by solving an inhomogenous Schrodinger equation.

Input file switches controlling how the integrations are performed: the RSEDEF token

To preserve backward compatibility, Questaal uses legacy algorithms, but it is generally ore accurate to use more modern replacements. Do this by adding this to the OPTIONS category : RSEDEF=0 in the input file. This does three things: (1) replaces the legacy Numerov solver with a Runge Kutta solver, (2) solves $\dot{\phi}$ with an inhomogenous equation, rather than finite differences of the homgenous one, and (3) safeguards against $ab$ becoming too small.

To set options individually, see the following table.