Higher-order finite difference method for accurate calculation of Wannier centers and position matrix elements

An accurate calculation of Wannier function centers and position matrix elements is crucial for calculating various physical quantities, such as electric polarization, orbital magnetization, and optical responses, from first principles. In practice, approximations can make convergence with respect to the ab initio k-grid significantly slow. Since the Fourier transform of the position operator is the gradient in momentum space, a finite difference technique is used to calculate the position matrix elements for Wannier functions. However, so far, only the first-order finite difference has been considered therein. We present a higher-order finite-difference method and demonstrate its advantage over the conventional first-order finite difference method.