Advanced usage

We will describe here the advanced usage of the qat.fermion module.

Variational Quantum Eigensolver (VQE)

The VQE is one of the flagship algorithms using near-term quantum computers. We present in this section how to do a simple VQE using qat.fermion.

We will define manually a Hubbard Hamiltonian, using one- and two-body integrals, and use them to initialize an ElectronicStructureHamiltonian. We will then convert this Hamiltonian to its spin representation, and find its ground state energy via VQE. To do so, we will manually create a parametric circuit (ansatz), whose parameters we will optimize.


In the following section, we define the Hubbard model manually, but we could also use the Hamiltonian constructor for the Hubbard model.

Let us define the Hubbard model:

import numpy as np
from qat.fermion import ElectronicStructureHamiltonian
from qat.lang.AQASM import H, RX, RY, CNOT, Program
from qat.qpus import get_default_qpu
from qat.plugins import ScipyMinimizePlugin

nqbits = 2
u = 2.0
hpq = np.zeros((nqbits,) * 2)
hpq[0,0] = hpq[1,1] = - u/2
hpqrs = np.zeros((nqbits, ) * 4)
hpqrs[0, 1, 0, 1] = hpqrs[1, 0, 1, 0] = - u

hamiltonian = ElectronicStructureHamiltonian(hpq=hpq, hpqrs=hpqrs)

We then convert the fermionic Hamiltonian to a spin Hamiltonian.

hamiltonian_sp = hamiltonian.to_spin()

We compute the eigenenergies of the Hamiltonian matrix to find the exact energy.

exact_energy = min(np.linalg.eigh(hamiltonian_sp.get_matrix()))
>>> print(f"Exact_energy = {exact_energy}")
Exact_energy = -1.0

Let us write the ansatz we will optimize:

# Initialize a Program and the number of qubits
prog = Program()
qbits = prog.qalloc(nqbits)

# Initialize the parameters
theta = prog.new_var(float, "\\theta")
phi = prog.new_var(float, "\\phi")

# Write the ansatz

Finally, we optimize the parameters such that they minimize the expectation value of the Hamiltonian we defined previously.

# Choose the QPU and define the VQE stack
qpu = get_default_qpu()
stack = ScipyMinimizePlugin(x0=[0.5, 1.23], method="COBYLA") | qpu

# Submit the job to the stack
res = stack.submit(prog.to_circ().to_job(observable=hamiltonian_sp))

We find \(\theta\) and \(\phi\):

>>> print(f"=== VQE COBYLA RESULTS ===\nEnergy = {res.value}\ntheta, phi = {res.meta_data['parameters']}")
Energy = -0.9999999974320526
theta, phi = [3.1416940034276815, 1.2147911128415545]

We show in more details a VQE for fermions in our Jupyter notebooks


Many more ansatz are available in qat.fermion. See our list of available fermionic ansatze.

Quantum Chemistry tools for VQE

Various methods are available to quantum chemists for the study of atomic and molecular systems:

  • Unitary coupled-cluster (UCC) ansatz construction,

  • Cluster operator generation,

  • Initial parameters guess via Møller-Plesset perturbation theory,

  • Trotterization,

  • Active space selection.

More information is available in the following notebooks:

We have introduced in qat.fermion the classes MolecularHamiltonian and MoleculeInfo. These classes are helper classes, meaning they are an interface to lower-level functions in qat.fermion, meant to simplify the way you interact with the code. We will see here how to use these classes to study a \(LiH\) molecule. A more thorough study is available in qat.fermion Jupyter notebooks.

We will assume you already have enough information about the molecule itself. You are free to use whatever quantum chemistry library you prefer. For simplicity’s sake, we provide a basic function based on the package PySCF perform_pyscf_computation() which collects the information we need.

We start by inputting the geometry, the basis, the spin and the charge of the molecule.

from qat.fermion.chemistry.pyscf_tools import perform_pyscf_computation

geometry = [("Li", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 1.75))]
basis = "6-31g"
spin = 0
charge = 0

) = perform_pyscf_computation(geometry=geometry, basis=basis, spin=spin, charge=charge)


From there, we can use the MolecularHamiltonian class.

from qat.fermion.chemistry import MolecularHamiltonian

# Define the molecular hamiltonian
mol_h = MolecularHamiltonian(one_body_integrals, two_body_integrals, nuclear_repulsion)
>>> print(mol_h)
 - constant_coeff : 0.9071609330057144
 - integrals shape
    * one_body_integrals : (11, 11)
    * two_body_integrals : (11, 11, 11, 11)
>>> print(mol_h.nqbits)

This class is different from the spin, fermion and electronic-structure Hamiltonian classes we have seen so far. It is defined using interaction integrals. You can convert it to an ElectronicStructureHamiltonian using the method get_electronic_hamiltonian().

MolecularHamiltonian allows:

To illustrate this, we will compute the natural orbitals occupation numbers (NOONs) as well as their basis transformation matrix. We will then use the transform_basis() method to change the basis of our MolecularHamiltonian.

import numpy as np

# Compute NOONs and the basis
noons, basis_change = np.linalg.eigh(rdm1)

# The noons should be in decreasing order. This means we should flip the basis as well.
noons = list(reversed(noons))
basis_change = np.flip(basis_change, axis=1)

# Change the hamiltonian basis
mol_h_new_basis = mol_h.transform_basis(basis_change)

We can then proceed to the active space selection.

# Active space selection
mol_h_active, active_indices, occupied_indices = mol_h_new_basis.select_active_space(
noons=noons, n_electrons=n_electrons, threshold_1=0.02, threshold_2=0.002
>>> print(mol_h_active)
 - constant_coeff : 0.9071609330057144
 - integrals shape
    * one_body_integrals : (2, 2)
    * two_body_integrals : (2, 2, 2, 2)
>>> print(mol_h_active.nqbits)

We reduced the number of qubits from 22 to 4 qubits !



While the MolecularHamiltonian brings simplicity, we can add another layer of simplicity by using the MoleculeInfo class. When restricting the active space, MoleculeInfo takes care of updating its attributes such as the active and occupied orbitals indices, as well as the direct update of the NOONs and orbital energies.

More information is available in the source code documentation. The MoleculeInfo class is used in our VQE UCC ansatz resolution of the \(LiH\) molecule.

Fermionic ansatz circuits

We also provide some circuits useful in a context of VQE on fermionic systems:

Quantum phase estimation

The quantum phase space estimation algorithm allows for the estimation of the eigenvalue associated with a given eigenvector. Let us see how to use this algorithm with a simple example.

We will first define a Hubbard model, and find its exact eigenenergies via direct diagonalization.

import numpy as np
from qat.fermion.hamiltonians import make_hubbard_model

# Define the Hubbard model
U = 1.0
t = 0.2
t_mat = -t * np.array([[0.0, 1.0], [1.0, 0.0]])

# We get an ElectronicStructureHamiltonian
hamiltonian = make_hubbard_model(t_mat, U, mu=U / 2)

# Find its eigenenergies
eigvals = np.linalg.eigvalsh(hamiltonian.get_matrix())
>>> print(np.round(eigvals, decimals=4))
[-1.1403 -1.     -1.     -1.     -0.7    -0.7    -0.7    -0.7    -0.3
 -0.3    -0.3    -0.3     0.      0.      0.      0.1403]

Now we know the exact eigenenergies, let us find them once again using the quantum phase estimation algorithm.

from qat.qpus import get_default_qpu
from qat.fermion.phase_estimation import perform_phase_estimation

qpu = get_default_qpu()

nqbits_phase = 8
n_trotter_steps = 6
guess_energy = -0.1  # try an energy which is off from 0
size_interval = 5

energy, prob = perform_phase_estimation(
>>> print(f"E = {energy}")
E = -0.002343750000000089

We find an energy very close to some of the exact eigenenergies we computed earlier !

More information is available in the source code documentation and in the notebook on quantum phase estimation on the Hubbard molecule.

Trotterization tools

You can trotterize any Hamiltonian using our trotterization tools.

Here is an example.

import numpy as np

from qat.lang.AQASM import Program
from qat.fermion.trotterisation import make_trotterisation_routine
from qat.fermion.hamiltonians import make_hubbard_model

# Define the Hubbard model
U = 1.0
t = 0.2
t_mat = -t * np.array([[0.0, 1.0], [1.0, 0.0]])

# We get an ElectronicStructureHamiltonian
hamiltonian = make_hubbard_model(t_mat, U, mu=U / 2)

# Trotterize the Hamiltonian (with 1 trotter step)
qrout = make_trotterisation_routine(hamiltonian, n_trotter_steps=1, final_time=1)

# Define an empty Program and apply the QRoutine on it
prog = Program()
reg = prog.qalloc(hamiltonian.nbqbits)
prog.apply(qrout, reg)
circ = prog.to_circ()
>>> circ.display()

For more information, see the make_trotterisation_routine() function documentation. The trotterization is used for the UCC ansatz construction in our UCC VQE notebooks, see here for the \(H_2\) molecule, and here for the \(LiH\) molecule.

Quantum subspace expansion

The quantum subspace expansion is a method that can allow to reach a better precision at the cost of doing additional measurements. You will find more details directly in the apply_quantum_subspace_expansion() documentation. See how to use it in the quantum subspace expansion notebook.