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

## What You Need

• A Web browser and a Google Account

## 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.

In a browser, go to

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

## Importing Libraries

Execute these commands to install:
• pennylane for quantum computing
• openfermionpyscf to make quantum algorithms to model electron states
• qchem for quantum chemistry
• numpy to perform math on arrays
• pyplot to draw graphs
 ```!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 + HH + H2

## 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

Posted 8-2-24