High-order finite element method for atomic structure calculations

Abstract

We introduce featom, an open source code that implements a high-order finite element solver for the radial Schrödinger, Dirac, and Kohn-Sham equations. The formulation accommodates various mesh types, such as uniform or exponential, and the convergence can be systematically controlled by increasing the number and/or polynomial order of the finite element basis functions. The Dirac equation is solved using a squared Hamiltonian approach to eliminate spurious states. To address the slow convergence of the κ=±1 states due to divergent derivatives at the origin, we incorporate known asymptotic forms into the solutions. We achieve a high level of accuracy (10−8 Hartree) for total energies and eigenvalues of heavy atoms such as uranium in both Schrödinger and Dirac Kohn-Sham solutions. We provide detailed convergence studies and computational parameters required to attain commonly required accuracies. Finally, we compare our results with known analytic results as well as the results of other methods. In particular, we calculate benchmark results for atomic numbers (Z) from 1 to 92, verifying current benchmarks. We demonstrate significant speedup compared to the state-of-the-art shooting solver dftatom. An efficient, modular Fortran 2008 implementation, is provided under an open source, permissive license, including examples and tests, wherein particular emphasis is placed on the independence (no global variables), reusability, and generality of the individual routines. Program summaryProgram Title: featomCPC Library link to program files:https://doi.org/10.17632/962fzmm7f7.1Licensing provisions: MITProgramming language: Fortran with command-line interfacesNature of problem: Solution of the Schrödinger, Dirac, and Kohn–Sham equations of density functional theory for isolated atoms.Solution method: A high-order finite element method is used to assemble a generalized eigenproblem that is then solved for eigenvalues and orbitals. The Poisson equation is solved using the finite element method also. Self-consistent field equations are solved using Pulay mixing.Additional comments including restrictions and unusual features: Our solution methodology for the Dirac equation does not suffer from spurious states. We maintain high accuracy even for κ=±1 states. Unlike other methods, our approach is not limited to Coulombic or self-consistent potentials and can handle non-uniform meshes, including exponential meshes.Restrictions: Spherical symmetry.