From 8642b10fe44ca3771bebbdd6d57b8a639669fea6 Mon Sep 17 00:00:00 2001 From: Diego Prada-Gracia Date: Fri, 20 Dec 2024 20:22:37 -0600 Subject: [PATCH] Some tests added --- .../structure/align_principal_axes.ipynb | 116 +++----- .../user/tools/structure/get_angles.ipynb | 247 ++++++++++++++++++ docs/contents/user/tools/structure/index.md | 4 +- .../_private/digestion/argument/triplets.py | 21 ++ molsysmt/structure/get_angles.py | 5 +- ...st_align_principal_axes_molsysmt_MolSys.py | 37 +++ .../test_get_angles_from_molsysmt_MolSys.py | 32 +++ 7 files changed, 375 insertions(+), 87 deletions(-) create mode 100644 docs/contents/user/tools/structure/get_angles.ipynb create mode 100644 molsysmt/_private/digestion/argument/triplets.py create mode 100644 tests/structure/align_principal_axes/test_align_principal_axes_molsysmt_MolSys.py create mode 100644 tests/structure/get_angles/test_get_angles_from_molsysmt_MolSys.py diff --git a/docs/contents/user/tools/structure/align_principal_axes.ipynb b/docs/contents/user/tools/structure/align_principal_axes.ipynb index d624d565f..35b9992e2 100644 --- a/docs/contents/user/tools/structure/align_principal_axes.ipynb +++ b/docs/contents/user/tools/structure/align_principal_axes.ipynb @@ -20,7 +20,7 @@ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "acabf3aa42d247ada3e3978689397dd7", + "model_id": "54b3043c3363436e8c477185c720171e", "version_major": 2, "version_minor": 0 }, @@ -76,37 +76,37 @@ "text/html": [ "\n", - "\n", + "
\n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", " \n", " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", - " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", " \n", " \n", "
formn_atomsn_groupsn_componentsn_chainsn_moleculesn_entitiesn_lipidsn_structuresformn_atomsn_groupsn_componentsn_chainsn_moleculesn_entitiesn_lipidsn_structures
molsysmt.MolSys1341111111molsysmt.MolSys1341111111
\n" ], "text/plain": [ - "" + "" ] }, "execution_count": 5, @@ -130,20 +130,10 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "922d611e-686b-46f3-a603-bed5a9adffb4", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[ 0.09213962 -0.02624376 -0.9954002 ] 13.378581090770417\n", - "[0.56239734 0.82631277 0.03027277] 73.67866342245162\n", - "[-0.82171742 0.56259975 -0.09089556] 78.32952437971095\n" - ] - } - ], + "outputs": [], "source": [ "for ii in range(3):\n", " print(axes[0,ii], momenta[0,ii])" @@ -151,32 +141,17 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "095cb193-a067-4940-8a4d-2c0ca5d34d80", "metadata": {}, - "outputs": [ - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "877c7f96fcaa43188a75af4eefe83d19", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - "NGLWidget()" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "msm.view(molsys)" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "id": "b46c0c26-d4be-41a5-a9fb-c9b97dd5bad2", "metadata": {}, "outputs": [], @@ -186,7 +161,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "id": "5d29c2f5-b9ab-494e-bf26-a71de19c00c2", "metadata": {}, "outputs": [], @@ -196,20 +171,10 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "id": "2f4e8bdb-7f79-4c72-883b-536f7ba02ee0", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[1.00000000e+00 2.87221485e-16 8.58936055e-17] 13.37858109077042\n", - "[-2.87221485e-16 1.00000000e+00 -3.60822483e-15] 73.67866342245169\n", - "[-8.58936055e-17 3.60822483e-15 1.00000000e+00] 78.32952437971095\n" - ] - } - ], + "outputs": [], "source": [ "for ii in range(3):\n", " print(axes[0,ii], momenta[0,ii])" @@ -217,25 +182,10 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "id": "40ec81fa-54be-4eee-8073-d7f07af8beaf", "metadata": {}, - "outputs": [ - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "0fdae2ea503a4798ba16c3d433d71084", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - "NGLWidget()" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "msm.view(molsys_2)" ] @@ -243,7 +193,7 @@ { "cell_type": "code", "execution_count": null, - "id": "39800e8d-836c-4ee6-b6c4-c05c9d139149", + "id": "485d758d-833b-425f-ad77-552d37ad421e", "metadata": {}, "outputs": [], "source": [] @@ -265,7 +215,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.11.10" } }, "nbformat": 4, diff --git a/docs/contents/user/tools/structure/get_angles.ipynb b/docs/contents/user/tools/structure/get_angles.ipynb new file mode 100644 index 000000000..d239db947 --- /dev/null +++ b/docs/contents/user/tools/structure/get_angles.ipynb @@ -0,0 +1,247 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "984ae1bb553740ac8f53d7d8abbd52cb", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import molsysmt as msm\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Get angles" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## get_angles" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "molecular_system = msm.systems['pentalanine']['traj_pentalanine.h5']\n", + "molecular_system = msm.convert(molecular_system, to_form='molsysmt.MolSys')" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
indexidnametypegroup indexgroup idgroup namegroup typecomponent indexchain indexmolecule indexmolecule typeentity indexentity name
00H1H01ACEterminal capping000peptide0peptide 0
11CH3C01ACEterminal capping000peptide0peptide 0
22H2H01ACEterminal capping000peptide0peptide 0
\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "msm.info(molecular_system, element='atom', selection=[0,1,2])" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "5000" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "msm.get(molecular_system, n_structures=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "angles = msm.structure.get_angles(molecular_system, [0,1,2])" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(5000, 1)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "angles.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/GU6VOAAAACXBIWXMAAA9hAAAPYQGoP6dpAABgv0lEQVR4nO3dd3wUZf4H8M+mh5gEQiAFAgaIoIQSQLoQBOmgoocIKjbUs3KAntGfAmfB484Kop56ooLi3YGIjd57CaGDlAABEmpIJwnJ/P4IWXaTLTO7MzvP7H7er9e+XsnO7MyzszPPfOepJkmSJBAREREJxE/vBBARERHVxACFiIiIhMMAhYiIiITDAIWIiIiEwwCFiIiIhMMAhYiIiITDAIWIiIiEwwCFiIiIhBOgdwJcUVlZiTNnziA8PBwmk0nv5BAREZEMkiShoKAA8fHx8PNzXEZiyADlzJkzSEhI0DsZRERE5IKsrCw0btzY4TqGDFDCw8MBVH3BiIgInVNDREREcuTn5yMhIcF8H3fEkAFKdbVOREQEAxQiIiKDkdM8g41kiYiISDgMUIiIiEg4DFCIiIhIOAxQiIiISDgMUIiIiEg4DFCIiIhIOAxQiIiISDgMUIiIiEg4DFCIiIhIOAxQiIiISDgMUIiIiEg4DFCIiIhIOAxQiEhIB3Py8cW6Yyi7Wql3UohIB4aczZiIvN/AD9YBAColCU/0aq5zaojI01iCQkRC23M6X+8kEJEOGKAQERGRcBigEBERkXAYoBCR0CRJ0jsJRKQDBihEREQkHAYoREREJBzFAcratWsxbNgwxMfHw2QyYeHCheZl5eXl+Otf/4o2bdogLCwM8fHxeOihh3DmzBmrbaSmpsJkMlm9Ro0a5faXIZJDkiS8MG8n/r74oN5JISIiOxQHKEVFRWjXrh1mzpxZa1lxcTHS09Px2muvIT09HQsWLMAff/yB4cOH11p33LhxyM7ONr8+++wz174BkUL7zuTjp4wz+GT1Ub2TQkREdigeqG3QoEEYNGiQzWWRkZFYtmyZ1XszZsxA586dcfLkSTRp0sT8fp06dRAbG6t090RuK+XIpEREwtO8DUpeXh5MJhPq1q1r9f7cuXMRHR2N1q1bY9KkSSgoKLC7jdLSUuTn51u9iIiIyHtpOtT9lStX8PLLL2P06NGIiIgwvz9mzBgkJiYiNjYWe/fuRVpaGnbt2lWr9KXatGnTMHXqVC2TSj7EZNI7BaQEOxkT+SbNApTy8nKMGjUKlZWVmDVrltWycePGmf9OTk5GUlISOnXqhPT0dHTo0KHWttLS0jBhwgTz//n5+UhISNAq6URERKQzTQKU8vJyjBw5EpmZmVi5cqVV6YktHTp0QGBgIA4fPmwzQAkODkZwcLAWSSUfxAIUIiLxqR6gVAcnhw8fxqpVq1C/fn2nn9m3bx/Ky8sRFxendnKIiIjIgBQHKIWFhThy5Ij5/8zMTGRkZCAqKgrx8fG49957kZ6ejl9++QUVFRXIyckBAERFRSEoKAhHjx7F3LlzMXjwYERHR2P//v2YOHEiUlJS0KNHD/W+GRERERmW4gBl+/bt6NOnj/n/6rYhY8eOxZQpU7Bo0SIAQPv27a0+t2rVKqSmpiIoKAgrVqzAhx9+iMLCQiQkJGDIkCGYPHky/P393fgqREQkgv1n8lFSfhUdm0bpnRQyMMUBSmpqqsPJu5xN7JWQkIA1a9Yo3S15mY1HLuDbzScwdXhrNIwI0Ts5RKSiwR+tAwDs+L9+qH8D2w+SazTtZkxkz+gvtgAAyisq8cXYWz26bxP7GRN5xNn8UgYo5DJOFki6OnP5it5JINFxIBQin8QAhYiIiITDAMWASq9W6J0EQ/P1Cp684nL0+edqTOdszqQxicVf5AYGKAbzn21ZaPl/i/FTxmlFn9t5Mhe3v7saqw6e0yhlZBTfbDqOzAtFmMXZnIlIYAxQDOal+bsBAC/My1D0ubH/3opj54vwyOxtGqSKjKTCSU87IiIRMEDxEcVlrBayxVm3eCIi0gcDFLKr7GolKiu97wbOXsZEROJjgEI2lZRVoP3flmL4x+s13Y8Rw583ftmP299djcLSq3onhYgElHWpGKP+tQkrDpzVOymGxgCFbEo/mYvisgrsPZ2vd1I05UoNz5frM3HsfBH+tz1L/QR5gMlg/ZjYE4SM5uUFu7H52CU89vV2vZNiaAxQyOeodYM2au0Xb/jkKb7axOtiYZneSfAKDFDIJk9lLHo/y/to/klEJDwGKEQuYmNbIsd4jZA7GKAQERGpiBOSqoMBio9Qer34SjsFXxwHxWiNZMm4RLm8fPE69wYMUMjnqPVww9s8kfgem70Nd83a6JVjOnk7Big+QukDhK88ZbuTZdn6rCRJmLHiMFYeFHf8A6OVjvHhl9yx4uA57Mq6jD/OFWi6nxMXizBr9REUXCnXdD++JEDvBJCYPHUT87Z7z+pD5/Husj8AAMffGaJzaojIUwZ+sA4l5RU4fqFI76R4DZag+Ai22bLNnadzW4f0TF6J6xv0EF8pHSN9+Gp7j5LyqvnOtmZe0jkl3oMBCtnko3kMEREJggEKERGRh7z16358svqo3skwBAYogth45AL2ns7TOxkep3dlgzttbTjWAZFxeLJq096ejp4vxOfrMvH3xQc9lhYjY4CigqxLxZi58jDyil1rvX3mcglGf7EFQ2doO3MwVdEyrmD7DtKbr7YBccZTDf8dPbiUlFV4JA3OvPnLfvzlhwzhzxUGKCq4e9YG/HPpH3h5wW6XPn/msvgNK72V2tfn6cvF6m6QSIG0Bbtxx/trcaVcjBuhrxL7tg98sT4TP+48jSPnCvVOikMMUFRw4drMlZuOXdQ5JcaTW1yGKYv2Yd8Zz1VvaVnK8fEq1i3LJffpTfCHPKF8vzULR84VYsm+HN3SINLvpUcJgeilEpbKK8ROKwMUsslTp2123hXM3ngcQz6SV711pbwCqw6dE6aoVATL9p9FpsHGXsgrLkfPv6/C337er3dSvJKB7pEe46nq1+MXi1nRqxIGKCqylylUVkr4edcZfLD8DyzcedqzibrGW9pGTP5pHx75ahsm/CdD76QIYcORCxj3zXb0+edqvZOiyNytJ3D6cgn+vSFT76QQGdLivTn4dtNxvZOhKQYobrpQWOp0nf/uyMJz3+/EB8sPY/wPGawftmPcN9tx8qLjNhw/bM8CAPy+V78i7GoidOLJyLqs+DMipJtP+PatO3weUxbtc5pP5BaVcX4ZH/bUnB147ad9Vu1IftuTjQe/3CLrvmQEDFDclLZgj9N1Nh61bptyVYdMxQjzryzbfxbPfJfu0X364o3SF7+zkTz45VbM3ngcX204bnedXVmXkfLGMoz7ZrvN5SIEoYD+55re+/eES0Vl5r+fnpuOdYcvYNpv8roxi35fYIDipkM52k5AJYqsS8UeeVrLytW+F4wemfeV8gq8t/QQdrlQ4qGnnzJOY8hH65yWbCklyg1ULSsOnMXBnHxVt3nKwbUwe+Pxqv0ePGdzuS/cmPXyU8ZpbHbSIULv8/tycZnzlQyAAYqb5ESg5RWVmux7+uKDOHFRXuNIpW1QLFuiz99xCrdNX2W33cey/WfxlUHbEnjqCeKT1Ufx0cojuPPjDW5v6/c92dh01DM9xl6Yl4F9Z/Lx6kLnJYW+au/pPDz29XYM/GCd3knR3KT/7sLjX2932FPFm2OjA9n5eGFeBkb9a7PeSfEJDFBUZHnRbjx6ATNWHEZlpYTf9uTYXQ9w/YKetfoo7vlko4uflm/GysMAgIUZZ2wuH/fNdkz9eb9hRsLV4+lGrZK2rEvF+PPcdNz/uWczyMLSqx7dnyXRi6H/OKtNKao731qLc1ySJPxvxyksP3AWx2T2GtO7JEFtp3O1HbOq7GolcouUlX7YChblnjuid54I0DsB3mr051sAAAlRdTTdT/UYLCI47yUNs0R2Nv+K29vwtpuG3tSsTjFKo1e56dS7qkmy+lvMY5t1qRjpJ3MxtG08BnywFpkXirDh5dvRqG6o3knTneISlLVr12LYsGGIj4+HyWTCwoULzcvKy8vx17/+FW3atEFYWBji4+Px0EMP4cwZ6yfv0tJSPPfcc4iOjkZYWBiGDx+OU6dOuf1l9GZriOOTl2rXI9dcT437xaqD55C2YA97CCnkTgaq1n3+UE4B/vJDBo4LPJaJ2jca0Z/c9DL+hwzz346OkGhHL6+4XPi8Z+AH64ScA+e26avwwrwM/Hd7lnk8o5V22hbZYuu+Y0JVacyOE7m4qlETA09QHKAUFRWhXbt2mDlzZq1lxcXFSE9Px2uvvYb09HQsWLAAf/zxB4YPH2613vjx4/Hjjz9i3rx5WL9+PQoLCzF06FBUVIh9gtdUWSkh69L1Ij89RxB8ZPY2fL/1JL5cb8y2IN7OUanFiFkb8OPO03h09ja39lFUehVr/zjvtM2Tq6dphYBP91fKK3DnzPV4+7cDuqVBzRKpRbtsV6PWpPcvYbn/y8VlaPe3peg6bYVu6QGAL9Ydwxfrjjlcx91ZhLUsfdySecmlz9m777w8fzfu+WSjOShbvDcb034/YFX6JWqpUjXFVTyDBg3CoEGDbC6LjIzEsmXLrN6bMWMGOnfujJMnT6JJkybIy8vDl19+iW+//Rb9+vUDAMyZMwcJCQlYvnw5BgwY4MLX0EfN7sO22Dp31GqDYgvn9VHGrWOvUm5VdG1UXLn1+vY8/vV2bDp2EX9ObY6/DmylRtLMLhSWov3fluLO9vF48642qm7bHT/vOoNdp/Kw61QeXhl8sy5p0LsaQ287r/VMu2xjslRP3QALS6/izV+rgtQ/dUxAZJ1Aj+xXVBKABdcGBf18XSZeHXILnppTNYRDSkI9HVOmjOaNZPPy8mAymVC3bl0AwI4dO1BeXo7+/fub14mPj0dycjI2brTd4LO0tBT5+flWLxEUl+nXcFBrltmKo9k5jchbqxeq54L6futJh+u58nOeyi1BwZWrmLPZ8bY9TY8xhbRwUWH7LTk/4f92nMI7vx9UrWTX7mZqvO/K/q6UV+CX3Wdc7h5bYTGnzJWrFci8UITx83Zq1oDZGZGzTCO1FdS0keyVK1fw8ssvY/To0YiIiAAA5OTkICgoCPXqWUdxMTExyMmxPTrotGnTMHXqVC2TqiuBz2USQPrJXJy8WIy7UhrpnRTVqJWB55fUfmo3mnlbT+LlBXvwbJ8Wqm530n93AQBub9UQnROjVN22ZQziqJRE7sPA3xcfxFcbjuOWuAj89sJtitPjZ/GoXSlJeOjfW5B1qQRL9p1VvC1Pmu6hNjGW1bOzVh3xyD7VoFkJSnl5OUaNGoXKykrMmjXL6fqSJNl9Uk9LS0NeXp75lZWVpXZyNWPrK+lSIqF0ly48dKnxrTxdXK5lu6EdJy7h9z3Zbm9nxKyNGP9DBnaezFUhVd7jp4zTmPa7axn8f7Zn4e+L1StdcMfrP+0DAMzU6MZhhEG7fr7W9mZ/tmul434WeWqlBHPbwBKdGu7KOa32ncnDLIs2MZb5p7t5adlV63Zond683vQiO+96T0DRS5M1KUEpLy/HyJEjkZmZiZUrV5pLTwAgNjYWZWVlyM3NtSpFOXfuHLp3725ze8HBwQgODtYiqaqydU7uOyPI2CBOLhhJkrAz6zJaxoQjLNi100L/rF4eT8WH93yyCQCwfEJvVfZ54mIxGtdz3vWwelflFZXIybvisKt79Q1a5Go8e5n9VDdmQn7pf7sBAOEhAXio2424wcVzHhD/vNfit5XbtsTWejtP5sJkMqF9Ql3V0mP5FY3SVVvuDO6uKKvRUD7XRvsgQPxGsqqXoFQHJ4cPH8by5ctRv359q+UdO3ZEYGCgVWPa7Oxs7N27126AYmS2ihhdbSSrZS+K/+04hRGzNl4f+E2n+5Waeem03w7gz3N26D7q5WmVGi4rzUwe+GILbpu+Cmv+OG93neEzN+Chf291N2lu0bN30PTFh3C/zqOCanGTcDUoKSq9Knt0aleUlFXg7lkbcdfHGzTrllypYamY3MPqyuEXO1TQh+IApbCwEBkZGcjIyAAAZGZmIiMjAydPnsTVq1dx7733Yvv27Zg7dy4qKiqQk5ODnJwclJVVFTNGRkbisccew8SJE7FixQrs3LkTDzzwANq0aWPu1WNUWt/T5YzUau/CqBlR1zQ/vWocmoPVI57KuFryr5TjsEUjNNGewT9bewy/780x9zJQm2jft6bqbotzN5+wu86e03lYd/iCS9vPulSMOz/egF92W3eNXZB+CpuPXURJWYX5afZqRSUW7802z7Jqeex6TV9Vq0jak/ZYXFeVlRJenr8b31ybxv7kxWKnN1IRzwNXq656/n0lev9jtc15hey2kVWwK8sRiUvKtApQNNmssGx9XRHPSVcoLtfcvn07+vTpY/5/woQJAICxY8diypQpWLRoEQCgffv2Vp9btWoVUlNTAQDvv/8+AgICMHLkSJSUlKBv376YPXs2/P39XfwaYnD1utD6ZJKTEbhSF9njnZUouCJ+T6arFQ5KUDyUmYle1+uKVxfuxa6sy3j2u50Y2jYeQFW7mwn/2WVep23jSCx6tie+XJ+Jab8fRExEMLa8Yv0gcvpyCfaczkPHpsq6P7p6RB1N6rf28HnM21bVxi25USRGzNqIZg3CsHJiqnmd4rKrqBN0Pet09xSydw7qcZ+trgpYdfA8WsVG2F2vOs3zd5zCxP/usl7mwn7dvQ4tP69lCYqnCFzj6lGKA5TU1FTHReYyTo6QkBDMmDEDM2bMULp7YeSVlOOfSw+59FlP1/fvOnXZ6TqOkmRvkRbBiRZ5S81z0tPXvggNMUuvVuCZueno2SJatW0WXLGu1y6vqDS3u6m2+1RV6cTS/VVVnWfz9e/iWOjgvLU8pxddm3vq2PnrVR6frz2Gt347gI/uT8HwdvHaJdJNnshjikqv1gpORPCigGmSw9VfzJtjGU4W6KK//bwff5wtdOmzWg7UphbRG08p4eibZGRdxkP/3qrbeAlKuFMKM3/HaSw/cA5T3GhY6owa8wQpocmEeBZ/26oafOvaiLUTLIak14rDr+fku1tWS2lxA5MkYNgM9xp5frjiMH7KOO3y5ysqJZujJqefvOxGqtThjSWmemCA4qLdMkol3HGu4Iq5vt6XaDMLq/1lY/+9FWv/OI9HvnJvmHlPcbXxXZFF3b8eWWetUiwD5N+7NGq75Krl+8/izo834Nh55w9GL8/f7da+Fu06gywb84hZcnfk49kbj+OFeRkAlJ8PkiSh//tr0P2dlSivqFT0OFVcdlXz2bk9+YBna09+RrjAZGCA4iItT78r5RXo/NYKdHpzudVET3LOubP5pSi96rjNyc6TuRj3zXaryem84XyW2+jO1ndVq6eNlqYvPohfdjsfV8VZ8b4W5+7WzEuaVStIqPptd2VdVqW6zJXqTFtpMv/tZprkfvrxb7ZjV9Zl/EVG6Y1lQ9ETl4rx1q/7FZVwHcjOx23TV1m9Z1nKKLubsUYZZaUEHD1fhPMFpTYnZHXklteXIHnyEpcaZnuqZIQlMFUYoAig5ql4qej6wErFFkW1ci72ZfvPYtCH6xyuc/esjVi2/yye+Hb7te1K2HDE+bxC7iivqMSB7Hynmfnl4nKX5hP6dvMJ3Pz6YszfUXtWbK2eZjwd1J3Ju4KvNhxX/DkTnKfVlZus5Ue2HXdtojOLrTlcOuaLzbjz4w3mBqxaUXIU9GpblKdw9Nw3ftmPz9dl4tnv0l3e5/ELRU7zFXcoOZT5V8qxbP/14RtemLcT7/yufLJILUuo9Q4wvKWKngGKmgQ5Jywb9QH2nwqrR1usnsNFS3+ek45BH67D1xuPO123zz9Xm//OKy7H9uOXnN4MXlu4FwCEbLQHAKcUBF3/23EKL8/fLWt8EBEa4FZzlCW7m8rqdgXWAYp6N4HFe21Ps2FPRaWEuz7egIpK+yNgiybDhSqrykoJkiS59FmtjP33Vjw1Z4f5/72n84WbI0pLFZUSRn/ueOwevQMktTBAUVFB6VVcKCxVPLy5o8zb8v7jaj5of46vqiVnLlsX/VZUSubht91xNv8K1h0+D0mSsPxA1RPPF+sznX6u1KLotf8Ha3Dvp5vw2x7bN5Al+3Lwyo973E6r1pS0Z5j0312Yty2r1vgiNUmShH8u/cPNlLnOnfuyqxmoVtnuzFWHFW9/16k8ZKhU7WSLzLn5ZAvws5/d2yuVafbKb0hM+w1L9llff46+stYx8049G8EKcN/ffOwiNh5V56Hyvs8245lrJWtXnYyVpQdNJwv0RZ3eXK53Etz2655snMq1eOJ3oQHbvjP5GHqtlf9XD9/qclqqu6Uu3peDIW3jai1/8tsdtd6rnSCXd++QCSYcO1+IhhEhDodKd3X3zoryd52SN41Czf1rlce617ZDXqos96F1Lx6tPqEGVwK8AH/bn1mQfgoT/rML3ZrVt7kcAH5XWMKkhT2n8hBfN0TfRAhQYKnmDN6FpVfx6+5s/KnjOTz81TZMG9EG93duotr23cUSFBmyLhXXmt/BU0XramTC9jZRncnVXH6hwL262a83HjcHJwCw4YhrI5VaWrY/B0fOudatuzZ17my7si7j9nfXoHeNxoSeorQtghzlFZVOe2+4Qq2rRbPgSsdHY3t5idopCvS3nd1XD6zniapeV6WfzMWwmetx61vqPQAu3puj2fACLuXbcofRd2HTzlQ/6KUtEKs0mgGKE3M2n8Bt01fh1YXa/XAO6+4FiNiV+reDhpynckvwlJxSjxqulFei33trXEqP0kN49Hwh/vJDBo6cc5x5/bC9qj3ExSJtZot19bd3loE52uyof23GbdNXYb2M4e/VPTflbcyyvYfrGXXtTzq7ocxabXum4bKr7h8ET13iAtROuPzAVX0+qjmM/d9+2Y/+769V9iERDmINat8j/rNd24boSjBAseFKeYW5geJ7y6rq+L/fav2jHT2v3oRaosUgrmYichsLLt7nvLhYyxIqZ8l84Ist+HHnadz3mb6TyDkj9xjJ+VWqN7XjRC4A4PttyhsdOiqFqJlUl88xy226tgmXTF9se9To+z/f7FVD3cs11M1B2kRyUYPePO5mX46GipBz7bjTi6d6pm8RMECpIf9KOVq9ttjtURKNoPokrnnCu3pxqRlUPONGl8ialCYrO6+q0fDFojKUlFVgwPtrMWWR+42GRWGvd5BlNabSTHtr5iWM/2GnW+nKzivB7A2ZtQbRcuW0OpCdj19ljBlTTcAHY69Q86dzPW9xOyl2dXxzOVYcqD3rvJ7e/FV+t2lbAYvS4yVqMMxGsjVsvDYeyP5s+xOKqW31oXNW/xuk16Km7PXaUYOSw/vzrjM4dLYAh84WYMrw1pqlSUs1M58PVxy2ud6cLddnPd587JLVeDzOrPnjvCtJs3L3xxuRk38Fe8/k459/amdzHatGsg62VT1mR/QNXdHFQeNPNcg5n3ZlXcb+7Hx8vOoIHuuZiEd6JAKwHyw62q7o1b5GHIPjwxWH0ffmGKfrVTiYeNSSvTy8+uHH6X7crMvylnsIS1Bq8Xxj2I9XHa2xT813CcCikaxKJ7NRxoNQcngrBL4b5Bar2/ZlQbr1vCh7TjvuJSR35F5bikqtP/v91izkXBvp9H82BturprQx6yG5jSDdOHflnCF3frwBaQv24FRuCaZazIf0+15lQxKo4Up51cSRjo6zWtTKEbTOWuRe5o9/s92t/Uz9Wf2SWIGzKLcxQPFhcp90lOYNIsQploGlO090alz8WzPdHWXVtr/8oO6gdEqnqXc1eFu48zTeX249fkvNm+XHq2w3SrU8GV09zxw9nXp64DslpVTuqv5mc7ecxK97sjHJA4MaOgtyqzk76p6+CUuSpMnkl0UqzAEkJ0i/XKx+Dz89MEBxQo9SAbV3qfQ7eENA7ihDU3I01GjR/v1Wz45yWfPnrhrq3vm3Vlqs7OrNfLyMuWT+seR6o9TlFu0DTKgqAdh/Rl4VbM0krjt8HiM/21RrPXcuuVUHzzlfyQ5XDuHJS8Vu5RGXVS55c+T/Fu61CtBFeHixpWZ345krj6DL2yvw6Zqjdj7hmHpDItRm6xiuO3wez39/vQ3YwRyF3acFzfQZoNRQM8MQaShxe9xNo7cMi1wt80KRVSMzdw6PSEN8u8PZOSLBhQBFo3UdMZmqukIP/mideRA/JZ6eY7vxdXWm78oDiQiDmMlR/c20ztKulFmPSLpSRgCndw5UerXSarymd6/13nzn94Muba/YTvWno/PLnXz4wS+3YtEuxyNPK9HjnZX4rwDdjRmgKKDFAFY1uRJsyG145So1u1R7wp8+3Yh/b3A+pL5etnh4QKzjF+X9fkqrePSw+dglTYNGkR5I1E6Jp77Z0v3KArZvNh3XbCwhJeTMFO4ukc4vSzWrwU9fLsGLAnQ3ZoCiwFsKun454uopai/C1uPG8pmdos+dJ3Nx4qL2gRwArLXTc+RCoXVmJ1qWcN+/NrvdSl+JP84W4kq58watag6hrQd3qlpIPbW6Flu8YyurUmPeL2+gZ+8nQeMmBig11Zq3xKJITpcMXLL8U4yzqLJSwjQ7RZ93z9rosXQ89O+tstYT8anF00HlpSLnjeZqTufgjKICaQ9830dmb8O8rSddmvSs+rt4us2Zo3NT7ZRUN8i1zEcOeGI4BfEuP92ocX5pcYaK+hNxHJQaZqy83ntg4AdrPdrKXi9Kr5mCK+63RNeLJHk+zBMwPqrlxf/twpVyZTd2Eb/Wywv2WDXKlD3a7rWL4A+ljQsdKK+oxIXCUsRFhqq2TbU9PVe9ARHtsfwFRG0kq7WsS8X459JDPnE/URMDlBosnygUt4TWgtXAVJ65up1l6q/9tFfR5y5rMKmdEtWp+nztMXy29ijeGdHW7rpfrDum+v5FyJT3nXHc3bPmGCiq8+BBWLDTte9ScKUcM+11b3bBqH9txo4TufjPk91c+ryjq1Cto6nFw0bNtIlYgulJX6w7JmtkWNn5uxazeAv6GzFA0YGjc0GSgEOWgZGA5429UUMf/mqbzfc9MSCUHG/9VpVJTP3Ffp23kiGm5fL8tV87B9ui0VgssgmaAVqq2XbJXdXzGn298biq21WTnxY3u5r/i//Ta0punnKxSF6vNC0eVEVtfsYARQFPPASeKyi1e6P3FKOMCCubzIvvxpd/1TYdFuQc4U/XHLXbEJjkkZvvmgAUalR1+ese13qHeOIq9DPItS5K+zstrT4k71rfXKMXoDcfGwYoguk6bYX1GzrkH6IW96kl61KJR/dn63jKOcLZeVc070LuKaKfUSYTcM8nnmvgXc3RcZFQNeDXy/N34/m+SZrs31+LIhQH5GQtpVcrEBzgr31iajCZgI1HL9gfxVgQ9ubS8kbsxaOALvdtFfbp7CHJ60pMaiirqFRluHm5o5fqrbyiUla3Yk/afUrekOdqk3tmbzueizIXev9o7fnvd2L7iVzZPdbksMzHPHHpK83CUv62DGVXrX8LT7W/G/35Fmw4Yl1CsXz/WZvX0+s/7fX4mEa+hgGKAntOX1ZnQ4I9Tnp3eAL87ef9Noc3V2rwR+tc+pynA8C8knK0mbJEs+1XVkrYfvySKvOKaE2wS00RE9SfELImtUtQ9p/Jr9VdXWlAVFxWgdOXrUs5PVGNkX6tzVBNj3+zHa1eW4yfMqwbX3+z6QTu+9dmzdPljLeNBG6JVTwKnM0vRcGVcoSHBHpupzqce/YawVYzWhVQzczO0/Q4XuUyp4V3xX93ZOGv8/dotn1Neeh6UjqmjD1K2oi40oVV7TYogz9ah+AA6+deo7SRcNZr84V5GZ5JiEJGOb6uYAmKQmrOEilJEp781sn03Tqce+sOX7C77PO1x1BogCdn0RksxrMyX+suySrS4zifuVyCW99aLmtdhz36AAT4yw8gev59pex1q2lRuFdao3rG1d+gpKwCxWVVeY03lxKQfQxQdHQguwBL9p11uE5BaTm2Hb+k2hOZu9767YCwXdLIM4x4q7haUYk/z9nhkUEGP1px2OncMscvFFkPJ2CHvSoYy8Di41VHcLm4zO4EdY6c8nCDcbkqJQk3v74Yt7y+BOUVlV5dSqCEJ+aDEwmreBRS84lDznws93yyEWfzS/HOiDZIbhSpaPuSJCHzQpHT4MbL28jqztHR31Vj4rv8K/oOauetlu0/67FZh+WUGKT+czUA4JEeN9pdxwTAX8bF+Y8lh7Dz5GVZaatJxIbBgHWX727TVuK2pGgdU+PY9MUHMeGOmzyyLzXa0hkJAxQdVD8NyAkMqqeUX5hxGq3jlQUoszcex9Sf96NZdJjtdFzLSEUtPu0ss5hcdLae2qvPgTs/3mD1fnGpWL1vbFH6LDt9sfW8TSaT56teilwoXfCErzYct7tMgvw2IssPOC6JrbldI7lQWIofXRwd2BNmrT6KuLqemc7A1rADizLOeGTfemAVj0Jq9Mh4b9kfij/jSob+/rX9HLtQpPzDAjhXIG9kRfIwhefirNXWM1/r0S5EzBDcse+2nLS77Jfdrg3+poej5wv1ToLmTl7UL4+dty1Lt31rjQGKQmpkdN9sOoEr5RWKW92rXRXDqh3ydssPnMXL83d7tCrDE+0lXGlvohdHje69hZEbvYtMcYCydu1aDBs2DPHx8TCZTFi4cKHV8gULFmDAgAGIjo6GyWRCRkZGrW2kpqbCZDJZvUaNGuXqd/AotW7q0347oGjwJS2CCV5U+pn80z6bDd7KBW0TYFQbj17EvG1ZmO2gKoU851w+S0VJPsUBSlFREdq1a4eZM2faXd6jRw+88847Drczbtw4ZGdnm1+fffaZ0qQY2tebTiha31EwMXzmBvsLZWBJiufN25aFB77cUuv9if/dpUNqlDFij4pjFzxTzZB/pRzpLjZYteXQWQFmVFeRq4Mdis54V4QxKG4kO2jQIAwaNMju8gcffBAAcPz4cYfbqVOnDmJjY5XuXnciNih1ZYAm0t+Ji7VLUNQYkl9r7pa86dFIVo2B6x60EVDWNPjDdTiVK2bXXT2tP3wBPQXuiUNi0q0Nyty5cxEdHY3WrVtj0qRJKCiw/6RQWlqK/Px8q5dejFLaIEkS8j0w5gORUkatWpTTlsIIwYkex99WaaE3Meo5LTpduhmPGTMGiYmJiI2Nxd69e5GWloZdu3Zh2bJlNtefNm0apk6d6uFUGpuzYZuB68GWQWIuIYyY5V51mjdgXkyuUDKBpVEeBKudzVd/1nFJklApeX7GaZHoEqCMGzfO/HdycjKSkpLQqVMnpKeno0OHDrXWT0tLw4QJE8z/5+fnIyEhwSNprUmvU0XpTeGqjCJtRv3KOWpfYITJ85R4d+khvZNAXmSBgikSjJY3/bpH/W7f93yyERcKy7BiYm8E+vtmh1shvnWHDh0QGBiIw4cP21weHByMiIgIq5dudAxm3/n9oPOVXGC0pxURHT1fiA5v2C4BNKrdp/L0TgJp4LLGMyTbUz2vDsmTfvIyTl4qxuGz3j+OjD1CjCS7b98+lJeXIy4uTu+kOKVXI1ktG09mXvCt+R20MPCDtZrOIEykFm8e2Iu8i+IApbCwEEeOHDH/n5mZiYyMDERFRaFJkya4dOkSTp48iTNnqobfPXSoqpg4NjYWsbGxOHr0KObOnYvBgwcjOjoa+/fvx8SJE5GSkoIePXqo9LW0442lDX9frE3JjC/xpeBEMlr5OwmhRMHgcm//dkDDlBiTkjY83kJxFc/27duRkpKClJQUAMCECROQkpKC119/HQCwaNEipKSkYMiQIQCAUaNGISUlBZ9++ikAICgoCCtWrMCAAQPQsmVLPP/88+jfvz+WL18Of39/tb6XZrwpPqk5LTqRHJzNmlzxzWb5Yz9tMUB3e0+pfihW0obHWyguQUlNTXX4BPXwww/j4Ycftrs8ISEBa9asUbpbYRglb84r4ay4pI2MGjMwE8nhiyUAaiq76nvHT4hGsqS+z9Yedb4SOLQ6EZERGOXhWE0MUBSS031XBKXl8gKPJ77ZrnFKiIiAAg4cSQoJ0YvHSP6x5BAyPTSvhzu2HpdXh7vq0HmNU0JEomA1i/F4Y8cMuRigKDQ//ZTeSSAicsn4eRl6J4FINlbxEBH5iMX7cvROAilUPfZWoQ9WkTFAscDxHYiISCT5V8rx5fpMvLvsD72T4nGs4iEiIhLUnz7dpHcSdMMSFCIiIhIOAxQLxQqGYiYiIiLtMECxkJ1XoncSiIiICAxQiIiISEAMUKz48Ig4REREAmGAQkRERMJhgEJERETCYYBCREREwmGAQkRERMJhgEJERETCYYBCREREwmGAYsHEXsZERERCYIBCREREwmGAQkRERMJhgEJERETCYYBigU1QiIiIxMAAhYiIiITDAIWIiIiEwwCFiIiIhMMAhYiIiITDAIWIiIiEwwCFiIiIhMMAxYKJY90TEREJgQEKERERCYcBChEREQlHcYCydu1aDBs2DPHx8TCZTFi4cKHV8gULFmDAgAGIjo6GyWRCRkZGrW2UlpbiueeeQ3R0NMLCwjB8+HCcOnXK1e9AREREXkZxgFJUVIR27dph5syZdpf36NED77zzjt1tjB8/Hj/++CPmzZuH9evXo7CwEEOHDkVFRYXS5KiKLVCIiIjEEKD0A4MGDcKgQYPsLn/wwQcBAMePH7e5PC8vD19++SW+/fZb9OvXDwAwZ84cJCQkYPny5RgwYIDSJBEREZGX8XgblB07dqC8vBz9+/c3vxcfH4/k5GRs3LjR08khIiIiASkuQXFXTk4OgoKCUK9ePav3Y2JikJOTY/MzpaWlKC0tNf+fn5+vaRqJiIhIX8L04pEkye44JNOmTUNkZKT5lZCQoEkaOAwKERGRGDweoMTGxqKsrAy5ublW7587dw4xMTE2P5OWloa8vDzzKysrS5O0SZImmyUiIiKFPB6gdOzYEYGBgVi2bJn5vezsbOzduxfdu3e3+Zng4GBERERYvYiIiMh7KW6DUlhYiCNHjpj/z8zMREZGBqKiotCkSRNcunQJJ0+exJkzZwAAhw4dAlBVchIbG4vIyEg89thjmDhxIurXr4+oqChMmjQJbdq0Mffq0QureIiIiMSguARl+/btSElJQUpKCgBgwoQJSElJweuvvw4AWLRoEVJSUjBkyBAAwKhRo5CSkoJPP/3UvI33338fd911F0aOHIkePXqgTp06+Pnnn+Hv76/GdyIiIiKDM0mS8Vpe5OfnIzIyEnl5eapW95y4WITe/1it2vaIiIiM6vBbgxDor25LECX3b2F68RAREZE4vttyUtf9M0CxYOJg90RERACAyYv26bp/BihEREQkHAYoREREJBwGKERERCQcBigWOA4KERGRGBigEBERkXAYoBAREZFwGKAQERGRcBigEBERkXAYoBAREZFwGKAQERGRcBigEBERkXAYoBAREZFwGKAQERGRcBigEBERkXAYoFjgUPdERERiYIBiIcCPh4OIiEgEvCNbiI0M0TsJREREBAYoREREJCAGKERERCQcBihEREQkHAYoREREJBwGKERERCQcBihEREQkHAYoREREJBwGKERERCQcBihEREQkHAYoREREJBwGKERERCQcBihEREQkHAYoREREJBwGKERERCQcxQHK2rVrMWzYMMTHx8NkMmHhwoVWyyVJwpQpUxAfH4/Q0FCkpqZi3759VuukpqbCZDJZvUaNGuXWF1HLiA6N9E4CERGRz1McoBQVFaFdu3aYOXOmzeXTp0/He++9h5kzZ2Lbtm2IjY3FHXfcgYKCAqv1xo0bh+zsbPPrs88+c+0bqGz6PW31TgIRkdfrd3NDvZNAggtQ+oFBgwZh0KBBNpdJkoQPPvgAr776KkaMGAEA+PrrrxETE4PvvvsOTz75pHndOnXqIDY21sVkayfAn7VeRFRbUsMbcPhcod7J8Bov9L0Jyw+c0zsZJDBV78aZmZnIyclB//79ze8FBwejd+/e2Lhxo9W6c+fORXR0NFq3bo1JkybVKmGxVFpaivz8fKsXEZEnJUTV0TsJXiW5UYTeSSDBKS5BcSQnJwcAEBMTY/V+TEwMTpw4Yf5/zJgxSExMRGxsLPbu3Yu0tDTs2rULy5Yts7ndadOmYerUqWomlYhIEZPeCfAyJhOPKDmmaoBSreaJJ0mS1Xvjxo0z/52cnIykpCR06tQJ6enp6NChQ63tpaWlYcKECeb/8/PzkZCQoEHKiYhs4/2UyLNUreKpblNSXZJS7dy5c7VKVSx16NABgYGBOHz4sM3lwcHBiIiIsHoREXlSSpN6eieByKeoGqBUV9tYVtWUlZVhzZo16N69u93P7du3D+Xl5YiLi1MzOUREqrktKVrvJBD5FMVVPIWFhThy5Ij5/8zMTGRkZCAqKgpNmjTB+PHj8fbbbyMpKQlJSUl4++23UadOHYwePRoAcPToUcydOxeDBw9GdHQ09u/fj4kTJyIlJQU9evRQ75t5sYbhwZg8rDWe+S5d76QQ+QwTW6EQeZTiAGX79u3o06eP+f/qtiFjx47F7Nmz8dJLL6GkpARPP/00cnNz0aVLFyxduhTh4eEAgKCgIKxYsQIffvghCgsLkZCQgCFDhmDy5Mnw9/dX6Wt5t7TBrTCkbRye+U7vlBD5DrZBIfIsxQFKamoqJEmyu9xkMmHKlCmYMmWKzeUJCQlYs2aN0t0SEdnV+6YGWPPHeb2TQRqJCAlA/pWreieDPIyjkhkQi5rJG9zfWb2eeI3rhaq2LXv8WISiG3ZJ9k0MUDTwwxNd9U4CkfCeTm2hdxIU4T1SP45K7eXq1qy+CikhT2KAooEuGl8IzCiJrLl/+3LOnesuiFNo6C4hSvtSNlIXrxpyqntzPnmQb5j7eBe7y1i1amz8/YyHAQo5tGdKf0y/lzM8e4Pn+ybpnQSzwW3UnShUhRoAAEBwgP0s0Y/3NyKPYoBCdnVOjEJ4SKDeyVBNl8QovZOgq5BAcS739+9rr1pQ4Y7IUPnnN6tWbWtaX4xJFCfecZPD5fz9jEecHItqaR2v3ZD+HZs6H7bb0dOkETGDEkdwgBhjHvVsoWR0WJ5Atiz9Sy+9kwAACA0S45zyJh2a1NV1/951B/IyFZW2HzHV6HInZwvs2mdMLRreoHcSdKB9cQwvB9vYANh76V3IyTNLYAdzCmy+r0Y+KSezZX5sTL2SGuidBIeqn8qEuOHXSIOjNImQXF8l50bprMowwF/8X3DysFv0ToJQGKAYUFcVujHfWD/M6TqiNAr86pFbVdmOr7Tif7J3M72T4NDtrRrqnQSXeHqgtkn9HbepIGWMMO7OIz0SMbC1ug3I3TFlWGtd988AxUDGdmuKXa/3R4PwYLe39XzfJDzUrSnmORhUrrqKR+/GjKGBrFuWq1l0GGIiQvDj0/ZnD/dGrp6jRaXyh0/XMj7pZKNN2LO3i9PryhFRqoKdjSYcFRbkoZS4J21wK72TYNYuoa6u+2eAYiCB/n6IrKNOr5qw4AD87c5kh6UxYmQ74hKxEXH1vSKlST28dXeyvonxoEZ1XRuEa/Uh+fP31L/B/QcDpXrfZF1d17FpPaTo3HDRkS/HdtJt3wOTY/Hc7dqUkrx5l+euJb0fCEUiXg5Lurk7pRE+faCj+f/qJyNBHpDc5i3fw1uomRGP66VOtZatNN0UcwPWvthH05I8e+fmZw92tPq/c2IUvhyrTpWnFvTsnWUymfCsRgGKKF2pfQ0DFBX8qWNjvZOgCpOp6inE8n9A/4he7/3bI2iy7PLmNjghgf6aXYfxdUPRpH4dt45evTDXSj5DWL0pBE9eO0bLV7TEAEUFQTKK+tUYOXOQyqNvOiNKI1m1+EIJiqP2ABIk1FOpilBEamTstg6fycEyub56uLPj/XpJ8Kj3NabVcZQ8EDa0ig3XfB9GwwDFhkCF3dHqy2h8NWtMR6frOLJiYm90bOrZkVCrL3a9Mx1PZA6Wjr8zxKP786RVk1L1ToKZ3ueVLVqV1t3iZNDFsGB5JSUCHrJaYiNC9E6Coakxc7O3YIBiQ2yksgvsyd7NNUrJdc0beH7wLVGqeEg9desYoyeDK9o2jtRku9fbYmkXHigZBVX0G9gCBT3IqkfLVqv0wNFP5E51mZxhGUh9DFBsUHL9t4oNR1hwgHaJ8aCaxaMiPuE60u/mGPzfkJvtLr8pxnkm+N04+7PZGoHlT1azB4i3G925CaaoMNBV31YN0c4i2DHYZaC7eAU9qhY92xN7pw5wPnSCzjFZQpT2jWQ7+/hcYbYwQHGTKGMAqKFmVYoRv9ujPRJxo50W93KeoLo3VzI3i3gsf7LG9epg66t9ry/z8lttgL8fHu6R6PZ2vnz4Vix8pocKKVIm+gZ5pVtKJjgUnb+fCTeo+ICn5AxPbaksgL8pRrtS7FcGt8JfB4oz/okoGKC4yV5xqx71sP4qt2oV4XbWLDrM6ukpzElRuJ+fCaO7NNE4VcbRMFzd8zA00B8dmtRVZbRakWsqLINzT8Xprw2VV/oTYGfum/AQ7yjJtaXfLTGqb3N05yYY0iZO9e264olezc0l8QJfFh7HAEUjfXQYzlvuE5g91U/YPVpUDd72YNembqfJXdPvbWv1//KJvZHU8Ab0u7n28W0YUVVM/HD3RLw0sCV+ea4nfniiK5rWr4NvHnXci8IlPpiTBAX4YcHTPZA2yH5VGikXFRaECDdLRppF699OwtVYzlmwOmlAS3n7VxBNilpCHB/p2qCD3ogBiptEOsnVKsL/5tEu2PJKX3RRYc4fdzUMD7GKA+IiQ7FsQm+M6WIdPA1uE4uXrmViQQF+eDq1BZIbRaJLs/pY82If9LqpgdBP7GpxdA6o0RtKzQaaAl06ZpZtDRKiqm4UQ9rafspe8HR3tx8Kqk3qL+8GbMurg2/G9Hvaomvz69frXe3j0dHG8Pla0+oSCwnwQ38NSlH00LdVQ3z6QAe7y5U0mPZ2DFC8iL0MX844LQAQHFi1nr+fCTGCdxWsebOdNaaj4XqoaFEkL+JNX2Q1G2danve/PHcb/vNkN9zVvpHNz3Zo4noA8EiPG63+l3vu2vp9U5rUxchbEzC2243w9zPhng6N8cGoFHzxkH7DzutFyemv16XSpH4dDEwWo2pJdAxQ3GTU+4FlukOuBSZGekLxhdIQtckpYWt47YZtr8uuSCWGahjbzX41ZmRoIDonRunynZXssvpSiK8bioNvDMS7I9s53cb/nurmeuIckJvsYe3irf6XM8O13qeekhLqCXfoMxO1szZ6RsMAxYvYu3ycXVZrX+yD78Z1QWpL9dvNuNtg1V6mJEKAUh3YuUqL/Nbdm+l/n+qGx3om4l8PKn/6fsPOhGr3dGiMAa3FDH6TZHQ9d8y1463V+RtopwFtTZ1ujFJlVnRXe+DMuD/F6v+x3W/Epw90wJLxvdxKj+LTX1EgWPtHa9HQds+eAAeDfWrZm87bHiAYoNgQp3CgNmd+ea6nqtuzx97J6ayxa8OIENnda6eNaKMoTUEyM0x77GXkruTvao9I+7WbDW8FiLFqaVo/DK8NvUXxYIWA/fMsKixQ9d5EanG11PCxnu53Z1aLvWvEE93Kt7zS1/lKMvj7mTAwOQ4t3RywzdM3aNF6TnlXeMIAxab3RrZHn5YNVBu0K7mRNiNcytW6UQRWCzTEuVKeKi3pd7Oym1WKgjYIQ9vGYfo9bVHHy4pg5TCZTB6frkAuV29oSs8VpVS50TjZiBrXlScHqRSh1FQuI6VVZAxQbEiIqoOvHulsuEG7htrpbWCCSbXBnSJCxBgkyrI3idJSHVse65mIj+5vb3PZRzWKo10xc3QHjLw1Af4aP+GJ+ASl5dDsYUH++NwHG4NasndKqTkAmjdSeinaKpFS43p7vm+SClupIvdKm9RfnzYySjFAcZNIVX6P2il2liC5PcZCtT6tGuD+zgmYrMKQ4nLZevq2fOf+zu4PzJbasgHqBNnO0HslqReo+iuciFLpxJUinY+WtIpR9kwZgDsM1LhbLiUlO/aOrb+fCfumDlApReL5cFR7j+5PSSmgkuvQ3bZsrhC9l2Y1BihexNFIsmqNMmuCCdNGtNV9tFa9i1BHdLDd9dSZAIvfQc4v8s8/tXO6TrDMbuRqjNnhSgCkZbsAP5VHTxad0nYljqpgav4s424Tp12NHHfa6f6thLtnj71zW0n+JGdduSUecr+PrXQ/f3sLmZ/2HAYoOmqkYFIttckdG8UeJVUV7jYorGdzjAj9IpT9fxuAd2UEDrb4WRy3mt/AXo8AZxY9e70Rtq2fZdaYDniwa1PcnVI7Q+/mgcH4JElCzxbXS6HivGCkzOqnaVdjL1sTV1peU0o266wKraXMnkojOyUo2KvxaV3a+MmYDjaH0ndlvw/IHdXbzrbv7djY6n9b58wENwYL1AoDFDe5c5K7WwrxyuBWVhNY2e9mbL2kdXyEy43bqr9vgL8fPn2gIz4c1d5p6Uwbi0bC34/rqnifyY0i8eKAlvjgvvbm91wqQXHwGSXbqxMUUOsJ5DkHTx+Wx6c6WJMz7oNclj0fbD1hD24ThzfuSjbP4VI/7HrAZ2+UVHtCXZyyfmByLB7ufiNGdGiEu2wESt5shI3va6tzW/fm2gSL9npk6V0KaTRK26AMahOHj8fYHzFWCbUHoWwVG6Hq9rTCAMXAnujVHHMfV37Dt/Uk7YqBybG4s30jqwDEFsu6224uZsLP9GlhdWPzZN4qp1i9aX3786Asevb6zLjjbmuGhc/0wCcPdKi1VS0bk1pyZabezx7siGbRYU4bpNrqUm8ymWAymTBleGu8N9J5QOtt3rMIrB2xNwmgp6hSoqDRTxsVFoTOiVGqbrNjU2XbE7UnmiW5pWWhQca49StO5dq1azFs2DDEx8fDZDJh4cKFVsslScKUKVMQHx+P0NBQpKamYt++fVbrlJaW4rnnnkN0dDTCwsIwfPhwnDp1yq0vIoqBrWP1ToLHffJABzzQ1X5pkDv3XSNkCs5YjgHi52dC+4S6CA7QpruxnJuM5Xwzcg1oHYuVk1KddplPbhSJxvWMX4XjjJxGhpP634Sf7ASDnRM9N8+VvStIkyoODS7X2Y/cCpPJhLHdb1Rtm+NuS1StZ6NIXhzQEg85GB0Z0KdRrqsUp7SoqAjt2rXDzJkzbS6fPn063nvvPcycORPbtm1DbGws7rjjDhQUFJjXGT9+PH788UfMmzcP69evR2FhIYYOHYqKigrXv4kgGjnInNMGtVJ9fyLcwOMiQ/HmXfa7+soNULooeEJyJeix/Iizm+izfZQ1GLNX+nFvx8ZujdgZ4GeczESpunWc3yCWT3BvZFG1zXmsC94b2Q7NGzhvL/Ts7Ulol1DX5rLE6DCHbbNcGere7nKD1eXU7LlWXZ0qd5RcOap77LlaZVnNpUbjNf539ffxM9XuSh4WHIBXBnvPTOOKf/FBgwbhzTffxIgRI2otkyQJH3zwAV599VWMGDECycnJ+Prrr1FcXIzvvvsOAJCXl4cvv/wS7777Lvr164eUlBTMmTMHe/bswfLly93/RgJ7sndzhGs4NoFWvSUsR0t0pXhe7uX35cO31nrPXvWKu4HZupf6OBxWW61aCEe9cOR8g363NEQ7O/Pi2GKkypPwkEC7JQzVQty8gTjyRK9mij/TMykaIzo0dr6iDWNqtDlLjLZfLagPFc4eFTaR8Xp/bFVphFpnXhzQEq1kjl5ruw2K8i/co4V6wxa8N7J2/uLn5D5ggskw7Y9UfTzLzMxETk4O+vfvb34vODgYvXv3xsaNGwEAO3bsQHl5udU68fHxSE5ONq9TU2lpKfLz861eerrHIoNSeoLG1VW5/7kHTrS6dYLwrwc74qtHbpX1FFNzoqxKDa4GdzdpMpmsGpc63JyMn1irnyE4wB8/PatgqgSFQareg3nZK2GopmVG6uknzS41ekw5+qns5SuuPIPYu2Zbx7veUDL6hmDc2T7e+YouCAsOQEOVx+l4srftYLRhRAgWj++l+kSptkpFnr+9BeY81gWpLRuoso/qtl01BQX44f+G3IyJdiYsFHWsJFtUDVBycnIAADEx1j92TEyMeVlOTg6CgoJQr149u+vUNG3aNERGRppfCQn6dYe7JS4Cfx3oencso0SuNfVvHYs+MicTtMwQI0IDzV1Z1az7NOhhtFIzn3ClBKzmkOtytzB1eGsMbxePgcnat5lytfu0MyJktO4kwZW8wJXSm7p22lr84952DtuOObLllb74cJT7IyzL4agKZPWkVDRvEOa0JK5mlZwrA1faKrUNVTB1RcOIEPRMiq51nbvVRs/Ohx+/rRmesxih1qj3HU0quGv/AJLTzNfROmlpacjLyzO/srKyVEurEo3qhuK3F26ziu5FyCRFYzIB/7i3LVJbNsATvZohIaoONrx8O7b/3x0Aqp6+bH7Oxnv2qnJcqbdV8hlPXM/uzKQ7oHUMVk9KxacPuNaNcWz3G/HR/Slu9ahZMbE3xvdzPEz3q4NvtipxVNNSN2e+FZqNn+XNu5JtVgu5evNpEB6MycNau/RZUXpi3RgdhhUTU52WxFV76+5k9Lu5Ya3qtkkDlD10vvundmgVG4437czgrYSnY4dmDUSrWrRP1QAlNrbqaaxmSci5c+fMpSqxsbEoKytDbm6u3XVqCg4ORkREhNXLUywnd5MTjDhbRcuT0e44KDrkJX/qlIDZj3Q2VyE0qhtq/nv2I7XbmiilxlTxABATUbWd9jIzuLmP25lA0oUfdsb9KVYj0ioNum6MDtO1a2rzBjc4HdxrXK9mqt/Mvnu8C/ZNHeBWgCc6W0fMEzPn9rpJneoHUY3p0hRfjL21Vtumm2LC8aSTNkmW1W73dGyMxeN72e0R17ZxXbfTKoeSHOPnZ3tiWLt4fDKmo2bpUZuquVtiYiJiY2OxbNky83tlZWVYs2YNunfvDgDo2LEjAgMDrdbJzs7G3r17zeuIxN5cFn/pV1W/5+oTiFrknKCiFe8pmd3ZXl18t2b18ULfJMxycyCkdS/djr1TBzjscmgZ4CXUU95F1574uqF4b2R7q/d+f+E2q4HU7KZJxTYK7oivG4r/PtUNS//iudKM7i2iPTqLrh6a1pd/nqnZk++Lhzrh9xduU217alCj8b+cLQQonPfKln892BGfPtBBVqBX/ZDl6szYSh9o2jSOxIz7U5AQVccwVeSKr/LCwkIcOXLE/H9mZiYyMjIQFRWFJk2aYPz48Xj77beRlJSEpKQkvP3226hTpw5Gjx4NAIiMjMRjjz2GiRMnon79+oiKisKkSZPQpk0b9OvXT71vphJ7F8cL/ZLwdJ/mtRqg3RznuHRHjIJR4zOZTPiLnUZgSgQF+Nkc9l9pUKfWTeLmuAh0bV4fv+7Odunzepxft95o3T1czSBJSR2/I9E3BCPAz4Sc/CuqbM8Vcg7LT8/0wCerjyJtsPpDEjhjMlVdD87yMIfb8OIcTs413qzBDbLbXK15MRXnC0odDvKoFhGGo3CF4gBl+/bt6NOnj/n/CRMmAADGjh2L2bNn46WXXkJJSQmefvpp5ObmokuXLli6dCnCw68Xx77//vsICAjAyJEjUVJSgr59+2L27Nnw99euS6Eaat60LIOTX5/vifQTuU5HadX6NHmydzN8tuYYfny6O+6eVdUrytNP1WruTs0LS8uSpPAQMQZ9Uvq0qclYXSoc54SoUDzZq7nd9kpKNWsQhh+e6IrEtN9U2Z4rah4WWz9Vu4S6+PRBZUXwreOclEg6+JEt5/+JDnP/WIt6I/RUqhzVZtb8vesEBaBpfX1KAo0SRio+OqmpqQ6LlkwmE6ZMmYIpU6bYXSckJAQzZszAjBkzlO5eWK3jI9E6Xn7VhRZMJiBt0M1IG+Q9A/WIwNHF/Nbdydh58jIG+OAIwlq6t0OC/AnSZDChKm/64L72GP9DBj6xqBoclByL3/fa7kHodLsC5PSRzga8q5FdN7NoaOvnZ8Li8beh7Gql8+2o5CU3ekHqSU7pkKOxbQJVHnTRZHJ9PBMxw8javLsiV0Bqj+ooWvsSQIxM21PGdGmKMV3Uu5ECkJV72DvGPnToFakemv6ulEYY2jZO93lvAHWuXVfaH9ec00bNiePk3MRvaqiscbMeI+G++6d2mLxoHwpLr8r+zHsj29kswXzu9hbYcOQChms0bow30/8qNRDRb7xGrP91Vl+r13eyV1QtqyeXCejYtJ7bDXgdsZdnyxmC3WhecaM9xlcP34p+N8fgtaG3mN+rGZyofV33vqlBrW6surL4fuHBAXhG4TQOahM9HwWqeunULMFztfpqYv+WWPB0D9VHRTbaFAauYAmKwVn2PlGrQaEnRak8jbgjvW5qgC/WZ7rV7VVOnhAfGYr5f1beI8389OVC8ub/uTsW7jyteDwHI3iiV3O0io3AQ//eqvizfVo1RJ9W8gYYdEV4SCDO5pdavff1o51trqvmfblenUDkFpcjpUk95ytbyJjcX5gxTLTyaI9E/HtDptV7rnxjUdvTOKLVoIh6YYCiIzWeJEKD/LFiYm/4mUw2e6MYnZqZRK+bGuCHJ7qimYxSBstARKs5jmrv89pOXfjKHZvWQ8emym5WNdURKMA1ys1h1pgOeGFeBl7o2wI/78pGSpO6mu2rrkUwv+DpHvh20wm7Q7hb+ku/m7Bs31k80K2pIYMTpdffIz1urBWgqJIOIUuor18nB/42UNO5q/TAAMULiFasL+aFXKXmfChKiVI8rUU65j3RVf2NqkTUcOWmmHDzuCEDk+MUfVbubzj93rY4kJ2PXknXJ5lLjA7D68NucfCp6xKi6vhEyUk1Ua5Rtc17oivOXC7BhP/ssrncYQl6jQtITknwiwKUxnrfI7eG1Kjys9yED1QhGpYv/jTBAd719CUaV8+pkZ0SMHlYa7dK8pQGJ8/d3gL1PNSrR22qBWIKfzCtA6OuzerXmotJy3uI3m2VAAYopAElF6plUb63Pvko5cnqDU9VX1FtIj+gTOzfEjuuzZ2lNq1PudiIEAxvF4+Rna7fzNU41JOvlVg9d7v+N27AvWvX2UdFyRZYxUNkg8nO3+QZNW/eWvZYELlKUk9+KlYJTR52C6b+vN+lzyr97U0mEz66v2qm5f9sP+XSPoHaQU33FtE4+IY47TwkSUIdF6d8cHZIRQmeWYKigNpRpeX2Pro/BasmpeLHp8Wbj0hLljcHUS4KQN8qHlFumFOHV80z9bwgT4xGJ8av6nkfjmqPR3ok6p0Mp0ZfG89oSBv7bYlECU6q3dYiGnenNELaIMdd8ZXmZ3LmA/MElqAIYni7qkF8zuWLdQFozldzbQdE6cEytvuNGNQmFg3DQ/ROCquiDEKLn8lTkwU2qhuKg28MRPC13pB6nnGN64XiVG6J0/X8/Ex4/772qu33+3Fd8ffFB/HmXcmqbdMdLEEhssGV0pybYqp6Uw1u49qw93rfhG3tXq/gRIwQTVuMucQTEuhvvg71PAeX/aW3Lvvt1rw+Fj7TQ9GM81piCYqn+ULO6yJRM2y56Zr3RDdsOHIB/Vt7Zvp0X8JjI7bbkqJx/GKRzXFgagbeolRhisyIg25qgQGKjmIj9C8611t8ZAga1wtFaKC/UN1cLatZ5N4bo8KCMKydZ+bbEDWYE8G2V/vhp4zTGNJW2bgkWnM0kZzRffNoZ1RKtrv4GjW4VNw410OBV5OoOjh5qRg9kxrI/kzNlIlSjewMAxQZbkuKxrrDF/CgirOrAsDgNnHIyLqM3jfJP9G8jZ+fCWte7MNnKguT+rfEqoPn8WjPGzXflwjHvXmDMBw9X4T4yBCcybvi9vYahAfj8ducj7Bq5qGDkNwoEjNHp6BxvToAxGoU7swDXZtg8d4cXCgss7ncZDLBX4STyQfMe6IrFqSfMjfqdeTFAS3xzabjmCjAoGuuYIAiw+cPdcL+7Hy0b1xX1e2GBPrjb3fab4zUo0V9bDhyUdV9eoLSthTVT12VlQbKsTXUrMEN2Dt1gM+M/Dn7kc74Yt0xPNozEb3/sbrqzRp3b285M4a2NeaMtnpWyxjlMvBUqUR83VA8e3uSrHWf6dMCT6c21719m6vYSFaGkEB/dGhST9VxAeRoWt97i4TJMWfBSffm0Q6XG0lCVB1MvTOZ57uvkJmN3t85ASlN6qKbm9NTuMpIJVyOGDU4AViC4nHOzvlgi372Qf7GjB/VvBx0K0UQNHMa2jYOA5NjMUjhvC9Gp+VZ8GiPG/Hr7mzcruGsx7YY6b7RMjYcv+/N9ug+p41o69H9uYuNf9XHAEUwkaGB+Ps9bWAymXDyYrEq2/R0BF0vTP4cHsEOZmC+v3OCub6eqsRFhmhWTSBS9urJ+LBj0yjs+L9+qFdHjMGpRPLzsz2x+dhFjLo1AR8s/8OlbYTWGNyshRuTm3718K0Y/0MG/nGvsYIX0cRFhuqdBFkYoHiYnJbh993aBADwzyWHtE6Oqt4b2Q5bMy9hmIwb6OtDb8EP27LwlztuMr9nGUf1adlA1ycoPQpQGoQH67BXY9D696h/A4+9LW0aR6JNY/fGxOh7c1W3+1WTUnGpqAwJUa4/dPRp1RAZr99h6GoLEUSGBmL5hF7o995avZPiEAMUH6RV3eqIDo1rzbZpz6M9E/FoT/GHv/aEbx7tjM/XHcPbd7fROylEqurZItpcTZsYHaZKV2tRgxNBk2VXi4bhCA8JQMGVq3onxS4GKB6m5OIy2gmvpgdU7tItsl43NUAvH+5qbou3NFB0xGjf0WjpdZfSr+trx8cTjNkK08BEGLTICIGPO8XAarirfSMAwC1xEaK2l1WFEc4FAMI2WiYyMtEvf5agaOTBrk3x/daT6NOST8ZGdEt8BLa+0hf1woJQerVS7+SYiVq8TcoZ7ad0Jb1G+44kFgYoGrklPgK7Xu+P8BDrQ8wHQeNoeG0qgjKBAhQi8gylhd1qB2OP90zEF+sz8ZiGbfVEvx8xQNFQZB353W3JmgA1YWYCJcVn1BqVk0/i5GNeGXwz7unYGC1jwvVOim4YoHiYHvlsrYmiBL3jGqH6wgBJdJnQ303Qc5ZIK35+JtwcF6HpPkS+5AE2kvU4JfmsWiePEfN2UW+WogZ33qLZtW6og9t4/0i5RjuXYiM5+zp5FktQiEgYv4+/DRcKy9CorjFGuvQlM+/vgNcX7cPTqc31Tgr5CAYoHibCU5OopRNG4M3HLiJU/zZTwQH+PhOcGO1cujE6DN882lnvZLjElepjT81OTPYxQPEBBssHyYNMJhNmP3Irissq0DBc3CJ83izIHSKMP0XKMUDxQbxWlbGc7CzaC+dsSW2p7iy+NwQzWyHj84V8UvSvyJyEhCRSMb+/nwm7Xu+PSklCSI2ZWam2Gfen4Lnvd2J8vyS9k0I669Cknt5JIANjgCIyo1VSq2DX5P4or6hEmGBP4RzTRr6kmHAsHt9L72SQjlZM7I01h85jdJcmeieFHBD9DqNJN+OCggKMHz8eTZs2RWhoKLp3745t27aZlz/88MMwmUxWr65du2qRFDKYyNBAr6xGUYvoGQoRADRvcAMe7ZnoUyWOreO1HbPEF2nymPr4449j7969+PbbbxEfH485c+agX79+2L9/Pxo1qpqEbeDAgfjqq6/MnwkKCtIiKcJR0tgvJoI3aiIid2k5COS6l/rgXEEpWjQ03oivPtcGpaSkBPPnz8dPP/2EXr2qinmnTJmChQsX4pNPPsGbb74JAAgODkZsbKzau/cqIzsl4EB2Pnq2UHfCQR+sOSIi0kRCVB3dZ1/3VqoHKFevXkVFRQVCQqy7LIaGhmL9+vXm/1evXo2GDRuibt266N27N9566y00bGi7N0FpaSlKS0vN/+fn56udbCEF+vvhzbvauL2dmgGJL7ROJ+/Cc5ZIfaI/q6reBiU8PBzdunXDG2+8gTNnzqCiogJz5szBli1bkJ2dDQAYNGgQ5s6di5UrV+Ldd9/Ftm3bcPvtt1sFIZamTZuGyMhI8yshIUHtZBMRkZfiOCjGpEkj2W+//RaSJKFRo0YIDg7GRx99hNGjR8Pfv6rB1H333YchQ4YgOTkZw4YNw++//44//vgDv/76q83tpaWlIS8vz/zKysrSItlei9cmEZEyDGr0p0kj2ebNm2PNmjUoKipCfn4+4uLicN999yExMdHm+nFxcWjatCkOHz5sc3lwcDCCg72jwSjPeSICgHp1fKNjAJGrNB1sIiwsDGFhYcjNzcWSJUswffp0m+tdvHgRWVlZiIvz/hlM9cBGsUTiGdA6FmO6NEEKBzPTnJa9eEg7mgQoS5YsgSRJaNmyJY4cOYIXX3wRLVu2xCOPPILCwkJMmTIF99xzD+Li4nD8+HG88soriI6Oxt13361FcoSS1PAGnMot0TUNvFaJ9OfvZ8Jbd7vfCJ7IW2kSoOTl5SEtLQ2nTp1CVFQU7rnnHrz11lsIDAzE1atXsWfPHnzzzTe4fPky4uLi0KdPH/zwww8IDzdeP3Kl/n5vW/xzySE80LWpbmlgNRMREYlOkwBl5MiRGDlypM1loaGhWLJkiRa7NYSG4SGYfm87vZNBZCgMqsnTeMrpT5NePESkEVbPESnGHjnGxACFiIiIhMMAhYiE50uTzpEYfKHQ5aFuNwIAet2k7nQqahFrTnvyiMjQQL2TQKRI9+b1MbRtHFrGeH9DelIfuxnbNr5fEnq0iEb7hLp6J8UmBig+wHSt4cL0e9ti45ELuLtDI51TRC7zgac6W/z8TJg5uoPeySAfEh7i/bfHAH8/dGteX+9k2OX9vwCZjeyUgJGdOI8REZEzT6U2R0bWZQxvF693UnwWAxQfEBTApkZegyXVRB4RERKI78Z11TsZPo0Bihd7vm8S9py6jD4txWwARUREZA8DFC824Y6b9E4CEZHukhreoHcSyAUMUIiIyCv9/sJtOHmpGO0E7aVCjjFAISIir3RzXARujovQOxnkIraeJDIQE1vJEpGPYIBCREREwmGAQkRERMJhgEJkIJKvDiVLRD6HAQoREREJhwEKERERCYcBCpGBsBcPEfkKBihEREQkHAYoREREJBwGKERERCQcBihEREQkHAYoREREJBwGKEQGYmInHiLyEQxQiAxE4kCyROQjGKAQERGRcBigEBERkXAYoBAREZFwGKAQGQgbyRKRr2CAQkRERMJhgEJERETCYYBCREREwtEkQCkoKMD48ePRtGlThIaGonv37ti2bZt5uSRJmDJlCuLj4xEaGorU1FTs27dPi6QQERGRAWkSoDz++ONYtmwZvv32W+zZswf9+/dHv379cPr0aQDA9OnT8d5772HmzJnYtm0bYmNjcccdd6CgoECL5BAREZHBqB6glJSUYP78+Zg+fTp69eqFFi1aYMqUKUhMTMQnn3wCSZLwwQcf4NVXX8WIESOQnJyMr7/+GsXFxfjuu+/UTg6RV2EnHiLyFaoHKFevXkVFRQVCQkKs3g8NDcX69euRmZmJnJwc9O/f37wsODgYvXv3xsaNG21us7S0FPn5+VYvIl/Eke6JyFeoHqCEh4ejW7dueOONN3DmzBlUVFRgzpw52LJlC7Kzs5GTkwMAiImJsfpcTEyMeVlN06ZNQ2RkpPmVkJCgdrKJiIhIIJq0Qfn2228hSRIaNWqE4OBgfPTRRxg9ejT8/f3N65hqjDglSVKt96qlpaUhLy/P/MrKytIi2URERCQITQKU5s2bY82aNSgsLERWVha2bt2K8vJyJCYmIjY2FgBqlZacO3euVqlKteDgYERERFi9iIiIyHtpOg5KWFgY4uLikJubiyVLluDOO+80BynLli0zr1dWVoY1a9age/fuWiaHyPDYSJaIfEWAFhtdsmQJJElCy5YtceTIEbz44oto2bIlHnnkEZhMJowfPx5vv/02kpKSkJSUhLfffht16tTB6NGjtUgOERERGYwmAUpeXh7S0tJw6tQpREVF4Z577sFbb72FwMBAAMBLL72EkpISPP3008jNzUWXLl2wdOlShIeHa5EcIiIiMhiTJEmG67mYn5+PyMhI5OXlsT0K+YQbX/4VAPB0anO8NLCVzqkhInKNkvs35+IhIiIi4TBAISIiIuEwQCEiIiLhMEAhIiIi4TBAISIiIuEwQCEiIiLhMEAhIiIi4TBAITIQO/NpEhF5HQYoREREJBwGKERERCQcBihEREQkHAYoREREJBwGKEQGcktcpN5JICLyiAC9E0BEzi0efxt2Z+VhcJtYvZNCROQRDFCIDKBVbARaxTqempyIyJuwioeIiIiEwwCFiIiIhMMAhYiIiITDAIWIiIiEwwCFiIiIhMMAhYiIiITDAIWIiIiEwwCFiIiIhMMAhYiIiITDAIWIiIiEwwCFiIiIhMMAhYiIiITDAIWIiIiEY8jZjCVJAgDk5+frnBIiIiKSq/q+XX0fd8SQAUpBQQEAICEhQeeUEBERkVIFBQWIjIx0uI5JkhPGCKayshJnzpxBeHg4TCaTqtvOz89HQkICsrKyEBERoeq26ToeZ8/gcfYMHmfP4bH2DK2OsyRJKCgoQHx8PPz8HLcyMWQJip+fHxo3bqzpPiIiInjyewCPs2fwOHsGj7Pn8Fh7hhbH2VnJSTU2kiUiIiLhMEAhIiIi4TBAqSE4OBiTJ09GcHCw3knxajzOnsHj7Bk8zp7DY+0ZIhxnQzaSJSIiIu/GEhQiIiISDgMUIiIiEg4DFCIiIhIOAxQiIiISDgMUC7NmzUJiYiJCQkLQsWNHrFu3Tu8kCW3t2rUYNmwY4uPjYTKZsHDhQqvlkiRhypQpiI+PR2hoKFJTU7Fv3z6rdUpLS/Hcc88hOjoaYWFhGD58OE6dOmW1Tm5uLh588EFERkYiMjISDz74IC5fvqzxtxPHtGnTcOuttyI8PBwNGzbEXXfdhUOHDlmtw2Ptvk8++QRt27Y1D0zVrVs3/P777+blPMbamDZtGkwmE8aPH29+j8fafVOmTIHJZLJ6xcbGmpcb4hhLJEmSJM2bN08KDAyUPv/8c2n//v3SCy+8IIWFhUknTpzQO2nC+u2336RXX31Vmj9/vgRA+vHHH62Wv/POO1J4eLg0f/58ac+ePdJ9990nxcXFSfn5+eZ1nnrqKalRo0bSsmXLpPT0dKlPnz5Su3btpKtXr5rXGThwoJScnCxt3LhR2rhxo5ScnCwNHTrUU19TdwMGDJC++uorae/evVJGRoY0ZMgQqUmTJlJhYaF5HR5r9y1atEj69ddfpUOHDkmHDh2SXnnlFSkwMFDau3evJEk8xlrYunWrdOONN0pt27aVXnjhBfP7PNbumzx5stS6dWspOzvb/Dp37px5uRGOMQOUazp37iw99dRTVu+1atVKevnll3VKkbHUDFAqKyul2NhY6Z133jG/d+XKFSkyMlL69NNPJUmSpMuXL0uBgYHSvHnzzOucPn1a8vPzkxYvXixJkiTt379fAiBt3rzZvM6mTZskANLBgwc1/lZiOnfunARAWrNmjSRJPNZaqlevnvTFF1/wGGugoKBASkpKkpYtWyb17t3bHKDwWKtj8uTJUrt27WwuM8oxZhUPgLKyMuzYsQP9+/e3er9///7YuHGjTqkytszMTOTk5Fgd0+DgYPTu3dt8THfs2IHy8nKrdeLj45GcnGxeZ9OmTYiMjESXLl3M63Tt2hWRkZE++9vk5eUBAKKiogDwWGuhoqIC8+bNQ1FREbp168ZjrIFnnnkGQ4YMQb9+/aze57FWz+HDhxEfH4/ExESMGjUKx44dA2CcY2zIyQLVduHCBVRUVCAmJsbq/ZiYGOTk5OiUKmOrPm62jumJEyfM6wQFBaFevXq11qn+fE5ODho2bFhr+w0bNvTJ30aSJEyYMAE9e/ZEcnIyAB5rNe3ZswfdunXDlStXcMMNN+DHH3/ELbfcYs5seYzVMW/ePKSnp2Pbtm21lvF8VkeXLl3wzTff4KabbsLZs2fx5ptvonv37ti3b59hjjEDFAsmk8nqf0mSar1HyrhyTGuuY2t9X/1tnn32WezevRvr16+vtYzH2n0tW7ZERkYGLl++jPnz52Ps2LFYs2aNeTmPsfuysrLwwgsvYOnSpQgJCbG7Ho+1ewYNGmT+u02bNujWrRuaN2+Or7/+Gl27dgUg/jFmFQ+A6Oho+Pv714r4zp07VyvCJHmqW4s7OqaxsbEoKytDbm6uw3XOnj1ba/vnz5/3ud/mueeew6JFi7Bq1So0btzY/D6PtXqCgoLQokULdOrUCdOmTUO7du3w4Ycf8hiraMeOHTh37hw6duyIgIAABAQEYM2aNfjoo48QEBBgPg481uoKCwtDmzZtcPjwYcOczwxQUJUpdezYEcuWLbN6f9myZejevbtOqTK2xMRExMbGWh3TsrIyrFmzxnxMO3bsiMDAQKt1srOzsXfvXvM63bp1Q15eHrZu3WpeZ8uWLcjLy/OZ30aSJDz77LNYsGABVq5cicTERKvlPNbakSQJpaWlPMYq6tu3L/bs2YOMjAzzq1OnThgzZgwyMjLQrFkzHmsNlJaW4sCBA4iLizPO+ex2M1svUd3N+Msvv5T2798vjR8/XgoLC5OOHz+ud9KEVVBQIO3cuVPauXOnBEB67733pJ07d5q7Zr/zzjtSZGSktGDBAmnPnj3S/fffb7MbW+PGjaXly5dL6enp0u23326zG1vbtm2lTZs2SZs2bZLatGnjM10FJUmS/vznP0uRkZHS6tWrrboMFhcXm9fhsXZfWlqatHbtWikzM1PavXu39Morr0h+fn7S0qVLJUniMdaSZS8eSeKxVsPEiROl1atXS8eOHZM2b94sDR06VAoPDzff04xwjBmgWPj444+lpk2bSkFBQVKHDh3M3TjJtlWrVkkAar3Gjh0rSVJVV7bJkydLsbGxUnBwsNSrVy9pz549VtsoKSmRnn32WSkqKkoKDQ2Vhg4dKp08edJqnYsXL0pjxoyRwsPDpfDwcGnMmDFSbm6uh76l/mwdYwDSV199ZV6Hx9p9jz76qPn6b9CggdS3b19zcCJJPMZaqhmg8Fi7r3pck8DAQCk+Pl4aMWKEtG/fPvNyIxxjkyRJkvvlMERERETqYRsUIiIiEg4DFCIiIhIOAxQiIiISDgMUIiIiEg4DFCIiIhIOAxQiIiISDgMUIiIiEg4DFCIiIhIOAxQiIiISDgMUIiIiEg4DFCIiIhIOAxQiIiISzv8DgKZzjB14M14AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(angles[:,0])\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/docs/contents/user/tools/structure/index.md b/docs/contents/user/tools/structure/index.md index 031324f5d..324ba1e15 100644 --- a/docs/contents/user/tools/structure/index.md +++ b/docs/contents/user/tools/structure/index.md @@ -5,6 +5,7 @@ | [Align principal axes](align_principal_axes.ipynb) | Aligning the principal inertia or geometric axes of a molecular system over a reference coordinates axes| | [Center](center.ipynb) | Centering a molecular system | | [Flip](flip.ipynb) | Flip a molecular system over a plane| +| [Get angles](get_angles.ipynb) | Getting the angles between specific triplets of atoms of a molecular system | | [Get center](get_center.ipynb) | Getting the center of a molecular system | | [Get contacts](get_contacts.ipynb) | Getting the contact matrix of specific elements of a molecular system or two different molecular systems| | [Get dihedral_angles](get_dihedral_angles.ipynb) | Getting the dihedral angles of a molecular system | @@ -33,7 +34,8 @@ align_principal_axis.ipynb center.ipynb - flip.ipynb + flip.ipynb + get_angles.ipynb get_center.ipynb get_contacts.ipynb get_dihedral_angles.ipynb diff --git a/molsysmt/_private/digestion/argument/triplets.py b/molsysmt/_private/digestion/argument/triplets.py new file mode 100644 index 000000000..6478e6972 --- /dev/null +++ b/molsysmt/_private/digestion/argument/triplets.py @@ -0,0 +1,21 @@ +from ...exceptions import ArgumentError +import numpy as np + +def digest_triplets(triplets, caller=None): + + if triplets is None: + return None + + if isinstance(triplets, (list, tuple, np.ndarray)): + if not isinstance(triplets, np.ndarray): + triplets = np.array(triplets) + shape = triplets.shape + if len(shape) == 1: + if shape[0] == 3: + return triplets[np.newaxis, :] + if len(shape) == 2: + if shape[1] == 3: + return triplets + + raise ArgumentError('triplets', value=triplets, caller=caller, message=None) + diff --git a/molsysmt/structure/get_angles.py b/molsysmt/structure/get_angles.py index e0ddd6b27..3f9369311 100644 --- a/molsysmt/structure/get_angles.py +++ b/molsysmt/structure/get_angles.py @@ -4,9 +4,8 @@ from molsysmt import lib as msmlib import gc -#@digest() -def get_angles(molecular_system, triplets, - structure_indices='all', pbc=False): +@digest() +def get_angles(molecular_system, triplets, structure_indices='all', pbc=False, skip_digestion=False): from molsysmt.basic import get diff --git a/tests/structure/align_principal_axes/test_align_principal_axes_molsysmt_MolSys.py b/tests/structure/align_principal_axes/test_align_principal_axes_molsysmt_MolSys.py new file mode 100644 index 000000000..ca4a2be96 --- /dev/null +++ b/tests/structure/align_principal_axes/test_align_principal_axes_molsysmt_MolSys.py @@ -0,0 +1,37 @@ +""" +Unit and regression test for the get_structure_alignment module of the molsysmt package on molsysmt MolSys molecular +systems. +""" + +# Import package, test suite, and other packages as needed +import molsysmt as msm +from molsysmt import systems +import numpy as np + +# Distance between atoms in space and time + +def test_get_structure_alignment_molsysmt_MolSys_1(): + + crd = msm.systems['POPC']['popc.crd'] + psf = msm.systems['POPC']['popc.psf'] + molsys = msm.convert([crd, psf]) + + axes_1, momenta_1 = msm.structure.get_principal_axes(molsys) + + good_axes_1 = np.array([[[ 0.09213962, -0.02624376, -0.9954002 ], + [ 0.56239734, 0.82631277, 0.03027277], + [-0.82171742, 0.56259975, -0.09089556]]]) + + molsys_2 = msm.structure.align_principal_axes(molsys, axes=[[1,0,0],[0,1,0],[0,0,1]]) + + axes_2, momenta_2 = msm.structure.get_principal_axes(molsys_2) + + good_axes_2 = np.array([[[ 1.00000000e+00, 2.87221485e-16, 8.58936055e-17], + [-2.87221485e-16, 1.00000000e+00, -3.60822483e-15], + [-8.58936055e-17, 3.60822483e-15, 1.00000000e+00]]]) + + assert np.allclose(momenta_1, momenta_2) + assert np.allclose(axes_1, good_axes_1) + assert np.allclose(axes_2, good_axes_2) + + diff --git a/tests/structure/get_angles/test_get_angles_from_molsysmt_MolSys.py b/tests/structure/get_angles/test_get_angles_from_molsysmt_MolSys.py new file mode 100644 index 000000000..d170ccfc7 --- /dev/null +++ b/tests/structure/get_angles/test_get_angles_from_molsysmt_MolSys.py @@ -0,0 +1,32 @@ +""" +Unit and regression test for the get_dihedral_angles module of the molsysmt package on molsysmt MolSys molecular +systems. +""" + +# Import package, test suite, and other packages as needed +import molsysmt as msm +from molsysmt import systems +import numpy as np + +# Distance between atoms in space and time + +def test_get_angles_from_molsysmt_MolSys_1(): + molsys = msm.systems['pentalanine']['traj_pentalanine.h5'] + molsys = msm.convert(molsys) + angles = msm.structure.get_angles(molsys, [0,1,2]) + check1 = ((5000,1)==angles.shape) + good_angles = [113.06116981671607, 107.22517056178846, 116.89755264175076, 114.74587424352482, 107.95984387161585] + check2 = np.allclose(good_angles, msm.pyunitwizard.get_value(angles[:5,0], to_unit='degrees')) + assert check1 + assert check2 + +def test_get_angles_from_molsysmt_MolSys_2(): + molsys = msm.systems['pentalanine']['traj_pentalanine.h5'] + molsys = msm.convert(molsys) + angles = msm.structure.get_angles(molsys, [0,1,2], structure_indices=[0,1,2,3,4], pbc=True) + check1 = ((5,1)==angles.shape) + good_angles = [113.06116981671607, 107.22517056178846, 116.89755264175076, 114.74587424352482, 107.95984387161585] + check2 = np.allclose(good_angles, msm.pyunitwizard.get_value(angles[:,0], to_unit='degrees')) + assert check1 + assert check2 +