Efficient algorithm based on Liechtenstein method for computing exchange coupling constants using localized basis set

TL;DR

This paper introduces an efficient algorithm for calculating exchange coupling constants using a localized basis set and a simplified energy integration method, significantly reducing computational cost while maintaining accuracy.

Contribution

The authors adapt the Liechtenstein formula with a finite pole approximation for localized orbitals, improving computational efficiency for large-scale magnetic system simulations.

Findings

Number of poles needed for 0.05 meV accuracy is ~60

Compared distance dependence of J_ij with KKR formalism

Algorithm reduces computational cost by an order of magnitude

Abstract

For large-scale computation of the exchange coupling constants $J_{ij}$, we reconstruct the Liechtenstein formula for localized orbital representation and simplify the energy integrations by adopting the finite pole approximation of the Fermi function proposed by Ozaki [Phys. Rev. B 75, 035123 (2007)]. We calculate the exchange coupling constant $J_{\mathrm{1NN}}$ of the first-nearest-neighbor sites in body-centered-cubic Fe systems of various sizes to estimate the optimal computational parameters that yield appropriate values at the lowest computational cost. It is shown that the number of poles needed for a computational accuracy of 0.05 meV is determined as $\sim$ 60, whereas the number of necessary Matsubara poles needed to obtain similar accuracy, which was determined in previous studies, is on the order of 1000. Finally, we show $J_{ij}$ as a function of atomic distance, and compared it with one derived from Korringa-Kohn-Rostoker Green's function formalism. The distance profile of $J_{ij}$ derived by KKR formalism agrees well with that derived by our study, and this agreement supports the reliability of our newly derived formalism.

Tables

Table 1. Table 1: Numbers of real and reciprocal space grids, N r × N r × N r subscript 𝑁 𝑟 subscript 𝑁 𝑟 subscript 𝑁 𝑟 N_{r}\times N_{r}\times N_{r} and N k × N k × N k subscript 𝑁 𝑘 subscript 𝑁 𝑘 subscript 𝑁 𝑘 N_{k}\times N_{k}\times N_{k} , for SCF calculations, together with the number of atoms and length of the unit cell vector for each system.

System	$N_r\times N_r\times N_r$	$N_k\times N_k\times N_k$	$N_{\mathrm{atom}}$	$a$ [Å]
$1\times 1\times 1$	$32\times 32\times 32$	$24\times 24\times 24$	2	2.866
$2\times 2\times 2$	$64\times 64\times 64$	$12\times 12\times 12$	16	5.732
$3\times 3\times 3$	$96\times 96\times 96$	$8\times 8\times 8$	54	8.598
$4\times 4\times 4$	$128\times 128\times 128$	$6\times 6\times 6$	128	11.464

Table 2. Table 2: Optimized values of computational parameters N P subscript 𝑁 P N_{\mathrm{P}} and N k subscript 𝑁 𝑘 N_{k} with corresponding computational times t eig subscript 𝑡 eig t_{\mathrm{eig}} for eigenvalue calculation and t indiv subscript 𝑡 indiv t_{\mathrm{indiv}} for J i ​ j subscript 𝐽 𝑖 𝑗 J_{ij} calculation based on Eq. ( 19 ), and numbers of nodes and cores assigned in each calculation. All calculations are performed with MPI parallelization on 24 cores on 1 node of Intel(R) Xeon(R) E5-2680 v3 processors.

System	$N_P^{\mathrm{opt}}$	$N_k^{\mathrm{opt}}$	$t_{\mathrm{eig}}$ [s]	$t_{\mathrm{indiv}}$ [s]
$1\times 1\times 1$	60	29	1.5	1.9
$2\times 2\times 2$	60	15	30.1	1.9
$3\times 3\times 3$	60	8	193.2	1.3
$4\times 4\times 4$	60	5	906.2	1.1

Table 3. Table 3: Difference Δ ​ J 1 ​ N ​ N Δ subscript 𝐽 1 N N \Delta J_{\mathrm{1NN}} between the J 1 ​ N ​ N subscript 𝐽 1 N N J_{\mathrm{1NN}} values for periodic images and those for individual pairs in systems of various sizes, together with the computational times t periodic subscript 𝑡 periodic t_{\mathrm{periodic}} and t indiv subscript 𝑡 indiv t_{\mathrm{indiv}} of the periodic image and individual pair calculations, respectively. For each system, the numbers of poles and k 𝑘 k -points are the same as those for the individual pair calculations shown in Table 2 .

System	$\Delta J_{1NN}$ [meV]	$t_{\mathrm{periodic}}$ [s]	$t_{\mathrm{indiv}}$ [s]
$1\times 1\times 1$	117.95	0.2	1.9
$2\times 2\times 2$	0.30	1.5	1.9
$3\times 3\times 3$	0.41	4.6	1.3
$4\times 4\times 4$	0.53	8.0	1.1

Equations

$$J_{ij}=\frac{1}{4\pi }\int d\epsilon f\left(\beta (\epsilon -{\epsilon }_{\mathrm{F}})\right)\mathrm{Im}\mathrm{Tr}\left[{\hat{\Delta }}_i{\hat{T}}_{\uparrow }^{ij}{\hat{\Delta }}_j{\hat{T}}_{\downarrow }^{ji}\right],$$

$$J_{ij}$$

$${\hat{P}}_i\equiv {\hat{H}}_{i\uparrow }-{\hat{H}}_{i\downarrow },$$

$$J_{i0,j\mathbf{R}}$$

$${\hat{P}}_{i\mathbf{R}}\equiv {\hat{H}}_{i\mathbf{R}}^{(\mathrm{KS})}(\uparrow )-{\hat{H}}_{i\mathbf{R}}^{(\mathrm{KS})}(\downarrow ),$$

$$|\mathbf{k},n,\sigma ⟩=\frac{1}{\sqrt{N}}\sum _{\mathbf{R}}\sum _i\sum _{\mu \in i}|\mathbf{R},i,\mu ⟩e^{i\mathbf{k}\cdot \mathbf{R}}{C}_{i\mu ,n\sigma }(\mathbf{k}),$$

$$J_{i0,j\mathbf{R}}$$

$$G_{j{\mu }^{\prime },i\nu }^{+}(\uparrow ,\epsilon ,-\mathbf{R})$$

$$=\int d^3\left(\frac{ka}{2\pi }\right)e^{i\mathbf{k}\cdot \mathbf{R}}\sum _n\frac{C_{j{\mu }^{\prime },n\uparrow }(\mathbf{k})C_{i\mu ,n\uparrow }(\mathbf{k})}{\epsilon +i\eta -{\epsilon }_{n\uparrow }(\mathbf{k})}$$

$$G_{i\mu ,j{\nu }^{\prime }}^{+}(\downarrow ,\epsilon ,\mathbf{R})$$

$$=\int d^3\left(\frac{ka}{2\pi }\right)e^{-i\mathbf{k}\cdot \mathbf{R}}\sum _{n^{\prime }}\frac{C_{i\mu ,n^{\prime }\downarrow }(\mathbf{k})C_{j{\nu }^{\prime },n^{\prime }\downarrow }(\mathbf{k})}{\epsilon +i\eta -{\epsilon }_{n^{\prime }\downarrow }(\mathbf{k})}.$$

$${\int }_{C_1}\mathcal{F}(z)dz-{\int }_{C_2}\mathcal{F}(z)dz=2\pi i\sum _{z_p\in S}{R}_p[\mathcal{F}(z);z],$$

$$J_{i0,j\mathbf{R}}$$

$$\mathcal{F}(z)$$

$${\tilde{f}}_{N_{\mathrm{P}}}(z)=\frac{1}{2}+\sum _{p=1}^{N_{\mathrm{P}}}\frac{{\tilde{R}}_p}{z-{\tilde{z}}_p}+\sum _{p=1}^{N_{\mathrm{P}}}\frac{{\tilde{R}}_p}{z+{\tilde{z}}_p},$$

$${\mathbf{Av}}_p={\lambda }_p{\mathbf{Bv}}_p,(p=1,2,\dots ,N_{\mathrm{P}}),$$

$$A_{q{q}^{\prime }}$$

$$B_{q{q}^{\prime }}$$

$${\tilde{z}}_p=\frac{i}{{\lambda }_p},{\tilde{R}}_p=\frac{1}{4}{v}_{pp}^2{\tilde{z}}_p^2.$$

$$J_{i0,j\mathbf{R}}$$

$$J_{ij}\equiv \sum _R{J}_{i0,j\mathbf{R}}.$$

$$J_{ij}$$

$$\Delta J_{1NN}\equiv \left|J_{1NN}^{\mathrm{indiv}}-J_{1NN}^{\mathrm{periodic}}\right|.$$

$$\frac{2}{3k_{\mathrm{B}}}\left(\begin{array}[]{cc}J_{11}-J_{1\mathbf{0},1\mathbf{0}}&J_{12}\\ J_{21}&J_{22}-J_{2\mathbf{0},2\mathbf{0}}\end{array}\right).$$

$$\sum _{\mathbf{k},n,\sigma }|\mathbf{k},n,\sigma ⟩⟨\mathbf{k},n,\sigma |=1,$$

$$\sum _{i,\mu ,\mathbf{R}}\sum _{j,\nu ,{\mathbf{R}}^{\prime }}|i,\mu ,\mathbf{R}⟩[{\mathbf{S}}_{\mathrm{RS}}^{-1}{]}_{i\mu \mathbf{R},j\nu {\mathbf{R}}^{\prime }}⟨j,\nu ,{\mathbf{R}}^{\prime }|=1,$$

$$J_{i0,j\mathbf{R}}$$

$$=\frac{1}{4\pi }\sum _{\mathbf{k},{\mathbf{k}}^{\prime },{\mathbf{k}}^{\prime \prime },{\mathbf{k}}^{\prime \prime \prime }}\sum _{n,n^{\prime },n^{\prime \prime },n^{\prime \prime \prime }}\sum _{\sigma ,{\sigma }^{\prime },{\sigma }^{\prime \prime },{\sigma }^{\prime \prime \prime }}\int d\epsilon f\left(\beta (\epsilon -{\epsilon }_{\mathrm{F}})\right)$$

$$\times \mathrm{ImTr}[⟨{\mathbf{k}}^{\prime \prime \prime },n^{\prime \prime \prime },{\sigma }^{\prime \prime \prime }|{\hat{G}}_{0\uparrow }^{+}(\epsilon )|\mathbf{k},n,\sigma ⟩$$

$$\times ⟨\mathbf{k},n,\sigma |{\hat{P}}_{i,0}|{\mathbf{k}}^{\prime },n^{\prime }{\sigma }^{\prime }⟩$$

$$\times ⟨{\mathbf{k}}^{\prime },n^{\prime },{\sigma }^{\prime }|{\hat{G}}_{0\downarrow }^{+}(\epsilon )|{\mathbf{k}}^{\prime \prime },n^{\prime \prime },{\sigma }^{\prime \prime }⟩$$

$$\times ⟨{\mathbf{k}}^{\prime \prime },n^{\prime \prime },{\sigma }^{\prime \prime }|{\hat{P}}_{j,\mathbf{R}}|{\mathbf{k}}^{\prime \prime \prime },n^{\prime \prime \prime },{\sigma }^{\prime \prime \prime }⟩].$$

$$⟨{\mathbf{k}}^{\prime \prime \prime },n^{\prime \prime \prime },{\sigma }^{\prime \prime \prime }|{\hat{G}}_{0\uparrow }^{+}(\epsilon )|\mathbf{k},n,\sigma ⟩$$

Full text

I Introduction

Considering recent industrial demands, detailed study of the structural and electronic properties of magnetic materials has become an important issue in materials science. In particular, theoretical understanding of permanent magnets is one of the most important yet difficult topics, because of the fundamental complexity of permanent magnets. Recent studies of permanent magnets have revealed complicated material microstructure in Nd–Fe–B-type permanent magnets S_Sugimoto_2011 ; K_Hono_2012 ; S_Hirosawa_2017 . Specifically, grain boundaries are crucial to enhancing the coercivity of permanent magnets because of their pinning effect, which prevents the movement of magnetic domain walls S_Sugimoto_2011 ; K_Hono_2012 ; S_Hirosawa_2017 ; S_Li_2002 ; W_F_Li_2009 ; T_H_Kim_2012 ; H_Sepehri_Amin_2012 ; U_M_R_Seelam_2016 . Moreover, it was reported that both the crystallinity and composition ratio of the grain boundary phase change depending on the relative angles with respect to the $c$ axes of neighboring grains of the main phase T_T_Sasaki_2016 .

Recently, first principles computation techniques have been employed to investigate the magnetic properties of permanent magnets B_Balasubramanian_2016 ; A_Saengdeejing_2016 ; Y_Tatetsu_2016 ; Z_Torbatian_2016 ; N_Umetsu_2016 ; H_Akai_2018 ; Y_Gohda_2018 ; C_E_Patrick_2018 ; Y_Tatetsu_2018 ; C_E_Patrick_2019 ; A_M_Schonhobel_2019 ; A_L_Tedstone_2019 . It is, however, still difficult to understand the details of magnetic interaction in the microstructures of permanent magnets because of their structural complexity.

The Liechtenstein method is a powerful tool that uses second-order perturbation theory to treat the magnetism of complicated systems by examining the response of the total energy to rotation of the spins at two atomic sites through infinitesimal angles A_I_Liechtenstein_1987 . The original Liechtenstein formula is written as follows:

[TABLE]

where $T^{ij}_{\sigma}$ is the path integral operator between sites $i$ and $j$, $\hat{\Delta}_{i}\equiv\hat{t}_{i\uparrow}-\hat{t}_{i\downarrow}$ is the single-site scattering matrix $\hat{t}_{i\sigma}$ at site $i$, $f(x)=1/[\exp(x)+1]$ is the Fermi function with $\beta=(k_{\mathrm{B}}T)^{-1}$, and $\varepsilon_{\mathrm{F}}$ is the Fermi level. This formula was also rewritten within the Green’s function formalism as follows M_Pajda_2001 ; I_Turek_2003 :

[TABLE]

[TABLE]

where $\hat{G}_{\sigma}^{+}(\varepsilon)$ is the retarded Green’s function of the spin $\sigma$ in unperturbed states, and $\hat{H}_{i\sigma}$ is the on-site term of the Hamiltonian for spin $\sigma$ at site $i$.

Owing to use of the Liechtenstein formulae and the development of first-principles calculations of magnetic materials, the magnetic properties of various magnetic materials have begun to be revealed V_A_Dinh_2009 ; I_Galanakis_2011 ; M_Seike_2012 ; T_Fukushima_2015 . It is, however, difficult to understand the magnetic properties of various phases and their interfaces in permanent magnets because the structural details of subphases in permanent magnets are still under investigation T_T_Sasaki_2016 ; X_D_Xu_2018 , and the interfacial structures between multiple phases tend to become very complicated for large numbers of atoms. Therefore, quantitative understanding of the magnetism of an entire permanent magnet contains many related problems to be solved, and the formulation of the exchange coupling constants for large-scale computation is a significant problem associated with this challenging task.

In this paper, we derive the explicit forms of the Liechtenstein formula on the basis of non-orthogonal localized orbitals for the interaction between two individual sites located in different cells and for the interaction between periodic images. For the former derivation, we simplify the formalism by adopting the approximate form of the Fermi function proposed by Ozaki T_Ozaki_2007 . We implement this formalism in an MPI code that uses the overlap and Hamiltonian matrices of OpenMX T_Ozaki_2003 , a first-principles calculation code based on a linear combination of pseudoatomic orbitals approximation. To determine the most efficient computational conditions, we examine the calculated values of $J_{\mathrm{1NN}}$ of the first-nearest-neighbor (1NN) sites in body-centered cubic (bcc) Fe crystals in systems of various sizes and under different computational conditions. We find that it is necessary to take approximately 60 poles of the approximated Fermi functions, whereas the number of Matsubara poles needed for sufficient accuracy was reported to be approximately 1000 or more by Kvashnin et al. Y_O_Kvashnin_2015 . Finally, we calculate the dependence of $J_{ij}$ on the atomic distance $r_{ij}$, and compared it with a $J_{ij}$ profile of bcc Fe obtained by Akai-KKR, a first-principles calculation code based on Korringa-Kohn-Rostoker (KKR) Green’s function formalism H_Shiba_1971 ; H_Akai_1977H_Akai_1982; J_Korringa_1947 ; W_Kohn_1954 . It is shown that the $J_{ij}$ profile in this study agrees well with that derived by Akai-KKR. Considering the variation of calculated $J_{ij}$ among different previous studies M_Pajda_2001 ; H_Hang_2010 ; Y_O_Kvashnin_2015 ; H_Yoon_2018 , it is possible to say that the agreement between two formalism support the reliability of calculated values with each other.

II Modification of Liechtenstein formula

II.1 Non-orthogonal basis representation

Our aim in this section is to represent the Liechtenstein formula by a localized basis set with an efficient algorithm for large-scale computations. As a starting point, let us consider the exchange coupling constant between site $i$ in cell $\mathbf{0}$ and site $j$ in cell $\mathbf{R}$, which are shown in Fig. 1(a). With the explicit cell indices $\mathbf{0}$ and $\mathbf{R}$, Eqs. (2) and (3) can be rewritten as

[TABLE]

[TABLE]

where $\hat{H}^{(\mathrm{KS})}_{i\mathbf{R}}(\sigma)$ represents the on-site partial matrix of the Kohn–Sham Hamiltonian of spin $\sigma$ at site $i$ of cell $\mathbf{R}$. To treat Eq. (4) with a localized orbital basis set, we used the expansion of Bloch functions by the localized orbitals,

[TABLE]

where $\mu$ represents the index of the orbital belonging to site $i$, and $N$ represents the number of cells related to the periodic boundary condition. As a result, Eq. (4) can be written in the following explicit form:

[TABLE]

where $\varepsilon_{n\sigma}(\mathbf{k}),\mathbf{C}_{n\sigma}(\mathbf{k})$ represent the corresponding eigenvalue and eigenvector, respectively, indexed by $n$ and $\sigma$ for the Kohn–Sham equation at the wave number $\mathbf{k}$; and $[\hat{P}_{i}]_{\nu\mu}$ and $[\hat{P}_{j}]_{\nu^{\prime}\mu^{\prime}}$ represent the partial matrices of the potential difference operator at sites $i$ and $j$, respectively. We show that Eq. (7) is relevant not only for orthogonal basis sets but also for non-orthogonal basis sets; a detailed derivation is given in Appendix A.

There are two ways to simplify Eq. (7). One is the eigenfunction representation,

[TABLE]

which can be obtained by replacing the imaginary parts of the Green’s functions with delta functions around the energy eigenvalues M_J_Han_2004 . Direct implementation of Eq. (8), however, may be computationally costly because of the double integrals of two wave vectors.

The other representation of Eq. (7) is the Green’s function representation, which is written as follows:

[TABLE]

where

[TABLE]

Equation (9) involves the integral for the energy $\varepsilon$, and Eqs. (10) and (11) involve single integrals for the wave vector $\mathbf{k}$.

II.2 Contour integrals and finite pole approximation of Fermi function

Because the integrand of Eq. (9) involves only the terms of the retarded Green’s function, it is possible to reduce the computational cost by using the residue theorem instead of real axis integration. The fundamentals of the complex contour integration of the retarded Green’s function are illustrated in Fig. 2.

To estimate the computational cost of energy integration, it is necessary to evaluate the steepness of the integrand. Because the Green’s function, and thus the integrand, of Eq. (9) have poles just infinitesimally below the real axis (which is shown as $C_{1}$ in Fig. 2), direct integration on the real axis would require very fine energy grid points. Instead, the retarded Green’s function changes gradually at the half-circle complex contour ($C_{2}$ in Fig. 2), because the contour is far from all the poles of the retarded Green’s functions. From the residue theorem, the difference between integration on $C_{1}$ and on $C_{2}$ can be written as the summation of the residuals of the integrand:

[TABLE]

where $\mathcal{F}(z)$ is the integrand, $R_{p}[\mathcal{F}(z);z]$ is the residue of the integrand $\mathcal{F}(z)$ indexed by $p$ with the corresponding pole $z_{p}$, and $z_{p}\in S$ means that the pole $z_{p}$ is included in the closed area surrounded by $C_{1}$ and $C_{2}$ (pale blue area in Fig. 2). Given a sufficiently large radius $\varepsilon_{\mathrm{max}}$ of the contour $C_{2}$, it is possible to replace real axis integration with complex contour integration and summation of the residuals on the poles:

[TABLE]

with

[TABLE]

For the first term of Eq. (12), one can easily find that this term becomes zero at the limit $\varepsilon_{\mathrm{max}}\rightarrow+\infty$ (see Appendix B). By contrast, we have to consider all the poles of the Fermi function in the upper complex plane for the second term of Eq. (12). In this way, the Matsubara approximation, a standard approximation of the Fermi function, would result in slow convergence with respect to the number of poles. This is because the Matsubara approximation involves all the exact Fermi poles equally spaced on the imaginary axis, where the remaining part of the integrand decays as $z^{-2}$. Instead, we adopted the finite pole approximation of the Fermi function proposed by Ozaki T_Ozaki_2007 . That is, the approximated Fermi function can be written as in the summation of a finite number of fractions:

[TABLE]

where $z=\pm\tilde{z}_{p}$, and $\tilde{R}_{p}$ is the $p$-th positive/negative pole and the corresponding residue of the approximated Fermi function $\tilde{f}_{N_{\mathrm{P}}}(z)$. This approximated Fermi function is derived from a hypergeometric function, and the poles and residues are calculated from the $N_{\mathrm{P}}$-dimensional generalized eigenvalue problem

[TABLE]

where $N_{\mathrm{P}}$ represents a sufficiently large number for convergence of the approximated Fermi function, and the elements of the matrices $\mathbf{A}$ and $\mathbf{B}$ can be defined as

[TABLE]

where $q,q^{\prime}=1,2,\ldots,N_{\mathrm{P}}$. Using the eigenvalue $\lambda_{p}$ and eigenvector $v_{p}$, $\tilde{z}_{p}$ and $\tilde{R}_{p}$ can be written as

[TABLE]

In summary, Eq. (9) can be rewritten as

[TABLE]

II.3 Liechtenstein formula for periodic images

In addition to the exchange coupling constant of the individual site representation in Eq. (19), it is possible to consider the exchange coupling constant $J_{ij}$ between the periodic images of $i$ and $j$, which is shown in Fig. 1(b). This can be derived by summing up $J_{i\mathbf{0},j\mathbf{R}}$ for all the considered cell indices $\mathbf{R}$. That is,

[TABLE]

It is then possible to reduce the integration variable to only $\mathbf{k}$ from Eq. (8), as follows:

[TABLE]

which is derived in Appendix C. This expression is useful when the unit cell is very large, because the interaction between the distant images of $i$ and $j$ [thin dashed lines in Fig. 1(b)] can be negligibly small compared with the interaction of the nearest images of $i$ and $j$ [thick solid line in Fig. 1(b)].

Hereafter, we refer to calculations of $J_{i\mathbf{0},j\mathbf{R}}$ using Eq. (19) as individual pair calculations and those of $J_{ij}$ according to Fig. 1(b) as periodic image calculations to distinguish them.

III Computational models and methods

Figure 3 shows the five systems examined in this paper. Each system consists of $N_{\mathrm{cell}}\times N_{\mathrm{cell}}\times N_{\mathrm{cell}}$ conventional unit cells of a bcc Fe crystal.

For first-principles calculation based on the localized basis set for the systems, we performed a density functional calculation using the OpenMX code T_Ozaki_2003 . Unless otherwise specified, we adopted the Perdew–Burke–Ernzerhof exchange-correlation functional GGA-PBE within the generalized gradient approximation (GGA-PBE). For the pseudoatomic orbital basis sets, we adopted the $s2p2d2$ basis set for Fe, where the 3$p$, 3$d$ and 4$s$ states of Fe are treated explicitly as valence states. The cutoff radius was set to 6.0 times the Bohr radius for Fe. We adopted the fully relativistic pseudopotentials generated by the Morrison–Bylander–Kleinman scheme MBK . We used an electronic temperature of 300 K, and the convergence criterion for the total energy was chosen as $1.0\times 10^{-6}$ Ha.

The numbers of real and reciprocal space grids are determined according to the system size, as shown in Table 1. We set the numbers of real and reciprocal space grids proportional to the system size and the inverse of the system size, respectively, so as to obtain strict quantitative consistency between the calculation results for various system sizes.

IV Results

IV.1 Dependence of coupling constants on the computational parameters $N_{\mathrm{P}}$ and $N_{k}$

After a self-consistent field (SCF) calculation is performed for each system, it is possible to obtain the Hamiltonian and overlap matrices corresponding to the converged electron density as an output of first-principles calculation. The eigenvalues $\varepsilon_{n}(\mathbf{k})$ and vectors $\mathbf{C}(\mathbf{k})$ needed for calculations using Eqs. (10) and (11) are obtained by solving the generalized eigenvalue problem represented by the Hamiltonian and overlap matrices. It is then possible to calculate the exchange coupling constants using Eq. (19). Here, the computational accuracy of the $J_{ij}$ calculation is determined mainly by the number of $k$-points, which is $N_{k}\times N_{k}\times N_{k}$ for calculations using Eqs. (10) and (11), and the number of poles $N_{P}$ in Eq. (19). To test the efficiency of the derived formalism in this paper, we examined the exchange coupling constant $J_{\mathrm{1NN}}$ between the first-nearest-neighbor sites in bcc Fe for various $N_{k}$ and $N_{\mathrm{P}}$ values. Note that the number of $k$-points in $J_{ij}$ is not necessarily the same as that in the SCF calculations. We employed a fixed value of $N_{k}$ for the SCF calculations, as shown in Table 1, and examined different values of $N_{k}$ for the $J_{ij}$ calculations.

Figure 4 shows the calculated $J_{\mathrm{1NN}}$ values as functions of $N_{k}$ and $N_{\mathrm{P}}$ for the $1\times 1\times 1$ system. As shown in Fig. 4(a), the $J_{\mathrm{1NN}}$ value oscillates within a range of approximately 0.1 meV even for a large $N_{k}$. By contrast, $J_{\mathrm{1NN}}$ increases monotonically and converges within a range of 0.05 meV at around $N_{\mathrm{P}}=60$, as shown in Fig. 4(b). These features indicate that the computational cost of $J_{ij}$ calculation scales better for the number of poles $N_{k}$ than for the number of $k$-points, $N_{k}\times N_{k}\times N_{k}$.

This result must also be compared with the computational results obtained by adopting the Matsubara poles Y_O_Kvashnin_2015 . As we discussed in Sec. II, the Matsubara approximation would result in slow convergence on the order of $N_{\mathrm{P}}$, and the number of necessary Matsubara poles needed to obtain an accuracy of $10^{-5}$ Ry (0.136 meV) is 1024 at 300 K. By contrast, Fig. 4(b) shows that the calculated $J_{\mathrm{1NN}}$ value already converges to within 0.05 meV at $N_{\mathrm{P}}$. Thus, computation can be speeded up by a few ten times by straightforward adoption of the finite pole approximation of the Fermi function.

IV.2 Optimal computational parameters for $J_{ij}$ calculation

To test the efficiency of our formalism for larger systems, we explore the optimal values of the computational parameters $N_{k}$ and $N_{\mathrm{P}}$ for various system sizes that give appropriate exchange coupling constants within an acceptable error tolerance.

For this purpose, $J_{\mathrm{1NN}}$ was again examined, and defined the acceptable error from the converged value as 0.05 meV for $J_{\mathrm{1NN}}$.

Table 2 shows the set of smallest values $N^{\mathrm{opt}}_{\mathrm{P}}$ and $N^{\mathrm{opt}}_{k}$ that give the appropriate $J_{\mathrm{1NN}}$ with an error tolerance of 0.05 meV for each system. We also show the computational times $t_{\mathrm{eig}}$ for eigenvalue calculation and the $t_{\mathrm{indiv}}$ value for the Liechtenstein calculation based on Eq. (19), together with the numbers of nodes and cores. All calculations are performed by Intel(R) Xeon(R) E5-2680 v3 processors, and we parallelize only the calculations of the $k$-points using MPI and not the energy integrations. We adopt BLAS routines for the matrix multiplications and general eigenvalue problems.

Table 2 shows a few remarkable features. The optimal value of $N_{\mathrm{P}}^{\mathrm{opt}}$ is almost independent of the system size, whereas the optimal value of $N^{\mathrm{opt}}_{k}$ decreases with increasing system size, and is roughly in inverse proportion to the system size.

It is also possible to compare the computational times $t_{\mathrm{eig}}$ required for the eigenvalue calculation with $t_{\mathrm{indiv}}$ for various system sizes. When the system is small, $t_{\mathrm{eig}}$ is comparable $t_{\mathrm{indiv}}$, whereas $t_{\mathrm{eig}}$ increases quickly with increasing system size. In contrast, $t_{\mathrm{indiv}}$ remains to the very small values even for large systems. This means that it is possible to calculate $J_{ij}$ with small calculation costs even for large systems when the full set of eigenvalues and vectors for the necessary $k$ points is given.

IV.3 Comparison of $J_{ij}$ calculations for individual sites and periodic images

Next, we compared the $J_{\mathrm{1NN}}$ values determined by Eq. (19), which we call $J_{\mathrm{1NN}}$ for an individual pair, $J_{\mathrm{1NN}}^{\mathrm{indiv}}$, and those determined by Eq. (21), which we call $J_{\mathrm{1NN}}$ for periodic images, $J_{\mathrm{1NN}}^{\mathrm{periodic}}$. For quantitative comparison, we first defined the absolute difference $\Delta J_{\mathrm{1NN}}$ as

[TABLE]

Table 3 shows $\Delta J_{\mathrm{1NN}}$ for various system sizes, together with the computational times of the periodic image calculations, $t_{\mathrm{periodic}}$. For consistency, $N_{k}$ and the computational conditions are set to the same values as those used in the individual pair calculations summarized in Table 2. We can see in Table 3 that $\Delta J_{\mathrm{1NN}}$ is much larger than $J_{\mathrm{1NN}}^{\mathrm{indiv}}$ itself in the $1\times 1\times 1$ system. This is because in the periodic image calculation, all the interactions of neighboring sites in the $1\times 1\times 1$ system are counted, and the calculated value is approximately 8 times larger than the realistic value. The difference $\Delta J_{\mathrm{1NN}}$ is much smaller for larger systems, whereas it exceeds the error tolerance of 0.05 meV. Since the difference of 0.1 meV in $\Delta J_{\mathrm{1NN}}$ roughly corresponds to the deviation of 6 K in the Curie temperature, the differences of 0.3–0.5 meV seen in Tab. 3 would result in errors of 18–31 K in the Curie temperature. It is thus possible to conclude that the approximation of exchange coupling between individual sites to that of periodic images gives no advantage in terms of computational cost.

It is also necessary to note that the computational times $t_{\mathrm{eig}}$ for large systems are much larger than $t_{\mathrm{indiv}}$ and $t_{\mathrm{periodic}}$ for large systems. It is thus necessary to store certain sets of eigenvalues and vectors to speed up the computation of the coupling constants for number of pairs.

IV.4 Distance dependence of $J_{ij}$ and comparison with plane wave calculation

Another important problem related to $J_{ij}$ calculation is that, we do not know what the most reliable $J_{ij}$ is. Particularly, it is necessary to clarify whether $J_{ij}$ calculation derived in this study agrees quantitatively with those in other formulations. For this purpose, we calculated the distance dependence of $J_{ij}$, and compared it with the Liechtenstein calculation result obtained from Akai-KKR.

Figure 5 shows $J_{ij}$ as a function of atomic distance $r_{ij}$ in a bcc-Fe crystal. The red symbols represents $J_{ij}$ based on our formalism, and the blue symbols represents $J_{ij}$ from Akai-KKR. Here, we adopted $N_{k}=30$ and $N_{\mathrm{P}}=60$ for $1\times 1\times 1$ system for the formalism in this study, and $N_{k}=20$, non-relativistic muffin-tin potential and GGA-PBE approximation in Akai-KKR calculation. It is possible to see in Fig. 5 that the profiles of two results agree well with each other.

Moreover, it is possible to derive the Curie temperature from the calculated results. We calculated the Curie temperature $T_{\mathrm{C}}$ of bcc-Fe based on mean field approximation using computational results of $J_{i\mathbf{0},j\mathbf{R}}$ and $J_{ij}$ of system $1\times 1\times 1$ with Eqs. (19) and (21), as the maximum eigenvalue of the matrix below:

[TABLE]

(For details, see Appendix D.) As a result, the Curie temperature of bcc Fe is calculated as 1322 K. Considering the well known fact that the mean field approximation overestimates the Curie temperature of bcc Fe by $\sim$ 40%, it is possible to say that our computational result is in a correspondence with the experimental value of 1043 K.

V Summary

In this paper, we derived an explicit form of the Liechtenstein formula within the localized orbital basis representation and developed a computational code for the output of the first-principles calculation code OpenMX. In the derivation, we adopted the finite pole approximation of the Fermi function, which simplified and speeded up the energy integrations. To test the efficiency and computational speed of our implementation, we calculated the exchange coupling constant $J_{\mathrm{1NN}}$ of the first-nearest-neighbor pairs in bcc Fe crystals with various system sizes. Using a new formalism based on the finite pole approximation, we were able to calculate $J_{\mathrm{1NN}}$ with an error tolerance of 0.05 meV with small computational times even for large systems. To obtain an efficient formalism for calculating $J_{ij}$ for large systems, we compared the values of $J_{\mathrm{1NN}}$ based on the finite pole approximation and those of periodic image calculations for various system sizes. It is shown that the approximation to periodic images gains almost no computational speedup, while the it gives the deviation of $J_{ij}$ at about 0.5 meV. We also calculated the dependence of $J_{ij}$ on the atomic distance $r_{ij}$, and compared it with that obtained by Akai-KKR calculation. It is shown that the two profiles agree well with each other, indicating the transferability of computational results derived by the two different formalisms.

Acknowledgements.

The authors acknowledge Hisazumi Akai, Sonju Kou and Shotaro Doi for fruitful discussions and valuable comments. This work was supported in part by MEXT, Japan, as a social and scientific priority issue CDMSI to be tackled by using the post-K computer, and the Elements Strategy Initiative Project under the auspices of MEXT, as well as KAKENHI Grant No. 17K04978. Some of the calculations were performed using the supercomputers at ISSP, The University of Tokyo, and TSUBAME, Tokyo Institute of Technology, as well as the K computer, RIKEN (Project Nos. H30-Cb-0009, H31-Ca-0025, hp170269, hp180206, and hp190169). One of the authors (MM)’s work in ISSP, Univ. of Tokyo is supported by Toyota Motor Corporation.

Appendix A Detailed derivation of non-orthogonal Liechtenstein formula

In this section, we present a detailed derivation of Eqs. (8) and (9) from the general form in Eq. (4) by applying the completeness relationships for the Bloch functions $|\mathbf{k},n,\sigma\rangle$,

[TABLE]

and the completeness relationship for non-orthogonal basis sets,

[TABLE]

where $\mathbf{S}^{-1}_{\mathrm{RS}}$ is the overlap matrix of the real space basis.

Applying Eq. (23) to Eq. (4), we get

[TABLE]

Although the equation above has four sets of $(\mathbf{k},n)$, it is possible to drop two of them because the Green’s function becomes diagonal in the eigenfunction representation, that is,

[TABLE]

The remaining task is to express the eigenfunction representations of the potential difference operators $\hat{P}_{i,\mathbf{0}}$ and $\hat{P}_{j,\mathbf{R}}$. This can be done by applying Eq. (24) as follows:

[TABLE]

Here we have to be careful that the localized basis representation of the potential difference operator, $\langle j^{(1)},\nu^{(1)},\mathbf{R}^{\prime(1)}|\hat{P}_{i,\mathbf{0}}|i^{(2)},\mu^{(2)},\mathbf{R}^{(2)}\rangle$, becomes nonzero only at $j^{(1)}=i^{(2)}=i$ and $\mathbf{R}^{\prime(1)}=\mathbf{R}^{(2)}=\mathbf{0}$. Equation (28) reduces to

[TABLE]

One can then replace the Bloch functions with the localized orbital representation in Eq. (6) to obtain

[TABLE]

Combining Eqs. (29), (30), and (31), we obtain

[TABLE]

with $[\hat{P}_{i}]_{\nu\mu}\equiv\langle i,\nu,\mathbf{0}|\hat{P}_{i,\mathbf{0}}|i,\mu,\mathbf{0}\rangle$. Similarly,

[TABLE]

with $[\hat{P}_{j}]_{\nu^{\prime}\mu^{\prime}}\equiv\langle j,\nu^{\prime},\mathbf{R}|\hat{P}_{j,\mathbf{R}}j,\mu^{\prime},\mathbf{R}\rangle$. The resulting equation can be written in explicit form as

[TABLE]

Because we treat periodic systems, the matrix representations of $\hat{P}_{i,\mathbf{0}}$ and $\hat{P}_{j,\mathbf{R}}$ by localized basis sets do not depend on the cell indices $\mathbf{0}$ and $\mathbf{R}$. Finally, we obtain Eq. (7) by replacing the summations of $\mathbf{k}$ and $\mathbf{k}^{\prime}$ with the corresponding integrals.

Appendix B Complex contour integration at the limit $\varepsilon_{\mathrm{max}}\rightarrow+\infty$

In this subsection, we consider the behavior of the integrand $\mathcal{F}(z)$ in Eq. (12) on the complex contour $C_{2}$ at the limit $\varepsilon_{\mathrm{max}}\rightarrow+\infty$. From the conditions, we can straightforwardly derive $|z|\rightarrow+\infty$ on $C_{2}$, and the energy eigenvalues and poles become negligible compared with $|z|$:

[TABLE]

Then it is possible to take the limit of the fractions in the Green’s functions and Fermi functions as follows:

[TABLE]

[TABLE]

where

[TABLE]

It is straightforward to derive $|\delta|,|\delta^{\prime}|\rightarrow 0$ from Eq. (35), and we can neglect the first- and higher-order terms of $\delta$ and $\delta^{\prime}$ in Eqs. (36) and (37). It is then possible to find that

[TABLE]

Consequently, the integration of $\mathcal{F}(z)$ on $C_{2}$ yields a value of zero.

Appendix C Summation for periodic images

In this section, we reduce the double integration of $\mathbf{k}$ in Eq. 8 to a single integration by taking the summation of periodic images shown in Fig. 1(b). From the definition in Eq. (20),

[TABLE]

Considering the relationship

[TABLE]

we obtain Eq. (21).

Appendix D Curie temperature of periodic systems

Using the mean field approximation, Curie temperature of general periodic system can be obtained by finding the maximum value of $T$ which satisfies the following relationship:

[TABLE]

Using the definition in Eq. (21), the relationship can be rewritten as

[TABLE]

That is, because $J_{i\mathbf{0},i\mathbf{0}}$ derived by Liechtenstein formula need not to become zero but doesn’t have physical meaning, it is necessary to subtract its contribution from the periodic sum.