{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# VQEによる量子化学計算" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "このチュートリアルでは、Amazon Braket で PennyLane を使用して量子化学の重要な問題、すなわち分子の基底状態エネルギーを見つける方法を説明します。この問題は、変分量子固有値ソルバー (Variational Quantum Eigensolver; VQE) アルゴリズムを実装することにより、NISQ ハードウェアを利用して解くことができます。量子化学と VQE の詳細については、[Braket VQE ノートブック](../Hybrid_quantum_algorithms/vqe_Chemistry/vqe_Chemistry_braket.ipynb) や [PennyLane チュートリアル](https://pennylane.ai/qml/demos/tutorial_qubit_rotation.html) を参考にして下さい。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "注: このノートブックの実行には pennylane>=0.18 以上と amazon-braket-pennylane-plugin>=1.5.0 が必要です。\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 量子化学から量子回路へ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "まず最初のステップは、量子化学の問題を量子コンピュータで扱えるよう変換することです。PennyLane では ``qchem`` パッケージを使います。ローカルマシン上で実行している場合、``qchem`` パッケージは [これら](https://pennylane.readthedocs.io/en/stable/introduction/chemistry.html) の指示に従って別途インストールする必要があります。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pennylane as qml\n", "from pennylane import qchem\n", "from pennylane import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "入力化学データは、多くの場合、分子に関する詳細を含むジオメトリファイルの形式で提供されます。ここで、[h2.xyz](./qchem/h2.xyz) ファイルに保存された $\\mathrm {H} _2$ の原子構造を考えます。量子ビットハミルトニアンは ``qchem`` パッケージを使って構成されています。" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " (-0.24274280513140456) [Z2]\n", "+ (-0.24274280513140456) [Z3]\n", "+ (-0.04207897647782281) [I0]\n", "+ (0.1777128746513994) [Z1]\n", "+ (0.17771287465139946) [Z0]\n", "+ (0.12293305056183798) [Z0 Z2]\n", "+ (0.12293305056183798) [Z1 Z3]\n", "+ (0.16768319457718958) [Z0 Z3]\n", "+ (0.16768319457718958) [Z1 Z2]\n", "+ (0.17059738328801055) [Z0 Z1]\n", "+ (0.17627640804319586) [Z2 Z3]\n", "+ (-0.04475014401535161) [Y0 Y1 X2 X3]\n", "+ (-0.04475014401535161) [X0 X1 Y2 Y3]\n", "+ (0.04475014401535161) [Y0 X1 X2 Y3]\n", "+ (0.04475014401535161) [X0 Y1 Y2 X3]\n" ] } ], "source": [ "symbols, coordinates = qchem.read_structure('qchem/h2.xyz')\n", "h, qubits = qchem.molecular_hamiltonian(symbols, coordinates, name=\"h2\")\n", "print(h)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "VQE アルゴリズムでは、変分量子回路上の上記のハミルトニアンの期待値を測定することにより、$\\mathrm {H} _2$ 分子のエネルギーを計算します。我々の目的は、ハミルトニアンの期待値が最小になるように回路のパラメータを訓練し、それによって分子の基底状態エネルギーを見つけることです。\n", "\n", "このチュートリアルでは、トータルスピンも計算します。そのために、``qchem`` パッケージを使ってトータルスピン演算子 $S^2$ を構築します。" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " (0.375) [Z1]\n", "+ (0.375) [Z0]\n", "+ (0.375) [Z2]\n", "+ (0.375) [Z3]\n", "+ (0.75) [I0]\n", "+ (-0.375) [Z0 Z1]\n", "+ (-0.375) [Z2 Z3]\n", "+ (-0.125) [Z0 Z3]\n", "+ (-0.125) [Z1 Z2]\n", "+ (0.125) [Z0 Z2]\n", "+ (0.125) [Z1 Z3]\n", "+ (-0.125) [Y0 X1 X2 Y3]\n", "+ (-0.125) [X0 Y1 Y2 X3]\n", "+ (0.125) [Y0 X1 Y2 X3]\n", "+ (0.125) [Y0 Y1 X2 X3]\n", "+ (0.125) [Y0 Y1 Y2 Y3]\n", "+ (0.125) [X0 X1 X2 X3]\n", "+ (0.125) [X0 X1 Y2 Y3]\n", "+ (0.125) [X0 Y1 X2 Y3]\n" ] } ], "source": [ "electrons = 2 # Molecular hydrogen has two electrons\n", "\n", "S2 = qchem.spin2(electrons, qubits)\n", "print(S2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ansatz 回路の定義" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ここで、ハミルトニアンの基底状態を準備するために訓練される ansatz 回路を設定します。最初のステップは、ローカルの Braket デバイスを読み込むことです。\n", "\n", "このチュートリアルでは、[Delgado et al. (2020)](https://arxiv.org/abs/2106.13840) の [`AllSinglesDoubles`](https://pennylane.readthedocs.io/en/stable/code/api/pennylane.templates.subroutines.UCCSD.html) ansatz という化学インスパイアドな回路を使います。これを使用するには、量子化学の観点から追加の入力項目をいくつか定義する必要があります。" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Hartree-Fock state\n", "hf_state = qchem.hf_state(electrons, qubits)\n", "# generate single- and double-excitations\n", "singles, doubles = qchem.excitations(electrons, qubits)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "注: さまざまな ansatz とテンプレートが利用可能で、違うものを選ぶと、回路の深さや学習可能なパラメータ数が異なります。\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "この ansatz 回路は簡単に定義できます:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def circuit(params, wires):\n", " qml.templates.AllSinglesDoubles(params, wires, hf_state, singles, doubles)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "出力の測定はまだ定義されていないことに注意してください。これは次のステップで行います。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## エネルギーとトータルスピンの測定" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "回路を実行するために、デバイスをインスタンス化します。Braket Local Simulator を使います。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dev = qml.device(\"braket.local.qubit\", wires=qubits)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "先に説明したように、$\\mathrm {H} _2$ のエネルギーに対応する量子ビットハミルトニアンの期待値を最小化したいと考えています。このハミルトニアンとトータルスピン $\\hat {S} ^2$ 演算子の期待値は、以下を使用して定義できます。" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "wires = dev.wires.tolist()\n", "\n", "@qml.qnode(dev)\n", "def energy_expval(params):\n", " circuit(params, wires)\n", " return qml.expval(h)\n", "\n", "@qml.qnode(dev)\n", "def S2_expval(params):\n", " circuit(params, wires)\n", " return qml.expval(S2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ここで、`dev` は `shots` パラメータを持たないことに注意してください。これは、1回の評価でハミルトニアンの厳密な期待値を計算できることを意味します。\n", "\n", "次に、ランダムな値をいくつか初期化し、エネルギーとスピンを評価しましょう。準備された状態のトータルスピン$S$は、$S=-\\frac {1} {2} +\\sqrt {\\frac {1} {4} +\\langle\\hat {S} ^2\\rangle}$ を用いて、期待値 $\\langle \\hat {S}^2 \\rangle$ から得ることができます。$S$ を計算する関数はこのように定義することができます:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def spin(params):\n", " return -0.5 + np.sqrt(1 / 4 + S2_expval(params))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "np.random.seed(1967)\n", "params = np.random.normal(0, np.pi, len(singles) + len(doubles))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "エネルギーとトータルスピンは、" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Energy: -0.27304966436310224\n", "Spin: 0.11000908988780533\n" ] } ], "source": [ "print(\"Energy:\", energy_expval(params))\n", "print(\"Spin: \", spin(params))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ランダムなパラメータを選んだので、測定されたエネルギーは基底状態エネルギーに対応しておらず、準備状態はトータルスピン演算子の固有状態ではありません。ここで、最小エネルギーを見つけるためにパラメータをトレーニングする必要があります。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## エネルギー最小化" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "エネルギーは、オプティマイザーを選択し、標準の最適化ループを実行することで最小化できます。" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "opt = qml.GradientDescentOptimizer(stepsize=0.4)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "iterations = 40" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def run_vqe(energy_expval, spin, opt, initial_params, iterations):\n", " energies = []\n", " spins = []\n", " params = initial_params\n", "\n", " for i in range(iterations):\n", " params = opt.step(energy_expval, params)\n", " \n", " e = energy_expval(params)\n", " s = spin(params)\n", " \n", " energies.append(e)\n", " spins.append(s)\n", " \n", " if (i + 1) % 5 == 0:\n", " print(f\"Completed iteration {i + 1}\")\n", " print(\"Energy:\", e)\n", " print(\"Total spin:\", s)\n", " print(\"----------------\")\n", " \n", " print(f\"Optimized energy: {e} Ha\")\n", " print(f\"Corresponding total spin: {s}\")\n", " return energies, spins" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Completed iteration 5\n", "Energy: -0.6355679018250266\n", "Total spin: 0.11730536055750063\n", "----------------\n", "Completed iteration 10\n", "Energy: -0.9411439477019683\n", "Total spin: 0.05338891122286116\n", "----------------\n", "Completed iteration 15\n", "Energy: -1.0953476071202999\n", "Total spin: 0.012154011361553807\n", "----------------\n", "Completed iteration 20\n", "Energy: -1.1289360208519634\n", "Total spin: 0.0022417391335066705\n", "----------------\n", "Completed iteration 25\n", "Energy: -1.1349263900521251\n", "Total spin: 0.00039974850346669033\n", "----------------\n", "Completed iteration 30\n", "Energy: -1.1359701193564353\n", "Total spin: 7.088996239612566e-05\n", "----------------\n", "Completed iteration 35\n", "Energy: -1.1361513824252214\n", "Total spin: 1.2559325300198765e-05\n", "----------------\n", "Completed iteration 40\n", "Energy: -1.136182845906552\n", "Total spin: 2.224718140930726e-06\n", "----------------\n", "Optimized energy: -1.136182845906552 Ha\n", "Corresponding total spin: 2.224718140930726e-06\n" ] } ], "source": [ "energies, spins = run_vqe(energy_expval, spin, opt, params, iterations)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "水素分子の基底状態エネルギーの正確な値は ``-1.136189454088`` Hartrees (Ha) として理論的に計算されています。最適化されたエネルギーの誤差は、Hartree の $10^ {-5} $ 未満であることに注意してください。さらに、最適化された状態は、$\\mathrm {H} _2$分子の基底状態に予想される固有値$S=0$を持つトータルスピン演算子の固有状態です。したがって、上記の結果は非常に有望に見えます！反復回数を増やすと、理論値にさらに近づきます。\n", "\n", "最適化中に 2 つの量がどのように変化したかを可視化してみましょう。" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/BklEQVR4nO3deVxVdf7H8dfn3guIuygqgoqImuKCgppt4taiqVOu1ZROi5Y1jU1NNUs1NVb+amq0mTarSasZNa1J0zSXNCstw8IySnHBEDfAfRf4/P64V0NZRATOBT7Px+N0z7nn3Hvf94R8+J7l+xVVxRhjjCmMy+kAxhhj/JsVCmOMMUWyQmGMMaZIViiMMcYUyQqFMcaYInmcDlAWGjRooJGRkU7HMMaYCmPNmjWZqhpa0LpKWSgiIyNJTEx0OoYxxlQYIrK1sHV26MkYY0yRrFAYY4wpkhUKY4wxRbJCYYwxpkhWKIwxxhTJCoUxxpgiWaEwxhhTJCsUPsdO5vDais2s3JjpdBRjjPErlfKGu5LwuITXPttMx4i6XBLdwOk4poJp0KAB1huA8UepqalkZl7YH8BWKHw8bhfXdQnn9c+2kHHwOKG1gpyOZCoQ6w3A+Kv4+PgLfg879JTHsLim5OQqH3yb7nQUY4zxG1Yo8ohuWJPOzeoya00aNkSsMcZ4OVYoRCRERBaLSIrvsV4R29YWkW0i8q+yzjUsrikbdh3iu237y/qjTFWXkQELFsBbb8GxY06nMaZQTrYoHgaWqmorYKlvuTB/A1aUR6hrO4VRLcDFrDVp5fFxpqrZvh1uvBEiI6FhQ+jfH0aNgg4dYMkSp9MZUyAnC8VgYJpvfhrwq4I2EpE4oBGwqDxC1a4WwNUxjZmbtJ1jJ3PK4yNNVXHgAFxzDcydCxdfDM8+C8uWeVsVItCvH9x8s7elYYwfcfKqp0aqusM3vxNvMTiDiLiA54BfA33LK9iw+KZ8kLSdRcm7GNSpSXl9rKnMTpyAIUMgORnmz4crrzxz/XffwVNPwcSJsGgRfP01NGvmTFZjzlKmLQoRWSIi6wqYBufdTr1njgs6ezwO+EhVtxXjs8aISKKIJGZc4F9kPaLqE143mFmJdvjJlAJVuOMO76Gl117LXyQAqlWDJ56AxEQ4ehRuuAFOniz/rMYUoEwLhar2VdX2BUxzgF0iEgbge9xdwFv0AO4RkVTg78AtIjKxkM+aoqrxqhofGlrgaH7F5nIJQ7qE8/nGTLbvO3pB72UMjz7qPWH9xBMwenTR23bsCFOmwMqV8Nhj5RLPmHNx8hzFXGCUb34UMOfsDVT1JlVtpqqRwAPAW6pa1EnvUjM0rimq8D+7p8JciC++gAkT4Pbb4S9/Kd5rRo70bv/00/Dxx2Wbz5hicLJQTAT6iUgK3vMPEwFEJF5EXncwFwDN6lene4sQZiXaPRXmAlxyCbz7Lrz0kveEdXFNngwxMd6T2zt2nHt7Y8qQY4VCVbNUtY+qtvIdotrjez5RVW8vYPupqnpPeWYcFt+U1KwjJG7dW54fayoTERg2DAICzu911at7C8yhQ3DLLd7zHMY4xO7MLkL/Do2pEei2k9rGGe3awTPPeE+C2yEo4yArFEWoHuhhQMcw5n+3gyMnsp2OY6qiMWOgeXPvCXFrVRiHWKE4h2HxTTl8IoePvt/pdBRTFQUGwiOPeO+rmD/f6TSmirJCcQ7xzevRokENO/xknHPLLRAVZa0K4xgrFOcgIgyNi+CrLXv4OeuI03FMVRQQ4C0S334LH3zgdBpTBVmhKIbru4TjEphtHQUap9x0E7Ru7b0JLzfX6TSmirFCUQxhdYK5rFUo732TTm6uNf2NAzweb5H4/nt47z2n05gqxgpFMQ2LiyB931E+33hhY88aU2IjRkDbtt6uQOxchSlHViiKqV+7RjSoGcQbn29xOoqpqtxuuP9+WLcOPv/c6TSmCrFCUUzVAtyMvqQ5n27I4KedB5yOY6qqkSOhTh145RWnk5gqxArFefj1xc2pHujmtRXWqvB3CxcupE2bNkRHRzNxYv4Oh1esWEGXLl3weDzMnj37jHXTpk2jVatWtGrVimnTpuV7raNq1PBeLjt7tg1wZMqNFYrzULd6IMPjmzJ3bTo799sYx/4qJyeHu+++mwULFpCcnMz06dNJTk4+Y5tmzZoxdepUbrzxxjOe37NnD48//jhfffUVq1ev5vHHH2fvXj/r62vsWO9ASG++6XQSU0VYoThPt13Wgpxc5c0vrFXhr1avXk10dDRRUVEEBgYycuRI5sw5sxf7yMhIOnbsiMt15j+Bjz/+mH79+hESEkK9evXo168fCxcuLM/45xYTA1dcAa++apfKmnJhheI8NQ2pTv8OYfz3q585eMxGIPNH6enpNG3a9PRyREQE6enFG1fkfF47ZcoU4uPjiY+P50JHVTxvd90FmzfD4sXl+7mmSrJCUQJjroji4PFsZqy2G/CqsjFjxpCYmEhiYiIXOqriebvuOggNtZPaplxYoSiBjhF1uTgqhH9/sYWTOdb09zfh4eGkpf1SxLdt20Z4eHiZv7ZcBQXBbbfB3Lmw7ZxDyhtzQaxQlNDYK1qyY/8x5n233eko5ixdu3YlJSWFLVu2cOLECWbMmMGgQYOK9dqrrrqKRYsWsXfvXvbu3cuiRYu46qqryjhxCY0Z473x7rXXnE5iKjlHCoWIhIjIYhFJ8T3WK2S7HBFJ8k1zyztnURLahNK6UU1e/XSzDZXqZzweD//617+46qqraNu2LcOHDycmJoZHH32UuXO9P0Zff/01ERERzJo1i7FjxxITEwNASEgIjzzyCF27dqVr1648+uijhISEOPl1CteiBVx5JUydaie1TZkSJ37JicgzwB5VnSgiDwP1VPWhArY7pKo1z/f94+PjNTExsTSiFmlWYhp/mP0db93ajStal/MxauNX4uPjKY+fuXymT4cbb4RlyyAhofw/3/i94v5sisgaVY0vaJ1Th54GA6fuZJoG/MqhHBdkUGwTGtYKYsqKzU5HMVXV4MFQsya8/bbTSUwl5lShaKSqO3zzO4FGhWxXTUQSReRLEflVUW8oImN82yaW16WKQR43v7m0BZ9vzOSH7fvL5TOropycHLZv387PP/98ejI+1avDkCHeO7WPHnU6jamkyqxQiMgSEVlXwDQ473bqPfZV2PGv5r6m0I3AJBFpWdjnqeoUVY1X1fjyvFTxxu7NqBHo5jVrVZSJf/7znzRq1Ih+/foxYMAABgwYwLXXXut0LP9y881w4ID3CihjyoCnrN5YVfsWtk5EdolImKruEJEwYHch75Hue9wsIsuBzsCmsshbUnWCAxjZrRlTV6byh6svIrxusNORKpXJkyezfv166tev73QU/5WQAOHh3sNPI0Y4ncZUQk4depoLjPLNjwLmnL2BiNQTkSDffAPgUiD57O38wa2XtQDgTeuCvNQ1bdqUOnXqOB3Dv7nd3hHwFi6E3QX+zWXMBSmzFsU5TATeFZHbgK3AcAARiQfuVNXbgbbAqyKSi7egTVRVvywU4XWDGdgxjOmrf+a3vVtRp3qA05EqjaioKBISEhgwYABBQUGnn//973/vYCo/dPPN8MwzMGMG3Huv02lMJeNIi0JVs1S1j6q2UtW+qrrH93yir0igqitVtYOqdvI9vuFE1uIa27Mlh0/k8Mbndq6iNDVr1ox+/fpx4sQJDh48eHoyZ2nfHmJj7eonUyacalFUOm3DatO/Q2P+/UUqt17WgrrVA52OVCk89thjTkeoOG6+2TsC3k8/wUUXOZ3GVCLWhUcp+l2f1hw+kc1rn1mr4kKNHz8egIEDBzJo0KB8kynADTeAy2WtClPqrEVRito0rsWADmG8+UUqt10WRUgNa1WU1M033wzAAw884HCSCiQsDPr2hf/+FyZMABGnE5lKwloUpWx831YcPZnDqyv86ireCicuLg6Anj170qNHD+rVq0dISAg9evSgZ8+eDqfzYzfeCKmp8OWXTicxlYgVilIW3bAWgzs14a2VW8k8dNzpOBXe/PnzadmyJffeey/33HMP0dHRLFiwwOlY/uu666BaNW+rwphSYoWiDNzbpxXHs3N49VNrVVyo+++/n2XLlrF8+XI+/fRTli1bxn333ed0LP9VuzYMHAgzZ0J2ttNpTCVhhaIMRIXW5Fedw3lr1VZ2HzjmdJwKrVatWkRHR59ejoqKolatWg4mqgBuvBEyMmDpUqeTmErCCkUZubd3K7JzlZetVXFB4uPj6d+/P1OnTmXatGkMHDiQrl278v777/P+++87Hc8/XXMN1KkD//mP00lMJWGFooxENqjBkC7h/Oern9m531oVJXXs2DEaNWrEp59+yvLlywkNDeXo0aN8+OGHzJs3z+l4/ikoCIYOhf/9D44ccTqNqQTs8tgy9NverXj/m3ReWr6RJwa3dzpOhfTmm286HaFiuvFGeOMNmDcPhg93Oo2p4KxFUYaahlRnWHwEM1ansX2fjRVQEg8++CAHDhzg5MmT9OnTh9DQUN555x2nY/m/nj2991XY1U+mFFihKGN394pGUV5cttHpKBXSokWLqF27NvPmzSMyMpKNGzfy7LPPOh3L/7ndMHIkfPQR7N3rdBpTwVmhKGMR9aozomtT3k1MI22PHS8+X9m+Szznz5/PsGHDrMvx83HjjXDypHf0O2MugBWKcnB3r2gEsVZFCVx77bVcdNFFrFmzhj59+pCRkUG1atWcjlUxxMVB69Zgh+rMBbJCUQ7C6gRzY/dmzFqzja1Zh52OU6FMnDiRlStXkpiYSEBAANWrV2fOnHzjXJmCiMAtt8CKFbDFBtUyJWeFopyMS2iJxyVMXpridJQKJyQkBLfbDUCNGjVo3Lixw4kqEF/nitajrLkQVijKScPa1Rh1SSQffJvOxt028I4pJ82aQa9e8NZboOp0GlNBOVIoRCRERBaLSIrvsV4h2zUTkUUi8qOIJItIZDlHLVVjr4giOMDNP5ZYq8KUo1tugU2bYOVKp5OYCsqpG+4eBpaq6kQRedi3/FAB270FPKmqi0WkJpBbniFLW/2aQdx6WQv++clGxiXsJ6aJXcFTmG+++abI9V26dCmnJJXAkCFw993eVsWllzqdxlRAThWKwUCCb34asJyzCoWItAM8qroYQFUPlWO+MnP75VFMW5nKPxZv4PVRXZ2O47fuv//+QteJCJ988kk5pqngatWC66/39ig7ebK3G3JjzoNThaKRqu7wze8EGhWwTWtgn4i8D7QAlgAPq2pOQW8oImOAMQDNmjUr/cSlpE5wAGOuiOLvizbw7c976dyswKNuVd6yZcucjlC5jBrlvUx27lzr0sOctzIrFCKyBCjo8pQ/511QVRWRgs6yeYDLgc7Az8BMYDTwRkGfp6pTgCkA8fHxfn3WbvSlLfj3F6k8v3gDb9/W3ek4fm/dunUkJydz7NgvnSvecsstDiaqgHr1gvBw7+EnKxTmPJVZoVDVvoWtE5FdIhKmqjtEJAzYXcBm24AkVd3se80HwMUUUigqkppBHu7q2ZInP/qRLzdncXFUfacj+a3HH3+c5cuXk5ycTP/+/VmwYAGXXXaZFYrz5XZ7L5V99lnYtQsaFdSIN6ZgTl0eOxcY5ZsfBRR0B9XXQF0RCfUt9waSyyFbufj1xc1pWCuI5xdtQO2yxULNnj2bpUuX0rhxY958803Wrl3L/v37z/m6hQsX0qZNG6Kjo5k4cWK+9cePH2fEiBFER0fTvXt3UlNTAUhNTSU4OJjY2FhiY2O58847S/srOeeWWyAnx+6pMOfNqUIxEegnIilAX98yIhIvIq8D+M5FPAAsFZHvAQFecyhvqQsOdHNP72hWp+7hs5RMp+P4reDgYFwuFx6PhwMHDtCwYUPS0tKKfE1OTg533303CxYsIDk5menTp5OcfObfGG+88Qb16tVj48aN3HfffTz00C/XUrRs2ZKkpCSSkpJ45ZVXyuR7OaJtW+9VT6++CrkV+gJCU84cKRSqmqWqfVS1lar2VdU9vucTVfX2PNstVtWOqtpBVUer6gkn8paVEV2bEl43mOcWrbdWRSHi4+PZt28fd9xxB3FxcXTp0oUePXoU+ZrVq1cTHR1NVFQUgYGBjBw5Ml+3H3PmzGHUKG+jdujQoSxdurRq/D8YNw42boQlS5xOYioQuzPbQUEeN/f2iWbttv0s+bGg0zTmpZdeom7dutx5550sXryYadOmnXMwo/T0dJo2bXp6OSIigvT09EK38Xg81KlTh6ysLAC2bNlC586d6dmzJ5999lmhnzNlyhTi4+OJj48nIyOjpF+xfA0ZAqGh8OKLTicxFYgVCodd3yWCyPrVeW7RenJzq8BftOepT58+p+cjIyPp2LHjGc+VtrCwMH7++We+/fZbnn/+eW688UYOHDhQ4LZjxowhMTGRxMREQkNDC9zG7wQFwR13eEe+27rV6TSmgrBC4bAAt4vxfVvz086DfLRux7lfUEUcO3aMPXv2kJmZyd69e9mzZw979uwhNTU1X+vgbOHh4Wecx9i2bRvh4eGFbpOdnc3+/fupX78+QUFB1K/vvQotLi6Oli1bsmHDhlL+dg4bM8b7OGWKszlMhWGFwg8M7NSEVg1r8o/FG8ixVgUAr776KnFxcfz000906dKFuLg44uLiGDx4MPfcc0+Rr+3atSspKSls2bKFEydOMGPGDAYNGnTGNoMGDWLatGmA98qq3r17IyJkZGSQk+O9p3Pz5s2kpKQQFRVVNl/SKc2bw7XXwuuvw/HjTqcxFYGqVropLi5OK5oP16Zr84fm6ZykdKej+JUXXnihRK+bP3++tmrVSqOionTChAmqqvrII4/onDlzVFX16NGjOnToUG3ZsqV27dpVN23apKqqs2fP1nbt2mmnTp20c+fOOnfu3GJ9XoX7mVu4UBVU//tfp5OYMlbcn00gUQv5nSpaCa/0iI+P18TERKdjnJfcXOWqSSsAWDj+CtwucTiRfzhx4gSvvPIKK1Z4901CQgJjx44lICDA4WRnio+Pp0L9zOXmeke/CwuDIk7Ym4qvuD+bIrJGVeMLWmeHnvyEyyXc26cVKbsP8dH3dq7ilHHjxrFmzRrGjRt3ev6uu+5yOlbF53LBXXfB55/D2rVOpzF+zgqFH+nfIYzohjX55ycpVf4KqOzsbAC+/vprpk2bRu/evenduzdvvvkmX3/9tcPpKonf/AZq1oQC7lw3Ji8rFH7E7RJ+2zuaDbsOsWDdTqfjOKpbt24AuN1uNm3adPr5zZs3nx4W1VygkBDvOBUzZ8KPPzqdxvgxKxR+5tqOTWgZWoPJSzdU6VbFqXNnf//73+nVqxcJCQkkJCTQu3dvnnvuOYfTVSL33w/BwTBhgtNJjB9zajwKUwi371zF72YksfCHnfTvEOZ0JEdkZGTw/PPPAzB27NjTl6y63W6+/fZbevXq5WS8yiM01NuqeO45ePRRaNPG6UTGD1mLwg9d27EJUaE1eGFp1T1XkZOTw6FDhzh48CDZ2dmnL9PLzs7m4MGDTserXB54wHvH9pNPOp3E+ClrUfght0u4t3crxs9M4uMfdnJNFWxVhIWF8eijjzodo2po2NB7BdSkSfDII9CqldOJjJ+xFoWfGtipCVENajC5irYqKuP9PX7tD3+AwEBrVZgCWaHwU26X8Ns+0fy08yCLknc5HafcLV261OkIVUvjxnDnnd5xtStb31bmglmh8GMDOzahRRVtVYSEhDgdoep56CGoUcN7cttadCYPKxR+zON28dve0fy44wCLf6x6rQpTzho3hqee8g5qNH2602mMH3GkUIhIiIgsFpEU32O9ArbpJSJJeaZjIvIrB+I6alCnJkTWr87kJSl23N6UvTvvhK5d4b77YO9ep9MYP+FUi+JhYKmqtgKW+pbPoKrLVDVWVWOB3sARYFG5pvQDHreLe3q3InnHARZXwXMVppy53fDKK5CZCX/8o9NpjJ8oVqEQkedEJKYUP3cwMM03Pw341Tm2HwosUNUjpZihwvhVbBOa16/Oi8s2WqvClL0uXeDee+HVV2HVKqfTGD9Q3BbFj8AUEflKRO4UkToX+LmNVPVUF6k7gUbn2H4kUORBUxEZIyKJIpJYYcYvLiaP28WdPVuydtt+vtiY5XQcUxU88QRERMDYsXDihNNpjMOKVShU9XVVvRS4BYgEvhOR/4pIof0oiMgSEVlXwDT4rPdWoNA/k0UkDOgAfHyOjFNUNV5V4yvM+MXn4fou4TSqHcSLyzY6HcVUBbVqwb/+Bd9/7z1fYaq0Yp+jEBE3cJFvygTWAr8XkRkFba+qfVW1fQHTHGCXrwCcKgS7i/jo4cD/VPVkcbNWRkEeN3dcHsWqzVms2WonGU05GDzY273HSy95h001VVZxz1H8A1gP9AeeUtU4Vf0/VR0IdC7B584FRvnmRwFzitj2Bs5x2KmquKFbM+pVD+Dl5daqMOVk4kS48koYNw5WrnQ6jXFIcVsU3wGdVHWsqq4+a123EnzuRKCfiKQAfX3LiEi8iJz+00VEIoGmwKcl+IxKp0aQh99c2oIlP+7mxx0HnI5jqgK323tPRdOmMGQIbN/udCLjgOIWirVAGxHpkmdqKSIeVd1/vh+qqlmq2kdVW/kOUe3xPZ+oqrfn2S5VVcNVNfd8P6OyGtUjkhqBbl5avuncGxtTGkJCYM4cOHgQrrvO+2iqlOIWipeAL4EpwGvAKmAWsF5EriyjbKYAdaoH8OsezZn/3Xa2ZB52Oo6pKtq3h//8B9asgd69vfdZmCqjuIViO9DZd1VRHN7zEpuBfsAzZRXOFOy2y1rgcbt49VNrVZhyNHgwfPABrFsHl18OaWlOJzLlpLiForWq/nBqQVWTgYtUdXPZxDJFaVirGiPim/LeN9vYsf+o03FMVXLttfDxx95zFZddZj3NVhHFLRTJIvKyiPT0TS/5ngsCqvRlq04Zc0UUuQqvrdjidBRT1VxxBSxfDkePQo8e8O67TicyZay4hWIUsBEY75s2A6PxFgkbvNgBTUOqMzi2CdNX/0zWoeNOxzFVTefO3stlW7aEESPghhsgy3oNqKzOWSh8N9p9pKrPqep1vunvqnpEVXNV9VA55DQFGJfQkmPZOUxdmep0FFMVRUd7i8Xf/gazZ3tPeH/wgY1lUQmds1Coag6QWwr9O5lSFt2wFle1a8zUlakcPGZHAI0DPB74y19g9Wpo0MB7+Wy3bvDhh1YwKpHiHno6BHwvIm+IyAunprIMZopnXK+WHDyWzTtf/ux0FFOVde4M33zj7eojKwsGDYL4eJg503suw1RoxS0U7wOPACuANXkm47COEXW5vFUD3vh8M8dO5jgdx1RlAQFw222wfj38+9+wfz+MHAmNGsFvfgNLl0KO/YxWRMXtPXYa8C7wpapOOzWVbTRTXHf3iibz0AneTbTr2o0fCAjwFob1673Dqg4dCu+9B337/nJ46oUXvPdj5FqnCxVBcTsFHAgkAQt9y7EiMrcMc5nz0L1FCF2a1eXVTzdzMsf+4Z2ycOFC2rRpQ3R0NBMnTsy3/vjx44wYMYLo6Gi6d+9Oamrq6XVPP/000dHRtGnTho8/LrKHe1MYtxv69PG2Lnbt8l5GO2QIfPcd/O530KED1KkDl1wCd90FL7/sbXVs3gwn7ZybP/EUc7u/4u38bzmAqiaJSFQZZTLnSUQYlxDN7W8lMu+77VzXOcLpSI7Lycnh7rvvZvHixURERNC1a1cGDRpEu3btTm/zxhtvUK9ePTZu3MiMGTN46KGHmDlzJsnJycyYMYMffviB7du307dvXzZs2IDb7XbwG1VwwcEwbJh3Ati6FZYt857XWLsWZszwDsF6itvtHTipSRNo3Ng7NWoE9er9MtWt6x03o0YNqFnT+xgc7D3BbkpVcffoSVXdLyJ5n7M/Xf1I74sa0rpRTV5evonBncJxueTcL6rEVq9eTXR0NFFR3r9nRo4cyZw5c84oFHPmzOGvf/0rAEOHDuWee+5BVZkzZw4jR44kKCiIFi1aEB0dzerVq+nRo0eJ8yQkJOR7bvjw4YwbN44jR47Qv3//fOtHjx7N6NGjyczMZOjQofnW33XXXYwYMYK0tDRuvvnmfOvvv/9+Bg4cyPr16xk7dmy+9X/5y1/o27cvSUlJjB8/Pt/6p556iksuuYSVK1fypz/9Kd/6SZMmERsby5IlS5gwYUK+9a+++ipt2rThww8/5Lnnnsu3/u2336bp6NHMnDmTl196idDjxwk/epSwY8cIO3aM4V27Um3vXvauXo1r1y7qZGfne4+C5ADH3W6q16sHgYHsP3qU/ceOkSNCjgjZLhfqctG2fXtwu9mSlkbW3r3kipArggKewEC6dusGLhfrkpPJ2rMHBdS3Prh6de/Pgwhr1qxhz759pz9fgdq1a3Ox7+dl5ZdfcuDAmb09h4SE0K17dwBWrFjB4cOH0Ty/XxuGhhIfHw/A0k8+4fjxM++VCmvShM6xsQB8/PHHXPXPf3oP6ZWR4haKH0TkRsAtIq2AewHrnN6PuFzCXQktuW/mWj75aTd9251rdNnKLT09naZNm55ejoiI4Kuvvip0G4/HQ506dcjKyiI9PZ2LL774jNemp6fn+4wpU6YwZcoUACrb8LvlToSMatXIqFaNJN9T/V9+mWoNGjBn6lSmTp2KW5Wa2dne6eRJJj/+ONWys1nywQd8t2oVwTk5BObmEpibS1BODsMHDYKTJ9n27bfs/PlnPKq4fVOgCAQFQU4O7pwcquXm4lJFAJcqATk53laPKg327qXGkSMIeCdVAo4d83aQCDTbvZvGZ/0iDzx6FD7/HFRpnZlJ9lmH0gKPHIHD3k4922VmknPWSf6gI0fAV3w6ZmaSm+dcjpxav9s73luXrCxISSmF/wlFUNVzTkB14EngayDRN1+tOK91YoqLi9Oq6ER2jl7y9FK97sXPNTc31+k4jpo1a5bedtttp5ffeustvfvuu8/YJiYmRtPS0k4vR0VFaUZGht5999369ttvn37+1ltv1VmzZhX5eVX1Z874v+L+bAKJWsjv1OJe9XREVf+sql3V24Psn1X1WFkVL1MyAW4XY3tG8c3P+1i9ZY/TcRwVHh5OWp7eTbdt20Z4eHih22RnZ7N//37q169frNcaU5UU96qn1iIyRUQWicgnp6ayDmfO37C4ptSvEcjLVbwL8q5du5KSksKWLVs4ceIEM2bMYNCgQWdsM2jQIKZN817lPXv2bHr37o2IMGjQIGbMmMHx48fZsmULKSkpdOtWkoEcjakcinuOYhbwCvA63nNFF0xEQoCZQCSQCgxX1b0FbPcMMABvUVsM/M7XTDIFCA50c+tlLXj24/X8sH0/MU2qZs8rHo+Hf/3rX1x11VXk5ORw6623EhMTw6OPPkp8fDyDBg3itttu4+abbyY6OpqQkBBmzJgBQExMDMOHD6ddu3Z4PB5efPFFu+LJVGlSnN+5IrJGvQMWld4HewvAHlWdKCIPA/VU9aGztrkEeBa4wvfU58AfVXV5Ue8dHx+viYmJpRm3Qtl/9CSXTvyEXhc15J83dHY6TpXQoEEDIiMjC12fkZFBaGho+QU6D5atZCpKttTUVDKLMSKh7/d8fEHritui+FBExgH/A06f3lffWNclNBhI8M1Pw3uPxkNnbaNANSAQ78n+AGDXBXxmlVAnOICbLm7Gays2c3+/1kQ2qOF0pErvXP8Q4+Pj8dc/XixbyVSlbOczHsUf8F4Se6qfpwtN0UhVd/jmdwL5rudU1VXAMmCHb/pYVX8s6M1EZIyIJIpIol2qCLdd6h0udcpnNgihMebCFKtFoaotSvLmIrIEaFzAqj+f9f4qIvmOgYlINNAWOHWr8WIRuVxVPysg4xRgCngPPZUkb2XSsHY1hsZFMDtxG+P7tKJh7WpORzLGVFBFtihE5ME888POWvfUud5cVfuqavsCpjnALhEJ871XGLC7gLe4Dm9HhIfUO0DSAqDkt8dWMWOviCI7N5c3vrDhUp02ZswYpyMUyrKVTFXKVuTJbBH5RlW7nD1f0PJ5f7DIs0BWnpPZIar64FnbjADuAK7Ge45iITBJVT8s6r2r+snsvH47/VuW/bSbLx7uTZ3gAKfjGGP8VFEns891jkIKmS9o+XxNBPqJSArQ17eMiMSLyOu+bWYDm4DvgbXA2nMVCXOmO3tGceh4Nu98udXpKMaYCupc5yi0kPmCls+LqmYBfQp4PhG43TefA+TvzcwUW0yTOiS0CeXfn2/h1ktbEBxo9wMYY87PuVoUnUTkgIgcBDr65k8tdyiHfKYUjEuIJuvwCaavtuFSy9u5xsQoT7feeisNGzakffv2p5/bs2cP/fr1o1WrVvTr14+9e/Pd81ou0tLS6NWrF+3atSMmJobJkyf7Tb5jx47RrVs3OnXqRExMDI899hgAW7ZsoXv37kRHRzNixAhOnDhR7tlOycnJoXPnzlx77bVlkq3IQqGqblWtraq1VNXjmz+1bAe8K4huLULo3iKEV1dssuFSy9GpMTEWLFhAcnIy06dPJzk52bE8o0ePZuHChWc8N3HiRPr06UNKSgp9+vRxrJh5PB6ee+45kpOT+fLLL3nxxRdJTk72i3xBQUF88sknrF27lqSkJBYuXMiXX37JQw89xH333cfGjRupV68eb7zxRrlnO2Xy5Mm0bdv29HKpZyust8CKPFlPnvl9npKhzR+ap2+tSnU6SpWxcuVKvfLKK08vP/XUU/rUU085mEh1y5YtGhMTc3q5devWun37dlVV3b59u7Zu3dqpaGcYNGiQLlq0yO/yHT58WDt37qxffvml1q9fX0+ePKmq+f9fl6e0tDTt3bu3Ll26VAcMGKC5ubklysaF9h5rKr5LWtYnrnk9Xl62kRPZNuZUeShoTIyCxrVw0q5duwgLCwOgcePG7NrlfMcHqampfPvtt3Tv3t1v8uXk5BAbG0vDhg3p168fLVu2pG7dunh8o+k5+f92/PjxPPPMM7hc3l/nWVlZpZ7NCkUVISL8tnc02/cf4/1vtjkdx/ghEeGsUSzL3aFDhxgyZAiTJk2idu3aZ6xzMp/b7SYpKYlt27axevVqfvrpJ0dynG3evHk0bNiQuLhS7YovHysUVUjP1qF0jKjDi8s3cjLHWhVlrSKMa9GoUSN27PD2pLNjxw4aNmzoWJaTJ08yZMgQbrrpJq6//nq/ywdQt25devXqxapVq9i3bx/ZvuFZnfp/+8UXXzB37lwiIyMZOXIkn3zyCb/73e9KPZsViipERLi3dyvS9hxlTtJ2p+NUesUZE8NpecfkmDZtGoMHD3Ykh6py22230bZtW37/+9/7Vb6MjAz2+YYlPXr0KIsXL6Zt27b06tWL2bNnO5rt6aefZtu2baSmpjJjxgx69+7Nf/7zn9LPVtjJi4o82cnswuXm5uo1k1ZowrPLNDunag+XWh7mz5+vrVq10qioKJ0wYYKjWUaOHKmNGzdWj8ej4eHh+vrrr2tmZqb27t1bo6OjtU+fPpqVleVIts8++0wB7dChg3bq1Ek7deqk8+fP94t8a9eu1djYWO3QoYPGxMTo448/rqqqmzZt0q5du2rLli116NCheuzYsXLPlteyZct0wIABJc5GESezizUeRUVjXXgUbcH3O7jrP98weWQsg2P961CIMcYZF9KFh6mEroppTOtGNfnXJxvJza18fygYY0qXFYoqyOUS7undipTdh1j4w06n4xhj/JwViipqQIcwokJr8MLSFGtVGGOKZIWiinK7hHt6RfPTzoMs+dH5m6yMMf7LCkUVNqhTE5rXr84Ln6RQGS9qMMaUDisUVZjH7WJcQkvWpR9g+XobZ9yUjSeffJKYmBg6duxIbGwsX331FZMmTeLIkSNORzPFZIWiiruucwThdYOZvNRaFab0rVq1innz5vHNN9/w3XffsWTJEpo2bWqFooKxQlHFBXpc3N0rmqS0fSxOtnMVpnTt2LGDBg0aEBQUBECDBg2YPXs227dvp1evXvTq1QuARYsW0aNHD7p06cKwYcM4dOgQAJGRkTz44IN06NCBbt26sXHjRgBmzZpF+/bt6dSpE1dccYUzX64KcaRQiEiIiCwWkRTfY71Ctvs/EVnnm0aUd86qYnh8BFENavDMx+vJtj6gTCm68sorSUtLo3Xr1owbN45PP/2Ue++9lyZNmrBs2TKWLVtGZmYmEyZMYMmSJXzzzTfEx8fz/PPPn36POnXq8P3333PPPfcwfvx4AJ544gk+/vhj1q5dy9y5cx36dlWHUy2Kh4GlqtoKWOpbPoOIDAC6ALFAd+ABEal99nbmwnncLh68ug0bdx9i9hrrWdaUnpo1a7JmzRqmTJlCaGgoI0aMYOrUqWds8+WXX5KcnMyll15KbGws06ZNY+vWX8Z4v+GGG04/rlq1CoBLL72U0aNH89prr5GTY4NxlbVzjZldVgYDCb75acBy4KGztmkHrFDVbCBbRL4DrgbeLaeMVcpVMY3p3Kwu/1iygcGx4Ta2tik1brebhIQEEhIS6NChw+lO/k5RVfr168f06dMLfH3ersVPzb/yyit89dVXzJ8/n7i4ONasWUP9+vXL7ktUcU61KBqp6g7f/E6gUQHbrAWuFpHqItIA6AU0LWA7AERkjIgkikhiRoZdwXO+RIQ/XtOWXQeO8+8vtjgdx1QS69evJyUl5fRyUlISzZs3p1atWhw8eBCAiy++mC+++OL0+YfDhw+zYcOG06+ZOXPm6ccePXoAsGnTJrp3784TTzxBaGjoGd25m9JXZi0KEVkCNC5g1Z/zLqiqiki+y21UdZGIdAVWAhnAKqDQNqaqTgGmgLdTwAuIXmV1axFC37YNeWX5Jm7o1oyQGoFORzIV3KFDh/jtb3/Lvn378Hg8REdHM2XKFKZPn87VV199+lzF1KlTueGGGzh+/DgAEyZMoHXr1gDs3buXjh07EhQUdLrV8Yc//IGUFO+Ven369KFTp06OfceqwJHeY0VkPZCgqjtEJAxYrqptzvGa/wLvqOpH53p/6z225DbsOsjVk1bwm0tb8Mi17ZyOY6q4yMhIEhMTadCggdNRKj1/7D12LjDKNz8KmHP2BiLiFpH6vvmOQEdgUbklrKJaN6rF0LgI3l61lbQ9dp27Mca5QjER6CciKUBf3zIiEi8ir/u2CQA+E5FkvIeUfu07sW3K2H39WiMCzy/ecO6NjSlDqamp1prwA45c9aSqWUCfAp5PBG73zR/De+WTKWdhdYL5zaUteHXFJm6/vAUxTeo4HckY4yC7M9sU6K6EltSuFsDEBT85HcUY4zArFKZAdYIDuKdXNJ+lZPJ5SqbTcYwxDrJCYQp1c4/mhNcNZuLCH21wI2OqMCsUplDVAtzcf2Vr1qUfYN73O879AmNMpWSFwhRpcGw4FzWuxf8t+InDx+2iM2OqIisUpkhulzDhV+1J33eU5xbZ5bLGVEVWKMw5xUeG8OuLmzF15RbWpu1zOo4xppxZoTDF8uDVFxFaK4iH3/+ekzZmhTFVihUKUyy1qwXw+KD2/LjjAK9/Zr3LGlOVWKEwxXZ1+8ZcFdOISUs2kJp52Ok4xphyYoXCnJfHB7Un0O3izx98jxM9Dxtjyp8VCnNeGtepxoPXXMQXG7N475t0p+MYY8qBFQpz3m7q1oy45vWYMD+ZzEPHnY5jjCljVijMeXO5hInXd+Dw8WwmzEt2Oo4xpoxZoTAl0qpRLe5KiOaDpO18usHGKDemMrNCYUrs7l4taRlagz//73uOnLDuPYyprKxQmBIL8rh5+vqObNt7lKc/snErjKmsHCkUIjJMRH4QkVwRKXAwb992V4vIehHZKCIPl2dGUzzdWoRw+2UtePvLrfzv221OxzHGlAGnWhTrgOuBFYVtICJu4EXgGrxDot4gIjY0qh966JqL6NYihD++/z3J2w84HccYU8ocKRSq+qOqrj/HZt2Ajaq6WVVPADOAwWWfzpyvALeLF2/sQp3gAMa+k8i+IyecjmSMKUX+fI4iHEjLs7zN91yBRGSMiCSKSGJGhl2FU95CawXx8q/j2Ln/GONnJtmIeMZUImVWKERkiYisK2Aqk1aBqk5R1XhVjQ8NDS2LjzDn0KVZPR4bGMPy9RlMWpridBxjTCnxlNUbq2rfC3yLdKBpnuUI33PGj93UvRlJaft4YWkKnSLq0KdtI6cjGWMukD8fevoaaCUiLUQkEBgJzHU4kzkHEe+IeO3DazN+ZpL1MmtMJeDU5bHXicg2oAcwX0Q+9j3fREQ+AlDVbOAe4GPgR+BdVf3Bibzm/FQLcPPyTXG4XcLYt9fYzXjGVHBSGbuKjo+P18TERKdjVHkrNmQw6s3VDOzYhMkjYxERpyMZYwohImtUtcD72vz50JOp4K5oHcoDV7Zh7trtPL94g41fYUwFVWYns40BuKtnS37OOsI/P9lIrioPXNnGWhbGVDBWKEyZcrmEp6/vgMslvLhsE9m5ysNXX2TFwpgKxAqFKXMul/Dkr9rjcQmvfrqZnBzlzwPaWrEwpoKwQmHKhcslPDE4BrdLeP3zLWTnKo8NbGfFwpgKwAqFKTciwmMD2+F2CW98voXs3FyeGNQel8uKhTH+zAqFKVciwl8GtMXj9h2GyoUnf2XFwhh/ZoXClDsR4eGrL8LjO8Gdk5vLU9d1wOO2q7WN8UdWKIwjRIQHrmyD2+XihaUpbMk8zOSRnWlSN9jpaMaYs9ifcMYxIsLv+7Vm8shYkrcfoP8Ln7E4eZfTsYwxZ7FCYRw3ODacefdeTkS9YO54K5G/zv2B49k5TscyxvhYoTB+oUWDGrx31yWMviSSqStTGfLySrZYz7PG+AUrFMZvBHnc/HVQDK/dEk/anqNc+8JnzEmyIUiMcZoVCuN3+rVrxEe/u5y2YbX53Ywk7n93LRkHjzsdy5gqywqF8UvhdYOZMeZi7ukVzQdJ6fR8dhnPL97AoeM2toUx5c0KhfFbHreLB65qw+L7riChTSgvLE0h4dllvLUqlRPZuU7HM6bKcGqEu2Ei8oOI5IpIgQNl+Lb7t4jsFpF15ZnP+Jeo0Jq8dFMc/xt3CS1Da/LonB/o949Pmffddhvjwphy4FSLYh1wPbDiHNtNBa4u8zSmQujcrB4zxlzMm6O7Us3j5p7/fsvgF79g2frd5OZawTCmrDhyZ7aq/gics+dQVV0hIpHlkclUDCJCr4sackXrUN7/ZhvPL97Ab978mvC6wQyPb8qw+Ai7u9uYUmZdeJgKye0ShsU3ZVBsExb9sIuZX6fxjyUbmLx0Az1bhzKiazP6tG1IgPUfZcwFK7NCISJLgMYFrPqzqs4pg88bA4wBaNasWWm/vfFTQR43Azs1YWCnJvycdYRZa9J4NzGNO99ZQ4OagQyJi2BgxybENKltY18YU0Li5MlAEVkOPKCqiUVsEwnMU9X2xX3f+Ph4TUws9C1NJZedk8uKlAxmrE5j6U+7yclVGtQM4orWDejZOpTLW4USUiPQ6ZjG+BURWaOqBV5cZIeeTKXjcbvofVEjel/UiIyDx1mxIYNPN2Sw7KfdvP9NOiLQMbwOPVuHckXrUDpE1CHI43Y6tjF+y5EWhYhcB/wTCAX2AUmqepWINAFeV9X+vu2mAwlAA2AX8JiqvnGu97cWhSlITq7yffp+Pl2fwYqUDL79eS+5CgFuoXWjWnQIr0P78Dp0CK/DRWG1rHiYKqWoFoWjh57KihUKUxz7j5xk5aZM1m7bz7r0/Xyfvp/9R08C4HF5i0f78NpEhdYksn4NIhtUp3lIDYIDrYCYyscKhTHFoKps23uU731FY136fpK3HyDr8IkztmtcuxrN61cnsn4NmtWvTqPa1WhYK+j0Y93qAXbi3FQ4do7CmGIQEZqGVKdpSHX6dwg7/fz+oyf5OesIqVmH2Zp1mC2ZR9iadZilP+0m81D+zgoD3S5CawXRqHYQDWp6C0fd6oHUCQ6gbvUA6lUPpG5wAHWqB1C7WgDVA93UCPIQ5HFZgTF+yQqFMedQJziADhF16BBRJ9+6oydy2H3wGLsPHmfXgWPsOnDcu3zAu7w16whrt51g75GT5+yfyuMSagR5qBnkOV08ggPcVAtwUS3A7ZtcBHl+mQ/0uAh0uwhwe+dPPQa6hQC3C4/bhccl3skteFwu3C7vOrcL3C4XbhFcLu+9KW4R76NLcLkElwguAZcIIni39c1bUas6rFCcJSEhId9zw4cPZ9y4cRw5coT+/fvnWz969GhGjx5NZmYmQ4cOzbf+rrvuYsSIEaSlpXHzzTfnW3///fczcOBA1q9fz9ixY/Ot/8tf/kLfvn1JSkpi/Pjx+dY/9dRTXHLJJaxcuZI//elP+dZPmjSJ2NhYlixZwoQJE/Ktf/XVV2nTpg0ffvghzz33XL71b7/9Nk2bNmXmzJm8/PLL+dbPnj2bBg0aMHXqVKZOnZpv/UcffUT16tV56aWXePfdd/OtX758OQB///vfmTdv3hnrgoODWbBgAQB/+9vfWLp06Rnr69evz3vvvQfAH//4R1atWnXG+oiICN555x0Axo8fT1JS0hnrW7duzZQpUwAYM2YMGzZsOGN9bGwskyZNAuDXv/4127ZtO2N9jx49ePrpp2levwZDhgwhKyvrjPV9+vThkUceAeCq/tdy+CTkeqqR66lGjieYjnHduLxXXw6fyGHKv6eR6w5knzuQva5Act2B1G/UmHoNGnLkeDY/b9uOujzkujyoywMu5//5Ct4Co6rk5GQjKOQ5nF0tKAiPx01OdjbHjh0FBd9/AKV2rdp4PB6OHz/GkcOHTz9/6iEkJASPx82RI0c4dOhQntd6NQwNxe12c+jQIQ4dPJgvX1hYY1wuF/v37+fQwUN51njfp2nTpgDs2bOHw4fPHCjLJUJERAQAWVlZHD5y5IzXu91uwsPDAcjIyODo0aOn9wmAx+OhSZMmAOzatYvjx89sfQYGBtK4sfdWs507d3LihO8Qp2//VatWjUa+9enp6WSfPHnG64OrV6dhw4YAbEtLIycnh87tWvPunT3y7YcL5fxPmjFVhEtz8Jw8Cid/+YXVtlobRl/aAoBZj3+W7zXDOw5n3LghBf6Rogi/vmUUI268iV27M/nN7Xeg4kbFDS7v43VDhtAzoTc7du3m6Yn/h4oLFRf4Hvv3H0DH2M5s376DaW+/jSIgLhBBEfr260d0dCvSt29n7ofzQAQQ73sAvXr1JjwinLS0dJZ/+qlvPYCACD179qRBaChbf07j669//GUdgAiXtL+CunXrsjV1K2u3bfrly/nep1tsM2rUqMHmzVv4KX2r93vnacl06dySoGpBbNq4l43b8w9y1alrGzxuN+uPZvDz4R351ncI996e9ePB7Ww/vPOMfG63m5jwGAB+2JtG9uFdvv3uFRQURLuwtgB8l5FKzuHMPN8fqgcHc1HYRQAc27mJvUf2nfHZNV21uCisNQCHtm3g4JEzC12toLq0aVQLgH2p+zly9MgZ6+tUr396fdamfRw/cZyo0M75vmNpsJPZxhhjijyZbR3hGGOMKZIVCmOMMUWyQmGMMaZIViiMMcYUyQqFMcaYIlmhMMYYUyQrFMYYY4pkhcIYY0yRKuUNdyKSAWwt4csbAJmlGKc0WbaSsWwlY9lKpqJma66qoQWtqJSF4kKISGJhdyc6zbKVjGUrGctWMpUxmx16MsYYUyQrFMYYY4pkhSK/KU4HKIJlKxnLVjKWrWQqXTY7R2GMMaZI1qIwxhhTJCsUxhhjimSFwkdErhaR9SKyUUQedjpPXiKSKiLfi0iSiDg+IpOI/FtEdovIujzPhYjIYhFJ8T3W86NsfxWRdN/+SxKR/OPZln2upiKyTESSReQHEfmd73nH91sR2fxhv1UTkdUistaX7XHf8y1E5Cvfv9eZIhLoR9mmisiWPPsttryz5cnoFpFvRWSeb7lk+01Vq/wEuIFNQBQQCKwF2jmdK0++VKCB0zny5LkC6AKsy/PcM8DDvvmHgf/zo2x/BR5weJ+FAV1887WADUA7f9hvRWTzh/0mQE3ffADwFXAx8C4w0vf8K8BdfpRtKjDUyf2WJ+Pvgf8C83zLJdpv1qLw6gZsVNXNqnoCmAEMdjiT31LVFcCes54eDEzzzU8DflWemU4pJJvjVHWHqn7jmz8I/AiE4wf7rYhsjlOvU4OMB/gmBXoDs33PO7XfCsvmF0QkAhgAvO5bFkq436xQeIUDaXmWt+En/1B8FFgkImtEZIzTYQrRSFVPjV6/E2jkZJgC3CMi3/kOTTlyWOwUEYkEOuP9C9Sv9ttZ2cAP9pvv8EkSsBtYjLf1v09Vs32bOPbv9exsqnpqvz3p22//EJEgJ7IBk4AHgVzfcn1KuN+sUFQMl6lqF+Aa4G4RucLpQEVRb7vWb/6yAl4GWgKxwA7gOaeCiEhN4D1gvKoeyLvO6f1WQDa/2G+qmqOqsUAE3tb/RU7kKMjZ2USkPfBHvBm7AiHAQ+WdS0SuBXar6prSeD8rFF7pQNM8yxG+5/yCqqb7HncD/8P7j8Xf7BKRMADf426H85ymqrt8/6BzgddwaP+JSADeX8T/UdX3fU/7xX4rKJu/7LdTVHUfsAzoAdQVEY9vleP/XvNku9p3KE9V9TjwJs7st0uBQSKSivdQem9gMiXcb1YovL4GWvmuCAgERgJzHc4EgIjUEJFap+aBK4F1Rb/KEXOBUb75UcAcB7Oc4dQvYp/rcGD/+Y4PvwH8qKrP51nl+H4rLJuf7LdQEanrmw8G+uE9h7IMGOrbzKn9VlC2n/IUfsF7DqDc95uq/lFVI1Q1Eu/vs09U9SZKut+cPivvLxPQH+/VHpuAPzudJ0+uKLxXYa0FfvCHbMB0vIciTuI9znkb3uOfS4EUYAkQ4kfZ3ga+B77D+4s5zIFcl+E9rPQdkOSb+vvDfisimz/st47At74M64BHfc9HAauBjcAsIMiPsn3i22/rgHfwXRnl1AQk8MtVTyXab9aFhzHGmCLZoSdjjDFFskJhjDGmSFYojDHGFMkKhTHGmCJZoTDGGFMkKxTGXAAR+bOv59DvfD2FdheR8SJS3elsxpQWuzzWmBISkR7A80CCqh4XkQZ4ex9eCcSraqajAY0pJdaiMKbkwoBM9XbVgK8wDAWaAMtEZBmAiFwpIqtE5BsRmeXrU+nUOCPPiHeskdUiEu17fpiIrPONc7DCma9mzC+sRWFMCfl+4X8OVMd7V/VMVf3U179OvKpm+loZ7wPXqOphEXkI792wT/i2e01VnxSRW4DhqnqtiHyPt8+gdBGpq95+hIxxjLUojCkh9Y5FEAeMATKAmSIy+qzNLsY7CNAXvu6oRwHN86yfnuexh2/+C2CqiNyBd1AtYxzlOfcmxpjCqGoOsBxY7msJjDprE8E7TsENhb3F2fOqeqeIdMc76MwaEYlT1azSTW5M8VmLwpgSEpE2ItIqz1OxwFbgIN4hRQG+BC7Nc/6hhoi0zvOaEXkeV/m2aamqX6nqo3hbKnm7wDem3FmLwpiSqwn809fVdDbeHjnHADcAC0Vku6r28h2Omp5npLO/4O2pGKCeiHwHHPe9DuBZXwESvD3Lri2PL2NMYexktjEOyXvS2+ksxhTFDj0ZY4wpkrUojDHGFMlaFMYYY4pkhcIYY0yRrFAYY4wpkhUKY4wxRbJCYYwxpkj/D3rgE2GytWidAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n", "\n", "theory_energy = -1.136189454088\n", "theory_spin = 0\n", "\n", "plt.hlines(theory_energy, 0, 39, linestyles=\"dashed\", colors=\"black\")\n", "plt.plot(energies)\n", "plt.xlabel(\"Steps\")\n", "plt.ylabel(\"Energy\")\n", "\n", "axs = plt.gca()\n", "\n", "inset = inset_axes(axs, width=\"50%\", height=\"50%\", borderpad=1)\n", "inset.hlines(theory_spin, 0, 39, linestyles=\"dashed\", colors=\"black\")\n", "inset.plot(spins, \"r\")\n", "inset.set_xlabel(\"Steps\")\n", "inset.set_ylabel(\"Total spin\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ここまでで、Pennylane/Braket パイプラインを使用して、分子の基底状態エネルギーを効率的に見つける方法を学びました！" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## `shots>0` での VQE の実行\n", "\n", "上記の例では、ハミルトニアンの厳密な期待値計算を、それぞれ1回の評価で行うことが可能でした。これはシミュレータが状態にアクセス可能であるためで、QPU のように測定により期待値を計算するデバイスの場合は不可能です。\n", "\n", "電子ハミルトニアン ``h`` の期待値を有限回の測定から求めたいとします。このハミルトニアンは、パウリ作用素のテンソル積である15の個々のオブザーバブルから構成されます。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Number of Pauli terms in h:\", len(h.ops))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "期待値を測定する簡単なアプローチは、回路を15回実装し、毎回ハミルトニアン ``h`` の一部を形成するパウリ項の1つを測定することです。しかし、もっと効率的な方法があるかもしれません。パウリ項は、単一の回路で同時に測定できるグループ（PennyLane の [グループ化](https://pennylane.readthedocs.io/en/stable/code/qml_grouping.html) モジュールを参照）に分けることができます。各グループの要素は、量子ビットごとに交換可能なオブザーバブルとして知られています。ハミルトニアン ``h`` は5つのグループに分けることができます:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "groups, coeffs = qml.grouping.group_observables(h.ops, h.coeffs)\n", "print(\"Number of qubit-wise commuting groups:\", len(groups))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "実際には、これは15の別々の回路を実行する代わりに、5つを実行するだけで済むことを意味します。この節約は、ハミルトニアンのパウリ項の数が増えるにつれて、さらに顕著になります。例えば、より大きな分子または異なる化学的基底集合に切り替えると、量子ビット数と項数の両方が増加する可能性があります。\n", "\n", "幸い、PennyLane/Braket パイプラインには、デバイスの実行回数を最小限に抑えるためにオブザーバブルを事前にグループ化するための機能が組み込まれており、リモートデバイスを使用するときの実行時間とシミュレーション料金の両方を節約できます。このチュートリアルの残りの部分では、最適化されたオブザーバブルのグループ化を使用します。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "グループ化の恩恵を受けるためには、PennyLane が各ハミルトニアンを量子ビットごとの可換なグループに分割するよう指示する必要があります。そのため、グルーピングを計算し、ハミルトニアンと保存します。デバイスは、回路を生成する際にこの情報を使います。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h.compute_grouping()\n", "S2.compute_grouping()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "次に、非ゼロの `shots` を指定してデバイスをインスタンス化します。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dev = qml.device(\"braket.local.qubit\", wires=qubits, shots=1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "コスト関数を再定義します。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wires = dev.wires.tolist()\n", "\n", "@qml.qnode(dev)\n", "def energy_expval(params):\n", " circuit(params, wires)\n", " return qml.expval(h)\n", "\n", "@qml.qnode(dev)\n", "def S2_expval(params):\n", " circuit(params, wires)\n", " return qml.expval(S2)\n", "\n", "def spin(params):\n", " return -0.5 + np.sqrt(1 / 4 + S2_expval(params))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "最後に、VQE の実験を行い期待値を推定します。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "energies, spins = run_vqe(energy_expval, spin, opt, params, iterations)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "厳密な計算と非常に近い値が得られました！この実験では Local Simulator を使いましたが、同じことが QPU やショット数を非ゼロとした他のシミュレータでも実行できます。\n", "\n", "
\n", "次のステップは？ qchem フォルダには、水素分子の異なる原子間距離を表す追加の分子構造ファイルが含まれています。原子間距離の 1 つを選択し、基底状態のエネルギーを求めましょう。原子間距離によって基底状態のエネルギーはどのように変化するでしょう？ \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.13" } }, "nbformat": 4, "nbformat_minor": 4 }