VQE for a spin model using qse Operators

This notebook builds a spin model using the qse Operator / Operators classes, converts it to a Qiskit SparsePauliOp via to_qiskit() method, and runs VQE to find the ground state energy.

import numpy as np
import qse
from qse.operator.operator import Operator
from qse.operator.operators import Operators

1. Build the spin model.

We build a 1D transverse-field Ising model (TFIM) on \(N=10\) qubits:

\[H = -J \sum_{i=0}^{N-2} Z_i Z_{i+1} - h \sum_{i=0}^{N-1} X_i\]

N = 10
J = 1.0
h = 0.5
spacing = 1.0

qbits = qse.lattices.chain(spacing, N)

# ZZ interaction terms between neighbouring qubits
hamiltonian = qbits.compute_interaction_hamiltonian(
    distance_func=lambda d: -J * np.isclose(d, spacing),
    interaction="Z",
)

# Transverse field X terms on each qubit
for i in range(N):
    hamiltonian += qse.Operator("X", i, N, coef=-h)

2. Convert to Qiskit SparsePauliOp

Once we have the Hamiltonian \(H\), we can use the to_qiskit() method to transform it to a Qiskit Sparse Pauli Operator.

pauli_op = hamiltonian.to_qiskit()
print(pauli_op)

SparsePauliOp(['ZZIIIIIIII', 'IZZIIIIIII', 'IIZZIIIIII', 'IIIZZIIIII', 'IIIIZZIIII', 'IIIIIZZIII', 'IIIIIIZZII', 'IIIIIIIZZI', 'IIIIIIIIZZ', 'XIIIIIIIII', 'IXIIIIIIII', 'IIXIIIIIII', 'IIIXIIIIII', 'IIIIXIIIII', 'IIIIIXIIII', 'IIIIIIXIII', 'IIIIIIIXII', 'IIIIIIIIXI', 'IIIIIIIIIX'],
              coeffs=[-1. +0.j, -1. +0.j, -1. +0.j, -1. +0.j, -1. +0.j, -1. +0.j, -1. +0.j,
 -1. +0.j, -1. +0.j, -0.5+0.j, -0.5+0.j, -0.5+0.j, -0.5+0.j, -0.5+0.j,
 -0.5+0.j, -0.5+0.j, -0.5+0.j, -0.5+0.j, -0.5+0.j])

3. Set up and run VQE.

We set up and run the VQE in qiskit.

from qiskit.primitives import StatevectorEstimator
from qiskit.circuit.library import efficient_su2

import numpy as np
from scipy.optimize import minimize

ansatz = efficient_su2(num_qubits=N, reps=2, entanglement="linear")
estimator = StatevectorEstimator()
energies_list = []


def cost_function(params):
    result = estimator.run([(ansatz, pauli_op, params)]).result()
    energy = result[0].data.evs
    energies_list.append(float(energy))
    return float(energy)


# Random initial parameters
rng = np.random.default_rng(42)
x0 = rng.uniform(-np.pi, np.pi, ansatz.num_parameters)

result = minimize(
    cost_function, x0, method="COBYLA", options={"maxiter": 2000, "rhobeg": 0.5}
)

ground_state = np.linalg.eigh(pauli_op.to_matrix())[0][0]
print(f"VQE ground state energy : {result.fun:.6f}")
print(f"Exact ground state energy: {ground_state}")

VQE ground state energy : -9.746544
Exact ground state energy: -9.765503957927187