Possible discrepency of DSD-PBEP86-D3 functional energy due to P86 correlation from libxc

[ajz34]

It was found that DSD-PBEP86-D3 energy is somehow flawed. This is probably due to the P86 correlation from libxc. The following code reproduces P86 correlation energy of H2O molecule in cc-pVDZ with differenet backends libxc and xcfun.

from pyscf import gto, scf, dft, mp
ni_libxc = dft.numint.NumInt()
ni_xcfun = dft.numint.NumInt()
ni_xcfun.libxc = dft.xcfun

mol = gto.Mole(atom="N 0. 0. 0.; H .9 0. 0.; H 0. 1. 0.; H 0. 0. 1.1", basis="cc-pVDZ", verbose=0).build()
mf = dft.RKS(mol, xc="0.69*HF + 0.31*PBE, 0.44*P86").run()
dm = mf.make_rdm1()
ao = mf._numint.eval_ao(mol, mf.grids.coords, deriv=1)
rho = mf._numint.eval_rho(mol, ao, dm, xctype="GGA")

exc_libxc = ni_libxc.eval_xc(", P86", rho)[0]
exc_xcfun = ni_xcfun.eval_xc(", P86", rho)[0]
ep86_libxc = (rho[0] * exc_libxc * mf.grids.weights).sum();
ep86_xcfun = (rho[0] * exc_xcfun * mf.grids.weights).sum()
print(ep86_libxc)  # -0.35002424588085046
print(ep86_xcfun)  # -0.3500966288622522
print((ep86_libxc - ep86_xcfun) * 627.5 * 4.16)  # 0.19 kJ/mole

With matplotlib's plt.hist(exc_libxc - exc_xcfun), one can see many enrgy grids have diference at 1e-5 magnitude. For functional such as LYP, the error is as small as 1e-14 sometimes, which is close to machinary precision. So I guess something flaws for P86.

Further details of the DSD-PBEP86-D3 values will be presented in the coming comments.

At this stage, I'm not going to raise an issue to other softwares. Give me some time to figure out what's wrong. Suggestions are welcomed!

[ajz34]

A code to reproduce DSD-PBEP86-D3 is to be done in conventional SCF way, so not using dh module here.

def get_energy_bDH(mol, xc_scf, c_pt2, c_os, c_ss, libxc, d3_eng):
    mf = dft.KS(mol, xc=xc_scf)
    mf.conv_tol_grad = 1e-10
    mf.grids.atom_grid = (99, 590)
    mf._numint.libxc = libxc
    mf.run()
    if not mf.converged:
        raise ValueError("SCF not converged.")
    eng_scf = mf.e_tot
    # handle SS/OS partition
    mf_pt2 = mp.MP2(mf).run()
    t2 = mf_pt2.t2
    nocc = t2.shape[0]
    e = mf.mo_energy
    d2 = (+ e[:nocc, None, None, None] + e[None, :nocc, None, None]
          - e[None, None, nocc:, None] - e[None, None, None, nocc:])
    eng_os = np.einsum("ijab, ijab, ijab ->", t2, t2, d2)
    eng_ss = eng_os - np.einsum("ijab, ijba, ijab ->", t2, t2, d2)
    eng_pt2 = c_pt2 * (c_os * eng_os + c_ss * eng_ss)
    return eng_scf + d3_eng + eng_pt2, eng_scf + d3_eng, eng_pt2

# libxc:      -56.39727349708301
print(get_energy_bDH(mol, "0.69*HF + 0.31*PBE, 0.44*P86", 1, 0.52, 0.22, dft.libxc, -0.0004626305459352599)[0])
# xcfun:      -56.39730652290776
print(get_energy_bDH(mol, "0.69*HF + 0.31*PBE, 0.44*P86", 1, 0.52, 0.22, dft.xcfun, -0.0004626305459352599)[0])
# Gaussian:   -56.397301565299
# QChem:      -56.39730178

Input card for Gaussian:

#p DSDPBEP86(Full)/cc-pVDZ Int(Grid=99590) NoSymm

NH3 DSD-PBEP86-D3(BJ) Single Point Energy

0 1
N  0.  0.  0.
H  0.9 0.  0.
H  0.  1.  0.
H  0.  0.  1.1

! Dispersion Correction: -0.0004626305    a.u.
! SCF Done Energy      : -56.3001820336   a.u.
! PT2 Corrected Energy : -56.397301565299 a.u.

Input card for QChem:

$molecule
0 1
N  0.  0.  0.
H  0.9 0.  0.
H  0.  1.  0.
H  0.  0.  1.1
$end

$rem
JOBTYPE   sp
EXCHANGE  DSD-PBEP86-D3
BASIS     cc-pVDZ
SCF_CONVERGENCE 8
XC_GRID 000099000590
$end

[jeanwsr]

https://gitlab.com/libxc/libxc/-/merge_requests/311
https://gitlab.com/libxc/libxc/-/merge_requests/566

Are these two PRs from libxc related to this issue?