[BUG] Electron repulsion energy

[ohzde]

Describe the bug

I was trying to calculate the electron–electron repulsion energy for HCl analytically using electron_repulsion_integral() from gbasis, but the results do not make sense and do not match the Gaussian output.

In Gaussian, the total electron–electron repulsion term is 182.7069747387 Hartree, but with gbasis, the Coulomb term is -2698096.1551624318 Hartree and the exchange term is 741942.5687709596 Hartree.

ee_int = electron_repulsion_integral(self.basis)

# Exchange
exch = np.einsum("mn,mnrs->rs", self.dm, ee_int)
ex_total_analytical = -0.25 * np.trace(self.dm.dot(exch))

# Coulomb
coul = np.einsum("mn,mrns->rs", self.dm, ee_int)
coul_total_analytical = 0.5 * np.trace(self.dm.dot(coul))

System Information:

• OS: Linux
• Python version: 3.13.3
• NumPy version: 2.2.2
• SciPy version: 1.15.1

Additional context

This is the Gaussian log file, and this is the Gaussian formatted checkpoint.

[PaulWAyers]

@marco-2023 and ShuoyangW, do you have any insight into this? I think you've used gbasis in meanfield before.

In my experience, I always transformed to the MO basis first, and then contracted with the (now diagonal) DM. When you do this without first contracting, one needs to be careful with how the DM is defined: in the "chemist" convention one often needs to use $\mathbf{P} \cdot \mathbf{S}$ where $\mathbf{P}$ is the DM and $\mathbf{S}$ is the overlap matrix. That would be my guess but I have not actually tested this against the example. I'll do that if there @marco-2023 or @ShuoyangW doesn't immediately recognize the issue.

[ShuoyangW]

I agree with @PaulWAyers on transforming to MO basis. Usually, what I do is
overlap = overlap_integral(basis)
overlap_inv = np.linalg.pinv(overlap)
lowdin_transform = scipy.linalg.sqrtm(overlap_inv)
ee_int = electron_repulsion_integral(basis, transform=lowdin_transform)
I believe gbasis has some much smarter functions, but I haven't carefully look through the doc.

[FarnazH]

Thanks for your comments @PaulWAyers and @ShuoyangW. The code snippets above produced the correct energies for H2 and H2O molecules, so my guess was that this might have to do with higher angular momentum basis functions.

[PaulWAyers]

Also, to be clear, this is using the native (very inefficient) gbasis implementation of the 2-electron integrals or the link to libcint? I think this will depend on how your gbasis is compiled. (pip install doesn't include libcint options yet.)

[marco-2023]

Hi @FarnazH and @ohzde,

As I understand it, for reproducing the Gaussian integrals, you need to use as the transformation the molecular orbitals coefficient. The Lowdin transformation will give you a set of orthogonal orbitals, but not necessarily the result of the computation you are using as a reference. Here is a code snippet for obtaining the electron repulsion integrals. (Also, depending on the formalism used by Gaussian you may need to change the notation to “chemist”) Can you please give it a try?

import numpy as np
from iodata import load_one
from gbasis.wrappers import from_iodata
from gbasis.integrals.electron_repulsion import electron_repulsion_integral

# load molecule info and basis from fchk file
mol_data = load_one("HCl.fchk")
ao_basis = from_iodata(mol_data) 

# extract MO coefficients from iodata object
mo_coeffs = mol_data.mo.coeffs

ee_int = electron_repulsion_integral(ao_basis, transform=mo_coeffs.T, notation='physicist')

[ohzde]

this is the full code snippet, which produces reasonable results for H₂, H₂O, and LiH, but not for HBr, HCl, and LiF (these are the only systems I've tested so far).

import numpy as np
from chemtools.wrappers.molecule import Molecule
from gbasis.wrappers import from_iodata
from gbasis.integrals.electron_repulsion import electron_repulsion_integral

molecule = Molecule.from_file("file.fchk")
basis = from_iodata(molecule._iodata)
dm = molecule._iodata.one_rdms.get("scf")

ee_int = electron_repulsion_integral(basis)

# Exchange
exch = np.einsum("mn,mnrs->rs", dm, ee_int)
ex_total_analytical = -0.25 * np.trace(dm.dot(exch))

# Coulomb
coul = np.einsum("mn,mrns->rs", dm, ee_int)
coul_total_analytical = 0.5 * np.trace(dm.dot(coul))

[PaulWAyers]

@ohzde can you share the *.fchk files for the tests you are using?

For LiF, have you tried a minimal basis set (so no polarization functions)? If the issue is with higher angular momentum that would diagnose it.

[ohzde]

thank you @PaulWAyers for your response. I just tested the LiF with RHF/STO-3G, and now the coulomb and exchange terms make sense. So I guess the issue is with the higher angular momentum?

I included the previous six .fchk files as well as this new one for LiF with STO-3G here. everything is with aug-cc-pVTZ, except water and one of the LiFs, which are with STO-3G basis sets.

[PaulWAyers]

Yeah, probably a bug in the higher angular momentum. We need to redo the 2-electron integrals....we knew that it was needed but we were being lazy about it. @marco-2023 and @msricher were going to look at the libcint interface which is reliable.

This will not be a fast fix, I suspect, so my best advice is to generate the integrals with PySCF in the interim. It may be a simple-ish bug, but digging into the e-e repulsion recursion is likely painful.

I will try to dig up my notes on this to see what, if anything, can be done in the short run.

[ohzde]

Okay, thank you! I can make it work for now by generating the integrals with PySCF, or even by comparing my numerical values with Gaussian-printed terms.