Exchange interactions from a nonorthogonal basis set: from bulk ferromagnets to the magnetism in low-dimensional graphene systems

We present a computational method to determine the exchange constants in isotropic spin models. The method uses the Hamiltonian and overlap matrices computed from density functional schemes that are based on nonorthogonal basis sets. We demonstrate that the new method as implemented in the SIESTA code reproduces the Heisenberg interactions of simple metallic bulk ferromagnets as obtained from former well--established computational approaches. Then we address $sp$ magnetism in graphene nanostructures. For fluorinated graphene we obtain exchange interactions in fairly good agreement with previous calculations using maximally localized Wannier functions and we confirm the theoretical prediction of a 120$^\circ$ N\'eel state. Associated with the magnetic edge-states of a zigzag graphene nanoribbon we find rapidly decaying exchange interactions, however, with an unconventional distance dependence of $\exp(-\sqrt{r/\delta})$. We show that the stiffness constant derived from the exchange interactions is consistent with previous estimate based on total energy differences of twisted spin configurations. We highlight that our method is an efficient tool for the analysis of novel hybrid nano-structures where metallic and organic components are integrated to form exotic magnetic patterns.

We present a computational method to determine the exchange constants in isotropic spin models. The method uses the Hamiltonian and overlap matrices computed from density functional schemes that are based on nonorthogonal basis sets. We demonstrate that the new method as implemeted in the SIESTA code reproduces the Heisenberg interactions of simple metallic bulk ferromagnets as obtained from former well-established computational approaches. Then we address sp magnetism in graphene nanostructures. For fluorinated graphene we obtain exchange interactions in fairly good agreement with previous calculations using maximally localized Wannier functions and we confirm the theoretical prediction of a 120 • Néel state. Associated with the magnetic edge-states of a zigzag graphene nanoribbon we find rapidly decaying exchange interactions, however, with an unconventional distance dependence of exp(− r/δ). We show that the stiffness constant derived from the exchange interactions is consistent with previous estimate based on total energy differences of twisted spin configurations. We highlight that our method is an efficient tool for the analysis of novel hybrid nano-structures where metallic and organic components are integrated to form exotic magnetic patterns.

I. INTRODUCTION
Heisenberg-like spin Hamiltonians form a solid basis for describing the ground state and thermal behavior of a wide range of magnetic systems, either characterized by itinerant electrons or by local moments. The exchange constants entering the spin model can be derived from first principles. One of the most frequently used approaches is based on the magnetic force theorem originally introduced by Liechtenstein and co-workers [1], referred to as the LKAG formalism as what follows. Related methods have been developed since then to tackle correlated systems [2,3], relativistic effects [4,5] or both of them [2,6].
The use of one-electron Green's functions is an integral part of the formalism in Ref. [1]. Therefore, the Korringa-Kohn-Rostoker Green's function (KKR-GF) [7] and the tight-binding linear muffin-tin orbital (TB-LMTO) methods [8,9] have been particularly successful for calculating magnetic exchange interactions for bulk materials, surfaces, interfaces, films, superlattices, and even for finite metallic clusters [10][11][12]. Furthermore, the calculation of tensorial exchange interactions, including two-ion magnetic anisotropy parameters and Dzyaloshinskyi-Moriya interactions, has become available by extending the LKAG formula to relativistic systems [4,5]. This extension has opened the door to the analysis, design and tuning of complex magnetic states like domain walls [13], spin spirals [14,15] and magnetic skyrmions [16][17][18][19] in ultrathin films. This extension also enables us to study the recently discovered van der Waals ferromagnets [20][21][22], whose magnetic state is stabilized by the anisotropy barriers that overcome the thermal spin fluctuations standing behind the Mermin-Wagner theorem [23].
The KKR-GF or TB-LMTO methods are not fully adapted to describe open systems such as atoms or molecules deposited on surfaces or suspended in nanoscale junctions because they commonly make use of the atomic sphere approximation (ASA). Such systems are accurately treated by other methods like the program packages VASP [24], Quantum Espresso [25], SIESTA [26] or ADF [27]. The first two expand the eigenstates onto a plane-wave basis [24,25], so a transformation to maximally localized Wannier functions [28] is needed to determine an orthogonal tight-binding basis set requested by the LKAG formalism [29]. In contrast, SIESTA and ADF conveniently use a basis set of wave functions that are localized on each individual atom, possess orbital quantum numbers and are nonorthogonal, hence producing selfconsistent Hamiltonians and overlap matrices that are already written in the tight-binding language. In order to calculate exchange constants in the spirit of the magnetic force theorem orthogonal basis sets, a generalization of the LKAG formula to nonorthogonal bases is called for. This is the main goal and accomplishment of the present article. We therefore develop a formalism which enables us to use some popular density-functional codes to evaluate exchange interactions of isotropic Heisenberg models in a broad range of magnetic systems. We implement the new scheme in the SIESTA code and present results for conventional bulk metallic magnets as well as for graphene ribbons and fluorinated graphene sheets displaying sp magnetism that compare well with those available in the literature. We outline the main features of our approach in Section II, then we present and discuss our results in Section III. A short summary closes the article. All the algebraic details of the approach can be found in the appendices.

II. METHOD
The classical Heisenberg model describes a lattice of localized classical spins characterized by unit vectors e i , where i denotes lattice (or atomic) sites. Within a nonrelativistic theory, the spin-spin interactions are encapsulated in terms of isotropic exchange parameters J ij entering the spin Hamiltonian, The magnetic force theorem enables us to extract the exchange parameters from the effective single-particle HamiltonianĤ that results from ab initio calculations [1,4]. We use a collinear-spin reference frame and writê H in terms of a basis set of localized orbitals centered at lattice sites. Consequently, the tight-binding Hamiltonian matrix H and the overlap matrix S are diagonal in the spin indices. Collecting all the basis functions assigned to a given site i, the corresponding spin-dependent and site-indexed blocks of the Hamiltonian are denoted by H σ ij and, similarly, S ij for the overlap matrix which is independent on the spin index since in practice for both spin channels the same basis functions are considered. We then find that the exchange parameters can be derived from the expression where ε F is the Fermi energy, Tr L denotes the trace of matrices in orbital space, is the appropriate site off-diagonal block of the matrix of expansion coefficients of the resolvent operatorĜ(z) = zÎ −Ĥ −1 with z = ε + iδ. A detailed derivation of Eq. (2) is given in the appendices. These expressions are a generalization of the seminal work of Liechtenstein et al. [1] to the case of a nonorthogonal tight-binding basis set.
We implemented the above equations by using the selfconsistent Hamiltonian and overlap matrices provided by the SIESTA code [30]. This can be achieved with the assistance of the sisl tool [31]. We devote the next section to validate our approach by giving three examples [32] by comparing the results of our proposed methodology with previous calculations. Moreover, we give a detailed description of the sp magnetism in low-dimensional graphene systems in terms of exchange interactions and analyze their asymptotic behavior.

Bulk ferromagnets
In this section we present the exchange interactions of selected bulk ferromagnets using our proposed approach and compare our results with former ones obtained from the screened KKR (SKKR) method [7] in the framework of the atomic sphere approximation. We considered ferromagnetic bcc Fe, hcp Co as well as fcc Ni. For both the SIESTA and SKKR calculations we used the same approximations for the exchange-correlation density functional and the same geometrical parameters. In case of Fe the generalized gradient approximation (GGA) as parameterized by the PBE scheme [33], while for Co and Ni the local spin density approximation (LSDA) [34] were employed. The lattice constants 2.87 Å, 2.50 Å and 2.49 Å were chosen for bcc Fe, hcp Co and fcc Ni, respectively. In addition, for Co we considered the ratio of c/a = 1.633 of an ideal hcp structure. scheme. However, we found that choosing a k-space cutoff of 50 Å −1 and a real space mesh cutoff of at least 500 Ryd ensured reliable accuracy for the SIESTA results. We noticed, however, that the choice of the pseudo-potential parameters had a considerable impact on the results for the ground state obtained from SIESTA and, subsequently, also on the calculated exchange parameters. In our calculations we used the pseudo-potential generation scheme described in Ref. [35]. The calculated spin magnetic moments of the three bulk ferromagnets are summarized in Table I for the two self-consistent schemes. The data in this table show an almost perfect agreement between the spin moments obtained from the two ab initio methods for bcc Fe, and relative differences of about 3 % and 8 % for hcp Co and for fcc Ni, respectively.
Next we calculated the isotropic exchange parameters for the three bulk ferromagnets by using the relativistic torque method within the SKKR [4] and via the formula in Eq. (2) with the tight-binding Hamiltonian and overlap matrices obtained from SIESTA. It should be mentioned that in the latter case we needed a 100 k-point mesh in each direction of the full Brillouin zone to ensure ad- The exchange parameters J ij obtained from the two methods are plotted in Fig. 1. For a more extensive comparison we also included the corresponding values reported in Ref. [11] in terms of the TB-LMTO ap-proach. Note that the spin model considered in Ref. [11] misses the factor of 1/2 in Eq. (1), therefore, the exchange interactions presented there should be multiplied by a factor of 2 in order to compare with those calculated from Eq. (2). Apparently, the three methods provide isotropic exchange interactions in remarkably good agreement with each other for all three bulk ferromagnets. Considering mainly the large ferromagnetic nearest neighbor interactions, but also in case of some farther couplings, the SIESTA and TB-LMTO values compare more precisely than those and the SKKR values, which is not surprising as the former two methods rely on the tight-binding scheme. As can be seen in Fig. 1(b), the exchange interactions derived from the SIESTA and SKKR calculations also compare remarkably well with those obtained from the TB-LMTO method for fcc Co.
The Curie temperature T C of ferromagnetic materials is one of the measurable quantities closely related to the exchange interactions. While the transition temperature can accurately be obtained from Monte Carlo or spin-dynamics simulations, here we present theoretical estimates based on the mean field approach which is extracted from the spin model parameters J ij as with the Boltzmann constant k B . We calculated T MFA C summing up the exchange parameters up to a distance of r 0j = 10 Å for hcp Co and fcc Ni, while r 0j = 25 Å for bcc Fe, reducing the numerical error of the results below 20 K. The data obtained within the SIESTA and SKKR methods shown in Table II are in fairly good agreement with each other and with those reported in Ref. [11], also presented in Table II. The somewhat large deviation of T MFA C of Co within the TB-LMTO method from the very similar values obtained using the SIESTA and SKKR codes can mainly be attributed to the different crystal structures used in these calculations. The meanfield approximation is known to overestimate the exact transition temperatures, which might explain the higher values of T MFA C as compared with the experimental T C in case of Fe and Co. The considerably lower mean-field estimates for the Curie temperature with respect to the experimental value in case of Ni is most possibly the consequence of the highly itinerant nature of the magnetism of bulk Ni [36,37].

Fluorinated graphene
We turn now to sp magnetism in the context of graphene. First we present results for the exchange interactions in single-side fluorinated graphene C 2 F and compare them to earlier calculations by Rudenko et al. [29], who used a maximally localized Wannier function basis [28] which was mapped from a plane wave basis [42]. Wannier orbitals form an orthonormal basis representation, thus Eq. (A28) can simply be evaluated using Jij for a fluorinated graphene sample calculated using a Wannier basis [29] (triangles) and extracted from a nonorthogonal basis using the SIESTA code. For both calculations a rowwise antiferromagnetic spin configuration was used as reference. Circles and crosses represent couplings between spins of parallel and anti-parallel orientations, respectively. The inset shows the atomic arrangement of the carbon atoms in this system. Here, black circles label carbon atoms bound to a fluorine atom. Arrows in the white circles represent the orientation of local moments, displaying a row-wise antiferromagnetic state. Color-coded disks denote the relative magnitude and sign (red positive, blue negative) of the exchange couplings between the central site and other magnetic carbon sites. (b) Exchange interactions Jij calculated from a nonorthogonal basis set using the SIESTA code and a ferromagnetic state as reference. The inset demonstrates that in this case the exchange parameters respect the symmetry of the underlying lattice. the corresponding matrices with respect to this representation. However, Wannier orbitals are not necessarily localized to a single atom. Therefore, local degrees of freedom like the atomic spin can not be unambiguously described in terms of a Wannier basis. In our approach, every nonorthogonal orbital is explicitly localized to a given atom in the system. The nonorthogonality of these orbitals can be handled by using appropriate local projection operators, as discussed in Appendix B. Hence, these orbitals describe properly atomic degrees of freedom. In Ref. [29] it was found that a row-wise antiferromagnetic (AFM) spin alignment is preferred with respect to the ferromagnetic (FM) state. In our self-consistent calculations performed with SIESTA we also considered a row-wise AFM configuration. For better comparison, we used the exchange-correlation functional and geometry parameters of Ref. [29]. Note that the F atoms are placed above the C atoms in only one of the two sublattices of graphene (say, in sublattice A), forming thus a triangular lattice. We found that the carbon atoms at sublattice B have a total magnetic moment of 0.736 µ B , with a contribution of 0.65 µ B coming from their p z orbitals and that the carbon atoms at sublattice A have negligible magnetic moments. This is in good agreement with the results of Ref. [29], where considerable spin polarization was found only for the p z type Wannier orbitals associated with the carbon sites at the B sublattice, with a magnitude of 0.59 µ B .
Choosing the row-wise AFM configuration as a reference, we calculated the exchange parameters of C 2 F by using Eq. (2) withG σ ij -s evaluated on a k-mesh of 200 points in each direction of the two-dimensional Brillouin zone. The resulting exchange interactions are shown in Fig. 2(a). Note that we label the interactions between the moments of the same and opposite orientations with different symbols. Since a row-wise AFM spin configuration does not respect the point symmetry of the triangular lattice, these two sets of interactions significantly differ from each other: the first nearest neighbor interactions between moments with the same orientation are much stronger antiferromagnetic than those between opposite moments. Notably, the exchange interactions reported in Ref. [29], labeled by triangles in Fig. 2(a), show a good agreement with those calculated by SIESTA between mo- ments of the same orientation at positions along a row of the row-wise AFM pattern. Moreover, farther couplings decrease exponentially as shown in Fig. 3(a). This can be understood as a consequence of the quasiparticle gap that is manifested in the total density of states shown in Fig. 3(c).
As mentioned above, the chosen reference state and, consequently, the calculated exchange parameters do not respect the C 3v symmetry of the underlying lattice. For this reason, we also chose the ferromagnetic state shown as an alternative reference state, and we re-calculated the exchange parameters. The results are shown in Fig. 2(b). We found that the exchange parameters are now consistent with the C 3v symmetry of the lattice and that they are characterized by large antiferromagnetic first nearest neighbor and ferromagnetic next-nearest neighbor couplings. The decay of the exchange interactions shows the expected power law ∝ r −2 [43] as depicted in Fig. 3 and is corroborated by the absence of a gap in the density of states that we show in Fig. 3(d).
It should be mentioned that both the first nearest neighbor antiferromagnetic and second nearest neighbor ferromagnetic couplings are consistent with the 120 • Néel ground state suggested in Ref. [29]. We evaluated the Fourier transform J( q) of the exchange interactions obtained from the ferromagnetic reference state and plotted it in Fig. 4 along the high symmetry directions of the hexagonal Brillouin zone. Indeed we obtain a clear maximum at the K point, which indicates the 120 • Néel state as ground state.

Graphene ribbon
In this section we analyse the one-dimensional itinerant ferromagnetic state that arises at the edge of zigzag graphene ribbons passivated by hydrogen atoms [44][45][46][47][48][49]. This rather elusive magnetic state has provoked enormous expectations in the past fifteen years, because of their plausible potential for spintronics applications [50,51].
We considered a hydrogen-passivated 28-carbon-atom wide graphene ribbon that extends along the x direction, as depicted in Fig. 5. In the self-consistent calculations we used the LSDA [34], we set a mesh cutoff of 200 Ry for the real-space integrals and we selected 100 k-points in the one-dimensional (1D) Brillouin-zone. The magnetic configuration was set ferromagnetic along the x direction, that means the 1D unit cell consisted of 28 carbon atoms and 2 H atoms. The magnetic moments yielded by this choice are indicated in Fig. 5. We find that the magnetic moments of the A and B carbon atoms are opposite in sign, implying that the ribbon's ground state displays an antiferromagnetic alignment with respect to the two edges. A total magnetic moment of 0.3 µ B per 1D unit cell is associated with each edge, that we obtain by adding the individual magnetic moments of the ribbon atoms from one edge to the center of the ribbon. The two sublattices A and B contribute 0.36 µ B and −0.06 µ B to this magnetic moment, respectively. These findings are in good agreement with previous calculations [52,53]. We calculated the exchange parameters J ij in the ribbon using the LKAG formula Eq. (2) with the groundstate Hamiltonian delivered by SIESTA. We found that the leading interactions occur between those edge carbon atoms on sublattice A that have the largest magnetic moments. Fig. 5 demonstrates that these interactions are ferromagnetic and fairly short-ranged. The decay of the interactions is non-monotonic: a small oscillatory behavior can be observed, but the interactions remain ferromagnetic for all distances. Our calculations also indicate non-negligible antiferromagnetic couplings between the A atoms at the edge and the first nearest neighbor atoms at the B sublattice.
In order to make a conceptional connection to the calculation presented in [52], we introduce meta magnetic moments as large as 0.3 µ B associated with one half of the graphene ribbon. To calculate the interactions between these meta moments we evaluate Eq. (2) considering all atomic orbitals in one half of the 1D unit cell of the ribbon containing 14 carbon atoms and one hydrogen atom. In Fig. 6 we also plot these exchange interactions as a function of the distance along the x direction by orange triangles. The nearest-neighbor interaction between the meta moments is somewhat reduced, while the farther interactions are enhanced compared to the corresponding interactions between the edge atoms, as can be seen in Fig. 6. Consequently, the interactions between the meta moments show a fairly monotonic decay, with similar characteristics as in case of the edge atoms for distances above 15 Å. We found that in this region the interactions fit well to a function ∝ e − √ r/δ , also shown in Fig. 6. In order to better visualize this unconventional decay we plotted the exchange parameters on a logarithmic scale against the square root of the distance in the inset of Fig. 6. Both sets of interactions display a nearly linear dependence on this graph, though the aforementioned oscillations for the edge atoms clearly show up. The experimentally accessible magnon spectrum E(q) of a ferromagnetic system is related to the Fourier transform of the exchange constants [11], where M is the magnitude of the magnetic moment per periodic unit. The calculated curves for J(0) − J(q) for the two considered spin models are plotted in Fig. 7, clearly proving that the ferromagnetic state (q = 0) is the ground state of the half ribbon system. The low-energy magnon spectrum of the ribbon was estimated in Ref.
where the index j runs over all considered magnetic atoms or 1D unit cells in the system other than the atom indexed by 0 and the x component of the displacement vector between the sites 0 and j is denoted by x j . Considering only the edge atoms in the sum we obtain D edge = 2308 meVÅ 2 , while considering the half ribbon meta magnetic moments we get D half = 3406 meVÅ 2 .
Though these values are larger than the one reported in Ref. [52], they are reasonably similar in magnitude and support the observation that sp magnets might have higher spin stiffness than the conventional d ferromagnets. Based on our calculations we conclude that the high apparent value of the spin stiffness is caused by the almost monotonic, unconventional decay of magnetic correlations along the ribbon edge.

IV. SUMMARY
We presented a computational approach that determines the exchange parameters of isotropic spin models based on the magnetic force theorem, directly from ab initio calculations using a nonorthogonal basis set to expand the eigenstates of the system. We demonstrated that the new method accurately reproduces the Heisenberg interactions of simple metallic bulk ferromagnets delivered by well-established computational approaches. We studied the magnetism of two systems based on graphene. For fluorinated graphene we obtained exchange constants in fairly good agreement with previous calculations using maximally localized Wannier functions and we confirmed the theoretical prediction of a 120 • Néel state. The long-range behavior of the exchange interactions was found consistent with the electron spectrum of the system around the Fermi level. For zigzag graphene nanoribbons we found that the stiffness constant derived from the exchange constants is consistent with previous estimates based on total energy differences of twisted spin configurations. We also found an unconventional exp(− r/δ)-like decay of the interaction. Understanding this exotic behavior poses a challenge for further investigations. The SIESTA code can easily handle large nano-scale systems of high chemical complexity, therefore we are convinced that the presented method is a very efficient tool for the analysis and design of novel hybrid nano-structures hosting exotic magnetic patterns.

Appendix A: Exchange interactions in an orthogonal basis
In this section we give a detailed derivation of the LKAG formula for J ij [1,2] to be evaluated by density functional theoretical calculations performed in an orthogonal tight binding basis. Andersen's force theorem [54] states that if the system is in its ground state then the change of the total energy E tot due to a small variation in the external potential can be directly linked to small variations in the Kohn-Sham energy E KS calculated at fixed density without the need for further selfconsistent calculations. The force theorem thus provides us with a computationally inexpensive way to obtain response functions.
Neglecting relativistic effects and longitudinal spin fluctuations, the energy of a spin system is usually mapped to a classical Heisenberg model, Consider now a ferromagnetic ground state where all spins point in the same direction e 0 that has the energy If a single spin located at site i is excited to e i = e 0 , the energy of this single-spin excitation is given by The energy of a two-site excitation, e i = e 0 and e j = e 0 (i = j), can be expressed as The interaction energy between spins i and j, E int ij , is then defined by with δ e i/j = e i/j − e 0 . Calculating the energy cost of appropriate local perturbations we can thus extract the classical J ij parameters from ab initio Green's function methods. Applying Lloyd's formula [55] in the spirit of the force theorem the energy cost of a perturbation δV can be cast in terms of the Green's operator (resolvent) where ε F is the Fermi energy,Î is the identity operator andĤ corresponds to a Hamiltonian which selfconsistently determines E tot . Assuming now that δV i and δV j are operators that describe the local perturbations corresponding to spin rotations at sites i and j, respectively, and using the identitŷ where the scattering operatorT i is defined aŝ (A11) In the spirit of Eq. (A5), the interaction energy for the two-site perturbation is then given by Since we are interested in small perturbations around the ground state, we can safely use the Born approximation, T i ≈ δV i , and we also expand the logarithm as ln(1−x) ≈ −x, thus Eq. (A12) reduces to Note that so far we have not considered anything specific about the perturbation operatorsV i/j . Within the tight binding (TB) scheme a matrix representation ofĤ is used in an orthogonal basis of localized atomic-like wavefunctions centered at sites of the lattice. Thus, the basis functions are labeled by lattice sites n, composite angular momentum indices L = ( , m) and the spin index s = ±1/2 (or ↑ and ↓). As what follows, we will note matrices of the entire site-angular momentumspin space with bold-face letters, double and single underlines will denote block matrices in common angular momentum-spin space and in only angular momentum space, respectively: with the Pauli matrices σ. For simplicity, we shall assume that H nL,n L = 0 for n = n (A16) and spin dependence applies only to the site-diagonal blocks of the Hamiltonian, When the spin is aligned parallel to the z axis, the form of the local Hamiltonian is thus In case of a ferromagnetic (in general, collinear) magnetic configuration of the host with a magnetic orientation e 0 , where I denotes the unit matrix in spin space. The corresponding matrix representation of the Green's function is of the same form, According to Eq. (A17) the change of the Hamiltonian due to local spin rotations is given by elements defined as where H s ii denotes the angular momentum representation of the spin-dependent part of the Hamiltonian confined to the site i. In order to calculate the interaction energy of two spins in Eq. (A13) we evaluate the trace by substituting Eqs. (A21) and (A22): where Tr Ls denotes the trace of a matrix in both angular momentum and spin space. Using the algebraic properties of Pauli matrices the traces can easily be evaluated in spin space yielding where Tr L denotes trace in angular momentum space only. For infinitesimal rotations δ e i ⊥ e 0 , therefore, the second term will be neglected.
that indeed cancels the third contribution to (A24). Thus, the interaction of two spins can indeed be written as with Rewriting G c/s ij in terms of G ↑/↓ ij the above expression can be reduced to This expression is the celebrated LKAG formula [1,2]. It is important to note that if the magnetic orientations at site i and j are opposite in sign, as happens in case of collinear antiferromagnetic configurations, the J ij obtained from Eq. (A28) should be changed in sign as the spin channels at the two sites are reversed with respect to each other.

Appendix B: Some identities in a nonorthogonal basis
Here we review useful identities related to nonorthogonal bases, some of them discussed in Ref. [56]. Using these identities we then generalize Eq. (A28) to nonorthogonal bases.
A basis formed by states {|i } is not orthogonal if its elements have finite overlap In practice real valued basis functions are chosen, therefore the overlap matrix is symmetric. The inverse of the overlap matrix S defines the dual basis ĩ as whose elements are orthogonal to the original basis, The expansion of a general operatorÂ with respect to basis {|i } is defined aŝ while the matrix elements in the original basis can be expressed as Obviously, the expansion coefficientsÃ ij are the matrix elements of the operator in the dual basis, As what follows we shall denote the matrix of an op-eratorÂ with respect to the nonorthogonal basis with A, while the matrix in the dual basis will be denoted bỹ A. Note that these two matrices are connected by the overlap matrix S as The trace of an operator is calculated with the help of an orthogonal basis {|α }, Next we consider matrix elements and traces of operator products. The trace of a simple product gives which generalizes to The matrix element of a simple product is expressed as that can be generalized to Using the Taylor expansion of an operator function, In this section we discuss a pragmatic approximation to treat local spin rotations in a nonorthogonal basis leading to the generalization of the formula for J ij (A28) derived for orthogonal basis. Restricting our discussion to collinear magnetic systems, it is natural to choose a basis where the site and orbital degrees of freedom form the nonorthogonal part of the basis, while the basis functions are eigenvectors of the spin operator projected to the orientation of the magnetization. That is we consider the basis |σ ⊗ |iL with the property ( σ| ⊗ iL|) (|σ ⊗ |jL ) = σ |σ iL |jL (C1) = δ σ,σ S iL,jL , where i and j denote lattice sites, L and L stand for orbital degrees of freedom and σ, σ label the eigenvectors of the spin operator.
Let us define the local perturbation operator as whereĤ is a Hamiltonian whose matrix elements have been calculated by some self-consistent scheme,Ô G describes a global rotation of the spin degrees of freedom around direction n with angle φ andP i is a projector built up from all orbital degrees of freedom associated with site i: The identity operatorsÎ L andÎ S act on all orbital degrees of all atomic positions and in spin space, respectively. Note that this direct Hermitian projection does not project to a subspace with integer dimension [56]. The operator δV i has the convenient property that its matrix elements are only finite between orbitals located at site i, and are equal to the matrix elements of the Hamiltonian rotated globally in spin space relative to the reference Hamiltonian, δĤ =Ô † GĤÔ G −Ĥ. Since the global spin rotation and local projection act independently in the local perturbation, Eq. (C2), the evaluation of formula (B21) follows the steps as for the orthogonal basis. Thus, the expression of the Liechtenstein formula is readily generalized to nonorthogonal bases: with the actual expressions of the above matrices in the nonorthogonal basis.