deepmodeling / abacus-develop

An electronic structure package based on either plane wave basis or numerical atomic orbitals.
http://abacus.ustc.edu.cn
GNU Lesser General Public License v3.0
173 stars 132 forks source link

Unoccupied state energy mismatches with q-e results when using PBE functionals in LibXC in PW basis #4695

Open Flying-dragon-boxing opened 3 months ago

Flying-dragon-boxing commented 3 months ago

Describe the bug

When using LibXC for planewave basis, ABACUS and quantum espresso give similar total energies, but different unoccupied state energies, which is shown as below:

State energies:

# abacus
 STATE ENERGY(eV) AND OCCUPATIONS    NSPIN == 1
 1/1 kpoint (Cartesian) = 0.0000 0.0000 0.0000 (131155 pws)
       1       -25.4444        2.00000
       2       -13.2333        2.00000
       3       -9.26499        2.00000
       4       -7.20492        2.00000
       5      -0.871762        0.00000
# qe
   -25.4454 -13.2343  -9.2660  -7.2060  -0.9226

Energies

# abacus
----------------------------------------------------------
     Energy           Rydberg                 eV          
----------------------------------------------------------
 E_KohnSham     -34.2381003066       -465.8332528647      
 E_KS(sigma->0) -34.2381003066       -465.8332528647      
 E_Harris       -34.2377708751       -465.8287707194      
 E_band         -8.1065484418        -110.2952499219      
 E_one_elec     -69.2735108555       -942.5144680993      
 E_Hartree      36.0427508013        490.3867824917       
 E_xc           -8.4227089386        -114.5968341601      
 E_Ewald        7.4153686862         100.8912669031       
 E_entropy(-TS) -0.0000000000        -0.0000000000        
 E_descf        0.0000000000         0.0000000000         
 E_exx          0.0000000000         0.0000000000         
 E_Fermi        -0.2673299563        -3.6372106521        
----------------------------------------------------------
# qe
!    total energy              =     -34.23811682 Ry
     estimated scf accuracy    <          3.1E-10 Ry
     smearing contrib. (-TS)   =      -0.00000000 Ry
     internal energy E=F+TS    =     -34.23811682 Ry

     The total energy is F=E-TS. E is the sum of the following terms:
     one-electron contribution =     -69.27308544 Ry
     hartree contribution      =      36.04224887 Ry
     xc contribution           =      -8.42264894 Ry
     ewald contribution        =       7.41536869 Ry

Expected behavior

Same orbital energies

To Reproduce

Try the input given in the additional context.

Environment

OS: Archlinux latest Compiler: gcc 14.1.1

Additional Context

Input for q-e

&CONTROL
  calculation = 'scf',
  prefix = 'pwscf',
  outdir = './',
  pseudo_dir = '/Your/Pseudo/Dir',
  disk_io = 'low',
  max_seconds = 10000000,
  verbosity = 'low',
/

&SYSTEM
  ibrav = 0,
  celldm(1) = 1,
  nat = 3,
  ntyp = 2,
  nbnd = 5,
  ecutwfc = 50,
  occupations = 'smearing',
  smearing = 'gaussian',
  degauss = 0.015,
  input_dft = 'XC-000I-000I-101L-130L-000I-000I',
!  input_dft = 'pbe',
/

&ELECTRONS
  diagonalization = 'davidson',
  electron_maxstep = 100,
  mixing_mode = 'plain',
  mixing_beta = 0.4,
  conv_thr = 1.0D-9,
/

CELL_PARAMETERS alat
28.0 0.0 0.0
0.0 28.0 0.0
0.0 0.0 28.0

ATOMIC_SPECIES
  O  1.00 O_ONCV_PBE-1.0.upf
  H  1.00 H_ONCV_PBE-1.0.upf

ATOMIC_POSITIONS alat
  O 18.966810872100 14.638172355700 6.501255519660
  H 17.984104813700 16.155190833700 6.682561292200
  H 19.500958935300 14.612241615100 4.761862131010

K_POINTS automatic
  1 1 1 0 0 0

ABACUS INPUT

INPUT_PARAMETERS
calculation scf
suffix hse
ntype 2
ecutwfc 50.000000
scf_nmax 100
basis_type pw
dft_functional GGA_X_PBE+GGA_C_PBE
nbands 5
mixing_type broyden
mixing_beta 0.80000
symmetry 0
nspin 1

out_wfc_pw 1
# out_wfc_r true
out_chg 1

pseudo_mesh            1
pseudo_rcut            10

ABACUS STRU

ATOMIC_SPECIES
O 1.00 ./O_ONCV_PBE-1.0.upf
H 1.00 ./H_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1

LATTICE_VECTORS
28.0 0.0 0.0 
0.0 28.0 0.0 
0.0 0.0 28.0 

ATOMIC_POSITIONS
Cartesian

O
0.0
1
18.966810872100 14.638172355700 6.501255519660 0 0 0
H
0.0
2
17.984104813700 16.155190833700 6.682561292200 0 0 0
19.500958935300 14.612241615100 4.761862131010 0 0 0

ABACUS STRU

ATOMIC_SPECIES
O 1.00 ./O_ONCV_PBE-1.0.upf
H 1.00 ./H_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1

LATTICE_VECTORS
28.0 0.0 0.0 
0.0 28.0 0.0 
0.0 0.0 28.0 

ATOMIC_POSITIONS
Cartesian

O
0.0
1
18.966810872100 14.638172355700 6.501255519660 0 0 0
H
0.0
2
17.984104813700 16.155190833700 6.682561292200 0 0 0
19.500958935300 14.612241615100 4.761862131010 0 0 0

ABACUS KPT

K_POINTS
0
Gamma
1 1 1 0 0 0

Task list for Issue attackers (only for developers)

kirk0830 commented 3 months ago

is it because the virtual space is not accurately/exactly diagonalized? If so, it is a known feature of QE rather than a bug. Have you tried their in-built PBE implementation?

Flying-dragon-boxing commented 3 months ago

is it because the virtual space is not accurately/exactly diagonalized? If so, it is a known feature of QE rather than a bug. Have you tried their in-built PBE implementation?

Results of built-in PBE is closer to the libxc pbe in q-e than libxc pbe in abacus

# q-e with built-in pbe
   -25.4450 -13.2339  -9.2657  -7.2056  -0.9152

     the Fermi energy is    -6.0272 ev

!    total energy              =     -34.23811398 Ry
     estimated scf accuracy    <          2.8E-10 Ry
     smearing contrib. (-TS)   =      -0.00000000 Ry
     internal energy E=F+TS    =     -34.23811398 Ry

     The total energy is F=E-TS. E is the sum of the following terms:
     one-electron contribution =     -69.27325493 Ry
     hartree contribution      =      36.04245100 Ry
     xc contribution           =      -8.42267873 Ry
     ewald contribution        =       7.41536869 Ry
Flying-dragon-boxing commented 3 months ago

is it because the virtual space is not accurately/exactly diagonalized? If so, it is a known feature of QE rather than a bug. Have you tried their in-built PBE implementation?

Results of libxc pbe with q-e input diago_full_acc set to true, couldn't see any difference.


          k = 0.0000 0.0000 0.0000 (131155 PWs)   bands (ev):

   -25.4454 -13.2343  -9.2661  -7.2061  -0.9226

     the Fermi energy is    -6.0329 ev

!    total energy              =     -34.23811682 Ry
     estimated scf accuracy    <          9.6E-10 Ry
     smearing contrib. (-TS)   =      -0.00000000 Ry
     internal energy E=F+TS    =     -34.23811682 Ry

     The total energy is F=E-TS. E is the sum of the following terms:
     one-electron contribution =     -69.27314073 Ry
     hartree contribution      =      36.04231640 Ry
     xc contribution           =      -8.42266117 Ry
     ewald contribution        =       7.41536869 Ry

     convergence has been achieved in  11 iterations
Flying-dragon-boxing commented 3 months ago

I'm currently implementing hybrid functionals for planewave, and currently we need range-separated exchange functionals in libxc to support hse functionals. In my local version of abacus I suspect this is causing band gap differences in hse calculations between qe and abacus. I've tested on more bands and this problem persists. image

WHUweiqingzhou commented 3 months ago

@pxlxingliang could you reproduce this mismatch?

kirk0830 commented 3 months ago

I have tested with QE and ABACUS with following settings: QE:

&CONTROL
calculation = 'scf'
forc_conv_thr = 0.00038
outdir = './output/'
prefix = 'VLAB'
pseudo_dir = './'
restart_mode = 'from_scratch'
tefield = .FALSE.
verbosity = 'high'
wf_collect = .TRUE.
/

&SYSTEM
degauss = 0.002
ecutrho = 400
ecutwfc = 100
ibrav = 0
lda_plus_u = .FALSE.
nat = 3
noinv = .FALSE.
noncolin = .FALSE.
nosym = .FALSE.
nspin = 1
ntyp = 2
occupations = 'fixed'
nbnd = 20
/

&ELECTRONS
conv_thr = 1e-06
electron_maxstep = 100
mixing_beta = 0.5
mixing_mode = 'plain'
scf_must_converge = .TRUE.
startingwfc = 'random'
diagonalization = 'david'
diago_full_acc = .TRUE.
/

ATOMIC_SPECIES
O 15.9994 O_ONCV_PBE-1.0.upf
H 1.00797 H_ONCV_PBE-1.0.upf

K_POINTS {GAMMA}

CELL_PARAMETERS {angstrom}
20.0000000000 0.0000000000 0.0000000000
0.0000000000 20.0000000000 0.0000000000
0.0000000000 0.0000000000 20.0000000000

ATOMIC_POSITIONS {angstrom}
O 10.27885 10 10.78217
H 10 10.49942 10
H 11.246739999999999 10.02499 10.74303

ABACUS:

INPUT_PARAMETERS
calculation scf
basis_type pw
ecutwfc 100
nspin 1
nbands 20
pseudo_dir .
out_chg 1
printe 1
gamma_only 1
pseudo_mesh 1
pseudo_rcut 10
ks_solver dav
ATOMIC_SPECIES
H 1.008 H_ONCV_PBE-1.0.upf
O 8.0 O_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1.889

LATTICE_VECTORS
20 0 0
0 20 0
0 0 20

ATOMIC_POSITIONS
Cartesian
H
0.00
2
10 10.49942 10
11.246739999999999 10.02499 10.74303

O
0.00
1
10.27885 10 10.78217

yes, I indeed find the unoccupied states seem really mismatch with QE. QE:

          k = 0.0000 0.0000 0.0000 (455957 PWs)   bands (ev):

   -25.2700 -13.0845  -9.3018  -7.2115  -0.9610   0.0464   0.2535   0.3358
     0.3510   0.3834   0.3851   0.4992   0.5459   0.6468   0.6659   0.7405
     0.7466   0.7562   0.7592   0.7624

     occupation numbers 
     1.0000   1.0000   1.0000   1.0000   0.0000   0.0000   0.0000   0.0000
     0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000
     0.0000   0.0000   0.0000   0.0000

     highest occupied, lowest unoccupied level (ev):    -7.2115   -0.9610

ABACUS:

     1                 -25.2743                        2
     2                 -13.0875                        2
     3                 -9.30289                        2
     4                 -7.21301                        2
     5                -0.959999                        0
     6               -0.0341803                        0
     7                 0.215135                        0
     8                 0.280001                        0
     9                 0.284157                        0
    10                  0.30875                        0
    11                 0.312123                        0
    12                 0.422734                        0
    13                 0.484018                        0
    14                 0.591274                        0
    15                 0.602132                        0
    16                  0.66789                        0
    17                 0.677375                        0
    18                 0.680566                        0
    19                 0.681608                        0
    20                  0.68866                        0

But I also want to inform of you about the physical meaning of those unoccupied state. For QE I can plot each unoccupied states and find many of them shows really unphysical shape. LUMO looks not bad: image But things will get worse for higher states: image image image image image My points are:

  1. yes, I also reproduce the result you observed
  2. but if you choose water molecule to benchmark some unoccupied states, it is very dangerous because you cannot easily avoid those unphysical artifact, but those will definitely mismatch with some experimental results or real-world data.