benchmark against openfermion, unit issues

[overshiki]

Hi, Neat work! I just tried to do some benchmark of Quiqbox.jl against openfermion, but got some unit issues. For H2 molecule with bond length of 0.275 angstrom, the Hartree Fock energy calculated from openfermion is -0.46576 Ha, while the result I got from Quiqbox.jl is -2.53861 Ha. I think it must be a unit issues, but I'm not sure where it is. The piece of code I used for the calculation is as below:

using Quiqbox

function calc_H2()
    bond_length = 0.275
    nuc = ["H", "H"]
    nucCoords = [[0.0, 0.0, 0.0], [0.0, 0.0, bond_length]]
    bs = genBasisFunc.(nucCoords, ("STO-3G", "H") |> Ref) |> flatten

    resRHF = runHF(bs, nuc, nucCoords)
    println("Quiqbox:")
    println("bond length = ", bond_length)
    @show resRHF.Ehf
    println()
end 

using PyCall
py"""
from openfermion.chem import MolecularData
from openfermionpyscf import run_pyscf
def get_H2():
    bond_length = 0.275
    geometry = [
                ["H", [0.0, 0.0, 0.0]],
                ["H", [0.0, 0.0, bond_length]],
            ]

    basis = "sto3g"
    spin = 0

    molecule_of = MolecularData(
        geometry,
        basis,
        multiplicity=2 * spin + 1
    )
    molecule_of = run_pyscf(
        molecule_of,
        run_scf=1,
        run_ccsd=0,
        run_fci=0
    )
    print("openfermion:")
    print('At bond length of {} angstrom, molecular hydrogen has:'.format(
        bond_length))
    print('Hartree-Fock energy of {} Hartree.'.format(molecule_of.hf_energy))
""" 

calc_H2()
py"get_H2"()

The result is as below:

Quiqbox:
bond length = 0.275
resRHF.Ehf = -2.538611549467369

openfermion:
At bond length of 0.275 angstrom, molecular hydrogen has:
Hartree-Fock energy of -0.46576369202200985 Hartree.

For more detailed examples of openfermion, please refer to this

Thanks a lot in advance!

[frankwswang]

Hi @overshiki! Thanks for trying out the package! All physical values in Quiqbox are in the atomic unit, so for the bond length, you need to convert 0.275 angstrom to 0.275*1.8897259886 Bohr radius. Then you should get the correct result. If you have more questions regarding the usage of the package, you can check out the official doc. If there's anything unclear please feel free to open an issue to reach out. Also, if you want to test out the performance of the package, you might want to use the latest version on the main branch since I've been trying to improve the performance since the last release, especially for the TTFP issue that's very common for Julia packages.

[overshiki]

@frankwswang Thanks! After correcting the units, I can now get the same result as that of openfermion.

One interesting thing is, in openfermion, HF energy is, by default, the total energy calculated using HF method(including the nuclear repulsion), while in Quibox, it is treated separately. More specifically, see the result below:

Hartree-Fock Energy | Nuclear Repulsion   | Total Energy
-2.3900444821 Ha      1.9242809054 Ha       -0.4657635767 Ha

Quiqbox:
bond length = 0.519674646865
resRHF.Ehf = -2.3900444820905005

openfermion:
At bond length of 0.275 angstrom, molecular hydrogen has:
Hartree-Fock energy of -0.46576369202200985 Hartree.

Despite this little difference in the notation, the energy calculated using these two packages are quite close(the difference is about 1.15e-7Ha where the chemical accuracy is about 0.0016Ha). I also tried a larger system(LiH), the accuracy is still satisfactory. That is very good news to me, and I could try switching to this julia version of scf.

Thanks for your wonderful work!

[frankwswang]

Hi @overshiki! As you noticed in most packages such as pyscf the Hartree-Fock energy is defined with the nuclear repulsion included for convenience. However, Quiqbox is more oriented for basis set design, at least for now, so runHF will store the nuclear repulsion separately as it does not depend on the choice of basis set. If you always want to include it, you can use resRHF.Ehf + resRHF.Enn as the returned value for your wrapper function calc_H2 in the above example.

Thanks again for your appreciation! I think you opened the first-ever user issue for Quiqbox on GitHub! So you might just be the first user of this package outside the people I know. I'm excited and I'm glad I can help! I hope this package can further aid your work and if you have suggestions, you are also welcome to open issues.