TB2J: a python package for computing magnetic interaction parameters Xu Hea,b,c,∗, Nicole Helbigc , Matthieu J. Verstraeteb,c , Eric Bousqueta a Physique Théorique des Matériaux, Q-MAT, CESAM, Université de Liége, B-4000 Sart-Tilman, Belgium b CatalanInstitute of Nanoscience and Nanotechnology (ICN2), CSIC, BIST, Campus UAB, Bellaterra, Barcelona, 08193, Spain c Nanomat, Q-Mat, CESAM, and European Theoretical Spectroscopy Facility, Université de Liège, B-4000 Liège, 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 high- throughput workflows. 1. Introduction culations or energy mapping analysis [11]. This ap- proach necessitates the calculation of the total energies First-principles simulations of magnetic materials of different magnetic configurations. The parameters of have attracted strong interest in the last decade due the Heisenberg Hamiltonian are then fit to these ener- to an increase in precision provided by density func- gies under the supposition that the change in energy is tional theory (DFT) codes for strongly correlated mag- only related to the magnetic interactions. This requires netic atoms and also because of the developments in to have at least as many calculated magnetic configura- spintronics applications [1, 2, 3, 4, 5, 6]. Understand- tions as the number of parameters of the Hamiltonian ing the complex microscopic origin of the magnetic in- (though more magnetic configurations are usually nec- teractions often requires a reduction of the full many- essary to converge the fit), which often requires the use body electronic interactions to effective Hamiltonians, of large supercells. This simple method gives good re- the two most important ones being the Heisenberg and sults in some cases but the Heisenberg model itself can the Hubbard Hamiltonians [7]. The parameters of these break down if the chosen magnetic configurations devi- Hamiltonians are fit to DFT data to give access to an ate too much from the ground-state one. This can hap- easier understanding of the magnetic interactions and pen, e.g., if a magnetic phase, often the ferromagnetic to allow for simulations of larger systems and/or with (FM) one, closes the band gap of an insulating antifer- dynamics. Novel sensing, storing and computing tech- romagnetic (AFM) ground-state. Also, the method be- nologies have been proposed using spin waves [8] and comes unsuitable if the supercell needed to probe all of skyrmions [9, 10]. Theoretical models have consis- the pertinent magnetic configurations is too large to be tently opened new vistas and explained experimental handled by DFT. For example, if certain configurations findings in complex new chemistries, geometries, and show a delocalized picture for the electrons they must heterostructures. Reliable ab initio calculations of the be excluded from the fit, which typically leads to an in- effective Hamiltonian parameters is crucial for predic- crease in the size of the supercell needed to provide a tive and quantitative simulations. Fitting these param- sufficiently large number of configurations. eters is, however, often very cumbersome and necessi- Another method to determine the parameters is the tates a case-by-case construction [11]. Generalized Bloch Theorem (GBT) [12, 13] which For the Heisenberg Hamiltonian, one of the most adapts the boundary conditions on the spin orientation common fitting procedures is through total energy cal- to account for a spin spiral. From the change in total energy one can extract the magnetic interaction parame- ∗ Correspondingauthor ters. As yet another alternative, Density Functional Per- Email address: mailhexu@gmail.com (Xu He) turbation Theory (DFPT) has been extended to magnetic Preprint submitted to Elsevier September 8, 2020 field perturbations by Savrasov [14] and yields the sus- from non-relativistic effects have not been integrated ceptibility and the magnetic exchange as a sub-product. yet. This perturbation has been implemented in abinit [15] Here, we present a Python package, TB2J, that al- and quantum Espresso [16], and has also been com- lows for automatic and systematic calculations of the bined with atomic displacements within non-collinear parameters of a Heisenberg Hamiltonian through the formalisms [17]. For several magnetic ions in the unit Korotin approach [34]. The script can compute the cell, a local perturbation DFPT scheme is presented in isotropic exchanges, the anisotropic exchanges, and the Ref. [18]. DMI from the output of the Wannier90 code or from a A different procedure, which avoids the problems of LCAO Hamiltonian (Siesta and OpenMX through the the total energy mapping, employs Green’s functions by SISL package [52], and GPAW) following the scheme taking the local spin rotation as a perturbation, as pro- proposed in Refs. [19, 24, 30]. TB2J is designed with posed in the seminal work of Liechtenstein, Katsnel- the goal to minimize the number of inputs and actions son, Antropov and Gubanov (LKAG) [19]. The Green’s from the user. In most cases, only the paths of the Wan- function method allows to determine the Heisenberg nier or DFT related files and the species of magnetic magnetic interaction parameters from the ground-state atoms are mandatory. Several types of output files are solution of the system, regardless of whether it is FM generated which can then be used in spin dynamics and or AFM. This method also allows access to a band- Monte Carlo codes. TB2J’s API is implemented in an by-band decomposition of the different magnetic in- abstract manner that eases interfacing with new tight- teractions, and it is easier to automatize than the to- binding-like Hamiltonians or spin dynamics simulation. tal energy mapping [20]. The method was also ex- We first present the general methods used in TB2J tended to correlated systems in Refs. [21, 22]. By and then describe how these methods are implemented. taking into account relativistic effects (spin-orbit cou- Next, we exemplify the usage of TB2J by applying it to, pling), the Dzyaloshinskii-Moriya interaction (DMI) BCC Fe, HCP Co bulk, SrMnO3 , BiFeO3 , and La2 CuO4 [23, 24, 25, 26] and the magnetic anisotropy can be crystals. In the end we discuss the advantages and limi- calculated as well [21, 22]. The method has been ex- tations of this method. tended to orbital-spin and orbital-orbital magnetic in- Throughout this paper, vectors are denoted with the~ teractions [27, 28]. Higher order terms in the Hamil- notation while matrices are represented with bold char- tonian, like the four-spin interaction or the biquadratic acters. term, were also calculated through this method [29, 30]. The LKAG method of computing magnetic interactions has been initially implemented with the Korringa-Kohn- 2. Formalism and algorithms Rostoker Greens function (KKR-Green) [31, 32] and tight-binding linear muffin-tin methods [33], and has The main idea of the method is to perturb a localized proven to be quite efficient. Its extension to other basis spin, in both the Heisenberg spin model and the DFT sets popular for DFT calculations will strongly broaden electronic model. For the latter, the rigid spin-rotation its usage. perturbation is done within the single-particle Green’s Constructing the Green’s function has proven to be function method. The Heisenberg parameters are then cumbersome for non-localized basis sets (e.g. plane mapped to the electron expressions [19]. We will con- waves), but this difficulty is greatly reduced by the use sider the quantization axis along z in the following. of Wannier functions (WF) as done by Kokorin et al. in [34]. The WFs can be constructed from first prin- 2.1. Heisenberg model ciples with the widely used open source code Wan- The Heisenberg Hamiltonian contains four different nier90 [35, 36]. Wannier90 has been interfaced with a parts and reads as large number of DFT codes, including Abinit [37, 38], Quantum Espresso [39], Siesta [40], VASP [41, 42], E = − ∑ Ki (~Si ·~ei )2 Wien2K [43], Fleur [44], Octopus [45, 46], OpenMX  i [47], GPAW [48], ELK [49], and many others. Sev- ~ ~ − ∑ Jiiso j Si · S j eral codes, including exchange.x [34], nojij [50], and Jx i, j [51], exist which calculate magnetic exchange param- eters from Wannier functions or linear combination of +~Si Jani~ ij Sj atomic orbitals (LCAO) DFT results. However, to the   best of our knowledge, DMI and anisotropic exchange +~Di j · ~Si × ~S j , (1) 2 Convention J˜i j indices, respectively. Due to translation symmetry, the tight-binding Hamiltonian, H, and the overlap, S, matri- − 12 ∑ J˜i j~Si · ~S j 2Ji j ces can be parameterized as ij − ∑ J˜i j~Si · ~S j 2Ji j Him jm0 σ σ 0 (~R) = hψimσ (~r)| H |ψ jm0 σ 0 (~r + ~R)i(3) , 1 ˜~ ~ Sim jm0 σ σ 0 (~R) = ~ 2 ∑ Ji j Si · S j −2Ji j hψimσ (~r)| |ψ jm0 σ 0 (~r + R)i , (4) ij ∑ J˜i j~Si · ~S j −2Ji j where ~r is the position inside the unit cell, ~R is the lat- tice vector, and H denotes the total Hamiltonian. The ∑ J˜i j~Si · ~S j −Ji j ij overlap matrix S reduces to the identity matrix when the basis functions are orthonormal. TB2J can use both Table 1: The conversion of Ji j to other conventions, where J˜ is the non-orthogonal LCAO basis sets and orthogonal Wan- exchange parameter in that convention. The notation < i j > means a nier basis sets. Below we discuss only the results for pair of i j without counting it twice. The DMI parameters ~D can be an orthogonal basis set. The non-orthogonal basis set converted in the same way. is discussed in Ref. [50], where it is shown that the ex- pressions for the exchange parameters are the same as where the first term represents the single-ion anisotropy for an orthogonal basis. (SIA), the second is the isotropic exchange, and the In the following, we drop all orbital and spin indices third term is the symmetric anisotropic exchange, where for simplicity, Hence Hi j is a sub-matrix of H contain- Jani is a 3 × 3 symmetric tensor. The final term is the ing all the spin and orbital components for atoms i and DMI, which is antisymmetric. Importantly, the SIA is j. We note that in DFT it is possible to do non-collinear not accessible from Wannier90 as it requires separately calculations with or without SOC (the exchange correla- the spin-orbit coupling part of the Hamiltonian [53]. tion potential can induce a non-collinear ground state). However, it is readily accessible from constrained DFT If one were to perform a calculation with SOC and the calculations [54]. We note that there are several con- spin constrained to one direction, it would be catego- ventions for the Heisenberg Hamiltonian, here we take a rized as collinear. In the collinear case, H is diagonal in commonly used one in atomic spin dynamics: we use a the spin subspace. minus sign in the exchange terms, i.e. positive exchange The Green’s function in reciprocal space is defined as J values favor ferromagnetic alignment. Every pair i j  −1 is taken into account twice, Ji j and J ji are both in the G(~k, ε) = εS(~k) − H(~k) , (5) Hamiltonian. The spin vectors ~Si are normalized to 1, so that the parameters are in units of energy. The other ~ ~ ~ ~ where H(~k) = ∑~R H(~R)eik·R , and S(~k) = ∑~R S(~R)eik·R . commonly used conventions differ in a prefactor 1/2 or The Green’s function in real space is obtained using the a summation over different i j pairs only. The conver- following expression sion factors to other conventions are given in Table 1. For other conventions in which the spins are not nor- Z ~ ~ malized, the parameters need to be divided by |Si | S j G(~R, ε) = G(~k, ε)e−ik·R d~k. (6) BZ in addition. From the total energy due to the spin interactions, Eq. In the following, we drop the ~R in the atom pair la- (1), we obtain the following variation with respect to the beled by i, j, ~R, e.g. Gi j means Gi j (~R), and G ji means ~Si and S~ j G ji (−~R). From Eq. (7) all equations are given in real space. δ Ei j = − 2Jiiso ~ ~ j δ Si · δ S j For each atom i, the intra-atomic component of H is − 2δ ~Si Jani ~ (2) defined as Pi = Hii (~R = 0) and is of size 2Norb × 2Norb . ij δSj Each Pi,mm0 is a 2 × 2 matrix in spin, which can be de- − 2~Di j · (δ ~Si × δ ~S j ) composed into its scalar and vector parts Pimm0 = p0imm0 I +~pimm0 · ~σ , (7) 2.2. Tight-binding Hamiltonian and Green’s function = p0imm0 I + pimm0~eimm0 · ~σ , We start from a generalized tight-binding Hamilto- nian. The localized basis functions are denoted as where ~pimm0 is the vector part of P, (upper case P de- ψimσ (~r) with i, m, and σ being the site, orbital, and spin notes matrices in spin and orbitals, lower case p is used 3 for matrices in orbitals only) which has the x, y, and z By putting this into Eqn. 12, we get components pximm0 , pyimm0 , pzimm0 . ~eimm0 is the unit orien- tation vector of ~pimm0 , and ~σ = (σx , σz , σz ) are the three δ Ei j = −2[A00 ij − ∑ Auv ei · δ~e j i j ]δ~ u=x,y,z Pauli matrices. In condensed form, the spin and orbital matrix for site i is now Pi = p0i I +~pi · ~σ . We can decom- −2 ∑ δ eui [Auv vu v i j + Ai j ]δ e j (13) u,v∈x,y,z pose the Green’s function for each inter-site orbital pair Gim, jm0 in the same way − 2d~i j · (δ~ei × δ~e j ) ~ im, jm0 · ~σ , Gim, jm0 = G0im, jm0 I + G (8) in which the 4 × 4 matrix Ai j is defined as Z EF ~ im, jm0 form the G0 and ~Gi j matrices, where G0im, jm0 and G 1 n o ij Auv ij = − Tr pzi Guij pzj Gvji dε, (14) respectively. π −∞ where u, v ∈ {0, x, y, z}, and the component of d~i j , diuj = 2.3. Magnetic force theorem Re(A0u u0 i j − Ai j ). Comparing Eq. (13) to Eq. (2), we can According to the force theorem the total energy vari- ation due to a small perturbation from the ground state find that the values of the isotropic exchange J iso , the coincides with the change of the single-particle energies anisotropic exchange Jani , and the DMI ~D can be ex- at fixed ground-state potential pressed as 00 yy Z EF Z EF Jiiso xx zz j = Im(Ai j − Ai j − Ai j − Ai j ), (15) δE = εδ n(ε) dε = − δ N(ε) dε (9) −∞ −∞ Jiani,uv j = Im(Auv vu i j + Ai j ), (16) where n(ε) = − π1 Im Tr(G(ε)) is the density of states Duij = Re(A0u u0 i j − Ai j ), (17) and N(ε) = − π1 Im Tr(ε − H) is the integrated density which is the same with Ref. [30]. of states. The traces are taken over orbitals only, not If spin-orbit coupling (SOC) is neglected, A0u u0 i j = Ai j , spin. Thus, the first order variation of N due to δ H can i.e. the DMI term is zero and both p and G only have be written as δ N(ε) = π1 Im Tr(δ HG); the second order components along the spin quantization axis, say the variation is δ 2 N(ε) = π1 Im Tr(δ HGδ HG). z direction. Then, the x and y components of ~p and Now we use the spin rotation as a perturbation. For ~G vanish. Thus, ~pi = (0, 0, pz ), Gx = Gy = 0, G0 = i the rotation of the spin at site i, the change in the energy 1 ↑ ↓ z 1 ↑ ↓ 2 (G + G ), G = 2 (G − G ). The isotropic exchange up to the second order reads as parameter then reduces to Z EF 1 00 δ Ei1spin = − Im Tr(δ Hi G + δ Hi Gδ Hi G) dε. Jiiso zz j = Im(Ai j − Ai j ). (18) π −∞ (10) Similarly, for the rotation of the spins at sites i and j, Defining ∆i = p↑i − p↓i = 2|~pi | = 2pzi we obtain from Eq. the change in energy is given by (18) the LKAG [19] expression for isotropic exchange Z EF 1 n o Im Tr ∆i G↑i j ∆ j G↓ji dε. Z EF δ Ei2spin =− 1 Im Tr [ δ Hi G + δ Hi Gδ Hi G Jiiso j =− (19) j π 4π −∞ −∞ (11) +δ H j G + δ H j Gδ H j G 2.4. xyz average strategy +2δ Hi Gδ H j G ] dε. We note that when all spins are oriented along one direction (e.g. z), the components with u = z or v = z in The energy variation due to the two-spin interaction is the J tensor are non-zero if the variation with respect to then the rotation is only kept to first order. The xz, yz, zx, zy, δ Ei j =δ Ei2spin − δ Ei1spin − δ E 1spin = and zz components of the anisotropic exchange, and the j j Z EF (12) Dz in the DMI cannot be obtained from a single calcula- 2 tion with magnetization along the z direction. It should − Im Tr(δ Hi Gδ H j G) dε. π −∞ be noted that this does not affect the properties too much if the system stays close to the reference spin state, as The change of H due to the rotation of spin is δ~φ ×~p these terms are small and higher order in the spin ro- with the rotation axis along δ~φ and the angle δ~φ . tation angles. For example, in a AFM with small spin 4 canting system, the parameters obtained from the AFM 3. Implementation reference state can be used directly to get the canting angles. To determine the missing parameters, we ro- The TB2J package is implemented with three main tate the whole system and lattice, (equivalently one can modules, the tight-binding (TB) module, the exchange rotate the quantization axis) from z to the x and y direc- module, and the output module (see Fig. 1 for a tions and apply the back-rotation to the DMI vectors. A schematic view). The TB module provides the inter- weighted average over all three directions is taken: the face with external DFT and Wannier function codes and inaccessible components are given a weight of 0 and the calculates the Green’s function. The exchange module rest have the same weight. This small trick allows not uses the Green’s function to calculate the magnetic in- only to obtain the full set of Jani ~D vector components teraction parameters. The output module writes these , but also to reduce the numerical noise. The DMI is parameters to output files and provides a python API usually much smaller than the isotropic exchange, such from which the magnon band structure can be calcu- that its calculation is more numerically delicate. This is- lated. sue arises in particular when Wannier functions are used The TB module defines the classes of the TB model whose symmetry is not guaranteed, and where the dis- and the Green’s function. The internal TB module can entanglement procedure can introduce numerical errors. read the Wannier90 output data files to build the Hamil- A direct approach to calculate the Dz term is pro- tonian. The WF’s are assigned to the nearest atom with posed in Ref. [25], taking the spin-rotation perturbation the crystal periodicity taken into account. A warning is to higher order. This is not implemented in TB2J, and issued if a WF is far away from any atom. The code uti- would also be more sensitive to numerical noise. lizes the “duck type” feature: a class can be plugged in if the specific required methods are defined. The 2.5. Higher order terms methods include the calculation of H(~k) and the eigen- It has been shown by several authors that the pa- values and eigenvectors. Hence, external libraries can rameters from the Green’s function method cannot be be wrapped easily. An interface to the sisl library [52] mapped exactly onto a bilinear J tensor [29, 55, 30] has been implemented, which includes the tight-binding and higher order terms need to be considered, which in- model from Siesta and OpenMX. We also wrap directly clude multi-spin interactions and the higher-order two- the python based DFT code GPAW in LCAO mode in- spin terms. In the simplest case, following Ref. [30], the stead of reading the Hamiltonian from the output files. Hamiltonian can be written in a biquadratic form as Starting from the eigenvalues and eigenvectors, the Green’s functions can be calculated without invert- ing the εS(~k) − H(~k) matrix for every ε. Instead,  2 H Q = − ∑ Ji0 j~Si~S j − ∑ Bi j ~Si~S j , (20) i, j i, j the Green’s functions are calculated with G(~k, ε) = h i−1 where Ji0 j and Bi j are determined by the A matrices as Ψ(~k) εI − Diag(E(~k)) Ψ(~k)† , where E(~k) and Ψ(~k) are the eigenvalues and the eigenvector matrices, re- Ji0 j = A00 zz i j − 3Ai j , (21) spectively. The real-space Green’s function can then be Bi j zz = Ai j . (22) calculated by Fourier transforming G(k). Once the Green’s functions are calculated, the ex- When ~Si and S~ j are close to their reference values, we change module calculates the magnetic interaction pa- have rameters using Eqs. (15)-(17). The Hamiltonian and d2H Q the Green’s functions are first decomposed into their Ji j = − = Ji0 j + 2Bi j~Si~S j (23) {0, x, y, z} n components. oThen, the elements of the 4 × 4 ~ ~ d Si d S j matrix, Tr pi Guij p j Gvji for each ε are calculated and ' Ji0 j + 2Bi j~Siref~Sref 0 00 zz j = Ji j + 2Bi j = Ai j − Ai j , later integrated to obtain the Ai j matrix. A contour integration method is used for the EF dε R which is equivalent to the Ji j expression with only bi- linear terms considered, and it shows that the effective integration. The range of the integration is (Emin , EF ), bilinear J term depends on the orientation of Si and S j . where Emin is either below the lowest band energy or This method was proposed in order to improve the de- chosen such that the orbitals below Emin have only negli- scription of the system when the deviation from the ref- gible interaction with those near EF . By default, a semi- erence state is large [30]. From our limited experience, circle path is used, which is centered at (Emin + EF )/2 however, the parameters produced by Eqs. (21) and (22) on the real axis and has a radius of (EF − Emin )/2, go- are not always physical and their fit is more complex. ing through the upper half of the complex plane with 5 DFT TB2J Abint Quantum Espresso TB Exchange VASP Wien2K Collinear Fleur Wannier90: Magnon Isotropic J Band ... TBGreen Bi-quadratic H(k) LCAO H(R) Spin Dynamics Siesta Sisl: Non-collinear IO Eigen(k) OpenMX Isotropic J G(k) Bi-quadratic txt G(R) netcdf DMI -Multibinit xml -Tom's GPU ASD Anisotropic J User's format -Vampire Effective TB Hamiltonian Figure 1: Schematic diagram view of TB2J workflow. Im(ε) > 0. We show how to install and use the code to calculate the The magnetic interaction parameters are calculated magnetic interaction parameters using body centered from the Aij matrix. The full list of quantities that are cubic (BCC) Fe, and hexagonal close packed (HCP) Co, available from TB2J is given in Table 2. BiFeO3 , and La2 CuO4 as examples. Quantity Expression 4.1. Installation collinear case TB2J can be installed with a simple command if n o Ji j 1 R EF − 4π −∞ dε Im Tr ∆i G↑i j ∆ j G↓ji python and pip environments are pre-installed Ji0 j Im(A00 zz i j − 3Ai j ) pip install TB2J Bi j Azz ij Alternatively, one can download the package from non-collinear case https://github.com/mailhexu/TB2J and run yy zz Ji j Im(A00 xx i j − Ai j − Ai j − Ai j ) python setup . py install D~uij Re(A0u u0 i j − Ai j ) to install the package. Jani J ani,uv = J ani,vu = Im(Auv + Avu ) By default, TB2J only installs the hard (non optional) dependencies automatically. Hence, the sisl package, Table 2: Summary of magnetic interaction quantities that are calcu- which is used to read the Hamiltonian from the Siesta or lated within TB2J. All come from Ref. [30] except the first (Ref. [19]) and the DMI, which was first written in Ref. [24] 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. TB2J user instructions 4.2. Computing the magnetic interaction parameters In this section, we describe the practical usage of TB2J. The users can also refer to the online doc- In this section we describe the TB2J procedures to umentation at https://tb2j.readthedocs.io/en/ compute the magnetic interactions, from an electronic latest/, which is more detailed and regularly updated. structure calculation to the final output file. 6 4.2.1. Preparation of electronic structure and tight- tional symmetry and modifies the Hamiltonian accord- binding Hamiltonian ingly. The WF Hamiltonian and the center positions To obtain the magnetic interaction parameters, the need to be output by Wannier90 using the “write hr” first step is to do a converged DFT calculation of a and the “write xyz” input flags. magnetic crystal, at the collinear or non-collinear level. For DFT codes based on LCAO basis sets, such Preferably, this calculation treats the magnetic ground as Siesta, the Hamiltonian is already localized and state of the system, however, the magnetic interaction needs no further transformation. For calculating parameters can also be calculated for a different refer- the parameters of the Heisenberg Hamiltonian only ence state. It should be noted that the spin quantization the localized DFT Hamiltonian and the overlap ma- axis is presumed to be along the z axis throughout TB2J. trix need to be saved. For example, one can use the options “CDF.Save=True”, “SaveHS=True”, and The next step is to construct the tight-binding Hamil- “Write.DMHS.Netcdf=True” in Siesta (version 4.x) to tonian. For DFT codes using non-local basis sets, like enable the saving of these matrices. plane waves, the WFs can be constructed for any code which has an interface with Wannier90. All the spin- 4.2.2. Running TB2J polarized orbitals of all atoms contributing to the mag- TB2J has two Python executables: wann2J.py and netic interaction should be carefully selected to compute siesta2J.py, for calculating J from Wannier90 and accurately the magnetic interactions interaction parame- Siesta output, respectively. A similar script named ters. For example, in the case of a transition metal oxide openmx2J.py is in the TB2J OpenMX package. These the d orbitals of the transition metal cation should be in- scripts are designed to have a minimal user input where cluded, but also the oxygen 2p orbitals that are involved only the paths to the files containing the electron Hamil- in the superexchange or DMI mechanism through hy- tonian information and defining the magnetic atom bridization with d orbitals. This forms the minimal basis species need to be provided. of orbitals to be included in the construction of the WFs but other orbitals might be important as well. The num- With Wannier90. The executable script wann2J.py can ber of orbitals to be included is system dependent and be used with the Wannier90 output files. For a non- should be checked by the user (convergence of the cal- collinear calculation, the spin up and spin down channel culated magnetic interactions with respect to the num- Wannier functions need to be present. In addition, the ber of Wannier orbitals). For example, in the case of following parameters need to be specified: SrMnO3 the Mn-3d, O-2p orbitals are necessary. The options to build the Wannier functions should be en- • Whether the calculation is collinear or non- abled in the DFT codes. For example, in Abinit, the collinear, given by the –spinor option. TB2J as- “prtwant 2” and “w90iniprj 2” options (ABINIT ver- sumes a collinear calculation by default, and the sion 9.x) are needed for the calculation of maximally usage of –spinor specify that the calculation is localized Wannier functions (MLWFs). The input file non-collinear, where the Hamiltonian is in a spinor for Wannier90 also needs to be present in the execution form. directory. THe quality of the Wannier function Hamil- tonian is need to be checked. The magnetic interaction • the prefix to the paths of up and down Wan- parameters are often meV or µeV, which requires the nier functions, (–prefix up abinito w90 up – noise in the Hamiltonian to be lower. prefix down abinito w90 down, for non-collinear calculations, and –prefix spinor for non-collinear The rigid spin rotation on one site is performed by calculations). The filename of the Hamiltonian is rotating all the spins of the WFs associated to a given the prefix plus “ hr.dat”. atom, which requires the WFs to be centered on an atom or at least very close to it. TB2J uses the Wannier cen- • The posfile specifies the name of a file contain- ters to decide which atom each WF “belongs” to. As ing the atomic structure and cell parameters. ASE a result, using the MLWFs [56] might not always be file formats are also readable (the full list of which the best choice. Projected WFs or selectively localized can be found on https://wiki.fysik.dtu.dk/ WFs [57], which add a constraint on the Wannier cen- ase/ase/io/io.html). It should be noted that ters, can be used instead. The Wannier centers for one some ASE formats cannot be used, e.g. xyz, be- atom might be located closer to a periodic copy of that cause they do not contain the cell parameters which atom than the original atom. To solve this problem, are required by TB2J. We also recommend using TB2J shifts the Wannier centers by using the transna- formats with the cell matrix rather than only the 7 (a, b, c, α, β , γ), which will cause trouble if there • emin, emax: the integration lower and upper are anisotropic or DMI terms as they are not rota- bounds, relative to the Fermi energy. The value tionally invariant. of emin is automatically determined if not given, emax should be about 0 for metallic systems, • The Fermi energy in units of eV. whereas it lies in the band gap for insulating sys- • The type of magnetic elements (symbol from the tems. periodic table). • rcut: the cutoff or max distance between spin pairs. Here is an example for calculating the Js in the non- collinear case: The full list of options can be accessed using “–help”. As discussed in the previous section, the z compo- wann2J . py -- posfile abinit . in -- efermi nent of the DMI, and the xz, yz, zx, zy, zz components 5.8 -- elements Fe -- prefix_up abinito_w90_up -- prefix_down of the anisotropic exchanges are non-physical, and an abinito_w90_down xyz average is needed to get the full set of magnetic interaction parameters. In this case, scripts to rotate In the case of a non-collinear, the WF Hamiltonian is in the structure and merge the results are provided, they a single file and we need to specify the calculation type: are named TB2J rotate.py and TB2J merge.py. The wann2J . py -- posfile abinit . in -- efermi TB2J rotate.py reads the structure file and generates 5.8 -- elements Fe -- prefix_spinor three files containing the z → x, z → y and the non- abinito_w90 -- spinor rotated 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 for- With Siesta. Only a minimal set of parameters is mat of the output structure files is provided using the needed for Siesta: the filename of the input for the “–format” parameter. An example for using the rotate Siesta calculation. The other information needed, in- file is: cluding whether the calculation has SOC enabled and the Fermi energy, are found in the Siesta results: TB2J_rotate . py BiFeO3 . cif -- format cif siesta2J . py -- element Fe -- input - fname = ’ siesta . fdf ’ The user has to perform DFT single point energy cal- culations for these three structures in different directo- ries, keeping the spins along the z direction, and run With OpenMX. The interface to OpenMX is distributed TB2J on each of them. After producing the TB2J re- as a plugin to TB2J called TB2J OpenMX under the sults for the three rotated structures, we can merge the GPL license, which need to be installed separately, be- DMI results with the following command by providing cause code from OpenMX which is under the GPL li- the paths to the TB2J results of the three cases, e.g.: cense is used in the parser of OpenMX files. TB2J_merge . py BiFeO3_x BiFeO3_y BiFeO3_z pip install TB2J_OpenMX -- type structure In the DFT calculation, the ”HS.fileout on” options A new TB2J results directory is then made which con- should be enabled, so that the Hamiltonian and the over- tains the merged final results. lap matrices are written to a ”.scfout“ file. The neces- sary input are the path of the calculation, the prefix of the OpenMX files, and the magnetic elements: 4.2.3. Output files In the following we describe the output files which openmx2J . py -- path ./ -- prefix openmx -- TB2J produces. By running wann2J.py or siesta2J.py, a elements Fe directory with the name TB2J results will be generated, which contains the following output files: General options. There are several tunable parameters, which the user usually does not need to specify. Among • exchange.out: A human readable output file, which them, there are: summarizes the results. • nz: The number of steps in the path of the contour • Multibinit: A directory containing output which integration. 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 BiFeO3 with SOC enabled, calculated with spin along the z axis. The exchange.out file contains three sections: cell, is not necessarily the same. However, for localized d atoms and exchange. The cell section contains the lat- and f orbitals the magnetic moments should be close to tice parameter matrix. The atoms section contains the their DFT counterparts, for TB2J to yield correct results positions, charges (for verification) and magnetic mo- for the parameters. Large differences between the TB2J ments of the atoms: see Listing 1. and DFT values indicate that something may have gone wrong: either in in the contour integration EF dε used R Here, the charge and magnetic moment of each atom in TB2J, or in the construction of the Wannier functions are only integrated with the WFs attached to this atom. (incorrect wannierization process or too small WF ba- As such they can differ from the quantities coming from sis set). Often, it comes from excluding an orbital that the direct DFT output, as not all bands are used in the is important for the magnetic interaction in the studied construction of WFs. The WF charges should be inte- system. In the case of metallic system, the Fermi energy gers, for LCAO the values depend on the band energy might have to be slightly shifted with respect to the DFT cutoffs. In addition, the exclusion of very deep lying reference due to different numerical method used in the levels from the calculation of J can also lead to devia- integration of charge density. tions in the charges which might appear both for WFs and for LCAO. Another source of difference between Each pair of atoms is labeled by three parameters, the the TB2J charges and magnetic moments and the DFT index i, j and R, where i and j are the indices in the unit ones is the integration volume around the atoms, which cell. The vector ~R specifies the cell the atom j is trans- 9 lated to, i.e. the reduced positions of the two atoms are where the qpath option specifies the q-point path. If it ~ri and ~r j + ~R, respectively. By default, the interaction is not specified, an automatically generate qpoint-path is calculated within a supercell corresponding to the k- will be used from the information of the lattice struc- mesh. For example, with a 7 × 7 × 7 k-mesh, all i j pairs ture by using ASE[58]. The details of the q-point path will be produced for the spin labeled i in the center cell can be found on https://wiki.fysik.dtu.dk/ase/ of a 7 × 7 × 7 supercell. With the “rcut” flag, only the ase/dft/kpoints.html. parameters for i j pairs within a distance of rcut are cal- A file, named “exchange magnon.pdf”, containing culated. The exchange parameters are reported as fol- the magnon band structure will be generated. The low- lows: magnetic atom i connected with magnetic atom j, est energy eigenvector of J(~q) is printed, from which, R is the lattice vector between the unit cells containing together with the ~q vector, one determines the ground- i and j, the value of J for this pair of magnetic atoms in state spin configuration. meV, the vector connecting them, and the distance be- tween the pair of atoms. If SOC is enabled, the DMI 4.3. Examples and anisotropic J ani parameters are given in addition. The DMI vectors ~D and the anisotropic J ani are printed In order to demonstrate the usage of TB2J, we use as vectors and matrices, respectively. five different examples, body centered cubic (BCC) Fe, Apart from the main exchange.out file, TB2J deliv- hexagonal closed packed (HCP) Co, SrMnO3 , BiFeO3 , ers several other outputs, which provide the input for and La2 CuO4 . The BCC Fe and HCP Co systems are spin dynamics (SD) and Monte Carlo (MC) simula- among the most studied magnetic materials and are used tions. TB2J is interfaced with several SD and MC as standard benchmarks. The other three materials were codes. It has native support to the Multibinit code de- chosen to represent a wide range of properties: SrMnO3 livered as part of the Abinit code since version 9.0 [38]. as a prototype structure for superexchange, BiFeO3 as The TB2J results/Multibinit directory contains the tem- a multiferroic material [61], and La2 CuO4 as a layered plates of input files for this code. One can usually run perovskite with a canted spin structure. spin-dynamics with slight or no modification of these files. “Experimental” inputs are also generated for Vam- 4.3.1. BCC Fe pire [59] and Thomas Ostler’s GPU-ASD code [60]. As a first benchmark, we calculate the magnetic inter- The Output module of TB2J provides a versatile API, action parameters for the probably most studied mag- described in the online documentation (see Discussion netic structure, BCC Fe. The isotropic exchange is section), which makes it easy to generate data for any calculated with TB2J-Siesta and compared with results other code. from the TB-LMTO method in Ref. [62], and KKR method in Refs. [60] and [30]. We do not show TB2J- Wannier90 results because the localization procedure 4.2.4. Magnon band structure creates Wannier functions centered on bonds instead of Once the J have been obtained in real space from atoms. Then, the mapping to a Heisenberg model is not TB2J, the magnon dispersion curves Emagnon (~q) can be well defined and the resulting exchange parameters are obtained from the Fourier transform of J to the ~q space: unphysical. ~ The DFT calculations were performed with the Siesta J(~q) = ∑ J(~R)ei~q·R , (24) code [40] Max-1.0.13 version, with a double-zeta- ~R polarized numerical atomic orbital basis set. The In the simple FM case this reduces to diagonalizing GGA-PBE [63] functional and the norm-conserving 4 pseudopotentials from the pseudo-dojo [64] “standard” M (J(0) − J(~q)), where M is the one site magnetic mo- ment. The magnon band structure for more complicated dataset in the psml [65] format were used. A 9 × 9 × 9 magnetic configurations is more complex. Also, multi- k-point grid was used to sample the Brillouin zone. ple types of magnetic sites cannot be treated at the mo- The results in Fig. 2 (a) show good agreement (within ment. Both features will be added in a future version of a few percent for the largest exchange parameters) be- TB2J. tween our calculation and those in Refs. [62] and [60] With TB2J, the magnon band structure is obtained us- obtained with the GGA-PBE exchange correlation func- ing the following command: tional. The exchange parameters from the KKR method in Ref. [30] are, however, significantly smaller. This TB2J_magnon -- qpath GNPGHP -- show might be attributed to a different density functional in the underlying DFT calculation, which is not reported in 10 20 (a) TB2J-Siesta TB-LMTO 15 KKR1 KKR2 J (meV) 10 5 500 0 400 Energy (meV) 60 300 (b) J 0 TB2J-Siesta B: TB2J-Siesta 40 J 0: KKR2 200 B: KKR2 J (meV) 20 100 0 0 N P H N 20 0.5 1.0 1.5 2.0 2.5 3.0 rij/a Figure 3: The magnon band structure of BCC Fe as obtained from TB2J. Figure 2: Exchange parameters for BCC Fe. (a) J from TB2J-Siesta, TB-LMTO (Ref. [62]), and KKR (KKR1: Ref. [60], and KKR2: Ref. [30]). (b) The bilinear (J 0 ) and biquadratic (B) exchange of BCC Fe using TB2J-Siesta, and KKR-Green method (Ref. [30]). ri j is the distance between the Fe atom pairs, and a is the cubic cell parameter. Ref. [30]. The decomposition of the exchange into bi- 0.03 TB2J-Siesta linear and biquadratic terms is also calculated, see Fig. 2 0.02 KKR (b), and compared with the KKR method from Ref. [30]. 0.01 The J 0 and B from our calculation follow the same trend Dx (meV) but are again larger than those in Ref. [30]. 0.00 Fig. 3 shows an example of the magnon band struc- 0.01 ture of BCC Fe. As we can see, the ground state corre- 0.02 sponds to a q = Γ eigenvector, which identifies a ferro- 0.03 magnetic structure. The energy minimum and the cor- responding eigenvector are also printed to the standard 0.03 output: 0.02 0.01 Dy (meV) The energy minimum is at : q = [0. 0. 0.] 0.00 The ground state eigenvector is : Fe1 : [1.0 0.0 0.0] 0.01 0.02 0.03 4.3.2. HCP Co 1.0 1.5 2.0 2.5 3.0 As a second example, we calculate the DMI in HCP rij/a Co with TB2J-Siesta and compare with the results from the KKR method of Ref. [26]. A DZP LCAO basis set, Figure 4: The DMI parameters of HCP Co calculated with TB2J com- together with a 9 × 9 × 9 k-point mesh, the GGA-PBE pared with results calculated with the KKR-Green method reported in functional [63], and the standard accuracy of v0.4 norm- Ref. [26]. conserving pseudopotentials from the Pseudo-Dojo [64] are used. The SOC is enabled. The DMI in HCP Co is small (only of the order of 0.01 meV), which is a chal- 11 |S| 1NN 2NN 3NN 4.3.4. BiFeO3 BiFeO3 has a cubic-perovskite-like structure, in TB2J-W90-Abinit 2.81 -8.16 -0.36 -0.02 which the Bi atoms are at the corners, the O atoms TB2J-Siesta 2.85 -7.70 -0.02 0.11 are at the face centers and the Fe ions are in the cen- TE [70] -6.98 -0.36 -0.01 ter of the cube. The R3c state structure has two units Table 3: The exchange interaction parameters of SrMnO3 , calculated of BiFeO3 , and can be seen as the cubic structure plus with TB2J-Wannier90-Abinit (TB2J-W90-Abinit) and TB2J-Siesta, the “a− a− a− ” rotation [71] of the O octahedra plus the compared to the results of the total energy (TE) method in Ref. [70]. polar distortion along the [111] direction. The DMI is |S| and the exchange parameters are given in units of µB and meV, related to the local inversion symmetry breaking due to respectively. Ref. [70] uses a different convention and a |S|2 = (3/2)2 is multiplied to the values in it to convert the values to the convention these distortions [54]. of the present paper. They also use a slightly smaller U(Mn)=2.7 eV. The DFT results are from Siesta Max-1.0-13 ver- sion, with a DZP LCAO basis set for all atoms. We used the full-relativistic norm-conserving pseudopoten- lenge to the accuracy of its calculation independent of tials from the Pseudo-Dojo ”standard” dataset, with the the method used. The results in Fig. 4 show good agree- GGA-PBE [63] density functional, corrected by a Hub- ment between the two methods (within 10% for most bard term[69] with U(Fe) = 4, and a 7 × 7 × 7 k-point values). The Dz component is not accessible when the ~S mesh. The SOC is enabled for the calculation of the point along z. anisotropic and DMI terms. The calculated magnetic interaction parameters are given in Table 4.3.4. The first nearest neighbor ex- 4.3.3. SrMnO3 change and DMI are compared with both a DFT total In order to compare different DFT inputs, we calcu- energy method from Ref. [72] and experimental result late the exchange parameters for SrMnO3 starting from from Ref. [73]. In Ref. [72], the DFT+U [69] correc- either the Wannier90 or Siesta results. SrMnO3 has a tion is used with U(Fe) = 4. The results show that cubic perovskite structure with a lattice parameter of the isotropic exchange and DMI parameters from TB2J- 3.81 , and is a collinear antiferromagnet. The G-type Siesta and the total energy method are close. The first antiferromagnetic ground state with a magnetic wave- NN J and DMI also agree well with the experimental vector of (1/2, 1/2, 1/2) is taken as the reference con- value. figuration in our DFT calculations. It requires a 10-atom supercell expanded along the [111] direction. 4.3.5. La2 CuO4 The DFT calculation for constructing the Wannier In this section we show an example of calculating the functions is performed with Abinit v9, using a plane- magnetic interaction parameters for a spin-canting ma- wave basis set with a cutoff of 30 Ha. The projected terial, namely La2 CuO4 . It has a layered structure where augmented wave (PAW) [66] potentials are from the the Cu-O octahedra are separated by the La-O planes, “JTH” dataset [67] v1.2. The Mn 3d and O 2p orbitals see Fig. 5. There is an out-of-plane displacement for are considered in the construction of the maximally lo- each of the O ions in the Cu-O plane (xy plane) due to calized Wannier functions from which the tight-binding the a− b− rotation of the Cu-O octahedra. Hamiltonian is calculated. In the Siesta calculation, we The DFT calculations were performed with OpenMX use a DZP LCAO basis set for all atoms and the pseu- [47] using the numerical pseudo-atomic orbital (PAO) dodojo v0.4 norm-conserving pseudopotentials [64]. In basis set, with La s3p2d2f1, Cu s3p2d1, and O s3p2d1 both DFT calculations, we use the PBEsol [68] func- PAOs as provided in the OpenMX database [74] v2019. tional to approximate the exchange-correlation, with a A 6 × 6 × 4 k-point grid and the GGA-PBE [63] func- Hubbard U [69] of 3 eV on the Mn atom. The SOC is tional with a Hubbard U correction [75] (U(Cu)=8 eV, not enabled. A 9 × 9 × 9 k-point mesh is used for the J(Cu)=0.8 eV) are used. The primitive cell of 4 × integration of the Brillouin zone. La2 CuO4 with an antiferromagnetic (AFM) spin align- We show the calculated magnetic interaction param- ment in the xy plane is taken as the reference DFT state. eters with both Abinit+Wannier90 and Siesta in Table 3 The isotropic exchange and DMI parameters are calcu- and compare them with total energy results from Ref. lated. The interactions with the first nearest neighbors of [70]. We can see that all three methods yield similar one Cu atom are listed in Table 5. The results are com- qualitative results with a large and negative parameter pared with the LMTO results from Ref. [25], in which a for first nearest neighbours and small parameter for sec- U(Cu)=10 eV and J(Cu)=1 eV are used. The exchange ond and third neighbour interaction. interaction shows an AFM alignment. The sum of the ~D 12 Method nth NN (Ri j ) (J iso ) (Jxy ) (Jyz ) (Jxz ) (Dx , Dy , Dz ) |D| TB2J-Siesta 1 (1,0,0) -26.79 0.003 0.001 -0.008 (0.158, -0.093, 0.325) 0.373 Exp[73] -20.25 0.507 TE[72] -18.296 (0.131, -0.087, 0.363) 0.396 TB2J-Siesta 2 (1,-1,0) -0.742 0.004 -0.001 -0.008 (-0.003, -0.009, -0.003) 0.010 TE[72] -0.603 (-0.003, -0.006, -0.006) 0.009 TB2J-Siesta 3 (1,1,0) -0.023 -0.000 0.012 -0.013 (-0.036, 0.003, 0.034) 0.049 TE[72] -0.009 (-0.000, -0.006,-0.012) 0.013 Table 4: The magnetic interaction parameters (in meV) of BiFeO3 compared with total energy (TE) method and experimental (Exp) measurement. In Refs. [72] and [73], a different convention is used, a −|S|2 /2 = −(5/2)2 /2 is multiplied to adapt to the convention used in this work. has a non-zero component along x, which results in a net magnetic moment perpendicular to x. The angle of ro- tation can be estimated as | ∑ j,i ~Di j |/|2 ∑ j,i Ji j |, which is about 0.8 × 103 , close to 0.7 × 103 in Ref. [25]. 4.4. Spin dynamics 1e4 4 Cv (a.u.) 2 0 15 T (K) 10 (a.u.) 5 0 1.0 T (K) M/Ms (a.u.) 0.5 0.0 50 100 150 200 250 300 T (K) Figure 5: The structure of La2 CuO4 . La atoms are depicted in green, Figure 6: The volume specific heat, CV (top), the susceptibility, χ Cu atoms in blue, and O atoms in red. The figure is in the x − z plane (middle), and the ratio of magnetic moments, M/MS (bottom), as with x and z in the horizontal and vertical direction. functions of T in SrMnO3 as obtained with MULTIBINIT spin dy- namics. The units are atomic units. We use a simple example of calculating the Néel tem- perature (TN ) of SrMnO3 to show the usage of the TB2J results in spin dynamics with the MULTIBINIT code. In the TB2J output, the templates of the Multibnit in- put files (mb.files and mb.in) are already generated. To 13 TB2J-OpenMX LMTO [25] ~R0 j J ~D J ~D (-0.49,0.5,0) -18.60 (-0.023, 0.025,0.001 ) -14.576 (-0.020, 0.032, 0.005) (0.49, 0.5, 0) -18.60 (-0.023, -0.025, -0.001) -14.576 (-0.020, -0.032, -0.005) (0.49, -0.5, 0) -18.60 (-0.023, 0.025, 0.001) -14.576 (-0.020, 0.032, 0.005) (-0.49, -0.5, 0) -18.60 (-0.023, -0.025, -0.001) -14.576 (-0.020, -0.032, -0.005) (0, 1, 0) 5.82 (0,0,0) 2.071 (-0.002, 0,0) (0, -1, 0) 5.82 (0,0,0) 2.071 (-0.002, 0,0) (-0.98, 0, 0) 4.96 (0,0,0) 1.943 (0.007, 0,0) (0.98, 0, 0) 4.96 (0,0,0) 1.943 (0.007, 0,0) Total -52.84 (-0.092,0,0) -48.476 (-0.070, 0, 0) Table 5: The isotropic exchange and DMI parameters for La2 CuO4 . The neighbors of the Cu1 atom shown in Fig. 5 is shown. The ~R0 j is the vector connecting the neighbors with Cu1. are in the unit of lattice vectors. The units of J and ~D are meV and are in laboratory coordinates. The convention used in Ref [25] is different from the one in the present work by a minus sign. calculate TN , we simply adjust a few parameters, e.g. close to the ground state, and fitting to extreme config- the size of the supercell and the temperature range and urations, e.g. Stoner excitations, will produce worse dy- steps. Then, by running namics. Even up to the Curie temperature, perturbative calculations typically perform well. multibinit < mb . files The Green’s function method is most accurate around the spin dynamics at various temperatures are per- the magnetic phase for which the DFT calculation has formed. For the parameters, we employ the ones from been performed. This is in contrast with the total-energy the TB2J-Siesta calculation. A file named mb.out.varT methods in which the parameters are averaged over di- contains the magnetic volume specific heat, CV , the sus- verse magnetic configurations. However, experience ceptibility, χ, and the ratio of the magnetic moment to shows that the two methods yield similar results for the the saturated value, M/MS . We plot the results as func- parameters for many systems. tions of temperature, T , in Fig. 6, from which we can An important caveat is that the parameters from total- see that TN is about 240 K, which is close to the experi- energy fitting methods are supercell-dependent. For ex- mental value of about 233 K. [76]. ample, if the supercell only contains the second nearest neighbours (2NN), the magnetic interaction of longer range is wrapped into the 1NN and 2NN parameters, 5. Discussion and large cells have to be built to check the conver- gence of the fitting, which is not always numerically 5.1. Advantages of the Green’s function method feasible. In the case of TB2J, the building of super- The Green’s function method implemented in TB2J cells is avoided. This makes the method suitable for has several advantages over the total-energy mapping large scale calculations, where the exchange magnetic methods. It is more robust, and only requires a DFT interaction parameters can be calculated at any distance calculation for a single spin configuration, typically the in real space, provided the k-point sampling is dense ground-state one. Therefore, it avoids the problem of enough. converging DFT calculations for magnetic structures TB2J can be used either in conjunction with WFs or which are much higher in energy or even unstable. For with LCAO codes. For WFs, TB2J works almost as a example, a structure with an insulating ground state may black box. Indeed, the whole procedure - from DFT become metallic in some spin configurations (often the data to the calculation of magnetic interaction param- FM phase is metallic while the AFM is not), or electrons eters - requires little or no user input. In the case of become delocalized, which modifies or destroys the LCAO codes, TB2J also runs with almost no extra in- magnetic character. The convergence of non-collinear put. As a result, TB2J can be easily integrated into high- cases with SOC can be especially difficult in DFT, and throughput workflows for automatic high-throughput avoiding several such calculations is a strong advantage. calculations of the magnetic interactions, using pack- Spin dynamics at moderate temperatures typically stay ages like ASE [58], Abipy [38] or Aiida [77], since 14 TB2J can work as a python library. or when the non-bilinear magnetic interactions are non- The construction of MLWFs was a sophisticated task negligible. The rotation of the spins is described in a in the past, since input parameters would vary for each basis set, either Wannier functions, or numerical atomic structure and had to be tuned by hand. However, the sit- orbitals. In some cases, especially when the spins are uation has been largely changed thanks to the recent de- not localized at or close to the centers of the basis func- velopment of several methods, in particular the selected tions, the result can be basis set dependent. columns of density matrix (SCDM) algorithms [78, 79]. With these advances, fully automated building of WFs 5.3. Code availability has been demonstrated [80]. The code is freely available under the BSD 2 The output of TB2J is directly usable by spin simula- clause license and can be found at https://github. tion codes. Therefore, the full procedure from the DFT com/mailhexu/TB2J/. Documentation is provided at exchange calculation to the spin system simulation can https://tb2j.readthedocs.io/en/latest/. The be easily automated. interface to OpenMX is in a separate package, 5.2. Limitations TB2J OpenMX (https://github.com/mailhexu/ TB2J-OpenMX), under the GPLv3 license, which use In the limit where the method implemented in TB2J TB2J as a library. Contribution to the code is wel- applies (rigid spin rotation approximation, magnetic come. In particular, we are happy to integrate interfaces moments localized on the atoms, etc.) and that the limi- to additional codes, be it on the input side (other first- tations of the DFT for magnetic systems (selection of a principles or tight-binding codes), or the output side correct enough exchange-correlation functional, choice (e.g. outputs to atomic spin simulation packages). of the values of DFT+U parameters, etc [81, 82, 83]) the calculation of the DMI and anisotropic exchange param- eters through the Wannier basis set might still be prob- 6. Conclusions and perspectives lematic. The user should always verify that the elec- tronic band structure obtained with the Wannier func- In this paper we present the TB2J python package tions is in good agreement with the DFT one (choice of which calculates the magnetic interaction parameters good disentanglement energy window, choice of projec- from Wannier and LCAO Hamiltonians, using a Green’s tors, etc). However, even if that is the case, the Wannier- function method. The isotropic and anisotropic ex- ization process itself can introduce noise of the order of change, and the DMI can be systematically calculated magnitude of a few µeV. For small quantities, like the by TB2J. The code can use results from a large number DMI or anisotropic exchanges, which are of the same of first-principles DFT codes, either through the Wan- size as the noise, the parameters cannot be determined nier90 interface (for plane wave codes but also other ba- with sufficient resolution. The problem occurs, in par- sis sets) or directly from LCAO codes (SIESTA, GPAW, ticular, when strong disentanglement of bands is neces- and OPENMX). TB2J requires a minimal number of in- sary, i.e. when no gap is present between selected bands. put parameters, such that it can be easily integrated in It can be reduced in cases where the bands of interest are high-throughput workflows. One of the most appeal- not entangled too much with higher unoccupied bands ing features of TB2J is the requirement of only one unit by including more unoccupied bands in the Wanneriza- cell DFT calculation (or eventually 3 for numerical av- tion process. Therefore, a validation of the quality of eraging of x, y, z spin orientations) to evaluate the mag- Wannier function Hamiltonian should be performed to netic interactions at any distance between the magnetic make sure that the calculation is meaningful, by com- atoms. We hope this development will simplify the life paring the band structures from DFT and from Wannier of future users, in their calculations of the magnetic in- function based Hamiltonian. We have tested using the teractions from first-principles ingredients, and enable symmetry-adapted Wannier functions [84] which, how- larger scale and more accurate micromagnetics calcula- ever, does not fully solve the problem. It appears that tions, exploring novel physical phenomena. the calculation of very small parameters using Wannier Future versions will include a wider range of inter- functions poses a challenge which hopefully future de- faces with DFT and ASD codes, as well as expanded velopments in the Wannierization methods can address feature sets, including single ion anisotropy, magnon by reducing the noise introduced in the disentanglement band structures for more complex systems. A promis- procedure. ing avenue is the calculation of higher order parame- At last, the Heisenberg model might not be valid for ters (3 spin, 4 spin), which has been explored only very all circumstances, for instance, when spins are itinerant sparsely in the literature (e.g. Refs. [85, 86, 87, 88]). 15 Acknowledgements [19] Liechtenstein, A. I., Katsnelson, M., Antropov, V., and Gubanov, V., Journal of Magnetism and Magnetic Materials 67 (1987) 65. The authors thank Yajun Zhang, Alireza Sasani, Jorge [20] Steenbock, T. and Herrmann, C., Journal of Computational Pilo Gonzlez, and Zachary Romestan for the testing of Chemistry 39 (2018) 81. the code, Thomas Ostler, Bertrand Dupé and Phivos [21] Katsnelson, M. I. and Lichtenstein, A. I., Phys. Rev. B 61 (2000) 8906. Mavropoulos for explanations about the limits and in- [22] Katsnelson, M. I., Kvashnin, Y. O., Mazurenko, V. V., and Licht- tricacies of fitting the Heisenberg model. enstein, A. I., Phys. Rev. B 82 (2010) 100403. This work has been funded by the Commu- [23] Solovyev, I., Hamada, N., and Terakura, K., Phys. Rev. Lett. 76 nauté Française de Belgique (ARC AIMED G.A. (1996) 4825. [24] Antropov, V., Katsnelson, M., and Liechtenstein, A., Physica B: 15/19-09).XH thanks the support by the EU H2020- Condensed Matter 237 (1997) 336. NMBP-TO-IND-2018 project ”INTERSECT” (Grant [25] Mazurenko, V. V. and Anisimov, V. I., Phys. Rev. B 71 (2005) No. 814487). EB thanks the FRS-FNRS for sup- 184434. port, as does MJV for an “out” sabbatical grant to [26] Mankovsky, S. and Ebert, H., Phys. Rev. B 96 (2017) 104416. [27] Secchi, A., Lichtenstein, A. I., and Katsnelson, M. I., Annals of ICN2 Barcelona in 2018-2019. The authors acknowl- Physics 360 (2015) 61. edge the CECI supercomputer facilities funded by the [28] Secchi, A., Lichtenstein, A. I., and Katsnelson, M. I., Journal of F.R.S-FNRS (Grant No. 2.5020.1) and the Tier-1 super- Magnetism and Magnetic Materials 400 (2016) 112. computer of the Fédération Wallonie-Bruxelles funded [29] Lounis, S. and Dederichs, P. H., Phys. Rev. B 82 (2010) 180404. [30] Szilva, A. et al., Phys. Rev. Lett. 111 (2013) 127204. by the Walloon Region (Grant No. 1117545). Com- [31] Yavorsky, B. Y. and Mertig, I., Phys. Rev. B 74 (2006) 174402. puting time was also provided by PRACE-3IP DECI [32] Ebert, H. and Mankovsky, S., Phys. Rev. B 79 (2009) 045209. grants 2DSpin and Pylight on Beskow (G.A. 653838 of [33] Andersen, O. K. and Jepsen, O., Phys. Rev. Lett. 53 (1984) 2571. H2020). [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 . References [36] Pizzi, G. et al., Journal of Physics: Condensed Matter 32 (2020) 165902. [1] Hirohata, A. et al., Journal of Magnetism and Magnetic Materi- [37] Gonze, X. et al., Comput. Phys. Commun. 180 (2009) 2582 als 509 (2020) 166711. , 40 YEARS OF CPC: A celebratory issue focused on quality [2] Baltz, V. et al., Rev. Mod. Phys. 90 (2018) 015005. software for high performance, grid and novel computing archi- [3] Bhatti, S. et al., Materials Today 20 (2017) 530 . tectures. [4] Joshi, V. K., Engineering Science and Technology, an Interna- [38] Gonze, X. et al., Comput. Phys. Commun. 248 (2020) 107042. tional Journal 19 (2016) 1503 . [39] Giannozzi, P. et al., Journal of Physics: Condensed Matter 29 [5] Lu, J. W., Chen, E., Kabir, M., Stan, M. R., and Wolf, S. A., (2017) 465901. International Materials Reviews 61 (2016) 456. [40] Soler, J. M. et al., Journal of Physics: Condensed Matter 14 [6] Žutić, I., Fabian, J., and Das Sarma, S., Rev. Mod. Phys. 76 (2002) 2745. (2004) 323. [41] Kresse, G. and Furthmüller, J., Phys. Rev. B 54 (1996) 11169. [7] Fazekas, P., Lecture notes on electron correlation and mag- [42] Kresse, G. and Joubert, D., Phys. Rev. B 59 (1999) 1758. netism, Series in modern condensed matter physics 5, World [43] Blaha, P., Schwarz, K., Madsen, G. K., Kvasnicka, D., and Scientific, 1999. Luitz, J., An augmented plane wave+ local orbitals program [8] Chumak, A. V., Vasyuchka, V. I., Serga, A. A., and Hillebrands, for calculating crystal properties (2001). B., Nat. Phys. 11 (2015) 453. [44] Blügel, S. and Bihlmayer, G., Forschungszentrum Jülich GmbH [9] Kiselev, N. S., Bogdanov, A. N., Schfer, R., and Rßler, U. K., (2006) 85. Journal of Physics D: Applied Physics 44 (2011) 392001. [45] Andrade, X. et al., Physical Chemistry Chemical Physics 17 [10] Fert, A., Cros, V., and Sampaio, J., Nat. Nanotechnol. 8 (2013) (2015) 31371. 152. [46] Tancogne-Dejean, N. et al., The Journal of Chemical Physics [11] Xiang, H., Lee, C., Koo, H.-J., Gong, X., and Whangbo, M.-H., 152 (2020) 124119. Dalton Trans. 42 (2013) 823. [47] Ozaki, T. and Kino, H., Physical Review B 72 (2005) 045121. [12] Herring, C., Magnetism: a treatise on modern theory and mate- [48] Enkovaara, J. et al., Journal of physics: Condensed matter 22 rials. 4. Exchange interactions among itinerant electrons, Aca- (2010) 253202. demic Press, 1966. [49] The Elk Code, http://elk.sourceforge.net/. [13] Sandratskii, L. M., physica status solidi (b) 136 (1986) 167. [50] Oroszlány, L., Ferrer, J., Deák, A., Udvardi, L., and Szunyogh, [14] Savrasov, S. Y., Phys. Rev. Lett. 81 (1998) 2570. L., Phys. Rev. B 99 (2019) 224412. [15] Romero, A. H. et al., The Journal of Chemical Physics 152 [51] Yoon, H., Kim, T. J., Sim, J.-H., and Han, M. J., Comput. Phys. (2020) 124102. Commun. 247 (2020) 106927. [16] Cao, K., Lambert, H., Radaelli, P. G., and Giustino, F., Phys. [52] Papior, N. R., B., J. L., Frederiksen, T., and Wuhl, S. S., zeroth- Rev. B 97 (2018) 024420. i/sisl: v0.9.8, 2020. [17] Ricci, F., Prokhorenko, S., Torrent, M., Verstraete, M. J., and [53] Solovyev, I., Dederichs, P., and Mertig, I., Physical Review B Bousquet, E., Phys. Rev. B 99 (2019) 184404. 52 (1995) 13419. [18] Phillips, J. J. and Peralta, J. E., The Journal of Chemical Physics [54] Weingart, C., Spaldin, N., and Bousquet, E., Phys. Rev. B 86 138 (2013) 174115. (2012) 094413. 16 [55] Udvardi, L., Szunyogh, L., Palotás, 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ý, 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., Com- putational Materials Science 42 (2008) 337 . [67] Jollet, F., Torrent, M., and Holzwarth, N., Comput. Phys. Com- mun. 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 Crys- tallography and Crystal Chemistry 28 (1972) 3384. [72] Xu, C., Xu, B., Dupé, 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ügel, S., Phys. Rev. B 101 (2020) 024418. 17