For B3LYP, Gaussian and GAMESS use the correlation functionals VWN_3 and VWN_5, respectively, instead of Psi4’s default VWNrpa_3. Of course, total energies of identical geometries differ between the applications.
If the ./python/functional.py:build_b3lyp_superfunctional() is modified accordingly, Gaussian and GAMESS energies are reproduced with 7 decimal point agreement for a few small test molecules.
The B3LYP_5 builder is currently commented out, noted as broken, although it was to have used the RPA parameterization of VWN_5.
Does it make sense to proliferate two more B3LYP variants?
What am I missing?
As far as I knew, Psi4 was using the industry-standard vwn3 for its b3lyp. The b3lyp5 is commented out because although nothing looked wrong with the setup in functional.py, b3lyp3 and b3lyp5 energies were (wrongly) equal. I wasn’t aware that there was any proliferation of two more b3lyp variants. Guess will have to investigate. You say psi4’s b3lyp (w/vwn_3) matches gaussian’s to 7 and psi4’s b3lyp (w/vwn_5) matches gamess’s to 7. How many decimals does psi4’s b3lyp (w/vwnrpa_3) match gaussian’s to?
For water at 6-31G(2df,p), they differ by 0.03 Ha , about 20 kcal in an isodesmic reaction.
Psi4 -76.4293542
GAMESS -76.392202
For small molecules, GAMESS and Gaussian’s B3LYP @ (random basis set) values are often the same to 6 decimals. GAMESS includes a Gaussian compliant version named B3LYPv3. Nwchem doesn’t bother but permits a completely user-defined functional, viz.,
I noticed this when implementing G4 and G4(MP2) under Nwchem since the optimization+ZPE steps require B3LYP. A good exegesis on B3LYP variabilty, B3LYP and WAH – the confusion, is at the Lindqvist blog.
The minimal addition is to copy build_b3lyp_3_superfunctional() from the default b3lyp and change ‘VWNRPA_3’ to ‘VWN_3’.
Backstory: Curtiss used Gaussian-94 while authoring the G4 and G4(MP?) methods. Its B3LYP was apparently different from the g03 and g09 standard.
The Lindqvist blog shows that “B3LYP” is an arbitrary term, the sort where you squint your eyes shut, massage your temples, and dose with Advil. As I attempt to replicate G4,G4(MP2) in Psi4, the phrase “defect compatibility” seems more factual than confrontational.
Here are some values for water at 6-31G(2df,p) using the same geometry:
rOH = 0.9622
aHOH = 103.84
Sloppy table follows. Pardon my unfamiliarity with this forum’s markup.
A B3LYP-? functional that gets within 0.1 mHa (~0.6 kcal) of GAMESS or NWCHEM is satisfactory. I’ll do the tests and post a candidate patch for examination.
For VWN3 (same as Gaussian) I get
-76.46290368244 (Turbomole “B3LYP_Gaussian”)
-76.46290355576 (psi4 energy(‘b3lyp’))
-76.46290344335 (nwchem default vnw-1-rpa)
VNW-1-rpa in nwchem seems to be VWN3 in Turbomole/Gaussian
and VWN5
-76.42582290565 (Turbomole default )
-76.42582303608 (psi4 energy(‘b3lyp5’) with VWN5_C not VWN5RPA_C , needs to be changed!)
-76.42582292329 (nwchem VWN_5)
besides the VWN5_C change it seems all okay for PSI4. Better agreement will be reached with tighter grids, i suppose.
Note: Gaussian does mix cartesian/spherical basis functions for the pople series (you can enforce all-spherical treatment with 5D keyword i think).
my water
O -1.769142 -0.076181 0.000000
H -2.065745 0.837492 0.000000
H -0.809034 0.001317 0.000000
I have gone through this several times with grad
students, so let me present a few results for comparison.
I am using 6-311g(d) to avoid any issues of cartesian
vs. spherical in the comparison.
Using G09 and this input
#p b3lyp/6-311g(d) int=ultrafine test
Water
0,1
O
H 1 rOH
H 1 rOH 2 aHOH
rOH = 0.9622
aHOH = 103.84
gives
SCF Done: E(RB3LYP) = -76.4338100903 A.U. after 9 cycles
O 8 0.000000 0.000000 0.118690
H 1 0.000000 0.757396 -0.474758
$end
gives
FINAL R-B3LYPV1R ENERGY IS -76.4338100577 AFTER 12 ITERATIONS
which shows that GAMESS and Gaussian can give results that agree
to better than 1 microhartree if suitable input options are chosen.
Using a recent psi4 binary and this input
molecule water {
O 0.000000 0.000000 0.118690
H 0.000000 0.757396 -0.474758
H 0.000000 -0.757396 -0.474758
}
set basis 6-311g(d)
set {
scf_type pk
dft_spherical_points 590 # Often needed
dft_radial_points 99 # Often needed
}
energy(“b3lyp”)
gives
@RKS Final Energy: -76.43381003672118
which also agrees with the previous results to better than 1 microhartree.
So if the goal is to reproduce Gaussian B3LYP calculations with psi4, you
seem to be there already.
Spherical handling of the basis set and (very) fine grid size would seem to be the obligatory bits.
Those changes made, Nwchem ,too, reproduces your values.
I’d never noticed the JANS dft keyword in GAMESS. Knowing to look for ‘pruned grids’, it’s clearly important.
Thanks everyone for your helpful reference values and sorting out of some issues. @hokru, what basis set are your numbers at?
I’ve started to fix this issue at https://github.com/psi4/psi4/pull/339 . Basically, fixed up B3LYP5 so it gives the right answers and re-enabled it. Just still uncertain on the reasoning behind some psi4 internals.
We’ve nearly got this straightened out. If anyone with access to Gamess, Turbomole, Gaussian would run open-shell calcs at the geom and basis described here and report back energies, I’d be obliged. UKS, not ROKS.