Inconsistency between expectation values via SparseHamiltonian and get_restricted_hamiltonian?

[cvjjm]

Suppose I have an OpenFermion InteractionOperator and I want to measure its energy on a FQE wave function.

It appears there are two routs to do this:

1) Convert the InteractionOperator to a FermionOperator and then to a FQE SparseHamiltonian (btw, it would be nice if there was a more direct way to do this)

2) Get the tensors from the InteractionOperator and feed them into FQE's get_restricted_hamiltonian()

Naively I would have thought that this should yield the same result, but this does not seem to be the case, not even in the HartreeFock state as the following "minimal" example shows (just to be sure I am also constructing the InteractionOperator in two different ways):

import fqe
import numpy as np
import openfermion as of

from openfermion.chem.molecular_data import spinorb_from_spatial

def integrals_to_coefficients(one_body_integrals, two_body_integrals):
    return spinorb_from_spatial(one_body_integrals, two_body_integrals/2)

nact = 2
nalpha = 1
nbeta = 1

core_energy = 0.12345134
one_body_integrals = np.array([[-1.37560394, -0.26984563],
                               [-0.26984563,  1.09725617]])
two_body_integrals = np.array([[[[-0.56763992,  0.24989089],
                                 [ 0.24989089, -0.50890225]],
                                [[ 0.24989089, -0.50890225],
                                 [ 0.45668102, -0.21692285]]],
                               [[[ 0.24989089,  0.45668102],
                                 [-0.50890225, -0.21692285]],
                                [[-0.50890225, -0.21692285],
                                 [-0.21692285, -1.23529049]]]])

# make InteractionOperator manually
one_body_coefficients, two_body_coefficients = integrals_to_coefficients(one_body_integrals, two_body_integrals)
interaction_operator1 = of.InteractionOperator(core_energy, one_body_coefficients, two_body_coefficients)

# make InteractionOperator via MolecularData
molecule = of.MolecularData(
    geometry=[['?', [0, 0, 0]]],
    basis='cc-pvdz',
    multiplicity=nalpha-nbeta+1,
)
molecule.n_orbitals = nact
molecule.n_qubits = 2 * molecule.n_orbitals
molecule.nuclear_repulsion = core_energy
molecule.one_body_integrals = one_body_integrals
molecule.two_body_integrals = two_body_integrals
interaction_operator2 = molecule.get_molecular_hamiltonian()

wfn = fqe.Wavefunction([[nalpha+nbeta, nalpha-nbeta, nact]])
wfn.set_wfn(strategy="hartree-fock")

for interaction_operator in [interaction_operator1, interaction_operator2]:
    # btw: Is there a better way to transform an InteractionOperator into a FermionOperator?
    fermion_operator = of.FermionOperator()
    for key in interaction_operator:
        fermion_operator += of.FermionOperator(key, interaction_operator[key])
    fermion_operator = fqe.hamiltonians.sparse_hamiltonian.SparseHamiltonian(fermion_operator)
    print("SparseHamiltonian energy:           ", wfn.expectationValue(fermion_operator))

    fermion_operator = fqe.get_restricted_hamiltonian(
        ([interaction_operator.one_body_tensor,
          interaction_operator.two_body_tensor]),
        e_0=interaction_operator.constant)
    print("get_restricted_hamiltonian() energy:", wfn.expectationValue(fermion_operator))

The output I get is:

SparseHamiltonian energy:            (-3.19539646+0j)
get_restricted_hamiltonian() energy: (1.26736786+0j)
SparseHamiltonian energy:            (-3.19539646+0j)
get_restricted_hamiltonian() energy: (1.26736786+0j)

Am I doing something wrong?

[shiozaki]

I believe that fqe.get_restricted_hamiltonian assumes that the input tensors are spin-free form (i.e., spatial -- in this case (2,2) and (2,2,2,2) size). It appears that you are passing (4,4) and (4,4,4,4) sized arrays to this function. If you zero out the two body and pass one_body_integrals to get_restricted_hamiltonian, they'd agree.

Second, your two-body operator is not consistent. This is the print out - the coefficient in front of 3^3^11 should be identical to two_body_integrals[1, 1, 0, 0]. It appears that for this particular term, there is a factor of 2 missing. Also your hamiltonian does not have some elements, for example, 3^ 0^ 3 0 element, which I find odd.

Further, due to symmetry, I'd imagine that 3^ 1^ 1 3 and 3^ 1^ 3 1 should be the same magnitude but different sign. Something is odd.

...
0.124945445 [3^ 0^ 0 1] +
0.22834051 [3^ 0^ 0 3] +
-0.254451125 [3^ 0^ 2 1] +
-0.108461425 [3^ 0^ 2 3] +
-0.26984563 [3^ 1] +
0.124945445 [3^ 1^ 1 1] +
0.22834051 [3^ 1^ 1 3] +
-0.254451125 [3^ 1^ 3 1] +
-0.108461425 [3^ 1^ 3 3] +
-0.254451125 [3^ 2^ 0 1] +
-0.108461425 [3^ 2^ 0 3] +
-0.108461425 [3^ 2^ 2 1] +
-0.617645245 [3^ 2^ 2 3] +
1.09725617 [3^ 3] +
-0.254451125 [3^ 3^ 1 1] +
-0.108461425 [3^ 3^ 1 3] +
-0.108461425 [3^ 3^ 3 1] +
-0.617645245 [3^ 3^ 3 3]
SparseHamiltonian energy:            (-3.3188478+0j)
[[[[-0.56763992  0.24989089]
   [ 0.24989089 -0.50890225]]

  [[ 0.24989089 -0.50890225]
   [ 0.45668102 -0.21692285]]]

 [[[ 0.24989089  0.45668102]
   [-0.50890225 -0.21692285]]

  [[-0.50890225 -0.21692285]
   [-0.21692285 -1.23529049]]]]
get_restricted_hamiltonian() energy: (-1.2179846899999998+0j)

[cvjjm]

Well, the docstring of get_restricted_hamiltonian and the RestrictedHamiltonian suggests otherwise. Both talk about "tensors of Hamiltonian elements" and not integrals. OpenFermion uses the term "integrals" to refer to the spin-free/spatial form and "coefficients" to the tensors of numbers in front of the operators in a representation of a Hamiltonian as a polynomial of, e.g., fermionic operators.

I can confirm that if I set the two body integrals to zero and pass the integrals instead of the coefficients to get_restricted_hamiltonian(), the numbers do agree.

If get_restricted_hamiltonian() really expects integrals, then may I suggest to make the terminology in the documentation of that function consistent with the one used in OpenFermion? And maybe a dimension check could be added?

Regardless of whether I pass integrals or coefficient tensors to get_restricted_hamiltonian(), I however still cannot get the result to agree with the one obtained via the InteractionOperator -> FermionOperator -> SparseHamiltonian path, as soon as the two body integrals/tensors are non-zero.

You seem to say that my Hamiltonians are not well formed. But (at least in the second example) they are generated with standard OpenFermion tools from integrals that, while random, should have all the symmetries also shared by real RHF integrals.

To exclude that the discrepancy I observe is caused by my integrals, I have now added tests with integrals generated by openfermion-pyscf and such from the FQE unittest_data and still see the same discrepancy:

import fqe
import numpy as np
import openfermion as of
import openfermionpyscf
from fqe.unittest_data import build_lih_data
from openfermion.chem.molecular_data import spinorb_from_spatial
from pyscf import scf

def integrals_to_coefficients(one_body_integrals, two_body_integrals):
    one_body_coefficients, two_body_coefficients = spinorb_from_spatial(one_body_integrals, two_body_integrals)
    return one_body_coefficients, two_body_coefficients/2

# random integrals
core_energy = 0.
one_body_integrals = np.array([[-1.37560394, -0.26984563],
                               [-0.26984563,  1.09725617]])
two_body_integrals = np.array([[[[-0.56763992,  0.24989089],
                                 [ 0.24989089, -0.50890225]],
                                [[ 0.24989089, -0.50890225],
                                 [ 0.45668102, -0.21692285]]],
                               [[[ 0.24989089,  0.45668102],
                                 [-0.50890225, -0.21692285]],
                                [[-0.50890225, -0.21692285],
                                 [-0.21692285, -1.23529049]]]])

# make InteractionOperator manually from random integrals
one_body_integrals1, two_body_integrals1 = one_body_integrals, two_body_integrals
one_body_coefficients1, two_body_coefficients1 = integrals_to_coefficients(one_body_integrals1, two_body_integrals1)
interaction_operator1 = of.InteractionOperator(core_energy, one_body_coefficients1, two_body_coefficients1)

# make InteractionOperator from random integrals via MolecularData
one_body_integrals2, two_body_integrals2 = one_body_integrals, two_body_integrals
molecule = of.MolecularData(
    geometry=[['?', [0, 0, 0]]],
    basis='cc-pvdz',
    multiplicity=1,
)
molecule.n_orbitals = one_body_integrals.shape[0]
molecule.n_qubits = 2 * molecule.n_orbitals
molecule.nuclear_repulsion = core_energy
molecule.one_body_integrals = one_body_integrals2
molecule.two_body_integrals = two_body_integrals2
interaction_operator2 = molecule.get_molecular_hamiltonian()

# make InteractionOperator with openfermionpyscf
geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.4))]
interaction_operator3 = openfermionpyscf.generate_molecular_hamiltonian(geometry, 'sto-3g', 1, 0, 2, 2)
one_body_integrals3, two_body_integrals3 = None, None

# make InteractionOperator from FQE unit test data
one_body_integrals4, two_body_integrals4, _ = build_lih_data.build_lih_data('energy')
molecule = of.MolecularData(
    geometry=[['?', [0, 0, 0]]],
    basis='cc-pvdz',
    multiplicity=1,
)
molecule.n_orbitals = one_body_integrals.shape[0]
molecule.n_qubits = 2 * molecule.n_orbitals
molecule.nuclear_repulsion = core_energy
molecule.one_body_integrals = one_body_integrals4
molecule.two_body_integrals = two_body_integrals4
interaction_operator4 = molecule.get_molecular_hamiltonian()

for nact, nalpha, nbeta, interaction_operator, one_body_integrals, two_body_integrals, description in [
        (2, 1, 1, interaction_operator1, one_body_integrals1, two_body_integrals1, "InteractionOperator manually from random integrals"),
        (2, 1, 1, interaction_operator2, one_body_integrals2, two_body_integrals2, "InteractionOperator from random integrals via MolecularData"),
        (2, 1, 1, interaction_operator3, one_body_integrals3, two_body_integrals3, "InteractionOperator with openfermionpyscf"),
        (6, 2, 2, interaction_operator4, one_body_integrals4, two_body_integrals4, "InteractionOperator from FQE unit test data:"),
]:

    print(description)
    wfn = fqe.Wavefunction([[nalpha+nbeta, nalpha-nbeta, nact]])
    wfn.set_wfn(strategy="hartree-fock")

    # calculate energy via the InteractionOperator->FermionOperator->SparseHamiltonian route:
    # btw: Is there a better way to transform an InteractionOperator into a FermionOperator?
    fermion_operator = of.FermionOperator()
    for key in interaction_operator:
        fermion_operator += of.FermionOperator(key, interaction_operator[key])
    fermion_operator = fqe.hamiltonians.sparse_hamiltonian.SparseHamiltonian(fermion_operator)
    print("  InteractionOperator->FermionOperator->SparseHamiltonian energy: ", wfn.expectationValue(fermion_operator))

    # calculate energy with get_restricted_hamiltonian():
    if one_body_integrals is not None:
        fermion_operator = fqe.get_restricted_hamiltonian(
            ([
                one_body_integrals,
                two_body_integrals,
            ]),
            e_0=interaction_operator.constant)
        print("  get_restricted_hamiltonian() energy (called with integrals):    ", wfn.expectationValue(fermion_operator))

    fermion_operator = fqe.get_restricted_hamiltonian(
        ([
            interaction_operator.one_body_tensor,
            interaction_operator.two_body_tensor,
        ]),
        e_0=interaction_operator.constant)
    print("  get_restricted_hamiltonian() energy (called with coeffs):       ", wfn.expectationValue(fermion_operator))

I get

InteractionOperator manually from random integrals
  InteractionOperator->FermionOperator->SparseHamiltonian energy:  (-3.3188478+0j)
  get_restricted_hamiltonian() energy (called with integrals):     (0.3152385+0j)
  get_restricted_hamiltonian() energy (called with coeffs):        (1.14391652+0j)
InteractionOperator from random integrals via MolecularData
  InteractionOperator->FermionOperator->SparseHamiltonian energy:  (-3.3188478+0j)
  get_restricted_hamiltonian() energy (called with integrals):     (0.3152385+0j)
  get_restricted_hamiltonian() energy (called with coeffs):        (1.14391652+0j)
InteractionOperator with openfermionpyscf
  InteractionOperator->FermionOperator->SparseHamiltonian energy:  (-7.860538661020691+0j)
  get_restricted_hamiltonian() energy (called with coeffs):        (-7.4798403623090595+0j)
InteractionOperator from FQE unit test data:
  InteractionOperator->FermionOperator->SparseHamiltonian energy:  (-13.178585813739+0j)
  get_restricted_hamiltonian() energy (called with integrals):     (-8.857341498221992+0j)
  get_restricted_hamiltonian() energy (called with coeffs):        (-12.421481015593002+0j)

but still think that the energy from InteractionOperator->FermionOperator->SparseHamiltonian should agree with either the value returned by get_restricted_hamiltonian() when called with integrals or when called with coefficients.

[shiozaki]

Regarding the first point, I will include the docstring update in the next commit to clarify that what's expected is a Hamiltonian tensor in a spin-free form.

Regarding the second, get_restricted_hamiltonian is correctly working (see below); therefore, I can only speculate that there is some inconsistency in your code, perhaps with regards to interacting with OpenFermion and/or FQE.

Nick @ncrubin -- could you take a look?

import numpy
import fqe 
from fqe.wavefunction import Wavefunction
from fqe.unittest_data import build_lih_data

"""Checking total energy with LiH
"""
eref = -8.877719570384043
norb = 6 
nalpha = 2 
nbeta = 2 
nele = nalpha + nbeta
h1e, h2e, lih_ground = build_lih_data.build_lih_data('energy')

elec_hamil = fqe.get_restricted_hamiltonian((h1e, h2e))
wfn = Wavefunction([[nele, nalpha - nbeta, norb]])
wfn.set_wfn(strategy='from_data',
            raw_data={(nele, nalpha - nbeta): lih_ground})

ecalc = wfn.expectationValue(elec_hamil)
print(ecalc)
print(eref)

which returns

(-8.87771957038547+0j)
-8.877719570384043

[cvjjm]

I can also reproduce the problem with a more minimal example starting from your code:

import numpy
import fqe
from fqe.wavefunction import Wavefunction
from fqe.unittest_data import build_lih_data
import openfermion as of

"""Checking total energy with LiH
"""
eref = -8.877719570384043+0j
norb = 6
nalpha = 2
nbeta = 2
nele = nalpha + nbeta
h1e, h2e, lih_ground = build_lih_data.build_lih_data('energy')

elec_hamil = fqe.get_restricted_hamiltonian((h1e, h2e), e_0=0.0)
print("elec_hamil", elec_hamil)
wfn = Wavefunction([[nele, nalpha - nbeta, norb]])
wfn.set_wfn(strategy='from_data',
            raw_data={(nele, nalpha - nbeta): lih_ground})

ecalc = wfn.expectationValue(elec_hamil)
print(ecalc)
print(eref)

# now we try to reach the same result via the InteractionOperator->FermionOperator->SparseHamiltonian route:
molecule = of.MolecularData(
    geometry=[['?', [0, 0, 0]]],
    basis='cc-pvdz',
    multiplicity=1,
)
molecule.n_orbitals = h1e.shape[0]
molecule.n_qubits = 2 * molecule.n_orbitals
molecule.nuclear_repulsion = 0.0
molecule.one_body_integrals = h1e
molecule.two_body_integrals = h2e
interaction_operator = molecule.get_molecular_hamiltonian()
fermion_operator = of.FermionOperator()
for key in interaction_operator:
    fermion_operator += of.FermionOperator(key, interaction_operator[key])
fermion_operator = fqe.hamiltonians.sparse_hamiltonian.SparseHamiltonian(fermion_operator)
ecalc2 = wfn.expectationValue(fermion_operator)
print(ecalc2)

assert numpy.isclose(ecalc, ecalc2)

Here, I think, the correctness of the integrals is beyond doubt and I am only using standard OpenFermion functions to convert the given integrals into an InteractionOperator. I really don't see where my code can explain the result that I am getting:

(-8.87771957038547+0j)
(-8.877719570384043+0j)
(-13.129513263353799+0j)
Traceback (most recent call last):
  File "fqe_test2.py", line 45, in <module>
    assert numpy.isclose(ecalc, ecalc2)
AssertionError

[shiozaki]

Thank you, this example is helpful!

To summarize the situation for others: there seems some inconsistency in the conversion process of numpy.ndarray->InteractionOperator->FermionOperator->SparseHamiltonian, and I have observed in the initial reply that the SparseHamiltonian object does not seem to be properly formed. The first two steps are implemented in OpenFermion (that I am not familiar with) and the third step is implemented in FQE. I can double-check the third step, should we narrow down the problem to this step.

I will first let Nick @ncrubin review this example to get his insights.

[cvjjm]

Great!

There are all sorts of subtleties here concerning tensor index ordering and subsystem ordering (both qubits and fermions) but from reading the examples (and the great overlap in developers of the two packages), my impression was that OpenFermion and FQE should use all the same conventions, which is why I am puzzled by the discrepancy described above.

[ncrubin]

@cvjjm Unfortunately, the default convention used in OpenFermion is not consistent with FQE (which was designed for efficiency). I'll take a look at the examples and report back.

[cvjjm]

Just to add a little more to this: The following code computes the energy from the cirq wave function and the Hamiltonian obtained from OF:

qubit_op = of.transforms.jordan_wigner(fermion_operator)
hamiltonian_mat = of.get_sparse_operator(qubit_op).toarray()
cirq_state = fqe.to_cirq(wfn)
ecalc3 = cirq_state.conj().T.dot(hamiltonian_mat @ cirq_state)
print(ecalc3)

The result is also -13.129513263353799+0j, i.e., in line with result from fqe.hamiltonians.sparse_hamiltonian.SparseHamiltonian but again different from the get_restricted_hamiltonian() result.

My suspicion thus was that the problem is caused by some difference in how qubits or orbitals are ordered during the jordan wigner transformation. I have some code that does JW in 'sequential' ordering and the qubit ordering of the resulting operator can be changed with with openfermion.transforms.reorder, which allowed me to check all four combinations of sequential and interleaved, but none of the four results agrees with get_restricted_hamiltonian():

(-13.129513263353813+0j)
(-13.135429661233069+0j)
(-7.592300238551463+0j)
(-7.59127391423463+0j)

[cvjjm]

I think I found the solution to the mystery here.

FQE seems to use "chemists' ordering" of the two electron integrals and conversion from OpenFermion convention (which uses "physicist's ordering" can be done by means of numpy.einsum('ijlk', -0.5 * h2e) (notice also the 0.5).

One can thus go back from FQE convention to OF convention by means of numpy.einsum('ijlk', -2 * h2e). Passing two electron integrals treated in this way into the InteractionOperator->FermionOperator->SparseHamiltonian route yields energies that are consistent with the results from fqe.get_restricted_hamiltonian().

I wish this were documented somewhere, or even better that at least the packages from the same developers stick to the same conventions....

[ncrubin]

Hi @cvjjm , circling back to this in more detail today. Indeed we provide the openfermion functions for interfacing the two Hamiltonian types. You are right that ideally documentation on OpenFermion and FQE highlights (the very first page) integral formatting. I will raise this as an issue and get someone to address it (myself or @jjgoings) Indeed it would be nice if the ordering was the same for OpenFermion and FQE. Just as a note, OpenFermion isn't even physics notation as I've traditionally seen it. The OpenFermion convention is to use the pqrs ladder op indices as the tensor indices. This is not the case in most quantum chemistry codes.

Another common way to build hamiltonians for FQE is to use the build_hamiltonian function. We provide this as a high level interface that accepts FermionOperators and returns relevant Hamiltonians. It is semi-smart and can recognize if you provide a general Hamiltonian, a quadratic Hamiltonian, or diagonal coulomb. building a restricted hamiltonian using conversion from OpenFermion is the best option.

For example, augmenting your original code slightly

import numpy as np
import openfermion as of

from openfermion.chem.molecular_data import spinorb_from_spatial

def integrals_to_coefficients(one_body_integrals, two_body_integrals):
    return spinorb_from_spatial(one_body_integrals, two_body_integrals/2)

nact = 2
nalpha = 1
nbeta = 1

core_energy = 0.12345134
one_body_integrals = np.array([[-1.37560394, -0.26984563],
                               [-0.26984563,  1.09725617]])
two_body_integrals = np.array([[[[-0.56763992,  0.24989089],
                                 [ 0.24989089, -0.50890225]],
                                [[ 0.24989089, -0.50890225],
                                 [ 0.45668102, -0.21692285]]],
                               [[[ 0.24989089,  0.45668102],
                                 [-0.50890225, -0.21692285]],
                                [[-0.50890225, -0.21692285],
                                 [-0.21692285, -1.23529049]]]])

# make InteractionOperator manually
one_body_coefficients, two_body_coefficients = integrals_to_coefficients(one_body_integrals, two_body_integrals)
interaction_operator1 = of.InteractionOperator(core_energy, one_body_coefficients, two_body_coefficients)
fermion_operator1 = of.get_fermion_operator(interaction_operator1)
fqe_ham = fqe.build_hamiltonian(fermion_operator1, conserve_number=True,
                                norb=nact)

# make InteractionOperator via MolecularData
molecule = of.MolecularData(
    geometry=[['?', [0, 0, 0]]],
    basis='cc-pvdz',
    multiplicity=nalpha-nbeta+1,
)
molecule.n_orbitals = nact
molecule.n_qubits = 2 * molecule.n_orbitals
molecule.nuclear_repulsion = core_energy
molecule.one_body_integrals = one_body_integrals
molecule.two_body_integrals = two_body_integrals
interaction_operator2 = molecule.get_molecular_hamiltonian()

wfn = fqe.Wavefunction([[nalpha+nbeta, nalpha-nbeta, nact]])
wfn.set_wfn(strategy="hartree-fock")

for interaction_operator in [interaction_operator1, interaction_operator2]:
    # btw: Is there a better way to transform an InteractionOperator into a FermionOperator?
    fermion_operator = of.FermionOperator()
    for key in interaction_operator:
        fermion_operator += of.FermionOperator(key, interaction_operator[key])
    fermion_operator = fqe.hamiltonians.sparse_hamiltonian.SparseHamiltonian(fermion_operator)
    print("SparseHamiltonian energy:           ", wfn.expectationValue(fermion_operator))

    fermion_operator = fqe.get_restricted_hamiltonian(
        ([interaction_operator.one_body_tensor,
          interaction_operator.two_body_tensor]),
        e_0=interaction_operator.constant)
    print("get_restricted_hamiltonian() energy:", wfn.expectationValue(fermion_operator))

print("FQE Ham energy : ", wfn.expectationValue(fqe_ham).real)

[cvjjm]

Great! And thanks for making me aware of build_hamiltonian(). Feel free to close this issue or leave it open until the changes suggested above have been implemented.

[ncrubin]

I will leave it open until changes are made. This is a good first issue for @jjgoings