[Q-e-developers] B3LYP in QE

[Andrea Ferretti]

Dear all,

these days I've addressed a problem that I have been facing for some time, 
i.e. the validation of B3LYP results obtained by QE.
Basically, comparing the restuls of hybrid calculations done with QE (say
for a molecule) with other quantum chemistry (localized basis set) codes I
was able to get very a good agreement on the results (mostly eigenvalues)
for HF and PBE0 (besieds LDA, PBE, etc), but not as good for B3LYP.

To be more precise, I have looked at small molecules based on light
elements (such Cl2, CH4, benzene, etc), and compared with MOLPRO,
NWchem, and, to a smaller extent, Gaussian and ORCA (courtesy of Daniele
Varsano).

Common experience was that, when using a rather large basis set (such as 
aug-cc-pTVZ), the agreement of eigenvalues compared to pp-plane-waves
reached few tens of meV (from 0.01 to 0.05 eV in the worst cases, at 
least for frontier orbitals). This worked across several functionals 
(local, semilocal, hybrids but B3LYP)
Some numbers are listed below.

All in all, this analysis pointed out a problem localized in the
B3LYP.
Indeed, looking into the implementation I think that I've found at least a
couple of differences wrt the reference paper

P.J. Stephens,F.J. Devlin,C.F. Chabalowski,M.J. Frisch,
J.Phys.Chem 98, 11623 (1994)
see Eq. (2)

as well as the implementations of NWchem and libxc.

Basically:

* QE implements vwn lda-correlation according to vwn-1
   of the vwn paper,
   S.H.Vosko, L.Wilk, M.Nusair, Can.J.Phys. 58,1200(1980)
   while vwn-1-rpa has to be used for b3lyp

* the vwn-1-rpa contribution has to be included with a proper
   scaling factor and summed to the LYP correlation (both
   icorr and igcc)

after these changes have been implemented, the agreement of the QE-b3lyp 
results with other reference data is much better (of the same order of 
PBE0 or alike).

In doing so I've also implemented X3LYP, according to the ref:
X. Xu, W.A Goddard III, PNAS 101, 2673 (2004).

Results are also reported and in the same quantitative agreement with 
localized basis set codes found for other functionals.

=========

In the attachment you can find a tgz file with a patch to few QE files 
implementing the xc-functionals.

So far I have not done any commit to the trunk, since I wanted to first 
hear any comments/suggestions.

What do you think ?

take care
Andrea

=========

Cl2 molecule
#========================================================
# NWCHEM results, AVTZ basis
#========================================================
                   HOMO [eV]    LUMO [eV]      lowest eig [eV]
LDA                -7.48         -4.65           -23.39
VWN-RPA            -7.9572       -5.1218         -23.8799
PBE                -7.36         -4.44           -23.33
PBE0               -8.79         -3.36           -26.02
HF                -12.16         +0.51           -32.98

X3LYP              -8.5772       -3.5198         -25.5847
B3LYP              -8.52         -3.63           -25.44

#========================================================
# QE, plane-waves (this work)
#========================================================
                   HOMO [eV]    LUMO [eV]      lowest eig [eV]
LDA                -7.45         -4.66           -23.60
VWN-RPA            -7.9332       -5.1423         -24.0844
PBE                -7.39         -4.54           -23.52
PBE0               -8.82         -3.49           -26.18
HF                -12.27          >0             -33.11

X3LYP              -8.6243       -3.6022         -25.7891

B3LYP (FIX)        -8.58         -3.72           -25.65
B3LYP (WRONG)      -8.88         -4.03           -25.96

===========

CH4 molecule
#========================================================
# NWCHEM results, AVTZ basis
#========================================================
                    HOMO [eV]     LUMO [eV]      lowest eig [eV]
LDA                -9.4622       -3.2760         -16.9573
PBE                -9.4413       -3.6852         -17.0532
PBE0              -10.9809       +0.0368         -19.4375
HF                -14.8200       +0.8104         -25.6458

X3LYP             -10.8204       -0.1998         -19.0752
B3LYP             -10.7594       -0.2262         -18.9560


#========================================================
# QE, plane-waves (this work)
#========================================================
                    HOMO [eV]     LUMO [eV]      lowest eig [eV]
LDA                -9.4896       -0.5150         -16.9666
PBE                -9.4829       -0.5322         -16.9490
PBE0              -11.0220       -0.2528         -19.3512
HF                -14.8492        0.1920         -25.6397

X3LYP             -10.8440       -0.3028         -19.0213

B3LYP (FIX)       -10.7898       -0.3041         -18.9041
B3LYP (WRONG)     -11.1000       -0.4030         -19.2163



-- 
Andrea Ferretti, PhD
S3 Center, Istituto Nanoscienze, CNR
via Campi 213/A, 41125, Modena, Italy
Tel: +39 059 2055322;  Skype: andrea_ferretti
URL: http://www.nano.cnr.it
-------------- next part --------------
A non-text attachment was scrubbed...
Name: patch_b3lyp.tgz
Type: application/x-gtar-compressed
Size: 50581 bytes
Desc: 
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20151103/ca85a92a/attachment.bin>