Fourier interpolation for matrices in atomic basis

Introduction

In atomic basis framework, the matrix representation of operator in reciprocal space is calculated by Fourier transforming the real space matrices

$$\begin{array}{}\text{(1)}& \mathbf{A}(\mathbf{k})=\sum _{\mathbf{R}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}\mathbf{A}(\mathbf{R})\end{array}$$

where the matrix elements

$$\begin{array}{}\text{(2)}& \begin{array}{r}{A}_{ij}(\mathbf{R})=⟨{\phi }_{i,\mathbf{0}}|\hat{A}{\phi }_{j,\mathbf{R}}⟩=\int \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }{\phi }_{i,\mathbf{0}}(\mathbf{r})A(\mathbf{r},{\mathbf{r}}^{\prime }){\phi }_{j,\mathbf{R}}({\mathbf{r}}^{\prime })\end{array}\end{array}$$

and ${\phi }_{i,\mathbf{R}}$ is basis function located on atom $I$ in the unit cell marked by vector $\mathbf{R}$. For non-periodic $\mathbf{A}(\mathbf{R})$, as long as it decays in real-space, one can employ a cut-off radius $R_c$ (for isotropic system) where $\mathbf{A}(\left|\mathbf{R}\right|\ge R_c)$ is effectively zero. Then the summation can be reduced to

$$\begin{array}{}\text{(3)}& \mathbf{A}(\mathbf{k})=\sum _{\left|\mathbf{R}\right|<R_c}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}\mathbf{A}(\mathbf{R})\end{array}$$

Examples are overlap matrix and Hamiltonian matrix. For the latter, the batch-based integration in FHI-aims is discussed in a previous CN post.

Sometimes, rather than direct evaluation in real space using Eq. $\text{(2)}$, the real-space matrix is obtained by the inverse Fourier transform of its counterpart in reciprocal space computed on a finite set of $\mathbf{k}$ points. In this case, matrix $A(\mathbf{R})$ is periodic under the Born-von Karman (BvK) boundary condition. Note that this period is artificial due to finite $\mathbf{k}$ sampling. The unit cells within a period then forms the so-called BvK super cell. The summation of Eq. $\text{(1)}$ is then constrained within the BvK super cell

$$\begin{array}{}\text{(4)}& \mathbf{A}(\mathbf{k})=\sum _{\mathbf{R}\in \{{\mathbf{R}}^{\text{BvK}}\}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}\mathbf{A}(\mathbf{R})\end{array}$$

Note that the two equations become formally the same when the BvK cell encloses the $R_c$ sphere.

Fourier interpolation

For the periodic case, we can write explicitly the construction of $\mathbf{A}(\mathbf{R})$ by inverse transform using the matrices on mesh points $\{{\mathbf{k}}^{\prime }\}$, commensurate with the BvK cell in the Brillouin zone. By commensurate, it means that for BvK cell with period $(N_1,N_2,N_3)$, the mesh points are

$$\begin{array}{}\text{(5)}& \{{\mathbf{k}}^{\prime }=\sum _{i=1}^3\frac{n_i}{N_i}{\mathbf{b}}_i\mid n_i=0,1,\dots ,N_i-1\}\end{array}$$

where ${\mathbf{b}}_i$ are the reciprocal lattice vectors of the unit cell. Therefore,

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathbf{A}(\mathbf{k})=& \sum _{\mathbf{R}\in \{{\mathbf{R}}^{\text{BvK}}\}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}\mathbf{A}(\mathbf{R})\\ =& \sum _{\mathbf{R}\in \{{\mathbf{R}}^{\text{BvK}}\}}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot \mathbf{R}}\frac{1}{N_{{\mathbf{k}}^{\prime }}}\sum _{{\mathbf{k}}^{\prime }}{\mathrm{e}}^{-\mathrm{i}{\mathbf{k}}^{\prime }\cdot \mathbf{R}}\mathbf{A}({\mathbf{k}}^{\prime })\\ =& \sum _{{\mathbf{k}}^{\prime }}\mathbf{A}({\mathbf{k}}^{\prime })\frac{1}{N_{{\mathbf{k}}^{\prime }}}\sum _{\mathbf{R}\in \{{\mathbf{R}}^{\text{BvK}}\}}{\mathrm{e}}^{\mathrm{i}(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{R}}\\ =& \sum _{{\mathbf{k}}^{\prime }}F(\mathbf{k},{\mathbf{k}}^{\prime })\mathbf{A}({\mathbf{k}}^{\prime })\end{array}\end{array}$$

This can be viewed as an interpolation scheme for $\mathbf{A}(\mathbf{k})$ from $\{\mathbf{A}({\mathbf{k}}^{\prime })\}$ with interpolation coefficients

$$\begin{array}{}\text{(7)}& F(\mathbf{k},{\mathbf{k}}^{\prime })=\frac{1}{N_{{\mathbf{k}}^{\prime }}}\sum _{\mathbf{R}\in \{{\mathbf{R}}^{\text{BvK}}\}}{\mathrm{e}}^{\mathrm{i}(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{R}}\end{array}$$

and we name it as Fourier interpolation.

Choice of Born-von Karman cell

To compute the interpolation coefficients Eq. $\text{(7)}$, the BvK super cell has to be specified. Naively, it can be chosen the same for the whole matrix $\mathbf{A}$. To illustrate, we look at a 2-dimensional case with 3 by 3 BvK super cell as shown in figure 1. The unit cell is a square one, and contains 4 atoms with fractional coordinates

• I: (0.15, 0.10)
• J: (0.90, 0.70)
• K: (0.60, 0.40)
• L: (0.90, 0.30)

The BvK cell is chosen to have (0,0) cell at its center, marked by the blue square. The BvK super cell then contains unit cells $(n_1,n_2)$ with $n_1=0,\pm 1$ and $n_2=0,\pm 1$. The Fourier transform Eq. $\text{(4)}$ is then explicitly written as

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathbf{A}(\mathbf{k})=& \sum _{n_1=-1,0,1}\sum _{n_2=-1,0,1}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot (n_1{\mathbf{a}}_1+n_2{\mathbf{a}}_2)}\mathbf{A}(n_1{\mathbf{a}}_1+n_2{\mathbf{a}}_2)\\ \equiv & \sum _{n_1=-1,0,1}\sum _{n_2=-1,0,1}{\mathrm{e}}^{\mathrm{i}\mathbf{k}\cdot (n_1{\mathbf{a}}_1+n_2{\mathbf{a}}_2)}\mathbf{A}(n_1,n_2)\end{array}\end{array}$$

If we focus on block ${\mathbf{A}}_{IJ}$, whose row indices correspond to basis on atom $I$ and column indices on $J$, we can visualize the contributions $\mathbf{A}(n_1,n_2)$ by blue arrows in the figure.

Figure 1: 3 by 3 BvK cell centered at (0,0) unit cell. Blue arrows stand for the submatrices of $\mathbf{A}(\mathbf{R})$ with row index on atom $I$ and column on $J$. This choice of BvK cell, however, is problematic from physical consideration. Suppose that operator $\hat{A}$ is short-ranged, and a 3 by 3 BvK cell barely encloses the sphere of cut-off radius $R_c$. This indicates that ${\mathbf{A}}_{IJ}$ is non-zero only within the BvK cell. On the other hand, the atom $J$ in unit cell (1,-1) is farther from atom $I$ in the center cell than its periodic image in cell (-2,-1), which implies by the decay property of $\hat{A}$ that

$$\begin{array}{}\text{(9)}& \Vert {\mathbf{A}}_{IJ}(-2,-1)\Vert >\Vert {\mathbf{A}}_{IJ}(1,-1)\Vert >0\end{array}$$

Since cell (-2,-1) is outside of our choice of BvK cell, it contradicts the condition that ${\mathbf{A}}_{IJ}$ is non-zero only within the BvK cell. The same situation applies to (-1,1), (0,1), (1,1), and (1,0).

A more physical scheme is to choose the BvK cell where the nearest periodic images are included. For pair $IJ$, we can refine the BvK cell in figure 1 according to analysis above. The final BvK cell for $IJ$ is $(n_1,n_2)$ with $n_1=0,-1,-2$ and $n_2=0,-1,-2$, as illustrated in figure 2. On the other hand, for pair $KL$, the BvK cell in figure 1 does include the nearest images of $L$ and is hence an appropriate choice.

Figure 2: 3 by 3 BvK cell chosen for atom pairs $IJ$ (blue) and $KL$ (red). It is therefore indicated that this scheme will lead to different choices of BvK cells for each pair of atoms. We thus modify our interpolation formula Eq. $\text{(6)}$ to

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathbf{A}}_{IJ}(\mathbf{k})=& \sum _{{\mathbf{k}}^{\prime }}{\mathbf{A}}_{IJ}({\mathbf{k}}^{\prime })\frac{1}{N_{{\mathbf{k}}^{\prime }}}\sum _{\mathbf{R}\in \{{\mathbf{R}}_{IJ}^{\text{BvK}}\}}{\mathrm{e}}^{\mathrm{i}(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{R}}\\ =& \sum _{{\mathbf{k}}^{\prime }}{F}_{IJ}(\mathbf{k},{\mathbf{k}}^{\prime }){\mathbf{A}}_{IJ}({\mathbf{k}}^{\prime })\end{array}\end{array}$$

where

$$F_{IJ}(\mathbf{k},{\mathbf{k}}^{\prime })=\frac{1}{N_{{\mathbf{k}}^{\prime }}}\sum _{\mathbf{R}\in \{{\mathbf{R}}_{IJ}^{\text{BvK}}\}}{\mathrm{e}}^{\mathrm{i}(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{R}}$$

The modified scheme for choosing BvK cell is shown to give smoother band structure by Fourier interpolation in hybrid functional calculation in FHI-aims, given the same Brillouin zone sampling.

Summary

This note briefly discusses the Fourier interpolation technique to obtain matrices in reciprocal space at arbitrary wave vector from those at available $\mathbf{k}$ vectors. When the matrices are represented by atomic basis, the Born-von Karman cell used to compute the interpolation coefficient should be chosen for each atom pairs to make the interpolation of more physical meaning.

Discussion with Andrey Sobolev at MS1P is kindly appreciated.