First-principles calculation of optical responses based on nonorthogonal localized orbitals

Based on ab initio software packages using nonorthogonal localized orbitals, we develop a general scheme of calculating response functions. We test the performance of this method by calculating nonlinear optical responses of materials, like the shift current conductivity of monolayer WS2, and achieve good agreement with previous calculations. This method bears many similarities to Wannier interpolation, which requires a challenging optimization of Wannier functions due to the conflicting requirements of orthogonality and localization. Although computationally heavier compared to Wannier interpolation, our procedure avoids the construction of Wannier functions and thus enables automated high throughput calculations of linear and nonlinear responses related to electrical, magnetic and optical material properties.


Introduction
Wannier functions encode all the information of a band structure for a pre-selected energy window. Moreover, due to their highly localized nature, Wannier functions can be constructed with Bloch functions at just a few reciprocal k points [1,2]. After obtaining Wannier functions, information of Bloch functions at arbitrary k points can be obtained by Fourier transformation. This procedure is called Wannier interpolation [3]. Wannier interpolation captures the physics of Bloch electrons and is an efficient way to compute various physical observables of solids from first-principles calculations. Previous researches have proven its success on studying anomalous Hall effect [4], optical conductivity [5], orbital magnetization [6] and nonlinear optical effects including shift current [7,8] and second harmonic generation [7].
Current dominating algorithm of obtaining Wannier functions is maximally localized Wannier function (MLWF) theory, which mixes bands in a way that maximizes the localization of Wannier functions [3]. This method is universal and serves as a post processing tool for various ab initio calculation packages [9], regardless of their implementation schemes of density functional theory (DFT). However, since MLWF algorithm is essentially an optimization process, to avoid local minimum, trial Wannier functions and other parameters have to be carefully tested. This inevitable human intervention hampers its usage in high throughput materials discovery. Moreover, the optimization of Wannier functions usually breaks material symmetries that are essential to determine various material properties, including nonlinear optical responses and topological properties. One important fact is that many existing DFT packages already use localized orbitals to span the Hilbert space. These orbitals are usually optimized before the self consistent calculation and inherit the symmetry of atomic orbitals. Therefore, if the ideas of Wannier interpolation can be borrowed to use these Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. orbitals, the disadvantages of Wannier interpolation mentioned above can be eliminated. However, localized orbitals adopted by most software packages are nonorthogonal and do not directly fit into the framework of Wannier interpolation. This practice is rationalized by two practical reasons: (i) nonorthogonal localized orbitals (NoLO) can be made more localized than orthogonal orbitals; (ii) It is generally easier to construct NoLOs, especially when a huge amount of them is needed.
In this work, we develop a general scheme for calculating response functions using NoLOs. This method bears similarity to Wannier interpolation and allows calculation of derivatives of arbitrary orders of band energies and Bloch wave functions. Therefore, it is suitable for calculation of various linear and nonlinear response functions. As an example, we calculate dielectric constant and shift current conductivity of monolayer WS 2 , as a proof of calculations of linear response and nonlinear response functions. The obtained results are in excellent agreement with previous works. This validates our computational formalism and its implementations. Finally, we discuss the performance and symmetry properties of this method.

Method
and its derivative A k  lie at the heart of response functions.
However, due to the random phase generated by diagonalization, direct derivative calculation by finite difference is not possible. Wannier interpolation scheme solves this problem by fixing the phase by choosing a definite gauge. Here we extend this method to NoLOs. We label localized orbitals as Rnñ | , where n is an index of the orbital inside the unit cell labeled by a lattice vector R. n runs from 1 to N, where N is the number of localized orbitals in one unit cell. Bloch summations of these localized orbitals constitute a complete basis of Hamiltonian at a k point, Under this basis, the Hamiltonian matrix is where V and E can be obtained by solving a generalized eigenvalue problem: which is just the Schrödinger equation expressed in the nonorthogonal basis k n L y ñ | ( ) . v k n is a column vector and has N independent solutions constituting the columns of matrix V. Notice that V is not a unitary matrix. Assuming the localized orbitals are linearly independent, S is Hermitian and positive definite. This type of generalized eigenvalue problem is well behaved and the eigenvector can be normalized as where r is the position operator, k ¶ º ¶ a a and α is the Cartesian indices. The position matrix depends on the choice of origin, as it should, since diagonal elements of Berry connection A depend on the choice of origin. However, since only off-diagonal elements of A are needed in the expression of dielectric constant and shift current, the results are independent of the choice of origin. For NoLOs, the position operator matrix satisfies is a Hermitian matrix (with fixed R) in the usual Wannier interpolation.
One critical issue is that the calculation of V S V ¶ a † still suffers from arbitrary phases from diagonalization. Fortunately, since the three matrices are known from ab initio calculations, arbitrary derivatives of H, S and A (L) are known, which makes it possible to calculate V S V ¶ a † [11].
Multiplying equation (12) with v n † on the left, we have However, due to the intrinsic phase ambiguity of v n , it is necessary to introduce an extra gauge fixing condition is actually positive around Γ. Differentiating both equations (8) and We have assumed in the above equations that n is not degenerate. Degenerate eigenvalues can be treated in principle [12] but are not important in our calculations [7]. In the calculations, the procedure outlined above is done independently at every k point, owing to the property that response functions can be expressed in a way that a global smooth gauge is not needed [7].
Every ingredient needed for linear response are obtained at this point. However, we still need A k  for nonlinear responses. Differentiating equation (9) again, we have The corresponding result for orthonormal Wannier functions [7] is reproduced by choosing S=I. The However, it is not difficult to imagine how complex the final expression would become. Therefore, we introduce in the appendix an iteration procedure to calculate derivatives of v n of arbitrary orders for nondegenerate E n .
Another way to find Berry connection is to do orthogonalization by S −1/2 in equation (7) and use the usual Wannier interpolation developed previously, but derivatives of S −1/2 will be needed to calculate derivatives of the transformed Hamiltonian. This orthogonalization scheme will be developed and discussed in detail in section 4.
Finally, we remark that since the full spectrum is needed at every k point, this method includes a non-selfconsistent calculation at every k point, in contrast to the usual Wannier interpolation scheme. Henceforth, we will refer to the scheme developed here as NoLO-based method and the usual Wannier interpolation as MLWFbased interpolation.
3. Shift current of monolayer WS 2

Background and computation details
Shift current [13][14][15][16] is a second-order bulk photovoltaic effect arising from the difference of real space positions of Bloch electrons between valence band and conduction band where shift current conductivity s abb is given by [13][14][15][16] k . These quantities can be calculated using the method described in section 2. Then an numerical integration would produce the result of s abb . Since a δ function is present in the expression of s abb , a very fine k mesh is needed to achieve convergence.
Since monolayer WS 2 has the point group D 6h , there's only one independent shift current conductivity component σ yyy [7,17] 7 . Following the convention of [7], we choose a two dimensional (2D) definition of current. Two DFT packages are utilized in the calculations: the full-potential, all-electron FHI-AIMS package [18] and the pseudopotential-based OPENMX [19,20] package. For comparison, MLWF-based calculations are also carried out with VASP [21,22] and WANNIER90 [9] packages. In all the calculations, slab model is used to characterize monolayer WS 2 with a vacuum layer thicker than 15 Å. A k grid 12×12×1 is used to sample the Brillouin Zone in self consistent calculations and a much finer k grid 400×400×1 is used to perform the numerical integration in the expression of shift current conductivity. Interaction effects are captured in Perdew-Burke-Ernzerhof exchange-correlation functional [23]. To expand the Kohn-Sham wave functions, we choose the so-called 'tight' numerical settings in FHI-AIMS calculations, pseudo atomic basis 'W7.0-s3p2d2f1' and 'S7.0-s3p3d2f1' in OPENMX calculations and plane wave basis with an energy cut of 258.689 eV in VASP calculations. Spin-orbit interaction is not included and δ function is simulated using the following numerical approximation where the broadening factor ò is chosen to be 0.1 eV.

Results
The shift current conductivity of monolayer WS 2 is presented in figure 1(a), compared with MLWF-based method introduced in [7] with the same parameters. Despite using different software packages and DFT schemes, the calculated shift current conductivity curves are almost the same. Optical absorptions, represented by imaginary part of dielectric constant calculated by both MLMF-based interpolation and the method developed here, are plotted in figure 1(b). Both dielectric constant and shift current conductivity have two peaks at 2.75 eV and 3.05 eV respectively. While for dielectric constant, the peak at 3.05 eV is higher than that at 2.75 eV, an opposite feature is observed for shift current conductivity. This difference should be attributed to the difference of shift vectors. The contribution of these two peaks to dielectric constant can be decomposed to individual bands and is shown in figure 1(c). Sizes of red and blue dots represent contributions to the 2.75 eV peak and 3.05 eV absorption peak respectively. It is obvious that both peaks are mainly contributed by the highest valence band and lowest four conduction bands around Γ and K points.

Discussion of the method
MLWFs are known to break symmetry slightly, which is revealed by small avoid crossings in the interpolated band structure where they should have been direct crossings. This behavior results in a small but nonvanishing value for symmetry forbidden components of shift current conductivity even for well convergent MLWFs. This problem does not arise for NoLO-based method since symmetry is enforced in ab initio calculations by only sampling the irreducible Brillouin Zone in the calculation. Figure 2 . Therefore, NoLO-based method preserves symmetry properties quite well. Compared to MLWF-based interpolation, NoLO-based calculation is computationally heavier, since ab initio packages, especially all-electron full-potential packages, need to use many NoLOs to span the Hilbert  space, while MLWFs are usually only constructed for bands near the Fermi surface. Therefore, it is necessary to test the scaling behavior of this method with respect to the number of NoLOs. This scaling behavior is presented in figure 2(b). It is observed that computation time roughly scales as N O 2 ( ), where N is the number of NoLOs. This is quite unexpected since diagonalization scales as O(N 3 ). A closer analysis of the performance reveals that the computation is dominated by Fourier transformations in calculating H, S, A (L) and their derivatives. Therefore, for common bulk materials, we can safely assume the time complexity is N O 2 ( ). There is another source of performance difference between the MLWF-based method and the NoLO-based method. Due to the orthogonality of MLWFs, the Fourier transformation of S is trivial in the MLWF-based method and thus saves one third of the time compared to NoLO-based method, even if they use the same number of orbitals.

Alternative scheme by orthogonalization
The Hamiltonian in the basis of k Notice that the principle square root of k S ( ) is well defined since k S ( ) is positive definite. In this way, MWLFbased interpolation method can be directly applied to calculate nonlinear optical responses.
However, one key ingredient of MWLF-based interpolation is the derivative of k H O ( ) ( ) with respect to k. By equation (23), this issue reduces to the calculation of derivatives of arbitrary order of k S 1 2 -( ) with respect to k. This can be done as follows. By differentiating S −1/2 S −1/2 =S −1 , we obtain Orthogonalization method presented in this section has the advantage of being directly connected to previous Wannier interpolation methods. However, the pervasive usage of matrix inversion makes this method slightly slower and less numerical stable than the method presented in section 2.

Conclusion
Real space localized orbital-based ab initio packages have the potential of achieving N O( )computational resource scaling with respect to number of atoms. In addition, vacuum can be treated without extra computational cost in these packages, making them competitive tools in research for low dimensional materials. Here we demonstrate these packages can be more powerful by extending Wannier interpolation to NoLOs.
Although computationally heavier compared to MLWF-based interpolation, NoLO-based scheme developed in this work avoids human intervention and can be used in high throughput material discovery. The correctness of this scheme is proved by calculating shift current conductivity of monolayer WS 2 . This NoLObased method is quite general and can be used to calculate different kinds of linear responses and nonlinear responses in different research fields [26], including anomalous Hall effect (see also [27] 8 ) for a optical matrix based approach, nonlinear Hall effect, orbital magnetization and second harmonic generation.