Add Qiskit-Nature qml.qchem integration

[Cmurilochem]

trafficstars

Context: This PR aims to provide a new functionality to allow for conversion between Qiskit-Nature qubit SparsePauliOp molecular Hamiltonians and Pennylane qubit Operators by adding extra functions to qml.qchem.convert module. The converters provided in this PR were adapted from the ones developed in connection with the qc2 software, a joint project between the the Netherlands eScience Center, SURF and Vrije Universiteit Amsterdam.

Description of the Change: All the changes made were uniquely done in the qml.qchem.convert module. The following extra functions were added/modified:

• _qiskit_nature_to_pennylane: converts qiskit-nature qubit molecular hamiltonians to pennylane format.
• _pennylane_to_qiskit_nature: converts pennylane qubit operators to qiskit-nature format.
• _qiskit_nature_pennylane_equivalent: checks equivalence between operators.
• import_operator: has the additional format keyword "qiskit-nature", allowing easy use of the _qiskit_nature_to_pennylane converter.

Benefits: Adds the possibility of a closer integration with Qiskit-Nature, possibly opening the doors for a complete integration in the near future with an additional module similar to openfermion_obs.py.

Possible Drawbacks: These new functions take Qiskit SparsePauliOp objects as input. So, to be able to use, the users would need to get SparsePauliOp forms from somewhere else, e.g., via a Qiskit-Nature driver or from native dataformats like QCSchema.

Related GitHub Issues: #5397

Todo:

• [ ] Add unit tests
• [ ] Add a new entry to the doc/releases/changelog-dev.md file, summarizing the change, and including a link back to the PR.

[isaacdevlugt]

Hey @Cmurilochem, thank you for your PR! I'll be taking a look at it and will get back to you by next week at the latest 🚀

[isaacdevlugt]

Hey @Cmurilochem, I might not be understanding your PR as you're intending it, but can you accomplish what you want with qml.from_qiskit_op (https://docs.pennylane.ai/en/stable/code/api/pennylane.from_qiskit_op.html)? This translates SparsePauliOp to PennyLane operators. It also looks like you're trying to do the reverse (PennyLane operators to Qiskit), and we currently don't have functionality for that. We'd be interested in having it, though!

Let me know if any of this helps!

[Cmurilochem]

Hey @Cmurilochem, I might not be understanding your PR as you're intending it, but can you accomplish what you want with qml.from_qiskit_op (https://docs.pennylane.ai/en/stable/code/api/pennylane.from_qiskit_op.html)? This translates SparsePauliOp to PennyLane operators. It also looks like you're trying to do the reverse (PennyLane operators to Qiskit), and we currently don't have functionality for that. We'd be interested in having it, though!

Let me know if any of this helps!

Hi @isaacdevlugt. Thanks for your your comments and availability in reviewing this PR.

Yes, the idea is indeed to provide qubit Hamiltonian converters from Qiskit-Nature to Pennylane and back, as you currently have for OpenFermion in convert.py.

I saw indeed the recent implementation of qml.from_qiskit_op. I also noticed that we can use the wires argument to cope with the different qubit conventions between Qiskit and PennyLane.

Specifically for qchem, one thing that I realized (which is the main reason for this PR) is that qml.from_qiskit_op may not work as expected when you have specifically as input Qiskit-Nature qubit Hamiltonian operators mapped using standard methods. e.g., JordanWignerMapper. Or maybe am I missing something ?

So, as an example, I provide below a small input where I test this (that you should be able to reproduce from my branch):

# import qiskit/qiskit-nature related packages
from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit.primitives import Estimator

# import pennylane packages
import pennylane as qml
from pennylane import numpy as np

# set up qiskit-nature pyscf driver
driver = PySCFDriver(
    atom="H 0 0 0; H 0 0 0.737166",
    basis="sto3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
)

# run qiskit-nature driver
problem = driver.run()

# get qiskit-nature fermionic hamiltonian
qtn_ferm_ham = problem.hamiltonian.second_q_op()

# get qiskit-nature qubit hamiltonian
mapper = JordanWignerMapper()
qtn_qubit_ham = mapper.map(qtn_ferm_ham) # <====== `SparsePauliOp` operator. This is the point of entrance into PennyLane

# get pennylane qubit hamiltonian from `SparsePauliOp`
pl_qubit_ham = qml.qchem.import_operator(qtn_qubit_ham, format="qiskit-nature")
#pl_qubit_ham = qml.from_qiskit_op(qtn_qubit_ham)

# Lets calculate expectation values using these hamiltonians
# 1. Qiskit-Nature

reference_state = HartreeFock(
    num_spatial_orbitals=2,
    num_particles=(1, 1),
    qubit_mapper=mapper,
)

# print("Reference circuit:\n")
# print(reference_state.draw())

ansatz = UCCSD(
    num_spatial_orbitals=2,
    num_particles=(1, 1),
    qubit_mapper=mapper,
    initial_state=reference_state
)

# print("Ansatz circuit:\n")
# print(ansatz.decompose().draw())

# set initial circ params and estimator
circuit_parameters = [1.57086611, 1.57080593, 1.45869142]
estimator = Estimator()

# get qiskit-nature expectation value
result = estimator.run(
    circuits=ansatz,
    observables=qtn_qubit_ham,
    parameter_values=circuit_parameters
).result()
qtn_expct_value = result.values[0]

print("########################################## \n")
print(f"Qiskit-Nature qubit Hamiltonian:\n {qtn_qubit_ham}\n")
print(f"Qiskit-Nature expectation value (hartree): {qtn_expct_value}")
print("########################################## \n")

# 2. PennyLane

# set up reference state and excitations
pl_reference_state = qml.qchem.hf_state(2, 4)
singles, doubles = qml.qchem.excitations(2, 4)
s_wires, d_wires = qml.qchem.excitations_to_wires(singles, doubles)
circuit_parameters = np.array([0.0, 0.0, 0.2242627774366624])

# set device and qnode
dev = qml.device("default.qubit", wires=4)
@qml.qnode(dev)
def circuit(params):
    qml.UCCSD(
        params, wires=range(4), s_wires=s_wires,
        d_wires=d_wires, init_state=pl_reference_state
    )
    return qml.expval(pl_qubit_ham)

# get expectation value
pl_expct_value = circuit(circuit_parameters)

print("########################################## \n")
print(f"PennyLane qubit Hamiltonian:\n {pl_qubit_ham}\n")
print(f"Pennylane expectation value (hartree): {pl_expct_value}")
print("##########################################")

The output is:

##########################################

Qiskit-Nature qubit Hamiltonian:
 SparsePauliOp(['IIII', 'IIIZ', 'IIZI', 'IIZZ', 'IZII', 'IZIZ', 'ZIII', 'ZIIZ', 'YYYY', 'XXYY', 'YYXX', 'XXXX', 'IZZI', 'ZIZI', 'ZZII'],
              coeffs=[-0.81125571+0.j,  0.17184932+0.j, -0.2247445 +0.j,  0.12078809+0.j,
  0.17184932+0.j,  0.16882419+0.j, -0.2247445 +0.j,  0.16605111+0.j,
  0.04526302+0.j,  0.04526302+0.j,  0.04526302+0.j,  0.04526302+0.j,
  0.16605111+0.j,  0.17454347+0.j,  0.12078809+0.j])

Qiskit-Nature expectation value (hartree): -1.855155066723801
##########################################

##########################################

PennyLane qubit Hamiltonian:
 -0.8112557075915288 * I(0) + 0.17184931866629144 * Z(0) + -0.22474449848985895 * Z(2) + 0.12078809196857414 * (Z(0) @ Z(2)) + 0.17184931866629144 * Z(1) + 0.16882419223387474 * (Z(0) @ Z(1)) + -0.2247444984898589 * Z(3) + 0.16605110981082186 * (Z(0) @ Z(3)) + 0.045263017842247726 * (Y(0) @ Y(1) @ Y(2) @ Y(3)) + 0.045263017842247726 * (Y(0) @ X(1) @ Y(2) @ X(3)) + 0.045263017842247726 * (X(0) @ Y(1) @ X(2) @ Y(3)) + 0.045263017842247726 * (X(0) @ X(1) @ X(2) @ X(3)) + 0.16605110981082186 * (Z(1) @ Z(2)) + 0.1745434714459974 * (Z(2) @ Z(3)) + 0.12078809196857414 * (Z(1) @ Z(3))

Pennylane expectation value (hartree): -1.855155078551354
##########################################

where we should be able to have equivalent hamiltonians with the same expectation values (within numerical uncertainties). However, if I replace pl_qubit_ham = qml.qchem.import_operator(qtn_qubit_ham, format="qiskit-nature") by pl_qubit_ham = qml.from_qiskit_op(qtn_qubit_ham) in the above code, I am not able to reproduce the result from Qiskit-Nature. The output is:

##########################################

Qiskit-Nature qubit Hamiltonian:
 SparsePauliOp(['IIII', 'IIIZ', 'IIZI', 'IIZZ', 'IZII', 'IZIZ', 'ZIII', 'ZIIZ', 'YYYY', 'XXYY', 'YYXX', 'XXXX', 'IZZI', 'ZIZI', 'ZZII'],
              coeffs=[-0.81125571+0.j,  0.17184932+0.j, -0.2247445 +0.j,  0.12078809+0.j,
  0.17184932+0.j,  0.16882419+0.j, -0.2247445 +0.j,  0.16605111+0.j,
  0.04526302+0.j,  0.04526302+0.j,  0.04526302+0.j,  0.04526302+0.j,
  0.16605111+0.j,  0.17454347+0.j,  0.12078809+0.j])

Qiskit-Nature expectation value (hartree): -1.855155066723801
##########################################

/Users/murilo/Library/CloudStorage/OneDrive-NetherlandseScienceCenter/Documents/NLeSC_Projects/qc2nl_quantum_computing_for_quantum_chemistry/pennylane/venv_pen/lib/python3.11/site-packages/autoray/autoray.py:81: ComplexWarning: Casting complex values to real discards the imaginary part
  return func(*args, **kwargs)
##########################################

PennyLane qubit Hamiltonian:
 (-0.8112557075915288+0j) * I(0) + (0.17184931866629144+0j) * Z(0) + (-0.22474449848985895+0j) * Z(1) + (0.12078809196857414+0j) * (Z(0) @ Z(1)) + (0.17184931866629144+0j) * Z(2) + (0.16882419223387474+0j) * (Z(0) @ Z(2)) + (-0.2247444984898589+0j) * Z(3) + (0.16605110981082186+0j) * (Z(0) @ Z(3)) + (0.045263017842247726+0j) * (Y(0) @ Y(1) @ Y(2) @ Y(3)) + (0.045263017842247726+0j) * (Y(0) @ Y(1) @ X(2) @ X(3)) + (0.045263017842247726+0j) * (X(0) @ X(1) @ Y(2) @ Y(3)) + (0.045263017842247726+0j) * (X(0) @ X(1) @ X(2) @ X(3)) + (0.16605110981082186+0j) * (Z(1) @ Z(2)) + (0.1745434714459974+0j) * (Z(1) @ Z(3)) + (0.12078809196857414+0j) * (Z(2) @ Z(3))

Pennylane expectation value (hartree): -1.245149406955886
##########################################

This behaviour arises from the way Qiskit-Nature organizes qubits after Jordan-Wigner mappings, i.e., alpha blocks followed by beta blocks. In case of PennyLane, we always have alpha, beta, alpha, beta, ...

Please, let me know what do you think. I am also happy to contribute to a separate PR implementing only the PennyLane-to-Qiskit-Nature converter.

Thanks once again.

[isaacdevlugt]

Hey @Cmurilochem, thanks for your response!

Maybe I'm missing something, but using qml.from_qiskit_op on your qtn_qubit_ham object gives what I would expect. We can look at the result from qml.from_qiskit_op and qtn_qubit_ham term by term:

qtn_ferm_ham = problem.hamiltonian.second_q_op()

# get qiskit-nature qubit hamiltonian
mapper = JordanWignerMapper()
qtn_qubit_ham = mapper.map(qtn_ferm_ham)

pl_ham =  qml.from_qiskit_op(qtn_qubit_ham)
coeffs, ops = pl_ham.terms()

qiskit_list = qtn_qubit_ham.to_list()

for i in range(len(coeffs)):
    print(ops[i], qiskit_list[i][0])
    print(coeffs[i], qiskit_list[i][1])

I() IIII
(-0.8112557075915288+0j) (-0.8112557075915288+0j)
Z(0) IIIZ
(0.17184931866629144+0j) (0.17184931866629144+0j)
Z(1) IIZI
(-0.22474449848985895+0j) (-0.22474449848985895+0j)
Z(0) @ Z(1) IIZZ
(0.12078809196857414+0j) (0.12078809196857414+0j)
Z(2) IZII
(0.17184931866629144+0j) (0.17184931866629144+0j)
Z(0) @ Z(2) IZIZ
(0.16882419223387474+0j) (0.16882419223387474+0j)
Z(3) ZIII
(-0.2247444984898589+0j) (-0.2247444984898589+0j)
Z(0) @ Z(3) ZIIZ
(0.16605110981082186+0j) (0.16605110981082186+0j)
Y(0) @ Y(1) @ Y(2) @ Y(3) YYYY
(0.045263017842247726+0j) (0.045263017842247726+0j)
Y(0) @ Y(1) @ X(2) @ X(3) XXYY
(0.045263017842247726+0j) (0.045263017842247726+0j)
X(0) @ X(1) @ Y(2) @ Y(3) YYXX
(0.045263017842247726+0j) (0.045263017842247726+0j)
X(0) @ X(1) @ X(2) @ X(3) XXXX
(0.045263017842247726+0j) (0.045263017842247726+0j)
Z(1) @ Z(2) IZZI
(0.16605110981082186+0j) (0.16605110981082186+0j)
Z(1) @ Z(3) ZIZI
(0.1745434714459974+0j) (0.1745434714459974+0j)
Z(2) @ Z(3) ZZII
(0.12078809196857414+0j) (0.12078809196857414+0j)

This looks fine to me 😄. Let me know if I'm missing something.

[Cmurilochem]

pl_ham =  qml.from_qiskit_op(qtn_qubit_ham)
coeffs, ops = pl_ham.terms()

qiskit_list = qtn_qubit_ham.to_list()

for i in range(len(coeffs)):
    print(ops[i], qiskit_list[i][0])
    print(coeffs[i], qiskit_list[i][1])

Thanks @isaacdevlugt. I appreciate your effort and input.

But still, if I use the hamiltonian directly from qml.from_qiskit_op(qtn_qubit_ham), this will not be equivalent to the one from Qiskit-Nature, and consequently they will have different expectation values; see complete example here.

I was investigating this and realised that the only way around is to play with the wires argument; see here.

In your example, this would be equivalent to:

qtn_ferm_ham = problem.hamiltonian.second_q_op()

# get qiskit-nature qubit hamiltonian
mapper = JordanWignerMapper()
qtn_qubit_ham = mapper.map(qtn_ferm_ham)

pl_ham =  qml.from_qiskit_op(qtn_qubit_ham, wires=[0, 2, 1, 3])

This reflects the different ways Qiskit-Nature and Pennylane map alpha and beta electrons. So, as an example, if I use the Jordan-Wigner mapper of Qiskit-Nature and print the resulting Qiskit's reference circuit (via reference_state.draw()), I would have:

Qiskit reference circuit:

     ┌───┐
q_0: ┤ X ├
     └───┘
q_1: ─────
     ┌───┐
q_2: ┤ X ├
     └───┘
q_3: ─────

So, for H2 in a minimal basis, you see that q_0-q_1 represent the alpha electron/qubit "block" and q_2-q_3 the beta ones. In PennyLane, the wires are represented in a more usual way (from the chemistry viewpoint): alpha, beta, alpha, beta,.... So, in this case, q_1 in Pennylane would map to q_2 in Qiskit-Nature.

That is essentially the source of all this confusion!

As an additional side note, the wire map wires=[0, 2, 1, 3] shown above should have to be different in case I use other Hamiltonian qubit mapping, e.g., BravyiKitaevMapper, which makes our life quite difficult.

I am excited to hear your thoughts on this!

[isaacdevlugt]

Hey @Cmurilochem, thanks! I'll get back to you next week as I need to confirm a few things with a colleague. Appreciate your patience! 🙏

[isaacdevlugt]

Hey @Cmurilochem! Thanks again for your patience. Indeed, you've been correct all along 😅. Apologies for not understanding that sooner. After speaking with a few of my colleagues, we recognize that there's a need for better wire mapping among different conventions, but we need some time to define how that should look in PennyLane. In the mean time, if your implementation works, please use that!

Thanks again for your incredible work! 🎉

[Cmurilochem]

Hey @Cmurilochem! Thanks again for your patience. Indeed, you've been correct all along 😅. Apologies for not understanding that sooner. After speaking with a few of my colleagues, we recognize that there's a need for better wire mapping among different conventions, but we need some time to define how that should look in PennyLane. In the mean time, if your implementation works, please use that!

Thanks again for your incredible work! 🎉

Hi @isaacdevlugt. Thank you for your time and careful consideration of the problem. I am happy to at least have triggered some discussions and raised any point for improvement. Also, let me know if I can help in anything. Always glad to contribute. Best, Carlos

[isaacdevlugt]

Hey @Cmurilochem, we'll close this PR and track this feature via your feature request here: https://github.com/PennyLaneAI/pennylane/issues/5397

Thanks again for contributing! Please reach out to us with any further questions.

[Cmurilochem]

Hey @Cmurilochem, we'll close this PR and track this feature via your feature request here: #5397

Thanks again for contributing! Please reach out to us with any further questions.

Thanks @isaacdevlugt. I am happy to help; please, feel free to contact me in case you think we may contribute further in the future. Best, Carlos