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 on a 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. Usually 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 having an issue and looking for a quick fix, include the RSEDEF token in your ctrl file as explained below.

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 coupling 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 on a Radial Mesh

$g$ and $f$ must be integrated numerically on a radial mesh. Usually this is handled fairly well with internal defaults. The default legacy implementations has been modernized in several ways, which you can apply as a bundle by adding OPTIONS_RSEDEF=0 to your ctrl file.

The following explains the process in more detail, together with the modernizations available.

$g$ and $f$ 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{\dots}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 $a{\cdot}b$; as $i$ increases it smoothly turns into a logarithmic mesh where $r_{i+1}/r_i$ is constant.

To control stability, several overlapping considerations must be taken into account.

1. When $a{\cdot}i$ is much larger than unity, 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 $a{\cdot}b$ (i.e. $r_{i+1}-r_i=ab$). If $ab$ is too small, points are too closely spaced and the integrator becomes inaccurate.

3. In a linear method not only the partial wave $\phi(r)$, but its energy derivative $\dot{\phi}(r)$, must be integrated. This can be solved by finding $\phi(r)$ at multiple energies and differentiating by finite difference (what the legacy mode does), or by solving an inhomogeneous Dirac equation.

4. Core states with $\epsilon$ very negative decay exponentially with $r$ for large $r$ (think of the H atom). This is the classic instance of a stiff differential equation.

To address the stiffness issue, $\phi$ is integrated 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 in an iterative manner. (In the Questaal codes, usually $D$ is fixed and $\epsilon$ is varied. But there can be exceptions.)

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 should adjust for this. In the legacy code this issue was not properly taken into account and choosing $a$ too small causes problems.

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 inhomogeneous Dirac 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 more accurate to use some modern improvements to them. They are included as a bundle when you put RSEDEF=0 in the OPTIONS category of the input file.

It does three things: (1) replaces the legacy Numerov solver with a Runge Kutta solver, (2) solves $\dot{\phi}$ with an inhomogeneous equation, rather than finite differences of the homogeneous one, and (3) safeguards against $ab$ becoming too small, in effect setting OPTIONS_RBFAC=-1, as explained below.

To set options individually, see the following table.

• OPTIONS_RBFAC=1 is equivalent to the default legacy setting. It does not safeguard against $a{\cdot}b$ becoming too small.

• OPTIONS_RBFAC=fac > 1 reduces the total number of points, preserving the logarithmic spacing at large $r$ but increases the product $ab$ by fac. It is not advised to choose fac as a positive number less than 1.

• OPTIONS_RBFAC=fac < 0 adjusts $b$ by choosing a canonical default value for it, such that it scales as $1/a$ so as to keep the product $ab$ fixed. OPTIONS_RBFAC=−1 is intended to be select a good default value of $b$ for a wide range of $a$. The figure shows $r_i$ for Au with $a{=}0.01$. The blue points are for the default mesh; the red ones were generated with OPTIONS_RBFAC=-1. The have the same spacing for large $r$, but at small $r$ the blue points are more coarsely spaced, so $r_i$ initially increases more quickly with $i$. The default mesh has 1212 points, while the red mesh has 1102 points.