Ground energy from integrals not matching PySCF value

[mafaldaramoa]

I am trying to create a Hamiltonian in openfermion from some data (core energy and one-/two-body integrals), but the ground energy I obtain from diagonalizing the Hamiltonian does not match the FCI energy I get from PySCF. Below I leave the code I am using. The first print statement is showing the correct energy, while the second print statement is showing an incorrect value. I'm not sure what could be wrong.

import numpy as np

from pyscf.tools.fcidump import read
from pyscf import ao2mo, fci
from openfermion.chem.molecular_data import spinorb_from_spatial
from openfermion import InteractionOperator, get_sparse_operator, get_ground_state, get_fermion_operator

data = np.load('integrals.npz')

# Get core energy and one-/two-body tensors
ecore = data['ECORE']
h1 = np.array(data['H1'])
h2 = np.array(data['H2'])
norb = data["NORB"]

# Get number of alpha and beta electrons
n_a = int(data["NELEC"]/2)
n_b = int(data["NELEC"]/2)

# Get FCI energy from pyscf
cisolver = fci.direct_spin1.FCI()
ground_energy_pyscf, _ = cisolver.kernel(h1, h2, norb, (n_a,n_b), ecore=ecore, nroots=1, verbose=1)

print("Ground energy from PySCF: ", ground_energy_pyscf)

# Convert spatial integrals -> spin-orbital integrals
h1_spin, h2_spin = spinorb_from_spatial(h1, h2)

# Get the Hamiltonian and transform it to FermionOperator
h = InteractionOperator(ecore.item(), h1_spin, h2_spin)
h = get_fermion_operator(h)

# Diagonalize the Hamiltonian to get the ground energy
h_sparse = get_sparse_operator(h)
ground_energy_of, ground_state_of = get_ground_state(h_sparse)

print("Ground energy from OpenFermion: ", ground_energy_of)

[fdmalone]

Openfermion uses a different convention for the two electron integrals, if you do:

    h2 = 0.5 * np.asarray(
            h2.transpose(0, 2, 3, 1), order="C"
    )
    h1_spin, h2_spin = spinorb_from_spatial(h1, h2)
    ...

Then you should get the correct energy, I tested using a modified version of your example using random integrals as you're reading them from a file. The factor of 1/2 is from the definition of the Hamiltonian ~ 1/2 \sum_{pqrs}(pq|rs), and the ordering is openfermion expects h[p,q,r,s] = (ps|qr), where (pq|rs) is your usual chemist ERI tensor from pyscf.

[fdmalone]

here's some reference code: https://github.com/quantumlib/OpenFermion-PySCF/blob/master/openfermionpyscf/_pyscf_molecular_data.py#L71