get_sparse_operator does not handle DOCIHamiltonian objects as expected

[vincentelfving]

DOCIHamiltonian subclasses PolynomialTensor. That means, if a user applies get_sparse_operator to an instantiation of a DOCIHamiltonian, it will 'work' and automatically assume that it is a FermionOperator according to this check . (The reason is (I guess) that so far all such representations (QuadraticHamiltonian, InteractionOperator) which derive from PolynomialTensor were in the Fermionic context, not hard-core bosons.)

Currently, the workaround is to just not give a DOCIHamiltonian directly, but pass its property, DOCIHamiltonian.qubit_operator to the get_sparse_operator function and it behaves correctly. However, such unexpected behaviour may lead to hidden bugs because no error is raised and it returns a matrix with the right dimensions.

Suggested fix: add DOCIHamiltonian to the checks like this, which should occur before the PolynomialTensor check:

if isinstance(operator, DOCIHamiltonian):
    return qubit_operator_sparse(operator.qubit_operator, n_qubits)

(I have not checked if similar behaviour happens elsewhere apart from in get_sparse_operator)

Example snippet:

from openfermion import random_hermitian_matrix, get_sparse_operator, haar_random_vector, expectation
from openfermion.ops import DOCIHamiltonian
N = 4
seed = 42
hamiltonian = DOCIHamiltonian(constant=1.234,
                              hc=haar_random_vector(N, seed=seed).real,
                              hr1=random_hermitian_matrix(N, real=True, seed=seed+1),
                              hr2=random_hermitian_matrix(N, real=True, seed=seed+2))
sp1 = get_sparse_operator(hamiltonian)
sp2 = get_sparse_operator(hamiltonian.qubit_operator)
state = haar_random_vector(2**N, seed=seed+3)
print(expectation(sp1, state).real)
print(expectation(sp2, state).real)

0.25819950204147385 -1.0570826874100623

[obriente]

So, a DOCIHamiltonian has a representation in fermion form (where it has twice as many spin orbitals) and one in qubit form (which we have correct). It was my intention for get_sparse_operator to provide the first, not the second, and with this specified I wouldn't expect the two representations to have the same expectation value on the same state. But, there's clearly a bug here - I would expect 'print(expectation(sp1, state).real)' to raise an error here as sp1 should be a 2^8x2^8 matrix, not 2^4x2^4. I'll check it out as soon as I can.