Closest Wannier functions to a given set of localized orbitals
The original has been published
in Phys. Rev. B 110, 125115 (2024).
Abstract
A non-iterative method is presented to calculate the closest Wannier functions (CWFs) to a given set of localized guiding functions, such as atomic orbitals, hybrid atomic orbitals, and molecular orbitals, based on minimization of a distance measure function. It is shown that the minimization is directly achieved by a polar decomposition of a projection matrix via singular value decomposition, making iterative calculations and complications arising from the choice of the gauge irrelevant. The disentanglement of bands is inherently addressed by introducing a smoothly varying window function and a greater number of Bloch functions, even for isolated bands. In addition to atomic and hybrid atomic orbitals, we introduce embedded molecular orbitals in molecules and bulks as the guiding functions, and demonstrate that the Wannier interpolated bands accurately reproduce the targeted conventional bands of a wide variety of systems including Si, Cu, the TTF-TCNQ molecular crystal, and a topological insulator of BiSe. We further show the usefulness of the proposed method for calculating effective atomic charges. These numerical results not only establish our proposed method as an efficient alternative for calculating WFs, but also suggest that the concept of CWF can serve as a foundation for developing novel methods to analyze electronic structures and calculate physical properties.
I INTRODUCTION
Wannier functions (WFs) [1, 2, 3] play a central role in analyzing electronic structures of real materials and advancing electronic structure methods along with density functional theory (DFT) [4] and the other electronic structure theories [5, 6]. A widely adopted method for calculating WFs involves maximizing their localization, which can be reformulated as minimizing the spread function characterizing the variance of the WFs in real space [7, 8]. The concept of maximal localization leads to an elegant formulation, resulting in maximally localized Wannier functions (MLWFs), and the following compact representation of the Hamiltonian allows us to analyze the electronic structures of real materials using a localized orbital picture. However, it is also recognized that the minimization of the spread function often encounters local minima, particularly in large-scale systems with complicated electronic structures [3]. To overcome this challenge, methods aimed at automated high-throughput Wannierisation have been developed [9, 10, 11]. In addition, other methods for generating (nonorthogonal) WFs have been proposed based on projection methods [12, 13, 14], which produce compact atomic like orbitals. The orthogonalization of such nonorthogonal orbitals obtained by the projection can be achieved using the Löwdin orthogonalization procedure [15, 16]. These orthogonalized orbitals are utilized as an initial guess of the WFs in the minimization of the spread function to obtain MLWFs [7, 17]. Apart from the role of the initial guess, it is worth noting that the Löwdin orthogonalized orbitals possess an important variational property, which in particular minimizes the sum of the squared distances between the orthogonal and original nonorthogonal orbitals in the Hilbert space [18]. As long as the original nonorthogonal orbitals are localized, the localization of the Löwdin orthogonalized orbitals in real space is guaranteed in a sense of minimization of the distance. In quantum chemistry, this property is leveraged to develop methods for generating localized orthogonal orbitals. These methods, which include the natural atomic and bond orbital method [19] and the intrinsic atomic orbital method (IAO) [20], involve obtaining orthogonalized orbitals through Löwdin type orthogonalization, giving heavy weight to occupied states via the density matrix. On the other hand, in solid state physics, the variational property inherent in the Löwdin orthogonalization has not yet been fully explored. This is because the Löwdin orthogonalization requires the calculation of for the overlap matrix , which causes numerical instability for nearly ill-conditioned matrices. Owing to this difficulty, advancements in methodologies along these lines appear to be slow. It should be mentioned that avoiding the calculation of has already been proposed based on singular value decomposition (SVD) to calculate an initial guess of WFs for the maximal localization [8]. However, the relationship between the WFs obtained through this procedure and the variational principle has not been sufficiently discussed so far.
Our study aims to develop an efficient and robust method for generating the closest Wannier functions (CWFs) to a given set of localized guiding functions, based on a variational principle and polar decomposition of the projection matrix. We will demonstrate that exploiting the variational property leads to a versatile method that eliminates the need for iterative optimization and avoids complications related to the gauge choice. The projection method based on the polar decomposition has already been employed as a workaround [21, 22] with successful applications [23, 24]. However, to the best of our knowledge, no study has thoroughly discussed the mathematical and physical significance of this method based on the closest principle. Our study establishes the theoretical foundation by formulating the closest principle as a variational problem, and clarifies the relationship between the Löwdin orthogonalization and the polar decomposition. We will also introduce a smoothly varying window function that enables us to flexibly control targeted states to be Wannierized in a similar way to that proposed in Ref. [10]. The present study offers a more rigorous justification for these practices. With the theoretical advancement, the concept of CWF can serve as a foundation for developing novel methods to analyze electronic structures and calculate physical properties as exemplified in the analysis of effective atomic charges in Sec. IV.
The structure of this paper is as follows: In Sec. II, we present the theory for calculating CWFs. Sec. III provides a detailed discussion on the implementation of the method. A series of benchmark calculations are presented in Sec. IV. Finally, in Sec. V, we summarize the theory of CWFs and suggest the potential role of CWFs as a foundation for developing novel methods.
II THEORY
Let us start by introducing Bloch functions , which can be obtained by solving the Kohn-Sham (KS) equation [25] within the density functional theory (DFT) [4], normalized as
(1) |
where and are the k-vector and band index, respectively, and is the number of primitive cells in the Born-von Karman (BvK) boundary condition. We consider the projection of a localized guiding function onto the Bloch functions weighted with a window function as
(2) |
with
(3) |
where and are the translational lattice vector and index of the localized orbital, respectively, is the eigenvalue of the KS equation, and . In the summation of Eq. (2), the numbers of k-points and the Bloch functions to be included per k-point are and , respectively. Throughout this paper, we do not consider the spin dependency on the formulation for the sake of simplicity, but the generalization is straightforward. The localized functions , such as atomic orbitals, hybrid orbitals, or molecular orbitals (MOs), are assumed to localize in the unit cell specified by . When the window function is taken to be unity, Eq. (2) is nothing but an expansion of with Bloch functions from the identity relation . To introduce an expansion of by the Bloch functions in a subspace, we use the following window function :
(4) |
with the definition of and , where with the Boltzmann constant of and a temperature of . The window function of Eq. (4) is obtained by subtracting 1 from the sum of the two functions of . It is also possible to introduce two temperatures and . The usefulness of the two temperature scheme will be demonstrated for the case of GaAs in Sec. IV. The last term in Eq. (4) is a small constant, e.g., , which is introduced to avoid the ill-conditioning of the matrix consisting of by Eq. (3) as will be discussed later. As shown in Fig. 1, the parameters of and () determine the energy range where the Bloch states are included in the expansion of Eq. (2) with a large weight, and the temperature provides the degree of smearing around and in . To control the degree of smearing, we introduce the temperature of so that one can intuitively understand the degree of smearing. It should be noted that the choice of the window function by Eq. (4) is not unique, and the other choices should give almost equivalent results as long as they are chosen properly. One may also consider that the window function has no special advantages over the outer and inner window scheme [8] and merely functions in a similar manner. However, the window function defined by Eq. (4) exhibits superior performance in the disentanglement of bands and calculation of effective charges. The advantages will be clearly demonstrated in Secs. IV A and IV B, respectively.
The function defined by Eq. (2) must be similar to the localized function , but they are generally non-orthogonal to each other. We now consider generating a set of closest Wannier functions (CWFs) to a set of localized functions in the sense that the sum of the squared distance between and in the Hilbert space is minimized. As well as Eq. (2), noting that the CWFs can be expressed [3] as
(5) |
and defining the residual function :
(6) |
the distance measure (DM) function we minimize is given by
(7) | |||||
with
(8) |
where the elements of the matrices and are given by and , respectively. The second line of Eq. (7) is obtained by considering the orthonormality of Eq. (1), and of Eq. (8) is regarded as the squared Frobenius norm of [26]. The same Bloch functions as in Eq. (2) are included in the summation of Eq. (5). The number of CWFs to be generated is per unit cell, being equivalent to the number of the guiding functions per cell, and should be smaller than or equal to to guarantee the linear independency of the subspace spanned by the CWFs, resulting in size of for the matrices and . Assuming , where is of size , the CWFs are ensured to form a set of orthonormal functions. Under this constraint, we consider the optimization of the DM function . Furthermore, the matrix is regarded as a part of a unitary matrix of size , and will be referred to as a partial unitary matrix in subsequent discussions.
The minimum of the DM function of Eq. (7) is obtained by choosing as which is calculated by the polar decomposition of , where and are unitary and hermitian, respectively [27]. Let us prove the statement below. The polar decomposition of is obtained via SVD of as
(9) |
where we dropped the dependency on for simplicity of the notation, and hereafter we will denote the dependency if necessary. and are the left and right singular matrices in size of and , respectively, and is the singular value diagonal matrix in size of . Note that and . It is worth mentioning that and are partial unitary matrices, and hold () and (), and that is a full unitary matrix holding (). We evaluate of Eq. (8) for each with both the matrix obtained by the polar decomposition and an arbitrary partial unitary matrix , and calculate the difference as
(10) | |||||
with
(11) |
The second line of Eq. (10) is derived by using the polar decomposition of and , and and in the third line are diagonal elements of and , respectively. The upper bound of the diagonal elements of the hermitian matrix is found to be unity, since the matrix consists of the sum of the product of partial unitary matrices of () and (). Another proof for the upper bound is also given based on the Cauchy-Schwarz inequality in the appendix. By considering the upper bound of the diagonal elements of and , the third line of Eq. (10) leads to the following inequality:
(12) |
The equality of Eq. (12) holds if . If some of is zero, the corresponding in Eq. (10) can be arbitrarily chosen. So, we see from Eq. (11) that is not uniquely determined. If all the singular values of are positive, all the s should be unity when . The uniqueness of the solution can be confirmed as follows: Since the matrix consists of the products of the partial unitary matrices as discussed above, the case that all the diagonal elements are unity happens when
(13) |
which is derived by equating and multiplying and from the left and right of the equation, respectively. By noting again that and in Eq. (13) are the partial unitary matrices, it is found that Eq. (13) holds if and only if . Thus, we have proven the statement that the minimum of the DM function of Eq. (7) is uniquely determined when as long as all the singular values of are positive, since of Eq. (7) consists of the sum of over k. The uniqueness of itself is related to that of the SVD for the matrix of Eq. (3). If has positive singular values which are non-degenerate, the SVD is uniquely determined except for the non-unique phases in and . Even if there are degenerate singular values, the matrix is invariant, since the freedom of the unitary transformation for the degenerate left and right singular vectors is canceled out as with . Further, noting that the non-unique phases in and are canceled out when the matrix is computed as , we conclude that is uniquely determined if has the positive singular values. The violation from the positive definiteness of the singular values can be avoided by the small constant of in Eq. (4). On the other hand, if in Eq. (4) is set to be 0 and a narrow window, including fewer than eigenstates for a given , is used with a small , the matrix has positive singular values, where . In this case, a subspace of the dimension is arbitrary chosen to form and of the dimension of in addition to the subspace of the dimension, which breaks the symmetry of CWFs. The symmetry breaking is due to the non-uniqueness of the SVD when some of the singular values are zero. The situation can be avoided by introducing a small constant in the window function of Eq. (4). The simple treatment restores the symmetry of CWFs, allowing us to calculate symmetry preserving CWFs for a wide variety of choices in the energy range of the window function. It is also noted that the non-unique phase of the Bloch functions is canceled out via the polar decomposition of and Eq. (5). We have the arbitrariness of phase both in calculating the Bloch functions and singular vectors and , and unitary freedom for and when singular values are degenerate. However they are all canceled out as discussed above. Thus, we see that the proposed method is free from complications arising from the choice of gauge.
Using Eqs. (8) and (9), the minimum of the DM function at is given by
(14) |
where is a singular value of the matrix . From Eq. (14), we find that the mean squared distance between and is related to the deviation of from unity.
The proposed method based on the polar decomposition of is closely related to the Löwdin orthogonalization [15] as shown below. The Fourier transform of the overlap integrals for is given by
(15) |
which is obtained by noting that . By writing Eq. (15) as in a matrix form, and using Eq. (9), we have the following relation:
(16) |
Comparing Eq. (16) with , one obtains a relation . If is invertible, the matrix is given with Eq. (9) by
(17) |
The matrix obtained by Eq. (17) is exactly equivalent to that by the Löwdin orthogonalization [15], and the closest property of the Löwdin orthogonalized orbitals to a given set of orbitals was proven in Ref. [18]. On the other hand, the proposed method does not require the calculation of , and can be applied in a numerical stable manner even to the case that the matrix is nearly ill-conditioned. The definition of by Eq. (3) with the window function allows us to calculate CWFs in various respects, e.g., heavily weighting to occupied states and disentangling of localized states embedded in the wider bands as demonstrated in Sec. IV. In this sense, the method we propose can be regarded as a generalization of the Löwdin orthogonalization. We also notice that the projection method in Ref. [8] is equivalent to the proposed method in the sense that in Eq. (17) is calculated via the polar decomposition, and that the WFs obtained through the procedure [8] can be regarded as CWFs under the constraints of how the matrix is constructed.
The disentanglement of bands can be easily performed by properly selecting the window function of Eq. (3). Since the Bloch functions per k-point are included in Eqs. (2) and (5), it should be emphasized that the proposed method always disentangles states to generate Wannier functions. A couple of examples will be shown in Sec. IV, including valence and low-lying conduction bands for the diamond Si and narrow -bands embedded in the wider -band for the face centered cubic (FCC) Cu.
Once the -dependent are obtained by the polar decomposition, tight-binding (TB) parameters are calculated as expectation values of the KS Hamiltonian by summing contributions over and as
(18) | |||||
where is the matrix element of . Like with MLWFs, the TB parameters can be used for the Wannier interpolation in calculations of band structures and physical quantities.
The computational procedure to generate the CWFs is summarized as follows:
-
1.
Determining , , and in Eq. (4) by checking the band dispersion of a system of interest. The parameters should be chosen properly so that the targeted eigenstates can be included.
-
2.
Choosing a set of localized orbitals {}. Atomic orbitals, hybrid orbitals, and MOs might be possible choices. Depending on the symmetry of the targeted eigenstates, a proper set of localized orbitals needs to be chosen, e.g., -orbitals need to be employed in the generation of CWFs for the -bands in FCC Cu as shown later on.
-
3.
Calculation of by Eq. (3). The calculation of the overlap depends on the implementation of the KS method.
-
4.
Performing the SVD of . A proper numerical library might be used.
-
5.
Calculation of . Using the left and right singular matrices, and , of , is calculated as .
- 6.
It should be stressed that the minimization of can be efficiently performed in a non-iterative manner by the steps 1 to 6 using the polar decomposition via the SVD. The computational cost of each step is estimated to be , , , or for the step 3, 4, 5, or 6, respectively, where is the number of basis functions to expand the KS orbitals. If the localized basis functions are used, the computational cost of the step 3 becomes .
For isolated systems, the Brillouin zone sampling is limited to only the point. No modification of the proposed method is needed to calculate the CWFs.
The theoretical framework we have discussed so far is general for any implementation of the DFT-KS method. However, the choice of the the localized functions in Eq. (3) may depend on the implementation. So, we will discuss how the localized functions can be properly generated in Sec. III.
III IMPLEMENTATION
III.1 General
We have implemented the proposed method into the OpenMX DFT software package [28, 29, 30] which is based on norm-conserving pseudopotentials (PPs) [31, 32] and optimized pseudo-atomic orbitals (PAOs) [33, 34] as basis set. All the benchmark calculations were performed with a computational condition of a production level. The basis functions used are listed in Table 1. In the abbreviation of basis functions such as H7.0-s2p2d1, H stands for the atomic symbol, 7.0 the cutoff radius (Bohr) in the generation by the confinement scheme [33, 34], and s2p2d1 means the employment of two, two, and one optimized radial functions for the -, -, and -orbitals, respectively. The radial functions were optimized by a variational optimization method [33]. These basis functions we used can be regarded as at least double zeta plus double polarization basis sets if we follow the terminology of Gaussian or Slater-type basis functions. Valence states included in the PPs are listed in Table 1. All the PPs and PAOs we used in the study were taken from the database (2019) in the OpenMX website [28], which were benchmarked by the delta gauge method [35]. Real space grid techniques were used for the numerical integrations and the solution of the Poisson equation using fast Fourier transform (FFT) with the energy cutoff of 250 to 500 Ryd [29]. We used a generalized gradient approximation (GGA) proposed by Perdew, Burke, and Ernzerhof to the exchange-correlation functional [36]. An electronic temperature of 300 K was used to count the number of electrons by the Fermi-Dirac function for all the systems we considered. For all the calculations, was used in the window function of Eq. (4), and was set to be equivalent to the number of basis functions in both the summation of Eqs. (2) and (5).
Element | Basis functions | Valence states in PP |
---|---|---|
H | H7.0-s2p2d1 | |
C | C6.0-s3p2d2 | , |
N | N6.0-s3p2d2 | , |
Na | Na9.0-s3p3d2f2 | , , |
Si | Si7.0-s2p2d2f1 | , |
Cl | Cl9.0-s3p3d2f2 | , |
S | S7.0-s3p2d2f1 | , |
Cu | Cu6.0-s3p3d3f1 | , , , |
Ga | Ga7.0-s3p2d2f1 | , , |
As | As7.0-s3p2d2f1 | , , |
Se | Se7.0-s3p2d2 | , |
Bi | Bi8.0-s3p2d2f1 | , , |
III.2 Choice of Guiding Functions
Depending on the targeted eigenstates in the window function of Eq. (4), one can choose either atomic orbitals, hybrid atomic orbitals, or MOs as the guiding functions in Eq. (3). In this subsection we discuss the three kinds of the localized guiding functions , i.e., atomic orbitals, hybrid atomic orbitals, and embedded MOs in molecules and bulks, and how they can be generated in our implementation.
In our implementation, the Bloch function is expanded by PAOs [33, 34] as
(19) |
where and are atomic and orbital indices, and is a linear combination of PAO (LCPAO) coefficient. Also note that , where is the position of atom . We use PAOs as the guiding functions of atomic orbitals, corresponding to valence orbitals or specific orbitals such as localized -orbitals in a narrow energy window.
The choice of the atomic orbitals gives a good guiding function, however, the resultant CWFs may not recover the symmetry of CWFs respecting the bond direction to the neighboring atoms unlike hybrid atomic orbitals due to the closest property of CWFs. In this case, it would be better to use the hybrid atomic orbitals rather than the atomic orbitals as explained below. Let us introduce a projection operator for the occupied space by
(20) |
where is the Fermi-Dirac function. Considering that are a set of non-orthogonal basis set [37], where the overlap integral is denoted by , the density matrix is calculated with the projection operator by
(21) | |||||
where is the dual orbital defined by
(22) |
holding the following orthonormality:
(23) |
By setting and in Eq. (21), and diagonalizing the diagonal block element consisting of , which is associated with selected atomic orbitals on the atomic site , we obtain the hybrid atomic orbitals respecting the symmetry of the bond direction to the neighboring atoms. We employ the hybrid atomic orbitals as the localized functions in our implementation. When the same atomic orbitals are chosen to form the diagonal block element as for the case of atomic orbitals , the resultant CWFs span the same subspace in the Hilbert space.
For some systems the electronic structures can be rather easily understood by employing MOs as building blocks. A molecular crystal is such a case, where the eigenstates near the chemical potential can be approximately constructed by a linear combination of the MOs of constituting molecules. Another example is the bond in molecules and bulks. The bond between two atoms embedded in the system might be analyzed by MOs associated with the two atoms. Here we show how embedded MOs in molecules and bulks can be calculated in a simple procedure. Let us start by noting that using Eq. (20) and assuming that the Bloch function is expressed by Eq. (19) the total number of electrons for a spin degenerate case is given by [37]
(24) | |||||
In the second line of Eq. (24) we used the following identity operator:
(25) |
Since with , it is enough to consider . We utilize the formula for to calculate embedded MOs in molecules and bulks, and rewrite it the sum of partial traces as
(26) |
where is the index of group including PAOs on grouped atoms, e.g., a set of PAOs on a molecule. The notation of stands for a subset of PAOs as
(27) |
where PAOs on all the atoms in the group are included in the subset. The notation of is the counterpart for the dual orbitals. Using the identity operator of Eq. (25), in Eq. (26) is given by
(28) |
with definition of block elements:
(29) | |||||
(30) |
We now introduce the embedded MOs orbitals in molecules and bulks by a right-singular vectors of [38]. Since the elements of are real in case of the collinear DFT, the right-singular vectors can be obtained with real components. However, the SVD of tends to produce the right-singular vectors with complex components, which results in the complex CWFs making analysis complicated. To obtain the real right-singular vectors, we perform the eigendecomposition of as
(31) |
where is the diagonal matrix consisting of the eigenvalues of , and is a unitary matrix consisting of the corresponding eigenvectors . By noting that the SVD of is given by , we see that the right-singular vectors can be obtained by diagonalizing as . The square roots of the eigenvalues of are singular values of , and can be approximately regarded as the occupation of the corresponding eigenvector . If is a real matrix, obtained by the diagonalization of is guaranteed to be a real unitary matrix. We further normalize the eigenvector by considering the overlap integrals of Eq. (30) as
(32) |
and calculate an expectation value of the KS Hamiltonian with as
(33) |
After are ordered based on the expectation values, we employ selected ones among , e.g., ones near the chemical potential, as the guiding functions in Eq. (3) by monitoring the expectation values and the corresponding singular values. They are what we call embedded MOs in molecules and bulks. The reason why the expectation values are used to select the embedded MOs in addition to the singular values is that the singular values are found to be almost similar to each other as long as they are associated with occupied states. As demonstrated in Sec. IV, work as the good guiding functions to capture the electronic structure of a molecular crystal.
IV NUMERICAL RESULTS
IV.1 Wannier Interpolated Bands
We present five numerical examples to demonstrate the broad applicability of the proposed method across various systems. In addition, a series of benchmark calculations for 30 systems are shown as supplementary information on a website [39]. In Figs. 2 (a) and (b), the interpolated bands of Si in the diamond structure, calculated by the TB Hamiltonian derived from CWFs, are compared with those by the conventional calculation. For comparison, the result with the projection-only WFs is also shown in Fig. 2 (c), calculated by the projection method in Ref. [8]. The hybrid atomic orbitals consisting of valence minimal orbitals, which are one - and a set of -, -, and -orbitals per Si, were used as the guiding functions . For the window function of Eq. (4) we use -15.0 and 0.0, relative to the chemical potential, in eV for and , respectively, which covers the energy range of valence bands. On the other hand, two of 0.01 and 3.0 eV were used to check how conduction bands are reproduced depending on the effect of smearing. As clearly seen, the two cases well reproduce the valence bands, and the conduction bands are also reproduced well in the case (b) of eV. We also note that the interpolated bands by the projection-only WFs reproduce well the reference as shown in Fig. 2 (c). Evaluating the DM function per CWF using Eq. (14), we obtain values of 0.571 and 0.444 for of 0.01 and 3.0 eV, respectively. For eV, the first four singular values at each point are close to 1, while the remaining four approach . Conversely, for eV, the eight singular values at each point distribute ranging from 1.2 to 0.1. These differences in the distribution of singular values explain why the value of the DM function for eV is larger than that for eV. One of the obtained CWFs is shown in Fig. 2 (d), which represents a -like CWF. Although the proposed method does not impose the localization of WFs, the decaying property of TB parameters exhibits a subexponential behavior as shown in Fig. 2 (e). In case of eV, the conduction bands are treated less importantly due to heavily weighting to the valence bands, resulting in the poor reproduction of the conduction bands. At first glance, one may consider that this is an undesirable feature in such a treatment. However, the feature enables us to generate the minimal atomic-like orthogonal orbitals, which well span the subspace by the valence bands, and to calculate an effective charge allocated to each atom using the minimal atomic-like orthogonal orbitals in an unbiased way. We will demonstrate the calculation of the effective atomic charges and stress the usefulness of the method later on.
The disentanglement of bands in metals is demonstrated for Cu in the FCC structure as shown in Figs. 3 (a) and (b). For comparison, the result with the projection-only WFs is also shown in Fig. 3 (c). By selecting the parameters and so that the bands are included, and using only five -orbitals as the guiding functions , the five -bands are reproduced as shown in Fig. 3 (a). We see that the -bands are properly disentangled from the dispersive band, and one of the obtained CWFs preserves the shape of -orbital as shown in Fig. 3 (d) in accordance with the subexponential decaying property of the TB parameters as shown in Fig. 3 (e). When -, -, and -orbitals are used as the guiding functions , a broad range of bands are reproduced as shown in Fig. 3 (b). Since our PP of Cu includes , , , and states, the original band structure has the deep - and -bands. We set the parameters in the window function so as to discard the - and -bands, and use the -, -, and -orbitals, not - and -orbitals, as the guiding functions . So, the - and -bands are not included in the interpolated bands (not shown). The examples show that the CWFs can be flexibly and easily constructed in accordance with the purpose of study without numerical difficulties. On the other hand, the interpolated bands by the projection-only Wannier functions (WFs) demonstrate good agreement with the reference, as shown in Fig. 3 (c). This agreement is achieved through the minimization of , which is originally called gauge invariant part of the spread function , as discussed in Ref. [8]. Without minimizing , however, the bands obtained by the projection-only WFs lose smoothness and exhibit slight deviations from the reference (not shown) as discussed in Ref. [3]. Our method achieves good agreement with the reference without any iterative calculations, as illustrated in Fig. 3 (a). This remarkable property can be attributed to the use of the smoothly varying window function of Eq. (4).
The interpolated bands for a molecular crystal, consisting of tetrathiafulvalene (TTF) and tetracyanoquinodimethane (TCNQ), are shown in Figs. 4 (a) and (b). By employing hybrid atomic orbitals as the guiding functions , the wide range of bands are reproduced as shown in Fig. 4 (a). On the other hand, as shown in Fig. 4 (b), only four bands near the chemical potential are reproduced using the embedded MOs in the bulk as the guiding functions . The MOs were calculated by grouping the TTF and TCNQ molecules separately as explained in Sec. III B. The expectation values calculated by Eq. (33) and the singular values of are listed in Table 2. It can be seen that the singular values largely change from to and from to for TTF and TCNQ, respectively. So, the 26th MO and the 37th MO for TTF and TCNQ, respectively, were used as the guiding functions . Since there are two TTF molecules and two TCNQ molecules in the unit cell, we have the four embedded MOs as , resulting in the four bands. In Figs. 4 (c) and (d) two of the resultant CWFs are shown, which localize in TTF and TCNQ molecules, respectively, as expected. The value of the DM function per CWF is found to be 0.022, which implies that the guiding MOs are very close to the resultant CWFs. The example demonstrates that the proposed method provides a direct way to calculate CWFs for targeted bands in molecular crystals.
TTF | TCNQ | ||||
---|---|---|---|---|---|
index | index | ||||
23 | -0.069 | 1.359 | 34 | 0.135 | 1.082 |
24 | 0.719 | 1.346 | 35 | 0.375 | 1.329 |
25 | 0.761 | 1.449 | 36 | 0.477 | 1.336 |
26 | 3.408 | 0.864 | 37 | 4.928 | 0.535 |
27 | 5.947 | 0.047 | 38 | 9.389 | 0.036 |
28 | 6.197 | 0.047 | 39 | 9.795 | 0.019 |
We extend the proposed method to the non-collinear DFT [42, 43] with spin-orbit interaction (SOI) [32]. The KS orbitals are expressed by two-component spinor, and the SOI is introduced by fully relativistic -dependent PPs [28]. The theoretical framework is readily extended to the non-collinear DFT with the SOI without any difficulty. However, it should be noted that the hybrid atomic orbitals and embedded MOs in molecules and bulks as the guiding functions can be complex functions in the case. In Figs. 5 (a) and (b), we show the Wannier interpolated bands of BiSe which is known to be a topological insulator when the SOI is included in the DFT calculation. For the cases without and with the SOI, corresponding to Figs. 5 (a) and (b), respectively, it is confirmed that the interpolated bands accurately reproduce the original ones regardless of existence of the band inversion, demonstrating a wide variety of applicability of the proposed method.
Finally, we show results for GaAs in Figs. 6 (a)-(d) to demonstrate how the window function of Eq. (4) with the two temperature scheme plays a significant role in reproducing the targeted bands. The interpolated bands calculated with the projection-only WFs calculated by the projection method in Ref. [8] reproduce well the reference ones as shown in Fig. 6 (a). When the four bands located just below the chemical potential, which mainly consists of Ga- and As- orbitals, are focused by discarding the other bands with = 0.001 eV, interpolated bands are obtained as shown in Fig. 6 (b). Other the other hand, including outer bands with = 1.0 eV leads to an erratic behavior as shown in Fig. 6 (c), where the CWFs attempt to reproduce the lower band located around -12 eV and low lying unoccupied bands depending on the k-point, resulting in the erratic behavior. Noting that the three bands in between -7 and 0 eV consists of mainly Ga- and As- orbitals and the low lying unoccupied bands are largely contributed by Ga- orbitals, the band around -12 eV needs to be completely discarded and the low lying unoccupied bands should be partly included when Ga-, , and As- orbitals are used as the guiding functions. The result in Fig. 6 (d) is the case, where = 0.001 and = 1.0 eV were used. We see that the interpolated bands fully reproduce the reference ones as expected, which clearly demonstrates the usefulness of the two temperature scheme.
IV.2 Effective Atomic Charges
From the construction, one can calculate CWFs closest to (hybrid) valence atomic orbitals, while fully respecting the occupied states. The feature of CWFs provides us a unique way to evaluate effective atomic charges based on CWFs being orthogonal atomic-like orbitals. Since the CWFs are a set of orthonormal functions, the effective atomic charge of atom is easily calculated without arbitrariness as
(34) | |||||
where the factor of 2 is due to the spin degeneracy, is the projection operator defined by Eq. (20), and is the number of valence electrons in the PP of atom . Though the summation over orbitals belonging to atom is taken into account in Eq. (34), one can analyze the orbital resolved charges as well. Table III shows effective atomic charges in a HCN molecule and the NaCl bulk, calculated by Eq. (34) together with effective charges calculated by MLWFs [3, 7], the projection-only WFs (POWF) [8], the Mulliken population analysis [45], and the Bader analysis [46, 47].
HCN | CWF | MLWF | POWF | MP | Bader | ||||||||||
Basis | H | C | N | H | C | N | H | C | N | H | C | N | H | C | N |
6.0-s1p1 | 0.077 | -0.052 | -0.025 | 0.665 | 0.303 | -0.969 | 0.262 | -0.215 | -0.047 | 0.384 | -0.164 | -0.221 | 1.000 | 1.746 | -2.746 |
6.0-s2p2 | 0.069 | 0.003 | -0.073 | 0.793 | 0.067 | -0.860 | 0.312 | -0.009 | -0.303 | -0.070 | 0.321 | -0.252 | 0.201 | 2.439 | -2.639 |
6.0-s2p2d1 | 0.066 | 0.009 | -0.076 | 0.793 | 0.102 | -0.907 | 0.307 | 0.010 | -0.317 | -0.008 | 0.393 | -0.385 | 0.183 | 2.523 | -2.706 |
6.0-s3p3d2 | 0.066 | 0.008 | -0.074 | 0.699 | 0.343 | -1.042 | 0.314 | 0.001 | -0.315 | 0.110 | 0.298 | -0.408 | 0.182 | 2.539 | -2.721 |
6.0-s3p3d2f2 | 0.066 | 0.008 | -0.074 | 0.721 | 0.344 | -1.065 | 0.308 | 0.007 | -0.315 | 0.167 | 0.297 | -0.464 | 0.182 | 2.532 | -2.715 |
NaCl | CWF | MLWF | POWF | MP | Bader | ||||||||||
Basis | Na | Cl | Na | Cl | Na | Cl | Na | Cl | Na | Cl | |||||
9.0-s2p1 | 0.391 | -0.391 | 0.760 | -0.760 | 0.627 | -0.627 | 0.716 | -0.716 | 0.858 | -0.858 | |||||
9.0-s3p2 | 0.422 | -0.422 | 0.628 | -0.628 | 0.698 | -0.698 | 0.595 | -0.595 | 0.865 | -0.865 | |||||
9.0-s3p3d2 | 0.421 | -0.421 | 0.876 | -0.876 | 0.697 | -0.697 | 0.158 | -0.158 | 0.853 | -0.853 | |||||
9.0-s3p3d2f2 | 0.421 | -0.421 | 0.877 | -0.877 | 0.697 | -0.697 | -0.062 | 0.062 | 0.850 | -0.850 |
Both the systems are known to be notorious cases because of the difficulty in calculating the effective charges [20]. We see that the effective charges by the CWFs quickly converge as a function of basis functions, while those by the Mulliken population are highly dependent on the choice of the basis function [48]. The non-ionic nature of the CN bond in the HCN molecule, obtained by the CWF method, is consistent with the picture given by the intrinsic atomic orbital method (IAO) [20]. On the other hand, it turns out that the effective charges calculated by MLWFs are relatively sensitive to the choice of basis functions compared to the cases of CWF, and exhibit large charge transfer in HCN and NaCl. The resultant ionic CN bond in the MLWF method differs from the picture derived by the IAO and CWF methods. The overestimation of the charge transfer in the MLWF method can be understood by comparing the WFs as shown in Fig. 7. It is confirmed that the CWFs (a)-(i) well preserve the shape of the atomic orbitals as expected. By contrast, the MLWFs (j), (k), (l), (o), and (p) deform significantly from the shape of atomic orbitals as a result of the maximally localization of the spread function. The MLWFs (j) and (k), originating from H- and C- orbitals, respectively, reduces the occupations by the deformation, while the occupations of the MLWFs (l), (o), and (p) increase due to the deformation. The bonds between H-C and C-N are contributed by the MLWFs (l) and (p), respectively, which originate from C-2 and N-2 orbitals, respectively. Thus, the MLWF (j), originating from H- orbital, does not largely contribute to the bonds between H-C, reducing the population on the H atom. Similarly, the population on the N atom increases, since the bond between C-N is formed by the MLWF (p), originating from N-. It is apparent that attributing the occupation of the MLWFs (l) and (p) to a single atom is not justified due to the fact that they are bond-centered. This implies that MLWFs should not be used to decompose electron population to each atom as demonstrated for the HCN molecule. On the other hand, CWFs are closest to atomic orbitals when the atomic orbitals are used as the guiding functions, and the population calculated by the CWFs can be naturally attributed to each atom, which allows us to calculate the atomic populations in a physically convincing manner.
For comparison, we also calculated the effective charges using the projection-only Wannier functions (POWF) [8] and found that the obtained values are relatively sensitive to the choice of basis functions and lie between the results obtained by the CWF and MLWF methods. Here, we analyze why the results from the POWF deviate from those by the CWF even when the inner window scheme [8] was used to focus on the valence states. In the HCN molecule, the projection of the -orbital for hydrogen and the -, -, -, -orbitals for carbon and nitrogen results in nine Wannier functions. Since there are five occupied orbitals, applying the inner window scheme to the energy range of these orbitals first selects a subspace spanned by them. Subsequently, four vectors are chosen from the subspace consisting of unoccupied orbitals so that they can be orthogonal to the former subspace and have maximal overlap with the subspace obtained by projection of guiding functions onto a subspace spanned by states included in the outer window. Once the latter subspace is selected, each subspace is given equal weight in constructing the projection matrix. As a result, the nine Wannier functions generated are not specifically designed to primarily accommodate the occupied orbitals. This statement can be verified using the CWF method: By setting the window function with eV, eV, and eV in the CWF method, and calculating the effective charges, contributions from low-lying unoccupied orbitals are included equally. We confirmed that this calculation results in effective charges nearly equivalent to those obtained by the projection-only method (not shown). It is evident that even when using the inner window scheme with the projection-only approach, properly computing effective charges remains a challenge. To calculate the effective charges assigned to atoms, it is reasonable to use atomic orbitals as projection orbitals. Typically, the number of atomic orbitals used exceeds the number of occupied orbitals, leading to similar challenges as observed in the case of NaCl. A simple solution to resolve this problem is to introduce a smoothly varying window function even in the inner window scheme when constructing the projection matrix. In fact, we confirmed that introducing a smoothly varying window function in the inner window scheme almost reproduces those by the CWF method (not shown). However, one should note that the use of the smoothly varying window function in the CWF is a much simpler method.
The Bader analysis shows erratic charge transfer in the CN bond, implying the limited applicability of the Bader analysis based on a real-space decomposition to the triple bond [20], while the effective atomic charges in NaCl are similar to those by the MLWF method.
It is considered from the comparison that the effective atomic charges obtained from the CWFs could serve as a valuable tool for analyzing electronic structures in a manner related to the approach discussed in Ref. [20]. Further investigation in this direction will be conducted in future studies.
V CONCLUSIONS
We presented a non-iterative method to calculate the closest Wannier functions (CWFs) to a given set of localized guiding functions such as atomic orbitals, hybrid atomic orbitals, and molecular orbitals (MOs). We defined the distance measure (DM) function by the sum of the squared distance between the projection functions and Wannier functions in the Hilbert space, and considered minimizing the DM function. It was shown that the minimization of the DM function is achieved by the polar decomposition of the projection matrix with a window function via the singular value decomposition (SVD) in a non-iterative manner. The CWFs can be uniquely constructed as long as the projection matrix is positive definite. It was also shown that the method is free from subtle choice of the gauge. The disentanglement of bands is naturally taken into account by introducing a smoothly varying window function, and including Bloch functions to generate CWFs, where . Even for isolated bands, the method always performs the disentanglement of bands without any iterative calculation and difficulty. The advantage of the smoothly varying window function over the inner window scheme [8] was clearly demonstrated in the disentanglement of -bands in copper and in the calculations of the effective charges. Also, it is noted that the implementation of the method is relatively simple, which requires mainly matrix-matrix product and SVD. We have implemented the proposed method into the OpenMX code, which is based on the density functional theory (DFT), numerical pseudo-atomic orbitals (PAOs), and norm-conserving pseudopotentials (PPs), and introduced three types of guiding functions, i.e., atomic orbitals, hybrid atomic orbitals, and embedded MOs in molecules and bulks. The first two are easily employed from the PAOs. For the last one, we developed a method to calculate embedded MOs in molecules and bulks, which focuses on a partial trace formula of the projection operator for the occupied space and applies SVD to the partial matrix for the projection operator. The interpolated bands by tight-binding (TB) models derived from the CWFs reproduce well the targeted conventional bands of a wide variety of systems including Si, Cu, the TTF-TCNQ molecular crystal, and a topological insulator of BiSe. These successful reproduction of targeted bands clearly demonstrates a wide variety of applicability of the proposed method. We further show the usefulness of the proposed method in calculating effective atomic charges, implying that the CWFs closest to atomic orbitals can be used as a measure to analyze electronic structures from one system to the others. Thus, we conclude that the proposed method is an alternative way in efficiently calculating WFs, and the concept of CWFs will provide a basis for development of novel methods of analyzing electronic structures and calculating physical properties [49, 50].
Acknowledgements.
The author wishes to express gratitude to Prof. Yoshiaki Sugimoto for inspiring the development of CWFs through collaborative research. Thanks are also extended to Mr. Ryotaro Koshoji for his insightful comments on the theoretical aspect. Part of the computation in this study was carried out using the computational facility of the Institute for Solid State Physics at the University of Tokyo.Appendix A: The upper bound of the diagonal elements of the matrix
A proof for the upper bound of the diagonal elements of the matrix is given in the appendix. Note that and in Eq. (11) are partial unitary matrices in size of . Since is hermitian, its eigenvalues are real. Let be an eigenvalue of and be the corresponding eigenvector. Then, the following equation holds:
(35) |
Defining and , and operating from the left side of Eq. (35), we have
(36) |
Noting , and using the Cauchy-Schwarz inequality , we obtain
(37) |
Since , combining Eq. (36) and Eq. (37) results in
(38) |
Considering , the upper bound of the eigenvalues of is found to be
(39) |
Noting that the matrix can be written by and as
(40) |
the diagonal elements of the matrix is given by
(41) |
where is the -th element in the vector . Since the upper bound of is unity, and forms a unitary matrix, we obtain the upper bound of the diagonal elements of the matrix as
(42) |
References
- [1] G.H. Wannier, The structure of electronic excitation levels in insulating crystals, Phys. Rev. 52, 191 (1937).
- [2] W. Kohn, Analytic properties of bloch waves and wannier functions, Phys. Rev. 115, 809 (1959).
- [3] N. Marzari, A.A. Mostofi, J.R. Yates, I. Souza, and D. Vanderbilt, Maximally localized wannier functions: Theory and applications Rev. Mod. Phys. 84, 1419 (2012).
- [4] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
- [5] D.R. Hamann and D. Vanderbilt, Maximally localized wannier functions for GW quasiparticles, Phys. Rev. B 79, 045109 (2009).
- [6] F. Lechermann, A. Georges, A. Poteryaev, S. Biermann, M. Posternak, A. Yamasaki, and O.K. Andersen, Dynamical mean-field theory using Wannier functions: A flexible route to electronic structure calculations of strongly correlated materials, Phys. Rev. B 74, 125120 (2006).
- [7] N. Marzari and D. Vanderbilt, Maximally localized generalized wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997).
- [8] I. Souza, N. Marzari, and D. Vanderbilt, Maximally localized wannier functions for entangled energy bands, Phys. Rev. B 65, 035109 (2001).
- [9] V. Vitale, G. Pizzi, A. Marrazzo, J.R. Yates, N. Marzari, and A.A. Mostofi, Automated high-throughput wannierisation, npj Computational Materials 6, 66 (2020).
- [10] A. Damle and L. Lin, Disentanglement via entanglement: A unified method for wannier localization, Multiscale Model. Simul. 16, 1392 (2018).
- [11] J. Qiao, G. Pizzi, and N. Marzari, Automated mixing of maximally localized wannier functions into target manifolds, npj Computational Materials 9, 206 (2023).
- [12] W. Ku, H. Rosner, W.E. Pickett, and R.T. Scalettar, Insulating ferromagnetism in LaBaCuO: An ab initio wannier function analysis, Phys. Rev. Lett. 89, 167204 (2002).
- [13] W.C. Lu, C.Z. Wang, T.L. Chan, K. Ruedenberg, and K.M. Ho, Representation of electronic structures in crystals in terms of highly localized quasiatomic minimal basis orbitals, Phys. Rev. B 70, 041101(R) (2004).
- [14] X. Qian, J. Li, L. Qi, C.-Z. Wang, T.-L. Chan, Y.-X. Yao, K.-M. Ho, and S. Yip, Quasiatomic orbitals for ab initio tight-binding analysis, Phys. Rev. B 78, 245112 (2008).
- [15] P.-O. Löwdin, On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals, J. Chem. Phys. 18, 365 (1950).
- [16] J. Des Cloizeaux, Energy bands and projection operators in a crystal: Analytic and asymptotic properties, Phys. Rev. 135, A685 (1964).
- [17] A.A. Mostofi, J.R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Wannier90: A tool for obtaining maximally-localised wannier functions, Comput. Phys. Commun. 178, 685 (2008).
- [18] B.C. Carlson and J.M. Keller, Orthogonalization procedures and the localization of wannier functions, Phys. Rev. 105, 102 (1957).
- [19] A.E. Reed, L.A. Curtiss, and F. Weinhold, Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint, Chem. Rev. 88, 899 (1988).
- [20] G. Knizia, Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts, J. Chem. Theory Comput. 9, 4834 (2013).
- [21] M. Engel, M. Marsman, C. Franchini, and G. Kresse, Electron-phonon interactions using the projector augmented-wave method and wannier functions, Phys. Rev. B 101, 184302 (2020).
- [22] A.A. Mostofi, J.R. Yates, G. Pizzi, Y.S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, An updated version of wannier90: A tool for obtaining maximally-localised wannier functions, Comput. Phys. Commun. 185, 2309 (2014).
- [23] T. Schäfer, F. Libisch, G. Kresse, A. Grüneis, Local embedding of coupled cluster theory into the random phase approximation using plane waves, J. Chem. Phys. 154, 011101 (2021).
- [24] M. Engel, H. Miranda, L. Chaput, A. Togo, C. Verdi, M. Marsman, and G. Kresse, Zero-point renormalization of the band gap of semiconductors and insulators using the projector augmented wave method, Phys. Rev. B 106, 094316 (2022).
- [25] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
- [26] Matrix Computations (3rd ed.). Baltimore: The Johns Hopkins University Press (1996), G. Golub and C.F Van Loan.
- [27] K. Fan and A. J. Hoffman, Some metric inequalities in the space of matrices, Proc. Amer. Math. Soc. 6, 111 (1955).
- [28] The code, OpenMX, pseudo-atomic basis functions, and pseudopotentials are available under terms of the GNU-GPL on a website (http://www.openmx-square.org/).
- [29] T. Ozaki and H. Kino, Efficient projector expansion for the ab initio LCAO method, Phys. Rev. B 72, 045121 (2005).
- [30] T.V.T. Duy and T. Ozaki, A three-dimensional domain decomposition method for large-scale DFT electronic structure calculations, Comput. Phys. Commun. 185, 777 (2014).
- [31] I. Morrison, D.M. Bylander, and L. Kleinman, Nonlocal hermitian norm-conserving vanderbilt pseudopotential, Phys. Rev. B 47, 6728 (1993).
- [32] G. Theurich and N.A. Hill, Self-consistent treatment of spin-orbit coupling in solids using relativistic fully separable ab initio pseudopotentials, Phys. Rev. B 64, 073106 (2001).
- [33] T. Ozaki, Variationally optimized atomic orbitals for large-scale electronic structures, Phys. Rev. B. 67, 155108, (2003).
- [34] T. Ozaki and H. Kino, Numerical atomic basis orbitals from H to Kr, Phys. Rev. B 69, 195113 (2004).
- [35] K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I.E. Castelli, S.J. Clark, A. Dal Corso, S. de Gironcoli, T. Deutsch, J.K. Dewhurst, I. Di Marco, C. Draxl, M. Dułak, O. Eriksson, J.A. Flores-Livas, K.F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E.K. Gross, A. Gulans, F. Gygi, D.R. Hamann, P.J. Hasnip, N.A. Holzwarth, D. Iu̧san, D.B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçü kbenli, Y.O. Kvashnin, I.L. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordströ m, T. Ozaki, L. Paulatto, C.J. Pickard, W. Poelmans, M.I. Probert, K. Refson, M. Richter, G.M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunströ m, A. Tkatchenko, M. Torrent, D. Vanderbilt, M.J. van Setten, V. Van Speybroeck, J.M. Wills, J.R. Yates, G.X. Zhang, and S. Cottenier, Reproducibility in density functional theory calculations of solids, Science 351, aad3000 (2016).
- [36] J.P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
- [37] T. Ozaki, M. Fukuda, and G. Jiang, Efficient O() divide-conquer method with localized single-particle natural orbitals, Phys. Rev. B 98, 245137 (2018).
- [38] As discussed in Ref. [37], strictly speaking the right eigenvectors obtained by the eigendecomposition of correspond to embedded MOs, since the eigenvalues are directly related to the occupations. However, the right eigenvectors can be complex even for the collinear case. Thus, we use the SVD which allows us to obtain real vectors at least within the collinear case. We confirmed from numerical experiments that the right singular vectors are qualitatively similar to the right eigenvectors.
- [39] A series of benchmark results are available on a website (http://www.openmx-square.org/cwf/).
- [40] Y.-T. Lee and T. Ozaki, OpenMX Viewer: A web-based crystalline and molecular graphical user interface program, J. Mol. Graph. Model. 89, 192 (2019); https://www.openmx-square.org/viewer/
- [41] T.J. Kistenmacher, T.E. Phillips, and D.O. Cowan, The crystal structure of the 1:1 radical cation-radical anion salt of 2,2’-bis-l,3-dithiole (TTF) and 7,7,8,8-tetracyanoquinodimethane (TCNQ), Acta Cryst. B30, 763 (1974).
- [42] U. von Barth and L. Hedin, A local exchange-correlation potential for the spin polarized case. i, J. Phys. C: Solid State Phys. 5, 1629 (1972).
- [43] J. Kubler, K.-H. Hock, J. Sticht, and A. R. Williams, Density functional theory of non-collinear magnetism, J. Phys. F: Met. Phys. 18, 469 (1988).
- [44] S. Nakajima, The crystal structure of BiTeSe, J. Phys. Chem. Solids 24, 479 (1963).
- [45] R.S. Mulliken, Electronic population analysis on LCAO–MO molecular wave functions. I, J. Chem. Phys. 23, 1833 (1955).
- [46] Atoms in Molecules: A Quantum Theory. Oxford University Press, New York (1990), R.F.W. Bader.
- [47] G. Henkelman, A. Arnaldsson, and H. Jónsson, A fast and robust algorithm for bader decomposition of charge density, Comput. Mater. Sci. 36, 354 (2006).
- [48] R.S. Mulliken, Iodine revisited, J. Chem. Phys. 55, 288 (1971).
- [49] R. Zhang, Y. Yasui, M. Fukuda, T. Ozaki, M. Ogura, T. Makino, D. Takeuchi, and Y. Sugimoto, Atomic observation on diamond (001) surfaces with non-contact atomic force microscopy, arXiv:2403.17427.
- [50] A. Inda, R. Oiwa, S. Hayami, H.M. Yamamoto, and H. Kusunose, Quantification of chirality based on electric toroidal monopole, J. Chem. Phys. 160, 184117 (2024).