{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction\n", "In this notebook, we illustrate how to implement the Variational Quantum Eigensolver (VQE) algorithm in Amazon Braket SDK to compute the potential energy surface (PES) for the Hydrogen molecule. Specifically, we illustrate the following features of Amazon Braket SDK:\n", "* `LocalSimulator` which allows one to simulate quantum circuits on their local machine\n", "* Construction of the ansatz circuit for VQE in Braket SDK\n", "* Computing expectation values of the individual terms in the Hamiltonian in Braket SDK" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem formulation\n", "In computational chemistry, we are typically interested in finding the ground-state energy (i.e., minimum energy) of a molecule for a given configuration of atomic positions.\n", "This is accomplished by finding the lowest eigenvalues and eigenstates of the molecular Hamiltonian.\n", "## Molecular Hamiltonian\n", "In [atomic units](https://en.wikipedia.org/wiki/Hartree_atomic_units), the full molecular Hamiltonian takes the form [1],\n", "\n", "$$\n", "H_{full} = - \\sum_i \\frac{\\nabla_i^2}{2} - \\sum_I \\frac{\\nabla_I^2}{2 M_I}\n", " - \\sum_{i,I} \\frac{Z_I}{|\\mathbf{r}_i - \\mathbf{R}_I|}\n", " + \\frac{1}{2} \\sum_{i \\neq j} \\frac{1}{|\\mathbf{r}_i - \\mathbf{r}_j|}\n", " + \\frac{1}{2} \\sum_{I \\neq J} \\frac{1}{|\\mathbf{R}_I - \\mathbf{R}_J|}\n", "\\tag{1}\n", "$$\n", "\n", "where $i$ and $I$ label the electrons and nuclei respectively.\n", "Here $\\mathbf{r}_i$ and $\\mathbf{R}_I$ denote the electron and nucleus positions, $Z_I$ denotes nuclear charge, and $M_I$ denotes nuclear mass respectively.\n", "The first two terms in Eq.(1) represent the electronic and nuclear kinetic energy, while the remaining terms denote the Coulombic interaction between the nuclei and electrons, amongst the electrons, and amongst the nuclei respectively.\n", "\n", "Since we are mostly interested in the electronic structure of the molecule, we invoke the [Born-Oppenheimer approximation](https://en.wikipedia.org/wiki/Born%E2%80%93Oppenheimer_approximation), which is based on the fact that nuclei are over a thousand times heavier than electrons.\n", "This approximation allows us to reduce the problem to finding the eigenvalues of the following electronic Hamiltonian for a given nuclear configuration [1], i.e., fixed set of nuclear positions $\\{\\mathbf{R}_I\\}$,\n", "\n", "$$\n", "H = - \\sum_i \\frac{\\nabla_i^2}{2}\n", " - \\sum_{i,I} \\frac{Z_I}{|\\mathbf{r}_i - \\mathbf{R}_I|}\n", " + \\frac{1}{2} \\sum_{i \\neq j} \\frac{1}{|\\mathbf{r}_i - \\mathbf{r}_j|}\n", "\\tag{2}\n", "$$\n", "\n", "## Second quantization\n", "In order to determine the eigenvalues and eigenstates of the electronic Hamiltonian of a molecule given in Eq. (2), we start with a set of single-electron basis functions (also called *spin orbitals*), $\\{\\phi_p(\\mathbf{x}_i)\\}$, where $\\mathbf{x}_i=(\\mathbf{r}_i, s_i)$ denotes the spatial position and spin of the $i$-th electron.\n", "The many-electron wavefunction is expressed as a [Slater determinant](https://en.wikipedia.org/wiki/Slater_determinant) of these basis functions, i.e., as an antisymmetrized product of the single-electron basis functions, so that the electrons automatically satisfy the [Pauli exclusion principle](https://en.wikipedia.org/wiki/Pauli_exclusion_principle).\n", "While the total number of single-electron basis functions considered $M$, is typically larger than the total number of electrons $N$, in the molecule, the $N$ electrons can only occupy $N$ ($< M$) of these orbitals in a Slater determinant, i.e., a Slater determinant contains only $N$ occupied orbitals [1].\n", "Since any Slater determinant of the basis functions is uniquely determined by which orbitals are occupied, a Slater determinant can be represented in an abstract vector space, called the *Fock space*, by an *occupation number* vector $|f\\rangle = |f_{M-1}, \\ldots, f_0 \\rangle$ where each $f_i$ is either 0 or 1 depending on whether the $i$-th orbitals is unoccupied or occupied respectively in the Slater determinant.\n", "\n", "Next, we introduce a set of operators $\\{a_p, a_p^\\dagger\\}$ , called the fermion *annihilation* and *creation* operators respectively, corresponding to each of the basis functions $\\{\\phi_p(\\mathbf{x}_i)\\}$, which satisfy the anticommutation relations,\n", "\n", "$$\n", "\\{a_p, a_q\\} = \\{a_p^\\dagger, a_q^\\dagger\\} = 0, \\qquad \\{a_p, a_q^\\dagger\\} = \\delta_{pq}\n", "\\tag{3}\n", "$$\n", "\n", "where the anticommutator of two operators $C$ and $D$ is defined as $\\{C,D\\} = CD + DC$ and $\\delta_{cd}$ denotes the Kronecker delta.\n", "\n", "Expressing the electronic Hamiltonian in Eq. (2) in terms of these fermionic operators we get [1],\n", "\n", "$$\n", "H = \\sum_{p,q} h_{pq} a_p^\\dagger a_q \n", " +\\frac{1}{2} \\sum_{p,q,r,s} h_{pqrs} a_p^\\dagger a_q^\\dagger a_r a_s\n", "\\tag{4}\n", "$$\n", "\n", "where,\n", "\n", "$$\n", "h_{pq} = \\int d\\mathbf{x}\\ \\phi_p^*(\\mathbf{x}) \\left( -\\frac{\\nabla^2}{2} - \\sum_I \\frac{Z_I}{|\\mathbf{r} - \\mathbf{R}_I|} \\right) \\phi_q(\\mathbf{x})\n", "\\tag{5}\n", "$$\n", "\n", "and\n", "\n", "$$\n", "h_{pqrs} = \\int d\\mathbf{x}_1 d\\mathbf{x}_2\\ \\frac{\\phi_p^*(\\mathbf{x}_1) \\phi_q^*(\\mathbf{x}_2)\\phi_s(\\mathbf{x}_1)\\phi_r(\\mathbf{x}_2)}{|\\mathbf{r}_1 - \\mathbf{r}_2|}\n", "\\tag{6}\n", "$$\n", "\n", "### Basis sets\n", "A variety of basis sets are used in computational chemistry depending on the chemical system being studied and the accuracy required, e.g.,\n", "* [STO-$n$G][r8] (also called *minimal basis sets*)\n", "* [Split-valence][r9]\n", "* [Correlation-consistent][r9]\n", "* [Polarization-consistent][r9]\n", "* [Plane-wave][r9]\n", "\n", "[r8]: https://en.wikipedia.org/wiki/STO-nG_basis_sets\n", "\n", "[r9]: https://en.wikipedia.org/wiki/Basis_set_(chemistry) \n", "\n", "## Hartree-Fock method\n", "The [Hartree-Fock (HF)](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) method is the starting point for most computational chemistry calculations of atoms and molecules [1].\n", "The HF method aims to find the single-most dominant Slater determinant that best approximates the system wavefunction by optimizing the spatial form of the spin orbitals in a self-consistent fashion (hence, also called *self-consistent field* or SCF method) to minimize the energy of the wavefunction.\n", "The Slater determinant generated by the optimized spin orbitals (also known as *canonical orbitals*) computed from the HF method is used as the reference state for post Hartree-Fock methods that try to correct for some of the approximations in the HF method.\n", "### Post Hartree-Fock methods\n", "There are a variety of post Hartree-Fock (post-HF) methods employed in classical computational chemistry to improve upon the accuracy of a simple HF computation to varying degrees [1], e.g.,\n", "* [Configuration interaction and Full Configuration Interaction](https://en.wikipedia.org/wiki/Configuration_interaction)\n", "* [Multi-configurational self-consistent field](https://en.wikipedia.org/wiki/Multi-configurational_self-consistent_field)\n", "* [Coupled cluster](https://en.wikipedia.org/wiki/Coupled_cluster)\n", "* [Perturbation theory](https://en.wikipedia.org/wiki/M%C3%B8ller%E2%80%93Plesset_perturbation_theory)\n", "\n", "## Mapping to quantum computer\n", "In order to represent the electronic Hamiltonian of the molecular system under study, we need to map it to operators that act on the qubits of a quantum computer.\n", "There are various encoding techniques to accomplish this, e.g., parity encoding, Bravyi-Kitaev encoding, etc., but the simplest one is the Jordan-Wigner encoding [2].\n", "### Jordan-Wigner encoding\n", "In the Jordan-Wigner (JW) encoding, we store the occupation number of an orbital $f_i$ in the $i$-th qubit $q_i$ as either $|0\\rangle$ or $|1\\rangle$ depending on whether the orbital is unoccupied or occupied respectively.\n", "Correspondingly, the fermionic annihilation and creation operators are mapped to the qubit operators as\n", "\n", "$$\n", "a_p = Q_p \\otimes Z_{p-1} \\otimes \\ldots \\otimes Z_0\n", "\\tag{7}\n", "$$\n", "\n", "and\n", "\n", "$$\n", "a_p^\\dagger = Q_p^\\dagger \\otimes Z_{p-1} \\otimes \\ldots \\otimes Z_0\n", "\\tag{8}\n", "$$\n", "\n", "where $Q = |0\\rangle \\langle 1|, Q^\\dagger = |1\\rangle \\langle 0|$, $Z=|0\\rangle \\langle 0| - |1\\rangle \\langle 1| $ denotes the Pauli-Z matrix and $\\otimes$ denotes the tensor product.\n", "This shows that in the JW encoding, while the occupation of an orbital is stored locally, i.e., in a single qubit, the parity is stored non-locally.\n", "This makes JW encoding not as efficient compared to other encoding techniques, e.g., Bravyi-Kitaev encoding, but for small systems with a few qubits, the efficiency gap is not significant.\n", "\n", "## Variational Quantum Eigensolver\n", "One of the most promising algorithms for performing quantum chemistry computations on the NISQ (noisy intermediate-scale quantum) era quantum computers we can currently use is the variational quantum eigensolver (VQE) [2].\n", "VQE is a hybrid quantum-classical algorithm using the quantum computer only for a classically intractable subroutine and has been shown to be robust against noise [2,4], and capable of finding ground-state energies of small molecules using low-depth quantum circuits.\n", "\n", "VQE is based on the Raleigh-Ritz variational principle: for a trial wavefunction $|\\psi(\\vec{\\theta})\\rangle$ parameterized by a vector parameter $\\vec{\\theta}$ for a system with Hamiltonian $H$, whose lowest energy eigenvalue is $E_0$ [2],\n", "\n", "$$\n", "\\langle \\psi(\\vec{\\theta}) | H | \\psi(\\vec{\\theta}) \\rangle \\geq E_0\n", "\\tag{9}\n", "$$\n", "\n", "Thus, VQE reduces the problem of finding the ground-state energy of the system Hamiltonian to finding the optimum value of the vector parameter $\\vec{\\theta}$.\n", "The parameterized wavefunction $|\\psi(\\vec{\\theta})\\rangle$ is prepared on a quantum computer by applying a series of parameterized unitary gates (also called a *parameterized circuit*) to the qubits initialized in some reference state $|\\psi_{ref}\\rangle$, e.g., a HF state, as [2] \n", "\n", "$$\n", "|\\psi(\\vec{\\theta})\\rangle = U(\\vec{\\theta}) | \\psi_{ref} \\rangle\n", "\\tag{10}\n", "$$\n", "\n", "As classical computers are unable to efficiently prepare, store and measure the expectation value of the Hamiltonian $H$ for the parameterized \n", "wavefunction $|\\psi(\\vec{\\theta})\\rangle$, we use the quantum computer for this subroutine, and then use the classical computer to update the parameters using an optimization algorithm.\n", "\n", "To complete the specification of VQE, one needs to specify the form of the parameterized wavefunction $|\\psi(\\vec{\\theta})\\rangle$, known as the *ansatz*. A popular chemically-inspired ansatz for VQE is the unitary coupled cluster ansatz truncated to singles and doubles excitations (UCCSD).\n", "### UCCSD ansatz\n", "The unitary coupled cluster (UCC) method is an extension of the coupled-cluster (CC) method, which is one of the most popular post-HF methods [2].\n", "In the UCC method, the parameterized trial function is given by\n", "\n", "$$\n", "U(\\vec{\\theta}) = e^{T - T^\\dagger}\n", "\\tag{11}\n", "$$\n", "\n", "where $T = \\sum_i T_i$ is the sum of various excitation levels, e.g.,\n", "\n", "$$\n", "T_1 = \\sum_{i \\in virt, \\alpha \\in occ} \\theta_{i \\alpha} a_i^\\dagger a_\\alpha, \\qquad\n", "T_2 = \\sum_{i, j \\in virt, \\alpha, \\beta \\in occ} \\theta_{i j \\alpha \\beta} a_i^\\dagger a_j^\\dagger a_\\alpha a_\\beta, \\ldots\n", "\\tag{12}\n", "$$\n", "\n", "denote the single and double excitation levels respectively.\n", "\n", "The UCC method possesses all the advantages of the CC method and is also fully variational as well as able to converge when using multi-reference states. As in CC, we truncate the expression for $T$ at a given excitation level -- usually single and double excitations leading to the UCCSD ansatz [2]." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Hydrogen molecule\n", "For simplicity, we use the minimal basis set STO-3G for $H_2$ molecule, which includes only the $\\{1s\\}$ orbital for each $H$ atom.\n", "Since each $H$ atom contributes one spin-orbital, and there are 2 possible spins for each orbital (up or down), this leads to a total of 4 orbitals for $H_2$ in the STO-3G basis set [2],\n", "$$\n", "| 1s_{A\\uparrow} \\rangle, | 1s_{A\\downarrow} \\rangle, | 1s_{B\\uparrow} \\rangle, | 1s_{B\\downarrow} \\rangle\n", "\\tag{13}\n", "$$\n", "where the subscripts A and B label the $H$ atom and $\\uparrow/\\downarrow$ label the electron spin respectively.\n", "\n", "In the molecular orbital basis for $H_2$, given by [2]\n", "$$\n", "| \\sigma_{g \\uparrow} \\rangle = \\frac{1}{\\sqrt{2}} ( | 1s_{A\\uparrow} \\rangle + | 1s_{B\\uparrow} \\rangle) \\\\\n", "| \\sigma_{g \\downarrow} \\rangle = \\frac{1}{\\sqrt{2}} ( | 1s_{A\\downarrow} \\rangle + | 1s_{B\\downarrow} \\rangle) \\\\\n", "| \\sigma_{u \\uparrow} \\rangle = \\frac{1}{\\sqrt{2}} ( | 1s_{A\\uparrow} \\rangle - | 1s_{B\\uparrow} \\rangle) \\\\\n", "| \\sigma_{u \\downarrow} \\rangle = \\frac{1}{\\sqrt{2}} ( | 1s_{A\\downarrow} \\rangle - | 1s_{B\\downarrow} \\rangle)\n", "\\tag{14}\n", "$$\n", "the Slater determinant can be expressed in the occupation number basis as\n", "$$\n", "| \\psi \\rangle = | f_{u\\downarrow}, f_{u\\uparrow}, f_{g\\downarrow}, f_{g\\uparrow} \\rangle = | f_3, f_2, f_1, f_0 \\rangle\n", "\\tag{15}\n", "$$\n", "where each $f_i \\in \\{0,1\\}$.\n", "\n", "The second-quantized electronic Hamiltonian in this basis is given by [2]\n", "$$\n", "H = H_1 + H_2 + H_3 + H_4\n", "\\tag{16}\n", "$$\n", "\n", "where\n", "\n", "$$\n", "H_1 = h_{00} a_0^\\dagger a_0 + h_{11} a_1^\\dagger a_1 + h_{22} a_2^\\dagger a_2 + h_{33} a_3^\\dagger a_3,\n", "\\tag{17}\n", "$$\n", "\n", "$$\n", "H_2 = h_{0110} a_0^\\dagger a_1^\\dagger a_1 a_0 + h_{2332} a_2^\\dagger a_3^\\dagger a_3 a_2 + h_{0330} a_0^\\dagger a_3^\\dagger a_3 a_0 + h_{1221} a_1^\\dagger a_2^\\dagger a_2 a_1,\n", "\\tag{18}\n", "$$\n", "\n", "$$\n", "H_3 = (h_{0220}-h_{0202}) a_0^\\dagger a_2^\\dagger a_2 a_0 + (h_{1331} - h_{1313}) a_1^\\dagger a_3^\\dagger a_3 a_1,\n", "\\tag{19}\n", "$$\n", "\n", "and\n", "\n", "$$\n", "H_4 = h_{0132} (a_0^\\dagger a_1^\\dagger a_3 a_2 + a_2^\\dagger a_3^\\dagger a_1 a_0) + h_{0312} (a_0^\\dagger a_3^\\dagger a_1 a_2 + a_2^\\dagger a_1^\\dagger a_3 a_0)\n", "\\tag{20}\n", "$$\n", "\n", "where the coefficients $h$ are given by Eqs. (5) and (6).\n", "\n", "Using the JW encoding, the electronic Hamiltonian for $H_2$ in Eq.(16) can be expressed as a 4-qubit operator [2]\n", "\n", "$$\n", "H_Q = H_{Q1} + H_{Q2} + H_{Q3}\n", "\\tag{21}\n", "$$\n", "\n", "where\n", "\n", "$$\n", "H_{Q1} = h_0 I + h_1 Z_0 + h_2 Z_1 + h_3 Z_2 + h_4 Z_3,\n", "\\tag{22}\n", "$$\n", "\n", "$$\n", "H_{Q2} = h_5 Z_0 Z_1 + h_6 Z_0 Z_2 + h_7 Z_1 Z_2 + h_8 Z_0 Z_3 + h_9 Z_1 Z_3 + h_{10} Z_2 Z_3,\n", "\\tag{23}\n", "$$\n", "\n", "and\n", "$$\n", "H_{Q3} = h_{11} Y_0 Y_1 X_2 X_3 + h_{12} X_0 Y_1 Y_2 X_3 + h_{13} Y_0 X_1 X_2 Y_3 + h_{14} X_0 X_1 Y_2 Y_3\n", "\\tag{24}\n", "$$\n", "\n", "In the JW encoding, the HF state, which is taken as the reference state for the UCCSD ansatz for VQE, is given by,\n", "\n", "$$\n", "| \\psi_{HF} \\rangle = | 0011 \\rangle\n", "\\tag{25}\n", "$$\n", "\n", "which represents the state in which the molecular orbitals $| \\sigma_{g \\uparrow} \\rangle$ and $| \\sigma_{g \\downarrow} \\rangle$ are occupied, while the molecular orbitals $| \\sigma_{u \\uparrow} \\rangle$ $| \\sigma_{u \\downarrow} \\rangle$ are unoccupied.\n", "The most general state for $H_2$ with the same charge and spin multiplicity as the HF state can be expressed as [2],\n", "\n", "$$\n", "| \\psi \\rangle = \\alpha | 0011 \\rangle + \\beta | 1100 \\rangle + \\gamma | 1001 \\rangle + \\delta | 0110 \\rangle\n", "\\tag{26}\n", "$$\n", "\n", "Finally, the UCCSD ansatz for $H_2$ in the STO-3G basis can be simplified to the following single-parameter unitary operator acting on the HF reference state given by Eq. (25) [2]\n", "\n", "$$\n", "U (\\theta) = e^{\\imath \\theta X_3 X_2 X_1 Y_0}\n", "\\tag{27}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Implementation details\n", "We use the following procedure to determine the ground-state energy of the $H_2$ molecule at various bond-lengths, i.e., its potential energy surface (PES).\n", "1. We use OpenFermion-PySCF to compute the coefficients of the different terms of the 4-qubit operator in Eq.(21) representing the electronic Hamiltonian for $H_2$ at the required bond lengths of $H_2$ molecule.\n", "1. Using Amazon Braket, we implement the UCCSD ansatz given by Eq.(27) through the quantum circuit shown below [2].\n", "1. We run the above quantum circuit using the local simulator of Amazon Braket for various values of the parameter in the UCCSD ansatz to determine its optimal value and using this value, compute the expectation value of the qubit operator in Eq.(21) to get the estimated ground-state energy of $H_2$ for each value of bond-length.\n", "\n", "Finally, we assess the accuracy of the ground-state energy of $H_2$ molecule for various bond-lengths computed by VQE to that computed using the classical chemistry post-HF method FCI (full configuration interaction), which is the most accurate computation that can be performed for a given basis set.\n", "\n", "
\n", "\n", "## Pre-requisites\n", "* Install [Braket Python SDK](https://github.com/aws/amazon-braket-sdk-python) by following the instruction provided in the repository." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Execution\n", "## Imports and setup" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Name: amazon-braket-sdk\r\n", "Version: 1.34.4.dev0\r\n", "Summary: An open source library for interacting with quantum computing devices on Amazon Braket\r\n", "Home-page: https://github.com/aws/amazon-braket-sdk-python\r\n", "Author: Amazon Web Services\r\n", "Author-email: \r\n", "License: Apache License 2.0\r\n", "Location: /home/ec2-user/anaconda3/envs/Braket/lib/python3.7/site-packages\r\n", "Requires: amazon-braket-default-simulator, amazon-braket-schemas, backoff, boltons, boto3, nest-asyncio, networkx, numpy, openpulse, openqasm3, oqpy, sympy\r\n", "Required-by: amazon-braket-algorithm-library, amazon-braket-pennylane-plugin, amazon-braket-strawberryfields-plugin\r\n" ] } ], "source": [ "# create a directory named \"data\" to store intermediate classical computation results from OpenFermion\n", "!mkdir -p \"data\"\n", "!pip show amazon-braket-sdk" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# import required libraries\n", "import os\n", "import time\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from openfermion import MolecularData\n", "from openfermion.transforms import get_fermion_operator, jordan_wigner\n", "from openfermionpyscf import run_pyscf\n", "from braket.circuits import Circuit, FreeParameter, observables\n", "from braket.devices import LocalSimulator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parameters defining the problem space" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INFO: Simulating ground-state of H2 molecule for 28 bond lengths\n" ] } ], "source": [ "# Initialize variable to determine total run-time of the program\n", "run_time = time.time()\n", "# Set parameters that won't change for all simulations\n", "# total charge\n", "tot_chg = 0\n", "# set the multiplicity\n", "# this parameter specifies how electrons are paired in orbitals\n", "# spin multiplicity is equal to the number of unpaired electrons plus one\n", "# e.g., spin multiplicity = 1 implies no unpaired electron, i.e., for every\n", "# spin-up electron in a spatial orbital there is a spin-down electron,\n", "# spin-multiplicity = 3 implies 2 unpaired electrons, and so on.\n", "# See https://en.wikipedia.org/wiki/Multiplicity_(chemistry) for more details\n", "spin_mult = 1\n", "# use the minimal basis set (STO-3G)\n", "basis_set = 'sto-3g'\n", "# number of active electrons and orbitals to consider\n", "# For H2 as we're consider all electrons and orbitals as active for STO-3G\n", "occ_ind = act_ind = None\n", "# bond lengths to simulate ground-state of H2 for (in Angstroms)\n", "bond_lengths = np.arange(start=0.24, stop=3.00, step=0.1)\n", "# no. of bond lengths to simulate H2 for\n", "n_configs = len(bond_lengths)\n", "print(f\"INFO: Simulating ground-state of H2 molecule for {n_configs} bond lengths\")\n", "# no. of VQE parameter values to scan\n", "n_theta = 24\n", "# no. of shots for measurements on quantum device\n", "n_shots = 2000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Classical pre-computation of electronic Hamiltonian of $H_2$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INFO: Computing Hamiltonian for bond-length 0.24 A...\n", "INFO: Computing Hamiltonian for bond-length 0.34 A...\n", "INFO: Computing Hamiltonian for bond-length 0.44 A...\n", "INFO: Computing Hamiltonian for bond-length 0.54 A...\n", "INFO: Computing Hamiltonian for bond-length 0.64 A...\n", "INFO: Computing Hamiltonian for bond-length 0.74 A...\n", "INFO: Computing Hamiltonian for bond-length 0.84 A...\n", "INFO: Computing Hamiltonian for bond-length 0.94 A...\n", "INFO: Computing Hamiltonian for bond-length 1.04 A...\n", "INFO: Computing Hamiltonian for bond-length 1.14 A...\n", "INFO: Computing Hamiltonian for bond-length 1.24 A...\n", "INFO: Computing Hamiltonian for bond-length 1.34 A...\n", "INFO: Computing Hamiltonian for bond-length 1.44 A...\n", "INFO: Computing Hamiltonian for bond-length 1.54 A...\n", "INFO: Computing Hamiltonian for bond-length 1.64 A...\n", "INFO: Computing Hamiltonian for bond-length 1.74 A...\n", "INFO: Computing Hamiltonian for bond-length 1.84 A...\n", "INFO: Computing Hamiltonian for bond-length 1.94 A...\n", "INFO: Computing Hamiltonian for bond-length 2.04 A...\n", "INFO: Computing Hamiltonian for bond-length 2.14 A...\n", "INFO: Computing Hamiltonian for bond-length 2.24 A...\n", "INFO: Computing Hamiltonian for bond-length 2.34 A...\n", "INFO: Computing Hamiltonian for bond-length 2.44 A...\n", "INFO: Computing Hamiltonian for bond-length 2.54 A...\n", "INFO: Computing Hamiltonian for bond-length 2.64 A...\n", "INFO: Computing Hamiltonian for bond-length 2.74 A...\n", "INFO: Computing Hamiltonian for bond-length 2.84 A...\n", "INFO: Computing Hamiltonian for bond-length 2.94 A...\n", "INFO: Computed Hamiltonians for 28 bond-lengths.\n" ] } ], "source": [ "# Dictionary to store molecular data and qubit Hamiltonians for each\n", "# molecular configuration(i.e., bond length)\n", "mol_configs = {}\n", "# Store molecular config data and Hamiltonian for each bond length\n", "for rr in bond_lengths:\n", " # round to 2 digits to get a reasonable-length number for bond-length\n", " r = round(rr,2)\n", " print(f\"INFO: Computing Hamiltonian for bond-length {r} A...\")\n", " geom = [('H', (0., 0., -r/2.)), ('H', (0., 0., r/2.))]\n", " # Make sure a directory named 'data' exists in the current folder,\n", " # else this statement will throw an error!\n", " h2_molecule = MolecularData(\n", " geometry=geom, basis=basis_set, multiplicity=spin_mult,\n", " description='bondlength_'+str(r)+'A', filename=\"\",\n", " data_directory=os.getcwd()+'/data')\n", " # Run PySCF to get molecular integrals, HF and FCI energies\n", " h2_molecule = run_pyscf(molecule=h2_molecule, run_scf=True,\n", " run_mp2=False, run_cisd=False,\n", " run_ccsd=False, run_fci=True, verbose=False)\n", " # Convert electronic Hamiltonian to qubit operator using JW encoding\n", " h2_qubit_hamiltonian = jordan_wigner(get_fermion_operator(\n", " h2_molecule.get_molecular_hamiltonian(occupied_indices=occ_ind,\n", " active_indices=act_ind)))\n", " # store molecular data and qubit operator for this config in an ordered list\n", " mol_configs[r] = [h2_molecule, h2_qubit_hamiltonian]\n", "print(f\"INFO: Computed Hamiltonians for {len(mol_configs)} bond-lengths.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Construction of UCCSD ansatz circuit in Braket" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "T : |0| 1 |2|3|4| 5 |6|7|8| 9 |\n", " \n", "q0 : -H--------------X-Rz(a_theta)-X-H-------------\n", " | | \n", "q1 : -H------------X-C-------------C-X-H-----------\n", " | | \n", "q2 : -X-H--------X-C-----------------C-X-H---------\n", " | | \n", "q3 : -X-Rx(1.57)-C---------------------C-Rx(-1.57)-\n", "\n", "T : |0| 1 |2|3|4| 5 |6|7|8| 9 |\n", "\n", "Unassigned parameters: [a_theta].\n", "Total number of qubits in the circuit: 4\n", "{'0011': 1.0}\n", "{'0011': 1.0}\n" ] } ], "source": [ "# Construct circuit for UCCSD ansatz parameterized by VQE parameter per McArdle et al. \n", "a_theta = FreeParameter(\"a_theta\")\n", "# Initialize HF state |0011>\n", "ansatz_uccsd = Circuit().x(2).x(3)\n", "# Perform intial rotations to measure in Y & X bases\n", "ansatz_uccsd.rx(3,np.pi/2.).h(range(3))\n", "# Entangle with CNOTs\n", "ansatz_uccsd.cnot(3,2).cnot(2,1).cnot(1,0)\n", "# Perform the rotation in Z-basis\n", "ansatz_uccsd.rz(0,a_theta)\n", "# Uncompute the rotations\n", "ansatz_uccsd.cnot(1,0).cnot(2,1).cnot(3,2)\n", "ansatz_uccsd.h(range(3))\n", "ansatz_uccsd.rx(3,-np.pi/2.)\n", "\n", "# initialize quantum device to run VQE over\n", "local_sim = LocalSimulator()\n", "# verify by running the circuit with no rotation and verifying we get HF state\n", "print(ansatz_uccsd)\n", "n_qubits = ansatz_uccsd.qubit_count\n", "print(f\"Total number of qubits in the circuit: {n_qubits}\")\n", "# We can fix the value of a_theta by calling the circuit with the value filled as follows.\n", "# print(local_sim.run(ansatz_uccsd(0)).result().state_vector)\n", "print(local_sim.run(ansatz_uccsd(0), shots=n_shots).result().measurement_probabilities)\n", "# We can also fix the value of a_theta by supplying the inputs argument to run\n", "print(local_sim.run(ansatz_uccsd, shots=n_shots, inputs={'a_theta': 0.0}).result().measurement_probabilities)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computation of expectation value of Hamiltonian operator for $H_2$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def calculate_observable_expectation(a_indices_gates, a_ckt, a_dev, a_shots):\n", " \"\"\" Calculates the expectation value of the given observable\n", " Parameters:\n", " a_indices_gates [sequence(tuple(int, str))]: List of tuple of qubit index & observable\n", " a_ckt [Circuit]: Braket circuit to which measurement gates are to be added\n", " a_dev [Braket Device]: Quantum device to run a_ckt on\n", " a_shots [int]: No. of shots for measuring in bases other than the computational basis\n", " Returns:\n", " float: The expectation value of the observable\n", " \"\"\"\n", " if not a_indices_gates:\n", " # this is the constant term of the Hamiltonian\n", " return 1\n", " factors = {}\n", " for ind, factor in a_indices_gates:\n", " # N.B.: Convert from OpenFermion's little-endian convention to Braket's big-endian convention\n", " qubit = n_qubits - 1 - ind\n", " if factor == \"X\":\n", " factors[qubit] = observables.X()\n", " elif factor == \"Y\":\n", " factors[qubit] = observables.Y()\n", " elif factor == \"Z\":\n", " factors[qubit] = observables.Z()\n", " qubits = sorted(factors)\n", " observable = observables.TensorProduct([factors[qubit] for qubit in qubits])\n", " # initialize measuring circuit and add expectation measurement\n", " measuring_ckt = Circuit().add(a_ckt).expectation(observable=observable, target=qubits)\n", " # compute expectation value\n", " return a_dev.run(measuring_ckt, shots=a_shots).result().values[0]\n", "\n", "def H_exp(a_qH, a_ckt, a_dev, a_shots=n_shots):\n", " \"\"\" Get expectation value of Hamiltonian for a given circuit result.\n", " Parameters:\n", " a_qH [OpenFermion QubitHamiltonian.terms]: Dictionary of OpenFermion QubitHamiltonian operator terms\n", " a_ckt [Braket Circuit]: Circuit to create the final state\n", " a_dev [Braket Device]: Quantum device to run a_ckt on\n", " a_shots [int]: No. of shots for measuring in bases other than the computational basis\n", " Returns:\n", " H_e [float]: Expectation value of a_qH for a_ket\n", " \"\"\"\n", " # initialize expectation value of Hamiltonian\n", " H_e = 0.\n", " # loop over each term in qubit operator\n", " for term in a_qH:\n", " # extract the real-valued coefficient for this term\n", " coeff = np.real(a_qH[term])\n", " # compute and add this term's contribution to the Hamiltonian\n", " H_e += coeff * calculate_observable_expectation(term, a_ckt, a_dev, a_shots)\n", " return H_e" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Results" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INFO: Computing ground-state energy for bond-length 0.24 A...\n", "min = -0.2294 Ha\n", "INFO: Computing ground-state energy for bond-length 0.34 A...\n", "min = -0.7463 Ha\n", "INFO: Computing ground-state energy for bond-length 0.44 A...\n", "min = -0.9747 Ha\n", "INFO: Computing ground-state energy for bond-length 0.54 A...\n", "min = -1.0775 Ha\n", "INFO: Computing ground-state energy for bond-length 0.64 A...\n", "min = -1.1273 Ha\n", "INFO: Computing ground-state energy for bond-length 0.74 A...\n", "min = -1.1400 Ha\n", "INFO: Computing ground-state energy for bond-length 0.84 A...\n", "min = -1.1312 Ha\n", "INFO: Computing ground-state energy for bond-length 0.94 A...\n", "min = -1.1117 Ha\n", "INFO: Computing ground-state energy for bond-length 1.04 A...\n", "min = -1.0907 Ha\n", "INFO: Computing ground-state energy for bond-length 1.14 A...\n", "min = -1.0732 Ha\n", "INFO: Computing ground-state energy for bond-length 1.24 A...\n", "min = -1.0496 Ha\n", "INFO: Computing ground-state energy for bond-length 1.34 A...\n", "min = -1.0321 Ha\n", "INFO: Computing ground-state energy for bond-length 1.44 A...\n", "min = -1.0083 Ha\n", "INFO: Computing ground-state energy for bond-length 1.54 A...\n", "min = -0.9944 Ha\n", "INFO: Computing ground-state energy for bond-length 1.64 A...\n", "min = -0.9820 Ha\n", "INFO: Computing ground-state energy for bond-length 1.74 A...\n", "min = -0.9722 Ha\n", "INFO: Computing ground-state energy for bond-length 1.84 A...\n", "min = -0.9553 Ha\n", "INFO: Computing ground-state energy for bond-length 1.94 A...\n", "min = -0.9538 Ha\n", "INFO: Computing ground-state energy for bond-length 2.04 A...\n", "min = -0.9455 Ha\n", "INFO: Computing ground-state energy for bond-length 2.14 A...\n", "min = -0.9435 Ha\n", "INFO: Computing ground-state energy for bond-length 2.24 A...\n", "min = -0.9404 Ha\n", "INFO: Computing ground-state energy for bond-length 2.34 A...\n", "min = -0.9340 Ha\n", "INFO: Computing ground-state energy for bond-length 2.44 A...\n", "min = -0.9369 Ha\n", "INFO: Computing ground-state energy for bond-length 2.54 A...\n", "min = -0.9354 Ha\n", "INFO: Computing ground-state energy for bond-length 2.64 A...\n", "min = -0.9286 Ha\n", "INFO: Computing ground-state energy for bond-length 2.74 A...\n", "min = -0.9317 Ha\n", "INFO: Computing ground-state energy for bond-length 2.84 A...\n", "min = -0.9339 Ha\n", "INFO: Computing ground-state energy for bond-length 2.94 A...\n", "min = -0.9346 Ha\n" ] } ], "source": [ "# Loop over all bond lengths\n", "for r in mol_configs:\n", " print(f\"INFO: Computing ground-state energy for bond-length {r} A...\")\n", " # intialize min_energy for this config\n", " mol_configs[r].append(np.inf)\n", " # Loop over all VQE parameter values\n", " for theta in np.linspace(start=-np.pi,stop=np.pi,num=n_theta,endpoint=False):\n", " # get expectation value of this config's Hamiltonian for this parameter value\n", " exp_H_theta = H_exp(mol_configs[r][1].terms, ansatz_uccsd(theta), local_sim)\n", " # if this expectation value is less than min found so far, update it\n", " if exp_H_theta < mol_configs[r][2]:\n", " # print(f\"DEBUG: New optimum value of theta = {theta}\")\n", " mol_configs[r][2] = exp_H_theta\n", " print(f\"min = {mol_configs[r][2]:.4f} Ha\")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total running time: 308.69 s\n" ] } ], "source": [ "# Plot ground-state energy for each configuration vs bond length\n", "plt.figure(figsize=(8,6), dpi=700)\n", "plt.plot(list(mol_configs.keys()), [val[0].hf_energy for val in mol_configs.values()], label='HF')\n", "plt.plot(list(mol_configs.keys()), [val[0].fci_energy for val in mol_configs.values()], label='FCI')\n", "plt.plot(list(mol_configs.keys()), [val[2] for val in mol_configs.values()], 'x-', label='VQE')\n", "plt.grid()\n", "plt.xlabel('Bond length of H2 [A]')\n", "plt.ylabel('Ground-state energy for H2 [Ha]')\n", "plt.legend()\n", "plt.title('VQE computed ground-state energy of H2 vs bond-length for STO-3G basis vs HF and FCI')\n", "plt.savefig('H2_PES_byVQE.png')\n", "run_time = time.time()-run_time\n", "print(f\"Total running time: {run_time:.2f} s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Discussion\n", "In the plot of the computed ground-state energy of the $H_2$ molecule against its bond length above, the curve labeled \"HF\" denotes the ground-state energy computed by using the classical Hartree-Fock method, i.e., classical pre-computation alone.\n", "The curve labeled \"FCI\" denotes the ground-state energy of $H_2$ for each bond-length computed using one of the most accurate classical post-Hartree-Fock method, viz., Full Configuration Interaction (FCI).\n", "As can be seen from the curve labeled \"VQE\" in the plot above, the computational results of the quantum computational method VQE agree extremely well with the FCI results over the entire range of bond-lengths simulated for the $H_2$ molecule.\n", "This shows that the quantum computational VQE method can be very competitive with the most accurate classical computational chemistry methods." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# References\n", "[1] Helgaker, T., Jorgensen, P., & Olsen, J. (2014). *Molecular\n", "electronic-structure theory*. John Wiley & Sons. \n", "\n", "[2] McArdle, S., Endo, S., Aspuru-Guzik, A., Benjamin, S., & Yuan, X. (2020). Quantum computational chemistry. *Reviews of Modern Physics, 92(1)*. \n", "\n", "[3] Cao, Y., Romero J., Olson, J. P., Degroote, M., Johnson, P. D., Kieferov, M., Kivlichan, I. D., Menke, T., Peropadre, B., Sawaya, N. P. D., Sim, S., Veis, L., & Aspuru-Guzik A. (2018). Quantum chemistry in the age of quantum computing. *Chemical Reviews, 119(19)*, 10856–10915. \n", "\n", "[4] Hempel, C., Maier, C., Romero, J., McClean, J., Monz, T., Shen, H., Jurcevic, P., Lanyon, B. P., Love, P., Babbush, R., Aspuru-Guzik, A., Blatt, R., & Roos C. F. (2018). Quantum Chemistry Calculations on a Trapped-Ion Quantum Simulator. *Phys.Rev. X, 8(3)*, 031022. \n", "\n", "[5] McClean, J. R., Sung, K. J., Kivlichan, I. D., Bonet-Monroig, X., Cao, Y., Dai, C., Fried, E. S., Gidney, C., Gimby, B., Gokhale, P., Häner, T., Hardikar, T., Havlíček, V., Higgott, O., Huang, C., Izaac, J., Jiang, Z., Kirby, W., Liu, X., ..., Babbush, R. (2017). OpenFermion: The electronic structure package for quantum computers. *arXiv:1710.07629*. \n", "\n", "[6] Sun, Q., Berkelbach, T. C., Blunt, N. S., Booth, G. H., Guo, S., Li, Z., Liu, J., McClain, J., Sayfutyarova, E. R., Sharma, S., Wouters, S., & Chan, G. K.-L. (2017). The Python-based simulations of chemistry framework (PySCF). *WIREs Compututational Molecular Science, 8(e1340).* \n", "\n", "[7] Sun. Q. (2015). Libcint: An efficient general integral library for Gaussian basis functions. *J. Comp. Chem. 36(1664)*.\n" ] } ], "metadata": { "kernelspec": { "display_name": "conda_braket", "language": "python", "name": "conda_braket" }, "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.15" } }, "nbformat": 4, "nbformat_minor": 4 }