quantum-espresso/PW/examples/example08

A simplified rotational invariant DFT+U method is presently implemented in
the pw.x code of Quantum ESPRESSO. The implemented functional is the
one proposed, among others, by S.L.Dudarev et al. in PRB, 57, 1505 (1998).

A discussion of the method, more details about the current implementation
and a description of a method to compute the Hubbard U parameter
are contained in Matteo Cococcioni's PhD thesis at SISSA and in the
paper "Linear-response approach to the calculation of the effective"
interaction parameters in the LDA+U method" by M. Cococcioni and
S. de Gironcoli, PRB 71, 035105 (2005). See also how to compute
the Hubbard parameters using the HP code whcih is part of Quantum ESPRESSO:
I. Timrov, N. Marzari, M. Cococcioni, PRB 98, 085127 (2018) and
PRB 103, 045141 (2021).

A classical example for DFT+U method is FeO that is incorrectly predicted
to be a metal by LDA and GGA while it is an insulating antiferromagnetic
material in real world.

In this example we use FeO in order to illustrate some of the input
variables involved in DFT+U calculations.

Computational parameters (as wfc and density cutoff, k-points grid etc.) 
are set so as to make the example reasonably fast and the results are
NOT meant to be converged in any sense.

The first run is just plain LDA calculation for FeO in the rhombohedral
(antiferromagnetic) cell. There are 2 types of Fe atoms in the input
because the desired magnetic structure is antiferromagnetic and opposite
starting_magnetization for the two types is suggested.  

The lda_plus_u flag is enabled in the input and a tiny value is set
for the Hubbard U of the two Iron types in order to force the code to
write out the occupation matrices for the localized Fe-3d states without
affecting the LDA result.

Looking at the output it is clear that the resulting solution is metallic:
the "correction for metal" energy term is clearly non zero and the Fermi
energy falls in the middle of the bands.

Coming to the occupation of the localized 3d states one can see that
they are completely filled for the majority spin [spin 1(up) for atom
3 and spin 2(down) for atom 4] while minority-spin components only are
partially filled and with FRACTIONAL occupations.

In the second run of the example a test value for the Hubbard U
parameter is adopted (for demonstration purposes) and the calculation 
is repeated. Note that for production calculations we suggest to compute
Hubbard U from first principles using the HP code of Quantum ESPRESSO. 
The DFT+U method is now active and disfavors fractional occupations.

In spite of that the system still, painfully, converges to a metallic
solution with similar fractional occupations as the LDA solution.  
This is due to the fact that DFT+U calculations can exhibit---even more
than spin polarized calculations do---several solutions and one is not
guaranteed to fall in the desired global minimum automatically.

We have to live with that and manage to explore several possibilities by 
suggesting to the system different starting points.
This can be done by setting the starting occupation matrices of the
system in a user defined way. 

This is done by exploiting the starting_ns_eigenvalue input variable as
in the third calculation of this example.

From literature or simple electron counting, one knows that in the
minority spin component one would like to occupy completely a single state 
leaving the other as empty as possible.

So, in the third run, by means of the starting_ns_eigenvalue variable,
one enforces the complete occupation of the fifth eigenvalue of the
minority spin components of each Fe atomic type.

Why the fifth eigenvalue? Because from the "standard DFT+U" run we know
that at the first iteration this is the one that is non-degenerate and 
if occupied completely could lead to an insulating result.

This calculation converges rather easily to the desired insulating solution.
In the output we can see that the "correction for metal" energy term is
essentially zero and Fermi energy falls in a gap.

A comment about energetics:
Plain LDA calculation has the lowest energy, as expected, since the +U
term is a positive defined penalty function added to it and energy can
only go up.
Notice however that the "standard DFT+U" calculation, the one with
fractional occupation of minority-spin levels, has a higher energy than
the "user defined ns" one, where one manages to completely fill the
desired level. This shows that this later one is indeed the ground state,
or at least, a better solution of the problem (still higher than plain
LDA, of course).

Looking at the output of these calculations one can notice that even in
the insulating solution obtained starting with user-defined ns matrices,
many of the minority spin occupations are still fractional while the DFT+U
method would like them to be either 0 or 1.

This is because the projector on localized 3d states used in the calculation
are based on atomic orbitals that are somehow different from the crystal
wavefunctions. So some "spurious" 3d level occupation comes from Oxygen
2s and 2p states that protrude toward Iron sites.

This is not wrong in general, the important thing is to be consistent
and use of the U parameter appropriate for the chosen projector.

See PP/examples/example06 for a calculation using Wannier functions.