Fermionic operators¶

A toy model¶

Our first example is a toy Hamiltonian inspired by the Hückel method, which is a method for describing molecules with alternating single and double bonds. Our toy model is a simplified version of the Hückel Hamiltonian and assumes only two orbitals and a single electron.

This Hamiltonian can be constructed with pre-defined values $\alpha = 0.01$ and $\beta = -0.02.$

h1 = 0.01 * (FermiC(0) * FermiA(0) + FermiC(1) * FermiA(1))
h2 = -0.02 * (FermiC(0) * FermiA(1) + FermiC(1) * FermiA(0))
h = h1 + h2
print(h)

0.01 * a⁺(0) a(0)
+ 0.01 * a⁺(1) a(1)
+ -0.02 * a⁺(0) a(1)
+ -0.02 * a⁺(1) a(0)

The fermionic Hamiltonian can be converted to the qubit Hamiltonian with:

h = jordan_wigner(h)

The matrix representation of the qubit Hamiltonian in the computational basis can be diagonalized to get its eigenpairs.

import numpy as np

val, vec = np.linalg.eigh(h.sparse_matrix().toarray())
print(f"eigenvalues:\n{val}")
print()
print(f"eigenvectors:\n{np.real(vec.T)}")

eigenvalues:
[-0.01  0.    0.02  0.03]

eigenvectors:
[[-0.         -0.70710678 -0.70710678 -0.        ]
 [ 1.          0.          0.          0.        ]
 [ 0.          0.          0.          1.        ]
 [ 0.         -0.70710678  0.70710678  0.        ]]

Hydrogen molecule¶

The second quantized molecular electronic Hamiltonian is usually constructed as

where $\alpha$ and $\beta$ denote the electron spin and $p, q, r, s$ are the orbital indices. The coefficients $c$ are integrals over molecular orbitals that are obtained from Hartree-Fock calculations. These integrals can be computed with PennyLane using the electron_integrals() function. We can build the molecular Hamiltonian for the hydrogen molecule as an example. We first define the atom types and the atomic coordinates.

import pennylane as qml
from jax import numpy as jnp

symbols = ["H", "H"]
geometry = jnp.array([[-0.67294, 0.0, 0.0], [0.67294, 0.0, 0.0]])

Then we compute the one- and two-electron integrals, which are the coefficients $c$ in the second quantized molecular Hamiltonian defined above. We also obtain the core constant, which is later used to calculate the contribution of the nuclear energy to the Hamiltonian.

mol = qml.qchem.Molecule(symbols, geometry)
core, one, two = qml.qchem.electron_integrals(mol)()

These integrals are computed over molecular orbitals. Each molecular orbital contains a pair of electrons with different spins. We have assumed that the spatial distribution of these electron pairs is the same to simplify the calculation of the integrals. However, to properly account for all electrons, we need to duplicate the integrals for electrons with the same spin. For example, the $pq$ integral, which is the integral over the orbital $p$ and the orbital $q,$ can be used for both spin-up and spin-down electrons. Then, if we have a $2 \times 2$ matrix of such integrals, it will become a $4 \times 4$ matrix. The code block below simply extends the integrals by duplicating terms to account for both spin-up and spin-down electrons.

for i in range(4):
    if i < 2:
        one = one.repeat(2, axis=i)
    two = two.repeat(2, axis=i)

We can now construct the fermionic Hamiltonian for the hydrogen molecule. The one-body terms, which are the first part in the Hamiltonian above, can be added first. We will use itertools to efficiently create all the combinations we need. Some of these combinations are not allowed because of spin restrictions and we need to exclude them. You can find more details about constructing a molecular Hamiltonian in reference 1.

import itertools

n = one.shape[0]

h = 0.0

for p, q in itertools.product(range(n), repeat=2):
    if p % 2 == q % 2:  # to account for spin-forbidden terms
        h += one[p, q] * FermiC(p) * FermiA(q)

The two-body terms can be added with:

for p, q, r, s in itertools.product(range(n), repeat=4):
    if p % 2 == s % 2 and q % 2 == r % 2:  # to account for spin-forbidden terms
        h += two[p, q, r, s] / 2 * FermiC(p) * FermiC(q) * FermiA(r) * FermiA(s)

We then simplify the Hamiltonian to remove terms with negligible coefficients and then map it to the qubit basis.

h.simplify()
h = jordan_wigner(h)

We also need to include the contribution of the nuclear energy.

h += np.sum(core * qml.Identity(0))

This gives us the qubit Hamiltonian which can be used as an input for quantum algorithms. We can also compute the ground-state energy by diagonalizing the matrix representation of the Hamiltonian in the computational basis.

np.linalg.eigh(h.sparse_matrix().toarray())[0].min()

-1.1368472302251442