{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "8bc4e00f", "metadata": { "execution": { "iopub.execute_input": "2025-12-15T08:10:40.108797Z", "iopub.status.busy": "2025-12-15T08:10:40.108692Z", "iopub.status.idle": "2025-12-15T08:10:40.113775Z", "shell.execute_reply": "2025-12-15T08:10:40.113120Z" } }, "outputs": [], "source": [ "# Colab 环境检测与依赖安装\n", "try:\n", " import google.colab # type: ignore\n", " IN_COLAB = True\n", "except ImportError:\n", " IN_COLAB = False\n", "\n", "if IN_COLAB:\n", " !pip install -q git+https://github.com/bud-primordium/AtomSCF.git@main\n", " !pip install -q git+https://github.com/bud-primordium/AtomPPGen.git@main\n" ] }, { "cell_type": "markdown", "id": "05baf8df", "metadata": {}, "source": [ "# AtomPPGen 教程 04:Kleinman-Bylander 可分离形式\n", "\n", "本 Notebook 聚焦如何把半局域通道转换成 Kleinman-Bylander (KB) 可分离势。你将会:\n", "- 在 Colab 或本地快速配置 AtomSCF + AtomPPGen 运行环境;\n", "- 理解从 $O(N_{\\mathrm{PW}}^2)$ 降到 $O(N_{\\mathrm{PW}})$ 的效率提升逻辑;\n", "- 推导并实操 $\\beta_l$ 与 $D_l$ 的构造方法;\n", "- 使用 `kb_transform` 获得局域势、投影子及诊断信息,并完成可视化。\n", "\n", "建议先完成 03-potential-inversion 教程,确保半局域势已稳定可用。\n" ] }, { "cell_type": "markdown", "id": "8575b40d", "metadata": {}, "source": [ "## 为什么需要 KB 可分离形式?\n", "\n", "半局域非局域势写作\n", "$$\n", "V_{\\mathrm{NL}}^{\\mathrm{SL}} = \\sum_{l,m} |Y_l^m\\rangle \\bigl[V_l(r) - V_{\\mathrm{loc}}(r)\\bigr] \\langle Y_l^m|\n", "$$\n", "在平面波基底中需要对 $(\\mathbf{G},\\mathbf{G}')$ 的所有组合计算积分,复杂度约为 $O(N_{\\mathrm{PW}}^2)$。随着平面波数量 $N_{\\mathrm{PW}} \\propto E_{\\mathrm{cut}}^{3/2}$ 增加,赝势应用成本迅速飙升,并且不同 $l$ 通道之间不能复用计算结果。\n", "\n", "KB 可分离形式将每个通道表示为 rank-1 投影:\n", "$$\n", "V_{\\mathrm{NL}}^{\\mathrm{KB}} = \\sum_l |\\beta_l\\rangle D_l \\langle \\beta_l|\n", "$$\n", "这样矩阵-向量乘法只需与波函数数目线性相关,复杂度降为 $O(N_{\\mathrm{PW}} N_{\\mathrm{proj}})$,通常 $N_{\\mathrm{proj}}$ 只等于 2-4。相比直接提升平面波截断能(方案一:牺牲效率换准确度),KB 方案(方案二)在保持半局域精度的同时把数值工作集中在少量投影子上,是主流平面波 DFT 代码的共同选择。\n" ] }, { "cell_type": "markdown", "id": "fa2461bd", "metadata": {}, "source": [ "## 投影子公式推导:$\\beta_l \\propto (V_l - V_{\\mathrm{loc}})\\phi_l$\n", "\n", "1. 以选定的局域通道 $V_{\\mathrm{loc}}$ 作为背景,把剩余通道的势差写成 $\\Delta V_l(r) = V_l(r) - V_{\\mathrm{loc}}(r)$。\n", "2. 对应角动量通道的径向规约波函数记为 $\\phi_l(r) = u_l(r)/r$,其中 $u_l$ 为 TM 伪轨道或 AE 轨道,满足 $\\int |u_l|^2 \\mathrm{d}r = 1$。\n", "3. 把半局域势作用在 $|\\phi_l\\rangle$ 上:\n", " $$ |\\chi_l\\rangle = \\Delta V_l |\\phi_l\\rangle $$\n", " 在径向表象下即 $\\chi_l(r) = \\Delta V_l(r)\\,\\phi_l(r)$。\n", "4. 为了得到正交归一的投影子,定义\n", " $$ \\beta_l(r) = \\frac{\\chi_l(r)}{\\sqrt{\\langle \\chi_l|\\chi_l \\rangle}} = \\frac{\\bigl[V_l(r) - V_{\\mathrm{loc}}(r)\\bigr] \\phi_l(r)}{\\sqrt{\\int |\\Delta V_l\\,\\phi_l|^2 w\\,\\mathrm{d}r}}. $$\n", " 因此 $\\beta_l \\propto (V_l - V_{\\mathrm{loc}})\\phi_l$,比例因子由归一化积分确定。\n", "5. 耦合系数采用\n", " $$ D_l = \\frac{\\langle \\chi_l|\\chi_l \\rangle}{\\langle \\phi_l|\\Delta V_l|\\phi_l \\rangle} = \\frac{W_l}{Z_l}, $$\n", " 其中 $W_l = \\int |\\Delta V_l\\,\\phi_l|^2 w\\,\\mathrm{d}r$, $Z_l = \\int \\Delta V_l\\,|\\phi_l|^2 w\\,\\mathrm{d}r$。当 $Z_l$ 接近零时意味着局域道选择不当,需要重新选择 $l^*$ 或增加 $r_c$。\n" ] }, { "cell_type": "markdown", "id": "d6c3ffa0", "metadata": {}, "source": [ "## 局域道选择的物理依据\n", "\n", "- 高角动量轨道(如 d、f)在核附近受到更强的离心排斥,势能曲线最硬、更接近\"屏障\"形状。把这类通道当作 $V_{\\mathrm{loc}}$,可让其余通道的 $\\Delta V_l$ 在截断半径内保持平滑,投影子更局域。\n", "- 价层中最软的 s/p 通道留给非局域投影能够减少幽灵态风险;一旦把 s 通道设为局域势,其它通道需在更硬的背景下构造 $\\beta_l$,往往导致 $Z_l \\to 0$ 或耦合系数巨大。\n", "- 方案 A(推荐):选择未占据或弱占据的高 $l$(Al 取 d,过渡金属可取 f),仅需使用较小的 $r_c$ 即可获得平滑投影。\n", "- 方案 B:当体系确实需要把 d 通道纳入可移植性测试时,可改用 p/s 通道做局域势,但必须增大对应通道的 $r_c$ 并重新验证 log-derivative 曲线。\n", "\n", "经验法则:局域势应该是所有通道里\"最排斥\"的一条,以便其它通道在其基础上构造出的 $\\Delta V_l$ 带来有限而稳定的耦合。\n" ] }, { "cell_type": "markdown", "id": "f1d5de1a", "metadata": {}, "source": [ "## `kb_transform` 使用要点\n", "\n", "`kb_transform` 负责把 `{V_l(r)}` 与波函数 `{u_l(r)}` 转换成 `KBResult`:\n", "\n", "```python\n", "from atomppgen.kb import kb_transform\n", "kb = kb_transform(\n", " invert_results=invert_dict, # {l: InvertResult}\n", " u_by_l={l: u_l for l in invert_dict},\n", " r=radial_grid,\n", " w=weights,\n", " loc_channel=2,\n", ")\n", "```\n", "\n", "- `invert_results`:来自 `invert_semilocal_potential` 的半局域势,键必须覆盖所有角动量。\n", "- `u_by_l`:与势相匹配的波函数,推荐直接使用 TM 伪轨道(`tm_result.u_ps`)。\n", "- `loc_channel`:局域通道索引,若不传则默认取 d 道。\n", "- 返回的 `KBResult` 携带 `V_loc`, `beta_l`, `D_l` 以及 `diagnostics`,方便在 Notebook 中打印或绘图。\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "446b6b79", "metadata": { "execution": { "iopub.execute_input": "2025-12-15T08:10:40.115555Z", "iopub.status.busy": "2025-12-15T08:10:40.115436Z", "iopub.status.idle": "2025-12-15T08:10:40.882780Z", "shell.execute_reply": "2025-12-15T08:10:40.882311Z" } }, "outputs": [], "source": [ "# 导入依赖与绘图配置\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import platform\n", "\n", "from atomppgen import (\n", " solve_ae_atom,\n", " tm_pseudize,\n", " invert_semilocal_potential,\n", " kb_transform,\n", ")\n", "\n", "# 中文字体配置(兼容多平台)\n", "if platform.system() == 'Darwin': # macOS\n", " plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'Heiti TC', 'STHeiti']\n", "elif platform.system() == 'Windows':\n", " plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei']\n", "else: # Linux / Colab\n", " plt.rcParams['font.sans-serif'] = ['DejaVu Sans', 'WenQuanYi Micro Hei']\n", "plt.rcParams['axes.unicode_minus'] = False\n", "plt.rcParams['figure.figsize'] = (7.0, 4.5)\n" ] }, { "cell_type": "markdown", "id": "9ac2300a", "metadata": {}, "source": [ "## 实战:Al (Z=13) TM→KB 流程\n", "\n", "流程回顾:\n", "1. `solve_ae_atom` 得到 LDA 参考解;\n", "2. 对每个角动量调用 `tm_pseudize` 构造伪轨道;\n", "3. 用 `invert_semilocal_potential` 反演半局域势;\n", "4. 交给 `kb_transform` 得到 KB 可分离形式(局域通道取 d)。\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "b7b7e6f9", "metadata": { "execution": { "iopub.execute_input": "2025-12-15T08:10:40.884108Z", "iopub.status.busy": "2025-12-15T08:10:40.883987Z", "iopub.status.idle": "2025-12-15T08:11:02.957512Z", "shell.execute_reply": "2025-12-15T08:11:02.956293Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Local channel l* = 2, covered channels: [0, 1, 2]\n", "l = 0 (s): W_l = 1.590e+00, D_l = 4.990 Ha\n", "l = 1 (p): W_l = 7.530e-05, D_l = 2.219 Ha\n" ] } ], "source": [ "# Al LDA 示例:求解 → TM 伪化 → 势反演 → KB 转换\n", "Z = 13\n", "rc_map = {0: 2.1, 1: 2.2, 2: 2.4} # s/p/d 截断半径(Bohr)\n", "loc_channel = 2 # d 通道作为局域势\n", "continuity_orders = 4\n", "\n", "# 1) 全电子参考解(LDA,自旋无关)\n", "ae = solve_ae_atom(\n", " Z=Z,\n", " lmax=max(rc_map.keys()),\n", " grid_type='exp_transformed',\n", " spin_mode='LDA',\n", ")\n", "\n", "# 2&3) TM 伪化 + 半局域势反演\n", "invert_results = {}\n", "tm_results = {}\n", "for l, rc in rc_map.items():\n", " # 选取该通道能量最高的参考态(列表末尾;p/d 通道通常为近零散射态)\n", " u_val = ae.u_by_l[l][-1]\n", " eps_val = ae.eps_by_l[l][-1]\n", " tm_res = tm_pseudize(\n", " r=ae.r,\n", " w=ae.w,\n", " u_ae=u_val,\n", " eps=eps_val,\n", " l=l,\n", " rc=rc,\n", " continuity_orders=continuity_orders,\n", " )\n", " tm_results[l] = tm_res\n", " invert_results[l] = invert_semilocal_potential(tm_res, ae.r)\n", "\n", "# 4) KB 可分离形式\n", "ae_like_u = {l: tm_results[l].u_ps for l in invert_results}\n", "kb = kb_transform(\n", " invert_results=invert_results,\n", " u_by_l=ae_like_u,\n", " r=ae.r,\n", " w=ae.w,\n", " loc_channel=loc_channel,\n", ")\n", "\n", "channel_labels = {0: 's', 1: 'p', 2: 'd'}\n", "print(f\"Local channel l* = {kb.loc_channel}, covered channels: {sorted(invert_results.keys())}\")\n", "for l in sorted(kb.beta_l.keys()):\n", " norm = kb.diagnostics['projector_norms'][l]\n", " print(f\"l = {l} ({channel_labels[l]}): W_l = {norm:.3e}, D_l = {kb.D_l[l]:.3f} Ha\")\n" ] }, { "cell_type": "markdown", "id": "324b593e", "metadata": {}, "source": [ "## 可视化:局域势曲线 $V_{\\mathrm{loc}}(r)$\n", "\n", "局域势的形状决定了其它通道的 $\\Delta V_l$。以下示例展示 d 道(局域势)的径向曲线,并在 $r_c^{(d)}$ 处做标记,检查是否足够平滑且单调。\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "feeb1a4b", "metadata": { "execution": { "iopub.execute_input": "2025-12-15T08:11:02.965463Z", "iopub.status.busy": "2025-12-15T08:11:02.965314Z", "iopub.status.idle": "2025-12-15T08:11:03.039150Z", "shell.execute_reply": "2025-12-15T08:11:03.038692Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAAGuCAYAAABFrtg/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAT/FJREFUeJzt3Ql8E2X6wPGn98VNW5D7vgQXOVVARORakeW/3qigKLjijawreKIrKKjrLaioCN7iBSuKCAioCAqeICKCgsi2XG2hdzv/z/PCxDRNS9JOaNL8vp/PkMlkMpm8mZCn7/G8EZZlWQIAAICAiQzcoQEAAEDABQAAcAxQwwUAABBgBFwAAAABRsAFAAAQYARcAAAAAUbABQAAEGAEXAAAAAFGwAUAOKYCmW+bXN4IVgRcwBFjxoyRiIgIefDBB0uVyWmnnWaWsqxYscI894UXXij12MqVK6VGjRpywgknSFpammtfvfWVHlefs3379pD/vOz34r7ExcVJhw4d5L777pPi4mK/jqfPv+uuu+RYqMhnFwzHPppLL71UWrRocUxe64cffpCTTz75mB27ItfHsbymED6iq/oEgGCQlZUlb775pnTr1k2efvppmThxovlPt7JWrVolf/3rX6Vdu3by0UcfSf369WXjxo1+H+fMM8+Uzz//XI477jipLp544glT3ionJ0fWrFkjd955p6mhmDx5clWfHgLkjTfekC+++OKYHVu/N02aNAnI6wH+IOACROS1114z5XD//ffLoEGDZNmyZTJw4MBKlc3q1atNsNWpUydZsmSJ1KlTp8LHSklJMUt1ouVy0kknue4PGDBAtmzZIk899RQBFxzjfo0BVYkmRUBEnnvuOTn99NNNkNWyZUuZPXt2pcrl008/lWHDhplmxKVLlx412NqwYYMMGTLEND3WrVtXLrzwQvn999/LbFLUJqChQ4fKiy++aGrPtEnuL3/5iyxatKjEc+Lj403gpzVJuq77vvfee/Ljjz+aACcxMVFatWpljuPuv//9r/Tv39/UyCUkJEjHjh3lySefLNX8pTV4gwcPNsdp2LCh3HTTTVJQUFDhctNy0tpGtW/fPnPOU6ZMKbGP1obVrl1bpk6d6vUYWoN47rnnSqNGjUy5aO3G9ddfL7m5ua599NyfffZZGT9+vNSrV09q1qwpZ599tvzxxx8ljvXYY4+Z8tHzOOWUU0yTlb+fnTZNtWnTRu69917zOm3btpUDBw74VB7fffednHXWWVKrVi1JSkoy16h+nu70WNdcc415v7pPr169SlwHeXl5cscdd0j79u3N56TneMYZZ8g333wjvrI/b62l1SY7LQ8tl0cffbTEfkVFRaZZWN9jbGysNG/eXG655RZX2WtZ2J+be7OdPk+363dPP7POnTvL888/X+LY2uT573//W/71r3+Za03fi157ei2Xd2zP5kFfrg8gICwgzG3atEl78FpvvPGGuT916lQrJibG2r17t2uf/v37m6Usy5cvN8d4/vnnrdWrV1s1atSwTj75ZCsrK6vMffVWbd682ezftWtX680337Tefvttq3fv3lbLli2tAwcOmH30uPqcbdu2mftjxoyx6tata3Xu3Nl67bXXrP/+97/WCSecYMXFxVk7d+50PScyMtJq2rSp9eSTT1offvih1b17d6tmzZpWmzZtrBkzZlhLly61hgwZYt7vr7/+ap63ZMkS87wbbrjB+vjjj60PPvjAGjlypHn9FStWlHgPjRo1sqZNm2b2u+aaa8y2mTNnllve9nvR1y4oKDDLoUOHrJUrV1qpqanWxIkTXfued955VrNmzazi4mLXtpdeesmKiIhwlYUe68477zTr+pklJydbQ4cOtRYtWmTOd8qUKWafu+66y3UMvV+vXj1r7Nix5jwee+wxKzY21jrzzDNd+zz77LNmv+uvv96U3W233Wb28fez03PT8tV93n//fXP+5V1D9rG//fZbKykpyerWrZv1+uuvW++88441YMAAKzo62pyzKiwstE455RRzLTzxxBNm+yWXXGJFRUW5PquLL77YatCggTVnzhxTxvPmzbOaNGliztGm11Pz5s3L/Mzsc9PXmTRpkrkmJkyYYLY988wzrv30tfX8br31VnMd3X///VZiYqI1cOBA8xnu2LHDuvzyy83zPv/8c3NfjR8/vsTz9LPSa/mRRx5xHVvPTz+zESNGmNefO3euOZ/jjz/ePF7WsSt6fdjPAZxCwIWwd/PNN5v/yPPy8kxZaOChAce9997rd8ClPxwa0Oh627ZtfQq4LrroIvP6e/bsce2jP9YafNxzzz1lBlx6/+eff3Y9R3+EdJv+8Ls/Z9asWa599Idbt2nwYFu7dq3ZpsGC0oBp9OjRJc553759JYIp+z24H6eoqMgEYMOHDy/3mrLPy9vSt29f1+egNNBxLyulAeLpp5/u9cfxo48+MgHIwYMHS7ymBqPuwZT9Wu4uuOACEzzZ76Vx48bW2WefXWIfO6j057PTc/N8D954Xhf62vXr17cyMjJc+2hw2q5dOxOEqYULF5rnaDBm03M/6aSTrLvvvtvKz8+3Bg0aZIJydw899JB5Xnp6ul8B15VXXlli+9/+9jdTTvqaGiDqPu7fG/Xiiy+a7e+9916J8rD99NNPJoC+4447SjxPA3n9LmkwrvT8WrRoYYJM2y233GKOZZe/57Erc30QcMFpNCkirGlTxrx580xzUn5+vhw8eNA0MWlz2jPPPOP3iDntcK9NLu+++678/PPPcvXVVx/1OR9//LFpHtRmssLCQrPYTUPal6ws2qerdevWrvv2+qFDh0rsp01h7s/x7NeSnJxsbu1mrkmTJsncuXNNU9S3334rCxYskBkzZng9tvuIsMjISNMkZO+jZWe/H120rN3NmjVL1q1bZxZtJtNm3F9//dU0xdrNO9r01axZM/MZKW3y0ybayy67zGuZ6P7anKvNoNof7P333zf98v73v/+Ve+52+dn77Ny50zQLjhw5ssQ+5513XoU/O3uAgK8++eQTGT58uGlOtEVHR5smy/Xr15umV32vMTExptnR/XPQjuK33367eUz7D+p5p6enm3KeM2eOaTJWnmVyNNqU7U6/N1pOWl46GleNGjWqxD56vlFRUeb9eKPlpDGOPs/9etHPUt/jl19+6dq3Z8+e5lhHu+bL4s/1ATiNgAthTf/D1R9xDa60f429LF++3PSX+vDDD/06ngYL2kdqxIgRpl+N9o166aWXyn3Onj175OWXXzY/ju6L9sPZtWtXmc/TPizu9MdYeQaJ+n48aVBg8xyNqeXx97//3fzQa5Bw2223yd69e73mOPJ2Dvbrjx07tsT7cQ8OlfYp6tGjh1n69Olj+lPNnz/f/ADb/Xc0eNAfeR1Bqn239HE9dz0/b7T/mAaM2vdM+6vpMfUHW/sc+XLu9j52H6wGDRqU2Ef7/VT0s3MPnHyxf/9+01fJk70tMzPTvL6+Vy2nsixcuFC6dOkiqamp5vrUvn32teJvzqqmTZuWuK/HtMtLz9f9/Gz6WhrUZ2RkeD2mvgelaUHcy1CDVuVejr5e82Xx5/oAnMYoRYQ1/WHXzrh2DYpN//PVoElrYfRHyldak6AdcZX+5ay1MVdddZWpUfIMONw7iutf3trh3JN9rGNJA5yffvrJBKMaCOmPkdaOaFDqD+2orEGnP+9Fgy+7s7hNa7PuueceE8i++uqrcv7555f64bU98MADpqO71pbp56e1ler444/369y1Y7nSvGnegoNj8dnpOezevbvUdg1ANMDSx3XRYFgDDvegS2vA7GNoB/FzzjlH3nrrLdN5XwPsmTNn+v3HhP3+Gzdu7Lpvn58GVHaZ6Tb3nF4a5Ojz7JpUT/aAEg20vf1xoLWmTnHq+gAqgoALYUuDCK2J0L94+/bt6zV40uYX9xFn/tBmCw3ktOnqggsukM8++8z85e5Jmy915FTXrl1L1DxooKE/BFo7cSzpeepf/u5pMd555x1z608Tq/7o+ptM0w4UdASc+3F0dJ6mi9DHH3/88XLPXUe4uTd9bdq0yQSQWqvhK61909Frr7/+ulx00UWu7dpUfKw+u1NPPdVcn1qTZdeOabOspjDp3r27CTr19TV40uBYmx/t19faRR0pqNedNg3ryD69b3+GWutlr/tD37+OhnXPe6VNvhrIaQ2k0ho/95GlGiTreffr18/ctwNDPU8N/vQ9KA3KdOSse3Ptf/7zHzM61tfPzvPYgbo+gIog4ELY0mBI//rWPiZlZZ7XPln+1uy40x9GHZKv/Wn0R0h/HD1pss/evXub5KZaG6Y/pJqyQGskPvjggwq/dmXOWZvutDlRm9Q0JYDW9GkAmZ2d7djraKCitWdKf6x1eL+mFNDkrp59hTSA0MBHA6HyspTruWvwoX3OtElK0zjoMbWPlT/nrj/WWkOpr/mPf/zD1BJp3x8th2P12emxtf+dznCgiWA1zYIGH7/88ossXrzY7KO1r1oeWl7Tpk0zgaoGPBpEaICqtUoaCGpqhhtuuMH0UdRrWh9X/n6e06dPN380aDlrM6+WtQZUWl6aAuXiiy82563H1UBKayr1+tdAyq4ptmu0XnnlFVPzqwGQPm/cuHGmGV9rOTUA0qZsDVo1oPOV57HdA3cnrw+gQhzvhg+ECE2pYA8pL4uONNRRWDqizde0EJ50VJWmiNCRWDqS0HM0mvryyy+twYMHm1FyOjJLX09H6Nm8jVL0HFWmj7mfg+dz3M/T/bU9n7d9+3Yz0rB27dpm0dFyGzduNCMD7ZF93o7jy2hO9/NyXzRtgpazjvrbunVrqefoSD0tv/vuu6/UY+4jynJzc62rr77apEGIj4+3+vTpY1Jm6Ig9TU+QmZlZ6jk2byPcNIVD+/btTTqIHj16mBGg/n523o7rjbcy1ZF/mgZBj5uQkGCdccYZ1qefflpqBOm4ceOslJQU8x71PS9btsz1uKY70etc0yzoKD8dDajpLPS1NPWFP6MUdZTqiSeeaMpDj2mnUnG/1jXdiH5v9DPVlB46kjUnJ8e1zx9//GH16tXLHEM/K6WjKbWcNFWFbtdUJpqWxP68lJ6fnqc7z2vc27Gduj6AyorQfyoWqgHAsaHNaFoL8ttvv1Wr6Y1ChdZyai2VDiYpb05RAGWjSRFA0NI+TJrJXeddHD16NMEWgJBFWggAQWvz5s2mv43mX3rwwQer+nQAoMJoUgQAAAgwargAAAACjIALAAAgwAi4AAAAAqxajVLUrMk67YVOD+EtyzAAAIBTNLOWTrKu86yWN6dptQu4NNjynFwVAAAgkHbs2GGmAwubgMue+FTfuD33mNM1aDr/XkpKSpmRbHF2tmzpd6pZb7tqpUSWMcluuPOlLEE5cj2GFr7XlGO4XZOZmZmmosfbxOvVOuCymxE12ApUwJWbm2uOXWbAFR0tNaKiXOdBwFXxsoQz1yQox2OF65FyDNdrMsKHbkz8Dw0AABBg1aqGKyhERUlS/1Nd6wAAAARcDouMi5Nms2dzZQEAABeaFAEAAAKMgAsAACDACLgcpmkhfjyxm1l0HQAAgD5cAWDl5HBlAQAAF2q4AAAAwi3gysnJkbFjx5qsrQ0aNJBp06ZV9SkBAABUrybFf/7zn/Ltt9/K2rVr5ffff5dzzz1XWrZsKRdeeGFVnxoAAEDo13Bp7dacOXPkgQcekI4dO8oZZ5whN954ozz11FNVfWoAAADVo4br66+/loKCAjnllFNc2/r06SPTp08Xy7J8mqsIoUs/42Lr8K1l7osU6z9H1g9v9eU4frymn+fn9HH9OVdvBy22iiUzt1DicgokMiKo/n4KKZQj5RhMuB6dLcucgiIJBkEVcO3cuVOSk5MlNjbWta1Ro0Zm4sm9e/eax9zl5eWZxX3WbnuySl2cpsc0QUE5x9ZHEnr2dK1LAM7DF3qeepFl5hRKVm6BZOXprb0USHZ+keQWFkt+weHbvIJiyS0sktyCIskrLJbcgmIpKCqWomJLCostKSouPnJrSVGRvc3LY8WHgyYTHpkg6c/A6fDmksGUvQ4AQCAMbFtHZo9pEJBj+xNrBFXApYFVXFxciW32fW1u9KQ1X1OnTi21PT093RwrEAWbkZFhgpnyZh2PnznD3O7RAPBIEOjoeViWpGUVyO8ZebIzI092Z+bJ3uxC2XuoQPZnF8je7ALZl10oBUVEMgCA8FZQUCBpaWnl/m5XVFZWVmgGXPHx8aUCJbsGKyEhodT+kydPlokTJ5ao4WratKmkpKRIrVq1AhJwabOmHj8QH5w32fmF8s2ODPnhj0zZuOvwsn3vIcn3MZiKjBCpFR8jNeOjjyyH1xNjoyQuOkriYyL/vI3RbZESr0tMlERHRUhMZKRERUaY9Wh7PTLCte3P+5HmVl9Py0hbf7UB2Kyb25LrGrTuO1Jrqc8tsX+EnvfhfcXLcXxx5Nm+7etHS7VfjdrH4Fz1mtyTvkeSU5KP2TVZHVGOlGMw4Xp0uCz37JHU1NSA/B+pcUtIBlyNGzc2TYcajcbExJhtOlJR31D9+vVL7a+1X541YkoLNVA/PvrDH8jjayCy8Y9MWfLD/+SzrXvk6x0HvNZUaXDTpG6CNK+fJE3rJUhqzXhJqRknKTXiJLlmnCTXiJW6ibEmsArGvm/6JYjOPyiptRMIFCpZjjHRGixHU46UY5XjeqQcg/KajIoM2O+2P8cMqoCra9euEh0dLZ999pn079/fbFu1apX06tUrKIMGb3Q6n58HnmHW23y8VCITE3163oHsfHl57W+y4KudsjX9UInHGtWOly5NasvxjXSpJe0a1JTjasdLdBQ1GgAAhIKgCrgSExNlzJgxctNNN8ncuXNNm+ujjz4qs2fPllBStH+/z/tmZBfIo8u2yMtf/OYaSREbHSkD2qfIae1TpU/rZFODFSoBJwAACPKASz344IMyYcIE6d27t9SoUUOmTJki5513nlRHb2/YKXcv3Gg6uquOx9WSy/q0kGGdG5q+VgAAoHoIuoBLg6wXX3zRLNWVpl649e3vZcH6neZ+29QacuuZHaV/uxRqsgAAqIaCLuCq7g7lFcr4eV/Kpz/vNSP6bjijnVx1WmvTqQ8AAFRPBFzHUGFRsfxj/lcm2EqKjZKnR/eQPm1KJnMFAADVDwHXMfTv/26SVVv2SEJMlMy/orec2KzusXx5AABQRQi4nBYZKfGdO7vWbau2pMsLn2036/85vyvBFgAAYYSAy2GR8fHS8s03SmzLyS+SyW99Z9bHnNxchnZu6PTLAgCAIEZP7WPgpS9+lZ37c0wC05uHdjgWLwkAAIIIAdcxmAvxqRVbzfr1Z7SVpDgqFQEACDcEXA4rzsmRn08faBZdf3vD77L3UL40q5cof+/WxOmXAwAAIYDqFqdZlhTs2uVaf2Xtb2Z19MnNybUFAECYooYrgDbuypTvf8+U2KhIOZvaLQAAwhYBVwAt2bjb3A7smCp1k2ID+VIAACCIEXAF0Meb0swtaSAAAAhvBFwB9HP6QYmOjJDT2qcG8mUAAECQI+AKsO7N60rthJhAvwwAAAhijFJ0WkSExLZpLbszcvWOnNSqvuMvAQAAQgsBl8MiExKk1cKFctH0ZZKXmSu9W9Vz+iUAAECIoUkxAHbsy5HdmbkmHUS3ZnUD8RIAACCEEHAFwPe7Msxth+NqSnxMVCBeAgAAhBACLofpdD7J146RWR/PlC7145w+PAAACEH04XKaZUnN3TukpoikN9R/AQBAuKOGK4C0SREAAICAy2H7s/Nd6+0a1OAKAwAABFxO27w7y7VeM46EpwAAgIDLcVvTD3FdAQCAEmhSdNi2PQedPiQAAAhxjFJ02La92fK/hLpSNynWTPMDAABADZfDfjpQKJcOuVUK5r9lpvkBAAAg4HJQTn6R/H4gx6y3SmGEIgAAOIyAy0Hb9x7uMF8nMUbqaZMiAAAAAZezdh3IkdiiApmx9CHZds65Upyby0UGAADoNO+kXRm5EmEVS9O0XyU3TSdWLOYSAwAANCk6aXfG4f5bAAAA7ujD5aA/DtCECAAASiPgctAfGQRcAACgNAIuB+3OJOACAAClEXA5aM/BPCcPBwAAqgkCLj98/3uGrPh5v/ySXnq+xMKiYsnKLTTrEXXqSlTdus59SgAAIKQxl6IfXvriN3nty51y06AoubZBrRKPHcgpMLd50XHS5tPVEh1FLAsAAA4jKnDIgex8c1srPppgCwAAlEDA5ZD92YdruOoypQ8AAPBAwOWQfYcO13Alx4r8esloszC1DwAAUPThcrhJsV58tGSvW3d4I1P7AAAAaricb1KskxjDhQUAAEqgSdEhB44EXLUJuAAAQCgEXFu2bJH4+HgJJQfzDgdcNeOo4QIAAEEecO3cuVNGjBgheXmhlbX9UF6RuU2Ki6rqUwEAAEEmqAKud955R7p37x5ytVvqYN7hLPNJcYxDAAAAJQVVdLB48WKZOnWqdOjQQQYMGCCh5JAdcMVGSURCQlWfDgAACCJBFXDNnj3b3K5YsUJCjR1wJdSqKR02rK/q0wEAAEEkqAIuf2k/L/e+XpmZmea2uLjYLE6zLMt163l8u0kxMSYyIK9d3WgZeStHUI5VgeuRcgwmXI+hU5b+HLfKAq5p06aZxb05sV+/fn4dY/r06aYJ0lN6errk5uaK03KOHPPgoUOSlpZW4rHMnMOJT/OzMyUt7XAHepR/kWZkZJgvQmRkUHUlDCmUI+UYTLgeKcdwuyazsrKCP+CaMGGCjBo1ynW/YcOGfh9j8uTJMnHixBI1XE2bNpWUlBSpVauWOC0hfre5rZGUJKmpqSUeyy08XPvVpG4tyb/jNrPe6JGHJTIuzvHzqC5fgoiICPNZEXBRjlWN65FyDCZcj6FTlv4M8quygKtOnTpmqYy4uDizeNJCDUTB6odm37ofXyNnuw9Xzbho2b9y5eHzoPbmqOUZqM8qnFCOlGMw4XqkHMPpmoz045j80jkgt6BYig9XcEliLHm4AABASQRcDrA7zKvEGAIuAAAQAgHXaaed5hoRGGo5uCIjDzc7AgAABHXAFWrIMg8AAMpDwOWA7Hx7HsWQTmsGAAAChIDLAbkFhwOuuGiKEwAAlEaVjIMBV3xMlEQmJkrHHzc5cVgAAFBNUCXjgNzCw6n942MoTgAAUBoRgsM1XAAAAJ4IuByQZwdc0VFSnJcnO6+/wSy6DgAAQMDlUKZ5V5NiUZFkffihWXQdAACAgMsBNCkCAIDyEHA5ILeQPlwAAKBsBFwONinGMUoRAAB4QcDlZJNiNKMUAQBAaQRcjnaaJ+ACAAAEXAHuw0X8CgAASmNqHyfzcMVESURCgrRf/5W5r+sAAAAEXA7n4YqIiJCIxESuLAAA4EIbmAPoNA8AAMpDwOVwHq7i/HzZdctks+g6AAAAAZfTebgKCyXjnXfMousAAAAEXA5gah8AAFAeAi4nO82T+BQAAHhBwOVoWgiKEwAAlEaE4IC8QjLNAwCAshFwVVJxsSX5RUc6zUdTnAAAoDQihEqygy0VS8AFAAC8INN8JRW4BVwxUZESEZcgbT/71Nxnah8AAEDA5YD8I/23VKwGXBEREl2vHlcXAABwoUmxkgqKLHMbExUhkZERlT0cAACohmhSdKiGS5sTlU7nk3bffWY99ZZbJDI2trIvAQAAQhw1XJWUX1RUssN8YaHsf/kVszC1DwAAIOByQH6h3aRI7AoAALwjSnAoLYR2mAcAAPCGKMGhtBDk4AIAAGUh4HKo0zw1XAAAoCwEXA41KcZEkxICAAB4R8BVSdRwAQCAoyEPl0N9uOxRihHx8dJ66VLXOgAAAAGXUzVcR/JwRURGSmyTxlxZAADAhSZFh2q44uzEpwAAAB6o4XJ4ah9Lp/Z5+BGznnrD9RLB1D4AAIQ9qmUqKc+jSdEqLJR9zz1nFl0HAAAg4KqkgiKm9gEAAOUj4HK40zwAAIAnogSnpvZhLkUAAFAGAi6nJq+mhgsAAIRCwLV8+XLp2bOn1KhRQ0444QRZtGiRhM4oRab2AQAAQR5w/fzzzzJ8+HC56KKL5Pvvv5errrpK/v73v8t3330nIVHDFRVV1acCAACCVNDk4XrllVekR48ecsMNN5j7GnC9/fbb8tprr0mXLl0kWBUUlpy8WqfzabXwPdc6AABA0ARc559/vgwbNqzEtoiICDl06JCERg3Xn1P7xLVtW8VnBQAAgknQBFzt2rUrcX/Tpk2mT9c111wjwYypfQAAQMgEXO7S0tJk5MiR0r17d9Ovqyx5eXlmsWVmZprb4uJiszjNsizXrX38vILDt9GREWabTu2z9+lnzLb648cxtU8ZTFm5lSMqhnJ0BuVIOQYTrsfQKUt/jltlAde0adPMYlu8eLH069dP9u3bJ4MGDZLc3Fx58803TbNiWaZPny5Tp04ttT09Pd0832k5R4558NAhExSa9ezD23KzD5ptVk6OHHjySbOt6KzhEpGQ4Ph5VAd6kWZkZJgvQmRk0IzdCDmUI+UYTLgeKcdwuyazsrKCP+CaMGGCjBo1ynW/YcOGJmAZOHCgeQPLli2Txo0bl3uMyZMny8SJE0vUcDVt2lRSUlKkVq1ajp9zQvxuc1sjKUlSU1PNekTUNnNbv24ds604O1sOHNlfzyMyMdHx86guXwINpk0ZEXBRjlyP1QLfa8ox3K7JeD8Gx1VZwFWnTh2z2LRG6q9//avk5OTI6tWrpUmTJkc9RlxcnFk8aaEGomDt2ja9tY9vd5qPi4k6vM3tdQN1HtWFXY6UEeUYDLgeKcdgwvUYGmXpzzGDpg/XI488YvJvrVixwkSMe/bsMdt1XROhBvvk1UztAwAAyhI01S/aX0s7wJ988smm6s9egn2UIpNXAwCAkKnhWrdunYQiOy1EDJNXAwAApwKu7du3y0cffWQCpN27D3ci187i3bp1M4lLW7ZsKeEk70imeSavBgAAlW5S1L5VQ4YMkQ4dOsgzzzwjRUVF0rFjR+ncubNER0fLq6++Kn/5y1/MPitXrpRw8WcN15EO9XFx0uKN182i6wAAAD7VcGn6ht9++03GjRsnb731liQlJXndT0cavvfee3LbbbdJs2bNZP78+dW+hAuLrRJNihFRUZIQxHM/AgCAIA24Lr30Uhk8ePBR99MRheedd55ZlixZIuFUw6WZ5gEAACrcpOhLsFVYWCiff/65X8+pDgqLStZwmal95swxi64DAAD4nRZCk5Jqvy2tzYqKinItmoC0f//+YVeihUfmUYo6UsNlFRZK2swHzKLrAAAAfgdc119/vbRt21YWLlwoiYmJ8vrrr8v9998v9evXL1HDFQ50biY78Wn0kU7zAAAAlU4LsWnTJlmwYIG0aNFCevXqZYKuSZMmmXkP77rrLhOIhYuiIx3mVQxT+AAAAKdquGrWrGkmg1SdOnWSDRs2mPV+/frJ8uXLJZzYIxQVNVwAAMCxgOuMM86QW2+9VdLS0kyQ9cYbb8iBAwdk0aJFUrduXQnXgItM8wAAwLGA6+GHHzYBlgZaZ599tmlSrFevnlx77bVy5513SjgpPJISQpEWAgAAONaHSyeUXrx4sev+qlWr5LvvvjPbGzVqJOHE7jDvPkoRAACgQgGXZpkvjzYlah4u3U8zzIdbSgid1ici4s+pfZrNnetaBwAA8Cng0hGJdkBhp0Ow79vr9q3OsRgu7KSn0W4jFHVqn6TevarwrAAAQEgGXNu2bStxX4Mrnaj6/ffflyZNmki4ck3rQw4uAABQ2YCrefPmpbZFRkZK06ZNw6oJ8WgTVyuroED2v/66Wa973nkSERNTZecHAABCtNM8yp+4WgOu/93zb7Ne5//+j4ALAAD4nxYCZU9cDQAA4E2lIgX3jvThyB6lSB8uAABQ6SbFsWPHltqWnZ0tN910k9SoUaPE9ueee07ChWvianJwAQCAygZcOirR06hRo8p8LFzQpAgAABwLuJ5//nmfDhZuCmhSBAAATvXh0jkSMzMzxVcZGRlhMa+iXcMV5Zb4FAAAwJNPkYLOkaiJTnWC6tWrV0tBQYHXPl0ff/yxXHHFFdK1a9ewSIhqT14d49aHKyI2VprMesosug4AAOBTk+KVV14p5513njz++ONywQUXSHp6urRq1Urq1asnxcXFsnfvXpONXoOscePGyYwZM8xj1V3BkcSn7qMUI6KjpeZpp1XhWQEAgJBNfKoTVN9+++1m+eabb8ySlpZmHmvYsKH07NlT2rdvL+GkyDV5NU2KAADA4Uzz2ryoS7jzlhZCM81nLFxk1mufNZxM8wAAgKl9nOg0H+0xl+IfU6aY9VpDhxBwAQAApvZxItN8jFsfLgAAAE90PnKkSZFiBAAAZSNScCAtBHMpAgAARzvNf/LJJ7JkyRJZt26d7N6922xLTU2Vbt26ybBhw2TAgAESLgqPpIWIoYYLAAA4EXDp9D6aXysrK0tOOeUU6devnyQnJ0tUVJTs27dPNm7caHJwxcTEyJQpU+SSSy6R6q6AGi4AAOBUwNW/f3+T1HTWrFlmvTxr166Vxx57TObMmSMrVqyQ6ozJqwEAgGMBlwZaHTt2NOuaUb5ly5Zl7turVy+ZN2+e/PjjjxI2k1d7TO3T+OH/uNYBAAB86jRvB1t2QKVZ5o+mQ4cOYZmHS6f2qTV0qFl0HQAAwO9RitpBfufOnZSc++TV5OECAADl8LsKpnv37vL3v/9dTjrpJGnWrJnpJO/uueeek3DhmrzabZSiVVgoWUuXmvWaZ5xBLRcAAPA/4NJRiaNGjXLdt6zDQUc48paHy8rPl99vuNGst1//FQEXAADwP+DS9BDwHKXI1D4AAKCSfbgGDhxoEp76atmyZTJo0CAJlybFKBKfAgCAytZwPfLII3LNNdfI3r175fzzz5e+ffuakYv16tWT4uJis/27776T1atXy6uvviqNGzeWJ554Qqo7Os0DAADHAq7OnTubJKYrV640Obk0ANMgKyIiwtWPq2HDhjJ48GCZO3euyUQfTlP7MHk1AABwrA/XqaeeahYNsHbs2CFpaWlmuwZbmok+3DB5NQAACEgeLrVo0SLZvHmz9OjRwywzZ86UDz74QMKNa/JqOs0DAAAnAy6dJ1HTQuzatcu1LT4+Xi666CIzf2JlfPjhh9K1a1epUaOGaZb88ssvJSQmr3brNB8REyPHTZtmFl0HAADwO+B66KGHZMGCBTJmzBjXtvvvv990lp82bVqFS3Tr1q0ycuRIGTt2rGzcuFHOOeccGTZsmOzbty+k0kJokFXn7/9nFgIuAABQoYArPT1dWrVqVWq7Tmj9xx9/VLhUtU/Y+PHj5brrrjMZ7CdOnCiFhYWyatWqkMo0DwAA4MnvSKFfv35y9913S05Ojmubrv/73/+Wk08+WSrqtNNOM6Mflaaa0CmCDhw4ILVr15aQyjSvU/usWGEWXQcAAPA70/zjjz8uQ4YMMSMT27dvb0Ys/vTTTyb31rvvvlvpEtWmRU1DkZubK5dddpn079+/zH3z8vLMYsvMzHQFbLo4zZ7GSG/1+HYfLo237Ncrzs2Vnf+4yqy3/XKdRCYmOn4e1YGWl12OoByrGtcj5RhMuB5Dpyz9Oa7fAVfr1q1NHyvt4K6BVnR0tHTo0MFklo90oGktJSVF1q5da5YXX3zRNFM2atTI677Tp0+XqVOnem321IDNaTlHjnnw0CGTEiMvv+Dw/cwMSUs7XOiWW82fnkdEQoLj51Ed6EWakZFhvghOXDfhinKkHIMJ1yPlGG7XZFZWVuACLhUbG2sWDba0n1VRUZF5U/68Ge1g797JfvHixaa5slatWtKlSxezfPPNNzJ58mSTTNUbfUz7ernXcDVt2tQEbXocpyXE7za3NZKSJDU1VayIw+83JbmepKbWNevF2dly4Mj+eh7UcHmn14smzjVlRMBVYZSjMyhHyjGYcD2GTllqloaABVw7d+6Uv/3tbyYPlzYp6pvRde3ortnotanRFxMmTDDpJWyHDh0yUwPptEE2rTkrb4qguLg4s3jSQg1EwdqZ9fVWj2+PUoyNjvrz9dxeN1DnUV3Y5UgZUY7BgOuRcgwmXI+hUZb+HNPvV9c5FRs0aGBGFX711VeyYcMGE4Rp8DVp0iSfj1OnTh1p0aKFa9FkqldeeWWJfTQPlwZywYpRigAAICAB19KlS00urrp1DzehKZ3E+r777jP9uirqkksukd9//13uvPNO2bZtmzz11FMyb948ufHGGyVYMXk1AAAISMCVnJws+/fvL7VdJ7NOqEQHce0YrwGbLtp/a/bs2fLGG2+YCbGDld2kGB1FsyEAAHCwD9eFF15ossFr36pevXqZnv9r1qwxTY1aS1UZvXv3NscKFQVHhoNGR5bMNN/g9ttc6wAAAH4HXJr0VOdR1JonOy+VdkjTQEybA8PJnzVcJQOuehddVIVnBQAAQj7giomJMWkaHn74YdmyZYvZ1rZt2xJ9usKBBpuFTO0DAACcCrg0AWl5fvzxR9f66NGjJRzYwZbn5NVWUZFkf/mVWU/s0V0ioqKq5PwAAECIBVy+NhVq02K4BFxFbgGXe6d5Ky9Pfhszxqy3X/+VRDC1DwAAYc+ngEvTNPhj5cqVctJJJ5ls9NWVPY+iZ6d5AAAATwHJZzB8+HDZvfvwNDjVvcO8iiEtBAAAONYBlz16sTqzU0LobD9R1HABAIBykLGzkjVcMcyVCAAAjoKAy8EcXAAAAN4QcDmYZR4AAMCRxKfwaFL06DAfER0tqf+c5FoHAADwKSLIz8/3K8VDp06dqnVKCPe0EJ5NihGxsVL/8sur6KwAAEDINimmpqbKFVdcIcuWLfNpBOIXX3whDRs2lOqMaX0AAICjAdfs2bNl//79Jr9WkyZNZOLEibJ+/XoJZ4VHarjcp/Wxp/bJ+e47s+g6AACATwHX+eefLwsWLJD09HSZOXOmyTzfp08f6dChg9x9993y888/h11JFrhGKZYsQp3aZ/u555lF1wEAAPwapZiUlCSjRo2St99+W9LS0uTWW2+Vr776Sk444QTp3bu3PProo2FTooWMUgQAAIFOC1GzZk255JJL5N1335W33npLCgoK5MYbb5RwH6UIAADgqULRQlFRkSxdulQmTJggjRo1ktGjR5vJqj/55BMJF/YoRab1AQAAR+NzoihNDbFkyRJTm/Xee++Z+yNGjJBnn31WBg8eLNFhlnPKHqXo2WkeAADAk09R0gUXXCCLFy+W3NxcGTp0qDzxxBMm2EpISJBw5crDxVyKAADAiYBLO8g/8MADcs4550jdunV9eUq1V3Skhou5FAEAgCMBlyY8he9T+yRffbVrHQAAgIjA4cmrdWqflGuv4coCAAAu5DSoINJCAAAAX1HD5fDk1VZxseRv3WrWY1u3lgg61QMAEPYIuByevNrKzZVfzhph1tuv/0oiEhPD/iIDACDc0aTo8OTVAAAAngi4Kj15NQEXAAAoHwFXpSevpggBAED5iBYqPUqRGi4AAFA+Aq5KNylShAAAoHxEC5VsUozxSHwKAADgibQQDtdw6XQ+9caOda0DAAAQEVQyLYTnKEWd2qfBzf/kygIAAC40KVY68SlNigAAoHzUcDmdab64WAp2/WHWYxodx9Q+AACAgMvpTPM6tc/WM84w60ztAwAAFE2KFURaCAAA4CsCrkpnmqcPFwAAKB8BV6UzzVOEAACgfEQLFVRQRloIAAAATwRcDo9SBAAA8ES04PAoRQAAAE/k4arkKMUoz07z0dFSd9SFrnUAAAAigspOXu3RaT4yNlYa3nEHVxYAAAjuJsWCggLp0qWL3HXXXRLsoxRJCwEAAEIy4JoxY4Z8//33EswK7DxcHjVclmVJ4b59ZtF1AACAoGtS3LJlizz55JNy/PHHSzArcuXh8pjaJydHtpzSx6wztQ8AAAjKGq4rr7xSpk6dKsnJyRLMCkgLAQAAQrGG6/nnn5f8/Hy5/PLLZf78+UfdPy8vzyy2zMxMc1tcXGwWp9lNhKbZ8EhaiKgIq8RrlVoPwHlUB1o2Wo6B+JzCCeVIOQYTrkfKMdyuyWI/jhs0AVd6erpMnjxZli5dKhERvuW2mj59uqkN83as3Nxcx88x58gxDx46JPmFRWY948B+SYvIKdGk6H4eEQkJjp9HdaAXaUZGhvkiRJI8lnLkeqwW+F5TjuF2TWZlZQV/wDVt2jSz2Hr37i2XXnqpdO7c2edjaIA2ceLEEjVcTZs2lZSUFKlVq5bj55wQv9vc1khKkiNduKRharKk1k107VOcnS0HjqzreUQm/vkYSn4JNLA2ZUTAVWGUozMoR8oxmHA9hk5ZxsfHB3/ANWHCBBk1apTrfsuWLWXt2rUya9Ysc//gwYOyZs0aeeutt+Tbb7/1eoy4uDizeNJCDUTB2jVveuuavDo6quRrua0H6jyqCy1HyohyDBZcj5RjMOF6DI2y9OeYVRZw1alTxyy2bdu2lXj8ggsukJNOOkkmTZokwUa7cjGXIgAACLk+XC1atChVTacBWZMmTSTY2MGW17kUo6Ol9siRrnUAAAAigkoGXNFepvZpdN90riwAABD8AdeKFSskWNkpIRRT+wAAgJANuEIh6am3yat16KmdGkJTQvia4gIAAFRfDKGrxLQ+GktFRZae2mdzt+5mcc/JBQAAwhcBVyUmro4h5QMAAPABAVcF2Dm4oj1HKAIAAHhBwFUBhUdquOgwDwAAfEHAVQEFrhouig8AABwdEUMl0kJQwwUAAHxBwFWJxKeeKSEAAAC8IQ9XpZoUvXSaj4qSmkOGuNYBAAAIuBzuNB8ZFydNHnmYKwsAALjQJlaJtBA0KQIAAF8QcFWm0zx5uAAAgA8IuCrRaT7aS6b54uxs2dSho1l0HQAAgICrUqMUyTQPAACOjoCrAgqONCl6TlwNAADgDQFXBdBpHgAA+IOAqwKYSxEAAPiDgKsCmEsRAAD4g4CrAug0DwAA/EGm+UpNXu0lXo2KkqT+p7rWAQAACLgcnktRp/ZpNns2VxYAAHChSbESneZjvNVwAQAAeCBiqERaCKb2AQAAviDgqtTUPqWbFHU6nx9P7GYWpvYBAAAmZqAYKi46ynu8auXkUKwAAMCFGq5KoEkRAAD4goCrEug0DwAAfEHAVQnUcAEAAF8QcFVCTBl9uAAAANwRMVSCt1GKAAAAnhilWAmx0V7i1chISezZ07UOAABAwOVwk2JkfLw0n/ciVxYAAHChCqYSYunDBQAAfEDAVQkx0fThAgAAR0fA5XCTok7n89PJp5iFqX0AAICiD5fDAZcq2r+fqwsAALhQw+X0KEUAAAAPRAyVQKd5AADgCwKuSiDTPAAA8AUBVyXERDFKEQAAHB0BVyVQwwUAAHzBKMUATO0T37mzax0AAICAKwBT+7R88w2uLAAA4ELAVQn04QKAqlNUVCQFBQXV7iMoLi427ys3N1ciaSmp8rKMjY115HMg4KoE8nABwLFnWZbs3r1bDhw4UG3fnwYKWVlZEhHB4KyqLksNtlq2bGkCr2oTcN10003y0EMPldj22GOPyTXXXCOhkoerOCdHfjlzuFlv9d9FEpmQUAVnBgDVlx1spaamSmJiYrULSjRIKCwslOjo6Gr33kKtLDVY27Vrl/zxxx/SrFmzSn0eQRVwbdy4UaZNmybjxo1zbatRo4aE1ChFy5KCXbtc6wAAZ5sR7WCrfv361bJoCbiCqyxTUlJM0KXHiYmJqfC5BNUwuk2bNslf/vIXSU5Odi3x8fESrEgLAQDHlt1nS2u2gGPBbkrUYL8ygibgOnTokPz222/Stm1bCRV0mgeAqkFTG0LtWguaJsUff/zRVP098cQT8u6770q9evVMn65Ro0aV+Zy8vDyz2DIzM11trro4Tc/PPdjS++7b7NcusR6A86gOtGzszoygHKsa12PolKP9Gt7+/w1mp59+umkGffXVV0s99tprr8mVV15p+qbZrTr2e3N/j9px+84775RLL73U8fP78MMP5aWXXpIXX3yx1GPbt2+XVq1ayS+//CItWrTw+ZgrVqww77uq/5+3vJTlgAEDpH///jJlyhTp3bu3fPDBB9KgQYMyn29f157vxZ/3FlQBV1RUlDRt2lTee+89Wb16tVx22WWm2njkyJFenzN9+nSZOnVqqe3p6elmCKjTctyOGR0ZIWlpaaX2sXJySpxHBJ3mvdKLNCMjw1zEDHuuOMrRGZRj6JSjNinq62h/Gl1CxXnnnSc333yzHDx4sFRXGQ3C9HdO+xnpe9Lys5uvPGtX7PfudJled9118vbbb3s9tr3N3zK330NVfk5WGWVpB1B6nY4fP14mTZokzz//vNdj6Pnrvnv37i3Vh0tHPwZ9wKWd43WxLV68WPbv3y81a9Y097t06WL6dD355JNlBlyTJ0+WiRMnlqjh0oBNO7jVqlXL8XNOiN/tWo+NjjJ/rXgqzs4We6Cynkck/Qy80otXL35TRuSZqTDK0RmUY+iUo/4xrT9yGpzoEkoB1/XXXy/Lli2TESNGuLZrAKa1S9qy4/l+vHXQ1nJ1+n2//PLLpjtPp06dvD5uv56/Za6VKO7Pr0qeZanXqV2Wo0ePNjVd2jFeRyJ60n10Xx2k4Rks+9PPvMr6cE2YMEG+//5719KzZ09XsGXr0KGDKYCyxMXFmcDKfVFaMIFY3KNj7TDvdb+oKIlt09osuh6oc6kOi33Bs1COwXANcD2GTjnqa4Taoj/WQ4cONbVI7tsXLlwoderUkYEDB7q2Kfdbz+26aK2UBgnaxKej+bViYufOna7Ht27dKsOGDTOP6T6aYqmsc5s1a5acf/75JY6tv9G1a9eWRo0amYDM81zcl/JeSz311FPSpEkT8xs9duxYyc/PN49pLdNtt91mnqOBizZbzp492/VcbfbTlqwhQ4ZIUlKSdOvWzQSs9uN6LWgzqFbQ6LF1P+0Lbj+uTZrdu3c3sYXGE/PmzfNalnreev7PPfdcuZ9hWddj0AdceoFpIdvLjBkzTHuquw0bNkjHjh0lGMV5m0dRCzQhQVovWmQWcnABQODpD3d2fmGVLP70I9M+ydplxj07/uuvv26CHbs2yFeaPumtt94yAce6detMEPO3v/3N1DJq32YNPjTIW79+vTz99NMmONNaNE/79u2TtWvXSt++fV3btPVJW5205m3RokWugMsbX15Lm0yXLFki77zzjnm/L7zwgtmuv/u634IFC0z/MO1GdPXVV5ucV7b777/fbNe0UVoDpwGbe7+pe++91wRpX3zxhXmeBnBK14cPHy4XXHCBqdR54IEHTL9w7avljb5/PcdAqvp6viPOPPNMueeee+TBBx+Uc845x7xxjUZXrVolwYgRigAQHHIKiqTTHR9WyWtvvHuIJMb69lOqTYlXXHGFLF++XAYPHmy6wWgA4O/vnAZJ8+fPN8+1AyUNvLRLjQZJGtBpH2MNRLT2pn379qbWyVvz13fffWea27RDvk278mgwdPLJJ5v7mpBca+e80dc72ms9/vjjcvzxx5tl0KBB8s0335jtGkDNmTNHevTo4eompIMCtm3bJscdd5yrzC688EKzrgGTBn/6eg0bNjTbtO9Vnz59zLoGZhrw2bVqWmv4z3/+0/TBatOmjdxwww3mPL29F60l+/LLL00AHagRsEETcGm1n0brt99+u6uK8ZVXXjFNjcEo2lvSUwAAyqCDwLQWSn/rNODS2h3tM+Tv79yWLVtMLY+OrrPpyH4NdnQAmnYS13X3xOEajHijwUvdunVdQYYmldUBX/qbbOvVq1eZ5/LTTz+V+VrapKe0Oc+9dSvnyOAyDaY++eQTE2hpLZTW1CmtrbO5j4rU96jcsxN4Pm4/9sMPP5iKG31vNj2uBl7eaN5PDcy0Y7yuV+uAS5111llmCQXxMd4DLp3aZ9u555r1lm+8QbMiAARYQkyUqWmqqtf2hzYrXn755aYW6Y033ig39dHReEtLZGdD96ep072Jzg683Du6lzeHoO53tNfyfL51ZP9//etfMnfuXBkzZoxcdNFF8sgjj0jr1q1L7OutqdX99Twftx/ToPPiiy+WW2+9tUSm+bIyxfvbpBvyAVcoiYsu48OxLMn/eatrHQAQWPpD6muzXlXTmi0NALQp7qOPPjJ9i/ylIwq1s/aaNWtM/ym7mXHz5s2m37MGF1rzpAnFtbO5HdxkZ2eb5j532jSnGQLspjTtKN+4cWNT22T3odb+1GXR2quyXuvss88u933Mnj3bNCna+2mTnlP0vLSpVmvA7IDrmWeeMeV0yy23lNpft2vQFcjpomgXc7jTPAAAZdEalnPPPdd0LO/cubO0a9fO78LSpjPt13TNNdeYoEJTKGkzno4E1IBO+0lp0KQjDbX5UQM77dvkre9S165dTbClTZE2TbekNUOaD1P7W2nfqbL481qeNNjTTvk6slCbFrXmT4Me9ybFirrqqqvkq6++MuWsHfJ1NKj29/KW9kFpp3wti0DOYEDUUEHxflYjAwCgtBnx66+/rlRz4rPPPmsGmOkxtOO89g/ToEWb7zRo0f5hO3bsMJ3BNbGnjvbTwWmeNGXCKaecYoIrm3Yu19F9OspPAzgdEVkWf17Lk9ZuaVCkQeeNN95oUkBoB3q7L1dlNG/e3ARz2o/rhBNOMO9JR1+WVeafffaZSQ0RSBFWKM2NcBQ64kOrQzXTcSASn/7rzW/ktS93mvVhnRvKUxf/2anQPfHp5m6Ht7df/xWJT8vpM6CdNTV5rD95TEA5BgLXY+iUoyY+1VFsOqrOn6SToUR/lt37HQWajnDUHFQff/yxhGNZ5uXlmdpBbdLUQM2fa86fuINfugqihgsA4CQdDaij/bTWSUfX6a3e18Ue8ReoLPhaQ6VNguFo3rx5pkbOW7DlpNDoZRiE6MMFAHCSJgbV/kveamXsvFOB6lemnek1F6a3yaurs/z8fJMnTAcxBBoBl9M1XDrstFEj1zoAAL7QaXSqoklR6WhHe8RjOImNjTX96Y4FAq4ATO3TZln1awcHAAAVRx+uCopjlCIAAPARAVcF0YcLAAD4ioDL4T5cxTp89JxzzaLrAAAA9OFyuoaruFhyv//etQ4AAEANVwWRhwsAAPiKgOsYzRAPAEAg6fyKmi1++/btJp2E3voiPT1dLr74YklOTjbzIurciZp9XU2ePNlMMo3Ko0mxgpLiCLgAAMFh+fLl8vvvv8vAgQN9DrRsOr+gBmg6Efb+/ftlzJgxJgeYzomoEz6feOKJZt7G+vXrB+z8wwE1XBVUI45YFQAQHDRL/LXXXuv383RKn6VLl5parI4dO5qJrKdOnWrmV1QaZJ1++unUcjmAgKuCkgi4AAA+spv57rrrLklKSvIpONI5/tq3by+JiYnSt29f2bBhg9f9Nm/eLGvXrpXhw4eXekznYNTX9ba88MILUqdOHVm4cKGZmNmmjx06dKhEU+XTTz/NZ11JVNMEoIYrqm7dih4WAFBBxdnZZT8YFSWRcXG+7RsZKZHx8UfdNzIx0e9zXLdunXz77bdH3e+jjz6S8ePHy1NPPSUDBgwwcx0OGzZMfv31V4lzex/qgw8+kG7duklCQkKp4/Tp08f00fJGJ8fWY7kHajqtkL6m1mrZtNZLa8I0sNMAEBVDwOWHIstyrSfGRpX5BWz3+WcV/DgAABW1uVv3Mh9L6n+qNHPr/P1Tn75i5eR43TexZ09pPu/PSZx/HniGFO3fX2q/jj9u8vscJ06cKK1btz7qflqjpDVLl112mbk/ffp0iYqKkgMHDkiDBg1K7KtzAWpzYFkTU2tneF/deOON8sUXX5jA0D0wa968udlGwFVxNCn6ITuvyLVOkyIAwF/NmjXzab8tW7aYWiv3SZZnzJhRKthSaWlpUq9ePa/HWb16tWk29La8/PLLJfa97bbb5NFHH5VnnnlGTjjhhBKPadC2e/duH98lvKGGyw8Hcgpc60ztAwDBpf36r8p+MKpkq0S7T1eXvW9kybqINh8vFadojZMvdJSgNu/5qriMRNs9evQwNWDeuNd83XzzzfLggw/KrFmzZPTo0aX21do1VA4Blx/2HDycl8TuVOiNTuezY9x4s970madL9AMAAASOP32qArWvU7Tpzr2TvAZUbdu2NR3d+/XrV2Lfhg0byt69e70eJz4+Xlq0aFHuaz3wwAMm2Jo7d67Jx+XNvn37JDU1tULvBYcRcPlh78H8o+9UXCzZdts3U/sAACrg6quvlqFDh5q8WjpCUdM2FBQUSM+ePUvtq9uef/75CpWzdoafMmWKybelr7dnz55SNWCaBPWXX34p0cQJ/9GHyw/JNUuODAEAIBB0ZODjjz8ut99+u3Tq1EnWrFkj77//vqmx8jR48GD55ptv5ODBg36/znvvvWcCOe0flpKSUmKxacqJ4447Tjp37lzp9xXOCLj8cN3pbSQ5KUbmXtYjcJ8IAKDa0WY97ZN1tOY9d2PHjpVt27ZJdna2ySRfVsDTqlUr08z49ttv+/1aWpOm+3pbbK+88opceeWVPp83vCPg8sOwzg1l0bgTpF/bPyN/AAD8paMBa9So4XW54oor/D7eLbfcYkYXOk37bmnN2oQJExw/drihDxcAAMfYpZdeKiNGjPD6mGai95f29Xr22WfN5NW67pSZM2eaPl51SehdaQRcAAAcY7Vq1TKLN9qcV1hY6PcxtenPaZpwFc4g4AqACC/TKwAAgPBFwOUwzdfSYcN6pw8LAABCGJ3mAQAhx58s7EAwXGsEXACAkGFPjaOpEoBjIT8/35HpjWhSdFhxXp7svO46s97k0UclMo5kqQDgFP3R04mXdcJmlZiYWOZUa6HK7jSv8ylWt/cWamWpUyqlp6eb60yPURkEXE4rKpJDn6x0rQMAnKVzByo76KqOQYL+0EdGRhJwBUFZ6nObNWtW6c+CgAsAEFL0h0+nmtHJlHVamupGAwSdjLp+/frmxx5VW5axsbGOfA4EXACAkG1erGy/mmANErSvms6bSMBVfcqS0BkAACDACLgAAAACjIALAAAgwKKrY3KyzMzMgLUFZ2VlldsWXJydLQePjE7U84iswHxY4cCXsgTlyPUYWvheU47hdk1mHok3fEmOWq0CLi1U1bRpUwkKxx1X1WcAAACOQfxRu3btcveJsKrR/Agaye7atUtq1qwZkNwlGslqMLdjx44yZ3kHZXkscU1SjsGE65FyDLdr0rIsE2w1atToqDVo1aqGS99skyZNAv46+qERcFGWwYRrknIMJlyPlGM4XZO1j1KzZaPzDAAAQIARcAEAAAQYAZcf4uLi5M477zS3qBzK0hmUI+UYTLgeKcdgExdEv9vVqtM8AABAMKKGCwAAIMAIuAAAAAKMgAsAACDACLh8lJOTI2PHjjVJVRs0aCDTpk0L7CdTTWzbtk3OOussk6ekZcuWMn36dJOgVq1fv1569Ohhplw48cQT5Ysvvqjq0w0JV155pZx22mmu+7/88osMGDDAlGP79u1l0aJFVXp+wU67rWon2oYNG0rdunVl/Pjxkpubax6jLH23e/duOfvss6VOnTrSrFkz8922UY5Hl5+fL506dZIVK1b4XG5LliyRjh07SkJCgpx66qny008/+fGJhVdZLl++XHr27Ck1atSQE044oVRZzp8/X1q0aCGJiYnmN+p///tf4E9UO83j6K6++mqre/fu1saNG62PPvrIqlOnjvXyyy9TdOXIy8uzjj/+eOviiy+2fv75Z+uDDz6wUlJSrCeffNLKysqyGjRoYN16663W9u3brVtuucVKTk62MjMzKdNyrFq1yoqIiLD69+9v7hcVFVldunSxLrvsMlOOjz32mBUfH2/98ssvlGMZ7r//fis1NdVavny59fXXX1vt2rWz/vWvf1GWfho6dKg1YMAAa9OmTdbKlSut4447znr++ecpRx/k5uZa55xzjg5YM9ehL9/lX3/91UpMTLSeeuopa9u2bdbo0aOtjh07mueFs1wvZbllyxZTVv/5z39MWelvTkxMjPXtt9+ax9esWWPKdsGCBWbfwYMHmyXQCLh8kJ2dbT4c+8NUU6dOtfr16xfIzybk6X/CcXFxpvxs9957r9WnTx/rueees1q0aGEVFxeb7XrbqlUra86cOVV4xsEfwHbq1Mlcd3bAtWzZMvMfy8GDB1376Y/g7bffXoVnGrz0x0mD/vnz57u26frAgQMpSz/pdbdw4ULX/YkTJ1ojR46kHI/ihx9+sLp27WoW9yDhaN/lO+64w/W9V4cOHbKSkpKsjz/+2ApXP5RRlnfffbd16qmnlth30KBB5g98pcHqmDFjXI9pgKt/yG7dujWg50uTog++/vprKSgokFNOOcW1rU+fPrJu3TqfZggPVx06dJCFCxea6m+bznF56NAhWbNmjfTt29c156XeavnSrFg2bbLRptfTTz/dtU3LsVu3bpKUlFTi2qQcvfvhhx9kz549MmLECNe2iy66SJYuXUpZ+kmvxVdeeUXy8vJMc8yHH35oug5wTZZv5cqVptlw9erVJbYfrdz08X79+rke06awcO+KsbKMsjz//PPlwQcfLLHN/u3xVpbNmzeXxo0bB7wsCbh8sHPnTklOTpbY2FjXNp2oUvt97N27N5CfT0hLSUmRQYMGue7rf8zPPfecCRi0TPUCd6dl+vvvv1fBmQa/H3/8UWbNmiUPPfRQie2Uo3+2bt1qrstVq1aZHyvtV3jTTTeZPiCUpX+0D4z+0GkfmeOOO06KiorktttuoxyP4h//+If5HrsHVr58l7k+fS/Ldu3amf7Btk2bNpk+XfYfq1VVltVq8upA0cDKM0utfV870+Po9D9jrUnQAHXSpEly8cUXey1TyrM0rUXVjvJ33XWXpKam+nRtUo7eHTx40CxTpkyRhx9+2Gy74oorzPVJWfpn9OjRJmB988035Y8//pD77rtP0tPTKccKOtr1x/VZMWlpaTJy5Ejp3r27DB8+vErLkoDLBzpixB7F5F5bo9yby+CdjkrUEZ7avPj++++bv4bLKlPKs7Rnn33WNGnraDpv16Y2kVGOvomOjpbs7GwTbNkjPWfOnCmjRo2SSy+91ARjlOXRff755/Lpp5+amgL9PqusrCzzR9WwYcO4JivgaN/lsv7PrFWrVkVeLizs27fPtLJouekfBnYXlqr6/aFJ0Qda9ag1M/qjZ9OqR/3Q6tevH8jPJ+RpzYHWZr3++uuyYMECGThwoKtM9a9id1qmntW8EHn55Zdlw4YNJoWBDsHXmgRtytF1TVFCOfpOU0EoHXLv3tdQ//PVxyhL3+zYscN0s7CDLaVNtJoGRmthKUf/He3/RP7P9L9mq3///pKRkSHLli0r8dtSVWVJwOWDrl27mr+MP/vsM9c27QPSq1cvV8QM766//np55513TA4UuzpXnXTSSSZosAcd6K3+xazbUZJ2TNY+CDp4Qxftt6D9E3Rdy0vzmdmdQe1rk3Is+7scExNjysymZav59XQQB2Xpm9atW5vaGPfcRdrPUPu5nnzyyZRjBRztu6y3et+m++kfYnzXS9M/oP7617+aJkL9ndHr1Z1nWeofChpwBbwsAzoGshoZP368ycP1/fffm+G7mofrtddeq+rTCmqa60Qvsccff9xKT093Lfv27bMyMjJM3i0dpqv5ZfRWcyNpfi6U78477yyRh0tz8di5ezTfTEJCgllH2d/lNm3aWJ999pm1bt06k4fr5ptvpiz91LdvXzPUXofmawqYli1bWv/4xz8oRz945uEq77us+bg0PdETTzxhtul+nTt3Dvs8XN7K8r777jMpiT7//PMSvz3278vq1avN42+++aZJBTFs2DCzBBoBl4/0g7rkkktM3hNN2DljxozAfjLVwKRJk8yXwHNp3ry5KyDT/CmxsbFWt27dzI8f/Au41ObNm82Pn/4H0r59e2vx4sUU41Hymd1www1WvXr1rJo1a1rjxo0zyRMpS//oD9iFF15oyrFx48bWTTfd5Mq5xzXpf5DgS7lp3jP9A0Ef1/8DNKE0Spdljx49vP72uOfe0pyPTZs2NUHsiBEjrLS0NCvQIo6cKAAAAAKEPlwAAAABRsAFAAAQYARcAAAAAUbABQAAEGAEXAAAAAFGwAUAABBgBFwAAAABRsAFoFqaPHmyzJkzx6y3aNHCTMNlLzrh71lnnWUmX/aVPm/FihU+7Zufn2+mEXKf+gZAeCPgAlDtbN68Wd5++20ZM2aMa9tjjz0m6enpZlLbzz//XAoLC83E6oGgcwpeddVV8s9//jMgxwcQegi4AFQ7999/v1x++eVm0nlbjRo1JDk5WVJSUuT444+Xu+++Wz755BM5cOBAQM5Bg7n33ntPfvvtt4AcH0BoIeACEDK2b99umvbuuusuSUpKkmuvvbbUPllZWfLSSy/J+eefX+6xEhISJD4+XuLi4sz94uJiefDBB6Vt27bmsdNPP12+//77Es9ZuXKldOnSxTzvtNNOM+ejXnjhBendu7eMGDHCnNerr75qbocNG+Zq1gQQ3gi4AIScdevWybfffis33HBDqceWLVsmDRo0kGbNmpX5/NzcXNPEeMcdd5jgSk2dOlVmzpwpjz76qDl28+bNZciQIXLw4EHX82bNmiWPPPKIrF+/Xvbs2WP6idnWrl1r+m1t2LBBBg0aZLb17dtXlixZ4vC7BxCKCLgAhJyJEydK69atzeLp66+/lo4dO5baPmHCBKlTp47Url1bEhMTZd68eaaZUVmWJY8//rjcc889plZKa7mefvppiYqKMvvZbrnlFlPz1alTJ7niiivkm2++KdFvSwOwdu3aSf369c02rQ378ssvzfEBhLc/OzgAQIgor/ZKO8XXq1ev1Hbts3XOOeeY4Edrp+bOnSvXXXedtGzZUnr27Cn79u0zzYK2mJgY6dGjh/z444+ubR06dHCta/CWk5Pjuq/9w+zaMvdt2jl/7969Zh1A+CLgAhByNBgqj/bH8qQBj6aHUHaQpc19H3zwgVlXnjVRehwNmNxrsdy57+/tnLSGDAAUTYoAqpWGDRuaGiVfaJCkOblSU1Olbt26smbNGtdjGmhpXzFvzZO+0lozDbrsJkYA4YsaLgDVitZWaX8sT9r5XZsS1f79++WVV16RLVu2yIUXXmhGPmofL+1E36RJE9MP64EHHpC8vLyjjnYsz8aNG01Hej0+gPBGwAWgWunXr59JDaHBlHZ+t2kKCTuNhHaa147vmhxVO7bboxS1lkv7dWlgpiMMP/30U5O3q6I+++wz0wkfACIshs8AqGbGjRtnOtbffvvtVXYOWjumtWU6SlFTTAAIb/ThAlDt3HTTTTJ//vwSHd6PNU0nceaZZxJsATCo4QJQLWlOrDZt2pgpfo41nby6V69e8uGHH5okrABAwAUAABBgNCkCAAAEGAEXAABAgBFwAQAABBgBFwAAQIARcAEAAAQYARcAAECAEXABAAAEGAEXAABAgBFwAQAASGD9P9RQ3QDGsNwrAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 绘制局域势\n", "fig, ax = plt.subplots()\n", "ax.plot(kb.r, kb.V_loc, label='V_loc (d channel)')\n", "rc_loc = invert_results[loc_channel].rc\n", "ax.axvline(rc_loc, color='tab:red', linestyle='--', label=f'r_c (l={loc_channel})')\n", "ax.set_xlabel('r (Bohr)')\n", "ax.set_ylabel('V_loc(r) (Ha)')\n", "ax.set_title('Al Kleinman-Bylander local potential')\n", "ax.legend()\n", "ax.grid(True, alpha=0.3)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "7737c0ab", "metadata": {}, "source": [ "## 可视化:各通道投影子 $\\beta_l(r)$\n", "\n", "根据 $\\beta_l \\propto (V_l - V_{\\mathrm{loc}})\\phi_l$,不同角动量的投影子应该在局域势截断半径附近集中并快速衰减。下图比较 s、p 通道的 $\\beta_l(r)$,方便检查其局域性与归一化情况。\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "7dc45d47", "metadata": { "execution": { "iopub.execute_input": "2025-12-15T08:11:03.040438Z", "iopub.status.busy": "2025-12-15T08:11:03.040364Z", "iopub.status.idle": "2025-12-15T08:11:03.146370Z", "shell.execute_reply": "2025-12-15T08:11:03.146023Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAGuCAYAAAD7+IBiAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAX2lJREFUeJzt3Qd4k1XfBvA73S0thQ5mKZS9994oWxmKgKh8KoK+AqKgqLgRBcervuIAJ6AiKqLIHiJ77z3K3qMtlO6Z7/qfkJqWlq4kT57k/l3XY56kGScnkd4902A0Go0gIiIiIl1y07oARERERFR0DHNEREREOsYwR0RERKRjDHNEREREOsYwR0RERKRjDHNEREREOsYwR0RERKRjDHNEREREOsYwR0SUB66pTkR6wDBHpFOPPvooDAYDPvroo9t+1rlzZ3XkZc2aNeqxM2fOvO1n69atg7+/Pxo2bIirV6/i9OnT6r6Wh4eHB8LCwvDkk0/i+vXr0JK8BymTlNNaUlJS8Nxzz2H27NnQ2nvvvYdHHnkEemf+Hpm/c7b43PLy2GOPoUqVKur8yJEjqFy5MmJiYmz+ukT2wjBHpENxcXH4/fff0bRpU3z99ddWa0Fav349evfujZo1a2L16tUoU6ZM1s+eeOIJbN68WR3//PMP3nnnHSxZsgQPPvggtHTPPfeoMpUvX95qz3np0iV8+umnSE9Ph5YOHjyowtyUKVPgbGzxuRVE7dq10bdvX4waNcqur0tkSx42fXYisolff/1VXb7//vvo1q2bCld33313sZ5zw4YNKsjVrVsXK1asQKlSpbL9XFriWrdunXW9Y8eOcHNzUy2EJ0+eRNWqVaGF0NBQdTijV155RbXKVapUCc5Gy8/thRdeQLVq1fDSSy+hcePGmpSByJrYMkekQ99//z3uuusuFeAiIiLw1VdfFev5Nm7ciF69eqmu1b///vu2IJcX8/2kpTAvb731lmoNmTt3LqpXrw4fHx80b95cBVAz6XKTrtsZM2YgKChIBcejR4+qn82fPx8tW7ZUjwsODsbDDz+Mc+fOZXtszu46ee527drBz89PtS4+/fTTiI2NzVauLVu2qPoLCAhQoeKhhx7C+fPn1fNInYrHH388q3vO3AXdpUsX9byBgYGqhefAgQP5vg8Ju/fdd58qf8mSJVXZVq1adce63b9/PxYsWKC6si2fX+ph27ZtaNOmjToPDw9XrXeWkpKSVFCRn3l5eamW1g8++ACZmZlZ95H3NW7cuKz3M2bMmKzud6k/Ceu+vr7qOb755htcvHgR/fr1Q4kSJdT7kueztGnTJvUdkvqWckm4nzhxYrbXtJTzc8vZlW95mF27dk19JvJ5SZllKIG8rqWoqCgVgOW7KcfIkSORlpaW7T7Szdq9e/fb6o1It4xEpCuHDx+WPlXj3Llz1fWJEycaPT09jZcvX866T6dOndSRl9WrV6vnmDFjhnHDhg1Gf39/Y5s2bYxxcXG33ffUqVPqvq+//roxLS1NHcnJycajR48a27VrZ2zWrJkxMzMzz9d688031fOXLl3a+MUXXxgXLlxovOuuu4xeXl7GgwcPqvtIOeQ1qlatavzrr7/U9YyMDON3332nbn/wwQeNixcvNs6aNctYpUoVY/ny5Y2XLl3K9lgpp1i7dq3Rw8NDvcaCBQuMP//8s7FWrVrGli1bGlNTU9V9du3aZfT29lZ19McffxjnzZtnrF69urF27drqvclt8pyvvfaauq9YsWKF0d3d3di1a1fj/Pnzjb/++quxUaNG6r0dOHDgju9DnrdLly7GRYsWGVeuXGm89957jT4+Plllzo28dnh4eLbb5Pnc3NyMlStXNn7++efGv//+2zh48OBs3wd5vc6dOxtLlChh/PDDD1W5J0yYoOrk8ccfz3oueQ75DMaMGWNctmyZ+h6Yvxdly5Y1vvfee6qsPXv2VO9b3sOrr75qXLVqlXHo0KHqfps3b1bPtX//fqOvr6/x4YcfNi5fvlzdZ/jw4eo+M2fOzPY9kveQ2+cmz2V5yOPkdUeOHKl+Lp9Lw4YNjSEhIcbp06erMks55D1s374963116NBB3ef7779Xn7/Uu7x3eb+Wpk6davTz8zMmJSXl+RkQ6QXDHJHOvPjii8agoCBjSkqKun7mzBn1C/7dd98tdJh78sknjQEBAeq8Ro0adwxzuR1lypS5YyAxhzm575w5c7Juk1+g5cqVUyHN8he7+Re9SE9PV8/frVu3bM93/PhxFV7HjRuX7bHmckjAlPdirh/ze5Bf+j/++KO6PmDAAPX6EhDMNm7cqILi3r17bwseQkKrBBopl9n169fVZ3H//ffn+T6ioqLUbb/88kvWbTdu3DCOHTs2KwTmRsKnlNOS+fm//fbbrNvkM5OwMnr0aHVdwrLcZ/bs2dke+/bbb6vb9+3bp65LuMkZcMzfi5dffjnrtm3btqnbHnnkkazbrl69qm775JNP1HWp1x49emQL9RIqJcCPGjWqQGHOUnR0tArXEszMAfybb75R9//nn3+y3bd79+7Gu+++W51LkJT7SB1YftdCQ0Nve68S0uW+EliJ9I7drEQ6kpGRgR9//BEDBgxAamoq4uPjVXdep06dVFdYXl1aeZHJE9Jd99dff+H48eN3HBQ+YsQIbN++XR3SRfnbb7+p7jvpjjtx4sQdX0e6xAYNGpR1XbrhZHze1q1bs91PJnSYyaxDmU0r3Z+WZKxT27ZtsXbt2tteJzExUQ2qHzhwoBrPJxMY5JBuQSmruWtXupXl9b29vbMeK8956tQp1dWck9Tzrl27MHjwYLi7u2fdLt14ffr0ua0slu9DulabNGmC4cOHY+jQofjpp5/UbNmPP/4Y9erVy7POzp49i4oVK+b6M/nMzGTmcdmyZZGQkKCuS1mkjJb1LeS1zT/PrZyWpC7MzOPaLMdLhoSEqMsbN26oS+nWXLZsmarrQ4cOqe7hSZMmqfdpLldByXPI55ecnKwm+Xh6eqrbpVta6qNDhw5Zn6scMmZUur/l/w35Xkq3sny2lt81+Yxyku+EkC5wIr1jmCPSEZk9KjMtJbjJWC/zITNPZezR8uXLC/V8MsZJfvHK2K/Ro0fjhx9+yHM5jgoVKqixbnK0atVK/cJdvHixGi8nA/XvRB4r4cqSjK26cOFCtttkPJmZecmTcuXK3fZ8clvOMXDmx0ignTx5sgoBloeMbZNxX+ZxVZYzdfMjoUV6MgpaFsv3ISToDBs2DCtXrlShSuqjf//+agxYXqKjo9Vnm1c4tiTj9MxBXupAwpbclrOcwrKsOctpltvrylg5M8txbOLmzZvq/ZUuXRr169dXy7rIHwcSrAo70/qZZ55R4+D+/PPPbJ+RfGbyfcn5uY4fP16NiZO6lJ/Le8/5XZP6zmu8p9Qzkd5xNiuRjsjAehm4Lq1zluQXpgSy6dOnq4BWUNJ6Y26dkpmxMvlBJgtIK4y0gOVHwkCNGjXUYP07kV/EOV2+fDmrhSc3EgzM98tJQlluj5VJCRI0xo4diyFDhuRaXvNz5xakJJzm1lolv/jleQtTFksSSmSpk//9738qVM6bN08t7fLaa6/lOXlFWvTMLV+FIe9N6ltarSwDnTnI5lfWopBAJeFL/hCQiTkSBuU7mVcYzctnn32mvsOzZs1SfzTk/Azku/bzzz/n+lhpoZb3LmE253vP7ftnnrQj9Uykd2yZI9IJCR+LFi1SIaV9+/bZDul6kmAmYSRna1dBycxFCYnSvSVrx+WcAZgb+YUYGRmZ77IkEkqkK8yyO1TKeqflVGQGrHTx5fzlLV2h0pUq7zkn6XKUMCZdtOZWRDkaNGigwpM8Tki3tLRySle12e7du3Hvvfeq7lRzy465VUmeV7pKZUkY6c6zbJFauHBhrmUxk+eTblCZgSqBUMoiM3zr1KmDK1eu5Pk4qVOZXVtY0u0tZZRucEvmFlf5vlibtKRJiJPZruYAJ38YSBdrQbv+pdVSQvizzz6L//u//7vt5/KZySxmCcaWn610G0sIlFa6rl27qpm88pmYyWe8dOnS257P/P+JzNYl0juGOSKdkKAlASu3Fich673JL3Hpgi2qZs2a4Y033sCOHTtu6zqVYCFjkuSQNenmzJmDnj17qkAnLTP5kbFv0o0rv2jlcTKeSl4rLzLuSxbLlV/yEi7lF7KMN5MxUhLy5Bd/bqSLVbqbZQkTCWwSgGUMlYQLeX/i1VdfVTsAyO1SHmkpk/LJz+X5pYXPPE5L3q/5eSW49ujRQ3VNy1IrEkalFej111/P831IeJPnk10IJAxKl/jLL7+sWuikqzovsnRGzmU3CroYrwQ6WdLkww8/VO/7zTffxNtvv62+I3cap1dUUm/SlSxjMOX9SSuv1L+8bwnu+ZEuWfljRJaEeeCBB7K+Z+ZDQrMsSSJLisjnIy138jpS77IEizxOgrJ8HvL5yALX8v+BfGckYObWlSrjNeUPGKkrIt3TegYGERVM/fr1jfXq1bvjfWQWZ8WKFY3t27cv8NIkOclsTVmmxGAwqOUfcpvNKrNnZRanzDSVZSgKMpv1hx9+MFaqVEktySHLhuzZsyfrPnea2SjLgLRq1UotJSKvKctRnDt37o6PleU4ZFarLJdRqlQpNdNSZmVaWr9+vZotKfeRWbNPPPGE8cqVK1k/Hz9+fNaSKrIci1i3bp0quzxGZgH3798/a3mVO70PuS4zd2XJD3lsgwYNspbsyIvMdJXnsix3Xs8vMzUfffTRrOuJiYlqRqrUt8z8rVmzpvGjjz7KNhM352MsvxdyaVn23L4rcpt8tubZp7IsSXBwsFoSRep706ZNxmHDhqnlVWRm651ms5rP8zrM5ZHldx577DE1O1W+R3Xr1lVLjFiKj483Pv300+q7Ip+RLJHywgsv3DabtV+/fsaBAwfe8TMg0guD/EfrQElEzku6FGXxWFv9UyNjrGScn8x8dbadIGQWpmx3JS1eZD0yiUha+WRmdqNGjVi1pHvsZiUiXZJwKN2jMjZMxmk540B2Gecn76+o4yAp73qVrnsGOXIWnM1KRLok67DJODRZR0xmieZcjsIZSNiQMWEyxi7nDGYqGpkcI+Mo9+zZwyokp8FuViIiIiIdc74/ZYmIiIhcCMMcERERkY4xzBERERHpGCdA5EFWLZftb2SWXM59CImIiIhsMUtfFmLPbT/rO2GYy4MEuUqVKlnr8yEiIiIqENm6LiwsrGB3ZpjLm3l/wTNnzqgNnsl2LaCy56gs9uqMS0s4EtY169nZ8DvNenY2N27cUAtamzNIQbFlLg/mrtWSJUuqg2z3j7Fs7C51zDBnW6xr+2A92w/rmvXsjN9pUdjhXWwKISIiItIxhjkiIiIiHWOYIyIiItIxjpkjIiKibDIyMpCWlpY1jkvOZXwzxzYXj6enJ9zd3WFtDHNERESUtc7Z5cuX1axKy9sk0Mn6Z1x3tfhkhYxy5cpZtS4Z5oiIiEgxB7kyZcrAz89PBQ4Jc+np6fDw8GCYKwapx8TERFy9elVdL1++PKyFYY6IiIhU16o5yAUHB2cLIQxz1uHr66suJdBJPVury5UTIIiIiChrjJy0yJHtmOvXXN/WwDBHREREWTguTn/1yzBHREREpGO6DHOpqamoW7cu1qxZk+d9du3ahebNm8PHxwdNmjTB1q1b7VpGIiIisq3Tp0+rli65LAiZlfviiy8iKCgIpUuXxtixY9VYwTtZvnw5hg4desf7REVFoVGjRkhKSoIWdBfmUlJS8PDDD+Pw4cN53ic+Ph69e/dGz549cfToUXV57733qmnVeiKDTpPTMnAzOQ0JKenqPD0jU91OREREhfPxxx9j7ty5+Pvvv7Fs2TL8+eef+OCDD/K8v4xrGzNmDF5//fU7Pm9ISIjKHVOmTNHkI9HVbNZDhw6pIJcf+aBkxsikSZNUYp88eTJ+++03dfuwYcOghbSMTFxPSEVUfCqiE1IQHS/nKYhOSEV0/K3rCamISUhBYkoGktJMR265zc0AlPbzQukSXggq4YWwUr6oGloCVUP9UatcAKqGlOCYByIiohw+//xzvPXWW2jatKm6/vbbb+O1117DhAkTkJs5c+agevXqqFmzJvLz5JNPolmzZnjppZdQokQJ2JOuwty6devQpUsXFdL8/f3zvN+WLVvQvn37rEAjl23btlVdrdYIc9I6Fp+SjrjkdMQmpanLuOQ0XE80hbVrcSkqqKmwdiu0XU+03qyVTCNMITAhVV3fluPnEvCahpdG66pB6FGvHCoFcWYSERE5tzVr1qiMkJsZM2agR48eOHPmDDp06JB1e7t27XDu3DlcunQp13Xfpk2bhqeffjrr+r59+9T13bt3q9a40aNHq25bERERgVq1aqnGo8cffxz2pKsw95///KdA9zt//jwaNGiQ7bYKFSrg4MGDhX7N7v9bB3j6qZa19Awj0jIzkZyWiaKQFrWgEt4I8fdCsBwlvG9dynXvW5de8Pf2hK+nO3y83NSll4ebaqGTMmRkGpGSnqmCY8yt1rxzMYk4cS0eJ68l4PClm4hJSMXfh6+o453Fh1G3fEnc26g8BjWvhBB/7yKVnYiIXIsM6TH1EN1aZy7TfjNd5XdfYV+rXbt2uHbtWq4/CwgIUEFMVKxYMVs2EBcuXLgtzMXExGDbtm2YPXt21m3SO9i5c2f88ssvOHHiBPr166da+bp27ap+Lg1JK1asYJizBtk/zts7e2iR63camChj8eQwu3nzprq8HJsCN+/cF/Xz8XRDgI8nSvp4IMDHA6X8vFRQC7kVzOTSfF0upWvUTRJdEXlZLC4Y6u8FlLn9PqnpmTh4MRY7zlzHmqPXsPVUDA5duqmOT1YeUy11IzpEoEHFQDgCGYxq3iqGWNfOgN9p1rXev7vmIzE1HfXeXKFJWQ5O7A4/r/zbm8xjyI1Go9qhwnKx45zMGcDLyyvrcXIuZGeGnOPRJfzJXqpVqlTJ+pm04snkCdmOKywsTE2OsPx5/fr11Ti8O41tN9ev1HfO331F/V2oq5a5gpIZrBLoLElQM6+8nBsZtDhx4sTbbn++UxiaVC0LT3cD3N0M8HAzwM/LHf5e7vBwL0gwkw8mCRmJSYhKhF1U9AEq1vJHv1r+uJFUCWtP3MCCA1E4eDkBi/ZdUkeHqoF4onUF1C6jbResfHFjY2PVF5sbOLOunQG/06xrvZLB/vL9lVY486EV9foFmKJpLmN6ejrWrl2L/v3753q/qVOnokaNGlmTJM15wDLg5Xy/0vUqM14tZ7tKl6qMsZOxd9JtK7NcZZyc+bFyf9kS7U51Jz+Teo6OjlZh0ZL8PiwKpwxz0oQqH4IlaUK1bFrNSQY/jhs3LlvLXKVKldCwajm0r18ZeiWNdzUrV8CIu6Ba7L7bcBoL9l7E+pOx6hjYLAwv9qipunm1IF9oaUoPDQ1lmGNdOwV+p1nXeiWNILLqg7RwyRHg7q5ayERaWjo8Pe0XGQrazSrlNF+2atVKjWXLjYxvM/e4SVds1apV1fmVK1fUZXh4eNZzmclWW/L/s+XtkhUeeOAB/PHHH1i5cqWawSpBUcbOWbb05XwuS/IzabyQVkRpfLJkfnxhOWWYa926Nd555x3V2mPeJHjjxo1444038nyMdMPm7JoVbgaD04SMBmGl8b8HS+OZu2tg6qpI/LXnIubuPI/lBy/jtXvqYmDzME1mwRpu1bGz1LMjY12znp0Nv9PWI/8GS31aHiW8Zcy2UbWSSQhxtN0hLCc6+vr6qkkIeZFxc9JIs2HDBlSrVk3dJucS5Mxj5yzJGLrr169nPb+EQWmVk9mqL7/8sjqGDx+uljh55pln1P3k/rLn6p3qyVy3uf3eK+rvQaf57SkL9pnHvElqlr8uZF2Ys2fPqkvpDx84cKDWxXQI1UL98emDTTDv6TZqcsTN5HS8OG8fnvpxp1omhYiIyBk9/fTTammS7du3Y+fOnWppklGjRuV638aNG6sge+TIkawwuGTJEtWLJ5MfpBVw8+bNqkXQcgk187In9uQ0YU666WQ9GFGyZEksWrQIixcvVn3kS5cuVed3Ws4kLw72R4hVNaschIXPtMfLvWqrMYErDl1Bz0/XY+eZGK2LRkREZHXjx49XM1C7d++uZqAOGDAAL7zwQq73lfAmy5pJ652Q1rT58+erbloJerIhgTyHtNCZbdq0Cb169bL7J2cwcjuBXElzamBgINbuP4WO9avA2cl4uud+2YPIq/Eq2E3qVx8Ptgy3+evKeISrV6+qZml2s7KunQG/06xrPY+ZO3XqlOqqtBzLlbU0iQN2s9ra7Nmz8f3332PVqlX53lfG6suyaLK12J0aj/KqZ3Hjxg01iUImQkjDlMu1zFHx1KsQiPmj2qFX/XJIyzDi5T/2Y9KiQ8iUFYqJiIhc0KBBg9RyJJGRkfne94svvlATIYrSC1hcDHP5cKW/QUp4e+DLh5vi+W6mbUu+23AKL8zdqxYrJiIicjWenp747LPP1M5TdyLLjMjwLpkcoQWnnM1KRSdN6DLbNSzIFy/M3Yc/dl9QW5Z98XBT+HjmvngyERGRs+rRo4c67kSWGdmzZw+0wpY5ytV9TcLw9dBm8PZww6ojVzFy9i61uwQRERE5Foa5/LhSP2sOd9cpi5mPt1SB7p8jVzH21z1qb1giIiJyHAxzdEdtqgXjq6HN1AzXxfsv4aV5++645xwRERHZF8Mc5atzrTL4bEgTtTft7zvPY+qq46w1IiIiB8EwRwXSs355vNO/vjr/5O9jan9XIiIi0h7DXD4MrjxoLochLcMxooNp3ztZsmTnGdOedURERKQdhjkqlJd71UHXOmXVzFbZy/VqXDJrkIiINCG7LRgMBnVZGDL2W7bimjlzZr73Xb58OYYOHZrv/WRf1o4dO2oyrpxhjgpFxs19+mBj1CobgKj4FIyZsxvpXFSYiIh0tOXec889V6AtutLS0jBmzBi8/vrr+d63bt26qFixImbMmAF7Y5jLh4ttQ1fwnSIeaYoSXu7YcjJGjaEjIiJydBcuXFAtcgsWLECpUqXyvf+cOXNQvXp11Kxp2hkpP0899RQ++OADu7fOMcxRkVQL9cd7Axqq8y9Wn8DqI1dZk0REpJk1a9aoLtfcDnN36q5duxAeHo5t27YhMDAw3+ecNm0aBg8enHX9sccew7PPPovevXvDx8cH9evXx7p167J+3qlTJ8TFxamy2BPDHBVZn0YV8H9tKmdNiJBuVyIichLSupSaoM1RhJatdu3a4dq1a7keQ4YMUffp06ePCnahoaH5Pl9MTIwKfe3bt892+/Tp09GmTRscOHAAPXv2RL9+/RAbG6t+JsGxbdu2WLFiBeyJe7Pmg72sd/bqPXWw7VQMjlyOw4Q/9qstwOTLTEREOpeWCEyuoH4Petr7tV+5CHiVKNRDPD09ERISYrUi7N+/Xz1nRIRpFQezpk2bZo2hky7VefPm4eeff8bTTz+tbmvQoEG21jp7YMscFYu3hzs+GdxY7RCx8tAVtagwERGRvW3YsEGNg8vtkLBVWFevXkXp0qVva6CwbKlzc3NT4e7o0aNZt0mgvHz5MuyJLXNUbHXKl8S4brXw/rIjmLjwEFpXDUalID/WLBGRnnn6qRYyGcyfnp4ODw8P+/W8yGsXUvPmzbFnz55cf1bUFjuZ+ZqT1EPO+1je5u7uDntjmMsHewwL5smOVfH34StqIeEXf9+Hn0e0YncrEZHefwFKV6eMX3NLlxTj0L8UfXx8UKVKFas9X7ly5XD9+nUVZi1D7O7du7POMzIy1KQKGYtnOdauTJkysCd2s5LV1p/7aGAj+Hi6YfPJaHa3EhGRrjVu3FgFuSNHjmS7XSY3fPnllzhx4gTGjRunAt1DDz2UbfFg6Xq1J4Y5spoqISXw7N2mtXjeXXIY0ZzdSkREOhUQEKBmpspYPEvSCrd06VI0bNgQW7ZsweLFi1WroNmmTZvQq1cvu5aVYS4f3Ju1cIZ3iEDtcgG4kZiGdxcfLvo3k4iIKB/SrWo0GovUvSpbgMm6cXfy5JNP4pdffsl2m0yKWLhwIRISErB161Y0atQo62eylImMn+vSpYtdPzuGObIqT3c3TLm/gRpW8cfuC9gQGcUaJiIiXRo0aBDOnTuHyMjIAt1/6tSpmDBhgprlak8Mc2R1TcJLY2hr02LCr83fj5T0DNYyERHpjqenJz777DNMmjQp3/sePnwYZ8+exaOPPgp742xWsonxPWph6YHLOB2diJkbT+OpTtVY00REpDs9evRQhzBvC5abOnXq2H2xYDO2zOXDgWdhO7QAH0+81LO2Ov/sn+O4GpesdZGIiIicEsMc2cz9TSqiUVgg4lPS8eGyf1fHJiIixyUTCkhf9cswRzbj5mbAG33qqfO5O89j77kbrG0iIgceHyYSExO1LopTS7xVv+b6tgaOmSObala5NO5rUhF/7r6AiQsPYt7TbbkzBBGRA5JtqGQfU9mTVPj5+al/rzXZzssJGY1GFeSkfqWerbntF8Mc2dzLvWpj+cHL2HX2Bpbsv4x7GpZnrRMROSDZwkqYA505hMj+o7LcBsNc8UmQM9eztTDMkc2VLemDER2q4tNVkfhw+RF0r1dWrUdHRESORcJa+fLl1d6iaWlp6jYJctHR0QgODrb7+mnOxtPT06otcmYMc/ngXyHWMaJjVczeekYtVfLL9nNZ69AREZHjkcBhDh0S5iSEyJZVDHOOiRGb7MLf2wNj7q6hzj/9OxIJKemseSIiIitgmCO7ebBFOCoH+yEqPgXfrj/FmiciIrIChrl8cN6O9Xh5uOGF7rXU+dfrTqhQR0RERMXDMEd2dU+D8mgYFoiE1Ax8ufoEa5+IiKiYGObI7gsJm1vnZELElZvc5ouIiKg4GObywfURra9DjRC1mHBKeiamrz1pg1cgIiJyHQxzpMlyL+O61VTnc7adxdW4VH4KRERERcQwR5poWy0YLSOCkJphxKztl/kpEBERFRHDHGnWOje2q6l17q8DUbhwI4mfBBERUREwzOWDS5PYTptqwWhTNQjpmUZ8uYYzW4mIiIqCYY409eytXSHm7jiPczGJ/DSIiIgKiWGONCXj5lqEB7B1joiIqIgY5gowtotsa1jL8upy3s7zuBzLdeeIiIgKg2GONNckLADNK5dGakYmvl3PdeeIiIgKg2GOHMLIztXU5eytZxGTwHXniIiICophLh/sZLWPTjVDUL9iSSSlZWDmxlN2elUiIiL9Y5gjhxmbOKpzdXU+c9NpxCWnaV0kIiIiXWCYI4fRo145VAstgZvJ6fhpy1mti0NERKQLDHP5YT+r3bi5GTDyVuvcdxtOIjktw34vTkREpFMMc+RQ+jaugLDSvoiKT8Wv289pXRwiIiKHxzBHDsXT3Q1PdTLNbP1q7QmkpmdqXSQiIiKHxjBHDmdgszCEBnjjYmwyFuy9qHVxiIiIHBrDXH64A4Td+Xi6Y1i7CHX+9boTMBqN9i8EERGRTjDMkUN6qFU4/L09cOxKPNYcvaZ1cYiIiBwWwxw5pEBfTwxpWUmdT197QuviEBEROSyGuXxwZRLtDGsfAQ83A7aeisGeczc0LAkREZHjYpgjh1U+0Bf9GlfMGjtHREREt2OYI4f2ZMeq6nLZgcs4HZWgdXGIiIgcDsNcPjiZVVu1ygWgS61QZBqBbzec1Lg0REREjodhjhyeeRHhuTvOIyo+ReviEBERORSGOXJ4rSKC0CgsECnpmfhh8xmti0NERORQdBfmkpKSMGzYMAQEBKBs2bKYPHlynvcdMGAADAZDtmPRokWFej3OZtWefG7m1rkfNp9GYmq61kUiIiJyGB7QmfHjx2Pfvn3Ytm0bLly4gIEDByIiIgJDhgy57b6HDh3Cjz/+iJ49e2bdFhgYaOcSkzX0qFcOlYP9cCY6UXW3Ptq2CiuWiIhIby1z0ir33Xff4b///S/q1KmDrl27YuzYsZg2bdpt901LS8Px48fRtGlThISEZB2enp6alJ2Kx93NgOEdTDNbv1l/EukZmaxSIiIivYW5PXv2qJDWtm3brNvatWuH7du337Z/Z2RkJDIzM1Gtmql7rjhdfOQYBjYLQ3AJL5y/noSlBy5rXRwiIiKHoKtu1vPnz6vWNS8vr6zbKlSogOTkZERHR6ufmR0+fBilSpXC6NGjsXz5coSFhWHixIno1q1brs+dkpKiDrObN2+qSwmEcpBtSN1KEC9IHXu5GzC0dTj+t+o4vlp7Ar3rl2XYtlFdU9Gxnu2Hdc16djaZRfz3WVdhTkKbt7d3ttvM16UL1tKRI0fU/du0aYMxY8ZgwYIF6N27NzZv3ozmzZvf9txTpkxRYS+nmJgYXC2hq2rS3Rc3NjZWhQw3t/wbintVL4Hpa91w4OJNLNl5Ai3CS9qlnK5Y18R6dnT8TrOenU1sbGyRHqerlOLj46MCmiVza5qvr2+221955RU8++yz8Pf3V9cbNGiArVu34uuvv841zE2YMAHjxo3L1jJXqVIlBAUFoUyZMjZ6RyT/GEtXdmhoaIEChnwSg5rfwA9bzuC3/ddxT/PqrEQb1TUVDevZfljXrGdn42XR8+i0Ya5ixYqqO1XGzZknMsiMVgl5wcHB2e4rv7TMQc6sdu3aaoZrbqSFL2ern3B3c+MvPhuTz0rCRUEDxoiOVTF721msj4zC4ctxqFeBM5RtVddUNKxn+2Fds56dSVH/bdbVv+iNGzeGh4cHNm3alHXb+vXr0bJly9vGTg0fPhxDhw7Ndtvu3bvVLFjSt0pBfrinQXl1/vU6bvFFRESuTVdhzs/PD48++iief/55HDx4EKtXr8bUqVMxatQo9fOoqKisbte+ffvil19+wU8//YTTp0+rxYVlvJyMnyP9e7KjaZmSRfsu4VxMotbFISIi0oyuwpz46KOPULduXbRq1UotFCxj4wYNGqR+JmOB5syZkxXmvvrqK0yaNEl1r/7xxx9YunSpGgdXGFyZxDHVrxiIDjVCkJFpxHcbTmldHCIiIs0YjDkXaKOsCRCyW8SBUxdQr0oF1ooNBzBfvXpVTTIp7FiBjcej8PC3W+Hj6YZNL9+NoBJFGzjqKopT18R6dkT8TrOenc2NGzdQunRpNau1ZMmCr9bAf9FJt9pWC0b9iiWRnJap9mwlIiJyRQxz+TCAO0A4Kpn08p9Oph0+Zm06jaTUDK2LREREZHcMc6RrPeuVQ3iQH64npuG3Hee0Lg4REZHdMcyRrnm4u6l158Q3608iPYNbVRERkWthmCPdG9gsDMElvHD+ehIW77+kdXGIiIjsimGOdM/H0x2Pta2izr9ae1LtPUpEROQqGObIKQxtUxm+nu44dOkmNhyP0ro4REREdsMwlw/OZdWHUn5eeLClaUHo6WtPaF0cIiIiu2GYI6cxvENVuLsZsPF4NPafj9W6OERERHbBMEdOo2IpX/RtZNqtY/o6ts4REZFrYJgjp/JUJ9MyJUv3X8LpqASti0NERGRzDHP5MHDQnK7ULlcSd9Uug0wj8OWa41oXh4iIyOYY5sjpjL6rurr8Y9cFnItJ1Lo4RERENsUwR06naXhptK8egvRMI2e2EhGR02OYy4eBi5Po0jO3Wufm7jiPy7HJWheHiIjIZhjmyCm1qhqMlhFBSM3IxFec2UpERE6MYY6c1pi7aqjLn7eexbW4FK2LQ0REZBMMc/ngbFb9alc9GI0rlUJKeia+XX9S6+IQERHZBMMcOS2DwYAxd5vGzv245QxiElK1LhIREZHVMcyRU+tSqwzqVyyJxNQMfL/hlNbFISIisjqGuXxwzWD9t86N7mIaOzdr02nEJqZpXSQiIiKrYpgjp9e9blnULheAuJR0zmwlIiKnwzBHTs/NzYBx3Wqq8xkbT3NmKxERORWGOXIJ3eqWRaNKpZCUlsE9W4mIyKkwzBVgzBU5x+c4vnstdT57y1lcuJGkdZGIiIisgmGOXGrdudZVTbtCfLYqUuviEBERWQXDHLlW61wPU+vc3J3ncSoqQesiERERFRvDHLmUZpWD0KVWKDIyjfh45TGti0NERFRsDHPkcl7oUUtt07Zw70XsOXdD6+IQEREVi0fxHg6cPn0aK1euxPbt23H58mV1W5kyZdC0aVP06tULERERxX0JIquqVyEQ9zcJw7xd5/HOokOY+582nOhCRESu1zK3Zs0a9OjRA7Vr18Y333yDjIwM1KlTB/Xr14eHhwd++eUXNGrUSN1n3bp10CtOZnVOMnbOx9MNO85cx9IDpj9CiIiIXKZl7qGHHsLZs2cxYsQI/PHHHyhRokSu90tOTsaCBQvw2muvITw8HD/99FNxy0tkFeUCffBUx2r4dFUk3lt6BHfXKQNvD3fWLhERuUbL3GOPPYYNGzbg0UcfzTPICR8fHwwaNEi1zP3f//1fccpJZHVPdaqKMgHeOBuTiB82nWENExGR64S57t272+UxjoBLBjsvPy8PNRlCfPZPJK4npGpdJCIiItuHuaSkJFy4cOG22w8ePFj4VyfS2ICmYahTviRuJqfjvyuOal0cIiIi24a533//HTVq1MA999yDhg0bYuvWrVk/Gzp0aOFfnUhj7m4GvNmnrjr/edtZ7OVSJURE5Mxh7p133sHOnTuxZ88ezJgxA0888QR+/vln9TOj0QhnxL1ZnV/rqsG4v0lFyFf4tfkH1ILCRERETjmbNS0tDWXLllXnzZo1UxMb7rvvPhw/fpyhh3RtQu86WHn4CvZfiMXsrWfwf22qaF0kIiIi67fMyWLA+/bty7oeFBSkFgw+fPhwttuJ9CY0wBsv3poM8eGyo7gal6x1kYiIiKwf5iZNmpTVMmfm5eWFOXPmYO3atYV5KiKH81CrymgYFoi4lHRMXHBI6+IQERFZP8x9+eWX6NevH0aOHIn58+cjLi4u62ft2rWDM+IOEK41GWLyfQ3U5eL9l7B43yWti0RERGTdMCeTHTZv3oxhw4bh0KFD6N+/P+666y7VYiczW511EgS5jvoVAzGqczV1/vpfBxAVn6J1kYiIiKy7zpzM7mzevDleeeUVrFq1Sm3XJXuw/vDDD2jVqlVhn47I4Yy+qwZqlwtATEIq3vjrgNbFISIism6Yu379OmJiYtT5tWvXsGLFClSrVg1ffPEFtm3bVtinI3I4Xh5u+O/ARvBwM2DJ/stYuPei1kUiIiKyTpj79ttv1ZIk0jI3bdo0tSyJtM4NGTJE/YzIqbpbu1RX56/8uR/nYhK1LhIREVHx15mbOnWq2rZLtvQKDw/HqVOnEBoaitjYWHTq1AnDhw8vzNMRObTRd1XH+shr2HX2BkbP2Y25T7VRrXZERESOpFC/mTw8PODr66vWl6tevboKciIwMNBpFw121vdF+fN0d8PUIU0Q6Ouptvn6YNkRVhsREek7zLm7uyM52bSYquW6cvHx8dYvGZEDCCvthw8faKjOv91wCn8fuqJ1kYiIiIoe5v7++294e3tntcaZJSYm4uuvvy7MUxHpRvd65fB4O9P2XmN/24MT1/jHCxER6TTM5exOvXz5ctY2Xy1atIAzYicriQm96qBpeCnEJadjxKwdiE1MY8UQEZFDKNZo7u7du1uvJEQOTCY+fDW0OSoE+uBkVAJG/bwL6RmZWheLiIioeGGOOz6QKwkN8MY3jzaHr6c7NhyPwkvz9iMzk7ueEBGRjsOcS8z0dIG3SAVXr0KgmuEq+7fO23Uek5cc5h81RESkKS6aRVRI3eqWxfsD/p3h+vk/x1mHRESkGYY5oiJ4oFkYXrunjjr/aOUxfLzyGFvoiIhIf2FO1p0jclXDO1TFSz1rq/OpqyIxZekRBjoiItJXmNu9ezecHYfM0Z083bka3uxTV51/ve4kxv66B8lpGaw0IiLSTzfrs88+q/ZoJXJVj7eLwPsDGqhJEfP3XMTD325FVHyK1sUiIiIXUeww991331mnJEQ6NrhFOGY93hIlfTyw88x19PlsA7aejNa6WERE5AKKHebGjx+PMWPGYM2aNThx4gTOnj2b7dA7l1h+hayifY0Q/DmqHaqGlMCl2GQM+WaLmhjBxYWJiMiWDMZirvzr5uaWa/CRp5XrGRn6HD908+ZNtX3ZucvXEFY2ROviOK3MzExcvXpVbQln+V3Ss4SUdLzx10G1Dp2oW74k3rmvPpqGl9a0XM5Y146I9cy6djb8TtvPjRs3ULp0acTGxqJkyZIFfpxHcV+Y4+WIsivh7YGPBjVCx5ohKtQdunQTA6ZtwuDmlfBc15ooF+jDKiMiIqsp9p/nlStXvuNhbUlJSRg2bBgCAgJQtmxZTJ48Oc/77tq1C82bN4ePjw+aNGmCrVu3Fvr12MlKRdWvcUWser6TWpNO2r9/2X4OHT9cjbcWHMSVm8msWCIisgoPazQJfvjhh9i7d68KWjl7bf/55x9Yk4zR27dvH7Zt24YLFy5g4MCBiIiIwJAhQ7LdLz4+Hr1798bw4cMxb948TJ8+Hffeey9OnjypgiCRPYT4e+O/AxthcItK+HDZUWw7HYOZm07jpy1n0L1eWTzSqjJaVw2Gmxv/bCAiIo3C3OOPP66C1aBBg1CqVCnYkoRFmT27dOlS1KlTRx1jx47FtGnTbgtzc+fOha+vLyZNmqTG7kkL3m+//aZul5Y9IntqUSUIvz7VGhuPR6sFhiXULdl/WR3lSvqgR72y6FGvHJpWLg0fTy7GTUREdgxzy5cvx/r169GsWTPY2p49e5CWloa2bdtm3dauXTtMmTIla8KF2ZYtW9C+ffus2+RSHiddrYUJc5zMStYi30GZ8SrHkcs3Vevc/N0XcflmMmZtPqMOLw83NK5UCs0rl0bt8iVRu1wAIkJKwNO9mCMiEmOAywfhe3onDAfjgYQrQNwVIDUeSEsE0pKAzAzAzePW4W669PABPH1uXfre8TLTwwfpbt63DrnuDaOHD4zuvjB6eGcNWjDCovXefGowWN6quqVz1N6dKhaONlg85noMkJbAiSasa6fA77T9yMQHTcJchQoVYC/nz59HSEgIvLy8sr1+cnIyoqOj1c8s79ugQYPbynrw4MFcnzslJUUdlrNZRWZaivoik21I3UoQd7U6rlnGH2/3rYdXe9XGxhPRWH7wClYfvYqo+FRsOxWjDsusEurvjbIlfVA+0AfBJbzUJAt/b3d16e3hphYsdjMYVHetG4zwjz+NkGtbUCZqK0Ju7EWJ1Cg1QDbQhu9Jnl/+z/z3/07XVU7rArgQ1jXr2Zn4pBjtF+Ys148bOXIkHnvsMbz33nuoVatWtqAlwsPDYS0S2ry95S/8f5mvSxdsQe6b835m0ro3ceLE226PO74FV93bWaH0lBsJcfKXiAQ6V10uo34QUL9DWYxrXwbnbqRg94V4HL6cgBPRSepITM3E1bgUdey/kPdfbVUMl9DffSP6uW1EhNuV235+LjMUx4xhuGQMwhVjaVxFacQZfZEEbyTDC5lGN7gZMuGBDLjDdOmDVHgjDT6G1FvnqbfO09R1ddz6mfq5weJ2i5/Jc+Qlv3a17G12ORVrZSUiIoeSarRjmKtSpUpW96V5wkOfPn1uu5+115mTWakS0iyZW9NkfFxB7pvzfmYTJkzAuHHjsrXMVapUCeVSTiKwzH1Wew90e5iT70loaKjLhjlLZcsCzWv9e13+/4pOSMXl2GTVHXvxRjJuJKYiITUD8SnpSEhKRe24Deh+fS6qpxzIelwaPHHcpx6O+TXFiRKNcdW3GtI9/JGemoIAfz/VbRvobkCou5s693Q33Lq8/Vy6fs3XPdzM1/+9j5f5MR63fi73cze1Frrqd/ratWv8TrOunQa/0/bjduMG8H6wfcJcbmvLyexRGc+Wnp6uWudkwV1rq1ixoupOldfx9PRUt8mMVgluwcHBt9330qVL2W6T+8rtuZFWu5wtecLtwg6GDBuTMCdBjmEud2VK+qqjoeWNMr5t7xxgwydA9PFbFekOVLsLaDgYnrV7o45XCdSxfAgXDbYbfqdZ186G32n7KOrvwSKFOcv14yRYvfjii/jiiy+yWuGkMEOHDlXLgVhT48aN4eHhgU2bNqFTp07qNpl80bJly9u23WrdujXeeeedrIkRcrlx40a88cYbhXvR8ztMo7EdbJA1ubBT64BlE4Art1rifAKB5k8ArZ4CAjiCiIjI1RS7X0uC3MKFC9Uha87FxMRgwYIFKji9++67sCY/Pz88+uijeP7559VEhtWrV2Pq1KkYNWqU+nlUVFRWt+sDDzyAuLg4vP7662qMn1wmJiaqdekKw5B4DYg5adX3QVQk8VeBX4cCs/qYgpyEuG5vA2MPAl3fZJAjInJRxQ5zs2fPxowZM9CjRw+1GK90r/bq1Qvff/89Zs6cCWv76KOPULduXbRq1UqtLffKK6+oNe6EjLuaM2eOOpc9zRYtWoTFixejRo0aam06Off39y/8i57dYu23QVQ4h/4CvmwNHF5g6k5tMQJ4ZjfQ7lnAm4tgExG5smIvTSLdl7KtVk5BQUGqJczaJIz98MMP6sitLJYk8O3evbv4L3puC9Dk4eI/D1FhpSYCi8eZxseJsg2A+6YB5bIvu0NERK6r2C1zd911F1566SXVpWkm3a1yW9euXeEU2DJHWrh+GviuuynIGdyADi8AI/5hkCMiIuu2zMmYtQ4dOqhZojVr1lStY8eOHUO5cuWwbt06OIWoY0BCFFDi30WJiWzq5Bpg7mNA0nWgRCgwcBZQhesdEhGRDcKcdLHKZIQlS5bg6NGj6jZZPFjGzeVcQFiPjME1gfhI4NxWoPY9WheHXMH+34E/nwIy04EKTYHBPwKBYVqXioiInDXMCVnzrV+/fnBKYc2BI5HA2c0Mc2R7W78Clr5oOq8/AOj3pWlvVCIiIluGObP9+/dj5cqVaiaphDuZXap3xooS5uZw3BzZ3toPgNW3lvNp+STQ831ZtJE1T0RE1g9zstXVc889h/nz56NUqVJqcV5pnZOlQmQ3Blk0WNafW758OVq0aAFdq3Sr/Bf3AGlJgGfu24ERFcv6j/8Ncl1eBTqO50LVRERUIEX6s3/kyJHYt28fvv76a3z22Wf47bff1MK9slCvzGSV48EHH1SBTvcCwwH/ckBmGnBhl9alIWe0+Qtg1UTTede3gE4vMsgREZFtW+ZkAd5Vq1ap7bVEx44d1WLB0lon222JMWPGqC21dE+28QpvDRyabxo3xxmFZE07ZwHLXzGdd54AtB/L+iUiItu3zF2/fh1hYf/OrpOdH2QxX1mOxEzGy1muPadr4W1Ml1xvjqwpciWw6FZ4a/cc0Okl1i8RERVakUdXy7g4S+7u7rfd5jSkZU6c2wZkZmhdGnIGl/aZ1pEzZgCNhpi6V6UVmIiIyF6zWTdt2qRmrZplZGRg69atOH36tLoeGxsLp1G2PuDlD6TEAlcPA+Xqa10i0rObF4GfBwGp8UBER6DPVAY5IiKyf5jr27fvbbcNHjw423WDs7Q0uHsAYS2Ak6tN4+YY5qio0lOAX4cCcZeA0NrAoB8BD/0vrk1ERNopUr9oamoqMjMz8z2ktc5MrjtFVyvHzVFxyILAF3YAPoHAkDmAbynWJxER2T/MtWrVCrNmzcoW1vKSkpKC77//Xj3GOcbNbdW6JKRXO2eaDhiAAd8DQVW1LhEREblqN6vsw/rSSy+pdeT69OmD9u3bo06dOggKClItcNHR0Wo3iA0bNqj7yn0WL14MXZOdIAzuQOw54MY5oFQlrUtEepvwsGS86fyu14AaXbUuERERuXKYK1u2LGbOnIlTp07hm2++wdSpU1V4M7fUyVpzzZo1Q48ePbB7926Eh4dD97z9gfINgYu7Ta1zDHNUUKkJwO/DgIxUoGZPoMPzrDsiInKMvVkjIiIwefJkdRiNRtUiJ4KDg51n8kPO9eYkzMkkiAYPaF0a0oslLwLRkUBABaDfl5y5SkREVmW1heEkvIWEhKjDKYOc4CQIKqz9vwN7fgIMbsCAb4ASwaxDIiKyKidd5ddGKt2aBHHlIJB0Q+vSkKO7eQlYPM503nE8UKW91iUiIiInxDBXGAFlgdIRAIzA+e02+1DICRiNwMJngeRYoHxjoOOLWpeIiIicFMNcYXGfViqIPT8DkcsBdy/gvummhaeJiIgcNcylpaXh3LlzOHr0KGJiYuDUOG6O8hN7AVj2sum8yytAmTqsMyIicrwwFxcXh2nTpqFTp05qj9YqVaqoteZCQ0NRuXJljBgxAtu3b3feljlZxT89VevSkEN2r44BUm6a1iZs84zWJSIiIidXpDD38ccfq/A2Y8YMdO3aFfPnz8eePXtw7NgxbN68GW+++SbS09PRvXt39OzZE5GRkXAaITUA3yAgPRm4tFfr0pCj2T8XOP434O4N9J/G7lUiIrK5Ig3kkRa3devWoV69ern+vGXLlhg2bBimT5+uAt/69etRo0YNOAVZdkVa544uNq03V6mF1iUiR5F0HVj+ium804tAaE2tS0RERC6gSC1zc+bMUUFOdnyQVjnpcs2Nt7c3/vOf/6hg51Q4bo5ys+ptIOEaEFILaDuGdURERI4/AcLd3R1DhgzBtWvX4FKywtxm0xgponPbgR0zTPVw7yeAhxfrhIiI9DGbtUWLFmqPVpdSvhHg4QMkxQBRTjQekIomIx1YNNa0/mDjh4Eq7ViTRESknzD3zDPP4JVXXlFLk7gMD2+gYjPT+bktWpeGtLbta+DKfsC3NNDtba1LQ0RELqbYYW7w4MFqQoSMoXvkkUfw7bffYufOnUhNdfJlOzhujkRCFLDmPVNddH0LKBHCeiEiIrsq9rL00sW6d+9etTSJXE6ZMgWnT5+Gh4cHatWqhX379sG5d4LYrHVJSEv/vAOkxALlGgJNhvKzICIi/YU5WSBYjr59+2bdJrNbJdw5bZATYbIkiQGIOQnEXTHt20qu5fJ+YNcs03mv9wE3d61LRERELsgme7MGBASgQ4cOGDVqFJyWbymg7K119jhuzvXILOZlEwBjJlDvPqByW61LRERELqpIYe7s2bOFuv+FCxfg1OPmzrCr1eUcXgicXm+a1cxJD0REpLcwJ8uRPPXUU3fcezU2NhbffPMN6tevj3nz5sG5x81t0rokZE+yJ++K10znbZ8BSoWz/omISF9j5g4dOoR3330X3bp1g4+PD5o1a4YKFSqo8+vXr6ufHzx4EE2bNsUHH3yA3r17wylVbvfv2KnkWMAnUOsSkT3snAHcOAP4lwXaPcc6JyIi/bXMBQcH4+OPP8alS5fw+eefq31Xo6KiEBlpWkD34YcfVsuTbN682XmDnChZHgiqaho3dZbrzbmE5JvA2vdN551fBrz9tS4RERG5uGLNZvX19cUDDzygDpclA99lRuuZjUDNHlqXhmxt8+dAYjQQXJ1LkRARkfPNZh07dixmzJihWuVSUlLgEiq3N12e3qh1ScjWZAmaTZ+bzu9+A3D3ZJ0TEZH+15mz1KVLF7W23NKlS9WYOYPBoHaGaNiwoTr69OkDp2Peh/PibiAlnt1uzmzdB0Bagmkrtzr/rqtIRETkNGFOFg62XDw4OTkZBw4cUAHvn3/+cc4wJzMZAysBseeA89uAandpXSKyhegTwM6ZpvOuEwGDgfVMRET67WZdsmSJ2vUhKCgId999N5YtW6Zuf/vtt3HPPffgvffew7Vr19Ts1ubNm2PYsGH45JNP4LTMs1rZ1erc23ZlpgPVuwERHbQuDRERUfHC3AsvvID7778fv/32G5o0aYL+/ftj4MCBahmS8PBwLFiwAI0bN8axY8fgEsxdrWe43pxTkqVnDv5h2r6t65tal4aIiKj43axnzpzBs88+iypVqqBr166oXbu2WkRYliuR28Vzzz2HV199FXPnzoXTM7fMXdgBpCUBnr5al4isac17pkvZtqtcA9YtERHpv2VOQty2bduyrsu6ckajEe3a3Qo1AEaOHIkNGzbAJchac/7lgIxU4PwOrUtD1nRpL3BkkalVrtNLrFsiInKOMDd+/Hg88cQTaoycbOnl7u6ugpu00JklJiYiISHBmmV1XDIYnl2tzt0q1+ABoMy/328iIiJdd7M+9thjCAgIUJMaJNBJmJMgJ9t3yVGnTh11e5s2t/YudZXFgw/MA85IayRbcJyCLDdzdAlgcGOrHBEROd/SJAMGDFBHfHw89u7diz179qjjhx9+UGvMybIksl+r3Me8ztx9990Hp188+Nx200bsHl5al4is1io3EAipwfokIiLnXGfO399fjZWzHC+XkZGBI0eOZAU86YL98ssvnTvMhdYC/IJNWz1Ji054K61LRMVxfidwbJmpVa7ji6xLIiJyjUWDzaTbVXZ+kEMmR7jMuDnpaj280NTVyjCnb2ummC4bPgiEVNe6NERERPbZm9XlcZ9W5yBd5cdXAgZ3oNN4rUtDRER0Rwxz1mSe0XpuK5CRbtWnJjtae2usXKMhpmVniIiIHBjDnDWVqQv4BAKp8cDlvVZ9arKTi3uA43/fGiv3PKudiIgcHsOcVWvTHQhvazrn1l76tOFj02X9B9gqR0REusAwZ6uu1tMusvuFM7l2FDi0wHTefqzWpSEiIioQhjlrq9Lh35Y5jpvTlw3/A2AEat8LlK2rdWmIiIgKhGHO2mQjdp9SQMpN4NIeqz892cj1M8C+X03n7cexmomISDcY5qxeo+5AlVu7QZxaa/WnJxvZ9BlgzACqdgbCmrGaiYhINxjmbCGik+ny1DqbPD1ZWdwVYNcPpvMOL7B6iYhIVxjmbCGio+ny7BYgPcUmL0FWtOULICMFCGv5b6sqERGRTjDM2Wqf1hJlgPRk4Px2m7wEWUnSdWD7d6bzDs+btmUjIiLSEV2FuZMnT6JLly7w8fFBrVq1sGjRojvePyQkBAaDIdsRHx9v+4JKIDC3zrGr1bFt+8a0yHPZ+kDNHlqXhoiIyHnDXGZmJvr374+IiAgcPXoUzzzzDAYOHIhTp07lev8rV64gOjpaBcBr165lHf7+/vYpMMOc40tNALZMM513GMdWOSIi0iXdhLm1a9fixIkT+Oyzz1C5cmWMHj0abdq0wYwZM3K9/+HDhxEaGqrCn7TQmQ+7MYc56WaV0ECOZ/dPQFIMUDoCqNtf69IQERE5d5jbsmULmjZtihIlSmTd1q5dO2zdujXX+x86dAg1atSAZkpXAQLDgcx04Oxm7cpBuZMFnTd/YTpvO9q0pAwREZEO6SbMnT9/HhUrVsx2W4UKFXDhwoU8W+aka7Z3796oVKkS+vXrp1r27Ibj5hzb4QXAjTOAXzDQ6CGtS0NERFRkHtCJ5ORkeHt7Z7tNriclJeV6/yNHjiAqKgpvvvkmypQpg/fee09Nnjh48CACAgJuu39KSoo6zG7evKkuJRDKUSRV2sNtz08wnloHY1Gfw8lJ3RqNxqLXcVEYjTBsmgqZt2psPhxGDx8pCJydJnXtgljPrGtnw++0/RT132eHDXOTJ09Wh1mrVq1uG/Mm4cvX1zfXxy9ZsgTp6elZP//xxx9VC93ChQvx0EO3t8RMmTIFEydOvO12mTSRmppapPfgFlAXZeTk0l5cOxcJo3dgkZ7H2b+4sbGxKmS4udmnodjz4jYEX9wNo7s3rkb0h/HqVbgCLeraFbGeWdfOht9p+5F/o50qzI0cOTJb6Pr555+xbNmybPeRLtacXa9mnp6e6rBsxatSpQouXryY6/0nTJiAcePGZWuZk/AnkyhKlSpVxHdRBsaQmjBEHUNoYiRQqXcRn8e5/5GQJWOknu0VMAyrfjSdNHkYoZVrw1VoUdeuiPXMunY2/E7bj5eXl3OFOQlQliGqdevWqqUuISEhaxLE+vXr0blz51wfX7NmTYwfPx4jRoxQ1+Vxx44dQ506dXK9v4S9nN24Qn7pFesXn8xqjToGt9PrgTr3Fv15nJgEjGLXc0FdPQJErpBXhaHNaBhcLNTYta5dGOuZde1s+J22j6L+26ybf9EltIWHh6v15c6cOYNp06Zh+/btGDZsmPp5WlqaGiNn7m++9957MWnSJKxbt06tS/fYY4+pZUp69uxp34Kblyg5uda+r0u52/SZ6VKCdXA11hIREemem57S6vz58xEZGal2f/j000/xxx9/qDXnxMaNG1X30dmzZ7PGwA0ePBgPPvggmjRpoiZQyHg5d3c7L0FRpYNqBcK1w0DcZfu+NmUn9b/vV9N522dZO0RE5BQctps1r65T6VrNq+VOBnabSZfphx9+qA5N+QUBFRoDF3cDJ1YDjYdoWx5XtnU6kJkGVGoNVGqhdWmIiIhcq2VO16rdZbo88Y/WJXFdKXHA9u9N5+3GaF0aIiIiq2GYs2eYO7naJdYzc0i7fgRSYoHg6kDNXlqXhoiIyGoY5uwhrCXg5Q8kXAOu7LfLS5KFjDRgy5em8zaydRe/9kRE5Dz4W80ePLxuTYRgV6smDs4HYs8BJUKBRhyzSEREzoVhzl44bk4bMilm06em85ZPAZ4+GhWEiIjINhjm7B3mzm4BUhPs9rIu79Ra4PJ+wNMPaPGEy1cHERE5H4Y5e5EFakuFAxmpwOmNdntZl7dxqqkKmjxiWiaGiIjIyTDM2YvBwK5We7t8ADixCjC4Aa1H2v3liYiI7IFhzp44bk6jrbv6AkERdn5xIiIi+2CYs/c+rdJKFHUUiD1v15d2ObEXgAO/m865SDARETkxhjl78i0NVGxuOuduEHbYuisdqNweqNjMxi9GRESkHYY5rbpaj6+y+0u7jOSbwM6ZpvO2o7UuDRERkU0xzNlb9bv/3dorI93uL+8Sdv0ApNwEQmoCNXpoXRoiIiKbYpizN+ny8w0CkmOBc1vt/vKusXXXNNM5t+4iIiIXwDBn9xp3B6p3NZ1HLrf7y7vE1l03z5u27mo4WOvSEBER2RzDnBZq3ur6O7ZCk5d37q27bi0SzK27iIjIRTDMaTUJQpYouXYYuHFWkyI4pVPrgMv7AA9fbt1FREQug2FOC7KtVKVWpvNj7Gq1+iLB3LqLiIhcCMOcVmp0N11GsqvVKq4eBo6vlH3TgDbcuouIiFwHw5zW4+akazA1UbNiOI3Nn5su6/QBgqpqXRoiIiK7YZjTSpm6QMkwID0ZOL1es2I4hbjLwL7fTOdtx2hdGiIiIrtimNOKwQDUvNXVynFzxbPtayAjFajUGqjUwhqfDhERkW4wzGnJvDuBjJuTZTWo8FLige3fmc7bPsMaJCIil8Mwp6WIjoCHDxB7DrhyUNOi6Nae2UDyDSCoGlCrl9alISIisjuGOS15+ZnWnBNHFmtaFF2SvW03f2E6bzPKtLsGERGRi2GY01rte0yXRxZqXRL9kTq7cQbwCwYaDdG6NERERJpgmNNazV6m3SAu7weun9a6NPohYww33tq6q8VwUysnERGRC2KY01qJYKByO9M5u1oL7uxm4OIuwN0baDHCVp8OERGRw2OYc6iuVo6bKzBzq1yjBwH/UNt8LkRERDrAMOdIYU5amxKitC6N47tyCDi21LR1FxcJJiIiF8cw5whKhQPlGwHGTODoEq1L4/g2fmq6rNsXCKmudWmIiIg0xTDnKGr3MV0eXqR1SRzbjbPA/rmm83bPaV0aIiIizTHMOVpX68nVQEqc1qVxXJs+B4wZQNXOQMWmWpeGiIhIcwxzjqJMHdMuBrLHKPdqzZ2MJ9z1g+m8/Vh7fjpEREQOi2HOURgMQL37TOcH/tC6NI5p61dAehJQoQkQ0Unr0hARETkEhjlHUv9+0+XxlUByrNalcSzS9bzt639b5ST8EhEREcOcQylTFwipZepq5Zpz2e2cBSTfAIKrA7Xv1egDIiIicjxsmXMk0tpUf4DpnF2t/0pPATZ//u8MVjd3TT4eIiIiR8Qw56hdrTKrNTFG69I4hn2/AXGXgIAKQMNBWpeGiIjIoTDMOZqQGkDZBkBmOnB4gdal0V5GOrDhE9N5m5GAh7fWJSIiInIoDHOO3Dp3YJ7WJdHewT+AmBOAbxDQ7HGtS0NERORwGOYckXmJktMbgLgrcFmZGcC6D03nbUcD3v5al4iIiMjhMMw5oqAIoGIz016trtw6d/BPIOoY4FMKaDFC69IQERE5JIY5R9VoiOlyz89wSZmZ/7bKtRkN+JTUukREREQOiWHOUckSJe5ewJX9wOX9cDky+ePaEcA7EGj1pNalISIiclgMc47KLwio2dN0vmcOXK5Vbu0HpvPWTwM+gVqXiIiIyGExzDmyxg+ZLvf9CmSkwWUcXQxcPQh4lwRa/0fr0hARETk0hjlHVr0rUCIUSIwCjv8Nl2A0AmvfN523egrwLa11iYiIiBwaw5wjc/cEGgxyrYkQMlZOxgh6+QOtR2pdGiIiIofHMKeXrtajS4H4a3D63R7+eeffGawybpCIiIjuiGHO0ZWrD1RoCmSmAXt+glOTsYGyrpx0rbYZpXVpiIiIdIFhTg9aPGG63DHDNNPTGaWnAGveM523H8d15YiIiAqIYU4P6t1vWm/txhng5D9wSrtmAbFngYDyQEvu9kBERFRQDHN64OUHNL61I8T27+FsDGmJMKz/yHSl43jA01frIhEREekGw5xeNHvcdHlsKRB7Ac7Eb/8PMCRcBUpXAZoM1bo4REREusIwpxdlagOV2wHGTGDnDDiN+Ksosfsr03mXVwEPL61LREREpCsMc3piHku243sgLQnOwLBmMtzSEmGs0ASo/4DWxSEiItIdhjk9qd0HCAwHEqOBvU6wX+uVg8DuH9Wpsfu7gBu/jkRERIXF35564u5h2nhebP5S38uUyLZdy1+FwZiJ5Ko9gPA2WpeIiIhIlxjm9KbpUNMG9NGRQOQK6JbsNXtyNYzuXohr9YLWpSEiItIthjm98Q4Amj1qOt/0GXQpI021yiktn0SGdB0TERFRkTDM6VGr/wBunsCZDcCZzdCdLdOAqKOAXzCMHdgqR0RE5HJhLjIyEj4+Pvne76effkKVKlXg5+eHPn364MqVK3AKgWFAk4dN52tvbYGlF7Hn/922q9vbgE+g1iUiIiLSNd2FufPnz6Nv375ISUm54/22bt2KESNG4OOPP8a+ffuQmpqK//u//4PT6PA84OYBnFyjr9a5ZS8DaQmmCQ+NHtK6NERERLqnqzA3f/58NGvWrECtcl9++SUGDx6M+++/H9WrV8fXX3+NlStX4uTJk3AKpcKBJo/oq3Xu2Arg8ELA4A7c8xGXIiEiInK1MLd06VJMnDgRn3zySb733bJlCzp06JB1vXLlyqhYsaJqsXPO1rlNcGgpccDi503nbUYCZetpXSIiIiKn4AEd+eor07ZPa9asKVB3rIQ3SxUqVMCFC7nvayrdtpZdtzdv3lSXmZmZ6nBIJcNgaPIIDDtnwrj8NRifWAEYHDOfG1a+AUPsWRgDK8HY8cWsNfKkbo1Go+PWsRNhXbOenQ2/06xnZ5NZxN+FugpzhZGcnAxvb+9st8n1pKTct8GaMmWKavXL6dq1a2q8naNyqzccIXt/g9vFnYjdNBPJNe6Fo/G6sBlBsgUZgOsdJyH1RiKAxKwvbmxsrAp0btwBwqZY1/bBerYf1jXr2dnExsY6V5ibPHmyOiy7WC27TfMj4+ok0FmSljdfX99c7z9hwgSMGzcuW8tcpUqVEBoailKlSsFxlQE6jAVWv4vAHZ+gZMshgGfu71ETKXEwzHldnRqbD0OpJv1u+8fYYDCoemaYsy3WtX2wnu2Hdc16djZeXl7OFeZGjhyJhx76d7ZjuXLlCvV46WK9dOlSttukizVn16tlq13OljwhAcPhQ0bbZ4Cds2CIPQ/D1i+BjuPhMFa8AsSeU3vKGrq9DUMudSlhThf17ARY16xnZ8PvNOvZmRT196DD/vaU1jBZI858FGQGq6XWrVtj/fr1WddPnTqlwpzc7nSkJa7rm6bzdf8FYhxkxu6+ucDun+SfW6D/F6bdK4iIiMiqHDbMFVZaWhqioqKyBg8+9dRTmDNnDubNm6eWIxk1ahR69eqFiIgIOKUGA4GIjkB6MrBorGkjey1FnwAWPWc6l5ZCKRsRERFZndOEuY0bN6pxV2fPnlXX27Vrp9aaGzt2LOrVqwdPT0/MmjULTstgAO79H+DhY1qqZO8c7cqSlgz8/jiQGg9Ubgd0ekm7shARETk5XYa5zp07q9mPud0mXbJmw4YNU+FOZrD+9ddfKuw5teBq/wanpS8B10/bvwzyuSx4Bri0F/ANAgZ8C7g77NBMIiIi3dNlmKM7aDsGCGsJpNwE5g0HMtLsW13rPwL2/2ZazHjgTKBkBfu+PhERkYthmHM20gomrWHegcD57WrJErs59BfwzyTTee8Pgaqd7PfaRERELophzhmVrgz0+Z/pfMMnwP7fbf+aJ/4xtQSKVk8DzYfZ/jWJiIiIYc5p1b/f1OUq5o8Ezm2z3WvJvrBzHgIyUoE6fYHu79jutYiIiCgbtsw5s65vATV7ARkpwE8PABd3W/81jq0AfrwfSE8CqncDBnzHCQ9ERER2xDDnzNzcTePnwtsAKbHAD/2B8zus9/y7fgTmPPhvkBv8I+BRtK1IiIiIqGgY5pydtz/w8FwgrAWQfAOYeQ9w4I/iryO38FlgwWjAmAE0HAwMmeNYe8ISERG5CIY5VyDbaA39E6jZ07RDhCzou/h5IDWh8M91ZjPwVUdg50zTNl1dXgX6TwfcPW1RciIiIsoHw5wrBboHfwbajDZd3/4t8HlL096p6an5P/7yAeC3R4EZPYGoo0CJUODh34FOL8rOwDYvPhEREeWOS/O72hi6Hu8C1buadmmIPQf8NQr4+y2g3v1ARAcguAbgWxpISwRunAEu7AQOL8w+eaLJUKD7JNP9iIiISFMMc66oWhdg9A5g21fA5i+A+CumcznyYnAD6vYDOrwAlKtvz9ISERHRHTDMuSpPH6Dds0DrUUDkCuD436a16KS1TiZKePoB/mWB8o2AKu2Buv0Bfyff25aIiEiHGOZcnWz/Vbu36TAzGgGDQctSERERUQFx5DrdjkGOiIhINxjmiIiIiHSMYY6IiIhIxxjmiIiIiHSMYY6IiIhIxxjmiIiIiHSMYY6IiIhIxxjmiIiIiHSMYY6IiIhIxxjmiIiIiHSMYY6IiIhIxxjmiIiIiHSMYY6IiIhIxzy0LoCjMhqN6vLmzZtwc2PmtZXMzEzExcXBx8eH9WxjrGv7YD3bD+ua9exsbt68mS2DFBTDXB6io6PVZeXKlYv72RAREREVKoMEBgYW+P4Mc3kICgpSl2fPni1UhVLh/wqpVKkSzp07h5IlS7L6bIh1bR+sZ/thXbOenU1sbCzCw8OzMkhBMczlwdy1KkGOIcP2pI5Zz/bBumY9Oxt+p1nPzqaww7s4GIyIiIhIxxjmiIiIiHSMYS4P3t7eePPNN9Ul2Q7r2X5Y16xnZ8PvNOvZ2XgXMXsYjIWd/0pEREREDoMtc0REREQ6xjBHREREpGMMc0REREQ6xjCXi6SkJAwbNgwBAQEoW7YsJk+ebP9PxoWkpqaibt26WLNmjdZFcUqnTp1Cnz591JqJERERmDJlitoGiawvMjISPXr0UP921KxZE7Nnz2Y129hTTz2Fzp07s55tZOHChTAYDNmOBx54gPVtZTJ9QSY+lCtXDqVLl8aTTz6J5OTkAj+eiwbnYvz48di3bx+2bduGCxcuYODAgeqX4JAhQ6z52RGAlJQUPPLIIzh8+DDrw0ZBWYJckyZNsGvXLhw/fhxDhw5FqVKl8PTTT7POrUgCstR1q1at1L8fR44cwUMPPYSwsDB06tSJdW0DGzZswDfffIOOHTuyfm3k0KFDuOeeezBz5sys27jKg/V9+OGHmD59On799VcV5gYNGoS33noL7733XoEezzCXS6vcd999h6VLl6JOnTrqGDt2LKZNm8YwZ4N/JB5++GFrPy1Z2Lp1qwpw27dvh6+vL6pVq4bnnntOtRgxzFnXlStXUL9+fXz++eeqZU7+AJQWo/nz5zPM2egPFWmVa9++vS2enm6RP7Sl5yQkJIR1YsM/BP/73//ik08+yWplfuONNzBjxowCPwe7WXPYs2cP0tLS0LZt26zb2rVrp34ZchUX61q3bh26dOmi/rom26hdu7bqJpEgZybdJAkJCaxyKytfvjx+//13FeSEtOyvXbuWezvbiAwXkBbnu+66y1YvQbf+6K5RowbrwoYOHjyIqKgo9O3bN+s2aej4+++/C/wcDHM5nD9/Xv0F4uXllXVbhQoVVN91dHS0NT43uuU///kPPv74Y5QoUYJ1YiOhoaHo1q1btm7t77//nr8AbUx++Ul3a9WqVfHMM8/Y+uVcjnRhS5eU/PtBtq/rTZs2oUGDBmocqLQYSYMHWc+JEyfUv9Xr169Xf6BIq/7zzz+vWp8Lit2sOUhoyzkewHxdumCJ9CojI0P9tSd/lLzwwgtaF8epzZ07F1evXsU777yD/fv3c4C+FUkPiXSvyniiMmXKWPOpKQcZMx4XF6c2fZcxc5cuXcLIkSNx8+ZN/O9//2N9WUl8fLw6Xnnllax6HT58uPo3u6D1zDCXg4+Pz20zSKQ1Q1h2VRHpbUyGzNCWLtclS5aoLkGyncaNG6tL+Qf68ccfx8mTJ1X3NhXft99+q1qGZLYf2VbFihVVcDMPHRDSWiR/FH700Udwd3fnR2AFHh4eSExMVMHNPGZOJkTIBCppfZYwne9zWKMgzvbllZYL+cfC09Mz668TCXnBwcFaF4+o0OSvO5nB+ueff2LevHm4++67WYs2mgCxefNm9O/fP9uYxdOnT6vxMNKNQsX3888/Y/fu3WrGn5A/vtPT09UMbZlFHB4ezmq2IssgZ/5OS51fv36dkyKsRJYjEbVq1bqtniWPFOTfDo6Zy+UvaknJMkbATPqxW7Zsyb+sSZeeffZZNaNy0aJFuPfee7UujtOS0DZgwABcvnw567YdO3bAz8+Pfwha0Zw5c9QMS5msJoeMvW3evLk6l/HNZD2rV69WE3ikhdlMgrSMK+fsVuvmDmk8kuWjzOQ7LkG6oPXMlrkc5B/eRx99VA0+nDVrlhr3MnXqVHz11VdW/OiI7Lc0yRdffKGWy2jUqJFqIRLSPWJu2SDraNGihTpGjBihuqBkseaXXnpJhemCdJNQ4VoxzKRFTnpOqlSpwiq0MmnEkDAnXdqTJk3K+k5PmDCBdW1F8h2W4RiybFRQUJAKdlLHsnxUQYdn8F+YXMg/xLKujsxGk4WCZVCiLOBHpDeyVIYYPXq0aqo3HzJjiqxLApt0ZctMePklKP8Qy0zWt99+m1VNuiQrDSxbtkx19UnrkQzXkGAna6+SdX322Weq50QOWW5HFhovzL8dBiMXTyMiIiLSLbbMEREREekYwxwRERGRjjHMEREREekYwxwRERGRjjHMEREREekYwxwRERGRjjHMEREREekYwxwRUT5kNfbvvvtOnctOA7Iqu/koWbIk+vTpg/Pnzxe4HuVxa9asKdB9ZWNzWbBV9n4lIsoNwxwR0R0cPXpU7ewg2/xZrtZ+7do1td3f5s2b1UbvjzzyiE3qUXaUkN0kxo8fz8+JiHLFMEdEdAfvv/8+nnjiCXh4/LuVtb+/v9oAW7ZGq1evntp2Z+3atbhx44ZN6lKC4oIFC3D27Fl+VkR0G4Y5InJJp0+fVt2db731ltqDUvZRzSkuLg6zZ8/G4MGD7/hcvr6+arN3b29vdT0zM1Pt8VyjRg31M9lr8cCBA9kes27dOjRo0EA9rnPnzqo8YubMmWpf6L59+6py/fLLL+qyV69eWV29RESWGOaIyKVt374d+/btw3PPPXfbz/755x+ULVsW4eHheT4+OTlZdbu+8cYbKriJiRMn4sMPP8TUqVPVc1euXBk9evRAfHx81uOmT5+OTz/9FLt27UJUVJQal2e2bds2NU5u9+7d6Natm7qtffv2WLFihZXfPRE5A4Y5InJp48aNQ7Vq1dSR0549e1CnTp3bbh85ciRKlSqFwMBA+Pn54ccff1Rdr8JoNOLzzz/HpEmTVGuatM59/fXXcHd3V/cze/nll1WLXd26dTF8+HDs3bs32zg5CXc1a9ZEcHCwuk1a8Xbs2KGen4jI0r+DQIiIXNCdWt1kgkNQUNBtt8sYuQceeEAFK2lVmzVrFsaMGYOIiAi0aNECMTExqqvUzNPTE82bN8eRI0eybqtdu3bWuQTDpKSkrOsyHs/cymd5m0y0iI6OVudERGYMc0Tk0iRo3YmMf8tJwpQsUSLMAU66QJctW6bORc4WNHkeCWOWrW+WLO+fW5mkZY+IKDfsZiUiykO5cuVUS1hBSACTNefKlCmD0qVLY8uWLVk/kxAnY/Ny67ItKGntk0Bn7nYlIjJjyxwRUR6klU3Gv+UkExmke1Vcv34dc+bMQWRkJIYMGaJmyMqYOpkQERYWpsa9/fe//0VKSkq+s2Lv5NChQ2pShDw/EZElhjkiojx06NBBLU8iQU0mMpjJMibmpUxkAoRMYpCFhWWSgnk2q7TOyTg6CX0yE3Xjxo1qXbqi2rRpk5pQQUSUk8HIqVFERHkaMWKEmiTx+uuva1ZL0qonrXwym1WWOSEissQxc0REd/D888/jp59+yjZ5wd5kSZN77rmHQY6IcsWWOSKifMiab9WrV1fbetlbamoqWrZsieXLl6sFjImIcmKYIyIiItIxdrMSERER6RjDHBEREZGOMcwRERER6RjDHBEREZGOMcwRERER6RjDHBEREZGOMcwRERER6RjDHBEREZGOMcwRERERQb/+Hz/otEtuqph2AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 绘制投影子 β_l(r)\n", "fig, ax = plt.subplots()\n", "for l, beta in sorted(kb.beta_l.items()):\n", " ax.plot(kb.r, beta, label=f'l={l} ({channel_labels[l]})')\n", "ax.set_xlim(0.0, 6.0)\n", "ax.set_xlabel('r (Bohr)')\n", "ax.set_ylabel(r'$\\beta_l(r)$ (Bohr$^{-3/2}$)')\n", "ax.set_title('Al KB projectors (normalized)')\n", "ax.legend()\n", "ax.grid(True, alpha=0.3)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "298eb853", "metadata": {}, "source": [ "## 耦合系数 $D_l$ 与诊断\n", "\n", "$D_l = W_l / Z_l$ 决定了投影子的权重;$|D_l|$ 越大意味着势差 $\\Delta V_l$ 与波函数 $\\phi_l$ 的重叠越弱。可以结合 `diagnostics` 中的 `projector_norms` / `coupling_strengths` 来定位异常通道。\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "17925174", "metadata": { "execution": { "iopub.execute_input": "2025-12-15T08:11:03.147877Z", "iopub.status.busy": "2025-12-15T08:11:03.147774Z", "iopub.status.idle": "2025-12-15T08:11:03.150007Z", "shell.execute_reply": "2025-12-15T08:11:03.149704Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "D_l (Hartree) and diagnostics:\n", "l = 0 (s): D_l = 4.989624 Ha, W_l = 1.590e+00\n", "l = 1 (p): D_l = 2.218599 Ha, W_l = 7.530e-05\n", "\n", "Local potential range:\n", "V_loc^max = 0.000 Ha, V_loc^min = -5.545 Ha\n" ] } ], "source": [ "# 打印耦合系数与诊断信息\n", "print('D_l (Hartree) and diagnostics:')\n", "for l in sorted(kb.D_l.keys()):\n", " D_val = kb.D_l[l]\n", " W_val = kb.diagnostics['projector_norms'][l]\n", " print(f\"l = {l} ({channel_labels[l]}): D_l = {D_val:.6f} Ha, W_l = {W_val:.3e}\")\n", "\n", "print('\\nLocal potential range:')\n", "print(f\"V_loc^max = {kb.diagnostics['loc_potential_max']:.3f} Ha, \"\n", " f\"V_loc^min = {kb.diagnostics['loc_potential_min']:.3f} Ha\")\n" ] }, { "cell_type": "markdown", "id": "5863dda9", "metadata": {}, "source": [ "## 小结\n", "\n", "- KB 可分离形式将半局域势的矩阵复杂度从 $O(N_{\\mathrm{PW}}^2)$ 降至 $O(N_{\\mathrm{PW}})$,是平面波 DFT 代码的标准做法。\n", "- 投影子 $\\beta_l \\propto (V_l - V_{\\mathrm{loc}})\\phi_l$,耦合系数 $D_l = W_l/Z_l$;局域道选择影响投影子的局域性与幽灵态风险。\n", "- `kb_transform` 返回的 `KBResult` 包含 `V_loc`、`beta_l`、`D_l` 及诊断信息,可直接用于验证或导出。\n", "\n", "**下一步**:进入 05-validation 教程,使用 `run_full_validation` 对范数守恒、对数导数进行定量检验。\n" ] } ], "metadata": { "kernelspec": { "display_name": ".venv", "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.12.2" }, "nbformat": 4, "nbformat_minor": 5 }, "nbformat": 4, "nbformat_minor": 5 }