How to calcuate the Hubbard Parameter U/J for Quantum Espresso? [SOLVED]

One of the major drawbacks of the conventional DFT methods like LDA & GGA is the underestimation of band-gap.
This occurs because of the error in self-interaction energy.

This self-interaction can be described as the spurious interaction of an electron with itself. It has two main consequences:

  • Electrons are over-delocalized;
  • Band gaps in semiconductors and insulators are predicted to be much lower than their real counterpart.

Referencehttps://docs.quantumwise.com/tutorials/nio_ldau/nio_ldau.html

As a result the occupied states are pushed upwards in energy due to over delocalization.
Hartree-Fock has the same problem but for the unoccupied states, so they are over-delocalised and forced up in energy, which increases the band-gap.

Reference: https://www.researchgate.net/post/Underestimation_of_Band_gap_by_DFT2

So how do you improve your calculations to approximate the experimental band-gap more closely?
The answer is to use the Hubbard correction.

The mean-field Hubbard correction, popularly called DFT+U, is a semi-empirical correction which tries to improve on the above mentioned deficiencies.

This is particularly useful for d and f-block elements, with localised atomic-like states.
The Hubbard U introduces a potential which favours localising the states to which it is applied, thus re-localising them (compared to standard Kohn-Sham DFT) and lowering their energy. Interestingly, the Hubbard U term actually has a derivative discontinuity. Determining the appropriate U can be tricky, but once you’ve done that it is a simple calculation to perform and computationally straightforward.

Note: There exist other better ways to approximate the band-gap (ex: Hybrid functionals, GW approximation, etc.) But these are relatively computationally expensive with respect to the DFT+U.

Now that we know what Hubbard correction (DFT+U) is, we need to know how do we use it.

In this post I will be discussing on how to implement DFT+U in Quantum ESPRESSO.

Quantum Espresso,  let’s you define a Hubbard parameter, U.
First of all you need to set lda_plus_u to TRUE.

This will enable DFT+U calculations.

Now there are two types of DFT+U calculations that Quantum Espresso supports.
One requires just one parameter U, the other requires two parameters: U and J.

You can choose between the two using lda_plus_u_kind.

lda_plus_u_kind=0 for simplified version of Cococcioni and de Gironcoli, PRB 71, 035105 (2005), using Hubbard_U

lda_plus_u_kind=1 for rotationally invariant scheme of Liechtenstein et al., using Hubbard_U and Hubbard_J

Finally you just need to specify the value of the correction(Hubbard parameter) using

Hubbard_U(i)

where i stands for the ith atomic species.
Make sure that this needs to be specified in eV in contrast to most other settings.
For example if you have two atomic species Ti and O necessarily in this order, then you specify the Hubbard parameter for Ti by Hubbard_U(1)=8, and for Oxygen by Hubbard_U(2)=4.

If using the lda_plus_u_kind=0 (simplified version), then the above is enough.

If using the lda_plus_u_kind=1 (rotationally invariant version), then you also need to specify

Hubbard_J0(i)

J0 parameter (eV) for species i used in DFT+U+J calculations.

Now that we know how to provide the Hubbard parameter, we need to know how to determine this parameter.

This gets a little tricky.

There are some ways, like the linear response method, etc. which could be used for systems with single type of atoms. It is explained in detail here: http://hjkgrp.mit.edu/content/calculating-hubbard-u.

However, this is not what most of the people in the field do.
Let’s say your study involves the study of change in band-gap of a system with different amount of vacancies created.
Since your study only concerns with the band-gap, so what you can do is, that you can try different Hubbard parameters for your system until you get the desired value of your band-gap for the pure system and then continue your study.

However, you should make sure that this value doesn’t get too high. Why?

Well, as you change the Hubbard parameter, various other results get affected. Maybe the lattice parameters will now no longer be close to experimental value, or maybe you would change the phonon frequencies, and so on.

So you can choose a Hubbard parameter depending on the kind of property you wish to study, but only as long as it doesn’t drastically change some other property.

What I’ve said till now is evident from the following research papers.

https://www.nature.com/articles/srep36875

https://www.researchgate.net/publication/51554133_DFT_U_Calculations_of_Crystal_Lattice_Electronic_Structure_and_Phase_Stability_under_Pressure_of_TiO2_Polymorphs

Moreover, you should always search for the available literature on your system and see what Hubbard value are they using. Although it differs from software to software, but still it would give you some idea. Besides many people have already studied a lot of basic materials and provided the Hubbard parameters in order to perfectly simulate the experimental values.

You could also make an account on materialsproject.org and see, what settings they used. They use VASP for computations.

I hope this answers your queries, and doubts.

If there are any more doubts, leave them in the comments section down below.

References:
https://www.researchgate.net/post/How_to_calculate_Hubbard_U_parameter_for_different_atoms_in_a_given_compound

https://www.researchgate.net/post/How_to_calculate_U_and_or_J_values_of_DFT_U

https://www.nature.com/articles/srep36875

http://hjkgrp.mit.edu/content/calculating-hubbard-u

https://www.researchgate.net/profile/ANJI_KAPAKAYALA

https://www.researchgate.net/publication/51554133_DFT_U_Calculations_of_Crystal_Lattice_Electronic_Structure_and_Phase_Stability_under_Pressure_of_TiO2_Polymorphs

http://pubs.rsc.org/en/Content/ArticleHtml/2012/RA/c2ra21349d

https://arxiv.org/pdf/1703.02496.pdf

 

[wpedon id="7041" align="center"]

4 thoughts on “How to calcuate the Hubbard Parameter U/J for Quantum Espresso? [SOLVED]

  1. Dear Manas Sharma,
    I am learning Burai, I found that I cannot set U values for any elements. Would you explain why?
    Your response will be highly appreciated.

    Yunlong Gao

    1. It’s not BURAI ‘s fault. Rather its quantum espresso that doesn’t suppor5 all the elements. Although there is a quivk fix. You just need to edit two files in quantum espresso, although I’m afraid I don’t think it would be possible if youre using Burai on windows. On linux I’m not sure. But if you have your own installation of quantum espresso then the files are easy to find.

  2. Dear Manas,
    I am unable to understand what to do, I am getting a very weird error while running hp.x that says WARNING: The Fermi energy shift is too big!
    This may happen in two cases:
    1. The DOS at the Fermi level is too small:
    DOS(E_Fermi) = -0.2495E-07
    This means that most likely the system has a gap, and hence it should NOT be treated as a metal
    (otherwise numerical instabilities will appear)
    2. Numerical instabilities due to too low cutoff for hard pseudopotentials.

    my code was
    &CONTROL
    calculation = ‘scf’ ,
    restart_mode = ‘from_scratch’ ,
    outdir = ‘./tmp’ ,
    pseudo_dir = ‘.’ ,
    prefix = ‘SrTaON’ ,
    verbosity = ‘high’ ,
    /
    &SYSTEM
    ibrav = 0,
    A = 3.9575 ,
    nat = 10,
    ntyp = 4,
    ecutwfc = 55 ,
    ecutrho = 550 ,
    input_dft = ‘pbe’,
    occupations = ‘smearing’,
    smearing=’mv’,
    degauss = 0.003
    Hubbard_U(4) = 1d-08,
    Hubbard_U(1) = 1d-08,
    lda_plus_u_kind = 0,
    U_projection_type = ‘ortho-atomic’,
    lda_plus_u = .TRUE.
    /
    &ELECTRONS
    electron_maxstep = 200
    conv_thr = 1.0D-12 ,
    mixing_beta = 0.4 ,
    diagonalization = ‘david’ ,
    startingpot = ‘atomic’
    startingwfc = ‘atomic+random’
    /
    &IONS
    ion_dynamics = ‘bfgs’ ,
    /
    &CELL
    /
    CELL_PARAMETERS {alat}
    1.000000000000000 0.000000000000000 0.000000000000000
    0.000000000000000 1.453287359557151 0.000000000000000
    0.000000000000000 0.000000000000000 1.453294845517233
    ATOMIC_SPECIES
    Sr 87.62000 Sr.pbesol-spn-kjpaw_psl.1.0.0.UPF
    N 14.00650 N.pbesol-n-kjpaw_psl.1.0.0.UPF
    O 15.99900 O.pbesol-n-kjpaw_psl.1.0.0.UPF
    Ta 180.94700 Ta.pbesol-spfn-kjpaw_psl.1.0.0.UPF
    ATOMIC_POSITIONS (crystal)
    Ta 0.0000000000 0.7356262105 0.0146174163
    Ta 0.0000000000 0.2643737895 0.5146174163
    Sr 0.5000000000 0.7619621319 0.5186764475
    Sr 0.5000000000 0.2380378681 0.0186764475
    O 0.5000000000 0.7737076687 0.9873815451
    O 0.0000000000 0.9802628043 0.7192719655
    O 0.0000000000 0.0197371957 0.2192719655
    O 0.5000000000 0.2262923313 0.4873815451
    N 0.0000000000 0.4872206367 0.7600936256
    N 0.0000000000 0.5127793633 0.2600936256
    K_POINTS automatic
    4 4 4 0 0 0

    for post processing i used

    &inputhp
    outdir = ‘./tmp’ ,
    prefix = ‘SrTaON’ ,
    nq1 = 2, nq2 = 2, nq3 = 2,
    iverbosity = 2
    perturb_only_atom(4) = .true.
    start_q = 1
    last_q = 2
    /
    Please let me know how can I rectify this error

    Regards,
    Varun

  3. Will you please help how to calculate the hubbard parameter of Ca2Sb2O7?

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.