LRZ QPU Tequila Tutorial

code
Author

Marlene Hermelink

Published

June 18, 2025

This tutorial shows how to initialize the MQP backend in Tequila for the LRZ QPU. It also provides a use case example for computing the groundstate energy of a HCB Hamiltonian using the MQP backend.

Initializing the MQP backend

The Munich Quantum Valley provides a Qiskit-backend, which can be used as a sampling backend for Tequila. If this backend is installed and configured with a valid API-Key, it is possible to run Tequila Quantum Circuits on the LRZ Quantum Computer.

Follow these steps:

  1. install the mqp qiskit provider pip install mqp-qiskit-provider
  2. successful installation can be checked with
import tequila as tq
tq.show_available_simulators()
  1. initialize the MQP backend as follows:
from mqp.qiskit_provider import MQPProvider

provider = MQPProvider('API_KEY')  # replace 'API_KEY' with your actual API key
[device] = provider.backends('AQT20')

device should now contain the initialized MQPBackend.

  1. This device can be passed to the tequila simulate function:
tq.simulate(circuit, backend="mqp", device=device, samples=200)

Note that the samples parameter is required, as the LRZ backend does not support statevector simulation and will throw an error if it’s omitted. The LRZ QPU currently does not support more than 200 samples in one job.

If the job was successfully submitted, it will show up on the Munich Quantum Portal, though it might take some time until it is scheduled and executed.

Local Simulation with AQT Simulator

Because of the long wait times and limited resources for jobs on the LRZ Quantum Computer, it is useful to test circuits locally before sending them to the LRZ QPU. The Qiskit AQT Provider offers a backend that represents a noiseless simulation of the AQT Hardware, which can be used for this purpose.

The AQT Backend is initialized like this:

# Select an execution backend.
# Any token (even invalid) gives access to the offline simulation backends.
provider = AQTProvider("INVALID_TOKEN")
device = provider.get_backend("offline_simulator_no_noise")

Then pass device to the tq.simulate function as before.

Use Case Example: HCB Groupings

The following demonstrates how to use Tequilas MQP Backend to compute the ground state energy of a Hardcore Boson Hamiltonian (HCB).

Background: VQE

The variational quantum eigensolver (VQE) is a technique for finding the smallest eigenvalue of the electronic Hamiltonian, which describes the groundstate energy of a molecule (Peruzzo). The algorithm uses the variational principle

\[argmin_{\theta} \langle U(\theta)|H|U(\theta)\rangle\]

where \(\theta\) somehow parametrizes the quantum circuit \(U(\theta)\) and \(H\) is the Hamiltonian. The expectation value of \(H\) is an upper bound for the smallest eigenvalue \(E_0\). To see why this works: because \(H\) is Hermitian, any state \(|\psi\rangle\) can be written as a linear combination of orthogonal eigenstates of the Hamiltonian:

\[|\psi\rangle = \sum_k c_k|E_k\rangle\]

so

\[\langle\psi|H|\psi\rangle = \sum_k |c_k|^2\langle E_k|H|E_k\rangle = \sum_k |c_k|^2 E_k \geq E_0\]

By finding the \(\theta\) that minimizes \(\langle U(\theta)|H|U(\theta)\rangle\) we can find an approximation for \(E_0\). However, we still need to find a suitable ansatz for \(U(\theta)\) and suitable approximation for \(H\).

The HCB Hamiltonian

The electronic hamiltonian of a molecule in second quantization is defined as

\[H = \sum_{kl} h_{kl} a^\dagger_k a_l + \frac{1}{2}\sum_{klmn} g_{klmn}a_k^\dagger a_l^\dagger a_n a_m\]

where \(a_k^\dagger\) and \(a_k\) are the creation and annihilation operators.

In order to evaluate the expectation value of this hamiltonian on a quantum computer, it has to be expressed as a sum of Paulistrings. Typically the Jordan-Wigner encoding is used to achieve this. But there is a problem: in Jordan-Wigner the creation and annihilation operators are encoded as follows:

\(a_k = (\otimes^{k-1}_{l=1} Z)\otimes \sigma^+\): annihilation operator

\(a_k^\dagger = (\otimes^{k-1}_{l=1} Z)\otimes \sigma^-\): creation operator

with \(\sigma^{±} = \frac{1}{2}(X±iY)\)

to encode the antisymmetry of the wavefunction inside the creation and annihilation operators. However, this encoding is rather nonlocal (i.e. there are many qubits affected by a single action), as there is a \(Z\) placed on every qubit before \(k\) (Santos and Kottmann).

The HCB-Hamiltonian is an approximation that improves the locality and reduces the number of gates in the circuit by ‘spin-pairing electrons into quasi-bosonic particles’ (Bincoletto and Kottmann) and restricting the number of unpaired electrons to zero. In this model, electrons can only be created and annihilated in pairs (so only double excitations are possible and every spatial orbital is either filled with two or zero electrons). This makes these operations symmetric. It also has the additional benefit of cutting the number of necessary qubits in half because now there is one qubit per spatial orbital instead of one qubit per spin-orbital.

But there is another advantage.

Usually, when sampling the Expectation Value of a Hamiltonian \(H = \sum_i c_iP_i\) that consists of multiple Pauli-Strings \(P_i\) with Tequila, every Paulistring is diagonalized into a Pauli-\(Z\)-String with a basis changing unitary operation, and then the expectation values for each diagonalized Paulistring are measured individually. This is necessary because measurements on the QPU are done by reading out the qubits, which means that we are measuring in the \(Z\)-basis.

As a consequence the number of jobs that are submitted to the backend is equal to the number of Paulistrings in the Hamiltonian.

However, because we are not the only users of the LRZ QPU, it might take a long time until our jobs are actually scheduled, so we would like to reduce the number of jobs. We want to find circuits that diagonalize as many Paulistrings simultaneously as possible.

Luckily, the HCB-Hamiltonian naturally decomposes into three commuting groups of Paulistrings: an \(X\), a \(Y\) and a \(Z\) group, depending on whether the Pauli-String contains only \(X\), \(Y\) or \(Z\) operators (Bincoletto and Kottmann). With this grouping, we can reduce the number of jobs to three.

The SPA Ansatz

It is important to choose an ansatz for \(U(\theta)\) that is capable of expressing the groundstate adequately. Separable Pair Approximation is one possible ansatz. It ‘leverages valence-bond structure’ (Bincoletto and Kottmann) – i.e. we derive the SPA circuit from the chemical graph of the molecule (where we have already made some assumptions about the bonds).

Implementation in Tequila

HCB Grouping

The following method make_hcb_grouping(H) performs the grouping for a HCB-Hamiltonian and simultaneously computes the unitary operators that are needed to diagonalize each group. (Bincoletto and Kottmann)

The three HCB-Groups are then sampled by appending their diagonalizing circuits to the original circuit, before sending it to the backend.

# H needs to be a HCB Hamiltonian
def make_hcb_grouping(H):
    H1 = tq.QubitHamiltonian()
    H2 = tq.QubitHamiltonian()
    H3 = tq.QubitHamiltonian()
 
    # z group is already in the right basis
    U1 = tq.QCircuit()
    # we can diagonalize X like this: H X H = Z 
    U2 = tq.gates.H([i for i in H.qubits])
    # we can diagonalize Y like this: Rx(pi/2) Y Rx(-pi/2) = Z
    U3 = tq.gates.Rx(angle=-np.pi/2, target=[i for i in H.qubits])
    for p in H.paulistrings:
        q = p.naked().qubits
        # hcb z group
        if p.is_all_z():
          H1 += tq.QubitHamiltonian().from_paulistrings(p)
        else:
            # hcb x group
            if (p.naked()[q[0]] == "X"):
                for k, v in p.items():
                    p._data[k] = "Z"
                H2 += tq.QubitHamiltonian().from_paulistrings(p)
            # hcb y group
            else:
              for k, v in p.items():
                  p._data[k] = "Z"
              H3 += tq.QubitHamiltonian().from_paulistrings(p)
    
    hamiltonians = [H1, H2, H3]
    circuits = [U1, U2, U3]
    result = [(H, U) for H, U, G in zip(hamiltonians, circuits)]

    return result

The following steps instruct how to create the SPA ansatz and the HCB Hamiltonian and how to use the HCB grouping:

  1. SPA ansatz and HCB Hamiltonian
mol = tq.Molecule(geometry="h 0.0 0.0 0.0\nh 0.0 0.0 1.0", basis_set="sto-3g", transformation="JordanWigner")
    
mol = mol.use_native_orbitals()

# initial guess for the chemical graph - SPA single bond between the two H atoms with bonding and antibonding MO
guess = np.array([[1.0, 1.0], [-1.0, 1.0]])
    
# HCB-SPA is necessary here because 'SPA' alone would encode one spin-orbital per qubit which wouldn't match the HCB hamiltonian
U_HCB  = mol.make_ansatz(name="HCB-SPA", edges=[(0, 1)])    # creating the SPA ansatz
opt = tq.chemistry.optimize_orbitals(mol, circuit=U_HCB, initial_guess=guess.T,silent=True, use_hcb=True)
H_HCB = opt.molecule.make_hardcore_boson_hamiltonian()      # creating the HCB Hamiltonian
  1. Then make the HCB grouping with the method from before:
hcbs_groups = make_hcb_grouping(H_HCB)
  1. Finally calculate the expectation value for each group and add the resulting energies:
E = tq.ExpectationValue(H=H_HCB, U=U_HCB, optimize_measurements=False)

result = tq.minimize(E, silent=True)      # find the parameter of the SPA circuit that minimizes E 
exact_energy = result.energy
v = result.variables
   
groups = ["Z", "X", "Y"]

result = 0
total_exact_energy = 0
for j, (H, U) in enumerate(hcbs_groups):

    Ehcb = tq.ExpectationValue(H=H, U=U_HCB + U)
    result_sampl_g = tq.simulate(Ehcb, variables=v, backend=backend, device=device, samples=200)
   `
    exact_result = tq.minimize(Ehcb, silent=True)
    exact_energy = exact_result.energy

    result += result_sampl_g
    total_exact_energy += exact_energy

The above code will send only three jobs to the backend.

Symmetry-Based Error Correction

After measuring the expectation values for each group on the LRZ QPU for different molecules, we noticed that the variance and error of the calculated groundstate energy were very high for the \(Z\)-group.

Figure 1 shows a histogram of counts for the \(Z\)-HCB group of the \(H_4\) HCB-Hamiltonian. The different bitstring counts for the ground state calculation that was closest to the exact energy and the bitstring counts for the energy that were the most inaccurate are shown.

Figure 1: Counts for Z-HCB Group

The \(Z\)-Group was the most error-prone. Here, the bitstring that appears most often in the best result was sampled not even once in the worst result. For comparision, here are the histrograms for the \(X\) and \(Y\) groups:

Counts for X-HCB Group Counts for Y-HCB Group

Luckily, the \(Z\)-Group is already in the \(Z\)-basis. This allows us to apply a symmetry-based error mitigation: As described above, in the HCB Hamiltonian, one qubit encodes one spatial orbital and if a \(1\) is measured in a qubit, that means that there are \(2\) electrons in the corresponding spatial orbital. So the number of 1s that can appear in a single bitstring must be equal to \(\frac{a}{2}\) where \(a\) is the number of non-core electrons in the molecule, i.e. for \(H_4\) we expect to measure two 1s.

If a bitstring contains more 1s, we know there must have been a severe error in the computation and we discard this sample. E.g. in the worst result for the \(Z\)-group, the most frequent bistring was ‘1111’, which is physically impossible and should be discarded. The following code applies this error mitigation to a dictionary of counts and recalculates the ground state energy. (Kottmann et al.)

import tequila as tq
from typing import Dict,Union


def count_ones(binstring: str) -> int:
    sum_ones = 0
    for char in binstring:
        if char not in '01':
            raise ValueError(f"Invalid character '{char}' in binary string '{binstring}'")
        if char == '1':
            sum_ones += 1
    return sum_ones


# discard bitstrings with wrong number of 1s
def correct_data(number_hs, counts : Dict) -> Dict:
    elecron_pairs = number_hs // 2
    for key, val in counts.items():
        if count_ones(key) != electron_pairs:
            counts[key] = 0
    return counts

# method copied from tequilas sample_all_z_hamiltonian
def recalculate_energy(hamiltonian, counts: Dict):
        # all qubits in the hamiltonian are also in the circuit
        abstract_qubits_H = hamiltonian.qubits
        assert len(abstract_qubits_H) != 0  
        # assert that the Hamiltonian was mapped before
        read_out_map = {q: i for i, q in enumerate(abstract_qubits_H)}
        bitstring_counts = {}
        for k, v in counts.items():
            bitstring_key = tq.BitString.from_binary(k[::-1])
            bitstring_counts[bitstring_key] = v

        # compute energy
        E = 0.0
        for paulistring in hamiltonian.paulistrings:
            n_samples = 0
            Etmp = 0.0
            if paulistring._data.keys() == 0:
                Etmp += paulistring.coeff
                E += Etmp 
                continue
            for key, count in bitstring_counts.items():
                # gives us the qubits the paulistring matrices are on
                mapped_ps_support = [read_out_map[i] for i in paulistring._data.keys()]
                # count all measurements that resulted in |1> for those qubits
                parity = [k for i, k in enumerate(key.array) if i in mapped_ps_support].count(1)
                # evaluate the PauliString
                sign = (-1) ** parity
                Etmp += sign * count
                n_samples += count
            E += (Etmp / n_samples) * paulistring.coeff
        return E


# hamiltonian string needs to be all z
def correct_energy(number_hs:int, hamiltonian:str, counts:Dict)-> float:
    hamiltonian = tq.QubitHamiltonian.from_string(hamiltonian)
    corrected_counts = correct_data(number_hs, counts)
    corrected_energy = recalculate_energy(hamiltonian, corrected_counts)
    return corrected_energy


 if __name__ == "__main__":
    # example usage for H2
    counts = {'01': 190, '00': 1, '10': 5, '11': 4}
    hamiltonian = "-0.0077+0.2743Z(0)-0.2607Z(1)+0.5233Z(0)Z(1)"
    corrected_energy = correct_energy(2, hamiltonian, counts)
    print(f"Corrected energy: {corrected_energy}")

Figure 2 compares the maximal error of the corrected and the uncorrected energies.

Figure 2: Maximal Energy Error for corrected and uncorrected bitstrings

Clearly the symmetry-based error mitigation significantly reduced the error. However, for larger molecule the error remained several Hartree high.

References

Bincoletto, Davide, and Jakob S. Kottmann. State Specific Measurement Protocols for the Variational Quantum Eigensolver. 2025. https://arxiv.org/abs/2504.03019.
Kottmann, Jakob S., Sumner Alperin-Lea, Teresa Tamayo-Mendoza, et al. Tequila: A platform for rapid development of quantum algorithms. Version 1.8.5. Released November 2020. https://github.com/tequilahub/tequila.
Peruzzo, McClean, A. A Variational Eigenvalue Solver on a Photonic Quantum Processor. 2014. https://doi.org/10.1038/ncomms5213.
Santos, Francisco Javier Del Arco, and Jakob S. Kottmann. A Hybrid Qubit Encoding: Splitting Fock Space into Fermionic and Bosonic Subspaces. 2024. https://arxiv.org/abs/2411.14096.