[BUG] Big difference in PySCF vs dqc vibrational frequencies

[sofroniewn]

Describe the bug When I compute vibrational frequencies using PySCF and dqc using the same level of theory I see very different results, and I'm not sure how to interpret the dqc frequencies.

To Reproduce To compute the dqc frequencies I do

import dqc

moldesc = 'C          -1.28652        0.25524        0.08013;C           0.07762        0.16535       -0.63495;O           0.22083        0.47033       -1.80613;C           1.24790       -0.33747        0.23567;H          -2.04623        0.61618       -0.60539;H          -1.22080        0.93293        0.92681;H          -1.57890       -0.72292        0.45183;H           1.03268       -1.33368        0.61210;H           2.16265       -0.36812       -0.34714;H           1.39077        0.32217        1.08708'

mol = dqc.Mol(moldesc=moldesc, basis="sto-3g")
mol.atompos.requires_grad_()
hf = dqc.HF(mol).run()
vib = dqc.vibration(hf, freq_unit='cm^-1') # calculate vibrational frequencies
freq_dqc = vib[0].detach().numpy()
print(freq_dqc)

# Which gives
array([ 1.70294616e+04,  1.69718735e+04,  1.68860090e+04,  1.68255299e+04,
        1.68225380e+04,  1.67753119e+04,  1.44644965e+04,  9.52891028e+03,
        9.07347876e+03,  1.88365866e-04,  1.02303786e-04, -6.48131738e-05,
       -5.43764168e+03, -6.09925459e+03, -6.22499857e+03, -6.99690325e+03,
       -9.01539897e+03, -1.08179418e+04, -1.08671386e+04, -1.08868149e+04,
       -1.11379683e+04, -1.12761626e+04, -1.14500718e+04, -1.19746118e+04,
       -1.25462017e+04, -1.26196271e+04, -1.27329229e+04, -1.27626848e+04,
       -1.29649564e+04, -1.29682415e+04])

To compute the PySCF frequencies I do

from pyscf import gto
from pyscf.hessian import thermo

moldesc = 'C          -1.28652        0.25524        0.08013;C           0.07762        0.16535       -0.63495;O           0.22083        0.47033       -1.80613;C           1.24790       -0.33747        0.23567;H          -2.04623        0.61618       -0.60539;H          -1.22080        0.93293        0.92681;H          -1.57890       -0.72292        0.45183;H           1.03268       -1.33368        0.61210;H           2.16265       -0.36812       -0.34714;H           1.39077        0.32217        1.08708'

mol = gto.Mole(atom=moldesc, basis='sto-3g')
mol.build()
mf = mol.RHF().run()
hessian = mf.Hessian().kernel()
freq_pyscf = thermo.harmonic_analysis(mf.mol, hessian)['freq_wavenumber']
print(freq_pyscf)

# Which gives
array([ 117.92712198,  163.62556757,  394.58567756,  531.09155165,
        557.37323521,  905.59174999, 1100.86423523, 1103.22104045,
       1263.96843255, 1299.84707843, 1432.34457908, 1705.964129  ,
       1706.20007991, 1806.21035873, 1810.62541007, 1810.74231001,
       1814.3241799 , 2129.55729236, 3564.58835085, 3565.97676807,
       3740.70588451, 3742.69670576, 3761.55937183, 3763.21894889])

The molecule I'm using is acetone which I have done geometry optimization on in PySCF, also using HF 'sto-3g'. The values in wavenumber are roughly in line for what I'd expect for acetone based on the literature. The dqc numbers are over a wide range, both positive and negative, despite using the same level of theory.

Expected behavior I expected the dqc numbers to be more inline with the PySCF numbers

Desktop (please complete the following information):

• DQC version: 0.1.0
• Python version: 3.9
• PyTorch version: 1.10.2
• OS: macOS

Additional context I am pretty new to PySCF, dqc and computational chemistry overall so I apologies if my problem here is a real beginners mistake and I am not understanding things correctly. I hope you also don't mind all the issues I am creating - I'm really impressed with dqc and am hoping to use it in one of my projects!

[mfkasim1]

This might be due to the difference in DQC and PySCF default units. In DQC, the default length unit is Bohr while in PySCF, the default unit is Angstrom. I think a good solution is to add a method for conversions.

[mfkasim1]

I have added dqc.utils.convert_length to convert it from/to other units, so you can use it with the latest version from github. With the converter, you can try the following (with the correct unit) to give you the correct frequencies

import dqc

moldesc = 'C          -1.28652        0.25524        0.08013;C           0.07762        0.16535       -0.63495;O           0.22083        0.47033       -1.80613;C           1.24790       -0.33747        0.23567;H          -2.04623        0.61618       -0.60539;H          -1.22080        0.93293        0.92681;H          -1.57890       -0.72292        0.45183;H           1.03268       -1.33368        0.61210;H           2.16265       -0.36812       -0.34714;H           1.39077        0.32217        1.08708'

atomzs, atomposs = dqc.parse_moldesc(moldesc)
atomposs = dqc.utils.convert_length(atomposs, from_unit="angstrom")
atomposs = atomposs.requires_grad_()

mol = dqc.Mol(moldesc=(atomzs, atomposs), basis="sto-3g")
hf = dqc.HF(mol).run()
vib = dqc.vibration(hf, freq_unit='cm^-1') # calculate vibrational frequencies
freq_dqc = vib[0].detach().numpy()
print(freq_dqc)

gives

[ 3.76326944e+03  3.76089076e+03  3.74274651e+03  3.74075565e+03
  3.56602803e+03  3.56301397e+03  2.12954923e+03  1.81435037e+03
  1.81065087e+03  1.81048472e+03  1.80623622e+03  1.70622020e+03
  1.70595952e+03  1.43231587e+03  1.29985771e+03  1.26397892e+03
  1.10323001e+03  1.10087823e+03  9.05593228e+02  5.57370295e+02
  5.31093572e+02  3.94587377e+02  1.63628519e+02  1.17929438e+02
  3.30398963e+00  1.32308125e+00  3.68292417e-01  1.94050796e-03
 -1.38886749e-03 -1.56565250e-03]

It gives you all the frequencies (even though not all of them is meaningful), but the first 24 frequencies agree with PySCF

[mfkasim1]

I'm closing this issue, because it's been resolved. Please feel free to reopen if you still have the problem.

[sofroniewn]

I have added dqc.utils.convert_length to convert it from/to other units, so you can use it with the latest version from github. With the converter, you can try the following (with the correct unit) to give you the correct frequencies

Amazing thanks! I've just tried this out and it worked perfectly - I now get less than 1.5cm^-1 difference between PySCF and dqc for the 24 frequencies, which is great agreement.

Thanks for the really quick turnaround here!!