Hamiltonians for Quantum Chemistry
Let's begin with a brief overview of the role Hamiltonians play in VQE.
The Hamiltonian in VQE Overview
Dr. Victoria Lipinska walks us through Hamiltonians and how to map them for use in quantum computing.
References
The following articles are referenced in the above video.
- Quantum Algorithms for Fermionic Simulations, Ortiz, et al.
- Simulating Chemistry using Quantum Computers, Kassal et al.
- A Comparison of the Bravyi–Kitaev and Jordan–Wigner Transformations for the Quantum Simulation of Quantum Chemistry, Tranter, et al.
- Quantum Chemistry in the Age of Quantum Computing, Cao, et al.
- Quantum computational chemistry, McArdle, et al.
- The Bravyi-Kitaev transformation for quantum computation of electronic structure, Seeley, et al., McArdle, et al.
Preparing Hamiltonians for Quantum Chemistry
A good first step in applying quantum computing to a chemistry problem is defining a Hamiltonian for the system of interest. Here, we will restrict the discussion to quantum chemistry Hamiltonians, as those Hamiltonians require some mapping specific to systems of identical fermions.
As someone working in quantum chemistry, you probably already have your favorite software for modeling molecules, which can generate a Hamiltonian that describes your system of interest. Here, we will use code built solely on PySCF, numpy, and Qiskit. But the process of Hamiltonian preparation transfers to prepackaged solutions as well. The only difference between this approach and other software will be minor syntax differences; some of these are addressed in the "Third-party software" subsection to facilitate integration of existing workflows.
Generating a quantum chemistry Hamiltonian for use on IBM Quantum systems involves the following steps:
- Define your molecule (geometry, spin, active space, etc.)
- Generate the fermionic Hamiltonian (creation and annihilation operators)
- Map from the fermionic Hamiltonian to a bosonic operator (in this context, using Pauli operators)
- If using third-party software: Handle any syntax mismatches between the generating software and Qiskit
The fermionic Hamiltonian is written in terms of fermionic operators, and in particular, takes into account that electrons are indistinguishable fermions. That means they obey completely different statistics from distinguishable, bosonic qubits. Hence the mapping process.
Those of you already familiar with these processes can likely skip this section.
Goal:
The end goal is to obtain a Hamiltonian of the form:
H = [(1, "XX"), (1, "YY"), (1, "ZZ")]
print(H)Output:
[(1, 'XX'), (1, 'YY'), (1, 'ZZ')]
Or
from qiskit.quantum_info import SparsePauliOp
H = SparsePauliOp(["XX", "YY", "ZZ"], coeffs=[1.0 + 0.0j, 1.0 + 0.0j, 1.0 + 0.0j])
print(H)Output:
SparsePauliOp(['XX', 'YY', 'ZZ'],
coeffs=[1.+0.j, 1.+0.j, 1.+0.j])
We'll start by importing some packages:
import numpy as np
from pyscf import ao2mo, gto, mcscf, scf- Define your molecule
Here we will specify attributes of the molecule of interest. In this example, we've chosen diatomic hydrogen (because the resulting Hamiltonians are short enough to display).
The Python-based Simulations of Chemistry Framework (PySCF) has a wide collection of electronic structure modules that can be used to, among other things, generate molecular Hamiltonians suitable for quantum computation. The PySCF Quickstart guide is an excellent resource for a full description of all the variables and functionality. We will give only the most cursory overview, since this will already be familiar to many of you. To understand these better please visit PySCF. Briefly:
distance can be used for diatomic molecules, or simply specify Cartesian coordinates for each atom. Distances are in units of Angstrom.
gto generates gaussian-type orbitals.
basis refers to the functions used to model molecular orbitals. Here 'sto-6g' is a common minimal basis, named for fitting Slater-Type Orbitals using 6 primitive Gaussian orbitals.
spin an integer value indicating the number of unpaired electrons (equal to ). Note that some software uses multiplicity instead ().
charge the charge of the molecule.
symmetry - the point symmetry group of the molecule, either specified with a string or automatically detected by setting "symmetry = True". Here "Dooh" is the appropriate symmetry group for diatomic molecules with two of the same atom species.
distance = 0.735
a = distance / 2
mol = gto.Mole()
mol.build(
verbose=0,
atom=[
["H", (0, 0, -a)],
["H", (0, 0, a)],
],
basis="sto-6g",
spin=0,
charge=0,
symmetry="Dooh",
)Output:
<pyscf.gto.mole.Mole at 0x7f45c63899d0>
Keep in mind that one can describe total energy (which includes nuclear repulsion energy as well as electronic), total electronic orbital energy, or the energy of some subset of electronic orbitals (with the complementary subset frozen). In the specific case of , note the different energies below, and note that the total energy less the nuclear repulsion energy does in fact yield the electronic energy:
print(
mol.energy_nuc(),
mol.energy_elec()[0],
mol.energy_tot(),
mol.energy_tot() - mol.energy_nuc(),
)Output:
0.7199689944489797 -1.8455976628764188 -1.125628668427439 -1.8455976628764188
mol.energy_elec()Output:
(-1.8455976628764188, 0.67565609246518)
active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)- Generate fermionic Hamiltonian
scf refers to a wide range of self-consistent field methods.
rhf as in mf = scf.RHF(mol) in mf is a solver that uses the Restricted Hartree Fock calculation. The kernel of this (E, below) is the total energy, including nuclear repulsion and molecular orbitals.
mcscf is a multi-configuration self-consistent fields package.
ao2mo is a transformation from atomic orbitals to molecular orbitals.
We also use the following variables:
ncas: number of orbitals in the complete active space
nelecas: number of electrons in the complete active space
mf = scf.RHF(mol)
E1 = mf.kernel()
mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))
mo = mx.sort_mo(active_space, base=0)
E2 = mx.kernel(mo)[:2]We want a Hamiltonian, and this is often separated into energy of an electronic core (ecore, not involved in minimization), single-electron operators (h1e), and two-electron energies (h2e). These are explicitly extracted below in the last two lines.
h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)These Hamiltonians are currently fermionic (creation and annihilation) operators, applicable to systems of (indistinguishable) fermions, and correspondingly subject to antisymmetry under exchange. This results in different statics than would apply to a distinguishable or bosonic system. To run calculations on IBM Quantum Systems, we require a bosonic operator describing the energy. The result of such a mapping is conventionally written in terms of Pauli operators, since they are both Hermitian and unitary. There are several mappings one can use. One of the simplest is the Jordan Wigner transformation.
- Mapping the Hamiltonian
It should be noted that there are many tools available for mapping a chemical Hamiltonian to one suitable for running on a quantum computer. Here, we implement the Jordan Wigner mapping directly using only PySCF, numpy, and Qiskit. We comment below on syntax considerations for other solutions.
The Cholesky function helps us obtain a low-rank decomposition of the two-electron terms in the Hamiltonian.
def cholesky(V, eps):
# see https://arxiv.org/pdf/1711.02242.pdf section B2
# see https://arxiv.org/abs/1808.02625
# see https://arxiv.org/abs/2104.08957
no = V.shape[0]
chmax, ng = 20 * no, 0
W = V.reshape(no**2, no**2)
L = np.zeros((no**2, chmax))
Dmax = np.diagonal(W).copy()
nu_max = np.argmax(Dmax)
vmax = Dmax[nu_max]
while vmax > eps:
L[:, ng] = W[:, nu_max]
if ng > 0:
L[:, ng] -= np.dot(L[:, 0:ng], (L.T)[0:ng, nu_max])
L[:, ng] /= np.sqrt(vmax)
Dmax[: no**2] -= L[: no**2, ng] ** 2
ng += 1
nu_max = np.argmax(Dmax)
vmax = Dmax[nu_max]
L = L[:, :ng].reshape((no, no, ng))
print(
"accuracy of Cholesky decomposition ",
np.abs(np.einsum("prg,qsg->prqs", L, L) - V).max(),
)
return L, ngThe functions identity and creators_destructors replace creation and annihilation operators in the fermionic Hamiltonian with Pauli operators; creators_destructors uses the Jordan-Wigner mapping.
def identity(n):
return SparsePauliOp.from_list([("I" * n, 1)])
def creators_destructors(n, mapping="jordan_wigner"):
c_list = []
if mapping == "jordan_wigner":
for p in range(n):
if p == 0:
ell, r = "I" * (n - 1), ""
elif p == n - 1:
ell, r = "", "Z" * (n - 1)
else:
ell, r = "I" * (n - p - 1), "Z" * p
cp = SparsePauliOp.from_list([(ell + "X" + r, 0.5), (ell + "Y" + r, 0.5j)])
c_list.append(cp)
else:
raise ValueError("Unsupported mapping.")
d_list = [cp.adjoint() for cp in c_list]
return c_list, d_listFinally, build_hamiltonian uses the cholesky, identity, and creators_destructors functions to create the final Hamiltonian suitable for running on a quantum computer.
def build_hamiltonian(ecore: float, h1e: np.ndarray, h2e: np.ndarray) -> SparsePauliOp:
ncas, _ = h1e.shape
C, D = creators_destructors(2 * ncas, mapping="jordan_wigner")
Exc = []
for p in range(ncas):
Excp = [C[p] @ D[p] + C[ncas + p] @ D[ncas + p]]
for r in range(p + 1, ncas):
Excp.append(
C[p] @ D[r]
+ C[ncas + p] @ D[ncas + r]
+ C[r] @ D[p]
+ C[ncas + r] @ D[ncas + p]
)
Exc.append(Excp)
# low-rank decomposition of the Hamiltonian
Lop, ng = cholesky(h2e, 1e-6)
t1e = h1e - 0.5 * np.einsum("pxxr->pr", h2e)
H = ecore * identity(2 * ncas)
# one-body term
for p in range(ncas):
for r in range(p, ncas):
H += t1e[p, r] * Exc[p][r - p]
# two-body term
for g in range(ng):
Lg = 0 * identity(2 * ncas)
for p in range(ncas):
for r in range(p, ncas):
Lg += Lop[p, r, g] * Exc[p][r - p]
H += 0.5 * Lg @ Lg
return H.chop().simplify()Finally, we use build_hamiltonian to construct our qubit Hamiltonian from Pauli operators by using the Jordan-Wigner transformation. This also gives us the accuracy of the Cholesky decomposition we used.
H = build_hamiltonian(ecore, h1e, h2e)
print(H)Output:
accuracy of Cholesky decomposition 2.220446049250313e-16
SparsePauliOp(['IIII', 'IIIZ', 'IZII', 'IIZI', 'ZIII', 'IZIZ', 'IIZZ', 'ZIIZ', 'IZZI', 'ZZII', 'ZIZI', 'YYYY', 'XXYY', 'YYXX', 'XXXX'],
coeffs=[-0.09820182+0.j, -0.1740751 +0.j, -0.1740751 +0.j, 0.2242933 +0.j,
0.2242933 +0.j, 0.16891402+0.j, 0.1210099 +0.j, 0.16631441+0.j,
0.16631441+0.j, 0.1210099 +0.j, 0.17504456+0.j, 0.04530451+0.j,
0.04530451+0.j, 0.04530451+0.j, 0.04530451+0.j])
Let's just briefly note two important points to consider when constructing the fermionic operators for a molecule. As the molecule type changes, the symmetry will change. Similarly, the number of orbitals with various symmetries, like the cylindrically symmetric "A1", will change. These changes are evident even with the simple extension to LiH, as seen here:
distance = 1.56
mol = gto.Mole()
mol.build(
verbose=0,
atom=[["Li", (0, 0, 0)], ["H", (0, 0, distance)]],
basis="sto-6g",
spin=0,
charge=0,
symmetry="Coov",
)
mf = scf.RHF(mol)
E1 = mf.kernel()
# %% ----------------------------------------------------------------------------------------------
mx = mcscf.CASCI(mf, ncas=5, nelecas=(1, 1))
cas_space_symmetry = {"A1": 3, "E1x": 1, "E1y": 1}
mo = mcscf.sort_mo_by_irrep(mx, mf.mo_coeff, cas_space_symmetry)
E2 = mx.kernel(mo)[:2]
h1e, ecore = mx.get_h1eff()
h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)It is also worth noting that one can quickly lose intuition for the final resulting Hamiltonian. The Hamiltonian for LiH (using the Jordan-Wigner mapper) already consists of 276 terms.
len(build_hamiltonian(ecore, h1e, h2e))Output:
accuracy of Cholesky decomposition 1.1102230246251565e-16
276
When in doubt regarding symmetries, one can also generate some symmetry information for the molecule by setting symmetry = True and verbose = 4:
distance = 1.56
mol = gto.Mole()
mol.build(
verbose=4,
atom=[["Li", (0, 0, 0)], ["H", (0, 0, distance)]],
basis="sto-6g",
spin=0,
charge=0,
symmetry=True,
)Output:
System: uname_result(system='Linux', node='IBM-R912JTRT', release='5.10.102.1-microsoft-standard-WSL2', version='#1 SMP Wed Mar 2 00:30:59 UTC 2022', machine='x86_64', processor='x86_64') Threads 16
Python 3.7.6 (default, Jan 8 2020, 19:59:22)
[GCC 7.3.0]
numpy 1.21.6 scipy 1.7.3
Date: Sun Oct 15 04:24:52 2023
PySCF version 2.1.1
PySCF path /home/porter284/anaconda3/lib/python3.7/site-packages/pyscf
[CONFIG] conf_file None
[INPUT] verbose = 4
[INPUT] num. atoms = 2
[INPUT] num. electrons = 4
[INPUT] charge = 0
[INPUT] spin (= nelec alpha-beta = 2S) = 0
[INPUT] symmetry True subgroup None
[INPUT] Mole.unit = angstrom
[INPUT] Symbol X Y Z unit X Y Z unit Magmom
[INPUT] 1 Li 0.000000000000 0.000000000000 0.000000000000 AA 0.000000000000 0.000000000000 0.000000000000 Bohr 0.0
[INPUT] 2 H 0.000000000000 0.000000000000 1.560000000000 AA 0.000000000000 0.000000000000 2.947972754321 Bohr 0.0
nuclear repulsion = 1.01764848253846
point group symmetry = Coov
symmetry origin: [0. 0. 0.73699319]
symmetry axis x: [1. 0. 0.]
symmetry axis y: [0. 1. 0.]
symmetry axis z: [0. 0. 1.]
num. orbitals of irrep A1 = 4
num. orbitals of irrep E1x = 1
num. orbitals of irrep E1y = 1
number of shells = 4
number of NR pGTOs = 36
number of NR cGTOs = 6
basis = sto-6g
ecp = {}
CPU time: 69.98
<pyscf.gto.mole.Mole at 0x7f5e9c7252d0>
Among other useful information, this returns both point group symmetry = Coov and also the number of orbitals in each irreducible representation.
point group symmetry = Coov
num. orbitals of irrep A1 = 4
num. orbitals of irrep E1x = 1
num. orbitals of irrep E1y = 1
number of shells = 4
This does not necessarily tell you how many orbitals you want included in your active space, but it helps you see what orbitals are present and their symmetries.
Specifying symmetry and orbitals is often helpful, but you can also specify the number of orbitals you want to include. Consider the case of ethene, below. Using verbose = 4, we can print the symmetries of the various orbitals:
# Replace these variables with correct distances:
a = 1
b = 1
c = 1
# Build
mol = gto.Mole()
mol.build(
verbose=4,
atom=[
["C", (0, 0, a)],
["C", (0, 0, -a)],
["H", (0, c, b)],
["H", (0, -c, b)],
["H", (0, c, -b)],
["H", (0, -c, -b)],
],
basis="sto-6g",
spin=0,
charge=0,
symmetry=True,
)We obtain:
num. orbitals of irrep Ag = 4
num. orbitals of irrep B2g = 2
num. orbitals of irrep B3g = 1
num. orbitals of irrep B1u = 4
num. orbitals of irrep B2u = 1
num. orbitals of irrep B3u = 2
But rather than specifying all the orbitals by symmetry, we can simply write:
active_space = range(mol.nelectron // 2 - 2, mol.nelectron // 2 + 2)In this approach, we take several orbitals near the filling level (valence and unoccupied). Here, 5 orbitals have been selected for inclusion in the active space (the 6th though 10th).
print(
mol.nelectron // 2 - 2,
mol.nelectron // 2 + 2,
)Output:
6 10
- Third-Party Software
There are several software packages developed for quantum chemistry, some offering multiple mappers and tools for restricting active spaces. The steps described above are general and apply to 3rd party software as well. But this other software may return Hamiltonians in a format that is not accepted by Qiskit. This is the case, for example, for OpenFermion and for Tangelo. Both of those return Hamiltonians of the form:
H = -0.042 [] + -0.045 [X0 X1 Y2 Y3] + ... + 0.178 [Z0] + ... + 0.176 [Z2 Z3] + -0.243 [Z3]
Note in particular that the gates are numbered, and identity operators are not shown. This is in contrast to Hamiltonians used in Qiskit, which would write the term [Z2 Z3] as "IIZZ" (qubits 0 and 1 being acting on by the identity operator, qubits 2 and 3 being acted on by the Z operator).
To accommodate any existing workstreams you have, the code block below converts from one syntax to the other. The function convert_openfermion_to_qiskit takes as its arguments a Hamiltonian generated in OpenFermion or Tangelo (and already mapped onto Pauli operators using any available mapper), and the number of qubits needed for the molecule.
from openfermion import QubitOperator
from qiskit.quantum_info import SparsePauliOp
def convert_openfermion_to_qiskit(
openfermion_operator: QubitOperator, num_qubits: int
) -> SparsePauliOp:
terms = openfermion_operator.terms
labels = []
coefficients = []
for term, constant in terms.items():
# Default set to identity
operator = list("I" * num_qubits)
# Iterate through PauliSum and replace I with Pauli
for index, pauli in term:
operator[index] = pauli
label = "".join(operator)
labels.append(label)
coefficients.append(constant)
return SparsePauliOp(labels, coefficients)You should now have an arsenal of tools for obtaining the Hamiltonian you need to perform quantum chemistry calculations on IBM Quantum Systems.