A Wigner Matrix Based Convolution Algorithm For Matrix Elements in the LCAO Method

TL;DR

The work introduces the Wigner Matrix Convolution Algorithm (WMCA), a method to compute matrix elements in LCAO frameworks with arbitrary angular momentum by recasting the full crystal potential into a site-centered multipole expansion, turning each angular channel into a two-center integral evaluated via Wigner-$\mathcal{D}$ matrices. WMCA not only reproduces the Slater–Koster results for the traditional two-center approximation but also extends to higher angular momentum, enabling a tight-binding–like efficiency for full-potential calculations. By decomposing the potential into multipoles and applying convolution theorems, the approach promises near $O(n)$ scaling and applicability to ab-initio all-electron LCAO, with silicon as a test case using a model empirical potential. The results show rapid convergence with angular momentum cutoff and good agreement with grid references for valence bands, suggesting a path toward transferable, efficient ab-initio LCAO methods and potential extensions to total-energy calculations in future work.

Abstract

The linear combination of atomic orbitals (LCAO) method uses a small basis set in exchange for expensive matrix element calculations. The most efficient approximation for the matrix element calculations is the two-center approximation (2CA) in tight binding (TB). In the 2CA, a variety of matrix elements are neglected with only "two-center integrals" (2CI) remaining. The 2CI are calculated efficiently by rotating to symmetrical coordinates where the integral is parameterized. This makes TB fast in exchange for diminished transferability. An ideal electronic structure method has both the efficiency of TB and the transferability of ab-initio methods. In this work, I expand the full crystal potential into multipoles where the resulting matrix elements are transformed into the form of 2CI between high angular momentum functions. The usual Slater-Koster formulae for TB are limited to $l\leq3$; to enable efficient evaluation of the full crystal potential 2CI, I derive a Wigner matrix based convolution algorithm (WMCA) that works for arbitrary angular momentum. Given a suitable method for generating a local ab-initio Kohn-Sham potential, the algorithm for calculating matrix elements is applicable to fully ab-initio LCAO methods (this is the subject of forthcoming work). In this paper, I apply the WMCA to silicon using a model crystal potential.

Figures (10)

• Figure 1: A $p_y$ and $p_z$ orbital are centered at the origin and the vector $\vec{\Delta}$. The azimuthal and polar angles are $\phi$ and $\theta$ and the internuclear distance is $\Delta = |\vec{\Delta}|$. Since there is no "roll" angle associated with a coordinate vector, the Euler angles for a space-fixed $ZYZ$ rotation are (i) $\alpha=\phi$ around the $z$ axis and (ii) $\beta=\theta$ around the $y$ axis.
• Figure 2: The overlap integral between a $p_y$ and $p_z$ orbital with internuclear vector $\vec{\Delta}$ (cf. Fig. \ref{['fig:fig_1']}) can be decomposed into three azimuthally symmetric integrals by rotating the coordinate system so that $\vec{\Delta}=\Delta \vec{z}'$. The azimuthal symmetry conserves magnetic quantum number $M$ in the new coordinate system so that only three integrals survive, with $M=M' = -1 = p_{y'}$, $M=M' = 0 = p_{z'}$, and $M=M' = 1 = p_{x'}$. $p_{x'}$ and $p_{y'}$ are perpendicular to the vector and $p_{z'}$ is parallel to it. In the conventional notation, $M$ is denoted by $\sigma = 0$, $\pi = \pm1$, and $\delta = \pm 2$.
• Figure 3: Band structures of (a) silicon in the fcc-diamond phase mattheiss1981electronic and (b) tetragonal LaFeAsO papaconstantopoulos2010tight calculated in SK 2CA and (c) of fcc thorium using the NRL method durgavich2016extensionmehl1996applicationscohen1994tight. The parameters were taken from the references and the band structures were calculated using the WMCA. The silicon data points in (a) are the pseudopotential results chelikowsky1974electronic that the band structure was fit to. In (b) and (c), the data points are the digitized tight binding band structure calculated with the same parameters in the refs. papaconstantopoulos2010tightdurgavich2016extension. My calculation in (c) is projected onto the $d$ (red) and $f$ (blue) orbitals, showing the $f$ character of the electrons near the Fermi level ($E=0$ eV) in thorium.
• Figure 4: Comparison between the band structure of silicon calculated using the minimal silicon SZV-MOLOPT-SR-GTH basis vandevondele2007gaussian on a grid (black lines) and using a plane wave basis (magenta lines); both calculations used the empirical pseudopotential of Wang et al.wang1994electronic.
• Figure 5: Convergence of the multipole expansion of the silicon crystal potential with increasing angular momentum cutoff, $L_\textrm{cut}$. The conventional 8-atom unit cell is used ($a=5.43~\textrm{\AA}$). Radial functions with amplitude $\max |V_{LM}(r)| \geq 0.1 \times \max|V_{00}(r)|$ up to $(L,M)=(24,24)$ are shown in (a). The radial functions are offset vertically for clarity; the zero of each curve is represented by a dashed line. The truncated potential evaluated on a grid in the $x-y$ plane (top) is compared to the exact potential (lower left) in (b-d). The error, $|V(\vec{r}) - \sum_{LM} V_{LM}(r)X_{LM}(\Omega_{\vec{r}})|$, is shown on the lower right in (b-e). Note that the error is shown in log scale.