High-order finite element method for atomic structure calculations

This paper introduces \texttt{featom}, an open-source, high-order finite element Fortran code that achieves high-accuracy solutions for radial Schrödinger, Dirac, and Kohn-Sham equations in heavy atoms through systematic convergence control and specialized techniques for handling singularities, while demonstrating significant speed improvements over existing methods.

The Gist

Imagine you are trying to build a perfect model of a single atom, like a tiny solar system where electrons orbit a nucleus. To do this, scientists use complex mathematical equations (the Schrödinger and Dirac equations) to predict exactly where the electrons are and how much energy they have.

For decades, solving these equations has been like trying to navigate a ship through a stormy sea using a map that keeps changing. Some methods are fast but inaccurate; others are accurate but take forever to compute.

This paper introduces a new tool called featom (Finite Element Atom). Think of featom as a high-tech, modular LEGO set that allows scientists to build atomic models with incredible precision and speed.

Here is a breakdown of how it works, using simple analogies:

1. The Problem: The "Spooky" Ghosts

In the world of quantum physics, there's a tricky equation called the Dirac equation used for heavy atoms (like Uranium). When scientists try to solve this with old methods, the math sometimes creates "ghosts." These are fake solutions (called spurious states) that look like real electrons but don't actually exist. It's like your GPS telling you there's a bridge where there is only a cliff.

The featom Solution:
Instead of trying to fix the bridge, the authors decided to square the map.

• Analogy: Imagine you have a number that could be positive or negative (like $+5$ or $-5$). If you square it, you get $25$ either way. By "squaring" the mathematical operator in their code, they force all the numbers to be positive. This eliminates the "ghosts" naturally because the math no longer allows for those fake, negative-energy solutions. It's a clever trick that makes the problem stable and easy to solve.

2. The Problem: The "Sharp Corner" at the Center

At the very center of an atom (the nucleus), the math gets messy. The equations behave like a sharp, jagged spike. If you try to draw a smooth curve (a polynomial) through a sharp spike, you get a terrible approximation. It's like trying to draw a smooth circle around a jagged rock; no matter how many lines you add, it never looks quite right near the rock.

The featom Solution:
The authors realized that the "spike" follows a known pattern. So, instead of trying to draw the jagged rock itself, they peeled off the jagged layer first.

• Analogy: Imagine you are trying to measure the shape of a mountain that has a sharp peak. Instead of measuring the whole mountain, you measure the "slope" of the mountain after you've mathematically removed the sharp tip. This leaves you with a smooth, gentle hill that is very easy to measure with standard tools.
• By doing this, their code can handle even the heaviest atoms (like Uranium) with extreme accuracy, reaching a precision of 1 part in 100 million.

3. The Method: The "Smart LEGO" Blocks

Most computer programs use a grid of points to solve these equations, kind of like a low-resolution photo. If you want a clearer picture, you have to add millions of points, which slows the computer down.

featom uses High-Order Finite Elements.

• Analogy: Imagine building a wall.
  • Old Method: You use thousands of tiny, flat bricks (low resolution). To make the wall smooth, you need millions of bricks.
  • featom Method: You use a few large, flexible, curved LEGO blocks. Each block can bend and twist to fit the shape of the wall perfectly.
• Because these "blocks" are so smart (using high-order polynomials), you need far fewer of them to get a perfect picture. This makes the calculation much faster.

4. The Results: Speed and Accuracy

The authors tested featom against the current best software (dftatom).

• Accuracy: It achieved the same ultra-high precision (10⁻⁸ Hartree) as the best existing tools. This is accurate enough to verify the most difficult benchmarks in physics.
• Speed: It was significantly faster. On a modern laptop, it solved a complex Uranium atom calculation in 28 milliseconds (for non-relativistic) compared to 166 milliseconds for the old method. That's like running a race in 10 seconds instead of 60.

Why Does This Matter?

This isn't just about solving a math puzzle.

• Modularity: The code is written like a set of interchangeable tools. If a scientist wants to change the rules of the game (the physics), they can swap out a "module" without breaking the whole machine.
• Open Source: The code is free for everyone to use, modify, and improve.
• Future Impact: Because it's so fast and accurate, it can be used as a building block for even larger simulations, helping us design new materials, understand nuclear reactions, or create better batteries.

In a nutshell:
The authors built a super-smart, flexible calculator that strips away the messy, jagged parts of atomic math before solving it. This allows them to see the atom with crystal-clear vision, much faster than anyone else, and without any "ghosts" confusing the results.

Technical Summary

1. Problem Statement

The paper addresses the computational challenges associated with solving the electronic structure equations (Schrödinger, Dirac, and Kohn-Sham) for isolated atoms, which are fundamental to Density Functional Theory (DFT) and materials science. Specific difficulties include:

• Spurious States in Relativistic Calculations: Solving the Dirac equation directly often leads to unphysical "spurious" eigenvalues due to the unbounded nature of the Dirac Hamiltonian and inconsistencies in discretization.
• Convergence Issues at the Origin: For heavy atoms (high atomic number $Z$), the derivatives of the wavefunctions diverge at the origin ($r=0$), particularly for states with relativistic quantum numbers $\kappa = \pm 1$. This non-polynomial behavior causes slow convergence in standard polynomial basis sets.
• Efficiency vs. Accuracy: Existing high-accuracy solvers (like the shooting method used in dftatom) can be computationally expensive or require careful tuning of grid parameters to maintain robustness.

2. Methodology

The authors introduce featom, an open-source Fortran 2008 code that implements a high-order finite element (FE) solver. The methodology relies on several key innovations:

A. High-Order Spectral Element (SE) Basis

Instead of low-order finite differences or standard finite elements, the code uses a high-order $C^0$ spectral element basis defined over Lobatto nodes.

• Benefits: This allows for exponential convergence with respect to the polynomial order ($p$), meaning high accuracy can be achieved with relatively few degrees of freedom.
• Flexibility: It supports various mesh types (uniform, exponential, or custom nodal distributions) and handles both singular (Coulombic) and non-singular potentials.

B. Squared Hamiltonian Formulation (for Dirac Equation)

To eliminate spurious states without modifying the Hamiltonian or imposing complex boundary constraints (like kinetic balance):

• The authors solve the square of the Dirac Hamiltonian ($(H + c^2)^2$) rather than the Hamiltonian itself.
• Rationale: The square of the operator is bounded from below, making it amenable to standard variational methods (like FE) without generating spurious states. The eigenfunctions remain identical to the original problem, while the eigenvalues are simply squared (recoverable via square root).
• Advantage: This approach preserves the correct non-relativistic limit exactly and stabilizes the numerics by naturally introducing second-derivative terms.

C. Asymptotic Transformation

To address the divergence of derivatives at the origin for $\kappa = \pm 1$ states:

• Instead of solving directly for the wavefunction components $P(r)$ and $Q(r)$, the code solves for transformed variables $\tilde{P}(r) = P(r)/r^\alpha$ and $\tilde{Q}(r) = Q(r)/r^\alpha$.
• The exponent $\alpha$ is chosen based on the known asymptotic behavior of the wavefunctions as $r \to 0$.
• Result: This transformation removes the non-polynomial behavior and derivative divergences at the origin, allowing the polynomial basis to converge rapidly for all quantum numbers.

D. Numerical Implementation

• Quadrature: Uses Gauss-Jacobi quadrature for integrals involving non-integer exponents (Dirac/Poisson) and Gauss-Legendre for standard cases. Gauss-Lobatto quadrature is used for the overlap matrix to ensure it is diagonal, converting the generalized eigenvalue problem into a standard one.
• SCF Loop: Implements self-consistent field (SCF) iterations with adaptive linear mixing and Periodic Pulay mixing to accelerate convergence.
• Initial Guess: Uses a Thomas-Fermi approximation for the initial density and potential.

3. Key Contributions

1. Robust Spurious State Elimination: Demonstrates that squaring the Hamiltonian is a simple, robust, and exact method to remove spurious states in Dirac calculations, avoiding the complexity of mixed-basis or kinetic-balance approaches.
2. Handling Singular Potentials: Successfully resolves the convergence bottleneck at the origin for heavy atoms by incorporating asymptotic forms into the basis functions, enabling high accuracy for $\kappa = \pm 1$ states.
3. Open-Source Modular Code: Provides featom, a modular, portable, and well-documented Fortran 2008 library under an MIT license. It is designed for extensibility and interoperability within the modern Fortran ecosystem (using fpm).
4. Comprehensive Benchmarking: Validates the method against analytic solutions and the state-of-the-art dftatom solver for atomic numbers $Z=1$ to $92$.

4. Results

• Accuracy: The code achieves total energy and eigenvalue accuracies of $10^{-8}$ Hartree (atomic units) for heavy atoms like Uranium ($Z=92$) in both Schrödinger and Dirac-Kohn-Sham formulations.
• Convergence:
  • Exponential Convergence: Error decreases exponentially with increasing polynomial order ($p$) until reaching machine precision ($\sim 10^{-9}$ to $10^{-10}$).
  • Mesh Independence: The method achieves target accuracy with fewer elements ($N_e$) compared to lower-order methods.
• Performance:
  • Speedup: In benchmarks on an Apple M1 Max processor, featom showed significant speedups over dftatom for non-relativistic (Schrödinger) calculations (28 ms vs. 166 ms).
  • Relativistic: For Dirac calculations, featom was competitive (360 ms vs. 276 ms for dftatom), though the overhead of the squared Hamiltonian formulation is noted.
• Validation: Results for Uranium and other elements match known analytic results and existing benchmarks, verifying the correctness of the implementation.

5. Significance

This work provides a general, high-accuracy, and efficient framework for atomic structure calculations. By combining high-order spectral elements with the squared Hamiltonian approach and asymptotic transformations, featom overcomes long-standing numerical instabilities in relativistic quantum chemistry.

Its significance lies in:

• Reliability: Offering a "black-box" solution that does not require manual tuning of grid parameters for different elements.
• Scalability: The exponential convergence makes it highly efficient for high-precision requirements.
• Accessibility: As an open-source tool, it lowers the barrier for researchers to perform high-fidelity atomic calculations, serving as a reliable reference for validating other electronic structure methods and as a component for larger-scale DFT codes.