Question on Hamiltonian Simulation

When working on this problem, I noticed my error was high. I looked around my code and noticed that the hamiltonian I had loaded wasn’t actually unitary, and if I put in a hamiltonian which is clearly unitary my error goes away. I can’t see anything wrong with my loading mechanism, does anyone have any idea what could be wrong?

from qiskit import IBMQ, Aer, assemble, transpile, execute
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.providers.ibmq import least_busy
from qiskit.quantum_info.operators import Operator, Pauli
from qiskit.algorithms import NumPyEigensolver
from numpy import linalg as LA
from scipy.linalg import expm

def is_unitary(m):
    return np.allclose(np.eye(len(m)), m.dot(m.T.conj()))

coefficients = []
bases = []
with open('LiH-Hamiltonian.txt') as f:
    for line in f:
        x = line.strip().split(" ")
        if (len(line) > 1):
            # print(x)
            c = float(x[1])
            c = (c if x[0] == '+' else -c)
            coefficients.append(c)
            bases.append(x[3])
coefficients[:10], bases[:10]
# print(list(zip(coefficients, bases)))
numBits = len(bases[0])
hamiltonian = np.zeros((2**numBits, 2**numBits), dtype=np.cdouble)
numTerms = len(bases)
for c, lbl in zip(coefficients[0:numTerms], bases[0:numTerms]):
    op = Operator(Pauli(label=lbl))
    # print(op.data)
    hamiltonian += c*op.data

print(hamiltonian)
print(is_unitary(hamiltonian))

# should give you the identity matrix if the hamiltonian is unitary, whose sum is just the length
print(np.sum(hamiltonian.dot(hamiltonian.T.conj()))) #  - np.eye(len(hamiltonian))))
print(np.sum(np.eye(len(hamiltonian))))
2 Likes

I suspect that nothing is wrong with your code. The Hamiltonian is not a unitary operator (a Hamiltonian generally is hermitian but not unitary). What is unitary is the time evolution operator U = exp(-i H t), this is the operator that you want to built a quantum circuit for.
Does that help? Or did I just misinterpret your question?