TB2J: A python package for computing magnetic interaction parameters

TB2J: a python package for computing magnetic interaction parameters Xu Hea,b,c,∗, Nicole Helbigc, Matthieu J. Verstraeteb,c, Eric Bousqueta aPhysique Th´ eorique des Mat´ eriaux, Q-MAT, CESAM, Universit´ e de Li´ ege, B-4000 Sart-Tilman, Belgium bCatalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC, BIST, Campus UAB, Bellaterra, Barcelona, 08193, Spain cNanomat, Q-Mat, CESAM, and European Theoretical Spectroscopy Facility, Universit´ e de Li` ege, B-4000 Li` ege, Belgium Abstract We present TB2J, a Python package for the automatic computation of magnetic interactions, including exchange and Dzyaloshinskii-Moriya, between atoms of magnetic crystals from the results of density functional calculations. The program is based on the Green’s function method with the local rigid spin rotation treated as a perturbation. As input, the package uses the output of either Wannier90, which is interfaced with many density functional theory packages, or of codes based on localised orbitals. A minimal user input is needed, which allows for easy integration into highthroughput workflows. 1. Introduction First-principles simulations of magnetic materials have attracted strong interest in the last decade due to an increase in precision provided by density functional theory (DFT) codes for strongly correlated magnetic atoms and also because of the developments in spintronics applications [1,2,3,4,5,6]. Understanding the complex microscopic origin of the magnetic interactions often requires a reduction of the full manybody electronic interactions to effective Hamiltonians, the two most important ones being the Heisenberg and the Hubbard Hamiltonians [7]. The parameters of these Hamiltonians are fit to DFT data to give access to an easier understanding of the magnetic interactions and to allow for simulations of larger systems and/or with dynamics. Novel sensing, storing and computing technologies have been proposed using spin waves [8] and skyrmions [9,10]. Theoretical models have consistently opened new vistas and explained experimental findings in complex new chemistries, geometries, and heterostructures. Reliable ab initio calculations of the effective Hamiltonian parameters is crucial for predictive and quantitative simulations. Fitting these parameters is, however, often very cumbersome and necessitates a case-by-case construction [11]. For the Heisenberg Hamiltonian, one of the most common fitting procedures is through total energy cal- ∗Corresponding author Email address: [email protected] (Xu He) culations or energy mapping analysis [11]. This approach necessitates the calculation of the total energies of different magnetic configurations. The parameters of the Heisenberg Hamiltonian are then fit to these energies under the supposition that the change in energy is only related to the magnetic interactions. This requires to have at least as many calculated magnetic configurations as the number of parameters of the Hamiltonian (though more magnetic configurations are usually necessary to converge the fit), which often requires the use of large supercells. This simple method gives good results in some cases but the Heisenberg model itself can break down if the chosen magnetic configurations deviate too much from the ground-state one. This can happen, e.g., if a magnetic phase, often the ferromagnetic (FM) one, closes the band gap of an insulating antiferromagnetic (AFM) ground-state. Also, the method becomes unsuitable if the supercell needed to probe all of the pertinent magnetic configurations is too large to be handled by DFT. For example, if certain configurations show a delocalized picture for the electrons they must be excluded from the fit, which typically leads to an increase in the size of the supercell needed to provide a sufficiently large number of configurations. Another method to determine the parameters is the Generalized Bloch Theorem (GBT) [12,13] which adapts the boundary conditions on the spin orientation to account for a spin spiral. From the change in total energy one can extract the magnetic interaction parameters. As yet another alternative, Density Functional Perturbation Theory (DFPT) has been extended to magnetic Preprint submitted to Elsevier September 8, 2020 arXiv:2009.01910v2 [cond-mat.mtrl-sci] 7 Sep 2020 field perturbations by Savrasov [14] and yields the susceptibility and the magnetic exchange as a sub-product. This perturbation has been implemented in abinit [15] and quantum Espresso [16], and has also been combined with atomic displacements within non-collinear formalisms [17]. For several magnetic ions in the unit cell, a local perturbation DFPT scheme is presented in Ref. [18]. A different procedure, which avoids the problems of the total energy mapping, employs Green’s functions by taking the local spin rotation as a perturbation, as proposed in the seminal work of Liechtenstein, Katsnelson, Antropov and Gubanov (LKAG) [19]. The Green’s function method allows to determine the Heisenberg magnetic interaction parameters from the ground-state solution of the system, regardless of whether it is FM or AFM. This method also allows access to a bandby-band decomposition of the different magnetic interactions, and it is easier to automatize than the total energy mapping [20]. The method was also extended to correlated systems in Refs. [21,22]. By taking into account relativistic effects (spin-orbit coupling), the Dzyaloshinskii-Moriya interaction (DMI) [23,24,25,26] and the magnetic anisotropy can be calculated as well [21,22]. The method has been extended to orbital-spin and orbital-orbital magnetic interactions [27,28]. Higher order terms in the Hamiltonian, like the four-spin interaction or the biquadratic term, were also calculated through this method [29,30]. The LKAG method of computing magnetic interactions has been initially implemented with the Korringa-KohnRostoker Greens function (KKR-Green) [31,32] and tight-binding linear muffin-tin methods [33], and has proven to be quite efficient. Its extension to other basis sets popular for DFT calculations will strongly broaden its usage. Constructing the Green’s function has proven to be cumbersome for non-localized basis sets (e.g. plane waves), but this difficulty is greatly reduced by the use of Wannier functions (WF) as done by Kokorin et al. in [34]. The WFs can be constructed from first principles with the widely used open source code Wannier90 [35,36]. Wannier90 has been interfaced with a large number of DFT codes, including Abinit [37,38], Quantum Espresso [39], Siesta [40], VASP [41,42], Wien2K [43], Fleur [44], Octopus [45,46], OpenMX [47], GPAW [48], ELK [49], and many others. Several codes, including exchange.x [34], nojij [50], and Jx [51], exist which calculate magnetic exchange parameters from Wannier functions or linear combination of atomic orbitals (LCAO) DFT results. However, to the best of our knowledge, DMI and anisotropic exchange from non-relativistic effects have not been integrated yet. Here, we present a Python package, TB2J, that allows for automatic and systematic calculations of the parameters of a Heisenberg Hamiltonian through the Korotin approach [34]. The script can compute the isotropic exchanges, the anisotropic exchanges, and the DMI from the output of the Wannier90 code or from a LCAO Hamiltonian (Siesta and OpenMX through the SISL package [52], and GPAW) following the scheme proposed in Refs. [19,24,30]. TB2J is designed with the goal to minimize the number of inputs and actions from the user. In most cases, only the paths of the Wannier or DFT related files and the species of magnetic atoms are mandatory. Several types of output files are generated which can then be used in spin dynamics and Monte Carlo codes. TB2J’s API is implemented in an abstract manner that eases interfacing with new tightbinding-like Hamiltonians or spin dynamics simulation. We first present the general methods used in TB2J and then describe how these methods are implemented. Next, we exemplify the usage of TB2J by applying it to, BCC Fe, HCP Co bulk, SrMnO3, BiFeO3, and La2CuO4 crystals. In the end we discuss the advantages and limitations of this method. Throughout this paper, vectors are denoted with the~ notation while matrices are represented with bold characters. 2. Formalism and algorithms The main idea of the method is to perturb a localized spin, in both the Heisenberg spin model and the DFT electronic model. For the latter, the rigid spin-rotation perturbation is done within the single-particle Green’s function method. The Heisenberg parameters are then mapped to the electron expressions [19]. We will consider the quantization axis along zin the following. 2.1. Heisenberg model The Heisenberg Hamiltonian contains four different parts and reads as E=−∑ i Ki(~ Si·~ei)2 −∑ i,jJiso i j ~ Si·~ Sj +~ SiJani i j ~ Sj +~ Di j ·~ Si×~ Sj,(1) 2 Convention ˜ Ji j −1 2∑ i j ˜ Ji j~ Si·~ Sj2Ji j −∑ <i j> ˜ Ji j~ Si·~ Sj2Ji j 1 2∑ i j ˜ Ji j~ Si·~ Sj−2Ji j ∑ <i j> ˜ Ji j~ Si·~ Sj−2Ji j ∑ i j ˜ Ji j~ Si·~ Sj−Ji j Table 1: The conversion of Ji j to other conventions, where ˜ Jis the exchange parameter in that convention. The notation <i j >means a pair of i j without counting it twice. The DMI parameters ~ Dcan be converted in the same way. where the first term represents the single-ion anisotropy (SIA), the second is the isotropic exchange, and the third term is the symmetric anisotropic exchange, where Jani is a 3 ×3 symmetric tensor. The final term is the DMI, which is antisymmetric. Importantly, the SIA is not accessible from Wannier90 as it requires separately the spin-orbit coupling part of the Hamiltonian [53]. However, it is readily accessible from constrained DFT calculations [54]. We note that there are several conventions for the Heisenberg Hamiltonian, here we take a commonly used one in atomic spin dynamics: we use a minus sign in the exchange terms, i.e. positive exchange Jvalues favor ferromagnetic alignment. Every pair i j is taken into account twice, Ji j and Jji are both in the Hamiltonian. The spin vectors ~ Siare normalized to 1, so that the parameters are in units of energy. The other commonly used conventions differ in a prefactor 1/2 or a summation over different i j pairs only. The conversion factors to other conventions are given in Table 1. For other conventions in which the spins are not normalized, the parameters need to be divided by |Si|Sj in addition. From the total energy due to the spin interactions, Eq. (1), we obtain the following variation with respect to the ~ Siand ~ Sj δEi j =−2Jiso i j δ~ Si·δ~ Sj −2δ~ SiJani i j δ~ Sj −2~ Di j ·(δ~ Si×δ~ Sj) (2) 2.2. Tight-binding Hamiltonian and Green’s function We start from a generalized tight-binding Hamiltonian. The localized basis functions are denoted as ψimσ(~r)with i,m, and σbeing the site, orbital, and spin indices, respectively. Due to translation symmetry, the tight-binding Hamiltonian, H, and the overlap, S, matrices can be parameterized as Him jm0σσ0(~ R) = hψimσ(~r)|H|ψjm0σ0(~r+~ R)i,(3) Sim jm0σσ0(~ R) = hψimσ(~r)| |ψjm0σ0(~r+~ R)i,(4) where~ris the position inside the unit cell, ~ Ris the lattice vector, and Hdenotes the total Hamiltonian. The overlap matrix Sreduces to the identity matrix when the basis functions are orthonormal. TB2J can use both non-orthogonal LCAO basis sets and orthogonal Wannier basis sets. Below we discuss only the results for an orthogonal basis set. The non-orthogonal basis set is discussed in Ref. [50], where it is shown that the expressions for the exchange parameters are the same as for an orthogonal basis. In the following, we drop all orbital and spin indices for simplicity, Hence Hi j is a sub-matrix of Hcontaining all the spin and orbital components for atoms iand j. We note that in DFT it is possible to do non-collinear calculations with or without SOC (the exchange correlation potential can induce a non-collinear ground state). If one were to perform a calculation with SOC and the spin constrained to one direction, it would be categorized as collinear. In the collinear case, His diagonal in the spin subspace. The Green’s function in reciprocal space is defined as G( ~ k,ε) = εS( ~ k)−H( ~ k)−1,(5) where H( ~ k) = ∑~ RH(~ R)ei ~ k·~ R, and S( ~ k) = ∑~ RS(~ R)ei ~ k·~ R. The Green’s function in real space is obtained using the following expression G(~ R,ε) = ZBZ G( ~ k,ε)e−i ~ k·~ Rd ~ k.(6) In the following, we drop the ~ Rin the atom pair labeled by i,j,~ R, e.g. Gi j means Gi j(~ R), and Gji means Gji(−~ R). From Eq. (7) all equations are given in real space. For each atom i, the intra-atomic component of His defined as Pi=Hii(~ R=0)and is of size 2Norb ×2Norb. Each Pi,mm0is a 2 ×2 matrix in spin, which can be decomposed into its scalar and vector parts Pimm0=p0 imm0I+~pimm0·~ σ,(7) =p0 imm0I+pimm0~eimm0·~ σ, where ~pimm0is the vector part of P, (upper case Pdenotes matrices in spin and orbitals, lower case pis used 3 for matrices in orbitals only) which has the x,y, and z components px imm0,py imm0,pz imm0.~eimm0is the unit orientation vector of ~pimm0, and ~ σ= (σx,σz,σz)are the three Pauli matrices. In condensed form, the spin and orbital matrix for site iis now Pi=p0 iI+~ pi·~ σ. We can decompose the Green’s function for each inter-site orbital pair Gim,jm0in the same way Gim,jm0=G0 im,jm0I+~ Gim,jm0·~ σ,(8) where G0 im,jm0and ~ Gim,jm0form the G0 i j and ~ Gi j matrices, respectively. 2.3. Magnetic force theorem According to the force theorem the total energy variation due to a small perturbation from the ground state coincides with the change of the single-particle energies at fixed ground-state potential δE=ZEF −∞ εδn(ε)dε=−ZEF −∞ δN(ε)dε(9) where n(ε) = −1 πImTr(G(ε)) is the density of states and N(ε) = −1 πImTr(ε−H)is the integrated density of states. The traces are taken over orbitals only, not spin. Thus, the first order variation of Ndue to δHcan be written as δN(ε) = 1 πImTr(δHG); the second order variation is δ2N(ε) = 1 πImTr(δHGδHG). Now we use the spin rotation as a perturbation. For the rotation of the spin at site i, the change in the energy up to the second order reads as δE1spin i=−1 πZEF −∞ ImTr(δHiG+δHiGδHiG)dε. (10) Similarly, for the rotation of the spins at sites iand j, the change in energy is given by δE2spin i j =−1 πZEF −∞ ImTr[δHiG+δHiGδHiG +δHjG+δHjGδHjG +2δHiGδHjG]dε. (11) The energy variation due to the two-spin interaction is then δEi j =δE2spin i j −δE1spin i−δE1spin j= −2 πZEF −∞ ImTr(δHiGδHjG)dε. (12) The change of Hdue to the rotation of spin is ~ δφ ×~ p with the rotation axis along ~ δφ and the angle  ~ δφ. By putting this into Eqn. 12, we get δEi j =−2[A00 i j −∑ u=x,y,z Auv i j ]δ~ei·δ~ej −2∑ u,v∈x,y,z δeu i[Auv i j +Avu i j ]δev j −2~ di j ·(δ~ei×δ~ej) (13) in which the 4 ×4 matrix Ai j is defined as Auv i j =−1 πZEF −∞ Trnpz iGu i jpz jGv jiodε,(14) where u,v∈ {0,x,y,z}, and the component of ~ di j,du i j = Re(A0u i j −Au0 i j ). Comparing Eq. (13) to Eq. (2), we can find that the values of the isotropic exchange Jiso , the anisotropic exchange Jani, and the DMI ~ Dcan be expressed as Jiso i j =Im(A00 i j −Axx i j −Ayy i j −Azz i j ),(15) Jani,uv i j =Im(Auv i j +Avu i j ),(16) Du i j =Re(A0u i j −Au0 i j ),(17) which is the same with Ref. [30]. If spin-orbit coupling (SOC) is neglected, A0u i j =Au0 i j , i.e. the DMI term is zero and both pand Gonly have components along the spin quantization axis, say the zdirection. Then, the xand ycomponents of ~ pand ~ Gvanish. Thus, ~ pi= (0,0,pz i),Gx=Gy=0, G0= 1 2(G↑+G↓),Gz=1 2(G↑−G↓). The isotropic exchange parameter then reduces to Jiso i j =Im(A00 i j −Azz i j ).(18) Defining ∆i=p↑ i−p↓ i=2|~ pi|=2pz iwe obtain from Eq. (18) the LKAG [19] expression for isotropic exchange Jiso i j =−1 4πZEF −∞ ImTrn∆iG↑ i j∆jG↓ jiodε.(19) 2.4. xyz average strategy We note that when all spins are oriented along one direction (e.g. z), the components with u=zor v=zin the Jtensor are non-zero if the variation with respect to the rotation is only kept to first order. The xz,yz,zx,zy, and zz components of the anisotropic exchange, and the Dzin the DMI cannot be obtained from a single calculation with magnetization along the zdirection. It should be noted that this does not affect the properties too much if the system stays close to the reference spin state, as these terms are small and higher order in the spin rotation angles. For example, in a AFM with small spin 4 canting system, the parameters obtained from the AFM reference state can be used directly to get the canting angles. To determine the missing parameters, we rotate the whole system and lattice, (equivalently one can rotate the quantization axis) from zto the xand ydirections and apply the back-rotation to the DMI vectors. A weighted average over all three directions is taken: the inaccessible components are given a weight of 0 and the rest have the same weight. This small trick allows not only to obtain the full set of Jani ~ Dvector components , but also to reduce the numerical noise. The DMI is usually much smaller than the isotropic exchange, such that its calculation is more numerically delicate. This issue arises in particular when Wannier functions are used whose symmetry is not guaranteed, and where the disentanglement procedure can introduce numerical errors. A direct approach to calculate the Dzterm is proposed in Ref. [25], taking the spin-rotation perturbation to higher order. This is not implemented in TB2J, and would also be more sensitive to numerical noise. 2.5. Higher order terms It has been shown by several authors that the parameters from the Green’s function method cannot be mapped exactly onto a bilinear Jtensor [29,55,30] and higher order terms need to be considered, which include multi-spin interactions and the higher-order twospin terms. In the simplest case, following Ref. [30], the Hamiltonian can be written in a biquadratic form as HQ=−∑ i,j J0 i j~ Si~ Sj−∑ i,j Bi j ~ Si~ Sj2,(20) where J0 i j and Bi j are determined by the Amatrices as J0 i j =A00 i j −3Azz i j ,(21) Bi j =Azz i j .(22) When ~ Siand ~ Sjare close to their reference values, we have Ji j =−d2HQ d~ Sid~ Sj =J0 i j +2Bi j~ Si~ Sj(23) ≃J0 i j +2Bi j~ Sref i~ Sref j=J0 i j +2Bi j =A00 i j −Azz i j , which is equivalent to the Ji j expression with only bilinear terms considered, and it shows that the effective bilinear Jterm depends on the orientation of Siand Sj. This method was proposed in order to improve the description of the system when the deviation from the reference state is large [30]. From our limited experience, however, the parameters produced by Eqs. (21) and (22) are not always physical and their fit is more complex. 3. Implementation The TB2J package is implemented with three main modules, the tight-binding (TB) module, the exchange module, and the output module (see Fig. 1for a schematic view). The TB module provides the interface with external DFT and Wannier function codes and calculates the Green’s function. The exchange module uses the Green’s function to calculate the magnetic interaction parameters. The output module writes these parameters to output files and provides a python API from which the magnon band structure can be calculated. The TB module defines the classes of the TB model and the Green’s function. The internal TB module can read the Wannier90 output data files to build the Hamiltonian. The WF’s are assigned to the nearest atom with the crystal periodicity taken into account. A warning is issued if a WF is far away from any atom. The code utilizes the “duck type” feature: a class can be plugged in if the specific required methods are defined. The methods include the calculation of H( ~ k)and the eigenvalues and eigenvectors. Hence, external libraries can be wrapped easily. An interface to the sisl library [52] has been implemented, which includes the tight-binding model from Siesta and OpenMX. We also wrap directly the python based DFT code GPAW in LCAO mode instead of reading the Hamiltonian from the output files. Starting from the eigenvalues and eigenvectors, the Green’s functions can be calculated without inverting the εS( ~ k)−H( ~ k)matrix for every ε. Instead, the Green’s functions are calculated with G( ~ k,ε) = Ψ( ~ k)hεI−Diag(E( ~ k))i−1Ψ( ~ k)†, where E( ~ k)and Ψ( ~ k) are the eigenvalues and the eigenvector matrices, respectively. The real-space Green’s function can then be calculated by Fourier transforming G(k). Once the Green’s functions are calculated, the exchange module calculates the magnetic interaction parameters using Eqs. (15)-(17). The Hamiltonian and the Green’s functions are first decomposed into their {0,x,y,z}components. Then, the elements of the 4 ×4 matrix, TrnpiGu i jpjGv jiofor each εare calculated and later integrated to obtain the Ai j matrix. A contour integration method is used for the REFdε integration. The range of the integration is (Emin,EF), where Emin is either below the lowest band energy or chosen such that the orbitals below Emin have only negligible interaction with those near EF. By default, a semicircle path is used, which is centered at (Emin +EF)/2 on the real axis and has a radius of (EF−Emin)/2, going through the upper half of the complex plane with 5 Wannier90: Sisl: DFT Abint Quantum Espresso VASP Wien2K Fleur ... TB2J TBGreen H(k) H(R) Eigen(k) G(k) G(R) TB Exchange Magnon Band txt netcdf xml User's format LCAO Siesta OpenMX Effective TB Hamiltonian Spin Dynamics -Multibinit -Tom's GPU ASD -Vampire Collinear Isotropic J Bi-quadratic Non-collinear Isotropic J Bi-quadratic DMI Anisotropic J IO Figure 1: Schematic diagram view of TB2J workflow. Im(ε)>0. The magnetic interaction parameters are calculated from the Aij matrix. The full list of quantities that are available from TB2J is given in Table 2. Quantity Expression collinear case Ji j −1 4πREF −∞dεImTrn∆iG↑ i j∆jG↓ jio J0 i j Im(A00 i j −3Azz i j ) Bi j Azz i j non-collinear case Ji j Im(A00 i j −Axx i j −Ayy i j −Azz i j ) ~ Du i j Re(A0u i j −Au0 i j ) Jani Jani,uv =Jani,vu =Im(Auv +Avu) Table 2: Summary of magnetic interaction quantities that are calculated within TB2J. All come from Ref. [30] except the first (Ref. [19]) and the DMI, which was first written in Ref. [24] 4. TB2J user instructions In this section, we describe the practical usage of TB2J. The users can also refer to the online documentation at https://tb2j.readthedocs.io/en/ latest/, which is more detailed and regularly updated. We show how to install and use the code to calculate the magnetic interaction parameters using body centered cubic (BCC) Fe, and hexagonal close packed (HCP) Co, BiFeO3, and La2CuO4as examples. 4.1. Installation TB2J can be installed with a simple command if python and pip environments are pre-installed pip install TB2J Alternatively, one can download the package from https://github.com/mailhexu/TB2J and run python setup . py install to install the package. By default, TB2J only installs the hard (non optional) dependencies automatically. Hence, the sisl package, which is used to read the Hamiltonian from the Siesta or OpenMX output, needs to be installed separately using pip. Also, the GPAW-TB2J interface is through python directly, which requires the GPAW python package. 4.2. Computing the magnetic interaction parameters In this section we describe the TB2J procedures to compute the magnetic interactions, from an electronic structure calculation to the final output file. 6 4.2.1. Preparation of electronic structure and tightbinding Hamiltonian To obtain the magnetic interaction parameters, the first step is to do a converged DFT calculation of a magnetic crystal, at the collinear or non-collinear level. Preferably, this calculation treats the magnetic ground state of the system, however, the magnetic interaction parameters can also be calculated for a different reference state. It should be noted that the spin quantization axis is presumed to be along the zaxis throughout TB2J. The next step is to construct the tight-binding Hamiltonian. For DFT codes using non-local basis sets, like plane waves, the WFs can be constructed for any code which has an interface with Wannier90. All the spinpolarized orbitals of all atoms contributing to the magnetic interaction should be carefully selected to compute accurately the magnetic interactions interaction parameters. For example, in the case of a transition metal oxide the dorbitals of the transition metal cation should be included, but also the oxygen 2porbitals that are involved in the superexchange or DMI mechanism through hybridization with dorbitals. This forms the minimal basis of orbitals to be included in the construction of the WFs but other orbitals might be important as well. The number of orbitals to be included is system dependent and should be checked by the user (convergence of the calculated magnetic interactions with respect to the number of Wannier orbitals). For example, in the case of SrMnO3the Mn-3d, O-2porbitals are necessary. The options to build the Wannier functions should be enabled in the DFT codes. For example, in Abinit, the “prtwant 2” and “w90iniprj 2” options (ABINIT version 9.x) are needed for the calculation of maximally localized Wannier functions (MLWFs). The input file for Wannier90 also needs to be present in the execution directory. THe quality of the Wannier function Hamiltonian is need to be checked. The magnetic interaction parameters are often meV or µeV, which requires the noise in the Hamiltonian to be lower. The rigid spin rotation on one site is performed by rotating all the spins of the WFs associated to a given atom, which requires the WFs to be centered on an atom or at least very close to it. TB2J uses the Wannier centers to decide which atom each WF “belongs” to. As a result, using the MLWFs [56] might not always be the best choice. Projected WFs or selectively localized WFs [57], which add a constraint on the Wannier centers, can be used instead. The Wannier centers for one atom might be located closer to a periodic copy of that atom than the original atom. To solve this problem, TB2J shifts the Wannier centers by using the transnational symmetry and modifies the Hamiltonian accordingly. The WF Hamiltonian and the center positions need to be output by Wannier90 using the “write hr” and the “write xyz” input flags. For DFT codes based on LCAO basis sets, such as Siesta, the Hamiltonian is already localized and needs no further transformation. For calculating the parameters of the Heisenberg Hamiltonian only the localized DFT Hamiltonian and the overlap matrix need to be saved. For example, one can use the options “CDF.Save=True”, “SaveHS=True”, and “Write.DMHS.Netcdf=True” in Siesta (version 4.x) to enable the saving of these matrices. 4.2.2. Running TB2J TB2J has two Python executables: wann2J.py and siesta2J.py, for calculating Jfrom Wannier90 and Siesta output, respectively. A similar script named openmx2J.py is in the TB2J OpenMX package. These scripts are designed to have a minimal user input where only the paths to the files containing the electron Hamiltonian information and defining the magnetic atom species need to be provided. With Wannier90. The executable script wann2J.py can be used with the Wannier90 output files. For a noncollinear calculation, the spin up and spin down channel Wannier functions need to be present. In addition, the following parameters need to be specified: •Whether the calculation is collinear or noncollinear, given by the –spinor option. TB2J assumes a collinear calculation by default, and the usage of –spinor specify that the calculation is non-collinear, where the Hamiltonian is in a spinor form. •the prefix to the paths of up and down Wannier functions, (–prefix up abinito w90 up – prefix down abinito w90 down, for non-collinear calculations, and –prefix spinor for non-collinear calculations). The filename of the Hamiltonian is the prefix plus “ hr.dat”. •The posfile specifies the name of a file containing the atomic structure and cell parameters. ASE file formats are also readable (the full list of which can be found on https://wiki.fysik.dtu.dk/ ase/ase/io/io.html). It should be noted that some ASE formats cannot be used, e.g. xyz, because they do not contain the cell parameters which are required by TB2J. We also recommend using formats with the cell matrix rather than only the 7 (a,b,c,α,β,γ), which will cause trouble if there are anisotropic or DMI terms as they are not rotationally invariant. •The Fermi energy in units of eV. •The type of magnetic elements (symbol from the periodic table). Here is an example for calculating the Js in the noncollinear case: wann2J . py -- posfile abinit . in -- efermi 5.8 -- elements Fe -- prefix_up abinito_w9 0_up -- prefix_down abinito_w90_down In the case of a non-collinear, the WF Hamiltonian is in a single file and we need to specify the calculation type: wann2J . py -- posfile abinit . in -- efermi 5.8 -- elements Fe -- prefix_spinor abinito_w90 -- spinor With Siesta. Only a minimal set of parameters is needed for Siesta: the filename of the input for the Siesta calculation. The other information needed, including whether the calculation has SOC enabled and the Fermi energy, are found in the Siesta results: siesta2J . py -- element Fe -- input - fname = ’ siesta .fdf ’ With OpenMX. The interface to OpenMX is distributed as a plugin to TB2J called TB2J OpenMX under the GPL license, which need to be installed separately, because code from OpenMX which is under the GPL license is used in the parser of OpenMX files. pip install TB2J_OpenMX In the DFT calculation, the ”HS.fileout on” options should be enabled, so that the Hamiltonian and the overlap matrices are written to a ”.scfout“ file. The necessary input are the path of the calculation, the prefix of the OpenMX files, and the magnetic elements: openmx2J .py --path ./ -- prefix openmx -- elements Fe General options. There are several tunable parameters, which the user usually does not need to specify. Among them, there are: •nz: The number of steps in the path of the contour integration. •emin, emax: the integration lower and upper bounds, relative to the Fermi energy. The value of emin is automatically determined if not given, emax should be about 0 for metallic systems, whereas it lies in the band gap for insulating systems. •rcut: the cutoff or max distance between spin pairs. The full list of options can be accessed using “–help”. As discussed in the previous section, the zcomponent of the DMI, and the xz,yz,zx,zy,zz components of the anisotropic exchanges are non-physical, and an xyz average is needed to get the full set of magnetic interaction parameters. In this case, scripts to rotate the structure and merge the results are provided, they are named TB2J rotate.py and TB2J merge.py. The TB2J rotate.py reads the structure file and generates three files containing the z→x,z→yand the nonrotated structures. The output files are named atoms x, atoms y, atoms z. A large number of output file formats is supported thanks to the ASE library [58] and the format of the output structure files is provided using the “–format” parameter. An example for using the rotate file is: TB2J_rotate . py BiFeO3 . cif -- format cif The user has to perform DFT single point energy calculations for these three structures in different directories, keeping the spins along the zdirection, and run TB2J on each of them. After producing the TB2J results for the three rotated structures, we can merge the DMI results with the following command by providing the paths to the TB2J results of the three cases, e.g.: TB2J_merge . py BiFeO3_x BiFeO3_y BiFeO3_z --type structure A new TB2J results directory is then made which contains the merged final results. 4.2.3. Output files In the following we describe the output files which TB2J produces. By running wann2J.py or siesta2J.py, a directory with the name TB2J results will be generated, which contains the following output files: •exchange.out: A human readable output file, which summarizes the results. •Multibinit: A directory containing output which can be read directly by the Multibinit code [38]. 8 ========================================================================================== Information : Exchange parameters generated by TB2J 0.2.8. ========================================================================================== Cell ( Angstrom ) : 0.030 3.950 3.950 3.950 0.030 3.950 3.950 3.950 0.030 ========================================================================================== Atoms : (Note : charge and magmoms only count the wannier functions .) Atom_number x y z w_charge M(x) M(y) M(z) Bi1 0.2413 0.2413 0.2413 2.1878 -0.0010 -0.0005 -0.0045 Bi2 4.2060 4.2060 4.2060 2.1878 0.0005 0.0010 0.0045 Fe1 2.0165 2.0165 2.0165 6.1722 -0.0027 -0.0021 4.1151 Fe2 5.9812 5.9812 5.9812 6.1722 0.0021 0.0027 -4.1151 O1 5.5238 2.1558 3.9388 4.8807 -0.0005 0.0011 -0.0568 O2 6.0903 5.5389 3.9539 4.8807 -0.0011 0.0005 0.0568 O3 3.9388 5.5238 2.1558 4.8806 -0.0016 -0.0027 -0.0559 O4 3.9539 6.0903 5.5389 4.8806 0.0001 -0.0031 0.0562 O5 2.1558 3.9388 5.5238 4.8806 0.0031 -0.0001 -0.0562 O6 5.5389 3.9539 6.0903 4.8806 0.0027 0.0016 0.0559 Total 46.0038 0.0015 -0.0015 -0.0000 ========================================================================================== Exchange : i j R J_iso ( meV ) vector distance ( A) ---------------------------------------------------------------------------------------- Fe2 Fe1 ( 0, 1, 1) -26.7976 ( 3.934 , 0.015 , 0.015) 3.934 J_iso : -26.7976 [ Experimental !] Jprime : -34.444 , B: -3.810 [ Experimental !] DMI : ( 0.1590 -0.0996 0.0358) [ Experimental !] J_ani : [[ -0.026 0.002 -0.01 ] [ 0.002 -0.027 -0.05 ] [ -0.01 -0.05 -7.62 ]] Listing 1: An example of the output sections for BiFeO3with SOC enabled, calculated with spin along the zaxis. The exchange.out file contains three sections: cell, atoms and exchange. The cell section contains the lattice parameter matrix. The atoms section contains the positions, charges (for verification) and magnetic moments of the atoms: see Listing 1. Here, the charge and magnetic moment of each atom are only integrated with the WFs attached to this atom. As such they can differ from the quantities coming from the direct DFT output, as not all bands are used in the construction of WFs. The WF charges should be integers, for LCAO the values depend on the band energy cutoffs. In addition, the exclusion of very deep lying levels from the calculation of Jcan also lead to deviations in the charges which might appear both for WFs and for LCAO. Another source of difference between the TB2J charges and magnetic moments and the DFT ones is the integration volume around the atoms, which is not necessarily the same. However, for localized d and forbitals the magnetic moments should be close to their DFT counterparts, for TB2J to yield correct results for the parameters. Large differences between the TB2J and DFT values indicate that something may have gone wrong: either in in the contour integration REFdεused in TB2J, or in the construction of the Wannier functions (incorrect wannierization process or too small WF basis set). Often, it comes from excluding an orbital that is important for the magnetic interaction in the studied system. In the case of metallic system, the Fermi energy might have to be slightly shifted with respect to the DFT reference due to different numerical method used in the integration of charge density. Each pair of atoms is labeled by three parameters, the index i,jand R, where iand jare the indices in the unit cell. The vector ~ Rspecifies the cell the atom jis trans9 Acknowledgements The authors thank Yajun Zhang, Alireza Sasani, Jorge Pilo Gonzlez, and Zachary Romestan for the testing of the code, Thomas Ostler, Bertrand Dup´ e and Phivos Mavropoulos for explanations about the limits and intricacies of fitting the Heisenberg model. This work has been funded by the Communaut´ e Franc¸aise de Belgique (ARC AIMED G.A. 15/19-09).XH thanks the support by the EU H2020NMBP-TO-IND-2018 project ”INTERSECT” (Grant No. 814487). EB thanks the FRS-FNRS for support, as does MJV for an “out” sabbatical grant to ICN2 Barcelona in 2018-2019. The authors acknowledge the CECI supercomputer facilities funded by the F.R.S-FNRS (Grant No. 2.5020.1) and the Tier-1 supercomputer of the F´ ed´ eration Wallonie-Bruxelles funded by the Walloon Region (Grant No. 1117545). Computing time was also provided by PRACE-3IP DECI grants 2DSpin and Pylight on Beskow (G.A. 653838 of H2020). References [1] Hirohata, A. et al., Journal of Magnetism and Magnetic Materials 509 (2020) 166711. [2] Baltz, V. et al., Rev. Mod. Phys. 90 (2018) 015005. [3] Bhatti, S. et al., Materials Today 20 (2017) 530 . [4] Joshi, V. K., Engineering Science and Technology, an International Journal 19 (2016) 1503 . [5] Lu, J. W., Chen, E., Kabir, M., Stan, M. R., and Wolf, S. A., International Materials Reviews 61 (2016) 456. [6] ˇ Zuti´ c, I., Fabian, J., and Das Sarma, S., Rev. Mod. Phys. 76 (2004) 323. [7] Fazekas, P., Lecture notes on electron correlation and magnetism, Series in modern condensed matter physics 5, World Scientific, 1999. [8] Chumak, A. V., Vasyuchka, V. I., Serga, A. A., and Hillebrands, B., Nat. Phys. 11 (2015) 453. [9] Kiselev, N. S., Bogdanov, A. N., Schfer, R., and Rßler, U. K., Journal of Physics D: Applied Physics 44 (2011) 392001. [10] Fert, A., Cros, V., and Sampaio, J., Nat. Nanotechnol. 8(2013) 152. [11] Xiang, H., Lee, C., Koo, H.-J., Gong, X., and Whangbo, M.-H., Dalton Trans. 42 (2013) 823. [12] Herring, C., Magnetism: a treatise on modern theory and materials. 4. Exchange interactions among itinerant electrons, Academic Press, 1966. [13] Sandratskii, L. M., physica status solidi (b) 136 (1986) 167. [14] Savrasov, S. Y., Phys. Rev. Lett. 81 (1998) 2570. [15] Romero, A. H. et al., The Journal of Chemical Physics 152 (2020) 124102. [16] Cao, K., Lambert, H., Radaelli, P. G., and Giustino, F., Phys. Rev. B 97 (2018) 024420. [17] Ricci, F., Prokhorenko, S., Torrent, M., Verstraete, M. J., and Bousquet, E., Phys. Rev. B 99 (2019) 184404. [18] Phillips, J. J. and Peralta, J. E., The Journal of Chemical Physics 138 (2013) 174115. [19] Liechtenstein, A. I., Katsnelson, M., Antropov, V., and Gubanov, V., Journal of Magnetism and Magnetic Materials 67 (1987) 65. [20] Steenbock, T. and Herrmann, C., Journal of Computational Chemistry 39 (2018) 81. [21] Katsnelson, M. I. and Lichtenstein, A. I., Phys. Rev. B 61 (2000) 8906. [22] Katsnelson, M. I., Kvashnin, Y. O., Mazurenko, V. V., and Lichtenstein, A. I., Phys. Rev. B 82 (2010) 100403. [23] Solovyev, I., Hamada, N., and Terakura, K., Phys. Rev. Lett. 76 (1996) 4825. [24] Antropov, V., Katsnelson, M., and Liechtenstein, A., Physica B: Condensed Matter 237 (1997) 336. [25] Mazurenko, V. V. and Anisimov, V. I., Phys. Rev. B 71 (2005) 184434. [26] Mankovsky, S. and Ebert, H., Phys. Rev. B 96 (2017) 104416. [27] Secchi, A., Lichtenstein, A. I., and Katsnelson, M. I., Annals of Physics 360 (2015) 61. [28] Secchi, A., Lichtenstein, A. I., and Katsnelson, M. I., Journal of Magnetism and Magnetic Materials 400 (2016) 112. [29] Lounis, S. and Dederichs, P. H., Phys. Rev. B 82 (2010) 180404. [30] Szilva, A. et al., Phys. Rev. Lett. 111 (2013) 127204. [31] Yavorsky, B. Y. and Mertig, I., Phys. Rev. B 74 (2006) 174402. [32] Ebert, H. and Mankovsky, S., Phys. Rev. B 79 (2009) 045209. [33] Andersen, O. K. and Jepsen, O., Phys. Rev. Lett. 53 (1984) 2571. [34] Korotin, D. M., Mazurenko, V. V., Anisimov, V. I., and Streltsov, S. V., Phys. Rev. B 91 (2015) 224405. [35] Mostofi, A. A. et al., Comput. Phys. Commun. 178 (2008) 685 . [36] Pizzi, G. et al., Journal of Physics: Condensed Matter 32 (2020) 165902. [37] Gonze, X. et al., Comput. Phys. Commun. 180 (2009) 2582 , 40 YEARS OF CPC: A celebratory issue focused on quality software for high performance, grid and novel computing architectures. [38] Gonze, X. et al., Comput. Phys. Commun. 248 (2020) 107042. [39] Giannozzi, P. et al., Journal of Physics: Condensed Matter 29 (2017) 465901. [40] Soler, J. M. et al., Journal of Physics: Condensed Matter 14 (2002) 2745. [41] Kresse, G. and Furthm¨ uller, J., Phys. Rev. B 54 (1996) 11169. [42] Kresse, G. and Joubert, D., Phys. Rev. B 59 (1999) 1758. [43] Blaha, P., Schwarz, K., Madsen, G. K., Kvasnicka, D., and Luitz, J., An augmented plane wave+ local orbitals program for calculating crystal properties (2001). [44] Bl¨ ugel, S. and Bihlmayer, G., Forschungszentrum J¨ ulich GmbH (2006) 85. [45] Andrade, X. et al., Physical Chemistry Chemical Physics 17 (2015) 31371. [46] Tancogne-Dejean, N. et al., The Journal of Chemical Physics 152 (2020) 124119. [47] Ozaki, T. and Kino, H., Physical Review B 72 (2005) 045121. [48] Enkovaara, J. et al., Journal of physics: Condensed matter 22 (2010) 253202. [49] The Elk Code, http://elk.sourceforge.net/. [50] Oroszl´ any, L., Ferrer, J., De´ ak, A., Udvardi, L., and Szunyogh, L., Phys. Rev. B 99 (2019) 224412. [51] Yoon, H., Kim, T. J., Sim, J.-H., and Han, M. J., Comput. Phys. Commun. 247 (2020) 106927. [52] Papior, N. R., B., J. L., Frederiksen, T., and Wuhl, S. S., zerothi/sisl: v0.9.8, 2020. [53] Solovyev, I., Dederichs, P., and Mertig, I., Physical Review B 52 (1995) 13419. [54] Weingart, C., Spaldin, N., and Bousquet, E., Phys. Rev. B 86 (2012) 094413. 16 [55] Udvardi, L., Szunyogh, L., Palot´ as, K., and Weinberger, P., Phys. Rev. B 68 (2003) 104436. [56] Marzari, N. and Vanderbilt, D., Phys. Rev. B 56 (1997) 12847. [57] Wang, R., Lazar, E. A., Park, H., Millis, A. J., and Marianetti, C. A., Phys. Rev. B 90 (2014) 165125. [58] Larsen, A. H. et al., Journal of Physics: Condensed Matter 29 (2017) 273002. [59] Evans, R. F. et al., Journal of Physics: Condensed Matter 26 (2014) 103202. [60] Di Gennaro, M., Miranda, A. L., Ostler, T. A., Romero, A. H., and Verstraete, M. J., Phys. Rev. B 97 (2018) 214417. [61] Catalan, G. and Scott, J. F., Adv. Matter. 21 (2009) 2463. [62] Pajda, M., Kudrnovsk´ y, J., Turek, I., Drchal, V., and Bruno, P., Phys. Rev. B 64 (2001) 174402. [63] Perdew, J. P., Burke, K., and Ernzerhof, M., Phys. Rev. Lett. 77 (1996) 3865. [64] Van Setten, M. et al., Comput. Phys. Commun. 226 (2018) 39. [65] Garc´ ıa, A., Verstraete, M. J., Pouillon, Y., and Junquera, J., Comput. Phys. Commun. 227 (2018) 51. [66] Torrent, M., Jollet, F., Bottin, F., Zrah, G., and Gonze, X., Computational Materials Science 42 (2008) 337 . [67] Jollet, F., Torrent, M., and Holzwarth, N., Comput. Phys. Commun. 185 (2014) 1246 . [68] Perdew, J. P. et al., Phys. Rev. Lett. 100 (2008) 136406. [69] Anisimov, V. I., Zaanen, J., and Andersen, O. K., Phys. Rev. B 44 (1991) 943. [70] Lee, J. H. and Rabe, K. M., Phys. Rev. B 84 (2011) 104440. [71] Glazer, A., Acta Crystallographica Section B: Structural Crystallography and Crystal Chemistry 28 (1972) 3384. [72] Xu, C., Xu, B., Dup´ e, B., and Bellaiche, L., Phys. Rev. B 99 (2019) 104420. [73] Matsuda, M. et al., Phys. Rev. Lett. 109 (2012) 067205. [74] Ozaki, T. and Kino, H., Phys. Rev. B 69 (2004) 195113. [75] Liechtenstein, A., Anisimov, V., and Zaanen, J., Phys. Rev. B 52 (1995) R5467. [76] Chmaissem, O. et al., Phys. Rev. B 64 (2001) 134412. [77] Pizzi, G., Cepellotti, A., Sabatini, R., Marzari, N., and Kozinsky, B., Computational Materials Science 111 (2016) 218. [78] Damle, A., Lin, L., and Ying, L., Journal of chemical theory and computation 11 (2015) 1463. [79] Damle, A. and Lin, L., Multiscale Modeling & Simulation 16 (2018) 1392. [80] Vitale, V. et al., npj Computational Materials 6(2020) 1. [81] Varignon, J., Bibes, M., and Zunger, A., Nature Commun. 10 (2019) 1658. [82] Himmetoglu, B., Floris, A., de Gironcoli, S., and Cococcioni, M., International Journal of Quantum Chemistry 114 (2014) 14. [83] Bousquet, E. and Spaldin, N., Phys. Rev. B 82 (2010) 220402. [84] Sakuma, R., Phys. Rev. B 87 (2013) 235109. [85] Fedorova, N. S., Ederer, C., Spaldin, N. A., and Scaramucci, A., Phys. Rev. B 91 (2015) 165122. [86] Ostler, T. A., Barton, C., Thomson, T., and Hrkac, G., Phys. Rev. B 95 (2017) 064415. [87] Mankovsky, S., Polesya, S., and Ebert, H., Phys. Rev. B 101 (2020) 174401. [88] Hoffmann, M. and Bl¨ ugel, S., Phys. Rev. B 101 (2020) 024418. 17