VQE in different spin sectors

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\).


../_images/energy_spectra_h2_sto3g.png

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.

H, qubits = qchem.molecular_hamiltonian(name, geometry, mapping="jordan_wigner")

print("Number of qubits = ", qubits)
print("Hamiltonian is ", H)

Out:

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

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\),

\[\hat{S}^2 = \frac{3}{4} N_e + \sum_{\alpha, \beta, \gamma, \delta} \langle \alpha, \beta \vert \hat{s}_1 \cdot \hat{s}_2 \vert \gamma, \delta \rangle ~ \hat{c}_\alpha^\dagger \hat{c}_\beta^\dagger \hat{c}_\gamma \hat{c}_\delta.\]

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.

electrons = 2
S2 = qchem.spin2(electrons, qubits, mapping="jordan_wigner")
print(S2)

Out:

(0.75) [I0]
+ (0.375) [Z1]
+ (-0.375) [Z0 Z1]
+ (0.125) [Z0 Z2]
+ (0.375) [Z0]
+ (-0.125) [Z0 Z3]
+ (-0.125) [Z1 Z2]
+ (0.125) [Z1 Z3]
+ (0.375) [Z2]
+ (0.375) [Z3]
+ (-0.375) [Z2 Z3]
+ (0.125) [Y0 X1 Y2 X3]
+ (0.125) [Y0 Y1 X2 X3]
+ (0.125) [Y0 Y1 Y2 Y3]
+ (-0.125) [Y0 X1 X2 Y3]
+ (-0.125) [X0 Y1 Y2 X3]
+ (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],

\[\hat{U}(\vec{\theta}) = \prod_{p > r} \mathrm{exp} \Big\{\theta_{pr}(\hat{c}_p^\dagger \hat{c}_r-\mathrm{H.c.}) \Big\} \prod_{p > q > r > s} \mathrm{exp} \Big\{\theta_{pqrs} (\hat{c}_p^\dagger \hat{c}_q^\dagger \hat{c}_r \hat{c}_s-\mathrm{H.c.}) \Big\}.\]

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,

\[E(\vec{\theta}) = \langle \mathrm{HF} \vert \hat{U}^\dagger(\vec{\theta}) \hat{H} \hat{U}(\vec{\theta}) \vert \mathrm{HF} \rangle,\]

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,

\[S = -\frac{1}{2} + \sqrt{\frac{1}{4} + \langle \hat{S}^2 \rangle}.\]

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
prev_energy = cost_fn(params)
for n in range(max_iterations):
    params = opt.step(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

    prev_energy = energy

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
prev_energy = cost_fn(params)
for n in range(max_iterations):
    params = opt.step(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

    prev_energy = energy

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: ( 4 minutes 51.140 seconds)

Gallery generated by Sphinx-Gallery