Tutorial 15: First principles calculation of \(U\) and \(J\)

Author:: Davide Sarpa
Date:: Aug 2024

Introduction

The goal of the tutorial is to provide a working example on how it is possible to compute the \(U\) and \(J\) parameters from first principles. We will work on Hematite +--+ antiferromagnetic configuration as you should already be familiar with it, if not, refer to tutorial 9.

The reason behind computing the parameters via first principles is because they directly correct the spurious localised self-interaction error (\(U\)) and static correlation error (\(J\)) and hence the physics of the system. While choosing and empirical \(U\) and \(J\) might give a better description of a specific property of the material, it does not guarantee that these errors are consistently corrected.

Fig. 10 Primitive rhombohedral cell (left), conventional hexagonal cell (right). Fe atoms with spin up and down are in green and pink, respectively. O atoms are in orange.

Theoretical background

We start by defining the response function \(\chi\), which describes how the occupation of localised orbitals changes with respect to a shift in the potential acting on these orbitals: The linear response method determines the Hubbard \(U\) parameter by comparing the response of the system to a perturbation in standard DFT and DFT+\(U\) frameworks.

We define the response function \(\chi\) as:

\[\chi = \frac{dn^{I\sigma}}{d\alpha}\]

where \(n\) is the occupation matrix of the localised orbitals and \(\alpha\) is a potential shift applied to these orbitals.

We compute two response functions:

\(\chi_0\): the bare Kohn-Sham (KS) response (without \(U\))
\(\chi\): the interacting response (with \(U\))

These are related by:

\[U = \chi^{-1} - \chi_0^{-1}\]

which allows us to compute \(U\).

The perturbation is applied by shifting the potential of the localised orbitals:

\[V_{\text{ext}}^{p} = V_{\text{ext}} + \alpha \sum_{m,m'}\lvert\varphi_{m'}^{(I)}\rangle\langle\varphi_m^{(I)}\rvert\]

This is the conventional linear response and its done in a supercell as the perturbation should not interact with its periodic images. Another approach to compute \(U\) and \(J\) is known as minimum tracking method [Moynihan2017] [Linscott2018].

Minimum Tracking Method

The minimum tracking method is based on a reformulation of the response matrices based on the ground state density of the perturbed system [Moore2024] . We can redefine the interacting and noninteracting response matrices as (in practice we’ll be using simpler yet equivalent formulae)

\[\chi_{IJ} = \frac{dn^I}{dv_\text{ext}^J},\]

\[(\chi_0)_{IJ} = \left[\frac{dn}{dv_\text{KS}}\left(\frac{dv_\text{KS}}{dv_\text{ext}}\right)^{-1}\right]_{IJ}\]

This allows us to work around the practical issues from the conventional linear response. This approach can also be extended to include the \(J\) exchange term In practice this is done by modifying the perturbation by including an additional term (spin-splitting):

\[V_{\text{ext}}^{p} = V_{\text{ext}} + \beta \sum_{m,m'}\lvert\varphi_{m'}^{(I\uparrow)}\rangle\langle\varphi_m^{(I\uparrow)}\rvert-\lvert\varphi_{m'}^{(I\downarrow)}\rangle\langle\varphi_m^{(I\downarrow)}\rvert\]

Setting up the calculations

We will configure a set (9 total) of bulk hematite single-point calculations to compute \(U\) and \(J\) for the Fe \(3d\) orbitals. We apply distinct labels to Fe atoms, enabling us to assign different parameters to spin-up and spin-down Fe atoms. We will be using a 4x4x1 supercell generated from the conventional cell.

Tutorial files

All the files needed for the simulations can be downloaded from

Practical calculation

The step by step approach to compute \(U\) and \(J\) is:

add hubbard_calculating_u : T in the input file,

choose an atom for the atom type we want to compute \(U\) or \(J\) for, and label it differently. In our case you can see from the input file that we have labelled this single atom Fe1U. It does not matter whether we choose a spin up or spin down atom for an AFM material.

apply the perturbation to this atom only and perform single-points calculations,

compute \(U\) and \(J\) with the following formulas:

\[U = \frac{1}{2} \frac{\delta v^\uparrow_{\text{Hxc+local}} + \delta v^\downarrow_{\text{Hxc+local}}}{\delta(n^\uparrow + n^\downarrow)}\]

\[J = -\frac{1}{2} \frac{\delta v^\uparrow_{\text{Hxc+local}} - \delta v^\downarrow_{\text{Hxc+local}}}{\delta(n^\uparrow - n^\downarrow)}\]

where \(\delta v^\uparrow_{\text{Hxc}}\) and \(\delta v^\downarrow_{\text{Hxc}}\) represent the derivative of the Hxc+local potential with respect to the applied potential (either \(\alpha\) to compute \(U\) or \(\beta\) to compute \(J\)) and \(\delta(n^\uparrow + n^\downarrow)\) and \(\delta(n^\uparrow - n^\downarrow)\) represent the derivative of the total occupation \(n^\uparrow + n^\downarrow\) with respect to \(\alpha\) and of \(n^\uparrow - n^\downarrow\) with respect to \(\beta\).

How and where to apply the perturbation

Looking at the input file provided you can see we activated the hubbard_calculating_u functionality and in the Hubbard block we have

%BLOCK HUBBARD
Fe1  2 0.0 0.0 -10.0 0.0 0.0
Fe1U 2 0.0 0.0 -10.0 0.0 0.0
Fe2  2 0.0 0.0 -10.0 0.0 0.0
%ENDBLOCK HUBBARD

where the columns of the hubbard block are described as follows:

Species Label

The species to apply the DFT+\(U\) correction to.
Angular Momentum: \(l\)

The angular momentum of the projectors which the Hubbard correction is applied to. In this example \(l=2\) which corresponds to d orbitals
Hubbard \(U\) value

The value of the Hubbard \(U\) for this sub-shell, in electron-volts. We are computing it so we can choose 0 as its value
Hund’s exchange \(J\) value

The value of the Hund’s exchange \(J\) for this sub-shell, in electron-volts. We are computing it so we can choose 0 as its value
Effective Charge \(\mathbf{Z}\) and Projectors type The default projectors are NGWFs. For other possibility, refer to the DFT+ \(U\) documentation
The \(\alpha\) prefactor

The perturbation term needed to compute \(U\)
The spin-splitting factor \(\beta\)

The perturbation term needed to compute \(J\).

To compute \(U\) you need to change the \(\alpha\) value while keeping \(\beta\) equal to 0. To compute \(J\) you need to change the \(\beta\) value while keeping \(\alpha\) equal to 0.

We have provided you only 1 input file – the one corresponding to 0 for both \(\alpha\) and \(\beta\), you need to generate the remaining 8 files.

The \(\alpha\) and \(\beta\) values you need to use for the \(U\) calculation are = -0.2, -0.1, 0.0, 0.1, 0.2.

Why these values? We want to apply a big enough perturbation to see an effect and to be able to compute derivatives but also remain in the linear regime. It is not necessary to use 5 datapoints to obtain a good value but it’s highly recommended.

Evaluating the outputs

In order to compute \(U\) and \(J\) we need the values of the \(v^\uparrow_{\text{Hxc}}\) and \(v^\downarrow_{\text{Hxc}}\),which can be found in the following block:

################################################################################
DFT+U information on Hubbard site     72 of species Fe1U and spin      1
The average Hxc+local potential is           -100.04043423 eV.
The average Hubbard potential is               -0.10000000 eV.
################################################################################
DFT+U information on Hubbard site     72 of species Fe1U and spin      2
The average Hxc+local potential is            -96.03296381 eV.
The average Hubbard potential is               -0.10000000 eV.
################################################################################

Note that we are looking only at the values for Fe1U atom which is the only atom we have applied the perturbation to. There are multiple instances of this block and we are only interested in the last one.

Next, we need to look at occupation of the Hubbard manifold \(n^\uparrow + n^\downarrow\) and \(n^\uparrow - n^\downarrow\),which can be found in the following block:

################################################################################
DFT+U information on atom      1 of Hubbard species Fe1U
################################################################################
Occupancy matrix of Hubbard site     72 and spin      1 is
m_l =   -2          -1           0           1           2
0.98583311  0.01105739  0.00017283  0.00149346 -0.00039754
0.01106973  0.98239066 -0.00021203  0.00037893  0.00244851
0.00017266 -0.00021405  0.99296562  0.00030517  0.00069962
0.00149451  0.00037878  0.00029134  0.98210951 -0.01203475
-0.00039830  0.00244943  0.00069122 -0.01204334  0.98340592
WARNING: DFT+U ENERGY of Hubbard site     72 and spin      1 is negative.
################################################################################
Occupancy matrix of Hubbard site     72 and spin      2 is
m_l =   -2          -1           0           1           2
0.32009924 -0.06393836 -0.00012245 -0.01033413 -0.00070413
-0.06400973  0.33409081 -0.00029354  0.00034179 -0.01142806
-0.00012106 -0.00027777  0.19025018 -0.00114325  0.00745246
-0.01034138  0.00034159 -0.00106271  0.33014982  0.06774687
-0.00070499 -0.01143070  0.00762074  0.06779446  0.29199808
WARNING: DFT+U ENERGY of Hubbard site     72 and spin      2 is negative.
################################################################################
Total occupancy of Hubbard site     72 is         6.39329292 e
Local magnetic moment of Hubbard site     72 is   3.46011669 mu_B
DFT+U energy of Hubbard site     72 is           -0.02349492 Ha
################################################################################

The total occupancy of Hubbard site is the \(n^\uparrow + n^\downarrow\), while the local magnetic moment of Hubbard site is the \(n^\uparrow - n^\downarrow\). We now have all the data we need to compute \(U\) and \(J\).

Step by step to compute \(U\) :

Calculate the slope of \(v^\uparrow_{\text{Hxc}}\) and \(v^\downarrow_{\text{Hxc}}\) with respect to \(\alpha\), these are the \(\delta v^\uparrow_{\text{Hxc}}\) and \(\delta v^\downarrow_{\text{Hxc}}\) that appear in the formula to compute \(U\)
Calculate the slope of the \(n^\uparrow + n^\downarrow\) with respect to \(\alpha\) this is the denominator appearing in the formula to compute \(U\)
Compute \(U\) using the formula provided above.

To compute \(J\) follow similar procedure but the derivatives are with respect to \(\beta\).

IMPORTANT: The actual \(\beta\) values in the calculations are half of the one specified in the input file.

To compute the slope, we first plot the Hxc+local for spin 1 and spin 2 as well as the occupation number against the values of \(\beta\), the same should be done with values of \(\beta\) to compute \(J\)

You can see from the plots that while the changes of the occupation numbers are perfectly linear at all \(\alpha\) values, this is not the case for the Hxc+local potential where a degree of non-linearity is present at a value of \(\alpha=0\), this is VERY important as if we were to include this data point in our calculation of \(U\), we would obtain a wrong value as our perturbation goes beyond the linear response regime.

If you discard the non-linear data point, you should obtain the following values.

\(U\) = 5.158 eV
\(J\) = 0.604 eV

What to do next

The tutorial is now complete, but you could still move forward. What can you do next?

Compute \(U\) for oxygen p states as this is commonly done in transition metal oxides, it’s usually large. For more information [Moore2024]

[Linscott2018]

E.B. Linscott, D. J. Cole, M. C. Payne, D. D. O’Regan, Phys. Rev. B 98, 235157 (2018). https://doi.org/10.1103/PhysRevB.98.235157

[Moore2024] (1,2)

G. C. Moore, M. K. Horton, E. Linscott, A. M. Ganose, Ma. Siron, D. D. O’Regan, K. A. Persson Phys. Rev. Materials 8, 014409 (2024). https://doi.org/10.1103/PhysRevMaterials.8.014409

[Moynihan2017]

G. Moynihan, G. Teobaldi, and D. D. O’Regan, A self-consistent ground-state formulation of the first-principles Hubbard U parameter validated on one-electron self- interaction error (2017), arXiv:1704.08076