An efficient method for hybrid density functional calculation with spin-orbit coupling
In first-principles calculations, hybrid functional is often used to improve accuracy from local exchange correlation functionals. A drawback is that evaluating the hybrid functional needs significantly more computing effort. When spin-orbit coupling (SOC) is taken into account, the non-collinear spin structure increases computing effort by at least eight times. As a result, hybrid functional calculations with SOC are intractable in most cases. In this paper, we present an approximate solution to this problem by developing an efficient method based on a mixed linear combination of atomic orbital (LCAO) scheme. We demonstrate the power of this method using several examples and we show that the results compare very well with those of direct hybrid functional calculations with SOC, yet the method only requires a computing effort similar to that without SOC. The presented technique provides a good balance between computing efficiency and accuracy, and it can be extended to magnetic materials.
pacs:63.20.dk, 73.22.-f, 71.70.Ej
Density functional theory (DFT) is a powerful method for predicting properties of materials such as crystal structures, electronic bands, phonon dispersions and other physical quantities. Practically, using appropriate exchange correlation (XC) functional is very important for accuracy, especially for predicting band gaps of materials. It is well known that the local density approximation (LDA) and general gradient approximation (GGA) XC functionals tend to severely underestimate band gaps Kohn-gap (). Consequently, hybrid functionals (HF) such as PBE0 PBE01 (); PBE02 (); PBE03 (); PBE04 (), HSE03 and HSE06 HSE06 (); HSE2 (); HSE3 (); HSE4 () were proposed and they often predict very good band gap values comparable to experiments. A drawback of HF is that it needs very significant computing resources, generally several orders of magnitude more compared to that of LDA or GGA.
In recent years, materials with strong spin-orbit coupling (SOC) have attracted great attention, including topological insulators Hasan10 (); Qi11 () BiSe Bi2Se3 (), silicene Liu1 (); Liu2 (), germanene Liu1 (); Liu2 (), stanene Liu1 (); Liu2 (); Xu2013 (), BiH Song2014 (); Liu2014 (), ZrTe ZrTe-Weng (), BiBr Zhou2014 (), ZrSiO WHM-family (), photoelectric materials PbI PbI2 () and BiOCl BiOCl (), two-dimensional group-VI transition metal dichalcogenides MoS, MoSe, WS, and WSe LiuRev2015 (), III-V direct band-gap semiconductors with heavy elements GaSb and InSb k35 (), etc. DFT calculations including SOC involve non-collinear spin which requires at least eight times more computing time as compared to that without SOC, due to the scaling for solving the Kohn-Sham DFT equations (KS-DFT). Since many of these SOC materials are semiconductors, HF calculations are desired to more accurately predict their band gaps and electronic structures. Unfortunately, HF+SOC calculations are numerically intractable thus rarely used - unless the unit cell is extremely small, due to the huge computational demand. It is the purpose of this paper to report a practical solution to this problem.
In particular, we propose an efficient approximate technique for HF+SOC calculations based on a mixed linear combination of atomic orbital (LCAO) scheme. The mixed LCAO Hamiltonian is constructed by two parts: an SOC-free part whose parameters are obtained from HF calculations without SOC, and an SOC part whose parameters are obtained from GGA+SOC calculations (DFT at the GGA level with SOC). Applying this approach to several non-magnetic materials, the results are demonstrated to be very close to those of direct HF+SOC calculation and much more accurate than the GGA+SOC calculation. Importantly, the required computing time of the mixed LCAO technique is comparable to that of HF calculation without SOC.
In the rest of the work, the DFT calculations are performed using the projector augmented wave method implemented in VASP Kresse1996 (). The Perdew-Burke-Ernzerhof (PBE) parametrization of GGA functional Perdew1996 (); KJ99 () and Heyd-Scuseria-Ernzerh hybrid functional (HSE06) HSE06 (); HSE2 (); HSE3 (); HSE4 () are used in the DFT calculations, and the VASP2WANNIER90 interface wannier2008 (); wannier2012 (); vasp2wannier90 () is used to obtain the LCAO parameters from the DFT results. Since numerical calculations are for the purpose of demonstrating the mixed LCAO technique, structure optimization is omitted.
Ii The method
WANNIER90 wannier2008 (); wannier2012 () is used to construct LCAO or Wannier-bases Hamiltonian from DFT calculations, and the resulting LCAO Hamiltonian can reproduce the original energy dispersion very well. We start by constructing an LCAO Hamiltonian to treat HF+SOC using DFT calculations.
For a given system, the required computing effort is most demanding for HF+SOC, followed by HF without SOC and next followed by GGA+SOC. Clearly and as explained in the Introduction, if HF+SOC were computationally affordable in general, the work of this paper would not be necessary. That is not the case. In the following we utilize HF without SOC and GGA+SOC to construct a mixed LCAO Hamiltonian which we show to be a very good approximation to . In particular, has two terms, which is obtained from HF without SOC, and which is obtained from GGA+SOC,
Clearly, constructing only consumes a time that is comparable to HF without SOC, thus much more efficient than that of a full HF+SOC calculation. The fact that compare very well with direct HF+SOC calculations (see below), suggests that the mixed LCAO scheme provides a viable approximation for the complicated HF+SOC analysis.
On the technical side, while can be constructed directly from DFT calculation of HF without SOC, is obtained from DFT of GGA+SOC involving a procedure for separating out the SOC contributions. The latter procedure and an associated technical detail are discussed in the following two subsections.
ii.1 Separating out the SOC contribution
The Hamiltonian of SOC systems can be divided into a non-SOC term plus the SOC term :
In LCAO representation, involves on-site energy and hopping integral between different atomic orbitals, and comes from SOC effects.
In the spin-up and spin-down bases and , the non-SOC term can be written as a diagonal matrix:
For simplicity, we consider non-magnetic systems in the rest of this work, but extension to magnetic system can be readily made without fundamental difficulty. For non-magnetic materials, .
For the SOC term , its original operator form is:
where is the reduced Planck constant, is the bare mass of electron, is the velocity of light, is the potential energy, the momentum, and the vector of Pauli matrices representing the spin degrees of freedom. For clarity we define a constant and a vector operator . can then be rewritten in the following matrix form:
in which and .
Then, from Eq. 6, we can separate the total Hamiltonian to obtain and as the following:
Hence, after obtaining the LCAO Hamiltonian from the corresponding DFT calculation, its SOC part can be separated out using Eq.8.
ii.2 Mixing the Hamiltonian
With the obtained non-SOC part and SOC part , the mixed LCAO Hamiltonian that approximates HF+SOC is determined by Eq.1. Hereinafter we use the HSE functional for HF, and PBE functional for GGA. Then Eq.1 becomes
The “mixing” procedure appears to be a simple addition. However it should be noted that only when and are constructed under the same bases can they be added directly. We achieve this by constructing and in the same bases , and details are presented in the appendix A. This way, we finally constructed the mixed Hamiltonian to treat HSE+SOC.
Iii Results, analysis and discussion
Having constructed to efficiently treat HSE+SOC, in this section we demonstrate its accuracy using several material systems. Predicting band gap is important, which is one of the reasons to use HSE in the first placeHSE06 (); HSE2 (); HSE3 (); HSE4 (). We calculated band gaps for eight semiconductor materials having heavy elements thus large SOC, including two-dimensional (2D) mono-layers of PbI, WSe, BiH, BiI and BiBr; 3D crystals BiOCl, GaSb, and InSb k35 (). The SOC effect is important for these materials, especially for their band gaps.
For the eight materials, we performed (very time-consuming) direct HSE+SOC calculations and using the results, we constructed an LCAO Hamiltonian : this would not be possible without the full direct HSE+SOC calculation. Then, we constructed following the procedure in the last section which does not require full HSE+SOC calculation. The three sets of results are compared: direct numerical data from full HSE+SOC calculations and from , as well as from .
First, band structures of the eight compounds calculated by our method via of Eq.9 is plotted in Fig.1, together with those from the direct HSE+SOC calculation and its fitting . The full HSE+SOC data are presented in black circles and the bands from are in thick red lines: these are used as benchmarks to compare to our results by which are presented in thin blue lines. We can see that the bands calculated by our method via of Eq.9 (thin blue lines) are qualitatively consistent to the benchmark results for all cases. In particular, band dispersions by and the benchmark agree well for PbI, WSe, GaSb, and InSb; and the agreement is somewhat reduced for BiOCl, BiH, BiI, and BiBr. In the latter cases, although the heavy element Bi gives rise to some quantitative difference, there is no qualitative discrepancy for the full range of the Brillouin zone.
Second, quantitatively we compare the band gaps of and the benchmark in Table 1. The band gaps obtained by are close to those of the benchmark for five of the eight compounds PbI, WSe, BiOCl, GaSb, and InSb [see the column labeled by column in Table 1]. But for three compounds BiH, BiI and BiBr - especially BiI which has a very small band gap, the discrepancy is large. As a supplement, we also show the calculated band gaps by PBE+SOC in Table 1: they are not only quantitatively quite different from the benchmark results, two of them are even qualitatively wrong (GaSb, InSb).
As for the three compounds with band gaps of showing large discrepancy to the benchmark, BiH, BiI and BiBr, they are all relevant to topological insulators which exhibit band inversion near the Fermi level when SOC is considered and have normal bands with SOC not considered. Concretely, 2D monolayers of BiH and BiBr are topological insulatorsZhou2014 (). As for BiI monolayer, although it is not a topological insulator, it lies near the transition point between normal and topological insulator which makes it also sensitive to SOC. The impact of topological properties on the accuracy of our method needs to be analyzed.
Recall that we used to approximate which can be written as
where and were obtained directly from using Eq.7 and Eq.8, respectively. Therefore, any discrepancy between and has two sources: (i) the discrepancy between the non-SOC part (obtained from HSE without SOC) and ; (ii) the discrepancy between the SOC part and . We analyze these terms in the next two subsections.
iii.1 HSE without SOC
In this subsection we analyze the source of discrepancy in the non-SOC part (obtained from HSE without SOC) and . Since DFT is based on ground state electronic density, any discrepancy should be due to the difference between densities calculated by the two approaches HSE-WF (). To illustrate this, we proceed by making use of a model Hamiltonian as follows:
where the first and second terms are the non-SOC part , the third to fifth terms — denoted as , , and respectively, constitute the SOC part . In Eq.11, and are respectively annihilation and creation operators for the conduction ( valence) bandAdd1 () eigen-state () of whose eigenenergy is ().
The Hamiltonian in Eq.11 can be analyzed by perturbation theory. We take the non-SOC as the unperturbed Hamiltonian and the SOC term as the perturbation. Note that can be divided into two types according to their SOC effects: and are the “type-I” terms, is the “type-II” term. The type-I term couples conduction bands with other conduction bands as well as valence bands with other valence bands. We distinguish two situations. (i) If the unperturbed band gap is large relative to the type-I SOC effect, the type-I term only splits bands and no band crossing occurs. This is the “weak type-I SOC” which does not make any significant difference between the unperturbed and perturbed charge densities (see Appendix B for details). (ii) If the unperturbed band gap is small relative to the type-I SOC effect, some split bands near the band gap will cross the Fermi level so that the charge density is altered. This is the “strong type-I SOC” [cf. Fig. 2(b)].
As for the type-II SOC , it mixes valence bands with conduction bands, for instance the first-order perturbed valence band state is
in which and is a normalization constant. As a result, the type-II SOC term alters charge density via the sum of all occupied states, i.e. all valence band states in the presence of band gap [cf. Fig. 2(c)].
|valence band||conduction band|
|WSe Liu2013 ()||W-d&d||W-d|
|BiH Liu2014 ()||Bi-p||Bi-p|
|BiI Zhou2014 ()||Bi-p&Bi-p||Bi-p&Bi-p|
Let us first understand why for the five normal band insulators PbI, BiOCl, WSe, GaSb and InSb, agrees with the benchmark very well. For these compounds, the properties of their lowest conduction band differ significantly from those of highest valence band, especially in orbital compositions and positions (cf. Table.2). This makes the type-II SOC hopping composed mostly of the SOC interactions between different atoms, and hence are very small. WSe is an exception because its band-edge orbitals locate at the same atom W. However the W-d&d orbitals have opposite mirror symmetry compared to W-d orbital, which still makes small. Furthermore, the gaps of the normal band insulators are large compared to , i.e. . Consequently the perturbed valence band eigenfunction can at most slightly mix with the conduction band since (see Eq. 12), which means type-II SOC has almost no effect. This is shown in table 3 which illustrates band gaps of and differ by less than 2.3%. In addition, their band gaps are also large enough compared to the type-I SOC hopping and , hence this is the weak type-I SOC case which makes no change of charge density if type-II SOC is not considered. We conclude that for normal band insulators, both type-I and type-II SOC effects contribute very little change to the charge density from to . This is why that in DFT calculations of these compounds, of HSE without SOC () and HSE+SOC () are very close to each other since charge densities are very close.
Next, we analyze the three compounds where has significant discrepancy to the benchmark . As discussed above, these compounds have inverted bands (BiH, BiBr) or near-inverted bands (BiI), the orbital compositions of their conduction bands are similar to those of valence bands around the Fermi energy (cf. Table.2). Hence is not small in general and can be close to or even larger than . According to Eq.12, this makes not small and the type-II SOC significantly perturb . In addition, the large SOC strength of Bi can lead to a strong type-I SOC that make conduction and valence bands cross the Fermi energy [see Fig. 2(b)]. Due to the strong type-I and type-II SOC effects, the charge density changes significantly from to , which makes and differ from each other.
Quantitatively, we compare the band gaps of and in Table. 4, from which we observe that the gaps of are very close to those of for PbI, BiOCl, WSe, GaSb and InSb, but not so close for BiI and BiBr. As a comparison, we also show the gaps of and in Table 4. It should be pointed out that because the PBE gaps without SOC are smaller than HSE gaps, i.e. are smaller in PBE, will mix more if and are assumed to be the same in HSE and PBE. Hence the gap differences of PBE are worse than those of HSE for GaSb, BiI and BiBr, especially for BiBr, as shown in table 4.
iii.2 PBE with SOC
Having understood the discrepancy in the non-SOC part (obtained from HSE without SOC) and , in this subsection we analyze the discrepancy between the SOC part and . For different XC functionals such as PBE and HSE, there are two major reasons that make different in PBE and HSE.
First, according to Eq. 4, different potential energy makes different SOC parameters. For HSE, its exchange potential is composed of a short-ranged part and a long-ranged part , where is produced by mixing the non-local Fock potential (i.e. the exact exchange potential) and the PBE exchange potential in short range, while is solely the PBE exchange potential in long range. Since the unscreened Fock potential is generally larger than the PBE exchange potential, the resulting with replaced partially by the unscreened , should also be larger than in general. This is demonstrated by the SOC hopping parameters of HSE and PBE in Fig. 3, in which the HSE ones are consistently larger than the PBE ones. Taking BiH as an example: it has a hexagonal structure like graphane and has Dirac cones at and points in the Brillouin zone without SOC. With SOC, a topological band gap opens at point and this gap depends only on the strength of SOC. Hence, the larger gap of HSE+SOC than that of our method shown in Fig1(f) means that the SOC strength of HSE+SOC is larger than that of our method, i.e. is larger than . For other topological insulators similar to BiH, i.e. the ones which have Dirac points or node lines before including SOC, such as ZrTeZrTe-Weng () and ZrSiO familiesWHM-family () which have band inversion before considering SOC just as BiH does, they share similar source of discrepancy to BiH in our methodAdd2 (). Second, if and of PBE and HSE are assumed to be the same, using the perturbation theory, for near-inverse band insulator (BiI) or inverse band topological insulator (BiBr), different gaps between PBE and HSE make mix at different ratios with , resulting in different charge densities which gives further differences of and between PBE and HSE.
To quantitatively understand the difference between and , we compare the gaps of two Hamiltonians: one is , and the other is with its SOC part replaced by , i.e. in Table 1. We can see from Table 1 (especially the last column) that the gap differences induced by the difference between and are relatively small for band insulators PbI, BiOCl, WSe, GaSb and InSb, but large for BiH, BiI and BiBr.
According to these analyses, for normal band insulators, the differences between and are as small as the differences between and . This is why our method works so wy our method works so well for the five compounds PbI, WSe, BiOCl, GaSb and InSb (see Table 1). But for the three near-inverse and inverse band insulators BiH, BiI and BiBr, the differences between and are much larger [cf. the column in table 4 and the column in table 1]. Furthermore, by comparing and in Table 1, we conclude that the error of our method is dominated by the difference between and , and this error are large for the near-inverse and inverse band insulators. Note, however, although the relative deviation of our with respect to is large for BiI and BiBr, their absolute deviations are actually not large, as illustrated by the not-so-large differences of the SOC hopping parameters between PBE and HSE shown in Fig.3.
We therefore conclude that our method is a very good approximation for normal band insulators and, for near-inverse or inverse band insulators which have very small band gaps, our method is still reasonable in that it can provide qualitatively correct results.
The purpose of this work is to develop a reasonably accurate, qualitatively correct and computationally efficient method to perform HSE+SOC calculations within DFT. We have so far demonstrated the accuracy of our method and understood the source of discrepancy when dealing with near-inverse and inverse band insulators.
Concerning computational efficiency: our method for HSE+SOC calculation takes essentially the same time as an HSE calculation without SOC. Taking WSe for example, in our calculations, one PBE+SOC electronic step takes 1.3 minute using 32 CPU cores, one HSE (without SOC) electronic step takes 5.6 minute using 128 CPU cores, and one full HSE+SOC electronic step takes 50.5 minute using 128 CPU cores. Therefore our technique is nearly an order of magnitude faster than the full HSE+SOC calculation. We checked that for all the cases we investigated, our method is faster by several times to more than an order of magnitude that the full direct HSE+SOC approach.
So far we analyzed non-magnetic compounds, but the method can be easily extended to magnetic materials for which and are not equal anymore but can be obtained in an additional spin-colinear DFT calculation. Namely, our method can be extended to magnetic materials by performing a spin-colinear PBE calculation and a PBE+SOC calculation to extract ; then performing a spin-colinear HSE calculation to obtain . The results are added together to obtain for the magnetic material.
The method developed in this work is not only suitable for LCAO, but also useful for accelerating HSE+SOC DFT calculations. Usually, two initializations are used to save computing time during HSE+SOC calculations: (a) using charge density and wave function from a PBE+SOC calculation as initialization, (b) using charge density and wave function from a HSE calculation without SOC as initialization. However, both will actually not accelerate calculation significantly. This is because for (a), of PBE+SOC differ quite a lot from that of HSE+SOC; and for (b), it lacks . Then, naturally, our method provides a better starting charge density and wave functions for full HSE+SOC DFT calculation because is closer to than or .
In summary, we have developed an efficient mixed LCAO technique to perform HSE+SOC DFT calculations. The LCAO Hamiltonian is obtained by mixing a non-SOC part and an SOC part. The non-SOC part is constructed by SOC-free HSE, and the SOC part by PBE+SOC. As a result, the mixed LCAO technique requires a computing time comparable to that of a non-SOC HSE calculation, thus saving about one order magnitude in computing time compared to a full direct HSE+SOC calculation.
Applying the method to eight non-magnetic compounds demonstrates that the mixed LCAO Hamiltonian can well approximate that of the full HSE+SOC. In particular, the method works very well for normal band insulators, and it is also reasonable to give qualitatively correct results for near-inverse and inverse band insulators having very small band gaps. We find that the errors in our method came from the difference between and more than the difference between and in most cases. Our method can be easily extended to other hybrid functionals and magnetic materials, and it can also be used to provide good initial conditions for full direct HSE+SOC calculation.
The authors thank Qing Shi of Mcgill University for the data of GaSb and InSb. The work is supported by the MOST Project of China (Grant No. 2014CB920903), the NSF of China (Grant Nos. 11734003, 11574029), the National Key R&D Program of China (Grant No. 2016YFA0300600) (MW and YY); by the National Key R&D Program of China with Grant No. 2017YFB0701600 and NSFC with Grant No. 11304014(GBL); and Natural Science and Engineering Research Council (NSERC) of Canada (HG). HG thanks Compute Canada for computational facilities where part of this work was carried out.
Appendix A LCAO representation and interpolation
where is localized trial orbitals serving as the initial guess of the Wannier functions, and is the number of bands considered or the dimension of the Bloch manifold at . Thus obtained are only determined by the local orbitals , for can be Fourier transformed to in which is lattice vector. To show this, first apply the Bloch theorem to eq. 13 to get
and then do Fourier transformation to get
where is the number of unit cells (also the number of points). In the above equation, we have used the completeness relation
Under the bases of , the Hamiltonian matrix is with matrix elements
where the overlap matrix is defined by
To orthonormalize the bases, we construct with properties as follows
in which is a Hermitian matrix with the property and hence can be denoted as in form. The existence of is guaranteed by the Hermiticity of . As a result, the Hamiltonian matrix under the orthonormalized bases is
which can be constructed from the data and generated by VASP and VASP2WANNIER90 respectively.
Because are determined only by (eq. 18), is determined only by (eq. 21), is determined only by , and are determined only by and (eq. 22), it can be concluded that the final bases are only determined by the initial local orbitals and irrelevant to the Bloch states . This is crucial for us to add directly and constructed from independent DFT calculations, because and have the same bases as long as the same initial local orbitals are used.
The obtained ) from eq. 23 is defined at only a finite number () of points. If the Hamiltonian at an arbitrary wave vector different from the points is required, interpolation has to be done. The interpolation can be achieved through two steps. First, do a Fourier transformation for () to get the Hamiltonian element in real space
in which is the Fourier transformation of
and is an orthonormalized local orbital with property . Then, construct the Hamiltonian at arbitrary wave vector using the -indepnedent quantities as follows :
We call this procedure LCAO interpolation.
Appendix B The weak type-I SOC effect
To interpret the effect of the type-I SOC, we consider the Hamiltonian eq. 11 in the absence of type-II SOC hoppings
in which represent valence bands, represent conduction bands, and are the number of valence and conduction bands respectively, and is the total number of bands considered. Note that the spin index is incorporated into the band index and here for simplicity. Choosing the eigen states of as bases, is a diagonal matrix with
and is a block-diagonal matrix with and being and matrices respectively. Hence, the charge density of is
where is the eigen state of and the second equality is due to the fact that is related to by a unitary transformation with and (For simplicity we use the usual normalization convention here). This is guaranteed by the absence of type-II SOC hoppings .
If the type-I SOC effect is weak, will not make conduction or valence bands cross the Fermi energy like Fig. 2(b). This makes the charge density of