scf total energy when U+V are non-zero is different between QE and QE+SIRIUS

[gsavva]

The total energies calculated using QE and QE-SIRIUS when the Hubbard U & V values are non-zero are different:

QE        !    total energy              =    -373.24928416 Ry
QE-SIRIUS !    total energy              =    -370.86624286 Ry

I'm using the input file of example10 (LiCoO2, non-magnetic insulator) as found in q-e/HP/examples https://gitlab.com/QEF/q-e/-/tree/develop/HP/examples.

The Hubbard U & V parameters are both changed from 0.0001 to 3.0

HUBBARD {ortho-atomic}
U Co-3d 3.0
V Co-3d O-2p 1 2 3.0

The input file, pseudopotentials, execution script and output files are available here: non_zero_U_V.zip

The version of QE used is 7.1
QE-SIRIUS was installed via spack so the branch used is q-e-sirius@develop-ristretto and sirius@develop.
The output files contain the exact commit hashes of the executables used.

For comparison, if the U & V values are not changed and thus example10 runs as provided (U, V = 0.0001), then QE and QE-SIRIUS calculate total energies as below:

QE        !    total energy              =    -373.36904960 Ry
QE-SIRIUS !    total energy              =    -373.36894085 Ry

[toxa81]

This small difference is ok. The GGA XC implementation is slighlty different between QE and LibXC. Whas is important is that forces and stress are similar.

[gsavva]

I initially started this comparison testing of forces & stresses (QE vs QE_SIRIUS) using a non-zero Hubbard U for the Hubbard atoms to evaluate (a) the total forces and (b) the Hubbard contribution to forces.
However, for some systems I had convergence issues (already communicated to @simonpintarelli ) and the total energies were different (as in the example reported above).
I'm now working on forces & stresses comparison without changing the U value (set to 0.0001 in the examples)

[gsavva]

Updating the issue with more results for different values of Hubbard U and V.
Using example10 as above, found in q-e/HP/examples.

U = 1e-4
V = 1e-4
QE     !    total energy              =    -373.36904960 Ry
SIRIUS !    total energy              =    -373.36894085 Ry

U = 3.0
V = 3.0
QE     !    total energy              =    -373.24928416 Ry
SIRIUS !    total energy              =    -370.86624283 Ry

U = 6.0
V = 0.0001
QE     !    total energy              =    -373.13921250 Ry
SIRIUS !    total energy              =    -368.92175896 Ry

U = 6.5
V = 0.0001
QE     !    total energy              =    -373.12143594 Ry
SIRIUS !    total energy              =    -368.61994388 Ry

U = 6.5
V = 0.7
QE     !    total energy              =    -373.12144331 Ry
SIRIUS !    total energy              =    -368.61994274 Ry
 

[mtaillefumier]

A factor 2 was missing in the calculation of the kinetic energy (the one-electron contribution). q-e-sirius now works fine for U \neq 0 V = 0.
U+V is still giving a different total energy though.

[mtaillefumier]

Should be fixed by PR #788. The mistake was a factor 2.0 during the calculation of the hubbard U contribution to the hubbard potential. The generalized version of Hubbard (reduce to dudarev version) is on by default.

[gsavva]

Testing results (@mtaillefumier):

HP/examples/example10 (LiCoO2, non-magnetic insulator): OK

U = 6.5
V = 0.7
QE     !    total energy              =    -373.12144331 Ry
SIRIUS !    total energy              =    -373.12142824 Ry

HP/examples/example01 (same as example10, without V): OK

U = 7.8
QE     !    total energy              =    -373.07581518 Ry
SIRIUS !    total energy              =    -373.07580009 Ry

HP/examples/example06 (Ni2MnGa, ferromagnetic metal): OK

U Mn-3d 5.0
U Ni-3d 8.3
QE     !    total energy              =    -585.95912167 Ry
SIRIUS !    total energy              =    -585.95908712 Ry
QE          total magnetization       =     5.18 Bohr mag/cell
SIRIUS      total magnetization       =    17.50 Bohr mag/cell

**Note that the energies were correct even before the fix (#788). The difference in magnetization is discussed in #786
HP/examples/example08 (NiO2, 2D system, nonmagnetic insulator): OK

U Ni-3d 8.65
QE     !    total energy              =    -167.47501933 Ry
SIRIUS !    total energy              =    -167.47572923 Ry

[toxa81]

New round of tests should be done with electronic-structure/q-e-sirius#60

[gsavva]

@toxa81 I'm checking this...

[gsavva]

Reporting below results after repeating some tests.
Codes used:

• QE: Version 7.2 Release
• SIRIUS: c871046
• q-e-sirius (QE-S): 5b06de4da8cef668c78ea60ffd9f381c044e9f01

The values in the tables are abs(total energy difference) in Ry (0 value means perfect agreement in total energy)
Used NC-PP

System	Type	U=1e-8, V=1e-8	U=5.0, V=1e-8	U=1e-8, V=5.0	U=5.0, V=5.0	source
LiCoO2	non-magnetic insulator	0	0	1e-8	6e-8	HP/examples/example10
NiO	antiferromagnetic insulator	0	5.44e-6	1.0e-3	0.51e-3	link
CTi	non-magnetic metal	0	3e-8	1e-8	6e-8	link

I'll add more soon
@toxa81

[toxa81]

Something is wrong with NiO. In CI/CD it passes dft+u+v with lower threshold for total energy.

[gsavva]

In ReFrame tests, for example this workflow, lines 359-376, there are no DFT+U +V tests that run (not yet, at least). I'm checking them locally

[gsavva]

UPDATE for NiO

Codes used:

• QE: Version 7.2 Release
• SIRIUS: 4bf993a
• q-e-sirius (QE-S): 694a67335623d128173e029e50b809157f92697b

Reporting abs(total energy difference) in Ry (0 value means perfect agreement in total energy).
Used NC-PP, and cutoffs as below:
ecutwfc = 200.0
ecutrho = 800.0

System	Type	HUBBARD	U=1e-8, V=1e-8	U=5.0, V=1e-8	U=1e-8, V=5.0	U=5.0, V=5.0	source
NiO	antiferromagnetic insulator	atomic	0	3e-8	7.2e-4	4.5e-4	link
NiO	antiferromagnetic insulator	ortho-atomic	0	5e-8	1.0e-3	5.5e-4	link

The biggest difference is still on the order of mRy

[toxa81]

I did "complex occupation matrix" experiment. For this change in QE:

diff --git a/PW/src/new_nsg.f90 b/PW/src/new_nsg.f90
index 2ee23d7ae..87b7d8cae 100644
--- a/PW/src/new_nsg.f90
+++ b/PW/src/new_nsg.f90
@@ -170,8 +170,8 @@ SUBROUTINE new_nsg()
                           ELSE
                              nrg(m2,m1,viz,na1,current_spin) = &
                              nrg(m2,m1,viz,na1,current_spin) + &
-                             DBLE( wg(ibnd,ik) * ( proj%k(off1,ibnd) * &
-                             CONJG(proj%k(off2,ibnd) * phase) ) )
+                             wg(ibnd,ik) * ( proj%k(off1,ibnd) * &
+                             CONJG(proj%k(off2,ibnd) * phase) )
                           ENDIF
                           !
                        ENDDO ! m2

I run NiO-ldapuv-ortho test with U=8.0 and V=4.3 using pure patched QE and QE/Sirius. Results are:
QE total energy: -750.95990817 Ry
QE/sirius total energy: -750.95990852 Ry

Forces looks fine. Hubbard stress has a minor difference. @mtaillefumier I would say it's a border case of being a bug.

[toxa81]

@gsavva please have a look at this PR (#1009). I run +U+V test with SIRIUS_CONFIG=cfg.json OMP_NUM_THREADS=4 pw.x -i pw.in where cfg.json contains

{
"settings" : {
 "real_occupation_matrix" : true
}
}

and got pretty good agreement for total energy, forces and stress

[gsavva]

@anton, indeed, I get almost perfect agreement in total energy with the change in PR #1009
Agreement in forces is very good (~1% difference)
Stresses almost identical (difference on the last digit)
The values below are abs(total energy difference) in Ry

System	Type	HUBBARD	U=1e-8, V=1e-8	U=5.0, V=1e-8	U=1e-8, V=5.0	U=5.0, V=5.0	source
NiO	antiferromagnetic insulator	ortho-atomic	0	3e-8	1e-8 (was 1.0e-3)	14e-8 (was 5.5e-4)	link

[gsavva]

A few more abs(total energy difference) comparisons for scf with U+V!=0

System	Type	PP	Hubbard Values	diff (Ry)
NW	non-magnetic metal	US	U=4.0 V=1.0	0
NW	non-magnetic metal	US+PAW	U=4.0 V=1.0	4.03e-4
FeO	FM metal	US	U=5.0 V=1.0	0

And almost identical stresses (difference only on the last -8th- digit)

For PAW-PP, the difference is expected since after extensive testing of DFT scf (without Hubbard), it was shown that the total energies differ by ~1e-4 Ry even for simple systems (unaries and oxides).

[gsavva]

@toxa81 I think we can finally close this...