This page was generated from docs/tutorials/vqe_with_estimator.ipynb.

Variational Quantum Eigensolver with Estimator Primitive


The Variational Quantum Eigensolver (VQE) is an optimization routine for finding the ground state energy (i.e. lowest eigenvalue) of a Hamiltonian and is a considered to be a viable candidate for NISQ hardware. In this tutorial, we will go over how to use Qiskit Runtime to submit variational jobs using the estimator. Specifically, we will be looking at calculating the ground state energy of the \(H_2\) molecule using the estimator primitive.

Set up your local development environment

This tutorial requires a Qiskit Runtime service instance. If you haven’t done so already, follow these steps to set one up.

Connect to the Qiskit Runtime service

First we have to create our service instance and specify our backend. In this example we will be working with a simulator for the sake of speed, but by simply changing this one line (for choice of backend) this could run on hardware as well.

from qiskit_ibm_runtime import QiskitRuntimeService, Estimator, Session

service = QiskitRuntimeService()

backend = "ibmq_qasm_simulator"

Generate Molecular Hamiltonians

We will be working with Qiskit Nature to generate and handle molecular Hamiltonians. First though, we need to install these packages. This can be done through pip, using the commands shown below. For more information on setting up and working with Qiskit Nature, see this guide.

pip install qiskit-nature
pip install 'qiskit-nature[pyscf]'

Now we need to generate the Hamiltonians that we wish to find the ground state energy of. For this task we will be utilizing Qiskit Nature. First we have to specify how we are converting the fermionic operators of the electronic molecular Hamiltonian to qubit operators. Here we will use the Parity Mapper, which maps annihilation to Pauli operators via \(\hat{a} \rightarrow \frac{1}{2} \left ( X_p Z_{p-1} + i Y_p \right ) X_{p+1} ... X_{N}\)

from qiskit_nature.converters.second_quantization import QubitConverter
from qiskit_nature.mappers.second_quantization import ParityMapper

qubit_converter = QubitConverter(ParityMapper(), two_qubit_reduction=True)

At our specified bond length we will create an electronic structure problem and generate the second quantized (i.e. fermionic) operators. We then compute the true ground state energy to compare against.

from qiskit_nature.drivers import Molecule
from qiskit_nature.drivers.second_quantization import ElectronicStructureDriverType, ElectronicStructureMoleculeDriver
from qiskit_nature.problems.second_quantization import ElectronicStructureProblem
from qiskit.algorithms import NumPyMinimumEigensolver
from qiskit.algorithms import MinimumEigensolverResult
from qiskit.opflow import TaperedPauliSumOp
import numpy as np

dist = 0.72

ops = []

# Generate the H2 molecule
molecule = Molecule(geometry=[['H', [0., 0., 0.]], ['H', [0., 0., dist]]])
driver = ElectronicStructureMoleculeDriver(molecule, basis='sto3g', \

# Create the optimization problem to solve
es_problem = ElectronicStructureProblem(driver)
second_q_ops = es_problem.second_q_ops()
hamiltonian = qubit_converter.convert(second_q_ops[0], num_particles=es_problem.num_particles)
# Generate the real solution
sol = MinimumEigensolverResult()
sol.eigenvalue = NumPyMinimumEigensolver().compute_minimum_eigenvalue(hamiltonian).eigenvalue
real_solution = es_problem.interpret(sol).total_energies[0]

# Convert the hamiltonian to an opflow op
h = hamiltonian.primitive
h.coeffs = np.real(hamiltonian.coeffs)
ops.append(TaperedPauliSumOp(h, hamiltonian.z2_symmetries))

Now we can use this Hamiltonian as the observable in our estimator factory. The VQE routine is formalized as \(\min_\theta \langle \Psi | \hat{H} | \Psi \rangle\), so we just need to minimize the expectation values of the Hamiltonian.

from qiskit.algorithms.optimizers import NELDER_MEAD
from qiskit.circuit.library import RealAmplitudes

circuit = RealAmplitudes(num_qubits=2, reps=2)
convergence = []

with Session(service=service, backend=backend):
    estimator = Estimator()

    def evaluate_expectation(x):
        x = list(x)
        results = estimator.run(circuits=circuit, observables=ops, parameter_values=[x]).result().values[0]
        return np.real(results)

    def callback(fx):

    initial_point = np.random.uniform(-np.pi, np.pi, len(circuit.parameters))
    optimizer = NELDER_MEAD(80, callback=callback)
    result = optimizer.minimize(evaluate_expectation, x0=initial_point)

    sol = MinimumEigensolverResult()
    sol.eigenvalue = result.fun
    vqe_sol = es_problem.interpret(sol).total_energies[0]

Now we’ve solved the electronic Hamiltonian, we have to add the nuclear energies.

vqe_interpret = []
for i in range(len(convergence)):
    sol = MinimumEigensolverResult()
    sol.eigenvalue = convergence[i]
    sol = es_problem.interpret(sol).total_energies[0]

Create the convergence plot

import matplotlib.pyplot as plt

plt.rcParams["font.size"] = 14

# plot loss and reference value
plt.figure(figsize=(12, 6))
plt.plot(vqe_interpret, label="Estimator VQE")
plt.axhline(y=real_solution, color="tab:red", ls="--", label="Target")

plt.ylabel("Energy [H]")
plt.title("VQE energy")
import qiskit_ibm_runtime
from qiskit.tools.jupyter import *


Version Information

Qiskit SoftwareVersion
System information
Python version3.9.6
Python compilerClang 10.0.0
Python builddefault, Aug 18 2021 12:38:10
Memory (Gb)32.0
Tue Jul 19 12:56:46 2022 MDT

This code is a part of Qiskit

© Copyright IBM 2017, 2022.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.