ML 170: Modeling Chemical Reactions with ML and Quantum Computing (10 pts)

What You Need

Purpose

To combine a simulated quantum computer and Machine Learning to solve some quantum chemistry problems.

Background

Consider a Hydrogen molecule, H2.

The diagram below shows the energy of the molecule as the distance between the two H atoms changes. I got that diagram from this page.

Normally the molecule sits at its equilibrium distance, where the energy is at its minimum.

Squeezing the atoms together causes them to repel one another, like compressing a spring.

Pulling the atoms apart requires inputting energy, like sending a rocket out of the Earth's gravity. When the distance becomes large, the atoms break apart, becoming separate H atoms.

The goal of this project is to model that energy curve, which determines how much energy is consumed or released in various chemical reactions.

Using Google Colab

In a browser, go to
https://colab.research.google.com/
If you see a blue "Sign In" button at the top right, click it and log into a Google account.

From the menu, click File, "New notebook".

Importing Libraries

Execute these commands to install: 
!pip install pennylane
!pip install openfermionpyscf

import pennylane as qml
from pennylane import qchem
from pennylane import numpy as np
import matplotlib.pyplot as plt
As shown below, the commands run without errors, so the versions are OK.

Defining a Hartree-Fock State

We'll model this molecule with four qubits, two for each atom. The first qubit represents the probability that the electron is in its lowest-energy state, and the second qubit represents a higher energy state.

We begin with both H atoms completely in their lowest-energy state, which is denoted as:

|1100>
Execute the code below to prepare the qubits we need. 
# Hartree-Fock state
hf = qml.qchem.hf_state(electrons=2, orbitals=4)
print(hf)
As shown below, the qubits are initialized.

Preparing Variables

Execute the commands below to prepare variables we'll use to save the predicted energy curve. 
# atomic symbols defining the molecule
symbols = ['H', 'H']

# list to store energies
energies = []

# set up a loop to change bond length
r_range = np.arange(0.5, 5.0, 0.25)

# keeps track of points in the potential energy surface
pes_point = 0
As shown below, the code runs without errors, and produces no output.

Defining the Quantum Device

Execute these commands to prepare a quantum computing device to simulate the H2 molecule. 
def circuit(parameters):
    # Prepare the HF state: |1100>
    qml.BasisState(hf, wires=range(qubits))
    qml.DoubleExcitation(parameters[0], wires=[0, 1, 2, 3])
    qml.SingleExcitation(parameters[1], wires=[0, 2])
    qml.SingleExcitation(parameters[2], wires=[1, 3])

    return qml.expval(H)  # we are interested in minimizing this expectation value

# Change only the z coordinate of one atom
coordinates = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])

# Construct the Molecule object
molecule = qchem.Molecule(symbols, coordinates)

# Obtain the qubit Hamiltonian
H, qubits = qchem.molecular_hamiltonian(molecule, method='openfermion')

# initialize the gate parameters
params = np.zeros(3, requires_grad=True)

drawer = qml.draw(circuit)
print(drawer(params))
As shown below, a diagram of the quantum circuit appears. The four inputs enter from the left, and pass through the gates, including some entanglements.

Calculating the Minimum Energy

Execute these commands to define a function that will find the energy of the hydrogen molecule with the two atoms a distance r from one another. That calculation involves machine learning, with up to 50 steps of gradient descent, stopping when the energy change becomes very small (less than 1e-6). 
# Set a fixed seed for reproducibility
np.random.seed(42)

def find_energy(r):
    global pes_point, params_old
    # Change only the z coordinate of one atom
    coordinates = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, r]])

    # Construct the Molecule object
    molecule = qchem.Molecule(symbols, coordinates)

    # Obtain the qubit Hamiltonian
    H, qubits = qchem.molecular_hamiltonian(molecule, method='openfermion')

    # define the device, optimizer and circuit
    dev = qml.device("default.qubit", wires=qubits)
    opt = qml.GradientDescentOptimizer(stepsize=0.4)

    @qml.qnode(dev, interface='autograd')
    def circuit(parameters):
        # Prepare the HF state: |1100>
        qml.BasisState(hf, wires=range(qubits))
        qml.DoubleExcitation(parameters[0], wires=[0, 1, 2, 3])
        qml.SingleExcitation(parameters[1], wires=[0, 2])
        qml.SingleExcitation(parameters[2], wires=[1, 3])

        return qml.expval(H)  # we are interested in minimizing this expectation value

    # initialize the gate parameters
    params = np.zeros(3, requires_grad=True)

    # initialize with converged parameters from previous point
    if pes_point > 0:
        params = params_old

    prev_energy = 0.0
    for n in range(50):
        # perform optimization step
        params, energy = opt.step_and_cost(circuit, params)

        if np.abs(energy - prev_energy) < 1e-6:
            break
        prev_energy = energy

    # store the converged parameters
    params_old = params
    pes_point = pes_point + 1

    return energy
As shown below, code executes without errors, and produces no output.

Calculating the Energy Surface

Execute these commands to step through various distances, calculating the energy of the H2 molecule for each distance. 
energies = []

for r in r_range:
    energies.append(find_energy(r))

fig, ax = plt.subplots()
ax.plot(r_range, energies)

ax.set(
    xlabel="Bond length (Bohr)",
    ylabel="Total energy (Hartree)",
    title="Potential energy surface for H$_2$ dissociation",
)
ax.grid()
plt.show()

# equilibrium energy
e_eq = min(energies)
# energy when atoms are far apart
e_dis = energies[-1]

# Bond dissociation energy
bond_energy = e_dis - e_eq

# Equilibrium bond length
idx = energies.index(e_eq)
bond_length = r_range[idx]

print()
print(f"The equilibrium bond length is {bond_length:.1f} Bohrs")
print(f"The bond dissociation energy is {bond_energy:.3f} Hartrees")

ML 170.1: Dissociation Energy (10 pts)

The flag is covered by a green rectangle in the image below.

Hydrogen Exchange Reaction

Consider this reaction, in which a hydrogen atom moves from one partner to another:
H2 + H → H + H2
The animation below shows the process, from this page.

Defining a Hartree-Fock State

Execute these commands to prepare a quantum computing device to simulate the H + H2 system. 
symbols = ["H", "H", "H"]
multiplicity = 2

energies = []
pes_point = 0

# get all the singles and doubles excitations, and Hartree-Fock state
electrons = 3
orbitals = 6
singles, doubles = qchem.excitations(electrons, orbitals)
hf = qml.qchem.hf_state(electrons, orbitals)

print(hf)
As shown below, the three qubits are initialized.

Defining the Quantum Device

Execute these commands to prepare a quantum computing device to simulate the H + H2 system. 
from pennylane.templates import AllSinglesDoubles

def circuit(parameters):
    AllSinglesDoubles(parameters, range(qubits), hf, singles, doubles)
    return qml.expval(H)  # we are interested in minimizing this expectation value

coordinates = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, r], [0.0, 0.0, 4.0]])

# We now specify the multiplicity
molecule = qchem.Molecule(symbols, coordinates, mult=multiplicity)

H, qubits = qchem.molecular_hamiltonian(molecule, method='openfermion')
    
# initialize the gate parameters
params = np.zeros(len(singles) + len(doubles), requires_grad=True)

drawer = qml.draw(circuit)
print(drawer(params))
As shown below, a diagram of the quantum circuit appears. The six inputs enter from the left, and pass through AllSinglesDoubles template, which conveniently performs the necessary calculations for this problem.

Calculating the Minimum Energy

Execute these commands to define a function that will find the energy of the H + H2 system with the middle atom a distance r from the left one. That calculation involves machine learning, with up to 60 steps of gradient descent, stopping when the energy change becomes very small (less than 1e-6). 
# Set a fixed seed for reproducibility
np.random.seed(42)

def find_energy(r):
    global pes_point, params_old
    
    coordinates = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, r], [0.0, 0.0, 4.0]])

    # We now specify the multiplicity
    molecule = qchem.Molecule(symbols, coordinates, mult=multiplicity)

    H, qubits = qchem.molecular_hamiltonian(molecule, method='openfermion')

    dev = qml.device("default.qubit", wires=qubits)
    opt = qml.GradientDescentOptimizer(stepsize=1.5)

    @qml.qnode(dev, interface='autograd')
    def circuit(parameters):
        AllSinglesDoubles(parameters, range(qubits), hf, singles, doubles)
        return qml.expval(H)  # we are interested in minimizing this expectation value

    params = np.zeros(len(singles) + len(doubles), requires_grad=True)

    if pes_point > 0:
        params = params_old

    prev_energy = 0.0

    for n in range(60):
        params, energy = opt.step_and_cost(circuit, params)
        if np.abs(energy - prev_energy) < 1e-6:
            break
        prev_energy = energy

    # store the converged parameters
    params_old = params
    pes_point = pes_point + 1
    
    return energy
As shown below, code executes without errors, and produces no output.

Calculating the Energy Surface

Execute these commands to step through various distances, calculating the energy of the H + H2 system for each distance. 
energies = []
r_range = np.arange(1.0, 3.0, 0.1)

for r in r_range:
    energies.append(find_energy(r))

fig, ax = plt.subplots()
ax.plot(r_range, energies)

ax.set(
    xlabel="Bond length (Bohr)",
    ylabel="Total energy (Hartree)",
    title="Energy surface for Hydrogen exchange ",
)
ax.grid()
plt.show()

# Energy of the reactants and products - two minima on the PES
e_eq1 = min(energies)
e_eq2 = min([x for x in energies if x != e_eq1])

idx1 = energies.index(e_eq1)
idx2 = energies.index(e_eq2)

# Transition state is the local maximum between reactant and products
idx_min = min(idx1, idx2)
idx_max = max(idx1, idx2)

# Transition state energy
energy_ts = max(energies[idx_min:idx_max])

# Activation energy
activation_energy = energy_ts - e_eq1

print()
print(f"The activation energy is {activation_energy:.3f} Hartrees")

ML 170.2: Activation Energy (10 pts)

The flag is covered by a green rectangle in the image below.

References

How to Build a Quantum Artificial Intelligence Model -- With Python Code Examples

Quantum computing simulation of the hydrogen molecular ground-state energies with limited resources

Modelling chemical reactions on a quantum computer

Posted 8-2-24