Comparative study of DFT+$U$ functionals for non-collinear magnetism
Siheon Ryee Myung Joon Han







Department of Physics, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Korea

mj.han@kaist.ac.kr
Department of Physics, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 34141, Korea





We performed comparative analysis for DFT+U functionals to better understand their applicability to non-collinear magnetism. Taking LiNiPO_4 and Sr_2IrO_4 as examples, we investigated the results out of two formalisms based on charge-only density and spin density functional plus U calculations. Our results show that the ground state spin order in terms of tilting angle is strongly dependent on Hund J. In particular, the opposite behavior of canting angles as a function of J is found for LiNiPO_4. The dependence on the other physical parameters such as Hubbard U and Slater parameterization F^4/F^2 is investigated. We also discuss the formal aspects of these functional dependences as well as parameter dependences. The current study provides 
useful information and important intuition for the first-principles calculation of non-collinear magnetic materials.




75.70.Cn, 73.20.-r, 75.47.Lx, 71.15.Mb

Comparative study of DFT+U functionals for non-collinear magnetism
    Myung Joon Han
    December 30, 2023
==================================================================








§ INTRODUCTION
  

DFT+U (density functional theory + U) method  <cit.> is a powerful tool to describe the correlated electron materials. Due to its computation efficiency and the straightforwardness of physical interpretation, DFT+U has become one of the most widely-used schemes to study complex magnetic phases in solids. On the other hand, an obvious drawback of DFT+U as a first-principles-based method is that the result is dependent on the interaction parameters such as Hubbard U and Hund J. The choice of double counting energy and density-functional type poses another important issue for DFT+U and related techniques such as DFT+DMFT (DFT + dynamical mean-field theory). 


A series of recent studies <cit.> have made notable progress in this regard. Careful investigations of the J dependence and the exchange-correlation (XC) functional dependence of widely-used DFT+U formulations show that adopting spin-density XC functional, such as LSDA (local spin-density approximation) and SGGA (spin-polarized generalized gradient approximation), can cause undesirable results. For clarity, we hereafter denote this type of formulation by `SDFT+U' (spin-density functional theory + U). The unphysical feature within SDFT+U originates from the intrinsic exchange energy within spin-density XC functionals. In order to avoid such an artifact, the use of charge density (i.e., spin-unpolarized) XC functional (denoted by `CDFT+U' to include both charge-only LDA+U and GGA+U simultaneously) is desirable <cit.>. Although many detailed features have been formally analyzed recently <cit.>, one important question still remains: Since all of these studies only dealt with collinear spin density, the case for non-collinear spin orders is basically unexplored to the best of our knowledge. While a previous study <cit.> provides useful information regarding the non-collinear magnetism within SDFT+U, more analysis at a deeper level and its comparison with CDFT+U results is still elusive. Not only because the non-collinear spin order is quite often realized in real materials <cit.>, but it also affects the other physical property such as electrical polarization in magneto-electric materials <cit.>, a systematic investigation is strongly requested.


In this paper, we performed a systematic comparative study of DFT+U formalisms for the case of non-collinear magnetic orders. With LiNiPO_4 (LNPO) and Sr_2IrO_4 (SIO) as our examples, we investigated the magnetic ground state property predicted by SDFT+U and CDFT+U, and its dependence on physical as well as numerical parameters. We show that the two formalisms give the opposite trend of spin canting angle as a function of J while the U dependence is weak in LNPO. The difference is mainly attributed to the different treatment of spin-off-diagonal potential terms. 
For the case of SIO, the canting angle exhibits a stronger dependence on the choice of U while it exhibits the same increasing trend as a function of J in both formalisms. The dependence on the Slater parameterization and the local projector is also investigated.
Our results clearly show that the current DFT+U schemes have limitations in describing non-collinear magnetism, and further methodological development is requested. As a step toward this direction, we present and discuss our analysis of the two functionals in non-collinear spin space.




 





§ COMPUTATION METHOD


All calculations were performed using rotationally invariant version of DFT+U method <cit.> as implemented in `OpenMX' software package <cit.>. Our formalism is based on the non-orthogonal LCPAO (linear combination of localized pseudoatomic orbitals) <cit.>. The effect of SOC (spin-orbit coupling) was treated within the fully relativistic pseudopotentials <cit.>. Experimental crystal structure was used for LNPO (Pnma) and SIO (I4_1/acd) with lattice parameters of a=10.02, b=5.83, c=4.66 Å <cit.>,  and a=b=5.48, c=25.80 Å <cit.>, respectively (see Fig. <ref>). We adopted 3 × 5 × 5 (7 × 7 × 3) 𝐤-points in the first Brillouin zone of LNPO (SIO) and the energy cutoff of 500 Ry in real space grid for numerical integrations. An energy criterion of 10^-9 Hartree was used for the self-consistency. The radial cutoff of PAOs were set to 8.0, 6.0, 7.0, 10.0, 7.0 and 5.0 a.u. for Li, Ni, P, Sr, Ir, and O, respectively. For XC functional, L(S)DA <cit.> parameterized by Perdew and Zunger <cit.> was used. Unless otherwise specified, the `dual' projector <cit.> was used for the on-site density matrix (DM), and the U value was set to 5 eV for LNPO and 2 eV for SIO following the previous DFT+U calculations <cit.>. 

The total energy correction by DFT+U, E^U, can be expressed as:

    E^U = E^int - E^dc,

where E^int represents the Hubbard-type on-site interaction energy <cit.>.
For the double counting term (E^dc), we adopted the so-called FLL (fully localized limit) functional. For the case of SDFT+U it can be written as <cit.>

    E_SDFT+U^dc,FLL = 1/2UN(N-1) - 1/2JN(N/2-1) -1/4JM⃗·M⃗.

It is worth to be noted that the widely-used Dudarev's approach <cit.> is equivalent to the J=0 version of SDFT+U employing FLL double counting <cit.>.
For CDFT+U, it reads <cit.> 

    E_CDFT+U^dc,FLL = 1/2UN(N-1) - 1/2JN(N/2-1),

where the d-electron occupation is N=Tr[𝐧], and n^σσ'_m_1m_2 are the elements of on-site DM, 𝐧, for given orbitals {m_i} and spins σ,σ' (σ,σ' = ↑ or ↓). Note that the correlated orbtials are responsible for the magnetization in CDFT+U, and the rest space is affected by hybridization with them. The spin moment is given by M⃗ = Tr[σ⃗𝐧] where σ⃗ is Pauli matrices. For clarity we use `sFLL' and `cFLL' in referring to `SDFT+U' and `CDFT+U' (expressed with the FLL double counting functional), respectively, throughout the rest of this paper. For more formulation details, see Ref. <cit.> and references therein.







§ RESULT AND DISCUSSION



 §.§ LiNiPO_4



The Li orthophosphates LiMPO_4 (M = Mn, Fe, Co, Ni) are magneto-electric materials <cit.>. A strong magneto-electric effect is observed in antiferromagnetic (AF) phase of LiMPO_4. In the present study, we focus on the non-collinear magnetic ground state of LNPO. Neutron diffraction shows that `C_zA_x'-type (C-type AF order along c-axis and A-type along a-axis with small cantings) spin alignment is stabilized at low temperature <cit.> (see Fig. <ref>(b)).







To investigate the predictive capability of current DFT+U methods for the non-collinear magnetic ground state, we performed the constrained spin-moment total energy calculations <cit.> in which the spins are restricted to lie within the ac-plane (as known from experiment) with an angle θ (as shown in Fig. <ref>(b)). Note that no constraint is imposed on orbital moment. The results are summarized in Fig. <ref>  where the zero energy points correspond to the uncanted (i.e., collinear) spin order. Since the stable spin configuration is sensitive to J as reported by Bousquet and Spaldin <cit.>, five different J values in the range of 0≤ J≤ 1 eV have been considered.


First of all, we note the different J dependence found in between sFLL (Fig. <ref>(a)) and cFLL (Fig. <ref>(b)). In sFLL, the larger canting angle is favored for the larger J with the greater stabilization energy: At J=0 eV, the canting angle is θ = 1.88^∘ and it gradually increases to θ =6.76^∘ at J=1 eV. In comparison to the experimental value of θ≃ 7.8^∘ <cit.>, the best agreement is expected at J>1 eV which might be too large to be realistic. In the entire range of J we considered, the moment size change is negligible; μ_s≃ 1.8 μ_B/Ni. 







In cFLL, on the other hand, the canting angle is reduced as J increases, and it eventually goes to θ=0 at J ≃ 0.7 eV. For J > 0.7 eV, θ becomes negative (e.g., θ≃ -2.9^∘ at J=1 eV), which means the canting angles are opposite along a axis. Contrary to sFLL, the size of spin moment is gradually increased as J increases as discussed in Ref. <cit.>; from μ_s =1.57 μ_B/Ni at J=0 eV to μ_s =1.73 μ_B/Ni at J=1 eV. The experimental value of θ≃ 7.8^∘ is achieved at J≃ 0 eV which is also unrealistic.


The contribution from orbital moment is notably smaller than the spin (see black lines in Fig. <ref>(a)); μ_o∼ 0.26 μ_B/Ni in the entire region of J within sFLL while it increases from 0.50 μ_B/Ni (J=0 eV) to 0.65 μ_B/Ni  (J=1 eV) within cFLL.


The U dependence of canting angle is less pronounced as reported in Ref. <cit.>. Our results of U=4 and 5 eV are presented in Fig. <ref>(b). The size of spin moment for U=4 eV is only slightly reduced from the U=5 eV result by ∼ 0.01 μ_B/Ni.


We also checked the effect of using different Slater parameterization in generating Coulomb interaction tensor. This parameterization has not attracted much attention in DFT+U methodology. As studied in Ref. <cit.>, however, due to the screening in solids, F^4/F^2 can be enhanced and thus be different from the conventional value of F^4/F^2=0.625 which corresponds to the atomic environment. To further elucidate this point within the context of non-collinear formalism, we performed the calculations with F^4/F^2=0.9 and 
the result is compared with that of F^4/F^2=0.625 in Fig. <ref>(c). 
It is found that the canting-angle change caused by different parameter choice is not significant; At J=1 eV, Δθ∼ +0.8^∘ for cFLL and ∼ +0.3^∘ for sFLL. The moment size change is also quite small, ∼ 0.01μ_B, in both cFLL and sFLL.



As a summary of current subsection, we note that, with the current DFT+U formalisms, it is difficult to determine or predict the non-collinear spin angle. The J dependence is noteworthy.
Considering the cRPA (constrained random phase approximation) estimation of J=0.7 – 0.9 eV for (isovalent) NiO <cit.>, sFLL and cFLL gives θ≃ +4^∘ – +6^∘ and ≃ 0^∘ – -2^∘, respectively. Further, the origin of opposite behavior of canting angles is far from clear. This issue together with the physical meaning of atomic Hund interaction certainly deserves further investigations. 







 §.§ Sr_2IrO_4


The second example is SIO which is known as a `relativistic Mott' insulator <cit.> exhibiting a canted AF order with a canting angle ϕ∼ 12-13^∘ within a IrO_2 plane (see Fig. <ref>(c)) <cit.>. 
Its low energy behavior can be described by so-called J_eff=1/2 Kramers doublet; namely, |↑̃⟩ = sinθ_t  |0,↑⟩ - cosθ_t  |+1,↓⟩ and |↓̃⟩ = sinθ_t  |0,↓⟩ - cosθ_t  |-1,↑⟩ where θ_t is parameterized via tan(2θ_t)=2√(2)λ/(λ-2Δ) <cit.>. Here, λ and Δ refers to the SOC strength and the tetragonal splitting within t_2g manifold, respectively, and |0⟩=|d_xy⟩ and |±1⟩=-1/√(2)(i|d_xz⟩± |d_yz⟩). 
The cooperation of isotropic superexchange and Dzyaloshinskii-Moriya interaction leads to the canted AF state carrying net ferromagnetic moment of M ∼ 0.1 - 0.2 μ_B/Ir <cit.> along a-direction (Fig. <ref>(c)). 







The calculation results of canting angle are summarized in Fig. <ref>(a). 
Here we present the angle ϕ of total moment μ_tot (the sum of spin and orbital moment at each Ir site) as a function of J (see Fig. <ref>(c)). On the contrary to the case of LNPO, sFLL and cFLL exhibit the same increasing trend of canting angle as J increases while cFLL gives a better agreement with the experimental value of ϕ∼ 12-13^∘ <cit.> at realistic J (see Fig. <ref>(a)). Note that the canting angle is sensitive to the structural distortion which is fixed to the experimental value in our calculations <cit.>. 


The choice of U and J can significantly affect the net moment size as shown in Fig. <ref>(b). The cFLL result at J=0.6 eV is in a good agreement with experimental value of M ∼ 0.1 - 0.2 μ_B/Ir <cit.>. The net moment enhancement as a function of J is originated from both the enhanced total moment μ_tot and the enlarged canting angle. See the inset of Fig. <ref>(b) which shows that, for U=2 eV, the total moment μ_tot is increased from μ_tot=0.35 (J=0 eV) to 0.5 μ_B/Ir (J=0.6 eV) within cFLL while μ_tot≃ 0.47 μ_B/Ir within sFLL in the entire range of J . 
 
 
The U dependence of canting angle is also presented in Fig. <ref>(a). The larger value of U by 0.5 eV results in the reduced canting angle by ∼ 1^∘  in the entire J range we considered for both sFLL and cFLL. The U dependence of moment size can be found in Fig. <ref>(b). For μ_tot (not shown), it is increased by ∼ 0.12μ_B/Ir (∼ 0.07 μ_B/Ir) within cFLL (sFLL) when we use U=2.5 eV instead of 2.0 eV. For the case of sFLL, this is a comparable size of change caused by varying J.




The results of F^4/F^2 dependence are presented in Fig. <ref>(c). In both formalisms of sFLL and cFLL, it is found that the canting angle is reduced by ∼ 0.5^∘ when using F^4/F^2=0.9 at J=0.6 eV. It is worth to be noted that increasing F^4/F^2 ratio leads to the change of off-diagonal elements in the Coulomb interaction matrix by about ± 0.01J (see Sec. <ref>). Considering that this is a small change in terms of interaction strengths, the canting angle is sensitive to the parameterization of Slater integrals especially in sFLL for which J dependence is much weaker.






 §.§ Further analysis and discussion
 

Our results of two material examples clearly show that standard DFT+U formalisms can give quite different solutions. It also implies that there is no a priori rule in choosing a functional for the description of non-collinear magnetism. For LNPO, sFLL seems to give a better agreement with experiment whereas cFLL is the case for SIO. The effect of U on the canting angle is material dependent in both formalisms; it is basically negligible in LNPO but not in SIO.

Thus our result raises a question regarding the predictive power of current DFT+U methodology for non-collinear magnetism. Since there is no good intuitive reason available for these behaviors, one needs to be careful when there is no guide from outside, e.g., from experiments. Further, as can be noticed in the cases of LaMnO_3 <cit.>, Li orthophosphates, and difluorites <cit.>, reproducing the correct electronic and spin structures does not always guarantee the physically relevant U and J range. In particular, within SDFT+U scheme, unphysically large J values can be required <cit.>. It can also be true in the non-collinear case (as in the collinear case) that some redundant contributions reside in LSDA/SGGA on top of atomic on-site interactions <cit.>. However, the usual sense of `interactions' does not seem directly responsible for non-collinearity, and therefore our intuition of choosing those parameters is not quite helpful any longer in this situation. 


At this moment, it is quite challenging to fully understand the physical meaning of the observed differences in between sFLL and cFLL, and to construct the revised DFT+U functional compatible with non-collinear spin density. As a step towards this direction, we further analyze below the origin of this behavior within the currently available standard formulation. In the non-collinear spin DFT, the XC potential is given by <cit.>

    𝐕_XC   = 
    (
    [ V^↑↑_XC V^↑↓_XC; V^↓↑_XC V^↓↓_XC ]) 
       =
    (
    [  V^0_XC+Δ V_XCcosβ Δ V_XCexp(-iα)sinβ;  Δ V_XCexp(iα)sinβ  V^0_XC-Δ V_XCcosβ ]),

where α and β are the azimuthal and polar Euler angles, respectively. V^0_XC = (V^↑_XC+V^↓_XC)/2 and Δ V_XC = (V^↑_XC-V^↓_XC)/2 is the charge potential and the spin up/down potential difference, respectively, obtained from the diagonalization of non-collinear spin-density matrix:

    𝐧= 
    (
    [ n^↑↑ n^↑↓; n^↓↑ n^↓↓ ]).

It is obvious from Eq. (<ref>) that α=β=0 is reduced back to the usual collinear spin-density functional.





Note that, in cFLL (i.e., with spin-unpolarized charge density in XC functional), Δ V_XC=0. As a result, the off-diagonal double counting potential should vanish for cFLL; V_m_1m_2^dc,σσ = ∂E^dc,FLL_CDFT+U/∂ n^σσ_m_1m_2=0 (where σ denotes the opposite spin to σ). For sFLL, on the other hand,  V^dc,σσ_m_1m_2=∂E^dc,FLL_SDFT+U/∂ n^σσ_m_1m_2 = -JN^σσδ_m_1 m_2 (where N^σσ=Tr[n^σσ]), which is not canceled out in general by Δ V_XC term in V^σσ_XC. As reported in Ref. <cit.>, enforcing V^dc,σσ_m_1m_2=0 within sFLL affects the spin canting. However, this treatment requires a careful physical interpretation because the contribution from spin-density XC functional should have been subtracted to avoid double counting. On the contrary, in cFLL, the off-diagonal spin potential and energy contain only E^int-related terms due to Δ V_XC=0 (see Eq. (<ref>)).


Since the off-diagonal potential is affected by Δ V_XC, it is useful to compare LSDA result with SGGA at J=0. SGGA has a stronger tendency toward the magnetic solution than LSDA, and the larger XC potential is expected in SGGA than LSDA <cit.>; |Δ V^SGGA_XC| > |Δ V^LSDA_XC|. It again affects the non-collinear spin potential through Eq. (<ref>). 
We found that the most stable canting angle of LNPO is θ =5.41^∘ and θ = 2.14^∘ within LSDA and SGGA (PBE parameterization <cit.>), respectively, at U=0 eV. For the case of SIO,  the difference is small; ϕ =7.81^∘ and 7.94^∘ in LSDA and SGGA, respectively. Note that the enhanced |Δ V_XC| gives the less (more) canted solution for LNPO (SIO). Thus, this feature is consistent with the result of cFLL (Fig. <ref>(b) and Fig. <ref>(b))) in which the spin-off-diagonal potential <cit.>,

    V_m_1m_2^U,σσ   = - ∑_m_3,m_4⟨ m_1,m_3|V_ee|m_4,m_2 ⟩ n^σσ_m_3m_4,

has the terms that are proportional to the elements of Coulomb interaction tensor, 

    ⟨ m_1,m_3|V_ee|m_4,m_2 ⟩
       = ∑_{m_i'}[S_m_1m_1'S_m_3m_3'{∑_k=0α_k(m_1',m_3',m_2',m_4')F^k}
    S^-1_m_2'm_4S^-1_m_4'm_2].

Note that the magnitude of matrix elements increases with J (see Eq. (<ref>)). Here F^k and α_k are Slater integrals and Racah-Wigner numbers, respectively, and S is a transformation matrix from spherical harmonics to the pre-defined local basis <cit.> (see, e.g., Ref. <cit.> for more discussion on Coulomb interaction in transition metal oxides). 


Another important issue is the on-site electron density. How to define the correlated subspace has been an important question in DFT+U formalism <cit.>. As Eq. (<ref>) clearly shows, it can also affect the non-collinear spin ground state through the off-diagonal spin potential. Since the choice of on-site density has always some degree of ambiguity, it is important to carefully check this dependence. In the previous subsections, we compared two different local projectors, previously named as `full' and `dual' projector <cit.>. The former takes the smaller d occupation compared to the latter by construction <cit.>. The calculation result of LNPO shows that the J dependence of the canting angle is weaker when the `full' projector is used in both sFLL and cFLL. It is found that within sFLL θ = 1.88^∘ (`dual') and 2.20^∘ (`full') at J=0 eV which is in reasonable agreement with θ = 1.6^∘ obtained by using PAW (projector-augmented-wave) formalism <cit.>. In cFLL, the difference between the two projectors is more pronounced as the off-diagonal potential is solely originated from interaction terms. In particular, for large J > 0.7 eV, the `dual' projector gives the negative θ whereas the `full' still favors the positive θ. It is noted that such technical detail of choosing local projector can lead to a qualitatively different magnetic solution. 


The U dependence of canting angle is weaker than the case of J for LNPO as reported in Ref. <cit.> while it was not negligible in SIO. We could not find a meaningful change in DM within our J range. Instead, in order to understand this dependence, we once again investigate the explicit form of ⟨ m_1,m_3|V_ee|m_4,m_2 ⟩ in Eq. (<ref>), which is given in general by

    ⟨ m_1,m_2|V_ee|m_2,m_1 ⟩ = 
    (
    [ U+1.14J   0.89J   0.89J   0.54J   0.54J;   0.89J U+1.14J   0.43J   0.77J   0.77J;   0.89J   0.43J U+1.14J   0.77J   0.77J;   0.54J   0.77J   0.77J U+1.14J   0.77J;   0.54J   0.77J   0.77J   0.77J U+1.14J ]),

where the diagonalized DM is assumed with real harmonics basis  (|d_z^2⟩, |d_x^2-y^2⟩, |d_xy⟩, |d_xz⟩, and |d_yz⟩). The conventional Slater parameterizations of F^0=U, (F^2+F^4)/14=J, and F^4/F^2=0.625 were used to generate Eq. (<ref>) <cit.>. One can notice that U is included only in the diagonal part of Eq. (<ref>). As for the case of LNPO, the weak U dependence of spin canting angle implies the dominant contribution of the inter-orbital exchange interactions (i.e., the off-diagonal elements of Eq. (<ref>)) in determining the non-collinear ground state. Note that changing F^4/F^2 does not affect the diagonal elements of interaction matrix (Eq. (<ref>)). Thus the canting angle change caused by different F^4/F^2 ratio shown in Fig. <ref>(c) is also indicative of the role of inter-orbital exchange interactions. On the other hand, in the case of SIO where both U and J dependences are noticeable, the situation becomes more complicated, and not well understood within this simple formal analysis. As discussed within the collinear SDFT+U case <cit.>, the J dependence of Coulomb interaction matrix (Eq. (<ref>)) gives rise to the nontrivial effect on the electronic structure. Our results of LNPO and SIO demonstrate that this is also the case for non-collinear magnetism.







§ CONCLUDING REMARKS


We performed a comparative study of DFT+U formalisms for the case of non-collinear magnetic orders by examining the results of LNPO and SIO. The two formulations of cFLL and sFLL are found to give different J dependences for the ground state magnetic order, and the difference is attributed to the spin-density XC functional which largely affects the canting angle especially through spin-off-diagonal potentials. Our result demonstrates that the inter-orbital exchange energy is a crucial component for non-collinear simulation and the care needs to be paid especially when SDFT+U functional is adopted. Both sFLL and cFLL fail to reproduce the experimental canting angle in the reasonable range of interaction parameters for LNPO while cFLL gives a fairly good description for SIO. Our result poses a challenging issue for the capability and the accuracy of current DFT+U and the related methods to describe the non-collinear magnetism in real materials. One promising aspect is that both formalisms give the qualitatively consistent results of C_zA_x spin order for LNPO with experiments. Since our investigation is based on the fixed experimental lattice parameters, the discrepancy found in this study might become less significant when the lattice parameters are optimized. Further investigations at a more fundamental level are requested for the better first-principles description of non-collinear magnetic materials.







§ ACKNOWLEDGEMENTS

This work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2018R1A2B2005204). 









§ REFERENCES





