numpy_1ddft_cn.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy 一维DFT\n",
    "\n",
    "## 理论背景\n",
    "\n",
    "- 宏观世界：经典力学\n",
    "- 纳米尺度：量子力学\n",
    "\n",
    "量子力学用薛定谔方程描述物理体系：\n",
    "\n",
    "$$\n",
    "\\hat{H} \\Psi=[\\hat{T}+\\hat{V}+\\hat{U}] \\Psi=\\left[\\sum_{i=1}^N\\left(-\\frac{\\hbar^2}{2 m_i} \\nabla_i^2\\right)+\\sum_{i=1}^N V\\left(\\mathbf{r}_i\\right)+\\sum_{i<j}^N U\\left(\\mathbf{r}_i, \\mathbf{r}_j\\right)\\right] \\Psi=E \\Psi \\tag{1}\n",
    "$$\n",
    "\n",
    "式1计算复杂度随电子数量N快速上升。密度泛函理论（DFT）大幅降低了计算量（将多电子体系的薛定谔方程转化成单电子形式KS方程，将精度损失全部归入交换关联泛函），使得通过数值计算分析纳米材料物化性质成为可能。\n",
    "\n",
    "$$\n",
    "\\left[-\\frac{\\hbar^2}{2 m} \\nabla^2+V_{\\mathrm{s}}(\\mathbf{r})\\right] \\varphi_i(\\mathbf{r})=\\varepsilon_i \\varphi_i(\\mathbf{r}) \\tag{2}\n",
    "$$\n",
    "\n",
    "$$\n",
    "V_{\\mathrm{s}}(\\mathbf{r})=V(\\mathbf{r})+\\int \\frac{n\\left(\\mathbf{r}^{\\prime}\\right)}{\\left|\\mathbf{r}-\\mathbf{r}^{\\prime}\\right|} \\mathrm{d}^3 \\mathbf{r}^{\\prime}+V_{\\mathrm{XC}}[n(\\mathbf{r})] \\tag{3}\n",
    "$$\n",
    "\n",
    "$$\n",
    "n(\\mathbf{r})=\\sum_{i=1}^N\\left|\\varphi_i(\\mathbf{r})\\right|^2 \\tag{4}\n",
    "$$\n",
    "\n",
    "式2到式4所示KS方程聚焦体系能量和电子密度的关系，二者一一对应（原则上，其他物化性质在求出这两项后可另行导出）。在量子力学中，能量与哈密顿量（算符）对应，而电子密度与波函数对应。**问题归结为如何正确写出作用在体系波函数上的哈密顿量（式2括号中内容）并求解KS方程**（更多细节请阅读 [wiki](https://en.wikipedia.org/wiki/Density_functional_theory) 或相关书籍）。\n",
    "\n",
    "## 代码概述\n",
    "\n",
    "- 理论部分见[此PDF](./Larsen_Lyon_Theory.pdf)\n",
    "- 目标: 针对一维谐振子，考虑动能、静电势、LDA，忽略关联项；写出这种情况下体系的KS方程并求解\n",
    "\n",
    "根据KS方程通用形式，哈密顿量由 动能和势能 组成：\n",
    "\n",
    "$$\\hat{H}=-\\frac{1}{2}\\frac{d^2}{dx^2}+v_s(x)$$\n",
    "\n",
    "势能包含 外部势（谐振子为$x^2$）、电子库伦势（Hartree）和电子交互关联势（局域密度近似LDA）\n",
    "\n",
    "$$v_s(x)=x^2+v_{Hartree}(x)+v_{LDA}(x)$$\n",
    "\n",
    "如何代码实现并求解上述方程？\n",
    "\n",
    "## 逐个击破\n",
    "\n",
    "顺序：\n",
    "- 动能\n",
    "- 外部势（一维谐振子，就是$x^2$）\n",
    "- 电子库伦势 Hartree\n",
    "- 电子交互关联式 LDA\n",
    "\n",
    "注意，本项目并不写for循环自洽求解，而是调用第三方库直接求解已经写好的哈密顿量的本征值和本征函数。问题只有一个：如何写出哈密顿量？\n",
    "\n",
    "### 先数值表示动能（算符）\n",
    "\n",
    "如何数值表示二阶微分？\n",
    "\n",
    "二阶微分基于一阶微分，一阶导的数值表示形式：\n",
    "\n",
    "$$(\\frac{dy}{dx})_{i}=\\frac{y_{i+1}-{y_{i}}}{h}$$\n",
    "\n",
    "然而我们没有网格上的具体值（求解KS方程之前未知），要想办法将算符与网格值剥离（变成向量相乘的形式），类似这样：\n",
    "\n",
    "$$(\\frac{dy}{dx})_{i}=D_{ij} y_{j}$$\n",
    "\n",
    "为此，我们可以定义一阶导算符：\n",
    "\n",
    "$$D_{ij}=\\frac{\\delta_{i+1,j}-\\delta_{i,j}}{h}$$\n",
    "\n",
    "\n",
    "- 端点处的导数不好定义，暂时忽略这个不算太大的问题\n",
    "- $\\delta_{ij}$ 是 Kronecker delta\n",
    "- 使用了爱因斯坦求和约定"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d_i1j\n",
      " [[0. 1. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 1.]\n",
      " [0. 0. 0. 0. 0.]]\n",
      "d_ij\n",
      " [[1. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 1.]]\n",
      "D\n",
      " [[-0.4  0.4  0.   0.   0. ]\n",
      " [ 0.  -0.4  0.4  0.   0. ]\n",
      " [ 0.   0.  -0.4  0.4  0. ]\n",
      " [ 0.   0.   0.  -0.4  0.4]\n",
      " [ 0.   0.   0.   0.  -0.4]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "n_grid_demo=5\n",
    "x=np.linspace(-5,5,n_grid_demo)\n",
    "h=x[1]-x[0]\n",
    "\n",
    "d_i1j = np.diagflat(np.ones(n_grid_demo-1),1)\n",
    "print(f'd_i1j\\n', d_i1j)\n",
    "d_ij = np.eye(n_grid_demo)\n",
    "print(f'd_ij\\n', d_ij)\n",
    "D=(d_i1j - d_ij)/h\n",
    "print(f'D\\n', D)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "同样的我们想把二阶微分算符剥离出来：\n",
    "\n",
    "$$(\\frac{d^2y}{dx^2})_{i}=D^2_{ij} y_{j}$$\n",
    "\n",
    "二阶算符就应该表示成:\n",
    "\n",
    "$$D^2_{ij}=\\frac{\\delta_{i+1,j}-2\\delta_{i,j}+\\delta_{i-1,j}}{h^2}$$\n",
    "\n",
    "它和一阶算符之间存在一些运算关系:\n",
    "\n",
    "$$D^2_{ij}=-D_{ik}D_{jk}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d_in1j\n",
      " [[0. 0. 0. 0. 0.]\n",
      " [1. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 1. 0.]]\n",
      "D2ij\n",
      " [[-0.32  0.16  0.    0.    0.  ]\n",
      " [ 0.16 -0.32  0.16  0.    0.  ]\n",
      " [ 0.    0.16 -0.32  0.16  0.  ]\n",
      " [ 0.    0.    0.16 -0.32  0.16]\n",
      " [ 0.    0.    0.    0.16 -0.32]]\n",
      "D2\n",
      "[[-0.32  0.16  0.    0.    0.  ]\n",
      " [ 0.16 -0.32  0.16  0.    0.  ]\n",
      " [ 0.    0.16 -0.32  0.16  0.  ]\n",
      " [ 0.    0.    0.16 -0.32  0.16]\n",
      " [ 0.    0.    0.    0.16 -0.16]]\n",
      "[[-0.32  0.16  0.    0.    0.  ]\n",
      " [ 0.16 -0.32  0.16  0.    0.  ]\n",
      " [ 0.    0.16 -0.32  0.16  0.  ]\n",
      " [ 0.    0.    0.16 -0.32  0.16]\n",
      " [ 0.    0.    0.    0.16 -0.32]]\n"
     ]
    }
   ],
   "source": [
    "d_in1j = np.diagflat(np.ones(n_grid_demo-1),-1)\n",
    "print(f'd_in1j\\n', d_in1j)\n",
    "D2ij = (d_i1j - 2*d_ij + d_in1j)/ h**2\n",
    "print(f'D2ij\\n', D2ij)\n",
    "\n",
    "D2=D.dot(-D.T)\n",
    "print('D2')\n",
    "print(D2)\n",
    "D2[-1,-1]=D2[0,0] # 由于周期性边界条件，最后一个元素需要特殊处理\n",
    "print(D2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "采用最经典的正弦函数例子，验证上述两个算符是否正确"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x119aef4d0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADEKElEQVR4nOydd3wc5Z3/37O72lXvvVuSbbl3y71hY2MCOJ0EQuAI3CWBuxx3JPHvEpIjuSOF43JHSIcD0oAk9GIw7kWWe5eLqtV777s7vz9mZyXZkiXZuzuzu8/79drXrqXZma+k9TOf51slWZZlBAKBQCAQCHwIg9YGCAQCgUAgELgaIXAEAoFAIBD4HELgCAQCgUAg8DmEwBEIBAKBQOBzCIEjEAgEAoHA5xACRyAQCAQCgc8hBI5AIBAIBAKfQwgcgUAgEAgEPodJawO0wG63U11dTVhYGJIkaW2OQCAQCASCcSDLMh0dHSQnJ2MwXN9H45cCp7q6mrS0NK3NEAgEAoFAcANUVFSQmpp63WP8UuCEhYUByi8oPDxcY2sEAoFAIBCMh/b2dtLS0pz38evhlwJHDUuFh4cLgSMQCAQCgZcxnvQSkWQsEAgEAoHA5xACRyAQCAQCgc8hBI5AIBAIBAKfQwgcgUAgEAgEPocQOAKBQCAQCHwOIXAEAoFAIBD4HELgCAQCgUAg8DmEwBEIBAKBQOBzCIEjEAgEAoHA53CrwNm7dy933HEHycnJSJLEm2++OeZ7du/ezfz587FYLOTk5PDiiy9ec8xzzz1HZmYmgYGB5OXlcfjwYdcbLxAIBAKBwGtxq8Dp6upizpw5PPfcc+M6vrS0lNtvv521a9dy8uRJvvGNb/CVr3yFDz/80HnMq6++ymOPPcb3vvc9jh8/zpw5c9i4cSP19fXu+jEEAoFAIBB4GZIsy7JHLiRJvPHGG2zZsmXUY771rW/x3nvvcfbsWefX7r77blpbW9m2bRsAeXl5LFq0iJ///OcA2O120tLSePTRR/n2t789Llva29uJiIigra1NzKISCAQCgcBLmMj9W1fDNvPz81m/fv2wr23cuJFvfOMbAPT393Ps2DG2bt3q/L7BYGD9+vXk5+ePet6+vj76+vqc/25vb3et4TdAe+8AV5q6qWjupqKlm64+G3PTIlmYGUVYYIDW5gl0iizLNPQ0UNRaRG1XLRnhGUyLnkZwQLDWpmlHVyPUnoH682AJg8RZEDcNAgK1tkwXyLLMleZuTlW2caGmnegQM5MTwpiSEEpieOC4hhYK/JOG7gaO1x9HRiYtNI2U0BQiLBFe85nRlcCpra0lISFh2NcSEhJob2+np6eHlpYWbDbbiMdcuHBh1PM+9dRT/Pu//7tbbJ4IPf02/nConBcOlFLT1jviMQYJZiRHkDcpmjvmJDMnLdKzRgp0R0tvCy+ff5njdce53HqZjv6OYd83SAayIrKYGTuTZcnL2Ji5EYPkw/UDdjuc+Quc/ZsibDqqrz1GMkLcVEieD0u/BgkzPG+nhlhtdl47Wsm2c7WcrmyltXtgxOPCAk3MTo3ga2tyWJ4T62ErBXqjs7+TXRW7OFp3lKO1R7nSceWaY0ICQsgIz+AzUz7DluwtBBj1uyHXlcBxF1u3buWxxx5z/ru9vZ20tDSPXb93wMYfC67wy93FNHYOepJiQsykRQeTFh1MgFHiaFkLV5q7OVPVxpmqNn63v5QvLE7jW5tyiQw2e8xegT7osfbwx8I/8vyZ5+kc6HR+3SgZSQtLIykkieK2Yuq76ylqLaKotYg3i97kxXMvsnXxVubGz9XOeHdRth8+/DeoOTn869HZiojpa4ea09DTrHh06s/DqT/B/Ptg7b9BaLwmZnuSg8WNPPnOeS7UDgphs9HAtKQwpidH0Nrdz6W6DsqauunotXKgqIkDRU2smhLHtzflMj1ZhO39DVmW+bDsQ3585Mc09jQ6vy4hkRudS6ApkMqOShp6Guga6OJ803mezH+S353+HQ/Pfpg7c+4kwKA/oaMrgZOYmEhdXd2wr9XV1REeHk5QUBBGoxGj0TjiMYmJiaOe12KxYLFY3GLzWLx5oor/fL+Q+g5F2KRGBfHouhw+MTuZEMu1v/6ath4OlzbzcWE975yq5s+HK/joXB3/dvs0PjkvxWtcg4Ibx2a38VbxWzx34jnqe5Tk+dzoXL40/UtMjZpKZkQmFuPg57mhu4GzjWc52XCS1y6+xvmm83zpgy+xedJm/nnBP5MYMvr/Da+hqRi2PwEX3lX+bQ6DZY9C1hpImK6EplRkGdqrFe/OqT/B+bfg2Itw5q+w8jFY8jUICNLip3ArFc3d/Of7hXxwthaAyOAA/mF1NsuyY5iaGIbFZBx2fJ/VRnF9F68dreAPh8rZe6mBfZcb+OTcFB7fNJWkCN/7HQmu5Ur7Ff6j4D84WH0QgNTQVDZkbGBh4kLmxs8l3DwoeHutvVR3VnOw+iDPn32e6q5qvp//fX575rf8/ey/5xPZn9CV0NFdkvH777/PmTNnnF/74he/SHNz87Ak48WLF/Pss88CSpJxeno6jzzyiK6SjGVZ5pntl3h2ZxEAKZFBPLIuh0/PT8VsGl/4oKCkie+8eZbL9crufVl2DM98bi6JESK3wFfpHujmsd2PcaD6AADJIck8Ov9RNk/aPK6wU2NPIz8/8XNev/w6MjJBpiC2Lt7KJyd/0t2mu48zf4U3vwq2fpAMsOABWLMVQuPG9/7yfPjw/0H1ceXfsVPg3tch0nNeXHfz1skqvvnX0/RZ7RgkuHdJBo9tmDJuz295Uxc//fAi756uASA21MyLDyxmZkqEO80WaIjNbuP5s8/z61O/pt/ej9lg5iuzv8KDMx/EbBz7c9Nr7eUvl/7C82eep6m3CYC/3vFXpkZPdavdE7l/u1XgdHZ2UlSk3ODnzZvHM888w9q1a4mOjiY9PZ2tW7dSVVXFyy+/DChl4jNnzuTrX/86f/d3f8fOnTv5x3/8R9577z02btwIKGXiX/7yl/n1r3/N4sWL+dnPfsZrr73GhQsXrsnNGQ13C5x+q51v/+00r5+oAuDra7P5p1umjFvYXH2u3+0v4X93XKZ3wE5GTDB/fmgJyZFid+VrtPS28PUdX+dM4xmCTEE8MvcR7s69e1yLzdWcazrHjw//mBP1JwD45qJv8qXpX3K1ye7nyO/gvX8FZMhaC5uegvhpEz+P3a7k7Gz/LnTUQHgq3PcmxE52tcUe508FV/i3N88gy7AkK5rv3zmD3MQbW9dOV7by7b+d4XxNO6EWE7/50gKWidwcn8Nmt/HEwSd4u/htAJYkLeE7S75DRnjGhM/VY+3htYuvUd5ezhNLn3C1qdegG4Gze/du1q5de83Xv/zlL/Piiy9y//33U1ZWxu7du4e955//+Z85f/48qampfPe73+X+++8f9v6f//zn/PSnP6W2tpa5c+fyv//7v+Tl5Y3bLncKnLaeAf7h98fIL2nCaJD4z0/O5POL0m/6vGWNXXzphQIqmntIiw7izw8tITXKjytnfIyazhoe3v4wZe1lRFoiee6W55gdN/umzinLMv997L/5v3P/B8DX536dv5/9994R5pRl2Pc07Pyh8u9FD8FtPwHDTSZPt1XB77dA4yUIjlE8Oclzb9Zazfjt3hL+4/1CAL60JIN/v3MGBsPN/X07egf4+98f42BxE2ajgZ/dPZfNs5JcYa5AB1jtVv5t/7/xfun7GCUj31v6PbbkbPGOdQEdCRy94i6BU9XawwP/d5hLdZ2EmI384t4FrJ4yTjf6OKhu7eELvz1EeVM3KZFBvPLwEtKihcjxdopbi3l4+8PUd9eTGJLIrzf8mqyILJecW5Zlfn361zx3Umm2ef+M+3lswWP6XsxkGT76DuQrva5Y9U1Y+//AVTZ3NcEfPqUkKpvD4IuvQuZy15zbQ8iyzM8+vsz/7LgMwN+vzuLbm3Jd9nfts9r451dP8v6ZWiQJfnDXTO5dMvHdvUBfWO1Wtu7byraybZgkEz9Z/RM2ZGzQ2qwJMZH7tw/Xknqe5/eVcqmuk/gwC6/9w1KXihuA5MggXn14KZNiQ6hq7eHu3xyivKnLpdcQeJaytjK+vO3L1HfXkxWRxe9v+73LxA0ouW//MOcfeHzh4wC8eO5Ffnjoh+h2XyPL8N6/DIqbjU/Bun9znbgBCImBL78DGSugv0MROyW7XXd+D/CjbRec4uZfb53iUnEDYDEZefYL87knLx1Zhu+8eZZXj1xbMizwHgbsA3xz7zcVcWMw8fSap71O3EwUIXBcyLdvy+WevHTe+PpyZiS7JzkvMSKQVx5eQlacInK++NsCWrr63XItgXvpHujmG7u+QVtfGzNjZvLSppfcVvF034z7+N7S7yEh8dql13jp3Etuuc5Nc+z/4OjzSjLxXb9Qeti4g8BwuPevMOU2sPbCXx6AVu+4gf/1WCW/3lMCwBOfmM4j6ya7xSNnNEj8cMtMvromG4DvvnWOUxWtLr+OwP3Issx3D3yX7eXbCTAE8N9r/ptb0m/R2iy3IwSOCzGbDPzHJ2eR4uYE4IRwReRkxgRT1drD4389pd8duWBE1AWnuK2YuKA4nr3lWSIDI916zc9M+Qz/L+//AfCz4z/jWN0xt15vwlQdgw++pby+5Xsw7x73Xi8gCD77IiTNVfrmvPZlsPaN9S5NuVDbznfeVKpMv7F+Mn+3YpJbrydJEo/fOpUN0xPot9r56h+O0dSp79+R4Fr+cukvvFfyHibJxM/W/ow1aWu0NskjCIHjpcSHBfLcPfMxmwx8XFjP8/tLtTZJMAFePv8yH5V/hMlg4pk1zxAb5JlKlc9P/Ty3Z92OTbbx+J7HhzX10pSuJkVg2Poh9xOw/J88c92AQPjcyxAYqZSRbxtfqwkt6Ogd4Kt/OE7vgJ1VU+L4x3WeqQAzGCT+63NzyIoNobqtl3985QRWm90j1xbcPBeaL/Djwz8G4BsLvsGq1FUaW+Q5hMDxYmYkR/DdT0wH4MfbLnBSuI+9gsM1h3nm2DMAfGvRtzzacViSJJ5Y8gTZEdk09DTwrb3fwma3eez6I2K3wetfgbYKiM6CLb9wbc7NWERlwKd/B0hw9AU4+WfPXXucyLLMt/92htLGLpIiAvnZ5+fedLXURAgPDODXX1pAsNnIgaImnv7okseuLbhxuga6eHzP4/Tb+1mdupr7pt+ntUkeRQgcL+fevHQ2z0pkwCbzyJ+O09Yz8swZgT6o7arl8b2PY5ft3Jl9J5+f+nmP2xAcEMwza54hyBTE4drDzgorzdjzYyjeCaYg+NzvIVCD5nKTN8BqR3js3W8oXZB1xEsHy3jvTA0mg8Rz98wnOsTzo1smJ4Tx08/MAeBXe4r54EyNx20QjB9ZlvnBoR9Q1l5GQnACP1z+Q31XT7oBIXC8HEmS+NGnZ5MWHURlSw/f+utpkY+jU+yynW/v+zbNvc3kRufy3SXf1WzByYrM4t+XKQNof3vmt+yt3KuJHRTvUgQOwB3/A4kztbEDFIGTs15JOn7tPujv1s6WIZysaHX2uvm326cxPz1KM1tun53Ew6uUKr9v/vU09e0jDw0WaM8bRW/wXsl7GCUjP139U7fn+OkRIXB8gPDAAJ774nwCjBLbztXyp8PeUQ3ib7xV9BbH6o4RZArimTXPEGjSduTGbZNu4+6pdwPw/YPfp7O/c4x3uJiBHnj3n5XXC+6HOZ73Zg3DYIBP/RbCU6C5BPb+VFt7gAGbnW/99TQDNpnbZyVx/7JMrU3imxunMic1go4+Kz98r1BrcwQjUNJawlMFTwHwyLxHmBc/T2OLtEEIHB9hdmok39qUC8CPP7ggKh10RktvizPv5mtzvkZamD7mIP3ron8lPSydhp4GfnHqF569+P6fQUsphCXBhh949tqjERwNmx3C5uCz0HBRU3NeOljGxboOooID+OGWmboIMZiMSrWoQYK3T1Wz91KD1iYJhiDLMv9R8B/02npZmrSUv5v5d1qbpBlC4PgQDyyfxLSkcNp7rTz9kbYLs2A4zxx7hta+VqZETeGe6W4uf54AFqPFWTr+p8I/canFQ8mjTcWwXxF8bPxPpS+NXpi6GaZsAvuA0nRQo5BvbVsv/71d+Xt8+7ZcojTIuxmNmSkRfNnhTfruW2fpHdA4UV3g5MPyDzlcexizwcwTS58Y15BeX8V/f3IfxGiQePKuGQC8cqRCNOXSCUdrj/Jm0ZsAfHfJdwkwBGhr0FUsT1nOhowN2GQbPzz0Q+yym0uAZRne/1elJDx7HczQ2aRzSYLbfqwkPZftg9OvaWLGD987T1e/jfnpkXx2gT48fkP5l1unkhgeSHlTN7/YXay1OQKU5qFPH3kagAdnPUhqWKrGFmmLEDg+xqLMaD45LwVZhifePofdLhKOtWTANsAPDinhl89M+YxHS8InwjcXfZMgUxAn6k84Jwy7jXNvKFVTRgtsftqzJeHjJSoTVv2r8vqjf4OeFo9e/kBRI++ersEgwQ+2zPRoSfh4CbWY+N4dSpuKX+0uprjBwzlcgmv43ZnfUdddR0poil+HplSEwPFBtt6WS4jZyKmKVv56vFJrc/yal86/RElbCdGB0Xxj/je0NmdUEkMS+eqcrwLw38f+m7a+NvdcqLcdtm1VXq/4Z4jJds91XMGyf4TYKdDVMDjV3AP0WW18962zANy3NNNtY19cwaaZiaydGke/zc533zwrKjg1pLy9nBfPvQjA44se17yIQQ8IgeODxIcH8o31UwAl4Vj0xtGGqs4qfnXqVwD868J/JcKi3xsVwL3T7yU7Ipvm3maePfGsey6y+ynorIWoSYrA0TMmM9z+X8rrI88royQ8wO/2lVLS0EVsqIV/3jDFI9e8USRJ4sm7ZmIxGThY3MTbp6q1NskvkWWZHx3+EQP2AZYnL2dd2jqtTdIFQuD4KF9elkl2XAhNXf3OREWBZ/nlyV/SZ+tjUeIiPpH1Ca3NGZMAQwD/tuTfAHjt4msUNrm4BLi5BAp+rby+/WllTILembQKZn0OkOHD77g94biuvZdndypTwv/f5lwigvSVrzUSadHBPLI2B4Bntl9iQIxx8Dh7Kvewv2o/JoOJby3+li6q7fSAEDg+itlk4Pt3KgnHvz9UTomIj3uUsrYy3il5B4B/nv/PXrPgLEpcxG2ZtyEj84uTLi4b3/MTkG1KM72c9a49tztZ/30lX+jKQSjZ7dZL/XJ3Mb0DduanR/LJeSluvZYreXDlJGJDzZQ3dfO3YyIs7kkG7AP85MhPALhv+n1MinDvAFZvQggcH2bl5DjW5cZjs8v8fFeR1ub4Fb849Qvssp01qWuYFTdLa3MmxNfmfg2DZGB35W7ONZ5zzUkbL8PpV5XXa/+fa87pKSJSYOEDyutd/+k2L05tW6+zSedjG6Z6jSgGCDab+IfVSj7VszuL6LOKsnFP8W7xu1R0VBAdGM3fz/57rc3RFULg+DjfWK9MHH7rZDVljV0aW+MfXG65zLbSbQB8fd7XNbZm4mRGZDpDai5r/rf7RyDbYcptkLLANef0JCv+GUyBUHkYij52yyV+sbuIfqudxZnRLM+Jccs13Mm9SzJICLdQ1drDa0cqtDbHL7Darfz2zG8BuH/G/QQHBGtskb4QAsfHmZ0aydqpcdjsMs/uFF4cT/CLk79ARmZDxgZyo3O1NueG+PvZf49RMrK3ci+nG07f3MnqC+Hs35TXa7fevHFaEJYIi76ivN71Hy734lS39vDKYUUUfGPDZK/y3qgEBhiduTjP7iwSzf88wHsl71HRUUGUJUqTwb16RwgcP+CfHBVVb56sorxJeHHcybmmc3x85WMkJL4+1/u8Nyrp4enckX0H4AIvzu4fATJMuwOS5ty8cVqx/BsQEALVJ+DSNpee+rldRfTb7ORNimZZdqxLz+1JPrcojZTIIOo7+vjDoXKtzfFprHYrvz6tJO3fP1N4b0ZCCBw/YG5aJKunKF6cnwsvjlt57sRzAGzO2kx2pI57vIyDh2c/jEkycaDqACfrT97YSWrPwPk3AQnWeFnuzdWExkHew8rrXf8BdtdUC1W2dPPaUcV7o/ey8LGwmIz84y2KF+eXu4vp6rNqbJHv8n7p+07vjTo0VzAcIXD8hH9y5OK8fqKKK03dGlvjm5ysP8m+qn0YJaOzaZ43kxaWxp05dwLceEXV7h8pzzM+CQnTXWSZhiz7RzCHKsLtwrsuOeVzu4oZsMksy45hSZb35d5czafmp5IRE0xTVz8v5ZdpbY5PYrVb+c3p3wDw5RlfFt6bURACx0+Ynx7FKocX5zlRUeUWVBFwZ/adZIRnaGyNa1C9OPk1+RyvOz6xN1efUESAZIA133aPgZ4mOBqWOMTr7qdu2otT0dzNX3zEe6MSYDQ4ixt+vadEeHHcwAelH1DeXk6kJZIv5H5Ba3N0ixA4fsQ/3aIsOn87XklFs/DiuJILzRfIr8nHKBl5ePbDWpvjMlJCU9gyeQtwA7k4+/9beZ75GYib6lrDtGTp18ESAfXn4dIHN3WqX+wuxmqXWZETy6LMaBcZqD13zklhUmwIbT0D/FX0xXEpNrtNeG/GiRA4fsSCjChWTo7Fapf55R4x/deV/P787wHYkLHB5yb4PjTrIYySkYKaAi42Xxzfm1rKoVBpdKj7kQwTJSgKFt6vvD70yxs+TXNXP687ZsU9ui7HBYbpB6NB4u+WZwLwwoFSbGLor8vYVraNsvYyIiwRwnszBkLg+Blfd5Rxvn68ktbufo2t8Q0auht4v/R9AL40/UsaW+N6kkOT2ZCxAYA/FP5hfG86/Bul703WGt/IvbmaxQ+DZISyfVBzY2X0fz58hT6rnVkpESye5DveG5VPL0glIiiA8qZudhTWaW2OTyDLsnMzde+0ewkJCNHYIn0jBI6fkTcpmulJ4fQO2HlFNONyCa9cfAWr3crcuLnMjputtTlu4Z5p9wBK342mnqbrH9zXAcdfVl4v8d5S+esSkQoztiivb8CL02+187IjAffvVmR6Zd+bsQg2m/jC4nQAnt9fqrE1vsGphlOcazqH2WDmc1M/p7U5ukcIHD9DkiTud7iOXz5YhlUMxrspeq29vHbxNcA3vTcqc+LmMCt2FgP2AV679Nr1Dz7xR+hrh5jJ3jVzaqIs+ZryfPav0DExD8UHZ2uoa+8jLszC7bOS3WCcPvjysgxMBomC0mbOVrVpbY7X88fCPwJKG4roQN/z+rkaIXD8kDvnJBMdYqa6rZft54Xr+GZ4p+QdWvtaSQlNYV36Oq3NcRuSJHHvtHsBePXCq/TbRglv2m1Q8Cvl9ZJ/AIMPLzGpCyF1Mdj64ejz436bLMtOj8aXlmRgNvnu7ygpIojbZycBwotzs9R21bK9fDsw6FEVXB/f/Z8lGJXAACNfdLiO/+9gmbbGeDF22c4fzis5KV/M/SImg0lji9zLhswNxAfF09TbxLayUTr5XtoGLaUQGAlz/CABcqnDi3PkeRjoHddbjl9p4XRlG2aTgS/mpbvROH3w4ApluvU7p6qpbRvf70hwLa9dfA2bbGNBwgKvHQHjaYTA8VPuXaK4jg+XNnOuWriOb4QDVQcoaSshJCCET03+lNbmuJ0AQwBfmKaIlj+c/wPySPOY8h2l5AsfALMfJEDm3gERadDdCGfGCN05eOFAGQBb5iYTG2pxo3H6YHZqJIszo7HaZWfekWBi9Fp7+culvwA4PamCsfGIwHnuuefIzMwkMDCQvLw8Dh8+POqxa9asQZKkax63336785j777//mu9v2rTJEz+Kz5AYEchtsxTX8YuOBVcwMdRqhk9N/hSh5lCNrfEMn5n8GSxGC4XNhRyrOzb8mzWnoHw/GEyw6CFtDPQ0RpNSUQVKsvEYQzirWnvYdrYWgAeWT3K3dbrh7xxenD8WXKG7XzT+mygflH5Aa18rSSFJrElbo7U5XoPbBc6rr77KY489xve+9z2OHz/OnDlz2LhxI/X19SMe//rrr1NTU+N8nD17FqPRyGc/+9lhx23atGnYcX/+85/d/aP4HPcvywTgrVPVNHX2aWuMl3Gp5RL5NfkYJINfxcMjAyOdQzivKRlXq4mmb4GIFM8apiXz71OGcNafh5Ld1z305fwybHZlLMO0pHDP2KcDNkxPID06mLaeAf52vEprc7wKWZadycVfyP2Cz4fCXYnbBc4zzzzDQw89xAMPPMD06dP51a9+RXBwMC+88MKIx0dHR5OYmOh8bN++neDg4GsEjsViGXZcVFSUu38Un2N+eiRzUiPot9r58+ErWpvjVaiVU7ek30JKqB/dzBl0ke+8spOKDkergc4GOPNX5bWal+IvBEXCPIfIVROsR6C738qfC5T/Z3/nR94bUBr/fdmxofrjofKRw5uCETlad5SLLRcJMgX5RSjclbhV4PT393Ps2DHWrx8sFTUYDKxfv578/PxxneP555/n7rvvJiRkeDx/9+7dxMfHM3XqVL761a/S1DR6b46+vj7a29uHPQTDS8Z/f6icAVEyPi56rD28X6I09vvslM+OcbTvkR2ZzdKkpcjIzrwATv0Z7AOQskB5+BtqmOryR9BePeIh756qob3XSkZMMOty4z1onD749PwUzCYDF2o7OFUp8v7Gy58K/wTAJ7I+QYQlQmNrvAu3CpzGxkZsNhsJCQnDvp6QkEBtbe2Y7z98+DBnz57lK1/5yrCvb9q0iZdffpkdO3bw4x//mD179nDbbbdhs9lGPM9TTz1FRESE85GWlnbjP5SPsXlWErGhFura+/hYlIyPi4/KPqJjoIOU0BTykvK0NkcT1CZjbxe9zYCtf7Cx3/wva2iVhsROhvRlSvfmk38c8ZBXHUM1716UjsHge439xiIy2MzmmYkAvHpEeIzHQ11XHTsrdgKiNPxG0HUV1fPPP8+sWbNYvHjxsK/ffffd3HnnncyaNYstW7bw7rvvcuTIEXbv3j3iebZu3UpbW5vzUVEhOviqWExGPrtQmZ302lHxexkPr19+HYBPT/40BknX/4XcxurU1UQHRtPU28S+k89D02UlD2WmH7vQ59+nPB///TVTxovqOzhW3oLRIPHpBf4V0hzK3Y72FG+frBZTxsfB28VvY5ftLEhYQHZkttbmeB1uXZ1jY2MxGo3U1Q33DNTV1ZGYmHjd93Z1dfHKK6/w4IMPjnmdrKwsYmNjKSoqGvH7FouF8PDwYQ/BIJ9doAicPZcaRJ+KMShpLeF4/XGMkpG7cu7S2hzNCDAGcGf2nQC8fsGR4D/zU2AJ09AqjZl+lzJlvLUcyvYO+9arjrEo63LjiQ8L1MI6XZA3KZpJsSF09dt49/TIoTyBgl2280bRGwB8MueTGlvjnbhV4JjNZhYsWMCOHTucX7Pb7ezYsYOlS5de971/+ctf6Ovr4957x675r6yspKmpiaSkpJu22R/JigtlcWY0dhn+5phuLBiZv13+GwArU1cSH+x/eRRD+eRkZdHdN9BMndHov+EpFXMwzHbkZKkhO5S5U687Koc+v9C/w+OSJPH5Rcrv4M+Hhcf4ehyrO0ZFRwUhASHOYbeCieF2//pjjz3Gb3/7W1566SUKCwv56le/SldXFw888AAA9913H1u3br3mfc8//zxbtmwhJiZm2Nc7Ozt5/PHHOXToEGVlZezYsYO77rqLnJwcNm7c6O4fx2cZGqay20WFw0j02/p5p/gdQOkH4+9kRWQxPygZuyTxdkKmMrrA31HDVIXvQJdS+LCjsI6mrn7iwyysmRqnoXH64NPzUzEZJE5WtHKhVhR8jMYblxXvzabMTQQHBGtsjXfidoHz+c9/nqeffponnniCuXPncvLkSbZt2+ZMPL5y5Qo1NTXD3nPx4kX2798/YnjKaDRy+vRp7rzzTqZMmcKDDz7IggUL2LdvHxaL73cFdRe3z04ixGykvKmbw2XNWpujS3ZW7KSlr4X4oHiWpyzX2hxd8Mn2VgDeCA3CjhDGJM1RHrZ+OP0qMJhc/JkFqZiM/pmzNZS4MAvrpynr/yvCizMiHf0dzrlTqqdUMHE88r/tkUceoby8nL6+PgoKCsjLG6w82b17Ny+++OKw46dOnYosy2zYcK1bLigoiA8//JD6+nr6+/spKyvjN7/5zTWVWoKJEWw2ccccZarxa0fEojMSf7ukhKe2TN4imm0BVJ/k1qqLhNjtVAy0c7T2qNYW6QNnsvHLVLd0s+dSAwCf8/Pw1FDuXqz8Lt44UUXvwMjVr/7MB6Uf0GvrJTsim9mxs7U2x2sR2wmBk885YuPvn62hvXdAY2v0RUVHBYdqDgEi4c/Jid8TLMtsNiubi9eLXtfYIJ0w8zNgCoKGQvbv3oYsw5KsaDJj/WA21zhZOTmOlMgg2noG+PDc2C1D/I03i94EFO+NJPlfSwFXIQSOwMm8tEhy4kPpHbDzzilR4TAUNR6+NGkpqWGpGlujA/q74bTS5O9TMxSPxfay7bT1iQZuBEXCjC3Ky7NKTxw1sVagYDRIzrw/EaYazuWWy5xpPINJMvGJrE9obY5XIwSOwIkkSc4qj9eOimoqFZvdxltFbwHw6Smf1tganVD4NvS1QWQ6M2Z9iSlRU+i39/N+6ftaW6YPHGGqddZ9JAQOcNtMUeF5NZ9bmIZBgvySJsqburQ2RzeopeGr01YTExQzxtGC6yEEjmAYn5yfgskgcaqilYu1HVqbowsO1x6mvqeeCEsEa9PWam2OPjjhGLQ570tIRqNzRo7aBNHvSV9KnTmNEKmPb6VdJDDAqLVFuiM5MojlObEAvHlCeIwBBmwDvFv8LoCYO+UChMARDCM21MIt05T+Lq+KZGMA3i1RFpxbM27FbDRrbI0OaK+Gsv3K6zl3A8qcHJPBxIXmC1xuuayhcfqgo8/Kn3qWAbDBtkdja/TLlrlKV+e3TlaJAZzA7srdtPS1EBcUx7LkZVqb4/UIgSO4BrXa4+1T1Vj9fABnr7WXHVeURpUiHu7g7N8AGdKXQqTSej/CEsHKlJUAIkwFfHSujr9ZlWamodUHoUMk0o7ExpmJBAYYKGns4kyVyN96u+htAO7MvlNUaroAIXAE17BqShxRwQE0dvaRXzL6lHZ/YHflbroGukgOSWZu/FytzdEHp19TnmcNn6R+e9btALxf8j522b+F8Zsnq6iU46kJm42E7BCFgqsJtZjYMF0Z2+PvYaq2vjb2Vyue0Tuy79DYGt9ACBzBNQQYDWyepSRFvnXSvxed90reA2Bz1ma/Haw5jIaLUHsaDCaYvmXYt1anriYkIITqrmpONZzSxj4d0NDRx4GiRgDM8z6vfFEVhYJr2DJX6b/l7x7j7eXbsdqtTImaIgZrugixYgtG5C5HbPzDs7V+24irtbeV/ZXKjkqEpxyoN+qc9RAyvMIj0BTILem3AIPC0B9593Q1dhnmpkUSk3c3SEaoOQmNIjdpJIZ6jA8W+6/HWA3tbp60WWNLfAchcAQjsjAjiqSIQDr6rOy+WK+1OZrwUflHWGUrudG5YkcFIMtwRul9c3V4SuX2SUqY6sOyDxmw+2ezSNXredfcZAiJhRxF9AkvzsgEGA3cPlvxGL95skpja7ShrqvO2Qn8tkm3aWyN7yAEjmBEDAaJO+cMuo79EdULod60/Z6Kw9BaDuZQmDryLnNx0mJiAmNo7WslvzrfwwZqT3lTFycrWjFIOG/azPqc8nzmNUUkCq5hyxCPcU+//3mMt5VtQ0Zmfvx8kkOTtTbHZxACRzAqdzpi4x8X1tPhZ6MbqjqrOF5/HAlJ7KhUVO9N7ifAPPJ0Y5PBxKZJmwD/DFOp3pvlObHEhwUqX5x6GwQEQ0sZVB3TzjgdsyAjitSoILr6bWwvrNPaHI+jhqfEWuNahMARjMr0pHCy40Lot9r58Jx/LTrvlygLzuLExSSEiEGu2AbgnKOJ3+yRw1MqqsdrV8Uuuge63W2ZbpBl2RliUXPYALCEQq7DCyjCVCMiSdJgT5wT/hWmKmsr43zTeYySkVszb9XaHJ9CCBzBqEiS5Fyo/SlMJcvyYHgqS4SnACjeBd1NEBIHk9Zc99CZsTNJC0ujx9rDropdHjFPD5yrbqekoQuLycDGGVeJYjVMde51sFk9b5wXsGWe4jHec6mB5q5+ja3xHB+UfgDA0uSlRAdGa2yNbyEEjuC6qHk4B4oaaezs09gaz3Cx5SLFbcWYDWbWZ6zX2hx9cMbheZj5aTBevwGZJEnOShB/ClO95fDerJ+WQFhgwPBvZq+F4BjoaoCS3Z43zgvIiQ9jRnI4VrvMe2dqtDbHI8iyLKqn3IgQOILrkhkbwpzUCGx2mff9ZNFRw1Or01YTZg7T2Bod0N8FFxxCZZTqqavZnKUs1gerD9Lc2+wuy3SDzS47vZx3zR0hSdQYADMcs4XOiDDVaKhhqrf9pJrqfPN5ytrLsBgtrEtfp7U5PocQOIIxudM5L8b3w1SyLPNR+UeASPhzculDGOiGqEmQsmBcb8mKyGJa9DRsso2Pyj5ys4HaU1DaRF17H+GBJlZPjRv5IFUcXngfBno9Z5wXoVaeHS1vob7d939HH5Qo4ak1aWsICQjR2BrfQwgcwZh8YnYSkgTHyluoaPbtpNHzzeep6qwiyBTEipQVWpujD86/pTzP+CRI0rjfprrct5dvd4dVuuKDM8qsqU0zE7GYRpkcnroIwlOgvwNK/Cc3aSIkRwYxLz0SWYZt53x7fpfNbnPm34jwlHsQAkcwJgnhgSyZpHSt9fUw1fYy5Wa8MmUlQaYgja3RAf3dcNnhgZl+14TeuiFzAwBH647S1OO7HWptdtl5M1ZHnIyIwQDT7lReq6JRcA2bZyq/Q19fa042nKS+p56wgDCxmXITQuAIxsVmh+v4g7O+u6uSZdnpbVBvzn5P0cdKeCoyA5LmTOitKaEpTI+Zjl22s7Nip5sM1J5j5S00dPQRFmhiWXbs9Q9WReKF98HqP5VCE2HTTGX45uHSZho6fLew4ePyjwFYm74Ws9GssTW+iRA4gnGxcUYCkgQnK1qpaevR2hy3cLHlIlc6rmAxWliVskprc/SB6mmYfteEwlMqGzIUoagu5r7IB2cVT8OG6QmYTWMsqWl5EJoIfW1QuscD1nkfadHBzE6NwC7Dhz4appJlmY+vKP8n1qeLSk13IQSOYFzEhwWyMCMKgG0+6sVRk2FXpqwkOGDkTr1+xUAvXNqmvJ5geEpFFTgFNQW09ra6yDD9YLfLzv8Pt828TnhKxWCAaXcor8+/6T7DvBw11KeKR1/jbONZartqCTYFsyxlmdbm+CxC4AjGzaaZvhumGhaeyhDhKQCKd0J/J4Snjrt66moywjOYGjUVm2zzyaZ/pypbqWnrJcRsZOXkMcJTKs4w1XtKh2jBNdzmCFMdKmmmyQf7b22/oqw1q1JXYTFaNLbGdxECRzBu1Nj4kTLfi41fbr1MWXsZZoOZ1WmrtTZHHzjDU3feUHhKRW2W6IvVVKrYv2VaAoEBo1RPXU3GMgiOhZ4WKNvnRuu8l4yYEGYkh2Ozy2w/71tjYmRZdhYziEai7kUIHMG4SYkMYk5qBLIMH533LS+OevNdnrJc9KMAsPbBRaWE9UbDUyq3ZijzdfJr8uno77hZy3SDLA82v1Q9DuPCYBwSphLVVKOhhql8ravxxZaLVHZWYjFaWJmyUmtzfBohcAQTwhmmOuNjAqdMhKeGUbJHSYQNTYTUxTd1qqzILLIjsrHareyu2O0S8/TAuep2Klt6CAowsmZq/MTerIrGwnfFbKpRUEXjweImWnxoNpVzM5W8XOT6uRkhcAQTQl108kt8Z9Epbi2muK0Yk8HEmrQ1WpujD4aGpww3v0yoZfe+FKZSvTdrc+MIMo8zPKWSuQKCoqG7Ea4cdIN13k9WXCi5iWFKmKrQd8JUakWhCE+5HyFwBBMiMzbE5xYddTTDsuRlYvYUKImvF95VXt9keEpFLYU9UHWAroEul5xTS2RZdubfbBpP9dTVGAMg1zGpXoSpRsVZTeUjYaqS1hJK2kowGUwi188DCIEjmDBqOayvlIurXgU1V8TvKd0Lva0QEgfpS11yyilRU8gIz6Df3s++Su9PrL1Y10FpYxdmk4F1uRMMT6lM36I8F74DdpvLbPMlNs9SPMb7ixpp6/H+ijN1rVmatJRwc7jG1vg+QuAIJsxt6qJzuZGOXu9edMrby7ncchmTJMJTTlSPwrQ7lIRYFyBJkjO/SfWYeTPvO3LQVk+JI9RiurGTTFoFgRHQWQdXDrnQOt8hJz6MyfGhDNhkdl7wfo+x2txP5Pp5BiFwBBNmcnwoWXEh9Nvs7LxQr7U5N8WuK0pvloWJC4mwRGhsjQ6w25T+LDA4N8lFqIv6/qr99Fi9uxv2trM3UD11NSYzTHUMWVR/54JruHVGAoDXl4tXtFdwofkCRsnI2rS1WpvjF3hE4Dz33HNkZmYSGBhIXl4ehw8fHvXYF198EUmShj0CAwOHHSPLMk888QRJSUkEBQWxfv16Ll++7O4fQ+BAkiTnwu7t1VRq8zmx4DioPKokvloilERYFzItehrJIcn0WHs4VO29Hovypi4u1XViMkjckptwcydTBc7F90CWb944H2TDdGWt2XOxgT6r94by1OZ+CxMXEhkYqa0xfoLbBc6rr77KY489xve+9z2OHz/OnDlz2LhxI/X1o+/8w8PDqampcT7Ky8uHff8nP/kJ//u//8uvfvUrCgoKCAkJYePGjfT29rr7xxE4UPNw9lxqoHfAOxed5t5mTjacBITAcXLR4UmYvEFJhHUhkiSxNl35PXtzV2PVk5CXFU1E8E3+jrLXgdECLWVQX3jzxvkgs1MiiA+z0NVv42Cx906l33lFGTi7IV2EpzyF2wXOM888w0MPPcQDDzzA9OnT+dWvfkVwcDAvvPDCqO+RJInExETnIyFhcJckyzI/+9nP+M53vsNdd93F7Nmzefnll6murubNN990948jcDAjOZykiEB6BmwcLG7U2pwbYk/FHuyynWnR00gKvYFKGF/kwvvKc+5mt5xezXPaU7kHm5cm1qoCZ/20m/TeAFhCIWuN8vqiCFONhMEgsX66d4epGnsaOd1wGkDk+nkQtwqc/v5+jh07xvr1g/X+BoOB9evXk5+fP+r7Ojs7ycjIIC0tjbvuuotz5845v1daWkptbe2wc0ZERJCXlzfqOfv6+mhvbx/2ENwckiQ5F/jt570zD2dnhbKjEt4bB42XoekyGAIgxz27zAUJCwgLCKO5t5kzjWfccg130tLVz5GyZsBFAgcGxaQqLgXXsMEhcD4+X4fd7n2hvL2Ve5GRmREzg4QQF31uBGPiVoHT2NiIzWYb5oEBSEhIoLZ25NyNqVOn8sILL/DWW2/xhz/8AbvdzrJly6isrARwvm8i53zqqaeIiIhwPtLS0m72RxMwuOjsKPS+RWdoHogaNvF7LjpusJNWQqB7SlgDDAGsSFVye7wxTLXrYj12GXITw0iLdlEX2imblOfq49DuG/1eXM2y7BhCzEbqO/o4XdWmtTkTRv2sC++NZ9FdFdXSpUu57777mDt3LqtXr+b1118nLi6OX//61zd8zq1bt9LW1uZ8VFRUuNBi/yUvK5pQi8krF5386nx6bb0khyQzNWqq1uboA9WDMNU94SmVdWnrAO8UOGqI5NbpLtyFhyVCykLl9aUPXHdeH8JiMrJ6ahwA271sDt6wzZTwFnsUtwqc2NhYjEYjdXXD46Z1dXUkJo6vvDIgIIB58+ZRVFQE4HzfRM5psVgIDw8f9hDcPBaTkdVTlEXnYy+LjQ/dUUk3MSnbZ+hsgIoC5fXU29x6qeUpyzEZTJS2lVLWVubWa7mS3gEbey41ADhzQlyGCFONyQYvzcMpqClwbqamRE3R2hy/wq0Cx2w2s2DBAnbs2OH8mt1uZ8eOHSxdOr4OqTabjTNnzpCUpCSBTpo0icTExGHnbG9vp6CgYNznFLiO9dOVLq4fe9HYBpvdxt7KvYAITzm5tA2QIWkORKS69VJh5jAWJSwC8Krhm/klTXT320gItzArxcU9k6Y6xjaU7oE+35m47krWTo3HaJC4VNdJeZP3jPtQP+NiM+V53B6ieuyxx/jtb3/LSy+9RGFhIV/96lfp6urigQceAOC+++5j69atzuOffPJJPvroI0pKSjh+/Dj33nsv5eXlfOUrXwGU5NZvfOMb/PCHP+Ttt9/mzJkz3HfffSQnJ7NlyxZ3/ziCq1AXnQu1HVQ0d2ttzrg41XCK5t5mwsxhLEhYoLU5+uCiIzSi3mjdjJqL4E1hqqHVUy6/UcVNhahJYOuH4p2uPbePEBlsZnFmNOA9Xhy7bB8mcASexe0C5/Of/zxPP/00TzzxBHPnzuXkyZNs27bNmSR85coVamoGE+taWlp46KGHmDZtGps3b6a9vZ2DBw8yffp05zHf/OY3efTRR3n44YdZtGgRnZ2dbNu27ZqGgAL3ExlsZmFGFOA9Xhz1proyZSUBBtf2evFK+rsHb6puKg+/GjUX4WTDSZp7mz1yzZvBbpfZ4fh8b3B1eApAkgaHb4ow1aiov/uPvETgnGk8Q1NvE6EBoSxMWKi1OX6HR5KMH3nkEcrLy+nr66OgoIC8vDzn93bv3s2LL77o/Pd///d/O4+tra3lvffeY968ecPOJ0kSTz75JLW1tfT29vLxxx8zZYqIbWqFs4TTCwSOLMuD3YtFeEqhZDdYeyAiHRJmeuSSSaFJTIuehl22O8OFeuZMVRt17X2EmI0szY5xz0XU5O7LH4LN6p5reDnqWnO0rJnmrn6NrRkb1XuzImUFAS5unCkYG91VUQm8j1sc/UAKSpp1P/G3tK2U8vZypVw52bWjCLwWtcFc7mbFk+AhnGGqK/oPU6niffXUOCwm1wwgvYa0PAiKhp4WuDJ6nzB/Ji06mNzEMOwyXjEHT4SntEUIHMFNMyk2hJz4UKx22VllolfU5n6LkxYTag7V2BodYLfBxW3KazdXT12NGqbKr8mn16rvMStqzodbwlMqRhNM2ai8vijCVKNxq7OaSt/l4hXtFRS1FmGSTKxIEZspLRACR+AShnYa1TN7KvYAsDZVhKeAweGagRGQsdyjl86NziUxJJEeaw8FNQUevfZEqGju5kJtB0aDxNqp8e692NDp4mL45oiowzf3XW7U9fDN3ZW7AaV7d4TFxVV3gnEhBI7AJaht63ddrGfAZtfYmpFp7W3ldKMyD2ZV6iqNrdEJlxzemxzXD9ccC0mSWJO6BtB3NZXqvVmUGUVksNm9F8teB0YztJZD4yX3XstLmZEcTlyYhe5+GwUl+k1QF92LtUcIHIFLmJsWSWyomY5eq24Xnf3V+7HLdiZHTRbDNVUuf6Q8q6ERD6Mu/vuq9iHr1GOx66KS63FLrgdmCFlCIdMRzrj0ofuv54UYDBLrHJ40vebhtPW1cbzuOACr01ZrbI3/IgSOwCUMdd/rddFRq3VWpQjvDQBtlVB3FiQD5Kwf+3g3sDBxIUGmIOq767nYclETG65HV9+gYF+b6+bwlMpkh9hUxafgGtS/xa6L9boUxgeqDmCTbeRE5pAWJmYfaoUQOAKXsc6x6Oy+qD+BY7VbOVB1ABA7KifqDTR1EQRHa2KCxWghL1FpG6HHcvEDRY302+ykRweTHRfimYtOuVV5vpIPvd41481TrJgcS4BRorypm5JG/XU13lvl2EyJULimCIEjcBkrJsdiMkiUNHZRprNF51TDKdr724mwRDA7drbW5uiDSw6BM/lWTc1YmboS0KfAUcNT63LjPddmPzoLYiaD3Sq6Go9CqMVE3iSlH9EunXmMbXabczMlBI62CIEjcBlhgQEscrRS36UzL45681yevByjwU19TLyJgV5l7hFoln+jot4ETjecpqW3RVNbhiLLMrsuKG0P1jgmWXsM9W9ySYSpRkMNU+ktJH6m8Qytfa2EmcOYEzdHa3P8GiFwBC5lnU4XHWf+jdhRKZTth4FuCE/xWPfi0UgMSWRK1BRkZA5UH9DUlqGcr2mntr2XoAAjS7Lc1L14NFSvWtF2sOuzKlFr1LXmcGkzHb36aTC6r2ofoGymTAaTxtb4N0LgCFzK2lxlp1tQ0kx3vz7azdd01lDUWoRBMoiGWyqXHRU6kzd4tHvxaKjCU09hqt0XFe/N8pwYAgM87PVLXwrmMOhqgOoTnr22lzApNoRJsSFY7TL7LzdqbY6TfZWKwFFDrwLtEAJH4FKy40JJiw6i32bnQFGT1uYAgzfNuXFzRcMtUBrIqSXIk7UNT6moAudA1QGsdn0IY9UL6bHqqaGYzJDtaEZ5WZSLj4beKjfru+spbC5EQmJ5smcbZwquRQgcgUuRpMFycb3k4agVDWJH5aDxktJIzmiBLH1UlM2OnU2EJYL2/nZON5zW2hxauvo5cUXJB1rj7u7Fo+HMwxECZzTWOcvFG7DbtS8X31+1H4CZsTOJCfJwWFNwDULgCFyOU+Bc0L5HxdAxACL/xoF6w8xcAWYPlT6PgdFgdO549RCm2nOpAbsMuYlhpEQGaWNEzgblueYkdOh77pJWLJ4UTYjZSGNnH2eqtC+pF+EpfSEEjsDlLM2OwWIyUNPWy8W6Dk1tOVJ7hD5bH4khiUyOnKypLbrhsj7Kw6/GmYdTpb3AUb2PmoSnVMISIGmu8rroY+3s0DFmk4EVk2MB7cNUA7YB8muUKfCimag+EAJH4HICA4wsy1Z7VGg7XXxo92KP9THRM71tSgM5GGwopxOWJy/HIBm43HKZ2i7tPBY2u8yeS8rndp2WAgdEmGocqCM0tA6JH68/TtdAFzGBMUyLmaapLQIFIXAEbsEZG9dwVyXLslPgiO7FDop3KQ3kYiYrDeV0RGRgpLMJo5ZhqhNXWmjtHiAiKIB5aZGa2QEMJoEX7wJrv7a26JQ1jsrN05Vt1Hf0amaH+pldkbICgyRurXpA/BUEbkFNzDx2pYW2bm16VBS1FlHTVYPFaGFR4iJNbNAdGg/XHAs1TKXmMmiBGupYNSUOk1HjJTJ5HoTEQX/HoOdNMIz4sEBmpSjVkWppvxao/W9E/o1+EAJH4BbSooOZHB+KzS6z97I2i47aLl0d6Oj32O1webvyWmf5NyqqwCmoLaDP1qeJDbsuquEpD3cvHgmDYTDZWAzfHJW1jk7TamjR01R0VFDaVopRMrIseZkmNgiuRQgcgdsYOvFXC/ZXKyWbK1PEjgqAujPQVQ/mUKWRnA6ZEjWFhOAEeqw9HK096vHr17b1UljTjiTB6ika59+oTHYInKId2tqhY1Y7BM7+y41YbZ7v/Kx6HOfFzyPMHObx6wtGRggcgdtQ5/fs0aBHRfdAN8frjgOI7sUqqvdm0mqlkZwOkSTJ+ffSYmzDXocHYHZqJNEhOvkdZa0ByQANhdBWqbU1umROaiThgSbaegY4Ven5cnExPVyfCIEjcBuLMqMJtZho6urnbLVnF53DtYcZsA+QGppKeli6R6+tW1QPQM4t2toxBstTlH44aojRk6ghjtVTdBCeUgmOhpSFymtRLj4iJqOBlVO0CVP1Wnud3kaxmdIXQuAI3EaA0eAsF9/r4UVH7Si6PGW5KA8HpTy8Qml4qHeBk5eUh1EyUtJWQnVntceua7XZ2XdZhwIHIGe98iwEzqis1kjgHKs7Rp+tj4TgBHIiczx6bcH1EQJH4FZWa5D8J8uyU+CIHZWDkj0g25Ty8KhMra25LuHmcObEzQEGhaonOFXZSnuvlYigAOak6mxm2WSHwCnZAzb9TM7WE6rAOV3ZSlOn5xLUh641YjOlL4TAEbiVVZOVRef4lVbaejyzMJe3l1PVWUWAIYDFiYs9ck3do+78VU+AztEiTLXHUT21YnKs9uXhV5M0D4JjoK8dKo9obY0uSQgPJDcxDFmG/UWemy6u5oqpn1mBftDZ/2KBr5EWHUxWXAg2u8xBDy066oIzP2E+wQHBHrmmrpHlQYEz2bsETkFtAQMe8ljoMv9GxWCA7HXKazVZXHANav+tPR7qh1PVWeUsD89LyvPINQXjRwgcgdtRbxie6oejNtxakSzCUwA0XID2KjAFQoZ37DKnRU8jOjCaroEuTjacdPv1mjr7OO0Y1qhLgQMiD2ccDF1rPFG5qXoYZ8fNJtwc7vbrCSaGEDgCt+NM/rvY4Pbp4qKiYQTUG2LmCgjwjoaHBsngbJjmiTDV/qJGZMf08ITwQLdf74ZQPTi1p6GjTltbdMqCjCjHdPF+zte0u/16B6sPAsocNYH+EAJH4HaWZCnTxavbeimq73TrtYZWNGRHZrv1Wl6Dl+XfqDjzcDzQD8cZnpqqU+8NQGj84HTx4p2amqJXzCYDy3OU6eK73dxgdMA+wKGaQ4DYTOkVIXAEbicwwEhellIu7u5qKlHRcBX9XVCu7DK9TeAsS16GhMSF5gs0dLvvc2O3y+y9pOSH6TY8pSLCVGPiqcrNU/Wn6BroIsoSJaaH6xQhcAQeYdVkZVflSYEjAEr3ga0fItMhxrt6dEQHRjM9ZjowGApwB+dr2mns7CPYbGRhRrTbruMSVIFTvAPsNm1t0SmeqtxUPYvLUpaJ6eE6xSN/leeee47MzEwCAwPJy8vj8OHDox7729/+lpUrVxIVFUVUVBTr16+/5vj7778fSZKGPTZt2uTuH0NwE6hjGwpKm+npd8/CXNlRSVl7mahoGMrQ8JQXerTUMJU7++GoontZdixmk85vVKmLwBIBPS1QfUJra3RJWnQw2R6o3FRzw0T+jX5x+//mV199lccee4zvfe97HD9+nDlz5rBx40bq60eOj+7evZsvfOEL7Nq1i/z8fNLS0rj11lupqqoadtymTZuoqalxPv785z+7+0cR3ATZcaGkRAbRb7VzqLTJLddQF5w5cXPEwDsVp8DZoK0dN4g6KPVg9UFsbvJYeEX+jYrRBFmrldciTDUqznJxN3mMG3saKWwuBBDTw3WM2wXOM888w0MPPcQDDzzA9OnT+dWvfkVwcDAvvPDCiMf/8Y9/5Gtf+xpz584lNzeX3/3ud9jtdnbsGD5J12KxkJiY6HxERUW5+0cR3ASSJLFqSDWVO1Cnh4vwlIOmYmgpBUMATPLOieozY2cSZg6jvb+ds01nXX7+9t4Bjpe3ALB6shcIHBgyXVwInNFQc6l2u6lyUw2ZTo+ZTkxQjMvPL3ANbhU4/f39HDt2jPXrB5MbDQYD69evJz8/f1zn6O7uZmBggOjo4bHx3bt3Ex8fz9SpU/nqV79KU9PoXoG+vj7a29uHPQSeZ/UUJQ/HHXOpBmwDHK5RQpnLUsSOChgcrpm+BCze6dEyGUwsTVoKuKdc/GBRE1a7TFZsCOkxXtIUMtsxS6zyKHQ3a2uLTlk8KRqLyUBtu3sqN52z7kR4Ste4VeA0NjZis9lISEgY9vWEhARqa2vHdY5vfetbJCcnDxNJmzZt4uWXX2bHjh38+Mc/Zs+ePdx2223YbCO7sJ966ikiIiKcj7S0tBv/oQQ3zLKcWIwGiZLGLiqau1167lMNp+i2dhMdGM20aFHRAAyWEut8uOZYqB45dwgcNYSxSu/VU0OJSIG4aYAMpXu0tkaXuLNy02a3kV+tbNCFt1jf6Dqj7kc/+hGvvPIKb7zxBoGBg8237r77bu68805mzZrFli1bePfddzly5Ai7d+8e8Txbt26lra3N+aioqHCf0f1d0OW5OSjeRHhgAAvSlVDibhcvOqrLeEnSElHRAGDthzKlo7OzQZyXsjRZ8eCcbTpLW1+by84ry7LTm6j78vCrUf+mRTuuf5wfo1Zu7rvs2vW4sLmQ1r5WQgNCmRU3y6Xn9lpaK8Bu19qKa3DrnSA2Nhaj0Uhd3fCum3V1dSQmJl73vU8//TQ/+tGP+Oijj5g9e/Z1j83KyiI2NpaioqIRv2+xWAgPDx/2cAv5v4AfZ8KeH7vn/D7AKkeYap+bBI5I+HNQeRj6OyEkDhK8exFODEkkOyIbu2ynoKbAZecta+qmqrWHAKNEXpbOy8OvRhU4xbuUWWOCa1g5Wa3cbKJ3wHUJ6qonMS8pjwBDgMvO67XYbfCr5fBfU5S8Px3hVoFjNptZsGDBsARhNWF46dKlo77vJz/5CT/4wQ/Ytm0bCxcuHPM6lZWVNDU1kZSU5BK7b5jIdKXniOgyOipqKCC/uIkBm2sUf0tvC+ebzgODu32/R/0MZq1VBjV6Oerf1ZX9cPY5ZqMtzIgm2Gxy2Xk9QsYyMFqgvRIaL2ttjS6ZkhBKQriF3gE7R8taXHZesZm6iuoT0NumeI0jM7S2ZhhuX/kee+wxfvvb3/LSSy9RWFjIV7/6Vbq6unjggQcAuO+++9i6davz+B//+Md897vf5YUXXiAzM5Pa2lpqa2vp7FQSxTo7O3n88cc5dOgQZWVl7Nixg7vuuoucnBw2btzo7h/n+kxaCZIRmoqgpVxbW3TKjOQIooID6Oizcqqi1SXnLKgpQEZmctRk4oPjXXJOr0cVOF4enlJRbyb51fkuq4pRw1MrHV5Fr8IcDBkOMS82VCMiSZLTi7PPRYN+O/s7Od1wGhACx4lzM7VKaWOgI9wucD7/+c/z9NNP88QTTzB37lxOnjzJtm3bnInHV65coaamxnn8L3/5S/r7+/nMZz5DUlKS8/H0008DYDQaOX36NHfeeSdTpkzhwQcfZMGCBezbtw+LxeLuH+f6BEYojbgASnZpa4tOMRok56yYvS6KjTs7iiaJBQeAriaoPqm8zl6rqSmuYkHCAgIMAVR3VVPWXnbT5+u32skvViovV3lLefjVOMNUIg9nNFZOdu1ac6T2CFbZSnpYOqlhqS45p9ej482UR+TWI488wiOPPDLi965ODC4rK7vuuYKCgvjwww9dZJkbyF4HFYeUP/qC+7W2RpesmhzHu6dr2He5gcc2TLmpc8myLFzGV1OyC5AhYSaEXT/XzVsIDghmfsJ8CmoKOFh9kEkRk27qfCeutNDVbyMmxMz0JDfl5Lmb7HWw/Qko2w/WPjBpvMHTISsnxyFJUFjTTn1HL/FhNzcpXt1MiVC4g952qHBMGtChwPH+4LzeUP/IJbvFrJhRWOHYVZ2qaKWt++ZmxRS3FlPfXY/FaGF+wnxXmOf9FDu8hz7ivVEZGqa6WdTKmhWTYzEYvG+EBaAI2JB4GOiGCtclX/sS0SFmZiZHALDfBV4c9bMn+t84KNsHsg2isyAqU2trrkEIHFeTPE8JVfW2iVkxo5AcGUROfCh2GQ4W39yio3pvFiQsINB0c7szn0CWde0yvhlUgXO49jADtpsTxmpOxkpvDU+BMlvMGaYSeTijsdJF5eIVHRVc6biCSTKxKHGRK0zzfpybKX2uNULguBqjCSY5ZsWIRWdUXBUbP1gjwlPDaLgAHdVgCoR03/qdTImaQnRgND3WHk42nLzh87R09XO6Sumno34OvRbRD2dMhiYa2+03nqCuem9mx80m1BzqEtu8Hp1vpoTAcQdiVzUmamLn3ks3Piumz9bHsdpjgBA4TtTPXMZyCPAtj5ZBMjj/zjdTLn6guBFZhtzEMBLCvfx3pIYha09Dp3tmvHk7CzKiCDYbaezsp7D2xsf0qAJHrDUOWsqguVipHM7U56w7IXDcgbroVBxWkrAE15CXFU2AUaKqtYeyphsb23C87ji9tl7iguLIicxxsYVeis53VDeLKwSOszzc2703AKHxkOho5FiyW1NT9IrZZGCpY2zD3ks35jG22q3OJpNC4DhQw1NpiyFQn4n6QuC4g6hMiM5Wkq/UdvmCYQSbTSzMULrH3miPCnVHtTR5KZLkpYmirmSgF8oc85q8fP7UaCxJWgJAYVMhzb0THzQpy7IzF8Or82+GIjzGYzKYh3Nja83ZxrN0DHQQbg5nesx0V5rmvXjBZkoIHHchFp0xWemcLn5juyq1ZFNUNDi4kg/WHghLgrhcra1xC3HBcUyJmoKMfENjG4obOqlp68VsMrB4kpeNZxiNoWuNGNswImoH9aNlLXT3Wyf8fnUztSRpCUaD0aW2eSU26+CgVyFw/BAhcMZEzcPJL26c8NiGhu4GLrVcAmBJ8hKX2+aVDN1R+bBHSxW0NxKmUsV03qRoAgN85EaVtgRMQdBZC/XntbZGl0yKDSElMoh+m52Ckol7/kSvratQxzMERiiVwzpFCBx3kbkCDCZoLoHmUq2t0SXTk8KJDjHT1W/jePnEZsUcqjkEwLToaUQH+shO/GbRecmmq3DOpao6OOEE9cHycB/Iv1EJCFTWGxAbqlGQJMk56HfvBMNU7f3tnGk8A4gGf06c4xnWgI49WkLguIvAcEhdrLwWYxtGxGCQWJFzYz0qREXDVXTWQ90ZQFIWHR9mfsJ8LEYL9T31FLUWjft9fVYbhxy7dzVk4TMIj/GYrMhR/uYTbfh3pOYINtlGZngmyaHJ7jDN+/CC/BsQAse9iEVnTG4k+U+WZfJrBhOMBQxW0CTNhhAf8k6MgMVoYWHCQmBiXY2PlbfQM2AjLszC1IQwd5mnDWrlZvlBJdlccA3LsmOQJLhc30lt2/h/RyI8dRW9bVB5RHmdpe9u6ULguBPn2Ia9SlKW4BrUSpbTVW20dPWP6z2XWy/T2NNIoDGQefH6jf96FKfLWN8LjqtQha0qdMeDunNfkRPre1V3cblKcrm1V0k2F1xDVIiZ2SmOsQ1F4/fiCIFzFaWO8QwxORCVobU110UIHHeSPBcCI6GvDaqPa22NLkmMCGRKQiiyDAcd053HQt21L0hYgNlodqd53oEs++z8qdFQy8WP1R2j3zY+Yaze1Hwq/0ZFGhKaFCHxUVHn4O0fp8e4or2Cys5KTJKJhYkL3Wma9+BFmykhcNyJwQhZ6tgGseiMhjM2XjS+RUeEp66i4YJSQWMKVCpq/IApUVOICYyhx9rDqYZTYx7f0tXPGcd4BjXvy+dQbzhirRmVwbWmcVxjG9S1ZnbcbEICQtxqm9dQ4j3FDELguBt10RG7qlEZOgxvrKqYflu/czyDEDgO1BtaxjKfG88wGpIkOdsDjCcPRx3PMDUhjHhvH88wGqoHp/Y0dN385GxfZH5GJEEBytiGC7UdYx4/tJmoAGgpVyqDJeNg5Z6OEQLH3aghg8oj0Df2fyh/RB3bUNnSQ/kYYxtO1J+g19ZLbFAskyMne8hCnaOKZy9wGbuSpUmOPJxxCBxn/o0vhqdUwhIgYabyWoxtGBGLyUheltJWYiyPsc1uo6BWaSYpBI4Dda1JXaTb8QxDEQLH3URlQtQksFuhbL/W1uiSYLOJ+elRAOwbI/lvaEdRn0sUvRGs/YPjGfwk/0ZFvemcazpHW1/bqMcNHc/g0wIHRB7OOBicLn79teZc0zk6+jsIM4cxI2aGJ0zTP16W6ycEjifIFrHxsVg5zuQ/NSYuKhocVB6GgS4IiYN4/1qE44PjyYnMQUZ2Nn4cibKmbqpaezAbDeT5yniG0XCuNbvF2IZRUNeaw6XN9A7YRj1O3UzlJeZhMpg8YpuusdsGxzN4ibdYCBxPIPJwxmSFY1d1sLgJ6yhjG1p6WyhsKgQGq2j8nmEdRf3vv7P6ObhemEoVzfMzIgk2+/iNKn0ZGM3QXgmNl7W2RpdMjg8lIdxCn9XOset0UBfFDFdRcwp6WsASDikLtLZmXPjfiqgFk1aBZIDGS9BWqbU1umRWSgQRQQF09Fo5XTVyuKGgpgAZmZzIHOKCfawT7Y3iJ+MZRkO9+RyqOTRqgrrPTQ+/HuZgSHeIf7GhGhFJkliec/2xDd0D3c7qPDXXy+9RP0+ZK8HoHRsFIXA8QVAkJM9XXosw1YgYDRLLsmOA0Vupix3VVXQ3K0PvwOfHM4zGwoSFmAwmqjqrqOiouOb7VpudfEd/JZ8tD78aUS4+JoMh8ZHXmqN1R7HaraSEppAWnuZJ0/SLl+XfgBA4niNbhKnGYsV1Fh1ZlsX8qasp3QvIShfbcP+ckRMcEMzcuLnAyGGqU5WtdPRZiQwOYKaji63Po641ZfvBNqCtLTpF9eCcq26nqbPvmu+L8vCr6O+CK448Ny/JvwEhcDyHMw9nN9hHzjHxd1Y6mnAdv9JCZ9/w0Rbl7eXUdNUQYAhgQYJ3xH/djp+Wh1/N9cY2qOGp5dmxGA1+UnWXOAeCoqG/AyqPam2NLokPCyQ3UZlHdmCEDurqeAYRnnJQfhDsAxCRBjHZWlszboTA8RSpiyAgBLqbHFOfBVeTHhNMenQwVrtMQcnwRUddcObFzyPIFKSFefrDC13G7kC9CR2uOYzVPlwY+0X/m6sxGES5+DgYrXKztquWkrYSJCTykvK0ME1/qGtN1hplLIiXIASOpzCZBzs/itj4qKwY0tV4KCL/5iqaS6C1HAwBkLFca2s0ZXrMdMLN4XQMdHCu6Zzz6x29A5yoaAX8KP9GRbSmGBO1cnP/VR3U1ZYDM2JmEGHxk7DmWJR452ZKCBxPIvJwxmSl40Y0dNrvgH2AI7VHAOEydqKWh6ctBkuotrZojNFgdO60VU8fwKGSZmx2mcyYYNKig7UyTxvUsGXVUehp1dQUvbI4MxqzyUB1Wy8ljV3Or4v8m6voqIX684AEk9ZobMzEEALHk6iLTnk+DPRoa4tOWZYdi0GCovpOatqU39HZxrN0DXQRYYkgNzpXYwt1gtqK38/zb1Sc5eLVgw3/1NCDX4WnVCLTICYHZLvooD4KQWYjCzOUDupqKNMu250eHCFwHKhrTdJsCInR1JSJIgSOJ4mbCmFJYOuDK2PPz/FHIoIDmJUaCQwuOupNKy8xD6PBqJVp+sFuc1RQ4XUuY3ehevZON5yma0DZjatjP9QJ0n6HaDA6JleHxC+3XKa5t5kgUxBz4uZoaZp+KPbeYgYhcDyJJIkeFeNADVOpi47Iv7mK6hPQ2waBEZA8T2trdEFqWCqpoalYZStHa49S3dpDSUMXBgmWZnvXrtNliDycMVFzsw6VKB3U1fDUgoQFmI1mLU3TB7I86MHxws2UEDieRuThjIm6qzpY3Eh7bwenG04DQuA4UW9Yk1aB8Gg5GVournr/5qRFEhEUoKVZ2pG5AiQjNBdD6xWtrdElM5IjiAwOoLPPyqnK1sHNlMj1U6gvhM5aMAVCmveNxxECx9Oo5Zu1Z6Dr+tNs/ZX56VEEm400dvbz5oW92GQb6WHppISmaG2aPnDm36zR0grd4RQ41fnOJHW/q54aSmDE4Mwg9TMjGIbRILE8W/mM7LpUzfG644DYTDlRPzcZyyAgUFNTbgSPCJznnnuOzMxMAgMDycvL4/Dhw9c9/i9/+Qu5ubkEBgYya9Ys3n///WHfl2WZJ554gqSkJIKCgli/fj2XL3vJYLnQeEiYqbwWi86ImE2DU58/dOSaiOGaDvo6oaJAee2FMXF3sjhxMQbJQElbCftLiwA/FzggwlTjQPUY7yg5TK+tl7igOHIiczS2Sid4eTNRtwucV199lccee4zvfe97HD9+nDlz5rBx40bq6+tHPP7gwYN84Qtf4MEHH+TEiRNs2bKFLVu2cPbsWecxP/nJT/jf//1ffvWrX1FQUEBISAgbN26kt7fX3T+OaxBNuMZEbaV+sU2ZtSR2VA7UjqKR6RCdpbU1uiLCEsGMmBkAtEvnCTYbmZcepbFVGqPemEr3iA7qo6CK4NIuZa1ZkrQEyYua2bkNaz+UHVBee2H+DXhA4DzzzDM89NBDPPDAA0yfPp1f/epXBAcH88ILL4x4/P/8z/+wadMmHn/8caZNm8YPfvAD5s+fz89//nNA8d787Gc/4zvf+Q533XUXs2fP5uWXX6a6upo333zT3T+Oa3DuqnYrSVyCa1g5OQ7J1EafVINBMrA4abHWJumDoTsqsQhfg+rpM4UUsSQrBrPJz6PwqQvBHKZ0UK89rbU1uiQtOpiMmGAMwUoUQGymHFQehoEuCImD+BlaW3NDuPV/f39/P8eOHWP9+vWDFzQYWL9+Pfn5I5dJ5+fnDzseYOPGjc7jS0tLqa2tHXZMREQEeXl5o56zr6+P9vb2YQ9NSV8GRjO0V0JTsba26JQpCaFExZQCkBEylXBzuMYW6QQvrmjwBOrNyRhS5JxO79cYAwY7qIuQ+KgszrJgCKwGRDjcifp5mbRaGf/hhbjV6sbGRmw2GwkJCcO+npCQQG1t7Yjvqa2tve7x6vNEzvnUU08RERHhfKSlpd3Qz+MyzMGQ5phxIsJUIyJJEjGx5QCEytM1tkYnDOsoulpra3TJ1MiZyPYADKZO0pPatDZHH4iQ+JjExl1BkmSM1iTigv20b9LV+MCsO++UZRNk69attLW1OR8VFRVamySS/8bALttpl84D0NCQrrE1OsHZUXQOBEdraopeOVPZha1byU2q7j2lsTU6QV1rRAf1UWlDWWt62rKcHdT9mp4WqFYqyrw1wRjcLHBiY2MxGo3U1dUN+3pdXR2JiYkjvicxMfG6x6vPEzmnxWIhPDx82ENz1A9N2T6wWa9/rB9yueUy3bY2ZLuZoisxtHT1a22S9vjAjsrd7LvciLVLqYBRW+77PbFTICxZdFAfBVmWOd6gVPZauyY7eyj5NaX7lDEfsVMgwnvbc7hV4JjNZhYsWMCOHTucX7Pb7ezYsYOlS0dO5Fq6dOmw4wG2b9/uPH7SpEkkJiYOO6a9vZ2CgoJRz6lLkuZAYCT0tUPVMa2t0R1qR9FA62Rk2cTB4iaNLdKYoR1FvXhH5W4OFDVi65oMwLG6Y/TbhDBGkoTH+DpUdFRQ1VmFARO27qxhg379FmcxwxpNzbhZ3B6ieuyxx/jtb3/LSy+9RGFhIV/96lfp6urigQceAOC+++5j69atzuP/6Z/+iW3btvFf//VfXLhwge9///scPXqURx55BFByM77xjW/wwx/+kLfffpszZ85w3333kZyczJYtW9z947gOgxGyHHkUIvnvGtSOormRSqMyv190Gi4M6Siap7U1uqSlq5+z1W3Y+xKItsTQa+vlZP1Jrc3SB848nN1aWqFL1M1UTsQMkM0cKGpE9vfqVh/ZTLld4Hz+85/n6aef5oknnmDu3LmcPHmSbdu2OZOEr1y5Qk1NjfP4ZcuW8ac//Ynf/OY3zJkzh7/+9a+8+eabzJw503nMN7/5TR599FEefvhhFi1aRGdnJ9u2bSMw0Ms6LYpheCPSZ+vjWJ3i1botexUA+4satDRJe9Sdt5d2FPUEB4obkWWYmhDOspTBsQ0ChnRQPy06qF+F+hlZl7GCoAClg/qF2g6NrdKQlnJoLlHGfKgVeF6KR5KMH3nkEcrLy+nr66OgoIC8vMEd6O7du3nxxReHHf/Zz36Wixcv0tfXx9mzZ9m8efOw70uSxJNPPkltbS29vb18/PHHTJkyxRM/imtR3caVR6DPj/9DXcXJ+pP02fqIC4rjrunzCTBKVDT3UN7UpbVp2uHlHUU9wQGHl295TuywsQ0CRAf1UbDZbRyuUfJvVqYuIy9LSd736zwcda1JXQSBOshXvQn8oopKt0RlQtQksFuhbL/W1ugG9aa0JGkJoYEBzm60+/x10fGBjqLuRpZl5+dj5eRYZy+T803nae1t1dAyHSHKxa/hXNM5OgY6CDOHMSNmhrOr8T5/Don7UDGDEDhaI2Lj13Cw+iAw2LRNXXQO+Oui4wMdRd1NeVM3lS09BBgl8rKiiQ+OJzsiGxmZgtoCrc3TB6KD+jWom6m8xDyMBiMrJys9cA6XNtE7YNPSNG2w25WxHuD1CcYgBI72iOqGYbT0tnCh+QIw2FFUHYZ3sLgJm90PF+biIRUNXtpR1N2oO25lEr0JQISprmZYB/Uira3RBWr+jfpZmZIQSnyYhd4BO8fLW7Q0TRtqTyk9cMxhg5PovRixWmrNpFUgGaDxIrRVaW2N5hTUFCAjkxOZ4+woOjslgrBAE209A5yp8sPutCL/Zkz2X1aS0FdOHpwert60DtUcElUxoHRQT3eMIRAbKroHujnVoDSDVDdTkiT5d5hK/VxMWqmM+fByhMDRmqAoSJ6nvBZhqmt2VAAmo8E5V0i9kfkNPS1QrUw59oWYuDuw2uzOPkkrJg+22V+YsBCTwURVZxUVHTroXq4HROWmk6N1R7HaraSEppAWNji+R/UY+2WisY/0v1ERAkcPiEUHUBJF1XDC0qThTRvVG5ffJRqX7nV0FJ0K4claW6NLTle10dFrJSIogFkpEc6vBwcEMzduLiDCVE5UkVy6D2wD2tqiMc61JnkpkiQ5v656cM5Wt/lXB/X+brji6P6dvU5bW1yEEDh6QF10SnYrSV5+Snl7OTVdNQQYAliQMDz+u9Kx6By/0kJXnx+NtvChigZ3oe60l2XHYDRIw77nzMMR/XAUEudAUDT0d/h9B/XRNlPx4YFMTQhDlpXeSn7DlYNg64fwVIjJ0doalyAEjh5IXQwBIdDVAHVntbZGM9Sb0Lz4eQQHBA/7XkZMMKlRQQzYZApK/Whsg8i/GRNV4KwYkn+jot68Dtccxmr3I2E8GgbDYAd1P87Dqe2qpbitGAmJvKRrO4P7ZZjKuZlao4z38AGEwNEDJjNkLlde+3GYaqjL+GokSXImkPpNmKq5FFrKwGAa/HwIhtHZZ+X4FaXaZWVO3DXfnx4znXBzOB0DHZxt9N/NwzBESNw5iHVGzAwiLBHXfH/FkLXGbxLUi31vMyUEjl5QY55+uquy2q0cqT0CXOsyVlF7VPjNrkq9AaXlgSVMW1t0SkFJE1a7THp0MOkxwdd832gwOnfoIkzlwNlB/Sj0+mFVItffTAHkTYrGbDRQ1dpDWVO3J03Tho46qD8HSELgCNyA+qG6kg8DPdraogFnG8/SOdBJhCWC3OjcEY9Zlh2DJMHl+k5q23o9bKEGFO9Unn1owXE1+64TnlJxlotXH/KITbonMh2is0G2KcnGfoZdtjs9OKMJnGCzifkZkYCfVG6qFbxJsyEkRlNTXIkQOHohbiqEJYG1VxE5fobavVjtKDoSkcFmZjuqZHx+urjdplRQgUgwvg7q50BNQh8J1SN4uuE0XQN+PM9sKNn+G6a63HKZ5t5mgkxBziq7kVjpT5WbPrqZEgJHL0hDXIN+GKYay2WsMhgb9/FdVfUJJXwQGDHYJ0kwjJq2HorqOzFIsCx7dIGTGpZKWlgaVnkwDOr3+PFao26mFiYsJOA6zezUcvH84iasNh+ubpXlQQ+Oj22mhMDRE2oejp/tqjr6OzjTeAYYh8BxJJIeKGrE7stjG5wdRVfBKB4tf0fNxZqdGklE8PW7ri5LXgaIfjhOJq0EyQjNxdB6RWtrPMp4N1MzUyKICAqgo8/KqUofzlWqL4TOWjAFQtoSra1xKULg6Am1e2TtGej0cQ/FEI7UHsEm20gPSyclNOW6x87PiCQowEhjZz8Xajs8ZKEGiPLwMRk6PXws1DCVunv3ewIjBmcN+ZEXp8/Wx/H648Cg6B0No0FieY7aQd2Hw1TqWpOxHAICtbXFxQiBoydC4yBhlvLaj8Y2jHdHBWAxGcnLigZgf5GPisC+Tqg4rLz2MZexq7DbZed0+RXXyb9RWZS0CINkoKy9jNquWneb5x34YR7O8brj9Nn6iA+KJysia8zjVY+xz641MJh/44NrjRA4eiN7jfLsR4uOs6JhlPLwq/H55L/yA2AfgMgMiB57EfZHCmvbaerqJ9hsZF561JjHh5vDmRk7ExBhKifOfjh7/KaDuvq3X5K8ZNh4htFQvYMnrrTS0euDoy2sfVB2QHntg95iIXD0xtDkPz9oMFXdWU1ZexlGyciipEXjeo+66BwubaZ3wOZO87RBjGcYEzVksCQrBrNpfMuYKqCFwHGQuhDMYdDTDDUntbbGI4w0zPd6pEUHkxkTjNUuc6ik2Z2maUNFAVh7ICQeEmZobY3LEQJHb2QsA6MFOqqh8ZLW1rgd9WYzK3YW4ebwcb1ncnwoCeEW+qx2jpa1uNM8bfDRkk1Xsn8C4SkVNefiUM0h7LJ/eCyuizFASTYGv/AYN/U0caH5AgBLksafTDvoMfbBMNXQzZSPjGcYihA4eiMgCDIcuwv1RufDHKhW3KNjJfwNRZIkljtubPt8LTbeVgmNF0EaMjNIMIzeARuHS5Xd9HgSjFVmxc0iJCCElr4W543O7/GjDuoFNQUATI2aSmzQ+D83Pj0ixsc3U0Lg6BE/6VFhs9vG7Cg6Gs5F55KPLTrq3zxlAQSNnVvijxwubabPaicpIpCc+NBxvy/AEMDixMWAqKZyogqcK4eU5HYfRv2bT3StWeqYUl/a2EVFsw+NbehuhppTymu1gtfHEAJHj6iLTtl+sPZra4sbOdd0jo7+DsICwpwJoONFrW44X9NOY2efO8zTBmdFwy3a2qFj1FDBysmx40oUHYrqKRQCx0F0lpLMbh9Qktt9FFmWJ1StOZSwwADmp0cCPubFKdkFyBA3DcKTtLbGLQiBo0cSZkJIHAx0KUlgPopzPENSHiaDaULvjQuzMD1JydnxmR4VdttgLoQqcgXXMNj/5trp4WOhCpwT9SfoHvCh3fiNIklDwlS+GxIvbi2mvqcei9HC/Pj5E36/T+bhqH/vHN/dTAmBo0cMhiFhKt9ddNQd1bKU8effDGXlFCVMtddXFp2ak9DTApbwwSZsgmHUt/dyobYDSZpYgrFKWlgaKaEpWO1WjtYddYOFXogqcIp2aGuHG1Fz/RYmLCTQNPFmdmpI/EBRo2+MbZBlv6jWFAJHr6iqutg3F53O/k5ONSjx3/H2v7maVUP64ci+UFKvitlJq8A4MY+Wv6B6b2anRBAVYp7w+yVJEmGqq5m0Sklqb7rss2MbbjQ8pTI7NZLwQBPtvVZOV/nA2IaGi9BepVTsZizX2hq3IQSOXlE9ODWnoMtHQjBDOFx7GJtsIyM8g9Sw1Bs6x8LMKAIDDDR09HGxzgfGNhSL8NRYDObfTDw8pbI8WVnQD1T5bs7JhAiKhJSFymsfLGzos/U5vXXq336iGA3S4KBfXyhsUDdTGcuUyl0fRQgcvRKW4NNjG5wVDTfovQFlbMOSLGVWzN5LXh6m6usYzLcSAmdE7HZ5QvOnRmNR0iKMkpGy9jKqO6tdZZ534/QY+15I/FjdMWU8Q3A82ZHZN3wetbDBJ/Jw1MiAD+ffgBA4+kaNjfpgbNyZfzOB/jcj4TNjG8r2g90KUZMgepLW1uiS8zXKeIaQcY5nGI1wczizYpXNg+hq7EAV1SW7lWR3H2LoWjPRqruhOMc2VLTS7s1jGwZ6B8cz+PhmSggcPTO0usEXckwcVHRUcKXjCibJxKLE8Y1nGI3VjkTjAm8f2+AsD/ftBedmUEXs0uzxj2cYDZGHcxXJ88ESAb2tUH1Sa2tcivo3vtnNVFp0MJNiQ7DZZfKLm1xhmjZcyVfGM4QmQvx0ra1xK0Lg6Jn0pWAKgs5aqD+vtTUuQ91RzY6bTah5/I3aRiI7LpSkiED6rXYKSr14VowQOGPiivwbFTXZ9FDNIWw+5rG4IYwmyFqlvPahMFVDdwOXWi4hIU1oPMNoDHY19uIw1dC1xgfHMwxFCBw9ExAImY6kOB9adFy1owKlKmawq7GXLjot5dBUBJJxcDaQYBjd/Vbn3LGbyb9RmRk7kzBzGO397ZxrOnfT5/MJnB5j3wmJq8M1p8dMJyrw5juD+0RI3I82U24VOM3Nzdxzzz2Eh4cTGRnJgw8+SGfn6O3Am5ubefTRR5k6dSpBQUGkp6fzj//4j7S1DS/LkyTpmscrr7zizh9FO7J9K/nParc6Z8IsT3FNeaLXLzrq3zZtMQRGaGuLTikoaabfZic1KohJsSE3fT6TweTc0YswlQP1hldxGHrbtbXFRbhyMwWwJCsak0GivKmb8qYul5zTo3TUQt1Z5bUP979RcavAueeeezh37hzbt2/n3XffZe/evTz88MOjHl9dXU11dTVPP/00Z8+e5cUXX2Tbtm08+OCD1xz7f//3f9TU1DgfW7ZsceNPoiHqolN+EAZ6tLXFBZxtPEvnQCcRlgimRU9zyTlX5MQiSXCxroO69l6XnNOj+NGO6kbZOyQ8dTOJokNRw1Qi0dhBVCZEZ4Nsg7J9Wltz09hlu8uKGVSUsQ2KJ8grN1RqRW7SHAi5eU+o3nGbwCksLGTbtm387ne/Iy8vjxUrVvDss8/yyiuvUF09cmnmzJkz+dvf/sYdd9xBdnY269at4z/+4z945513sFqtw46NjIwkMTHR+QgMnHh3Sq8gbiqEp4C1VxE5Xo66o1qStASjweiSc0aFmJmdong+vK5c3GaF0j3KayFwRkW9maxyQXhKRW1RcKrhFB39PtBHyRX40NiGi80Xae5tJtgUzJy4OS47r1fn4agVuX4y685tAic/P5/IyEgWLlzo/Nr69esxGAwUFIx/vlJbWxvh4eGYTMM7u379618nNjaWxYsX88ILL1y3k21fXx/t7e3DHl6DJA26En1g0VGbq7lqR6XitWGq6uPQ26aEppLnaW2NLqlu7aGovhODBMuyXSdwUsNSyQjPwCbbOFx72GXn9Wp8aGyDOp5hceJiAowBLjvvqinKWnOwqIkBbxrbYLf73aw7twmc2tpa4uPjh33NZDIRHR1NbW3tuM7R2NjID37wg2vCWk8++SSvvfYa27dv59Of/jRf+9rXePbZZ0c9z1NPPUVERITzkZaWNvEfSEt8JA+ntbeVM41nANcLHHXR2V/UiN3uRSX16o0kay24yKPla6g75TlpkUQEu+5GBYOfQ9HV2MGklWAwQUspNBVrbc1NcbOz7kZjVkoE0SFmOvqsnLjS6tJzu5W6M9DVAAEhkJantTUeYcIC59vf/vaISb5DHxcuXLhpw9rb27n99tuZPn063//+94d977vf/S7Lly9n3rx5fOtb3+Kb3/wmP/3pT0c919atW2lra3M+Kioqbto+j5K1BpCUUvH2Gq2tuWHya/KRkcmJzCExJNGl556XHkmI2UhzVz/nqr3IQ1f0sfKcs15bO3TMHkfYcfWUmy8Pvxq1df/B6oO+Mc/sZrGEKe0pwKs3VN0D3RyvPw64fjNlMAxWbu65VO/Sc7sV56y7lWCa+Bw3b2TCAudf/uVfKCwsvO4jKyuLxMRE6uuH//GtVivNzc0kJl7/5tbR0cGmTZsICwvjjTfeICDg+ru2vLw8Kisr6evrG/H7FouF8PDwYQ+vIjgaUuYrr7140VF3yStSVrj83AFGA8tyvGy6eHczVB1TXvt4y/QbxWqzO8OO7hA4ixIXEWAIoKqzirL2Mpef3ytRP4uq+PZCjtQewWq3khKaQnpYusvPr34W93hTzp+f5d/ADQicuLg4cnNzr/swm80sXbqU1tZWjh075nzvzp07sdvt5OWN7h5rb2/n1ltvxWw28/bbb48refjkyZNERUVhsVgm+uN4D9nevejIsuyMibuqPPxqnIvORS9ZdIp3AjLEz4DwZK2t0SUnK1rp6LUSGRzA7NRIl58/OCCYBQkLABGmcqJ6E0v3gnXkTaPe2V+1H1A2U66quhuKmvN3tqqdhg4v+B31dcCVQ8prP8m/ATfm4EybNo1Nmzbx0EMPcfjwYQ4cOMAjjzzC3XffTXKysphXVVWRm5vL4cNKgp8qbrq6unj++edpb2+ntraW2tpabDal2+g777zD7373O86ePUtRURG//OUv+c///E8effRRd/0o+kBddIp3euWsmEstl2jsaSTIFMT8+PluuYYqcI5daaGtxwtmxRT5x8C7m0HdIa+cHIfR4J6uq6pHcX/1frec3+tImAmhCTDQrbT19zJkWR4mcNxBXJiFGclKJMArqqlK94F9QGkFEHPjA0e9Dbf2wfnjH/9Ibm4ut9xyC5s3b2bFihX85je/cX5/YGCAixcv0t3dDcDx48cpKCjgzJkz5OTkkJSU5HyoeTMBAQE899xzLF26lLlz5/LrX/+aZ555hu9973vu/FG0J2WBUmnT2zoY1vAi1AVnceJizEb3xH/TooPJilNmxRws0nk1ld0u8m/GgTvzb1TUPJyjtUfptXphHyVXI0mDn0kv9Bhf6bhCZWclJoOJxYmL3XYd9TPpFa0pirYrzznrfX48w1BMYx9y40RHR/OnP/1p1O9nZmYOS+xbs2bNmIl+mzZtYtOmTS6z0WswmhTX4rk3lEUnzX3/cd2BGp5ydcLf1ayZEk9JQyl7LjVw26wkt17rpqg7C131SkVD+s3PyPFFGjv7OF2pdDF3Zf+bq8mOzCYhOIG67jqO1R1zWwjVq8i5BU7+UfEy3vpDra2ZEOpmakH8AoIDgt12ndVT4vjF7mL2XlYqNw1u8jDeNLI8ZDO1QVtbPIyYReVNqLuqy9u1tWOCdA10caLuBOA+l7HK6qnKrmr3xQZ9V8WoC86klWDy4dyxm2C/I7l4elI48eHua+QpSdJgmKpKhKkApW2BWrnZVqW1NRNC/Ru6W6jOz4gi1GKiuaufs9VtY79BK5qKoPUKGM1+N+tOCBxvQhU41SegS+chmCEcrjmMVbaSFpZGerjrKxqGkjcpGovJQG17L5fqRp97pjnO/BsRnhoNZ3hqqvvCUyrqzVD1NPo9wdFKWBy8avhmr7WXo7VHAfdvpgKMBpbnxAA6L2xQN8QZy8B883PcvAkhcLyJsERImAXIXlUu7qyeSna/6z8wwMiSLMeio9ceFb3tUOGoaBAJxiNit8vO3AZ35t+o5CXlYZSMlLaVUtXpXR4Lt+GFeTjH647Ta+slPjienMgct19vlTeUi/txrp8QON6Gl/Wo8ERFw9WsmarzRad0L9itEJ2lPATXcK66naaufkItJudwQ3cSbg53zisS5eIOnJWbu5WZaV7AviplSKi7ysOvZpWjXPxERas+KzcHeqDc8XkWAkegeyY7ksSKdiiVODqnvL2cqs4qAgwBLEpc5JFrqjv+I6UtdPXpcGEuFuGpsVC9b8uyYzCbPLNMqQJcCBwHKfMhMBL62rymclP1FntqM5UWHUy2nis3y/Yrg5rDUyEuV2trPI4QON5G6mIwh0F3I9Sc1NqaMVEXnPnx891a0TCUSbEhpEUH0W+zk1/c5JFrjpthFQ1C4IyGJ/NvVNQ8nILaAgZsOtyNexqDccjwTf17jKs6qyhtK8UoGclL8tyspdVTlJmLuvQYO9eaW/yqPFxFCBxvw2SGrNXKay+Y+OupioahSJLEGr0uOkMrGjI9s8v0Ntp6BjjuGGKohgA8QW50LtGB0XQNdHGy4aTHrqtrvCgPR/W8zYmbQ7jZc+N4Vg8JieuucvPykP43fogQON6Ic9HRd7l4n63PWdHg6d4iaphq96V6fS066o3CDysaxsvBokZsdpnsuBDSoj3j9QMwSAZnIrwIUzlQc/68oHJTi80UDFZu1rTprHKzuQSai5Xp8Oqm2M8QAscbUQVO5RHoadHWlutwpPaIs6JhcuRkj157aXYMAUaJiuYeShu7PHrt6+LnO6rxMNi9ON7j11ZvjqIfjgMvqdwcsA1QUFMAeC7/RmVo5ebuizqq3FQ9/Gl5Shd8P0QIHG8kMk1JGJPtULxLa2tGZV+lUtGwMmWlRyoahhJiMbEoMxrQUZiqv0tJ+gMhcEZBlmVN8m9UliUvQ0LiYstF6rt1dLPSkslqg9GPtLXjOpyoP0G3tZvowGhyoz2fTLvW8VndpSuBI3L9hMDxVpxhKn3m4ciyzN7KvQCsTNWme+aaIV2NdUHpXrD1QUS6X1Y0jIfCmg5q2noJCjCSNyna49ePCoxiVuwsQHhxnEy+VXku+li3g37VQanLk5djkDx/W1uXmwDA0bIW2nt1kKBu7VPWGxACR+CFDO2Ho6ccEwdl7WXOgXdLk5ZqYoMa4jhU0kRPvw4WZnUHPOVWv6xoGA/qDnh5TgyBAUZNbFiRqoQ4VIHu96QuVkIcPS1QeVRra0ZE9RZrNUcsPUYZ9Gu1y84RI5pSflCZBh+aAImztLZGM4TA8VYylkNAMHTWQu1pra25BnXBWZiw0GPl4VczJSGUlMgg+qx2DhZrvOjIMlxyCBx1Ryy4hl0XFIGzZqrn829UVqWuAiC/Op9+W79mdugGowmyHRuqyx9qa8sIVHdWU9RahEEyeDz/ZijrHJ9Z9TOsKUNz/fx4MyUEjrdiskDWGuX1Jf0tOnurlN2verPQAkmSWJurhKl2ar3o1BdCeyWYAiHTvwbejZeWrn6OX1GS5tfmaidwpkVPIzYolm5rN8fqvKPBnduZslF51mEejuppmxs3lwiLdsm06md218UG7HaNveqXtinPfr6ZEgLHm1EXHZ0JnK6BLueNYWWKtjfzdbmDuypNy8XVne+kVWDWxqOld/ZebsAuQ25iGCmRQZrZYZAMzs+tCFM5yFkPSFB7BtqrtbZmGFrn+qksyowmxGyksbOPc9Xt2hnSVDxYHq42avRThMDxZlR1XnVMVz0qDlUfwmq3kh6WTmZEpqa2LM2KxWIyUN3Wy8W6Du0MEeGpMVG9bFp6b1TUm6VINHYQEjs4XVxHXpweaw+Haw8D2nqLAcwmAysmxwIae4zVDW/GMgj0XMNDPSIEjjcTnuxIIJMHY646QB14p/WOCiDIbGRZttKjQrNFp6cFKpQeHULgjIzNPlgevlbD/BuVpUlLMRlMlLWXcaX9itbm6ANnmEo/a82R2iP02fpIDEn0eK+tkVA/u5qWi6ve4skbtbNBJwiB4+1M2aQ8qzFXjZFl2ZlgvCpF2x2VytAwlSYU7wTZppSGR2VoY4POOVnRQmv3ABFBAcxPj9TaHELNoSyIVzwWIkzlQBXnxbuUMmQdoP5tVqWs8nivrZFQvY+nKltp6tTgd9TXAWWOLtzqvcGPEQLH21FVevFO0MGAwIstF6nvqSfIFMSCxAVamwMMLjrHylto7dagKsYZntrg+Wt7Cap3bdWUOExGfSxLqgdSCBwHibOVsuOBLijXfpSFLMvsqdwDaB+eUkkID2R6UjiyrFGD0eJdYB+A6CyIzfH89XWGPlYSwY2TMh+CY6CvHa4c0toa580gLykPi9GisTUKqVHBTEkIxa7FomO3D84MEy7jUdl5Qfm7rMv1fPfi0VAFztG6o3QPdGtsjQ4wGAZFug7CVJdbL1PbVYvFaGFx0mKtzXGybkg1lccR4alhCIHj7RiMkONYdHQQpho6nkFPrNUqTFV9HLqbwBIO6Us8e20voaath8KadiTJs9PDx2JS+CRSQ1MZsA9wqEb7zYMuUMNUOqjcVDdTixMXE2TSruruatTWFHsu1mO12T13Ybt90Fs8ReT6gRA4voFOelS09rZyulFpOqgXl7HKLY5W6nsuNWDzZI8K9W+SvRaMAZ67rhehjtKYmxZJTKg+vH6g9FFSP8ciTOUgay0YApQy5KZiTU1x5vrpbK2ZmxZFZHAA7b1WTlS0eu7CNSehqx7MoUojWIEQOD5B9jqQjNB4CZpLNTPjQPUB7LKdyVGTSQxJ1MyOkZifHklEUAAt3QOcrPDgBPZLwmU8Fmr+zTodVE9djXrz3Fe1T9s+SnohMBwyHKNXNNxQtfW1cbLhJKA/gWM0SKyeokGDUfXvkbVGaQQrEALHJwiKhHTtFx014U9v4SkAk9HAKseis6PQQ4tOR62yqwK/Hnh3PfqsNg4UKT2c9ND/5moWJi4kyBREfXc9F1suam2OPtBBmOpAlbKZyonMITk0WTM7RkOTyk317zFFbKZUhMDxFZxdjbXJwxmwDziboq1NW6uJDWOxztNjG9REzOR5EJbgmWt6GQUlzXT324gPszAjWX9NySxGC3mJeYAIUzlRvZHlB5SyZA1QR8HoodfWSKyeEofRIHGhtoOKZg8kqHfUKfl+IHptDUEIHF9BFThl+6Gv0+OXP153nI7+DqIDo5kVq8/ptaunxCNJcKG2g+rWHvdf8OL7yvPUze6/lpeyo7AOUBqk6aGPyUisSlNCIHsq9mhsiU6InQzR2WDrh6IdHr+8zW5zbqZWp672+PXHQ2SwmYUZUQB87PiMuxW1UjNpLoTpKz1AS4TA8RVip0BkhrLolHp+Id5dsRtQ4uFGg9Hj1x8P0SFm5qVFAh7w4vR3Kz0pAKbe5t5reSmyLPOxI1y4Ybp+PVzqTfR042kae/QzEkUzJGnwM33xA49f/nTjadr62gg3hzMnbo7Hrz9e1M+0RwSOCE+NiBA4voIkDXau9PCiI8syuyqUm/matDUevfZEuWWahxadkt1g7YGIdEiY6d5reSnna9qpau0hMMDA8pxYrc0ZlfjgeGbGKH9DVcj7PapX8vKHYLN69NK7rihrzYqUFZgMJo9eeyKoa01BSTNtPW5swmrtG9xMiWKGYQiB40tMHTK2wW7z2GWLWouo6qzCbDCzNGmpx657I9zq2FUdLGqis8+NC7MzPHWbIj4F1/DxecV7s3JyHEFmfXr9VNamK3llqpD3e9LyICjaMWfNsz2C1L+B+jfRK5NiQ8iJD8U6ZM6aWyjdB/0dEJqo5PsJnAiB40tkrABLBHQ1QOVRj11W3dUuSV5CcECwx657I+TEh5IZE0y/zc5edy06dvtgsrcIT43K9sJaADZM0294SkVNnD9UfUh0NQYwmgbDIR70GJe0lVDWXobJYGJF8gqPXfdGWa96jM+70WN88T3leeptSrdpgRPx2/AlTObBVurqh94DOHdUOq2eGookSc7Y+HZ3LTpVRxWRaYmATP0vwlpQ09bD2Sqle/G6aforD7+anMgcUkJT6Lf3k1+dr7U5+kAV7xfeAw/1CFLDU3mJeYSaQz1yzZtBXWt2XaxnwB1dje32QYGZe7vrz+/luFXgNDc3c8899xAeHk5kZCQPPvggnZ3Xr/BZs2YNkiQNe/zDP/zDsGOuXLnC7bffTnBwMPHx8Tz++ONYrZ6NA+uWXEds/ML7HrlcQ3cDZxrPAPqtaLiaDdOVKoOdF9y06KjhqcnrRffiUVB3tPPTo4jVUffi0ZAkySngRZjKQfYtYDRDSyk0eKZHkDdtpkDpzh0baqaj18qR0mbXX6D6BHTUKN2LJ+mr4aEecKvAueeeezh37hzbt2/n3XffZe/evTz88MNjvu+hhx6ipqbG+fjJT37i/J7NZuP222+nv7+fgwcP8tJLL/Hiiy/yxBNPuPNH8R5yNiit1JsuQ8Mlt19Obe43K3YWccH6mSN0PRZkRBEdYqatZ4AjZW5YdNQdlSgPH5XtXlA9dTXr0tcBSj8cmwdz3HSLJRQmOTY1F92/oWrsaeR0gzIKRu/FDCpGg+Rs+rfdHYUNqqc+Z73oXjwCbhM4hYWFbNu2jd/97nfk5eWxYsUKnn32WV555RWqq6uv+97g4GASExOdj/DwwQZgH330EefPn+cPf/gDc+fO5bbbbuMHP/gBzz33HP39/e76cbyHwPBBJe+BMJWaf+MtCw5ctei4OkzVVAwNF8BgEt2LR6Gjd4D8YqXcer0X5N+ozIufR7g5nJa+FueYAL/HWS7ufoGzp2IPMjIzYmaQEOI9nxv1M779fJ3rx32onnoRnhoRtwmc/Px8IiMjWbhwofNr69evx2AwUFBQcN33/vGPfyQ2NpaZM2eydetWursHk/ry8/OZNWsWCQmDH/CNGzfS3t7OuXPnRjxfX18f7e3twx4+jYfCVN0D3c4py94kcIBheTguXXRU703GcmWEhuAa9l5qZMAmk+WoMvEWTAaTc+6RKBd3oAqcyqNKN1034m3hKZUVk2OxmAxUtvRwsc6FnZ+biqGhUJlDqOZeCobhNoFTW1tLfPzw5EGTyUR0dDS1tbWjvu+LX/wif/jDH9i1axdbt27l97//Pffee++w8w4VN4Dz36Od96mnniIiIsL5SEtLu9EfyztQQyOVR9y66ByqOUSfrY+U0BQmR05223Xcwcohi86FWhcuOiLhb0y2n3dUT3lReEplaB6OGL4JhCc7SpNlpSeOmxi6mdJ7efjVBJtNrHD0eXJpNZXqNctcAUFRrjuvDzFhgfPtb3/7miTgqx8XLly4YYMefvhhNm7cyKxZs7jnnnt4+eWXeeONNyguLr7hc27dupW2tjbno6Ki4obP5RUMXXQuua+Ec+iOSq9t9kfDLYtOdzNcOai8FuXhIzJgszu7SK/3QoGzPGU5AYYAytvLKW0r1docfTDV/R7j/Op8r91MweBnfbsrB/2K8NSYTFjg/Mu//AuFhYXXfWRlZZGYmEh9/fA/ptVqpbm5mcTE8c/KyMtTBt0VFRUBkJiYSF3d8BuS+u/RzmuxWAgPDx/28HnUD72bFh2b3eYcPuht4SkVZ5jKVcl/lz8C2Q4JsyAy3TXn9DGOlDXT3mslOsTM/HTv23WGBISwOGkxIKqpnKgCp2SXMqLEDeys2Al452YK4BZHzt+pilbq23tv/oRdjYMNFsVmalQmLHDi4uLIzc297sNsNrN06VJaW1s5duyY8707d+7Ebrc7Rct4OHnyJABJSUkALF26lDNnzgwTT9u3byc8PJzp06dP9MfxXaY6BE7JbrcM3zzVcIrm3mbCAsKYnzDf5ef3BLdMS0CS4HRlG7VtLlh0hnYvFoyI2r14XW48RoP33agA1qUp1VRC4DhImKGMJLH2KiLHxVjtVudmSq1k8zbiwwOZ45iD55IN1aVtymYqcbbYTF0Ht+XgTJs2jU2bNvHQQw9x+PBhDhw4wCOPPMLdd99NcnIyAFVVVeTm5nL48GEAiouL+cEPfsCxY8coKyvj7bff5r777mPVqlXMnj0bgFtvvZXp06fzpS99iVOnTvHhhx/yne98h69//etYLKJMzkn8NIjKBFsfFO90+em3lyvTa1enrSbA4J29XuLCLM7hmze96Az0wOWPlddC4IyILMuD3Yu9MDyl4hy+2SCGbwLKKBK1sKHwXZef/mT9SVr7WomwRDAv3ntHEWycoXzmt50dPQd13Ijw1Lhwax+cP/7xj+Tm5nLLLbewefNmVqxYwW9+8xvn9wcGBrh48aKzSspsNvPxxx9z6623kpuby7/8y7/w6U9/mnfeecf5HqPRyLvvvovRaGTp0qXce++93HfffTz55JPu/FG8D0ka9OJccG25uCzL7LiyA4D1Gd5dCq02/bvpcvGiHTDQpexkxTyYETlf005Fcw8Wk4GVk/U7XHMsEkISmBkzExmZnVdcv3nwSqbdoTxffB9srh0sqXrKVqWs0vVwzbG4baYShcgvbqK1+yZamvR3D25aRa+t6+LWT0t0dDR/+tOfRv1+ZmbmsEqEtLQ09uzZM+Z5MzIyeP99z3Tq9Wpyb4dDzw1O/DW65s99rukcNV01BJmCWJ683CXn1IoN0xP48bYL5Bc30tYzQETQDXqjCt9WnqfdIYZrjoK6c10zNY5gs/feqAA2ZG7gbNNZtpdv53NTP6e1OdqTvhRC4pQRJaV7IecWl5x26GbK26qnrmZSbAi5iWFcqO1g+/k6PrvwBqt5S3aBtUfZTCXOcq2RPoaYReXLpOVBcIwy8bd8v8tOq4anVqasJNAU6LLzakFOfCiT40MZsMnsuNEwlbVvsDx8+l2uM86HkGWZ987UALB5VpLG1tw8G9KVviNHao/Q2tuqrTF6wGCE3E8or8+/5bLTFjYXUtVZRaAx0Os3UwCbZioe45sKUxU6Ihq5m8VmagyEwPFljKbBGO25N11ySlmW+bhcyTXZkOkbzaVuc9xw3z9zg4tOyR7oa4ewJEhd5ELLfIfL9Z2UNHRhNhqcXaS9mbTwNHKjc7HJNpFsrDL9TuX5wnvgolEWH5V9BMDK1JUEBwS75Jxaooap9l1upKP3BkJ51r7B/BuxmRoTIXB8nelblOfCd1yy6FxqucSVjitYjBZWpfjGcLfbHQJn7+WGG1t0Ch071txPgEH8lxqJ9x3em5WTYwkL9M6k9KvZkKEI/I/KP9LYEp2QuVJpONfdCOUHb/p0siw7vcW3Ztx60+fTA1MSQsmKDaF/SD+oCVGyB/raIDQR0pa43kAfQ6zGvs6kVUMWnQM3fbqPryjem2XJy3xiRwWORScuhH7rDSw6toHBJG6xoxoV1SV/mw+Ep1TUBPtDNYdo7/fx8S/jwRgwWNjggjDVsM1Uqm9spiRJurkw1fk3ledpd4jN1DgQvyFfxxjg0jCVMzyV4RvhKVAWnc0O1/EHEw1Tle1XcpyCYyFjmRus836KGzq5UNuBySCxwYuGa45FVkQWOZE5WO1W9lSMXRzhF6hhqsJ3wG6/qVN9WKaMfliRssJnNlMwGKbafbGBnv4JeNWt/XDBUYY/Y4vrDfNBhMDxB6Z/Unm+yTBVSVsJRa1FmAwmVqetdpFx+uC2WcquatfFerr6rON/o1o9lXu7kmgpuAZ1p7osJ5aIYN8IT6mIMNVVZK0BSzh01iqz8G6QoeEpX9pMAcxMCSc1KoieARt7Lk3AY1y6F3rbICReqVoTjIkQOP7ApFUQGAFd9XAl/4ZPo3pv8pLyCDf71riL6UnhZMYE02e1s+viOBcdu22wokGEp0blg7OO6qmZ4x/R4i2oYaqDVQfpGujS2BodYLLAlE3K65sIU11uvUxZexlmg9nZWNFXkCSJTTOU/wsfTCRMNSw8JTZT40EIHH/AZHZJCaczPJXuWzsqUBYdNT9k3GGqK4eUvh+BkYqIFFzDlaZuzla1YzRI3DrD9wTO5MjJZIZn0m/vd44T8HtUsV/4DtzgxHXVe7MsZRmh5lBXWaYb1LVmZ2E9fdZxeNVtA4PhKbGZGjdC4PgL6n+K82/fUGy8oqOCwuZCjJLRa+fBjIWah7PzQv34YuNqeGrqZiXXSXANqvcmb1I00SFmja1xPZIkOUMo6k3Z78m5BQJCoO0KVJ+4oVNsL/Ot6qmrmZcWSUK4hY4+KweKxjHuo2zfkFw/7+8H5CmEwPEXstaCJUKJjatTaCfAjnKlm+jChIVEBXrfFOjxMDQ2vnusMJXdrohFGEysFFzDBz5YPXU1aphqf9V+ugfcM03bqwgIgskOL+8NeIyLW4spbiv2yVw/FYNhSJhqPB5jtUBk2idc1pHeHxACx18wmQcH4t3AoqMmUXr77KnrIUmSsyfO+2PFxquOQUc1mMMU8Si4hqrWHk5WtCJJg4MGfZFp0dNICU2hx9rDgeqbb8XgEzjDVG9POEylrjXLkpf5XK7fUDY5PMbbC+vot17Hq26zDglPbXG/YT6EEDj+hDNM9daEwlQV7RWcaTyDQTL4tMCBobHxOnoHrhOmOve68jxlIwR497gKd6FWTy3KiCY+zHd/R5IkOUMpauddv2fyrWAKhOYSqDk1obf6avXU1SyeFE1cmIXW7gH2FzWMfmD5fuhugqBopZmiYNwIgeNPZK9TPA4dNVB5eNxve79UaQ2+JGkJsUHeOwV6PMxJjSA5IpCufht7L42y6NhtcPZvyutZn/WccV7G26eqAdg8y/eSi6/m1kxF4Oyp3CPCVACWUJh6m/L6zF/G/bbStlIut1zGJJlYm+bbnlGjQeITs5UN1Vsnq0c/UISnbhghcPwJk2Vw0Rln0z9Zlp0CZ/OkzW4yTD8MraZ653TNyAeV7YPOOqVDdLZvJlzfLGWNXZyqaMUgwe2zk7U2x+3MiJlBelg6PdYedlbs1NocfaCK/7N/G3f/rW2l2wDIS84jwhLhLst0w11zUwD46Fwd3f0j9N8S4ambQggcf2Pmp5Tnc6+Pa9G51HKJkrYSzAYzt6Tf4mbj9MFdc5Ub8vbztXSO1PRP3ZFO36LkNgmuQd2RLs+JJS7MorE17keSJG7PUjqGv1fynsbW6ISc9Ur/rY6acc2mkmWZ90qV393tk253t3W6YE5qBBkxwfQM2Pi4cITChtI9SiuKoGjRiuIGEALH38i+RfE8dNYp/3nGQF1wVqet9sl+FCMxKyWCrNgQegfsfHTuqmTjgV4472juJ8JTIyLLMm+dqgJgi2OH6g+oHs786nyae5s1tkYHmCyDeX/jCFOdbTxLeXs5gcZAn21FcTWSJHHnHGVD9fbJqmsPOP2a8jzzU6IVxQ0gBI6/YTLDDIcXR/3PMwp22c4HpR8A/hGeUpEkyek6fuPEVYtO0XZlmm94imiXPgrnqtspaejCYjJwqw9XT11NZkQmM2JmYJNtzjlKfo+6CTj/Flj7rnuouplam76WkIAQd1umG1SP8Z5LDbR29w9+o79rsFP67M9rYJn3IwSOPzL7c8pz4TvQP3pC5In6E9R21RIaEMrKVP/K3lcXnQNFjdR39A5+Q92Jzvy0mOY7Cm86ROH66QmEBfrXrlPdCLxf8r7GluiEjOUQlgS9rVC0Y9TDrHarczP1iaxPeMg4fZATH8a0pHAGbPLw0Q0XP4CBLojKhNRFmtnnzYgV2h9Jy4PIdOjvhIujL8TqgnNL+i1YjL6fRzGUzNgQ5qZFYpfh3VOOZOPeNrioJEGK8NTI2Owy75xW8m/umuP7ycVXc9uk25CQONlwksqOSq3N0R6DUdkMAJwZ3WNcUFNAc28zUZYolib7n2dU3VC9NTRMdfpV5XnW50CSNLDK+xE1Z/6IJCkuz70/VcJUsz5zzSED9gGnm31zlv+Ep4byyXkpnKxo5a2TVfzdiklQ+C7Y+iB2KiTO0to8XVJQ2kRbdz+TY8wsyQijt7d37Df5EGGGMG5NuZWTDSfZVrKNr8z5itYmac+sz0D+zxWPRF8HWMKuOeTdEqVSaGPmRgIM/uX1A7hjTjI/+uACBaXN1LT1kGTsHPR4qR53wYQRAsdfmfU5ReAU74CuRggZ3t8mvzqf1r5WogOjWZy4WCMjteX22Uk8+e55TlW2UdLQSZYanpr1WbGjGgFZlrlUUsEzG+OJDA6gquKK1iZpwhfiv8CmiE1I/RLV1dUkJSUh+fPnJWkuxORAUxFceA/m3D3s290D3ey4otzM1Uo0fyMlMohFmVEcKWvh3VM1PGT5GGQbJM+D2Mlam+e1CIHjr8RNURaempNw7g1Y/NCwb6vhqU2ZmzAZ/PNjEhtqYeXkWHZfbOCjw6f5B7XqbASPlwAqq6uZFCEREZVAWlwkIX6Wf6Nis9sobSvF3menqaUJSZJISvLdWVxjIknKpmD3U0oO21UCZ3fFbnqsPaSGpjInbo42NuqAO+emcKSshbdPVfNQoCM8JZKLbwr/vHMJFGZ/XhE4p18dJnB6rD3OHdVtk27TyDh9sGVuCrsvNtB/6m8g25Vkv+hJWpulO2w2G3UNzURExWAJjSA6ItSvvRYRtgjaA9oJMgfR2tpKfHw8RqNRa7O0Y+ZnFIFTvAs6GyA0zvktZ++brNv9+jOzeWYi33/7HB3VF8ByFKQh+UuCG0IkGfszMz8NkgEqjygzYxzsvLKTHmsPKaEpfr2jAtgwPYGgACMre3crXxDJxSMyMDBA74ANyWQmMjjAr29UABFmpQtvr6EXGZmBgQGNLdKY2Bwl3CLbFI+xg+beZg5UKQNK/TU8pRLj8BjfZXA0RcxaA6Hxmtrk7QiB48+EJSj/iQBODzbieqNIWYDuzL7T729UIRYT9+X0Ms9QhA0jzPik1ibpks6+AXqtdpAkIoP8MzQ1lFBzKAbJgE220W/rH/sN/sAsR7LsqT85v/Rh2YfYZBvTY6YzKUJ4Rj81L4Utxv0A2GeJ5OKbRQgcf0eN8Z55DWSZyo5KCmoKkJC4K+cubW3TCfdY9gKwj3kM+Piw0RtlZ2E9sgxmk5HAAD8OxTgwSAbnLCUxfNPB7M+BIQCqT0DtWWBwrIW/9b4ZjY2RVUwy1NEtWzgYsERrc7weIXD8ndzbwRSkVDhUHuWt4rcAyEvKIyXUf9rsj4ptgLSKtwH4Q/9qdl4YYV6MwNmgLCLI+8JTsizz8MMPEx0djSRJnDx50iXnjbREAtBr7aWzv9Ml5/RqQmIHh/2e+D3FrcWcajiFUTL6fa6fiuWc0ivoI/sCXjklxn3cLELg+DuWMOe8GNuxF3mz6E0APjX5UxoapSMubUPqaqAzIJrd9jm8eqRCa4t0x7nqNi7VdSBJEB7ofXUL27Zt48UXX+Tdd9+lpqaGmTNnuuS8QaYgzEYzMjJ7Ksee++YXzL9PeT79Kq9fVMLiq1JXESs8o0pXecf4nL/ZVvHRuTpaukR482YQAkcAC74MwKGid6ntqiXcHO43w+7G5PjvAbDOvBsrJnZfrKemrUdjo/TFaw7RFxRgxGT0viWluLiYpKQkli1bRmJiIiaTa0SaJEmEm8MB2F6+3SXn9Hqy10F4Cv09Lbx9+XUAPjNFtF0AlHldfW0QmUFz/DL6bfbhnY0FE8b7ViOB60lfCrFTeCNIyZ24Pet2vxvNMCLt1cpwTSBy+d+xeFI0dhn+clS04FfpHbA5B5IGm70v9+b+++/n0Ucf5cqVK0iSRGZmpkvPH2YOQ0KipLWE803nXXpur8RghLlfZGdIMK22HuKD41mWvExrq/TB8ZeU5/lf4nOLMwB49WglsixraJR3IwSOACSJ1jmfY2dIMACfzBGVQgCc/JPS+yZ9KcRO5guL0wB49UgFdrtYdAA+PFdLe6+V+LBALKbhAkeWZbr7rZo8xntT+J//+R+efPJJUlNTqamp4ciRIy79/ZgMJgJNgQC87vBY+D1z7+FvYcq08E+mrvPbRqLDaLgIV/KV3jdz7+WuucmYjQYKa9o5V92utXVei1s/Wc3NzTz66KO88847GAwGPv3pT/M///M/hIaGjnh8WVkZkyaNXCr42muv8dnPKj1IRkpi/POf/8zdd999zdcF4+O9sAgGJIlpff1M6xdxX+x2OPEH5fW8LwFw28wkvvfWOapae9hf1MiqKXHXOYF/oOYkbZyRcM30ip4BG9Of+FADq+D8kxsJNo+9vEVERBAWFobRaCQxMdEttgQHKBuH90ve518X/qtT8PgrFQEmDgUFIckyn+wUFWYAHH9ZeZ6yCcKTiARunZHAu6dreO1oBTNTIrS0zmtxqwfnnnvu4dy5c2zfvp13332XvXv38vDDD496fFpaGjU1NcMe//7v/05oaCi33TY8y/7//u//hh23ZcsWd/4oPo0sy7x+RbkRfbKjE469pLFFOqD8ALSUgjkMZmwBIDDAyCfnKZVlrxzxzzlLQylv6uJgcROSBJtmukcc+AJmo5n4kHg6BjpELg7wxmWlz9bSnl5Szr4FdpvGFmmMtU/xFoMzHxLg84sUj/GbJ6roHfDz39EN4jYPTmFhIdu2bePIkSMsXLgQgGeffZbNmzfz9NNPk5ycfM17RtpFvfHGG3zuc5+7xusTGRnpth2Xv3G+6TyXWi5hlkxs7upS5sXc+gMwh2htmnaoO6qZnxr2e7h7cTov5Zez/XwdjZ19xIb6b66Smou0IieWxIggSq+qag0KMHL+yY0aWKZcWy9ISKxPX8+pllO8UfQGd2TfobVJmmG1W3mrSGlF8ek+O7RVQMluyLlFW8O0pPAd6GmG8BTIWe/88vLsWFIig6hq7eHDc7XcNVe07ZgobvPg5OfnExkZ6RQ3AOvXr8dgMFBQUDCucxw7doyTJ0/y4IMPXvO9r3/968TGxrJ48WJeeOGF68bc+/r6aG9vH/YQDKJ2Lr4lYwMRERnQ1w7n3tTWKC3paYVCpfeNs6zVwbSkcOakRTJgk3n9uP8mG1ttdv5yTAlP3b0ofcRj/n979x0W1ZU+cPw7M/QqIII0QSwoRhBUbInYiRpTTDGxYbKmWDaWXWOSzbpppmiKGoMmUXEtP5MYSzRKYmLXBIiIFQsqiiBFkd5n5vfHpYS10KZzPs8zD5eZO/e+XHE49z3nvEcmk2FjYaaXh6HV4hniMwQZMuIz4rmW33Kzf4euHyKrJAtnK2cGdXhMevL4Or3GpHfVg4t7TJAGYVeRy2U8GeoFiIkNTaW1Bk5GRgZt2tRdR8PMzAxnZ2cyMjIadIxVq1bRpUsX+vWrO8r+nXfe4bvvvmPPnj2MHTuWadOmsWzZsnse54MPPsDR0bHm4e3t3fgfyEQVlhey8/JOAB7v+HjtH/SEFtxNdep7qCwF1y7gGXrHy+OqUseb4lNb7AyHAxeyycwvw8nGnKFdxXo59XG1caW/Z3+g9oaiJaoeaD3Gfwzm1d0x536Colt6jEqPci7DlYOATGrg/I8nQ72QyeBw8k2u3CzSfXxGrtENnPnz5yOTye77OHfuXLMDKykpYePGjXfN3rz11lv079+fHj168NprrzFv3jwWLVp0z2O9/vrr5OXl1TxSU0WxtmrbL22nqKIIXwdfwtqGQfAEkJtBaixkJek7PN1TqyHua2k7dDJ3jJwFHgnywMZCweXsIuJTbus4QMOwMVbKQjwR4nXH7Cnh7sZ2lFaG3nJxS4tcnyqzKJODadKyJ090fALaBkHbYFCWQ0K0XmPTm+qu8A5DodWdmVBvZxsGdZZuINYeTdFhYKah0Q2cuXPnkpSUdN9H+/btcXd3Jyurbln7yspKcnJyGjR2ZvPmzRQXFzNp0qR69w0LC+P69euUlZXd9XVLS0scHBzqPARQqVX837n/A2B8l/HIZXJpAc5OEdIOLXGw8aW9cPM8WNhB8HN33cXO0oxHuktjyDbEXtVldAbhys0i9p6X/m8/F3b37iljMmvWLFJSUrR+noHeA2lj04ac0hxiUmK0fj5Ds+XiFlRqFSFtQmoX1gx7Wfoa9w0oW9iK68oKOL5B2v7L4OL/FdnPF4DNx65TUNrCrlEzNbqB4+rqSkBAwH0fFhYW9O3bl9zcXI4dO1bz3r1796JSqQgLC6v3PKtWrWLMmDG4utY/FTcxMREnJycsLVvugM+mOJx2mKv5V7E3t2eM/5jaF0KnSF8TN0BZgX6C05fYFdLXHhPA6t5TMyf0kQpx/XTyBhl5pbqIzGCsPZqCWg2DOrvi73r3kg/Cnczl5jwb8CwA68+ub1Hdm2XKMjad3wTAuIC/lPPo9gTYtoGCdKmSb0tyZhsUZYHdX24q7+LBjq3p0MaOwrJKNh8TY3EaQ2tjcLp06UJERARTp04lLi6OI0eOMGPGDMaNG1czgyotLY2AgADi4uLqvDc5OZmDBw/yt7/97Y7j7tixg2+++YbTp0+TnJxMVFQUCxcuZObMmdr6UUzWhiTp7uHxjo/X1OoApHLqLh2lwcbVtWBagpvJcPEXQAa9713OAOABL0d6+zlTqVLz399TdBKeIcgvreD7P6Uu3ucH3L1mlXBvT3Z8EiuFFUk5SRzLPFb/G0zErsu7yCnNwd3WnWHthtW+YGYJvaqGIVTfXLQEajX8XjVutNdUUJjfc1eZTMbkqizO2qMposhoI2i1Ds6GDRsICAhgyJAhjBw5kgEDBvDVV1/VvF5RUcH58+cpLq5b7Gn16tV4eXkxfPjwO45pbm7O8uXL6du3L8HBwaxcuZJPP/2UBQsWaPNHMTmXcy9zNP0oMmR176gA5HLoO13a/uPLllOnovoDttMIcPGvd/cXqv7Ab4i9RnF5pTYjMxjfxadSVK6kYxs7BnQQCyQ2ViurVjXTxNcntYybB7VazX/PSmNNxgeMv7Nycc/nQWEB1+Phegtp9F09CjdOgJm19PPX44kenthbmZFyq5gDF7J1EKBp0GoDx9nZmY0bN1JQUEBeXh6rV6+uU8/G19cXtVpNeHh4nfctXLiQa9euIZffGV5ERATHjx+noKCAwsJCEhMTeemll+66r3BvG89JhaXCvcPxtr/LrLKgcWDtDLnX4NxOHUenByW5tcW2qscF1GNoFzd8nG3IK6nghwTTXxRPqVITXTXQ8fkBfgY3FdtYjO8yHoC91/aSWmD6Ex5+v/E7ybnJWJtZ80SnJ+7cwa4NdJMGYBMbpdvg9OX35dLX4GfB1qXe3W0tzWpmb64+ckWbkZkU0SpogfLL8/nxklTnpfrD9g7m1tCrqovw6Bc6ikyPjq+HiiJpanj78Aa9RSGX8Xx/XwDWHL5i8qnjPWczuX67BCcb85qKzkLj+bfyp79Hf9Soawb5m7J1Z6U6N090fKJmdfU7VN9UnNkK+Td0FJme3LoE53dJ232mNfhtk/r6IpfBoYs3Sc5qYWMjm0g0cFqgrRe3UlJZQodWHejt3vveO/b6W1XqOA5S4+69n7FTKSFupbTd5+W7Tg2/l6d6emNvZcblm0XsO59V/xuMWPWd43NhPlgZUKVgYzShq1TzZMvFLRSWF+o5Gu25nHuZw2mHkSFjfMA9bqYAPILBpx+oKiH+G53Fpxd/fAmopYHFrTs2+G3ezjYM7eIGUJNJFe5PNHBaGKVKWWdq+H27Gezd4IGnpe3qlKopOr9b6oqzdqr9eRvI1tKMZ3tLU6VXHTbd1PHptDziruRgJpcxsY+vvsMxev08+uHn6EdRRRHbL5nu7KF1SVL2ZrDPYLwd6imw2qcqi3NsDVSY6MzE4pzaqeHV4xwbIbIqY/zDsTTySsSU8fqIBk4Lsy91H2mFaThaOjKq/aj639C3KoWa9CPcNtGaL9WDi0MjwcLmvrvezeR+vijkMo5eusWZ9DzNxmYg1hxJAWBU97a4O7bs1bA1QS6TM6GLlMXZkLQBpQkO5L9depsdl3YAMLHrxPrf0HkUOHpD8S2pmrgpOrYGKkvA/QHwfbDRb+/b3oXObvaUVCj5Viz4Wy/RwGlB1Go1X52UZrE93elprM2s63+TWyC0HwRqlWlO47x+DFIOgUxRO+aokTxbWfNw1Wraqw+naDA4w5CVX8qOE+kATOkvpoZryuj2o3GwcCC1IJX9qfv1HY7GfXf+O8qUZQS6BBLSJqT+NyjMoPdUafvoUtObvVlZDrFVs4j7zmhUV3g1mUxWM3vz60NXxCrj9RANnBZkf+p+knKSsDazbtgdVbW+M6SvCf+FUhPLUOz/QPoaNA4cvZp8mOoPnR9PpJGZb1rp9RUHLlOuVNGznRPB3q30HY7JsDG34ZnOzwCw4uQKkyr8V1pZWtMVPrHrxIbPuAuNBKtWcPOCNODYlJz+AQozwM4dAu8ym6yBHuvhiWcra7ILymqWTBHuTjRwWgi1Wk3UCWkK5nMBz+Fk5dTwN3cYAq4BUF5YewdiCq7/Ccl7pOzNQ/9o1qF6+DjRy9eJCqWaL/clayhA/cvML2V91XIUrw5t+IBIYxIZGVmzjp65uTlubm4MGzaM1atXo1KptHruSV0nYWNmw7mcc+y9tler59Klb89/y63SW3jYejDc9856Zvdk5Vh7Q3XgI9PJ4qiUcPgzaTvsRTCzaPKhLMzkTB/UAYAVBy6JLM59iAZOC3Hw+sGa7M3kwHuve3JXMhk89E9p++gyKDGRBSb3fyh9DXoWnNs3+3Czh3YC4P/iUknLLWn28QxB1P5LlFeq6OXrZNKF/SIiIrhx4wYpKSns3r2bQYMG8eqrrzJ69GgqK7VXxLGVVauaUg1fnvgSlVq7DSpdKK4oZvXp1QC8HPQy5vJ7V+m9q7CXTC+Lc+p7aY07q1ZN7gr/qydDvfBsZU1WQRn/FyeyOPciGjgtgFqt5ssTXwLSOjCNyt5UC3wC2nSFsjypkWPsUuP/kr2Zq5FD9uvQmr7tXShXqvhi70WNHFOfbuSVsLHqw3P20E4mXdjP0tISd3d3PD09CQkJ4Y033mD79u3s3r2b6OhorZ57cuBk7MztuHD7Ar9e/VWr59KFjec2klOag4+9T03V5kaxcoB+JpTFUVbUdoUPmHXfNe4aysJMzrRBUrV1kcW5N9HAaQEOpR3i7K2zWJtZExkY2bSDyOUw6E1p+48VUGjk5cIPaDZ7U23ucCmL892f17l6q0hjx9WHL/dJ2Zvefs709a+/2uod1GooL9LPQwPjWQYPHkxQUBBbtmxp9rHux9HSsaYuTtSJKKPO4hSUF7Dm9BoAXgl+5c5lGRqqtwllcY6vh9sp0qKi9axx1xhPhXrj4WhFZn4Z38abfkXspmjib59gLNRqNVGJ0tibZzo/g7OVc9MPFjAKPEIgPUHqT45YqKEodSw1HpJ/1Wj2plpPX2cGdnLlwIVslvx2kU+fDtbo8XUlPbek5kOzydmbimJY6KHhyBrojXSwsG32YQICAjh58qQGArq/iV0nsuHsBpJzk/nl6i9E+N57dWlDtu7sOvLL8/F39Odh34ebfqDqLM7e96QsTuDjIDfC4pIVpXBwkbT94FyN/E5WszCT88qgDry17TRR+y/xTC9vUYDzf4gMjok7nHaY07dOY6Wwanr2pppMBoP/JW3HfwN5Rrr+UnW6OFiz2Ztqc4ZJWZxtx9NIzjLOKrXL9yVTrlTRt71L07I3JkKtVuuka87BwoGJgdLMxqjEKKOsi5NbmluzqOa04Gkomtsg6f2SVHzz5gU4rd0smtb8uRry08DBC3pO0fjhn+7pRVtHKzLyS/nuT5HF+V8ig2PC/jpz6pnOz+BirYE/VP6DpZLq147CocUw+rPmH1OXUuPg0m9S9ubB5s2cupcg71YM6+rGnrOZfP7rBb54rgE1QAzI9dvFNR+Ws6saa01ibiNlUvTBvPEFG+8mKSkJPz/d1P6Z0GUC686u43LeZX5O+ZmR7Ufq5LyaEn0mmqKKIjo7dWZou6HNP6CVgzSjau+7Uhan2xPGlcUpK4TDn0rbA/8JZpYaP4WlmYJp4f68tf0MX+67xNM9RRbnr0QGx4TtvLyTUzdPSWNvukVq5qAyGQx5S9pO+C/kGNHyBGo17FkgbQc/C87a+8NVncXZefIGSTfytXYebVj2WzIVSjX9O7jQ268ZXZoymZSS18dDA1mXvXv3curUKcaOHdvsYzWEvYU9k7tKMxyjTkRRoTSeUvw3S26y8dxGAKYHT0cu09Cflt4vSlmcWxfh+DrNHFNX4lZCUTY4+UHwfdbhaqane3nXZHG+PnhZa+cxRqKBY6IKywv55M9PAHip+0u0ttbgFN92/cB/iLQwXvVUa2Nw8jsp82RmDQPna/VUXdo6MOqBtgAs+vm8Vs+lSSdSc/numJS9mdOc7I2RKSsrIyMjg7S0NBISEli4cCGPPvooo0ePZtKkSTqLY3yX8ThbOZOSn1LT3WMMVpxYQUllCd1cuhHuHa65A1s5wEPzpO1f35bWcjIGJblwZKm0PegNUDRyqnwjWJopmP9wAADL9ydz/Xax1s5lbEQDx0RFnYjiVuktfB18G1e1uKGqx+Kc3ATX/tD88TWtNA/2VGWeBv4TWtWz8J8GzB7WCTO5jL3nsvjlTIbWz9dcSpWaf207jVoNT/TwJLRdM7I3RiYmJoa2bdvi6+tLREQE+/btY+nSpWzfvh2FQncpfzsLO+b2lAa+rzy5khuFN3R27qY6lX2K785/B8Ds0NmaH7PU+0VoEwglOVJ3lTH47R0ozQXXLtBN+xnAMUEe9PZzprRCxcJdSVo/n7EQDRwTlHw7mQ1J0oq183vPx0LR9KqZ9+QZAiFVd7Y7ZknrrBiy/R9CYSY4+9dWStWyDm3sePEhaRDzf348Q1GZ9grGacLGuGucSsvD3sqM10d20Xc4OhMdHY1arUatVlNRUUFWVhZ79uxhypQpyOW6/4h8pP0jhLQJoaSyhI/jP9b5+RujUlXJu3+8ixo1o9uPpnfb3po/icIMRlbNRPpzDaQlaP4cmpQaJw0uBhi1WCfjhmQyGW+PCUQhl7HrVAaHL97U+jmNgWjgmBi1Ws2HcR+iVCsZ7D2Y/p79tXeyoW+DTWvIToLfDbj4X+YZiF0pbY/8WCuD/e5l5uCOeDlZk55XypLfDLf4383CMhbFnAPgH8M742qvu2sk1CWTyXizz5soZAp+vfYrh9MO6zuke9p0bhNJOUnYW9jXZJ60wrc/dH8GUMOuf4CWl9BoMmWFdMOHGoIngO8AnZ26S1sHJvZpB8CCH09TXmmg10iHRAPHxPxy9RdiM2KxVFjyz17/1O7JbJxhRFUtnAMfQ44BDnBTq2HXP0GthC6PQAcNzO5oBGsLBe8+2g2AVYevcDbdMAccf7j7HPmllQR6ODCh6kNS0J9OTp1qlnBYGLuQMmWZniO6U2ZRJsuOSzc2s0Nna3ac390Mewcs7CHtmOEOOP79C8g6A9bOUrw6NntYJ1xsLbiUXcTaoyk6P7+hEQ0cE1JcUcyieCmV+3y35/Gyb/rq2A3W/WnwGwiVpfDTXI1UkNWoU9/D1SPSwOIRH+glhEEBbRj5gDtKlZo3t51CpTKsa/RnSg6bj10H4N3HuqGQm+6SDMZkWvA02li3IbUgtWZtJ0PyUfxHFFcWE+QaxNiOOphpZu8Og16Xtn/9j+ENOM65Avs/krZHvA+2uq8f5WhtzmsR0oDjz3+9QFZ+qc5jMCSigWNCliQsIbM4E087T57v9rxuTiqTSbVwFJZwaS+c/kE3522IktvwS9Vg6If+oZOBxffy79GB2FooOH4tl00GVFa9UqniX9tOAzCulzchPk1Yp0zQCltz25os7DcnvyE133B+bw5eP8ieq3tQyBS81ectzU0Lr0/vF6U18UpypEaOoVCrpRu8yhLwfVBaAkZPngz1Isi7FUXlShb8eAa1od106pBo4JiI3679VlOH4s2wN7Eys9LdyV38a1cbj5lvGKuNq9Wwbbo0sNilA/Sbqddw3B2tmDu8MwAf7k4ymDurJb9d5FxGAa1szJlXdecnGI4RviPo07YP5apy5h+abxC1cQrLC1kYK3VNT+w6kc7OnXV3coU5jFwsbSeshTPbdHfu+zmzRSogqrCQbvj0uDCtXC7jvUe7YSaXsft0Rs2CuS2RaOCYgLTCNN46Ik2BjgyM5EGvB3UfRP+/Q+vOUmGrbdP1Pwjw9+Vw/ifpA2fsKp0OLL6XSX3b0c3TgfzSSmZsPE6FUr/XaN/5LJbtTQbg7TGBONtqYbad0CwymYwFfRdgb2HPyZsn+eTYJ3qNR6VW8ebhN0krTMPD1oNXgl7RfRC+/aH/LGl7+wy4dUn3MfxVzhXYOUfafnAutO6o33iAB7wc+ecIqeH59o6zBjv2T9tEA8fIVSgrmHdgHgXlBXRv3Z2/h/xdP4GYWcLjK6SuqvM/1S4wpw+p8fBrVcXiEQvBI1h/sfyFmULOknE9sLM0Iy4lhw92ndNbLGm5Jcz+NhGAiX3a8Wiwp95iEe7Py96LhQOkjMmGpA3EpMToLZbVp1ezN3Uv5nJzPgn/BBsNLYnRaIPfAp++UF4A30+WFrXUh/Ii+HaCVPPGMxQGzNZPHHcx9cH2DOrsSnmlihkbEwy+TIU2iAaOkVuSsISTN09ib2HPxwM/xlyuvYqZ9fIMqV2bav9COK+HD+LiHPg+UqqyHPgE9Pqb7mO4D39XOz55OgiA1UeusD1R9wuWlleqmLYhgdziCrp7OfKv0S2n5o2xCvcO54VuLwCw4MgCruTpfomUo2lHa2ZNvRH2Bt1ad9N5DDUUZvDkarBxgYxT8PPruo9BrZYySJmnwbYNPLPeIDLF1eRyGZ88HYy7gxWXbxZVFfFsWeNxRAPHiB1IPcDas2sBeLf/u3jaGcBdeI/x0GuqtL1lKtxM1t25VSrY+hLkX5cK+j2yRK994fcyItCd6YP8AZj/wynOZeg2fbxwVxInUnNxtDZn+XMhWJqJxfmMwYweM+jp1pPiymLm7J9DcYXuSvKnFaYx79A8VGoVYzuO5clOT+rs3Pfk4AFPfA3IpMJ6pzbr9vxHl0ljb+Rm8PRaKR4D42xrwdJne6CQy9h6PI3vq2ZLthSigWOkkm4l8cbhNwBpFeIhPkP0HNFfjFgopY/L8mHTc1Cqoz/gBxfBxV/AzEr6wLFy0M15m2DOsM482LE1JRVKXlp3jLwS3Qwe3XEineiq+hifPh2Et7OeuhiERjOTm7Fo4CJaW7cmOTeZd/54B5Va++O4SitLmb1vNnlleQS6BPJ6mB6yJffSYYg0QxLgx79L2RxduLSvths84kNpfT4D1dvPuWZduX9vP82xqwY2vV6LRAPHCJ25dYYXfnmB/PJ8glyDmB1qOP2+AJhZwFNrwd4Dbp6HLS9CpZYLlR34WOoWA3j4I3B/QLvnayaFXMbScT3wbGXN1VvFzNiYQEm5Uqvn/PlMBnO/OwHAtHB/hnRx0+r5jEVkZCQymQyZTIa5uTlubm4MGzaM1atXo6oaLJ+Tk8PMmTPp3Lkz1tbW+Pj48Pe//528vDydxtraujUfP/Qxcpmcny7/xNu/v41Spb3fmzJlGa8dfI2knCScLJ34LPwzLBWG0w0DQPjr0tTsiiJY+wikH9fu+bLPw+YpoFZJq4QbWDf43bwy0J+BnVwprVAxaVUc8SnabeTklVQwa9NxbuSVaPU89RENHCNzKvsUU3+eSkF5AUGuQUQNjdLOWlPNZe8Gz6yTZjFd2A0bntJOJketlha22/e+9P3gf0FopObPowVOthasnBiKlbmcQxdvMv6bP8gt1s6aXtuOpzFtQwLlShUPd3NvUSuFN0RERAQ3btwgJSWF3bt3M2jQIF599VVGjx5NZWUl6enppKens3jxYk6fPk10dDQxMTG88MILOo+1l3sv3u3/LnKZnC0Xt0jTx1WazwAWlBfw8p6XawYVLxq4iLZ2bTV+nmaTK6TxL549pRIVax+V1oPShmuxsHqEdB6PHjDqU4PsBv9fcrmMFRNC6d/BhaJyJZNXxxF7+ZZWzpWWW8JTK46yLTGdmRuP63Xcj2jgGJET2Sd4cc+LFFQUENImhJXDVmJvYa/vsO7Nqyc89y1Y2MGVA7BmJORrcHVktRp+fhMOVU2dHf5+bT0eI9HN05ENfwvDwcqMhGu5PLXid43f9WyIvcrs7xJRqtSMDfFi2bM9MFOI//p/ZWlpibu7O56enoSEhPDGG2+wfft2du/eTXR0NN26deOHH37gkUcewd/fn8GDB/P++++zY8cOKit1PztljP8YFj20CDO5GTEpMczeN1ujyzncLLnJlJgp/Jn5J7bmtqwYuoKwtmEaO77GWbeCSdvApx+U5cG6xyHliGbPkbQT/jtGatx49YLxP4C5DuuNNZO1hYJVk3vxYMfWFJcriVwTz9FLml2U83RaHo8vP8KFzELcHCx5+9FAza8u3wjiU85IHLx+kJf2vERhRSGhbqFEDY3C1txW32HVz38wRP4kzTLIPAWrhkP2heYfV1khVQ79Y7n0/cjF0E83q4RrWmg7Z75/uR9uDpZczCpk7JdHSc4q1Mixvzp4iTe3nkatlqaDL3qyu84aN2q1muKKYr08NHHXOHjwYIKCgtiyZctdX8/Ly8PBwQEzM7Nmn6sphvsOZ+mgpVgqLDlw/QDTf51OUUVRs4+bmp/KxF0TOX/7PC5WLqwZsUY7q4RrmqU9TNgM7cOhvBDWj4WLv2rm2PGr4LuJ0pI0nSJg0o96WYqhuazMFXw9qScDO7lSUqHk+eh49p7L1Mix953P4pmVv5NVUEZnN3u2TutPoIejRo7dVDK1lvJH77//Pj/99BOJiYlYWFiQm5tb73vUajULFizg66+/Jjc3l/79+xMVFUXHjrWFk6r7wnfs2IFcLmfs2LEsWbIEOzu7BseWn5+Po6NjzQeUNq08sZLM4kymPjC1SendrOIsPo7/mJ9Tfgagt3tvlg1epr/6E02Vc0X6wMm5BNZO0nTyLo+CvAl/bFOOSCsKZ50FZDBmKYRM0njIunb9djGTVsdxObsIJxtz/jMmkEe6eyBvwtpQqTnFvP9TEjFnMgB4eaA/r0V01trdVGlpKVeuXMHPzw8rK+mutriimLCN+rnrj30utsH/RyIjI8nNzWXbtm13vDZu3DhOnjzJ2bNn6zx/8+ZNQkNDmTBhAu+///49j32366Jp8RnxzPhtBsWVxbjZuDG351wifCMa/W9dqarkx0s/siRhCTmlOXjZefHVsK/wdtDfEidNUlEqNUYu/iJ9Hzwehr4Ndq6NP1bJbWl9qdgo6fuQyVK3lEI/jVpNKa1QMm1DAnvPZQHweA9P3hzVhdZ2jR9fpVSpWf/HVd7ZeRalSk3/Di5ETQjFwUo7JUsa8/dba7dy5eXlPPXUU7zySsMrXX788ccsXbqUFStWEBsbi62tLSNGjKC0tLaI0/jx4zlz5gx79uxh586dHDx4kBdffFEbP0KzFZYXsubMGr6/8D0jt47kvT/eI6Moo0HvVaqUbEzayKPbHuXnlJ+Ry+RM6DKB5UOWG1/jBsDZD17YU9tP/n0kfB0u3WE1tI1dkAE/TIXokVLjxtoZnlpjEo0bAC8nGza/3I8gL0duF1fw6qZExiw/zOGLDU8jF5dXsvjn8wz59AAxZzJQyGXMi+jM/IcD9JoqNlZqtfqO65afn8+oUaPo2rUr//nPf/QT2F/0cu/FN8O/wdPOk8ziTOYdnEdkTCRnb52t/81IP+O+a/t48scnWXB0ATmlOQQ4B7Bu5Drja9yA1G30zPrasXiJG2BZKMR+BcoGdieWF0ld30uCahs34a9LpSeMvHEDUiYnakIIk/u2QyaDrcfTGLx4PxtjrzV4MWClSs32xDSGf3aABT+eQalS80SIJ2sie2utcdNYWsvgVIuOjmbWrFn1ZnDUajUeHh7MnTuXf/xDmvaXl5eHm5sb0dHRjBs3jqSkJLp27Up8fDw9e/YEICYmhpEjR3L9+nU8PBpWh0CXGZxjmcf4MvFL4jKkQW/mcnPGdhzLw34P42nniauNa81Cdfnl+RzPPM6fmX9y6PohLuVJJcgfaP0Ab/V5iy4uJlCQrbwYjiyRllIoL5Ce8+knrRXlFgiO3nWzOgUZ0oDBa39Awn+r3iODnlOkaqY2znr5MbSptELJqsNXiNp/icKq6qMPdmzN1Afb09XD4Y67rLJKJUk3Cjh29TZfH7xMRtU6V/38XVjwSCCd3bU/TutumQq1Wk1JpX5mUVibWTe4QXe/DE737t3x8fFh586dABQUFDBixAhsbGzYuXNnvVkZXWRwas5VWcraM2tZdXoVJZUlyJAR4RdBb/fedGvdDf9W/jWFQEsqS7icd5lLuZf44cIPJGQlAOBo6cjUB6YyLmCc4c2WaorUeNg1F25IswdxewAeGAveYdIgYXPr2n0rSiDnMlw5CIc+hSIpu0GbrjBkAXSO0H38OnAiNZc3tp7iTNVyDkHerXike1vC/Fzo0tb+ji7tvOIKDiVns+TXi1ys6kp3tDZn5uAOvDDAT+s3Uo35+20wDZzLly/j7+/P8ePHCQ4Ornl+4MCBBAcHs2TJElavXs3cuXO5fbt2McfKykqsrKz4/vvvefzxx+967LKyMsrKagfg5efn4+3trZMGTrX4jHiWJy7nWOaxOs9byC3wsPPAXGFO8u1k1NT+c9ib2/NqyKs82elJFHITK8ZWdAuOfAZxX0v92tUUltLinY5e0nTM3Kt13+cZKo238QzRbbx6cKuwjGV7k9kQe5UKZe3vhbOtBR3b2OHZypqLWYWcy8iv87q3szVvjuzKiEA3nWVtdPmHXNPu1cDZu3cvQ4YMYfXq1UyZMoX8/HxGjBiBpaUlu3btwsam/kyqPq5LRlEGnx37jF1XdtV53lJhScdWHcktyyWtMK3OZ42lwpIJXSbw/APP42BhuPWjmkSlhGNr4Ld3pSUVqsnNwL27VC/r1iXI+5/V2p18YdCb0G2sNFPLhFUqVaz74yqf/HKh5qYKwM7SjJ6+TlibK7iWU0xqTjH5pbWvO1iZMfXB9kT298VeR1mbxjRwDCbXlpEhdd24udWtzeHm5lbzWkZGBm3atKnzupmZGc7OzjX73M0HH3zA22+/reGIG6eXey/WjFhDXEYc686uIzk3mYyiDMpV5aTkp9Ts186hHT3dehLqFsoAzwE4WTnpL2htsnWB4e9Bn2nS3VLKIenuSVkmdT9lVafXZVJmx7s3+A2ELmOaNm7HCLnYWfKfMYFM6e/LF3uTiU/J4WpOMTlF5cReqVvHwtnWgu5ejgzo0JoJfdphZW7aH8iaVlZWRkZGBkqlkszMTGJiYvjggw8YPXo0kyZNIj8/n+HDh1NcXMz69evJz88nP1+643V1dUWhMJzr7W7rzkcPfcT4LuP57dpvnLl5hjO3zlBYUcjpW6dr9nOydMK/lT+BLoFM6DoBd1t3PUatRXKFVKum62NwYhOkxkpZ4cIMSE+ou69VK2jdCYLGSV3fCsPoatE2M4WcKf39GPlAW7YdTyP2Sg7xV3IoKKtk//nsO/Z3c7Dk2d4+PD/Az2C6o+6mUQ2c+fPn89FHH913n6SkJAICApoVlKa9/vrrzJkzp+b76gyOrslkMsLahtVMt6xQVZBRlEFaYRpFFUV0b90dV5smDIQzZg4eMGqxtK1SQu416W4q9yo4t5cyNgZckVgX2rnYsugpaf2qknIll7ILuZBZQNrtEtq72tHdyxEvp4Z3yQh3iomJoW3btpiZmeHk5ERQUBBLly5l8uTJyOVyEhISiI2NBaBDhw513nvlyhV8fX31EPX9dXftTnfX7oC0CvjV/KtcuH2hpmHjYm18s4CaxbZ11UzLGdK4v7xUqaFTWQouHcGlg9Tl3YL/H7k5WPHSQH9eGuiPUqUm6UY+f6bkoAa8nWzwcbHBy8kaGwuDyY3cV6OinDt3LpGRkffdp3379k0KxN1dunvIzMykbdva2UaZmZk1XVbu7u5kZWXVeV9lZSU5OTk1778bS0tLLC0Nrz/ZXG6Ot7033vZGOJBPG+QKaTCys5++IzFY1hYKunk60s1Tv9MvTUl0dDTR0dH33Sc8PNyoFyqUy+T4Ofrh5yj+bwFSI6aVj/QQ7kohlxn9Z02jGjiurq64umonw+Dn54e7uzu//fZbTYMmPz+f2NjYmplYffv2JTc3l2PHjhEaGgpI/eQqlYqwMAMuQiUIgiAIgk5pbTDDtWvXSExM5Nq1ayiVShITE0lMTKSwsLaAWUBAAFu3bgWk7ptZs2bx3nvv8eOPP3Lq1CkmTZqEh4cHjz32GABdunQhIiKCqVOnEhcXx5EjR5gxYwbjxo1r8AwqQRAEQRBMn9Y60v7973+zdu3amu979OgBwL59+wgPDwfg/PnzdRarmzdvHkVFRbz44ovk5uYyYMAAYmJi6sw+2LBhAzNmzGDIkCE1hf6WLl2qrR9DEARBEAQjpPVp4oZIl3VwBKElMOZp4tokrosgaJZBVDIWBEEQBEHQF9HAEQRBY1Qqlb5DMCjiegiC/hjHZHZBEAyahYUFcrmc9PR0XF1dsbCwaNF1edRqNeXl5WRnZyOXy7GwsNB3SILQ4ogGjiAIzSaXy/Hz8+PGjRukp6frOxyDYWNjg4+PD/IWUn1bEAyJaOAIgqARFhYW+Pj4UFlZiVKp1Hc4eqdQKDAzM2vRmSxB0CfRwBEEQWNkMhnm5uaYmxvu+jSCILQMIm8qCIIgCILJEQ0cQRAEQRBMjmjgCIIgCIJgclrkGJzq4s35+fl6jkQQBEEQhIaq/rvdkEUYWmQDp6CgAABvb289RyIIgiAIQmMVFBTg6Oh4331a5FpUKpWK9PR07O3txRROpBaxt7c3qampYm0uLRLXWTfEddYNcZ11Q1znutRqNQUFBXh4eNRbX6pFZnDkcjleXl76DsPgODg4iP9AOiCus26I66wb4jrrhrjOterL3FQTg4wFQRAEQTA5ooEjCIIgCILJEQ0cAUtLSxYsWIClpaW+QzFp4jrrhrjOuiGus26I69x0LXKQsSAIgiAIpk1kcARBEARBMDmigSMIgiAIgskRDRxBEARBEEyOaOAIgiAIgmByRANHuKuysjKCg4ORyWQkJibqOxyTkpKSwgsvvICfnx/W1tb4+/uzYMECysvL9R2aSVi+fDm+vr5YWVkRFhZGXFycvkMyKR988AG9evXC3t6eNm3a8Nhjj3H+/Hl9h2XyPvzwQ2QyGbNmzdJ3KEZDNHCEu5o3bx4eHh76DsMknTt3DpVKxcqVKzlz5gyfffYZK1as4I033tB3aEbv22+/Zc6cOSxYsICEhASCgoIYMWIEWVlZ+g7NZBw4cIDp06fzxx9/sGfPHioqKhg+fDhFRUX6Ds1kxcfHs3LlSrp3767vUIyKmCYu3GH37t3MmTOHH374gcDAQI4fP05wcLC+wzJpixYtIioqisuXL+s7FKMWFhZGr169+OKLLwBp3Tlvb29mzpzJ/Pnz9RydacrOzqZNmzYcOHCAhx56SN/hmJzCwkJCQkL48ssvee+99wgODubzzz/Xd1hGQWRwhDoyMzOZOnUq69atw8bGRt/htBh5eXk4OzvrOwyjVl5ezrFjxxg6dGjNc3K5nKFDh/L777/rMTLTlpeXByB+f7Vk+vTpjBo1qs7vtdAwLXKxTeHu1Go1kZGRvPzyy/Ts2ZOUlBR9h9QiJCcns2zZMhYvXqzvUIzazZs3USqVuLm51Xnezc2Nc+fO6Skq06ZSqZg1axb9+/enW7du+g7H5GzatImEhATi4+P1HYpREhmcFmD+/PnIZLL7Ps6dO8eyZcsoKCjg9ddf13fIRqmh1/mv0tLSiIiI4KmnnmLq1Kl6ilwQmmb69OmcPn2aTZs26TsUk5Oamsqrr77Khg0bsLKy0nc4RkmMwWkBsrOzuXXr1n33ad++PU8//TQ7duxAJpPVPK9UKlEoFIwfP561a9dqO1Sj1tDrbGFhAUB6ejrh4eH06dOH6Oho5HJxv9Ec5eXl2NjYsHnzZh577LGa5ydPnkxubi7bt2/XX3AmaMaMGWzfvp2DBw/i5+en73BMzrZt23j88cdRKBQ1zymVSmQyGXK5nLKysjqvCXcSDRyhxrVr18jPz6/5Pj09nREjRrB582bCwsLw8vLSY3SmJS0tjUGDBhEaGsr69evFB5WGhIWF0bt3b5YtWwZIXSg+Pj7MmDFDDDLWELVazcyZM9m6dSv79++nY8eO+g7JJBUUFHD16tU6z02ZMoWAgABee+010SXYAGIMjlDDx8enzvd2dnYA+Pv7i8aNBqWlpREeHk67du1YvHgx2dnZNa+5u7vrMTLjN2fOHCZPnkzPnj3p3bs3n3/+OUVFRUyZMkXfoZmM6dOns3HjRrZv3469vT0ZGRkAODo6Ym1trefoTIe9vf0djRhbW1tcXFxE46aBRANHEHRsz549JCcnk5ycfEfDUSRUm+eZZ54hOzubf//732RkZBAcHExMTMwdA4+FpouKigIgPDy8zvNr1qwhMjJS9wEJwj2ILipBEARBEEyOGNUoCIIgCILJEQ0cQRAEQRBMjmjgCIIgCIJgckQDRxAEQRAEkyMaOIIgCIIgmBzRwBEEQRAEweSIBo4gCIIgCCZHNHAEQRAEQTA5ooEjCIIgCILJEQ0cQRAEQRBMjmjgCIIgCIJgckQDRxAEQRAEk/P/cc2VyxN8A+YAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "n_grid=100 # 为画图较平滑，取较大的n_grid\n",
    "x=np.linspace(-5,5,n_grid)\n",
    "h=x[1]-x[0]\n",
    "D=(-np.eye(n_grid)+np.diagflat(np.ones(n_grid-1),1))/h\n",
    "D2=D.dot(-D.T)\n",
    "D2[-1,-1]=D2[0,0]\n",
    "\n",
    "y=np.sin(x)\n",
    "plt.plot(x,y, label=\"f\")\n",
    "plt.plot(x[:-1],D.dot(y)[:-1], label=\"D\")\n",
    "plt.plot(x[1:-1],D2.dot(y)[1:-1], label=\"D2\")\n",
    "plt.legend()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到结果符合三角函数的性质，将$D$和$D^2$两个算符作用到函数上，确实得到了它的一阶而二阶导数\n",
    "\n",
    "我们可以先不管势能项，试试看只有动能项的时候（自由电子）KS方程求解出来会是什么样子\n",
    "\n",
    "注：在量子力学中，所有可观测物理量 $A$ 都有一个对应的算符 $\\hat{A}$，代表作用于波函数的一种运算。\n",
    "\n",
    "$$\\hat{A} \\psi = a \\psi$$\n",
    "\n",
    "$a$ 即算符本征值，$\\psi$ 是算符 $A$ 的本征态。对这个本征态波函数 $\\psi$ 测量可观测物理量 $A$，得到的值就是 $a$\n",
    "\n",
    "类似的，如果想要求出体系的总能量，需要求解哈密顿算符的本征值。求出本征值的同时，本征态也会求出，对应特定状态下的波函数（我们预先并不知道体系的波函数长什么样子）。\n",
    "\n",
    "- 不考虑电子相互作用时，体系的哈密顿量:\n",
    "\n",
    "$$\\hat{H} = \\hat{T} = - \\frac{1}{2} \\frac{d^2}{dx^2}$$\n",
    "\n",
    "- 直接调用现成的线性代数求解器:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n"
     ]
    }
   ],
   "source": [
    "eig_non, psi_non=np.linalg.eigh(-D2/2)\n",
    "print(len(eig_non))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 求出本征值和对应的本征态后，绘图观察"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(5):\n",
    "    plt.plot(x,psi_non[:,i], label=f\"{eig_non[i]:.4f}\")\n",
    "    plt.legend(loc=1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "所绘制的曲线反映的就是波函数（能量的本征态）的信息，而图例中写的则是对应的能量（能量的本征值）大小\n",
    "\n",
    "后面的求解原理类似，只不过逐渐追加哈密顿算符的组分，都直接调用成熟的线性代数求解器完成本征值问题求解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 外势\n",
    "X=np.diagflat(x*x)\n",
    "print(X)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 引入交换关联势和库伦势\n",
    "\n",
    "多电子体系最麻烦的地方在于电子之间的相互作用，上面我们考虑了体系整体动能、外部势能，但是还没有考虑到电子之间的相互作用能。**我们需要先得到电子的信息——电荷密度分布，然后才能导出库仑势和交换关联势**\n",
    "\n",
    "先定义几个常用的函数\n",
    "- 积分\n",
    "- 波函数必须满足归一化要求:\n",
    "\n",
    "$$\\int \\lvert \\psi \\rvert ^2 dx = 1$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 积分\n",
    "def integral(x,y,axis=0):\n",
    "    dx=x[1]-x[0]\n",
    "    return np.sum(y*dx, axis=axis)\n",
    "\n",
    "# 归一化的波函数\n",
    "def get_np(psi, x):\n",
    "    # normalization\n",
    "    I=integral(x,psi**2,axis=0)\n",
    "    normed_psi=psi/np.sqrt(I)[None, :]\n",
    "    return normed_psi"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来可以求出电荷密度\n",
    "\n",
    "- 定义 $f_n$ 为占据数, 电荷密度 $n(x)$ 可以写成:\n",
    "\n",
    "$$n(x)=\\sum_n f_n \\lvert \\psi(x) \\rvert ^2 $$\n",
    "\n",
    "- 注意:\n",
    "    - DFT框架下，求解的都是最低能量态（基态）\n",
    "    - 每个态最多只能容纳两个电子，自旋分别向上和向下\n",
    "    - 电子从低能级到高能级依次填充"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_nx(num_electron, psi, x):\n",
    "    normed_psi=get_np(psi,x) # 调用上面的函数将波函数归一化\n",
    "\n",
    "    # occupation num\n",
    "    fn=[2 for _ in range(num_electron//2)]\n",
    "    if num_electron % 2:\n",
    "        fn.append(1)\n",
    "\n",
    "    # density\n",
    "    nx=np.zeros_like(normed_psi[:,0])\n",
    "    for ne, psi in zip(fn,normed_psi.T):\n",
    "        nx += ne*(psi**2)\n",
    "    return nx"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "准备工作做好后，试试看随便给定n个电子，体系的电荷密度长什么样？以及外部势对此有什么影响？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of electrons\n",
    "num_electron=17\n",
    "\n",
    "# 没有外势\n",
    "plt.plot(get_nx(num_electron,psi_non, x), label=\"non\")\n",
    "# 外势为x^2\n",
    "plt.plot(get_nx(num_electron,psi_harm, x), label=\"harm\")\n",
    "# 外势为无限深势阱\n",
    "plt.plot(get_nx(num_electron,psi_well, x), label=\"well\")\n",
    "# 打上标签\n",
    "plt.legend(loc=1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "有了电荷密度，就可以根据理论公式，写出电荷密度泛函项\n",
    "\n",
    "#### 库仑势的贡献（简化版）\n",
    "\n",
    "$$ E_{Ha}=\\frac{1}{2}\\iint \\frac{n(x)n(x')}{\\sqrt{(x-x')^2+\\varepsilon}}dxdx'$$\n",
    "\n",
    "其中 $\\varepsilon$ 是一个小正数\n",
    "\n",
    "- 库伦势:\n",
    "$$ v_{Ha}=\\int \\frac{n(x')}{\\sqrt{(x-x')^2+\\varepsilon}}dx'$$\n",
    "\n",
    "- 统一写成矩阵形式\n",
    "$$E_{Ha} = \\frac{1}{2} \\frac{n_in_jh^2}{\\sqrt{(x_{i}-x_{j})^2+\\varepsilon}}$$\n",
    "$$v_{Ha, i} = \\frac{n_jh}{\\sqrt{(x_{i}-x_{j})^2+\\varepsilon}}$$\n",
    "\n",
    "#### 交换能的贡献（忽略关联势）\n",
    "\n",
    "$$ E_X^{LDA}[n] = -\\frac{3}{4} \\left(\\frac{3}{\\pi}\\right)^{1/3} \\int n^{4/3} dx$$\n",
    "\n",
    "- 交换势是能量的导数:\n",
    "\n",
    "$$ v_X^{LDA}[n] = \\frac{\\partial E_X^{LDA}}{\\partial n} = - \\left(\\frac{3}{\\pi}\\right)^{1/3} n^{1/3}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_exchange(nx,x):\n",
    "    energy=-3./4.*(3./np.pi)**(1./3.)*integral(x,nx**(4./3.))\n",
    "    potential=-(3./np.pi)**(1./3.)*nx**(1./3.)\n",
    "    return energy, potential\n",
    "\n",
    "def get_hatree(nx,x, eps=1e-1):\n",
    "    h=x[1]-x[0]\n",
    "    energy=np.sum(nx[None,:]*nx[:,None]*h**2/np.sqrt((x[None,:]-x[:,None])**2+eps)/2)\n",
    "    potential=np.sum(nx[None,:]*h/np.sqrt((x[None,:]-x[:,None])**2+eps),axis=-1)\n",
    "    return energy, potential"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 激动人心的最终boss\n",
    "\n",
    "将上述所有项目加到一起，开始自洽求解电荷密度、本征态等\n",
    "\n",
    "流程：\n",
    "1. 随便猜一个电荷密度分布\n",
    "2. 计算泛函项（交换关联势和库仑势）\n",
    "3. 把它们加入哈密顿量\n",
    "4. 求出哈密顿量的本征值和对应的本征态\n",
    "5. 根据本征态积分出电荷密度，与猜测的对比，如果差别很大，重复这个流程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_log(i,log):\n",
    "    print(f\"step: {i:<5} energy: {log['energy'][-1]:<10.4f} energy_diff: {log['energy_diff'][-1]:.10f}\")\n",
    "\n",
    "max_iter=1000\n",
    "energy_tolerance=1e-5\n",
    "log={\"energy\":[float(\"inf\")], \"energy_diff\":[float(\"inf\")]}\n",
    "\n",
    "# 1\n",
    "nx=np.zeros(n_grid)\n",
    "\n",
    "for i in range(max_iter):\n",
    "    # 2\n",
    "    ex_energy, ex_potential=get_exchange(nx,x)\n",
    "    ha_energy, ha_potential=get_hatree(nx,x)\n",
    "    \n",
    "    # 3\n",
    "    H=-D2/2+np.diagflat(ex_potential+ha_potential+x*x)\n",
    "    # 4\n",
    "    energy, psi= np.linalg.eigh(H)\n",
    "    \n",
    "    # 5\n",
    "    ## 记录\n",
    "    log[\"energy\"].append(energy[0])\n",
    "    energy_diff=energy[0]-log[\"energy\"][-2]\n",
    "    log[\"energy_diff\"].append(energy_diff)\n",
    "    print_log(i,log)\n",
    "    ## 检查是否收敛\n",
    "    if abs(energy_diff) < energy_tolerance:\n",
    "        print(\"converged!\")\n",
    "        break\n",
    "    ## 更新电荷密度分布\n",
    "    nx=get_nx(num_electron,psi,x)\n",
    "    ## 也可以用beta参数混合新旧电荷密度，加速收敛\n",
    "    # beta=0.1\n",
    "    # nx=beta*get_nx(num_electron,psi,x) + (1-beta)*nx\n",
    "    \n",
    "else:\n",
    "    print(\"not converged\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "求解完成！让我们看看这个简单体系的本征值和本征态长什么样？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(5):\n",
    "    plt.plot(x,psi[:,i], label=f\"{energy[i]:.4f}\")\n",
    "    plt.legend(loc=1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "看起来还不错，那么电荷密度分布呢？与前面只有动能项的那个比较一下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['font.sans-serif'] = 'SimHei'\n",
    "plt.plot(get_nx(num_electron,psi_non,x), label=\"动\")\n",
    "plt.plot(get_nx(num_electron,psi_harm,x), label=\"动+外势\")\n",
    "plt.plot(nx, label=\"动+外势+交换+库伦\")\n",
    "plt.legend()"
   ]
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}