{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# **Use case: Quantum chemistry**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## **Practice session: implement VQE for $H_2$**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the lecture we gave an example of the variational quantum eigensolver (VQE) algorithm: computing the ground energy of molecular hydrogen, $H_2$. This is one of the smallest systems one can imagine to start experimenting with VQE. It has also been implemented on an actual quantum computer as reported in [O'Malley etal 2015, *Scalable Quantum Simulation of Molecular Energies*, arXiv:1512.06860v2](https://arxiv.org/abs/1512.06860). The following figure given in the article summarizes the entire setup:\n", "\n", "![Fig. 1. The variational quantum eigensolver circuit and hardware. Source: arXiv:1512.06860v2](img/omalley-etal-2016-vqe.png)\n", "\n", "Fig. 1. The variational quantum eigensolver circuit and hardware. Source: arXiv:1512.06860v2. \n", "\n", "In this session, we are going to implement the VQE algorithm for $H_2$ from scratch using the circuit in Fig. 1 and the data given in the paper." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Synopsis**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we begin, let us summarize the task ahead of us, and what exactly we are going to do. Remember that VQE starts with classical preparations. This involves all the classical computational chemistry computations we summarized in the lecture. As a result of these preparations we have the qubit Hamiltonian and the ansatz circuit that represents the parameterized unitary. We assume that all this work has been completed: the result is the ansatz circuit and the Hamiltonian shown in the figure above.\n", "\n", "Let us summarize the computation shown in the figure.\n", "\n", "The algorithm consists of two parts. Firstly, there is a computation carried out on a quantum computer. This is divided into three parts in the figure, (1) Prepare initial state, (2) Apply Parameterized Ansatz, (3) Measure Expectation values. The entire quantum computation is performed using the circuit shown at the bottom part, 'Software'. Note the circuit contains only two qubits.\n", "\n", "Secondly, there is the classical part, comprising (1) Calculate Energy and (2) Classical Optimizer Suggests New Parameters $\\theta$. The qubit Hamiltonian of $H_2$ is also shown on the right. \n", "\n", "Note that the top part of the figure also shows a rough outline of physical pulses with timings that are needed to run the circuit.\n", "\n", "We will implement both the classical and quantum parts in Python using Qiskit. We will implement the algorithm from scratch *without* using the more advanced Qiskit libraries. We will only rely circuits and measurements. In principle, we will show how to implement the algorithm in *any language that contains circuits and measurements*.\n", "\n", "The rationale of showing this simple, low level implementation is to give a taste of what is going on behind the scenes of the more advanced (Qiskit) libraries." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### **Implementation**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us start by creating a routine for the circuit in the figure." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcUAAACoCAYAAACR3PW0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAd40lEQVR4nO3deXhU5d3G8W8WskEiWcqWIBIIEBICAUUR2YQIUgFRREBtRV7AsLQCtr4VC1UqCGJpLRJ5i4rVoiWghCIUUEgEg60gKImEsIpsQoiBsJpl3j9miASyTiY55yT357q49Dpz5sw9Dw/P7zxnzuJms9lsiIiICO5GBxARETELFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHFUUREREHT6MDiDllZGSUu87ChQuZNGlSmeu0a9fOVZFExEKsOoZopihOe+2114yOICIWZsYxREVRRETEQUVRRETEQUVRnLZixQqjI4iIhZlxDFFRFBERcVBRFKcNGzbM6AgiYmFmHEN0SYaL7N0Euadq/nP9G0Hbu2v+c42ktq45auuao7Y2BxVFF8k9BTlHjU5RN6ita47auuaorc1Bh0/FaRMnTjQ6gohYmBnHEBVFcVp5d6IQESmLGccQHT4Vp/Xs2ZNPP/3U6Bi1UlYu7PsevsuGsxfBZoMGPtA8CFo1gmaBRiesPU6ehf3fw9FsOHcJ3NwgwBfCgiCiMTQKMDph7WXGMURFUZx2+vRpoyPUOge+h43pkHGi5Nf/e9D+3xbBcHd7iGluH8Sl8tKPwSfpcLCcbhzRGPpFQdumNZOrLjHjGKKiWIOmJfRmz7fb8PCoh7u7B00CWzKq73R6dXzI6Gi1jtXa+sd8WP0lbN1XsfW/PQNvbYEOYTD8dvD3qd58ZbFaW1+8Aiu3w47DFVt/3/f2P13D4YFbwadetcYrk9Xa2opUFGvYI/1+zyP9nqOgIJ+k1IXMWTaK1qGxhIa0NjpapbVv397oCGWySltfzoPFm+GQEzvNu4/CiRyY2A8C67s8WoVZpa1zL8Frn9gPmVbWfw/C8R8gvi/U93Z9toqySltXhBnHEJ1oYxAPD0/uvX0sBYX5HDi+y+g4Tlm5cqXRESrEzG1daIO3PnWuIF6VdR4SNsGVPNflcpaZ2zqvAF7f7FxBvOroD/C3ZCgodFksp5m5rSvKjGOI5Yvi2bNnefLJJ2nUqBF+fn50796drVu3Gh2rXHn5P7ImNQGAsJA2BqdxzowZM4yOUCFmbuste2HvybLX+fMj9j9lOXUO/rXTdbmcZea2/vfXcOyHstepSFsfzrL/Fmk0M7d1RZlxDLH04VObzcaQIUPYs2cP8+fPp1mzZvz1r38lLi6O1NRUYmNjjY54g2WfvEhiynwuXcnFw6MeUx9aQnizGACOZe1nzrJR/HniZ3h61GPZJ7MBGNX3WSMjlyoxMZEXXnjB6BilMntbX7gCa3a5bntb98Edre1nTda0stp63X/f4OMd7xSteyL7IB1a9uB3o/5RY/lOn4NN37hue+vT4LZwYw5Zm71fV4YZxxBLzxTXrFlDSkoKS5cu5Re/+AX9+vUjMTGRsLAwpk+fbnS8Eo3qO51Vs3JY8YcsurYbyK79m4peCw1pzZ1R9/PP5Hkcy9pPavoqhvf+jYFprc3sbf2fA/ZDeq60NdO126uostr63q5jeCU+mVfik5n+yPv4eNXn8QF/rNF8n+0Dmwu3V1AI2/a7cIOVYPZ+bXWmLYqFhYXMnz+fiIgIfHx86NixIykpKbRt25Zx48YBkJSURHBwMAMGDCh6n5eXFyNGjGDjxo1cuHDBqPjl8vcLZOpDS/hvxlpS05KKlg/v/Ru2pSfx0rJHmDD4L3h6GHiqWy1h1rb+4pDrt/nlYWN/7yqtrcH+b3rOe4/wxIDZNA1qWWOZbDbYXg1tXR1/f5Vh1n5tdaYtik888QSzZs1i/PjxrFu3juHDhzNy5EgOHjxIly5dAEhLSyMqKgq36y7Uio6OJj8/n4yMDCOiV1iAXxAP9pjKm/9+lsJC+0jm6VGP6JY9uHgll8gWdxicsGwpKSlGR6gws7X1lfyqnfBRmh8Lqme7lVFSWwO8s/F5WjbpwF0dhtZonuwLcP6K67f7wwXIvez67VaG2fp1ZZlxDDFlUVy2bBlvv/02q1ev5umnn6ZPnz5Mnz6dbt26kZ+fX1QUs7OzCQy88dYeQUFBRa8DfP/999xzzz34+fnRsWNHdu40wRkJDkN7/JrscyfYuOPvABw+mc43h1Pp2Ko36794y+B0ZUtPN8HZBpVgprY+mWOfwVSHEznVs93KuL6tv9z3CTsyNzD25/NqPMvxck6uqQoztrXGkKpxs9mq65+m8zp06EBYWBjr1q0rtvyZZ55hwYIF5Obm4u3tTUREBNHR0Xz44YfF1vv444+Ji4tjw4YNxMXF8cADDxAWFsbcuXN55513mDt3LpmZmXh4eJSb5fpZaGnmP7mZjq16V/g7lsRmszEtoRfxg/9Ms5DWTHntLuaN/4SGDX5W6nu+OpDM06/3qdLnlmTKlCnlrrNgwYJy11uwYIGrIhWpDW0d1r4PDz67qdiy8s56LM1T152vsumteHZ/8rqTyYpzRVtnnzvJ04v7MHvMOpoE3VKh97iyrdveOYoBE4o3kqva+l8L7ufgjqSSV66k2tCvr2WmMaQyZc50M8WjR4+SlpbGQw/deIeGI0eOEBUVhbe3/crZ4ODgotngta4uCwoKIjc3l48++ogZM2bg6+vLuHHjKCgo4PPPP6/eL+KE1amLaNP8NiLCOlPfJ4BH42bw+r+mGh2rVjK6rQvzf6y2bRfkVcOxwip49+NZXLh8lpf/+TjTEnozLaE3f14xvsY+v6AOtbXR/bo2MN1M8fPPP6dbt2589NFHDBw4sGj5pUuXaNWqFQMHDmTJkiUAjBkzhtWrV3Pq1KliM7oZM2YwZ84ccnJy2Lt3L/fddx/Hjx8ven3QoEEMHjyYsWPHuiz39veNeRZawzC4dYTrt1uR32MjIyPZs2dPmeu0a9fOVZGK1Ia2zr0Mv6/gdctXZzXXz1JK86s4CG/kXK7r1Ya2Pv4DzFtbsXUr29bPDYYQf+dyXa82tPW1zDyGlMV0M8WQkBAAMjOLn1s+b948Tpw4QefOnYuWDRkyhKysLNavX1+0LC8vj/fff59+/fpRv359Lly4QEBA8dvcBwQEcP78+Wr8FnXD888/b3QEy/L3gYZ+rt+uGxCqJ2gU0/gm8Cz/l5JK8/WC4Aau325dYsYxxHQX74eHhxMTE8Ps2bMJCgoiNDSUFStWsHatfVfv6kk2YJ/x9ejRg9GjRzNv3jyaNm3KwoULOXLkCO+99x4A9evXJzc3t9hnnDt3jgYN1Juravjw4UZHsLSON0OKi0+QbtcMvHUGfjEe7vYbp+/81rXb1RNKqs6MY4jpZoru7u4kJiYSFRVFfHw8o0ePJiQkhIkTJ+Lp6UlMTEzRum5ubqxevZrBgwczdepUBg0axKlTp9iwYUNR8YyIiCArK4szZ84UvS8tLc2UN6K1msjISKMjWFr3CGtssza4S21tSmYcQ0w3UwRo06YNmzdvLrbsscceIzIyEl9f32LLGzZsyOLFi1m8eHGJ2/L39+fnP/85s2bN4qWXXuLdd9/Fzc2NO+4w9/U7Uvs1CoBurV13Z5TWjaF9qGu25QpvrP0d6Yc/I+qW7owZOMfQLOGNICrU/gxFV4htATcHu2ZbYi6mmymWZvv27cUOnVZGQkIC6enpBAYG8uqrr7Jy5coKXY4htVvW2WM8u+Refr3wTqYs6kHC6vJPIXe1IZ1dc/9Mb08YeQe4m+Rw3qGTaVy4fI4/TfiUcxfPcPiksdejubnB8K7g51X1bfn7wIO3Vn07rpR19hiLkp5iz5H/GNqfawNTzhSvd/78eTIzM5kwYYJT72/cuDEbN250caqq++pAMn94eyjhTTtyMvsQrUNjef7xVUbHqrDevXsbHaFKdmRupG/nR+nR4UG86vkwZ9kjHDqxm5ZNO9RYBp96ML4PLNxY+l1XyjsT0tMDxvQy10kfaYe2cmubewDoHBHH7kNbuKVJlKGZbvKDsb3h9U32OwqVpLy29q0H4/pAAwMf6lySHZkb6dImjsYNW/Dy+E2G9efKMuMYYomi2KBBAwoKXHznZBPo0LInbZt35aWx65n73i8MP8RUWQkJCUZHqJDSdj6+OpDM5KEL8apnH+E83D1xd6/5IwhNboJf3QNLt1b+7isN/eAX3V13CYar5F7MZs2211m5ZQHnL+XQq6M5Tqho+TP7A5n/vtX+HMrKaBQAj98FzQw+u7ek/lzfpyGThy7E1/unPSOj+nNlmHEMsczh09roRPZBmgaFA3D67HeE3GSiH4QqID4+3ugIFXJ15+OV+GRiwnsxeehr2Gw2ruRdLBpEDh7/mrMXsmjR2JgTsBoFwNT+cG+MfTZSHk8PuKsN/O995iuIAP5+Qfyy/wu8Ep/M6AF/xN/PgOdZleLmYPjNz6FPJHhVYFrg7QlxUfCbgcYXRCi5P1/bl8H4/lxRZhxDLDFTrK2+PZlOiyZRFBQW4OZmvf2T5ORkoyMUk33uJC/+o/hVyEH+TXh8wB9v2PnYf2wX4c06AnDuYjYLV03iuUeX13jma3l6QP8O9sF657eQeRK+y7Y/QBggqD40D7IXwVtbQn1vQ+OWKbrlXWz4Yil3Rg3mqwPJ9L9ttNGRivH2tP+ee0807DgM+7+Ho9k/zR5/5m9/LmVEY+h8i/0wt1lcvzOdc/50UV8G8/TnijDbGAIqioY6/H067Vt0Iy//CjnnT3Hm3AmCA5oaHcuyggKa8Ep88g3LU9OSbtj52JG5gS4RcRQU5PPSe48y9r6XCQpoUsOJS+blCbe3sv+Bn37nmnG/YZEqrWWTaDw96jEtoTftW3Qz/PfE0vh62WfcdzkeXH+1racPNi5Tea7fmb7alwFT9merUVE00LVPw/7btN0GJqndStr52HdsBw/1eprkr/5J5ndfsOSjZwAYc+8c2t/SzeDEtYPVfiO3iuv78879n/BQr6cBSPk6Uf25ilQUxWnl3bPQLEra+ejR4UHc3d25O3Ykd8eONCqaSKVd359TvlqOu7v9CIjV+rMZxxAVRRfxN+hkB6M+F2D58uWG3KbJFd95SFjlcxvZ1kapi/3aKM5+Z2f6sis+1xWMGkPKoqLoIm3vNjpBzZs5c6YhHboutrVR1NY1py62tVFjSFmsd8qjiIhINVFRFBERcVBRFKctWrTI6AgiYmFmHENUFMVpUVHmvPZMRKzBjGOIiqI4rVevXkZHEBELM+MYoqIoIiLioKIoTrvtttuMjiAiFmbGMURFUZz2xRdfGB1BRCzMjGOIiqKIiIiDiqKIiIiDiqI4bcWKFUZHEBELM+MYoqIoIiLioKIoThs2bJjREUTEwsw4hugpGWI5ezdB7qma/1z/RnXzSQZSM9SvzUFFUSwn9xTkHDU6hYhrqV+bgw6fitMmTpxodAQRsTAzjiEqiuK0SZMmGR1BRCzMjGOIiqI4rWfPnkZHEBELM+MYoqIoTjt9+rTREUTEwsw4huhEG6mVpiX0Zs+32/DwqIe7uwdNAlsyqu90enV8yOhoIk5Tv65+KoritPbt2xsdoUyP9Ps9j/R7joKCfJJSFzJn2Shah8YSGtLa6GgiTqtN/dqMY4gOn4rTVq5caXSECvHw8OTe28dSUJjPgeO7jI4j4hK1oV+bcQxRURSnzZgxw+gIFZKX/yNrUhMACAtpY3AaEdeoDf3ajGOI5Yvi2bNnefLJJ2nUqBF+fn50796drVu3Gh2rTkhMTDQ6QpmWffIi9/++Ifc968tb659j6kNLCG8WA8CL/xjJZ2mritadufR+du7fZFDS0tlsP/1/oa309aTqrm1fm4nbujb066vMOIZYuijabDaGDBnChx9+yPz581m9ejUhISHExcWxc+dOo+OJwUb1nc6qWTms+EMWXdsNZNc1g8OEIX/h7xtmcvFyLlt2f4C/XxCxrc1zr6tCG2zbD3M/+mnZC6vg43T4Md+wWLVSfgEk74EXV/+0bPa/YMteKCg0LldprNyvrcDSJ9qsWbOGlJQU1q5dy7333gvYr3uJiopi+vTprF271uCEYgb+foFMfWgJv3ypFalpSdwZPYTABo14oMcUXkv6FQdPfMXccR8bHbNIoQ3+kQo7DoPbNctzLsKaXZB+DOLvBi9L/+s1h/wC+Fsy7D1ZvK2zcmHldvvy0T3Aw4TTB6v1a6sw4V+1XWFhIfPnzyciIgIfHx86duxISkoKbdu2Zdy4cQAkJSURHBzMgAEDit7n5eXFiBEj2LhxIxcuXDAqfp2QkpJidIQKC/AL4sEeU3nz389SWGjf/e9/2+Mcy9rH/d1/RYBfkMEJf/KfA/aCCFDSUbxDp2Hd1zWZqPb6ON1e+KB4W1/9/7Sj8Onemk5VcVbq1yUx4xhi2qL4xBNPMGvWLMaPH8+6desYPnw4I0eO5ODBg3Tp0gWAtLQ0oqKicHNzK/be6Oho8vPzycjIMCJ6nZGenm50hEoZ2uPXZJ87wcYdfy9a1iy4talOZbfZICWj+KylJNv26zBqVRUUwtZ95a/3aYa5f8+1Qr8ujRnHEFMegFm2bBlvv/02ycnJ9OrVC4A+ffrw5Zdf8sEHHxQVxezs7BKvcwkKCip6HWDmzJkkJiaSkZHB8uXLK/UMr+sLbl0xZcqUctdZsGBBuestWLDAVZGKzH9yMx1b9S5znVfik29YVt8ngA9eyHb6c1NSkrltZB+n318R3vUDeXJx+Rkv50GLyG6c3P95teapzYJC2/PY3PIH5R8uQsNGt5Cb9W215qlt/dpMY4itEmdOmXKmOGfOHAYMGFBUEK9q3bo19erVo0OHDoD9i5ZUtK5fFhERwV/+8he6du1afaFFXMDd3aNa1pUbqa2lJG62ypTQGnD06FGaN2/OG2+8wRNPPFHstZEjR5KRkVF0Zukdd9yBt7f3Dcelly9fzsMPP8z27duLZpUAvXv3ZtKkSaZ82rPZVOTQc2RkJHv27ClznXbt2rkqUpHt7xvz3LmGYXDriOr9jEKb/SzTnItlr+fhDi88APW9qzdPbfZjPvx+JVwp5zC0n5e9rT2ruS7Wtn5t5jGkLKabKR49au8VTZo0Kbb80qVLpKSkFCtyUVFRfPPNNzdMjdPS0vD09Kzxxqxrnn/+eaMj1DrubnBXOddhuwFdblFBrCovT7i9Vfnr3RlR/QWxrjLjGGK6ohgSEgJAZmZmseXz5s3jxIkTdO7cuWjZkCFDyMrKYv369UXL8vLyeP/99+nXrx/169evmdB11PDhw42OUCv1bAstf1bya25AYH0Y1KkmE9VeAzpA44DSXw8NhH5RNZenrjHjGGK6E23Cw8OJiYlh9uzZBAUFERoayooVK4quObx2pjho0CB69OjB6NGjmTdvHk2bNmXhwoUcOXKE9957z6ivUGdU5NCHVJ6Xp/06xHVfQ+p+uJJnX+7hbp8hDuoE/r5GJqw9/Lzh1/fAv3bBF4fs1y2C/e+gazjc1wl86hmZsHYz4xhiuqLo7u5OYmIi48ePJz4+nuDgYH75y18yceJEpk+fTkxMTNG6bm5urF69mmeeeYapU6dy/vx5YmNj2bBhQ7HiKWI1Xp4wpDPcGwPHfrBfqtH4Jh0yrQ5+3vDw7TA4Fk7kgJsbNG2oYlhXma4oArRp04bNmzcXW/bYY48RGRmJr2/xXeSGDRuyePFiFi9eXOr28vLyKCgooLCwkLy8PC5fvoy3t3edvdyiLss6e4zlyS/TJ3Ykr6+egru7B23CbiV+sOsvHXEFL8/SD6WKa/l6QXgjo1NUXdbZY/wp8X+4cPms6fu3GZnuN8XSXH8maWWMHTsWX19ftmzZwqhRo/D19eXbb6v3mqO6oHfv3kZHqLQdmRvp0iaOxg1b8PL4TSyYsIWc86c4dGK30dFEXGJH5kb6dn7UEv3bjGOIKWeK1zt//jyZmZlMmDDBqfcvXbqUpUuXujaUkJCQYHSEMn11IJk/vD2U8KYdOZl9iNahsdT3acjkoQvx9W5QtJ6Hu6euQxPLKal/P//4Kr46kMzkoQvxqucDmLt/m3EMscRMsUGDBhQUFDB58mSjo8g14uPjjY5Qpg4te9K2eVdeiU8mJrwXk4e+xpW8i8UK4sHjX3P2QhYtGpvvCeAiZSmpf9tstmJ93Oz924xjiCVmimJOycnJRkco04nsgzQNCgfg9NnvyDl/mvBmHYteP3cxm4WrJvHco8uNiihSruxzJ3nxH8Wvrg/yb8LjA/5YrH+H3BTK/mO7ivq4Ffq3GccQFUWptb49mU6LJlEUFBbg5ubOjswNdImIA6CgIJ+X3nuUsfe9TFBAk3K2JGKcoIAmJd7zNDUtqVj/Bor6uPq38yxx+FTEGYe/T+eWxlHk5V8h5/wpdu7/hDZhtwKQ8nUimd99wZKPnmFaQm++ObzN4LQilXN9/z5z7gT7ju2gTdit6t9VYLp7n4o5uOqxW2a692nKV8vp1dH5O2jUxL1Ppe5yxb1PnenjRt77tCLq/L1PxTqWLzfvbxUlqUpBFLECq/VxM44h+k1RnDZz5kxD7l3ob9AF1kZ9rtQNdbFfGzWGlEVFUSyn7d1GJxBxPfVrc9DhUxEREQcVRXHaokWLjI4gIhZmxjFERVGcFhWlB82JiPPMOIaoKIrTevXqZXQEEbEwM44hKooiIiIOKooiIiIOuiRDSlSRu0jMnDmzxu82ISLWYNUxRLd5ExERcdDhUxEREQcVRREREQcVRREREQcVRREREQcVRREREQcVRREREQcVRREREQcVRSd999139O3bl8jISKKjo/nd735ndCQREakiFUUneXp6MnfuXPbs2cOXX35JamoqSUlJRscSEZEq0G3enNS0aVOaNm0KgJeXFzExMRw5csTgVCIiUhWaKbrAmTNnWLVqFXFxcUZHERGRKlBRrKIrV64wbNgwnnrqKdPd2FZERCpHNwSvgoKCAh5++GFuvvlm/vSnPxkdR0REqkhFsQrGjBlDYWEhb775Jm5ubkbHERGRKqr1h0+PHz/OyJEjCQwMpEGDBvTv35/09PQqb/ezzz7jzTffZPv27cTGxtKpUydeffXVote1ryEiYj21eqZ46dIlYmNjKSwsZPbs2fj5+TF79mwyMjLYtWsXYWFh1fK5NpuNN5evpV2rFnS/NbpaPkNERFyvVl+SsWTJEjIzM9m9ezdRUVEAdOvWjZYtWzJ79mwWLVpULZ+beego+w4fo0Pb8GrZvoiIVA/LHz7dvXs3Dz74ICEhIfj4+BAREcH06dMBSEpKIjY2tqggAgQGBjJo0CA+/PDDasljs9n4eOsOGgY0oHOHNtXyGSIiUj0sPVPcsWMHPXv2pEWLFsyfP5+bb76ZQ4cOkZqaCkBaWhr9+/e/4X3R0dG8++67nDlzhuDg4DI/43/n/p/T+Z6b/4bT7xUREdd46ZlxFV7X0kVx2rRp+Pv78/nnnxMQEFC0fMyYMQBkZ2cTGBh4w/uCgoKKXi+vKIqISN1h2aJ48eJFtmzZwuTJk4sVxOuVdKlEZS6fqMwext6D3/FW4joe6N+Drp0iK/w+ERExB8sWxR9++IHCwkJCQ0NLXScoKIjs7Owbll9ddnXGWBZnDp9+sH4LH6zfUun3iYiI61VmcmPZE20CAwNxd3fn2LFjpa4TFRVV4jWJaWlpNGnSRIdORUSkGEtfp9inTx/27NlDZmZmiYdQX331VZ566inS09OJjLQfzszJyaFly5aMGDGChIQEl+Sw2WwseieJ3AsXeXrcw3h6eLhkuyIiUrMsXRSvPfv0t7/9LS1atODIkSNs2bKFJUuWcPHiRTp16oSbm1uxi/e/+eYbdu3aRfPmzV2SQ78liojUDpY9fArQpUsXtm3bRkREBFOmTGHgwIG8+OKLNGvWDAA/Pz82b95Mp06dGDNmDMOGDcPX15fk5GSXFUSAs+fO0zgkUNcliohYnKVnimZSWFiIu7ul9zFEROo8FUUREREHTW1EREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQcVBRFREQc/h/sXv1hVsN+6AAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "from qiskit import Aer, ClassicalRegister, QuantumCircuit, QuantumRegister, assemble\n", "from qiskit.circuit import Parameter\n", "\n", "\n", "##########################################################\n", "#\n", "# Routines\n", "#\n", "##########################################################\n", "def create_ansatz_circuit(qreg, creg, parameter):\n", " \"\"\"\n", " Create H2 ansatz circuit given in O'Malley etal 2015.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " parameter - parameter for parameterized circuit\n", "\n", " returns\n", " ansatz_circuit - a parameterized circuit\n", " \"\"\"\n", " ansatz_circuit = QuantumCircuit(qreg, creg)\n", "\n", " # H2 circuit\n", " # 1 - create HF reference state\n", " ansatz_circuit.rx(np.pi, qreg[0])\n", " ansatz_circuit.barrier()\n", " # 2 - ansatz: U(theta)\n", " ansatz_circuit.rx(-np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(np.pi / 2, qreg[1])\n", " # 3\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 4 - parameterized Z-rotation\n", " ansatz_circuit.rz(parameter, qreg[0])\n", " # 5\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 6\n", " ansatz_circuit.rx(np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(-np.pi / 2, qreg[1])\n", " ansatz_circuit.barrier()\n", "\n", " return ansatz_circuit\n", "\n", "\n", "##########################################################\n", "#\n", "# Main program\n", "#\n", "##########################################################\n", "\n", "# variables\n", "theta = Parameter(\"θ\")\n", "q = QuantumRegister(2)\n", "c = ClassicalRegister(2)\n", "\n", "ansatz_circuit = create_ansatz_circuit(q, c, theta)\n", "\n", "ansatz_circuit.draw(\"mpl\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have added barriers for better visualization. \n", "\n", "The new feature here is a *parameter* that is inserted in the `rz()` rotation in L40. It means every time when we update the parameter, we do not need to create a new circuit from scratch. This saves a lot of time in general. VQE circuits can be quite large. Imagine a circuit containing $10$k gates and $300$ parameters. Creating fixed circuits each time one parameter needs to be updated will take much longer than updating a single value in the circuit. For this reason, when implementing VQE, it's always best to use *parameterized circuits* if they are available in the programming language.\n", "\n", "At this point we have accomplished the first two quantum components, (1) Prepare initial state, (2) Apply Parameterized Ansatz." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### **Measure Expectation values**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let us focus on (3) Measure Expectation values. The qubit Hamiltonian we need to measure is as follows\n", "\n", "$$ \n", "H = g_0 \\mathbb{1} + g_1 Z_0 + g_2 Z_1 + g_3 Z_0 Z_1 + g_4 Y_0 Y_1 + g_5 X_0 X_1\n", "$$\n", "\n", "where $g_i$ are (real valued) coefficients that have been computed classically. They are given in the Table 1 in the appendix of O'Malley etal 2015. The values of $g_i$ are functions of hydrogen-hydrogen bond length $R$. We will take the coefficients at bond length $R = 0.75$ (in units $10^{-10}$ m) where the energy is lowest (the actual bond length where the energy is lowest is $0.74$, which is close).\n", "\n", "There is a problem though. We can perform measurements in the $Z$ basis using the Qiskit `measure()` routine. However, we cannot directly measure the Pauli $X$ or $Y$ operators of the terms $X_0 X_1$ and $Y_0 Y_1$. In order to measure the latter, we have to use the following trick: \n", "\n", "- to measure $X$, perform basis transformation from the $X$ and $Z$\n", "- then measure the qubits as usual (in the $Z$ basis) with `measure()`.\n", "\n", "The same applies to $Y$ measurements, in this case we need to perform basis transformation from $Y$ to $Z$.\n", "\n", "In order to compute the average value $\\langle H \\rangle$ of the (total) Hamiltonian we will use the fact that the total average is equal to the sum of averages of its terms,\n", "\n", "$$\n", "\\langle H\\rangle = g_0 \\langle \\mathbb{1}\\rangle + g_1 \\langle Z_0\\rangle + g_2 \\langle Z_1\\rangle + g_3 \\langle Z_0 Z_1\\rangle + g_4 \\langle Y_0 Y_1\\rangle + g_5 \\langle X_0 X_1\\rangle.\n", "$$\n", "\n", "So we need to compute the expectation value of each term, then multiply each with the respective coefficient $g_i$, and finally add averages of all terms. Note that $\\langle\\mathbb{1}\\rangle = 1$, so we can just add the coefficient $g_0$.\n", "\n", "Pseudo code wise, it means: \n", "\n", "- Compute avg of H:\n", " - Loop over terms, compute avg of each term\n", "\n", "Now computing the average of a term like $X_0 X_1$ means that we need to add the respective measurement mini-circuit at the end of the ansatz. This is how the part '(3) Measure Expectation values' is realized in the circuit above in the figure, see the final light grey box.\n", "\n", "In terms of implementation, it means that *instead of looping over terms we can loop over measurement mini-circuits*. And in order to be more efficient, let us create these mini-circuits prior to running the algorithm.\n", "\n", "So let us add a routine that takes in a term (e.g. $X_0 X_1$) and creates the corresponding mini-circuit, which we can add at the end of the main ansatz circuit." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " ┌─┐ \n", "q1_0: ┤M├───\n", " └╥┘┌─┐\n", "q1_1: ─╫─┤M├\n", " ║ └╥┘\n", "c1: 2/═╩══╩═\n", " 0 1 \n", " ┌─┐ \n", "q1_0: ┤M├───\n", " └╥┘┌─┐\n", "q1_1: ─╫─┤M├\n", " ║ └╥┘\n", "c1: 2/═╩══╩═\n", " 0 1 \n", " ┌─┐ \n", "q1_0: ┤M├───\n", " └╥┘┌─┐\n", "q1_1: ─╫─┤M├\n", " ║ └╥┘\n", "c1: 2/═╩══╩═\n", " 0 1 \n", " ┌───┐┌─┐ \n", "q1_0: ┤ H ├┤M├───\n", " ├───┤└╥┘┌─┐\n", "q1_1: ┤ H ├─╫─┤M├\n", " └───┘ ║ └╥┘\n", "c1: 2/══════╩══╩═\n", " 0 1 \n", " ┌─────┐┌───┐┌─┐ \n", "q1_0: ┤ Sdg ├┤ H ├┤M├───\n", " ├─────┤├───┤└╥┘┌─┐\n", "q1_1: ┤ Sdg ├┤ H ├─╫─┤M├\n", " └─────┘└───┘ ║ └╥┘\n", "c1: 2/═════════════╩══╩═\n", " 0 1 \n" ] }, { "data": { "text/plain": [ "[None, None, None, None, None]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "from qiskit import Aer, ClassicalRegister, QuantumCircuit, QuantumRegister, assemble\n", "from qiskit.circuit import Parameter\n", "\n", "\n", "##########################################################\n", "#\n", "# Routines\n", "#\n", "##########################################################\n", "def create_ansatz_circuit(qreg, creg, parameter):\n", " \"\"\"\n", " Create H2 ansatz circuit given in O'Malley etal 2015.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " parameter - parameter for parameterized circuit\n", "\n", " returns\n", " ansatz_circuit - a parameterized circuit\n", " \"\"\"\n", " ansatz_circuit = QuantumCircuit(qreg, creg)\n", "\n", " # H2 circuit\n", " # 1 - create HF reference state\n", " ansatz_circuit.rx(np.pi, qreg[0])\n", " ansatz_circuit.barrier()\n", " # 2 - ansatz: U(theta)\n", " ansatz_circuit.rx(-np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(np.pi / 2, qreg[1])\n", " # 3\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 4 - parameterized Z-rotation\n", " ansatz_circuit.rz(parameter, qreg[0])\n", " # 5\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 6\n", " ansatz_circuit.rx(np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(-np.pi / 2, qreg[1])\n", " ansatz_circuit.barrier()\n", "\n", " return ansatz_circuit\n", "\n", "\n", "def paulis_to_measure_circ(qreg, creg, hamlist):\n", " \"\"\"\n", " Generate measure circuits from hamiltonian pauli list\n", " for measurements in different bases.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " hamlist - total hamiltonian, i.e. pauli strings\n", "\n", " returns\n", " circuits - list of circuits that can be used to average over\n", " \"\"\"\n", " circuits = []\n", "\n", " for elem in hamlist:\n", " minicirc = QuantumCircuit(qreg, creg)\n", "\n", " for qubitno, qubitop in elem:\n", " if qubitop == \"Z\":\n", " pass\n", " elif qubitop == \"X\":\n", " minicirc.h(qreg[qubitno])\n", " elif qubitop == \"Y\":\n", " minicirc.sdg(qreg[qubitno])\n", " minicirc.h(qreg[qubitno])\n", " else:\n", " assert False, \"Error: INVALID qubit operation\"\n", "\n", " minicirc.measure(qreg, creg)\n", " circuits.append(minicirc)\n", "\n", " return circuits\n", "\n", "\n", "##########################################################\n", "#\n", "# Main program\n", "#\n", "##########################################################\n", "\n", "# variables\n", "theta = Parameter(\"θ\")\n", "q = QuantumRegister(2)\n", "c = ClassicalRegister(2)\n", "hamiltonian = (\n", " ((0, \"Z\"),),\n", " ((1, \"Z\"),),\n", " ((0, \"Z\"), (1, \"Z\")),\n", " ((0, \"X\"), (1, \"X\")),\n", " ((0, \"Y\"), (1, \"Y\")),\n", ")\n", "\n", "ansatz_circuit = create_ansatz_circuit(q, c, theta)\n", "ham_mini_circuits = paulis_to_measure_circ(q, c, hamiltonian)\n", "\n", "# print result\n", "[print(circ) for circ in ham_mini_circuits]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The new routine `paulis_to_measure_circ()` takes in the total qubit Hamiltonian and outputs circuits which correspond to the terms of the Hamiltonian. We also needed to represent the Hamiltonian, to this end we use a tuple of tuples in L96. The smallest element has the form `(qubit_number, Pauli_matrix)`, where `qubit_number` is the index of the qubit for the Pauli observable.\n", "\n", "Let's summarize the computation of the average of the Hamiltonian. We need to compute the averages of terms. We do this by looping over terms, which are represented as mini-circuits. For each mini-circuit, we will take the ansatz, add to it the mini-circuit, and run the resulting circuit, collect measurement statistics and compute the average.\n", "\n", "### **Scan the parameter interval**\n", "\n", "We will not bother with classical optimization in this small implementation. The circuit has only one parameter—the angle $\\theta$ of the $R_z$ rotation—which ranges between $[-\\pi, \\pi]$. Instead of searching for the minimum value of $\\langle H\\rangle$, we will scan this interval and record the minimum energy.\n", "\n", "Let us add a loop to do this." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "H_min = 0.2252\n", "\n" ] } ], "source": [ "import numpy as np\n", "from qiskit import Aer, ClassicalRegister, QuantumCircuit, QuantumRegister, assemble\n", "from qiskit.circuit import Parameter\n", "\n", "\n", "##########################################################\n", "#\n", "# Routines\n", "#\n", "##########################################################\n", "def create_ansatz_circuit(qreg, creg, parameter):\n", " \"\"\"\n", " Create H2 ansatz circuit given in O'Malley etal 2015.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " parameter - parameter for parameterized circuit\n", "\n", " returns\n", " ansatz_circuit - a parameterized circuit\n", " \"\"\"\n", " ansatz_circuit = QuantumCircuit(qreg, creg)\n", "\n", " # H2 circuit\n", " # 1 - create HF reference state\n", " ansatz_circuit.rx(np.pi, qreg[0])\n", " ansatz_circuit.barrier()\n", " # 2 - ansatz: U(theta)\n", " ansatz_circuit.rx(-np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(np.pi / 2, qreg[1])\n", " # 3\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 4 - parameterized Z-rotation\n", " ansatz_circuit.rz(parameter, qreg[0])\n", " # 5\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 6\n", " ansatz_circuit.rx(np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(-np.pi / 2, qreg[1])\n", " ansatz_circuit.barrier()\n", "\n", " return ansatz_circuit\n", "\n", "\n", "def paulis_to_measure_circ(qreg, creg, hamlist):\n", " \"\"\"\n", " Generate measure circuits from hamiltonian pauli list\n", " for measurements in different bases.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " hamlist - total hamiltonian, i.e. pauli strings\n", "\n", " returns\n", " circuits - list of circuits that can be used to average over\n", " \"\"\"\n", " circuits = []\n", "\n", " for elem in hamlist:\n", " minicirc = QuantumCircuit(qreg, creg)\n", "\n", " for qubitno, qubitop in elem:\n", " if qubitop == \"Z\":\n", " pass\n", " elif qubitop == \"X\":\n", " minicirc.h(qreg[qubitno])\n", " elif qubitop == \"Y\":\n", " minicirc.sdg(qreg[qubitno])\n", " minicirc.h(qreg[qubitno])\n", " else:\n", " assert False, \"Error: INVALID qubit operation\"\n", "\n", " minicirc.measure(qreg, creg)\n", " circuits.append(minicirc)\n", "\n", " return circuits\n", "\n", "\n", "def get_ham_avg(params, ansatz, coeffs, ham_mini_circuits, obs, shots):\n", " \"\"\"\n", " Compute average of the total qubit Hamiltonian.\n", "\n", " args\n", " params - parameter value for the ansatz circuit\n", " ansatz - ansatz circuit\n", " ham_mini_ciruits - pre-built circuits to measure Hamiltonian terms\n", " obs - data for observables\n", " shots - number of shots\n", "\n", " returns\n", " avg - average of Hamiltonian\n", " \"\"\"\n", " avg: float = 0.0\n", "\n", " # TODO\n", "\n", " return avg\n", "\n", "\n", "##########################################################\n", "#\n", "# Main program\n", "#\n", "##########################################################\n", "\n", "# variables\n", "theta = Parameter(\"θ\")\n", "q = QuantumRegister(2)\n", "c = ClassicalRegister(2)\n", "hamiltonian = (\n", " ((0, \"Z\"),),\n", " ((1, \"Z\"),),\n", " ((0, \"Z\"), (1, \"Z\")),\n", " ((0, \"X\"), (1, \"X\")),\n", " ((0, \"Y\"), (1, \"Y\")),\n", ")\n", "shots: int = 1000\n", "no_of_points: int = 40\n", "param_range = np.arange(-np.pi, np.pi, np.pi / no_of_points)\n", "# g_i coeffs for R = 0.75\n", "g_coeff = [0.3435, -0.4347, 0.5716, 0.091, 0.091]\n", "g_coeff_id = 0.2252\n", "# graph\n", "graph = []\n", "#\n", "obs = []\n", "\n", "# start\n", "ansatz_circuit = create_ansatz_circuit(q, c, theta)\n", "ham_mini_circuits = paulis_to_measure_circ(q, c, hamiltonian)\n", "\n", "# record min val of hamiltonian\n", "h_min: float = 100.0\n", "\n", "# scan theta domain [-pi, pi]\n", "for x in param_range:\n", " h_avg = get_ham_avg(\n", " {theta: x}, ansatz_circuit, g_coeff, ham_mini_circuits, obs, shots\n", " )\n", " h_avg += g_coeff_id # add coeff of id term\n", " graph.append(h_avg) # collect data for plotting\n", " if h_avg < h_min: # record min value\n", " h_min = h_avg\n", "\n", "# result\n", "print(\"\\nH_min =\", h_min)\n", "print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The loop in lines L141-147 walks through the interval $[-\\pi, \\pi]$ and finds the minimum energy. Of course, the answer is wrong—the value $0.2252$ is just the value of the constant term—since we have not implemented the actual computation of the average.\n", "\n", "\n", "### **Mapping outcome states to values**\n", "\n", "To actually compute the average of a Hamiltonian term, we need to know how the outcomes of the measurement $00, 01, 10, 11$ that Qiskit uses correspond to the actual real valued outcomes $+1, -1$ for a particular observable like $Z_0$ or $Z_0 Z_1$ we measure. Without discussing this at length here, we are going to introduce dictionaries that record the mapping for each observable. Of course, there are other ways to implement this.\n", "\n", "Let us fill in the generic routine for computing the average." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "H_min = 0.2252\n", "\n" ] } ], "source": [ "import numpy as np\n", "from qiskit import Aer, ClassicalRegister, QuantumCircuit, QuantumRegister, assemble\n", "from qiskit.circuit import Parameter\n", "\n", "\n", "##########################################################\n", "#\n", "# Routines\n", "#\n", "##########################################################\n", "def create_ansatz_circuit(qreg, creg, parameter):\n", " \"\"\"\n", " Create H2 ansatz circuit given in O'Malley etal 2015.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " parameter - parameter for parameterized circuit\n", "\n", " returns\n", " ansatz_circuit - a parameterized circuit\n", " \"\"\"\n", " ansatz_circuit = QuantumCircuit(qreg, creg)\n", "\n", " # H2 circuit\n", " # 1 - create HF reference state\n", " ansatz_circuit.rx(np.pi, qreg[0])\n", " ansatz_circuit.barrier()\n", " # 2 - ansatz: U(theta)\n", " ansatz_circuit.rx(-np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(np.pi / 2, qreg[1])\n", " # 3\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 4 - parameterized Z-rotation\n", " ansatz_circuit.rz(parameter, qreg[0])\n", " # 5\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 6\n", " ansatz_circuit.rx(np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(-np.pi / 2, qreg[1])\n", " ansatz_circuit.barrier()\n", "\n", " return ansatz_circuit\n", "\n", "\n", "def paulis_to_measure_circ(qreg, creg, hamlist):\n", " \"\"\"\n", " Generate measure circuits from hamiltonian pauli list\n", " for measurements in different bases.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " hamlist - total hamiltonian, i.e. pauli strings\n", "\n", " returns\n", " circuits - list of circuits that can be used to average over\n", " \"\"\"\n", " circuits = []\n", "\n", " for elem in hamlist:\n", " minicirc = QuantumCircuit(qreg, creg)\n", "\n", " for qubitno, qubitop in elem:\n", " if qubitop == \"Z\":\n", " pass\n", " elif qubitop == \"X\":\n", " minicirc.h(qreg[qubitno])\n", " elif qubitop == \"Y\":\n", " minicirc.sdg(qreg[qubitno])\n", " minicirc.h(qreg[qubitno])\n", " else:\n", " assert False, \"Error: INVALID qubit operation\"\n", "\n", " minicirc.measure(qreg, creg)\n", " circuits.append(minicirc)\n", "\n", " return circuits\n", "\n", "\n", "def get_paulistr_avg(counts, obs, shots):\n", " \"\"\"\n", " Calculate average of a Pauli string (e.g. XXY).\n", "\n", " args\n", " counts - result of simulation: {'00': 45, '01': 34, . . .}\n", " obs - data for calculating probabilities\n", " shots - number of repetitions\n", "\n", " returns\n", " avg - the average value\n", " \"\"\"\n", " avg: float = 0.0\n", "\n", " # TODO\n", "\n", " return avg\n", "\n", "\n", "def get_ham_avg(params, ansatz, coeffs, ham_mini_circuits, obs, shots):\n", " \"\"\"\n", " Compute average of the total qubit Hamiltonian.\n", "\n", " args\n", " params - parameter value for the ansatz circuit\n", " ansatz - ansatz circuit\n", " ham_mini_ciruits - pre-built circuits to measure Hamiltonian terms\n", " obs - data for observables\n", " shots - number of shots\n", "\n", " returns\n", " avg - average of Hamiltonian\n", " \"\"\"\n", " avg: float = 0.0\n", " sim = Aer.get_backend(\"aer_simulator\")\n", "\n", " for coeff, measure_circuit, obsval in zip(coeffs, ham_mini_circuits, obs):\n", " total_circuit = ansatz.compose(measure_circuit)\n", " total_circuit_bind = total_circuit.bind_parameters(params)\n", " qobj = assemble(total_circuit_bind)\n", " counts = sim.run(qobj, shots=shots).result().get_counts()\n", " paulistring_avg = get_paulistr_avg(counts, obsval, shots)\n", " avg += coeff * paulistring_avg\n", "\n", " return avg\n", "\n", "\n", "##########################################################\n", "#\n", "# Main program\n", "#\n", "##########################################################\n", "\n", "# variables\n", "theta = Parameter(\"θ\")\n", "q = QuantumRegister(2)\n", "c = ClassicalRegister(2)\n", "hamiltonian = (\n", " ((0, \"Z\"),),\n", " ((1, \"Z\"),),\n", " ((0, \"Z\"), (1, \"Z\")),\n", " ((0, \"X\"), (1, \"X\")),\n", " ((0, \"Y\"), (1, \"Y\")),\n", ")\n", "shots: int = 1000\n", "no_of_points: int = 40\n", "param_range = np.arange(-np.pi, np.pi, np.pi / no_of_points)\n", "# g_i coeffs for R = 0.75\n", "g_coeff = [0.3435, -0.4347, 0.5716, 0.091, 0.091]\n", "g_coeff_id = 0.2252\n", "# graph\n", "graph = []\n", "# map outcome states to outcome values\n", "obs = [\n", " {\"00\": 1, \"10\": 1, \"01\": -1, \"11\": -1}, # 1 x Z\n", " {\"00\": 1, \"10\": -1, \"01\": 1, \"11\": -1}, # Z x 1\n", " {\"00\": 1, \"10\": -1, \"01\": -1, \"11\": 1}, # Z x Z\n", " {\"00\": 1, \"10\": -1, \"01\": -1, \"11\": 1}, # Y x Y\n", " {\"00\": 1, \"10\": -1, \"01\": -1, \"11\": 1},\n", "] # X x X\n", "\n", "# start\n", "ansatz_circuit = create_ansatz_circuit(q, c, theta)\n", "ham_mini_circuits = paulis_to_measure_circ(q, c, hamiltonian)\n", "\n", "# record min val of hamiltonian\n", "h_min: float = 100.0\n", "\n", "# scan theta domain [-pi, pi]\n", "for x in param_range:\n", " h_avg = get_ham_avg(\n", " {theta: x}, ansatz_circuit, g_coeff, ham_mini_circuits, obs, shots\n", " )\n", " h_avg += g_coeff_id # add coeff of id term\n", " graph.append(h_avg) # collect data for plotting\n", " if h_avg < h_min: # record min value\n", " h_min = h_avg\n", "\n", "# result\n", "print(\"\\nH_min =\", h_min)\n", "print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The routine `get_ham_avg()` in L105, and more specifically the loop in L122 implements the idea we outlined above. It iterates over all mini-circuits, combining each mini-circuit with the ansatz circuit, running it, collecting the statistics and finally computing the average for the term. \n", "\n", "The answer is still wrong because the very last step, computing the average of a single Pauli string in `get_paulistr_avg()` is still not implemented. Let us fill in that too." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "H_min = -1.1556946\n", "\n" ] } ], "source": [ "import numpy as np\n", "from qiskit import Aer, ClassicalRegister, QuantumCircuit, QuantumRegister, assemble\n", "from qiskit.circuit import Parameter\n", "\n", "\n", "##########################################################\n", "#\n", "# Routines\n", "#\n", "##########################################################\n", "def create_ansatz_circuit(qreg, creg, parameter):\n", " \"\"\"\n", " Create H2 ansatz circuit given in O'Malley etal 2015.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " parameter - parameter for parameterized circuit\n", "\n", " returns\n", " ansatz_circuit - a parameterized circuit\n", " \"\"\"\n", " ansatz_circuit = QuantumCircuit(qreg, creg)\n", "\n", " # H2 circuit\n", " # 1 - create HF reference state\n", " ansatz_circuit.rx(np.pi, qreg[0])\n", " ansatz_circuit.barrier()\n", " # 2 - ansatz: U(theta)\n", " ansatz_circuit.rx(-np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(np.pi / 2, qreg[1])\n", " # 3\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 4 - parameterized Z-rotation\n", " ansatz_circuit.rz(parameter, qreg[0])\n", " # 5\n", " ansatz_circuit.cx(qreg[1], qreg[0])\n", " # 6\n", " ansatz_circuit.rx(np.pi / 2, qreg[0])\n", " ansatz_circuit.ry(-np.pi / 2, qreg[1])\n", " ansatz_circuit.barrier()\n", "\n", " return ansatz_circuit\n", "\n", "\n", "def paulis_to_measure_circ(qreg, creg, hamlist):\n", " \"\"\"\n", " Generate measure circuits from hamiltonian pauli list\n", " for measurements in different bases.\n", "\n", " args\n", " qreg - qubits\n", " creg - classical bits\n", " hamlist - total hamiltonian, i.e. pauli strings\n", "\n", " returns\n", " circuits - list of circuits that can be used to average over\n", " \"\"\"\n", " circuits = []\n", "\n", " for elem in hamlist:\n", " minicirc = QuantumCircuit(qreg, creg)\n", "\n", " for qubitno, qubitop in elem:\n", " if qubitop == \"Z\":\n", " pass\n", " elif qubitop == \"X\":\n", " minicirc.h(qreg[qubitno])\n", " elif qubitop == \"Y\":\n", " minicirc.sdg(qreg[qubitno])\n", " minicirc.h(qreg[qubitno])\n", " else:\n", " assert False, \"Error: INVALID qubit operation\"\n", "\n", " minicirc.measure(qreg, creg)\n", " circuits.append(minicirc)\n", "\n", " return circuits\n", "\n", "\n", "def get_paulistr_avg(counts, obs, shots):\n", " \"\"\"\n", " Calculate average of a Pauli string (e.g. IX, XX or YY).\n", "\n", " args\n", " counts - result of simulation: {'00': 45, '01': 34, . . .}\n", " obs - data for calculating probabilities\n", " shots - number of repetitions\n", "\n", " returns\n", " avg - the average value\n", " \"\"\"\n", " avg: float = 0.0\n", " fshots = float(shots)\n", "\n", " for k, v in counts.items():\n", " avg += obs[k] * (v / fshots)\n", "\n", " return avg\n", "\n", "\n", "def get_ham_avg(params, ansatz, coeffs, ham_mini_circuits, obs, shots):\n", " \"\"\"\n", " Compute average of the total qubit Hamiltonian.\n", "\n", " args\n", " params - parameter value for the ansatz circuit\n", " ansatz - ansatz circuit\n", " ham_mini_ciruits - pre-built circuits to measure Hamiltonian terms\n", " obs - data for observables\n", " shots - number of shots\n", "\n", " returns\n", " avg - average of Hamiltonian\n", " \"\"\"\n", " avg: float = 0.0\n", " sim = Aer.get_backend(\"aer_simulator\")\n", "\n", " for coeff, measure_circuit, obsval in zip(coeffs, ham_mini_circuits, obs):\n", " total_circuit = ansatz.compose(measure_circuit)\n", " total_circuit_bind = total_circuit.bind_parameters(params)\n", " qobj = assemble(total_circuit_bind)\n", " counts = sim.run(qobj, shots=shots).result().get_counts()\n", " paulistring_avg = get_paulistr_avg(counts, obsval, shots)\n", " avg += coeff * paulistring_avg\n", "\n", " return avg\n", "\n", "\n", "##########################################################\n", "#\n", "# Main program\n", "#\n", "##########################################################\n", "\n", "# variables\n", "theta = Parameter(\"θ\")\n", "q = QuantumRegister(2)\n", "c = ClassicalRegister(2)\n", "hamiltonian = (\n", " ((0, \"Z\"),),\n", " ((1, \"Z\"),),\n", " ((0, \"Z\"), (1, \"Z\")),\n", " ((0, \"X\"), (1, \"X\")),\n", " ((0, \"Y\"), (1, \"Y\")),\n", ")\n", "shots: int = 1000\n", "no_of_points: int = 40\n", "param_range = np.arange(-np.pi, np.pi, np.pi / no_of_points)\n", "# g_i coeffs for R = 0.75\n", "g_coeff = [0.3435, -0.4347, 0.5716, 0.091, 0.091]\n", "g_coeff_id = 0.2252\n", "# graph\n", "graph = []\n", "# map outcome states to outcome values\n", "obs = [\n", " {\"00\": 1, \"10\": 1, \"01\": -1, \"11\": -1}, # 1 x Z\n", " {\"00\": 1, \"10\": -1, \"01\": 1, \"11\": -1}, # Z x 1\n", " {\"00\": 1, \"10\": -1, \"01\": -1, \"11\": 1}, # Z x Z\n", " {\"00\": 1, \"10\": -1, \"01\": -1, \"11\": 1}, # Y x Y\n", " {\"00\": 1, \"10\": -1, \"01\": -1, \"11\": 1},\n", "] # X x X\n", "\n", "# start\n", "ansatz_circuit = create_ansatz_circuit(q, c, theta)\n", "ham_mini_circuits = paulis_to_measure_circ(q, c, hamiltonian)\n", "\n", "# record min val of hamiltonian\n", "h_min: float = 100.0\n", "\n", "# scan theta domain [-pi, pi]\n", "for x in param_range:\n", " h_avg = get_ham_avg(\n", " {theta: x}, ansatz_circuit, g_coeff, ham_mini_circuits, obs, shots\n", " )\n", " h_avg += g_coeff_id # add coeff of id term\n", " graph.append(h_avg) # collect data for plotting\n", " if h_avg < h_min: # record min value\n", " h_min = h_avg\n", "\n", "# result\n", "print(\"\\nH_min =\", h_min)\n", "print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This result is now close to the result obtained in `arXiv:1512.06860v2` where the circuit was executed on a real, physical quantum computer—compare it with the result in the paper in Fig. 3 at bond length $R = 0.75$.\n", "\n", "\n", "### **Further problems**\n", "\n", "The article `arXiv:1512.06860v2` contains several figures that visualize the results in more detail. Using the above code, you can reproduce them with a little more work. This will give a better understanding of how VQE works.\n", "\n", "> **Exercise 1**. Record the average of each Hamiltonian term at $R = 0.75$ and plot the averages. The resulting graph should be similar to Fig. 2a in `arXiv:1512.06860v2`.\n", "\n", "> **Exercise 2**. Extend the program to find the minimum energies at all bond lengths given in Table 1 in the Appendix. Plot the resulting graph. It should look like Fig. 3a in `arXiv:1512.06860v2`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## **Summary**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this session, we have given an implementation of the VQE algorithm for the $H_2$ molecule based on the circuit and the data in `arXiv:1512.06860v2`." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.12" } }, "nbformat": 4, "nbformat_minor": 4 }