[Pw_forum] example25 and LDA+U forces

[stefano fabris]

Dear Kevin,

first a warning and a suggestion. LDA+U calculations are trickier  
than standard LDA ones, so my suggestion is: experience with your  
system at the LDA/GGA levels first, and get a clear understand of  
where/why these levels of approximations fails. Only then, move to  
LDA/GGA+U

>
> I am new to pwscf.  I have been looking at example25, and I am  
> confused because the README does not seem to correspond to the  
> reference output files, or my results.

I have just run the example with the current CVS version and I do not  
notice relevant differences between the calculated results and those  
reported in the reference output files. If this is not the case for  
you, please provide more details on which version you are using,  
which compiler, ...

> For example, how does one conclude from feo_user_ns.out that the  
> result is an insulator?

Let's have a look at the eigenvalues of feo_user_ns.out. Here are  
those at gamma (in eV):

spin-up
     -8.7192  -7.4687   1.4531   3.6648   3.6648   5.4898   5.4898    
6.8743
      7.8273   7.8797   7.8797   8.4598   8.4598   9.8922  11.5963   
12.5867
     12.5867  13.4550  13.4550  20.0157

spin-down
     -8.7192  -7.4687   1.4531   3.6648   3.6648   5.4898   5.4898    
6.8743
      7.8273   7.8797   7.8797   8.4598   8.4598   9.8922  11.5963   
12.5867
     12.5867  13.4550  13.4550  20.0157

The calculated Fermi level is at 10.6806 eV. You can notice that  
there is a gap of ~1.7eV between the highest occupied level (9.8922)  
and the lowest unoccupied one (11.5963) hence one conclude that the  
system is insulating. You can check that there is a gap at each  
calculated k point. Consistently, the reported correction for metals  
is zero (correction for metals     =     0.00000000 ryd) and the  
calculated occupations are either close to 0 or to 1.

On the contrary let's have a look at the eigenvalues at the plain LDA  
level (feo_LDA.out), for example at the kpoint k =-0.2500-0.2500-0.2500:

     -8.4146  -8.3806   4.4378   4.8298   5.5355   5.5355   6.5045    
6.5045
      7.8487   7.8487   8.2987   9.9328   9.9328  10.9602  10.9724   
10.9724
     12.4961  12.4961  13.9943  14.2633

The Fermi level (10.9769 eV) crosses the band formed by the d states  
(10.9602  10.9724  10.9724), hence, from the point of view of the  
band structure, the oxide is formally described as metallic. As a  
result, the correction for metal is not zero (correction for  
metals     =    -0.00318382 ryd) and the corresponding calculated  
occupations are fractional, ~0.44 e.

> The "correction for metals" energy term is larger in magnitude than  
> say, feo_LDA.out.
> Also, in feo_wannier.out, I don't see how the d occupations are  
> close to 0 or 1.

Please report your results: I get the "correction for metals"  
different from zero only in the LDA calculation (feo_LDA.out), and  
the calculated occupations are always  between 0 and 1.

>
> My other question is, is the correct calculation of the forces in  
> LDA+U implemented?

Yes, forces are implemented, but only for the atomic-like projectors  
(forces with the wannier-like projectors are not implemented).  
Concerning accuracy, note that the parameter U is considered to be  
weakly dependent on atomic displacements, therefore the \partial U/  
\partial tau term in the functional derivative does not contribute to  
the forces.

Hope this can help you. All the best,

Stefano