High-order finite element method for atomic structure calculations

Substituting (11) into (12) and integrating, we obtain

$$\lim_{r\to\infty}r^{2}V_{H}^{\prime}(r)=-N,$$

(13)

from which it follows that the asymptotic behavior of $V_{H}^{\prime}(r)$ is

$$V_{H}^{\prime}(r)\sim-\frac{N}{r^{2}},\qquad r\to\infty.$$

(14)

Integrating (14) and requiring $V_{H}\to 0$ as $r\to\infty$ then gives the corresponding asymptotic behavior for $V_{H}(r)$,

$$V_{H}(r)\sim\frac{N}{r},\qquad r\to\infty.$$

(15)

For small $r$, the asymptotic behavior can be obtained by expanding $n(r)$ about $r=0$: $n(r)=\sum_{j=0}^{\infty}{c_{j}r^{j}}$. Substituting into Poisson equation (11) gives

$$(r^{2}V_{H}^{\prime})^{\prime}=-4\pi\sum_{j=0}^{\infty}{c_{j}r^{j+2}}.$$

(16)

Integrating and requiring $V_{H}(0)$ to be finite then gives

$$V_{H}^{\prime}(r)=-4\pi\sum_{j=0}^{\infty}{c_{j}\frac{r^{j+1}}{j+3}},$$

(17)

with linear leading term, so that we have

$$V_{H}^{\prime}(r)\sim r,\qquad r\to 0.$$

(18)

Integrating (17) then gives

$$V_{H}(r)=-4\pi\sum_{j=0}^{\infty}{c_{j}\frac{r^{j+2}}{(j+2)(j+3)}}+C,$$

(19)

with leading constant term $C=V_{H}(0)$ determined by Coulomb’s law:

$$V_{H}(0)=4\pi\int_{0}^{\infty}{rn(r)\differential r}.$$

(20)

Finally, from (18) we have that

$$V_{H}^{\prime}(0)=0.$$

(21)

So $V_{H}\to 0$ as $r\to\infty$ and $V_{H}^{\prime}(0)=0$. This asymptotic behavior provides the initial values and derivatives for numerical integration in both inward and outward directions.

The known asymptotic behavior of $P_{nl}$ as $r\to 0$ is [33]

$$P_{nl}(r)\sim r^{l+1},$$

(42)

where $P_{nl}$, being a solution of a homogeneous system of equations, is only unique up to an arbitrary multiplicative constant. We use (40) as our starting point but from now on drop the $nl$ index from $P_{nl}$ for simplicity:

$$-\frac{1}{2}P^{\prime\prime}(r)+\left(V(r)+\frac{l(l+1)}{2r^{2}}\right)P(r)=EP(r).$$

(43)

In order to facilitate rapid convergence and application of desired boundary conditions in a finite element basis, we generalize (40) to solve for $\tilde{P}=\frac{P(r)}{r^{\alpha}}$ for any chosen real exponent $\alpha\geq 0$ by substituting $P=r^{\alpha}\tilde{P}$ to obtain

$$-\frac{1}{2}{1\over r^{2\alpha}}\left(r^{2\alpha}\tilde{P}^{\prime}(r)\right)^{\prime}+\left(V(r)+\frac{l(l+1)-\alpha(\alpha-1)}{2r^{2}}\right)\tilde{P}(r)=E\tilde{P}(r).$$

(44)

We note that:

• •

  For $\alpha=0$ we get $\tilde{P}(r)=\frac{P(r)}{r^{0}}=P(r)$ and (44) reduces to (40).

• •

  For $\alpha=1$ we get $\tilde{P}(r)=\frac{P(r)}{r^{1}}=R(r)$ and (44) reduces to (4).

• •

  Finally, for $\alpha=l+1$, (which corresponds to the known asymptotic (42)) we obtain $\tilde{P}(r)=\frac{P(r)}{r^{l+1}}$, which tends to a nonzero value at the origin for all $l$ and (44) becomes

  $$-\frac{1}{2}{1\over r^{2(l+1)}}\left(r^{2(l+1)}\tilde{P}^{\prime}(r)\right)^{\prime}+V(r)\tilde{P}(r)=E\tilde{P}(r).$$

  (45)

To obtain the weak form, we multiply both sides of (44) by a test function $v(r)$ and integrate from $0$ to $\infty$. In addition, to facilitate the construction of a symmetric bilinear form, we multipy by a factor $r^{2\alpha}$ to get

	$\displaystyle\int_{0}^{\infty}\left[-\frac{1}{2}\left(r^{2\alpha}\tilde{P}^{\prime}(r)\right)^{\prime}v(r)+\left(V(r)+\frac{l(l+1)-\alpha(\alpha-1)}{2r^{2}}\right)\tilde{P}(r)v(r)r^{2\alpha}\right]\differential r$		(46)
	$\displaystyle=E\int_{0}^{\infty}\tilde{P}(r)v(r)r^{2\alpha}\differential r.$

We can now integrate by parts to obtain

	$\displaystyle\int_{0}^{\infty}\left[\frac{1}{2}r^{2\alpha}\tilde{P}^{\prime}(r)v^{\prime}(r)+\left(V(r)+\frac{l(l+1)-\alpha(\alpha-1)}{2r^{2}}\right)\tilde{P}(r)v(r)r^{2\alpha}\right]\differential r$		(47)
	$\displaystyle-{\textstyle\frac{1}{2}}\left[r^{2\alpha}\tilde{P}^{\prime}(r)v(r)\right]_{0}^{\infty}=E\int_{0}^{\infty}\tilde{P}(r)v(r)r^{2\alpha}\differential r.$

Setting the boundary term to zero,

$$\left[r^{2\alpha}\tilde{P}^{\prime}(r)v(r)\right]_{0}^{\infty}=0,$$

(48)

we then obtain the desired symmetric weak formulation

	$\displaystyle\int_{0}^{\infty}\left[{\textstyle\frac{1}{2}}\tilde{P}^{\prime}(r)v^{\prime}(r)+\left(V(r)+\frac{l(l+1)-\alpha(\alpha-1)}{2r^{2}}\right)\tilde{P}(r)v(r)\right]r^{2\alpha}\differential r$		(49)
	$\displaystyle=E\int_{0}^{\infty}\tilde{P}(r)v(r)r^{2\alpha}\differential r.$

As discussed below, by virtue of our choice of $\alpha$ and boundary conditions on $v(r)$, the vanishing of the boundary term (48) imposes no natural boundary conditions on $\tilde{P}$.

To discretize (49), we introduce finite element basis functions $\phi_{i}(r)$ to form trial and test functions

$\displaystyle\tilde{P}(r)$

$\displaystyle=\sum_{j=1}^{N}c_{j}\phi_{j}(r)\qquad\mathrm{and}\qquad v(r)=\phi_{i}(r),$

(50)

and substitute into (49). In so doing, we obtain a generalized eigenvalue problem

$$\sum_{j=1}^{N}H_{ij}c_{j}=E\sum_{j=1}^{N}S_{ij}c_{j},$$

(51)

where $H_{ij}$ is the Hamiltonian matrix,

$$H_{ij}=\int_{0}^{\infty}\left[{\textstyle\frac{1}{2}}\phi_{i}^{\prime}(r)\phi_{j}^{\prime}(r)+\phi_{i}(r)\left(V+\frac{l(l+1)-\alpha(\alpha-1)}{2r^{2}}\right)\phi_{j}(r)\right]r^{2\alpha}\differential r,$$

(52)

and $S_{ij}$ is the overlap matrix,

$$S_{ij}=\int_{0}^{\infty}\phi_{i}(r)\phi_{j}(r)r^{2\alpha}\differential r.$$

(53)

While, by virtue of the weak formulation, the above discretization can employ any basis in Sobolev space $H^{1}$ satisfying the required boundary conditions (including standard $C^{0}$ FE bases, $C^{1}$ Hermite, and $C^{k}$ B-splines), we employ a high-order $C^{0}$ spectral element (SE) basis [34,35] in the present work. FE bases consist of local piecewise polynomials. SE bases are a form of FE basis constructed to enable well-conditioned, high polynomial-order discretization via definition over Lobatto nodes rather than uniformly spaced nodes within each finite element (mesh interval). In the present work, we employ polynomial orders up to $p=31$.

In the electronic structure context, SE bases have a number of desirable properties, including:

1. 1.

   Polynomial completeness and associated systematic convergence with increasing number and/or polynomial order of basis functions.

2. 2.

   Well conditioned to high polynomial order by virtue of definition over Lobatto nodes.

3. 3.

   Exponential convergence with respect to polynomial order, enabling high accuracy with a small basis.

4. 4.

   Applicable to general nonsingular as well as singular Coulombic potentials, whether fixed or self-consistent.

5. 5.

   Applicable to bound as well as excited states.

6. 6.

   Can use uniform as well as nonuniform meshes.

7. 7.

   Basis function evaluation, differentiation, and integration are inexpensive, enabling fast and accurate matrix element computations.

8. 8.

   Dirichlet boundary conditions are readily imposed by virtue of cardinality, i.e.,

   $$\phi_{i}(x_{j})=\delta_{ij}$$

   where $x_{j}$ is a nodal point of the basis.

9. 9.

   Derivative boundary conditions are readily imposed via the weak formulation.

10. 10.

    Diagonal overlap matrix using Lobatto quadrature, thus enabling solution of a standard rather than generalized eigenvalue problem, reducing both storage requirements and time to solution.

Properties (1)-(4), (7), and (10) are advantageous relative to Slater-type and Gaussian-type bases, which can be nontrivial to converge; can become ill-conditioned as the number of basis functions is increased, limiting accuracies attainable in practice; and produce a non-diagonal overlap matrix, requiring solution of a generalized rather than standard matrix eigenvalue problem.

Properties (2) and (10) are advantageous relative to B-spline bases, which are generally employed at lower polynomial orders in practice and give rise to a non-diagonal overlap matrix, thus requiring solution of a generalized rather than standard matrix eigenvalue problem.

Since we seek bound states, which vanish as $r\to\infty$, we impose a homogeneous Dirichlet boundary condition on $v(r)$ and $\tilde{P}(r)$ at $r=\infty$ by employing a finite element basis $\{\phi_{i}\}$ satisfying this condition.

For $\alpha=0$, $\tilde{P}(r)=P(r)=0$ at $r=0$ for all quantum numbers $l$ according to (42). Thus we impose a homogeneous Dirichlet boundary condition on $v(r)$ and $\tilde{P}(r)$ at $r=0$ by employing a finite element basis $\{\phi_{i}\}$ satisfying this condition, whereupon the boundary term (48) as a whole vanishes, consistent with the weak formulation (49).

For $\alpha=1$, $\tilde{P}(r)=\frac{P(r)}{r}=R(r)=0$ at $r=0$ for all quantum numbers $\ell>0$ according to (42). However, $\tilde{P}(r)\neq 0$ at $r=0$ for $\ell=0$. Thus a homogeneous Dirichlet boundary condition cannot be imposed for $\ell=0$. However, for $\alpha>0$, the $r^{2\alpha}$ factor in (48) vanishes at $r=0$, whereupon the boundary term as a whole vanishes, consistent with the weak formulation (49), regardless of boundary condition on $v(r)$ and $\tilde{P}(r)$ at $r=0$. Hence for $\alpha>0$, we impose a homogeneous Dirichlet boundary condition on $v(r)$ and $\tilde{P}(r)$ at $r=\infty$ only, by employing a finite element basis $\{\phi_{i}\}$ satisfying this condition. This is sufficient due to the singularity of the associated Sturm-Liouville problem.

Similarly, for $\alpha=\ell+1$, $\tilde{P}(r)=\frac{P(r)}{r^{\ell+1}}\neq 0$ at $r=0$ for all quantum numbers $\ell$ according to (42) and, since $\alpha>0$, we again impose a homogeneous Dirichlet boundary condition on $v(r)$ and $\tilde{P}(r)$ at $r=\infty$ only.

In practice, the featom code defaults to $\alpha=0$ with homogeneous Dirichlet boundary conditions at $r=0$ and $r=\infty$.

Note that, while $\alpha=0$ yields a rapidly convergent formulation in the context of the Schrödinger equation, it does not suffice in the context of the Dirac equation, as discussed in Section 3.2.4.

The above eigenvalue formulation of the radial Dirac equation can be solved using the finite element method. However, due to the fact that the operator is not bounded from below, one obtains spurious eigenvalues in the spectrum [10]. To eliminate spurious states, we work with the square of the Hamiltonian [22, 15], which is bounded from below, rather than the Hamiltonian itself. Since the eigenfunctions of the square of an operator are the same as those of the operator itself, and the eigenvalues are the squares, the approach is straightforward and enables direct and efficient solution by the finite element method.

Let us derive the equations for the squared radial Dirac Hamiltonian. First we shift the energy by $c^{2}$ to obtain the relativistic energy, making the Hamiltonian more symmetric, using (58) to get

$$H+\mathbb{I}c^{2}=\begin{pmatrix}V(r)+c^{2}&c\left(-{\derivative{r}}+{\kappa\over r}\right)\\ c\left({\derivative{r}}+{\kappa\over r}\right)&V(r)-c^{2}\\ \end{pmatrix}.$$

(61)

Then we square the Hamiltonian using (59) to obtain the following equations to solve for $P(r)$ and $Q(r)$:

$$(H+\mathbb{I}c^{2})^{2}\begin{pmatrix}P(r)\\ Q(r)\\ \end{pmatrix}=(E+c^{2})^{2}\begin{pmatrix}P(r)\\ Q(r)\\ \end{pmatrix},$$

(62)

where $\mathbb{I}$ is the $2\times 2$ identity matrix.

As can be seen, the eigenvectors $P(r)$ and $Q(r)$ are the same as before but the eigenvalues are now equal to $(E+c^{2})^{2}$ and the original non-relativistic energies $E$ can be obtained by taking the square root of these new eigenvalues and subtracting $c^{2}$.

We follow a similar approach as for the Schrödinger equation: instead of solving for $P(r)$ and $Q(r)$, we solve for $\tilde{P}(r)={P(r)\over r^{\alpha}}$ and $\tilde{Q}(r)={Q(r)\over r^{\alpha}}$, which introduces a parameter $\alpha$ that can be chosen to facilitate rapid convergence and the application of the desired boundary conditions in a finite element basis. We then multiply the eigenvectors by $r^{\alpha}$ to obtain $P(r)$ and $Q(r)$.

Now we can write the finite element formulation as

	$\displaystyle A$	$\displaystyle=\int_{0}^{\infty}\!\!\begin{pmatrix}\phi_{i}^{a}(r)&\phi_{i}^{b}(r)\\ \end{pmatrix}r^{\alpha}(H+\mathbb{I}c^{2})^{2}r^{\alpha}\begin{pmatrix}\phi_{j}^{a}(r)\\ \phi_{j}^{b}(r)\end{pmatrix}\differential r,$		(63a)
	$\displaystyle S$	$\displaystyle=\int_{0}^{\infty}\begin{pmatrix}\phi_{i}^{a}(r)&\phi_{i}^{b}(r)\\ \end{pmatrix}r^{\alpha}r^{\alpha}\begin{pmatrix}\phi_{j}^{a}(r)\\ \phi_{j}^{b}(r)\end{pmatrix}\differential r,$		(63b)
	$\displaystyle Ax$	$\displaystyle=(E+c^{2})^{2}Sx.$		(63c)

This is a generalized eigenvalue problem, with eigenvectors $x$ (coefficients of $(\tilde{P}(r),\tilde{Q}(r))$), eigenvalues $(E+c^{2})^{2}$, and $2\times 2$ block matrices $A$ and $S$. Note that the basis functions from both sides were multiplied by $r^{\alpha}$ due to the substitutions $P(r)=r^{\alpha}\tilde{P}(r)$ and $Q(r)=r^{\alpha}\tilde{Q}(r)$. We can denote the middle factor in the integral for $A$ as $H^{\prime}=r^{\alpha}\left(H+\mathbb{I}c^{2}\right)^{2}r^{\alpha}$ and compute it using (61) with rearrangement of $r^{\alpha}$ factors as follows:

	$\displaystyle H^{\prime}$	$\displaystyle=r^{\alpha}\left(H+\mathbb{I}c^{2}\right)^{2}r^{\alpha}$		(64a)
		$\displaystyle=r^{2\alpha}\begin{pmatrix}V(r)+c^{2}&c\left(-{\derivative{r}}+{\kappa-\alpha\over r}\right)\\ c\left({\derivative{r}}+{\kappa+\alpha\over r}\right)&V(r)-c^{2}\\ \end{pmatrix}^{2}.$		(64b)