I have tried to use SCAN0 in psi4, since it’s in the new libXC, but am getting somewhat different answers than Q-Chem… it looks like I’m off in the 3rd decimal place. I then tried PBE0 and found I’m off in the 5th decimal place… which still seems kind of high. I was wondering what the cause of these discrepancies could be?
Here’s one of my scan0 inputs:
molecule mol{
0 1
O 0.0 0.0 0.0
H 1.0 0.0 0.0
H 0.0 1.0 0.0
}
set {
basis 6-31G
scf_type direct
dft_spherical_points 590
dft_radial_points 99
}
Looks like you’re using a non-standard integration grid. Are you specifying that same grid in your Q-Chem inputs? The $rem block specification should look like
XC_GRID 000099000590
to request the (99, 590) grid you’ve used in the Psi4 input provided. Although, since Psi4 doesn’t prune the grids while Q-Chem does, it’s unlikely that you can compare more closely than around 10^-6 anyways.
This was psi4 1.2 downloaded from github master branch ~1 month ago. I played around with the grids in Q-Chem and Psi4 assuming that that was likely the source of error, but this didn’t really fix it. Q-Chem uses the pruned SG grids by default but I tried using the unpruned grids in my comparison (see below). I also noticed grid convergence is slow in both codes… I’m not sure how concerned I should be, but off in the 3rd decimal place seems scary.
Q-Chem input:
$molecule
0 1
O 0.0 0.0 0.0
H 1.0 0.0 0.0
H 0.0 1.0 0.0
$end
Third decimal place seems very high. Typically we see errors in the 5th decimal or so (this seemed high to me, but is somewhat typical). One note is that absolute energies are hard to compare while relative are much cleaner. We use a water interaction energy test case and a ionization energy to compare between codes. Dom may be able to point you in the right direction.
SCAN is ridiculously sensitive to the grid. 99 radial points is way too little. Try the calculation with 500 radial grid points, then you should get microhartree level agreement for absolute energies.
Relative energies may possibly be reproducible with fewer radial points.
so you actually need over 600 radial points to get the converged value, and god knows how many angular points…
The reason why the results disagree between codes is that the radial grids are slightly different, so the only way to get a reasonable comparison is to do it at the radial grid limit. The angular grid doesn’t matter if both codes use Lebedev grids.
It is quite a sad issue, since SCAN/SCAN-D3 is otherwise really good. Turbomole e.g. explicitly recommends a very high value for the radial grid to get stable gradients for weakly bound complexes (effectively something like 300-400 radial points for lighter elements).
no special sensitivity for the axial grid (beyond being a meta-GGA)