On site Coulomb interaction: L(S)DA+U

Supported only by VASP.4.6 and on. THIS FEATURE IS IN LATE BETA STAGE (BUGS ARE POSSIBLE).

The L(S)DA often fails to describe systems with localized (strongly correlated) $d$ and $f$ electrons (this manifests itself primarily in the form of unrealistic one-electron energies). In some cases this can be remedied by introducing a strong intra-atomic interaction in a (screened) Hartree-Fock like manner, as an on site replacement of the L(S)DA. This approach is commonly known as the L(S)DA+U method.

VASP allows one to choose between two different approaches to L(S)DA+U:

• The rotationally invariant version introduced by Liechtenstein et al. [64], which is of the form
  
  $$E_{\rm HF}=\frac{1}{2} \sum_{\{\gamma\}} (U_{\gamma_1\gamma_... ...mma_2}){ \hat n}_{\gamma_1\gamma_2}{\hat n}_{\gamma_3\gamma_4}$$
  
  
  and is determined by the PAW on site occupancies
  
  $${\hat n}_{\gamma_1\gamma_2} = \langle \Psi^{s_2} \mid m_2 \rangle \langle m_1 \mid \Psi^{s_1} \rangle$$
  
  
  and the (unscreened) on site electron-electron interaction
  
  $$U_{\gamma_1\gamma_3\gamma_2\gamma_4}= \langle m_1 m_3 \mid \... ...e\vert} \mid m_2 m_4 \rangle \delta_{s_1 s_2} \delta_{s_3 s_4}$$
  
  
  ($\vert m \rangle$ are the spherical harmonics)

  The unscreened e-e interaction $U_{\gamma_1\gamma_3\gamma_2\gamma_4}$ can be written in terms of Slater's integrals $F^0$, $F^2$, $F^4$, and $F^6$ (f-electrons). Using values for the Slater integrals calculated from atomic wave functions, however, would lead to a large overestimation of the true e-e interaction, since in solids the Coulomb interaction is screened (especially $F^0$).

  In practice these integrals are therefore often treated as parameters, i.e., adjusted to reach agreement with experiment in some sense: equilibrium volume, magnetic moment, band gap, structure. They are normally specified in terms of the effective on site Coulomb- and exchange parameters, $U$ and $J$. ($U$ and $J$ are sometimes extracted from constrained-LSDA calculations.)

  These translate into values for the Slater integrals in the following way (as implemented in VASP at the moment):

  • $p$-electrons: $F^0=U$, $F^2=5J$
  • $d$-electrons: $F^0=U$, $F^2=\frac{14}{1+0.625}J$, and $F^4=0.625 F^2$
  • $f$-electrons: $F^0=U$, $F^2=\frac{6435}{286+195 \cdot 0.668+250 \cdot 0.494}J$, $F^4=0.668 F^2$, and $F^6=0.494 F^2$

  The essence of the L(S)DA+U method consists of the assumption that one may now write the total energy as:
  

  $$E_{\rm tot}(n,\hat n)=E_{\rm DFT}(n)+E_{\rm HF}(\hat n)-E_{\rm dc}(\hat n)$$
  
  
  where the Hartree-Fock like interaction replaces the L(S)DA on site due to the fact that one subtracts a double counting energy ($E_{\rm dc}$) which supposedly equals the on site L(S)DA contribution to the total energy.

  Currently VASP allows for the choice between two different definitions for the double counting energy:

  LSDA+U	$E_{\rm dc}(\hat n) = \frac{U}{2} {\hat n}_{\rm tot}({\hat n}_{\rm tot}-1) - \frac{J}{2} \sum_\sigma {\hat n}^\sigma_{\rm tot}({\hat n}^\sigma_{\rm tot}-1)$
  LDA+U	$E_{\rm dc}(\hat n) = \frac{U}{2} {\hat n}_{\rm tot}({\hat n}_{\rm tot}-1) - \frac{J}{4} {\hat n}_{\rm tot}({\hat n}_{\rm tot}-2)$

• The simplified (rotationally invariant) approach to the LSDA+U , due to Dudarev et al. [65], is of the following form:
  
  $$E_{\rm LSDA+U}=E_{\rm LSDA}+\frac{(U-J)}{2}\sum_\sigma \left... ...n_{m_1,m_2}^{\sigma} \hat n_{m_2,m_1}^{\sigma} \right) \right]$$
  
  

  This can be understood as adding a penalty functional to the LSDA total energy expression that forces the on site occupancy matrix in the direction of idempotency, i.e., $\hat n^{\sigma} = \hat n^{\sigma} \hat n^{\sigma}$. (Real matrices are only idempotent when their eigenvalues are either 1 or 0, which for an occupancy matrix translates to either fully occupied or fully unoccupied levels.)

  Note: in Dudarev's approach the parameters $U$ and $J$ do not enter seperately, only the difference $(U-J)$ is meaningfull.

The L(S)DA+U in VASP is switched on by means of the following tags

• LDAU = .TRUE. Switches on the L(S)DA+U.
• LDAUTYPE = 1|2|4 Type of L(S)DA+U (Default: LDAUTYPE = 2)

  1

  Rotationally invariant LSDA+U according to Liechtenstein et al.

  4

  Idem 1., but LDA+U instead of LSDA+U (i.e. no LSDA exchange splitting)

  2

  Dudarev's approach to LSDA+U (Default)

• LDAUL = L .. $l$-quantum number for which the on site interaction is added
  (-1: no on site terms added, 1: p, 2: d, 3: f, Default: LDAUL = 2)
• LDAUU = U .. Effective on site Coulomb interaction parameter
• LDAUJ = J .. Effective on site Exchange interaction parameter
• LDAUPRINT = 0|1|2 Controls verbosity of the L(S)DA+U module
  (0: silent, 1: Write occupancy matrix to OUTCAR, 2: idem 1., plus potential matrix dumped to stdout, Default: LDAUPRINT = 0)

NB: LDAUL, LDAUU, and LDAUJ must be specified for all atomic species!

It is important to be aware of the fact that when using the L(S)DA+U, in general the total energy will depend on the parameters $U$ and $J$. It is therefore not meaningful to compare the total energies resulting from calculations with different $U$ and/or $J$ [c.q. $(U-J)$ in case of Dudarev's approach].

Note on bandstructure calculation: The CHGCAR file also contains only information up to LMAXMIX for the on-site PAW occupancy matrices. When the CHGCAR file is read and kept fixed in the course of the calculations ( ICHARG=11), the results will be necessarily not identical to a selfconsistent run. The deviations can be (or actually are) large for L(S)DA+U calculations. For the calculation of band structures within the L(S)DA+U approach, it is hence strictly required to increase LMAXMIX to 4 (d elements) and 6 (f elements). (see Sec. 6.56).