How to construct and load Hamiltonians in PennyLane

The Hamiltonian H is a core part of many NISQ algorithms. Depending on the problem that we are interested in - quantum simulation, calculating orbital energies or solving a combinatorial optimization problem, we need to construct H in the appropriate way. In this how-to, we will highlight how PennyLane can be used to build to the Hamiltonians for combinatorial optimization problems, quantum many-body physics and quantum chemistry.

In general, we can write the Hamiltonian as a linear combination of M local operators

with c_m\in \mathbb{R} . To work with Hamiltonians in PennyLane, we make use of the qml.Hamiltonian object (check the docs here), an abstract class that can be used directly with other PennyLane functions to work with observables. For example, consider the 4-qubit Hamiltonian

We can directly define this Hamiltonian in PennyLane by summing qml.Observable objects and multiplying them by scalars:

import pennylane as qml
import pennylane.numpy as np
hamiltonian = 0.5 * qml.PauliX(0) @ qml.PauliX(1) - 0.25 * qml.PauliZ(0) @ qml.PauliY(2) \
+ qml.PauliZ(3)

Note that the matrix multiplication operator @ works as the tensor product \otimes between the Pauli operators. A more general approach for constructing a Hamiltonian in PennyLane is from a list of qml.Observable objects and a corresponding list of coefficients:

coeffs = [0.5, -0.25, 1.]
obs = [qml.PauliX(0) @ qml.PauliX(1), qml.PauliZ(0) @ qml.PauliY(2), qml.PauliZ(3)]
hamiltonian = qml.Hamiltonian(coeffs, obs)

The returned qml.Hamiltonian object also has user-friendly string representation:

>>> print(hamiltonian)
  (1.0) [Z3]
+ (-0.25) [Z0 Y2]
+ (0.5) [X0 X1]

Combinatorial optimization

Classical combinatorial optimization problems are defined over a set of N bits b \in \{0,1\}^N and require minimizing the objective function

Mathematically, there exists a mapping from the binary cost function C(b) to the Ising Hamiltonian on a graph with a polynomial number of steps [Lucas, 2014]. PennyLane offers several functions that perform this mapping to the Ising Hamiltonian for its users.

For example, let us consider the following graph:

and construct the appropriate Hamiltonian for the MaxCut problem on this graph. With NetworkX, we can create an abstract graph class by calling nx.Graph on a tuple of edges:

import networkx as nx
graph = nx.Graph(((0, 1), (2, 0), (3, 2), (2, 1)))

Creating a VQE-ready Hamiltonian is now a piece of cake; we simply call maxcut from the qml.qaoa module on the graph we just created:

from pennylane.qaoa import maxcut
hamiltonian_maxcut, hamiltonian_mixer = maxcut(graph)

This returns a qml.Hamiltonian object corresponding to the cost function.:

>>> print(type(hamiltonian_maxcut))
<class 'pennylane.vqe.vqe.Hamiltonian'>
>>> print(hamiltonian_maxcut)
  (-2.0) [I0]
+ (0.5) [Z0 Z1]
+ (0.5) [Z0 Z2]
+ (0.5) [Z1 Z2]
+ (0.5) [Z2 Z3]

The second output that is returned, hamiltonian_mixer is the commonly used mixer Hamiltonian for a QAOA optimization., Analogously, we can create the Hamiltonians for a number of other combinatorial problems, detailed here.

Quantum many-body Hamiltonians

In quantum many-body physics we are often interested in studying the ground state properties of a quantum spin system. Depending on the system, the spins are located at specific points on a 1D, 2D or 3D lattice. This spatial layout, combined with the different interaction between the spins governs the microscopic and macroscopic physics of the system. By combining NetworkX and PennyLane, we can quickly construct Hamiltonians for quantum many-body systems with different lattice structures.

For instance, we can look at the Heisenberg spin model defined on a rectangular lattice. The Hamiltonian of is this system is given by

with the sum over the edges of the following graph:

After we relabel the nodes to integers, we only have to loop over the edges and add the interactions of the model with the correct coefficients:

graph = nx.generators.lattice.grid_2d_graph(3,3)
graph = nx.relabel.convert_node_labels_to_integers(graph)
obs = []
coeffs = []
for edge in graph.edges():
    coeffs.extend([1.0, 1.0, 1.0])
    obs.extend([qml.PauliX(edge[0]) @ qml.PauliX(edge[1]),
                        qml.PauliY(edge[0]) @ qml.PauliY(edge[1]),
                        qml.PauliZ(edge[0]) @ qml.PauliZ(edge[1])])
hamiltonian_heisenberg = qml.Hamiltonian(coeffs, obs)

Quantum Chemistry

Note: to run the code below, make sure that the pennylane-qchem package is installed with pip!

In quantum chemistry we want to calculate properties of different kinds of molecules. One such property is called the electronic structure of the molecule. Characterizing the electronic structure of a molecule is a notoriously difficult computational problem where quantum computers can potentially outshine classical methods. It should come as no surprise that the physics of this problem can be described with a Hamiltonian.

To create a molecular Hamiltonian in PennyLane, we can either directly define the types and positions types of the atoms in the molecule

from pennylane import qchem
symbols, coordinates = (['H', 'H'], np.array([0., 0., -0.66140414, 0., 0., 0.66140414]))

or load .xyz files from disk (more formats are supported if OpenBabel is installed)

symbols, coordinates = qchem.read_structure('h2.xyz')

Then, we can use the symbols and coordinates to call the function molecular_hamiltonian() to create the Hamiltonian and the qubits required to perform

h, qubits = qchem.molecular_hamiltonian(
        symbols,
        coordinates,
        charge=0,
        mult=1,
        basis='sto-3g',
        active_electrons=2,
        active_orbitals=2
)

which gives the following qml.Hamiltonian object:

>>> print(h)
  (-0.24274280706347268) [Z2]
+ (-0.24274280706347268) [Z3]
+ (-0.04207897085891747) [I0]
+ (0.1777128752691193) [Z0]
+ (0.17771287526911933) [Z1]
+ (0.12293305078371519) [Z0 Z2]
+ (0.12293305078371519) [Z1 Z3]
+ (0.16768319474683535) [Z0 Z3]
+ (0.16768319474683535) [Z1 Z2]
+ (0.17059738347065292) [Z0 Z1]
+ (0.17627640822374238) [Z2 Z3]
+ (-0.04475014396312016) [Y0 Y1 X2 X3]
+ (-0.04475014396312016) [X0 X1 Y2 Y3]
+ (0.04475014396312016) [Y0 X1 X2 Y3]
+ (0.04475014396312016) [X0 Y1 Y2 X3]

To understand all the arguments of this function, refer to the documentation here, or check out this quantum chemistry tutorial.