Discrepancy between vibrational frequencies computed using xtb and xtb-python (but not with tblite) #99
Description
Overview
I'm getting different values for the vibrational frequencies when I use the standalone xtb executable with the --hess flag and when I use the Python API provided with python-xtb and the ASE Vibrations module. Interestingly, when I use the Python API of tblite, I get values that are consistent with the standalone xtb executable, suggesting that the root issue may reside with xtb-python.
Versions: xtb 6.4.1 from conda-forge (although 6.6.0 gives the same results), python-xtb 22.1 from pip, tblite 0.3.0 from pip, the master branch of ASE, Python 3.8.16, Ubuntu.
Reproducing the Discrepancy
Internal xTB vibrations
Input:
xtb c2h6.xyz --hess
c2h6.xyz.txt. Please drop the .txt extension.
Output:
projected vibrational frequencies (cm⁻¹)
eigval : -0.00 -0.00 -0.00 0.00 0.00 0.00
eigval : 341.39 873.74 873.75 1054.99 1203.59 1203.61
eigval : 1414.48 1447.57 1514.23 1514.30 1516.92 1517.01
eigval : 2971.79 2971.81 2989.26 2989.29 2989.29 2996.97
xtb-python vibrations
Input:
from ase.build import molecule
from ase.vibrations import Vibrations
from xtb.ase.calculator import XTB
atoms = molecule('C2H6') # this is identical to c2h6.xyz
atoms.calc = XTB(method='GFN2-xTB')
vib = Vibrations(atoms)
vib.run()
vib.summary()Output:
The values with <--x should be dropped due to the 3N-6 rule for vibrations.
---------------------
# meV cm^-1
---------------------
0 2.3i 18.8i <--x
1 0.0i 0.0i <--x
2 0.0i 0.0i <--x
3 50.6 408.0 <--x
4 50.9 410.7 <--x
5 87.8 708.0 <--x
6 90.4 729.0
7 108.3 873.6
8 130.7 1053.8
9 130.8 1055.1
10 131.5 1060.6
11 144.0 1161.1
12 150.8 1216.1
13 158.1 1274.9
14 175.3 1414.2
15 179.5 1447.4
16 188.0 1516.7
17 190.8 1538.6
18 370.6 2989.0
19 370.6 2989.1
20 370.8 2990.8
21 371.6 2996.9
22 371.6 2997.5
23 372.0 3000.4
---------------------
Note how the lowest real vibrational mode in the two codes do not match (341.39 vs. 729.0 cm^-1).
tblite vibrations
Input:
from ase.build import molecule
from ase.vibrations import Vibrations
from tblite.ase import TBLite
atoms = molecule('C2H6')
atoms.calc = TBLite(method='GFN2-xTB')
vib = Vibrations(atoms)
vib.run()
vib.summary()Output:
---------------------
# meV cm^-1
---------------------
0 0.0i 0.1i <--x
1 0.0i 0.1i <--x
2 0.0i 0.1i <--x
3 13.0 104.9 <--x
4 13.0 105.0 <--x
5 21.8 175.5 <--x
6 42.2 340.7
7 108.3 873.6
8 108.3 873.7
9 130.8 1055.1
10 150.1 1210.4
11 150.1 1210.4
12 175.3 1414.3
13 179.4 1447.2
14 187.8 1514.4
15 187.8 1514.7
16 188.0 1516.6
17 188.1 1516.9
18 368.5 2972.3
19 368.6 2972.6
20 370.6 2988.9
21 370.6 2989.2
22 370.6 2989.3
23 371.6 2997.1
---------------------
Note how this is largely in agreement with xtb but not with xtb-python.
Additional Debugging
Energies
As a sanity check, I have also computed the single-point energy with each code and get the following:
xtb: -7.336177090885 Eh
xtb-python: -7.3361770263 Eh
tblite: -7.3361770297 Eh
So, there is nothing obviously wrong with the energies.
Forces
I also compared forces. Again, this is all without any geometry optimization.
xtb-python:
array([[ 2.81356077e-15, -1.57075849e-05, 1.31289683e-01],
[ 1.95063714e-14, 1.57075849e-05, -1.31289683e-01],
[ 1.75860275e-15, -1.18139575e-01, -1.26627069e-01],
[ 1.02323065e-01, 5.90773242e-02, -1.26632632e-01],
[-1.02323065e-01, 5.90773242e-02, -1.26632632e-01],
[-7.42620433e-15, 1.18139575e-01, 1.26627069e-01],
[ 1.02323065e-01, -5.90773242e-02, 1.26632632e-01],
[-1.02323065e-01, -5.90773242e-02, 1.26632632e-01]])
tblite:
array([[-7.43730571e-15, -1.52828580e-05, 1.31291022e-01],
[-1.10273807e-14, 1.52828580e-05, -1.31291022e-01],
[-3.56890669e-15, -1.18139858e-01, -1.26627303e-01],
[ 1.02323311e-01, 5.90774659e-02, -1.26632866e-01],
[-1.02323311e-01, 5.90774659e-02, -1.26632866e-01],
[ 1.65071646e-15, 1.18139858e-01, 1.26627303e-01],
[ 1.02323311e-01, -5.90774659e-02, 1.26632866e-01],
[-1.02323311e-01, -5.90774659e-02, 1.26632866e-01]])
Again, there is nothing obviously different here.