{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# This cell is added by sphinx-gallery\n# It can be customized to whatever you like\n%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Modelling chemical reactions on a quantum computer\n==================================================\n\n*Author: PennyLane dev team. Posted: 23 July 2021. Last updated: 23 July\n2021*\n\nThe term \"chemical reaction\" is another name for the transformation of\nmolecules -- the breaking and forming of bonds. They are characterized\nby an energy barrier that determines the likelihood that a reaction\ntakes place. The energy landscapes formed by these barriers are the key\nto understanding how chemical reactions occur, at the deepest possible\nlevel.\n\n![An example chemical\nreaction.](/demonstrations/vqe_bond_dissociation/reaction.png){width=\"50%\"}\n\nIn this tutorial, you will learn how to use PennyLane to simulate\nchemical reactions by constructing potential energy surfaces for\nmolecular transformations. In the process, you will learn how quantum\ncomputers can be used to calculate equilibrium bond lengths, activation\nenergy barriers, and reaction rates. As illustrative examples, we use\ntools implemented in PennyLane to study diatomic bond dissociation and\nreactions involving the exchange of hydrogen atoms.\n\nPotential Energy Surfaces\n-------------------------\n\n[Potential energy surfaces\n(PES)](https://en.wikipedia.org/wiki/Potential_energy_surface) describe\nthe energy of molecules for different positions of its atoms. The\nconcept originates from the fact that the electrons are much lighter\nthan protons and neutrons, so they will adjust instantaneously to the\nnew positions of the nuclei. This leads to a separation of the nuclear\nand electronic parts of the Schr\u00f6dinger equation, meaning we only need\nto solve the electronic equation:\n\n$$H(R)|\\Psi \\rangle = E|\\Psi\\rangle.$$\n\nFrom this perspective arises the concept of the electronic energy of a\nmolecule, $E(R)$, as a function of nuclear coordinates $R$. The energy\n$E(R)$ is the expectation value of the molecular Hamiltonian,\n$E(R)=\\langle \\Psi_0|H(R)|\\Psi_0\\rangle$, taken with respect to the\nground state $|\\Psi_0(R)\\rangle$. The potential energy surface is\nprecisely this function $E(R)$, which connects energies to different\ngeometries of the molecule. It gives us a visual tool to understand\nchemical reactions by associating stable molecules (reactants and\nproducts) with local minima, transition states with peaks, and by\nidentifying the possible routes for a chemical reaction to occur.\n\nTo build the potential energy surface, we compute the energy for fixed\npositions of the nuclei, and subsequently adjust the positions of the\nnuclei in incremental steps, computing the energies at each new\nconfiguration. The obtained set of energies corresponds to a grid of\nnuclear positions and the plot of $E(R)$ gives rise to the potential\nenergy surface.\n\n![Illustration of a potential energy surface for a diatomic\nmolecule.](/demonstrations/vqe_bond_dissociation/pes.png){width=\"75%\"}\n\nBond dissociation in a Hydrogen molecule\n----------------------------------------\n\nWe now construct a potential energy surface and use it to compute\nequilibrium bond lengths and the bond dissociation energy. We begin with\nthe simplest of molecules: $H_2$. The formation or breaking of the $H-H$\nbond is also the most elementary of all reactions:\n\n$$H_2 \\rightarrow H + H.$$\n\nUsing a minimal [basis\nset](https://en.wikipedia.org/wiki/STO-nG_basis_sets), this molecular\nsystem can be described by two electrons in four spin-orbitals. When\nmapped to a qubit representation, we need a total of four qubits. The\n*Hartree-Fock (HF) state* is represented as $|1100\\rangle$, where the\ntwo lowest-energy orbitals are occupied, and the remaining two are\nunoccupied.\n\nWe design a quantum circuit consisting of \\~.pennylane.SingleExcitation\nand \\~.pennylane.DoubleExcitation gates applied to the Hartree-Fock\nstate. This circuit will be optimized to prepare ground states for\ndifferent configurations of the molecule.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import pennylane as qml\nfrom pennylane import qchem\n\n# Hartree-Fock state\nhf = qml.qchem.hf_state(electrons=2, orbitals=4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To construct the potential energy surface, we vary the location of the\nnuclei and calculate the energy for each resulting geometry of the\nmolecule. We keep an $H$ atom fixed at the origin and change only the\ncoordinate of the other atom in a single direction. The potential energy\nsurface is then a one-dimensional function depending only on the bond\nlength, i.e., the separation between the atoms. For each value of the\nbond length, we construct the corresponding Hamiltonian, then optimize\nthe circuit using gradient descent to obtain the ground-state energy. We\nvary the bond length in the range $0.5$ to $5.0$\n[Bohrs](https://en.wikipedia.org/wiki/Bohr_radius) in steps of $0.25$\nBohr. This covers the point where the $H-H$ bond is formed, the\nequilibrium bond length, and the point where the bond is broken, which\noccurs when the atoms are far away from each other.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n\n# atomic symbols defining the molecule\nsymbols = ['H', 'H']\n\n# list to store energies\nenergies = []\n\n# set up a loop to change bond length\nr_range = np.arange(0.5, 5.0, 0.25)\n\n# keeps track of points in the potential energy surface\npes_point = 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We build the Hamiltonian using the\n\\~.pennylane\\_qchem.qchem.molecular\\_hamiltonian function, and use\nstandard Pennylane techniques to optimize the circuit.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for r in r_range:\n # Change only the z coordinate of one atom\n coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, r])\n\n # Obtain the qubit Hamiltonian \n H, qubits = qchem.molecular_hamiltonian(symbols, coordinates)\n\n # define the device, optimizer and circuit\n dev = qml.device(\"default.qubit\", wires=qubits)\n opt = qml.GradientDescentOptimizer(stepsize=0.4)\n\n @qml.qnode(dev)\n def circuit(parameters):\n # Prepare the HF state: |1100>\n qml.BasisState(hf, wires=range(qubits))\n qml.DoubleExcitation(parameters[0], wires=[0, 1, 2, 3])\n qml.SingleExcitation(parameters[1], wires=[0, 2])\n qml.SingleExcitation(parameters[2], wires=[1, 3])\n\n return qml.expval(H) # we are interested in minimizing this expectation value\n\n # initialize the gate parameters\n params = np.zeros(3)\n\n # initialize with converged parameters from previous point\n if pes_point > 0:\n params = params_old\n\n prev_energy = 0.0\n for n in range(50):\n # perform optimization step\n params, energy = opt.step_and_cost(circuit, params)\n\n if np.abs(energy - prev_energy) < 1e-6:\n break\n prev_energy = energy\n\n # store the converged parameters\n params_old = params\n pes_point = pes_point + 1\n\n energies.append(energy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's plot the results \ud83d\udcc8\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n\nfig, ax = plt.subplots()\nax.plot(r_range, energies)\n\nax.set(\n xlabel=\"Bond length (Bohr)\",\n ylabel=\"Total energy (Hartree)\",\n title=\"Potential energy surface for H$_2$ dissociation\",\n)\nax.grid()\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the potential energy surface for the dissociation of a hydrogen\nmolecule into two hydrogen atoms. It is a numerical calculation of the\nsame type of plot that was illustrated in the beginning. In a diatomic\nmolecule such as $H_2$, it can be used to obtain the equilibrium bond\nlength --- the distance between the two atoms that minimizes the total\nelectronic energy. This is simply the minimum of the curve. We can also\nobtain the bond dissociation energy, which is the difference in the\nenergy of the system when the atoms are far apart and the energy at\nequilibrium. At sufficiently large separations, the atoms no longer form\na molecule, and the system is called \"dissociated\".\n\nLet's use our results to compute the equilibrium bond length and the\nbond dissociation energy:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# equilibrium energy\ne_eq = min(energies)\n# energy when atoms are far apart\ne_dis = energies[-1]\n\n# Bond dissociation energy\nbond_energy = e_dis - e_eq\n\n# Equilibrium bond length\nidx = energies.index(e_eq)\nbond_length = r_range[idx]\n\nprint(f\"The equilibrium bond length is {bond_length:.1f} Bohrs\")\nprint(f\"The bond dissociation energy is {bond_energy:.6f} Hartrees\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These estimates can be improved by using bigger basis sets and\nextrapolating to the complete basis set limit \\[\\#motta2020\\]\\_. The\ncalculations are of course are subject to the grid size of interatomic\ndistances considered. The finer the grid size, the better the estimates.\n\n> **note**\n>\n> Did you notice a trick we used to speed up the calculations? The\n> converged gate parameters for a particular geometry on the PES are\n> used as the initial guess for the calculation at the adjacent\n> geometry. With a better guess, the algorithm converges faster and we\n> save considerable time.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hydrogen Exchange Reaction\n==========================\n\nAfter studying a simple diatomic bond dissociation, we move to a\nslightly more complicated case: a hydrogen exchange reaction.\n\n$$H_2 + H \\rightarrow H + H_2.$$\n\nThis reaction has a barrier, the [transition\nstate](https://en.wikipedia.org/wiki/Transition_state), that must be\ncrossed for the exchange of an atom to be complete. In this case, the\ntransition state corresponds to a specific linear arrangement of the\natoms where one $H-H$ bond is partially broken and the other $H-H$ bond\nis partially formed. The molecular movie \u269b\ufe0f\ud83c\udfa5 below is an illustration of\nthe reaction trajectory. It depicts how the distance between the\nhydrogen atoms changes as one bond is broken and another one is formed.\nThe path along which the reaction proceeds is known as the [reaction\ncoordinate](https://en.wikipedia.org/wiki/Reaction_coordinate).\n\n![](/demonstrations/vqe_bond_dissociation/h3_mol_movie.gif)\n\n> width\n>\n> : 50%\n>\n> align\n>\n> : center\n>\nIn a minimal basis like STO-3G, this system consists of three electrons\nin six spin molecular orbitals. This translates into a six-qubit\nproblem, for which the Hartree-Fock state is $|111000\\rangle$. As there\nis an unpaired electron, the spin multiplicity is equal to two and needs\nto be specified, since it differs from the default value of one.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"symbols = [\"H\", \"H\", \"H\"]\nmultiplicity = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To build a potential energy surface for the hydrogen exchange, we fix\nthe positions of the outermost atoms, and change only the placement of\nthe middle atom. For this circuit, we employ all single and double\nexcitation gates, which can be conveniently done with the\n\\~.pennylane.templates.subroutines.AllSinglesDoubles template. The rest\nof the procedure follows as before.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from pennylane.templates import AllSinglesDoubles\n\nenergies = []\npes_point = 0\n\n# get all the singles and doubles excitations, and Hartree-Fock state\nelectrons = 3\norbitals = 6\nsingles, doubles = qchem.excitations(electrons, orbitals)\nhf = qml.qchem.hf_state(electrons, orbitals)\n\n\n# loop to change reaction coordinate\nr_range = np.arange(1.0, 3.0, 0.1)\nfor r in r_range:\n\n coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, r, 0.0, 0.0, 4.0])\n\n # We now specify the multiplicity\n H, qubits = qchem.molecular_hamiltonian(symbols, coordinates, mult=multiplicity)\n\n dev = qml.device(\"default.qubit\", wires=qubits)\n opt = qml.GradientDescentOptimizer(stepsize=1.5)\n\n @qml.qnode(dev)\n def circuit(parameters):\n AllSinglesDoubles(parameters, range(qubits), hf, singles, doubles)\n return qml.expval(H) # we are interested in minimizing this expectation value\n\n params = np.zeros(len(singles) + len(doubles))\n\n if pes_point > 0:\n params = params_old\n\n prev_energy = 0.0\n\n for n in range(60):\n params, energy = opt.step_and_cost(circuit, params)\n if np.abs(energy - prev_energy) < 1e-6:\n break\n prev_energy = energy\n\n # store the converged parameters\n params_old = params\n pes_point = pes_point + 1\n\n energies.append(energy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Once the calculation is complete, we can plot the resulting potential\nenergy surface.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\nax.plot(r_range, energies)\n\nax.set(\n xlabel=\"Distance (Bohr)\",\n ylabel=\"Total energy (Hartree)\",\n)\nax.grid()\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The two minima in the curve represent the energy of the reactants and\nproducts. The transition state is represented by the local maximum.\nThese are the configurations illustrated in the animation above.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Activation energy barriers and reaction rates\n=============================================\n\nThe potential energy surfaces we computed so far can be leveraged for\nother important tasks, such as computing activation energy barriers and\nreaction rates. The activation energy barrier ( $E_{a}$) is defined as\nthe difference between the energy of the reactants and the energy of the\ntransition state.\n\n$$E_{a} = E_{TS} - E_{R}.$$\n\nThis can be computed from the potential energy surface:\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Energy of the reactants and products - two minima on the PES\ne_eq1 = min(energies)\ne_eq2 = min([x for x in energies if x != e_eq1])\n\nidx1 = energies.index(e_eq1)\nidx2 = energies.index(e_eq2)\n\n# Transition state is the local maximum between reactant and products\nidx_min = min(idx1, idx2)\nidx_max = max(idx1, idx2)\n\n# Transition state energy\nenergy_ts = max(energies[idx_min:idx_max])\n\n# Activation energy\nactivation_energy = energy_ts - e_eq1\n\nprint(f\"The activation energy is {activation_energy:.6f} Hartrees\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The reaction rate constant ($k$) has an exponential dependence on the\nactivation energy barrier, as shown in the [Arrhenius\nequation](https://en.wikipedia.org/wiki/Arrhenius_equation) (Arrr!\n\ud83c\udff4\u200d\u2620\ufe0f):\n\n$$k = Ae^{-{E_{a}}/k_BT},$$\n\nwhere $k_B$ is the Boltzmann constant, $T$ is the temperature, and $A$\nis a pre-exponential factor that can be determined empirically for each\nreaction. Crucially, the rate at which a chemical reaction occurs\ndepends exponentially on the activation energy computed from the PES ---\nthis is a good reminder of the importance of performing highly-accurate\ncalculations in quantum chemistry!\n\nFor example, let's calculate the ratio of reaction rates when the\ntemperature is doubled. We have\n\n$$\\frac{k_2}{k_1}=\\frac{Ae^{-{E_{a}}/k_B(2T)}}{Ae^{-{E_{a}}/k_BT}}=e^{E_a/2k_B T}.$$\n\nWe choose $T=300$ Kelvin, which is essentially room temperature. This\nmakes doubling the temperature roughly equivalent to the temperature\ninside a pizza oven. We have\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# convert to joules\nactivation_energy *= 4.36e-18\n# Boltzmann constant in Joules/Kelvin\nk_B = 1.38e-23\n# Temperature\nT = 300\n\nratio = np.exp(activation_energy / (2 * k_B * T))\n\nprint(f\"Ratio of reaction rates is {ratio:.0f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Doubling the temperature can increase the rate by a factor of almost two\nmillion! For a similar reason, changing the activation energy can lead\nto drastic changes in the reaction rates, which means we have to be\ncareful to compute it very accurately.\n\nConclusion\n==========\n\nWe can learn how atoms combine to form different molecules by performing\nexperiments; this is the approach many of us learn as children by\nplaying with chemistry sets. However, a deeper quantitative\nunderstanding of chemical reactions can be achieved by performing\ntheoretical simulations of the mechanisms for forming and breaking\nbonds. This tutorial described how simple chemical reactions can be\nsimulated using quantum algorithms that reconstruct potential energy\nsurfaces, allowing us to identify reactants and products as minima of\nthe energy, and transition states as local maxima. These results can\nthen be used to calculate activation energies and reaction rates. The\ngoal (and challenge!) for quantum computing is to improve both hardware\nand algorithms to reach the regime where providing accurate simulations\nbecomes intractable for existing methods. If successful, this quest will\nallow us to understand the properties of quantum systems in ways that\nhave so far been out of reach.\n\nReferences\n==========\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}