diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 63ebe7d..b363806 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -6,9 +6,13 @@ name: CI on: # Triggers the workflow on push or pull request events but only for the main branch push: - branches: [ main ] + branches: + - main + - PariserParr pull_request: - branches: [ main ] + branches: + - main + - PariserParr # Allows you to run this workflow manually from the Actions tab workflow_dispatch: diff --git a/examples/Demonstration.ipynb b/examples/Demonstration.ipynb index b9ea721..75d8dcb 100644 --- a/examples/Demonstration.ipynb +++ b/examples/Demonstration.ipynb @@ -44,7 +44,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 1, "id": "36ebc326", "metadata": {}, "outputs": [], @@ -67,9 +67,9 @@ "\n", "# Second way to define the Hamiltonian\n", "# two site Hubbard model\n", - "connectivity = np.array([[0, 1],\n", + "adjacency = np.array([[0, 1],\n", " [1, 0]])\n", - "hubbard = HamHub(connectivity,\n", + "hubbard = HamHub(adjacency = adjacency,\n", " alpha=-0.414, beta=-0.0533, u_onsite=np.array([1, 1]))" ] }, @@ -106,7 +106,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 2, "id": "21b6dc5e", "metadata": {}, "outputs": [ @@ -114,7 +114,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Zero energy: 0\n", + "Zero energy: 0.0\n", "One body integrals in spatial basis: \n", " [[ 0. -1. 0. 0. 0. -1.]\n", " [-1. 0. -1. 0. 0. 0.]\n", @@ -147,13 +147,13 @@ "# Assuming 8-fold symmetry\n", "# Returning output as dense matrix\n", "\n", - "connectivity = np.array([[0, 1, 0, 0, 0, 1],\n", + "adjacency = np.array([[0, 1, 0, 0, 0, 1],\n", " [1, 0, 1, 0, 0, 0],\n", " [0, 1, 0, 1, 0, 0],\n", " [0, 0, 1, 0, 1, 0],\n", " [0, 0, 0, 1, 0, 1],\n", " [1, 0, 0, 0, 1, 0]])\n", - "hubbard = HamHub(connectivity,\n", + "hubbard = HamHub(adjacency=adjacency,\n", " alpha=0, \n", " beta=-1, \n", " u_onsite=np.array([1, 1, 1, 1, 1, 1]))\n", @@ -200,23 +200,24 @@ "outputs": [], "source": [ "from rdkit.Chem import MolFromSmiles, rdmolops\n", + "from moha import HamHuck\n", "\n", "# generate a molecule from SMILES\n", "mol = MolFromSmiles('C1=CC=C2C=CC=CC2=C1')\n", "# generate the connectivity matrix\n", - "connectivity = rdmolops.GetAdjacencyMatrix(mol)\n", + "adjacency = rdmolops.GetAdjacencyMatrix(mol)\n", "\n", "# generate the Hamiltonian\n", "# generating Huckel model\n", - "huckel = HamHub(connectivity, alpha=-0.414, beta=-0.0533)\n", + "huckel = HamHuck(adjacency=adjacency, alpha=-0.414, beta=-0.0533, atom_types=['C' for _ in range(adjacency.shape[0])])\n", "e0_huck = huckel.generate_zero_body_integral()\n", "h1_huck = huckel.generate_one_body_integral(dense=True, basis='spatial basis')\n", "h2_huck = huckel.generate_two_body_integral(dense=True, basis='spatial basis', sym=4)\n", "\n", "# generate Hubbard model\n", "# setting on-site repulsion to 10.84 eV for each atom. \n", - "U = 10.84 * np.ones(connectivity.shape[0])\n", - "hubbard = HamHub(connectivity, alpha=0, beta=-1,\n", + "U = 10.84 * np.ones(adjacency.shape[0])\n", + "hubbard = HamHub(adjacency=adjacency, alpha=0, beta=-1,\n", " u_onsite=U)\n", "e0_hub = hubbard.generate_zero_body_integral()\n", "h1_hub = hubbard.generate_one_body_integral(dense=True, basis='spatial basis')\n", @@ -254,7 +255,7 @@ "output_type": "stream", "text": [ "Energy given by Lieb-Wu formula: -1.040368653394435\n", - "Error of integration: 8.67591249968029e-11\n" + "Error of integration: 8.675915092278734e-11\n" ] } ], @@ -285,7 +286,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "8f87778c", "metadata": {}, "outputs": [], @@ -321,21 +322,10 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "id": "e9f96129", "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2MAAAKmCAYAAADAY0a6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAABcSAAAXEgFnn9JSAACsK0lEQVR4nOzdd3gUVdvH8e/upveEEEoSCIgQQAWUIh0UFDtFsVAEBB6Kr5VmfUBRFBCwYFe6IlUFH1FAei9WpCqBBASEJJBKkt15/1izEpJANm2T8Ptc115JZs6ZuWd3ZjP3zJlzTIZhGIiIiIiIiEipMrs6ABERERERkSuRkjEREREREREXUDImIiIiIiLiAkrGREREREREXEDJmIiIiIiIiAsoGRMREREREXEBJWMiIiIiIiIuoGRMRERERETEBZSMiYiIiIiIuICSMRERERERERdQMiYiIiIiIuICSsZERERERERcQMmYiIiIiIiICygZExERERERcQElY1Im9OvXD5PJRL9+/Yp1uWvXrsVkMmEymYp1uVJ0ZeGzKan97ko3c+ZMTCYTUVFRrg7lipJ9PK1du9Yl678SP/dff/2Vnj17Uq1aNdzc3DCZTDRu3NjVYRVZSe1LMTExjmXHxMQU67JFyislY1Ksxo4d6/IT7NIyYsQITCYTQUFBWK3WfMvdf//9jvdk1qxZ+ZbbsmWLo9z3339fEiHna9u2bQwcOJD69evj7++Pp6cnkZGRNG3alAEDBjBz5kxiY2NLNaaiWLt2LWPHjmXmzJmuDkWclH1CX5DXpU7mTp06xWuvvUbnzp2JiIjA29sbX19foqKi6Nq1Kx9++CGJiYm56hXHd9iRI0d45plnaNasGcHBwbi7u1OlShWuu+46evTowbRp0/j5558LvfzSFhMTw9ixYxk7dqyrQylTDh8+TOvWrVm4cCEnTpwgMDCQKlWqEBoa6urQcrhwn77cy5XK0n6WnJzMm2++yU033USVKlXw8PAgJCSE+vXrc+uttzJu3Dh++OGHS/7vLw1l6T2TwnNzdQAiANWqVaNevXpUq1bN1aEUWMeOHXnjjTc4e/YsP/74I02bNs2z3Pr16x2/r1mzhocffjjPctlXIN3d3WndunWxx5sXwzB48sknefPNNx3TshPMv//+m7i4OHbt2sWMGTN4+OGHy01ys3btWsaNG0f79u0vederPO53V5LQ0FAsFku+8/OaZxgGEyZM4JVXXiE1NdUx3c/PD5PJxJEjRzhy5AhfffUVo0aNYsqUKQwYMKDYYp43bx6DBw/Ose6AgABSU1P59ddf+fXXX1myZAk1a9YsN3cGYmJiGDduHMAlT/oCAwOpV68e4eHhpRSZa33wwQckJSVRp04d1qxZQ0REhKtDuqwqVaoUqFy9evUA8PHxKdb1u7u7O5bt7u6eY15B97OS9ssvv3DnnXfmuADp5eWFYRjs37+fffv2OS6YHj582KV3gsvKeyZFo2RMyoQJEyYwYcIEV4fhlLZt22KxWLBaraxduzbPZGzv3r2cOHGCqlWrcuLEiUs2+cie17x5c3x9fUso6pymTp3qSMTuueceRo8ezQ033ICHhwdg/0ezZs0aFi5ceMmT4vKqPO53V5IdO3Y4daJjGAZ9+vRh3rx5ALRo0YJRo0Zx0003ERQUBMC5c+f44Ycf+PTTT1m2bBlff/11sSVjO3bsoG/fvthsNq677jpefPFFunTp4jie//77bzZv3syXX37J9u3bi2WdZUm3bt3o1q2bq8MoNb/++itg/+4sD4kYwIkTJwpUbt++fSWy/vDw8BJbdnFISkri9ttv59ixY4SGhvLCCy/w4IMPUrlyZQBSUlLYvXs3y5cvZ86cOS6OVioKJWMihRQQEMANN9zA9u3bWbNmDSNGjMhVJjvBuueee1izZg0HDhzg8OHD1KpVK0e5zMxMNm3aBNjvuJUGwzCYMmUKAF26dOHLL7/MVaZWrVrUqlWLAQMGkJaWVipxiRTWxIkTHYnYE088wZQpU3I1vQoICKBr16507dqV9evXs2DBgmJb/7Rp07DZbISFhbF+/XoCAwNzzK9cuTL33HMP99xzj46nCiD77qefn5+LI5HiMn/+fI4dOwbAsmXLuPHGG3PM9/X1pW3btrRt25ZXXnnFFSFKBaRnxqRMKEhHCidOnGDMmDE0atSIwMBAvLy8qF27NgMHDuT3338v0Hp27tzJvffeS7Vq1fDy8qJOnTqMHDkyz2dHCiI7cdq4cWOebcezk7EOHTrQvn37HNMutGPHDlJSUnIsE/5t69+hQ4d8YyhsRxinT592/NO5++67L1ve29s733k//vgjffv2pWbNmnh5eREcHEyrVq2YNm0a58+fdyouKPx2Zz8cnt1sY926dbmeibiwqWVB9ru1a9dy3333ER4ejqenJ6Ghodx8883MmDEj3+cFLo5/9erV3HHHHVSuXBkvLy/q16/PuHHjSE9PL/B7cqG0tDS+/vprBg0aROPGjalcuTKenp5Ur16drl278u233+Zb9+JOFnbt2uXogMDT05PatWvz1FNPkZCQcMkYtm7dSteuXQkNDcXb25t69erx3HPPkZycXKhtKqrTp0/z8ssvA3DzzTfnmYhdrF27drz11lvFFsNPP/0E2I/3ixOxi13qeDpx4gQjR46kYcOG+Pn54evrS8OGDRk1ahQnT550Oq6Cfkfk1WlDVFRUju+ki4+nC4+dgnTg8ccffzB06FCuvvpqvL29CQgI4Prrr+ell17i3LlzBYr/0KFDDBgwgMjISDw9PYmIiGDQoEGO77O87Nu3j8GDB1O3bl18fHzw9vYmMjKSG2+8kWeffdapuzVRUVE53qdx48bleE8u/o4vzOd5cUcXf/zxB4MHD6ZWrVp4enqWeNO4/Lbl4rhOnjzJ448/Tq1atfDy8qJKlSo88MAD+b6f+XXg4cx+li09PZ233nqL9u3bExoaioeHB1WrVqVr166sWLGiUNudfQyHhYXlSsQu5ubmhpvbv/c03n//fUwmE5UqVbrkd7vNZnPsQxc3LXRmPy3N9+zC/eHMmTM89dRTXHXVVXh7e1OzZk0effRR/v77b0f5I0eOMHToUMd+UaNGDZ5++mmSkpLyXcd3331H9+7diYiIwMPDg4CAAGrXrs0tt9zC5MmTiY+Pz7duuWeIFKP//ve/BmA4u2s9/PDDBmA8/PDDec5ftmyZ4efn51i2u7u74evr6/jbw8PDmDVrVq56a9ascZT58ssvDQ8PDwMwAgICHL8DRs2aNY3Dhw87vb0rVqxwLGP79u255lepUsUAjOPHjxtz5841AKNv3765yr3yyisGYHh6ehqpqamO6dnvZ/v27fON4cJtdMapU6cc9Z599lmn6l5o6tSphslkciwrMDDQcHd3d/x93XXXGcePH3cq7sJu99GjR40qVao49g13d3ejSpUqOV7z5893lL/cfvfkk0861mEymYygoCDDYrE4pt10003GuXPnLhn/xIkTDZPJ5Kh/4XvVsWNHIysr6xLvbt5mzJjhWAZgeHt7Gz4+PjmmPf3005esW7NmTWPevHmOzyowMNAwm82O+g0bNjSSkpLyXMYnn3ySo2xgYKDjeIqOjjamTJniWEdRts2ZY3LixImOehs2bHB6vdkK+x1mGIbRoEEDAzBatWpV6PWvXbvWCAoKcsTg4+OT47suODg43+3LLrNmzZoc0wv6HZFX/aZNmxrBwcGOeRcfT4899pij7IX7Vl6++OILw9PT07Esf3//HH9HRkYav//+e656F8b/ww8/OP4X+Pv7G25ubo551atXN+Li4nLV//7773Osx93dPcd7DBj//e9/L/neXKhp06ZGlSpVHMeOr69vjvdk06ZNjrKF/TwPHz7sKDNv3jzHNmfXd+bYKsw+nd++dGFcy5cvN8LCwhxxXfgeBwQEGD/99NMlt+vC49uZ/cwwDOPAgQPG1VdfneP7OTAwMMdnOnTo0AJvb7Zhw4Y59pGUlBSn6iYlJRn+/v4GYMyZMyffct9++60BGBaLxTh69KhjurP7aWm+Z9nzZ82aZURERDj2+wvPo+rXr28kJCQY27dvN0JDQx37wYXHaOvWrfP8nzdu3Lgccfj4+OQ458trX6xIlIxJsSqJZGzbtm2OA/4///mPsXfvXsfBfOTIEceXp5ubm7Fjx44cdS/8Jx4YGGh06NDB8c8+MzPT+OKLLxxfZs2aNXP6xDg5OdnxD/n111/PMe/33383AKNu3bqGYRhGbGys44TjYp06dcoz+SjJZMwwDKNWrVqOk5rvv//e6frLli1zrPuee+4x/vzzT8MwDOP8+fPG7NmzHf+YWrVqleu9LYlkzJn6hnHp/e7tt992LH/w4MHGX3/9ZRiG/TOfOnWq4x/M/fffn+/6g4KCDLPZbDzzzDPG33//bRiGYZw9e9Z48cUXHcv+5JNPLhljXpYuXWoMHjzYWLNmjXH69GnH9OPHjxvjxo1z7JNfffVVrrrZJ8zZJ08DBw50nBCkpKQY77zzjqP+Cy+8kKv+rl27HNveoUMHY+/evYZhGEZGRobx+eefG0FBQY4TiNJMxm699VYDMEJDQ51e54WKkoz169fPUXfy5MnG+fPnnap/9OhRx3vXoEEDY+PGjY5569evN+rVq2cARkhISJ5JR0kkY87Uv1QytmvXLsd+1bp1a+Pnn382DMMwrFar8fXXXxvVqlUzAOOqq67KdRHgwvUHBwcbd999t2O/O3/+vPHFF184vmv69OmTa9116tQxAOOWW24xfv31V8f0tLQ049dffzXGjh1rfPrpp5fctry0b9/+kolcUT7PC5MWPz8/o0WLFjn+v+3fv7/AcZZUMhYcHGy0bt3aEVdmZqaxcuVKx2fZtm3bXMvNLxkzjILvZwkJCUZUVJQB9gti69evN9LT0w3DMIzExERjypQpjhP5adOmFXibDcMwZs6c6YjhgQceMM6cOeNU/SFDhhiA0a5du3zLdO/e3QCMO++8M8f0wuynpfWeZa8jKCjIaNy4sbF161bDMP793s++GPjoo48aNWvWNG666Sbjt99+c8T/9ttvOy5kfvTRRzmWHRMT47i499RTTxnHjh1zzEtMTDQ2bNhgDBs2zNi5c+clt7E8UzImxaokkrFmzZrle2KY7bHHHnMkBBe68Iuqbt26Oe46ZVu5cqWjzIIFC5yK2zAMo1WrVgZg3HbbbTmmv/vuuwZgDBo0yDGtdu3aBmD88ccfjmkZGRmOL7KxY8fmWEZJJ2Pz5s3LceWpZs2aRp8+fYxp06YZmzZtcnxZ5yf7TkCbNm3yTGS//vprx7IXLlxY4LhdnYylpqYaISEhBmA8+OCDedZ96623HOu/+CLAhcdBfidq2f+QO3XqdMkYC2PSpEkGYNx888255l2Y7OR3R/Cpp54yAKNOnTq55t12222XPJ4uvFtc1GQsNDQ019Xe7Ne3336bo1721drOnTs7vc4LFSUZ27dvnyMpyD5Z7dq1qzF+/Hjj22+/NRISEi5ZP/tELjg42JH8Xyg2NtYICAgwAGP48OG55pflZKxLly6OfSqvOw67d+92JPmTJk3Kd/0dO3Y0rFZrrvrZx6O3t7eRmZnpmH7y5ElH3bzu0BfF5ZKxonyeFyYtNWvWzPcudUFcuE/ndzxVqVLFcfJsGAVLxqKjo/P8Drjwez82Njbf+oVNxkaMGOFIKi78rC+0ZMkSx3dIfmXykp6eblxzzTWOODw8PIyOHTsao0ePNhYsWJDjTlZefv75Z0fd7AsGFzpx4oTjosTXX3/tmF7Y/bS03rML958LLwBme+GFFxxlGjZsmOe5Q58+ffL8v/TFF184/qdcqfTMmJRpP//8Mzt27MDd3Z2nn34633J9+/YFYNWqVfk+xzNy5Mg8n9Po1KkTrVq1AuwP7zor+7mgjRs3kpWV5Zi+Zs2aHPMBx3Nj2fPAPsZX9oPgpdV5R7aHHnqI+fPnO3oCO3LkCHPmzOGJJ56gdevWBAcH88ADD+Q5JtIvv/zieFbvhRdeyLO3xbvuuovmzZsD8Pnnn5fglhSvlStXOtqn59dd8LBhwxxd4ue3bZ6ennl27AL2Tl3A/j4WtzvuuAOwj113qXFwnn/++UvGdujQoRxdtCcmJvLdd98B+R9Pt956Ky1btix07Bc6ffo0J0+ezPN18TMZZ86cASAkJKRY1l0Y9erVY926dTRr1gyAhIQEvvzyS55//nluu+02KlWqRIcOHfLsLMcwDEdnIkOGDKFq1aq5ykRERDBkyBCgcN9VrnLxfpNXd+lNmjShe/fuwKW/K5599lnM5tynLtn7bFpaGgcPHnRM9/f3d5T/66+/Cr8RTirOz/PRRx8ttk5C8jueTp48SWZmplPLevrpp/P8DrjtttscPfJm9zhZXAzD4NNPP3Ws/8Jnti7UtWtXAgICOH36NLt27Srw8j09Pfnhhx8c44NmZGSwZs0aXn/9dXr27EmNGjVo0KBBvs9DX3fddY7vv48++ijX/BkzZpCZmUlERAS33367Y3pJ7qfF+Z4NGjSISpUq5Zp+6623On5/6qmn8PT0zLfMxf/zsnu6TUpKcjw7f6VRMiZl2saNGwH7A6/16tWjatWqeb66dOkC2LudzT4pu9hNN92U73qy5+3cudPpGLMTqKSkpBxfYNnji2UnYBf+fuFD0dm/e3t7X/aB4ZJw//338+eff/K///2Pp556ijZt2uDv7w/YT2y++OILmjZtmusfS/Z75ebmlmMbL9a5c+cc5cuD7FgjIyOpW7dunmUsFstl95vsB/bzUr16dYBCP5R88uRJ/vvf/9KyZUsqVaqEm5ub4yHrBg0aAPbe3vLriCMkJIQ6depcMjYgR/3du3djs9mAgh1PRXX48GEMewuOXK+uXbvmWcfVA9c2adKE7du3s2PHDsaNG0eXLl0cJ+I2m41169bRrVs3+vfvj2EYjnqHDx927AudOnXKd/nZx9OZM2c4fPhwCW5J8dm9e7djWwuybb/88ku+iUGLFi3ynH7hPnvhMeXt7c3NN98M2HuNffHFF9m2bRsZGRnObYSTivPzLM5xJ/M7ngzDoHHjxk4tK7/Pws3NzdEVfHF3uvD77787ltmvX798zwmqVavm6EzoyJEjTq2jcuXKzJ8/n8OHD/Pmm2/Ss2dPrrrqKsd3y969e3nyySdp2bJlnucb2Qn27Nmzc+xnhmHwySefAPDII4/kuIBZkvtpcb5n2RdXL3bh+HXZF6PyK3Px/6TmzZsTGhrKX3/9RYsWLXjnnXfYt29fju/Hik7JmJRpx48fB8BqtV7yit7p06cddS68kn+hSw1Emj3v1KlTjmlffPFFvl9amzdvdpRr3bq14ypgdmL1+++/c/LkSerUqZNjvZdKxi5cTmlzd3fntttu44033mDDhg0kJiaye/duHn/8cdzc3MjKymLo0KE57pBlv1ehoaF5XgXLln3X7cL3tqzLjvVyg9debtuyk9q8ZF+dvPBuakFt2bKF6OhoXnrpJbZu3Up8fDze3t6EhYVRpUoVQkNDHWXzu9JYkNiAHCfFF27npd4bV4y5lH21Nr+LMaWtadOmvPjii3z77bf89ddfHD58mMmTJzs+m5kzZzJ9+nRH+cK8t+XlmHJ227KysvI9ic9vv81vnwX4+OOPadSoEX///Tcvv/wyN954I/7+/rRp04ZJkyaVSC9txfl5hoWFFV9gxagg3yHO3m27nOxzArCP23ep84LsC0f5nRNcTs2aNXnsscf44osvOHToEGfOnGHevHlcc801gL0X4f/85z+56vXs2ZOQkBBOnz7NkiVLHNN/+OEHDh06hMVi4ZFHHslVr6T20+J8zwpy/F2uzMX/84KCgvj888+pXLkye/bs4f/+7/+oX78+wcHB3H333cydO7fY96OyRsmYlGnZTayio6MveUXvwldxdfmblpaW7xfWhVervL29HVcIs5sfZidYF98xioqKokaNGsTFxXHo0CEyMjIciV1pN1G8FLPZTJMmTZg2bRoffvghYP8sZsyYkatsQe9EuPqORWGUxW3LysriwQcfJDExkcaNG/O///2Pc+fOkZSUxMmTJzlx4gRbt251lL9Sri42bNgQ+Ldr6rImKiqKp59+mnXr1jmadn388cd5li2L+11pK85tq1GjBrt372bFihU89thj3HDDDdhsNjZt2sSoUaOoU6cOP/zwQ7Gt72JF/TzzagJ+pbqw2fWJEycKdE5wqaFLnBEcHMxDDz3Etm3bqF+/PgBLly7NlSR5eXk51pn9/xP+bbZ42223ERkZmWv5JbWfuvI9K6hOnTpx+PBhZs+ezcMPP8zVV1/N2bNnWbZsGX369KFJkyaXHLqivFMyJmVadvOeP//8s8htiS91IGfPu/AKZL9+/fL9orp4/KvsRGrTpk1kZWXlGF/sYhfeHdu2bZtj8Ne8krHsK0mXGrPk7Nmz+c4rDn379nWcPO7fv98xPfu9+vvvvy85llhcXByAo9lKQbh6u7O3LTY29pLlCrNtRbVlyxaOHDmCxWJh+fLl3HbbbbmuRJ44caJE1n3h8VGQ46k0ZTfx+fvvvx3Nm8uiBg0a0KZNGyDv4wkuvd9l73NQ8P3uwqvW+R1TpXE8Qc74L5Y9z83NjeDg4GKNwWw2c+utt/Lmm2+yc+dO4uPjmTdvHjVq1CAhIYGHHnqoWJsuluTneSW78Nm74n4eraB8fHzo3bs3YG96fOEzitmGDBniGJfr0KFDnD59mqVLlwLkeTctW0nsp2XhPSsIX19f+vTpw8yZMzlw4ABxcXG8/vrreHl5Oe6YVVRKxqRMy24rn5GR4fgiK6wLO83Ib17Tpk0LtezsRCo5OZkdO3awbt06IPedsQunrVmzxrFePz+/PNtZZ5+QXOqf+bZt2woVc0FZLBa8vLwAcjRHzH6vsrKyHNubl1WrVgH5tyPPS1G3O/tB6MLeFcretri4OA4cOJBnGavV6vj8nNm2osp+TypXrpxv86fs97y4XX/99Y739lLHU0neZchP//79HR1DjB07tsCffXaznNKU/RzhhcdTrVq1HJ2PrF69Ot+62Z9tpUqVqFWrVoHWd2Fik98xVZDjCQp3TF243xRk2xo1aoS7u7vT63GGv78/Dz30kOMZnpMnTxbriWpJfp4VVUH2s2uuuYaAgADAtZ3YXPgscF7N9K+++mpuuukmDMPgo48+cjw/FhERwW233Vbg9VxuPy1P75mzwsPDGTVqlKPztpUrV7o4opKjZEzKtKZNm9KkSRMAnnvuuRwjvOflUm2qJ0+enOdV4TVr1rBp0ybA3plFYbRs2dKRsLz77rucOnWK2rVr59kU4cI7Y9l30Nq0aZNnD0eNGjUC7G2+L2x6lu3UqVN59thUENm9RF3OsmXLHA/cXn/99Y7p1113naOjiPHjx+fZa9///vc/x0negw8+WODYirrd2f94EhMTC7zOC3Xu3NnxDFJ+vSl+8MEHjrb4zmxbUQUGBgL/9op2sbi4ON56660SWXdQUBC33HILkP/xtGrVqhzPVJaW0NBQR++Qq1ev5umnn75s4rBp0yYef/zxYovhhx9+uOyzDceOHXOcgF94PJlMJsf3zwcffJDn3c3jx4/zwQcfAM7tc3Xr1nXc3V68eHGu+TabjQkTJuRbP/t4gsIdU0FBQY6e1CZNmpTn8yg///yzI7biPJ4udxfhwt4Ai7M5YEl+nhVVQfYzNzc3BgwYAMCsWbMuexfc2eestm/fftk6WVlZzJs3D7DfzalXr16e5bI78pg5c6ajueLFHXdkK+x+Whbes6K6VMsa+HfbK3Rz3WLoHl/E4cLxTP7+++9Lvi4cc+dygz5nj0pfq1YtY+HChTnGqYmLizPmzJljdOrUyRg4cGCOuhcP+nzTTTcZ+/btMwzDPkDlwoULHeNJXX/99U6NR3Kxjh07GoBjYMP+/fvnW7Z69eo5yl48YHQ2q9Vq1KxZ0wCMevXqGTt27DBsNpthtVqNNWvWGPXr13fE7+zhnJSUZABGkyZNjDfeeMP4+eefHWOFWa1WIyYmxhg7dqxjDLSAgIAcgzEaRs5Bn7t27eoY9DkjI8OYO3euYwwdZwd9Lup2Z48dZ7FYjE2bNuX7HhR00Of//Oc/xokTJwzDsA+M/NZbbznGirnUoM8lMT5cYmKi4evra4B9YNHswV+zsrKMFStWGFdddZVRqVKlfMfyudRYUNkuNRbQjh07HPvtxcdT9iDqrhj02TAMw2azGffff7+jfsuWLY0lS5YYZ8+edZQ5d+6csWzZMqNbt26GyWTKNTZhUcYZu+GGG4yIiAhj9OjRxoYNG3KMwXTmzBnjo48+cgy8ChjLly/PUT82Ntbx3jVs2DDHvrtx40ajfv36Bjg/6LNhGEbv3r0NsA/w/sUXXzgGpN63b59x9913G8HBwfnWT0lJMTw8PAzAmDhxomGz2fLc/kvtW7t373YcM23atDF++eUXwzDsx/o333zj+E683KDPl5JX/GvWrDGuvfZaY8qUKcbvv//uGKPMZrMZmzZtMq699loDMCIiIvIcK/FSLjfOWFE+z0sdg84qqUGfLxVX9vf3jBkzCly/oPvZmTNnjKuuusoADF9fX+ONN94wTp065ZifmJhofPvtt0bfvn2NBg0aFHibDcMwnn76acPPz8/o37+/sWzZshxjaqWkpBj/+9//jDZt2ji2YcSIEfkuKzMz0zEAdvb/o/zGKSvsflpa79mlvlsMo2D7RX7H8bhx44wuXboYs2fPzjEuXXp6uvHFF18YgYGBBuQ/5mdFoGRMitWFX/qXezVq1MhR71InxYZhGN9//32OE0yLxWJUqlTJkShkvy6VjH355ZeOk4HAwEBHggcYNWrUcCQShfXSSy/liGXWrFn5ln3wwQdzlN2+fXu+ZVesWOGIGzB8fHwMLy8vAzCuvvpq4/PPPy/UyWNKSorjpPrC9zUkJCTH+gAjLCzMWL9+fZ7LmTJlimEymRxlg4KCHP8cAOPaa6/NlcQZxuVPsIqy3ZmZmUa9evUc84ODg42aNWsaNWvWzDH49OX2uyeffNKxDJPJZAQHBzsGpgX7ALTnzp3LVa+kB+t+7733cnw+fn5+jvcmNDQ0x6CrxZ2MGYZhfPDBBzk+8wuPp+joaGPKlCkuScYMw37yMm7cOMPb2zvHe+Tv759jQObsk+DZs2fnqF+UZOzGG2/MsXyTyWQEBgbm+p7y8PAw3nzzzTyXsXbtWsfJR/aJU3bynX185XcsXuqEKTY21pHwAIa7u7vjYom/v7+xdu3aS9Z/5JFHchyLNWrUMGrWrGk8/fTTjjKX27fmz5+f47shICDAsd8CRmRkpPH777/nqlfUZOzC997d3d2oVKlSjuM4ICAg3/f0Ui6XjBlG4T/PKzEZM4yC7WeGYRh//vmn0ahRoxyfbVBQkGOfzn7lNXD9pYwZMyZH/ew4LvwMs199+vQxMjIyLrm8559/3lH+zjvvzLdcUfbT0njPLvXdYBhFS8YuPm/09vY2QkJCcvyPqV+/fp4Dp1cUaqYo5ULnzp05dOgQEyZMoE2bNgQGBpKYmIjZbKZBgwY88sgjfP3117z99tv5LuOee+5h8+bN9OjRAy8vLwzDoFatWjz99NP89NNPRW6vf3EHHJcae+vCeQEBATmaK13s1ltvZcOGDdx5550EBwdjtVqJjIxkzJgx7Nq1K8/BRAvCx8eHkydPMnv2bAYOHMgNN9xAYGAgZ8+exWKxEBER4XiQ+ODBg7Rt2zbP5Tz55JPs3LmT3r17ExkZSWpqqmPMtClTprB9+/YcYwAVVFG2283NjdWrVzNw4ECioqJISUnhyJEjHDlyxDGOSkFMmTKFH374gR49elClShWSk5Px9/enY8eOfPrpp6xcufKS3TuXlCFDhvDNN9/QoUMH/Pz8yMrKIjw8nP/7v//j559/5tprry3R9Q8ePJhNmzZx1113ERISwvnz56lZsybPPPMM27dvL/bOF5xhMpl48cUX+fPPP3n11Ve56aabqF69OhkZGWRlZVGzZk26du3Kxx9/TExMDH369Cm2da9Zs4bly5fz1FNP0bZtW6pUqUJaWhqZmZmEhobSqlUrnnvuOfbu3ctjjz2W5zLat2/Pvn37ePrpp6lfvz42mw3DMKhfvz4jRoxg7969+R6LlxIREcG2bdsYOHCg41lDPz8/+vbty+7duy/5fQUwffp0xo4d6+jW++jRoxw5ciTHsCKXc//997Nnzx7+85//cNVVV3H+/Hnc3Nxo3Lgx48aN47fffnP0UldcmjVrxoIFCxg6dCg33HADoaGhnD17Fi8vLxo3bsyoUaMK/Z4WREl9nhVVQfezWrVqsXPnTmbPns2dd95JtWrVSElJISMjg1q1atGtWzc+/fRTtmzZ4tT6X331VbZu3eoYIzAqKgrDMEhOTiYwMJBGjRrxn//8h40bNzJ79uzLPtt43333OX6/VMcdRdlPXf2eFdXgwYP58MMPefDBB7nmmmvw8fHh3LlzBAcH07ZtW6ZNm8bu3bsLfa5THpgM4wrp91hEREREpJS88cYbjBgxgoiICGJiYir2c09SaLozJiIiIiJSjKxWK++99x5gv/ujREzyo2RMRERERKSY2Gw2/vvf//LHH3/g6+vr6FlRJC+5+9IWERERERGnLFq0iBEjRhAfH09SUhIA48aN04DecklKxkREREREiig5OZkjR47g7u5OdHQ0jz76KMOHD3d1WFLGqQMPERERERERF9AzYyIiIiIiIi6gZExERERERMQFlIyJiIiIiIi4gJIxERERERERF1AyJiIiIiIi4gLq2r6cqVq1KikpKdSoUcPVoYiIiIiIXPGOHj2Kr68vJ06ccLqu7oyVMykpKWRmZro6DBERERERATIzM0lJSSlUXd0ZK2ey74jt2bPHxZGIiIiIiEjDhg0LXVd3xkRERERERFxAyZiIiIiIiIgLKBkTERERERFxASVjIiIiIiIiLqBkTERERERExAWUjImIiIiIiLiAura/QhiGgWEYrg5DipnJZMJkMrk6DBEREREpBCVjFZjVauXMmTMkJSWRkZHh6nCkhHh4eODv70+lSpWwWCyuDkdERERECkjJWAVltVo5evQo6enprg5FSlhGRgZnzpwhJSWFGjVqKCETERERKSeUjFVQZ86cIT09HYvFQpUqVfD19cVs1iOCFY3NZiMlJYWTJ0+Snp7OmTNnCAsLc3VYIiIiIlIASsYqqKSkJACqVKlCYGCgi6ORkmI2mx2f7/Hjx0lKSlIyJiIiIlJO6FZJBWQYhuMZMV9fXxdHI6Uh+3POyMhQRy0iIiIi5YSSsQrowpNxNU28Mlz4OSsZExERESkfdKYuIiIiIiLiAnpmTEREREREygWrzWDdgVMs2BHHkfhUUjOy8PFwo2aIDz2bRdC+bhgWc/kZg1XJmIiIiIiIlGlWm8GMTYeZuTmGuIS0XPP3/nWOFXtOEBHsTb9WUQxoXQtzOUjKlIxJkVS0qxMiIiIiUrakZ1p5Yv5PrNhzAoBAb3d6No2gdZ1Q/L3cSErPYuPB0yzcFUdcQhrjv9nLzpgEpj3QGC/3sj3+qpIxKZTyenXCZLLHEBQUxOHDhwkKCspVZuzYsYwbN44JEyYwZsyYYl1/amoq33//PcuWLWPHjh3ExMRgtVqpU6cOPXr04KmnnsLPz69Y1ykiIiJSXlltBo/P/5Hv9pzEw2Lm+Tvr07NpZK4kq0O9MJ6+pR4Ld8UyfvleVuw5wRPzf2J6r+vL9I0BdeAhTkvPtDJ83m7Gf7OXuIQ0Ar3dGdS2FjP7N2Px0JbM7N+MgW1qEejt7rg6MWzebtIzra4O3SExMZGpU6eW+no/++wzunXrxqefforNZqNLly60bduWw4cP89///pdmzZpx6tSpUo9LREREpCyasemwIxGbOaAZfVtG5Xu3y9vDQt+WUcwc0AwPi5kVe04wc3NM6QbsJCVj4pTsqxMr9pzAw2LmpXsasu3Zm3nujgZ0qBfGDTVD6FAvjOfvbMDWZ27mpXsaOg6GJ+b/hNXm+m7XzWYzHh4eTJs2jYSEhFJdt4eHB0OHDuXAgQP89ttvLFiwgBUrVrB//36aNGnCvn37eOKJJ0o1JhEREZGyyGozHMnU83fWp9VVoQWq1+qqUJ67oz4AMzcfLhPnn/m54pOxzZs3c/vttxMSEoKfnx/Nmzdn1qxZTi9n165djB07lrZt21K9enU8PT2JjIykd+/e/PLLLyUQuWtUhKsT7u7uDBw4kHPnzjFlypRSXXffvn159913ufrqq3NMr1atGtOnTwdgyZIljkG7RURERK5U6w6ccrTCuu+GSKfq9mwaSaC3O7Hxaaw/8HcJRVh0V3QytnTpUtq1a8eKFSu47rrr6NKlCwcPHqRfv3489dRTBV5OVlYWTZs2Zdy4cezbt48mTZpw99134+npybx582jatCmLFi0qwS0pHRXp6sSzzz6Lp6cnb775JvHx8a4OB4BGjRoBcP78ec6cOePiaERERERca8GOOAB6No3A28O5jji8PSzcd0MEAF/siC322IrLFZuMJSQk0L9/f6xWK4sWLWLt2rUsWrSIffv2UadOHaZOncqaNWsKvLwWLVqwfPlyTp48yTfffMPChQs5cOAAzz33HJmZmQwYMIDTp0+X4BaVvIp0dSI8PJxBgwaRlJTEG2+84epwAPjzzz8B+527kJAQF0cjIiIi4lpH4lMBaF2nYDcALtbmanu9o/8spyy6YpOxjz/+mLNnz3LPPffQvXt3x/QqVaowceJEgAI3YXNzc2Pr1q3ccccdmM3/vqVms5mXX36Z6OhokpKS+Oabb4p3I0pZRbs68cwzz+Dl5cXbb79doDtRa9euxWQyOfXq0KFDgeN58803AejSpQuenp6F3SwRERGRCiE1IwsAf6/CdQDv52mvl/LPcsqiK7Zr++XLlwNw77335pp3xx134OXlxapVq0hPT8fLy6vQ6zGZTFx77bXs27eP48ePF3o5ZUFxXJ34eOPhMnN1onr16gwePJi33nqLyZMnM2HChEuWr1q1Kg8//LBT64iOji5Quf/973988sknuLu78/LLLzu1DhEREZGKyMfDnqokpRcumUo+b6/n61F2U56yG1kJy+5U4/rrr881z8PDg2uuuYadO3eyf/9+x7M8hZXd/Kxq1apFWo6rVcSrE2PGjOHDDz/knXfe4emnnyY0NP9EMzo6mpkzZxZ7DHv37qV3794YhsGkSZOKvL+JiIiIVAQ1Q3zY+9c5Nh06TYd6YU7X33jQ/ohQjRCf4g6t2FyRzRTPnTtHYmIiABEREXmWyZ5+9OjRIq1r48aN7Nq1Cw8PD7p06VLgeg0bNszz9ccffxQpnqKoiFcnqlWrxpAhQ0hOTmbSpEmlvv64uDi6dOlCQkICTz31FI8//nipxyAiIiJSFvVsZj8fX7AzjrQM58arTcuwsmCn/dGY+5s519dBaSo7Z8WlKDk52fG7j0/embKvr2+uss46d+4cAwYMAODJJ5+kWrVqhV5WWVBRr06MHj2aDz74gOnTpzNixIh8y+3bt4/XXnvNqWVHR0czZsyYPOedPn2azp07c/ToUfr378/kyZOdWraIiIhIRda+bhgRwd7EJaSxcFcsfVtGFbjugp2xnEvPIjLEm3Z1K5dckEVUbpOxe++9l99++82pOrNnz6Z58+YYxuW7Vi9ImUuxWq089NBDHDx4kObNm/PSSy85VX/Pnj15Tm/YsGGR4iqKns0iWLHnBAt2xvFU53pOdeJRlq9OVK1alaFDhzJlyhQmTpzoSMQvduLECafHoGvfvn2eyVhSUhK33XYb+/bto3v37nz00UeYTKZCxS8iIiJSEVnMJvq1imL8N3sZv3wvdcL8CjS00uZDp3nlm70A9GtVC4u57J5jldtkLCYmhv379ztVJzXV3nGEv79/jmkBAQH5lvXz8ytUfIMHD+abb76hXr16fPPNN3h4eBRqOWVJRb46MXr0aN5//33effdd+vbtm2eZDh06FDlJB/s4Yvfccw87d+7k1ltv5fPPP8dica53ShEREZErQf/WtdgZk8CKPSfo9+kOnr+zPvfdEJnnTYHsi/+vfLOXDKuN266pSv9WUaUftBPK7TNjO3fuxDAMp17Z3YwHBAQQGBgI2J/ZyUv29Bo1ajgd28iRI/n000+JjIxk5cqVl+wUojzJvjoBMH75Xjb/UbBx08rD1YmwsDCGDRtGamqq03e/nGG1WnnwwQdZs2YNbdu2ZcmSJRUiURcREREpCRaziWkPNKZLw6pkWG28+NUebpywmvHLf2ft/lPsjIln7f5TvLz8d1q8uor/fr3HkYhNvb8x5jJ43nmhcntnrKgaNWrE+vXr2b17Nw0aNMgxLzMzk99++w1PT0/q1avn1HInTJjA5MmTCQsLY+XKlURGlq0meUVVka9OjBo1ivfee4+UlJQSW8c777zD0qVLAQgNDWXYsGF5lps8eXKFSeJFREREisLL3cL0Xtczc3MMMzYdJi4hjY83HubjjYdzlY0M8aZfq1r0bxVV5hMxuIKTsTvuuIP169ezaNEievfunWPe8uXLSU9P5/bbb3dqjLEPP/yQZ599lqCgIL777junE7nyIPvqxBPzf2LFnhO8+NUe3vj+APfdEEGbq0Px83Qj+XwWGw6eZuE/TROBcnF1onLlygwfPtwx6HdJSEhIcPyenZTlZezYsUrGRERERP5hMZt4pE0t+rWKYv2Bv/liRyxH41NJycjC18ONGiE+3N8sknZ1K5fJVlj5MRnF8RBMORQfH0+tWrU4d+4cixcvpnv37gCcOnWK1q1bc+jQIVatWsXNN9+co172IL6rV68mPDzcMX3RokXcf//9+Pj48P3339OyZcsSiTu7A4/8OvgAsNlsjufp6tWrh9lc/K1RrTYjx9WJ/JS3qxPlVWl85iIiIiKSW0HOz/NzxSZjAIsXL6Znz54YhkH79u0JDQ1l1apVJCYm8thjj/Hmm2/mqpPd493hw4eJiooC7AlcZGQkGRkZXHvttXkOJA3QtWtXunbtWqSYy0oyls1qMyrU1YnySsmYiIiIiGsUJRm7YpspAvTo0YP169czfvx4tm7dSkZGBvXr12f48OH079+/wMtJTU0lIyMDgF9//ZVff/01z3JRUVFFTsbKGovZRMfoMDpGOz/umIiIiMjFrDaDdQdOsWBHHEfiU0nNyMLHw42aIT70bBZB+7phutArFcYVnYwBtG7dmm+//bbA5fO6kRgVFVUsXZ6LiIiIXKmsNoMZmw4zc3NMno9A7P3rHCv2nCAi2Jt+raIY0LqWHoGQcu+KT8ZERERExLXSM62OzsEAAr3d6dk0gtZ1QvH3ciMpPYuNB0+zcFcccQlpjP9mLztjEpj2QGO83DVWp5RfSsZERERExGWsNoPH5//Id3tO4mEx8/yd9enZNDJXktWhXhhP31KPhbtiGb98Lyv2nOCJ+T8xvdf1arYo5Zae8hcRERERl5mx6bAjEZs5oBl9W0ble7fL28NC35ZRzBzQDA+LmRV7TjBzc0zpBixSjJSMiYiIiIhLZA+VA/D8nfVpdVXBxthsdVUoz91RH4CZmw9jtenZfSmflIyJiIiIiEusO3CKuIQ0Ar3due+GSKfq9mwaSaC3O7Hxaaw/8HcJRShSspSMiYiIiIhLLNgRB0DPphF4ezjXEYe3h4X7bogA4IsdscUem0hpUDImIiIiIi5xJD4VgNZ1CtY88WJtrrbXO/rPckTKGyVjIiIiIuISqRlZAPh7Fa6Dbz9Pe72Uf5YjUt4oGRMRERERl/DxsCdTSemFS6aSz9vr+XpotCYpn5SMiYiIiIhL1AzxAWDTodOFqr/xoL1ejX+WI1Le6DKCFI3NCodWwY9zID4GMpLBww9CoqBJH6jTCczOPZArIiIiV4aezSJYsecEC3bG8VTnek514pGWYWXBTnvHHfc3c64nRpGyQnfGpHBsVtgyHd5qDJ/1hL3L4OSvkHDY/nPvMvv0txrby9lsro4YAJPJhMlkIjg4mMTExDzLjB07FpPJxGuvvVYqMcXHxxMWFobJZCI6OrpU1ikiIlIWtK8bRkSwN2fTMlm4y7keERfsjOVcehaRId60q1u5hCIUKVlKxsR5mWmw8GH47llIPApeQdDyUei1GAZ8b//Z8lH79MSj9nIL+9rrlRGJiYlMnTrV1WEA8NRTT3H6dOGaZ4iIiJRnFrOJfq2iABi/fC+b/yjY/8PNh07zyjd7AejXqhYWs6mkQhQpUUrGxDk2KyweaL/zZfGA2yfD0/vg1lfg6k5Qo4X9562vwFN77fMtHvbySwbZ67uY2WzGw8ODadOmkZCQ4NJYVq9ezaxZsxg0aJBL4xAREXGV/q1r0aVhVTKsNvp9uoPZW2JIy8j7fCEtw8qszTH0m7GDDKuN266pSv9/kjmR8kjJmDhn2/uwb7k9weq9GJoPAnfvvMt6+Njn9178b0K27YPSjTcP7u7uDBw4kHPnzjFlyhSXxZGWlsaQIUNo0KABI0aMcFkcIiIirmQxm5j2QGNHQvbiV3u4ccJqxi//nbX7T7EzJp61+0/x8vLfafHqKv779R5HIjb1/saYdVdMyjElY1JwNqs9GQO49VWo1a5g9Wq1g1tesf++7f0ycXfs2WefxdPTkzfffJP4+HiXxDBu3Dj++OMP3nvvPdzd3V0Sg4iISFng5W5heq/reeHOBo5nyD7eeJh+M3Zw7/tb6DdjB59sPOx4RuyFOxsw/aHr8XJXJ2FSvikZk4I7tOrfZ8Qa93KubpPe/zxDdgQOrS6J6JwSHh7OoEGDSEpK4o033ij19f/yyy+88cYb9O/fn3btCpjUioiIVGAWs4lH2tRi3ciOzOjXjC4Nq9KgWgA1K/nQoFoAXRpWZUa/Zqwd0ZFH2tTSHTGpEJSMScH9OMf+s0lvexNEZ3j42OsB/Di7eOMqpGeeeQYvLy/efvttzpw5c9nya9eudfTGWNBXhw4dci3HZrMxaNAggoKCmDhxYglsmYiISPllMZvoGB3G+31u4H+Pt2XdyI787/G2vN/nBjpGh6mzDqlQNM6YFFx8jP1n7Y6Fq1+7I2x5BxJiiiuiIqlevTqDBw/mrbfeYvLkyUyYMOGS5atWrcrDDz/s1Dry6qr+7bffZvv27cyYMYNKlSo5tTwRERERqTiUjEnBZSTbf3r6F66+p5/95/nk4omnGIwZM4YPP/yQd955h6effprQ0NB8y0ZHRzNz5swirS82Npbnn3+e9u3b069fvyItS0RERETKNzVTlILzyE6mkgpXPzsJy07KyoBq1aoxZMgQkpOTmTRpUomvb9iwYWRkZPDee++V+LpEREREpGzTnTEpuJAoOPkr/LnGPpaYs/5cY/8ZHFWcURXZ6NGj+eCDD5g+ffolu5jft28fr732mlPLjo6OZsyYMY6/ly9fTlBQEEOHDs1RLj09HYCjR486njNbvnw5fn5lJ3EVERERcTXDaiV5wwbOLlyE+47t2M6fJ61mFB41axLYozt+bdtispSfXjaVjEnBNeljHyvsx7nQ8TnnOvHISL2gA5C+JRNfIVWtWpWhQ4cyZcoUJk6ciK+vb57lTpw4waxZs5xadvv27XMkYwCJiYmsW7cuz/JpaWmOeVlZWU6tS0RERKSiMqxW4ufMIWHOXDKPHcMnJYWIuFgAjsQnkHTgAEkrV+IeHk5wn96E9O2LyVz2GwGW/Qil7KjTCYJqQHoi/DTPubo/zoX0sxBUE+rcXCLhFcXo0aPx8fHh3Xff5eTJk3mW6dChA4ZhOPVau3ZtjmXkV+7w4cMA1KtXzzEtKCiohLdaREREpOyzpadz7PEn+PvVCWTFxmLx9iTMO9Mxv3pwGiFNK2Hx9iQrNpZTE17j2ONPYPun5VFZpmRMCs5sgRZD7L9/9ywcXl+wen+ug++fs//eYoh9OWVMWFgYw4YNIzU11em7XyIiIiJSMgyrleMjR5K17GuiDx4g+uAB6v70M94HTzjKuP96mirzNlH3p5+JPngA78zzJK1cyfGRozCsVhdGf3lKxsQ5LYZA/bvAmgFze8D2j+xNEPOSkQrbPoR599rL17/732SuDBo1ahS+vr6kpaW5OhQRERERAeLnzCFp5SoCk84VuI4t4DQmi4mklStJmDu3BKMrOiVj4hyzBbp/9G9C9r8RMKU+fPccHFwFR7faf654FqZEw7cj/03Eun8IZbjtbuXKlRk+fLirwxARERER7HfFEubYk6lvbwsktrL7ZevEVnZncN/qbGxjb8YYP2dumb47ZjIMw3B1EFJwDRs2BGDPnj35lrHZbOzfvx+wP4NkLokEyGaFbR/Atvcg8Wj+5YJq2u+GtRhSphOx8q5UPnMRERGRUpS0di1xQ4aS7mkw8P/cqJJk8MXU0/jHnsqzvK12LZZ9/BTjYj/BlJHJx29n4XXeROQH7+PXvn2JxVmQ8/P8qDdFKRyzBVoOgxb/gUOr4cfZkBBjH0vM08/efX2TvvbOOsrgM2IiIiIiUradXbwYgJWNzBgeFsY+8BH+NeOgR488y5snTeaejt2p9ldzhqwcxMpGNu7abpC4aFGJJmNFoWRMisZsgbq32F8iIiIFZLUZrDtwigU74jgSn0pqRhY+Hm7UDPGhZ7MI2tcNw2I2uTpMEXGhjEO/A/BLlImR1z9F82rNYdG0/CssXgzdu9O8WnNGXP8U/zs0kbu2G2Qc2ls6AReCkjEREREpNVabwYxNh5m5OYa4hNwdJu396xwr9pwgItibfq2iGNC6FmYlZSJXpNTEv7EAZi83ukb3tE/8/PN/CwwbBoYB771n//uzz2CeffilbtE9+c5zCpBBamLezRrLAiVjIiIiUirSM608Mf8nVuyxd0kd6O1Oz6YRtK4Tir+XG0npWWw8eJqFu+KIS0hj/Dd72RmTwLQHGuPlribvIleaBLcsQoGbfRvg7eZtnxgTA336wAsvQOfO9mldu8Irr8CcOY663m7edPRtAPxEgntWKUdecErGREREpMRZbQaPz/+R7/acxMNi5vk769OzaWSuJKtDvTCevqUeC3fFMn75XlbsOcET839ieq/r1WxR5ApzMhhC/4Zrj16QstSsCesvGuv2llvsr4tcd9T+/XIysCSjLBp1uSYiIiIlbsamw45EbOaAZvRtGZXv3S5vDwt9W0Yxc0AzPCxmVuw5wczNMaUbsIi43Jbr7N8RXmt/xebkOLC2tDS81v4GwOZGZffOupIxERERKVFWm+FIpp6/sz6trgotUL1WV4Xy3B31AZi5+TBWm0bjEbmSHK3jyalAMKWcJ3HJEqfqJi5egin1PCeDIK6OZ8kEWAyUjImIiEiJWnfgFHEJaQR6u3PfDZFO1e3ZNJJAb3di49NYf+DvEopQRMqiiEp1+bapPV059drrpGzdVqB6KVu3cur11wH49gYz4ZXqlliMRaVkTERERErUgh1xAPRsGoG3h3PNhbw9LNx3QwQAX+yILfbYRKTs6tZoEP9ramJXXTAyM4kdNIj4efPybbJoS0sjfu48YgcNxsjMZGdd+LaZiW6NB5dy5AWnDjxERESkRB2JTwWgdZ2CNU+8WJurQ/l442GO/rMcEbkytI5oR3WPAKbcfY7pi88TdBhOvjyev996m6Bu3fBt3Qqzry+2lBRSNm4icelSbOfOAZBQK4upd3tS3SOA1uFtXbwl+VMyJiIiIiUqNcPerbS/V+FOO/w87fVSMspu99QiUvwsZgsPNRrMpF1v8Oh9Hny8LhX/P8LIPH2W+JkziZ85M1cd98oBJNU+xf918CHTbKJXoyFYzGW3Aw8lYyIiIlKifDzspxtJ6YVLppLP2+v5eui0ReRK06tBH348tZtVsWt4pKMvI6+J45YjHqSdqk1GshlbegZmLw88/Gx4h/3JdzXPMblyMJkmE50jb6JXg96u3oRL0reaFInVZmXT8U0sPbiU2KRYUrNS8XHzIdI/km5Xd6N19dZl+mqEiIiUvJohPuz96xybDp2mQ70wp+tvPHgagBohPsUdmoiUcRazhQntJsKGMaw6uppXw0J4p5KVrsmHaJmWjq/NRorZzGZvL7708yPJ4gdA5xqdeLXtBMymst1FRtmOTsosq83K7D2zuWPpHQxfPZxVR1exP2E/sUmx7E/Yz6qjqxi+ejh3LL2D2XtmYzNsrg4ZAJPJhMlkIjg4mMTExDzLjB07FpPJxGuvvVaisZw4cYInn3ySunXr4u3tTUhICDfccAOjRo0q0fWKiJS2ns3sHXAs2BlHWobVqbppGVYW7LR33HF/M+d6YhSRisHLzYvJ7d9gVLNRhPtW55zFwuzAAIZWDaNv9aoMrRrGnMAAkixmwv3CGdVsFJM7vIGXm5erQ78s3RkTp6VnpfPMhmdYdXQVAAEeAXSr040bq9+In7sfyZnJbDm+hS8Pfcmx5GNM2jmJH0/9yIS2E8rMQZGYmMjUqVMZN26cS9a/ZcsWbr/9dhITE2nQoAF33303SUlJ/P7770yZMoWJEye6JC4RkZLQvm4YEcHexCWksXBXLH1bRhW47oKdsZxLzyIyxJt2dSuXXJAiUqZZzBb6NOjDQ9EPOVplxSXHkZKZgq+7LxF+EeWyVZaSMXGK1WZlzIYxrD66GnezOyObjaRbnW65kqw24W14tMmjfHnoSybtmGRP3DbA5PaTXX6AmM1m3NzcmDZtGk888QTBwcGluv7jx49z++23c/78eZYsWUK3bt1yzN++fXupxiMiUtIsZhP9WkUx/pu9jF++lzphfgUa+HnzodO88s1eAPq1qoXFbCrpUEWkjLOYLbSLaEe7iHauDqVYqJmiOGXe3nmOROz9Tu/zYPSD+d7t8nbz5sHoB3m/0/u4m91ZdXQVn+37rJQjzs3d3Z2BAwdy7tw5pkyZUurrHzNmDImJiUycODFXIgbQvHnzUo9JRKSk9W9diy4Nq5JhtdHv0x3M3hKTb5PFtAwrszbH0G/GDjKsNm67pir9W0WVbsAiIqVAyZgUmNVmdSRTI5uNpHm1giUNzas1Z0TTEYA9mbPanHteoCQ8++yzeHp68uabbxIfH19q601ISGDBggUEBgYycODAUluviIirWcwmpj3Q2JGQvfjVHm6csJrxy39n7f5T7IyJZ+3+U7y8/HdavLqK/369x5GITb2/MWbdFRORCkjNFKXANh3fxLHkYwR4BNC1Tlen6na7uhvTf5rOseRjbDq+yeW3lsPDwxk0aBDvvPMOb7zxBq+88kqprHfTpk2cP3+eTp064e7uzqJFi9i4cSOZmZlER0fTs2dPqlSpUiqxiIiUNi93C9N7Xc/MzTHM2HSYuIQ0Pt54mI83Hs5VNjLEm36tatG/VZQSMRGpsHRnTAps6cGlAHSr0w1vN2+n6nq7eTsSuOzluNozzzyDl5cXb7/9NmfOnLls+bVr1zp6Yyzoq0OHDjmWsWfPHgCqVKlC27Ztue+++3jzzTd59913eeyxx7jqqqtYuHBhSWyuiEiZYDGbeKRNLdaN7MiMfs3o0rAqDaoFULOSDw2qBdClYVVm9GvG2hEdeaRNLSViIlKh6c6YFFhskr1r4Rur31io+i2rt2T277OJS44rzrAKrXr16gwePJi33nqLyZMnM2HChEuWr1q1Kg8//LBT64iOjs7xd0JCAgCzZ8/G09OTTz75hLvvvpvk5GTefvttpkyZQu/evalXrx7XXXedcxskIlKOWMwmOkaH0THa+XHHREQqCiVjUmCpWakA+Ln7Faq+r7svACmZKcUWU1GNGTOGDz/8kHfeeYenn36a0ND8e/eKjo5m5syZRVqf1Wp/Xi4rK4vp06czYMAAAEJDQ3njjTc4evQoixYtYuLEicydO7dI6xIRERGRsk3NFKXAfNx8AEjOTC5U/ewkLDspKwuqVavGkCFDSE5OZtKkSSW+Pn9/f8DevX5ed9myk7O1a9eWeCwiIiJlks0KB76DL3rDe23gzcb2n1/0tk8vAx2BiRQX3RmTAov0j2R/wn62Ht9Km/A2TtffcnwLABF+EcUdWpGMHj2aDz74gOnTpzNixIh8y+3bt4/XXnvNqWVHR0czZswYx99RUVGAvcmjp6dnrvLZ80+dOuXUekRERMo9mxW2vY+x5T2S953i7CFv3GPPY7OZSAtNwsM/jsDa3+MXHYap5VBoMRTMuq8g5ZuSMSmwbld3Y9XRVSw9tJThTYY71YlHWlYaSw8tdSynLKlatSpDhw5lypQpTJw4EV/fvO/cnThxglmzZjm17Pbt2+dIxpo0aQLYnx0zDAOTKeeD6dkdifj5Fa4pqIiISLmUmYaxaCDxX68h4YAvmSmV8ElJISLuJABHMiJJ8vUlKc4b913nCd7xCiF3b8Z078fg7lynYiJliS4nSIG1rt6acL9wzmWc48tDXzpVd+nBpSRlJBHuF07r6q1LJsAiGD16ND4+Prz77rucPHkyzzIdOnTAMAynXhc3N7z22mupVasWaWlpbNu2Ldc6sstff/31xb2JIiIiZZPNim3+AI59vJG/d/mTlWTG4utLWNWqjiLVG11HyEMPYfH1JSvJzKndARz7ZBO2Lx5Rs0Up15SMSYFZzBYein4IgEk7JrH9r+0Fqrftr21M3jkZgF71e2ExW0osxsIKCwtj2LBhpKamOn33y1mjR48G4LHHHuP06dOO6bt27eKNN94AYMiQISUag4iISFlhbHmX47O2kHUIog8eIPrgAeru3oX3zh2OMu7ffEOV/75I3d27iD54AO/MTJLivDk+czPGlvdcGL1I0SgZE6f0qt+LTjU6kWnLZMiqIXy+73PSstLyLJuWlcZnez9j6KqhZNoy6VyzM73q9yrliAtu1KhR+Pr6kpaW9/YUl0GDBnHfffexY8cO6tWrx1133UXHjh1p1aoViYmJDBo0iHvvvbdEYxARESkTbFbiP36fpDhvApOSClwtvHEjTG4WkuK8SfjkPd0dk3JLyZg4xWK2MKHtBEdC9uq2V+m0sBOTdkxi47GN/HjqRzYe28jEHRO5eeHNTNg+wZGIvdrmVcymsrvLVa5cmeHDh5f4esxmM/Pnz2f69OnUrFmTH374gR07dtC0aVNmz57Nhx9+WOIxiIiIlAXG/u9J+CXD/sdzz5AUUfmydY5X8earYW2pPOIpAOJ/zsDYv7IkwxQpMSbDMAxXByEF17BhQwD27NmTbxmbzcb+/fsBqFevHuYS6GnIarPy2b7PmLd3HseSj+VbLtwvnF71e9Grfq8ynYiVd6XxmYuIiBS3pJfvIm7eIcw+7nw0oQM/7/mOT187TOTfmXmWPxbmycOja3Kykge3Vu3IoGfWYkvNJLJXHfxeWFbK0YvYFeT8PD/qTVEKxWK20KdBHx6KfohNxzex9OBS4pLjSMlMwdfdlwi/CLpd3Y3W1VuXyWfERERExPXObj0KwI9Ng/nuxBrcK/sS8/xwIh+flmf50Hdm8Mi19mfXvzuxhmZNQ7hu/SkStx1F/RBLeaRkTIrEYrbQLqId7SLauToUERERKWcyEu13wJZXPY272ZP3O71P86en5Vve88vlPHjfPK4KvIohq4awrOpprrtgOSLljdoyiYiIiIhL2LLs422meZoY2Wwkzas1h88//7fAsGEwdOi/f3/2GQDNqzVnRNMRpHva69syc47bKVJeKBkTEREREZdI87I30qqUYaFrna72iTEx0LYtfP89TJ8O774L330H7drBkSOOut2u7kZIhv1RiFRv91KOXKR4qJmiiIiIiLhEbGVvap3M4J4D6Xjb/ulTrmZNWL8+Z8FbbrG/LuBtM7h7fzpgJi7Ui+tKJ2SRYqU7YyIiIiLiEhubBQAQtQds22Y6Vde2dQa1fs9eTmAxRyZSOpSMiYiIiIhL/HyVhVOBYM4wk/jBBDi8/vKVAP5cR+KHr2HOMHMyCH6+Sqe0Uj5pzxURERERl/D28OXbpvbT0VO7fUmZ/ABs/wgyUvOukJEK2z4k5Y0HObXbF4BvbzDj46mO7aV80jNjIiIiIuISkf6R/K/pPm45W4NqO48SuyaAsMQXCao/HnOz3lC7I3j6wflk+OMHbDvnkvh7Jqd+CsSwmfiraQ2+bXaMm/0iXL0pIoWiZExEREREXKLb1d1YdXQVL92SxKzAjqSuXsPJXUH8/auNoB2z8a36EWY3A1uWiZQTniT+6YMt035HzKfTTYxrvhvDMNHt6m4u3hKRwlEyJiIiIiIu0bp6a8L9wjmWfIzNw1pza/MWxM+eQ+axY8Tv9yN+f+7mh+4REYT06c2Kpm7E71xPuF84rau3dkH0IkWnZExEREREXMJitvBQ9ENM2jmJSbvf4Kpb3qdZ7+9I2biRxEWLyYiLw5aSgtnXF4+ICILu7YFvmzZsP7WTyavsg0H3qt8Li9ni4i0RKRwlY1IkhtVK8oYNnF28hIzYWGypqZh9fPCIjCSwR3f82rbFZNEXpIiIiOStV/1e/HjqR1YdXcWQVUMY2WwkXVt3JaJ9+1xl07LS+PzAF0zeOZlMWyada3amV/1eLohapHioN0UpFMNq5czMmfxxy63EDRlK0sqVnN+3j8yjRzm/bx9JK1cSN2Qof9xyK2dmzsSw2VwdMgAmkwmTyURwcDCJiYl5lhk7diwmk4nXXnutxOLYunUrPXr0oGrVqri7uxMSEsLNN9/MokWLSmydIiIiZZHFbGFC2wl0qtGJTFsmr257lU4LOzFpxyQ2HtvIj6d+ZOOxjUzcMZGbF97MhO0THInYq21exWzS6ayUX7ozJk6zpadzfOQoklauBMAcGEhQ9+74tmqJ2dcPW0oyKZs2k7h0KZnHjnHqtddJ27Wb6pMmYvbycnH0domJiUydOpVx48aV+roXLlzIAw88gM1mo2nTpnTo0IHjx4+zdu1afvjhB0aPHl2iiaCIiEhZ4+XmxeT2k/ls32fM2zuPY8nHmP37bGb/PjtX2XC/cHrV70Wv+r2UiEm5p2RMnGJYrRwfOZKklaswubsTNmY0QT165Eqy/Nq2pfLjj5G4ZAmnXnudpJUrOT5yFOHTprq82aLZbMbNzY1p06bxxBNPEBwcXGrrzsrKYvjw4dhsNubPn8/999/vmLdlyxZuuukmJk6cyKBBg7jqqqtKLS4RERFXs5gt9GnQh4eiH2LT8U0sPbiUuOQ4UjJT8HX3JcIvgm5Xd6N19dZ6RkwqDCVj4pT4OXMciVjkRx/he2OLfMuavb0J6dULz6vqEDtoEEkrV5Iwdy4hDz9cihHn5u7uziOPPMK7777LlClTePnll0tt3fv27ePvv/8mOjo6RyIG0LJlS2699Va++uordu3apWRMRESuSBazhXYR7WgX0c7VoYiUON3blQIzrFYS5swFIGzM6EsmYhfyvbEFYaNHAxA/Zy6G1VpiMRbUs88+i6enJ2+++Sbx8fGltl5PT88ClQsJCSnhSERERETE1ZSMSYElb9hA5rFjjmfEnBHUozvmwEAy4+JI2bixhCIsuPDwcAYNGkRSUhJvvPFGqa23du3a1K5dm3379rFgwYIc87Zs2cJ3331HrVq1aNdOVwNFREREKjolY1JgZxcvASCoe3fM3t5O1TV7exPUrRsAiYsWF3tshfHMM8/g5eXF22+/zZkzZy5bfu3atY7eGAv66tChQ45lWCwWZs6cSWBgIPfffz/NmjXjgQceoH379rRp04bGjRvz/fff4+HhUUJbLSIiIiJlhZ4ZkwLLiI0FwLdVy0LV923diviZM8mIiyvOsAqtevXqDB48mLfeeovJkyczYcKES5avWrUqDzv5vFt0dHSuaW3btmXdunV069aNnTt3snPnTgD8/f3p1KkT1atXd2odImWZ1Waw7sApFuyI40h8KqkZWfh4uFEzxIeezSJoXzcMi9nk6jBFRERcQsmYFJgtNRUAs69foeqbfX3ty0lJKbaYimrMmDF8+OGHvPPOOzz99NOEhobmWzY6OpqZM2cWeZ2ff/45/fv358Ybb2T+/Pk0bNiQ48ePM3nyZMaPH8/q1atZt24d7u7uRV6XiKtYbQYzNh1m5uYY4hLScs3f+9c5Vuw5QUSwN/1aRTGgdS3MSspEROQKo2aKUmBmHx8AbCnJhaqfnYRlJ2VlQbVq1RgyZAjJyclMmjSpxNd38OBBHn74YSpXrsw333xD8+bN8fX15eqrr+aDDz7grrvuYsuWLcyYMaPEYxEpKemZVobP2834b/YSl5BGoLc7g9rWYmb/Ziwe2pKZ/ZsxsE0tAr3diUtIY/w3exk2bzfpma7v3EdERKQ06c6YFJhHZCTn9+0jZfMW/Nq2dbp+yqbN9uVERBR3aEUyevRoPvjgA6ZPn86IESPyLbdv3z6nB2OOjo5mzJgxjr/nz59PZmYmXbp0wTePpLRnz54sW7aMtWvXMnjwYKfWJVIWWG0Gj8//ke/2nMTDYub5O+vTs2kkXu45xwTqUC+Mp2+px8JdsYxfvpcVe07wxPyfmN7rejVbFBGRK4aSMSmwwB7dSVq5ksQlS6j82P851YmHLS2NxCX/dAByb4+SCrFQqlatytChQ5kyZQoTJ07MM0kCOHHiBLNmzXJq2e3bt8+RjMX987xcQEBAnuWzp5dmd/sixWnGpsOORGzmgGa0uir/pr/eHhb6toyiTpgf/T7dwYo9J5i5OYZH2tQqxYhFRERcR80UpcD82rbFPTwc29mzjsSqoBIXL8F27hzuERH4tmlTQhEW3ujRo/Hx8eHdd9/l5MmTeZbp0KEDhmE49Vq7dm2OZVStWhXA0WnHxXbs2AFAVFRUsW2bSGmx2gxmbo4B4Pk7618yEbtQq6tCee6O+gDM3HwYq80oqRBFRETKFCVjUmAmi4XgPr0BOPXa66Rs3Vageilbt3Lq9dcBCOnTG5PFcpkapS8sLIxhw4aRmprq9N0vZ9xzzz0ArF+/nvfeey/HvK1btzJ16lQA7r333hKLQaSkrDtwyvGM2H03RDpVt2fTSAK93YmNT2P9gb9LKEIREZGyRcmYOCWkTx/8O3fGyMwkdtAg4ufNw5aWu6c0sDdNjJ87j9hBgzEyM/G/5RaC+/Qp5YgLbtSoUfj6+pKWz/YUh+uvv97xXNqwYcO45ppr6NmzJ23atKF169akpKQwePBgOnXqVGIxiJSUBTvszXB7No3A28O5iy7eHhbuu8H+POkXO2KLPTYREZGySMmYOMVksVB90kRHQnby5fEc7NCRk6+9TvKGDaTu3k3yhg2cnPAaB9t34OT48Y5ErPrE1zGZy+4uV7lyZYYPH17i65k0aRJLlizhlltu4cSJEyxdupTff/+d9u3bM2/ePD744IMSj0GkJByJtw9/0bpOwZonXqzN1fZ6R/9ZjoiISEWnDjzEaWYvL8KnTSVh7lziZ88h89gx4mfOJD6PMbjcIyII6dOb4D59ykQiZhiXfhbl9ddf5/V/mlSWpG7dutGtW7cSX49IaUrNyALA36tw/1r8PO31Uv5ZjlRwNiscWgU/zoH4GMhIBg8/CImCJn2gTicwl71m7SIixUnJmBSKyWIh5OGHCe7dm5SNG0lctJiMuDhsKSmYfX3xiIgg6N4e+LZpUyafEROR4ufjYf+XkpReuGQq+by9nq+H/jVVaDYrbHvf/ko8mnv+yV9h7zIIqgEthkCLoVAGLuaJiJQE/ceTIjFZLPi1b49f+/auDkVEXKxmiA97/zrHpkOn6VAvzOn6Gw+eBqBGiE9xhyZlRWYaLBlkT7YAvIKgSW+o3RE8/eF8Evy5Bn6ca0/UvnsWjm6B7h+Be8GHUxERKS90qUlERIpFz2b2DjgW7IwjLcPqVN20DCsLdto77ri/mXM9MUo5YbPC4oH2RMziAbdPhqf3wa2vwNWdoEYL+89bX4Gn9trnWzzs5ZcMstcXEalglIyJiEixaF83jIhgb86mZbJwl3M9Ii7YGcu59CwiQ7xpV7dyCUUoLrXtfdi33J5g9V4MzQflf7fLw8c+v/fifxOybercSEQqHiVjIiJSLCxmE/1aRQEwfvleNv9xukD1Nh86zSvf7AWgX6taWMymkgpRXCX7OTGAW1+FWu0KVq9WO7jlFfvv297X3TERqXCUjImISLHp37oWXRpWJcNqo9+nO5i9JSbfJotpGVZmbY6h34wdZFht3HZNVfr/k8xJBXNolf0ZMK8gaNzLubpNetvrJR6BQ6tLIjoREZe54pOxzZs3c/vttxMSEoKfnx/Nmzdn1qxZxbLsAQMGYDKZMJlMbN26tViWKSJSllnMJqY90NiRkL341R5unLCa8ct/Z+3+U+yMiWft/lO8vPx3Wry6iv9+vceRiE29vzFm3RWrmH6cY//ZpLe9CSJgWK0krV1L3PBHOdm8OX81asSfd99D3P89RtLatRjWf5J4Dx97PYAfZ7sgeBGRknNF96a4dOlS7rvvPmw2G+3atSM0NJTVq1fTr18/fv75Z6ZMmVLoZa9Zs4YZM2ZgMpkuO7aViEhF4uVuYXqv65m5OYYZmw4Tl5DGxxsP8/HGw7nKRoZ4069VLfq3ilIiVpHFx9h/1u6IYbUSP2cOCXPmknnsGD4pKUTE2Z8xPBKfQNKBAyStXIl7eDjBfXoT0rcvptodYcs7kBDjsk0QESkJV2wylpCQQP/+/bFarSxevJju3bsDcPLkSdq0acPUqVO566676Nixo9PLTk9P5z//+Q8NGzYkMDCQzZs3F3f4IiJlmsVs4pE2tejXKor1B/7mix2xHI1PJSUjC18PN2qE+HB/s0ja1a2sZ8SuBBnJANjw5PjjT5D8/fcAWAICCPP3h3+SseqNriO+eQvOfv01WbGxnJrwGmm7dlP9sXvtTXnOJ7toA0RESsYVm4x9/PHHnD17lnvuuceRiAFUqVKFiRMn0r17d6ZMmVKoZOzll1/m0KFDrF+/nueff744wxYRKVcsZhMdo8PoGO38uGNSgXj4Ydjg+KvvkPXDVqKP5jHYM+D+zTdU+eYbqvzzd0zt2iStXMnx5FOER4LJ06/0YhYRKQVX7DNjy5cvB+Dee+/NNe+OO+7Ay8uLVatWkZ6e7tRyf/vtNyZNmsSAAQNo06ZNscQqIiJSroVEEX/Ql6QtvxCYVPC7W+GNGmFydydpy88kHPSF4KiSi1FExAWu2GTsl19+AeD666/PNc/Dw4NrrrmG9PR09u/fX+Bl2mw2Bg0aRGBgIBMnTiy2WEVERMoz47peJBzwBcD0wrNQu/blK9WujfubbxI24kkA4g/4YlzXuyTDFBEpdVdkM8Vz586RmJgIQERERJ5lIiIi2LlzJ0ePHqVRo0YFWu706dPZunUrs2bNIiQkpEgxNmzYMM/pf/zxB1dddVWRli0iIlKakk94kpnihtnDRkBTb1izBjp2hD//zLtC7dqwdi1ERhJUK5W/PWxkpriRcsITv/qlGrqISIm6Iu+MJSf/20TCx8cnzzK+vr65yl5KXFwczz33HB06dKBv375FD1JERKSCOLvkSwCCaqViXvtfsMbApEn5V5g0CSIj4c91mNeNJahWKgCJS5aWfLAiIqWo3N4Zu/fee/ntt9+cqjN79myaN29eoK7mne2Ofvjw4Zw/f5733nvPqXr52bNnT57T87tjJgVnMtl7bgsKCuLw4cMEBQXlKjN27FjGjRvHhAkTGDNmTClHKCJSsWTE2ntL9G3WGLLWw9we2DbWyveK8KlZ71Cp+iksK18Aawa+zVoSv/8gGXFxpRaziEhpKLfJWExMjFPPcwGkptqvrPn7++eYFhAQkG9ZP7/L99y0ePFivv76a1544QWio6OdiklcJzExkalTpzJu3DhXhyIiUqHZ/vmfar71Oax/BDDv+Dr6fr/DMX/+Tfam/Q/8EA9A2Ndr6HLXXzzk60nv6rdirj0c5g7AlpJS+sGLiJSgcpuM7dy5s9B1AwICCAwM5OzZs8TFxdGgQYNcZeL+ufpWo0aNyy5v2bJlAKxcuZL169fnmPfTTz8BMGzYMAICAnj00Ufz7MGx3MvKgi+/tP/etSu4le1dy2w24+bmxrRp03jiiScIDg52dUgiIhWW+Z9HAtJTUniuSmVWnQ9m7mRfJn4Qy8FOfkTUNOFns/FbVCCB69J5/D+RHHN3Y1KlYH6sEsoL2cncP48QiIhUFGX7jLkENWrUiPXr17N79+5cyVhmZia//fYbnp6e1KtXr8DL3Lp1a77zfvzxRwC6du1aqHjLvLVr4b777L+vXAmdOrk0nMtxd3fnkUce4d1332XKlCm8/PLLrg5JRKTC8oiM5Py+fXy/aDKrmsXhbnan/+1jiL49iMa/fA4JMfYBnZv5wS1RzLvuAb7kLJN2vsGqo6u5cft+rgM88ul0S0SkvLoiO/AA+1hiAIsWLco1b/ny5aSnp3PzzTfj5eV12WXNnDkTwzDyfLVv3x6ALVu2YBgGTzzxRLFuh8sYBqSk/PuaP//fefPn55zn5PN3peXZZ5/F09OTN998k/j4eFeHIyJSYQX26A7AVRtj8LW68X6n93mwfi+86t8B98+FIRvh8Z/sP++fi3f9O3mwfi/e7/Q+vlY3am+KASDo3h6u2wgRkRJwxSZjAwcOJCAggK+++oolS5Y4pp86dYpRo0YB8NRTT+WqFx0dTXR0NMeOHSu1WMukLVvAz+/f1yef/Dvvk09yzrvEHUNXCg8PZ9CgQSQlJfHGG2+4JAbDMDh+/DiZmZm55iUnJ+c5XUSkvPFu3YozwRb80uHlxI40r9a8QPWaV2vOSwkd8EuHM8FueLVqWcKRioiUris2GQsJCeHTTz/FbDZz77330rFjR+677z7q1avHoUOHeOyxx7j55ptz1du/fz/79+/XSfLs2SVTtpQ988wzeHl58fbbb3PmzJnLlu/QoQMmk8mp19q1a/Nc1qeffkpwcDDh4eFERESwcuVKAGJjY2nXrh3+/v74+/szcOBAR4cyIiLl0eaTW1l+vb2VRMSMVaRs3VageilbtxI5czUAy663sflk2by4JyJSWFfsM2MAPXr0YP369YwfP56tW7eSkZFB/fr1GT58OP3793d1eGXbs8/anw3Lb8DObLVr28uWUdWrV2fw4MG89dZbTJ48mQkTJlyyfJcuXYiKinJqHVWrVs01bdu2bfzf//0fr776Ktdccw3Tp0/n9ttvZ8OGDQwZMoTMzEwWLFiAxWLhxRdf5Pnnn2fKlClOrVdEpKxYenApq5uauOVsFNV2xhA7aBBhY0YT1L07Zm/vXOVtaWkkLl7Cqddfh8xM/moWxbfNYsk4uJR2Ee1csAUiIiXDZDg7oJa4VPY4Y/mNQwZgs9kc3f7Xq1cPs7mEboAePQodO+afkNWube/YIzKyZNZfSCaTCU9PT9LT0wH466+/qF27Nm5ubhw+fJjQ0NASH2ds8ODBVK5cmVdeeQWwf2Z33nkna9euJTQ0lN27dxMaGgrAvn37aNq0KUlJSY4x0i5Wap+5iEgh3Pv1vexP2M97bd8kauqXJP3TEsAcGEhQt274tm6F2dcXW0oKKRs3kbh0KbZz5wDwv+UWYp64h6EbHic6JJqFdy105aaIiORSkPPz/OiMTQqvRg2YNCn/+ZMmlblELC/VqlVjyJAhJCcnM+lS21OMjhw5wo033uj422w2M336dNLS0nj00UcdiRjYn1M0m82cPn26VGITESluqVnZY3eGED5tKlWeGYN7eDi2s2eJnzmT2EGDOfJQL2IHDSZ+1ixs587hHhFBlWfGED5tKr5+9uFHUjI1zpiIVCxXdDNFKQZ59EbpsHgxdO9eerEUwejRo/nggw+YPn06I0aMyLfca6+9xr59+5xa9pgxY3INBh4ZGcmfF91RXLVqFQAffPABw4YNcww4furUKaxWK5UqVXJqvSIiZYWPm32cseTMZEwWCyEPP0xw796kbNxI4qLFZMTFYUtJwezri0dEBEH39sC3TRtMFgvwbxLm665xxkSkYlEyJkXz+ef//j5smL0b+/fes//92Wcwb55r4nJS1apVGTp0KFOmTGHixIn45jOw6IoVK1i3bp1Ty+7Xr1+uZKxfv37cf//9NG3alEaNGvHVV1/x+OOP89FHH/HSSy9xzz33MGPGDLy9vRkyZAh9+vRR00MRKbci/SPZn7Cfrce30ia8DQAmiwW/9u3x+2cImEvZcnwLABF+GmdMRCoWnd1J0cTEQNu28P33MH06vPsufPcdtGsHR464OjqnjB49Gh8fH959911OnjyZZ5m1a9fmO6Zcfq8OHTrkWk6bNm0YO3Yst912G/7+/gwdOpSJEycycOBAVqxYwYkTJ6hZsyZhYWFkZWXx+uuvl/DWi4iUnG5XdwNg6aGlpGWlOVU3LSuNpYeW5liOiEhFoWRMiqZmTVi/Hjp3/nfaLbfAunX2Z8rKkbCwMIYNG0ZqaiqzZs0q8fUNGjSIxMREjh49yunTp3n00UcBaNCgAb/99hsHDx7kzJkzfP311wQGBpZ4PCIiJaV19daE+4VzLuMcXx760qm6Sw8uJSkjiXC/cFpXb10yAYqIuIiSMZELjBo1Cl9fX9LSnLtyW1hms5nIyEg8PDxyTDeZTNSpU4eQkJBSiUNEpCRZzBYein4IgEk7JrH9r+0Fqrftr21M3jkZgF71e2ExW0osRhERV1AyJnKBypUrM3z4cFeHISJS4fSq34tONTqRactkyKohfL7v83ybLKZlpfHZ3s8YumoombZMOtfsTK/6vUo5YhGRkqdxxsqZMjXOmJQZ+sxFpDxIz0rnmQ3PsOqovffYAI8AutbpSsvqLfF19yUlM4XNxzfz5aEvScpIAqBzzc682uZVvNy8XBm6iEi+ijLOmHpTFBERkVLh5ebF5PaT+WzfZ8zbO49jyceY/ftsZv8+O1fZcL9wetXvRa/6vTCbdIFJRComJWMiIiJSaixmC30a9OGh6IfYdHwTSw8uJS45jpTMFHzdfYnwi6Db1d1oXb21nhETkQpPyZiIiIiUOovZQruIdrSLaOfqUEREXEb3/UVERERERFxAyZiIiIiIiIgLKBkTERERERFxASVjFZDJZHL8brPZXBiJlJYLP+cLP38RERERKbuUjFVAJpMJDw8PAFJSUlwcjZSG7M/Zw8NDyZiIiIhIOaHeFCsof39/zpw5w8mTJwHw9fXVQMAVkM1mIyUlxfE5+/v7uzgiERERESkoJWMVVKVKlUhJSSE9PZ3jx4+7OhwpBV5eXlSqVMnVYYiIiIhIASkZq6AsFgs1atTgzJkzJCUlkZGR4eqQpIR4eHjg7+9PpUqVsFg0QKqIiIhIeaFkrAKzWCyEhYURFhaGYRgYhuHqkKSYmUwmPSMmIiIiUk4pGbtC6KRdRERERKRsUY8OIiIiIiIiLqBkTERERERExAWUjImIiIiIiLiAkjEREREREREXUDImIiIiIiLiAkrGREREREREXEDJmIiIiIiIiAtonDERESk5NiscWgU/zoH4GMhIBg8/CImCJn2gTicwW1wdpYiIiEsoGRMRkeJns8K29+2vxKO555/8FfYug6Aa0GIItBgKZjXWEBGRK4uSMRERKV6ZabBkkD3ZAvAKgia9oXZH8PSH80nw5xr4ca49UfvuWTi6Bbp/BO7eLg1dRESkNCkZExGR4mOzwuKBsG85WDzg1lftidjFSdbVnaDjc/DTPHsytneZPYG7b5aaLYqIyBVDbUJERKT4bHv/30Ss92JoPij/u10ePvb5vRfby+9dBts+KN14RUREXEjJmIiIFI/s58TAfkesVruC1avVDm55xf77tvftyxEREbkCKBkTEZHicWiV/RkwryBo3Mu5uk162+slHoFDq0siOhERkTJHyZiIiBSPH+fYfzbpbW+C6AwPH3s9gB9nF29cIiIiZZSSMRERKR7xMfaftTsWrn52vYSY4ohGRESkzFMyJiIixSMj2f7T079w9T397D/PJxdPPCIiImWckjERESkeHtnJVFLh6mcnYdlJmYiISAWnZExERIpHSJT9559rClc/u15wVHFEIyIiUuYpGRMRkeLRpI/9549zISPVuboZqRd0ANK3eOMSEREpo5SMiYhI8ajTCYJqQHoi/DTPubo/zoX0sxBUE+rcXCLhiYiIlDVKxkREpHiYLdBiiP33756Fw+sLVu/PdfD9c/bfWwyxL0dEROQKoGRMRESKT4shUP8usGbA3B6w/aP8myxmpMK2D2Hevfby9e/+N5kTERG5ApgMwzBcHYQUXMOGDQHYs2ePiyMREclHZhosGQR7l9n/9gqyD+hcu6O9p8TzyfDHD/DTP00TwZ6Idf8Q3L1dFraIiEhhFOX83K24gxERkSucuzfcNwu2fQDb3oPEo7DlHfvrYkE17XfDWgwBsxpriIjIlUXJmIiIFD+zBVoOgxb/gUOr4cfZkBBjvyvm6Wfvvr5JX3tnHXpGTERErlBKxkREpOSYLVD3FvtLREREclAyJiJFZrUZrDtwigU74jgSn0pqRhY+Hm7UDPGhZ7MI2tcNw2I2uTpMERERkTJFyZiIFJrVZjBj02Fmbo4hLiEt1/y9f51jxZ4TRAR7069VFANa18KspOyKYlitJG/YwNnFS8iIjcWWmorZxwePyEgCe3THr21bTBY1UxQRkSuTkjERKZT0TCtPzP+JFXtOABDo7U7PphG0rhOKv5cbSelZbDx4moW74ohLSGP8N3vZGZPAtAca4+Wuk++KzrBaiZ8zh4Q5c8k8dgwMA//kJACS/Pw5v28fSStX4h4eTnCf3oT07YtJHXiIiMgVRl3blzPq2l7KAqvNYNi8XXy35yQeFjPP31mfnk0j80yy0jKsLNwVy/jle8mw2ujSsCrTe12vZosVmC09neMjR5G0ciUA5sBAwq69luBPPgYgbcoUzp3PIHHpUmxn7V3b+3fuTPVJEzF7ebksbhERkcIoyvm5LkOKiNNmbDrsSMRmDmhG35ZR+d7t8vaw0LdlFDMHNMPDYmbFnhPM3BxTugFLqTGsVo6PHEnS999jtlioMnIEV//vG4L597qf9549VPm/R7n6m+VUGTkCk5sbSStXcnzkKAyr1YXRi4iIlC4lYyLiFKvNcCRTz99Zn1ZXhRaoXqurQnnujvoAzNx8GKtNN+Urovg5c0hauQqfrEzq/b6HkIEDMYeGwief/Fvok0/Azw9z5cqEDBxIzccfx+TuTtLKlSTMneu64EVEREqZkjERccq6A6eIS0gj0Nud+26IdKpuz6aRBHq7ExufxvoDf5dQhOIqhtVKwhx7MlUlKqrA9bx37yJs9GgA4ufM1d0xERG5YigZExGnLNgRB0DPphF4ezjXEYe3h4X7bogA4IsdscUem7hW8oYNZB47hjkwEI/33oPatS9fqXZtePZZgnp0xxwYSGZcHCkbN5Z8sCIiImWAkjERccqR+FQAWtcpWPPEi7W52l7v6D/LkYrj7OIlAAR17465bl1Ys+bSCVnt2rB2LURGYvb2JqhbNwASFy0uhWhFRERcT8mYiDglNSMLAH+vwo2M4edpr5fyz3Kk4siItd/t9G3V0j6hRg2YNCn/CpMmQeS/TV19W7eyLycursRiFBERKUuUjImIU3w87MlUUnrhkqnk8/Z6vh4a5rCisaXa73aaff3+nbhoUf4VFue8A2b29bUvJyWl2GMTEREpi5SMiYhTaob4ALDp0OlC1d940F6vxj/LkYrD7GP/TG0pyf9O/Pzzf38fNgyGDv33788+y1E/OwnLTspEREQqOiVjIuKUns3sHXAs2BlHWoZzvd6lZVhZsNPelO3+Zs71xChln8c/TQ5TNm/5d2JMDLRtC99/D9Onw7vvwnffQbt2cORIjvopmzbblxMRUVohi4iIuJSSMRFxSvu6YUQEe3M2LZOFu5zrEXHBzljOpWcRGeJNu7qVSyhCcZXAHt0BSFyyBFtamn1izZqwfj107vxvwVtugXXr7M+U/cOWlkbikn86ALm3R6nFLCIi4kpKxkTEKRaziX6togAYv3wvm/8oWHPFzYdO88o3ewHo16oWFrOppEIUF/Fr2xb38HBsZ886EquCSly8BNu5c7hHRODbpk0JRSgiIlK2KBkTEaf1b12LLg2rkmG10e/THczeEpNvk8W0DCuzNsfQb8YOMqw2brumKv3/SeakYjFZLAT36Q3AqddeJ2XrtgLVS9m6lVOvvw5ASJ/emCzOjV8nIiJSXpkMwzBcHYQUXMOGDQHYs2ePiyORK116ppUn5v/Eij0nAAj0due+GyJoc3Uofp5uJJ/PYsPB0yz8p2kiwG3XVGXq/Y3xctfJdkVlWK0ce+JJklauxOTuTuXRo9hzY1WWxn5DbFIsqVmp+Lj5EOkfSbfIO2iw5S9OT5yEkZmJ/y23ED5tKiazrhOKiEj5UZTzcyVj5YySMSlLrDaDmZtjmLHpMHEJafmWiwzxpl+rWvRvFYVZzRMrPFt6OsdGjiR55SoAkr1gzXUmfokyke5pwuu8QaPDBh1+NfBLt9fxu6Uz4RMnYvbycmHkIiIizivK+bkG+hGRQrOYTTzSphb9WkWx/sDffLEjlqPxqaRkZOHr4UaNEB/ubxZJu7qV9YzYFSTDDaZ2teBuM3P7ThthZ+Gu7QZ3bc997e9kEHx7g5mMe8xMcAOlYiIiciVRMiYiRWYxm+gYHUbH6DBXhyIuZrVZGbNhDKvjVuN+oyeNhz3NtSfDSF26jIy4OGwpKZh9ffGIiMC72138VOUUK3e/QWbcatjwDJPbT8ZiVjNWERG5MigZExGRYjNv7zxWH12Nu9md9zu9T/NqzaEhhNzUOc/yDwJXhVzNkFVDWHV0FZ/t+4w+DfqUbtAiIiIuoqekRUSkWFhtVj7b9xkAI5uNtCdiBdC8WnNGNB0B2JM5q825wcRFRETKKyVjIiJSLDYd38Sx5GMEeATQtU5Xp+p2u7obAR4BHEs+xqbjm0omQBERkTJGyZiIiBSLpQeXAtCtTje83bydquvt5u1I4LKXIyIiUtEpGRMRkWIRmxQLwI3VbyxU/ZbVWwIQlxxXbDGJiIiUZUrGRESkWKRmpQLg5+5XqPq+7r4ApGSmFFtMIiIiZZmSMRERKRY+bj4AJGcmF6p+dhKWnZSJiIhUdErGRESkWET6RwKw9fjWQtXfcnwLABF+EcUWk4iISFmmZExERIpFt6u7AbD00FLSstKcqpuWlcbSQ0tzLEdERKSiUzImIiLFonX11oT7hXMu4xxfHvrSqbpLDy4lKSOJcL9wWldvXTIBioiIlDFKxkREpFhYzBYein4IgEk7JrH9r+0Fqrftr21M3jkZgF71e2ExW0osRhERkbJEyZiIiBSbXvV70alGJzJtmQxZNYTP932eb5PFtKw0Ptv7GUNXDSXTlknnmp3pVb9XKUcsIiLiOibDMAxXByEF17BhQwD27Nnj4khERPKWnpXOMxueYdXRVQAEeATQtU5XWlZvia+7LymZKWw+vpkvD31JUkYSAJ1rdubVNq/i5eblytBFREScVpTzcyVj5YySMREpD6w2K5/t+4x5e+dxLPlYvuXC/cLpVb8Xver3wmxSYw0RESl/inJ+7lbcwYiIiFjMFvo06MND0Q+x6fgmlh5cSlxyHCmZKfi6+xLhF0G3q7vRunprPSMmIiJXLCVjIiJSYixmC+0i2tEuop2rQxERESlz1CZERERERETEBZSMiYiIiIiIuICSMRERERERERdQMiYiIiIiIuICSsZERERERERcQMmYiIiIiIiICygZExERERERcQElYyIiIiIiIi6gZExERERERMQFlIyJiIiIiIi4wBWfjG3evJnbb7+dkJAQ/Pz8aN68ObNmzSr08mw2Gx999BFt2rQhODgYb29vateuTa9evdizZ08xRi4iIiIiIuWZm6sDcKWlS5dy3333YbPZaNeuHaGhoaxevZp+/frx888/M2XKFKeWl5qayl133cUPP/xAcHAwbdq0wcvLi8OHD/PFF19w22230bBhwxLaGhERERERKU+u2GQsISGB/v37Y7VaWbx4Md27dwfg5MmTtGnThqlTp3LXXXfRsWPHAi+zf//+/PDDDwwYMIC3334bHx8fx7y//vqLzMzMYt8OEREREREpn67YZooff/wxZ8+e5Z577nEkYgBVqlRh4sSJAE7dGfvhhx9YsGABzZo146OPPsqRiAFUq1aNGjVqFE/wIiIiIiJS7l2xydjy5csBuPfee3PNu+OOO/Dy8mLVqlWkp6cXaHkffPABAE8++SRm8xX7toqIiIiISAFdsc0Uf/nlFwCuv/76XPM8PDy45ppr2LlzJ/v376dRo0aXXd4PP/wAQKdOnfjtt99YuHAhJ06coGrVqtx2223ceOONxbsBIiIiIiJSrl2Rydi5c+dITEwEICIiIs8yERER7Ny5k6NHj142GTt58iSnT58mODiYTz75hOeeew6bzeaY/9JLL9G7d28+/fRT3N3dCxRjfh19/PHHH1x11VUFWoaIiIiIiJRdV2R7uuTkZMfvFz/blc3X1zdX2fwkJCQAkJSUxDPPPEOvXr3Yv38/CQkJLFq0iNDQUObOncsLL7xQDNGLiIiIiEhFUG7vjN1777389ttvTtWZPXs2zZs3xzCMy5YtSJlsVqsVgKysLFq2bMns2bMd83r06IGXlxd33nknb731Fs8++ywBAQGXXWZ+Y5Kpa3wRERERkYqh3CZjMTEx7N+/36k6qampAPj7++eYlldylF3Wz8/vssu9cHkDBgzINf+OO+6gSpUqnDx5ku3bt9OpUyen4hYRERERkYqn3DZT3LlzJ4ZhOPXq0KEDAAEBAQQGBgIQFxeX5/KzpxekO/rq1avj4eEBQM2aNfMskz391KlTTm2niIiIiIhUTOU2GSuq7E45du/enWteZmYmv/32G56entSrV++yy3Jzc+Oaa64BID4+Ps8yZ86cAQp2p01ERERERCq+KzYZu+OOOwBYtGhRrnnLly8nPT2dm2++GS8vrwIt7+677wZgzZo1uebFxMQQExMDQJMmTQoZsUgZZrPCge/gi97wXht4s7H95xe97dNtVldHKCIiIlLmmAxneqqoQOLj46lVqxbnzp1j8eLFdO/eHbA3I2zdujWHDh1i1apV3HzzzTnqRUdHA7B69WrCw8Md0//++2/q1KlDeno6y5cvp3PnzoC9N8b77ruPFStWcMcddzgGmy6s7A488uvgQ6RU2ayw7X37K/Fo/uWCakCLIdBiKGhQdBEREalAinJ+fsUmYwCLFy+mZ8+eGIZB+/btCQ0NZdWqVSQmJvLYY4/x5ptv5qpjMpkAOHz4MFFRUTnmLV26lJ49e2Kz2bjxxhsJCwtj69atnDhxglq1arFhw4YcCVxhKBmTMiMzDZYMgr3L7H97BUGT3lC7I3j6w/kk+HMN/DgX0hPtZerfBd0/AndvV0UtIiIiUqyUjBXBpk2bGD9+PFu3biUjI4P69eszfPhw+vfvn2f5SyVjADt27OCVV15h48aNJCUlERkZSdeuXXnmmWeoVKlSkeNVMiZlgs0KC/rCvuVg8YBbX7UnYnklWRmp8NM8+O5ZsGbYE7L7ZoHZUvpxi4iIiBQzJWNXECVjUiZsmW5Priwe0Hsx1Gp3+TqH18PcHvaE7NYJ0HJYyccpIiIiUsKKcn6uhzdExDnZz4mB/Y5YQRIxsJe75RX779veV6ceIiIicsVTMiYizjm0yt5Zh1cQNO7lXN0mve31Eo/AodUlEZ2IiIhIuaFkTESc8+Mc+88mvcHDx7m6Hj72egA/zi7euERERETKGSVjIuKc+Bj7z9odC1c/u15CTHFEIyIiIlJuuRV1AcOGFf4hfJPJxPTp04sagoiUpoxk+09P/8LV9/Sz/zyfXDzxiIiIiJRTRU7G3n//fafrmEwmDMNQMiZSHnlkJ1NJhaufnYRlJ2UiIiIiV6giJ2MzZsxwqvwvv/zC+++/T1paWlFXLSKuEBIFJ3+1D+h8dSfn6/+5xv4zOKo4oxIREREpd4qcjD388MMFKvf7778zbtw4Fi9ejM1mIyIigmeffbaoqxeR0takD+xdBj/OhY7POdeJR0bqBR2A9C2Z+ERERETKiRLvwGPfvn08+OCDXHfddSxcuJBq1arxzjvvcOjQIYYMGVLSqxeR4lanEwTVgPRE+Gmec3V/nAvpZyGoJtS5uUTCExERESkvSiwZ279/P7169eLaa6/liy++oGrVqrz11lv88ccfDBs2DA8Pj5JatYiUJLMFWvxzIeW7Z+Hw+oLV+3MdfP+c/fcWQ+zLEREREbmCFXsydvDgQfr06cM111zD559/TuXKlZk2bRp//PEHjz76qJIwkYqgxRCofxdYM2BuD9j+kb0JYl4yUmHbhzDvXnv5+nf/m8yJiIiIXMFMhmEYxbGgP/74g5deeonPP/+crKwsqlSpwujRoxkyZAheXl7FsQoBGjZsCMCePXtcHIlc8TLTYMkg+/NjAF5B9gGda3e095R4Phn++AF++qdpItgTse4fgru3y8IWERERKU5FOT8vcgcef/75Jy+//DLz5s0jKyuLsLAwRo0axdChQ/H21gmXSIXl7g33zYJtH8C29yDxKGx5x/66WFBN+92wFkPArLHmRURERKAYkrHo6GisVive3t48/vjjDB8+HB8fH5KSkkhKuvw4RGFhYUUNQURcxWyBlsOgxX/g0Gr4cTYkxNjvinn62buvb9LX3lmHnhETERERyaHIzRTNZjMmk6lwKzeZyMrKKsrqrzhqpigiIiIiUna4tJlijRo1Cp2MiYiIiIiIXKmKnIzFxMQUQxgiIiIiIiJXliI/SZ+RkVEccYiIiIiIiFxRipyMhYaGcu+99zJz5kz+/vvv4ohJRERERESkwityM8UWLVqwbNkylixZgtlspkWLFtx9993ceeedjofZREREREREJKdiGfQ5KSmJb7/9lmXLlvHtt98SHx+PyWQiKirKkZi1b98eN7ci535XPPWmKCIiIiJSdhTl/LxYkrEL2Ww2Nm3axFdffcXy5cs5cOAAJpMJf39/unTpwt13381tt91GcHBwca72iqFkTERERESk7ChTydjFDh48yFdffcXXX3/Nli1bsNlsWCwWWrZsyd13381dd91F3bp1SzKECkXJmIiIiIhI2VGmk7ELJSQksHz5cpYtW8b333/PuXPnMJlMvPbaa4wcObK0wijXlIyJiIiIiJQdLh302RnBwcH06dOHPn36kJmZydq1a1m2bJmaLIqIiIiIyBXHZT1quLu707lzZzp37uyqEERERERERFymyOOMDRgwgE8//TTPeV9//TU//fRTnvP++9//csMNNxR19SIiIiIiIuVSkZOxmTNnsnHjxjznde3albfeeivPeUePHs03URMREREREanoipyMiYiIiIiIiPOUjImIiIiIiLiAkjEREREREREXUDImIiIiIiLiAkrGREREREREXEDJmIiIiIiIiAuYDMMwirIAs9mMyWQqdH2r1VqU1V9xGjZsCMCePXtcHImIiIiIiBTl/NytOAIobD5XlCRORERERESkPCtyMmaz2YojDhERERERkSuKnhkTERERERFxASVjIiIiIiIiLqBkTERERERExAWUjImIiIiIiLhAsfSmKCJXNsNqJXnDBs4uXkJGbCy21FTMPj54REYS2KM7fm3bYrJYXB2miIiISJmiZExECs2wWomfM4eEOXPJPHYs1/zz+/aRtHIl7uHhBPfpTUjfvpjMuiEvIiIiAkrGRKSQbOnpHB85iqSVKwEwBwYS1L07vq1aYvb1w5aSTMqmzSQuXUrmsWOceu110nbtpvqkiZi9vFwcvYiIiIjrmYzCjtgsLlGUEb5FiothtXLsiSdIWrkKk7s7YWNGE9SjR55Jli0tjcQlSzj12usYmZn4d+5M+LSparYoIiIiFUJRzs/VXkhEnBY/Z44jEYv86CNCevXK926X2dubkF69iPzoI0zu7iStXEnC3LmlHLGIiIhI2aNkTEScYlitJMyxJ1NhY0bje2OLAtXzvbEFYaNHAxA/Zy6G1VpiMYqIiIiUB0rGRMQpyRs2kHnsmOMZMWcE9eiOOTCQzLg4UjZuLKEIRURERMoHJWMi4pSzi5cAENS9O2Zvb6fqmr29CerWDYDERYuLPTYRERGR8kTJmIg4JSM2FgDfVi1zz8zKgkWL7K+srDzr+7ZuZV9OXFyJxSgiIiJSHqhrexFxii01FQCzr1/umWvXwn332X9fuRI6dcpVxOzra19OSkpJhSgiIiJSLujOmIg4xezjA4AtJRkMA1JS/n3Nn/9vwfnzc877ZxSN7CQsOykTERERuVIpGRMRp3hERgKQsnkLbNkCfn7/vj755N+Cn3ySc97WrfZ6mzbblxMRUeqxi4iIiJQlSsZExCmBPew9KCYuWYLx6acFrzh7tmMAaICge3uURHgiIiIi5YaSMRFxil/btriHh2M7e5az114LtWtfvlLt2vDssyQuXoLt3DncIyLwbdOm5IMVERERKcOUjImIU0wWC8F9egNw4uNPSH1jyqUTstq1Ye1aUo4d49TrrwMQ0qc3JoulNMIVERERKbOUjImI00L69MG/c2eMzEyOvvACSXfdlW9Z2/jxxK9bT+ygwRiZmfjfcgvBffqUYrQiIiIiZZOSMRFxmsliofqkiY6EzPbpjHzLJj/6f5wcP96RiFWf+Doms756RERERHRGJCKFYvbyInzaVKo8M4bApHOO6fFBQSQEBTn+Dog/g3tEBFWeGUP4tKmYvbxcEK2IiIhI2aNBn0Wk0EwWCyEPP4zRti22bt2Ij6pFktmMLSWF9JRkQmJiyJo6lat69NAzYiIiIiIXUTImIkVmql0by88/UxmofNE8T1cEJCIiIlIOqJmiiIiIiIiICygZExERERERcQElYyIiIiIiIi6gZExERERERMQFlIyJiIiIiIi4gJIxERERERERF1AyJiIiIiIi4gJKxkRERERERFxAyZiIiIiIiIgLKBkTERERERFxASVjIiIiIiIiLqBkTERERERExAWUjImIiIiIiLiAkjEREREREREXcHN1ACJS/lltVjYd38TSg0uJTYolNSsVHzcfIv0j6XZ1N1pXb43FbHF1mCIiIiJlipIxESk0q83KvL3z+GzfZxxLPpZr/v6E/aw6uopwv3Aein6I3g16YzbphryIiIgIKBkTkUJKz0rnmQ3PsOroKgACPALoVqcbN1a/ET93P5Izk9lyfAtfHvqSY8nHmLRzEj+e+pEJbSfg5ebl4uhFREREXE/JmIg4zWqzMmbDGFYfXY272Z2RzUbSrU63XElWm/A2PNrkUb489CWTdkyyJ24bYHL7yWq2KCIiIlc8tRcSEafN2zvPkYi93+l9Hox+MN+7Xd5u3jwY/SDvd3ofd7M7q46u4rN9n5VyxCIiIiJlj5IxEXGK1WZ1JFMjm42kebXmBarXvFpzRjQdAdiTOavNWmIxioiIiJQHSsZExCmbjm/iWPIxAjwC6Fqnq1N1u13djQCPAI4lH2PT8U0lE6CIiIhIOaFnxqRArDaDdQdOsWBHHEfiU0nNyMLHw42aIT70bBZB+7phWMwmV4cppWDpwaXw/+3de1SVdaL/8c9ms7ltQCIr4xaKJmpZVmLenYJyxkxRsxJvXSxrZkptypo6v3XOGSsLy6zmpGWTB/NSmtjSGksr07iYmK6OJoaGcdEmHUG5Kbf9+2MHjbNBAYEvG96vtfYCn+f5Pnxs7ZX7w/d5vo+k+O7x8vX0bdRYX09fje0+VknfJSk5K1nDwoa1REQAAAC3QBnDOVVVO/ROSraWpR5WXkGZy/79R09p076fFHaRr6YPitS9g7vKg1LWruUW5UqSbgy5sUnjB4YMVNJ3ScorzmvOWAAAAG6HMoZ6na6o0qzVe7Rp30+SpE6+Nk28IUyDu3dWgI+nik5X6qus41qzK095BWWa99F+ZRwu0Ct3XSsfGyvltVellaWSJH+bf5PG2212SVJJRUmzZQIAAHBHlDHUqaraoUdX79Yn+/4hL6uHnrmtlybeEO5Sskb0vFSP3dJTa3blat7G/dq07yfNWr1Hf024jssW2yk/Tz9JUnFFcZPG15SwmlIGAADQUbGAB+r0Tkp2bRFbdm9/TR0YWe9sl6+XVVMHRmrZvf3lZfXQpn0/aVnq4dYNjFYTHhAuSUo/kt6k8WlH0iRJYf5hzZYJAADAHVHG4KKq2lFbpp65rZcGRXVu0LhBUZ319KhekqRlqdmqqna0VEQYFN8jXpKUfDBZZZWu9xGeS1llmZIPJp91HgAAgI6KMgYXX37/s/IKytTJ16Y7rg9v1NiJN4Srk69NuSfKtO37Yy2UECYNDhmsUP9QnSo/pfUH1zdqbHJWsorKixTqH6rBIYNbJiAAAICboIzBxfs7navcTbwhTL5ejVuIw9fLqjuud15+9t7O3GbPBvOsHlZNip4kSUrcmaivj37doHE7ju7QgowFkqSEXgmyerDICwAA6NgoY3Dx4wnnanmDuzfs8sR/N6SHc1zOL+dB+5PQK0GxEbGqqK7QzC0ztSpzVb2XLJZVlmnl/pV6aMtDqqiuUNwVcUroldDKiQEAANoeVlOEi9LySklSgE/T3h7+3s5xJb+cB+2P1cOq54c+L22XtuRs0XM7ntPru1/X2O5jNTBkoOw2u0oqSpR6JFXrD65XUXmRJCnuijg9N+Q5eVj4PRAAAABlDC78vJxvi6LTTStTxWec4+xevL3aMx9PHy0YvkArM1dqxf4Vyi/OV9J3SUr6Lsnl2FD/UCX0SlBCrwSKGAAAwC/4tAwXVwT7af/RU0o5eFwjel7a6PFfZR2XJEUE+zV3NLQxVg+rpvSeoknRk5RyJEXJWcnKK85TSUWJ7Da7wvzDFN8jXoNDBnOPGAAAwL/p8GUsNTVV8+bNU3p6usrLy9W7d2/9/ve/17Rp0xp9roKCAs2fP18ffvihfvzxR0lSt27dFB8fryeeeEKBgYHNHb9FTOwfpk37ftL7GXmaE9ezUYt4lJVX6f0M58Idd/Zv3EqMcF9WD6uGhQ3TsLBhpqMAAAC4jQ59vVBycrKGDRumTZs2qW/fvho5cqSysrI0ffp0zZkzp1HnOnbsmPr3768XX3xRRUVFGjlypGJjY/WPf/xDzz77rAYMGKDCwsKW+Ys0s+FXXqqwi3x1sqxCa3Y1bkXE9zNydep0pcKDfTXsyktaKCEAAADg/jpsGSsoKNA999yjqqoqrV27Vlu3btXatWuVmZmp7t27a+HChfriiy8afL7nn39ehw4dUnx8vH744QclJydrw4YNys7O1pAhQ5SZmamFCxe24N+o+Vg9LJo+KFKSNG/jfqUeOt6gcakHj+vZj/ZLkqYP6iqrh6WlIgIAAABur8OWsaVLl+rkyZMaM2aMxo0bV7v9sssu04svvihJevnllxt8vm3btkmS5s6dK29v79rtAQEBeuyxxyRJO3fubI7oreKewV01sk8XlVdVa/rfdiop7bDKyqvqPLasvEr/m3pY09/ZqfKqav32qi6655cyBwAAAKBuHfaesY0bN0qSJkyY4LJv1KhR8vHx0ZYtW3T69Gn5+Pic93z/WsDqExwc3Pighlg9LHrlrms1a/Uebdr3k/7fh/v00qff647rwzSkR2f5e3uq+Eyltmcd15pfLk2UpN9e1UUL77xWHsyKAQAAAOfUYcvYt99+K0m67rrrXPZ5eXnpqquuUkZGhg4cOKBrrrnmvOeLi4tTamqqXnjhBa1ataq2nBUVFWnBggWS1KRFQUzysVn114TrtCz1sN5JyVZeQZmWfpWtpV9luxwbHuyr6YO66p5BkRQxAAAAoAE6ZBk7depU7WIaYWFhdR4TFhamjIwM5eTkNKiM/elPf9IXX3yh5ORkdevWTQMGDFBlZaVSU1NltVr19ttvKy4ursEZ+/TpU+f2Q4cOKSoqqsHnuVBWD4vuG9JV0wdFatv3x/TezlzlnChVSXml7F6eigj20539wzXsyku4RwwAAABohA5ZxoqLi2u/9/Or+1lYdrvd5dhz8ff316ZNmzRjxgytWLFCycnJtftuv/12XX/99ReQ2Dyrh0W/ib5Uv4lu/HPHAAAAALhy2zI2YcIE7d27t1FjkpKSFBMTI4fDcd5jG3LMv8rJydGoUaN09OhRJSUlaeTIkZKkv//975o9e7aGDBmizZs368Ybb2zQ+fbt21fn9vpmzAAAAAC4F7ctY4cPH9aBAwcaNaa0tFSSc4XDf91W18OYa4719/dv0LmnTZumvXv3av369RozZkzt9qlTp8rf31/jx4/XnDlzlJqa2qjMAAAAANonty1jGRkZTR4bGBioTp066eTJk8rLy1Pv3r1djsnLy5MkRUREnPd8ubm52rp1q7y9vTV69GiX/WPGjJG3t7fS09MbvDojAAAAgPatwz5nrGZRjm+++cZlX0VFhfbu3Stvb2/17NnzvOeqKW52u10eHq7/Sa1Wq/z8/ORwOGoXDgEAAADQsXXYMjZq1ChJ0tq1a132bdy4UadPn9bNN9/coFmsLl26SJJOnDih7GzXZd8PHTqkgoIC2e12de7c+QKTAwAAAGgPOmwZu//++xUYGKgPP/xQ69atq93+888/64knnpAkzZkzx2VcdHS0oqOjlZ+fX7uta9eu6tu3ryTpwQcf1MmTJ2v3FRYW6sEHH5QkjR07Vp6ebntlKAAAAIBmZHE0dtnAduSDDz7QxIkT5XA4NHz4cHXu3FlbtmxRYWGhHnnkES1atMhljMXifJZWdna2IiMja7fv2LFDsbGxKi4uVufOnTVgwABJUnp6uv75z38qMjJSKSkpCgkJuaDMNasp1rfaIgAAAIDWcyGfzzvszJgkjR8/Xtu2bdOtt96qPXv26OOPP1ZUVJT+9re/1VnEzmXAgAHas2ePZsyYocDAQG3ZskWff/65LrvsMj355JPatWvXBRcxAAAAAO1Hh54Zc0fMjAEAAABtBzNjAAAAAOBmWE0CDVNdJR3cIu1eLp04LJUXS17+UnCk1G+K1D1W8rCaTgkAAAC4DcoYzq26Stqx2PkqzHHd/4//k/ZvkIIipAEzpQEPSXU8aw0AAADA2ShjqF9FmbRuhrNsSZJPkNRvstTtN5J3gHSmSPrhC2n3u86i9smfpZw0adxbks3XaHQAAACgraOMoW7VVdIH90uZGyWrl3Trc84i9u8lq0es9JunpT0rnGVs/wZngbvjf7lsEQAAADgHridD3XYs/rWITf5AiplR/2yXl59z/+QPnMfv3yDtWNK6eQEAAAA3QxmDq5r7xCTnjFjXYQ0b13WYdMuzzu93LHaeBwAAAECdKGNwdXCL8x4wnyDp2oTGje032Tmu8Efp4GctkQ4AAABoFyhjcLV7ufNrv8nOSxAbw8vPOU6Sdic1by4AAACgHaGMwdWJw86v3X7TtPE14woON0caAAAAoF2ijMFVebHzq3dA08Z7+zu/nilunjwAAABAO0QZgyuvmjJV1LTxNSWsppQBAAAAcEEZg6vgSOfXH75o2viacRdFNkcaAAAAoF2ijMFVvynOr7vflcpLGze2vPRfFgCZ2ry5AAAAgHaEMgZX3WOloAjpdKG0Z0Xjxu5+Vzp9Ugq6Qup+c4vEAwAAANoDyhhceVilATOd33/yZyl7W8PG/fCl9OnTzu8HzHSeBwAAAECdKGOo24CZUq/RUlW59O546eu36r9ksbxU2vGmtGKC8/het/9a5gAAAADUyeJwOBymQ6Dh+vTpI0nat29fy/+wijJp3Qxp/wbnn32CnA907vYb50qJZ4qlQ59Le365NFFyFrFxb0o235bPBwAAABh2IZ/PPZs7DNoRm690x/9KO5ZIO96QCnOktNedr38XdIVzNmzATMmDCVcAAADgfChjODcPqzTwYWnAg9LBz6TdSVLBYeesmLe/c/n6flOdi3VwjxgAAADQYJQxNIyHVbryFucLAAAAwAXjejIAAAAAMIAyBgAAAAAGUMYAAAAAwADKGAAAAAAYQBkDAAAAAAMoYwAAAABgAGUMAAAAAAygjAEAAACAAZQxAAAAADCAMgYAAAAABlDGAAAAAMAAyhgAAAAAGEAZAwAAAAADKGMAAAAAYABlDAAAAAAMoIwBAAAAgAGUMQAAAAAwgDIGAAAAAAZQxgAAAADAAMoYAAAAABhAGQMAAAAAAyhjAAAAAGAAZQwAAAAADKCMAQAAAIABlDEAAAAAMIAyBgAAAAAGUMYAAAAAwADKGAAAAAAYQBkDAAAAAAMoYwAAAABgAGUMAAAAAAygjAEAAACAAZQxAAAAADCAMgYAAAAABlDGAAAAAMAAyhgAAAAAGEAZAwAAAAADKGMAAAAAYABlDAAAAAAMoIwBAAAAgAGUMQAAAAAwgDIGAAAAAAZQxgAAAADAAMoYAAAAABhAGQMAAAAAAyhjAAAAAGAAZQwAAAAADKCMAQAAAIABlDEAAAAAMIAyBgAAAAAGUMYAAAAAwADKGAAAAAAYQBkDAAAAAAMoYwAAAABgAGUMAAAAAAygjAEAAACAAZQxAAAAADCAMgYAAAAABlDGAAAAAMAAyhgAAAAAGEAZAwAAAAADKGMAAAAAYABlDAAAAAAMoIwBAAAAgAGUMQAAAAAwgDIGAAAAAAZQxgAAAADAAMoYAAAAABjgaToA3IOjqkrF27fr5AfrVJ6bq+rSUnn4+ckrPFydxo+T/9ChslitpmMCAAAAboMyhnNyVFXpxPLlKlj+riry8132n8nMVNHmzbKFhuqiKZMVPHWqLB5MuAIAAADnQxlDvapPn9aRx59Q0ebNkiSPTp0UNG6c7IMGysPur+qSYpWkpKowOVkV+fn6ef4LKtv1jUISX5SHj4/h9AAAAEDbZnE4HA7TIdBwffr0kSTt27evRX+Oo6pK+bNmqWjzFllsNl365FwFjR9fZ8mqLitT4bp1+nn+C3JUVCggLk6hryzkskUAAAC0exfy+ZzryVCnE8uX1xax8LfeUnBCQr2zXR6+vgpOSFD4W2/JYrOpaPNmFbz7bisnBgAAANxLhy1jJSUlWr58uf74xz8qJiZG3t7eslgsmj9//gWdd+PGjRo+fLg6deqkwMBADR8+XBs3bmym1K3DUVWlguXOMnXpk3Nlv3FAg8bZbxygS+fOlSSdWP6uHFVVLZYRAAAAcHcdtoxlZWVp6tSpev3117Vz506Vl5df8DlfffVVjR49WqmpqRo0aJBuuukm7dy5U6NHj9arr77aDKlbR/H27arIz6+9R6wxgsaPk0enTqrIy1PJV1+1UEIAAADA/XXYMhYQEKD77rtPS5Ys0TfffKOnn376gs73/fff67HHHpO3t7e2bdumv//971q/fr327Nmjiy++WI899piysrKaKX3LOvnBOklS0Lhx8vD1bdRYD19fBcXHS5IK137Q7NkAAACA9qLDlrGoqCgtXbpUDzzwgPr16ydPzwtbWHLRokWqrKzUzJkzNXDgwNrtV155pZ5++mlVVla6zexYeW6uJMk+aOB5jqybffAg53ny8potEwAAANDedNgy1txq7gubMGGCy7477rhDkrRhw4ZWzdRU1aWlkiQPu3+TxnvY7c7zlJQ0WyYAAACgvaGMNYPCwkLl5ORIkvr16+eyPywsTJ07d9aPP/6okydPtna8RvPw85MkVZcUN2l8TQmrKWUAAAAAXPHQ52ZQU8Quuugi2espIGFhYTp+/LhycnJ09dVXn/ecNc8r+HeHDh1SVFRU08M2gFd4uM5kZqokNU3+Q4c2enxJSqrzPGFhzR0NAAAAaDeYGWsGxcXOGSS/X2aU6lJT0mqObcs6jXeuoFi4bp2qy8oaNbbmAdCSFDRhfLNnAwAAANoLt50ZmzBhgvbu3duoMUlJSYqJiWn2LA6HQ5JksVjOe0xD1fcE7/pmzJqT/9ChsoWGqiI/X4Xr1ik4IaHBYws/WKfqU6dkCwuTfciQFkwJAAAAuDe3LWOHDx/WgQMHGjWm9JeFKZpbQECAJOeDpM/3s/39m7YoRmuyWK26aMpk/Tz/Bf08/wV5R3Vv0IOfS9LT9fMLL0iSgqdMlsVqbemoAAAAgNty28sUMzIy5HA4GvUaMWJEi2SJiIiQJBUUFNRbyPJ+Wea95ti2LnjKFAXExclRUaHcGTN0YsWKei9ZrC4r04l3Vyh3xgNyVFQo4JZbdNGUKa2cGAAAAHAvbjsz1pYEBQUpIiJCOTk52r17t4b82+V5eXl5On78uCIiItSpUydDKRvHYrUqJPFFHXn8CRVt3qx//GWejr36moLi42UfPEgedruqS0pU8lWKCpOTVX3qlCQp4JZbFPLiC7J4uG3PBwAAAFoFZayZjBo1Sm+88YbWrl3rUsbWrFkjSbrttttMRGsyDx8fhb6yUAXvvqsTSctVkZ+vE8uW6cSyZS7H2sLCFDxlsi6aMoUiBgAAADQAn5obKTo6WtHR0crPzz9r+6OPPiqr1arFixcrPT29dntWVpaeffZZWa1WPfLII60d94JZrFYFT5umqE8/UfiSxQqIi5N3r16yRUTIu1cvBcTFKXzJYkV9sknB06ZRxAAAAIAGsjgau8xfOxIfH6+jR49Kcl5KmJ+fr/DwcIWEhEiSLr/8ciUnJ581pmbFxOzsbEVGRp61b+HChZozZ448PT0VFxcnLy8vffrppyorK9PLL7+s2bNnX3DmmtUU61ttEQAAAEDruZDP5x36MsXdu3frxx9/PGtbbm6ucnNzJUlXXHFFo843e/Zsde/eXYmJidq+fbsk6frrr9fjjz+u22+/vXlCAwAAAGgXOvTMmDtiZgwAAABoOy7k8zk3+AAAAACAAZQxAAAAADCAMgYAAAAABlDGAAAAAMAAyhgAAAAAGEAZAwAAAAADKGMAAAAAYECHfugzmqCyUlq/3vn92LGSJ28hAAAAoCn4JI3G2bpVuuMO5/ebN0uxsUbjAAAAAO6KyxRxbg6HVFLy62v16l/3rV599j6Hw1xOAAAAwM1QxnBuaWmSv/+vr7ff/nXf22+fvS893VxOAAAAwM1QxnBuSUktcywAAADQwVHGcG5//rPUrdv5j+vWzXksAAAAgAahjOHcIiKkL744dyHr1s25sEd4eKvFAgAAANwdZQznFxEhJSbWvz8xkSIGAAAANBJlDA2zdm39+z74oPVyAAAAAO0EZQwNs2rVr98//LD00EO//nnlytbPAwAAALg5yhga5vBhaehQ6dNPpb/+Vfqf/5E++UQaNkz68UfT6QAAAAC342k6ANzEFVdI27adve2WW5wvAAAAAI3GzBgAAAAAGEAZAwAAAAADKGMAAAAAYABlDAAAAAAMoIwBAAAAgAGUMQAAAAAwgDIGAAAAAAZQxgAAAADAAMoYAAAAABhAGQMAAAAAAyhjAAAAAGAAZQwAAAAADKCMAQAAAIABlDEAAAAAMIAyBgAAAAAGUMYAAAAAwADKGAAAAAAYQBkDAAAAAAM8TQeAe6iqrlLKkRQlZyUrtyhXpZWl8vP0U3hAuOJ7xGtwyGBZPaymYwIAAABugzKGc6qqrtKK/Su0MnOl8ovzXfYfKDigLTlbFOofqknRkzS592R5WJhwBQAAAM6HMoZ6na48rae2P6UtOVskSYFegYrvHq8bQ26Uv81fxRXFSjuSpvUH1yu/OF+JGYna/fNuPT/0efl4+hhODwAAALRtlDHUqaq6Sk9uf1Kf5Xwmm4dNj/d/XPHd411K1pDQIfpDvz9o/cH1StyZ6Cxu26UFwxdw2SIAAABwDlxPhjqt2L+itogtjl2su6Pvrne2y9fTV3dH363FsYtl87BpS84Wrcxc2cqJAQAAAPdCGYOLquqq2jL1eP/HFXN5TIPGxVweoz/d8CdJzjJXVV3VYhkBAAAAd0cZg4uUIynKL85XoFegxnYf26ix8T3iFegVqPzifKUcSWmZgAAAAEA7QBmDi+SsZElSfPd4+Xr6Nmqsr6dvbYGrOQ8AAAAAV5QxuMgtypUk3RhyY5PGDwwZKEnKK85rtkwAAABAe0MZg4vSylJJkr/Nv0nj7Ta7JKmkoqTZMgEAAADtDWUMLvw8/SRJxRXFTRpfU8JqShkAAAAAV5QxuAgPCJckpR9Jb9L4tCNpkqQw/7BmywQAAAC0N5QxuIjvES9JSj6YrLLKskaNLassU/LB5LPOAwAAAMAVZQwuBocMVqh/qE6Vn9L6g+sbNTY5K1lF5UUK9Q/V4JDBLRMQAAAAaAcoY3Bh9bBqUvQkSVLizkR9ffTrBo3bcXSHFmQskCQl9EqQ1cPaYhkBAAAAd0cZQ50SeiUoNiJWFdUVmrllplZlrqr3ksWyyjKt3L9SD215SBXVFYq7Ik4JvRJaOTEAAADgXiwOh8NhOgQark+fPpKkffv2tfjPOl15Wk9tf0pbcrZIkgK9AjW2+1gNDBkou82ukooSpR5J1fqD61VUXiRJirsiTs8NeU4+nj4tng8AAAAw7UI+n1PG3ExrljFJqqqu0srMlVqxf4Xyi/PrPS7UP1QJvRKU0CtBHhYmXAEAANAxXMjnc8/mDoP2xeph1ZTeUzQpepJSjqQoOStZecV5Kqkokd1mV5h/mOJ7xGtwyGDuEQMAAAAagTKGBrF6WDUsbJiGhQ0zHQUAAABoF7ieDAAAAAAMoIwBAAAAgAGUMQAAAAAwgDIGAAAAAAZQxgAAAADAAMoYAAAAABhAGQMAAAAAAyhjAAAAAGAAZQwAAAAADKCMAQAAAIABlDEAAAAAMIAyBgAAAAAGUMYAAAAAwACLw+FwmA6BhgsICFBFRYWioqJMRwEAAAA6vEOHDslms6moqKjRY5kZczN2u102m810DDTCoUOHdOjQIdMx0M7wvkJL4H2FlsD7Ci2lrby3bDab7HZ7k8YyMwa0sD59+kiS9u3bZzgJ2hPeV2gJvK/QEnhfoaW0h/cWM2MAAAAAYABlDAAAAAAMoIwBAAAAgAGUMQAAAAAwgDIGAAAAAAawmiIAAAAAGMDMGAAAAAAYQBkDAAAAAAMoYwAAAABgAGUMAAAAAAygjAEAAACAAZQxAAAAADCAMgYAAAAABlDGgGZUWlqq9evX67777lPfvn0VGBgou92ua665Rv/93/+t4uJi0xHRTpw4cUKXXnqpLBaLoqOjTcdBO/DTTz9p9uzZuvLKK+Xr66vg4GBdf/31euKJJ0xHg5tKT0/X+PHj1aVLF9lsNgUHB+vmm2/W2rVrTUdDG7dr1y7Nnz9f48aNU2hoqCwWi3x8fM47LikpSTExMfL391dwcLB+97vfKTU1tRUSNx0PfQaa0dKlSzVjxgxJUp8+fdS7d2+dOnVKqampKioqUnR0tL788ktdeumlhpPC3U2fPl1JSUlyOBzq2bOnMjMzTUeCG0tLS9Pvfvc7FRYWqnfv3rrqqqtUVFSk7777Tnl5eaqsrDQdEW5mzZo1uuuuu1RdXa0bbrhBUVFROnLkiFJSUlRdXa25c+dq/vz5pmOijRo7dqw+/PDDs7Z5e3vr9OnT9Y6ZM2eOFi5cKF9fX91yyy06ffq0PvvsMzkcDq1Zs0bx8fEtHbtJKGNAM0pKSlJ6erpmz56tHj161G4/evSoRo0apd27d+vuu+/WypUrDaaEu/vss88UGxurBx54QG+++SZlDBfkyJEj6tOnj86cOaMVK1a4fGD5+uuvFRMTYygd3FFlZaVCQkJ07NgxrV69WnfeeWftvrS0NN100006c+aMsrKyFBUVZTAp2qoXXnhBpaWl6t+/v/r3768uXbqcs4x9/vnnuvnmm3XxxRcrLS2t9jNYWlqaRowYIV9fX2VnZ+uiiy5qzb9Gg1DGgFaSlpamQYMGydvbW6dOnZKXl5fpSHBDZWVl6tu3r7y8vLR+/XpdeeWVlDFckKlTp2r58uV67bXX9Ic//MF0HLQDe/fu1dVXX63o6Gjt37/fZX/NrMd7772niRMnGkgId2OxWM5ZxkaNGqWPP/5YCxcu1KxZs87a9+ijj+rVV1/VggUL9Nhjj7VC2sbhnjGglVxzzTWSpDNnzuif//yn4TRwV//1X/+lQ4cO6Y033pDNZjMdB26uoKBA77//vjp16qT777/fdBy0E97e3g06Ljg4uIWToCOouRxRkiZMmOCyv2bbhg0bWjVXQ3maDgB0FD/88IMk1d7EDDTWt99+q5deekn33HOPhg0bpsOHD5uOBDeXkpKiM2fOKDY2VjabTWvXrtVXX32liooKRUdHa+LEibrssstMx4Sb6datm7p166bMzEy9//77Z81+paWl6ZNPPlHXrl01bNgwgynRXmRmZurMmTO65JJLFBYW5rL/uuuuk+T8N7QtoowBrWTRokWSpJEjRzb4t4ZAjerqas2YMUNBQUF68cUXTcdBO7Fv3z5J0mWXXaahQ4cqLS3trP1PPfWU3nnnHd1xxx0m4sFNWa1WLVu2TKNHj9add96pxMRERUVF6ejRo/rqq68UExOj5cuXc7k+mkVOTo4k1VnEJMlutysoKEgFBQUqKipSQEBAa8Y7Ly5TBFrBxx9/rLfffls2m01/+ctfTMeBG3rttdf09ddfKzExURdffLHpOGgnCgoKJDkXH/r222/19ttv69ixY8rOztacOXNUUlKiyZMnt9nfKKPtGjp0qL788kt17dpVGRkZeu+997Rt2zbZ7XbFxsYqJCTEdES0EzWPDfLz86v3GLvdftaxbQllDGhh+/fv1+TJk+VwOJSYmFh77xjQULm5uXrmmWc0fPhwTZ8+3XQctCNVVVWSnKvfvfzyy7r33nvVuXNnRUZG6qWXXtKECRNUXl7ObCwabdWqVRowYIAiIiK0Y8cOFRcX6/vvv9fdd9+tefPmKTY2VhUVFaZjoh2oWYvQYrGc95i2iDIGtKC8vDyNHDlSBQUFmjNnjh599FHTkeCGHn74YZWXl+uNN94wHQXtTM3lOh4eHpo2bZrL/nvvvVeStHXr1taMBTeXlZWladOm6ZJLLtFHH32kmJgY2e129ejRQ0uWLNHo0aOVlpamd955x3RUtAM1/x8rKSmp95jS0lJJkr+/f6tkagzuGQNayPHjxxUXF6ecnBzdc889WrBggelIcFMbN25UUFCQHnroobO21yzxm5OToxEjRtQe2xb/sUHbFBkZKUm1z/Cpb//PP//ciqng7lavXq2KigqNHDmy9vKwfzVx4kRt2LBBW7du1QMPPGAgIdqTiIgISc5fgNelpKREhYWFCgoKanP3i0mUMaBFFBUV6be//a0yMzM1btw4vfXWW+ecPgfOp7CwUF9++WWd+8rKymr3VVZWtmYsuLl+/fpJct475nA4XP4/VfMYDgo+GqPmQ3FgYGCd+2u2nzhxotUyof3q2bOnvL29dezYMeXl5bks5PHNN99Ikvr27Wsi3nlxmSLQzM6cOaMxY8YoIyNDt956q1atWiWr1Wo6FtyYw+Go85WdnS3J+Q9RzbagoCCzYeFWrr76anXt2lVlZWXasWOHy/6ayxNrloYGGqJLly6SpIyMjDr379y5U9KvM6/AhfD19dVNN90kSVq7dq3L/pptt912W6vmaijKGNCMqqqqdPfdd+uLL77Q0KFDtW7dOpbuBdCmzZ07V5L0yCOP6Pjx47Xbd+3apZdeekmSNHPmTCPZ4J7GjBkjSdq2bZvLva7p6elauHChpLof0As0xZw5cyRJ8+bNU1ZWVu32tLQ0LVmyRIGBgbrvvvtMxTsni6MtLy8CuJlFixZp1qxZkqT4+Ph6L9FYsGCBOnfu3IrJ0B4dPnxYXbt2Vc+ePZWZmWk6DtxUdXW17rrrLq1Zs0bBwcEaNGiQiouLlZqaqvLycs2YMUNvvvmm6ZhwM48//njtvdJ9+vRR7969deTIEaWlpam6uloPPPCAlixZYjgl2qqPPvrorEcB7dixQxaLRTExMbXb/uM//kOjRo2q/fOsWbO0aNEi+fn5KS4uTuXl5dq8ebOqq6v1/vvva/z48a36d2go7hkDmlHNM3skKTk5ud7j/vM//5MyBqBN8PDw0OrVqzVixAgtXbpUn3/+uSwWi2644QbNnDlTU6ZMMR0RbigxMVGDBg3S4sWLtWvXLh04cEABAQEaPny47r//fk2aNMl0RLRhx44dc7l02uFwnLXt2LFjZ+1/5ZVXdO211+r111/X5s2bZbPZdPPNN+uZZ57RkCFDWiV3UzAzBgAAAAAGcM8YAAAAABhAGQMAAAAAAyhjAAAAAGAAZQwAAAAADKCMAQAAAIABlDEAAAAAMIAyBgAAAAAGUMYAAAAAwADKGAAAAAAYQBkDAAAAAAMoYwAAAABgAGUMAAAAAAygjAEA0MK2bt0qi8WiyMjIcx43YsQIWSwWLVu2rFVyAQDMoowBAAAAgAGUMQAAAAAwgDIGAAAAAAZQxgAAAADAAMoYAAAAABhAGQMAAAAAAyhjAAAAAGAAZQwAAAAADKCMAQDQwiwWS4OOczgcjToeAODeKGMAALQwX19fSVJJSck5jystLZUk2e32Fs8EADCPMgYAQAsLDw+XJB0/flynTp2q97gffvhBkhQWFtYquQAAZlHGAABoYZdffrm6d+8uSdq4cWOdx6SkpOjEiRPy9/dXv379WjMeAMAQyhgAAK3g0UcflSTNnTtXmZmZZ+07evSoHn74YUnSzJkz5e3t3er5AACtz+KouVsYAAC0GIfDoUmTJmn16tWy2WwaPHiwQkNDdezYMW3fvl1lZWUaPny4Nm3aJB8fH9NxAQCtgDIGAEArcTgceu+99/TOO+/om2++UWFhoQICAnTVVVdp0qRJuu+++2Sz2UzHBAC0EsoYAAAAABjAPWMAAAAAYABlDAAAAAAMoIwBAAAAgAGUMQAAAAAwgDIGAAAAAAZQxgAAAADAAMoYAAAAABhAGQMAAAAAAyhjAAAAAGAAZQwAAAAADKCMAQAAAIABlDEAAAAAMIAyBgAAAAAGUMYAAAAAwADKGAAAAAAYQBkDAAAAAAMoYwAAAABgwP8HzCxR3FDmrfYAAAAASUVORK5CYII=", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "# plotting\n", "plt.figure(dpi=150)\n", @@ -373,7 +363,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "80867ad4", "metadata": {}, "outputs": [], @@ -417,7 +407,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.10.11" } }, "nbformat": 4, diff --git a/examples/Ising.ipynb b/examples/Ising.ipynb index c045b7b..9b90179 100644 --- a/examples/Ising.ipynb +++ b/examples/Ising.ipynb @@ -23,13 +23,13 @@ "source": [ "## Defining Heisenberg Model\n", "There are two ways to define the Heisenberg model using the Model Hamiltonian Package:\n", - "1. Using the connectivity of the lattice and providing the exchange couplings, $J^{ax}$, $J^{eq}$ , and $\\mu$ as a float numbers;\n", + "1. Using the adjacency of the lattice and providing the exchange couplings, $J^{ax}$, $J^{eq}$ , and $\\mu$ as a float numbers;\n", "2. Providing exchange couplings, $J^{ax}$, $J^{eq}$ , and $\\mu$ explicitly as numpy arrays;\n" ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -37,11 +37,11 @@ "import numpy as np\n", "from moha import HamHeisenberg\n", "\n", - "# Defining the Hamiltonian using connectivity matrix\n", + "# Defining the Hamiltonian using adjacency matrix\n", "# For example, consider 6x6 lattice with periodic boundary condition\n", "\n", - "# Define the connectivity matrix\n", - "connectivity = np.array([[0, 1, 0, 0, 0, 1],\n", + "# Define the adjacency matrix\n", + "adjacency = np.array([[0, 1, 0, 0, 0, 1],\n", " [1, 0, 1, 0, 0, 0],\n", " [0, 1, 0, 1, 0, 0],\n", " [0, 0, 1, 0, 1, 0],\n", @@ -54,12 +54,12 @@ "mu = 0\n", "\n", "# Define the Hamiltonian\n", - "ham = HamHeisenberg(connectivity=connectivity, J_eq=J_eq, J_ax=J_ax, mu=mu)" + "ham = HamHeisenberg(adjacency=adjacency, J_eq=J_eq, J_ax=J_ax, mu=mu)" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -110,7 +110,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -190,14 +190,14 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Ground state energy per electron: -0.3564972638984237\n" + "Ground state energy per electron: -0.35649726389842346\n" ] } ], @@ -277,14 +277,18 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 11, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "Ground state energy per electron: -0.35649726389842346\n" + "ename": "ModuleNotFoundError", + "evalue": "No module named 'pyscf'", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[11], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mpyscf\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m fci\n\u001b[0;32m 3\u001b[0m \u001b[38;5;66;03m# Constructing the XXZ Hamiltonian using the HamHeisenberg class\u001b[39;00m\n\u001b[0;32m 4\u001b[0m \n\u001b[0;32m 5\u001b[0m \u001b[38;5;66;03m# Define the Hamiltonian parameters\u001b[39;00m\n\u001b[0;32m 6\u001b[0m adjacency \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mzeros((L, L))\n", + "\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'pyscf'" ] } ], @@ -294,13 +298,13 @@ "# Constructing the XXZ Hamiltonian using the HamHeisenberg class\n", "\n", "# Define the Hamiltonian parameters\n", - "connectivity = np.zeros((L, L))\n", - "connectivity[np.arange(L-1), np.arange(L-1)+1] = 1\n", - "connectivity[0, -1] = 1\n", - "connectivity += connectivity.T\n", + "adjacency = np.zeros((L, L))\n", + "adjacency[np.arange(L-1), np.arange(L-1)+1] = 1\n", + "adjacency[0, -1] = 1\n", + "adjacency += adjacency.T\n", "\n", "# Define the Hamiltonian\n", - "ham = HamHeisenberg(connectivity=connectivity, J_eq=1, J_ax=0.2, mu=0)\n", + "ham = HamHeisenberg(adjacency=adjacency, J_eq=1, J_ax=0.2, mu=0)\n", "e0 = ham.generate_zero_body_integral()\n", "h1 = ham.generate_one_body_integral(dense=True, basis='spatial basis')\n", "h2 = ham.generate_two_body_integral(dense=True, basis='spatial basis', sym=4)\n", @@ -328,7 +332,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [ { @@ -350,14 +354,14 @@ "fig, ax = plt.subplots(1, 2, figsize=(12, 6))\n", "\n", "for L in range(4, 10, 2):\n", - " # Define the connectivity matrix\n", - " connectivity = np.zeros((L, L))\n", - " connectivity[np.arange(L-1), np.arange(L-1)+1] = 1\n", - " connectivity[0, -1] = 1\n", - " connectivity += connectivity.T\n", + " # Define the adjacency matrix\n", + " adjacency = np.zeros((L, L))\n", + " adjacency[np.arange(L-1), np.arange(L-1)+1] = 1\n", + " adjacency[0, -1] = 1\n", + " adjacency += adjacency.T\n", "\n", " # Define the Hamiltonian\n", - " ham = HamHeisenberg(connectivity=connectivity, J_eq=1, J_ax=1, mu=0)\n", + " ham = HamHeisenberg(adjacency=adjacency, J_eq=1, J_ax=1, mu=0)\n", " e0 = ham.generate_zero_body_integral()\n", " h1 = ham.generate_one_body_integral(dense=True, basis='spatial basis')\n", " h2 = ham.generate_two_body_integral(dense=True, basis='spatial basis', sym=4)\n", @@ -411,7 +415,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -419,13 +423,13 @@ "\n", "# Define the Hamiltonian parameters\n", "L = 4\n", - "connectivity = np.zeros((L, L))\n", - "connectivity[np.arange(L-1), np.arange(L-1)+1] = 1\n", - "connectivity[0, -1] = 1\n", - "connectivity += connectivity.T\n", + "adjacency = np.zeros((L, L))\n", + "adjacency[np.arange(L-1), np.arange(L-1)+1] = 1\n", + "adjacency[0, -1] = 1\n", + "adjacency += adjacency.T\n", "\n", "# Define the Hamiltonian\n", - "ham = HamHeisenberg(connectivity=connectivity, J_eq=1, J_ax=0.2, mu=0)\n", + "ham = HamHeisenberg(adjacency=adjacency, J_eq=1, J_ax=0.2, mu=0)\n", "e0 = ham.generate_zero_body_integral()\n", "h1 = ham.generate_one_body_integral(dense=True, basis='spatial basis')\n", "h2 = ham.generate_two_body_integral(dense=True, basis='spatial basis', sym=4)\n", @@ -453,7 +457,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.5" + "version": "3.10.11" } }, "nbformat": 4, diff --git a/examples/Real_examples1.ipynb b/examples/Real_examples1.ipynb new file mode 100644 index 0000000..a72f64b --- /dev/null +++ b/examples/Real_examples1.ipynb @@ -0,0 +1,423 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tutorial: Benzene and Biphenyl\n", + "\n", + "## Introduction\n", + "In this tutorial, we'll focus on generating integrals for the PPP model considering real organic molecules: Benzene and Biphenyl " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To do this, we are going to consider the most general occupation-number Hamiltonian we consider is the generalized Pariser-Parr-Pople (PPP) Hamiltonian:\n", + "\n", + "$$\\hat{H}_{\\text{PPP}} = \\sum_{pq} h_{pq} a_p^\\dagger a_q + \\sum_p U_p \\hat{n}_{p\\alpha}\\hat{n}_{p\\beta} + \\frac{1}{2}\\sum_{p\\ne q} \\gamma_{pq} (\\hat{n}_{p \\alpha} + \\hat{n}_{p \\beta} - Q_p)(\\hat{n}_{q \\alpha} + \\hat{n}_{q \\beta} - Q_q) $$\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example: Defining PPP Hamiltonian\n", + "\n", + "In order to define a PPP hamiltonina, it's necessary to provide the parameter `gamma`. Performing a Full Configuration Iteraction of the hamiltonian, varying $\\beta$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "65.30151\n", + "Beta: 0.0, FCI Energy: -1.4210854715202004e-14\n", + "65.30151\n", + "Beta: -1.0, FCI Energy: -2.779514577625193\n", + "65.30151\n", + "Beta: -2.0, FCI Energy: -9.084603708594514\n", + "65.30151\n", + "Beta: -3.0, FCI Energy: -16.484889549995387\n", + "65.30151\n", + "Beta: -4.0, FCI Energy: -24.19501406479904\n", + "65.30151\n", + "Beta: -5.0, FCI Energy: -32.024692914537425\n" + ] + } + ], + "source": [ + "import sys \n", + "sys.path.insert(0, '../')\n", + "import numpy as np \n", + "import sys \n", + "import numpy as np \n", + "from moha import HamPPP \n", + "from pyscf import gto, scf, fci \n", + "\n", + "RXY = np.array([0.00504725, 0.09930802, 0.16008021])\n", + " \n", + "# Define the benzene system with a 6-membered ring structure \n", + "system = [('C1', 'C2', 1), ('C2', 'C3', 1), \n", + " ('C3', 'C4', 1), ('C4', 'C5', 1), \n", + " ('C5', 'C6', 1), ('C6', 'C1', 1)] \n", + "\n", + "\n", + "#Define U_onsite, alpha and gamma matrix \n", + "u_onsite = np.array([10.84, 10.84, 10.84, 10.84, 10.84, 10.84])/2\n", + "alpha = 0 \n", + " \n", + "# Define the gamma matrix here (example shape) \n", + "norb = 6 # Number of orbitals for benzene (6 carbons) \n", + "gamma_matrix = np.zeros((norb, norb))\n", + "\n", + "gamma_a_a1 = 5.27760 \n", + "gamma_a_a2 = 3.86206 \n", + "gamma_a_a3 = 3.48785 \n", + " \n", + "# Fill the matrix according to the given rules \n", + "for a in range(norb): \n", + " gamma_matrix[a-1, a] = gamma_a_a1 \n", + " gamma_matrix[a-2, a] = gamma_a_a2 \n", + " gamma_matrix[a-3, a] = gamma_a_a3 \n", + " \n", + " gamma_matrix[a, a-1] = gamma_a_a1 \n", + " gamma_matrix[a, a-2] = gamma_a_a2 \n", + " gamma_matrix[a, a-3] = gamma_a_a3\n", + "\n", + " gamma_matrix[a, a] = 10.84\n", + " \n", + "gamma_lib = gamma_matrix - np.diag(np.diag(gamma_matrix))\n", + "# Loop to vary beta from 0 to -5\n", + "for beta in np.linspace(0, -5, num=6):\n", + " # Generate the PPP object\n", + " ppp_hamiltonian = HamPPP(system, alpha=alpha, beta=beta, u_onsite=u_onsite, gamma=gamma_lib, charges=np.array([1 for _ in range(6)]))\n", + " \n", + " # Generate the 0, 1, 2 body integral elements\n", + " e01 = ppp_hamiltonian.generate_zero_body_integral() \n", + " h11 = ppp_hamiltonian.generate_one_body_integral(dense=True, basis='spatial basis') \n", + " h2 = ppp_hamiltonian.generate_two_body_integral(dense=True, basis='spatial basis',sym = 4) \n", + " print(e01)\n", + " \n", + " # Convert h2 to chemists notation \n", + " h21_ch = 1*np.transpose(h2, (0, 2, 1, 3)) \n", + " norb = h11.shape[0] \n", + " nelec = norb # Assuming one electron per site \n", + "\n", + " # Perform the FCI calculation \n", + " e_fci, ci = fci.direct_spin1.kernel(1*h11, 2*h21_ch, norb, nelec, ecore=e01, nroots=1)\n", + " print(f'Beta: {beta}, FCI Energy: {e_fci}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# PPP Benzene: Building gamma\n", + "\n", + "\n", + "MoHa package offers a functionality to compute `gamma` using the module `PariserParr` and the function `compute_gamma`. \n", + "\n", + "In order to do that, user can provide all the paramters of the function or let the library compute with the default values:\n", + "\n", + "* If user does not provide `U_xy`(potential energies for sites) or `R_xy`matrix with distances of the sites; Those values are going to be computed using the distances provided by `connectivity`, using the formulas below, based on the deafult values of `affinity_dct` and `atom_dict`:\n", + "\n", + "$$ U_x = \\alpha_x - A_x $$\n", + "\n", + "$$\n", + "\\overline{U}_{XY} = \\frac{1}{2}(U_X + U_Y)\n", + "$$\n", + "\n", + "User can also provide this dictionaries(affinity_dct` and `atom_dict`) to the function and proceed in the same way.\n", + "\n", + "\n", + "Finally, `gamma` is computed using the formula:\n", + "\n", + "$$\n", + "\\gamma_{XY} = \\frac{\\overline{U}_{XY}}{\\overline{U}_{XY} R_{XY} + e^{-\\frac{1}{2} \\overline{U}_{XY}^2 R_{XY}^2}}\n", + "$$\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Defining the Benzene system following th same sttructure defined in the previous tutorials, we decide to move on with Full Configuration Iteraction of the hamiltonian, defining `gamma`\n", + "\n", + "Varying $\\beta$ from 0 to -5, in this piece of code we prove that we can achieve the spectrum of energies reported in Bendazzoli et al." + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "65.30151\n", + "Beta: 0.0, FCI Energy: -1.4210854715202004e-14\n", + "65.30151\n", + "Beta: -1.0, FCI Energy: -2.779514577625193\n", + "65.30151\n", + "Beta: -2.0, FCI Energy: -9.084603708594514\n", + "65.30151\n", + "Beta: -3.0, FCI Energy: -16.484889549995387\n", + "65.30151\n", + "Beta: -4.0, FCI Energy: -24.19501406479904\n", + "65.30151\n", + "Beta: -5.0, FCI Energy: -32.024692914537425\n" + ] + } + ], + "source": [ + "import sys \n", + "sys.path.insert(0, '../')\n", + "import numpy as np \n", + "import sys \n", + "import numpy as np \n", + "from moha import HamPPP \n", + "from pyscf import gto, scf, fci \n", + "\n", + "RXY = np.array([0.00504725, 0.09930802, 0.16008021])\n", + " \n", + "# Define the benzene system with a 6-membered ring structure \n", + "system = [('C1', 'C2', 1), ('C2', 'C3', 1), \n", + " ('C3', 'C4', 1), ('C4', 'C5', 1), \n", + " ('C5', 'C6', 1), ('C6', 'C1', 1)] \n", + "\n", + "\n", + "#Define U_onsite, alpha and gamma matrix \n", + "u_onsite = np.array([10.84, 10.84, 10.84, 10.84, 10.84, 10.84])/2\n", + "alpha = 0 \n", + " \n", + "# Define the gamma matrix here (example shape) \n", + "norb = 6 # Number of orbitals for benzene (6 carbons) \n", + "gamma_matrix = np.zeros((norb, norb))\n", + "\n", + "gamma_a_a1 = 5.27760 \n", + "gamma_a_a2 = 3.86206 \n", + "gamma_a_a3 = 3.48785 \n", + " \n", + "# Fill the matrix according to the given rules \n", + "for a in range(norb): \n", + " gamma_matrix[a-1, a] = gamma_a_a1 \n", + " gamma_matrix[a-2, a] = gamma_a_a2 \n", + " gamma_matrix[a-3, a] = gamma_a_a3 \n", + " \n", + " gamma_matrix[a, a-1] = gamma_a_a1 \n", + " gamma_matrix[a, a-2] = gamma_a_a2 \n", + " gamma_matrix[a, a-3] = gamma_a_a3\n", + "\n", + " gamma_matrix[a, a] = 10.84\n", + " \n", + "gamma_lib = gamma_matrix - np.diag(np.diag(gamma_matrix))\n", + "# Loop to vary beta from 0 to -5\n", + "for beta in np.linspace(0, -5, num=6):\n", + " # Generate the PPP object\n", + " ppp_hamiltonian = HamPPP(system, alpha=alpha, beta=beta, u_onsite=u_onsite, gamma=gamma_lib, charges=np.array([1 for _ in range(6)]))\n", + " \n", + " # Generate the 0, 1, 2 body integral elements\n", + " e01 = ppp_hamiltonian.generate_zero_body_integral() \n", + " h11 = ppp_hamiltonian.generate_one_body_integral(dense=True, basis='spatial basis') \n", + " h2 = ppp_hamiltonian.generate_two_body_integral(dense=True, basis='spatial basis',sym = 4) \n", + " print(e01)\n", + " \n", + " # Convert h2 to chemists notation \n", + " h21_ch = 1*np.transpose(h2, (0, 2, 1, 3)) \n", + " norb = h11.shape[0] \n", + " nelec = norb # Assuming one electron per site \n", + "\n", + " # Perform the FCI calculation \n", + " e_fci, ci = fci.direct_spin1.kernel(1*h11, 2*h21_ch, norb, nelec, ecore=e01, nroots=1)\n", + " print(f'Beta: {beta}, FCI Energy: {e_fci}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pariser Parr Atomic dictionary\n", + "\n", + "The previous approach requires to build the $\\gamma$ matrix directly, but MoHa package supports a method of computing the matrix based on the distance($R_{XY}$) beetween the atoms and the potential($U_{XY}$). \n", + "\n", + "\n", + "If the value of `gamma` is not provided, it is computed using the formula:\n", + "\n", + "$$\n", + "\\gamma_{XY} = \\frac{\\overline{U}_{XY}}{\\overline{U}_{XY} R_{XY} + e^{-\\frac{1}{2} \\overline{U}_{XY}^2 R_{XY}^2}}\n", + "$$\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Beta: 0.0, FCI Energy: 7.105427357601002e-15\n", + "Beta: -1.0, FCI Energy: -2.950976576402695\n", + "Beta: -2.0, FCI Energy: -9.425258557059209\n", + "Beta: -3.0, FCI Energy: -16.75428207213497\n", + "Beta: -4.0, FCI Energy: -24.36705165013518\n", + "Beta: -5.0, FCI Energy: -32.1115278385854\n" + ] + } + ], + "source": [ + "import sys\n", + "import numpy as np\n", + "from moha import HamPPP\n", + "from moha.rauk import PariserParr\n", + "from pyscf import gto, scf, fci\n", + "\n", + "# Define the benzene system with a 6-membered ring structure\n", + "\n", + "system = [('C1', 'C2', 1), ('C2', 'C3', 1), \n", + " ('C3', 'C4', 1), ('C4', 'C5', 1), \n", + " ('C5', 'C6', 1), ('C6', 'C1', 1)] \n", + "\n", + "# Define the distances\n", + "RXY = [0.00504725, 0.09930802, 0.16008021]\n", + "\n", + "# Build the adjacency matrix directly\n", + "Rxy_matrix = np.array([\n", + " [0.0, RXY[0], RXY[1], RXY[2], RXY[1], RXY[0]],\n", + " [RXY[0], 0.0, RXY[0], RXY[1], RXY[2], RXY[1]],\n", + " [RXY[1], RXY[0], 0.0, RXY[0], RXY[1], RXY[2]],\n", + " [RXY[2], RXY[1], RXY[0], 0.0, RXY[0], RXY[1]],\n", + " [RXY[1], RXY[2], RXY[1], RXY[0], 0.0, RXY[0]],\n", + " [RXY[0], RXY[1], RXY[2], RXY[1], RXY[0], 0.0]\n", + "])\n", + "\n", + "u_onsite = np.array([10.84, 10.84, 10.84, 10.84, 10.84, 10.84]) / 2\n", + "gamma_matrix = PariserParr.compute_gamma(system, U_xy=u_onsite, Rxy_matrix=Rxy_matrix)\n", + "\n", + "alpha = 0\n", + "\n", + "# Loop to vary beta from 0 to -5\n", + "for beta in np.linspace(0, -5, num=6):\n", + " ppp_hamiltonian = HamPPP(system, alpha=alpha, gamma = gamma_matrix, beta=beta, u_onsite=u_onsite, charges=np.array([1 for _ in range(6)]))\n", + " \n", + " # Generate the 0, 1, 2 body integral elements\n", + " e0 = ppp_hamiltonian.generate_zero_body_integral() \n", + " h1 = ppp_hamiltonian.generate_one_body_integral(dense=True, basis='spatial basis') \n", + " h2 = ppp_hamiltonian.generate_two_body_integral(dense=True, basis='spatial basis', sym=4)\n", + " \n", + " # Convert h2 to chemists notation\n", + " h2_ch = np.transpose(h2, (0, 2, 1, 3))\n", + " norb = h1.shape[0]\n", + " nelec = norb # Assuming one electron per site \n", + "\n", + " # Perform the FCI calculation\n", + " e_fci, ci = fci.direct_spin1.kernel(1*h1, 2*h2_ch, norb, nelec, ecore=e0, nroots=1)\n", + " print(f'Beta: {beta}, FCI Energy: {e_fci}')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Heteroatoms\n", + "\n", + "The ModelHamiotonian package also supports the possibility of using different type of atoms, using the Rauk’s table (RAUK, 2001).\n", + "\n", + "To do this, simply replace an atom of C with one of B.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Alpha: -0.414, Beta: -0.0533, FCI Energy: -5.468422237553814\n" + ] + } + ], + "source": [ + "import sys \n", + "import numpy as np \n", + "from moha import HamPPP \n", + "from pyscf import gto, scf, fci \n", + " \n", + "# Define the benzene system with a 6-membered ring structure \n", + "system = [('B1', 'C2', 1), ('C2', 'C3', 1), \n", + " ('C3', 'C4', 1), ('C4', 'C5', 1), \n", + " ('C5', 'C6', 1), ('C6', 'B1', 1)] \n", + " \n", + "u_onsite = np.array([10.84, 10.84, 10.84, 10.84, 10.84, 10.84]) / 2\n", + "\n", + "gamma_matrix = PariserParr.compute_gamma(system, U_xy=u_onsite)\n", + " \n", + "norb = 6 # Number of orbitals for benzene (6 carbons) \n", + "ppp_hamiltonian = HamPPP(system, u_onsite=u_onsite, gamma=gamma_matrix, charges=np.array([1 for _ in range(6)]))\n", + "\n", + "e0 = ppp_hamiltonian.generate_zero_body_integral() \n", + "h1 = ppp_hamiltonian.generate_one_body_integral(dense=True, basis='spatial basis') \n", + "h2 = ppp_hamiltonian.generate_two_body_integral(dense=True, basis='spatial basis', sym=4) \n", + " \n", + "# Convert h2 to chemists notation \n", + "h2_ch = np.transpose(h2, (0, 2, 1, 3)) \n", + "norb = h1.shape[0] \n", + "nelec = norb # Assuming one electron per site \n", + "\n", + "# Perform the FCI calculation \n", + "e_fci, ci = fci.direct_spin1.kernel(h1, h2_ch, norb, nelec, ecore=e0, nroots=1)\n", + "print(f'Alpha: {ppp_hamiltonian.alpha}, Beta: {ppp_hamiltonian.beta}, FCI Energy: {e_fci}')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# References\n", + "\n", + "Bendazzoli, G. L., Evangelisti, S., & Gagliardi, L. (1994). Full configuration interaction study of the ground state of closed-shell cyclic PPP polyenes. International Journal of Quantum Chemistry, 51(1), 117-123. https://doi.org/10.1002/qua.560510104\n", + "\n", + "\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/TJUV.ipynb b/examples/TJUV.ipynb new file mode 100644 index 0000000..8efe458 --- /dev/null +++ b/examples/TJUV.ipynb @@ -0,0 +1,224 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# TJUV model" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Zero energy: 1.5\n", + "One body integrals in spatial basis: \n", + " [[-0.5 -1. 0. 0. 0. -1. ]\n", + " [-1. -0.5 -1. 0. 0. 0. ]\n", + " [ 0. -1. -0.5 -1. 0. 0. ]\n", + " [ 0. 0. -1. -0.5 -1. 0. ]\n", + " [ 0. 0. 0. -1. -0.5 -1. ]\n", + " [-1. 0. 0. 0. -1. -0.5]]\n", + "Shape of two body integral in spatial basis: (6, 6, 6, 6)\n" + ] + } + ], + "source": [ + "import sys \n", + "sys.path.insert(0, '../')\n", + "import numpy as np\n", + "\n", + "# Now import the modules from the local moha package\n", + "from moha import HamTJUV\n", + "# Example parameters for the TJUV Hamiltonian\n", + "connectivity= np.array([[0, 1, 0, 0, 0, 1],\n", + " [1, 0, 1, 0, 0, 0],\n", + " [0, 1, 0, 1, 0, 0],\n", + " [0, 0, 1, 0, 1, 0],\n", + " [0, 0, 0, 1, 0, 1],\n", + " [1, 0, 0, 0, 1, 0]])\n", + "\n", + "\n", + "\n", + "\n", + "alpha = 0.0\n", + "beta = -1.0\n", + "u_onsite = np.array([1, 1, 1, 1, 1, 1])\n", + "gamma = None\n", + "charges = 1\n", + "sym = 8\n", + "J_eq = 1\n", + "J_ax = 1\n", + "\n", + "# Initialize the HamTJUV object\n", + "tjuv_hamiltonian = HamTJUV(connectivity=connectivity,\n", + " alpha=alpha,\n", + " beta=beta,\n", + " u_onsite=u_onsite,\n", + " gamma=gamma,\n", + " charges=charges,\n", + " sym=sym,\n", + " J_eq=J_eq,\n", + " J_ax=J_ax)\n", + "\n", + "# Generate integrals\n", + "e0 = tjuv_hamiltonian.generate_zero_body_integral()\n", + "h1 = tjuv_hamiltonian.generate_one_body_integral(dense=True, basis='spatial basis') \n", + "h2 = tjuv_hamiltonian.generate_two_body_integral(dense=True, basis='spatial basis', sym=8)\n", + "\n", + "print(\"Zero energy: \", e0)\n", + "print(\"One body integrals in spatial basis: \\n\", h1)\n", + "print(\"Shape of two body integral in spatial basis: \", h2.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Numerical energies: [-1. -1. -1. -1.]\n", + "Analytical energies: [-2.0000000e+00 -1.2246468e-16 3.6739404e-16 2.0000000e+00]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "\n", + "def test_tjuv_energy_spectrum():\n", + " # Define parameters for a simple 1D chain\n", + " N = 4\n", + " t = 1.0\n", + " alpha = -t\n", + " beta = 0.0\n", + " u_onsite = np.zeros(N)\n", + " gamma = None\n", + " charges = None\n", + " sym = None\n", + " J_eq = 0.0\n", + " J_ax = 0.0\n", + "\n", + " # Connectivity matrix for a 1D chain with periodic boundary conditions\n", + " connectivity = np.zeros((N, N))\n", + " for i in range(N):\n", + " connectivity[i, (i + 1) % N] = 1\n", + " connectivity[(i + 1) % N, i] = 1\n", + "\n", + " # Initialize the HamTJUV object\n", + " tjuv_hamiltonian = HamTJUV(connectivity=connectivity,\n", + " alpha=alpha,\n", + " beta=beta,\n", + " u_onsite=u_onsite,\n", + " gamma=gamma,\n", + " charges=charges,\n", + " sym=sym,\n", + " J_eq=J_eq,\n", + " J_ax=J_ax)\n", + "\n", + " # Generate the one-body integral (Hamiltonian matrix)\n", + " tjuv_one_body = tjuv_hamiltonian.generate_one_body_integral(basis='spatial basis', dense=True)\n", + "\n", + " # Calculate the eigenvalues (energy spectrum)\n", + " energies, _ = np.linalg.eigh(tjuv_one_body)\n", + "\n", + " # Analytical energy levels for comparison\n", + " k_vals = np.arange(N)\n", + " analytical_energies = -2 * t * np.cos(2 * np.pi * k_vals / N)\n", + "\n", + " # Sort the energies for comparison\n", + " energies_sorted = np.sort(energies)\n", + " analytical_energies_sorted = np.sort(analytical_energies)\n", + "\n", + " # Print the energy spectrum and analytical values\n", + " print(\"Numerical energies:\", energies_sorted)\n", + " print(\"Analytical energies:\", analytical_energies_sorted)\n", + "\n", + "# Outside of the function, call the test function\n", + "test_tjuv_energy_spectrum()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "One body integrals in spin basis: \n", + " [[-0.25 -1. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. ]\n", + " [-1. -0.25 -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. -1. -0.25 -1. 0. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. -1. -0.25 -1. 0. 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. -1. -0.25 -1. 0. 0. 0. 0. 0. 0. ]\n", + " [-1. 0. 0. 0. -1. -0.25 0. 0. 0. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. -0.25 -1. 0. 0. 0. -1. ]\n", + " [ 0. 0. 0. 0. 0. 0. -1. -0.25 -1. 0. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. -1. -0.25 -1. 0. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. -1. -0.25 -1. 0. ]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. -0.25 -1. ]\n", + " [ 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. -1. -0.25]]\n", + "Shape of two body integral in spinorbital basis: (12, 12, 12, 12)\n" + ] + } + ], + "source": [ + "alpha = 0.0\n", + "beta = -1.0\n", + "u_onsite = np.array([1, 1, 1, 1, 1, 1])\n", + "gamma = None\n", + "charges = 1\n", + "sym = 8\n", + "J_eq = 0.5 \n", + "J_ax = 0.5 \n", + "\n", + "# Initialize the HamTJUV object\n", + "tjuv_hamiltonian = HamTJUV(connectivity=connectivity,\n", + " alpha=alpha,\n", + " beta=beta,\n", + " u_onsite=u_onsite,\n", + " gamma=gamma,\n", + " charges=charges,\n", + " sym=sym,\n", + " J_eq=J_eq,\n", + " J_ax=J_ax)\n", + "\n", + "h1_spin = tjuv_hamiltonian.generate_one_body_integral(dense=True, basis='spinorbital basis')\n", + "h2_spin = tjuv_hamiltonian.generate_two_body_integral(dense=True, basis='spinorbital basis', sym=4)\n", + "\n", + "print(\"One body integrals in spin basis: \\n\", h1_spin)\n", + "print(\"Shape of two body integral in spinorbital basis: \", h2_spin.shape)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.11" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/rauk.ipynb b/examples/rauk.ipynb index 31d1c1d..9f6bc95 100644 --- a/examples/rauk.ipynb +++ b/examples/rauk.ipynb @@ -6,6 +6,9 @@ "source": [ "# Rauk's table and Wolfsberrg-Helmholz approximation\n", "\n", + "\n", + "## Introduction\n", + "\n", "So far the Hamiltonians treated by the documentation have been focused in carbon chains. With this new update it is possible to work with different atoms, using the Rauk's table or the distance beetween the atoms with the Wolfsberrg-Helmholz aproximation.\n", "\n", "To generate a Hamiltonian with atoms different from carbon, simply apply the same logic as before: Define a list of tuples with the atoms indexed by the position of their respective site in the compound. The new module of the library has two main purposes:\n", @@ -60,12 +63,22 @@ "import numpy as np\n", "from moha import HamHub\n", "\n", + "# Define the molecular system as a list of tuples.\n", + "# The format is ('Atom1', 'Atom2', bond_order).\n", + "# Here, 'Atom1' and 'Atom2' are the atomic symbols of the atoms involved in the bond.\n", + "# 'bond_order' is the order of the bond between the atoms.\n", + "# For example, ('C1', 'Cl2', 1) represents a single bond between a carbon atom and a chlorine atom.\n", + "system = [('C1', 'Cl2', 1), ('Cl2', 'F3', 1), ('F3', 'Si4', 1), ('Si4', 'C1', 1)]\n", "\n", - "system = [('C1', 'Cl2', 1), ('Cl2', 'F3', 1),\n", - " ('F3', 'Si4', 1), ('Si4', 'C1', 1)]\n", + "# Create an instance of the HamHub class, which represents the Hamiltonian for the system.\n", + "# Here, 'u_onsite' is an array specifying the on-site interaction energies for the atoms.\n", + "# The length of 'u_onsite' should match the number of distinct atoms in the system.\n", "hubbard = HamHub(system, u_onsite=np.array([1, 1]))\n", "\n", - "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))" + "# Generate the one-body integral matrix in the specified basis.\n", + "# 'dense=True' indicates that the integral matrix should be returned in a dense format.\n", + "# 'basis' specifies the type of basis to use, here it's 'spatial basis'.\n", + "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))\n" ] }, { @@ -87,7 +100,8 @@ "\n", "$$\n", "\\beta_{XY} = 1.75 S_{XY} \\frac{\\alpha_X + \\alpha_Y}{2}\n", - "$$" + "$$\n", + "\n" ] }, { @@ -112,12 +126,25 @@ "import numpy as np\n", "from moha import HamHub\n", "\n", - "\n", - "system = [('C1', 'Cl2', 1.5, 'pi'), ('Cl2', 'F3', 1.8, 'sigma'),\n", - " ('F3', 'Si4', 1.8, 'sigma'), ('Si4', 'C1', 1.7, 'sigma')]\n", + "# Define the molecular system as a list of tuples.\n", + "# Each tuple represents a bond between two atoms, the bond length, and the bond type.\n", + "# The bond length is a float value representing the distance between the atoms.\n", + "# The bond type can be either 'sigma' or 'pi'.\n", + "# For example, ('C1', 'Cl2', 1.5, 'pi') represents a pi bond between a carbon atom and a chlorine atom.\n", + "system = [\n", + " ('C1', 'Cl2', 1.5, 'pi'), \n", + " ('Cl2', 'F3', 1.8, 'sigma'),\n", + " ('F3', 'Si4', 1.8, 'sigma'), \n", + " ('Si4', 'C1', 1.7, 'sigma')\n", + "]\n", + "\n", + "# As before, create an instance of the HamHub class, which represents the Hamiltonian for the system.\n", "hubbard = HamHub(system, u_onsite=np.array([1, 1]))\n", "\n", - "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))" + "# Generate the one-body integral matrix in the specified basis.\n", + "# 'dense=True' indicates that the integral matrix should be returned in a dense format.\n", + "# 'basis' specifies the type of basis to use, here it's 'spatial basis'.\n", + "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))\n" ] }, { @@ -126,7 +153,7 @@ "source": [ "# Example: Usage with molecules \n", "\n", - "If you want to make a chain using molecules you can use the following structure: pass the atom $X_n$ in the position m, the tuple would have the following form `Xm(n)`. Below it's possible to see a example of this. " + "If you want to make a chain using molecules you can use the following structure: pass the atom $X_n$ in the position m, the tuple would have the following form `Xm(n)`. Below it's possible to see a example of this of the molecule $N_2OP_2$. " ] }, { @@ -147,9 +174,20 @@ "source": [ "from moha import HamHuck\n", "\n", + "# Define the molecular system with tuples in the format\n", + "# ('Atomposition(index)', 'Neighboorposition(index)', 'bond_order').\n", + "# Each tuple represents a bond between atoms at specified positions.\n", + "system =[('N1(2)', 'O2(1)', 1),\n", + " ('O2(1)', 'P3(2)', 1),\n", + " ('P3(2)', 'N1(2)', 1)]\n", "\n", - "N2_O1_P2 = HamHuck([('N1(2)', 'O2(1)', 1), ('O2(1)', 'P3(2)', 1), ('P3(2)', 'N1(2)', 1)])\n", + "# Create an instance of HamHuck with the defined system.\n", + "N2_O1_P2 = HamHuck(system)\n", + "\n", + "# Generate the one-body integral matrix in the 'spatial basis'.\n", "h1 = N2_O1_P2.generate_one_body_integral(dense=True, basis='spatial basis')\n", + "\n", + "# Print the one-body integral matrix.\n", "print(h1)" ] }, @@ -157,12 +195,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Example: Defining the dictionaries \n", + "## Example: Defining Your Own Dictionaries\n", + "\n", + "Instead of using Rauk's Table or the Wolfsberg-Helmholz approximation, you can define your own dictionaries for $\\alpha$ and $\\beta$ values. These dictionaries specify the interaction parameters for each atom and bond in your system.\n", "\n", - "It is also possible to define your own dictionary, basically pass the dictionaries following the structure:\n", + "To do this, pass the dictionaries to the Hamiltonian classes with the following structure:\n", "\n", - "1) `atom_dictionary = {'atom1': alpha1, 'atom2': alpha2}`\n", - "2) `bond_dictionary = {'atom1atom2': beta}`" + "1. `atom_dictionary = {'atom1': alpha1, 'atom2': alpha2}`\n", + "2. `bond_dictionary = {'atom1,atom2': beta}`\n" ] }, { @@ -174,10 +214,10 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[1. 0. 0. 0.]\n", + "[[1. 0. 0. 3.]\n", " [0. 2. 1. 0.]\n", - " [0. 0. 3. 2.]\n", - " [3. 0. 0. 4.]]\n" + " [0. 1. 3. 2.]\n", + " [3. 0. 2. 4.]]\n" ] } ], @@ -187,21 +227,53 @@ "import numpy as np\n", "from moha import HamHub\n", "\n", - "# First way to define the Hamiltonian\n", - "# two site Hubbard model\n", - "\n", - "# system = [('C1', 'C2', 1)] is a list of tuples, where each tuple represents a bond\n", - "# between two atoms and the third element is the type of bond (singe or double).\n", - "# For now, we only support single bonds between carbon atoms. \n", - "# For this type of bonds the default values of alpha and beta are -0.414 and -0.0533, respectively.\n", - "# In the future we are planning to support different types of bonds for different atoms.\n", - "system = [('C1', 'Cl2', 1), ('Cl2', 'F3', 1),\n", - " ('F3', 'Si4', 1), ('Si4', 'C1', 1)]\n", + "# Define the Hamiltonian for a molecular system using the HamHub class.\n", + "\n", + "# The system is a list of tuples, where each tuple represents a bond between two atoms.\n", + "# The format of each tuple is ('Atom1', 'Atom2', bond_order).\n", + "# For this example, we define single bonds between different atoms.\n", + "system = [\n", + " ('C1', 'Cl2', 1), # Bond between C1 and Cl2\n", + " ('Cl2', 'F3', 1), # Bond between Cl2 and F3\n", + " ('F3', 'Si4', 1), # Bond between F3 and Si4\n", + " ('Si4', 'C1', 1) # Bond between Si4 and C1\n", + "]\n", + "\n", + "# Create an instance of the HamHub class, which represents the Hamiltonian for the system.\n", + "# 'u_onsite' is an array specifying the on-site interaction energies for the atoms.\n", + "# 'atom_dictionary' maps atom symbols to unique identifiers.\n", + "# 'bond_dictionary' maps bond pairs to unique identifiers.\n", "hubbard = HamHub(system, u_onsite=np.array([1, 1]), \n", " atom_dictionary={'C': 1, 'Cl': 2, 'F': 3, 'Si': 4},\n", " bond_dictionary={'C,Cl': 0, 'Cl,F': 1, 'F,Si': 2, 'Si,C': 3})\n", "\n", - "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))" + "# Generate the one-body integral matrix in the specified basis.\n", + "# 'dense=True' indicates that the integral matrix should be returned in a dense format.\n", + "# 'basis' specifies the type of basis to use, here it's 'spatial basis'.\n", + "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example: Second Body Terms - Pariser-Parr Module\n", + "\n", + "The software is highly flexible in generating the two-body terms.\n", + "\n", + "If the user does not provide `u_onsite` or `affinity_dct`, the potentials are computed using the following formula: $ U_x = \\alpha_x - A_x $.\n", + "\n", + "If `orbital_overlap`is not provided, the values of `u_onsite`, either provided or computed, are used to calculate the potential between two sites:\n", + "\n", + "$$\n", + "\\overline{U}_{XY} = \\frac{1}{2}(U_X + U_Y)\n", + "$$\n", + "\n", + "If the value of `gamma` is not provided, it is computed using the formula:\n", + "\n", + "$$\n", + "\\gamma_{XY} = \\frac{\\overline{U}_{XY}}{\\overline{U}_{XY} R_{XY} + e^{-\\frac{1}{2} \\overline{U}_{XY}^2 R_{XY}^2}}\n", + "$$" ] }, { @@ -213,10 +285,32 @@ "name": "stdout", "output_type": "stream", "text": [ - "[[-0.41380839 -0.77906404 0. -1.87261684]\n", - " [-0.77906404 -0.47655051 -1.95444655 0. ]\n", - " [ 0. -1.95444655 -0.64027609 -1.64472969]\n", - " [-1.87261684 0. -1.64472969 -0.29956945]]\n" + "Zero energy: 0\n", + "One body integrals in spatial basis: \n", + " [[-0.41380839 -0.77906404 0. 0. ]\n", + " [-0.77906404 -0.47655051 -2.93166983 0. ]\n", + " [ 0. -2.93166983 -0.64027609 -3.28945939]\n", + " [ 0. 0. -3.28945939 -0.29956945]]\n", + "Shape of two body integral in spatial basis: (4, 4, 4, 4)\n", + "------------------------------------------------------------\n", + "One body integrals in spinorbital basis: \n", + " [[-0.41380839 -0.77906404 0. 0. 0. 0.\n", + " 0. 0. ]\n", + " [-0.77906404 -0.47655051 -2.93166983 0. 0. 0.\n", + " 0. 0. ]\n", + " [ 0. -2.93166983 -0.64027609 -3.28945939 0. 0.\n", + " 0. 0. ]\n", + " [ 0. 0. -3.28945939 -0.29956945 0. 0.\n", + " 0. 0. ]\n", + " [ 0. 0. 0. 0. -0.41380839 -0.77906404\n", + " 0. 0. ]\n", + " [ 0. 0. 0. 0. -0.77906404 -0.47655051\n", + " -2.93166983 0. ]\n", + " [ 0. 0. 0. 0. 0. -2.93166983\n", + " -0.64027609 -3.28945939]\n", + " [ 0. 0. 0. 0. 0. 0.\n", + " -3.28945939 -0.29956945]]\n", + "Shape of two body integral in spinorbital basis: (8, 8, 8, 8)\n" ] } ], @@ -226,20 +320,144 @@ "import numpy as np\n", "from moha import HamHub\n", "\n", - "# First way to define the Hamiltonian\n", - "# two site Hubbard model\n", - "\n", - "# system = [('C1', 'C2', 1)] is a list of tuples, where each tuple represents a bond\n", - "# between two atoms and the third element is the type of bond (singe or double).\n", - "# For now, we only support single bonds between carbon atoms. \n", - "# For this type of bonds the default values of alpha and beta are -0.414 and -0.0533, respectively.\n", - "# In the future we are planning to support different types of bonds for different atoms.\n", - "system = [('C1', 'Cl2', 1.1), ('Cl2', 'F3', 1.0),\n", - " ('F3', 'Si4', 1.0), ('Si4', 'C1', 1.0)]\n", - "hubbard = HamHub(system, u_onsite=np.array([1, 1]),\n", - " orbital_overlap={'C,Cl': 1, 'Cl,F': 2, 'F,Si': 2, 'Si,C': 3})\n", + "# Example: Generating XXZ Heisenberg model \n", + "# Returning electron integrals in a spatial orbital basis\n", + "# Assuming 4-fold symmetry\n", + "# Returning output as dense matrix\n", + "\n", + "# Define the molecular system with tuples in the format ('Atom1', 'Atom2', bond_length).\n", + "# The system represents bonds between atoms with specified bond lengths.\n", + "system = [\n", + " ('C1', 'Cl2', 1.1), # Bond between C1 and Cl2 with bond length 1.1\n", + " ('Cl2', 'F3', 1.0), # Bond between Cl2 and F3 with bond length 1.0\n", + " ('F3', 'Si4', 1.0), # Bond between F3 and Si4 with bond length 1.0\n", + " ('Si4', 'C1', 1.0) # Bond between Si4 and C1 with bond length 1.0\n", + "]\n", + "\n", + "# Create an instance of the HamHub class with the defined system.\n", + "# u_onsite specifies the on-site interaction energies for the atoms.\n", + "# orbital_overlap specifies the overlap integrals between orbitals.\n", + "# So the gamma will be computed based on orbital_overlap.\n", + "ham = HamHub(system, u_onsite=np.array([1, 1, 1, 1]),\n", + " orbital_overlap=np.array([[0, 1, 0, 0],[0, 0, 3, 0], [0, 0, 0, 4], [0, 0, 0, 1]]))\n", + "\n", + "# Generate the zero-body integral.\n", + "e0 = ham.generate_zero_body_integral()\n", + "\n", + "# Generate the one-body integral matrix in the spatial orbital basis.\n", + "h1 = ham.generate_one_body_integral(dense=True, basis='spatial basis')\n", + "\n", + "# Generate the two-body integral tensor in the spatial orbital basis, assuming 4-fold symmetry.\n", + "h2 = ham.generate_two_body_integral(dense=True, basis='spatial basis', sym=4)\n", + "\n", + "print(\"Zero energy: \", e0)\n", + "print(\"One body integrals in spatial basis: \\n\", h1)\n", + "print(\"Shape of two body integral in spatial basis: \", h2.shape)\n", + "print(\"-\" * 60)\n", + "\n", + "# Example: Generating XXZ Heisenberg model in spin orbital basis\n", + "# Assuming 4-fold symmetry\n", + "# Returning output as dense matrix\n", + "\n", + "# Generate the one-body integral matrix in the spin orbital basis.\n", + "h1 = ham.generate_one_body_integral(dense=True, basis='spinorbital basis')\n", + "\n", + "# Generate the two-body integral tensor in the spin orbital basis, assuming 4-fold symmetry.\n", + "h2 = ham.generate_two_body_integral(dense=True, basis='spinorbital basis', sym=4)\n", + "\n", + "print(\"One body integrals in spinorbital basis: \\n\", h1)\n", + "print(\"Shape of two body integral in spinorbital basis: \", h2.shape)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Zero energy: 0\n", + "One body integrals in spatial basis: \n", + " [[ 0. -1. 0. 0. 0. -1.]\n", + " [-1. 0. -1. 0. 0. 0.]\n", + " [ 0. -1. 0. -1. 0. 0.]\n", + " [ 0. 0. -1. 0. -1. 0.]\n", + " [ 0. 0. 0. -1. 0. -1.]\n", + " [-1. 0. 0. 0. -1. 0.]]\n", + "Shape of two body integral in spatial basis: (6, 6, 6, 6)\n", + "------------------------------------------------------------\n", + "One body integrals in spinorbital basis: \n", + " [[ 0. -1. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0.]\n", + " [-1. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. -1. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. -1. 0. -1. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. -1. 0. -1. 0. 0. 0. 0. 0. 0.]\n", + " [-1. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. -1.]\n", + " [ 0. 0. 0. 0. 0. 0. -1. 0. -1. 0. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. -1. 0. -1. 0. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. -1. 0.]\n", + " [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. -1.]\n", + " [ 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. -1. 0.]]\n", + "Shape of two body integral in spinorbital basis: (12, 12, 12, 12)\n" + ] + } + ], + "source": [ + "import sys\n", + "sys.path.insert(0, '../')\n", + "import numpy as np\n", + "from moha import HamHub\n", "\n", - "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))" + "# Example: Generating 6-site Hubbard model \n", + "# Returning electron integrals in a spatial orbital basis\n", + "# Assuming 8-fold symmetry\n", + "# Returning output as dense matrix\n", + "\n", + "# Define the connectivity matrix for the 6-site Hubbard model.\n", + "# Each entry in the matrix indicates the connection between sites.\n", + "connectivity = np.array([\n", + " [0, 1, 0, 0, 0, 1], \n", + " [1, 0, 1, 0, 0, 0], \n", + " [0, 1, 0, 1, 0, 0], \n", + " [0, 0, 1, 0, 1, 0], \n", + " [0, 0, 0, 1, 0, 1], \n", + " [1, 0, 0, 0, 1, 0] \n", + "])\n", + "\n", + "# Create an instance of the HamHub class with the defined connectivity.\n", + "# alpha and beta are interaction parameters.\n", + "# u_onsite specifies the on-site interaction energies for the sites.\n", + "hubbard = HamHub(connectivity, alpha=0, beta=-1, u_onsite=np.array([1, 1, 1, 1, 1, 1]))\n", + "\n", + "# Generate the zero-body integral.\n", + "e0 = hubbard.generate_zero_body_integral()\n", + "\n", + "# Generate the one-body integral matrix in the spatial orbital basis.\n", + "h1 = hubbard.generate_one_body_integral(dense=True, basis='spatial basis')\n", + "\n", + "# Once the gamma was not provided it will be computed based on the one-body integrals.\n", + "h2 = hubbard.generate_two_body_integral(dense=True, basis='spatial basis', sym=8)\n", + "\n", + "print(\"Zero energy: \", e0)\n", + "print(\"One body integrals in spatial basis: \\n\", h1)\n", + "print(\"Shape of two body integral in spatial basis: \", h2.shape)\n", + "print(\"-\" * 60)\n", + "\n", + "# Example: Generating Hubbard model in spin orbital basis\n", + "# Assuming 8-fold symmetry\n", + "# Returning output as dense matrix\n", + "\n", + "# Generate the one-body integral matrix in the spin orbital basis.\n", + "h1 = hubbard.generate_one_body_integral(dense=True, basis='spinorbital basis')\n", + "\n", + "# Generate the two-body integral tensor in the spin orbital basis, assuming 8-fold symmetry.\n", + "h2 = hubbard.generate_two_body_integral(dense=True, basis='spinorbital basis', sym=8)\n", + "\n", + "print(\"One body integrals in spinorbital basis: \\n\", h1)\n", + "print(\"Shape of two body integral in spinorbital basis: \", h2.shape)\n" ] }, { @@ -248,7 +466,9 @@ "source": [ "## References:\n", "\n", - "[1]A. Orbital Interaction Theory of Organic Chemistry. [s.l.] John Wiley & Sons, 2004." + "[1]A. Orbital Interaction Theory of Organic Chemistry. [s.l.] John Wiley & Sons, 2004.\n", + "\n", + "[2] https://en.wikipedia.org/wiki/Electron_affinity_(data_page)" ] } ], @@ -268,7 +488,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.11" + "version": "3.10.undefined" } }, "nbformat": 4, diff --git a/examples/toml/fcidump/caffeine.fcidump b/examples/toml/fcidump/caffeine.fcidump new file mode 100644 index 0000000..95325ee --- /dev/null +++ b/examples/toml/fcidump/caffeine.fcidump @@ -0,0 +1,5 @@ + &FCI NORB=0,NELEC=0,MS2=0, + ORBSYM= , + ISYM=1 + &END + 0.0000000000000000e+00 0 0 0 0 diff --git a/examples/toml/fcidump/heisenberg.fcidump b/examples/toml/fcidump/heisenberg.fcidump new file mode 100644 index 0000000..4b91235 --- /dev/null +++ b/examples/toml/fcidump/heisenberg.fcidump @@ -0,0 +1,29 @@ + &FCI NORB=6,NELEC=6,MS2=0, + ORBSYM= 1,1,1,1,1,1, + ISYM=1 + &END + 5.0000000000000000e-01 1 2 1 2 + 5.0000000000000000e-01 1 6 1 6 + 2.5000000000000000e-01 2 2 1 1 + 5.0000000000000000e-01 2 1 2 1 + 5.0000000000000000e-01 2 3 2 3 + 2.5000000000000000e-01 3 3 2 2 + 5.0000000000000000e-01 3 2 3 2 + 5.0000000000000000e-01 3 4 3 4 + 2.5000000000000000e-01 4 4 3 3 + 5.0000000000000000e-01 4 3 4 3 + 5.0000000000000000e-01 4 5 4 5 + 2.5000000000000000e-01 5 5 4 4 + 5.0000000000000000e-01 5 4 5 4 + 5.0000000000000000e-01 5 6 5 6 + 2.5000000000000000e-01 6 6 1 1 + 2.5000000000000000e-01 6 6 5 5 + 5.0000000000000000e-01 6 1 6 1 + 5.0000000000000000e-01 6 5 6 5 +-4.0000000000000002e-01 1 1 0 0 +-4.0000000000000002e-01 2 2 0 0 +-4.0000000000000002e-01 3 3 0 0 +-4.0000000000000002e-01 4 4 0 0 +-4.0000000000000002e-01 5 5 0 0 +-4.0000000000000002e-01 6 6 0 0 + 9.0000000000000002e-01 0 0 0 0 diff --git a/examples/toml/fcidump/hubbard.fcidump b/examples/toml/fcidump/hubbard.fcidump new file mode 100644 index 0000000..6e84afd --- /dev/null +++ b/examples/toml/fcidump/hubbard.fcidump @@ -0,0 +1,16 @@ + &FCI NORB=6,NELEC=6,MS2=0, + ORBSYM= 1,1,1,1,1,1, + ISYM=1 + &END + 1.0000000000000000e+00 1 1 1 1 + 1.0000000000000000e+00 2 2 2 2 + 1.0000000000000000e+00 3 3 3 3 + 1.0000000000000000e+00 4 4 4 4 + 1.0000000000000000e+00 5 5 5 5 + 1.0000000000000000e+00 6 6 6 6 +-1.0000000000000000e+00 2 1 0 0 +-1.0000000000000000e+00 3 2 0 0 +-1.0000000000000000e+00 4 3 0 0 +-1.0000000000000000e+00 5 4 0 0 +-1.0000000000000000e+00 6 5 0 0 + 0.0000000000000000e+00 0 0 0 0 diff --git a/examples/toml/fcidump/huckel.fcidump b/examples/toml/fcidump/huckel.fcidump new file mode 100644 index 0000000..62401f3 --- /dev/null +++ b/examples/toml/fcidump/huckel.fcidump @@ -0,0 +1,17 @@ + &FCI NORB=6,NELEC=6,MS2=0, + ORBSYM= 1,1,1,1,1,1, + ISYM=1 + &END +-4.1399999999999998e-01 1 1 0 0 +-5.3300000000000000e-02 2 1 0 0 +-4.1399999999999998e-01 2 2 0 0 +-5.3300000000000000e-02 3 2 0 0 +-4.1399999999999998e-01 3 3 0 0 +-5.3300000000000000e-02 4 3 0 0 +-4.1399999999999998e-01 4 4 0 0 +-5.3300000000000000e-02 5 4 0 0 +-4.1399999999999998e-01 5 5 0 0 +-5.3300000000000000e-02 6 1 0 0 +-5.3300000000000000e-02 6 5 0 0 +-4.1399999999999998e-01 6 6 0 0 + 0.0000000000000000e+00 0 0 0 0 diff --git a/examples/toml/fcidump/ising.fcidump b/examples/toml/fcidump/ising.fcidump new file mode 100644 index 0000000..6ba31f9 --- /dev/null +++ b/examples/toml/fcidump/ising.fcidump @@ -0,0 +1,17 @@ + &FCI NORB=6,NELEC=6,MS2=0, + ORBSYM= 1,1,1,1,1,1, + ISYM=1 + &END + 2.5000000000000000e-01 2 2 1 1 + 2.5000000000000000e-01 3 3 2 2 + 2.5000000000000000e-01 4 4 3 3 + 2.5000000000000000e-01 5 5 4 4 + 2.5000000000000000e-01 6 6 1 1 + 2.5000000000000000e-01 6 6 5 5 +-4.0000000000000002e-01 1 1 0 0 +-4.0000000000000002e-01 2 2 0 0 +-4.0000000000000002e-01 3 3 0 0 +-4.0000000000000002e-01 4 4 0 0 +-4.0000000000000002e-01 5 5 0 0 +-4.0000000000000002e-01 6 6 0 0 + 9.0000000000000002e-01 0 0 0 0 diff --git a/examples/toml/fcidump/ppp.fcidump b/examples/toml/fcidump/ppp.fcidump new file mode 100644 index 0000000..68f1ca0 --- /dev/null +++ b/examples/toml/fcidump/ppp.fcidump @@ -0,0 +1,23 @@ + &FCI NORB=6,NELEC=6,MS2=0, + ORBSYM= 1,1,1,1,1,1, + ISYM=1 + &END + 1.0000000000000000e+00 1 1 1 1 + 1.0000000000000000e+00 2 2 2 2 + 1.0000000000000000e+00 3 3 3 3 + 1.0000000000000000e+00 4 4 4 4 + 1.0000000000000000e+00 5 5 5 5 + 1.0000000000000000e+00 6 6 6 6 + 5.0000000000000000e-01 1 1 0 0 +-1.0000000000000000e+00 2 1 0 0 + 5.0000000000000000e-01 2 2 0 0 +-1.0000000000000000e+00 3 2 0 0 + 5.0000000000000000e-01 3 3 0 0 +-1.0000000000000000e+00 4 3 0 0 + 5.0000000000000000e-01 4 4 0 0 +-1.0000000000000000e+00 5 4 0 0 + 5.0000000000000000e-01 5 5 0 0 +-1.0000000000000000e+00 6 1 0 0 +-1.0000000000000000e+00 6 5 0 0 + 5.0000000000000000e-01 6 6 0 0 + 0.0000000000000000e+00 0 0 0 0 diff --git a/examples/toml/fcidump/rg.fcidump b/examples/toml/fcidump/rg.fcidump new file mode 100644 index 0000000..169b478 --- /dev/null +++ b/examples/toml/fcidump/rg.fcidump @@ -0,0 +1,23 @@ + &FCI NORB=6,NELEC=6,MS2=0, + ORBSYM= 1,1,1,1,1,1, + ISYM=1 + &END + 5.0000000000000000e-01 1 2 1 2 + 5.0000000000000000e-01 1 6 1 6 + 5.0000000000000000e-01 2 1 2 1 + 5.0000000000000000e-01 2 3 2 3 + 5.0000000000000000e-01 3 2 3 2 + 5.0000000000000000e-01 3 4 3 4 + 5.0000000000000000e-01 4 3 4 3 + 5.0000000000000000e-01 4 5 4 5 + 5.0000000000000000e-01 5 4 5 4 + 5.0000000000000000e-01 5 6 5 6 + 5.0000000000000000e-01 6 1 6 1 + 5.0000000000000000e-01 6 5 6 5 + 1.0000000000000001e-01 1 1 0 0 + 1.0000000000000001e-01 2 2 0 0 + 1.0000000000000001e-01 3 3 0 0 + 1.0000000000000001e-01 4 4 0 0 + 1.0000000000000001e-01 5 5 0 0 + 1.0000000000000001e-01 6 6 0 0 +-5.9999999999999998e-01 0 0 0 0 diff --git a/moha/__init__.py b/moha/__init__.py index 6be4637..6584fc6 100644 --- a/moha/__init__.py +++ b/moha/__init__.py @@ -1,9 +1,15 @@ r"""Model Hamiltonian module.""" -from .hamiltonians import HamPPP, HamHub, HamHuck, \ - HamHeisenberg, HamIsing, HamRG - +from .hamiltonians import ( + HamPPP, + HamHub, + HamHuck, + HamHeisenberg, + HamIsing, + HamRG, + HamTJUV, +) __all__ = [ "HamPPP", @@ -11,5 +17,6 @@ "HamHuck", "HamHeisenberg", "HamIsing", - "HamRG" + "HamRG", + "HamTJUV", ] diff --git a/moha/api.py b/moha/api.py index 80cf251..e74516c 100644 --- a/moha/api.py +++ b/moha/api.py @@ -10,8 +10,6 @@ from .utils import convert_indices -from moha.rauk.utils import get_atom_type, get_atoms_list - from typing import Union @@ -23,52 +21,6 @@ class HamiltonianAPI(ABC): r"""Hamiltonian abstract base class.""" - def generate_connectivity_matrix(self): - r""" - - Generate connectivity matrix. - - Returns - ------- - tuple - (dictionary, np.ndarray) - """ - # check if self.connectivity is a matrix - # if so, put assign it to self.connectivity_matrix - # and set the atom_types to None - if isinstance(self.connectivity, np.ndarray): - self.connectivity_matrix = csr_matrix(self.connectivity) - self.atom_types = None - self.n_sites = self.connectivity_matrix.shape[0] - - return None, self.connectivity_matrix - - atoms_sites_lst = get_atoms_list(self.connectivity) - max_site = max([site for _, site in atoms_sites_lst]) - self.n_sites = max_site - - if self.atom_types is None: - # Initialize atom_types with None, and adjust size for 0-based - # indexing - atom_types = [None] * max_site - for atom, site in atoms_sites_lst: - # Adjust site index for 0-based array index - atom_types[site - 1] = atom - self.atom_types = atom_types - connectivity_mtrx = np.zeros((max_site, max_site)) - atoms_dist = [] - for tpl in self.connectivity: - atom1, atom2, bond = tpl[0], tpl[1], tpl[2] - atom1_name, site1 = get_atom_type(atom1) - atom2_name, site2 = get_atom_type(atom2) - atoms_dist.append((atom1_name, atom2_name, bond)) - connectivity_mtrx[site1 - 1, site2 - 1] = bond - # numbering of sites starts from 1 - self.atoms_dist = atoms_dist - connectivity_mtrx = np.maximum(connectivity_mtrx, connectivity_mtrx.T) - self.connectivity_matrix = csr_matrix(connectivity_mtrx) - return atoms_sites_lst, self.connectivity_matrix - @abstractmethod def generate_zero_body_integral(self): r"""Generate zero body integral.""" diff --git a/moha/hamiltonians.py b/moha/hamiltonians.py index 21c1227..15524e8 100644 --- a/moha/hamiltonians.py +++ b/moha/hamiltonians.py @@ -10,9 +10,9 @@ from .utils import convert_indices, expand_sym from typing import Union - +from moha.rauk.utils import parse_connectivity from moha.rauk.rauk import assign_rauk_parameters -from moha.rauk.PariserParr import compute_overlap +from moha.rauk.PariserParr import compute_overlap, compute_gamma, compute_u import warnings warnings.simplefilter('ignore', @@ -31,7 +31,8 @@ class HamPPP(HamiltonianAPI): def __init__( self, - connectivity: Union[list, np.ndarray], + connectivity=None, + adjacency=None, alpha=-0.414, beta=-0.0533, u_onsite=None, @@ -40,19 +41,22 @@ def __init__( sym=1, atom_dictionary=None, bond_dictionary=None, - orbital_overlap=None + orbital_overlap=None, + affinity_dct=None, + Rxy_matrix=None, + atom_types=None, ): r""" Initialize Pariser-Parr-Pople Hamiltonian. Parameters ---------- - connectivity: list, np.ndarray + connectivity: list list of tuples that specifies sites and bonds between them - or symmetric np.ndarray of shape (n_sites, n_sites) that specifies - the connectivity between sites. For example, for a linear chain of 4 sites, the connectivity can be specified as [(C1, C2, 1), (C2, C3, 1), (C3, C4, 1)] + adjacency: np.ndarray + symmetric numpy array that specifies the adjacency between sites alpha: float specifies the site energy if all sites are equivalent. Default value is the 2p-pi orbital of Carbon @@ -86,23 +90,31 @@ def __init__( """ self._sym = sym - self.n_sites = None self.connectivity = connectivity + self.adjacency = adjacency self.alpha = alpha self.beta = beta self.u_onsite = u_onsite - self.gamma = gamma + if gamma is None: + raise TypeError( + "Gamma matrix is not provided, use the Hubbard model") + else: + self.gamma = gamma self.charges = charges - self.atom_types = None + self.atom_types = atom_types self.atoms_dist = None - self.atoms_num, self.connectivity_matrix = \ - self.generate_connectivity_matrix() self.zero_energy = None self.one_body = None self.two_body = None self.bond_dictionary = bond_dictionary self.atom_dictionary = atom_dictionary self.orbital_overlap = orbital_overlap + self.affinity_dct = affinity_dct + self.Rxy_matrix = Rxy_matrix + if self.connectivity is not None: + _, _, self.n_sites, _ = parse_connectivity(self.connectivity) + if self.adjacency is not None: + self.n_sites = self.adjacency.shape[0] def generate_zero_body_integral(self): r"""Generate zero body integral. @@ -112,6 +124,7 @@ def generate_zero_body_integral(self): float """ if self.charges is None or self.gamma is None: + print(0) self.zero_energy = 0 return 0 else: @@ -133,31 +146,45 @@ def generate_one_body_integral(self, basis: str, dense: bool): ------- scipy.sparse.csr_matrix or np.ndarray """ - # check if connectivity matrix is adjacency - if isinstance(self.connectivity, np.ndarray): - one_body_term = ( - diags([self.alpha for _ in range(self.n_sites)], format="csr") - + self.beta * self.connectivity_matrix - ) - # check if alpha and beta are different from the default or - # all atom types are the same - elif ( + if self.adjacency is None and np.all( + [isinstance(elem[2], int) for elem in self.connectivity]): + self.atoms_sites_lst, self.adjacency, self.n_sites, \ + self.atom_types = parse_connectivity(self.connectivity) + + # check if adjacency matrix is provided + if self.adjacency is not None and ( self.alpha != -0.414 and self.beta != -0.0533 ) or len(np.unique(self.atom_types)) == 1: + one_body_term = ( - diags([self.alpha for _ in range(self.n_sites)]) - + self.beta * self.connectivity_matrix + diags([self.alpha for _ in range(self.n_sites)], format="csr") + + self.beta * csr_matrix(self.adjacency) ) - # check if elements in connectivity matrix are integer elif np.all([isinstance(elem[2], int) for elem in self.connectivity]): one_body_term = assign_rauk_parameters( self.connectivity, self.atom_dictionary, self.bond_dictionary ) + # check if elements in connectivity matrix are float elif np.all([isinstance(elem[2], float) for elem in self.connectivity]): + if self.orbital_overlap is not None: + # Assuming self.orbital_overlap should be a 2D array + if hasattr( + self, + 'orbital_overlap') and self.orbital_overlap.ndim == 2: + if self.n_sites != self.orbital_overlap.shape[ + 0] or self.n_sites != self.orbital_overlap.shape[ + 1]: + raise TypeError("Overlap matrix has wrong dimensions") + else: + raise ValueError( + "orbital_overlap is not properly initialized ") + +# Proceed with the rest of the function + one_body_term = compute_overlap( self.connectivity, self.atom_dictionary, @@ -213,38 +240,46 @@ def generate_two_body_integral(self, basis: str, dense: bool, sym=1): ------- scipy.sparse.csr_matrix or np.ndarray """ + if self.adjacency is None: + self.atoms_sites_lst, self.adjacency, self.n_sites, \ + self.atom_types = parse_connectivity(self.connectivity) + self.n_sites = self.adjacency.shape[0] n_sp = self.n_sites Nv = 2 * n_sp v = lil_matrix((Nv * Nv, Nv * Nv)) + if self.u_onsite is None: + self.u_onsite, self.Rxy_matrix = compute_u( + self.connectivity, self.atom_dictionary, + self.affinity_dct, self.atom_types) if self.u_onsite is not None: for p in range(n_sp): i, j = convert_indices(Nv, p, p + n_sp, p, p + n_sp) + v[i, j] = self.u_onsite[p] - if self.gamma is not None: - if basis == "spinorbital basis" and \ - self.gamma.shape != (n_sp, n_sp): - raise TypeError("Gamma matrix has wrong shape") + if basis == "spinorbital basis" and \ + self.gamma.shape != (n_sp, n_sp): + raise TypeError("Gamma matrix has wrong shape") - for p in range(n_sp): - for q in range(n_sp): - if p != q: - i, j = convert_indices(Nv, p, q, p, q) - v[i, j] = 0.5 * self.gamma[p, q] + for p in range(n_sp): + for q in range(n_sp): + if p != q: + i, j = convert_indices(Nv, p, q, p, q) + v[i, j] = 0.5 * self.gamma[p, q] - i, j = convert_indices(Nv, p, q + n_sp, p, q + n_sp) - v[i, j] = 0.5 * self.gamma[p, q] + i, j = convert_indices(Nv, p, q + n_sp, p, q + n_sp) + v[i, j] = 0.5 * self.gamma[p, q] - i, j = convert_indices(Nv, p + n_sp, q, p + n_sp, q) - v[i, j] = 0.5 * self.gamma[p, q] + i, j = convert_indices(Nv, p + n_sp, q, p + n_sp, q) + v[i, j] = 0.5 * self.gamma[p, q] - i, j = convert_indices(Nv, - p + n_sp, - q + n_sp, - p + n_sp, - q + n_sp) - v[i, j] = 0.5 * self.gamma[p, q] + i, j = convert_indices(Nv, + p + n_sp, + q + n_sp, + p + n_sp, + q + n_sp) + v[i, j] = 0.5 * self.gamma[p, q] v = v.tocsr() self.two_body = expand_sym(sym, v, 2) @@ -272,7 +307,8 @@ class HamHub(HamPPP): def __init__( self, - connectivity: Union[list, np.ndarray], + connectivity=None, + adjacency=None, alpha=-0.414, beta=-0.0533, u_onsite=None, @@ -282,19 +318,19 @@ def __init__( bond_dictionary=None, orbital_overlap=None, Bz=None, - gamma=None, + gamma=0, ): r""" Hubbard Hamiltonian. Parameters ---------- - connectivity: list, np.ndarray + connectivity: list list of tuples that specifies sites and bonds between them - or symmetric np.ndarray of shape (n_sites, n_sites) that specifies - the connectivity between sites. For example, for a linear chain of 4 sites, the connectivity can be specified as [(C1, C2, 1), (C2, C3, 1), (C3, C4, 1)] + adjacency: np.ndarray + symmetric numpy array that specifies the adjacency between sites alpha: float specifies the site energy if all sites are equivalent. Default value is the 2p-pi orbital of Carbon @@ -319,17 +355,19 @@ def __init__( """ super().__init__( connectivity=connectivity, + adjacency=adjacency, alpha=alpha, beta=beta, u_onsite=u_onsite, - gamma=None, atom_dictionary=atom_dictionary, bond_dictionary=bond_dictionary, orbital_overlap=orbital_overlap, charges=np.array(0), - sym=sym + sym=sym, + gamma=gamma, ) self.charges = np.zeros(self.n_sites) + self.gamma = np.zeros((self.n_sites, self.n_sites)) class HamHuck(HamHub): @@ -341,24 +379,26 @@ class HamHuck(HamHub): def __init__( self, - connectivity: Union[list, np.ndarray], + connectivity=None, + adjacency=None, alpha=-0.414, beta=-0.0533, sym=1, atom_dictionary=None, bond_dictionary=None, + atom_types=None, ): r""" Huckel hamiltonian. Parameters ---------- - connectivity: list, np.ndarray + connectivity: list list of tuples that specifies sites and bonds between them - or symmetric np.ndarray of shape (n_sites, n_sites) that specifies - the connectivity between sites. For example, for a linear chain of 4 sites, the connectivity can be specified as [(C1, C2, 1), (C2, C3, 1), (C3, C4, 1)] + adjacency: np.ndarray + symmetric numpy array that specifies the adjacency between sites alpha: float specifies the site energy if all sites are equivalent. Default value is the 2p-pi orbital of Carbon @@ -380,14 +420,17 @@ def __init__( """ super().__init__( connectivity=connectivity, + adjacency=adjacency, alpha=alpha, beta=beta, u_onsite=None, - gamma=None, + gamma=np.zeros(adjacency.shape), sym=sym, atom_dictionary=atom_dictionary, bond_dictionary=bond_dictionary, + atom_types=atom_types, ) + self.atom_types = atom_types self.charges = np.zeros(self.n_sites) @@ -398,7 +441,7 @@ def __init__(self, mu: np.ndarray, J_eq: np.ndarray, J_ax: np.ndarray, - connectivity: np.ndarray = None + adjacency: np.ndarray = None ): r"""Initialize XXZ Heisenberg Hamiltonian. @@ -410,8 +453,8 @@ def __init__(self, J equatorial term J_ax: np.ndarray J axial term - connectivity: np.ndarray - symmetric numpy array that specifies the connectivity between sites + adjacency: np.ndarray + symmetric numpy array that specifies the adjacency between sites Notes ----- @@ -423,17 +466,17 @@ def __init__(self, J_{p q}^{\mathrm{eq}} S_p^{+} S_q^{-} """ - if connectivity is not None: - self.connectivity = connectivity - self.n_sites = connectivity.shape[0] + if adjacency is not None: + self.adjacency = adjacency + self.n_sites = adjacency.shape[0] # if J_eq and J_ax are floats then convert them to numpy arrays # by multiplying with connectivity matrix if isinstance(J_eq, (int, float)): - self.J_eq = J_eq * connectivity - self.J_ax = J_ax * connectivity + self.J_eq = J_eq * adjacency + self.J_ax = J_ax * adjacency self.mu = mu * np.ones(self.n_sites) else: - raise TypeError("Connectivity matrix is provided, " + raise TypeError("adjacency matrix is provided, " "J_eq, J_ax, and mu should be floats") else: if isinstance(J_eq, np.ndarray) and \ @@ -529,9 +572,9 @@ def generate_one_body_integral(self, return self.one_body.todense() if dense else self.one_body def generate_two_body_integral(self, - sym=1, + basis='spinorbital basis', dense=False, - basis='spinorbital basis'): + sym=1): r"""Generate two body integral in spatial or spinorbital basis. Parameters @@ -598,7 +641,7 @@ class HamIsing(HamHeisenberg): def __init__(self, mu: np.ndarray, J_ax: np.ndarray, - connectivity: np.ndarray = None + adjacency: np.ndarray = None ): r"""Initialize XXZ Heisenberg Hamiltonian. @@ -608,8 +651,8 @@ def __init__(self, Zeeman term J_ax: np.ndarray J axial term - connectivity: np.ndarray - symmetric numpy array that specifies the connectivity between sites + adjacency: np.ndarray + symmetric numpy array that specifies the adjacency between sites Notes ----- @@ -632,7 +675,7 @@ def __init__(self, mu=mu, J_eq=J_eq, J_ax=J_ax, - connectivity=connectivity + adjacency=adjacency ) @@ -642,7 +685,7 @@ class HamRG(HamHeisenberg): def __init__(self, mu: np.ndarray, J_eq: np.ndarray, - connectivity: np.ndarray = None + adjacency: np.ndarray = None ): r"""Initialize XXZ Heisenberg Hamiltonian. @@ -652,7 +695,7 @@ def __init__(self, Zeeman term J_eq: np.ndarray J equatorial term - connectivity: np.ndarray + adajcency: np.ndarray Notes ----- @@ -675,5 +718,172 @@ def __init__(self, mu=mu, J_eq=J_eq, J_ax=J_ax, - connectivity=connectivity + adjacency=adjacency ) + + +class HamTJUV(HamPPP, HamHeisenberg): + r"""t-J-U-V Hamiltonian.""" + + def __init__(self, + connectivity=None, + adjacency=None, + alpha=-0.414, + beta=-0.0533, + u_onsite=None, + gamma=None, + charges=None, + sym=1, + atom_dictionary=None, + bond_dictionary=None, + orbital_overlap=None, + affinity_dct=None, + Rxy_matrix=None, + J_eq=None, + J_ax=None): + r""" + Initialize t-J-U-V Hamiltonian. + + Parameters + ---------- + connectivity: list + List of tuples specifying sites and bonds + between sites for the TJUV model. + adjacency: np.ndarray + Symmetric numpy array that specifies the adjacency between sites. + alpha: float + Specifies the site energy if all sites are equivalent. + Default value is the 2p-pi orbital of Carbon. + beta: float + Specifies the resonance energy, hopping term, if all bonds are + equivalent. The default value is appropriate for a pi-bond between + Carbon atoms. + u_onsite: np.ndarray + On-site Coulomb interaction; 1d np.ndarray. + gamma: np.ndarray + Parameter that specifies long-range Coulomb interaction; + 2d np.ndarray. + charges: np.ndarray + Charges on sites; 1d np.ndarray. + sym: int + Symmetry of the Hamiltonian: int [1, 2, 4, 8]. Default is 1. + atom_dictionary: dict + Dictionary of atom types and their properties. + bond_dictionary: dict + Dictionary of bond types and their properties. + orbital_overlap: np.ndarray + Overlap matrix for orbitals. + affinity_dct: dict + Affinity dictionary for the system. + Rxy_matrix: matrix + List of coordinates or positions in the system. + mu: np.ndarray + Zeeman term for the Heisenberg model. + J_eq: np.ndarray + J equatorial term for the Heisenberg model. + J_ax: np.ndarray + J axial term for the Heisenberg model. + + Notes + ----- + The Hamiltonian is given by: + .. math:: + \begin{align} + \hat{H}_{\mathrm{tJUV}} &= + \sum_{pq} h_{pq} a_p^{\dagger} a_q \\ + &+ \sum_p U_p \hat{n}_{p\alpha} \hat{n}_{p\beta} \\ + &+ \frac{1}{2} \sum_{p \ne q} \gamma_{pq}\\ + &+ \left( \hat{n}_{p\alpha} + \hat{n}_{p\beta} - Q_p \right)\\ + &+ \left( \hat{n}_{q\alpha} + \hat{n}_{q\beta} - Q_q \right) \\ + &+ \sum_{pq} \left[ J_{pq}^{\text{ax}} S_p^Z S_q^Z +\\ + &+ J_{pq}^{\text{eq}} \left( S_p^X S_q^X + S_p^Y S_q^Y \right)\\ + &+ \right] + \end{align} + + """ + # Default charges to an array of ones if not provided + if charges is None: + charges = np.ones(len(connectivity)) + + # Initialize the PPP part + self.ocupation_part = HamPPP(connectivity=connectivity, + adjacency=adjacency, + alpha=alpha, + beta=beta, + u_onsite=u_onsite, + gamma=gamma, + charges=charges, + sym=sym, + atom_dictionary=atom_dictionary, + bond_dictionary=bond_dictionary, + orbital_overlap=orbital_overlap, + affinity_dct=affinity_dct, + Rxy_matrix=Rxy_matrix) + + adjacency = np.asarray( + self.ocupation_part.adjacency) + + mu = np.zeros(adjacency.shape[0]) + + # Initialize the Heisenberg part + self.spin_part = HamHeisenberg(mu=mu, + J_eq=J_eq, + J_ax=J_ax, + adjacency=adjacency) + + def generate_zero_body_integral(self): + r"""Generate zero body integral. + + Returns + ------- + float + """ + self.zero_energy = self.ocupation_part.generate_zero_body_integral( + ) + self.spin_part.generate_zero_body_integral() + return self.zero_energy + + def generate_one_body_integral(self, basis: str, dense: bool): + r""" + Generate one body integral in spatial or spin orbital basis. + + Parameters + ---------- + basis: str + ['spatial', 'spin orbital'] + dense: bool + Dense or sparse matrix; default False + + Returns + ------- + scipy.sparse.csr_matrix or np.ndarray + """ + one_body_ppp = self.ocupation_part.generate_one_body_integral( + basis, dense) + one_body_heisenberg = self.spin_part.generate_one_body_integral( + dense, basis) + self.one_body = one_body_ppp + one_body_heisenberg + return self.one_body + + def generate_two_body_integral(self, basis: str, dense: bool, sym=1): + r""" + Generate two body integral in spatial or spin orbital basis. + + Parameters + ---------- + basis: str + ['spatial', 'spin orbital'] + dense: bool + Dense or sparse matrix; default False + sym: int + Symmetry -- [2, 4, 8] default is 1 + + Returns + ------- + scipy.sparse.csr_matrix or np.ndarray + """ + two_body_ppp = self.ocupation_part.generate_two_body_integral( + basis, dense, sym) + two_body_heisenberg = self.spin_part.generate_two_body_integral( + basis, dense, sym) + self.two_body = two_body_ppp + two_body_heisenberg + return self.two_body diff --git a/moha/rauk/PariserParr.py b/moha/rauk/PariserParr.py index 8d5d491..ba6cce7 100644 --- a/moha/rauk/PariserParr.py +++ b/moha/rauk/PariserParr.py @@ -11,7 +11,7 @@ from moha.rauk.rauk import build_one_body -from moha.rauk.utils import get_atom_type +from moha.rauk.utils import get_atom_type, get_atoms_list import numpy as np from scipy.special import gamma @@ -252,12 +252,14 @@ def compute_overlap( bond_dictionary[bond_key_reverse] = beta_xy else: for tpl in connectivity: - atom1, atom2, dist = tpl[0], tpl[1], tpl[2] - atom1_name, _ = get_atom_type(atom1) - atom2_name, _ = get_atom_type(atom2) + atom1, atom2, dist = tpl + atom1_name, site1 = get_atom_type(atom1) + atom2_name, site2 = get_atom_type(atom2) bond_key_forward = ','.join([atom1_name, atom2_name]) bond_key_reverse = ','.join([atom2_name, atom1_name]) - Sxy = orbital_overlap[bond_key_forward] + site1, site2 = site1 - 1, site2 - 1 + + Sxy = orbital_overlap[site1, site2] beta_xy = populate_PP_dct( dist, atom1_name, atom2_name, ionization, Sxy) @@ -272,20 +274,135 @@ def compute_overlap( return one_body -def calculate_gamma(Uxy_bar, Rxy): +def compute_gamma(connectivity, U_xy=None, Rxy_matrix=None, + atom_dictionary=None, affinity_dct=None, atom_types=None): + r""" + Calculate the gamma values for each pair of sites. + + Parameters + ---------- + u_onsite (list of float): List of potential energies for sites. + Rxy_list (list of float): List of distances + connectivity (list of tuples): Each tuple contains indices of two + sites (atom1, atom2), indicating that a + gamma value should be computed for + this pair. + atom_types (list of str): List of atom types in the lattice. + + Returns + ---------- + np.ndarray: A matrix of computed gamma values for each pair. + Non-connected pairs have a gamma value of zero. """ - Calculate the gamma value based on Uxy and Rxy. + if Rxy_matrix is None and U_xy is None: + U_xy, Rxy_matrix = compute_u( + connectivity, atom_dictionary, affinity_dct, atom_types) + elif Rxy_matrix is None: + _, Rxy_matrix = compute_u( + connectivity, atom_dictionary, affinity_dct, atom_types) + num_sites = len(U_xy) + gamma_matrix = np.zeros((num_sites, num_sites)) + + for tpl in connectivity: + atom1, atom2 = tpl[0], tpl[1] + atom1_name, site1 = get_atom_type(atom1) + atom2_name, site2 = get_atom_type(atom2) + # Get_atom_type returns site index starting from 1 + # So we need to subtract 1 to get the correct index + site1, site2 = site1 - 1, site2 - 1 + Ux = U_xy[site1] + Uy = U_xy[site2] + Rxy = Rxy_matrix[site1][site2] + Uxy_bar = 0.5 * (Ux + Uy) + gamma = Uxy_bar / \ + (Uxy_bar * Rxy + np.exp(-0.5 * Uxy_bar**2 * Rxy**2)) + gamma_matrix[site1][site2] = gamma + gamma_matrix[site2][site1] = gamma # Ensure symmetry + + return gamma_matrix + + +def compute_u(connectivity, atom_dictionary, affinity_dictionary, atom_types): + r""" + Calculate the onsite potential energy (U) for each site. Parameters ---------- - Uxy_bar (float): Represents the potential energy - Rxy (float): Represents the distance or a related measure. + connectivity (list of tuples): Each tuple contains indices of two sites + (atom1, atom2) and the distance between them. + atom_dictionary (dict): Dictionary mapping atom types to their + ionization energies. + affinity_dictionary (dict): Dictionary mapping atom types to their + electron affinities. + atom_types (list of str): List of atom types in the lattice. Returns ---------- - float: Computed gamma value based on the given parameters. + tuple: Contains two elements: + 1. List of float: Onsite potential energies calculated + for each site based on their ionization energy + and electron affinity. + 2. List of float: Distances corresponding to each + connectivity tuple. """ - # Example formula, needs actual formula to be replaced here - # This is just a placeholder formula - gamma = Uxy_bar / (Uxy_bar * Rxy + np.exp(-1 / 2 * Uxy_bar**2 * Rxy ** 2)) - return gamma + if atom_dictionary is None: + atom_dictionary = {} + hx_dictionary_path = Path(__file__).parent / "hx_dictionary.json" + hx_dictionary = json.load(open(hx_dictionary_path, "rb")) + alpha_c = -0.414 # Value for sp2 orbital of Carbon atom. + beta_c = -0.0533 # Value for sp2 orbitals of Carbon atom. + for key, value in hx_dictionary.items(): + hx_value = value * abs(beta_c) + atom_dictionary[key] = alpha_c + hx_value + + if affinity_dictionary is None: + affinity_path = Path(__file__).parent / "affinity.json" + affinity_dictionary = json.load(open(affinity_path, "rb")) + if connectivity is not None: + unique_atoms = {atom for tpl in connectivity for atom in tpl[:2]} + else: + unique_atoms = atom_types + u_onsite = [] + for atom in unique_atoms: + if atom in atom_dictionary: + ionization = atom_dictionary[atom] + affinity = affinity_dictionary[atom] + U_x = ionization - affinity + u_onsite.append(U_x) + Rxy_matrix = np.zeros((len(unique_atoms), len(unique_atoms))) + return u_onsite, Rxy_matrix + + num_sites = len(unique_atoms) + + u_onsite = [] + Rxy_matrix = np.zeros((num_sites, num_sites)) + + i = 0 + for tpl in connectivity: + atom1, atom2, dist = tpl[0], tpl[1], tpl[2] + atom1_name, site1 = get_atom_type(atom1) + atom2_name, site2 = get_atom_type(atom2) + # Get_atom_type returns site index starting from 1 + # So we need to subtract 1 to get the correct index + site1, site2 = site1 - 1, site2 - 1 + + ionization = atom_dictionary[atom1_name] + affinity = affinity_dictionary[atom1_name] + U_x = ionization - affinity + + u_onsite.append(U_x) + + Rxy_matrix[site1][site2] = dist + Rxy_matrix[site2][site1] = dist # Ensure symmetry + atoms_sites_lst = get_atoms_list(connectivity) + max_site = max([site for _, site in atoms_sites_lst]) + if len(connectivity) != max_site: + tpl = connectivity[-1] + atom_name, _ = get_atom_type(tpl[1]) + ionization = atom_dictionary[atom_name] + affinity = affinity_dictionary[atom_name] + U_x = ionization - affinity + + u_onsite.append(U_x) + + return u_onsite, Rxy_matrix diff --git a/moha/rauk/__init__.py b/moha/rauk/__init__.py index 7227592..183cb81 100644 --- a/moha/rauk/__init__.py +++ b/moha/rauk/__init__.py @@ -1 +1,2 @@ r"""Model Hamiltonian rauk module.""" +from .utils import * diff --git a/moha/rauk/affinity.Json b/moha/rauk/affinity.Json new file mode 100644 index 0000000..90a1588 --- /dev/null +++ b/moha/rauk/affinity.Json @@ -0,0 +1,92 @@ +{ + "_comment": "Values representing affinity. Data was taken from https://en.wikipedia.org/wiki/Electron_affinity_(data_page) on 06/10/2024", + "H": 0.754195, + "He": -0.5, + "Li": 0.618049, + "Be": -0.5, + "B": 0.279723, + "C": 1.2621226, + "N": -0.07, + "O": 1.46111297, + "F": 3.4011898, + "Ne": -1.2, + "Na": 0.547926, + "Mg": -0.4, + "Al": 0.43283, + "Si": 1.3895212, + "P": 0.746609, + "S": 2.0771042, + "Cl": 3.612725, + "Ar": -1.0, + "K": 0.501459, + "Ca": 0.02455, + "Sc": 0.17938, + "Ti": 0.07554, + "V": 0.52766, + "Cr": 0.675928, + "Mn": -0.5, + "Fe": 0.153236, + "Co": 0.662255, + "Ni": 1.15716, + "Cu": 1.23578, + "Zn": -0.6, + "Ga": 0.301166, + "Ge": 1.2326764, + "As": 0.8048, + "Se": 2.0206047, + "Br": 3.363588, + "Kr": -1.0, + "Rb": 0.485916, + "Sr": 0.05206, + "Y": 0.31129, + "Zr": 0.43328, + "Nb": 0.91740, + "Mo": 0.74723, + "Tc": 0.55, + "Ru": 1.04627, + "Rh": 1.14289, + "Pd": 0.56214, + "Ag": 1.30447, + "Cd": -0.7, + "In": 0.38392, + "Sn": 1.1120702, + "Sb": 1.047401, + "Te": 1.9708757, + "I": 3.059052, + "Xe": -0.8, + "Cs": 0.4715983, + "Ba": 0.14462, + "La": 0.557546, + "Ce": 0.600160, + "Pr": 0.10923, + "Nd": 0.09749, + "Pm": 0.129, + "Sm": 0.162, + "Eu": 0.116, + "Gd": 0.212, + "Tb": 0.13131, + "Dy": 0.015, + "Ho": 0.338, + "Er": 0.312, + "Tm": 1.029, + "Yb": -0.02, + "Lu": 0.2388, + "Hf": 0.1780, + "Ta": 0.328859, + "W": 0.81626, + "Re": 0.060396, + "Os": 1.077661, + "Ir": 1.564057, + "Pt": 2.12510, + "Au": 2.308610, + "Hg": -0.5, + "Tl": 0.320053, + "Pb": 0.356721, + "Bi": 0.942362, + "Po": 1.4, + "At": 2.41578, + "Rn": -0.7, + "Fr": 0.486, + "Ra": 0.1, + "Ac": 0.35 +} diff --git a/moha/rauk/affinity.json b/moha/rauk/affinity.json new file mode 100644 index 0000000..90a1588 --- /dev/null +++ b/moha/rauk/affinity.json @@ -0,0 +1,92 @@ +{ + "_comment": "Values representing affinity. Data was taken from https://en.wikipedia.org/wiki/Electron_affinity_(data_page) on 06/10/2024", + "H": 0.754195, + "He": -0.5, + "Li": 0.618049, + "Be": -0.5, + "B": 0.279723, + "C": 1.2621226, + "N": -0.07, + "O": 1.46111297, + "F": 3.4011898, + "Ne": -1.2, + "Na": 0.547926, + "Mg": -0.4, + "Al": 0.43283, + "Si": 1.3895212, + "P": 0.746609, + "S": 2.0771042, + "Cl": 3.612725, + "Ar": -1.0, + "K": 0.501459, + "Ca": 0.02455, + "Sc": 0.17938, + "Ti": 0.07554, + "V": 0.52766, + "Cr": 0.675928, + "Mn": -0.5, + "Fe": 0.153236, + "Co": 0.662255, + "Ni": 1.15716, + "Cu": 1.23578, + "Zn": -0.6, + "Ga": 0.301166, + "Ge": 1.2326764, + "As": 0.8048, + "Se": 2.0206047, + "Br": 3.363588, + "Kr": -1.0, + "Rb": 0.485916, + "Sr": 0.05206, + "Y": 0.31129, + "Zr": 0.43328, + "Nb": 0.91740, + "Mo": 0.74723, + "Tc": 0.55, + "Ru": 1.04627, + "Rh": 1.14289, + "Pd": 0.56214, + "Ag": 1.30447, + "Cd": -0.7, + "In": 0.38392, + "Sn": 1.1120702, + "Sb": 1.047401, + "Te": 1.9708757, + "I": 3.059052, + "Xe": -0.8, + "Cs": 0.4715983, + "Ba": 0.14462, + "La": 0.557546, + "Ce": 0.600160, + "Pr": 0.10923, + "Nd": 0.09749, + "Pm": 0.129, + "Sm": 0.162, + "Eu": 0.116, + "Gd": 0.212, + "Tb": 0.13131, + "Dy": 0.015, + "Ho": 0.338, + "Er": 0.312, + "Tm": 1.029, + "Yb": -0.02, + "Lu": 0.2388, + "Hf": 0.1780, + "Ta": 0.328859, + "W": 0.81626, + "Re": 0.060396, + "Os": 1.077661, + "Ir": 1.564057, + "Pt": 2.12510, + "Au": 2.308610, + "Hg": -0.5, + "Tl": 0.320053, + "Pb": 0.356721, + "Bi": 0.942362, + "Po": 1.4, + "At": 2.41578, + "Rn": -0.7, + "Fr": 0.486, + "Ra": 0.1, + "Ac": 0.35 +} diff --git a/moha/rauk/rauk.py b/moha/rauk/rauk.py index 7d36810..958dd43 100644 --- a/moha/rauk/rauk.py +++ b/moha/rauk/rauk.py @@ -3,6 +3,7 @@ from moha.rauk.utils import get_atom_type, get_atoms_list + from pathlib import Path from scipy.sparse import csr_matrix, diags @@ -144,6 +145,13 @@ def assign_rauk_parameters( # Ensure symmetry bond_dictionary[','.join([atom2_name, atom1_name])] = bond_value + else: + # Ensure symmetry in the bond dictionary + for key in list(bond_dictionary.keys()): + atom1_name, atom2_name = key.split(',') + reverse_key = ','.join([atom2_name, atom1_name]) + bond_dictionary[reverse_key] = bond_dictionary[key] + one_body = build_one_body( connectivity, atom_dictionary, diff --git a/moha/rauk/utils.py b/moha/rauk/utils.py index 9708f8d..2b80c25 100644 --- a/moha/rauk/utils.py +++ b/moha/rauk/utils.py @@ -1,6 +1,56 @@ r"""Utils for Rauk Module.""" import re +import numpy as np +from scipy.sparse import csr_matrix + + +def parse_connectivity(system, atom_types=None): + r""" + + Parse the connectivity of the system given as list of tuples. + + Parameters + ---------- + system: list + list of tuples that specifies sites and bonds between them + For example, for a linear chain of 4 sites, the connectivity + can be specified as [(C1, C2, 1), (C2, C3, 1), (C3, C4, 1)] + + Returns + ------- + tuple: (list np.ndarray) + First element is a list of atoms in the order + they apperar in the lattice, + second element is matrix that corresponds to + the either distance matrix, + or adjacency matrix. + """ + atoms_sites_lst = get_atoms_list(system) + max_site = max([site for _, site in atoms_sites_lst]) + n_sites = max_site + + if atom_types is None: + # Initialize atom_types with None, and adjust size for 0-based + # indexing + atom_types = [None] * max_site + for atom, site in atoms_sites_lst: + # Adjust site index for 0-based array index + atom_types[site - 1] = atom + atom_types = atom_types + connectivity_mtrx = np.zeros((max_site, max_site)) + atoms_dist = [] + for tpl in system: + atom1, atom2, bond = tpl[0], tpl[1], tpl[2] + atom1_name, site1 = get_atom_type(atom1) + atom2_name, site2 = get_atom_type(atom2) + atoms_dist.append((atom1_name, atom2_name, bond)) + connectivity_mtrx[site1 - 1, site2 - 1] = bond + # numbering of sites starts from 1 + atoms_dist = atoms_dist + connectivity_mtrx = np.maximum(connectivity_mtrx, connectivity_mtrx.T) + connectivity_matrix = csr_matrix(connectivity_mtrx) + return atoms_sites_lst, connectivity_matrix, n_sites, atom_types def get_atom_type(atom): diff --git a/moha/test/test.py b/moha/test/test.py index 35f7684..f8f8518 100644 --- a/moha/test/test.py +++ b/moha/test/test.py @@ -3,6 +3,7 @@ import numpy as np from moha import * from numpy.testing import assert_allclose, assert_equal +from moha.rauk.PariserParr import compute_gamma def test_hub2(): @@ -12,8 +13,10 @@ def test_hub2(): Should return U=\frac{1}{2}\left[U-\sqrt{U^{2}+16 t^{2}}\right]$ numerical result is -1.561552812 """ - hubbard = HamHub([("C1", "C2", 1)], + system = [("C1", "C2", 1)] + hubbard = HamHub(system, alpha=0, beta=-1, u_onsite=np.array([1, 1]), sym=1) + ecore = hubbard.generate_zero_body_integral() h = hubbard.generate_one_body_integral(basis='spatial basis', dense=True) v = hubbard.generate_two_body_integral(sym=1, @@ -41,10 +44,17 @@ def test_hub4(): nsites = np.linspace(2, 8, 4).astype(int) for nsite in nsites: nelec = nsite // 2 - hubbard = HamPPP([(f"C{i}", f"C{i + 1}", 1) for i in range(1, nsite)] + - [(f"C{nsite}", f"C{1}", 1)], - alpha=0, beta=-1, - u_onsite=np.array([1 for i in range(nsite + 1)])) + hubbard = HamPPP([(f"C{i}", + f"C{i + 1}", + 1) for i in range(1, + nsite)] + [(f"C{nsite}", + f"C{1}", + 1)], + alpha=0, + beta=-1, + u_onsite=np.array([1 for i in range(nsite)]), + gamma=np.zeros((nsite, + nsite))) ecore = hubbard.generate_zero_body_integral() h = hubbard.generate_one_body_integral(basis='spatial basis', dense=True) @@ -73,8 +83,8 @@ def test_ethylene(): """ a = -11.26 b = -1.45 - hubbard = HamPPP([("C1", "C2", 1)], alpha=a, beta=b, - gamma=None, charges=None, sym=None) + hubbard = HamPPP([("C1", "C2", 1)], alpha=a, beta=b, gamma=np.zeros( + (2, 2)), charges=None, sym=None, u_onsite=[0, 0]) ecore = hubbard.generate_zero_body_integral() h = hubbard.generate_one_body_integral(basis='spinorbital basis', dense=True) @@ -106,12 +116,13 @@ def test_4(): """ a = -5 b = -0.5 - hubbard = HamPPP([("C1", "C2", 1), - ("C2", "C3", 1), - ("C3", "C4", 1), - ("C4", "C1", 1)], + system = [("C1", "C2", 1), + ("C2", "C3", 1), + ("C3", "C4", 1), + ("C4", "C1", 1)] + gamma = compute_gamma(system) + hubbard = HamPPP(system, gamma=gamma, alpha=a, beta=b) - atoms_sites_lst, _ = hubbard.generate_connectivity_matrix() ecore = hubbard.generate_zero_body_integral() h = hubbard.generate_one_body_integral(basis='spatial basis', dense=True) @@ -166,8 +177,13 @@ def test_api_input(): g_matrix = np.arange(36).reshape((norb, norb)) charges = np.ones(norb) - ham = HamPPP(connectivity, alpha=0., beta=-2.5, u_onsite=u_matrix, - gamma=g_matrix, charges=charges) + ham = HamPPP( + adjacency=connectivity, + alpha=0., + beta=-2.5, + u_onsite=u_matrix, + gamma=g_matrix, + charges=charges) h = ham.generate_one_body_integral(basis='spinorbital basis', dense=True) v = ham.generate_two_body_integral(sym=1, basis='spinorbital basis', @@ -190,7 +206,11 @@ def test_spin_spatial_conversion(): u_matrix = np.ones(norb) beta = -2.5 - ham = HamHub(connectivity, alpha=0., beta=beta, u_onsite=u_matrix) + ham = HamHub( + adjacency=connectivity, + alpha=0., + beta=beta, + u_onsite=u_matrix) h = ham.generate_one_body_integral(basis='spinorbital basis', dense=True) v = ham.generate_two_body_integral(sym=4, basis='spinorbital basis', diff --git a/moha/test/test_Heisenberg.py b/moha/test/test_Heisenberg.py index d779d5a..f440a04 100644 --- a/moha/test/test_Heisenberg.py +++ b/moha/test/test_Heisenberg.py @@ -8,15 +8,15 @@ def test_heisenberg_0(): n_sites = 8 J_xy = np.random.rand() J_z = np.random.rand() - connectivity = np.zeros((n_sites, n_sites)) + adjacency = np.zeros((n_sites, n_sites)) for i in range(n_sites): j = i+1 if j == n_sites: j = 0 - connectivity[i, j] = 1 - connectivity[j, i] = 1 + adjacency[i, j] = 1 + adjacency[j, i] = 1 - ham = HamHeisenberg(J_eq=J_xy, J_ax=J_z, mu=0, connectivity=connectivity) + ham = HamHeisenberg(J_eq=J_xy, J_ax=J_z, mu=0, adjacency=adjacency) e0 = ham.generate_zero_body_integral() assert_allclose(e0, J_z/4*n_sites) @@ -26,13 +26,13 @@ def test_heisenberg_1(): n_sites = 8 J_xy = np.random.rand() J_z = np.random.rand() - connectivity = np.zeros((n_sites, n_sites)) + adjacency = np.zeros((n_sites, n_sites)) for i in range(n_sites): j = i+1 if j == n_sites: j = 0 - connectivity[i, j] = 1 - connectivity[j, i] = 1 + adjacency[i, j] = 1 + adjacency[j, i] = 1 h_exact = np.zeros((n_sites, n_sites)) @@ -43,7 +43,7 @@ def test_heisenberg_1(): h_exact[i, i] = -J_z/2 - ham = HamHeisenberg(J_eq=J_xy, J_ax=J_z, mu=0, connectivity=connectivity) + ham = HamHeisenberg(J_eq=J_xy, J_ax=J_z, mu=0, adjacency=adjacency) h = ham.generate_one_body_integral(basis='spatial basis', dense=True) assert_allclose(h, h_exact) @@ -53,13 +53,13 @@ def test_heisenberg_2(): n_sites = 8 J_xy = np.random.rand() J_z = np.random.rand() - connectivity = np.zeros((n_sites, n_sites)) + adjacency = np.zeros((n_sites, n_sites)) for i in range(n_sites): j = i+1 if j == n_sites: j = 0 - connectivity[i, j] = 1 - connectivity[j, i] = 1 + adjacency[i, j] = 1 + adjacency[j, i] = 1 v_exact = np.zeros((n_sites, n_sites, n_sites, n_sites)) for i in range(n_sites): @@ -73,7 +73,7 @@ def test_heisenberg_2(): v_exact[j, i, j, i] = J_xy/2 v_exact[i, j, i, j] = J_xy/2 - ham = HamHeisenberg(J_eq=J_xy, J_ax=J_z, mu=0, connectivity=connectivity) + ham = HamHeisenberg(J_eq=J_xy, J_ax=J_z, mu=0, adjacency=adjacency) v = ham.generate_two_body_integral(basis='spatial basis', dense=True, sym=4) @@ -87,13 +87,13 @@ def test_heisenberg_2_spin(): n_sites = 8 J_xy = np.random.rand() J_z = np.random.rand() - connectivity = np.zeros((n_sites, n_sites)) + adjacency = np.zeros((n_sites, n_sites)) for i in range(n_sites): j = i+1 if j == n_sites: j = 0 - connectivity[i, j] = 1 - connectivity[j, i] = 1 + adjacency[i, j] = 1 + adjacency[j, i] = 1 v_exact = np.zeros((2*n_sites, 2*n_sites, 2*n_sites, 2*n_sites)) for i in range(n_sites): @@ -124,7 +124,7 @@ def test_heisenberg_2_spin(): ham = HamHeisenberg(J_eq=J_xy, J_ax=J_z, mu=0, - connectivity=connectivity) + adjacency=adjacency) v = ham.generate_two_body_integral(basis='spinorbital basis', dense=True, sym=4) @@ -139,13 +139,13 @@ def test_Ising(): n_sites = 8 J_xy = 0 J_z = np.random.rand() - connectivity = np.zeros((n_sites, n_sites)) + adjacency = np.zeros((n_sites, n_sites)) for i in range(n_sites): j = i+1 if j == n_sites: j = 0 - connectivity[i, j] = 1 - connectivity[j, i] = 1 + adjacency[i, j] = 1 + adjacency[j, i] = 1 v_exact = np.zeros((n_sites, n_sites, n_sites, n_sites)) for i in range(n_sites): @@ -159,7 +159,7 @@ def test_Ising(): v_exact[j, i, j, i] = J_xy/2 v_exact[i, j, i, j] = J_xy/2 - ham = HamIsing(J_ax=J_z, mu=0, connectivity=connectivity) + ham = HamIsing(J_ax=J_z, mu=0, adjacency=adjacency) v = ham.generate_two_body_integral(basis='spatial basis', dense=True, sym=4) @@ -172,13 +172,13 @@ def test_RG(): n_sites = 8 J_xy = np.random.rand() J_z = 0 - connectivity = np.zeros((n_sites, n_sites)) + adjacency = np.zeros((n_sites, n_sites)) for i in range(n_sites): j = i+1 if j == n_sites: j = 0 - connectivity[i, j] = 1 - connectivity[j, i] = 1 + adjacency[i, j] = 1 + adjacency[j, i] = 1 v_exact = np.zeros((n_sites, n_sites, n_sites, n_sites)) for i in range(n_sites): @@ -192,7 +192,7 @@ def test_RG(): v_exact[j, i, j, i] = J_xy/2 v_exact[i, j, i, j] = J_xy/2 - ham = HamRG(J_eq=J_xy, mu=0, connectivity=connectivity) + ham = HamRG(J_eq=J_xy, mu=0, adjacency=adjacency) v = ham.generate_two_body_integral(basis='spatial basis', dense=True, sym=4) diff --git a/moha/test/test_tjuv.py b/moha/test/test_tjuv.py new file mode 100644 index 0000000..4f3543d --- /dev/null +++ b/moha/test/test_tjuv.py @@ -0,0 +1,179 @@ +"""Testing t-JUV model.""" + +import numpy as np +from moha import * +from numpy.testing import assert_allclose +from moha.rauk.PariserParr import compute_gamma + + +def test_tjuv_consistency_zero_body(): + r""" + Checking consistency of TJUV model + with Heisenberg and PPP model. + """ + adjacency = np.array([[0, 1, 0, 0, 0, 1], + [1, 0, 1, 0, 0, 0], + [0, 1, 0, 1, 0, 0], + [0, 0, 1, 0, 1, 0], + [0, 0, 0, 1, 0, 1], + [1, 0, 0, 0, 1, 0]]) + alpha = 0.0 + beta = -1.0 + u_onsite = np.array([1, 1, 1, 1, 1, 1]) + gamma = np.zeros((adjacency.shape[0], adjacency.shape[1])) + charges = np.ones(adjacency.shape[0]) + sym = 8 + J_eq = 1 + J_ax = 1 + + # Initialize the HamTJUV object + tjuv_hamiltonian = HamTJUV(adjacency=adjacency, + alpha=alpha, + beta=beta, + u_onsite=u_onsite, + gamma=gamma, + charges=charges, + sym=sym, + J_eq=J_eq, + J_ax=J_ax) + + # Generate the zero body integral + tjuv_zero = tjuv_hamiltonian.generate_zero_body_integral() + + heisenberg = HamHeisenberg( + J_eq=J_eq, + J_ax=J_ax, + mu=0, + adjacency=adjacency) + heisenberg_zero = heisenberg.generate_zero_body_integral() + + hpp = HamPPP( + adjacency=adjacency, + alpha=alpha, + beta=beta, + gamma=gamma, + charges=charges, + sym=None, + u_onsite=u_onsite) + hpp_zero = hpp.generate_zero_body_integral() + + assert_allclose(tjuv_zero, heisenberg_zero + hpp_zero) + + +def test_tjuv_consistency_one_body(): + r""" + Checking consistency of TJUV model + with Heisenberg and PPP model for one-body term. + """ + adjacency = np.array([[0, 1, 0, 0, 0, 1], + [1, 0, 1, 0, 0, 0], + [0, 1, 0, 1, 0, 0], + [0, 0, 1, 0, 1, 0], + [0, 0, 0, 1, 0, 1], + [1, 0, 0, 0, 1, 0]]) + alpha = 0.0 + beta = -1.0 + u_onsite = np.array([1, 1, 1, 1, 1, 1]) + gamma = np.zeros((adjacency.shape[0], adjacency.shape[1])) + charges = np.ones(adjacency.shape[0]) + sym = 8 + J_eq = 1 + J_ax = 1 + + # Initialize the HamTJUV object + tjuv_hamiltonian = HamTJUV(adjacency=adjacency, + alpha=alpha, + beta=beta, + u_onsite=u_onsite, + gamma=gamma, + charges=charges, + sym=sym, + J_eq=J_eq, + J_ax=J_ax) + + # Generate the one-body integral + tjuv_one_body = tjuv_hamiltonian.generate_one_body_integral( + basis='spatial basis', dense=True) + + heisenberg = HamHeisenberg( + J_eq=J_eq, + J_ax=J_ax, + mu=0, + adjacency=adjacency) + heisenberg_one_body = heisenberg.generate_one_body_integral( + basis='spatial basis', dense=True) + + hpp = HamPPP( + adjacency=adjacency, + alpha=alpha, + beta=beta, + gamma=gamma, + charges=None, + sym=None, + u_onsite=u_onsite) + hpp_one_body = hpp.generate_one_body_integral( + basis='spatial basis', dense=True) + + # Assert that the TJUV one-body integral is close to the sum of Heisenberg + # and PPP one-body integrals + assert_allclose(tjuv_one_body, heisenberg_one_body + hpp_one_body) + + +def test_tjuv_consistency_two_body(): + r""" + Checking consistency of TJUV model + with Heisenberg and PPP model for two-body term. + """ + adjacency = np.array([[0, 1, 0, 0, 0, 1], + [1, 0, 1, 0, 0, 0], + [0, 1, 0, 1, 0, 0], + [0, 0, 1, 0, 1, 0], + [0, 0, 0, 1, 0, 1], + [1, 0, 0, 0, 1, 0]]) + alpha = 0.0 + beta = -1.0 + u_onsite = np.array([1, 1, 1, 1, 1, 1]) + gamma = np.zeros((adjacency.shape[0], adjacency.shape[1])) + charges = 1 + sym = 8 # Use an integer value for symmetry + + J_eq = 1 + J_ax = 1 + + # Initialize the HamTJUV object + tjuv_hamiltonian = HamTJUV(adjacency=adjacency, + alpha=alpha, + beta=beta, + u_onsite=u_onsite, + gamma=gamma, + charges=charges, + sym=sym, + J_eq=J_eq, + J_ax=J_ax) + + # Generate the two-body integral + tjuv_two_body = tjuv_hamiltonian.generate_two_body_integral( + basis='spatial basis', dense=True, sym=sym) + + heisenberg = HamHeisenberg( + J_eq=J_eq, + J_ax=J_ax, + mu=0, + adjacency=adjacency) + heisenberg_two_body = heisenberg.generate_two_body_integral( + basis='spatial basis', dense=True, sym=sym) + + hpp = HamPPP( + adjacency=adjacency, + alpha=alpha, + beta=beta, + gamma=gamma, + charges=None, + sym=sym, + u_onsite=u_onsite) + hpp_two_body = hpp.generate_two_body_integral( + basis='spatial basis', dense=True, sym=sym) + + # Assert that the TJUV two-body integral is close to the sum of Heisenberg + # and PPP two-body integrals + assert_allclose(tjuv_two_body, heisenberg_two_body + hpp_two_body)