diff --git a/docs/optim.ipynb b/docs/optim.ipynb index bcd4086..eaea392 100644 --- a/docs/optim.ipynb +++ b/docs/optim.ipynb @@ -72,7 +72,7 @@ "import seaborn as sns\n", "from tqdm.notebook import tqdm\n", "\n", - "from mess import Hamiltonian, basisset\n", + "from mess import Hamiltonian, basisset, molecule\n", "from mess.structure import nuclear_energy\n", "\n", "sns.set_theme(style=\"whitegrid\")" @@ -84,9 +84,7 @@ "metadata": {}, "outputs": [], "source": [ - "from mess.interop import from_pyquante\n", - "\n", - "mol = from_pyquante(\"c6h6\")\n", + "mol = molecule(\"water\")\n", "basis = basisset(mol, \"6-31g\")\n", "H = Hamiltonian(basis, xc_method=\"pbe\")\n", "optimiser = optax.adam(learning_rate=0.1)" @@ -162,46 +160,6 @@ " bar.set_description(f\"Total energy: {e:0.06f} (Hartree)\")" ] }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/ubuntu/miniforge3/envs/jax/lib/python3.10/site-packages/pyscf/gto/mole.py:1280: UserWarning: Function mol.dumps drops attribute spin because it is not JSON-serializable\n", - " warnings.warn(msg)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "converged SCF energy = -231.891123588896\n" - ] - }, - { - "data": { - "text/plain": [ - "-231.89112358889585" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "from mess.interop import to_pyscf\n", - "from pyscf import dft\n", - "\n", - "scfmol = to_pyscf(mol, \"6-31g\")\n", - "s = dft.RKS(scfmol, xc=\"pbe\")\n", - "s.kernel()" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -209,31 +167,6 @@ "Next we can look at how the variation in the total energy through the optimisation process." ] }, - { - "cell_type": "code", - "execution_count": 26, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAHACAYAAACGW+2YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABtJUlEQVR4nO3deVhU5fsG8HsGhp0ZZBFEVAQEUVREBRElxaUojVxKrNyycsFUTFPLn2ZmGdnXSisrLbUy0tQsM6zUQMF9XxBZBAUUEHQYNllmfn+QYwQmzMJh4P5cl5fOOe+ceeZpkttz3nmPSKVSqUBEREREGhELXQARERGRIWOYIiIiItICwxQRERGRFhimiIiIiLTAMEVERESkBYYpIiIiIi0wTBERERFpgWGKiIiISAsMU0RERERaYJgiIiIi0oKx0AU0RUqlEuPGjUNZWRmUSiU6duyId955B1ZWVkKXRkRERE2MiPfmq1tRUZE6PL377rswNzfHnDlzhC2KiIiImpxmdZkvIyMDS5YsQVhYGLp06YLhw4fXOS41NRWTJ0+Gr68vgoKCEBUVhfLy8hpj7gUppVKJ0tJSiEQivddPREREhqdZXeZLTk5GbGwsevToAaVSibpOusnlckycOBGurq5Ys2YNcnJysHLlSpSVlWHJkiU1xk6aNAmJiYnw9PTEggULGuttEBERkQFpVpf5lEolxOLqk20LFy7EhQsXsHv37hpjPv/8c6xbtw4HDhyAjY0NAOCHH37AsmXLcODAATg6OtYYX1VVhVWrVsHW1hYvvfRSo7wPIiIiMhzN6szUvSD1X+Li4hAYGKgOUgAQGhqKpUuXIj4+HqNGjaox3sjICCNHjsTcuXM1DlOnT5+GSqWCRCLR6PlERETU+CoqKiASidCzZ8//HNeswlR9pKWlYfTo0TW2SaVSODg4IC0tDQBQUFAAALC1tYVKpcLevXvRqVMnjV9TpVJBpVLVmpdFREREhq/FhanCwkJIpdJa22UyGeRyOQAgPz8fCxYsQEVFBQDAw8MDixcv1vg1JRIJVCoVPDw8ND5GXUpLS5Geng5XV1eYm5vr9NhUjT3WL/ZXv9hf/WJ/9U/oHqekpNTrC2gtLkzVR6dOnbBjxw6dHlMkEsHCwkKnx7zH3Nxcb8emauyxfrG/+sX+6hf7q39C9bi+3+RvVksj1IdUKoVCoai1XS6XQyaTCVARERERGbIWF6bc3NzUc6PuUSgUyMvLg5ubm0BVERERkaFqcWEqODgYCQkJKCwsVG+LiYmBWCxGUFCQgJURERGRIWpWc6ZKS0sRGxsLAMjKykJRURFiYmIAAP7+/rC1tUV4eDi++eYbREREYOrUqcjJyUFUVBTCw8NrrTFFRERE9DDNKkzl5+dj9uzZNbbde7x582YEBARAJpNh06ZNWL58OSIiImBpaYkxY8YgMjJSiJKJiIjIwDWrMOXi4oKkpKSHjnN3d8fGjRv1XxARERE1ey1uzhQRERGRLjFMEREREWmBYYqIiIhICwxTRERERFpgmCIiIiLSAsMUERERkRYYpgxYQWEZDpyTQ1FSLnQpRERELRbDlAHbfzILsRcU+OHPVKFLISIiarEYpgyYZ3sbAEDChZuoqFQKWwwREVELxTBlwHw62sLKXIzi0kqcvpIrdDlEREQtEsOUAROLRfBoYwYASL1+R9hiiIiIWiiGKQNnL5UAADJziwSuhIiIqGVimDJw9tLqe1Vn5jFMERERCYFhysDZWVeHqay8IiiVKoGrISIiankYpgycrbUxjMQi3C2vQr68TOhyiIiIWhyGKQNnJBahdStzAEBWnkLgaoiIiFoehqlmoK2DJQAgi5PQiYiIGh3DVDPQxt4CACehExERCYFhqhlwtq8+M8XlEYiIiBofw1QzcC9MZfHMFBERUaNjmGoGnP++zJd3uxRl5ZUCV0NERNSyMEw1A1JLE9hYmQIAknlbGSIiokbFMNVM+HVuDQA4euGmwJUQERG1LAxTzURfHycAwJELN1DFldCJiIgaDcNUM9HTszUszIyRU1CCX+PThC6HiIioxWCYaibMTI0xaXhXAMDmPYm4catY4IqIiIhaBoapZuTRgA7o7mGPu+VV+DD6FKqqlEKXRERE1OwxTDUjYrEIrzzjC3NTY1y6WoBvfksUuiQiIqJmj2HqAZYsWYIBAwbAy8tL6FIaxMnOErPG+gIAth9IwZELN4QtiIiIqJljmHqAESNGYOfOnUKXoZH+PdpieP+OAIBV353ElWu3Ba6IiIio+WpWYSojIwNLlixBWFgYunTpguHDh9c5LjU1FZMnT4avry+CgoIQFRWF8vLyGmP69OkDe3v7xihbL6Y86QM/r9a4W16FtzYcQWauQuiSiIiImqVmFaaSk5MRGxuLDh06wN3dvc4xcrkcEydOREVFBdasWYPIyEhs3boVK1eubORq9cvYSIwFE3rD3UUGeVE53vgsnvfuIyIi0oNmFaZCQkIQGxuLjz/+GF27dq1zTHR0NIqLi7F27VoMGDAAY8aMwfz58xEdHY2cnJxGrli/LMwkWPZSIDo4WaOg8C5e/zQe2QxUREREOtWswpRY/PC3ExcXh8DAQNjY2Ki3hYaGQqlUIj4+Xo/VCUNmZYq3pwWhvZM1CgrLsPCTQ0i+zjlUREREumIsdAGNLS0tDaNHj66xTSqVwsHBAWlp+ls5XKVSoaSkRKfHLC0trfH7g5gYAYsn+uHtjSdxLacICz85hOkju6JfNyed1tMc1bfHpBn2V7/YX/1if/VP6B6rVCqIRKKHjmtxYaqwsBBSqbTWdplMBrlcrn68cOFCJCQkAACCg4MREBCA999/X+PXraioQGKiftZ9Sk9Pr9e4Z4Ol+PFQBVJu3MVHW88j9kQqQnvZwFTSrE5Q6kV9e0yaYX/1i/3VL/ZX/4TssYmJyUPHtLgwVV+6npAukUjg4eGh02OWlpYiPT0drq6uMDc3r9dzundVYtuBNPwUdxVn0kpw47YKM5/2gWc7G53W1lxo0mOqP/ZXv9hf/WJ/9U/oHqekpNRrXIsLU1KpFApF7WUC5HI5ZDKZ3l5XJBLBwsJCL8c2Nzdv0LFfeLI7Anza4n9bTiLndimWrj+Bpwd3wtghXpAY8yxVXRraY2oY9le/2F/9Yn/1T6ge1+cSH9DMJqDXh5ubW625UQqFAnl5eXBzcxOoqsbX1c0OH786CAP9XKBUqvDDH1cw7+M4ZNwsFLo0IiIig9LiwlRwcDASEhJQWHg/NMTExEAsFiMoKEjAyhqfpbkErz7XCwsm9Ia1hQRpWXJEro7Fzr9SUKVUCV0eERGRQWhWYaq0tBQxMTGIiYlBVlYWioqK1I8LCgoAAOHh4bC0tERERAQOHTqE7du3IyoqCuHh4XB0dBT4HQijf4+2WDs/BL29HVFRqcRXv1zEG5/F42Z+sdClERERNXnNas5Ufn4+Zs+eXWPbvcebN29GQEAAZDIZNm3ahOXLlyMiIgKWlpYYM2YMIiMjhSi5ybCVmmHJlAD8fvQaNvx8HhfT8jHrgwOY8mQ3DAtoX+/rxkRERC1NswpTLi4uSEpKeug4d3d3bNy4Uf8FGRiRSIRH+3ZAj072+DD6NC6m5WPttjM4cuEGXnnGF7ZSM6FLJCIianKa1WU+0g0nO0usmB6EF0Z0hcRYjBOJOZj5/n4cOpsldGlERERNDsMU1clILMLIgR5YHfkI3F1kUJRU4L3NJ7By03Hk3eZqv0RERPcwTNF/6uAkxapZwQgf6gWxWIT4c9mYHrUPP/yZhPKKKqHLIyIiEhzDFD2UsZEYzz3WGR9GPoIuHW1xt7wK3/52GTPfP4BjF29CpeIyCkRE1HIxTFG9dXSWYWVEf7z6XC/YSk1xI78Yy786ikWfxuPS1XyhyyMiIhIEwxQ1iEgkwkA/F3y2YDBGD/KAxFiMi2n5WLD2EJatP4Kr2fKHH4SIiKgZaVZLI1DjsTCTYNLwrhje3w3RfyThj2PXcCIxBycv56Bfd2eMHeKJjs76u9chERFRU8EwRVqxtzHHzKd9MXKgB76LuYyDZ7IQfzYb8WezEdDVCc8M8YRn+1ZCl0lERKQ3DFOkE20drPDa+N54Zogntv55BYfOZuHoxZs4evEm/LxaY8zgTvBxs+NK6kRE1OwwTJFOubaR4rXxvTFumBd+3J+Mv05l4lRSLk4l5cK1jRSP93PFwF7tYG7Kjx4RETUPnIBOetHO0RqR4/ywbsFgPNq3A0wkRki/UYhPt5/DxGV78fmOc5ysTkREzQJPD5BetbG3xMynfTFpeFfsO34Ne+KvIvtWMXbHX8Xu+KvwaGeDYf7tEdzTBZbmEqHLJSIiajCGKWoUVuYShAW7Y0R/N5xJzsPeI+k4dvEmUq7fQcr1O1j/80X07+GMof7t0ZVzq4iIyIAwTFGjEotF8PNqDT+v1pAX3cWBk9fx+9FruJ6jwP4T17H/xHU421tiYK92GOjngjb2lkKXTERE9J8YpkgwMitTPPWIB8KC3ZF07Tb+OHoNB89kIvtWMbbsvYwtey/Dq0MrDPRzwQDftpBZmQpdMhERUS0MUyQ4kUiEzh1s0bmDLV4M88Hh89n462QmzibnISnjNpIybuPLXRfQ09MBA/1c0NenDcz4bUAiImoi+BOJmhRzU2OE9G6PkN7tcbuwDAfPZOHAqUykXL+Dk5dzcfJyLkxNjNDH2xH9e7RFL+/WMDPhx5iIiITDn0LUZLWSmuHJYHc8GeyOrLwi/HUyE7GnMnEjvxiHzmbj0NlsBisiIhIcf/KQQWjrYIXnHuuMZx/1QvL1O4g/m41D57KRW1BSI1j19nbEAAYrIiJqRPxpQwZFJBLBs30reLZvhUnDuyAlszpYHTxbHazu3RfQ1MQIvTs7IqiHM/p4O3KOFRER6Q1/wpDBEolE6NSuFTq1a4WJT9wPVofOZiOnoATx57IRfy4bJhIj9OrcGv17OKNPFyfeyoaIiHSKP1WoWfh3sErNlFeHqbPZuJFfjMPnb+Dw+RswMRbDr3NrBPVoC/8ujrAw46rrRESkHYYpanZEIhE82tnAo50NJjzujbSs+8Eq+1Yxjly4iSMXbkJiLIafV2sE9XCGfxcn3s6GiIg0wjBFzZpIJIK7iw3cXWwwPtQb6TcK1ZcCs/KKcPTiTRy9eBPGRmL09HJAUHdnBPi0gRWDFRER1RPDFLUYIpEIHZ1l6Ogsw3OPdca1mwocOpuN+HNZuJ5ThOOXcnD8Ug6Mjc7A17M1grq3QXc3G6HLJiKiJo5hilokkUiEDm2k6NBGiuce64yMm4VI+Hu5hWs3FTiRmIMTiTkwEovQ0dEEg0ukCPbrAKmlidClExFRE8MwRQSgg5MUHZykGPdoZ1zPUajnWKXfKETKjbtI+ekSvvw5ET087DGwlwsCuznzW4FERASAYYqolnaO1ggf6oXwoV5IycjDL39dQGquEhk3i3D6Sh5OX8nDp9vPIdCnDQb1aocenexhZCQWumwiIhIIw9R/WLJkCQ4cOIDc3FwkJSUJXQ4JwNnBEsE+Ukz19sadYiViT2fhwMnruHGrGH+dysRfpzLRytoUwT1dMKiXC9zayiASiYQum4iIGhHD1H8YMWIEZs2ahaCgIKFLoSbA2cEK44Z5IXyoJ65cu40DJzMRdzoLtxV3sSsuFbviUtHeyRohvdohpE87tLI2E7pkIiJqBAYVpjIyMrBhwwacPXsWycnJcHNzw+7du2uNS01Nxdtvv43Tp0/D0tISYWFhmDNnDkxMGjZ5uE+fProqnZoRkUgErw628OpgiylP+uDU5RwcOJWJYxdv4tpNBTb+egnf/JYI/65OeLRvB/h6toaRmGeriIiaK4MKU8nJyYiNjUWPHj2gVCqhUqlqjZHL5Zg4cSJcXV2xZs0a5OTkYOXKlSgrK8OSJUsEqJqaM4mxGAE+bRDg0wZFpRWIP5uNP45lICnjtnrV9datzDE0oAOG9GkPextzoUsmIiIdM6gwFRISgiFDhgAAFi5ciAsXLtQaEx0djeLiYqxduxY2NjYAgKqqKixbtgxTp06Fo6MjAGDkyJHIzs6u9XwfHx9s2LBBf2+Cmi0rcwke7dsBj/btgPQbhdh7JB0HTmYi93Ypvou5jO/3XkYvb0c8FuiKXp0debaKiKiZMKgwJRY//BtTcXFxCAwMVAcpAAgNDcXSpUsRHx+PUaNGAQB27typrzKJ4NpGiqkju2PS8K5IOJeNvUcycDEtX70wqJOdBZ4IcsMQ//ZcbZ2IyMAZVJiqj7S0NIwePbrGNqlUCgcHB6SlpQlUFaBSqVBSUqLTY5aWltb4nXRPFz0O8LZDgLcdsvKKse9EJv46lY2b+SXY8PMFfBuTiEd82+DRgHZwaW2lq7INBj/D+sX+6hf7q39C91ilUtXrG9rNLkwVFhZCKpXW2i6TySCXyxt0rIULFyIhIQEAEBwcjICAALz//vsa1VVRUYHExESNnvsw6enpejku3aerHvdxBXq4tMa59BIcTSpCnrwSvx/LxO/HMuHmZIq+XlbwcDaDuIUtr8DPsH6xv/rF/uqfkD2uz5fXml2Y0qWVK1fq7FgSiQQeHh46Ox5QndTT09Ph6uoKc3NObNYHffW4Rzfg+eEqXLx6GzFHruHE5Tyk3byLtJt34eJgiRH9XdG/uxOMjZv3YqD8DOsX+6tf7K/+Cd3jlJSUeo1rdmFKKpVCoVDU2i6XyyGTyQSoqJpIJIKFhYVejm1ubq63Y1M1ffXY38cS/j4uyCkowa/xV7H3SDoy84rx2c6L2Lo/FU8OcMdjgR1gYda851XxM6xf7K9+sb/6J1SP67sIc7P7Z6+bm1utuVEKhQJ5eXlwc3MTqCqi/+Zoa4EXRnTFV4uHYfLwLrCVmiJfXoavd1/E5OW/Y+Pui8iXc14GEVFT1OzCVHBwMBISElBYWKjeFhMTA7FYzJXMqcmzNJdg1KBOWP/GUMwe64t2jlYoKavE9gMpeHHFn/hs+1nk3WaoIiJqSgzqMl9paSliY2MBAFlZWSgqKkJMTAwAwN/fH7a2tggPD8c333yDiIgITJ06FTk5OYiKikJ4eLh6jSmipk5ibIQh/h0Q0rs9TiTm4Mf9yUhML8CehHT8fjQDQ/07YExIJ7S25aUFIiKhGVSYys/Px+zZs2tsu/d48+bNCAgIgEwmw6ZNm7B8+XJERETA0tISY8aMQWRkpBAlE2lFLBbBv6sT+nRxxIXUfHz/exLOp97Cb4erQ9XgPu3x9OBOcLKzFLpUIqIWy6DClIuLC5KSkh46zt3dHRs3btR/QUSNRCQSoZuHPbp52ONC6i1E/5GEs8m38PvRDPx5/BpCerXD00M6wdm+5a1VRUQkNI3C1OXLl3Hy5Emkpqbi9u3bEIlEaNWqFdzc3ODn5wdvb29d10lEf/Nxt8fb7va4dDUfP/xxBaeScvHn8WvYf/I6hvq3R/hQL94DkIioEdU7TOXn52PLli346aefkJ2dDZVKBYlEAplMBpVKhcLCQlRUVEAkEqFNmzYYOXIkxo0bB3t7e33WT9Rideloh2UvByIpowDf/56Ek5dzsfdIBvafuI4ngjpiTEgnyKxMhS6TiKjZq1eYev/997FlyxZYWlriscceQ79+/dC1a9daE7pzcnJw8eJFxMfHY+vWrfjqq6/w/PPP49VXX9VL8UQEeHWwxZsvBeJiWj6++S0RF9Py8VNsKvYeycDIR9wR9oh7s1+niohISPUKUydOnMD777+PwYMH/+cCVo6OjnB0dERISAgWL16Mffv2Yf369TorlogerKubHd6dEYRTSbnYvCcRaVlybPk9Cbvjr+LpwZ0Q2q8jTCVGQpdJRNTs1CtM/fDDDw0+sEgkwpAhQzBkyJAGP5eINCMSidCrsyN6erZGwvlsfPvbZWTlFWHDzxexKzYV4x/3xkC/dhCLW9a9/4iI9KnZLdpJRNVLKvTv0RafzB+EWc/4wt7GHLfkZVj9/WlEfhiLs8l5QpdIRNRsaLw0QlVVFWJiYnD06FHk5+dj1qxZ8PLygkKhwOHDh+Hn58fJ50QCMzISY2hABzzi54JfDqZh674rSMuSY/G6BPT2dsTk4V3Q3kkqdJlERAZNozBVWFiIF198EefOnYOFhQVKS0vx/PPPAwAsLCzw9ttv46mnnsLcuXN1WiwRacZEYoTRIZ0wxL89ov9Iwm8J6TiRmINTl3MwrK8rnh3mhVZSM6HLJCIySBpd5lu1ahWSk5OxYcMG/Pnnn1CpVOp9RkZGePTRR9W3fSGipkNmZYqpI7vjk9dCENitDZQqIOZwOl5+909E/5GEsruVQpdIRGRwNApT+/btw/jx4xEUFFTnt/tcXV2RlZWldXFEpB9tHazw+iR/rIzoD8/2Nigrr8J3MZcx/b19iD2VWeMfSERE9N80ClMKhQIuLi4P3F9ZWYmqqiqNiyKixtHVzQ6rZgXjted7o7WtBW7Jy7Dqu5NY+MkhpGbeEbo8IiKDoFGYat++PS5evPjA/fHx8XB3d9e4KCJqPCKRCAN6tsVnr4Xg+dDOMDUxwqWrBYj8MBZrt52BvOiu0CUSETVpGoWpMWPGYPv27dizZ4/6coBIJEJ5eTlWr16NgwcPYuzYsTotlIj0y0RihLFDvLBuwWAE92wLlQrYeyQDU9/9E7viUlFZpRS6RCKiJkmjb/NNnDgRKSkpmDt3LqTS6q9Vz5s3D3fu3EFlZSXGjh2Lp59+WqeFElHjsLcxx/zne+Pxfh3xxU/nkZYlx/pdF7D3SDpeDOsGP6/WQpdIRNSkaBSmRCKRevmDmJgYXLt2DUqlEu3bt0doaCj69Omj6zqJqJF1dbPD/+Y8gj+PZWDznkRczynC0i8OI6CrE14M84GTnaXQJRIRNQkaL9oJAL1790bv3r11VQsRNTFGYhEe7euKoB5t8f3vl/Hroas4evEmTifl4pkhnhg1yAMSY97vj4haNq1uJ5OTk4Pdu3dj06ZNuHnzJoDqldHv3LnDb/MRNSNW5hK8FNYNa+YNQncPe5RXKvFtzGW8suoAzlzJFbo8IiJBaXRmSqVSYeXKlfjuu+9QWVkJkUgET09PODk5oaSkBCEhIZg1axYmTZqk43KJSEjtHK3x9rR+iDudhQ0/X0BWXjH+7/PDGODbFlOe7Ao7mbnQJRIRNTqNzkytX78emzdvxgsvvICvv/66xgJ/1tbWGDZsGH7//XedFUlETYdIJMIjfi74bMFgjBjgBrEIOHgmC9Pf24+f41JRxW/9EVELo1GY2rZtm/ree507d66138vLC+np6drWRkRNmKW5BC8/1Q3/m/MIvNq3QundSny56wIiP4zF5fQCocsjImo0GoWpGzduoGfPng/cb25ujqKiIo2LIiLD4e5ig6hXBmDm0z1gZS7B1exCzF9zEGu2noGipFzo8oiI9E6jMGVnZ4cbN248cP/FixfRpk0bjYsiIsMi/vtbf+sWDsZQ//YAgN+PZmDGe/sRd5r3+iOi5k2jMDV06FBER0fj+vXr6m33bnh86NAh7Ny5E4899phuKiQigyGzMsWssT3x3sz+aOdojTtFd/H+tyfx1oajyC0oEbo8IiK90ChMzZo1Cw4ODggLC8OCBQsgEonw5ZdfYty4cXjppZfg6emJadOm6bpWIjIQXTra4aO5j+DZRzvD2EiME4k5iHh/P3bFpaJKybNURNS8aBSmrK2tsXXrVrz44ovIycmBqakpjh8/DoVCgYiICGzZsgXm5vyKNFFLJjE2wrhhXvj41YHo6maHsvIqrN91AfM/jsPVbLnQ5RER6UyD15m6e/cufvjhB3h7e2PGjBmYMWOGPuoiomainaM13pkehD+OZeDrXy4i+fodzFkdi1EDPRDWv53Q5RERaa3BZ6ZMTU2xatUqXL16VR/1EFEzdG+C+qcLBiOouzOUShV+3J+MeWsPI+1mmdDlERFpRaPLfJ06dUJWVpauayGiZs5WaoaFE/tg8WR/2MnMkFNQis37b+GLXZdQUlYhdHlERBrRKExFRkYiOjoaCQkJuq6HiFqAAJ82+PS1EAzzdwEA7DuRhYj3D+DUZd7nj4gMj0b35vv2229hY2ODKVOmwMXFBS4uLjA1Na0xRiQS4bPPPtNJkY1NqVRi3LhxKCsrg1KpRMeOHfHOO+/AyspK6NKImg0LMwmmjPBGG+syxJwqRs7tUiz98jCG9GmPKWE+sDKXCF0iEVG9aBSmrly5AgBo06YNqqqqkJGRodOihCYWi7FhwwZ1eHr33Xexfv16zJkzR9jCiJqhjo5miJrZHdv/Sscvh9Lw5/FrOJWUi4ine8C/i5PQ5RERPZRGYWr//v26ruOhMjIysGHDBpw9exbJyclwc3PD7t27a41LTU3F22+/jdOnT8PS0hJhYWGYM2cOTExMGvR694KUUqlEaWkpLCwsdPI+iKg2MxMjvPRUNwT1cMZH0aeRfasYyzccxcBeLnj5qW6wtmjY/79ERI1JozlTx48fR0HBg29kWlBQgOPHj2tcVF2Sk5MRGxuLDh06wN3dvc4xcrkcEydOREVFBdasWYPIyEhs3boVK1eu1Og1J02ahMDAQFy9ehUvvviiNuUTUT106WiHj+cNwsiBHhCLgL9OZmJG1H4cPv/g21cREQlNozA1YcIExMfHP3D/kSNHMGHCBI2LqktISAhiY2Px8ccfo2vXrnWOiY6ORnFxMdauXYsBAwZgzJgxmD9/PqKjo5GTk6MeN3LkSAQEBNT6NWXKlBrH27hxIxISEuDj44MtW7bo9P0QUd1MJUZ4YURXvPfKALi0tsIdxV28s/EYVn17EkW8cTIRNUEaXeZ72E1Ly8vLYWRkpFFBDyIWPzz3xcXFITAwEDY2NuptoaGhWLp0KeLj4zFq1CgAwM6dO+v9ukZGRhg5ciTmzp2Ll156qcF136NSqVBSott7k5WWltb4nXSPPdav/+pvewczvDvNHz8eSMPPh9IRezoT51LyMH1UV/TwsGvsUg0SP7/6xf7qn9A9VqlU6nsP/5d6h6ns7Owaa0ulpaXVeSmvsLAQ0dHRcHZ2ru+hdSYtLQ2jR4+usU0qlcLBwQFpaWn1Ps69S5i2trZQqVTYu3cvOnXqpFVtFRUVSExM1OoYD5Kenq6X49J97LF+/Vd/e7YD7Ic6YOfh2yhQ3MU7m06hdydLDOspg4mxRifXWxx+fvWL/dU/IXtcnznX9Q5TO3bswNq1ayESiSASibBu3TqsW7eu1jiVSgUjIyMsW7asYdXqQGFhIaRSaa3tMpkMcnn97wWWn5+PBQsWoKKiehFBDw8PLF68WKvaJBIJPDw8tDrGv5WWliI9PR2urq68F6KesMf6Vd/+egN4JKAK3/2ejL1Hr+NEcjEy85WYMdoHXu1tGq1eQ8PPr36xv/ondI9TUlLqNa7eYSo0NBSdOnWCSqXCnDlzMH78ePTu3bvGGJFIBHNzc3h7e8Pe3r5hFTchnTp1wo4dO3R6TJFIpLdvBJqbm/PbhnrGHutXffprYQHMfMYP/X1d8FH0adwsKMWb649jdEgnjBvWGRKepXogfn71i/3VP6F6XJ9LfEADwpS7uzvc3d2hUqmwaNEi9O/f/4HfqhOKVCqFQqGotV0ul0MmkwlQERHpmq9na6yZH4Ivdp7DgZOZ2LYvGScScxA5zg8dnfn/ORE1vgb/U66iogLvvfeeIGtNPYybm1utuVEKhQJ5eXlwc3MTqCoi0jUrcwnmPtsLiyb2gdTSBFezCzH3wzj8uD8ZVcr//oIMEZGuNThMmZiYwN7evsGLYDaG4OBgJCQkoLCwUL0tJiYGYrEYQUFBAlZGRPrQr7sz1s4fBP8uTqisUmLTr5ew6JNDuJlfLHRpRNSCaDTJYOTIkdi1axfKyxtvzZfS0lLExMQgJiYGWVlZKCoqUj++9+278PBwWFpaIiIiAocOHcL27dsRFRWF8PBwODo6NlqtRNR4WlmbYfEL/pg91hfmpsZITC/ArA/+woGT14UujYhaCI3WmfLy8sK+ffswfPhwjBw5Em3btoWZmVmtccOGDdO6wHvy8/Mxe/bsGtvuPd68eTMCAgIgk8mwadMmLF++HBEREbC0tMSYMWMQGRmpszqIqOkRiUQY4t8B3Twc8MF3J5GYXoD/bTmFk4m5mD66Oyx502Qi0iONwtTcuXPVf/7oo4/qHCMSiXS6rpKLiwuSkpIeOs7d3R0bN27U2esSkeFwtLXAuzOCsG1/Mr7/PQmxpzORmJ6PV5/rhS4dudAnEemHRmFq8+bNuq6DiEgnjIzECB/qBd9ODlj13UnkFJRg0SeH8MwQL4QP9YSREZdQICLd0ihM+fv767oOIiKd6uxqi49fHYh1O6qXUIj+IwlnruTi1ed6wcnOUujyiKgZ4T/RiKjZsjCrXkJh3nO9YGFmjMsZtzk5nYh0TqMzUwCQl5eHH3/8EZcuXYJCoYBSqayxXyQSYdOmTVoXSESkrUf8XNDZ1Rb/23ISl67en5w+Y0x3WJhxcjoRaUejMHX58mVMmDABZWVl6NixI65cuQIPDw8UFhYiJycH7du3h5OTk65rJSLSmKOtBd6ZXnNy+pXrt/Ha+N7wcLERujwiMmAaXeb74IMPYGFhgZiYGHz99ddQqVR4/fXXERsbi9WrV0Mul2PevHm6rpWISCv3JqevnNEfDq3MceNWMeZ/fBC7D6VBpeLK6USkGY3C1KlTpzB27Fg4OztDLK4+xL2/iEJDQzFixAhERUXprkoiIh3y7miLj+YOREDX6pXTP995Hu9uOo6iksZbiJiImg+NwpRSqYS9vT2A6psLGxkZ4c6dO+r9Xl5euHjxok4KJCLSB2sLE7wx2R8vPeUDYyMRDp+/gdn/+wuXMwqELo2IDIxGYcrFxQWZmZnVBxCL4eLigsOHD6v3nzp1CtbW1rqpkIhIT0QiEZ4c4I6oVwagjZ0lcm+XYuHaQ9hxIBlK3jCZiOpJozDVv39/xMTEqB+PGzcO27Ztw6RJkzBx4kT89NNPGD58uM6KJCLSp07tWuHDuY9ggG9bVClV+Hr3Jby14QjkRXeFLo2IDIBGYWratGn44IMPUFFRAQCYOHEiZs2ahTt37kChUGDGjBmYM2eOLuskItIrCzMJ5j/fCzOf7gETYzFOXs7FrA/+woXUW0KXRkRNnEZLI8hkMshkMvVjkUiEGTNmYMaMGTorjIiosYlEIjza1xVeHWzx3ubjyMwtwhvrEjDx8S4YOdAdIpFI6BKJqAniCuhERP/i2kaK1XMewcBeLlAqVfh690W8u+k4iksrhC6NiJqgep+Z+vrrrxt0YJFIhEmTJjW0HiKiJsHM1Bhzx/mhi6stvvjpAg6fv4H0G4VYNLEPOjrLHn4AImox6h2m3nvvvVrbRCLRAxe6Y5giIkMnEokQ2q8j3F1ssHLzcdy4VYx5Hx/EjNHdMbhPe6HLI6Imot5hat++fTUey+VyjBo1CqtWrULPnj11XhgRUVPh2b4VPowciA+2nMSpy7n4MPo0EtML8PJT3WAiMRK6PCISWL3DVNu2bWs8trCwAADY2dnV2kdE1NxILU2wdEpfbN13BVv2XsbeIxlIybyDhRP6wMnOUujyiEhAnIBORFRPYrEI4UO98OZLgbC2MEFqphyRq2Nx/NJNoUsjIgExTBERNZCfV2t8OPcReLVvhaLSCry14Si27L3MVdOJWiiGKSIiDbRuZYF3I/rjiaCOAIDvf0/COxuPoaSMyycQtTRahykuYkdELZXEWIxpo7pjTnhPSIzFOHrxJl79KA6ZuQqhSyOiRlTvCegjRoyo8VipVAIAFi9eDHNz81rjRSIRfv75Zy3LIyJq+gb3aY92jtZ4d+MxZOYW4dWP4vDqs73g39VJ6NKIqBHUO0zZ2NjU2mZra6vLWoiIDJZn+1b4X+QjeG/zCVxMy8fyr47iucc645nBnhCLeQafqDmrd5j65ptv9FkHEZHBa2VthuVT+2HDzxfwa/xVfBdzGamZdxA5zg8WZhKhyyMiPeEEdCIiHbo3j2rWM74wNhLjyIWbmPfxQWTnFQldGhHpSb3C1I0bNzR+AW2eS0RkqIYGdMDKiCDYSs1wPUeBuR/G4lRSrtBlEZEe1CtMDR06FIsWLcK5c+fqfeBTp07htddew7BhwzQujojIkHl1sMWHkY/A29UWxWWVWLb+CHYfSnvgPU2JyDDVa87Uli1b8OGHH+KZZ56Bs7Mz+vbti65du8LFxQVSqRQqlQqFhYXIzMzEhQsXcOTIEeTk5CAgIADfffedvt8DEVGT1UpqhhXT++GTH89i3/Hr+HzneVy7qcDLI7vB2IgzLYiag3qFqe7du+Orr75CYmIitm/fjv3792PHjh0A7q8zde9fWm3atMGQIUMwevRoeHt766ls/QsJCYGZmRkkkupJox988AE8PDwEroqIDJHE2Aizx/ZEe0drbPz1En47nI6svCIsnNgH1hYmQpdHRFqq97f5AMDb2xuLFy/G4sWLkZOTg7S0NNy5cwdA9dIJbm5ucHR01Eedgvjiiy/g4uIidBlE1AyIRCKMGtQJLq2tseq7EziXcguvfhSHJVMC4NLaWujyiEgLDQpT/+To6NiowSkjIwMbNmzA2bNnkZycDDc3N+zevbvWuNTUVLz99ts4ffo0LC0tERYWhjlz5sDEhP/6IyLh+Xd1wnszB2D5V0dx41Yx5n0UhwUT+qCnV2uhSyMiDRnMBfvk5GTExsaiQ4cOcHd3r3OMXC7HxIkTUVFRgTVr1iAyMhJbt27FypUrNXrNmTNn4sknn8QHH3yAigreb4uIdKOjswz/m31/Yvqb64/g10NpQpdFRBrS+MxUYwsJCcGQIUMAAAsXLsSFCxdqjYmOjkZxcTHWrl2rXrG9qqoKy5Ytw9SpU9Vn0kaOHIns7Oxaz/fx8cGGDRsAVE+6d3JyQnFxMV577TWsX78e06dP19O7I6KWxsbaFCum98PabWex/8R1rNt5HtdyFHj5qW4w4sR0IoNiMGFKLH74Xy5xcXEIDAysceub0NBQLF26FPHx8Rg1ahQAYOfOnQ89lpNT9T21LC0tMWbMGPzwww+aFU5E9AASYyPMCa+emL5pzyXsSUhH7u1SvDa+N8xNDeavZ6IWr1n935qWlobRo0fX2CaVSuHg4IC0tPqfQi8pKYFSqYSVlRUqKyvx+++/w8vLS6vaVCoVSkpKtDrGv5WWltb4nXSPPdYv9rdaaN+2sJNKsObH8ziRmIPX1sRhwfO+sJWaaXVc9le/2F/9E7rHKpVKvWrBf2lWYaqwsBBSqbTWdplMBrlcXu/j5OfnY+bMmVAqlaiqqkLPnj0xbdo0rWqrqKhAYmKiVsd4kPT0dL0cl+5jj/WL/QWsRcCEQXbYEpuP9BsKLPgkAc8NtIejjfb39GN/9Yv91T8he1yfL7BpFKb27NmDwYMHw9TUVJOnN3nt2rXDrl27dHpMiUSi83WqSktLkZ6eDldXV5ibm+v02FSNPdYv9rcmbwDdu5Zg5TenkX2rBBv35WPuuO7o7m6n0fHYX/1if/VP6B6npKTUa5xGYWru3LmwsrLCsGHD8OSTT6Jv376aHEbnpFIpFApFre1yuRwymUyAiu4TiUSwsLDQy7HNzc31dmyqxh7rF/t7X0cLC6ya/QhWfH0MF9PysXLzacx8ugeG+HfQ+Jjsr36xv/onVI/rc4kP0HBphC1btmDEiBE4cOAAJk+ejIEDB2LVqlW4cuWKJofTGTc3t1pzoxQKBfLy8uDm5iZQVUREDWNtYYLlUwMR3LMtqpQqfPTDGXwbk8h7+hE1URqFKT8/PyxduhQHDx7Ep59+Cj8/P3z33XcICwtDWFgYvvrqK+TmNv7d0YODg5GQkIDCwkL1tpiYGIjFYgQFBTV6PUREmpIYG+HVZ3vhmSGeAIAf/riC1d+fQkWlUuDKiOjftJqAbmxsjEGDBmHQoEEoLi7GH3/8gZ07d+L999/HBx98AH9/fzz11FMIDQ3VegXy0tJSxMbGAgCysrJQVFSEmJgYAIC/vz9sbW0RHh6Ob775BhEREZg6dSpycnIQFRWF8PDwZnWbGyJqGcRiEcaHeqN1Kwt8uv0sDpzMRL68DG9M9oeFmfYT04lIN3S2MlxycjLOnz+PK1euQKVSwc3NDXfu3MGCBQswdOhQnDhxQqvj5+fnY/bs2Zg9ezaOHTuGGzduqB8nJycDqP7W3qZNm2BkZISIiAh88MEHGDNmDBYuXKiLt0hEJIhH+3bA0il9YW5qhHMpt7Dok3jcLiwTuiwi+ptWZ6auXr2KX375Bbt378b169fRqlUrDB8+HE899RS6du0KADh//jzeeOMNvPnmm3XeS6++XFxckJSU9NBx7u7u2Lhxo8avQ0TUFPl1bo13ZvTHsi+PIC1bjvlrDuKtqYFwtrcSujSiFk+jM1ObNm3CmDFj8Pjjj2PDhg3o0qULPvvsMxw8eBBvvPGGOkgBQLdu3TB58uQGLZpJRES1ebjYIOqVAXCys0BOQQleW3MQKdfvCF0WUYunUZh69913YWJigmXLluHQoUP48MMPMXDgQBgZGdU53sfHBzNmzNCqUCIiAtrYWyLqlQFwayuDvKgcr392CGeuNP4XfojoPo0u8/3xxx9o165dvcd36tQJnTp10uSliIjoX1pZm+HdGUFY8fUxnEu5hWXrj2DuuF4Y0LOt0KURtUganZlqSJAiIiLdszCT4M2X+iKohzMqq1R4/7sT+OUgp1MQCUGjM1OLFi36z/0ikQimpqZwcnKCv78/evbsqVFxRET0YBJjI8x/vjdsrM7j1/ir+OKn87itKMP4UO96r9xMRNrTKEwdPXoUZWVlKCgoAAD1rVru3UzY1tYWSqUSd+7cgUgkQv/+/fHxxx/z3kVERDpmJBZh6shuaCU1xbe/Xca2fckoKavEy091E7o0ohZDo8t8X375JUxMTDBz5kwcPXpU/evIkSOYOXMmzMzM8P333+P48eOYMWMGDh48iI8++kjXtRMREaqvBowd4oUZo7tDJAJ+jb+Kj344jaoqrpZO1Bg0ClNvvfUWgoODMXPmzBo3ELaxscHMmTPRv39/LF++HNbW1njllVfwxBNPYO/evTormoiIagvt1xFzx/lBLBZh/4nrWP3DeVRW8X5+RPqmUZg6e/YsOnfu/MD9nTt3xunTp9WPe/XqhVu3bmnyUkRE1AADe7XDool9YGwkxvHEXHwfewtl5VVCl0XUrGkUpqytrREfH//A/QcPHoSV1f1VeUtKSmo8JiIi/enr0wZvvtgXpiZGSL15F+9sOoWi0gqhyyJqtjQKU8888wz27duHWbNm4fDhw8jKykJWVhYOHz6MWbNm4a+//sIzzzyjHh8bGwtvb2+dFU1ERP+th6cDFk/0g5lEhKRrd/DGZ/GQF90VuiyiZkmjb/PNnDkTZWVl2LRpE/74448a+4yMjDBp0iTMnDkTAHD37l2MGjUKXl5e2ldLRET15tneBpOGOOD7uDtIy5Jj4SeHsHxqP9jb8JvVRLqkUZgSiUSYP38+XnjhBRw+fBjZ2dkAAGdnZwQGBsLOzk491tTUFCNHjtRNtURE1CBOrUzw5pTeWLHpFDJzi7Dwk0N4Z3oQWttaCF0aUbPR4DBVWlqK5557Dk8//TTGjRuH4cOH66MuIiLSEWcHS7w3cwAWr0vAjfxiLPr0EFZMD4KTnaXQpRE1Cw2eM2Vubo7MzEyurktEZEBa21rgnRlBcLa3RO7tUiz6NB43bhULXRZRs6DRBPQBAwbg0KFDuq6FiIj0yN7GHO/MCEJbByvculOKRZ8eQnZekdBlERk8jcLUjBkzkJ6ejvnz5+PEiRPIycnBnTt3av0iIqKmxU5mjndnBKGdoxXy5WVY9OkhZOYqhC6LyKBpNAH9iSeeAACkpKRg9+7dDxyXmJioWVVERKQ3raRmeGd6fyxeF4+Mmwos+jQeK6b1Q3snqdClERkkjcJUREQE50wRERkwG2tTrJgehMXrEpB+oxCvfxaPFdOC0KENAxVRQ2kUpl555RVd10FERI1MZlUdqP5vXQLSsuXVZ6im90NHZ9nDn0xEahrNmfo3hUKBqire+4mIyNBILU3w9vR+8HCRQVFSjsXrEpBxo1DosogMisZh6vz585gyZQp69OiBgIAAHDt2DABQUFCA6dOn4+jRozorkoiI9MfawgTLpwXBw0WGwuJyLP48AddzOCmdqL40ClOnTp3Cs88+i4yMDDz55JNQKpXqfba2tigqKsIPP/ygsyKJiEi/rMwleGtqP3R0luKO4i4Wr4tH9i0um0BUHxqFqdWrV8Pd3R179uxBZGRkrf0BAQE4e/as1sUREVHjsbYwwfKp/dDeyRoFhXfxxmcJyCkoEbosoiZPozB1/vx5jBo1CiYmJnV+q8/R0RG3bt3SujgiImpcMitTvD2tn3phz9c/i0fe7VKhyyJq0jQKU8bGxjUu7f1bTk4OLCx4E00iIkPUytoMK6b3Qxt7S+QWlOCNdfHIlzNQET2IRmGqR48e2Lt3b537SkpKsGPHDvTp00erwoiISDh2MnOsmBaE1rYWuHGrGIvXJeC2okzosoiaJI3C1KxZs3DhwgW8/PLLiIuLAwAkJSVh27ZtGDVqFAoKCjBjxgydFkpERI3LoZU5VkzrB3sbc2TmFmHJ54ehKCkXuiyiJkfjM1NffPEFMjIysGDBAgDAypUr8X//939QKpX44osv0LlzZ50W2pjy8/MRFham/hUcHIynnnpK6LKIiBqdk50lVkzvB1upKdJvFGLZ+iMovVspdFlETYpGK6ADQGBgIPbu3YvExESkp6dDpVKhXbt28PHxMfhbzdjZ2WHXrl3qx6+++iq8vLwErIiISDjO9lZ46+V+WPjJISRl3MY7G49hyZQASIyNhC6NqEnQegV0b29vhIaG4vHHH0e3bt30FqQyMjKwZMkShIWFoUuXLhg+fHid41JTUzF58mT4+voiKCgIUVFRKC/X/LR0UVER9u/fj7CwMI2PQURk6Dq0keLNl/rCzMQIZ67k4f1vT6Kq6sFfRCJqSTQ+MwUAKSkpuH79OuRyeZ37dXlpLDk5GbGxsejRoweUSiVUKlWtMXK5HBMnToSrqyvWrFmDnJwcrFy5EmVlZViyZIlGr7t371707NkTjo6O2r4FIiKD5tXBFosnB+DN9Udw+PwNrN12Fq884wux2LCvRhBpS6Mwde3aNcyfPx/nzp2rM9QAgEgk0mmYCgkJwZAhQwAACxcuxIULF2qNiY6ORnFxMdauXQsbGxsAQFVVFZYtW4apU6eqA9HIkSORnZ1d6/k+Pj7YsGFDjW0//fQTnn76aZ29DyIiQ9bD0wGvje+FlZuO48/j12BpLsGUJ7sa/PQOIm1oFKaWLFmCK1eu4PXXX0fv3r0hlUp1XVctYvHDr0jGxcUhMDBQHaQAIDQ0FEuXLkV8fDxGjRoFANi5c2e9XjM7OxuJiYkYOnSoRjUTETVHgd2cMWtsT3wYfRq74lJhZSFB+FDOK6WWS6MwderUKUydOhXjx4/XdT1aSUtLw+jRo2tsk0qlcHBwQFpaWoOPt2vXLgwbNgzm5uZa16ZSqVBSotvbMpSWltb4nXSPPdYv9le/9NnfwK72uP24FzbtScJ3MZdhYqTCY33b6/x1mjJ+fvVP6B6rVKp6nXXVKEy1atUK1tbWmjxVrwoLC+s8SyaTyR44r+u/7Nq1C2+99ZYuSkNFRQUSExN1cqx/S09P18tx6T72WL/YX/3SV3872gCP+Fgj9oICX/+ahCJ5Hrq2b3l3v+DnV/+E7LGJiclDx2gUpsLDw/Hzzz/jueeeg5FR8/1qbExMjM6OJZFI4OHhobPjAdVJPT09Ha6urjo5e0a1scf6xf7qV2P0t3NnFUwtLuP3Y5nYefg2unh2RJeOtnp5raaGn1/9E7rHKSkp9RqnUZhydXWFUqlEWFgYRo8eDScnpzpD1bBhwzQ5vMakUikUCkWt7XK5HDKZrFFr+TeRSKS3+xWam5vzXoh6xh7rF/urX/ru74yn/VBYUokjF25i1ZazeG/mAHRoo/+5tE0FP7/6J1SP6/vFCo3CVGRkpPrP77333gML0NdlrQdxc3OrNTdKoVAgLy8Pbm5ujVoLEVFLYSQWYd7zvfF/6xKQmF6AN788jKhXguHQimdrqGXQKExt3rxZ13XoRHBwMNatW1dj7lRMTAzEYjGCgoIEro6IqPkylRjh/6YE4LU1B5GZW4Q31x/GezMHwMpcInRpRHqnUZjy9/fXdR0PVVpaitjYWABAVlYWioqK1HOa/P39YWtri/DwcHzzzTeIiIjA1KlTkZOTg6ioKISHh3PRTSIiPbO2MMGylwIxf81BXLupwIqvj2LZS4EwkTTfubVEgJYroJeXl+PixYvIz8+Hn58fbG31N+kwPz8fs2fPrrHt3uPNmzcjICAAMpkMmzZtwvLlyxEREQFLS0uMGTOmxmVJIiLSn9a2Fnjzpb5Y+MkhXEjNx/++P4XXnu/NVdKpWdM4TG3evBlr165VT/j+6quvEBgYiIKCAoSGhmL+/PkYM2aMzgp1cXFBUlLSQ8e5u7tj48aNOntdIiJqmI7OMrw+yR9vfnkY8WezsUF2AS+FdRO6LCK90ehGx9u3b8c777yDAQMGYMWKFTVuKWNra4u+fftiz549OiuSiIgMS49ODogc5wcA+DkuDb8eavjCyUSGQqMw9fXXX2Pw4MH44IMPMGjQoFr7u3btiuTkZK2LIyIiwxXc0wUTHvcGAHzx03mcvJwjcEVE+qFRmMrIyEBwcPAD99vY2ODOnTua1kRERM3EmJBOGNKnPZQq4L3NJ5B+o1Dokoh0TqMwJZVKcfv27QfuT0lJgYODg8ZFERFR8yASiTBjTA90c7dH6d1KvLXhCG4XlgldFpFOaRSmgoODsXXrVhQW1v4XRnJyMrZt24aQkBCtiyMiIsMnMRZj0aQ+aOtgibzbpXj766MoK68UuiwindEoTM2ZMwdVVVUYPnw4PvzwQ4hEIvz000+YN28eRo8eDVtbW8yYMUPXtRIRkYGytjDBkhf7wtpCgivX7uDD709DqVQ9/IlEBkCjMOXo6IgdO3ZgwIAB+O2336BSqbBr1y4cOHAATzzxBLZu3arXNaeIiMjwONtb4Y3JATA2EiH+XDa++a1xbzlGpC8arzNlZ2eHFStWYMWKFSgoKIBSqYStrS3EYo3yGRERtQBd3ezwyjM9sfr7U/hxfzI6OFljYK92QpdFpBWdJB9bW1vY29szSBER0UOF9G6Hpwd3AgCs2XoGydcf/IUmIkPA9ENERI3u+ce80dvbEeWVSrzz9THcVvAbfmS4GKaIiKjRicUizHuuF9o6WOGWvAzvbjyOikql0GURaYRhioiIBGFpLsHiF/xhaWaMxPQCfL7zXI3bkxEZCoYpIiISjEtra8x7vjdEImDvkQz8djhd6JKIGoxhioiIBNXb2xETHu8CAPhi53mcT70lcEVEDVOvpRGOHz+u0cH79Omj0fOIiKhlGT3IA1ez5Ig7k4WVm45jdeQjaN3KQuiyiOqlXmFq/PjxEIlE9T6oSqWCSCRCYiIXZCMioocTiUR4ZawvMvOKkJYlx3ubj2NlRH9IjI2ELo3ooeoVpjZv3qzvOoiIqIUzMzHGool9ELk6Fleu3cGGny9i2qjuQpdF9FD1ClP+/v76roOIiAhOdpaY+6wf3tpwFL/GX0XnDq24Qjo1eZyATkRETUqfLk4YO8QTALD2x7PIuFEocEVE/03je/PdvXsXe/fuxaVLl6BQKKBU1lxsTSQS4Z133tG6QCIiannGPdoZSRm3cSY5D+9uOob/zXkEFmYSocsiqpNGYSorKwsTJkxAVlYWpFIpFAoFZDIZFAoFqqqq0KpVK1hY8FsYRESkGSOxCPOe74U5//sLWXnF+PiHM1gwoXeDvgxF1Fg0uswXFRWFoqIibN26FTExMVCpVFi9ejVOnz6NefPmwczMDBs2bNB1rURE1ILIrEyxYGIfGBuJEH8uGz8fTBO6JKI6aRSmjhw5gnHjxqF79+4Qi+8fwsTEBC+++CL69u3LS3xERKS1zh1sMeVJHwDAxt0XkXL9jrAFEdVBozBVVlaGtm3bAgCsrKwgEomgUCjU+3v27ImTJ0/qpkIiImrRngjqiMBubVBZpULUtydQUlYhdElENWgUptq0aYOcnBwAgLGxMRwdHXHmzBn1/pSUFJiamuqkQCIiatlEIhFmPeMLh1bmuHGrGJ9t5w2RqWnRKEz17dsX+/btUz8eOXIkNm3ahMWLF+P111/Hli1bMGjQIJ0VSURELZuVhQnmPdcLYrEIf53KxP4T14UuiUhNo2/zvfzyyzh//jzKy8thYmKCadOmITc3F3v37oVYLMbw4cOxcOFCXddKREQtWJeOdnj2US98+9tlfLbjHLw6tIJLa2uhyyLSLEw5OzvD2dlZ/djU1BQrVqzAihUrdFYYERHRv40J8cS55Fs4l3ILUd+cwKpZwTCR8P59JCyNLvMtWrQIZ8+efeD+c+fOYdGiRRoXRUREVBcjsQhzn/WDzMoEV7ML8fXui0KXRKRZmNq5cyeuXbv2wP2ZmZn46aefNK2pSdi2bRtGjBiB0NBQvPrqqygrKxO6JCIiAmAnM8eccD8AwO5DV3Hqcq7AFVFLp5d78+Xm5sLMzEwfh24UKSkpWLduHb777jv89ttvkEql+Prrr4Uui4iI/tbb2xHD+3cEAHz0wykUFpcLXBG1ZPWeM/Xnn3/W+Abf1q1bkZCQUGucQqFAQkICfHx8dFPh3zIyMrBhwwacPXsWycnJcHNzw+7du2uNS01Nxdtvv43Tp0/D0tISYWFhmDNnDkxMTOr9WsnJyejWrRukUikAoH///li9ejWmT5+us/dDRETamfhEF5y5kofM3CJ8+uNZ3m6GBFPvMJWamoqYmBgA1Wt+nD17FhcuXKgxRiQSwcLCAn369NH5t/mSk5MRGxuLHj16QKlU1rnGiFwux8SJE+Hq6oo1a9YgJycHK1euRFlZGZYsWVLv1/Ly8sK7776LnJwc2NvbY+/evcjOztbl2yEiIi2ZmRjj1Wd7Yd7HcYg/l40DJzMR0rud0GVRC1TvMDV16lRMnToVANC5c2esWLECI0aM0Fth/xYSEoIhQ4YAABYuXFgryAFAdHQ0iouLsXbtWtjY2AAAqqqqsGzZMkydOhWOjo4AqtfFqisc+fj4YMOGDXBzc8Orr76K6dOnQyKRoG/fvjA21uiLj0REpEce7Www7u/lEj7feQ4+bnZobWshdFnUwmiUEC5fvqzrOh7qn/cAfJC4uDgEBgaqgxQAhIaGYunSpYiPj8eoUaMAVE+gf5iwsDCEhYUBAH777Te4ublpVvjfVCoVSkpKtDrGv5WWltb4nXSPPdYv9le/Wkp/Hw9oi2MXbuDKdTk++O4E/m9y9eKe+tZS+iskoXusUqnqdelYq9Mt169fR1xcnPosj7OzM4KDg9GunTCnWdPS0jB69Oga26RSKRwcHJCW1rC7jd+6dQv29vZQKBRYv349pkyZolVtFRUVSExM1OoYD5Kenq6X49J97LF+sb/61RL6+6ivOa7eKMSl9NvY9PNx9PVqvMU8W0J/hSZkj+sz51rjMLVy5Ups3rwZSqWyxnaxWIyJEydiwYIFmh5aY4WFhepJ4/8kk8kgl8sbdKz58+cjNzcX5eXlGDt2LB5//HGtapNIJPDw8NDqGP9WWlqK9PR0uLq6wtzcXKfHpmrssX6xv/rV0vpbAhts+OUyDpxTIDTYB056vtzX0vorBKF7nJKSUq9xGoWpr776Chs3bsSjjz6KF154Ae7u7gCqJ6lv3LgRGzduhKOjIyZNmqTJ4ZsEXS+FcG9yvj6Ym5vr7dhUjT3WL/ZXv1pKf58M9sSxS7dwPvUWNvyShOVT+zXK5b6W0l8hCdXj+n47VKN1prZu3YqQkBB89NFH6NGjB6ysrGBlZYUePXpg9erVGDRoEKKjozU5tFakUikUCkWt7XK5HDKZrNHrISKixiMWi/DKM74wNTHCuZRb2HskXeiSqIXQKExlZWWhf//+D9zfv39/ZGVlaVyUptzc3GrNjVIoFMjLy9N6AjkRETV9bewtMSHUGwDw9e6LyL2t2y/+ENVFozBlZ2f3n9/ou3z5MmxtbTUuSlPBwcFISEhAYWGheltMTAzEYjGCgoIavR4iImp8w/u7wdvVFqV3q7B265k61yUk0qV6h6njx4+joKAAAPDYY4/hxx9/xBdffFHj6/4lJSX44osv8OOPP2o9YfvfSktLERMTg5iYGGRlZaGoqEj9+F5d4eHhsLS0REREBA4dOoTt27cjKioK4eHh6jWmiIioeROLRZg11hcmxmKcvpKHfccffC9ZIl2o9wT0CRMmICoqCiNGjMDs2bORmJiI//3vf/j444/RunVrANX35KusrERAQABmzZql00Lz8/Mxe/bsGtvuPd68eTMCAgIgk8mwadMmLF++HBEREbC0tMSYMWMQGRmp01qIiKhpc2ltjece64yvd1/CV79cRJ8uTpBZmQpdFjVT9Q5T/zxNam5ujk2bNuHPP/+ssc5U//798cgjjyAkJETn90dycXFBUlLSQ8e5u7tj48aNOn1tIiIyPGHB7vjrVCauZhfi690XMSfcT+iSqJnSatHOIUOGqG/xQkRE1JQYGYkxY0wPvLbmIPYdv47Bfdqjm7u90GVRM9SgCei8GzcRERmSzh1s8WhfVwDAZ9vPoqJS+d9PINJAg85MzZ8/H/Pnz6/XWJFIhEuXLmlUFBERka5MfNwbh89n43pOEX6KTcHTgz2FLomamQaFqX79+sHV1VVPpRAREemelYUJpjzpg/9tOYXo35MwwLctnOwshS6LmpEGhamnnnoKI0aM0FctREREejHQzwV/HruGcym3sG7HOSx9sS+nrpDOaLRoJxERkSERiUSYPro7jI3EOHk5F8cv5QhdEjUjDFNERNQiuLS2Rlhw9a3F1v98gZPRSWcYpoiIqMV4ZognWlmb4satYvxyMO3hTyCqh3qHqcuXL3O+FBERGTQLMwkmPF59I+Qf/kzCbUWZwBVRc8AzU0RE1KKE9G4PDxcZSsoq8e1vl4Uuh5oBhikiImpRxGIRXn6qOwDgj2MZSM28I2xBZPAYpoiIqMXx7miL4J5toVIBX/1yscb9Z4kaimGKiIhapImPd4GxkRjnUm7h9JU8ocshA8YwRURELVJrWws8EdQRALBp9yUolTw7RZphmCIiohbr6cGdYGFmjLRsOQ6eyRK6HDJQDFNERNRiyaxMMWqQBwDg25hELuRJGmGYIiKiFi1sgDtaWZviZn4J9h5JF7ocMkAMU0RE1KKZmRpj3DAvAED0H0kovVspcEVkaBimiIioxRsa0AFt7C0hLyrHnvirQpdDBoZhioiIWjxjIzHGDvEEAOz4KwVlPDtFDcAwRUREBGCgnwva2FmisLgcexJ4dorqj2GKiIgIgJGRGM/w7BRpgGGKiIjob4N6ucDJzgLyonL8djhd6HLIQDBMERER/c3oH3Ondv6VgorKKoErIkPAMEVERPQPj/i1g53MDLcVdxF7KlPocsgAMEwRERH9g8RYjCcHuAMAdvyVynv20UMxTBEREf3Lo307wNzUGNdzFDh5OUfocqiJY5giIiL6F0tzCR4LdAVQ/c0+ov/S4sPUkiVLMGDAAHh5eTVoHxERNW9PDnCDkViEC6n5uHLtttDlUBPW4sPUiBEjsHPnzgbvIyKi5s3exhyP+LkAqP5mH9GDNLkwlZGRgSVLliAsLAxdunTB8OHD6xyXmpqKyZMnw9fXF0FBQYiKikJ5eXmDX69Pnz6wt7dv8D4iImr+nnqkeiJ6wrls3MwvFrgaaqqaXJhKTk5GbGwsOnToAHd39zrHyOVyTJw4ERUVFVizZg0iIyOxdetWrFy5spGrJSKi5qyjsww9PR2gVAF7EtKFLoeaKGOhC/i3kJAQDBkyBACwcOFCXLhwodaY6OhoFBcXY+3atbCxsQEAVFVVYdmyZZg6dSocHR0BACNHjkR2dnat5/v4+GDDhg36exNERNRsDO/vhtNX8vDnsQw891hnmEqMhC6JmpgmF6bE4oefLIuLi0NgYKA6SAFAaGgoli5divj4eIwaNQoAON+JiIi01svbEa1bmSP3dikOns7EEP8OQpdETUyTC1P1kZaWhtGjR9fYJpVK4eDggLS0NIGq+m8qlQolJSU6PWZpaWmN30n32GP9Yn/1i/3VnSG922LLHyn45WAq+vk4AGB/G4PQPVapVBCJRA8dZ5BhqrCwEFKptNZ2mUwGuVzeoGMtXLgQCQkJAIDg4GAEBATg/ffff+i+hqqoqEBiYqJGz32Y9PR0vRyX7mOP9Yv91S/2V3ttratgJAbSshXYF38WzrYm6n3sr/4J2WMTE5OHjjHIMKVL/zVpXZcT2iUSCTw8PHR2PKA6qaenp8PV1RXm5uY6PTZVY4/1i/3VL/ZXt/yvKHH4Qg4ybpticFBn9rcRCN3jlJT6LYlhkGFKKpVCoVDU2i6XyyGTyQSo6OFEIhEsLCz0cmxzc3O9HZuqscf6xf7qF/urG6H93HD4Qg4OnbuJl0b2wL0f7eyv/gnV4/pc4gOa4NII9eHm5lZrbpRCoUBeXh7c3NwEqoqIiJqzHp0c4NDKHMWlFThy/obQ5VATYpBhKjg4GAkJCSgsLFRvi4mJgVgsRlBQkICVERFRcyUWizCkT3sAwB/HMgSuhpqSJhemSktLERMTg5iYGGRlZaGoqEj9uKCgAAAQHh4OS0tLRERE4NChQ9i+fTuioqIQHh6uXmOKiIhI1wb/HabOpdzCLXmZwNVQU9Hk5kzl5+dj9uzZNbbde7x582YEBARAJpNh06ZNWL58OSIiImBpaYkxY8YgMjJSiJKJiKiFcLS1gI+7HS6k5uPQ2RvwchC6ImoKmlyYcnFxQVJS0kPHubu7Y+PGjfoviIiI6B9CerXDhdR8xJ25Ac8hNkKXQ01Ak7vMR0RE1JQF9XCGibEYWXnFyC6oELocagIYpoiIiBrAwkyCvt3aAADOXdXtnS3IMDFMERERNVBI73YAgPMZJaisUgpcDQmNYYqIiKiBfDs5QGZlgpK7SpxNzhe6HBIYwxQREVEDGRmJ0b+7EwAg9ky2wNWQ0BimiIiINBDs6wwAOHk5D0Ul5QJXQ0JimCIiItKAaxtrONpIUFmlwsGzPDvVkjFMERERaahHx+qb7x44cV3gSkhIDFNEREQa6uZqAZEISEwvQPatIqHLIYEwTBEREWnI2twIPTzsAAB/ncwUuBoSCsMUERGRFoJ9qxfw3H/iOlQqlcDVkBAYpoiIiLTQu3NrmJsaIaegBCmZd4QuhwTAMEVERKQFUxMj+Hq2BgCcuJQjcDUkBIYpIiIiLfXxdgQAHDiZiYrKKo2Pky8vxZqtZ5CUUVBje1l5Jb7+5SKOX7qpVZ2kH8ZCF0BERGToens7QmIsxo38Yrz99TEseymw3s9NuX4HienV4enAyetIvn4Hfx7LwIth3eDjbofConIcT8zBrrhU/Hn8GjYvfRRGRjwX0pQwTBEREWmpldQMCyf0wfKvjuLU5VykXL8Dj3Y2D33eb4fTsW77WSj/NW9dqQK++Ol8rfGFxeW4eDUf3T0cdFQ56QLDFBERkQ74d3XCIz1dEHs6E3sSrmLW2J51jjt1ORdR3xxHeaUSFZVKAEB3D3vIrEwhEgGdO9gi7nQmLmfcrvE8YyMxKquUWPblEYjFIliZSzB9dA/8cjANYrEIkeP8YGNtqvf3SbUxTBEREenIUP/2iD2diZOXc3DtZiFcWltDLBap92ffKsI7m47hbvn9eVX2NuZYPrVfjXHD+3fEqm9PIu5MFgDA0dYCk4Z3wXubT6D87wBWVl6F5V8dVT9n5ebjmDG6O5wdrGDMy4CNimGKiIhIR7xcW8FILEJB4V1EvH8AU0d2w5A+7VFcVoELqflY9d3JWs/p29WpRpACAJFIhFef64XZ4T0hFougUgHGRiK0sbfEjVvFAAC/zq1x6nKu+jkX0/IR8f4BPN7PFdNH99DvG6UaGKaIiIh0xMzEGO2drHE1uxAA8PnO89j5Vwpyb5eqx1hbmGDaqG4QQYR9J67hmSGedR5LLBbBRGxUY9uiiX2wftcFjBvmhc6utvj2t0QUl1UiNfMOkq/fAQDsSUjHc495Q2ppop83SbUwTBEREemQk52lOkwBqBGkAro6YdEkfxj9fSZqQM+2DTp2R2cZVkwPUj+eNLwrAOCrXy6qwxQA/BSbggmPd9GkfNIAL6oSERHp0PhQb7i2kdbY5tLaCn19nPDqc73UQUqXnnrEHZ3a2aBLR1sAwPYDKUjLkuv8dahuPDNFRESkQ+0crbFm3iBs+vUSth9IxsynfTEsoINeX9NWaob/zXkEALBs/RGcSMzB1j+vwLN9KwCAWAz06eIEawsTnLycgw5OUuQUFCOwm7Ne62opGKaIiIj0YHyoN556xB0yq8ZdrqCbux1OJOYg/lw24s9lq7fvP3EdphKjGksuLJkSgD5dnBq1vuaIl/mIiIj0QCwWNXqQAoD2TvcvMRobiTGwlwsA4Gp2Ya21q2JPZdV4fDEtH0u/PIxV355ESVmFVnVUVCqxZusZHDqb9fDBBo5hioiIqBlp72it/rNbWylefbYX2jpY1Tk261ZRjcdb/7yCU5dzEXs6E7/GX9Wqjn3Hr+H3oxl4b/MJrY5jCHiZj4iIqBlxaGWu/rONlRkAoEtHW2TlFdUam3L9Dn6Nv4qb+cUoLq3AqaT761b9djgdYlH12bXMXAX692hbr1vkAMCfxzLwbUyi+rFKpYJIVHPifWWVEtv+vIKenVujcwdbXEzLx+9HM6BUqdDF1Rah/TrWOu6hs1m4dacMdysq0bqVBZKv30FlpRJPD/as8b4bG8MUERFRMyISieBoa4GcghKE9G4HAOjeyQF/HLtW5/h1O87V2mZsJELe7VJs/PWSetv2Ayn45YOwh75+VZUSH/1wpsa2krJKWJpLamz7OS4VW35Pwpbfk/DzqiexcfdF9WXIv05mIrCbM0z+scxWlVL1wLNcpiZGmPKkz0Nr05cWH6aWLFmCAwcOIDc3F0lJSertSqUS48aNQ1lZGZRKJTp27Ih33nkHVlZ1nyolIiJqKt6d0R8ZNwvR29sRADDAty2USiUuphXg96MZ//lcawsJXFpbIzG9QKPXVtWxraCwrFaYOnz+hvrP6TcKcafobo39N/OL0b61mfpx/p1SPEh2XrFGtepKi58zNWLECOzcubPWdrFYjA0bNmDXrl345Zdf0KZNG6xfv16AComIiBrGoZW5OkgBgJFYhJDe7TFumJd6Wzd3ezz/WOdaz62oVMLJzqLO41ZWKfH93suY/b+/8H/rEvBR9GmoVHXFp5pmRO3H8Us3AQBrt53BnNV/1ZgMH38uG0Ul1RPebf6etH8u5RYWfnYEp9OKse9EJqas+OOBx79ZwDBVQ0ZGBpYsWYKwsDB06dIFw4cPr3NcamoqJk+eDF9fXwQFBSEqKgrl5eUNfr0+ffrA3t6+zn33zkIplUqUlpbWut5LRERkSGRW928xY2wkwuA+7WFlLkFPTwfMHusLAIgY0wNt7CzrfH7u7RJs+T0JaVlynEnOw5/Hr9VY4f2/rNx0HFVKFfYeyUBqZs0FRc9cyUPx398evDcv65vfEnE1W4FdR27ji12J/z5cDTfzS+oV6vSlyV3mS05ORmxsLHr06AGlUllnc+RyOSZOnAhXV1esWbMGOTk5WLlyJcrKyrBkyRKd1jNp0iQkJibC09MTCxYs0OmxiYiIGpPEuOa9/uxtzPH1kmEwNhLD2EiMft2dYW5qjL9OZdb5/Ju3SmptS7yaD0fbus9k/VN5pRJFJTVPevh6OuDMlTwk/eMslYeLDU4k5tTn7dw/dkUVCgrLYCcTZhJ6kwtTISEhGDJkCABg4cKFuHDhQq0x0dHRKC4uxtq1a2FjYwMAqKqqwrJlyzB16lQ4Olaf2hw5ciSys7NrPd/HxwcbNmyoVz0bN25EVVUVVq1ahS1btuCll17S8J0RERE1PWYm96OAhVn1vKYHnZk6l5JXa9uhs9l4xM8F5ZVKmBiLUV5RVedzW7cyR+7tmmFsUK921fOlFNXzpUwkRjWWdqiP1rYWyC0owc38Eoape8Tih195jIuLQ2BgoDpIAUBoaCiWLl2K+Ph4jBo1CgDqnAulCSMjI4wcORJz585lmCIiombhXnCqi+MD5kxtP5BSa9vRizcxbvEelN6thFgsfuDlttzbpZj7YVyNbV062sLb1VY9Gd3aQgIn+4ef5fonZztL5BaUIKegBF3d7Br0XF1pcmGqPtLS0jB69Oga26RSKRwcHJCWlqaT1ygoqP4Wg62tLVQqFfbu3YtOnTppfDyVSoWSktqnR7VRWlpa43fSPfZYv9hf/WJ/9ctQ+/vC8M74NSEDzwzu+MCfSyZiFXp3dkB+YRnullfBxtoUl67WXD3dw0WKlMxCAEBxWSUAQFmlrDHGxcESN/JLUKWsHbACfRxhbQa0a22Bw39vszA1hqNMgk7tZMi4qYBYJEJZed1nugDg2aEesLE2xY38Ijjbmer852xd62PVxSDDVGFhIaRSaa3tMpkMcnnD7pK9cOFCJCQkAACCg4MREBCA999/H/n5+ViwYAEqKv6eEOfhgcWLF2tcc0VFBRIT/3sCnabS09P1cly6jz3WL/ZXv9hf/TK0/raXAtMfs8PtnAzc/o+pScP9TAHcvx3OcXsVfj1+BwDQ1k6C54OlkBdbYPWumw88xnOPyGAkluGdrTWn3HRpZ45Hu0tw+fJlVJXeD0AiVTmSk5Pw3ABrANa4klWKLbH5AAArMzHuVqhQUVUdzKzNjeDpUAagDNMfs0NxwXUkaraaw38yMTF56BiDDFO6tHLlyjq3d+rUCTt27NDZ60gkEnh4eOjseED1v4bS09Ph6uoKc3PhVn5tzthj/WJ/9Yv91a+W1t9yST5+PX4KAOBoL4O3tzeUShWMd+egsqruS3tenp7V60v9K0y5tLGDt7c3AEBiLcf2hGMAgNZ2MvV2ADC1zgf+DlO2suqFSPH3a0kkxjXG6kNKSu3LmnUxyDAllUqhUChqbZfL5ZDJZAJU9HAikQgWFg27Dlxf5ubmejs2VWOP9Yv91S/2V79aSn9dne8HplbS++9ZZmWKfHlZnc8xt7CAhXntuVn2Npbq53dwvv8NQwszkxq9bOd0/7KhWCwG/nHJTZ8/V//5GvXR5NaZqg83N7dac6MUCgXy8vLg5uYmUFVERETN1z/vfffPOeaarMBobXn/0pn0H38u+nutqXskxvdjSundSg1eqXEYZJgKDg5GQkICCgsL1dtiYmIgFosRFBQkYGVERETNk7GR7oJNReX9M07/PPvz73Wo/olhqgFKS0sRExODmJgYZGVloaioSP343jfswsPDYWlpiYiICBw6dAjbt29HVFQUwsPD1WtMERERkW6Zm1bPDurhcf/OIQN7tXvo87q517zTSEfnml8ia/33Wa/enWv/DG9rV33mKti3bcOKbURNbs5Ufn4+Zs+eXWPbvcebN29GQEAAZDIZNm3ahOXLlyMiIgKWlpYYM2YMIiMjhSiZiIioRfj41YE4n3ILg3rfD1Dhw7zgZGeBjs4yvPpRXJ3Pe218bxy5cAOuzlLcvFWMHp0cauyPemUATiTmYmAvl1rPHfeIHQqrbDAkwA1/HL+m2zekI00uTLm4uCApKemh49zd3bFx40b9F0REREQAACc7Szj9a3V0U4kRHu3rCgDwat8KSddu13qejbUpHgusHtO5g22t/XYyczzat0Odr2llZoQ+3i7qs2JNUZO7zEdERERkSBimiIiIiLTAMEVERESkBYYpIiIiIi0wTBERERFpgWGKiIiISAsMU0RERERaYJgiIiIi0gLDFBEREZEWGKaIiIiItMAwRURERKQFhikiIiIiLTBMEREREWlBpFKpVEIX0dydOnUKKpUKJiYmOj2uSqVCRUUFJBIJRCKRTo9N1dhj/WJ/9Yv91S/2t7aCwjJUVCrVj1u3soA2rfl3j3MKStT7xGIRHGzMtSn3ocrLyyESieDn5/ef44z1WgUBgN7+JxOJRDoPaFQTe6xf7K9+sb/6xf7WZis10+nx/t1jR1sLnR6/Pq9fn5/hPDNFREREpAXOmSIiIiLSAsMUERERkRYYpoiIiIi0wDBFREREpAWGKSIiIiItMEwRERERaYFhioiIiEgLDFNEREREWmCYIiIiItICwxQRERGRFhimiIiIiLTAMEVERESkBYYpA5WamorJkyfD19cXQUFBiIqKQnl5udBlNXkZGRlYsmQJwsLC0KVLFwwfPrzOcdu2bcOjjz6Kbt264cknn8SBAwdqjVEoFHj99dfh7++Pnj17YtasWcjNzdX3W2jSfvvtN0yfPh3BwcHw9fVFWFgYfvzxR/z7fursr2ZiY2Px/PPPo2/fvvDx8cHgwYPx7rvvQqFQ1Bi3f/9+PPnkk+jWrRseffRRbN++vdaxysvL8d577yEoKAi+vr6YPHky0tLSGuutGITi4mIEBwfDy8sL58+fr7GPn+GG27FjB7y8vGr9WrVqVY1xhthbhikDJJfLMXHiRFRUVGDNmjWIjIzE1q1bsXLlSqFLa/KSk5MRGxuLDh06wN3dvc4xv/76K/7v//4PoaGh+PLLL+Hr64uZM2fizJkzNcbNmTMH8fHxePPNN7Fq1SpcvXoVL730EiorKxvhnTRNGzduhLm5ORYuXIjPPvsMwcHB+L//+z988skn6jHsr+bu3LmD7t27Y9myZdiwYQMmT56Mn376CbNnz1aPOXHiBGbOnAlfX198+eWXCA0NxRtvvIGYmJgax3r77bexbds2REZGYs2aNSgvL8ekSZNqBbOW7NNPP0VVVVWt7fwMa2f9+vX44Ycf1L+ee+459T6D7a2KDM66detUvr6+qtu3b6u3RUdHq7y9vVU3b94UrjADUFVVpf7zggULVE888UStMcOGDVPNnTu3xraxY8eqXnzxRfXjU6dOqTw9PVUHDx5Ub0tNTVV5eXmpfv31Vz1Ubhjy8/NrbVu8eLHKz89P3Xv2V7d++OEHlaenp/r//RdeeEE1duzYGmPmzp2rCg0NVT++ceOGytvbWxUdHa3edvv2bZWvr6/qiy++aJzCm7iUlBSVr6+v6vvvv1d5enqqzp07p97Hz7Bmtm/frvL09Kzz74l7DLW3PDNlgOLi4hAYGAgbGxv1ttDQUCiVSsTHxwtXmAEQi//7I3/9+nWkp6cjNDS0xvbHH38chw8fVl9KjYuLg1QqRVBQkHqMm5sbvL29ERcXp/vCDYStrW2tbd7e3igqKkJJSQn7qwf3/h6oqKhAeXk5jh49iscee6zGmMcffxypqanIzMwEABw6dAhKpbLGOBsbGwQFBbG/f3v77bcRHh6Ojh071tjOz7D+GHJvGaYMUFpaGtzc3Gpsk0qlcHBw4JwHLd3r37//AnV3d0dFRQWuX7+uHtexY0eIRKIa49zc3Pjf4F9OnjwJR0dHWFlZsb86UlVVhbt37+LixYv45JNPEBISAhcXF1y7dg0VFRW1/n64d0n7Xu/S0tJgZ2cHmUxWaxz7C8TExODKlSuIiIiotY+fYe0NHz4c3t7eGDx4MD7//HP1pVRD7q2xIK9KWiksLIRUKq21XSaTQS6XC1BR83Gvf//u773H9/YXFhbC2tq61vNlMhkuXLig5yoNx4kTJ7Bnzx4sWLAAAPurK4MGDUJOTg4AYMCAAfjggw8AaN9fqVTa4v8OKS0txcqVKxEZGQkrK6ta+/kZ1pyDgwNeeeUV9OjRAyKRCPv378eHH36InJwcLFmyxKB7yzBFRHpx8+ZNREZGIiAgABMmTBC6nGbliy++QGlpKVJSUvDZZ59h2rRp+Prrr4Uuq1n47LPPYGdnh9GjRwtdSrMzYMAADBgwQP24f//+MDU1xaZNmzBt2jQBK9MeL/MZIKlUWuc3buRyea3T9tQw9/r37/4WFhbW2C+VSlFUVFTr+fxvUK2wsBAvvfQSbGxssGbNGvVcNfZXNzp37oyePXvi6aefxqeffoqjR4/ijz/+0Lq/hYWFLbq/WVlZ+OqrrzBr1iwoFAoUFhaipKQEAFBSUoLi4mJ+hnUsNDQUVVVVSExMNOjeMkwZoLquCysUCuTl5dWaK0ENc69//+5vWloaJBIJ2rVrpx539erVWusnXb16tcX/NygrK8PUqVOhUCiwfv36Gqfj2V/d8/LygkQiwbVr19C+fXtIJJI6+wvc77+bmxtu3bpV65JeXfMxW5LMzExUVFTg5ZdfRp8+fdCnTx/1GZMJEyZg8uTJ/AzrkSH3lmHKAAUHByMhIUGd1oHqCZNisbjGtxuo4dq1awdXV9daa/Ls2bMHgYGBMDExAVD930Aul+Pw4cPqMVevXsWlS5cQHBzcqDU3JZWVlZgzZw7S0tKwfv16ODo61tjP/ure2bNnUVFRARcXF5iYmCAgIAB79+6tMWbPnj1wd3eHi4sLgOrLK2KxGL///rt6jFwux6FDh1p0f729vbF58+YavxYtWgQAWLZsGZYuXcrPsI7t2bMHRkZG6NKli0H3lnOmDFB4eDi++eYbREREYOrUqcjJyUFUVBTCw8Nr/fCimkpLSxEbGwug+pR+UVGR+n9cf39/2Nra4pVXXsG8efPQvn17BAQEYM+ePTh37hy+/fZb9XF69uyJ/v374/XXX8eCBQtgamqK1atXw8vLC8OGDRPkvTUFy5Ytw4EDB7Bw4UIUFRXVWGivS5cuMDExYX+1MHPmTPj4+MDLywtmZma4fPkyNmzYAC8vLwwZMgQAMH36dEyYMAFvvvkmQkNDcfToUezevRurV69WH8fJyQljxoxBVFQUxGIxHB0d8fnnn8Pa2hrh4eFCvT3BSaVSBAQE1Lmva9eu6Nq1KwDwM6yhKVOmICAgAF5eXgCAffv2YevWrZgwYQIcHBwAGG5vRap/nycjg5Camorly5fj9OnTsLS0RFhYGCIjI9XJneqWmZmJwYMH17lv8+bN6r9It23bhi+//BLZ2dno2LEj5s6di0GDBtUYr1Ao8O677+KPP/5AZWUl+vfvj8WLF7foQBsSEoKsrKw69+3bt099ZoT91cwXX3yBPXv24Nq1a1CpVGjbti2GDh2KKVOm1Pjm2b59+/Dhhx/i6tWrcHZ2xssvv4wxY8bUOFZ5eTlWr16NXbt2obi4GH5+fli8ePED7wzQUh09ehQTJkzAjz/+iG7duqm38zPccG+//TYOHjyImzdvQqlUwtXVFU8//TTGjx9fY5kDQ+wtwxQRERGRFjhnioiIiEgLDFNEREREWmCYIiIiItICwxQRERGRFhimiIiIiLTAMEVERESkBYYpIiIiIi0wTBERNZKjR4/Cy8sLR48eFboUItIhhikiMlg7duyAl5cXzp8/DwCIjY3FmjVrBK4K+O6777Bjxw6hyyCiRsIwRUTNRmxsLNauXSt0Gfj++++xc+fOWtv79OmDc+fOoU+fPgJURUT6wjBFRPQfVCoVysrKdHIssVgMU1NTiMX8q5eoOeH/0UTULCxcuBDfffcdAMDLy0v96x6lUomNGzfiiSeeQLdu3dCvXz8sWbIEcrm8xnFCQkIwdepUHDx4EKNGjUL37t0RHR0NANi+fTsmTJiAwMBA+Pj44PHHH8eWLVtqPT85ORnHjh1T1zB+/HgAD54z9dtvv6lfKyAgAPPmzUNOTk6t99ezZ0/k5ORgxowZ6NmzJ/r27Yv33nsPVVVVumkiEWnEWOgCiIh0YezYscjNzUV8fDyioqJq7V+yZAl27tyJUaNGYfz48cjMzMR3332HS5cu4fvvv4dEIlGPvXr1Kl599VWMHTsWzzzzDDp27Aig+vJdp06dEBISAmNjYxw4cADLli2DSqXCc889BwB4/fXXsXz5clhYWGDatGkAAHt7+wfWvWPHDixatAjdunXD3LlzkZ+fj82bN+PUqVP46aefIJVK1WOrqqowZcoUdO/eHa+99hoOHz6Mr776Cu3atcOzzz6rkz4SkQZUREQGavv27SpPT0/VuXPnVCqVSrVs2TKVp6dnrXHHjx9XeXp6qn7++eca2+Pi4mptHzRokMrT01MVFxdX6zilpaW1tr3wwguqwYMH19j2xBNPqJ5//vlaY48cOaLy9PRUHTlyRKVSqVTl5eWqwMBA1fDhw1VlZWXqcQcOHFB5enqqPvroI/W2BQsWqDw9PVVr166tccynnnpKNXLkyFqvRUSNh5f5iKjZi4mJgbW1NYKCglBQUKD+1bVrV1hYWNS67Obi4oIBAwbUOo6ZmZn6zwqFAgUFBfD398f169ehUCgaXNeFCxeQn5+PcePGwdTUVL194MCBcHNzw19//VXrOePGjavxuFevXsjMzGzwaxOR7vAyHxE1exkZGVAoFAgMDKxzf35+fo3HLi4udY47efIk1qxZgzNnzqC0tLTGPoVCAWtr6wbVlZ2dDQDqy4j/5ObmhpMnT9bYZmpqCltb2xrbZDJZrXlfRNS4GKaIqNlTKpWws7PDqlWr6tz/74DyzzNQ91y7dg2TJk2Cm5sbFi5ciDZt2kAikSA2NhYbN26EUqnUS+3/ZGRkpPfXIKKGY5giomZDJBLVub19+/Y4fPgw/Pz86gxK9bF//36Ul5fjs88+g7Ozs3p7XauZP6iOf7t3nKtXr9Y6a3b16tUar0NETRfnTBFRs2Fubg4AKCwsrLE9NDQUVVVV+PTTT2s9p7Kystb4utw7K6RSqdTbFAoFtm/fXmcd9Tmmj48P7OzsEB0djfLycvX22NhYpKamYuDAgQ89BhEJj2emiKjZ6Nq1KwDg7bffRv/+/WFkZIQnnngC/v7+GDt2LD7//HMkJiYiKCgIEokE6enpiImJwRtvvIHHHnvsP4997znTpk1DeHg4iouLsW3bNtjZ2SEvL69WHd9//z0+/fRTdOjQAba2tnXO15JIJJg3bx4WLVqE559/Hk888YR6aYS2bdti0qRJOusNEekPwxQRNRvDhg3D+PHj8euvv+Lnn3+GSqXCE088AQB466234OPjg+joaKxevRpGRkZo27YtnnzySfj5+T302G5ubvj444/x4Ycf4r333oO9vT3GjRsHW1tbvP766zXGRkREIDs7G+vXr0dxcTH8/f0fOPl91KhRMDMzw5dffolVq1bBwsICQ4YMwfz582usMUVETZdI9c9z1kRERETUIJwzRURERKQFhikiIiIiLTBMEREREWmBYYqIiIhICwxTRERERFpgmCIiIiLSAsMUERERkRYYpoiIiIi0wDBFREREpAWGKSIiIiItMEwRERERaYFhioiIiEgLDFNEREREWvh/mp0mBvsk41MAAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "history = jnp.stack(history)\n", - "history = history - jnp.min(history)\n", - "ax = sns.lineplot(history)\n", - "ax.set_yscale(\"log\")\n", - "ax.set_xlabel(\"Iteration\")\n", - "ax.set_ylabel(\"Total energy (Hartree)\");" - ] - }, { "cell_type": "code", "execution_count": 6, @@ -256,180 +189,6 @@ "ax.set_xlabel(\"Iteration\")\n", "ax.set_ylabel(\"Total energy (Hartree)\");" ] - }, - { - "cell_type": "code", - "execution_count": 15, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "'cart'" - ] - }, - "execution_count": 15, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "backend = \"pyscf_cart\"\n", - "backend.startswith(\"pyscf_\")\n", - "backend.split(\"_\")[1]" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/ubuntu/miniforge3/envs/jax/lib/python3.10/site-packages/pyscf/dft/libxc.py:771: UserWarning: Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, corresponding to the original definition by Stephens et al. (issue 1480) and the same as the B3LYP functional in Gaussian. To restore the VWN5 definition, you can put the setting \"B3LYP_WITH_VWN5 = True\" in pyscf_conf.py\n", - " warnings.warn('Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, '\n" - ] - } - ], - "source": [ - "from mess import molecule\n", - "from mess.interop import to_pyscf\n", - "import numpy as np\n", - "import numpy.linalg as npl\n", - "\n", - "np.set_printoptions(linewidth=120, precision=6)\n", - "\n", - "mol = molecule(\"water\")\n", - "scfmol = to_pyscf(mol, basis_name=\"def2-TZVPPD\")\n", - "# scfmol.cart = False\n", - "# scfmol = scfmol.build()\n", - "S = scfmol.intor(\"int1e_ovlp_cart\")\n", - "s, U = npl.eigh(S)" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([2.227773e-04, 4.519461e-04, 1.138834e-03, 5.175996e-04, 3.437612e-03, 8.825716e-04, 1.000282e-03, 4.020117e-03,\n", - " 3.256600e-03, 5.913186e-03, 5.233487e-03, 6.822549e-03, 6.770662e-03, 1.884973e-03, 1.011013e-02, 8.236594e-04,\n", - " 1.085325e-02, 8.021005e-04, 1.065478e-02, 7.274325e-03, 7.035361e-04, 2.655205e-02, 1.293354e-02, 1.585945e-02,\n", - " 6.655871e-03, 1.770198e-02, 3.971056e-02, 9.050910e-04, 5.530916e-03, 6.185895e-03, 1.594431e-02, 1.441710e-02,\n", - " 1.686322e-02, 2.529782e-02, 7.736134e-04, 1.433575e-02, 4.884687e-03, 5.009570e-02, 2.044539e-02, 2.933490e-02,\n", - " 4.542135e-02, 6.228932e-02, 1.673826e-02, 1.258803e-02, 1.978192e-02, 6.286715e-02, 2.073157e-02, 6.370569e-02,\n", - " 1.735051e-02, 5.812182e-03, 1.039793e-01, 3.430022e-02, 5.516339e-02, 3.182214e-02, 3.353706e-03, 1.008233e-01,\n", - " 4.393568e-02, 2.208200e-02, 4.042325e-02, 1.120331e-01, 2.074160e-02, 4.688071e-02, 2.036013e-01, 7.769408e-02,\n", - " 5.559280e-02, 5.492349e-02, 8.038952e-02, 4.374833e-01, 1.038853e-01, 2.356155e-02, 1.074352e-01, 5.213506e-02,\n", - " 4.636286e-01, 1.819283e-01, 2.744310e-01, 6.971762e-01, 1.009265e+00, 5.941142e-01, 5.952234e-01, 1.692097e+00,\n", - " 7.889425e+00])" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "np.diff(s)" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "70" - ] - }, - "execution_count": 16, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "S = scfmol.intor(\"int1e_ovlp_sph\")\n", - "s, U = npl.eigh(S)\n", - "np.sum(np.diff(s) > 1e-3)" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(82, 82)" - ] - }, - "execution_count": 12, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "U.shape" - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "True" - ] - }, - "execution_count": 21, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "sc = s[s > 1e-3]\n", - "Uc = U[:, s > 1e-3]\n", - "\n", - "X = Uc @ np.diag(np.power(sc, -0.5)) @ Uc.T\n", - "\n", - "sc = np.where(s < 1e-3, 1.0, s)\n", - "Uc = np.where(s < 1e-3, 0.0, U)\n", - "\n", - "X2 = Uc @ np.diag(np.power(sc, -0.5)) @ Uc.T\n", - "\n", - "np.allclose(X, X2)" - ] - }, - { - "cell_type": "code", - "execution_count": 23, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "(82, 82)" - ] - }, - "execution_count": 23, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "(Uc @ np.diag(np.power(sc, -0.5))).shape" - ] } ], "metadata": {