A brief overview of VQE

Author: PennyLane dev team. Last updated: 8 Apr 2021.

The Variational Quantum Eigensolver (VQE) [1], [2] is a flagship algorithm for quantum chemistry using near-term quantum computers. VQE is an application of the Ritz variational principle where a quantum computer is used to prepare a wave function ansatz of the molecule and estimate the expectation value of its electronic Hamiltonian while a classical optimizer is used to adjust the quantum circuit parameters in order to find the molecule’s ground state energy.

For example, if we use a minimal basis, the ground state wave function of the hydrogen molecule \(\vert \Psi \rangle = \alpha \vert 1100 \rangle + \beta \vert 0011 \rangle\) consists of only the Hartree-Fock component and a doubly-excited configuration where the two electrons occupy the highest-energy molecular orbitals. If we use a quantum computer to prepare the four-qubit entangled state \(\vert \Psi \rangle\), the ultimate goal of the VQE algorithm is to find the values of \(\alpha\) and \(\beta\) that minimize the expectation value of the electronic Hamiltonian.

The PennyLane library allows users to implement the full VQE algorithm using only a few lines of code. In this tutorial, we guide you through a calculation of the ground-state energy of the hydrogen molecule. Let’s get started! ⚛️

Building the electronic Hamiltonian

The first step is to import the required libraries and packages:

import pennylane as qml
from pennylane import qchem
from pennylane import numpy as np

The second step is to specify the molecule whose properties we aim to calculate. This is done by providing three pieces of information: the geometry and charge of the molecule, and the spin multiplicity of the electronic configuration.

The geometry of a molecule is given by the three-dimensional coordinates and symbols of all its atomic species. There are several databases such as the NIST Chemistry WebBook, ChemSpider and SMART-SNS that provide geometrical data for a large number of molecules. Here, we make use of a locally saved file in .xyz format that contains the geometry of the hydrogen molecule, and specify its name for later use:

geometry = 'h2.xyz'

Alternatively, you can download the file here: h2.xyz.

The charge determines the number of electrons that have been added or removed compared to the neutral molecule. In this example, as is the case in many quantum chemistry simulations, we will consider a neutral molecule:

charge = 0

It is also important to define how the electrons occupy the molecular orbitals to be optimized within the Hartree-Fock approximation. This is captured by the multiplicity parameter, which is related to the number of unpaired electrons in the Hartree-Fock state. For the neutral hydrogen molecule, the multiplicity is one:

multiplicity = 1

Finally, we need to specify the basis set used to approximate atomic orbitals. This is typically achieved by using a linear combination of Gaussian functions. In this example, we will use the minimal basis STO-3g where a set of 3 Gaussian functions are contracted to represent an atomic Slater-type orbital (STO):

basis_set = 'sto-3g'

At this stage, to compute the molecule’s Hamiltonian in the Pauli basis, several calculations need to be performed. With PennyLane, these can all be done in a single line by calling the function molecular_hamiltonian(). The charge, multiplicity, and basis set can also be specified as keyword arguments. Finally, the number of active electrons and active orbitals may be indicated, as well as the fermionic-to-qubit mapping, which can be either Jordan-Wigner (jordan_wigner) or Bravyi-Kitaev (bravyi_kitaev). The atomic symbols and their nuclear coordinates can be read directly from the geometry file. The outputs of the function are the qubit Hamiltonian of the molecule and the number of qubits needed to represent it:

symbols, coordinates = qchem.read_structure(geometry)

h, qubits = qchem.molecular_hamiltonian(

print('Number of qubits = ', qubits)
print('Hamiltonian is ', h)


Number of qubits =  4
Hamiltonian is    (-0.24274280513140462) [Z2]
+ (-0.24274280513140462) [Z3]
+ (-0.04207897647782276) [I0]
+ (0.1777128746513994) [Z1]
+ (0.17771287465139946) [Z0]
+ (0.12293305056183798) [Z0 Z2]
+ (0.12293305056183798) [Z1 Z3]
+ (0.1676831945771896) [Z0 Z3]
+ (0.1676831945771896) [Z1 Z2]
+ (0.17059738328801052) [Z0 Z1]
+ (0.17627640804319591) [Z2 Z3]
+ (-0.04475014401535161) [Y0 Y1 X2 X3]
+ (-0.04475014401535161) [X0 X1 Y2 Y3]
+ (0.04475014401535161) [Y0 X1 X2 Y3]
+ (0.04475014401535161) [X0 Y1 Y2 X3]

That’s it! From here on, we can use PennyLane as usual, employing its entire stack of algorithms and optimizers.

Implementing the VQE algorithm

PennyLane contains the ExpvalCost class, specifically built to implement the VQE algorithm. We begin by defining the device, in this case a simple qubit simulator:

dev = qml.device('default.qubit', wires=qubits)

In VQE, the goal is to train a quantum circuit to prepare the ground state of the input Hamiltonian. This requires a clever choice of circuit, which should be complex enough to prepare the ground state, but also sufficiently easy to optimize. In this example, we employ a variational circuit that is capable of preparing the normalized states of the form \(\alpha|1100\rangle + \beta|0011\rangle\) which encode the ground state wave function of the hydrogen molecule described with a minimal basis set. The circuit consists of single-qubit rotations on all wires, followed by three entangling CNOT gates, as shown in the figure below:


In the circuit, we apply single-qubit rotations, followed by CNOT gates:

def circuit(params, wires):
    qml.BasisState(np.array([1, 1, 0, 0], requires_grad=False), wires=wires)
    for i in wires:
        qml.Rot(*params[i], wires=i)
    qml.CNOT(wires=[2, 3])
    qml.CNOT(wires=[2, 0])
    qml.CNOT(wires=[3, 1])


The qubit register has been initialized to \(|1100\rangle\) which encodes the Hartree-Fock state of the hydrogen molecule described with a minimal basis.

The cost function for optimizing the circuit can be created using the ExpvalCost class, which is tailored for VQE optimization. It requires specifying the circuit, target Hamiltonian, and the device, and returns a cost function that can be evaluated with the circuit parameters:

cost_fn = qml.ExpvalCost(circuit, h, dev)

Wrapping up, we fix an optimizer and randomly initialize circuit parameters. For reliable results, we fix the seed of the random number generator, since in practice it may be necessary to re-initialize the circuit several times before convergence occurs.

opt = qml.GradientDescentOptimizer(stepsize=0.4)
params = np.random.normal(0, np.pi, (qubits, 3))



[[ 5.54193389  1.25713095  3.07479606]
 [ 7.03997361  5.86710646 -3.07020901]
 [ 2.98479079 -0.47550269 -0.32427159]
 [ 1.28993324  0.45252622  4.56873497]]

We carry out the optimization over a maximum of 200 steps, aiming to reach a convergence tolerance (difference in cost function for subsequent optimization steps) of \(\sim 10^{ -6}\).

max_iterations = 200
conv_tol = 1e-06

for n in range(max_iterations):
    params, prev_energy = opt.step_and_cost(cost_fn, params)
    energy = cost_fn(params)
    conv = np.abs(energy - prev_energy)

    if n % 20 == 0:
        print('Iteration = {:},  Energy = {:.8f} Ha'.format(n, energy))

    if conv <= conv_tol:

print('Final convergence parameter = {:.8f} Ha'.format(conv))
print('Final value of the ground-state energy = {:.8f} Ha'.format(energy))
print('Accuracy with respect to the FCI energy: {:.8f} Ha ({:.8f} kcal/mol)'.format(
    np.abs(energy - (-1.136189454088)), np.abs(energy - (-1.136189454088))*627.503
print('Final circuit parameters = \n', params)


Iteration = 0,  Energy = -0.88179557 Ha
Iteration = 20,  Energy = -1.13380513 Ha
Iteration = 40,  Energy = -1.13558756 Ha
Iteration = 60,  Energy = -1.13585794 Ha
Iteration = 80,  Energy = -1.13600617 Ha
Iteration = 100,  Energy = -1.13608848 Ha
Iteration = 120,  Energy = -1.13613394 Ha

Final convergence parameter = 0.00000099 Ha
Final value of the ground-state energy = -1.13615709 Ha
Accuracy with respect to the FCI energy: 0.00003237 Ha (0.02031093 kcal/mol)

Final circuit parameters =
 [[ 5.54193389e+00  1.30219523e-08  3.07479606e+00]
 [ 7.03997361e+00  6.28318530e+00 -3.07020901e+00]
 [ 2.98479079e+00 -2.09540998e-01 -4.16893297e-02]
 [ 1.28993324e+00  1.30907669e-12  4.56873497e+00]]

Success! 🎉🎉🎉 The ground-state energy of the hydrogen molecule has been estimated with chemical accuracy (< 1 kcal/mol) with respect to the exact value of -1.136189454088 Hartree (Ha) obtained from a full configuration-interaction (FCI) calculation. This is because, for the optimized values of the single-qubit rotation angles, the state prepared by the VQE ansatz is precisely the FCI ground-state of the \(H_2\) molecule \(|H_2\rangle_{gs} = 0.99 |1100\rangle - 0.10 |0011\rangle\).

What other molecules would you like to study using PennyLane?


[1]Alberto Peruzzo, Jarrod McClean et al., “A variational eigenvalue solver on a photonic quantum processor”. Nature Communications 5, 4213 (2014).
[2]Yudong Cao, Jonathan Romero, et al., “Quantum Chemistry in the Age of Quantum Computing”. Chem. Rev. 2019, 119, 19, 10856-10915.

Total running time of the script: ( 0 minutes 29.805 seconds)

Gallery generated by Sphinx-Gallery