Note

Click here to download the full example code

# VQE in different spin sectors¶

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

Quantum computers offer a promising avenue to perform first-principles simulations of the electronic structure of molecules and materials that are currently intractable using classical high-performance computers. In particular, the Variational Quantum Eigensolver (VQE) algorithm [1], [2] has proven to be a valuable quantum-classical approach to find the lowest-energy eigenstate of the electronic Hamiltonian by using Noisy Intermediate-Scale Quantum (NISQ) devices [3].

In the absence of spin-orbit coupling, the electronic Hamiltonian matrix is block diagonal in the total spin projection quantum number \(S_z\). This allows us to compute the energy spectrum of the Hamiltonian in a given sector of \(S_z\). For example, the figure below shows the energy spectra of the hydrogen molecule calculated in different sectors of the total-spin projection. The ground state with energy \(E_\mathrm{gs}=-1.136189\) Ha has total spin \(S=0\) and spin projection \(S_z=0\). The lowest-lying excited states, with energy \(E^*=-0.478453\) Ha and total spin \(S=1\), show a three-fold spin degeneracy related to the values of \(S_z=-1, 0, 1\).

Similarly, in the framework of VQE, if the quantum computer can be programmed to prepare many-body states in a specific sector of the total-spin projection \(S_z\), the variational optimization algorithm will allow us to estimate the energy of the lowest-energy state in this spin sector. More specifically, if we run a VQE simulation for the \(\mathrm{H}_2\) molecule in the subspace of states with \(S_z=0\), we will find the ground-state energy of the molecule. If the VQE simulation is performed in the subspace with \(S_z=1\), the optimized state will be in practice an excited state of the molecule, as shown in the figure above.

At the core of the VQE algorithm is the variational quantum circuit that is optimized to prepare the desired quantum states. The choice of the circuit is crucial for the success of the algorithm. The unitary coupled cluster ansatz [4] is a powerful quantum circuit that is believed to outperform the classical coupled cluster method [5], traditionally referred to as the gold standard of quantum chemistry.

We demonstrate how different functionalities implemented in PennyLane can be put together to run VQE simulations to find the lowest-energy states of the \(\mathrm{H}_2\) molecule in different sectors of the total spin \(S\). We also specify how to use the unitary coupled cluster ansatz, restricted to single and double excitations (UCCSD), as the variational circuit for the algorithm. These functionalities can be combined to estimate the energies of the ground and the lowest-lying excited states of the hydrogen molecule.

Let’s get started! ⚛️

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

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

```
import pennylane as qml
from pennylane import numpy as np
from pennylane import qchem
from pennylane.templates.subroutines import UCCSD
from functools import partial
```

The second step is to specify the molecule whose properties we aim to calculate. This is done by providing the name and the geometry of the molecule.

```
name = "h2"
geometry = "h2.xyz"
```

The geometry of the molecule can be given in any format recognized by Open Babel. Here, we used a locally saved file in xyz format specifying the three-dimensional coordinates and symbols of the atomic species.

In this example, we use a minimal basis set to model the hydrogen molecule.
In this approximation, the qubit Hamiltonian of the molecule in the Jordan-Wigner
representation is built using the `molecular_hamiltonian()`

function.

Out:

```
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]
```

The `molecular_hamiltonian()`

function allows us to define
an additional set of keyword arguments to provide the user with ample flexibility
to generate the Hamiltonian of more complicated systems. For more details take a look
at the tutorial Quantum Chemistry with PennyLane.

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

In the equation above, \(N_e\) is the number of active 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-particle spin operator
\(\hat{s}_1 \cdot \hat{s}_2\) in the basis of Hartree-Fock *spin* orbitals
[6]. The \(\mathrm{H}_2\) molecule has two electrons that populate,
within the minimal basis set approximation, four *spin* orbitals. As a reminder, the
variable `qubits`

output by the `molecular_hamiltonian()`

above stores the number of spin orbitals included the basis.

In order to build the spin operator \(\hat{S}^2\) we call the
`spin2()`

function.

Out:

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

The `spin2()`

function uses
`_spin2_matrix_elements()`

and
`observable()`

to compute the
matrix elements in the equation above and to build the many-body observable,
respectively.

Note

The `observable()`

function can be used to build any
many-body observable as long as we have access to the matrix elements of the
one- and/or two-particle operators in the basis of single-particle states.

## Implementing VQE with the UCCSD ansatz¶

PennyLane contains the `ExpvalCost`

class to implement the VQE algorithm.
We begin by defining the device, in this case a qubit simulator:

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

The next step is to define the variational quantum circuit used to prepare the state that minimizes the expectation value of the electronic Hamiltonian. In this example, we use the unitary coupled cluster ansatz truncated at the level of single and double excitations (UCCSD) [4].

The UCCSD method is a generalization of the traditional CCSD formalism used in quantum chemistry to perform post-Hartree-Fock electron correlation calculations. Within the first-order Trotter approximation [7], the UCCSD ground state of the molecule is built via the exponential ansatz [8],

In the latter equation, the indices \(r, s\) and \(p, q\) run respectively over the occupied and unoccupied molecular orbitals. The operator \(\hat{c}_p^\dagger \hat{c}_r\) creates a single excitation [5] since it annihilates an electron in the occupied orbital \(r\) and creates it in the unoccupied orbital \(p\). Similarly, the operator \(\hat{c}_p^\dagger \hat{c}_q^\dagger \hat{c}_r \hat{c}_s\) creates a double excitation.

The quantum circuits to exponentiate the excitation operators in the
Jordan-Wigner representation [8] are implemented by the
functions `SingleExcitationUnitary()`

and
`DoubleExcitationUnitary()`

contained
in the PennyLane templates library. Finally, the parameters \(\theta_{pr}\) and
\(\theta_{pqrs}\) have to be optimized to minimize the expectation value,

where \(\vert \mathrm{HF} \rangle\) is the Hartree-Fock state. The total number of
excitations determines the number of parameters \(\theta\) to be optimized by the
VQE algorithm as well as the depth of the UCCSD quantum circuit. For more details,
check the documentation of the
`SingleExcitationUnitary()`

and
`DoubleExcitationUnitary()`

functions.

Now, we demonstrate how to use PennyLane functionalities to build up the UCCSD
ansatz for VQE simulations. First, we use the `excitations()`

function to generate the whole set of single- and double-excitations for \(N_e\)
`electrons`

populating `qubits`

spin orbitals. Furthermore, we can define the selection rules
\(s_{z_p} - s_{z_r} = \Delta s_z\) and
\(s_{z_p} + s_{z_q} - s_{z_r} - s_{z_s}= \Delta s_z\) for the spin-projection of the
molecular orbitals involved in the single and double excitations
using the keyword argument `delta_sz`

. This allows us to prepare a
correlated state whose total-spin projection \(S_z\) is different from the one of
the Hartree-Fock state by the quantity `delta_sz`

. Therefore, we choose `delta_sz = 0`

to
prepare the ground state of the \(\mathrm{H}_2\) molecule.

```
singles, doubles = 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 indices representing the single and
double excitations. For the hydrogen molecule in a minimal basis set we have two single
and one double excitations. The latter means that in preparing the UCCSD ansatz the
`SingleExcitationUnitary()`

function has to be called
twice to exponentiate the two single-excitation operators while the
`DoubleExcitationUnitary()`

function is invoked once to
exponentiate the double excitation.

We use the function
`excitations_to_wires()`

to generate the set of wires that the UCCSD
circuit will act on. The inputs to this function are the indices stored in the
`singles`

and `doubles`

lists.

```
s_wires, d_wires = qchem.excitations_to_wires(singles, doubles)
print(s_wires)
print(d_wires)
```

Out:

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

Next, we need to define the reference state that the UCCSD unitary acts on, which is
just the Hartree-Fock state. This can be done straightforwardly with the
`hf_state()`

function. The output of this function is an array containing the occupation-number vector
representing the Hartree-Fock state.

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

Out:

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

Finally, we can use the `UCCSD()`

function to define
our VQE ansatz,

```
ansatz = partial(UCCSD, init_state=hf_state, s_wires=s_wires, d_wires=d_wires)
```

Next, we use the PennyLane class `ExpvalCost`

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

```
cost_fn = qml.ExpvalCost(ansatz, H, dev)
```

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:

```
S2_exp_value = qml.ExpvalCost(ansatz, S2, dev)
```

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))
```

Wrapping up, we fix an optimizer and randomly initialize the circuit parameters.

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

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 \(\sim 10^{-6}\). Furthermore, we track the value of the total spin \(S\) of the prepared state as it is optimized through the iterative procedure.

```
max_iterations = 100
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)
spin = total_spin(params)
if n % 4 == 0:
print("Iteration = {:}, E = {:.8f} Ha, S = {:.4f}".format(n, energy, spin))
if conv <= conv_tol:
break
print()
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.1361894507)), np.abs(energy - (-1.1361894507)) * 627.509474
)
)
```

Out:

```
Iteration = 0, E = -0.89301846 Ha, S = 0.4428
Iteration = 4, E = -1.03708792 Ha, S = 0.2303
Iteration = 8, E = -1.09698238 Ha, S = 0.1053
Iteration = 12, E = -1.12030623 Ha, S = 0.0458
Iteration = 16, E = -1.12965584 Ha, S = 0.0194
Iteration = 20, E = -1.13348679 Ha, S = 0.0081
Iteration = 24, E = -1.13506944 Ha, S = 0.0034
Iteration = 28, E = -1.13572502 Ha, S = 0.0014
Iteration = 32, E = -1.13599683 Ha, S = 0.0006
Iteration = 36, E = -1.13610956 Ha, S = 0.0002
Iteration = 40, E = -1.13615631 Ha, S = 0.0001
Iteration = 44, E = -1.13617571 Ha, S = 0.0000
Iteration = 48, E = -1.13618375 Ha, S = 0.0000
Final convergence parameter = 0.00000090 Ha
Final value of the ground-state energy = -1.13618578 Ha
Accuracy with respect to the FCI energy: 0.00000367 Ha (0.00230257 kcal/mol)
```

Success! 🎉🎉🎉 We have estimated the lowest-energy state with total-spin projection \(S_z=0\) for the hydrogen molecule, which is the ground state, with chemical accuracy. Notice also that the optimized UCCSD state is an eigenstate of the total spin operator \(\hat{S}^2\) with eigenvalue \(S=0\).

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

In the last part of the tutorial we want to demonstrate that VQE can also be used to find the lowest-energy excited states with total spin \(S=1\) and \(S_z \neq 0\). For the hydrogen molecule, this is the case for the states with energy \(E = -0.4784529844\) Ha and \(S_z=1\) and \(S_z=-1\).

Let’s consider the case of \(S_z=-1\) for which we can use the
`excitations()`

function with the keyword argument `delta_sz=1`

.

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

Out:

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

Notice that for the hydrogen molecule in a minimal basis set there are no double excitations but only a single excitation corresponding to the spin-flip transition between orbitals 1 and 2. And, that’s it!. From this point on the algorithm is the same as described above.

```
s_wires, d_wires = qchem.excitations_to_wires(singles, doubles)
hf_state = qchem.hf_state(electrons, qubits)
ansatz = partial(UCCSD, init_state=hf_state, s_wires=s_wires, d_wires=d_wires)
cost_fn = qml.ExpvalCost(ansatz, H, dev)
S2_exp_value = qml.ExpvalCost(ansatz, S2, dev)
```

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

```
np.random.seed(0)
params = np.random.normal(0, np.pi, len(singles) + len(doubles))
max_iterations = 100
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)
spin = total_spin(params)
if n % 4 == 0:
print("Iteration = {:}, E = {:.8f} Ha, S = {:.4f}".format(n, energy, spin))
if conv <= conv_tol:
break
print()
print("Final convergence parameter = {:.8f} Ha".format(conv))
print("Energy of the lowest-lying excited state = {:.8f} Ha".format(energy))
print(
"Accuracy with respect to the FCI energy: {:.8f} Ha ({:.8f} kcal/mol)".format(
np.abs(energy - (-0.4784529849)), np.abs(energy - (-0.4784529849)) * 627.509474
)
)
```

Out:

```
Iteration = 0, E = 0.37442032 Ha, S = 0.2839
Iteration = 4, E = 0.01440973 Ha, S = 0.6423
Iteration = 8, E = -0.33716938 Ha, S = 0.9068
Iteration = 12, E = -0.45310078 Ha, S = 0.9837
Iteration = 16, E = -0.47444690 Ha, S = 0.9974
Iteration = 20, E = -0.47783396 Ha, S = 0.9996
Iteration = 24, E = -0.47835772 Ha, S = 0.9999
Iteration = 28, E = -0.47843838 Ha, S = 1.0000
Iteration = 32, E = -0.47845080 Ha, S = 1.0000
Final convergence parameter = 0.00000084 Ha
Energy of the lowest-lying excited state = -0.47845164 Ha
Accuracy with respect to the FCI energy: 0.00000134 Ha (0.00084308 kcal/mol)
```

As expected, we have successfully estimated the lowest-energy state with total spin \(S=1\) and \(S_z=-1\) for the hydrogen molecule, which is an excited state.

Now, you can run a VQE simulation to find the degenerate excited state with spin quantum numbers \(S=1\) and \(S_z=1\). Give it a try!

## References¶

[1] | A. Peruzzo, J. McClean et al., “A variational eigenvalue solver on a photonic
quantum processor”. Nature Communications 5, 4213 (2014). |

[2] | Y. Cao, J. Romero, et al., “Quantum chemistry in the age of quantum computing”.
Chem. Rev. 2019, 119, 19, 10856-10915. |

[3] | A. Kandala, A. Mezzacapo et al., “Hardware-efficient variational quantum
eigensolver for small molecules and quantum magnets”. arXiv:1704.05018 |

[4] | (1, 2) J. Romero, R. Babbush, et al., “Strategies for quantum computing molecular
energies using the unitary coupled cluster ansatz”. arXiv:1701.02691 |

[5] | (1, 2) F. Jensen. “Introduction to computational chemistry”.
(John Wiley & Sons, 2016). |

[6] | A. Fetter, J. D. Walecka, “Quantum theory of many-particle systems”. Courier Corporation, 2012. |

[7] | M. Suzuki. “Decomposition formulas of exponential operators and Lie exponentials with some applications to quantum mechanics and statistical physics”. Journal of Mathematical Physics 26, 601 (1985). |

[8] | (1, 2) P. Kl. Barkoutsos, J. F. Gonthier, et al., “Quantum algorithms for electronic structure
calculations: particle/hole Hamiltonian and optimized wavefunction expansions”.
arXiv:1805.04340. |

**Total running time of the script:** ( 2 minutes 14.439 seconds)

## Contents

## Downloads

## Related tutorials