Accuracy of atoms.get_potential_energy

[nowozin]

Describe the bug

The computed potential energy for a simple CH4 molecule computed through the ASE interface seems wrong when compared against Psi4 and NWChem.

To Reproduce

I am running Ubuntu 20.04.2 LTS with ase 3.22.0 and xtb-python version 20.1. The python environment is provided by conda, using the following conda packages from the conda-forge channel: ase, psi4, xtb, xtb-python, nwchem.

I run the following xtb_test.py script:

import ase
import xtb
from ase.calculators.psi4 import Psi4
from ase.calculators.nwchem import NWChem
from xtb.ase.calculator import XTB

atoms = ase.Atoms(
    "CH4",
    [
        [-0.0034502 ,  0.01017081,  0.01938033],
        [-0.7954868 ,  0.5766599 , -0.5472012 ],
        [-0.39378393, -0.97992676,  0.2722862 ],
        [ 0.6344988 ,  0.4473651 ,  0.93568736],
        [ 0.59581804, -0.16517928, -0.8915708 ]
    ],
)

print("ase version: %s" % ase.__version__)
print("xtb version: %s" % xtb.__version__)

# From ANI-1 DFT reference data
energy_ref_Ha = -40.4805881714
print("reference: %.5f Ha" % energy_ref_Ha)

for name, calc in zip(
    ["psi4", "nwchem", "xtb"],
    [Psi4(), NWChem(), XTB(method="GFN2-xTB")]):

    atoms.calc = calc
    potential_energy_eV = atoms.get_potential_energy()
    # 27.211396641308 eV = 1 Ha
    potential_energy_Ha = potential_energy_eV / 27.211396641308

    print("'%s': %.5f Ha" % (name, potential_energy_Ha))

On my setup running this script produces the following output:

ase version: 3.22.0
xtb version: 20.1
reference: -40.48059 Ha
'psi4': -40.19151 Ha
'nwchem': -39.86620 Ha
'xtb': -4.15454 Ha

Expected behaviour

I expect that the potential energy computed for the XTB calculator approximately agrees with the energy computed by the other quantum chemistry codes, i.e. the Psi4 and NWChem calculators.

Additional context

The CH4 conformation is an out-of-equilibrium conformation from the ANI-1 DFT data set.

[awvwgk]

You can't compare absolute energies between different methods.

For context, you are comparing HF/aug-cc-pVTZ (Psi4 default), some DFT/3-21G (NWChem default), ωB97x/6–31G(d) (ANI-1) and GFN2-xTB/STO-NG. First since those are completely different methods, Hartree–Fock, GGA, RS-Hybrid and tight-binding, you can't compare absolute energies here. Second the basis set between the methods are different, for Hartree–Fock you are using an augmented triple-ζ basis, for your DFT methods you use a double-ζ basis sets and for tight-binding a valence minimal basis.

The most important point here is that while your DFT/WFT methods are using all-electron basis sets, the tight binding method is using a valence only basis set. Since the number of electrons differs between the calculations they can't be reasonably compared.

[nowozin]

Thanks Sebastian for the quick and informative answer, I understand the difference now.

I understand your argument for the potential energy, as in the above example code. In general I am also interested in force computation. For forces, does a variation of your argument still hold or does the "constant shift" drop out when forces are computed as gradients of the energy? I.e., can I expect a higher degree of agreement in forces if the dominant effects creating these forces is already accurately captured at the GFN2-xTB/STO-NG level of theory?

[awvwgk]

Energy derivatives or energy differences usually compare well between different methods.

If you are interested in the theoretical background, I can recommend Introduction to Computational Chemistry by Frank Jensen.

[nowozin]

Thanks a lot Sebastian, I very much appreciate your answer!

[awvwgk]

You're welcome.