Qiskit-Nature `qml.qchem` integration

[Cmurilochem]

Feature details

Release 0.35.0 adds Qiskit 1.0 integration via the implementation of from_qiskit and from_qiskit_op functions. Yet, in the context of quantum chemistry, and specifically Qiskit-Nature, circuits (activate spaces) and operators are defined in a slightly different, specific manner: in addition to the well-known right-to-left format of Qiskit operators, in Qiskit-Nature, alpha and beta electrons are defined in separate "blocks" and not sequentially as in PennyLane. So, from_qiskit and from_qiskit_op are not enough to convert, e.g., qubit Hamiltonian operators from Qiskit-Nature to qml Operator.

To this end and similarly to _openfermion_to_pennylane, I suggest the implementation of a _qiskit_nature_to_pennylane function in connection with qchem/convert.py.

Implementation

Based on the format differences between Qiskit-Nature and PennyLane discussed above, the new qml.qchem.convert._qiskit_nature_to_pennylane function would assume something like:

def _qiskit_nature_to_pennylane(qubit_operator, wires=None):
    """Convert Qiskit SparsePauliOp to 2-tuple of coeffs and PennyLane Paulis.

    This functions is useful to convert fermionic-to-qubit transformed
    qchem operators from Qiskit-Nature to Pennynale format.

    Args:
        qubit_operator (qiskit.quantum_info.SparsePauliOp): Qiskit operator
            representing the qubit electronic Hamiltonian from Qiskit-Nature.
        wires (Wires, list, tuple, dict): Custom wire mapping used to convert
            the qubit operator to an observable terms measurable in PennyLane.
            For types Wires/list/tuple, each item in the iterable represents a
            wire label corresponding to the qubit number equal to its index.
            For type dict, only int-keyed dict (for qubit-to-wire conversion)
            is accepted. If None, will use identity map (0->0, 1->1, ...).

    Returns:
        tuple[array[float], Iterable[pennylane.operation.Operator]]: coeffs
        and their corresponding PennyLane observables in the Pauli basis.

    **Example**

    >>> from qiskit.quantum_info import SparsePauliOp
    >>> qubit_op = SparsePauliOp.from_list([("XIIZI", 1), ("IYIIY", 2)])
    >>> qubit_op
    SparsePauliOp(['XIIZI', 'IYIIY'],
              coeffs=[1.+0.j, 2.+0.j])
    >>> _qiskit_nature_to_pennylane(qubit_op,wires=['w0','w1','w2','w3','w4'])
    (tensor([1., 2.], requires_grad=False),
    [PauliX(wires=['w2']) @ PauliZ(wires=['w3']),
    PauliY(wires=['w0']) @ PauliY(wires=['w4'])])

    If the new op-math is active, the list of operators will be cast as
    :class:`~.Prod` instances instead of :class:`~.Tensor` instances when
    appropriate.
    """
    n_wires = qubit_operator.num_qubits
    wires = _process_wires(wires, n_wires=n_wires)
    wire_map = {wires[x]: y for x, y in zip(range(n_wires), range(n_wires))}

    if qubit_operator.coeffs.size == 0:
        return np.array([0.0]), [qml.Identity(wires[0])]

    def _get_op(term, wires):
        """A function to translate Qiskit to Pennylane Pauli terms."""
        if len(term) == n_wires:

            # the Pauli term '...XYZ' in Qiskit is equivalent to [Z0 Y1 X2 ..]
            # in Pennylane. So, invert the string..
            term = term[::-1]
            # wires in Qiskit-Nature are grouped by separated alpha and beta
            # blocks, e.g., for H2 the circuit is represented by:
            #      ┌───┐
            # q_0: ┤ X ├
            #      └───┘
            # q_1: ─────
            #      ┌───┐
            # q_2: ┤ X ├
            #      └───┘
            # q_3: ─────
            #
            # However, in Pennylane they are represented by alpha-beta
            # sequences. So, organize the term accordingly...
            n = len(term)//2  # => valid for closed shell systems only ?
            term = ''.join([term[i::n] for i in range(n)])
            # this could also be done by using the `_process_wires` function.

            if active_new_opmath():
                return qml.prod(
                    qml.pauli.string_to_pauli_word(term, wire_map=wire_map))

            return qml.pauli.string_to_pauli_word(term, wire_map=wire_map)

        return qml.Identity(wires[0])

    coeffs, ops = zip(
        *[(coef, _get_op(term, wires))
          for term, coef in qubit_operator.to_list()]
    )

    return np.real(np.array(coeffs, requires_grad=False)), list(ops)

This could also be implemented in qml.qchem.convert.import_operator as

def import_operator(qubit_observable, format="openfermion",
                    wires=None, tol=1e010):
    r"""Convert an external operator to a PennyLane operator.

    The external format currently supported is openfermion and qiskit.

    Args:
        qubit_observable: external qubit operator that will be converted
        format (str): the format of the operator object to convert from
        wires (.Wires, list, tuple, dict): Custom wire mapping used to convert
            the external qubit operator to a PennyLane operator.
            For types ``Wires``/list/tuple, each item in the iterable
            represents a wire label for the corresponding qubit index.
            For type dict, only int-keyed dictionaries (for qubit-to-wire
            conversion) are accepted. If ``None``, the identity map
            (0->0, 1->1, ...) will be used.
        tol (float): Tolerance in `machine epsilon
            <https://numpy.org/doc/stable/reference/generated/numpy.real_if_close.html>`
            for the imaginary part of the coefficients in ``qubit_observable``.
            Coefficients with imaginary part less than 2.22e-16*tol are
            considered to be real.

    Returns:
        (.Operator): PennyLane operator representing any operator expressed as
        linear comb of Pauli words, e.g., :math:`\\sum_{k=0}^{N-1} c_k O_k`

    **Example**

    >>> from openfermion import QubitOperator
    >>> h_of = QubitOperator('X0 X1 Y2 Y3', -0.0548)
    + QubitOperator('Z0 Z1', 0.14297)
    >>> h_pl = import_operator(h_of, format='openfermion')
    >>> print(h_pl)
    (0.14297) [Z0 Z1]
    + (-0.0548) [X0 X1 Y2 Y3]

    >>> from qiskit.quantum_info import SparsePauliOp
    >>> h_qt = SparsePauliOp.from_list([("XXYY", -0.0548), ("ZZII", 0.14297)])
    >>> h_pl = import_operator(h_qt, format='qiskit')
    >>> print(h_pl)
    (0.14297) [Z1 Z3]
    + (-0.0548) [Y0 X1 Y2 Y3]

    If the new op-math is active, an arithmetic operator is returned instead.

    >>> qml.operation.enable_new_opmath()
    >>> h_pl = import_operator(h_of, format='openfermion')
    >>> print(h_pl)
    (-0.0548*(PauliX(wires=[0]) @ PauliX(wires=[1]) @ PauliY(wires=[2]) @
    PauliY(wires=[3]))) + (0.14297*(PauliZ(wires=[0]) @ PauliZ(wires=[1])))
    """
    if format not in ["openfermion", "qiskit"]:
        raise TypeError(f"Converter does not exist for {format} format.")

    if format == "openfermion":
        # dealing with openfermion `QubitOperator`
        coeffs = np.array([np.real_if_close(coef, tol=tol)
                           for coef in qubit_observable.terms.values()])
    elif format == "qiskit":
        # dealing with qiskit `SparsePauliOp`
        coeffs = np.array([np.real_if_close(coef, tol=tol)
                           for coef in qubit_observable.coeffs])

    if any(np.iscomplex(coeffs)):
        warnings.warn(
            f"The coefficients entering the QubitOperator"
            f" or SparsePauliOp must be real;"
            f" got complex coefficients in the operator"
            f" {list(coeffs[np.iscomplex(coeffs)])}"
        )

    if format == "openfermion":
        if active_new_opmath():
            return qml.dot(*_openfermion_to_pennylane(
                qubit_observable, wires=wires))

        return qml.Hamiltonian(*_openfermion_to_pennylane(
            qubit_observable, wires=wires))

    if format == "qiskit":
        if active_new_opmath():
            return qml.dot(*_qiskit_nature_to_pennylane(
                qubit_observable, wires=wires))

        return qml.Hamiltonian(*_qiskit_nature_to_pennylane(
            qubit_observable, wires=wires))

If this idea is of interest, I could assign myself to this task and create a separate PR to implement these functions in addition to _pennylane_to_qiskit_nature and _qiskit_nature_pennylane_equivalent.

How important would you say this feature is?

2: Somewhat important. Needed this quarter.

Additional information

No response

[trbromley]

Hi @Cmurilochem! This sounds like a really nice idea and we're definitely open to contributors.

I just want to make sure I understand as I'm not so familiar with Qiskit Nature. Will your proposed conversion functionality support Qiskit SparsePauliOp instances as input or Qiskit Nature objects like operators or hamiltonians?

[Cmurilochem]

Hi @Cmurilochem! This sounds like a really nice idea and we're definitely open to contributors.

I just want to make sure I understand as I'm not so familiar with Qiskit Nature. Will your proposed conversion functionality support Qiskit SparsePauliOp instances as input or Qiskit Nature objects like operators or hamiltonians?

Hi @trbromley. Thanks for your prompt feedback!

Yes, this would only support Qiskit SparsePauliOp objects as input like it is does currently in the case of openfermion hamiltonians. So, the user would have to get SparsePauliOp forms from somewhere else, e.g., via a Qiskit-Nature driver or from native dataformats like QCSchema.

I will give you an example with some expectation values calculations:

# 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 entrance point to PennyLane

# get pennylane qubit `Operator` hamiltonian from `SparsePauliOp`
pl_qubit_ham = qml.qchem.import_operator(qtn_qubit_ham, format="qiskit-nature") # <====== Assuming that `_qiskit_nature_to_pennylane` has been implemented and the extra option "qiskit-nature" is added to `import_operator`

# 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 of which would be:

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

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) [I0]
+ (-0.22474449848985895) [Z2]
+ (-0.2247444984898589) [Z3]
+ (0.17184931866629144) [Z0]
+ (0.17184931866629144) [Z1]
+ (0.12078809196857414) [Z0 Z2]
+ (0.12078809196857414) [Z1 Z3]
+ (0.16605110981082186) [Z0 Z3]
+ (0.16605110981082186) [Z1 Z2]
+ (0.16882419223387474) [Z0 Z1]
+ (0.1745434714459974) [Z2 Z3]
+ (0.045263017842247726) [Y0 Y1 Y2 Y3]
+ (0.045263017842247726) [Y0 X1 Y2 X3]
+ (0.045263017842247726) [X0 Y1 X2 Y3]
+ (0.045263017842247726) [X0 X1 X2 X3]

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

If that is OK, I will create a PR to implement this function into qml.qchem.convert.

[trbromley]

Nice, thanks @Cmurilochem for moving quickly on opening a PR! We'll take a look over the coming days and give feedback.

[Cmurilochem]

Nice, thanks @Cmurilochem for moving quickly on opening a PR! We'll take a look over the coming days and give feedback.

Thanks @trbromley. I might ask your help later on for devising units tests; tried this morning but they look complicated at first glance. Looking forward to receiving your feedback on the PR.