MGGA Total Energy has Unphysical Cell Convention Dependence
Hi GPAW developers,
I recently came across a large cell convention dependence of the total energy of a meta-GGA functional (specifically R2SCAN). As shown in the attached scripts (run_gpaw.sh calls run_gpaw.py for each cell and functional), I computed the total energies of silicon in the primitive and conventional cells, as well as 2x1x1 supercells of the primitive along the three lattice vectors. The deviation between primitive and conventional for PBE was roughly 0.001 meV/atom, while for R2SCAN it was 6 meV/atom. This is surprising because I used extremely dense k-point meshes and grids (the total energy for each cell is converged to well under 1 meV). In addition, the 2x1x1 and 1x2x1 supercells differ in energy by 1 meV, even though they should be identical by symmetry. This leads me to think there is a bug in the generation of the kinetic energy density and/or orbital-dependent potential.
I tried my best to find the source of the error but could not identify it; all I found out is that if the kinetic energy density is replaced with the single-orbital value |\nabla n|^2/8n
, the issue goes away, even if the pseudo-core kinetic energy is still included. This suggests the issue is probably not related to incorrect interpolation of core kinetic energy onto the FFT grid.
Also, to confirm this is not some weird issue with numerical errors getting amplified in meta-GGAs, I confirmed that there is no difference between the primitive cell and conventional cell energies in VASP (INCAR/POSCAR attached, same k-points as for GPAW); in fact the difference is <10 nano-eV/atom. For PBE the difference in VASP is 0.005 meV/atom, similar to the error found for PBE in GPAW.
Here are the energy/atom results for each functional and cell, in PW and FD mode, with the number labels for each cell below. This is all on the master branch of GPAW, with Python 3.7.
- 0=primitive
- 1=2x1x1 supercell of primitive
- 2=1x2x1 supercell of primitive
- 3=1x1x2 supercell of primitive
- 4=conventional
Struct, xc, mode: 0 PBE pw; Energy/atom: -5.399759899550386 eV
Struct, xc, mode: 1 PBE pw; Energy/atom: -5.399758400926492 eV
Struct, xc, mode: 2 PBE pw; Energy/atom: -5.3997584487273995 eV
Struct, xc, mode: 3 PBE pw; Energy/atom: -5.399759327676579 eV
Struct, xc, mode: 4 PBE pw; Energy/atom: -5.3997585409012 eV
Struct, xc, mode: 0 R2SCAN pw; Energy/atom: -8.655205185091724 eV
Struct, xc, mode: 1 R2SCAN pw; Energy/atom: -8.649115188686968 eV
Struct, xc, mode: 2 R2SCAN pw; Energy/atom: -8.650413415629014 eV
Struct, xc, mode: 3 R2SCAN pw; Energy/atom: -8.650429625007096 eV
Struct, xc, mode: 4 R2SCAN pw; Energy/atom: -8.649540306593847 eV
Struct, xc, mode: 0 PBE fd; Energy/atom: -5.400118175329186 eV
Struct, xc, mode: 1 PBE fd; Energy/atom: -5.400082264393628 eV
Struct, xc, mode: 2 PBE fd; Energy/atom: -5.40008231227739 eV
Struct, xc, mode: 3 PBE fd; Energy/atom: -5.400083191679324 eV
Struct, xc, mode: 4 PBE fd; Energy/atom: -5.400484771066662 eV
Struct, xc, mode: 0 R2SCAN fd; Energy/atom: -8.655621735570751 eV
Struct, xc, mode: 1 R2SCAN fd; Energy/atom: -8.649479367296736 eV
Struct, xc, mode: 2 R2SCAN fd; Energy/atom: -8.650778286925716 eV
Struct, xc, mode: 3 R2SCAN fd; Energy/atom: -8.65079449810575 eV
Struct, xc, mode: 4 R2SCAN fd; Energy/atom: -8.650347094996 eV
Note: The attached scripts are expensive to run due to the large k-point mesh, but the basic results can be reproduced on smaller k-point meshes. I am current running M06-L, and the problem is larger than for R2SCAN, with >20 meV/atom deviations between the primitive cell and 2x1x1 cell.
Please let me know if there is anything I can do to help, as I am happy to contribute code to help solve this bug if I can find out what is causing it.