All-electron fully relativistic Kohn-Sham theory for solids based on the Dirac-Coulomb Hamiltonian and Gaussian-type functions

We present the first full-potential method that solves the fully relativistic four-component Dirac-Kohn-Sham equation for materials in the solid state within the framework of atom-centered Gaussian-type orbitals (GTOs). Our GTO-based method treats one-, two-, and three-dimensional periodic systems on an equal footing, and allows for a seamless transition to the methodology commonly used in studies of molecules with heavy elements. The scalar relativistic effects as well as the spin-orbit interaction are handled variationally. The full description of the electron-nuclear potential in the core region of heavy nuclei is straightforward due to the local nature of the GTOs and does not pose any computational difficulties. We show how the time-reversal symmetry and a quaternion algebra-based formalism can be exploited to significantly reduce the increased methodological complexity and computational cost associated with multiple wave-function components coupled by the spin-orbit interaction. We provide a detailed description of how to employ the matrix form of the multipole expansion and an iterative renormalization procedure to evaluate the conditionally convergent infinite lattice sums arising in studies of periodic systems. Next, we investigate the problem of inverse variational collapse that arises if the Dirac operator containing a repulsive periodic potential is expressed in a basis that includes diffuse functions, and suggest a possible solution. Finally, we demonstrate the validity of the method on three-dimensional silver halide ($\mathrm{Ag}X$) crystals with large relativistic effects and two-dimensional honeycomb structures (silicene and germanene) exhibiting the spin-orbit-driven quantum spin Hall effect. Our results are well-converged with respect to the basis set limit using standard bases developed for molecular calculations and indicate that the common rule of removing basis functions with small exponents should not be applied when transferring the molecular basis to solids.

Figure 1

Band structure diagram (left) and the DOS (right) of AgI obtained at the 4c (full line) and 1c (dashed line) levels of theory with the PBE XC functional. The horizontal dashed black line separates the occupied and the vacant states, and was placed in the middle of the band gap. The path traversing high-symmetry $\mathit{k}$ points in the reciprocal-space unit cell was chosen according to Ref. [141]. The figure was generated using python matplotlib library [142].

Figure 2

Speedup on the wall-clock time of the evaluation of the Coulomb and XC contributions to the Fock matrix in real space for 1 SCF cycle of the 4c AgI as a function of the number of CPUs used. The reference calculation was performed using 128 CPUs (4 nodes) and the largest calculation used 2048 CPUs (64 nodes). The dashed line denotes a hypothetical (linear) scaling of wall-clock time with given computational resources. The figure was generated using python matplotlib library [142].