Note

Click here to download the full example code

# VQE in different spin sectors¶

*Author: Alain Delgado — Posted: 13 October 2020. Last updated: 25 June 2022.*

The Variational Quantum Eigensolver (VQE) algorithm is an approach for finding the lowest-energy state of a molecule using a quantum computer 1.

In the absence of spin-orbit coupling, the eigenstates of the molecular Hamiltonian can be calculated for specific values of the spin quantum numbers. This allows us to compute quantities such as the energy of the electronic states in different sectors of the total-spin projection \(S_z\). This is illustrated in the figure below for the energy spectrum of the hydrogen molecule. In this case, the ground state has total spin \(S=0\) while the lowest-lying excited states, with total spin \(S=1\), form a triplet related to the spin components \(S_z=-1, 0, 1\).

In this tutorial we demonstrate how to run VQE simulations to find the lowest-energy states of the hydrogen molecule in different spin sectors. First, we show how to build the electronic Hamiltonian and the total spin operator \(\hat{S}^2\). Next, we use excitation operations, implemented in PennyLane as Givens rotations 2, to prepare the trial states of the molecule. In order to probe the molecular states with different total spin, \(S=0\) and \(S=1\), we apply excitation operations which preserve and modify, respectively, the total-spin projection of the initial state encoded in the qubit register. Finally, we run the VQE algorithm to compute the energy of the states.

Let’s get started!

## Building the Hamiltonian and the total spin operator \(\hat{S}^2\)¶

First, we need to specify the structure of the molecule. This is done by providing a list with the symbols of the constituent atoms and a one-dimensional array with the corresponding nuclear coordinates in atomic units.

```
from pennylane import numpy as np
symbols = ["H", "H"]
coordinates = np.array([0.0, 0.0, -0.6614, 0.0, 0.0, 0.6614])
```

The geometry of the molecule can also be imported from an external file using
the `read_structure()`

function.

Now, we can build the electronic Hamiltonian. We use a minimal basis set approximation to represent
the molecular orbitals. Then,
the qubit Hamiltonian of the molecule is built using the
`molecular_hamiltonian()`

function.

```
import pennylane as qml
H, qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)
print("Number of qubits = ", qubits)
print("The Hamiltonian is ", H)
```

Out:

```
Number of qubits = 4
The Hamiltonian is (-0.2427450126094144) [Z2]
+ (-0.2427450126094144) [Z3]
+ (-0.042072551947439224) [I0]
+ (0.1777135822909176) [Z0]
+ (0.1777135822909176) [Z1]
+ (0.12293330449299361) [Z0 Z2]
+ (0.12293330449299361) [Z1 Z3]
+ (0.16768338855601356) [Z0 Z3]
+ (0.16768338855601356) [Z1 Z2]
+ (0.17059759276836803) [Z0 Z1]
+ (0.1762766139418181) [Z2 Z3]
+ (-0.044750084063019925) [Y0 Y1 X2 X3]
+ (-0.044750084063019925) [X0 X1 Y2 Y3]
+ (0.044750084063019925) [Y0 X1 X2 Y3]
+ (0.044750084063019925) [X0 Y1 Y2 X3]
```

The outputs of the function are the Hamiltonian and the number of qubits
required for the quantum simulations. For the \(\mathrm{H}_2\) molecule in a minimal
basis, we have four molecular **spin**-orbitals, which defines the number of qubits.

The `molecular_hamiltonian()`

function allows us to define
additional keyword arguments to simulate more complicated molecules. For more details
take a look at the tutorial Building molecular Hamiltonians.

We also want to build the total spin operator \(\hat{S}^2\),

In the equation above, \(N\) is the number of electrons, \(\hat{c}\) and \(\hat{c}^\dagger\) are respectively the electron annihilation and creation operators, and \(\langle \alpha, \beta \vert \hat{s}_1 \cdot \hat{s}_2 \vert \gamma, \delta \rangle\) is the matrix element of the two-body spin operator \(\hat{s}_1 \cdot \hat{s}_2\) in the basis of spin orbitals.

We use the `spin2()`

function to build the
\(\hat{S}^2\) observable.

```
electrons = 2
S2 = qml.qchem.spin2(electrons, qubits)
print(S2)
```

Out:

```
((0.375+0j)) [Z0]
+ ((0.375+0j)) [Z1]
+ ((0.375+0j)) [Z2]
+ ((0.375+0j)) [Z3]
+ ((0.75+0j)) [I0]
+ ((-0.375+0j)) [Z0 Z1]
+ ((-0.375+0j)) [Z2 Z3]
+ ((-0.125+0j)) [Z0 Z3]
+ ((-0.125+0j)) [Z1 Z2]
+ ((0.125+0j)) [Z0 Z2]
+ ((0.125+0j)) [Z1 Z3]
+ ((-0.125+0j)) [Y0 X1 X2 Y3]
+ ((-0.125+0j)) [X0 Y1 Y2 X3]
+ ((0.125+0j)) [Y0 X1 Y2 X3]
+ ((0.125+0j)) [Y0 Y1 X2 X3]
+ ((0.125+0j)) [Y0 Y1 Y2 Y3]
+ ((0.125+0j)) [X0 X1 X2 X3]
+ ((0.125+0j)) [X0 X1 Y2 Y3]
+ ((0.125+0j)) [X0 Y1 X2 Y3]
```

## Building the quantum circuit to find the ground state¶

We build the variational circuit for the ground state of the hydrogen molecule using the
`SingleExcitation`

and `DoubleExcitation`

operations
which act as Givens rotations on subspaces of two and four qubits, respectively
2. The total number of excitation operations is given by
all possible single and double excitations of the Hartree-Fock reference state. The indices
of the qubits they act on correspond to the spin-orbitals involved in each excitation.
For more details on how to use the excitation operations see the
tutorial Givens rotations for quantum chemistry.

First, we use the `hf_state()`

function to generate the vector representing the Hartree-Fock state
\(\vert 1100 \rangle\) of the \(\mathrm{H}_2\) molecule.

```
hf = qml.qchem.hf_state(electrons, qubits)
print(hf)
```

Out:

```
[1 1 0 0]
```

Next, we use the `excitations()`

function to generate all single and double excitations of the Hartree-Fock state.
This function allows us to define the keyword argument `delta_sz`

to specify the total-spin projection of the excitations with respect to the reference
state. This is illustrated in the figure below.

Therefore, for the ground state of the \(\mathrm{H}_2\) molecule we choose
`delta_sz = 0`

.

```
singles, doubles = qml.qchem.excitations(electrons, qubits, delta_sz=0)
print(singles)
print(doubles)
```

Out:

```
[[0, 2], [1, 3]]
[[0, 1, 2, 3]]
```

The output lists `singles`

and `doubles`

contain the qubit indices involved in the
single and double excitations. Even and odd indices correspond, respectively, to spin-up
and spin-down orbitals. For `delta_sz = 0`

we have two single excitations, one from qubit
0 to 2 and the other from qubit 1 to 3, and one double excitation from qubits 0, 1 to 2, 3.
In order to build the variational circuit, we apply the corresponding
excitation operations. This can be done straightforwardly using the
`AllSinglesDoubles`

template.

```
def circuit(params, wires):
qml.AllSinglesDoubles(params, wires, hf, singles, doubles)
```

The circuit above prepares trial states for the \(\mathrm{H}_2\) molecule of the form:

where the coefficients \(c\) are functions of the variational parameters \(\theta = (\theta_1, \theta_2, \theta_3)\) to be optimized by the VQE algorithm. Note that the prepared state \(\vert \Psi(\theta) \rangle\) is a superposition of the Hartree-Fock (HF) state with all possible single and double excitations that preserve the spin projection of the HF reference state.

## Running the VQE simulation¶

Now we proceed to optimize the variational parameters. First, we define the device, in this case a qubit simulator:

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

Next, we define the cost function as the following QNode, where we make use of
the `expval()`

function to compute the expectation value of the hamiltonian.
This requires specifying the circuit, the target Hamiltonian, and the device. It returns
a cost function that can be evaluated with the circuit parameters:

```
@qml.qnode(dev, interface="autograd")
def cost_fn(params):
circuit(params, wires=range(qubits))
return qml.expval(H)
```

As a reminder, we also built the total spin operator \(\hat{S}^2\) for which we can now define a function to compute its expectation value:

```
@qml.qnode(dev, interface="autograd")
def S2_exp_value(params):
circuit(params, wires=range(qubits))
return qml.expval(S2)
```

The total spin \(S\) of the trial state can be obtained from the expectation value \(\langle \hat{S}^2 \rangle\) as

We define a function to compute the total spin.

```
def total_spin(params):
return -0.5 + np.sqrt(1 / 4 + S2_exp_value(params))
```

Now, we proceed to minimize the cost function to find the ground state. We define the classical optimizer and initialize the circuit parameters.

```
opt = qml.GradientDescentOptimizer(stepsize=0.8)
np.random.seed(0) # for reproducibility
theta = np.random.normal(0, np.pi, len(singles) + len(doubles), requires_grad=True)
print(theta)
```

Out:

```
[5.54193389 1.25713095 3.07479606]
```

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

```
max_iterations = 100
conv_tol = 1e-06
for n in range(max_iterations):
theta, prev_energy = opt.step_and_cost(cost_fn, theta)
energy = cost_fn(theta)
spin = total_spin(theta)
conv = np.abs(energy - prev_energy)
if n % 4 == 0:
print(f"Step = {n}, Energy = {energy:.8f} Ha, S = {spin:.4f}")
if conv <= conv_tol:
break
print("\n" f"Final value of the ground-state energy = {energy:.8f} Ha")
print("\n" f"Optimal value of the circuit parameters = {theta}")
```

Out:

```
Step = 0, Energy = -0.09929557 Ha, S = 0.1014
Step = 4, Energy = -0.87153518 Ha, S = 0.0982
Step = 8, Energy = -1.11692841 Ha, S = 0.0087
Step = 12, Energy = -1.13529755 Ha, S = 0.0004
Step = 16, Energy = -1.13614888 Ha, S = 0.0000
Step = 20, Energy = -1.13618734 Ha, S = 0.0000
Final value of the ground-state energy = -1.13618832 Ha
Optimal value of the circuit parameters = [3.14350662 3.14087516 2.93185886]
```

As a result, we have estimated the lowest-energy state of the hydrogen molecule with total spin \(S = 0\) which corresponds to the ground state.

## Finding the lowest-lying excited state with \(S=1\)¶

In the last part of the tutorial, we will use VQE to find the lowest-lying
excited state of the hydrogen molecule with total spin \(S=1\).
In this case, we use the `excitations()`

function to generate
excitations whose total-spin projection differs by the quantity `delta_sz=1`

with respect to the Hartree-Fock state.

```
singles, doubles = qml.qchem.excitations(electrons, qubits, delta_sz=1)
print(singles)
print(doubles)
```

Out:

```
[[1, 2]]
[]
```

For the \(\mathrm{H}_2\) molecule in a minimal basis set there are no
double excitations, but only a spin-flip single excitation from qubit 1 to 2.
In this case, the circuit will contain only one `SingleExcitation`

operation. Additionally, as we want to probe the excited state of the hydrogen molecule,
we initialize the qubit register to the state \(\vert 0011 \rangle\).

```
def circuit(params, wires):
qml.AllSinglesDoubles(params, wires, np.flip(hf), singles, doubles)
```

This allows us to prepare trial states of the form

where the first term \(\vert 0101 \rangle\) encodes a spin-flip excitation with \(S_z=-1\) and the second term is a double excitation with \(S_z=0\). Since an eigenstate of the electronic Hamiltonian cannot contain a superposition of states with different total-spin projections, the double excitation coefficient should vanish as the VQE algorithm minimizes the cost function. The optimized state will correspond to the lowest-energy state with spin quantum numbers \(S=1, S_z=-1\).

Now, we define the new functions to compute the expectation values of the Hamiltonian and the total spin operator for the new circuit.

```
@qml.qnode(dev, interface="autograd")
def cost_fn(params):
circuit(params, wires=range(qubits))
return qml.expval(H)
@qml.qnode(dev, interface="autograd")
def S2_exp_value(params):
circuit(params, wires=range(qubits))
return qml.expval(S2)
```

Finally, we generate the new set of initial parameters, and proceed with the VQE algorithm to optimize the variational circuit.

```
np.random.seed(0)
theta = np.random.normal(0, np.pi, len(singles) + len(doubles), requires_grad=True)
max_iterations = 100
conv_tol = 1e-06
for n in range(max_iterations):
theta, prev_energy = opt.step_and_cost(cost_fn, theta)
energy = cost_fn(theta)
spin = total_spin(theta)
conv = np.abs(energy - prev_energy)
if n % 4 == 0:
print(f"Step = {n}, Energy = {energy:.8f} Ha, S = {spin:.4f}")
if conv <= conv_tol:
break
print("\n" f"Final value of the energy = {energy:.8f} Ha")
print("\n" f"Optimal value of the circuit parameters = {theta}")
```

Out:

```
Step = 0, Energy = 0.31463319 Ha, S = 0.3539
Step = 4, Energy = -0.38517129 Ha, S = 0.9391
Step = 8, Energy = -0.47698618 Ha, S = 0.9991
Step = 12, Energy = -0.47842743 Ha, S = 1.0000
Step = 16, Energy = -0.47844667 Ha, S = 1.0000
Final value of the energy = -0.47844667 Ha
Optimal value of the circuit parameters = [3.14259046]
```

As expected, the VQE algorithms has found the lowest-energy state with total spin \(S=1\) which is an excited state of the hydrogen molecule.

## Conclusion¶

In this tutorial we have used the standard VQE algorithm to find the ground and the lowest-lying excited states of the hydrogen molecule. We used the single and double excitation operations, implemented as Givens rotations, to prepare the trial states of a molecule. We also showed how to build the total-spin operator \(\hat{S}^2\) and used it to compute the total spin of the optimized states. By choosing the total-spin projection of the generated excitations, we were able to probe the lowest-energy eigenstates of the molecular Hamiltonian in different sectors of the spin quantum numbers.

## References¶

- 1
A. Peruzzo, J. McClean

*et al.*, “A variational eigenvalue solver on a photonic quantum processor”. Nature Communications 5, 4213 (2014).- 2(1,2)
J.M. Arrazola, O. Di Matteo, N. Quesada, S. Jahangiri, A. Delgado, N. Killoran. “Universal quantum circuits for quantum chemistry”. arXiv:2106.13839, (2021)