Quantinuum H-Series - Quantum Subspace Expansion

This tutorial illustrates the computation of molecular excited states through the application of the Quantum Subspace Expansion (QSE) technique on the Quantinuum H-series remote emulator. The QSE method leverages a precise calculation of the ground state energy for a specific molecular configuration to approximate the energies of corresponding excited states. Consider a stable molecule, denoted by a Hamiltonian \(H\), and let \(\left|\Psi_{0}\right\rangle\) represent the outcome of the Variational Quantum Eigensolver (VQE) algorithm employed to estimate the ground state energy of \(H\). More information on the QSE method can be found in the publication available at this link.

The QSE technique constructs a subspace of state vectors \(\left|\Psi_j^k\right\rangle\) formed by one-electron excitations of the ground state wavefunction:

\begin{equation} \left|\Psi_{j}^{k}\right\rangle = c_k^{\dagger}c_{j}\left|\Psi_0\right\rangle. \end{equation}

where \(c_k^{\dagger}, c_{j}\) are the fermionic creation and annihilation operators over spin orbitals \(k\) and \(j\), respectively. That is, these vectors are formed by reducing the occupation of spin orbital \(j\) by one, and increasing the occupation of spin orbital \(k\) by one. The vectors are not in general orthogonal to \(\Psi_{0}\) hence we will need to calculate an overlap matrix.

Within this subspace, we solve a generalized eigenvalue problem. Consider the operator $:nbsphinx-math:hat{H} $ with matrix elements given by

(109)\[(H)_{jk}^{lm} = \langle\Psi_j^l \left| \hat{H} \right| \Psi_k^m\rangle,\]

and define an overlap matrix \(S\) whose matrix elements are given by

(110)\[S_{jk}^{lm} = \langle \Psi_j^l \left|\Psi_k^m\right\rangle.\]

The generalized eigenvalue equation to be solved is

\begin{equation} HC=SCE, \end{equation}

where \(C\) is the matrix of eigenvectors, and \(E\) is the vector of eigenvalues. Crucially, the energy eigenvalues \(E\) provide an estimate of the excited state energies of \(H\) as well as a refined value of the ground state energy.

Notice that the solution to the generalized eigenvalue equation can be done on a classical computer, provided \(H\) and \(S\) have been calculated. The matrix elements of both of these matrices can be constructed using a quantum computer, in the following way. First, re-write the matrix elements in terms of \(\left | \Psi_{0}\right \rangle\):

\begin{align} (H)_{jk}^{lm} =& \langle\Psi_j^l \left| \hat{H} \right| \Psi_k^m\rangle = \langle \Psi_{0} | c_{j}^\dagger c_{l} \hat{H}c_{m}^{\dagger}c_{k}|\Psi_{0}\rangle\\ S_{jk}^{lm} &= \langle \Psi_j^l \left|\Psi_k^m\right\rangle = \langle \Psi_{0} | c_{j}^\dagger c_{l} c_{m}^{\dagger}c_{k}|\Psi_{0}\rangle. \end{align}

The matrix elements can be calculated using a quantum computer or simulator. By transforming the operators \(c_{j}^\dagger c_{l} c_{k}^{\dagger}c_{m}\) and \(c_{l}^\dagger c_{j} \hat{H} c_{m}^{\dagger}c_{k}\) to a set of Pauli quantum gates according to an appropriate scheme such as Jordan-Wigner or Bravyi-Kitaev, apply this gate set to the ground state wavefunction (constructed with the coefficients obtained from the VQE calculation) and perform a measurement to obtain the expected value. By measuring these expectation values, we obtain estimates of \(S_{jk}^{lm}\) and \(H_{jk}^{lm}\), respectively.

This tutorial will guide you through the process of performing QSE and will show how to collate the resulting eigen-energies and eigen-vectors. To demonstrate these concepts, we will utilize the methane (CH\(_4\)) chemical system. Note that to follow this tutorial, you must have been granted access to Quantinuum systems by your InQuanto administrator.

The steps below are such:

  • Define the system

  • Prepare the approximate ground state obtained by the VQE algorithm

  • Configure the Quantinuum backend

  • Establish the PMSV error mitigation

  • Define protocol for measuring QSE matrix elements

  • Build InQuanto computables

  • Submit and retreive InQuanto computables

  • The black-box approach using the AlgorithmQSE

After defining the Z-matrix for Methane, the Hamiltonian object is created using a restricted Hartree-Fock (RHF) InQuanto-PySCF driver.

The get_system function is responsible for computing the fermionic Hamiltonian operator, Fock space, and Hartree Fock state.

[ ]:

import warnings
warnings.filterwarnings('ignore')

from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHF

zmatrix = """
C
H             1    1.083000
H             1    1.083000      2  109.471000
H             1    1.083000      2  109.471000      3  120.000000
H             1    1.083000      2  109.471000      4  120.000000
"""

driver = ChemistryDriverPySCFMolecularRHF(
    zmatrix=zmatrix,
    charge=0,
    frozen=[0, 1, 2, 3, 7, 8],
    basis="STO-3G",
    point_group_symmetry=True,
)

hamiltonian, space, state = driver.get_system()

The qubit_encode function carries out qubit encoding, utilizing the mapping class associated with the current integral operator. The default mapping approach is the Jordan-Wigner method. Additionally, we employ the compress method, which eliminates Hamiltonian terms with coefficients less than abs_tol, chosen as \(10^{-6}\) here. This step is a means of reducing the number of measurement circuits/quantum computational resources needed (which means less time to run which is preferred for a tutorial), and the threshold we use is fairly reasonable compared to the scale of errors induced by noisy quantum simulation.

[ ]:

qubit_hamiltonian = hamiltonian.qubit_encode()
qubit_hamiltonian.compress(abs_tol=1e-6)
qubit_hamiltonian.df()
Coefficient Term Coefficient Type
0 -37.822198 <class 'numpy.float64'>
1 -0.496463 Z5 <class 'numpy.float64'>
2 -0.496463 Z4 <class 'numpy.float64'>
3 0.128963 Z4 Z5 <class 'numpy.float64'>
4 -0.496463 Z3 <class 'numpy.float64'>
5 0.110247 Z3 Z5 <class 'numpy.float64'>
6 0.123222 Z3 Z4 <class 'numpy.float64'>
7 -0.012975 X2 X3 Y4 Y5 <class 'numpy.float64'>
8 0.012975 X2 Y3 Y4 X5 <class 'numpy.float64'>
9 0.012975 Y2 X3 X4 Y5 <class 'numpy.float64'>
10 -0.012975 Y2 Y3 X4 X5 <class 'numpy.float64'>
11 -0.496463 Z2 <class 'numpy.float64'>
12 0.123222 Z2 Z5 <class 'numpy.float64'>
13 0.110247 Z2 Z4 <class 'numpy.float64'>
14 0.128963 Z2 Z3 <class 'numpy.float64'>
15 -0.067741 Z1 <class 'numpy.float64'>
16 0.106870 Z1 Z5 <class 'numpy.float64'>
17 0.118715 Z1 Z4 <class 'numpy.float64'>
18 0.106870 Z1 Z3 <class 'numpy.float64'>
19 0.118715 Z1 Z2 <class 'numpy.float64'>
20 -0.011845 X0 X1 Y4 Y5 <class 'numpy.float64'>
21 -0.011845 X0 X1 Y2 Y3 <class 'numpy.float64'>
22 0.011845 X0 Y1 Y4 X5 <class 'numpy.float64'>
23 0.011845 X0 Y1 Y2 X3 <class 'numpy.float64'>
24 0.011845 Y0 X1 X4 Y5 <class 'numpy.float64'>
25 0.011845 Y0 X1 X2 Y3 <class 'numpy.float64'>
26 -0.011845 Y0 Y1 X4 X5 <class 'numpy.float64'>
27 -0.011845 Y0 Y1 X2 X3 <class 'numpy.float64'>
28 -0.067741 Z0 <class 'numpy.float64'>
29 0.118715 Z0 Z5 <class 'numpy.float64'>
30 0.106870 Z0 Z4 <class 'numpy.float64'>
31 0.118715 Z0 Z3 <class 'numpy.float64'>
32 0.106870 Z0 Z2 <class 'numpy.float64'>
33 0.123178 Z0 Z1 <class 'numpy.float64'>

To construct our ansatz for the specified fermion space and fermion state, we have employed the Chemically Aware Unitary Coupled Cluster method with singles and doubles excitations (UCCSD). The circuit is synthesized using Jordan-Wigner encoding.

[ ]:

from inquanto.ansatzes import FermionSpaceAnsatzChemicallyAwareUCCSD

ansatz = FermionSpaceAnsatzChemicallyAwareUCCSD(space, state)

Next, a noiseless VQE simulation is executed on the QulacsBackend. The goal is to achieve a minimized expectation value and a set of circuit parameters corresponding to the ground state. These optimized parameters are subsequently employed in the Hamiltonian averaging experiment later on.

[ ]:

from inquanto.express import run_vqe
from pytket.extensions.qulacs import QulacsBackend

state_backend = QulacsBackend()
vqe = run_vqe(ansatz, hamiltonian, backend=state_backend, with_gradient=False)

parameters = vqe.final_parameters

# TIMER BLOCK-0 BEGINS AT 2024-02-19 10:02:38.479944
# TIMER BLOCK-0 ENDS - DURATION (s):  0.8142243 [0:00:00.814224]

To simulate the specific noise profiles of machines, we can load and apply them to our simulations using the QuantinuumBackend, which retrieves information from your Quantinuum account. The QuantinuumBackend offers a range of available emulators, such as H1-1E and H1-2E. These emulators are designed for specific devices and they run remotely on a server.

Parameters used:

device_name – Name of device, e.g. “H1-1”

label – Job labels used if Circuits have no name, defaults to “job”

group – string identifier of a collection of jobs, can be used for usage tracking.

The pytket-quantinuum extension allows the user to access the following quantum devices, emulators and syntax checkers.

[ ]:

from pytket.extensions.quantinuum import QuantinuumBackend

backend = QuantinuumBackend(device_name="H1-1E", group ="Default - UK")

To reduce errors and inaccuracies caused by quantum noise and imperfections in the Quantinuum device, we can employ noise mitigation techniques. In this case, we will define the Qubit Operator symmetries within the system, enabling us to utilize PMSV (Partition Measurement Symmetry Verification). More information about QubitOperator can be found in the provided link.

[ ]:

from inquanto.operators import QubitOperator
from inquanto.protocols import PMSV

stabilizers = [QubitOperator("Z0 Z1 Z2 Z3 Z4 Z5", 1)]

mitms_pmsv = PMSV(stabilizers)

To expand the subspace for computing low-lying excited states, it is essential to define the basis. In this context, the function generate_subspace_singlet_singles is employed, which sequentially provides spin-adapted singlet-single excitation operators to preserve spin symmetry which are then mapped to qubits using the Jordan-Wigner method.

[ ]:

from inquanto.mappings import QubitMappingJordanWigner

mapping = QubitMappingJordanWigner()

fermionic_excitations_singlets = list(space.generate_subspace_singlet_singles())
qubit_excitations_singlets = mapping.operator_map(fermionic_excitations_singlets)

Now that the Hamiltonian, Ansatz and QSE basis operators are defined, we can build the computable. The computable is the chemical quantity that we are interested in. In this example, we are interested in building the QSE computable using the QSEMatricesComputable function.

Parameters used:

state – Ansatz used to represent the ground state, it is used as a reference state for the excitations to generate the subspace.

hermitian_operator – A hermitian operator, typically a hamiltonian, to be expanded in the subspace.

expansion_operators – A List of excitation operators spanning the subspace.

[ ]:

from inquanto.computables.composite import QSEMatricesComputable

qse = QSEMatricesComputable(
    state=ansatz,
    hermitian_operator=qubit_hamiltonian,
    expansion_operators=qubit_excitations_singlets,
)

To compute the expectation value of a Hermitian operator through operator averaging on the system register, we employ the PauliAveraging protocol. This protocol effectively implements the procedure outlined in Operator Averaging.

[ ]:

from inquanto.protocols import PauliAveraging
from pytket.partition import PauliPartitionStrat

protocol = PauliAveraging(
    backend=backend,
    shots_per_circuit=6000,
    pauli_partition_strategy=PauliPartitionStrat.CommutingSets,
)

protocol.build_from(parameters=parameters, computable=qse, noise_mitigation=mitms_pmsv)

<inquanto.protocols.averaging._pauli_averaging.PauliAveraging at 0x16ac06990>

Despite having many Qubit Operators to measure, this protocol ends up generating 30 measurement circuits due to the commutation relations measurement partitioning strategy pytket uses.

[ ]:

protocol.n_circuit

30

To launch the circuits to the backend we have used the launch function. This method processes all the circuits associated with the expectation value calculations and returns a list of ResultHandle objects representing the handles for the results. We can pickle these ResultHandle objects so we can retrieve the results once they are ready. After our experiments have finished, we can obtain the results by utilizing the retrieve function, which retrieves distributions from the backend for the specified source.

[ ]:

handles = protocol.launch()

[ ]:

import pickle

with open("handles.pickle", "wb") as handle:
    pickle.dump(handles, handle, protocol=pickle.HIGHEST_PROTOCOL)

[ ]:

completion_check=backend.circuit_status(handles[-1])[0].value #n-1
print(completion_check)

Circuit is queued.

[ ]:

protocol.retrieve(handles)

<inquanto.protocols.averaging._pauli_averaging.PauliAveraging at 0x7fa080330990>

Once the results have been retrieved the Hamiltonian expectation value is evaluated as a classical post-processing step.

We now need to pass the output of the QSEMatricesComputable into pd_safe_eigh. This will solve the generalized eigenvalue equation and return the eigenvalues and eigenvectors.

[ ]:

H,S = qse.evaluate(protocol.get_evaluator())

[ ]:

from inquanto.core import pd_safe_eigh

e_vals_singlets, e_coeffs_singlets, _ = pd_safe_eigh(H,S)

print(e_vals_singlets)
print(e_coeffs_singlets)

[-39.729 -38.91  -38.871 -38.768 -38.273 -37.903]
[[-0.156+0.476j  0.004-0.008j  0.002-0.006j -0.018+0.023j  0.002-0.015j
  -0.   -0.002j]
 [-0.311-0.135j -0.24 +0.368j -0.567+0.08j   3.806+1.333j  0.888+0.622j
  -0.223-0.734j]
 [-0.388-0.069j -0.606-0.105j -0.781+0.008j  4.621+0.967j -1.663-0.86j
  -0.85 +0.603j]
 [-0.003-0.003j  0.576+0.164j -0.384+0.022j  0.06 +0.043j -0.061+0.026j
  -0.004-0.055j]
 [-0.077-0.165j -1.04 +0.671j -0.19 -0.518j  0.644+1.722j -1.634+2.515j
  -5.204+1.559j]
 [ 0.16 +0.204j  0.08 +1.131j  0.391+0.141j -2.197-2.655j  0.879+3.241j
  -2.243+2.55j ]
 [-0.006+0.003j  0.111+0.338j  0.363+0.46j   0.113+0.117j -0.12 -0.056j
  -0.007+0.04j ]
 [-0.098-0.106j -0.638+0.085j -0.395+0.104j  1.81 +1.216j -1.329+1.704j
   1.878-1.188j]
 [ 0.17 +0.018j  0.057+0.8j    0.203+0.411j -1.287-0.708j  1.636+4.126j
   3.24 -1.954j]]

We perform additional post-processing to report the operator expansion for each excited state.

[ ]:

from inquanto.operators import FermionOperator

e_vecs = []
for i in range(e_coeffs_singlets.shape[1]):
    operator = FermionOperator()
    for c, o in zip(e_coeffs_singlets[:, i], fermionic_excitations_singlets):
        operator += c * o
    operator.compress(abs_tol=1e-5)
    print(operator)
    e_vecs.append(operator)

(-0.15565783977198872+0.47621487975804033j, F0^ F0 ), (-0.15565783977198872+0.47621487975804033j, F1^ F1 ), (-0.3109942131939058-0.13521605338785778j, F0^ F2 ), (-0.3109942131939058-0.13521605338785778j, F1^ F3 ), (-0.38839307178519544-0.0693114600174391j, F0^ F4 ), (-0.38839307178519544-0.0693114600174391j, F1^ F5 ), (-0.003308852334800071-0.002955219215334559j, F2^ F0 ), (-0.003308852334800071-0.002955219215334559j, F3^ F1 ), (-0.07692730809458342-0.1649068223822534j, F2^ F2 ), (-0.07692730809458342-0.1649068223822534j, F3^ F3 ), (0.16028553390102782+0.2040034257419802j, F2^ F4 ), (0.16028553390102782+0.2040034257419802j, F3^ F5 ), (-0.006135282730581221+0.003137920667681536j, F4^ F0 ), (-0.006135282730581221+0.003137920667681536j, F5^ F1 ), (-0.09846811137359278-0.10620905658730515j, F4^ F2 ), (-0.09846811137359278-0.10620905658730515j, F5^ F3 ), (0.16993533787082993+0.01769895101511043j, F4^ F4 ), (0.16993533787082993+0.01769895101511043j, F5^ F5 )
(0.003529209585032674-0.00795687921061023j, F0^ F0 ), (0.003529209585032674-0.00795687921061023j, F1^ F1 ), (-0.24049475157327982+0.36783111046288186j, F0^ F2 ), (-0.24049475157327982+0.36783111046288186j, F1^ F3 ), (-0.6060191423196757-0.10482632989727117j, F0^ F4 ), (-0.6060191423196757-0.10482632989727117j, F1^ F5 ), (0.5763459636232429+0.16388842490954422j, F2^ F0 ), (0.5763459636232429+0.16388842490954422j, F3^ F1 ), (-1.0398193192980176+0.6711364842936922j, F2^ F2 ), (-1.0398193192980176+0.6711364842936922j, F3^ F3 ), (0.0799938237057527+1.1312076955513257j, F2^ F4 ), (0.0799938237057527+1.1312076955513257j, F3^ F5 ), (0.11138656902183633+0.3380893725102996j, F4^ F0 ), (0.11138656902183633+0.3380893725102996j, F5^ F1 ), (-0.6378073443599778+0.08473450806257327j, F4^ F2 ), (-0.6378073443599778+0.08473450806257327j, F5^ F3 ), (0.057037317243711416+0.7995538108062838j, F4^ F4 ), (0.057037317243711416+0.7995538108062838j, F5^ F5 )
(0.002165953897558483-0.005944635760931628j, F0^ F0 ), (0.002165953897558483-0.005944635760931628j, F1^ F1 ), (-0.5667818952627056+0.08033281708487372j, F0^ F2 ), (-0.5667818952627056+0.08033281708487372j, F1^ F3 ), (-0.7806938566797863+0.008375225559272712j, F0^ F4 ), (-0.7806938566797863+0.008375225559272712j, F1^ F5 ), (-0.38358059039903575+0.02164733189057025j, F2^ F0 ), (-0.38358059039903575+0.02164733189057025j, F3^ F1 ), (-0.18951924122622935-0.5175643220754953j, F2^ F2 ), (-0.18951924122622935-0.5175643220754953j, F3^ F3 ), (0.39067878123749633+0.14105366840747063j, F2^ F4 ), (0.39067878123749633+0.14105366840747063j, F3^ F5 ), (0.3633037576320181+0.45991596424530234j, F4^ F0 ), (0.3633037576320181+0.45991596424530234j, F5^ F1 ), (-0.3948958613818514+0.10360791672209882j, F4^ F2 ), (-0.3948958613818514+0.10360791672209882j, F5^ F3 ), (0.2025074399411661+0.41092059234737977j, F4^ F4 ), (0.2025074399411661+0.41092059234737977j, F5^ F5 )
(-0.017840047182897188+0.022744956740089296j, F0^ F0 ), (-0.017840047182897188+0.022744956740089296j, F1^ F1 ), (3.806126291783216+1.332658212611577j, F0^ F2 ), (3.806126291783216+1.332658212611577j, F1^ F3 ), (4.621131902658994+0.9671995436795001j, F0^ F4 ), (4.621131902658994+0.9671995436795001j, F1^ F5 ), (0.05993966785070938+0.04310122860320857j, F2^ F0 ), (0.05993966785070938+0.04310122860320857j, F3^ F1 ), (0.6441688080331356+1.7220627585285808j, F2^ F2 ), (0.6441688080331356+1.7220627585285808j, F3^ F3 ), (-2.1965755386805927-2.655483790470051j, F2^ F4 ), (-2.1965755386805927-2.655483790470051j, F3^ F5 ), (0.1127899167896759+0.11659167840761259j, F4^ F0 ), (0.1127899167896759+0.11659167840761259j, F5^ F1 ), (1.809508151433246+1.2160491220244634j, F4^ F2 ), (1.809508151433246+1.2160491220244634j, F5^ F3 ), (-1.2874605040396814-0.7077883650219732j, F4^ F4 ), (-1.2874605040396814-0.7077883650219732j, F5^ F5 )
(0.0021294584990190677-0.01537836611697584j, F0^ F0 ), (0.0021294584990190677-0.01537836611697584j, F1^ F1 ), (0.8881021827827966+0.6217291428756166j, F0^ F2 ), (0.8881021827827966+0.6217291428756166j, F1^ F3 ), (-1.6634091829941275-0.85957877147374j, F0^ F4 ), (-1.6634091829941275-0.85957877147374j, F1^ F5 ), (-0.06072464916897717+0.025946744324903372j, F2^ F0 ), (-0.06072464916897717+0.025946744324903372j, F3^ F1 ), (-1.6336871970371127+2.515156393240569j, F2^ F2 ), (-1.6336871970371127+2.515156393240569j, F3^ F3 ), (0.8787624912303986+3.2406507977554675j, F2^ F4 ), (0.8787624912303986+3.2406507977554675j, F3^ F5 ), (-0.12020930798769339-0.055773317529686006j, F4^ F0 ), (-0.12020930798769339-0.055773317529686006j, F5^ F1 ), (-1.3286485705307103+1.7035098710284213j, F4^ F2 ), (-1.3286485705307103+1.7035098710284213j, F5^ F3 ), (1.6355067063551598+4.126046236307394j, F4^ F4 ), (1.6355067063551598+4.126046236307394j, F5^ F5 )
(-0.00020756753326981142-0.0016959811491862975j, F0^ F0 ), (-0.00020756753326981142-0.0016959811491862975j, F1^ F1 ), (-0.22283094986115085-0.7342082125387235j, F0^ F2 ), (-0.22283094986115085-0.7342082125387235j, F1^ F3 ), (-0.8502303263381792+0.6025252787670932j, F0^ F4 ), (-0.8502303263381792+0.6025252787670932j, F1^ F5 ), (-0.004028844274401662-0.054943191515560284j, F2^ F0 ), (-0.004028844274401662-0.054943191515560284j, F3^ F1 ), (-5.203599573643404+1.5592992892443864j, F2^ F2 ), (-5.203599573643404+1.5592992892443864j, F3^ F3 ), (-2.2427076382898403+2.550269786999596j, F2^ F4 ), (-2.2427076382898403+2.550269786999596j, F3^ F5 ), (-0.006964622239363893+0.040135082617393666j, F4^ F0 ), (-0.006964622239363893+0.040135082617393666j, F5^ F1 ), (1.878027835548134-1.18788450746501j, F4^ F2 ), (1.878027835548134-1.18788450746501j, F5^ F3 ), (3.2403778890630868-1.953712459189611j, F4^ F4 ), (3.2403778890630868-1.953712459189611j, F5^ F5 )

Alternatively, the Quantum Subspace Expansion algorithm implemented in InQuanto, known as AlgorithmQSE, can be employed. This is a black-box method enabling direct retrieval of eigenvalues and eigenvectors, stored as final_value and final_states, respectively. The algorithm addresses numerical instabilities during matrix diagonalization by eliminating near-linear dependencies in the basis. However, it’s worth noting that the user has less control over the specific details of the process. Some examples that cannot be executed include launching, retrieving, and performing PMSV.

[ ]:

from inquanto.algorithms import AlgorithmQSE

algorithm = AlgorithmQSE(
    computable_qse_matrices=qse,
    parameters=parameters,
)

protocol_expression = PauliAveraging(
    backend=backend,
    shots_per_circuit=1000,
    pauli_partition_strategy=PauliPartitionStrat.CommutingSets,
)

algorithm.build(protocol=protocol_expression)

#algorithm.run()

#print(algorithm.generate_report())

<inquanto.algorithms.qse._algorithm_qse.AlgorithmQSE at 0x7fa0534b7850>