VQE_CN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 变分量子本征求解器\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概览\n",
    "\n",
    "目前普遍认为，量子计算在近期很有前景的一个应用是处理量子化学问题 [1-2]。**变分量子本征求解器** （VQE）作为这个研究方向的核心应用之一，为研究者们提供了可以在目前含噪的中等规模量子设备（NISQ device）上研究量子化学的可能 [1-4]。其核心任务是求解一个量子尺度上封闭物理系统的哈密顿量 $\\hat{H}$ 的基态能量及其对应的量子态。主要的实现方法是通过在量子设备上准备一个参数化的试探波函数 $|\\Psi(\\boldsymbol\\theta)\\rangle$ 然后结合经典机器学习中的优化算法（例如梯度下降法）去不断地调整、优化参数 $\\boldsymbol\\theta$ 使得期望值  $\\langle \\Psi(\\boldsymbol\\theta)|\\hat{H}|\\Psi(\\boldsymbol\\theta)\\rangle$ 最小化。这套方案的基本原理是基于 **Rayleigh-Ritz 变分原理**。 \n",
    "\n",
    "$$\n",
    "E_0 = \\min_{\\boldsymbol\\theta} \\langle \\Psi(\\boldsymbol\\theta)|\\hat{H}|\\Psi(\\boldsymbol\\theta)\\rangle.\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "其中 $E_0$ 表示该系统的基态能量。从数值分析的角度来看，该问题可以被理解为求解一个**离散化**哈密顿量 $H$（埃尔米特矩阵）的最小本征值 $\\lambda_{\\min}$ 和其对应的本征向量 $|\\Psi_0\\rangle$。具体的离散化过程是如何通过建立模型实现的，这属于量子化学的专业领域范畴。精确地解释该过程需要很长的篇幅，这超过了本教程所能处理的范围。我们会在下一节背景知识模块粗略的介绍一下相关知识，感兴趣的读者可以参考 `量子化学: 基本原理和从头计算法`系列丛书 [5]。通常来说，为了能在量子设备上处理量子化学问题，哈密顿量 $H$ 会被表示成为泡利算符 $\\{X,Y,Z\\}$ 的加权求和形式。\n",
    "\n",
    "$$\n",
    "H = \\sum_k c_k ~ \\bigg( \\bigotimes_{j=0}^{M-1} \\sigma_j^{(k)} \\bigg),\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "其中 $c_k$ 表示权重系数, $\\sigma_j^{(k)} \\in \\{I,X,Y,Z\\}$ 并且 $M$ 表示所需的量子比特个数。这样一种哈密顿量的表示形式被称为 **泡利字符串**。以下为一个2量子比特的具体例子，\n",
    "\n",
    "$$\n",
    "H= 0.12~Y_0 \\otimes I_1-0.04~X_0\\otimes Z_1.\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "在下一节，我们会补充一些关于电子结构问题的背景知识。本质上讨论的就是上述哈密顿量 $H$ 是如何计算的。对于熟悉相关背景的读者，或者主要关心如何在量桨上实现 VQE 的读者，请直接跳转至第三节分析氢分子（$H_2$）基态的具体例子。 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 背景： 电子结构问题\n",
    "\n",
    "这里，我们集中讨论下量子化学中的一个基本问题 -- **电子结构问题**。更准确的说，我们关心的是给定分子（molecule）的低位能量本征态。这些信息可以帮助我们预测化学反应的速率和分子的稳定结构等等 [6]。假设一个分子由 $N_n$ 个原子核和 $N_e$ 个电子组成，描述该分子系统总能量的哈密顿量算符 $\\hat{H}_{mol}$ 在一次量子化表示下可以写为，\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{H}_{\\text{mol}} & = -\\sum_{i}\\frac{\\nabla_{R_i}^2}{2M_i} - \\sum_{i} \\frac{\\nabla_{r_i}^2}{2} -\\sum_{i,j}\\frac{Z_i}{\\lvert R_i - r_j\\lvert} + \\sum_{i,j>i}\\frac{Z_iZ_j}{\\lvert R_i - R_j\\lvert} + \\sum_{i, j>i}\\frac{1}{\\lvert r_i - r_j\\lvert}, \n",
    "\\tag{4}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "其中 $R_i、M_i$ 和 $Z_i$ 分别表示第 $i$ 个原子核的位置、质量和原子序数（原子核内质子数），第 $i$ 个电子的位置则表示为 $r_i$。以上公式右边前两项分别代表原子核和电子的总动能。第三项表示带正电的质子和带负电的电子之间的库伦相互吸引作用。最后两项则表示原子核-原子核之间，电子-电子之间的相互排斥作用。这里，分子哈密顿量 $\\hat{H}_\\text{mol}$ 使用的是原子单位制能量 **哈特里能量**（Hartree），记为 Ha。$1$ 哈特里能量的大小为 $[\\hbar^2/(m_ee^2a_0^2)] = 27.2$ 电子伏或 $630$ 千卡/摩尔，其中 $m_e、e$ 和 $a_0$ 分别表示电子质量、基本电荷和玻尔半径。\n",
    "\n",
    "**注释1：** 在处理电子结构问题时，我们不考虑自旋-轨道耦合以及超精细结构。如果出于计算需要，可以作为微扰加入。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 玻恩-奥本海默近似\n",
    "\n",
    "由于原子核的质量要远大于电子，因而在同样的相互作用下电子的运动速度会比原子核快很多。所以，将原子核所处的位置看成固定 $R_i =$常数 是一种合理的近似。这种通过在时间尺度上将电子行为和原子核行为去耦合的近似处理思想被称为玻恩-奥本海默近似。作为近似的直接结果，公式（4）中原子核的动能项会被消去并且表示原子核-原子核相互排斥作用的项可以被认为是一个能量移位（这个项是与电子位置 $r_i$ 无关的）从而也可以作为常数项被忽略。经过这些步骤后，我们可以把哈密顿量近似为：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{H}_{\\text{electron}} & =  - \\sum_{i} \\frac{\\nabla_{r_i}^2}{2} -\\sum_{i,j}\\frac{Z_i}{\\lvert R_i - r_j\\lvert} + \\sum_{i, j>i}\\frac{1}{\\lvert r_i - r_j\\lvert} \n",
    "\\tag{5},\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "在经过以上近似后，分子中多电子结构的能级在理论上可以通过求解以下不含时薛定谔方程获得：\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\text{electron}} |\\Psi_n \\rangle = E_n |\\Psi_n \\rangle,\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "其中 $n$ 指代能级。值得注意的是，电子哈密顿量中电子-电子相互排斥作用的求和项数会随着电子数 $N_e$ 的增多至 $N_e(N_e-1)/2$ 项。这意味着对于一个含有16个电子的氧分子（$O_2$）我们需要计算多达120项的相互排斥作用项。 一般来说，这样的问题是无法从理论上精确求解的。正如狄拉克在 [Quantum mechanics of many-electron systems](https://royalsocietypublishing.org/doi/10.1098/rspa.1929.0094) [7] 所指出的那样，\n",
    "\n",
    "> *The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble.* \n",
    "> \n",
    "> -- Paul Dirac (1929)\n",
    "\n",
    "由于解析的方法太复杂，那么我们可以采用数值方法来处理。一个最简单的数值方法（离散化方法）就是把上述作用中无限维度希尔伯特空间离散化为等间距排开的立方体晶格点。在这样一个离散化的空间里，主要运算规则为复数域的线性代数。假设空间的每个轴都离散为等间距排开的 $k$ 个点，则 $N$-电子（为了方便去掉下标 $e$）的多体波函数可以写为 [2]：\n",
    "\n",
    "$$\n",
    "|\\Psi \\rangle = \\sum_{\\mathbf{x_1}, \\ldots, \\mathbf{x_N}} \\psi(\\mathbf{x_1}, \\ldots, \\mathbf{x_N}) \\mathcal{A}(|\\mathbf{x_1}, \\ldots, \\mathbf{x_N}\\rangle).\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "其中坐标 $|\\mathbf{x_j}\\rangle = |r_j\\rangle |\\sigma_j\\rangle$ 记录第 $j$ 个电子的空间位置信息和自旋，$|r_j\\rangle  = |x_j,y_j,z_j\\rangle$ 且 $j\\in \\{1,2,\\cdots,N\\}$, $x_j,y_j,z_j \\in \\{0,1,\\cdots,k-1\\}$ 同时 $\\sigma_j \\in \\{\\downarrow,\\uparrow\\}$ 表示自旋向下和向上。这样一种离散化方式共计需要 $k^{3N}\\times 2^{N}$ 个数据来表示波函数。在这里，$\\mathcal{A}$ 表示反对称化操作（根据泡利不相容原理）并且 $\\psi(\\mathbf{x_1}, \\mathbf{x_2}, \\ldots, \\mathbf{x_N})=\\langle\\mathbf{x_1}, \\mathbf{x_2}, \\ldots, \\mathbf{x_N}|\\Psi\\rangle$。 可以看出，经典计算机存储这样一个波函数需要的内存是随着电子个数呈指数增长的。这使得基于这种离散化的经典数值方法，无法模拟超过几十个电子的系统。那么，我们是不是能够通过量子设备来存储和准备这样一个波函数然后求解基态能量 $E_0$ 呢？在下一节中，我们将以最简单的分子系统 -- 氢分子（$H_2$）为例，讲解 VQE 算法。\n",
    "\n",
    "**注释2：** 关于量子化学和现有数值计算方法的综述也超过了本教程的处理范围，我们推荐感兴趣的读者去查阅以下经典教材 Helgaker 等人撰写的 *'Molecular Electronic-Structure Theory'* [6] 以及 Szabo & Ostlund 撰写的 *'Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory'* [8]。 如果需要弥补量子计算和量子化学之间知识空缺，请参考以下综述文章 [Quantum chemistry in the age of quantum computing](https://pubs.acs.org/doi/10.1021/acs.chemrev.8b00803) [1] 和  [Quantum computational chemistry](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.92.015003) [2] 。\n",
    "\n",
    "**注释3：** 对于量子化学中的能量计算，我们期望能够达到 **化学精度**（chemical accuracy）$1.6\\times10^{-3}$ Ha 或者 1 千卡/摩尔。\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 氢分子 $H_2$ 基态能量\n",
    "\n",
    "### 构造电子哈密顿量\n",
    "\n",
    "首先，让我们通过下面几行代码引入必要的 library 和 package。量桨的量子化学工具包是基于 `psi4` 和 `openfermion` 进行开发的，所以需要读者先行安装这两个语言包。在进入下面的教程之前，我们强烈建议您先阅读[哈密顿量的构造](./BuildingMolecule_CN.ipynb)教程，该教程介绍了如何使用量桨的量子化学工具包。\n",
    "\n",
    "**注意：关于环境设置，请参考 [README_CN.md](https://github.com/PaddlePaddle/Quantum/blob/master/README_CN.md).**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:13:45.528201Z",
     "start_time": "2021-04-30T09:13:43.385553Z"
    }
   },
   "outputs": [],
   "source": [
    "import paddle\n",
    "import paddle_quantum.qchem as qchem\n",
    "from paddle_quantum.utils import Hamiltonian\n",
    "from paddle_quantum.circuit import UAnsatz\n",
    "\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy\n",
    "from numpy import pi as PI\n",
    "from numpy import savez, zeros\n",
    "\n",
    "# 无视警告\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于具体需要分析的分子，我们需要其**几何构型** (geometry)、**基组**（basis set，例如 STO-3G 基于高斯函数）、**多重度**（multiplicity）以及**分子的净电荷数** (charge) 等多项信息来建模计算出该分子单体积分 (one-body integrations)，双体积分(two-body integrations) 以及哈密顿量等信息。接下来，通过量桨的量子化学工具包将分子的哈密顿量提取出来并储存为 paddle quantum 的 `Hamiltonian` 类，方便我们下一步的操作。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:13:45.545018Z",
     "start_time": "2021-04-30T09:13:45.531302Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FCI energy for H2_sto-3g_singlet (2 electrons) is -1.1372838344855134.\n",
      "\n",
      "The generated h2 Hamiltonian is \n",
      " -0.0970662686176252 I\n",
      "-0.04530261550868938 X0, X1, Y2, Y3\n",
      "0.04530261550868938 X0, Y1, Y2, X3\n",
      "0.04530261550868938 Y0, X1, X2, Y3\n",
      "-0.04530261550868938 Y0, Y1, X2, X3\n",
      "0.1714128263940238 Z0\n",
      "0.16868898168693292 Z0, Z1\n",
      "0.12062523481381847 Z0, Z2\n",
      "0.1659278503225078 Z0, Z3\n",
      "0.17141282639402383 Z1\n",
      "0.1659278503225078 Z1, Z2\n",
      "0.12062523481381847 Z1, Z3\n",
      "-0.22343153674664024 Z2\n",
      "0.17441287610651632 Z2, Z3\n",
      "-0.2234315367466403 Z3\n"
     ]
    }
   ],
   "source": [
    "geo = qchem.geometry(structure=[['H', [-0., 0., 0.0]], ['H', [-0., 0., 0.74]]])\n",
    "# geo = qchem.geometry(file='h2.xyz')\n",
    "\n",
    "# 将分子信息存储在 molecule 里，包括单体积分（one-body integrations），双体积分（two-body integrations），分子的哈密顿量等\n",
    "molecule = qchem.get_molecular_data(\n",
    "    geometry=geo,\n",
    "    basis='sto-3g',\n",
    "    charge=0,\n",
    "    multiplicity=1,\n",
    "    method=\"fci\",\n",
    "    if_save=True,\n",
    "    if_print=True\n",
    ")\n",
    "# 提取哈密顿量\n",
    "molecular_hamiltonian = qchem.spin_hamiltonian(molecule=molecule,\n",
    "                                               filename=None, \n",
    "                                               multiplicity=1, \n",
    "                                               mapping_method='jordan_wigner',)\n",
    "# 打印结果\n",
    "print(\"\\nThe generated h2 Hamiltonian is \\n\", molecular_hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**注释4：** 生成这个哈密顿量的几何构型中，两个氢原子间的原子间隔（interatomic distance）为 $d = 74$ pm。\n",
    "\n",
    "除了输入分子的几何结构外，我们还支持读取分子的几何构型文件 (`.xyz` 文件)，关于量子化学工具包更多的用法请参考[哈密顿量的构造](./BuildingMolecule_CN.ipynb)教程。如果你需要测试更多分子的几何构型，请移步至这个[数据库](http://smart.sns.it/molecules/index.html)。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 搭建量子神经网络（QNN）和试探波函数\n",
    "\n",
    "在实现VQE的过程中，我们首先需要设计量子神经网络QNN（也可以理解为参数化量子电路）来准备试探波函数 $|\\Psi(\\boldsymbol\\theta)\\rangle$。这里，我们提供一个预设好的的深度为 $D$ 层的 4-量子比特的量子电路模板，图中的虚线框内为一层：\n",
    "\n",
    "![Utheta.jpg](https://release-data.cdn.bcebos.com/PIC%2FUtheta.jpg)\n",
    "\n",
    "- 我们预设一些该参数化电路的参数，比如宽度为 $N = 4$ 量子位。\n",
    "\n",
    "- 初始化其中的变量参数，${\\bf{\\theta }}$ 代表我们量子神经网络中的参数组成的向量。\n",
    "\n",
    "接下来我们根据上图中的电路设计，通过 Paddle Quantum 的 `UAnsatz` 函数和内置的 `real_entangled_layer(theta, D)` 电路模板来高效搭建量子神经网络。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def U_theta(theta, Hamiltonian, N, D):\n",
    "    \"\"\"\n",
    "    Quantum Neural Network\n",
    "    \"\"\"\n",
    "    \n",
    "    # 按照量子比特数量/网络宽度初始化量子神经网络\n",
    "    cir = UAnsatz(N)\n",
    "    \n",
    "    # 内置的 {R_y + CNOT} 电路模板\n",
    "    cir.real_entangled_layer(theta[:D], D)\n",
    "    \n",
    "    # 铺上最后一列 R_y 旋转门\n",
    "    for i in range(N):\n",
    "        cir.ry(theta=theta[D][i][0], which_qubit=i)\n",
    "        \n",
    "    # 量子神经网络作用在默认的初始态 |0000> 上\n",
    "    fin_state = cir.run_state_vector()\n",
    "    \n",
    "    # 计算给定哈密顿量的期望值\n",
    "    expectation_val = cir.expecval(Hamiltonian)\n",
    "\n",
    "    return expectation_val, cir, fin_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 配置训练模型 - 损失函数\n",
    "\n",
    "现在我们已经有了数据和量子神经网络的架构，我们将进一步定义训练参数、模型和损失函数。通过作用量子神经网络 $U(\\theta)$ 在初始态 $|0..0\\rangle$ 上，我们将得到输出态 $\\left| {\\psi \\left( {\\bf{\\theta }} \\right)} \\right\\rangle $。进一步，在VQE模型中的损失函数一般由量子态 $\\left| {\\psi \\left( {\\bf{\\theta }} \\right)} \\right\\rangle$ 关于哈密顿量 $H$ 的期望值 (能量期望值 expectation value) 给出，\n",
    "\n",
    "$$\n",
    "\\min_{\\boldsymbol\\theta}  \\mathcal{L}(\\boldsymbol \\theta) = \\min_{\\boldsymbol\\theta} \\langle \\Psi(\\boldsymbol\\theta)|H |\\Psi(\\boldsymbol\\theta)\\rangle\n",
    "= \\min_{\\boldsymbol\\theta} \\sum_k c_k~\\langle \\Psi(\\boldsymbol\\theta)| \\bigotimes_j \\sigma_j^{(k)}|\\Psi(\\boldsymbol\\theta)\\rangle.\n",
    "\\tag{8}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class StateNet(paddle.nn.Layer):\n",
    "    \"\"\"\n",
    "    Construct the model net\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, shape, dtype=\"float64\"):\n",
    "        super(StateNet, self).__init__()\n",
    "        \n",
    "        # 初始化 theta 参数列表，并用 [0, 2*pi] 的均匀分布来填充初始值\n",
    "        self.theta = self.create_parameter(shape=shape, \n",
    "                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*PI),\n",
    "                                           dtype=dtype, is_bias=False)\n",
    "        \n",
    "    # 定义损失函数和前向传播机制\n",
    "    def forward(self, N, D):\n",
    "        \n",
    "        # 计算损失函数/期望值\n",
    "        loss, cir, fin_state = U_theta(self.theta, molecular_hamiltonian.pauli_str, N, D)\n",
    "\n",
    "        return loss, cir, fin_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 配置训练模型 - 模型参数\n",
    "\n",
    "在进行量子神经网络的训练之前，我们还需要进行一些训练的超参数设置，主要是学习速率（LR, learning rate）、迭代次数（ITR, iteration）和量子神经网络计算模块的深度（D, Depth）。这里我们设定学习速率为 0.5, 迭代次数为 50 次。读者不妨自行调整来直观感受下超参数调整对训练效果的影响。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:14:03.744957Z",
     "start_time": "2021-04-30T09:14:03.738881Z"
    }
   },
   "outputs": [],
   "source": [
    "ITR = 80  # 设置训练的总迭代次数\n",
    "LR = 0.4   # 设置学习速率\n",
    "D = 2      # 设置量子神经网络中重复计算模块的深度 Depth\n",
    "N = molecular_hamiltonian.n_qubits # 设置参与计算的量子比特数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 进行训练\n",
    "\n",
    "当训练模型的各项参数都设置完成后，我们将数据转化为 Paddle 中的张量，进而进行量子神经网络的训练。过程中我们用的是Adam Optimizer，也可以调用Paddle中提供的其他优化器。我们将训练过程中的结果存储在summary_data文件中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 20 loss: -1.0861\n",
      "iter: 20 Ground state energy: -1.0861 Ha\n",
      "iter: 40 loss: -1.1304\n",
      "iter: 40 Ground state energy: -1.1304 Ha\n",
      "iter: 60 loss: -1.1363\n",
      "iter: 60 Ground state energy: -1.1363 Ha\n",
      "iter: 80 loss: -1.1372\n",
      "iter: 80 Ground state energy: -1.1372 Ha\n",
      "\n",
      "训练后的电路：\n",
      "--Ry(1.572)----*--------------x----Ry(1.568)----*--------------x----Ry(0.012)--\n",
      "               |              |                 |              |               \n",
      "--Ry(4.724)----x----*---------|----Ry(4.722)----x----*---------|----Ry(1.619)--\n",
      "                    |         |                      |         |               \n",
      "--Ry(-0.03)---------x----*----|----Ry(4.954)---------x----*----|----Ry(-0.41)--\n",
      "                         |    |                           |    |               \n",
      "--Ry(4.288)--------------x----*----Ry(1.564)--------------x----*----Ry(3.137)--\n",
      "                                                                               \n"
     ]
    }
   ],
   "source": [
    "# 确定网络的参数维度\n",
    "net = StateNet(shape=[D + 1, N, 1])\n",
    "\n",
    "# 一般来说，我们利用Adam优化器来获得相对好的收敛，\n",
    "# 当然你可以改成SGD或者是RMS prop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# 记录优化结果\n",
    "summary_iter, summary_loss = [], []\n",
    "\n",
    "# 优化循环\n",
    "for itr in range(1, ITR + 1):\n",
    "\n",
    "    # 前向传播计算损失函数\n",
    "    loss, cir, fin_state = net(N, D)\n",
    "\n",
    "    # 在动态图机制下，反向传播极小化损失函数\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # 更新优化结果\n",
    "    summary_loss.append(loss.numpy())\n",
    "    summary_iter.append(itr)\n",
    "\n",
    "    # 打印结果\n",
    "    if itr % 20 == 0:\n",
    "        print(\"iter:\", itr, \"loss:\", \"%.4f\" % loss.numpy())\n",
    "        print(\"iter:\", itr, \"Ground state energy:\", \"%.4f Ha\" \n",
    "                                            % loss.numpy())\n",
    "    if itr == ITR:\n",
    "        print(\"\\n训练后的电路：\") \n",
    "        print(cir)\n",
    "\n",
    "# 储存训练结果到 output 文件夹\n",
    "os.makedirs(\"output\", exist_ok=True)\n",
    "savez(\"./output/summary_data\", iter = summary_iter, \n",
    "                               energy=summary_loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 测试效果\n",
    "我们现在已经完成了量子神经网络的训练，通过 VQE 得到的基态能量的估计值大致为 $E_0 \\approx -1.137$ Ha，这与通过 `Psi4` 在 sto-3g 基底下使用 FCI (full configuration-interaction) 方法计算得到的基态能量值 $E_0 = -1.13728$ Ha 是在化学精度 $\\varepsilon = 1.6 \\times 10^{-3}$ Ha 内相符合的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:14:21.341323Z",
     "start_time": "2021-04-30T09:14:20.710152Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGwCAYAAAC5ACFFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAA9hAAAPYQGoP6dpAABgnklEQVR4nO3deVxUVeMG8OcOy7CDArIosrjihoppqKWvUu5LmUsv5ZrWm5aaWZqllq9ppVaaabaI/dJXW9zK0swtF9xQ3MENxQVEQ0FA1jm/P44zMAIj4AwDw/P9fO6H4c65957LwPDMOeeeqwghBIiIiIioWCpzV4CIiIioMmNYIiIiIjKAYYmIiIjIAIYlIiIiIgMYloiIiIgMYFgiIiIiMoBhiYiIiMgAa3NXwBJoNBpcv34dzs7OUBTF3NUhIiKiUhBC4O7du/D19YVKVXL7EcOSEVy/fh1+fn7mrgYRERGVw5UrV1CnTp0Sn2dYMgJnZ2cA8oft4uJi5toQERFRaaSlpcHPz0/3f7wkDEtGoO16c3FxYVgiIiKqYh42hIYDvImIiIgMYFgiIiIiMoBhiYiIiMgAjlkiIrJwGo0GOTk55q4GUYWzsbGBlZXVI++HYYmIyILl5OQgPj4eGo3G3FUhMgs3Nzd4e3s/0jyIDEtERBZKCIHExERYWVnBz8/P4KR7RJZGCIHMzEwkJycDAHx8fMq9L4YlIiILlZeXh8zMTPj6+sLBwcHc1SGqcPb29gCA5ORk1KpVq9xdcvyYQURkofLz8wEAtra2Zq4JkfloPyjk5uaWex8MS0REFo73rKTqzBi//wxLRERERAYwLBEREREZwLBEREREZADDUmV2+zaQmAhwMjkiIiKzYViqzN54AxgzBrh82dw1ISIyi0WLFuH69etl2uaff/5BrVq1cOnSJd06IQQWLFiAwMBAODg4oH///khNTdU9P2TIEMyfP19vPzt37kRAQEC56v0o2xZXf8DwORRXfwD47rvvEBsbW656UAGGpcpMe7kvW5aIqBo6f/483n77bdSoUaNM282ePRv9+vXTCyuTJ0/GkiVLsGLFCuzevRvR0dGYOXOm7vl3330Xs2fP1gtQJenUqRNGjhxZZP2XX34JJyenR54tvbj6P+wcSqr/gQMH8NVXXz1SfYhhqXJjWCIiYxICyMoyzyJEmau7YcMGPPXUU7qJBUsjMzMT3377LUaNGqVbd+DAASxYsABr1qzBk08+idDQUIwePRq///67rkyzZs1Qr149/PDDDw/5EQocPXoUoaGhRZ47fPgwWrZs+UgzpRdX/9KcQ0n179evHzZu3Fju+pDEGbwrM7VafmVYIiJjyM4GBg40z7F/+gmwsyvTJhs2bMCwYcN0369fvx4jRozA7du3ceHCBdSvXx+JiYnw8PCAs7Mz1q1bh/T0dKjVajz++OO67ebNm4euXbuidevWunVeXl64deuW3vH69OmD1atXY+zYsSXW6dy5c7h7926JYelf//pXiduWt/6lPYfi6t+1a1fcuHEDJ0+eRLNmzUqsGxnGlqXKTNuylJ1t3noQEVWwW7duYf/+/ejdu7duXUxMDEJCQgAAx44dg5eXF7y9vREbG4usrCy0bNkSu3fv1gsy2dnZ2LRpE5555hm9/WdlZcHV1VVvXdu2bXHw4EFkG3jPjY6OhpWVla4eWvfu3cPp06f1wsyDylP/spxDcfVXq9V4+umn2br0iNiyVJmxG46IjEmtli085jp2Gfz2229o06YNvLy8dOuOHTumFzaKCx6XL1+Gr6+vbpsjR47g3r17mDRpEt566y3d+tzc3CKtQL6+vsjJyUFSUhL8/f2LrdeRI0eQn59f4r32DIWl8tS/LOdQUv379euHJUuW4J133imxbmQYw1Jlxm44IjImRSlzV5i5/P777+jZs6feupiYGPTp0weAftiIiYlBy5YtAcgWHrtC53j27Fk4OjoiJiZGb1+9evVChw4d9NZpx0ZlZmaWWK8jR47gmWeewfTp0/XWr169GgsXLkSTJk1K3LY89S/LOZRU/549e2LEiBG4desWPDw8Sqwflcwiu+EWL16MgIAA2NnZoV27djh48KDB8j/99BMaN24MOzs7NG/eXG/Qn1mxG46IqqmAgADEx8frvk9LS8OlS5d0424Kh40jR46gVatWAAAPDw/cvn1bbzsPDw/Ur19ft9jY2ODcuXMYMGCA3jFTUlIAAJ6eniXW68iRI+jcuTNatmypt6SkpKBFixYl3tW+vPUvyzmUVP/4+Hi4ubnBzc2txPMiwywuLK1ZswZvvPEGZsyYgSNHjiAkJATdunVDcnJyseX37duH559/HqNGjcLRo0fRv39/9O/fHydPnqzgmheDLUtEVE3169cPmzZt0l2Gn5iYCABwdnZGamoqLl26hJCQECQnJ2PPnj0IDw8HALRq1QqnT5/W7cfDwwOpqakQha7Gmz17Nnr27FmkFejkyZOoU6dOia0vFy9exJ07d4rtajty5Eixg761ylv/spxDSfXfuHEjevbsCWtrdiaVl8WFpQULFmD06NEYMWIEmjRpgqVLl8LBwQHfffddseU///xzdO/eHZMnT0ZwcDBmzZqF1q1b44svvijxGNnZ2UhLS9NbTIJjloiomgoLC4MQAgcOHAAA1K5dG/b29liwYAF27twJGxsb3Lt3D8888wzatWuHLl26AAC6deuGU6dO6VpnunTpgqysLMydOxfx8fH473//i19//RVLliwpcszdu3fj6aefLrFO0dHRUKlUui4zrdzcXJw8edLgeKXy1r8s51BS/Tdu3Ih+/fqVWDd6OIsKSzk5OYiOjtYldABQqVQIDw9HVFRUsdtERUXplQfkL2tJ5QFgzpw5cHV11S1+fn7GOYEHaVuW2A1HRNWMSqVC7969sWHDBgCAk5MTfvzxR2zfvh39+/dHbm4uevTogfbt22PTpk1QFAUA0Lx5c7Ru3Ro//vgjAHl5fWRkJJYsWYKmTZti//792LNnT5H37aysLKxfvx6jR48usU5HjhxBgwYN4OTkpLf+9OnTyM7ONhiWylv/0p5DSfWPj49HXFwcunfvXmLdqBSEBbl27ZoAIPbt26e3fvLkyaJt27bFbmNjYyNWrVqlt27x4sWiVq1aJR4nKytLpKam6pYrV64IACI1NfXRT6KwVauE6N1biMWLjbtfIqoW7t27J06fPi3u3btn7qqUy4YNG0RwcHCR9c8//7x4/vnnhUajKXa73377TQQHB4v8/PxSH+vLL78UTz31lN66HTt2CH9//zLVuTTbVlT9hRDis88+E08//XSp92OJDP0dpKamlur/t0W1LFUUtVoNFxcXvcUk2A1HRNXYU089hcuXL+P8+fN66+Pi4tCuXTtda8yDevXqhTFjxuDatWulPpaNjQ0WLVr0SPUtrYqs/8aNG9G3b99y15Ukixrt5eHhASsrK9y4cUNv/Y0bN+Dt7V3sNt7e3mUqX6F4NRwRVWP29vbIyMjQW5eXl4dTp04VGTf0oAkTJpTpWC+99FIZa1c+FV3/bdu2lWk/VDyLalmytbVFaGio3i+HRqPBtm3bEBYWVuw2YWFhRX6Ztm7dWmL5CsWr4YiI9FhbWyMrKwudOnUy+bECAgLKHFoetm1F1p+Mx6JalgDgjTfewLBhw9CmTRu0bdsWn332GTIyMjBixAgAwNChQ1G7dm3MmTMHADB+/Hh06tQJ8+fPR69evbB69WocPnwYy5YtM+dpSOyGIyIyG1OEJaqaLC4sDR48GDdv3sT06dORlJSEli1bYvPmzbop8xMSEvTuCN2+fXusWrUK7777Lt555x00aNAA69evrxw3HGQ3HBERkdlZXFgCgHHjxmHcuHHFPrdz584i6wYOHIiB5roTtyHshiMiIjI7ixqzZHHYDUdERGR2DEuVGbvhiIiIzI5hqTJjNxwREZHZMSxVZuyGIyIiMjuGpcqM3XBERFXKzJkzHzrhJFU9DEuVmbYbLj9fLkRE1URSUhLGjx+P+vXrw87ODl5eXujQoQOWLFmCzMxMc1fPZIYPH47+/fuXeTuGNNOyyKkDLIa2ZQmQXXH29uarCxFRBbl48SI6dOgANzc3fPjhh2jevDnUajVOnDiBZcuWoXbt2iXe7yw3Nxc2NjYVXGN6VPn5+VAURW8exMqkctaKpAfDEhFRNfDqq6/C2toahw8fxqBBgxAcHIygoCD069cPmzZtQp8+fXRlFUXBkiVL0LdvXzg6OmL27NkAgCVLlqBevXqwtbVFo0aN8H//93+6bS5dugRFURATE6Nbd+fOHSiKopuLb+fOnVAUBdu2bUObNm3g4OCA9u3bIy4uTq+uc+fOhZeXF5ydnTFq1ChkZWU99Px+/vlnNG/eHPb29nB3d0d4eDgyMjIwc+ZMrFixAhs2bICiKHr1efvtt9GwYUM4ODggKCgI7733HnJzcwEAkZGReP/993Hs2DHddpGRkbrzeumll+Dp6QkXFxd06dIFx44dM1i/K1euYNCgQXBzc0PNmjXRr18/XLp0Sfe8tvVr3rx58PHxgbu7O8aOHaurDwBkZ2fjzTffRO3ateHo6Ih27drpzXMYGRkJNzc3bNy4EU2aNIFarUZCQgISExPRq1cv2NvbIzAwEKtWrUJAQAA+++wzAMDIkSPRu3dvvfrm5uaiVq1a+Pbbbx/6sy83QY8sNTVVABCpqanG3/kzzwjRu7cQN24Yf99EZNHu3bsnTp8+Le7du/fAerloNAXrcnPlupycB/dRctns7NKVLYtbt24JRVHEnDlzSlUegKhVq5b47rvvxIULF8Tly5fF2rVrhY2NjVi8eLGIi4sT8+fPF1ZWVmL79u1CCCHi4+MFAHH06FHdfm7fvi0AiB07dgghhNixY4cAINq1ayd27twpTp06JZ544gnRvn173TZr1qwRarVafPPNNyI2NlZMmzZNODs7i5CQkBLre/36dWFtbS0WLFgg4uPjxfHjx8XixYvF3bt3xd27d8WgQYNE9+7dRWJiokhMTBTZ93/Is2bNEnv37hXx8fFi48aNwsvLS3z00UdCCCEyMzPFpEmTRNOmTXXbZWZmCiGECA8PF3369BGHDh0SZ8+eFZMmTRLu7u7in3/+KbZ+OTk5Ijg4WIwcOVIcP35cnD59Wvz73/8WjRo10tVl2LBhwsXFRbzyyivizJkz4tdffxUODg5i2bJluv289NJLon379uLvv/8W58+fF5988olQq9Xi7NmzQgghli9fLmxsbET79u3F3r17RWxsrMjIyBDh4eGiZcuWYv/+/SI6Olp06tRJ2Nvbi08//VQIIcTevXuFlZWVuH79uu5Ya9euFY6OjuLu3bvFnlNJfwdClP7/N8OSEZg0LA0eLMPSlSvG3zcRWbSS/kn07i2XO3cK1q1ZI9ctXKi/jwEDin5eW79ervvkE/2y//63XH/5csG6zZvLVuf9+/cLAGLt2rV6693d3YWjo6NwdHQUb731lm49ADFhwgS9su3btxejR4/WWzdw4EDRs2dPIUTZwtJff/2lK7Np0yYBQPfzDAsLE6+++qrecdq1a2cwLEVHRwsA4tKlS8U+P2zYMNGvX78St9f65JNPRGhoqO77GTNmFDnu7t27hYuLi8jKytJbX69ePfHVV18Vu9//+7//E40aNRKaQok3Oztb2Nvbiy1btujq6O/vL/Ly8nRlBg4cKAYPHiyEEOLy5cvCyspKXLt2TW/fXbt2FVOnThVCyLAEQMTExOieP3PmjAAgDh06pFt37tw5AUAXloQQokmTJrqgKIQQffr0EcOHDy/2fIQwTlhiN1xlx7mWiIhw8OBBxMTEoGnTpsh+4ArhNm3a6H1/5swZdOjQQW9dhw4dcObMmTIft0WLFrrHPj4+AIDk5GTdcdq1a6dXPiwsTPd49+7dcHJy0i0rV65ESEgIunbtiubNm2PgwIH4+uuvcfv27YfWY82aNejQoQO8vb3h5OSEd999FwkJCQa3OXbsGNLT0+Hu7q5Xj/j4eFy4cKHEbc6fPw9nZ2dd+Zo1ayIrK0tvm6ZNm8LKykrvZ6P9uZw4cQL5+flo2LCh3nF37dqltw9bW1u9n29cXBysra3RunVr3br69eujRo0aenV86aWXsHz5cgDAjRs38Mcff2DkyJEP+xE+Eg7wruw4fQARGdlPP8mv2s9iAPDss0DfvkCh/38AgB9+KFq2Vy+gWzfgwbG42iEjhct27Vq2utWvXx+KohQZGxQUFAQAsC/mQhdHR8cyHUM7iFgIoVtXeLxNYYUHiyuKAgDQaDSlOk6bNm30xkV5eXnBysoKW7duxb59+/Dnn39i0aJFmDZtGg4cOIDAwMBi9xMVFYWIiAi8//776NatG1xdXbF69WrMnz/f4PHT09Ph4+NT7D1R3dzcStwmNDQUK1euLPKcp6en7vGDg+gVRdH9XNLT02FlZYXo6Gi9QAUATk5Ousf29va6n2lZDB06FFOmTEFUVBT27duHwMBAPPHEE2XeT1kwLFV2nJiSiIzMzq7oOmtruZiibFm4u7vjqaeewhdffIHXXnutzEEIAIKDg7F3714MGzZMt27v3r1o0qQJgIJ/+omJiWjVqhUA6IWashznwIEDGDp0qG7d/v37dY/t7e1Rv379ItspioIOHTqgQ4cOmD59Ovz9/bFu3Tq88cYbsLW1Rf4DU8Xs27cP/v7+mDZtmm7d5cuX9coUt13r1q2RlJQEa2trBAQElOqcWrdujTVr1qBWrVpwcXEp1TYPatWqFfLz85GcnFymENOoUSPk5eXh6NGjCA0NBQCcP3++SMubu7s7+vfvj+XLlyMqKgojRowoVz3Lgt1wlR274Yiomvnyyy+Rl5eHNm3aYM2aNThz5gzi4uLwww8/IDY2tkhrxYMmT56MyMhILFmyBOfOncOCBQuwdu1avPnmmwBkiHn88ccxd+5cnDlzBrt27cK7775b5nqOHz8e3333HZYvX46zZ89ixowZOHXqlMFtDhw4gA8//BCHDx9GQkIC1q5di5s3byI4OBgAEBAQgOPHjyMuLg63bt1Cbm4uGjRogISEBKxevRoXLlzAwoULsW7dOr39BgQEID4+HjExMbh16xays7MRHh6OsLAw9O/fH3/++ScuXbqEffv2Ydq0aTh8+HCx9YuIiICHhwf69euH3bt3Iz4+Hjt37sTrr7+Oq1evlurn0rBhQ0RERGDo0KFYu3Yt4uPjcfDgQcyZMwebNm0qcbvGjRsjPDwcY8aMwcGDB3H06FGMGTOm2Baol156CStWrMCZM2f0QrHJGBzRRKVi0gHeb78tR0zu3m38fRORRTM0sLWyu379uhg3bpwIDAwUNjY2wsnJSbRt21Z88sknIiMjQ1cOgFi3bl2R7b/88ksRFBQkbGxsRMOGDcX333+v9/zp06dFWFiYsLe3Fy1bthR//vlnsQO8b9++rdvm6NGjAoCIj4/XrZs9e7bw8PAQTk5OYtiwYeKtt94yOMD79OnTolu3bsLT01Oo1WrRsGFDsWjRIt3zycnJ4qmnnhJOTk569Zk8ebJwd3cXTk5OYvDgweLTTz8Vrq6uuu2ysrLEgAEDhJubmwAgli9fLoQQIi0tTbz22mvC19dX2NjYCD8/PxERESESEhJKrGNiYqIYOnSo8PDwEGq1WgQFBYnRo0fr/scVNwh9/PjxolOnTrrvc3JyxPTp00VAQICwsbERPj4+4plnnhHHjx8XQsgB3oXrr3X9+nXRo0cPoVarhb+/v1i1apWoVauWWLp0qV45jUYj/P39dYP2DTHGAG9FiEKdtlQuaWlpcHV1RWpqarmbLUs0fTpw9CgwcSLQpYtx901EFi0rKwvx8fEIDAyEXXF9ZESV3NWrV+Hn54e//voLXQsNgEtPT0ft2rWxfPlyPPvsswb3YejvoLT/vzlmqbLjmCUiIqomtm/fjvT0dDRv3hyJiYl46623EBAQgCeffBKAHFx/69YtzJ8/H25ubiXO5G5sDEuVnXbMEq+GIyIiC5ebm4t33nkHFy9ehLOzM9q3b4+VK1fqrr5LSEhAYGAg6tSpg8jISFiX9QqCcmJYquzYskRERNVEt27d0K1btxKfDwgIgDlGD/FquMqOYYmIiMisGJYqO3bDEdEj4nU8VJ0Z4/efYamyY8sSEZWTdj6iHL5/UDWWmZkJoOis42XBMUuVHcMSEZWTtbU1HBwccPPmTdjY2Ohu80FUHQghkJmZieTkZLi5uT10MlNDGJYqO3bDEVE5KYoCHx8fxMfHF7k9BlF14ebmBm9v70faB8NSZceWJSJ6BLa2tmjQoAG74qhasrGxeaQWJS2GpcqO94YjokekUqk4gzfRI2AHdmXHliUiIiKzYliq7DhmiYiIyKwYlio7tiwRERGZFcNSZcewREREZFYMS5Udu+GIiIjMimGpsmPLEhERkVkxLFV2DEtERERmxbBU2RWeZ4k3wyQiIqpwDEuVnbZlCQByc81XDyIiomqKYamyKxyWOMibiIiowjEsVXZWVnIBOG6JiIjIDBiWqgIO8iYiIjIbhqWqQBuW2A1HRERU4RiWqoLCV8QRERFRhWJYqgrYskRERGQ2DEtVAccsERERmQ3DUlXA+8MRERGZDcNSVcCWJSIiIrNhWKoKGJaIiIjMhmGpKmA3HBERkdlYVFhKSUlBREQEXFxc4ObmhlGjRiE9Pd1g+ddeew2NGjWCvb096tati9dffx2pqakVWOtSYMsSERGR2VhUWIqIiMCpU6ewdetW/Pbbb/j7778xZsyYEstfv34d169fx7x583Dy5ElERkZi8+bNGDVqVAXWuhQYloiIiMxGEUIIc1fCGM6cOYMmTZrg0KFDaNOmDQBg8+bN6NmzJ65evQpfX99S7eenn37CCy+8gIyMDFhbW5dqm7S0NLi6uiI1NRUuLi7lPocSffMNsGED8NxzwLBhxt8/ERFRNVTa/98W07IUFRUFNzc3XVACgPDwcKhUKhw4cKDU+9H+wAwFpezsbKSlpektJsVJKYmIiMzGYsJSUlISatWqpbfO2toaNWvWRFJSUqn2cevWLcyaNctg1x0AzJkzB66urrrFz8+v3PUuFXbDERERmU2lD0tTpkyBoigGl9jY2Ec+TlpaGnr16oUmTZpg5syZBstOnToVqampuuXKlSuPfHyDeG84IiIisyndoBwzmjRpEoYPH26wTFBQELy9vZGcnKy3Pi8vDykpKfD29ja4/d27d9G9e3c4Oztj3bp1sLGxMVherVZDrQ0wFYHdcERERGZT6cOSp6cnPD09H1ouLCwMd+7cQXR0NEJDQwEA27dvh0ajQbt27UrcLi0tDd26dYNarcbGjRthZ2dntLobDbvhiIiIzKbSd8OVVnBwMLp3747Ro0fj4MGD2Lt3L8aNG4chQ4boroS7du0aGjdujIMHDwKQQenpp59GRkYGvv32W6SlpSEpKQlJSUnIz8835+noYzccERGR2VT6lqWyWLlyJcaNG4euXbtCpVJhwIABWLhwoe753NxcxMXFITMzEwBw5MgR3ZVy9evX19tXfHw8AgICKqzuBrEbjoiIyGwsKizVrFkTq1atKvH5gIAAFJ5WqnPnzqgS00yxG46IiMhsLKYbzqIxLBEREZkNw1JVwBvpEhERmQ3DUlXAliUiIiKzYViqCng1HBERkdkwLFUFha+GqwoD0omIiCwIw1JVoG1ZEgKoTPM/ERERVQMMS1WBtmUJ4CBvIiKiCsawVBVYWwOKIh9z3BIREVGFYliqChSFV8QRERGZCcNSVcGwREREZBYMS1UFJ6YkIiIyC4alqoItS0RERGbBsFRVMCwRERGZBcNSVcFuOCIiIrNgWKoq2LJERERkFgxLVQXDEhERkVkwLFUV7IYjIiIyC4alqoItS0RERGbBsFRVMCwRERGZBcNSVcFuOCIiIrNgWKoqtC1LDEtEREQVimGpqmA3HBERkVkwLFUVDEtERERmwbBUVXDMEhERkVkwLFUVbFkiIiIyC4alqoJhiYiIyCwYlqoKdsMRERGZBcNSVcGWJSIiIrNgWKoqGJaIiIjMgmGpqmA3HBERkVkwLFUVbFkiIiIyC4alqoJhiYiIyCwYlqoKbTccwxIREVGFYliqKrQtS3l5gEZj3roQERFVIwxLVYU2LAFsXSIiIqpADEtVhbYbDuAVcURERBWIYamqUBTA2lo+ZssSERFRhWFYqko4yJuIiKjCMSxVJZyYkoiIqMIxLFUlnGuJiIiowjEsVSUMS0RERBWOYakqYTccERFRhWNYqkrYskRERFThGJaqEoYlIiKiCsewVJVowxK74YiIiCoMw1JVwnmWiIiIKpxFhaWUlBRERETAxcUFbm5uGDVqFNLT00u1rRACPXr0gKIoWL9+vWkrWl7shiMiIqpwFhWWIiIicOrUKWzduhW//fYb/v77b4wZM6ZU23722WdQFMXENXxE7IYjIiKqcNbmroCxnDlzBps3b8ahQ4fQpk0bAMCiRYvQs2dPzJs3D76+viVuGxMTg/nz5+Pw4cPw8fGpqCqXHbvhiIiIKpzFtCxFRUXBzc1NF5QAIDw8HCqVCgcOHChxu8zMTPz73//G4sWL4e3tXapjZWdnIy0tTW+pEOyGIyIiqnAWE5aSkpJQq1YtvXXW1taoWbMmkpKSStxu4sSJaN++Pfr161fqY82ZMweurq66xc/Pr9z1LhN2wxEREVW4Sh+WpkyZAkVRDC6xsbHl2vfGjRuxfft2fPbZZ2XaburUqUhNTdUtV65cKdfxy4zdcERERBWu0o9ZmjRpEoYPH26wTFBQELy9vZGcnKy3Pi8vDykpKSV2r23fvh0XLlyAm5ub3voBAwbgiSeewM6dO4vdTq1WQ60NLhWJ3XBEREQVrtKHJU9PT3h6ej60XFhYGO7cuYPo6GiEhoYCkGFIo9GgXbt2xW4zZcoUvPTSS3rrmjdvjk8//RR9+vR59MobG8MSERFRhav0Yam0goOD0b17d4wePRpLly5Fbm4uxo0bhyFDhuiuhLt27Rq6du2K77//Hm3btoW3t3exrU5169ZFYGBgRZ/Cw/FGukRERBXO6GFJo9Fg165d2L17Ny5fvozMzEx4enqiVatWCA8PN+lg6JUrV2LcuHHo2rUrVCoVBgwYgIULF+qez83NRVxcHDIzM01WB5NiyxIREVGFM1pYunfvHubPn48lS5YgJSUFLVu2hK+vL+zt7XH+/HmsX78eo0ePxtNPP43p06fj8ccfN9ahdWrWrIlVq1aV+HxAQACEEAb38bDnzYpXwxEREVU4o4Wlhg0bIiwsDF9//TWeeuop2NjYFClz+fJlrFq1CkOGDMG0adMwevRoYx2+euDVcERERBXOaGHpzz//RHBwsMEy/v7+mDp1Kt58800kJCQY69DVB1uWiIiIKpzR5ll6WFAqzMbGBvXq1TPWoauPh41ZSk4GcnMrrj5ERETVgEmvhsvMzERCQgJyHvjn3qJFC1Me1nIZ6obbvx+YPRtwdQV69JBLzZoVWz8iIiILZJKwdPPmTYwYMQJ//PFHsc/n5+eb4rCWr3DLkhCAohQ8t3+//JqaCqxeDfz8M9CxI9CnD9CwYcXXlYiIyEKY5HYnEyZMwJ07d3DgwAHY29tj8+bNWLFiBRo0aICNGzea4pDVgzYsAUW7286elV/79QOCg4G8PGDnTmDSJOCtt2SIIiIiojIzScvS9u3bsWHDBrRp0wYqlQr+/v546qmn4OLigjlz5qBXr16mOKzlKxyWsrMLvs/MBK5elY+few5wcwPOnwd+/RX4+2/gzBkgKgro3r3Cq0xERFTVmaRlKSMjA7Vq1QIA1KhRAzdv3gQgbyVy5MgRUxyyerC2BlT3X7LC45bOn5fdcrVqyaAEAPXrAxMnAv37y+/j4yuypkRERBbDJGGpUaNGiIuLAwCEhITgq6++wrVr17B06VL4+PiY4pDVR3FXxGm74Iobm6S9bcvFi6atFxERkYUySTfc+PHjkZiYCACYMWMGunfvjpUrV8LW1haRkZGmOGT1oVYDWVn6cy2dOye/GgpLly4VHRRORERED2WSsPTCCy/oHoeGhuLy5cuIjY1F3bp14eHhYYpDVh/FtSzdb8UrNiz5+sptsrKApCSALXtERERlYpJuuAc5ODigdevWDErG8GBY+ucfuSgKUNxEn1ZWgL+/fMyuOCIiojIzasvSG2+8UapyCxYsMOZhqxftxJTabjhtF5y/P2BnV/w2gYGyXHw80KGD6etIRERkQYwalo4ePar3/Z49exAaGgp7e3vdOoVjZh7Ngy1LhgZ3a3GQNxERUbkZNSzt2LFD73tnZ2esWrUKQUFBxjxM9VaesKT9+XP6ACIiojKrkDFLZESF7w8nhOEr4bQCAuTXW7eAu3dNWj0iIiJLw7BU1RQes3T1qpy9W60G6tYteRsHB8DbWz5m6xIREVGZMCxVNYW74bRdcPXry6veDNGOW2JYIiIiKhOjjlk6fvy43vdCCMTGxiI9PV1vfYsWLYx52OqluLBkqAtOKyhI3h+OYYmIiKhMjBqWWrZsCUVRIITQrevduzcA6NYrioL8/HxjHrZ6KdwNV5rxSlq8Io6IiKhcjBqW4tlqYXralqWMjIJWorKEpStXgLw8eVNeIiIieiij/sf0184UTaajDUuxsTL0uLoCnp4P387TE3B0lCHrypWC8EREREQGGW2Ad0JCQpnKX7t2zViHrl603XCXLsmvDRuW7ua4isJB3kREROVgtLD02GOP4eWXX8ahQ4dKLJOamoqvv/4azZo1wy+//GKsQ1cv2pYl7biw0nTBaXHcEhERUZkZrRvu9OnTmD17Np566inY2dkhNDQUvr6+sLOzw+3bt3H69GmcOnUKrVu3xscff4yePXsa69DVizYsaTVoUPpt2bJERERUZkZrWXJ3d8eCBQuQmJiIL774Ag0aNMCtW7dw7v4VWxEREYiOjkZUVBSD0qPQdsNplaVlqfBtTwpdsUhEREQlM/olUfb29njuuefw3HPPGXvXBOi3LPn4AM7Opd/Wzw9QqeQtT/75B/DwMH79iIiILAxn8K5qCoelsrQqabf185OP2RVHRERUKgxLVc2jhCWA45aIiIjKiGGpqik8ZolhiYiIyOQYlqoaOzv51cqqYMB2WWi34fQBREREpWKSe15kZGTA0dHRFLum2rWB9u1lC9GD0wiUhrZlKTERyMoqCF9ERERULJO0LHl5eWHkyJHYs2ePKXZfvalUwNSpwJAh5dve1RWoWVNOHaCdBZyIiIhKZJKw9MMPPyAlJQVdunRBw4YNMXfuXFy/ft0Uh6Ly4LglIiKiUjNJWOrfvz/Wr1+Pa9eu4ZVXXsGqVavg7++P3r17Y+3atcjLyzPFYam0eNsTIiKiUjPpAG9PT0+88cYbOH78OBYsWIC//voLzz33HHx9fTF9+nRkZmaa8vBUErYsERERlZpJBnhr3bhxAytWrEBkZCQuX76M5557DqNGjcLVq1fx0UcfYf/+/fjzzz9NWQUqjvaKuEuXAI1GjoMiIiKiYpkkLK1duxbLly/Hli1b0KRJE7z66qt44YUX4ObmpivTvn17BAcHm+Lw9DC+voCNDZCdDdy8CXh5mbtGRERElZZJwtKIESMwZMgQ7N27F4899lixZXx9fTFt2jRTHJ4eRqWSgenyZeDqVYYlIiIiA0wSlhITE+Hg4GCwjL29PWbMmGGKw1Np1K5dEJZCQ81dGyIiokrLJGEpLy8PaWlpRdYrigK1Wg3b8kymSMZVp478evWqeetBRERUyZkkLLm5uUFRlBKfr1OnDoYPH44ZM2ZAxcHF5sGwREREVComCUuRkZGYNm0ahg8fjrZt2wIADh48iBUrVuDdd9/FzZs3MW/ePKjVarzzzjumqAI9DMMSERFRqZgkLK1YsQLz58/HoEGDdOv69OmD5s2b46uvvsK2bdtQt25dzJ49m2HJXGrXll/v3AEyMgDey4+IiKhYJukD27dvH1q1alVkfatWrRAVFQUA6NixIxISEkxxeCoNBwd5jziArUtEREQGmCQs+fn54dtvvy2y/ttvv4Wfnx8A4J9//kGNGjVMcXgqLXbFERERPZRJuuHmzZuHgQMH4o8//tDNs3T48GHExsbi559/BgAcOnQIgwcPNsXhqbTq1AGOHweuXTN3TYiIiCotk7Qs9e3bF3FxcejZsydSUlKQkpKCHj16IDY2Fr179wYA/Oc//8GCBQuMetyUlBRERETAxcUFbm5uGDVqFNLT0x+6XVRUFLp06QJHR0e4uLjgySefxL1794xat0qJLUtEREQPZfSWpdzcXHTv3h1Lly7FnDlzjL17gyIiIpCYmIitW7ciNzcXI0aMwJgxY7Bq1aoSt4mKikL37t0xdepULFq0CNbW1jh27Fj1mNKAYYmIiOihFCGEMPZOPT09sW/fPjRo0MDYuy7RmTNn0KRJExw6dAht2rQBAGzevBk9e/bE1atX4evrW+x2jz/+OJ566inMmjWr3MdOS0uDq6srUlNT4eLiUu79VLibN4GRIwFra+DnnwErK3PXiIiIqMKU9v+3SZpPXnjhhWIHeJtSVFQU3NzcdEEJAMLDw6FSqXDgwIFit0lOTsaBAwdQq1YttG/fHl5eXujUqRP27Nlj8FjZ2dlIS0vTW6okDw/A1hbIywNu3DB3bYiIiColk93u5LvvvsNff/2F0NBQOD4wh4+xxyoBQFJSEmrVqqW3ztraGjVr1kRSUlKx21y8eBEAMHPmTMybNw8tW7bE999/j65du+LkyZMltozNmTMH77//vnFPwBwURXbFXbwou+JKaH0jIiKqzkzSsnTy5Em0bt0azs7OOHv2LI4ePapbYmJiyrSvKVOmQFEUg0tsbGy56qnRaAAAL7/8MkaMGIFWrVrh008/RaNGjfDdd9+VuN3UqVORmpqqW65cuVKu41cKHLdERERkkElalnbs2GG0fU2aNAnDhw83WCYoKAje3t5ITk7WW5+Xl4eUlBR4e3sXu52Pjw8AoEmTJnrrg4ODDU6YqVaroVarS1H7KkA7kzfDEhERUbFMEpa0zp8/jwsXLuDJJ5+Evb09hBAGb7BbHE9PT3h6ej60XFhYGO7cuYPo6GiEhoYCALZv3w6NRoN27doVu01AQAB8fX0RFxent/7s2bPo0aNHmepZZbFliYiIyCCTdMP9888/6Nq1Kxo2bIiePXsiMTERADBq1ChMmjTJFIdEcHAwunfvjtGjR+PgwYPYu3cvxo0bhyFDhuiuhLt27RoaN26MgwcPAgAURcHkyZOxcOFC/Pzzzzh//jzee+89xMbGYtSoUSapZ6XDsERERGSQScLSxIkTYWNjg4SEBDg4OOjWDx48GJs3bzbFIQEAK1euROPGjdG1a1f07NkTHTt2xLJly3TP5+bmIi4uDpmZmbp1EyZMwNSpUzFx4kSEhIRg27Zt2Lp1K+rVq2eyelYq2kHdd+8CVfWqPiIiIhMyyTxL3t7e2LJlC0JCQuDs7Ixjx44hKCgIFy9eRIsWLUo1q3ZVUmXnWdIaOVLOufTxx0BwsLlrQ0REVCHMOs9SRkaGXouSVkpKiuUMjLYk7IojIiIqkUnC0hNPPIHvv/9e972iKNBoNPj444/xr3/9yxSHpEehDUtVeQoEIiIiEzHJ1XAff/wxunbtisOHDyMnJwdvvfUWTp06hZSUFOzdu9cUh6RHoQ1L166Ztx5ERESVkElalpo1a4azZ8+iY8eO6NevHzIyMvDss8/i6NGj1WfgdFXCbjgiIqISmWyeJVdXV0ybNs1Uuydj0oalpCR5nzhrk06/RUREVKWY7L/inTt3cPDgQSQnJ+tuK6I1dOhQUx2WyqNGDcDODsjKAhITAT8/c9eIiIio0jBJWPr1118RERGB9PR0uLi46M3arSgKw1JloygyIJ07J7viGJaIiIh0TDJmadKkSRg5ciTS09Nx584d3L59W7ekpKSY4pD0qDhuiYiIqFgmCUvXrl3D66+/XuxcS1RJMSwREREVyyRhqVu3bjh8+LApdk2mUru2/MrpA4iIiPSYZMxSr169MHnyZJw+fRrNmzeHjY2N3vN9+/Y1xWHpURRuWRJCjmMiIiIi09wbTqUqucFKURTk5+cb+5BmVeXvDQcAOTnAc8/JoPT99/IKOSIiIgtm1nvDaTSaEhdLC0oWw9YW8PKSj9kVR0REpGOSsERVFAd5ExERFWHUsNSzZ0+kpqbqvp87dy7u3Lmj+/6ff/5BkyZNjHlIMiaGJSIioiKMGpa2bNmC7Oxs3fcffvih3rxKeXl5iIuLM+YhyZgYloiIiIowalh6cKy4CcaOkyn5+MivSUnmrQcREVElwjFLVEA7wPvGDXlVHBERERk3LCmKoncfOO06qiI8PACVCsjLA3hbGiIiIgBGnpRSCIHhw4dDrVYDALKysvDKK6/A0dERAPTGM1ElZGUlA1Nysmxdcnc3d42IiIjMzqhhadiwYXrfv/DCC0XKDB061JiHJGPz9i4IS7xykYiIyLhhafny5cbcHZlD4XFLRERExAHe9ABtWOIVcURERAAYluhBbFkiIiLSw7BE+hiWiIiI9DAskT5tWLp1S04hQEREVM0xLJG+GjUAW1s5KeXNm+auDRERkdkxLJE+RQFq1ZKP2RVHRETEsETF4LglIiIiHYYlKophiYiISIdhiYry9pZfGZaIiIgYlqgYbFkiIiLSYViiojiLNxERkQ7DEhWlDUupqUBWlnnrQkREZGYMS1SUkxPg6CgfJyebty5ERERmxrBExeO4JSIiIgAMS1QSjlsiIiICwLBEJWHLEhEREQCGJSoJwxIREREAhiUqCcMSERERAIYlKknhWbyFMG9diIiIzIhhiYpXq5b8mpkJpKebty5ERERmxLBExVOrATc3+ZhdcUREVI0xLFHJOG6JiIiIYYkMYFgiIiJiWCIDCg/yJiIiqqYsKiylpKQgIiICLi4ucHNzw6hRo5D+kMHJSUlJePHFF+Ht7Q1HR0e0bt0av/zySwXVuJJjyxIREZFlhaWIiAicOnUKW7duxW+//Ya///4bY8aMMbjN0KFDERcXh40bN+LEiRN49tlnMWjQIBw9erSCal2J8ZYnRERElhOWzpw5g82bN+Obb75Bu3bt0LFjRyxatAirV6/G9evXS9xu3759eO2119C2bVsEBQXh3XffhZubG6Kjoyuw9pWUNiwlJ3OuJSIiqrYsJixFRUXBzc0Nbdq00a0LDw+HSqXCgQMHStyuffv2WLNmDVJSUqDRaLB69WpkZWWhc+fOJW6TnZ2NtLQ0vcUieXgAigLk5gK3b5u7NkRERGZhMWEpKSkJtbQTKd5nbW2NmjVrIslAN9KPP/6I3NxcuLu7Q61W4+WXX8a6detQv379EreZM2cOXF1ddYufn5/RzqNSsbYGPD3lY45bIiKiaqrSh6UpU6ZAURSDS2xsbLn3/9577+HOnTv466+/cPjwYbzxxhsYNGgQTpw4UeI2U6dORWpqqm65cuVKuY9f6XHcEhERVXPW5q7Aw0yaNAnDhw83WCYoKAje3t5ITk7WW5+Xl4eUlBR4ay+Bf8CFCxfwxRdf4OTJk2jatCkAICQkBLt378bixYuxdOnSYrdTq9VQq9VlP5mqyMsLOHGCLUtERFRtVfqw5OnpCU9tV5ABYWFhuHPnDqKjoxEaGgoA2L59OzQaDdq1a1fsNpmZmQAAlUq/gc3KygoajeYRa24hOH0AERFVc5W+G660goOD0b17d4wePRoHDx7E3r17MW7cOAwZMgS+vr4AgGvXrqFx48Y4ePAgAKBx48aoX78+Xn75ZRw8eBAXLlzA/PnzsXXrVvTv39+MZ1OJMCwREVE1ZzFhCQBWrlyJxo0bo2vXrujZsyc6duyIZcuW6Z7Pzc1FXFycrkXJxsYGv//+Ozw9PdGnTx+0aNEC33//PVasWIGePXua6zQqF87iTURE1ZwiBCfQeVRpaWlwdXVFamoqXFxczF0d40pJAYYNk1MIrF0rr5AjIiKyAKX9/21RLUtkAjVqADY2clLKW7fMXRsiIqIKx7BEhikKoJ2/il1xRERUDTEs0cNx3BIREVVjDEv0cPevJkRCgnnrQUREZAYMS/Rw9erJrxcumLceREREZsCwRA9XOCzx4kkiIqpmGJbo4fz8AFtb4N494Pp1c9eGiIioQjEs0cNZWQGBgfIxu+KIiKiaYVii0uG4JSIiqqYYlqh06teXX8+fN289iIiIKhjDEpUOB3kTEVE1xbBEpVO3rrwvXEYGJ6ckIqJqhWGJSsfaGggIkI85bomIiKoRhiUqPY5bIiKiaohhiUqPV8QREVE1xLBEpacNS+fPc5A3ERFVGwxLVHr+/nKCyrt3gVu3zF0bIiKiCsGwRKVnaysDE8BxS0REVG0wLFHZcNwSERFVMwxLVDa8Io6IiKoZhiUqGw7yJiKiaoZhicomMBBQqYDUVOD2bXPXhoiIyOQYlqhsbG0BPz/5mF1xRERUDTAsUdlpxy1xkDcREVUDDEtUdoXHLREREVk4hiUqO04fQERE1QjDEpVdUBCgKMA//wB37pi7NkRERCbFsERlZ2cH1K4tH7N1iYiILBzDEpVPWSan1GiAa9eA48eBvDzT1ouIiMjIrM1dAaqi6tUDdu4s2rKUmQnExgLx8cDly3K5ehXIyZHP9+oFvPJKhVeXiIiovBiWqHy0LUvnzgFRUcCpU8DJk8DFi8XP7G1rKwPTli3AgAGAp2fF1peIiKicGJaofAID5ddbt4APP9R/ztsbaNAACAgA/P2BunXlunfflV1xP/8M/Oc/FV5lIiKi8mBYovJxdARatgRiYmQgatq0YHF3L36b55+XYenPP4GBAwEPj4qsMRERUbkwLFH5ffABkJ0tr44rjWbN5HLyJPDLL8DLL5u2fkREREbAq+Go/BSl9EFJ6/nn5dctW4CUFOPXiYiIyMgYlqhiNW8OBAcDubly7BIREVElx7BEFUtRgH//Wz5m6xIREVUBDEtU8UJCgMaN5VQCa9eauzZEREQGMSxRxVMUYMgQ+fiPP3h/OSIiqtQYlsg8WrcGGjZk6xIREVV6DEtkHopScGXc778DqanmrQ8REVEJGJbIfEJD5W1TsrOBzZvNXRsiIqJiMSyR+SgK0LevfPzHH0B+vnnrQ0REVAyGJTKvjh0BV1fgn3+AAwfMXRsiIqIiGJbIvGxsgG7d5ONNm8xbFyIiomIwLJH5de8uu+SOHwcSEsxdGyIiIj0MS2R+np7A44/Lx2xdIiKiSsaiwtLs2bPRvn17ODg4wM3NrVTbCCEwffp0+Pj4wN7eHuHh4Th37pxpK0pF9eolv27fDmRkmLcuREREhVhUWMrJycHAgQPxn//8p9TbfPzxx1i4cCGWLl2KAwcOwNHREd26dUNWVpYJa0pFtGgB1KkDZGXJwERERFRJWFRYev/99zFx4kQ0b968VOWFEPjss8/w7rvvol+/fmjRogW+//57XL9+HevXry9xu+zsbKSlpekt9IgUpaB1adMmQAjz1oeIiOg+iwpLZRUfH4+kpCSEh4fr1rm6uqJdu3aIiooqcbs5c+bA1dVVt/j5+VVEdS1fly6AnR1w7Rpw7Ji5a0NERASgmoelpKQkAICXl5feei8vL91zxZk6dSpSU1N1y5UrV0xaz2rDwUEGJoADvYmIqNKo9GFpypQpUBTF4BIbG1uhdVKr1XBxcdFbyEh695ZfDxwAbt4s/34yMoD164FffgE0GqNUjYiIqidrc1fgYSZNmoThw4cbLBMUFFSufXt7ewMAbty4AR8fH936GzduoGXLluXaJz0iPz852Pv4cXkLlKFDy7Z9cjKwcSOwZYscLA4A588DkyYB1pX+152IiCqhSv/fw9PTE56enibZd2BgILy9vbFt2zZdOEpLS8OBAwfKdEUdGVmvXjIs/f478K9/yQD1MBcuAGvXAnv2FLQk+fkBiYlyXU4O8PbbgK2taetOREQWp9J3w5VFQkICYmJikJCQgPz8fMTExCAmJgbp6em6Mo0bN8a6desAAIqiYMKECfjvf/+LjRs34sSJExg6dCh8fX3Rv39/M50FoV07oEED2ZX27rsy8JQkLw9YuBCYMAH4+28ZlEJCgPffBxYvltvb2gIHDwKzZgHZ2RV2GkREZBkUISznGu3hw4djxYoVRdbv2LEDnTt3BiAD0vLly3Vde0IIzJgxA8uWLcOdO3fQsWNHfPnll2jYsGGpj5uWlgZXV1ekpqZy/JKxpKUBU6fK25/UqgXMnStn+n6wzIcfAqdOyakHOnUCnnkGeLBb9sQJ4IMPZLdc06bA9OlyMDkREVVrpf3/bVFhyVwYlkzk9m1gyhTg+nXAx0cGppo15XNXrsgAlJQkg89bbwGhoSXvKzYWmDlTtlY1aCBbnpydK+Q0iIiocirt/2+L6oYjC1OjBjB7tmxZSkyUXWppacCRI8Cbb8qg5O0NzJtnOCgBQOPGcl/OzsC5c3L8UnJyxZwHERFVaWxZMgK2LJlYUpIMNykpgJeXDDlCyC61d94ByvIzv3xZdsOlpACursB77wGNGpmu7kREVGmxZYksh7e3bBVydQVu3JBBqWtXOWC7rOHU3x+YP1+Oa0pNleOidu82Tb2JiMgisGXJCNiyVEEuXQK+/RZ47DGgTx85qLu8srKATz6RV8kBwAsvAIMGPdo+iYioSuEA7wrEsFRFaTTAd98BGzbI7//1L+DVV+X96YiIyOKxG47oYVQq4KWXgLFj5eMdO4Dhw2XrlYF7AxIRUfXCliUjYMuSBYiJkZNYakOSogBt2sh71bVqxe45IiILxG64CsSwZCGEAA4fBn77TU5PoOXnB0ycKOdnIiIii8GwVIEYlizQtWvApk3AX38B9+7Jm/COGPHoA8uJiKjS4JglokdRuzYwZowcAN6+vbwH3ddfy9urFLrXIBERWT6GJSJDnJzkLVdeflm2Lu3fD4wfD8TFmbtmRERUQRiWiB5GUeRA708+kfeoS06WM4pv2CDHORERkUVjWCIqrfr1gU8/BTp2BPLzgW++ARYtkl10RERksRiWiMrC0RF46y1g9GjZ4rR1q7y/3N275q4ZERGZCMMSUVkpCtC3LzBjBmBvD5w8CUyaJK+gIyIii8OwRFReoaFyHFOtWkBiogxMx46Zu1ZERGRkDEtEj8LfH1iwAAgOBjIygOnTgdWrgexsc9eMiIiMhGGJ6FG5ugL//S/QubO8Oe/KlXKOpr/+kt8TEVGVxhm8jYAzeBMAOY3A7t3AihVyegEACAgARo6U95czluxs4NQpeT+7Y8eA1FQgKEherdeggfxao4bxjkdEZKF4u5MKxLBEenJz5f3l1qyRXXMA0LIl0L+/DE2qcjToJicDf/8NHD0KnD798OkK3N2BZs2Axx+XNwS2syv7MYmILBzDUgViWKJi3b0rA9OmTQXhxsMDeOopIDxcDgw3RKMBoqOB33+XXwv/qXp4yADWqhXg6QlcuACcPw+cOwdcuaJf1sZGlg0LA9q2ld2GRETEsFSRGJbIoKQk4NdfgR07CuZjUhQZdNq2BZyd5fxNTk6Ag4O8rcqePcDmzQXdeQAQEiJbilq1Anx9S76hb1aWDE6HDgFRUfJKPS1FAerWLeiya9BAdhXa2prs9ImIKiuGpQrEsESlkpMj7y3355+ln2LAyUm2QnXvLm/uW1ZCAAkJMjRFRQEXLxYtY2Ulr+oLDJRLQIBc2AJFRBaOYakCMSxRmSUlyavlLl+W45oKL/fuyZafnj3lrVWM2eqTkiK76rRddufOAWlpxZetWVMGNA8P2dXn7i6/enjIIOXsLLv4iIiqKIalCsSwREYlRMldbKY41s2bcszTpUtyiY/X77ozxM5OhiZnZxmg/Pxkq5S/v+zu48ByIqrEKiwsrV0LLF0qx5+mpMiLdVq2NLzNqVNy7r7oaPnB+tNPgQkT9MssWSKXS5fk902bym169JDfX7okewyK8+OPwMCBsqdj7lw5/OPWLfke/sorwPjx+uVXrgQ+/lh+yHZ1lcf45BP5Qbo0GJbI4mRlyT/OxET5x/PgkpamP4i8OIoCeHvL4KTt2gsIAHx8yndFIBGRkZX2/7f1ox4oI0P2FAwaJO8tWhqZmXJamIEDgYkTiy9Tp44MOg0ayPfkFSuAfv1kGGvaVH6AffDD77JlMuRoA1V0tLzg6IcfZPl9++RcgVZWwLhxsszevcDQoTKw9ekjb+/1yivyXNauLd/PhKjKs7MDGjWSS3GEkH/8d+8WLCkpMmBpW6hSU+UfaWKiHKulZWsrW538/GSY8vGRA9a9vQEXl4prVSMiKiWjdcNpW3pK07JUWECAbFV6sGWpODVryjA0alTxz7dqBbRuDXz7bcn7GDsWOHMG2L5dfj9vnmzBunChoMyiRcBHHwFXr5buHNiyRFSM1NSC4KQNUZcvy4HuJXF0lJ+QmjaV80Q1bMgr9YjIZCqsZaki5OcDP/0kP8iGhRVfJjpaTmi8eLHhfaWmytClFRYGvPOOnMqmRw95pfbPP8uxtUT0CFxd5XQHISEF6zQaObj90iXg+vWClidtd19GhvxDjomR5a2tZXhq3Fi2RNWpIxdnZ/1jCSGbrFNS5FK4xSs9XX7NzpbbubnpLw4OsrlZpSr4qlLJ8sUNvrexkS1vanXBVwcHOfDd2ZktY0QWqFKHpRMnZJjJypJXUK9bBzRpUnzZb7+V9zJt377k/e3bVzBHoFaHDnLM0uDB8jh5ebI77mGhi4jKQaWSXW6+vkWfy8mR/eCnT8uBjadOyeBz5oxcCnN2lqHJykqW+eefynHzYhsbGZq0i4uLfPNydCxYHBxk2fx8uWg08mtenjyHBxchZCjTLmo1YG8v912jhlw4kJ7IpMoUllauBF5+ueD7P/4AnnjC2FUq0KiR/ICZmipbe4YNA3btKhqY7t0DVq0C3nuv5H2dPCnHPM2YATz9dMH606flgO/p04Fu3eQH3MmT5bglQ915FSUrS35Vqws+sOblycXKSv/KbUNlVSr93oyylNW+X9vaFozLzc+Xd/V4lLI5OfL/hI2NPBdAfp+TU7ayiiLP48Gy1tZyKWtZIQr+7xb+H5SbK8/FGGWL+7mXpWxZXvtH+T0p7vU0xu+J9ueuV9bKFrk+gVDVDoRtr15yZ0lJyDl2BprzF2GTmACr61eAW7egSUtHzsnzUCCgtiq49UuOnQs0NdxhU8MJVq5OgJMTNI7OyLF3haK2hfreHfmGcvs2cv65C82dNFhnpcMaeYBGA5GXj+w8KyA/H3YOKt1koTl2LtDYO8LawRbWIhfIyoLIykZ2hgw4dvduy/3m5iL3WjLyr96CtZIPa5Wm4PXUyBfBziq34PdEY4V8oSqxrFqVW/AaaVTIE1awUjSwUeUXvJ75sqVL7e4Exc0VsLFBniIXKyvARl3wwmVlyRdFLbKg5OcBioI8WCNPsYHKSpGvkZUVoCjIypVf1Taaor8nioCtkqu7UXR2vjWEygq21hqorBRAUZAvVMjNV8my1hrdL1t2vrX8PbHKl2UB5GsU5AprqKxVsg73W/dyhA00ihVsbBVYWSuASgWNYoUcjTVUttawVSu6FsEcjTU0UMHGSgMrGxWgKNAIBTm5ivy7t8qTP1wh5O+fYiX/5mz096uoFKjtC1oZczTW0Ail4O9TUeRrlCPrbmd7/2bZQiA3T5F/9yoNrK1Ewev54HuEoiA3X4V8jaJXBwEF2bmqoj/3XCHfI6xEwXuERujqoLYVUFTycV6+gjyNSr72tgWtnFk59/er/btXFFk2Xyl4j7h/wKwc+TujVqNgv3ly37q/Ze3rmQ0IKLC1EVApQr+sIvTfI3IUWdYWBb8nhd8j1AX1Le49At7eBW+SFaxMR+3bF2jXruD78syRVxa2tnK6GQAIDZUTEn/+OfDVV/rlfv5ZtsAPHVr8fk6fBrp2lYO7331X/7k5c2Tr0uTJ8vsWLeT74xNPyBvJ+/gY95zKauBA+fWHHwrmCFy7Fvi//5Oh77XXCsq+8IL8xf3224I7aWzaBHzzDdCpE/DmmwVlR42SFzQtXizH2gLAtm3AF1/I17jwz+nVV2X35IIFskcEkPeLnT9fjk+bNaug7MSJ8m4bH34ING8u1x06BMyeLVv+Pv64oOyUKfIKxOnTgccek+uOH5ehNzAQWLiwoOyMGTLwvv22vKAAAOLigLfekq/RsmUFZefMAQ4fluPgunaV6y5dkqG4Zk15sYDWggVykP8rrwC9esl1iYnyQ4GjI7B6dUHZxYvlz2jECODZZ+W6lBRg+HD5vrp+fUHZb76RXbvPPw/8+99yXWYmMGSIfLxuXcHf/Pffy++feUbecxeQIUP72q9eLesCyCs9//c/2U38n/8UHG/IELlNZGTBVZwbNwLLl8ufQeExgcOHyx6lr74qaODZskVe1dqhg3xdtMaMkef4+efyogxAfmD57DN5y7kZMwrKvvaa/Nl9/LF8rQE5D+ZHH8nhR3PmFJR98005Q8GsWQVjHI8cAT74QP6OLVgA+Wbs44N3P/XBmTNdMG2anMAcWVk4teMm3pnlAL9a2fjyvynypGvWxKz/qhETA0z6N9C5s9zvhXPAG2/Iv4nCH4A+/i9w4Ly82KNbN7nuSoIc1+jiKj8cai2cJ8/7pZfkhy4AuJks/47UavkehNxcICUFSz/X4M+dtngx7DwGtYgF0tORdisXL6zsAeTn49enF8lzs7JC5JkO2BjfHIMan8CLzY4BajWyVQ4YuGYAoFLhp1d3ws4mH8jOxv/2BOLHo/XRt95pjG60R/4Bp6Rg4O63AAA/tJoP1yR51cva6x3xf1f/hac9j+K1wN905/HC4SnI1tjg25CFqKVOBQBsSmqHbxKeRif3k3iz3jpd2VFHJiEtzwGLmy9FXfubAIBtya3wxaXeaFcjDu82+FFX9tVjryE52w0LmnyDBk6yDrtvNcP8i8+gpetFzGpU8MOceOI/uHLPAx82/h7NXS4DAA7dboTZ5wYh2OkKPm4SqSs75dQonMvwxfSGq/GY2zkAwPHUQLwX9wICHW5gYbOCP/wZZ4bi5F1/vF3/F3SseRoAEHe3Dt46MwI+dilY1qLgH8ecs0Nw+E4DTAjagK4exwEAlzK8MP7UGNS0vYsVLefryi44/xz2pgTjFf8/0MvrMAAgMasmXj4+Fo7WWVjd+hNd2cUX+2LbrRCM8PsLz/pEAQBScpwxPGYCrBQN1j82W1f2m0s98HtyGzxf+2/8u/YuAEBmnhpDjsjXc12b2boA/X1CONYlheEZ7yiMrPsXACBfo8LAw9MAAKtbfwxHa5nIfrzWCf+79iR61jqM/wT8oTvekEPTkC9UiGz5Gdxt5Z0ENiaGYfmVcHT1OIYJQRt1ZYcfmYyMPDt81WIxfO1SAABbbrTB0ss90KHmGUyp/7Ou7JiYCUjJccbnTZchyPEGAGDXrRb47GI/tHE7hxkNC95EXzs+FolZNfFx8HIEO8tBwVEpTfDR+QFo5nwZc4K/15V98+QYxGd6YVajH9DSNf7+D+0bwMsL5lCmsKSdTsVcNJriW9q//VYGOU/Pos+dOgV06SJbpWbPLvp8ZmbRoKptueAMVESVkJ0dUMcPcAfgC6BZMV165mBjI9/IfSDr1s4dGHT/02UqgMP3yy1ZUrDN1wA2AhjUFHjxfpLOAnDo/vMjAgFta4QLgHQAfRsCo/sX7KNXPpCTC7z/PqC5LVPzXzWBLZ5Aaxegj798M7O2BmY1APKs5CcHj/stIjtdgI1eQIgHMDBQbi8EMKcRkGkt3zw978myRzyBX+sCjd2BoXUKPvJ/1BC4bSOTeO1M+WZ9rAbwiz8Q5A4MdSuo7xdNgJt2cuxDYLpcd8YVWF0XqF0DGDqyoHtyRWMg8f4s9gEhcv0lN+COF1DTTn5i1JbN8AOSagItmgMBNeW6GzWBGzUAZ5uCT1kqFZAdAFjVBJo2A+o5yX3ccgGuuQJ2NvJTs7Z79I4HkO0oP1kE3v/0nu4CXHAAbKzlJ3ptE9BdT+De/bIN719Jes8BOO8EqDQFV5cKAWTWAjIcZYoPCpLrcmyAM/ayjJ+f3EZRgAx3IM1edu0GBNzfhxVw+n6XbkAAcD8sIdsDSLGX3bN16hT8IzthBwhF/o6q73/6uucGJNsCzk6yHvdb3mCrBlS2ch8Oyv36OgO2NrILuEaNgtfTVg0Im/sT1d5vvs92uV/WQX5C1e3XFsi3ll3T2jCRdf8WT7a2cr12njkbG8DaSn5S1JY143jAR74aLiVF3k3h+nX5yXz1avn74O0tF0C2+NSuXfDJMidHtvYA8hNyRIRcnJwKWpKmTpUDruvWlWMzV62Sn1C3bJH3IdU6f15eMPP77/KOEIWdPCmDUrdu8io6LSurgmAVGSmnCVi4sKAbbsIE+fd04EDpfgamvBqO3XAPL8tuOAvphivH70lJr2d5f09Kej0r0+9JRb32fI8ouSzfIwrKVuR7hClU2KSUkZGyW+JBM2YAM2fKx507y+AbGSm/L2lCyU6dgJ075eNRo2SXR2KiDKwtWsgumMJBCZBXsv3wg9zngz/QmTPlh60H+fsXTHYJyKkCli6V3QJubjJgffRR6bsZOXUAERFR1cPbnVQghiUiIqKqp7T/v3nPASIiIiIDGJaIiIiIDGBYIiIiIjKAYYmIiIjIAIYlIiIiIgMYloiIiIgMYFgiIiIiMoBhiYiIiMgAhiUiIiIiAxiWiIiIiAxgWCIiIiIygGGJiIiIyABrc1fAEmjvRZyWlmbmmhAREVFpaf9va/+Pl4RhyQju3r0LAPDz8zNzTYiIiKis7t69C1dX1xKfV8TD4hQ9lEajwfXr1+Hs7AxFUcq8fVpaGvz8/HDlyhW4uLiYoIaVQ3U4T56jZeA5Wgaeo2Uw5TkKIXD37l34+vpCpSp5ZBJbloxApVKhTp06j7wfFxcXi/1lL6w6nCfP0TLwHC0Dz9EymOocDbUoaXGANxEREZEBDEtEREREBjAsVQJqtRozZsyAWq02d1VMqjqcJ8/RMvAcLQPP0TJUhnPkAG8iIiIiA9iyRERERGQAwxIRERGRAQxLRERERAYwLBEREREZwLBUCSxevBgBAQGws7NDu3btcPDgQXNXqdz+/vtv9OnTB76+vlAUBevXr9d7XgiB6dOnw8fHB/b29ggPD8e5c+fMU9lymjNnDh577DE4OzujVq1a6N+/P+Li4vTKZGVlYezYsXB3d4eTkxMGDBiAGzdumKnGZbdkyRK0aNFCNwlcWFgY/vjjD93zVf38ijN37lwoioIJEybo1lX185w5cyYURdFbGjdurHu+qp+f1rVr1/DCCy/A3d0d9vb2aN68OQ4fPqx7vqq/7wQEBBR5HRVFwdixYwFYxuuYn5+P9957D4GBgbC3t0e9evUwa9YsvXu2mfV1FGRWq1evFra2tuK7774Tp06dEqNHjxZubm7ixo0b5q5aufz+++9i2rRpYu3atQKAWLdund7zc+fOFa6urmL9+vXi2LFjom/fviIwMFDcu3fPPBUuh27duonly5eLkydPipiYGNGzZ09Rt25dkZ6erivzyiuvCD8/P7Ft2zZx+PBh8fjjj4v27dubsdZls3HjRrFp0yZx9uxZERcXJ9555x1hY2MjTp48KYSo+uf3oIMHD4qAgADRokULMX78eN36qn6eM2bMEE2bNhWJiYm65ebNm7rnq/r5CSFESkqK8Pf3F8OHDxcHDhwQFy9eFFu2bBHnz5/Xlanq7zvJycl6r+HWrVsFALFjxw4hhGW8jrNnzxbu7u7it99+E/Hx8eKnn34STk5O4vPPP9eVMefryLBkZm3bthVjx47VfZ+fny98fX3FnDlzzFgr43gwLGk0GuHt7S0++eQT3bo7d+4ItVot/ve//5mhhsaRnJwsAIhdu3YJIeQ52djYiJ9++klX5syZMwKAiIqKMlc1H1mNGjXEN998Y3Hnd/fuXdGgQQOxdetW0alTJ11YsoTznDFjhggJCSn2OUs4PyGEePvtt0XHjh1LfN4S33fGjx8v6tWrJzQajcW8jr169RIjR47UW/fss8+KiIgIIYT5X0d2w5lRTk4OoqOjER4erlunUqkQHh6OqKgoM9bMNOLj45GUlKR3vq6urmjXrl2VPt/U1FQAQM2aNQEA0dHRyM3N1TvPxo0bo27dulXyPPPz87F69WpkZGQgLCzM4s5v7Nix6NWrl975AJbzOp47dw6+vr4ICgpCREQEEhISAFjO+W3cuBFt2rTBwIEDUatWLbRq1Qpff/217nlLe9/JycnBDz/8gJEjR0JRFIt5Hdu3b49t27bh7NmzAIBjx45hz5496NGjBwDzv468ka4Z3bp1C/n5+fDy8tJb7+XlhdjYWDPVynSSkpIAoNjz1T5X1Wg0GkyYMAEdOnRAs2bNAMjztLW1hZubm17ZqnaeJ06cQFhYGLKysuDk5IR169ahSZMmiImJsYjzA4DVq1fjyJEjOHToUJHnLOF1bNeuHSIjI9GoUSMkJibi/fffxxNPPIGTJ09axPkBwMWLF7FkyRK88cYbeOedd3Do0CG8/vrrsLW1xbBhwyzufWf9+vW4c+cOhg8fDsAyfk8BYMqUKUhLS0Pjxo1hZWWF/Px8zJ49GxEREQDM//+DYYnoEYwdOxYnT57Enj17zF0Vo2vUqBFiYmKQmpqKn3/+GcOGDcOuXbvMXS2juXLlCsaPH4+tW7fCzs7O3NUxCe2ncgBo0aIF2rVrB39/f/z444+wt7c3Y82MR6PRoE2bNvjwww8BAK1atcLJkyexdOlSDBs2zMy1M75vv/0WPXr0gK+vr7mrYlQ//vgjVq5ciVWrVqFp06aIiYnBhAkT4OvrWyleR3bDmZGHhwesrKyKXLVw48YNeHt7m6lWpqM9J0s533HjxuG3337Djh07UKdOHd16b29v5OTk4M6dO3rlq9p52traon79+ggNDcWcOXMQEhKCzz//3GLOLzo6GsnJyWjdujWsra1hbW2NXbt2YeHChbC2toaXl5dFnGdhbm5uaNiwIc6fP28xr6OPjw+aNGmity44OFjX3WhJ7zuXL1/GX3/9hZdeekm3zlJex8mTJ2PKlCkYMmQImjdvjhdffBETJ07EnDlzAJj/dWRYMiNbW1uEhoZi27ZtunUajQbbtm1DWFiYGWtmGoGBgfD29tY737S0NBw4cKBKna8QAuPGjcO6deuwfft2BAYG6j0fGhoKGxsbvfOMi4tDQkJClTrPB2k0GmRnZ1vM+XXt2hUnTpxATEyMbmnTpg0iIiJ0jy3hPAtLT0/HhQsX4OPjYzGvY4cOHYpM3XH27Fn4+/sDsJz3HQBYvnw5atWqhV69eunWWcrrmJmZCZVKP5JYWVlBo9EAqASvo8mHkJNBq1evFmq1WkRGRorTp0+LMWPGCDc3N5GUlGTuqpXL3bt3xdGjR8XRo0cFALFgwQJx9OhRcfnyZSGEvPTTzc1NbNiwQRw/flz069evSl3CK4QQ//nPf4Srq6vYuXOn3uW8mZmZujKvvPKKqFu3rti+fbs4fPiwCAsLE2FhYWasddlMmTJF7Nq1S8THx4vjx4+LKVOmCEVRxJ9//imEqPrnV5LCV8MJUfXPc9KkSWLnzp0iPj5e7N27V4SHhwsPDw+RnJwshKj65yeEnPbB2tpazJ49W5w7d06sXLlSODg4iB9++EFXxhLed/Lz80XdunXF22+/XeQ5S3gdhw0bJmrXrq2bOmDt2rXCw8NDvPXWW7oy5nwdGZYqgUWLFom6desKW1tb0bZtW7F//35zV6ncduzYIQAUWYYNGyaEkJd/vvfee8LLy0uo1WrRtWtXERcXZ95Kl1Fx5wdALF++XFfm3r174tVXXxU1atQQDg4O4plnnhGJiYnmq3QZjRw5Uvj7+wtbW1vh6ekpunbtqgtKQlT98yvJg2Gpqp/n4MGDhY+Pj7C1tRW1a9cWgwcP1pt/qKqfn9avv/4qmjVrJtRqtWjcuLFYtmyZ3vOW8L6zZcsWAaDYelvC65iWlibGjx8v6tatK+zs7ERQUJCYNm2ayM7O1pUx5+uoCFFoekwiIiIi0sMxS0REREQGMCwRERERGcCwRERERGQAwxIRERGRAQxLRERERAYwLBEREREZwLBEREREZADDEhEREZEBDEtEVOlcunQJiqIgJibG3FXRiY2NxeOPPw47Ozu0bNmy2DKdO3fGhAkTKrRepaEoCtavX2/uahBVWQxLRFTE8OHDoSgK5s6dq7d+/fr1UBTFTLUyrxkzZsDR0RFxcXF6N/MsbO3atZg1a5bu+4CAAHz22WcVVENg5syZxQa5xMRE9OjRo8LqQWRpGJaIqFh2dnb46KOPcPv2bXNXxWhycnLKve2FCxfQsWNH+Pv7w93dvdgyNWvWhLOzc7mPUZJHqTcAeHt7Q61WG6k2RNUPwxIRFSs8PBze3t6YM2dOiWWKa8n47LPPEBAQoPt++PDh6N+/Pz788EN4eXnBzc0NH3zwAfLy8jB58mTUrFkTderUwfLly4vsPzY2Fu3bt4ednR2aNWuGXbt26T1/8uRJ9OjRA05OTvDy8sKLL76IW7du6Z7v3Lkzxo0bhwkTJsDDwwPdunUr9jw0Gg0++OAD1KlTB2q1Gi1btsTmzZt1zyuKgujoaHzwwQdQFAUzZ84sdj+Fu+E6d+6My5cvY+LEiVAURa9Fbs+ePXjiiSdgb28PPz8/vP7668jIyNA9HxAQgFmzZmHo0KFwcXHBmDFjAABvv/02GjZsCAcHBwQFBeG9995Dbm4uACAyMhLvv/8+jh07pjteZGSkrv6Fu+FOnDiBLl26wN7eHu7u7hgzZgzS09OLvGbz5s2Dj48P3N3dMXbsWN2xiKobhiUiKpaVlRU+/PBDLFq0CFevXn2kfW3fvh3Xr1/H33//jQULFmDGjBno3bs3atSogQMHDuCVV17Byy+/XOQ4kydPxqRJk3D06FGEhYWhT58++OeffwAAd+7cQZcuXdCqVSscPnwYmzdvxo0bNzBo0CC9faxYsQK2trbYu3cvli5dWmz9Pv/8c8yfPx/z5s3D8ePH0a1bN/Tt2xfnzp0DILuxmjZtikmTJiExMRFvvvnmQ8957dq1qFOnDj744AMkJiYiMTERgGyh6t69OwYMGIDjx49jzZo12LNnD8aNG6e3/bx58xASEoKjR4/ivffeAwA4OzsjMjISp0+fxueff46vv/4an376KQBg8ODBmDRpEpo2bao73uDBg4vUKyMjA926dUONGjVw6NAh/PTTT/jrr7+KHH/Hjh24cOECduzYgRUrViAyMlIXvoiqHUFE9IBhw4aJfv36CSGEePzxx8XIkSOFEEKsW7dOFH7bmDFjhggJCdHb9tNPPxX+/v56+/L39xf5+fm6dY0aNRJPPPGE7vu8vDzh6Ogo/ve//wkhhIiPjxcAxNy5c3VlcnNzRZ06dcRHH30khBBi1qxZ4umnn9Y79pUrVwQAERcXJ4QQolOnTqJVq1YPPV9fX18xe/ZsvXWPPfaYePXVV3Xfh4SEiBkzZhjcT6dOncT48eN13/v7+4tPP/1Ur8yoUaPEmDFj9Nbt3r1bqFQqce/ePd12/fv3f2i9P/nkExEaGqr7vrjXQwghAIh169YJIYRYtmyZqFGjhkhPT9c9v2nTJqFSqURSUpIQouA1y8vL05UZOHCgGDx48EPrRGSJrM0b1Yiosvvoo4/QpUuXUrWmlKRp06ZQqQoasr28vNCsWTPd91ZWVnB3d0dycrLedmFhYbrH1tbWaNOmDc6cOQMAOHbsGHbs2AEnJ6cix7tw4QIaNmwIAAgNDTVYt7S0NFy/fh0dOnTQW9+hQwccO3aslGdYeseOHcPx48excuVK3TohBDQaDeLj4xEcHAwAaNOmTZFt16xZg4ULF+LChQtIT09HXl4eXFxcynT8M2fOICQkBI6Ojrp1HTp0gEajQVxcHLy8vADI18zKykpXxsfHBydOnCjTsYgsBcMSERn05JNPolu3bpg6dSqGDx+u95xKpYIQQm9dceNabGxs9L5XFKXYdRqNptT1Sk9PR58+ffDRRx8Vec7Hx0f3uHAoqAzS09Px8ssv4/XXXy/yXN26dXWPH6x3VFQUIiIi8P7776Nbt25wdXXF6tWrMX/+fJPU81FfHyJLwrBERA81d+5ctGzZEo0aNdJb7+npiaSkJAghdAOYjTk30v79+/Hkk08CAPLy8hAdHa0bW9O6dWv88ssvCAgIgLV1+d/KXFxc4Ovri71796JTp0669Xv37kXbtm0fqf62trbIz8/XW9e6dWucPn0a9evXL9O+9u3bB39/f0ybNk237vLlyw893oOCg4MRGRmJjIwMXSDbu3cvVCpVkdeXiCQO8Caih2revDkiIiKwcOFCvfWdO3fGzZs38fHHH+PChQtYvHgx/vjjD6Mdd/HixVi3bh1iY2MxduxY3L59GyNHjgQAjB07FikpKXj++edx6NAhXLhwAVu2bMGIESMeGhgeNHnyZHz00UdYs2YN4uLiMGXKFMTExGD8+PGPVP+AgAD8/fffuHbtmu4qvbfffhv79u3DuHHjEBMTg3PnzmHDhg1FBlg/qEGDBkhISMDq1atx4cIFLFy4EOvWrStyvPj4eMTExODWrVvIzs4usp+IiAjY2dlh2LBhOHnyJHbs2IHXXnsNL774oq4Ljoj0MSwRUal88MEHRbphgoOD8eWXX2Lx4sUICQnBwYMHH2ls04Pmzp2LuXPnIiQkBHv27MHGjRvh4eEBALrWoPz8fDz99NNo3rw5JkyYADc3N73xUaXx+uuv44033sCkSZPQvHlzbN68GRs3bkSDBg0eqf4ffPABLl26hHr16sHT0xMA0KJFC+zatQtnz57FE088gVatWmH69Onw9fU1uK++ffti4sSJGDduHFq2bIl9+/bprpLTGjBgALp3745//etf8PT0xP/+978i+3FwcMCWLVuQkpKCxx57DM899xy6du2KL7744pHOlciSKeLBAQdEREREpMOWJSIiIiIDGJaIiIiIDGBYIiIiIjKAYYmIiIjIAIYlIiIiIgMYloiIiIgMYFgiIiIiMoBhiYiIiMgAhiUiIiIiAxiWiIiIiAxgWCIiIiIy4P8BFN+rN7DOGpMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "result = numpy.load('./output/summary_data.npz')\n",
    "\n",
    "eig_val, eig_state = numpy.linalg.eig(\n",
    "                     molecular_hamiltonian.construct_h_matrix())\n",
    "min_eig_H = numpy.min(eig_val.real)\n",
    "min_loss = numpy.ones([len(result['iter'])]) * min_eig_H\n",
    "\n",
    "plt.figure(1)\n",
    "func1, = plt.plot(result['iter'], result['energy'], \n",
    "                  alpha=0.7, marker='', linestyle=\"-\", color='r')\n",
    "func_min, = plt.plot(result['iter'], min_loss, \n",
    "                  alpha=0.7, marker='', linestyle=\":\", color='b')\n",
    "plt.xlabel('Number of iteration')\n",
    "plt.ylabel('Energy (Ha)')\n",
    "\n",
    "plt.legend(handles=[\n",
    "    func1,\n",
    "    func_min\n",
    "],\n",
    "    labels=[\n",
    "        r'$\\left\\langle {\\psi \\left( {\\theta } \\right)} '\n",
    "        r'\\right|H\\left| {\\psi \\left( {\\theta } \\right)} \\right\\rangle $',\n",
    "        'Ground-state energy',\n",
    "    ], loc='best')\n",
    "plt.text(-15.5, -1.145, f'{min_eig_H:.5f}', fontsize=10, color='b')\n",
    "#plt.savefig(\"vqe.png\", bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 通过 VQE 确定原子间隔\n",
    "\n",
    "还记得在前面的注释中提到我们默认使用的两个氢原子间原子间隔为 $74$ pm 吗？VQE 的另一个用法便是通过在不同的原子间隔下多次运行然后观察运行结果的最小值是在什么原子间隔发生的，这个间隔即为估计得真实原子间隔。\n",
    "\n",
    "![vqe-fig-dist](figures/vqe-fig-distance.png)\n",
    "\n",
    "从上图可以看出，最小值确实发生在 $d = 74$ pm (1 pm = $1\\times 10^{-12}$ m) 附近，这是与[实验测得数据](https://cccbdb.nist.gov/exp2x.asp?casno=1333740&charge=0)相符合的 $d_{exp} (H_2) = 74.14$ pm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "[1] Cao, Yudong, et al. Quantum Chemistry in the Age of Quantum Computing. [Chemical reviews 119.19 (2019): 10856-10915.](https://pubs.acs.org/doi/10.1021/acs.chemrev.8b00803)\n",
    "\n",
    "[2] McArdle, Sam, et al. Quantum computational chemistry. [Reviews of Modern Physics 92.1 (2020): 015003.](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.92.015003)\n",
    "\n",
    "\n",
    "[3] Peruzzo, A. et al. A variational eigenvalue solver on a photonic quantum processor. [Nat. Commun. 5, 4213 (2014).](https://www.nature.com/articles/ncomms5213)\n",
    "\n",
    "[4] Moll, Nikolaj, et al. Quantum optimization using variational algorithms on near-term quantum devices. [Quantum Science and Technology 3.3 (2018): 030503.](https://iopscience.iop.org/article/10.1088/2058-9565/aab822)\n",
    "\n",
    "[5] 徐光宪, 黎乐民, 王德民. 量子化学: 基本原理和从头计算法(上)[M], 第二版. 北京: 科学出版社, 2012; \n",
    "\n",
    "[6] Helgaker, Trygve, Poul Jorgensen, and Jeppe Olsen. Molecular electronic-structure theory. John Wiley & Sons, 2014.\n",
    "\n",
    "[7] Dirac, Paul Adrien Maurice. Quantum mechanics of many-electron systems. [Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 123.792 (1929): 714-733.](https://royalsocietypublishing.org/doi/10.1098/rspa.1929.0094)\n",
    "\n",
    "[8] Szabo, Attila, and Neil S. Ostlund. Modern quantum chemistry: introduction to advanced electronic structure theory. Courier Corporation, 2012."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}