Skip to content

Instantly share code, notes, and snippets.

@ljmartin
Last active September 3, 2024 20:55
Show Gist options
  • Select an option

  • Save ljmartin/b253152f5824d3922e7c73737cdb0b0e to your computer and use it in GitHub Desktop.

Select an option

Save ljmartin/b253152f5824d3922e7c73737cdb0b0e to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{"cells":[{"id":"801fb9f5-7773-4f04-b9f1-31bbadddb774","cell_type":"markdown","source":["# Matching ESP contour figures\n","\n","From the paper: Converting SMILES to Stacking Interaction Energies"],"metadata":{},"attachments":{}},{"id":"0","cell_type":"code","execution_count":1,"outputs":[{"output_type":"execute_result","execution_count":1,"data":{"text/plain":["0"]},"metadata":{}}],"source":["from rdkit.Chem import rdEHTTools,AllChem\n","from rdkit import Chem\n","from scipy.spatial.distance import cdist\n","import matplotlib.pyplot as plt\n","import numpy as np\n","\n","mol = Chem.MolFromSmiles('c1cnccc1')\n","molH = Chem.AddHs(mol)\n","\n","# Calculate coordinates. Normally you'd want to do 3D coords\n","# to get accurate values, but I'm only interested in a planar molecule \n","# and it's just for fun. \n","AllChem.Compute2DCoords(molH)"],"metadata":{}},{"id":"882664b1-05e0-4ce4-b3f4-6dca96ae7ce7","cell_type":"code","execution_count":2,"outputs":[{"output_type":"execute_result","execution_count":2,"data":{"text/plain":["array([ 0.0772853 , 0.46514589, -0.89892606, 0.46514589, 0.0772853 ,\n"," 0.24763075, -0.08665155, -0.07808348, -0.07808348, -0.08665155,\n"," -0.10409697])"]},"metadata":{}}],"source":["# hint from: https://github.com/rdkit/rdkit/discussions/7196#discussioncomment-8624209\n","_, res = rdEHTTools.RunMol(molH)\n","static_chgs = res.GetAtomicCharges()\n","static_chgs"],"metadata":{}},{"id":"8731e18f-4328-492d-864d-304f8c4ee15b","cell_type":"code","execution_count":10,"outputs":[],"source":["# hints from:\n","# - for using meshgrid around molecule coordinates https://ljmartin.github.io/blog/04_meshing_ses.html\n","# - for calculating ESP - just using a simple approximation taken from ChimeraX:\n","# https://www.cgl.ucsf.edu/chimerax/docs/user/commands/coulombic.html\n","\n","xyz = molH.GetConformer(0).GetPositions()\n","xy = xyz[:,:2]\n","\n","lo = xy.mean(0)-5\n","hi = xy.mean(0)+5\n","\n","n = 100\n","grids = np.meshgrid(\n"," np.linspace(lo[0], hi[0], n), \n"," np.linspace(lo[1], hi[1], n))\n","pts = np.vstack([i.ravel() for i in grids] + [np.array([3.5]*n*n)]).T\n","\n","distances = cdist(pts, xyz)\n","\n","esp = (static_chgs/(4*distances)).sum(1).reshape(grids[0].shape)\n","#convert to hartree then kcal/mol:\n","esp /= 0.529177 #divide by bohr\n","esp *= 627.51 #hartree to kcal/mol"],"metadata":{}},{"id":"5be3b8d1-55c0-48e5-988f-847d6b2885b1","cell_type":"code","execution_count":15,"outputs":[{"output_type":"stream","name":"stderr","text":["[06:54:46] The new font size 0.8 is below the current minimum (6).\n"]},{"output_type":"execute_result","execution_count":15,"data":{"text/plain":["[]"]},"metadata":{}},{"output_type":"display_data","data":{"text/plain":["<Figure size 640x480 with 1 Axes>"],"image/png":["iVBORw0KGgoAAAANSUhEUgAAAYUAAAGFCAYAAAASI+9IAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABkU0lEQVR4nO3dd3hV6Xnv/e9auzdtdZCQkEAgQAjR+zDA9GaPxzO2xyUusZMc503snOQkseOTOO/Jm8ROr8dJnMRxG9dp8XQYht6LCiCBQB3Uu3bfe633D4FGwwiQQNLa5f5c11z2NaNyS3rW+q2nrOdRdF3XEUIIIQDV6AKEEELEDwkFIYQQYyQUhBBCjJFQEEIIMUZCQQghxBgJBSGEEGMkFIQQQowxT+aDNE3j6tWreDweFEWZ6ZqEEEJMM13XGR4eJj8/H1W9eX9gUqFw9epVCgsLp604IYQQxmhtbaWgoOCm/31SoeDxeK59sRWkpZmmp7JUUHcBfm2J0VUIkXz+9QIslWtrKoaGYhQW1ozdz29mUqFwfcgoLc0koTAVbgVM8vsSYtq5FZB70R253RSATDQLIYQYI6EghBBijISCEEKIMRIKM+V8LXx6mdFVCCHElEgoCCGEGCOhIIQQYoyEghBCiDESCkIIIcZIKAghhBgjoSCEEGKMhMJMKVsG36s1ugohhJgSCQUhhBBjJBSEEEKMkVAQQggxRkJBCCHEGAkFIYQQYyQUhBCJRVb1zSgJBSFE4imTHYhnioSCEEKIMRIKM0leYBNCJBgJBSGEEGMkFIQQQoyRUBBCJA4Zjp1xEgpCiMQiK49mlITCTJPJZiFEApFQEEIIMUZCQQiRGKTHPSskFIQQiUPmE2achMJskHkFIUSCkFAQQggxRkJBCBH/vlcrQ0ezREJhNskQkhAizkkozBZ5yhHizsjD1KySUBBCxD95qJo1EgqzSVYhCSHinISCECJ+yUPUrJNQMII0dCEmT4aOZpWEwmyTBi6EiGMSCkaR3oIQtybvJhhCQsEI0tCFEHFKQkEIEX+kl2AYCQWjyPJUIUQcklAQQsQXeVgylISCkaS3IMTEZOjIMBIK8UCCQYhRMpdgOAkFo8kFIISIIxIK8UJ6CyLVSS8hLkgoxAO5EESqk4eiuCGhEE/kwhCpTB6O4oKEQry4fkFIMIhUI20+rkgoxBN5UhKp5nogSNuPGxIK8UienEQqkUCIKxIK8UaGkUSqkDYelyQU4pE8OYlkJ8NGcUtCIZ7Jk5RIZhIIcUlCIV7JMJJIVvKSWlwzG12AuIWyZXBeQmG6jMRitEcidEQitEciDCoKYUUhrKqEAJOiYNU06vx+aoaH2eFy8VR6OnMtFnItFiyKYvSPkPjkISfuSSjEu+s7qX5anqwmK6LrVPr9nAiHabLZuGqz0eNwEHO58Hq9pKen4/F4sFqtqKo69o+u62iahlfT2KJpDPj9/P3QEAMDA4wMDZHm95MbClEQClEaibDF5aLQYkGRsJgcmUdICBIKiUKC4aZius4xn4+3NY06t5suj4e5ZWXk5+fjdrmosNvfvXFHo1hjMWyxGNZYDDOgahpmTUMDYiYT1ZEIqqqSnp7OwrlziZjNmMxmdF0nGo0SCAQ4MjjIj5qaiHR1Uezzsdrv51GnkwKr1chfRfyTQIh7iq7r+u0+aGhoCK/Xy+DgKtLSTLNRl7jR9WEkCQYAwprGG8PD7LJauej1kltSQlFREWlpaaiqCqEQGaEQmdEo2bEYeSYTuRYLTnXq02iarjMQi9ERjdKh6/SpKn02GyN2O6rJRCgUoq+vjwsXLmDr6GD14CBPW60sczhm4CdPUDKPYLihoRhebyWDg4OkpaXd9OMkFBKJBAMnfD5+DlRlZlJUXs68efNwuVzooRC5gQBF4TALLRbmmM2oMzysE9I0msNhLus6bXY7I04nKAr9/f1cuHABtamJ7QMDPOvxkG1O4U65DBvFBQmFZJWCwRDRdV4YGuJ5txt1yRJKS0vxer0ogQCFPh/LgRKbDZPBY/vDsRjnwmEuWiz0ut1oQFdXF+cqKym7epVfVlXKU633IIEQNyQUklmKBENQ0/j3oSFey85m/qpVFBUVYVZV5gwPszIaZandbngQ3MxwLMapcJg6u52Qx8Pg4CBnz57FXV/P/4hE2OJ2G13izJNAiCsSCsnufPJOPIc1je8ODfFidjaLNmygsLAQNRSidHiYrTYbaabEaYO6rnM5HOaoqtLj9RIIBjl37hz2s2f5zViMtS6X0SXODAmEuCOhkAqSMBheHx7m/6alMX/DBhYsWIAaCLB6ZIRNDkfCvyfQE42yJxrlSkYGgWCQ6upqcs+d439bLOQn06olCYS4tPrvz1P5W4HpDYVVf+fgzJfLprVQcReSaBipKRTi/8RihNesYcWKFZijUSqGh7nHbk/4MLhRVyTC25pGZ2YmAwMDnDp8mEdaW/n1tLTE/1klEOLOun8b/ZvEAqVU/tY09xRW/d0qTI6LAJz8Vfmjx4UEDwZN1/nXoSFeKS5m3ZYtpLndFPX28ojNhv0Olo8mkuZQiNfNZkJeL42NjbQdPMifRCKsSNTJaAmEuHM9EGAZsUBspkLBBLz7qrqEQxxI0GBoDYf5SiyGa+tWlixZgnNggMd1nXyLxejSZo2m6xwJBjmVloYvFuPEkSPsrK/nS2lpcTuJPiEJhLgyPgyum2wo3OGj2LKxb/buNxeGScDN814ZHubXMjMpefpplpSUsKq7m8+bTCkVCACqorDV4eAzwSBzAwG2338/9Q8+yCeDQTojEaPLmxwJhLix7t9qJwyEqbjL/vloOLy3EGGIBAmGmK7zJ4ODfG/lSnY+/ji5us7HhofZ7nDM+Mtm8cxrMvFZq5X13d0sKi5m0VNP8Tmnk6M+n9Gl3ZoEQtx4bxjc+d/jDoePJiJDSnEhjoeSBmMxfiMUwrNjB4sWLWJedzcfioOXzuJNeyTC82YzQbudIwcP8uTFi/zyLbr7hpFAiAuT7RlMdvhoGt+9v17Qu70GCYd36brOlSsR6uqCNDXF6Okx0d9vwe9XiURUIhEFVQWLRcNm00lLi5KZGWXOHI3SUhtLlthwuSaxHHh8jyGOguFqOMwXVZWlH/gAc7Ky2NjdzUa73eiy4lKexcIXNI2fDA1x786d7ElP5+rx43zN44mfHVmnGAi6rtPePtr+Gxtj9Paa6esz4/erhMPvbf9Wq0ZamkZmZoTcXI3SUitLl9pxu2U5/HjvHZ2Zvmt9GnsKN0rtYPD7NfbtG+bkSYX6ejft7S6czkzy8vLwer1YrVYsFgtmsxlVVccu9lgshq7rhMNhIpEIgUCA7u5uOjo6cDiGKS72sXx5gHvvtbNsmf3WN4k46TXUBYP8jsvF2ocfJt1m40m/n6JkWpc/Q3Rd5/VAgPrcXOrq6nDu3ctfp6VhNjoYJhEIwaDG/v0jnDgBFy86aW93Y7NlkJeXR0ZGBhaLBavV+r72r2kamqYRiUQIh8OEQiG6u7tpb2/Hbh9m/vwRli8Pcs89VlascMRPSM6iOw2DGV59NBWpEw69vVGef36EI0dctLVlUlKyhJycHNxu97X9+mPY7RFcrghOZwynU8Nu1zGbdcxm0DSIxRQiEQW/X8HvN+HzmRkZMaNpFlRVJRAIMDw8TGNjI4FAC2VlQzz+uMa2bW5UdYILxOBgqAkE+P3MTLY8/DDuWIxno9HU3hzuDhwIBDidnU1TczP+Xbv4F5fLuPcZbhEIAwNRXnjBx6FDDlpaMliwYAm5ubl4PB5MJhOa9m77dziiuFwTt/9odLT9+3wm/H4TIyMWNM2KoiiEQiGGhoZoampieLiZsrIhHnkkys6dHkym5A6Iu+0ZxFEoXJec4aBpOm+8McyLL9pobZ1LWdkKcnJysFgsmEwBcnNDzJsXobjYwty5ZszmO2u4w8Mara1hmpp0OjqsDA66UBTTWED09l5g8+Y+PvtZBwUFNzyFGxQM5wIBfi8ri62PPIInEOBTioIngbaoiCdVgQB7s7Joamkh+NZbfMvtnv0ewwSBoOs6b789ws9/bqGpaQ5Ll64gNzf32gFGo+0/Pz/CggVm5s61YLHcWc0+n0ZLS5jmZp3OTht9fXYUxcLIyAjNzc10dtayYUMfn/2sjeJi23T8tHFjuoaJ4jAUIJkmo4NBje9+d5hXXvGSk1PBggULcDgc2Gw+SkoClJebycszz1j3NhTSqKsLUldnpr3dja6b6evro7q6kqKiFj7/eY21a8ftqzPLwVAXDPLb6enc89hjuAMBPquqd3SWgXhXTTDIO5mZNDQ3E37rLb7l8czeJP0NgRAOa/zwh8O8/HIaXm85ixYtwul0YrGMUFISoKzMREHBzJ1KFw7rXLwYpK5O5coVN7pupa+vj5qaKvLymvjc52Js3pzYmw5O95xBnIbCdYkbDuGwxne/O8SLL2ZTUrJ+dLM2NUZR0RDr149eCLNfk05VVYDqahsjI16Gh4c5e/Ys6em1fPnLsGLFtTdkZykY2iMRPm+3s+UDH8ATDvNZRcElgTAtqgIB9mVnc6mhAfeuXXwzLW1mx9VvCINoVOe554b4yU8yKC5eT1FREaqqUVg4yPr1KoWFs388aSSiU1MToLLSxtBQGn6/f3TTQftZvvSlGx6OEsBMTSDHeShcl1jh8NZbw/zjP7qZP38TxcXFmExBVqzwsXmzDbs9Pm56ra1hDh7U6erKZHh4hMrKShYurOVrX3OQlWWe8WAYicX4VCzGiiefJMts5jOahkcCYVqdDAY5lJ1NdXU16w4f5je93pn5RjcEwv79I/z1XzvIy9tISUkJJlOIsrIRtmyx4XTGx9/46tUIBw5otLdn4PcHqK6uJi/vLP/7f9uYMye+X4ycqTC4LkFCYbz4nXPo7Izwh38YZmRkNRUVFVitMSoqhtm61X7HY6Qzrb09wu7d0NeXQU9PD5WVB3n22Q4+/elrT5YzEA5RXefzfj85jz9OflYWz/r95Mqk8ozYFQhQm5vLkcOH+URNDU95PNP7DcYFQm9vlD/6oyDd3RWsXr0aq1Vn+fIh7r3XjtUan+2/qyvC229rdHVl0dvby+nTh3jqqSv86q96427F0kyHwXUJGArXxVfv4cUXh/i3f5vD2rXbycjIoLCwh0cescbNk9Ht1NcHeestG5GIm/Pnz6NpR/jmN62jT03TGAwxXefXh4exPvwwxYWFPDEwQIktuSb84omu6/w8GORKdjYH9+3jf9TX88h0HdwzLhBee22Yf/zHTFav3klmZibz5vXw6KM23O7EaP+NjSHefNNMKOSlrq6OQOAw3/iG6f2LMQwwW2FwXQKHwnXGhkMopPEHfzBCR8d6Vq1ahdU6wgc/GKWw0PjGNFWxmM6ePUHOn8+iv3+QU6f28nu/18POne4pBYNf0+iIRGiPROiIRBjRdcKqShCoslqx7dzJwgUL2NHTw8pE3ekzgcR0nR8Gg/RmZXFgzx4219eTDlh1HZumkWU2M9diYa7FQq7ZfPtJ6XFhEInofP3rQzQ0rGXNmjVYrX4efzzMggWJF/SaprNvX5Dq6kwGB0c4cWIvX/pSF48+asxE9GyHwXVJEArXzX44tLWF+fKXdQoL76ewsJAFC7p59NH4HSqarLa2MC+/bCYUcnHkyGG2bDnLZz7jIjf3vb2GgWiU434/5zWNFpuNdrudXrsd3enE6/Xi9XrxeDxYraNrx1VVxWw2k56ezraeHtbIm8qzJqbrfD8UYiAri76+PnRdR9M0YrEYwWCQgYEBBgcH8Y2M4A4EyAkGmRcKURwOs8Zmo9zhGF3aei0QtKVLqakJ8Gd/ppKbu5Pi4mIKC3t44glb3A4VTVZ7e4SXXlIIBDycOHGcNWuq+P3fT5v4/Z5pZlQQjJdEoXDd7IRDdbWfr341nXXrHsTrdXDffYOUlyfPTS4c1vn5z4P09OTS3NxMfX09FksPCxf6MA0PcLkjG18ol4UlJWRmZuJwOLCPu8nrkQi2WAxbLIZF1zGN+2eFrrNM3lSedRFd551wmH5VJaaqRIGYohA2mQipKhGzGZPZPBYYfr+fQCBAa2srPVeukD+ni/m5/VwNzqW11UNmZgGlpaWkp7vYtq2f1auTp9cXiei89FKQ9vZcLl68CLzDP/yDC6t1ZobD4iEMrkvCUBhvZgJi374R/uIv8tm69QFcrhAf/ahGTk7yTZTqus6hQ0Gqqx0Egw5UVcXv96PrOk6nEwCLyUfGgEZmOEy2rpNnNpNlNmNXlLibqBO3puk6Q5pGVyRCh6bRazYzYDYz6LaCYiMWi+H3+3E4HJhMJiBKRoafhx/WycuL7xU7d+rw4QAnTmTT0tJKV9cu/vVfbdN61HA8hcF1SR4K101fOOzZM8Lf/m0x27bdh8s1xKc+ZcLlSozJtLvh82lcvhykudkE6BQUaJSU2Egb6Rr9gO/kGVqfmCGfa0fT4CrZNDZG6Ooyk54epaTERGGhNem3jAA4dy7I7t0ZdHb20NDwGt/5jgWP587vb/EYBOOlSCiMd+cBsX//CH/5l6OBkJ4+wC/9kjXh5w+mxdX20f+VYEgen7v2N82XvylAU1OI//7vNDo6+mhsfJX/+i/r5HYjHifew+C6FAyF66YWDmfPBvj935/D9u2P4PUO8OlPSyC8j4RDcpBAmNBoMHi5erWb7u5X+c//dN62p5QoQTBeCofCeLcOiKtXw3zhCw7uuedJvF4/v/zLZgmEm7keDCDhkGgkDG7r8uUQr7ySTkNDM07nW/zFX7x/+5D3ny6ZGGFwnQGH7MSj9x/8A6MBEQ5r/OZv6qxd+xBOZ5hPfcokgXAr128oV9tHbzISDPHvc+OCXALhlkpKbGzb1g8spLp6C//xH4f5whdGtw9JxF7B3UjyULhu/B9yNCAGX3HyaPHHychw8dGP+nC5UuRXcbfy894NBpBwiFfSO5iyNWvsdHd3ARU8/1IXf9O0F+d8lVQIgvGSf3nN+yxjuKaAda4HKCoqoiuniW/XdRtdVGLJz3v3ZjP+aVQY73PtEgh36E/3t3PCPkhvrJ1tW7bhPLicWKjU6LJmXco9Hkd9UbLOZbHm0TV0m7txzBu9cP50/7s3t6/dKxfTpNwYDNJrMI4MFd2R8dc9gKLkYanQiFXF2L52O7/Y/QvSHr/5+HsySrlQCO4K8ujGR/HFfNjLxr+p/O6FJAExRRIOxpEwmLIbg2D8tQ+gWlX88/zkm/OpaK7gbP1Z3IsT+8CeqUipUBi5NEKFu4KsrCz65vdhM91sc6/rjaRdAmIqZL5hdskw0ZS8Nwxu/Tuz59vp6e5h7dq1NLzVgL5QR0mBF/oghUJB13TsJ+ysenAVvZZe7NmT2c9ofMORgJiUG1cpgYTDdJMwmLSpBMGNrEut6NU6m5ZtYtfRXXi3ztBhRnEmZUJh+OQwO0t3olpULIvvZD8XCYgpkXCYfhIGk3I3QTCealMZSh9ioWkhua/m4l/rx2RPpPe07kxKhIIe0/Fe9lLyaAl97j4c9rvd9VECYtIkHO6OzBlMynQFwY3sJXb8p/xsWLWBV46+gndH8vcWUiIUhk4Pcf+S+4noEWwLp/uQEAmISZFwmBoJg9uaqSAYT1EVAjkB5pnmkVuZiz/kx2RL7t5C0oeCrut4GjwseGgB/Z5+HJaZ3BteAuK2xt/gxt/4JCBGSRjc0mwEwY1s8234u/2sq1jHq6deJX1L+qx8X6MkfSj4GnysmbMGFLDN6lGCEhC3Jb2HUTe+AChh8B5GBMF4iqrg8/qYp8zDW+NF36wn9ZkiSR8K6jmV0nWlDFmHsNqMOhVMAuKWUrH3IEFwS0YHwY3sxXa0Ko1l+cs4efkk7kXJ+95CUoeCFtbI9mfj8Xjozes1upxrJCBuaaLeAyRHQEgQ3FK8BcF4qlVlyDpESUkJpypPwSKjK5o5SR0Kw2eH2bxsM/6oH2tOPJ4dLAFxU+NvmDcGBCROSEgQ3FI8B8GNtBwNT8xD5nAmwVgwaV9mS+pQsLXZmLtxLj6nDwfxfvi4BMRN3XgjjeeQmGiDQAmC90ikIBjPPtdO6EqIJcVLOHTpEJ4lHqNLmhFJGwq6puP1e3E6nQRyA0aXM0USELc0mZCAmQ2Km+0OKwEwoUQNgvEUVSFgDTBv3rzR87uWGF3RzEjaUPA3+ymbV0YoGsKaFY9DR5MlAXFbN7sRj7txRyIa/oCOpoNJhbTbHNAeieiM+DVUBZwOBYtlgl3mJQBuKRmC4EZRb5SsaBaegeTsJUASh0KkOUJBUQEBSwCzmiw/pgTElIy7adfVRPnRjx7G77fgdvv50pfeITf35sFQfz6Tf/u3DTidET7xiTcpL0+WNjSzkjEIxrPOtUIvZCgZ9AR6Eux44slJ2pbuGHTgdrsZtg9jTsofUwJiKvr7NU6cmMfQkA2TSSMry8+v//pxbLaJz5nq77dx5EghaWkhHn1Um+VqE0uyB8F4JruJoB6kuLiY5uZm0pYm31kLyXi3BMAZcGK1WlHSk3OFwHtJQEyWouhomsLu3UvYsuUsGzcGjS4pId3uTIJkFrFEyMnJgUZgqdHVTL+kDAUtrJGupgNgzUjk+YQ7IQFxK4oCHk+I3l4H3/nONsrKXsPjuZNdc1NPKgfBeGF7GLfTjWUkOdtNUoZCsCtITk4O4VgY9SbDA6lh4oBI5XAwmTSefLKOn/yknLNnC3n11YV89KMtqGoq9CinToLg/TS7ht1uxx6azJksiScpQyHUHSInJ4eoOWp0KXFEjhu9buXKI3R0ONm9ezEvvbSRioomysqSb8LwTkkQ3JrFY0Hv0bFH7YQIGV3OtEvKUFB8Cs58JzE1hgm52N8vtQPCbg/xzDNHqanJo6kpnR//eCd/+If7sFhSt7cgQTB5JufoPcWKNSnfbE7KUDCFTVgsFmImCYXbS8WA0Fi+PMAjj+znu999jAMHFnLwYCU7dgwm9e6XN5IguDOqVSUWi+Fyuujz9WFJS665haQNBZPJRNQkw0dTkzoBYTarPP10K/v29dLYmMUPf7idJUteID8/uS7wG0kQ3D1FVdAUDbfbTdQXlVBIBIquoKoqpPIc811L/oDIylL5xCdO8rd/ex91dTm8/PIyfuVXLmI2J1fDkSCYfrqiYzab0aO60aVMu6QMBTRSahhg5iVnQKiqwo4dDRw/ns7u3Rt54YXNbN/eRllZ4r+7IEEws3R0VFVF1yQURMqbOCAgMUPC47HwyU+e5/jxCoaGHDz33Gq++tWDRpd1RyQIZpeiKJB8mZCcoaCbdGKxGIouvYWZlRwvypWWRnn66Rq++931HD26mIMHq3C5EmMISYLAGIquEAqFkvI9qKQMhZgphqZpKDEJhdkTfwERCMTQdRWn89btQFUVHnuskuPH8zl7toCf/ewBnnnm/Ps+zu/XgRhOp3GXjYRAfFBRCQaDmGzJt7oxOUPBEiMajWLSku8Plhjee6Oa7WGmaFSnutrBCy9sxu0O86Uv7b/t5xQUqHzyk/v54z/+GOfO5eLxhN/z3/3+GP/wD/fS32/nwx8+zJo1/ll5r+H9IQASBMbSIhpm1Yzf7x97ZyGZJGUoaC6NYDCIN+Y1uhQB3KoXAdMXEpqm09mp8bOfrWP37qV0d7tISwvx4INngMHbfv7mzcNs2XKZvXtLOXEi/z3/ra7OwzvvLGRoyMa5c3O4//46PvaxE8ydq077FhnSG4hvsUAMAH/Yj8kuoZAQrBlWBnsHyZibgYZsexxfbrzBTU9I9PXFeOutRfz0p+tob/ei62CxaCxb1o3FMrn3Vex2E5/+9HHq6+fS2vreBwqzOUp5eRfHj8+jp8fJT3+6hgMHSvjIR07xwAP15OTc2c1BegKJJ+YbDYWwJXybj0xMSRkKthwbnec7KV1SSiAWSLrX0JPL+2+AE90obxYUkUiM48fd/OhHO6mpmUMwaMZk0pg3r4+PfewUO3deJjvbwv79k1smsnixn0cf3cd3vvM4kci7N/qKijBf+9qr7Nu3gJ/+dA3NzdlcuZLOt761nf37y/jYx/awZcswVuutw0FCIPFpIxqxWIyQJYSF5HpxDZI0FMxpZvp8faiKSngwjC3TZnRJYkpu35vQNJ2uRjPZLR/mnXdKCAYtKArMmTPCAw+MDu3k5ppQlKldtGazykc/2s6hQ52cO/feIaTMTDMf+lALW7c2jg1RdXS4OXNmLrW1H+Heexv40If2sXKlwjcOdUzyZxOJxhww4/P5CLqCEgqJQlEUgrYguq4TG4hBptEVibvz7o1U13X8QybqDi7g4tFyhnvTAVBNQeaX11G4cw/9C0b49wsKXHi3h1FWpvH1r3+LWEynpOTWF7LHo/KVrzzPlStRLBaFpUvfffpXFIVv13UTW/YaBc79DLyznebq5cRGTDS+sJDn92axaP1hyu414UqPyUuUScgStTA4Mkg0Izm30UnKUADwuX0MDw9jtiXtj5hy/ENRWs5mUr3rQQa70tF1BdWkkbugifKdRylcPozFmga8e0Ti+4ZrFDhYM7nvF4lpmFWFfefff2M3mfOZWwLZhSdpq7vA2Xc20nFpASN92VS++UEaTg9S8cBbFFf04vRKG0wWekzHoTmobKzEucFpdDkzImlbq5an0dPTQ74nXyabE1wsptHd5ODMG/fTXl9INDzabF3pI5Tde5ZF6yvxZJlhwh1x72y4JhyMUbVrCQ5PN4vW92F3TTxXYLaaKK7wk124m4ZTKzi7dwUjfR6Gur0c/flTNFW2sfLBg8xd5MeUZHsqpaJQdwi36qZjqANrVnKe6pi0oeBe5Kb+lXpKSkoY9A1idiXtj5q0tJhGf4eDs3tW03imlHBw9CJ0pAUorjjJygcv4s7SUNXp/9u21RZw7p0tRMNmLhzuZPWjR5m3pAvbTdaluzNMlN93jgWra6l5exENZ9bjH3TSVltMZ8M8iiouUfHASTLyA5hMEg4Jqw8CgQAjzhHsipy8llBMDhO9Si+xWIzI1QjmxUn7oyaloE+j9sBc6o/fy0BHOgCqSSNv8RUq7j9J3uIuzFYTMP1j9rquEw5YUc1RtKCV3rY89v/gcQrLWql4cC9ZBdqE7yaoqoInS2fDUxeYX1HL2T2PcKWukEjIwqUTy+huzqdk3X7Kd7Rjcyky35CA7H473d3dhPJD2JFQSDjDOcP09fVhNyXnHy8ZBX1hWs4u5Ow7q+hry0XTVFSTRua8Xsp37qd4ZS9Wu8rEQ0XTQ1EUFm9sJLe4kTNv7KDtfBEhv43Lp0q4cmEuizfUsnTrCbxzzBOGg9miUrBUJbf4LVrPZVK9azu9V7IY7PJy5o0naKrqYsXOSopWXsbuSs4hiGQUGYqQrqazv3Y/rsddRpczY5I6FBzLHdQdq2Nr9laG/cOYDdyzRtyapml0Xu6k9mAnVy8+jX9wDgAOT4DlOw6yaH0DadkmZuuQDJNJJTMftn1iL1dqM6h5Zysdl/IJjrio2bOOqxcXsGTLEUo3dWC9yTOH1a6ycE0/OUUvcPnkQs7t3Yp/yEnflVxO/GIBbXUvsOyeecxdNHf0/A8R16JXooTDYXrUHizO5FuKel1S3yVt2Taa/c1sjm0m0hLBvDSpf9yEpGkaAx0DXDx6kbpDdUSCURyen+H0fpr55ZdZvuMcGXP9qCZjthOw2lWKVw0wd9HLXD6VT+2BzQx05NDblsXR5x+j8Uw7Kx/cxdxFIawTbHmgKApp2SZWPthI8coOzu5dTnN1Eejf5/LJM7SePcfiTYtZsnkJmfMyJRzilK7ruEZcXLlyBX+RHy/Ju4VO0t8lffN9XLlyhRwlB03XZBw3joT8IS6duMS5fecY6BgY+/fZ8ytZvl0hf4nl2oodY2+UiqLg8Fgo39FNfunLXDi8iQtHSgkHbLTXz6O37RMsWF1L+c5TZOTFJp5vMKlk5IXY/MxJiisOUXeolcZKCAfDnNt7jtZzrSzfvpxFGxbhcDtm/4cUtxTqCJFtzqa6thrPhz1GlzOjFF3Xb/v+/9DQEF6vl1V/twqTI7E2gIqFYrhecvHUY0/Rnd2No1AuOKOFAiHa69upfKOS3rZeYtEYiqqQmZ/J0q1LWbR+EVaHNW4DPBqJ0dvqomrXCq5eWEYoYENRdJxeP0u3nmHJ5gu4MiYOB7g2kR0M03CqgdqDtfS29aJrOqpJJXNeJqsfXk1eaR52l8yFxYvo6Shav8bPz/0cz2OJGQqxQIzK36pkcHCQtLS0m35c0ocCwOBbgzy96GlsmTaUdQrKDKxYmcgc7eOz8n1up1P9kdElAKM3w56WHmr21NBc3UwkFAHA7rZTurGUZfcuw5uTON3ycDBGc1UxNe+soaclBwBF0cku7GHZvbtZtH4Ys+XWvZzhvmFq99dy8dhF/IN+AMxWM/PL51O+s5w5C+fEXTimWrsO94fJaMxgzzt7aN3aii07MbfNkVAYJzIUIfvNbB5/9HG6c7ux50/PE9jtLo4899xp+T53q33kZvvwvGsmLzBN0xjuHab2QC2XTlwau/lZHVYKlhWw6uFVpM9Nx2xJvNFMTdMJjkD9sSLO79/MSJ8bXVcwW6PMLWli5YNnmLOwF7P15j9bLBpjsGuQM2+coa22jZAvBIAzzcnCtQspu7eMtJy0WZlvmMwNPx7a9e3a9HS25/DpMOZhMz8/+XMcTybuSIOEwg2GXh/iqdKncGQ6YB1Tfvq62cUSDxfI3ZrJCywcDHPp+KWxYZLrchfkUvFABYVlhVhsib+SQ9N0+q7o1B5cxaXjq4mELIBORt6rzFv6Iku3LiEz/9abcEXDUdrq2qh5u4aOSx1cvzQz8jJYunUppZtLsTnu7ik1UR5k7tbN2vRU23KoO0RWWxZ79+2lcV0j9rmJO6QnoXCD6EiUtFfTePKxJ+lO78ax4OaJn8wBcCduFRo3u8gioQgdlzuoequKrsYuopEoiqrgyfJQvrOcBasX4PIm31rvSDhGV0MWVbvX0NtqIRb9P0RC1bjSXZTvKGfh2oW4M9y3/Br+IT/N1c1U765mqGcIXdMxWUzkFOVQ8UAF+aX5WO23fr/hVjd/acfvulVI6LpO7GQMfVjn55U/x/3BW//d4p2EwgQG9w7ycObDzJs/j+CK4Nih2xNdQKl64UzFRBeZpml0X21m/4G/4PKpy4QDoweR2F12SjeVsmTL7Z+Yk0E4qNFcNcz5/fvpbGxH13VUVSV3QS7lO8qZv2L+bXtIA50DXDhygQuHLhAYCQBgsVlYuGZ0SCl7fjZ5fHLCz5X2Ozk3tuHxIeFv8JM7mMtrb71G933dWDMS+0VDCYUJ6DEd/ac6H3noI4RtGcypuG/sv8lFdHc0TWN4sI8TB19n357nGervQdc0bHYn80uXs/GBJ5k7fxF99ueNLnXWaJqGf8BP3eE6Lh69yEjfCABWu5W80jxWP7KazHmZt5xLiUai6M3rOfH2L2i6UE0o4AdFwePNpHzTDh55+FOkpWfL+w3TYHxARAI+fLVv09LUwK6hXXjvTZwFEDcjoTCBJdF/oqehCk/zbjZv2oS7cBNz568wuqyENzzYR83pAxzf/xpXmi+io6MoKvMXLGPTjg+wasNOLNbRsfA7GYpKdLqu09vWy9l3ztJU1fSe3tOiDYsou7eMJTlfvOXXyLZlcvb0AQ7veYnmy+fRtBgKCnMLFrLx3sepWL+DNG/y98Bmgw7UHnqOqL+bV/YcZuGHv4ZqMnPB/BtGl3ZXJBSuWRL9p7H/vyxnGQAHX/gzluQGyZ83n+J1H8XuTPynAKMMD/bx8o/+iXOVh4mEgyiKSnpmLpt3fJC1Wx7C4538W7qzuaJkts3RPk44GKCtoY5Dr/+czrZGouEQiqqSNbeAVVsfZOf2D+Nye2+6CELX9dHe2KE3OLr3Fwz0daFpMaxWO8tXb+XJj/8GHgmGu9Zcd4hI73ne3vMOJfd9jTnzy6jtrn3PxyRiQKR0KEwUBONFIyHe/s8v8OC96zA7Mli25ZOopsRbDhkPzhzdzY//4xtEoxEcTg+rNuzkngc+zJz84mlfXz+ZpbXXzUaATHW9/vUhyoB/hMrjezi85yU62hpHe1YolCxbzcZ7H6d8zTYslpuPX+u6Tk9nGwd2P0/lsT34RgYxm608+4WvsHrj/Xf1M6W6vq4mei/u4sKFOnrUpax94Avv+5hEDYjJhkLS3AnHBwFMHAbXmS021n3wjzn46v/mgft3UH/qFUo3fGiWXmlLLr3d7USjEVRFZeO9j/PQk5/BZp+ZE6mmNO8zMvMvWN3pPJTD6WbjvU+wcMlKTh16i2MHXsU/PMil2tNcab5IbdVR7n3oI8wtWIDZ/P7JaEVRyJlbyOPP/BpWq519b/yEaDRMX3f7BN9NTFbQP0hX3S56uju4cCXEg5/9/IQfN/7eUttdO3bvSZRwuJ2ED4Xb9QpuJiuvhHnrv8CxY8+xdesWms7uZUH5jhmoMHWkpWfNWCBMVbwvHFBVlTl5RTz69BdYvnorh995mZpT+wn4Rzh9dBcXzh5n044PsG7rI2Tnzpuw12WzO3B50me/+CQUCQe4fOJ5IsERjpyqY+fn/mVSPd333HO6370XJXJAJGQoTKVXcCslK+/nZOdlKisrWblSp+2inYLSTdNRYspweUbnY0xmS9wEQiJRFIWikjJy8+ZTtnIz+9/8GW3NF/GNDLL39R9zvvII2x58mtUb7x+brB/P4XBhMlvQIiEJiDsUjYa5ePRnKFqAvQeOsuHD38Bmn/o7CdfvQ4nee0ioULjTXsGtrHvoVzn04jepq6tjKdCmqBQs3jAtXzsVWK4dJqCo6oRDHWJyHE43K9fvYHHZGo7ufYWTh9+ku72F9rbLvPXyd5mTX0xRSdn7Pk81mVGuTeSbbzEPISYWjYa5cOSnmHQfe/YdpOzhr5E5p/iuvmaih0NChMJMhMF4Wz70e+z7yR+hXLjAEqAlFmH+0q3T/n1S0cjwAP/9o38mHA6yeNkaNu34ICaDzkZIBE5XGtsf+RhlqzZzeM9LnDm2h7KVm8ifv+i2nxtvG+fFu0g4wIVjP8Ws+Xln736Ktv4mecXTt0Q9UYeW4joUZjoMrlMUhe0f+z8c+PmfED13jvJyhUtnhihZ9YhcaHcpEg5RW32EgH8Eh9ODrmvM5FGaycBkMjF33gKe/PhvsnrTA2TnzrvlaiQxdf7hPprOvIwa8/P23gMsvPe3KShdP2PfL5F6D3H5GuSS6D+N/eKW5Syb0UC4TlEUtj3zh3TF5nPixAnwt1F75MdEwsEZ/97JTtdhEiufxQ1MZjMLFq+Qdw+mWW9HA62VzxP09fHm7n0suv8rMxoI442/n42/z8WTuOopzFbP4GYURWHTE/+Ts0ee5529L7N1yybqj/yA/LKHSc8pnPV6hBDTR9d1ms7tJTZ4ia7Odo6eucjmj/wN3qx5s15LPPcc4iIUjA6DG5Vvfpq27CJe3/VXbN24CsuFN+jvXExR2TZUVYY+hEg0/uFemqpex6oEqKmpoanXzH2f/ResdmN36r0xHOIhGAwNhXgLg/EKFq8jK+9fOPCzP2TR3E7KynRqDzUyb9mDpGcXGF2eEGISNF2jtfYgwd4LhHzD7D92EnfJY9z/2Cfiar5w7P7XbXyvwbBQGD9nEK8c7nQe/Ow/UHPwJ7z2+gts3rSO7rrX6LRmU1R+v+yZJESc0oHulvP0NB3Foka5VH+RusY+1jzx/5Kdf/uVXEZZlrPM8CGlWQ+FRAiD8RRFoWLbs4ysuI8jr/wFmZYLrFy5kuZTP0F1zaNw6T0SDkLECR3obrtAb/MJrEqQ/p4OTpyqImvpB3jgC59MiC3GjZ5vmLVQiOehoslwp+dy36f+itb6k+ze+y8U5zkpLQ3RcqoN7LnkFK8lPWe+7J8khAGi0TAdjVUMdZzDqkYY6u3l9JlKSF/Bll/6NxyudKNLnDKj5htmJRQSrXdwK4WL11Gw6NtcrtrDrgM/Jj/TTGlpKWq4h44LFuYsupfMuQuNLlOIlKDrOo01uwkPNKIqOv3d3ZyprMaUsZwVH/xLvEkw/zfbQ0ozGgrJFAbjKYrColX3s2jV/bRcOMbBU8+j+I6yelUFZvagqCYycouMLlOIpKYD9adfQQ12cKHuPA0tnXgKN7PuY/+MKy3b6PKm1Wz2GmYsFJI1EG40f8lG5i/ZiH+kn/0/+n02VWho+luYzI+TlplvdHlCJK2GqrdQgx0cO3aMWPa93PcriTFncDdmY5XStP8Gr7+lN1tvIscLpzuD+37pbzla2UB/Xw9Xz75GKDBidFlCJKW2+mPoI82cPn2aWPZWVu/8paQPhPHGvxU93aa1p5AqvYObsdpdbP/U33Lwh7/BYw/cQ0djJUVl9xhd1h3Z9+L7t2meSONlC8P9KmaLQs0RC76u936eb8TKQLdCOKzSVGti/0s2TKbb76a6/anQHdWdaCb7e76Zy5csDPYqRCMqNYctjHRM7esl6u95uOM8V1ub6VcWsmXnZ40uxxAz1WuYllBI9JVF08nhSie9aCt9fX1kWluMLmdCk7kReZ2Te3PbaVMxqWBWFVw29X2fp8ZMmE0KMRXsltH/PpldUu/2ZpkoJvt7vhmXTcWsKugqE/7+b+d2v+d4DA3fYDdmJUJtXR0bPvFto8sx3PiJ6OkIhrsOhVTvHUykZOUjXHj7j9mSlUU45MdqM+7wmZtd9Hd7M5pp8V5fsrjd7/lm7cfIsOhuO080GiVmyUnIpaYzYTpXKN1VKEggTCxz7gJO+0b/f1frOQoWzc4OjPD+i1huruJuTNR+Bv2x97Wz2QyJ4GALvb29eAvWztr3TATTtULpjkJBwuD2XHkr6e/vRzdfhBkMBQkBMdtubGM3hsRMBoRvsBuTFuDcufOUP/n5Gfs+iexuh5OmHAoSCJOzdOPT1Lz+B2zfnol/uA+nZ/r2xJcgEPHkxvY3kwHR0XSaUChEUM0iLTNvWr92Mrmb4aQphcLi6F+BRQJhMjJyi/DpaUQiEdobT1NS8cBdf83xF5sEgYhX19vmdPcgdF0nNNBCW1sbc5be/fWU7G4cTjrPFyf1eVNe2CuBMHm5pfdz9epVwoOtd/w19r1oG/sHRi84CQSRCK631evtdXw7vhN9HQ1YzVB3sYElax6ZrjKT3vV79uLoX03q46fUUyjNLp16RSmsuGwb515+kfnz5xOLRTGZJv/rll6BSCYT9R6m2nMYGeggEolgds/DbEmNJcvTZVnOMoL+oUl9bOq8AmgAtzeXEX8IRVHwDXZP+vOkVyCS1d30HCKBAQKBABb33JkqTxAnx3EmK0VRMNkyURSF4YHO206MjQ8DIZLZnfQcYqEhfD4fzsySGa8vlUkozDCTK5tAIIDJ13vTj5GhIpGqphIOWtRHd3c3WUsenbX6UpEMH80wizOLUChELDQ84X+XoSIh3j+sdKNYLIqJGD09PXizZPfhmSShMIN0XWekvYq0tDSszve+p3B9PFXCQIh33SwYVJOZqG5h8eLFXDrzphGlpQwJhRnUXHeEwjkedBTyS959q1nmDoS4uesPSuMnohXAlV1Kfn4+vZfeRtd1Y4tMYhIKM6jx1PMsXrwYk3MuFpvjfe8bCCFu7sZeQ/6itaCoFOVn0HjugJGlJTUJhRkSGBnAFGjF7XaTs2CtzB0IcQfGB4PF6kB1zqWkpISW6lcNrix5SSjMkAunXqGsbBkRzUzV/tEldBIGQkzd+GDIKqzA7XYTG2ogHPIbXFlyklCYIX2Nh5k7dy4dHYsACQQh7sb1Hnb1wSWM+FTKli6mvnKX0WUlJQmFGTA80IlN68dsNoOySgJBiGmS7jIxNFREXl4e3Rf3GV1OUpJQmAGXK3exbNlShoYdZHjmGF2OEEnFbluFw+EgPNBCNBo2upykI6EwAyp315CZmUk4VGB0KUIkHYulgEAASkuL+KffqTK6nKQjoTDN/u/v28nNbsDhcKCqS4wuR4ikoygqweAc8vLycFgP8K2vuIwuKalIKEwTXdf5+9+qQQ//JosX5REKaVitC4wuS4ikFI0uIC0tjfnz3iEc+C+C/kGjS0oasiHeXQoHfZw7/DPqD7/N9s1WlixZQmZmJv39JVgsMsEsxEyw2VYQDB7ngx98nKtXazj03Z9hSitiwdqnmb9kI4qiGF1iwpJQuEMjg91UvfOfhDpOsaJ8KRs+sRUw4fPl0de3BotlkdElCpG0VNVJNPrL+HwHycx08MTj8xkZGeFi3ffYfeBfyV/xQZZteBJVlcGQqZJQmKLAyACnd/0r0Z4zrFu7iuw1DzE8otLXV47dvgGTyYVJOghCzDhVdeFwPIyu6/T2thCJHWNFhYsKYly6dIi3v/0CBWs+xtJ1j0vPYQokFCZJ0zSq9z9H34VfsG5NBXPWPkBEs9DcthmTpZx0l8XoEoVISYqiYLMVYaOIjp5BihbvY9FiEyUlJVy4sJvd//EyFQ//NnMK5Xz5yZBQmISutjoqX/smSxdks+mRB4hoZtwF66k7tQ6zVd5WFiJeZLi9tF5+nC2PDdBSu5+yMjMLFvg5uffPqHeVsekDvy3nO9+GDLjdgqZpnHzzX7m0+094YNsqFi0uxZq5lKX3fIa5RStQFAkEIeLR4dfSKVn1MAWrnsHqnsuOHTtYkhvg7f/4Au1N1UaXF9ekp3ATgZEBDvzkqywr8lB6307CupOiNU9gd3qBiU+HEkIYz+s0jR3vuf2pTJZteZaOlhryoxq5ubkc2feXdDY9xMrtn5K5hglIT2ECHc1nOfCDL7J1dRGLF5diyypn2T2ffF8gSC9BzLRIJIimaUaXkXBuvDbnzl/BgvWfQHVks3PndlwjJ9j73FeJRiY+DzqVSSjcoOHsPi69/ec8fN9W0tKzyVn2KIVLN3Pj84QEgphpPT0XOXbsn2hrO2p0KQlrfI/e5nCzbMvHMXkXU1GxkvIiK2//16/Li283kOGjcc4dfRFf/Uvct/NeooqTRRuewWpzvudjZNhoapzOLB577B/Q9RgORyaqKk3uVnRdZ2Skk8uXd3H+/Av4fJ309l4kO3sJTmeW0eUllOvDSOMpwILyHXS35VEIOBwO9n7319n67F/jyZhrSJ3xRq7Qa84deYFw86ts376NiOpl2YYPo5om/vVIL2HyTCYLubllRpeREMJhP83N+6mp+TE9PRfQdQ1FUVBVE4FAv4TCHRqdW3jvMFFOwRKsDjcx7XV2bl3N3p/8Dlue/Rs86bKrsYQCUHv8ZYKNv+Cee7aiWXNYuu7JCSegpJfwfoFALwCxWIRQaNjgahJTNBqir+8Sp09/h/b2M0QiPkDB48mntPQxlix5Apdr4ptVODxCLBYBIBDom8WqE8NEvYWx/5Y1D/Oqp2g58yI7tqzmnR/9Nts++fe40rJnucr4kvKh0HT+IMMXXuDee+8has5iyU0C4TrpJbzX9RuSrmtoWtTgahKLrusMDrZw7tzPuXTpLUKh0bFtq9VDcfF2Vqz4GFlZi2/5NWKxCLquXfv/crbAzUzUWwBwebKYv+pDtFS+xI4tq3jnud/hgV/+FyxWhwFVxoeUDoXuq/U0H/pnHnxgBxE1jWUbnrppIEgvYWKqar72O9O5dOkN5sxZQXb2Usxmq9GlxS1d1wmFhqivf526ul8wMNCIrmuYTDays0tZs+bzzJlTgdXqvOnXiMUi9PTU0dCwG9CvDTOl9OV8U7fqLQC40rIpqPgg7Wf/m02rSnjnB7/HA5/9+5TdNyllW1EoMMzpl7/OQzs3E8XO0o1P33bNsvQS3i8/fw1OZw4+Xxe9vfW8/vr/pKzsKZYt+zBpaflGlxd3otEQTU37OXfuZ3R21gA6oJCRUcKKFc+ycOF9WK23Ph9gaOgqdXUvc/7884TDI+i6jts9h7y8NbPyMySqm/UWADzpuYRLdqKwhzKfjxOv/yMbH//yLFcYH1IyFHRd58BPv86mtWVYbC4Wrn8Gk+nmexdJL+HmsrJKWb/+f3Dq1L/j83USifioqfkJV66cZPnyj1BUtA27Pc3oMg0Xi0Xp66unpubHtLYeIRQaAhTs9nSWLn2SxYsfxeudf8un01BohNbWI9TUPEdv7yU0LQIouN1zWL36s7cdakplt+stAGTllRDw9VGqKnQfOkzjuX0sWL59liqMHykZCtX7n2NBjsbcuXPJXHQ/Nof7tp8jvYSJmUwWlix5nPz8NZw//zxNTQcYHGymp6eO/fv/jHnz1rFixbPk5a1NySGl6/MG9fVvcP7889fCACwWJwsW7KSs7MNkZy+7ZRjEYhE6Oio5e/antLYeQdOi6LpOWloBCxZsp6zsGemVTdKtegsABYvWc6GvjY0bN/D6m/+XOfNX4PRkzmKFxku5UBjsaWPg4itsfOg+TN7FZM6R09Gmg8eTx/r1v05JyUPU1r7I5cu7CYeHaWs7Rm9vPUVF97Jy5SfxeOalzFhtKDRMU9M+amp+wsBAI5oWRVUt5OQsZcWKT1BQsOGWQ0WapjEy0kF19Q9patqP398D6FitHhYs2MGyZU+Rnb0EVZUHlsmYTG8BYNHaJ7hw6Pts2rCaoy//Gfd96q9mobr4kVKhoOs6J/77z7l30zoiuo1Fk+ga7nvRJr2ESVJVlezsUrZu/V1KSh64Nox0gkCgj7q6l2ltPcyyZU+xePGjeDzJ+6JQJBLg6tVT1NT8iI6OqmurshQyMxdTWvoYS5c+ectJZACfr5v6+jeorX2R4eF2QMdkspGfv4blyz9CQcGmlAnX2WYyWchb9gCWC2+S72nlwqnXWLL2MaPLmjUpFQoXT79BUa6JtDQvc8sflc2wZoiqquTnryUzcxHNzQeorn6OgYFmfL4uTp/+D9rajlFW9jTFxduTakhJ0zQGBpqorn6OlpZDBIP9gILN5qWk5EGWL3/mtvMG0WiY1tYjnDv3Uzo6qtG0CIqi4vEUUFHxCYqLt+N0ptZwxnTyOk23HUICyMiZT29bIRUVMV5543ssXHEfFqt9lqo0VsqEQjQSovn493j84e2orvl40nNv+zkywXx37HYvpaWPU1i4mdraF7l48TWGh9vp6Kikq+sc8+dvpazsw+TlrcF0k7fHE8XgYBuXLr3B+fMvXnuJTMdicTFv3noqKj5Bbu7yWw7zaFqMjo4qzp9/gebm/dfeOVBwu+eyePGjLFv2FC5XjjzIzKLiFfdRf+i7rF5Ryuld/8bGx79kdEmzIrGvxCmofOd7rFxegq5aKCrfMenPk6Gju6MoCk5nFqtXf46Cgo2cP/8iDQ27icXCNDXtpavrHIsWPURZ2dO43XMTbkhkdGuKA5w79zO6u2vR9RiKopKRsYiVKz9FYeHmW66+0jQNv7+H2toXqK9/k5GRdgBU1UJx8XaWL3+G3NzyhA/NeDOZ3oLZbCWtYC0WM5x7czf+4U+lxKRzSrS0SDhIf8PbbH3kPuy5y+XkJQOoqok5c1aQlbWEkpIHOHfuZ1y5chK/v5vq6udoaHibiopPsmDBDlyuHKPLva1oNERHRxXV1c9x9erJsXkDr7eIJUueYMmSJ7Db02/5ZO/399HUtI+qqu+PzRuoqoW8vFWUl3+MefM2JNXwWryY7IQzQH7JGmqvVrNm1QrOvP3vbP3Q781wdcZLiVCoOfAjKsoWEdVNFCxaP6nPkaGjmWE2W5k/fwvZ2UtoaNjDuXM/Y3CwlZGRDo4d+ycaG9+houITcXtDfHdF0HM0Ne29tiIIzGYHixY9xNKlT5KdvfS2S0zb289QVfUDOjtriEYDjO51NI+ysg9TUvIgbvfthzfFzFOAzMK1WNUIp87swj/Sj9OdYXRZMyrpQ0HTNHovvc3mh+/Flr3spjufTkSGjmaO05k1Ntl89uyPaWh4h5GRdtrbT9PTU0dh4WYqKj5BVlbpLV8snC26ruPzdV+bN3iBkZFOdF3DYnGSn7+GiopPMmfOilvWGotF6Ou7RE3NT2huPnBt4ztwuXKvLdn9FC5XbsINoSW7OcUr6Gs9yaqK5Zw7+GPWP/JFo0uaUUkfCo01e1lQkIWOyrzF64wuR4yjqipudy7r13+RoqJ7OXv2p7S0HCIS8dPQ8Dbd3XWUlDxAeflHcTgyDZtkjcWiNDXto7b2RTo6KseGijIyFrBixbMUF2/H4bj506Ou6wSDA2Mb3w0NtQJgMlmZN2/DtYno8rjsGSWrya5CgtHegmdOGWbCnKh8B037taQO7qQPhZbqX7BjfQmqKw+zWYaE4pHJNDqOnp29lPb2U1RWfo/u7jqGh69QVfUDmpr2U17+0VlfjhmLRejsrOHcuZ/T2npkbJjH7c5j8eKHKSt7Gocj65Y3iEBggJaWg9TU/IT+/gZ0PYaqWsjOXkJFxScpKNh423cWhPHyFq7mUkcVi4vzaKjZy6KV9xld0oxJ6lAI+gdRg1dxOMqYUzz5zcJkPsEYFoud+fO3kpm5mIaGtzl79ieMjHQwMNDI4cN/Q0vLIZYvf5r8/PUzuhpn9PSzDmprX+Lixdfw+7uB0RVBCxfeT1nZhye1xPTq1VOcP/8Cra2Hx7a1drlyKS//KAsXPpDUL/AlG7PZismVx/z5Pg5WvSWhkKgunnyNsmVLiOoWPJlT2xtG5hOM43bnUl7+MebP38r588/T0LAHv7+blpZDtLefHrsxZ2YunvZw8Pm6aWk5TGXl9xgZ6UDXY5hMNubMWXFtRdA6LJab77WvaTH6+i5TV/ff1Ne/TiTiB3QcjkyKi3dQXv4RvN6ipB5+SFaZ85ajBDuIDjYQCQeT9mW2pA6FvuajrNpUis1biLzyk1hUVSU9fT6bNn1p7OW36/MNFy68wpUrJ1m27EMsWfLEtBxTOboi6DSVld+/tjXF6OFBaWkFlJV9mMWLH73lvAFAMDjIhQu/oLb2ZYaG2hg958BEQcEGysqeprBws+xTFEemMq8AkDlnAV0XoHRRMY1n91O65qEZrtAYSRsK0WiY2MgVbLYVZOYvNboccYdU1URBwUbmzKmguXn0HIKenouMjLRz6tS/c/nyblat+iUKCjZit3vv6HsMDLRw6tS3aW09Sjg8eqSo2z332ulnz15bEXTzm3koNLrxX1XV9+nra0DTIqiqmczMRSxf/hEWLNiBxeKUt5ETnKIoWNy5zJ07xLELRyUUEk3bxRMUFxUQiel4s+YZXY64C4qiYLU6Wbz4EfLy1nDhwi+oq3sZn6+Lvr569u37UwoLN1Ne/hHy8tZM+ebb2VlNY+M7aFoUk8lGYeEmVqz4OHPmrLjtk31HR/XYqqnRiWhwOLJYuvQDLFnyQdnSOsk4M4vR/B0Ee08bXcqMSdpQ6Gw8xYr8PEz2rCndJGRX1PjmdueyatVnWLjwfqqrf0hz8wGCwQGamvbR0VE1Nt+Qnl486XH7BQvuG9vqu6LiU5Pa0npwsOXaFuG7CAT6AR2bLY358+9h5cpP4fUWxsX7FWJ65eSXMtx6BIc5im+oF1fa3Q9dxpukDYVAXwOeJaVY3fG/ZYKYGpPJTEZGMVu3/i4LFuzg3Lmf09Z2jGCw/9r2GSdYtuxJFi9+DIcj/bZfz2p1cs89v4vZbL/t/EQwOMSlS29SV/cyfX2Xrv1bhfz8dZSXf5SCgk3yvkESs9gcaIqV4uJirlw+TenqB40uadolbShEfZ2YzWV4MguMLkXMkNEtM7aSk7Oc5uYDnD37EwYGmhgcbObEiX+hoWEPq1Z9mry8Ndhstz5dLy3t1kOM4bCPjo4qKiu/S3d3HbFYCEUx4fXOZ/nyZ1iwYKdsaZ2gpjLZDGCypZOVlUVr61mQUEgMQf8QVlMUgLTMPIOrETPN4Uhn6dIPkJ+/lrq6l7l48VX8/h46O2vYs+frFBXdQ3n5R8nJWT7lpaCaptHbe4Gamp/Q1LSPSMSPoig4HJljW1p7vfLgkaimsjnedRZXFg7HFUJDLTNUlbGSMhT6OhrIzc1F0xUs1puvKRfJJS0tnzVrPs+CBTuoqvoBbW3HCYeHuXx5N+3tlSxa9CDLl38El2vObcNh/JbWFy++js/Xxei8gYe8vNWsWvUZsrJKZagoBdldGUSsVmLBPqNLmRHJGQqdTeTl5IBZtg9INWazlZycZezY8XXa2o5SXf0cXV1nr23R/SNaWo6MvXdwsyGlcNg/dmDO6LzB6PsGublllJc/S1HRNgmDFOZKy2EY0EKDRpcyI5IyFEK+XmwZNlRTcr5xKG7PbLZSVLSN3NzlNDbupbr6ubEtM44f/2eamvZRUfFx5s5dNbbSKBz209V1lurq5+joqCIaDaAoKi7XXMrLP0pJyQM4ndnyvkGKsznT0HUdlSixaASTOblWmSVlKESCQ1gsFtQpboAny1GTy/VT38rKPkxBwYaxvYyCwX6uXj1JZ2cNhYWbKS19HEVRqK9/g5aWg0SjQQBsNi+LFz/CsmVPkZ5eJGEgALDanOi6jtPpxD/Shyd9jtElTaukDIVYaASz2TPlUBDJSVEUvN5C1q37NQoLN1NT82OuXj1FNBqgqWkv7e1nUBSFYHAAALPZTn7+OlaseJY5c1bI7rriPRRFQVdMuN1uAiP9EgqJYPScXAUU2XRMvMtstjJv3jrmzKng8uW3qKr6IQMDjWNhMBoeRaxY8SylpY/LvIG4OUXFbDYTi0aMrmTaJWUooGsoioIim4+JCZjNVhYvfoycnLJrW1QcBKCgYBMVFR8nPX2B7GIqbkNBVVU0bWrLWRNBcobCNbquG12CiFOqqpKZuZBt236fWOx/AqMnocm8gZgsRVEgCe8xSRkKqslGLBZD15Kvayeml6IoMmcgpk6PEgqFSLfffI+sRJWUfWTFbEPTNPRo2OhShBBJRgcUXSMYDGKVUEgMFpubSCQiPQUhxG1dP2xnsqLhAKqq4vP5cLhvffBSIkrKULB78wgGg8QifqNLEUIkmcDIAACRmILVlny7JiRlKHgz8+nv70ePTX7nQyGEmIyAbwAA5RZnbiSypAyFrLwSurq6UPUI0YgEgxBi+gSHu4lEIpjtyXfADiRpKDg9mfjDOoqiMNzfYXQ5Qog4N5XzFML+Pvx+P7Yk3TI9KUMBwOzIRtM0hvuuGl2KECKOTfU8hWhwgKGhITw5JTNUkbGSNhRsGQsYHh4mONRudClCiCQRCQdQtSANDQ0ULl5vdDkzImlDYc6CdXR0dBAL9pF87xwKIYzQ19mEoiiMhJSkPdUxaUOhYPF6GhubMKsaIwNdRpcjhEgCwz0N+P1+rN4io0uZMUkbClabExx5RCIRetrOG12OECLB6UBkpIPOzk5yFm42upwZk7ShAOAtWEdPTw/BwVajSxFCJLjhvqtYVI0LFy9RsmKn0eXMmKQOhcVrHuH8+VrMeoCRwW6jyxFCJLCulioCgQA4C5Jyz6PrkjoUPBlzCZmyCYVCdDaevu3Hb38qNOXlaUKI5KdpMSJDbbS0tJC37EGjy5lRSR0KAHOXPkBrayvhoRa0WNTocoQQCaiz5SxmVaf+cguLV0koJLTStY9SW9+MWdW5evn2vQUhhBhPBwbaKunp6cE+dx0ms8XokmZU0oeC2WzFW3Qv7e3tDHWeldPYhBBT0td+GYsS5vSZKlbu+IzR5cy4pA8FgIrtn+R05VmsaoyrDZVGlyOESBA60HX5MP39/eBZjDs91+iSZlxKhILd6cVVsIUrV64weOV0Uh62LYSYukF/7Jab4XW3nMeiBDl56jQr7/+1WazMOCkRCgBrHvgVzlSdx6xEaa09aHQ5Qog4p2kxepqO0tXVhe4tx5udnLui3ihlQsFqd5G97APU19cT6rs4dnrSRGRZqhCi+fwBLGqUE6eqWfvwrxtdzqxJmVAAWLHtE9Q19hEI+Giqfn3CjfKmsq+6ECI5+YZ6iQzUU1tbS075h3Em4VnMN5NSoaCqKisf/V8cP34CKz6uXDxmdElCiDijaTGaql5hZHiI+tZhyrd8xOiSZlVKhQLAnMJlmOZs5eLFi/i7ahgelB1UhRDvaqx5GwtBDh0+xoan/ghFUYwuaValXCgArH3o1zjX0MvQYD8tVb8gEg4aXZIQYpZNtPKos+Uc+kgzVVVV5FY8S3p2oUHVGSclQ0FVVbZ85E85cPgUSixI/YkX5aU2IVLc8EAHgy1HaGpqoiOQSdmmp4wuyRApGQoAnvQ5LH3gdzl06AhWRqg/9crYxLNsjCdEagn4Bmit+gUD/b2cqb3Ktmf+yOiSDJOyoQBQsGgNaUuf5tixY6ihDhqqdxldkhBiloVDfhpOPk8oMMz+I5Vs+/hfJP3+RreS0qEAULbxQ0SztnDmzBkYaeZy1S4501mIJHd9JCAUGKH+2I/Ro372HjjOxme+gSst2+DqjGU2uoB4sOa+z3H8dR+nT59m9Wqdy6dfQ9eTfzwxGu0FXsNq9REMZhONLsRuL0VV3UaXJmaJrutEo+1Eo3VYLG04HCFGRpZgtW5L+lU3Gx7q5NLx54mFR9iz/wgVT/wx6TnzjS7LcBIK12x49Dc4tevbnDhxgvXrwcRzDPg+QrrLbnRpMyIUuoTH8xrh8DBdXT1kZnaTnt6Kru8jGLQTDmcQieSiqvOwWPJQVU/S3ySSna5HiUZ7iUTaUNV2rNZerNZB3G6NaDTK4OAg7e0jFBf76eu7isXyYRTFanTZMyIababp9CuEA0O8s/84az70/5GVV2J0WXFBQmGctQ/+CjUHPezf/zpbt25hYOg7RKPPYjYn19uMgcARsrKO09nZzq5d9QwMPITD0YbHc5jc3DDFxcVkZGSQnt6K2VwFQCwGkYidSMRJNLoMm22DwT+FuB1N8xGNvoXV2ofJFMBiCWMyqei6TiAQwOfz0dTUQVPTVQYH5+HzrSAS8VBW9g733nsPgcB/EIs9g9mcY/SPMq0CgZMUzDlMb08nh06cY9Mz35QewjgSCjdYcc+zNGXM481d/8y2rRux2X7A4OB2HI4Ko0u7a5oWJBR6gdzcHi5erGfv3ggOx3Pk5IyeN6vr0NbWyYULJzCZzmG3N+BwXMXp9JGZmUl2djZer5c5cwbp6mrGbn8KRZEmFI8ikatYLM/j8YTo6uqir6+Pnp4ehoaCBAKZ+P3zCYUWYzY/hdO5GrPZhtc7+rkNDe8wOPhNHnhgIy7XcwwObsVuX2fsDzQNdD1MMPgyOTlXaWxspKq+lx2f+b/YnV6jS4srij6JBfpDQ0N4vV6+9l8t2J1ps1GX4fo6GnnzX/9fHrh/PvPnz6enJx+r9QlU1WF0aXckGKzHbn8Tmy3IsWPHOXduAx7P76Eot19roGlhwuEWwuFGNO0Qq1adY+PGDQwNOVGUj2MypUabSBSBQCWZmfsZGupj167TdHd/HKt1MTbbAszm3EkNA4bDLej677B9+xwWLlxId3cOVusHUdXEPLA+FGrAan0dpzPE8eMnGLEuYvMH/xcmU+o81AT9Q/zpZ+czODhIWtrNr1kJhVuIRkJ8+6t/w9aNDaxZs4ZIxIrPtw27fVXCjK9rmo9Q6DVycq7S1dXFvn01DAx8Fbd72x1/zZGR18jP/2d27twEuAgEnsRmK5q+osUd0XWdYPANcnIu0tTUxJ493ZjN/4jFcmfDP7oeYWjoz1mxoor169ehaTaGhzdjt69LoPYfIBx+naysFvr6+ti37wwXm77Ml/9us9GlzbrJhkLqxOQdMFtsYP7/qK4+QGPjX7NtWykFBQfo7z+Jpj2M1Rq/45C6HiUY3E9aWg12e5hjx85QVTUPh+M53O70u/rabvdjdHWV8OKL/4uHHlqO1/syfv8nMJtTeymf0YLBt8nNraeysoojR/JIS/s+inLn6+0VxYLX+0fU1h6mufkbbN1aTFGRzsDAaWKxB7HZFk5j9dNL12MEgwfxeCpxOCKcOlXFmTNZ2O3fx2pPB3xGlxi3pKdwG9/6iou8DDOa5md4+C9ZuPA4Gzaswev1MjCQhabdG1dPyaPjpkdwuaqx26NcvnyZ48c78Pu/jNu9c1q/Vyw2SCz2yzz1VDmQAXwBVU3O1VrxLhg8S3b2Hmpqqjl4cDXp6V+e1q+vaSGGh/+G+fP3s2nTGjIzM+nv96Jp27DZ4mfVjq5HCAaP4XJV4nBEaWxs5NixFoaH/x88nkdo74/yxW+kZiBIT2EatfdHyctw4vV+natXm3n++W9SWlrDypUrycjoY3DQQzi8Ert9pWETr6NLDY/i8VzC5dJobW3lxIl6eno+QlraZ3C7TdP+PU0mL7HYP7B796/ygQ9sY2joJ1itvzSpeQoxfSKRTrzet2lpaeHw4XS83i9N+/dQVRte71fp7v40zz//lyxaVM3q1avIyhpkaMhJKLQSu331XfVM7kY02k8kchS3+xJZWVGuXLnCiRN1dHZ+CK/37/B4UvcN5amSnsIkXO8tjBcInEPTvsWCBU2sXLmC7OxsgkHw+QrR9RXYbCUzfnPUtGFCoWqs1gt4PINEIhFaW1s5fbqF3t4Pkpb26VmZGPf5DlNS8k127ryX7u4iHI4Pzfj3FKNisSHM5u8RifTy0ku1WCw/RFVtM/59g8GLRKPforj4IqtWrSArK4twWMXnK0DXl2OzLUZRpv9BZLzR+bIqLJaLeDz9aFrsWvtvpKfnUdzuX8Zkeu/EuPQUpKcwYxyO5cA/0draTH3998jO3k1FxXzy80dwu1sJBiEYzCUcno/FsgizOeeuQ0LTgoTDzej65WtLRYdxuXT6+/s5ePA8TU1p+P0fJC3tQ6Snz96Tkcu1hYsXn8HtfpX166GnZxd2+4Oz9v1Tlab5gR8AQ7z55mlU9T9mJRAA7PZS4G+5evUKDQ3fJzNzL+XleRQUDOPxtBEMvkEwmDOu/c+ZhvYfIhJpQdMuYbNdxeUawuXSGRgY4MiRWhobHfh8H8Dj+Ru83vf/HlI5EKZCQmGSRoeQ3v/rstmKsNn+kHA4wp49b2K17iIj4xBLlhSSl9dPZuZVTKYTRCIQDrsJh9OIxbzoejqK4kZVHdee5q+v5ogSi/nRtAAwgKoOYjYPY7UOYLeHcLsVQqEQAwMDnD5dT1ubmaGhTdhsf43dvmhsrflsS0v7JU6f7sNqPcXKldDd7cDhuMeYYlKAroeJxb6Py+XjlVeO4vf/A3b73Fmvw2qdh9X6FSKR32Xfvt1YLG+Snn6YJUvmkZ+fT2ZmBybTKSIRjVAojWg0jVgsDU273v6dE7T/AJrmR9eHMJkGsFiGsVgGcDiCqKo61v5rai7T0gJDQ+uxWv8cu32pYe0/mUgoTMIXv+HjW1+59frs0ZUaTwBP4PeHOXjwMLp+CLe7Eo+nm7y8bPLz83G5XDgcNmy20ScZXddRFGXsPIfr/19RFGKxGMFgkFAoRHd3P62trfT16QwPLyAYXIPD8f9gsy2ImwshLe3LHDnyx9hsdSxapOH3r0RVPUaXlZSCwX1kZwd4/fVD9Pf/OU7nIkPrURQTXu/DwMOEQhEOHz6Krh/E7a7G7e5k7txM5s2bh8vlwum0Ybfbx9o58J72f52mae9r/729MUZGigkEVmO3/xp2+yJuMRIi7oCEwgxQVSsezw5gBwAjI1Fqauo4caISk6kVq7UTi6UXiyWIyRTEZIrAtb1ZdV0lGrURi9mJRt2EQjlEInkoyg6cztVYLLnY7WCP00U+aWlf58yZpykthVDoLA5H6q0Hnw12exMdHR00NT1BRsZKo8t5D0Wx4PFsA0bfhfH5Ypw7d5FTp85gMrVgtXZitfZiNl9v/2Hebf/KtbZvIxp1Ew7nEonkAduutf+52Gxgm+IomQwdTZ6EwiRd7y1MNIR0O4pixuEox+Eof8+/17TRfyKRiT/PbB79J5EoisLg4BoGBgawWBoACYXpFosN4XaPcPJkHS7Xbxldzm0pigmHYxkOx7L3/PtkbP/JQNYOimlnNj9EU1MTTmc3un6TK17csXD4LKDT0eHAai0wupy4194fNbqEhCKhMEXSwG7P6VxLW9tVLJbRF9zE9FLVXgKBAH5/sdGlJAwZOpo8CYUpkIY1OZFIO16vh2g0hsmUXNuOx4NYzIvNZsNq7TO6FJGEJBTEtAuF6pkzZw7RqH3GX2BKRaqag8lkwmbrNbqUuCc9+6mTUJiiL37DJw3tNiKRK3i9XiIRp9GlJCWzOQtFUbDZRowuJSFID39qJBTEtFPVPiwWC5o2O2/XphpV9aDrOlbr9beaxUTk4e3OSCjcIWlwN2cyDV4LhTh9mSLBKYqdaFTH7XYTjfYYXU5ck17C1Eko3AFpaLemqiFUVUXXZWfKmTD61rsJq9V6bTsUcSN5aLtzEgp3QRrexBQldu3GJc1rpug6mEwmdF3a4M3Iw9udkav2DkmDuzlF0a7tYZMYRzYmJkVC4SbkYe3uSCjcJWmA76dpFjRNA2JGl5K0FCVGOByelfMyEpE8tN05CYW7IA1vYprmQNM0VFW2uJgJo7uLju4gqqpuo8uJK/KQdvckFKaBNMT30jQn0WgUVQ0ZXUpS0vUwZrNKMBjEZJKtyW8kD2t3R0LhLkkDfL9YLPfa0IasjJkJsdgAAMGgLj2FceThbHpIKEwTaZDvslgKGRwcxGKRUJgJsdjonkehkPs9h9IIeUibDhIK00Aa4ntZrUV0d3djNofQdZlsnm6xWC+aphEOpxtdStyQh7LpI6EwjaRhjrJaC+nr82E2q0SjnUaXk3TM5m78fj/BYKHRpcQVeTibHhIK00Qa5LsURcXvzwYgEmkzuJrkY7H04/P5iESMPZc5XsjD2PSSUJhGsoPqu/z++YyMjGA2XzW6lKSi6zHs9kHa2tqwWlcZXU7ckIey6SOhMAMkGCAUWk1/fz9WqwwfTadwuBWLRaGtre99Z36novb+qATCNJNQmGbSQEe5XNu4fPkyDocfTZPfyXTRtHoikQjDw/NRlNS+fOXha2akdquaQaneYK3WeXR2WlAUCIXOGl1O0rDZmunr68PnW210KXFBHsKmn4TCDLjeUFM9GAYHV197X6He6FKSQiw2jMs1zIULF3A4njC6HEOl+rU1kyQUZog8wYDJ9CgNDQ243T1oWtDochJeKHQG0LlyxYLNVmR0OYa5Hghyjc0MCYUZlspPNE7nOurrh1FVnVDotNHlJDy7vZbu7m6Ghu41uhTDSSDMHAmFGZTqw0iKojAwsJ2enh5stnNGl5PQIpErpKUFqao6i9P5caPLMUyqXkuzSUJhhqX6E43D8Qmqqs6SluYnFGoyupyEpWkH8fl8XLlSiMWSa3Q5hpBho9khoTBLUvUJx2rNp6VlIUNDQ6jqQaPLSUix2CBebzt1dXXEYp8xuhxDSSDMPAmFWZDqw0jwK5w7d4709F4ikXaji0k44fA7aFqE8+djuN1bjS7HEPKS2uyRUJglqdygnc4KamvTGRkZRtd3G11OQolG+8jIaKK2tpZAIDV7Can7MGUMCYVZlNp7I/1PzpypJDOzj1CowehiEkYstotIJEhlZQS3+1Gjy5l1Mo8w+yQUDJCKweBwLOfChSJ6e3uxWN6UcxYmIRS6TFZWB1VVVYTDX07ZA3UkEGaXhMIsS+X5BZvtaxw8eAaPJ0QwuM/ocuKarsewWN6gp6eH6up83O5NRpc062QewRgSCgZI1WCwWHJpb3+Gixcvkpl5lkikw+iS4lYw+AZud5iDBytxOP7I6HJmXapdG/FEQsEgqRoMaWmf49ChEYaH+1HVF9H11Pr5JyMUqic7+xJVVVV0dn4OsznL6JJmlcwjGEtCwUCp2OgVRcVk+kv27TuJxxMkGPyF0SXFlVhsGKfzDbq6OjlxIp20tI8aXdKskkAwnoSCwVJxRZLVOo+urv/BiRMnyM1tJRg8YnRJcUHXY2jaj9D1EfbsqcXp/DOjS5pVEgjxQUIhTqRaMLjdT1JZWXFtfuE4wWBqb6+t6zrB4M/wen3s3XuUQOAvMJnSjC5r1kggxA8JhTiQqvMLHs8fsHdviI6Oq6Snv0443Gp0SYYJBl8hJ6eTI0eO0db2RRyOZUaXNGskEOKLhEKcSMVgUBQVp/OfefPNFgYGunG5XiQSuWp0WbMuEHiLnJwGTp06RU3No3g8qXeAjgRC/JBQiCOpGAyq6sRs/javvnqekZEeXK6fEw63GF3WrNB1nUDgVXJz66iurubkybV4vZ83uqxZJe8ixB8JhTiTisFgNqejqv/JL35xnoGBTtLSXiIUqjW6rBml6zGCwZ+Rk3OJU6dOceTIStLSfsfosmaVBEJ8klCIQ6kZDFmYzd/l1Vcb6OhoJTPzLQKBA0aXNSM0zU8k8h2ysq5w9OhRjh+/h7S03za6rFnT3h+VQIhjEgpxKhWDwWTyYrV+n9deC3DhQh25uWcIBn+MpoWMLm3ahELNmEz/idvdz+7d+6ip+RRe7xeNLmvWyKRy/JNQiGOpGAyqasft/hZ79y7m4MFDpKdfQVG+nfArk3RdIxDYQ0bGSwQCnbz44gFaW/8Et/tJo0ubNRIIiUFCIc6lYjAoioLX+3ucP/8bvPzyPiKRbtLTXyAQeAVdjxhd3pRFIleJxb5NTk4NdXXneemlTsLhH+B0rja6tFkjgZA4zEYXIG7v+oX0ra+4AMjLSI0/m9u9E5+vnBde+Crr17eyYoVGINCE378Vu31V3G8lrWl+QqFdZGU1MTQ0wKuvnqS19UOkpf1q3Nc+nSQQEov0FBJIKvYaLJYcXK5/5+jRp3jxxXcYHGwlJ2c/mvbvBIMX0XXd6BLfR9NC+P17sNm+TXp6PTU1VbzwQgs9Pf+G1/trEggirkkoJJhUDAaAtLSnCAR+zH//t5fXX9+Npl0hO/t1dP3b+P2VcXFoTyw2RCDwOlbrv5CdXUVb22V+9rPDHDv2LHb7f2G1Fhpd4qySQEhMij6JR62hoSG8Xi9f+68W7M7U2Y8lnqXaUNJ44XAb4fDfU1BwlvXrV5GdnY3fr+DzLcJiWYfZnDtrteh6jFCoDpPpDB5PN6DR0NDAmTNX6e9/irS0T6MoqfU3Gv/AIoEQP4L+If70s/MZHBwkLe3m9/HUaq1J5Ivf8PGtr7ho74+mXDBYrQVYrX9Jd3c7zz//r+Tk7GLlyoUUFoawWC7i89kIBIpR1SVYrUUoimlav7+m+QgGL2I21+N0tuPxwMDAAKdOXaC+HkZGPkxa2jN4van1dwHpHSSD1Gu1SSRVJ6Cvs1jysFj+mEAgxK5dL+J0vsacOV0sX76MnJxhLJY6IhGFYDCdSCQbTcvHYpmLyZSFqtpu+/V1XUfTholGu4nFrmI2d2K19uJy+XG7YWRkhEuXWrh4sZ2+vgoU5Q9wOlfh9c7CDx+HJBCSQ2rdRZJUKvcaAFTVhtf7LPAsXV09NDb+AofjCGlpTRQV5VBQUIDLdRWH493tucNhBU2zoWkWNM0MXJ/81VDVCKoaxmyOYL7269Q0DZ/PR1/fIJcvX6ary8Lw8BI07dO43dtwuVLv936dDBclF5lTSCLXewyQer2Gieh6DL+/kkjkJDbbRRyOVuz2AVwuhZycHJxOJw6HA6vVOrYiSNM0gsEggUCAwcFB+vr6CAQcBIO5BALFRCLLcLu3YbHkGfzTxQfpHSQOmVNIQeOHk1K11zCeophwudYCa8f+XTAIIyMDtLY2E4v1E432o+vDwOizka6bMZszMJnSMZtzsVrno6pWTCZwu435OeKVBEJySu27RpJK9eGk2zGb0zGb040uI2FJGCQ3uWMkqVSfhBYzQwIh+cmdIsmN7zWAhIO4MxIGqUPeaE4BX/yGL2XfhBZ35/rZByCBkCrksTGFyJCSmAoJg9Qkd4UUJENK4lYkDFKb3A1S1I3LV0HCIdXJS2gCJBRSnoSDkDAQ48nVLwAJh1QkYSAmIle9eA8Jh+QnYSBuRa52MSEJh+QiQSAmS65ycUsThQNIQCQKCQMxVXJli0kZf0OR3kN8u/EFRQkDMRVyRYspk95DfJJegZgOchWLO3az3gNIQMwWCQIx3eTKFdNCAmL2yPCQmElytYppJwExvSQExGySK1TMqBtvYOOPDAUJiYlMtJOtBIGYLXJFill1u5CA1AsKCQERT1Lr6hNxZ6KQmOgmmQxBcbOzLCQARDxJ/CtNJJWb3SAn6lFcF0+BcbtDjCQARLyLn6tJiFu41c30VoFxK7cKk7s5oU5u/CKRSSiIhHc3N+GbBYrc2EWqklAQKU1u/kK8l2p0AUIIIeKHhIIQQogxEgpCCCHGSCgIIYQYI6EghBBijISCEEKIMRIKQgghxkgoCCGEGCOhIIQQYoyEghBCiDESCkIIIcZIKAghhBgjoSCEEGLMpHZJ1XUdgFBgeEaLEUIIMTOu37+v389vRtFv9xFAW1sbhYWF01OZEEIIw7S2tlJQUHDT/z6pUNA0jatXr+LxeFAUZVoLFEIIMfN0XWd4eJj8/HxU9eYzB5MKBSGEEKlBJpqFEEKMkVAQQggxRkJBCCHEGAkFIYQQYyQUhBBCjJFQEEIIMUZCQQghxJj/H8ISac6bkHD+AAAAAElFTkSuQmCC"]},"metadata":{}}],"source":["# hints from:\n","# - plotting molecules in matplotlib, and matching coordinates of points in plot-space:\n","# https://ljmartin.github.io/blog/16_2d_molpics.html\n","\n","from rdkit.Chem import Draw, rdDepictor\n","from rdkit import Geometry\n","import PIL\n","import matplotlib as mpl\n","import io\n","\n","##\n","# Step 1: draw the molecule into a PNG. We'll later convert all the coordinates into PNG-coordinates.\n","## \n","\n","#draw to a canvas in PNG format:\n","d = Draw.MolDraw2DCairo(350, 200)\n","d.SetFontSize(0.8)\n","opt = d.drawOptions()\n","opt.bondLineWidth=4\n","opt.baseFontSize=0.8\n","opt.clearBackground=False\n","d.DrawMolecule(Chem.RemoveHs(molH))\n","d.FinishDrawing()\n","txt = d.GetDrawingText()\n","\n","#convert PNG to PIL object so we can plot:\n","image_stream = io.BytesIO(txt)\n","image = PIL.Image.open(image_stream)\n","\n","#plot the actual molecule, with a transparent background. \n","fig, ax = plt.subplots()\n","\n","# plot the actual RDKit molecule. \n","ax.imshow(image,zorder=10,alpha=0.8)\n","\n","##\n","# Plot the electrostatic potential. We'll need to convert the coordinates of \n","# the grid into PNG-space. This is all in https://ljmartin.github.io/blog/16_2d_molpics.html\n","# then we just use 'contourf', after defining a colourmap that mostly reproduces\n","# the source colours. \n","## \n","\n","#figure out mapping to canvas coordinates:\n","pt = d.GetDrawCoords(Geometry.Point2D(0,0))\n","interceptx, intercepty = pt.x, pt.y\n","rise, run = 1, 1\n","pt = d.GetDrawCoords(Geometry.Point2D(run,rise))\n","gradx = (pt.x - interceptx)/run\n","grady = (pt.y - intercepty)/rise\n","\n","canvas_coords = pts[:,:2] * np.array([gradx, grady]) + np.array([interceptx, intercepty])\n","xs = canvas_coords[:,0].reshape(grids[0].shape)\n","ys = canvas_coords[:,1].reshape(grids[0].shape) #both grids are same shape\n","\n","colors = [\"red\", \"yellow\", \"green\", \"cornflowerblue\", \"blue\"]\n","cmap= mpl.colors.ListedColormap(colors)\n","#unchanged ESP\n","#boundaries=[esp.min(), -0.004, -0.002, 0.0, 0.002, 0.004, esp.max()]\n","#converted ESP\n","boundaries=[esp.min(), -4.5, -1.5, 1.5, 4.5, esp.max()]\n","ax.contourf(xs, ys, esp,zorder=0, cmap=cmap,alpha=0.8, \n"," levels=boundaries)\n","\n","##\n","# Bonus! Plot a VdW surface. Most of this is from https://ljmartin.github.io/blog/16_2d_molpics.html,\n","# however this time I generate a concave hull so it's a smooth line. \n","# Matplotlib conveniently allows us to shade the inside of the VdW surface. \n","##\n","\n","#this function returns points around a circle, like an atom:\n","def circle(centre, radius):\n"," t = np.linspace(0, 2*np.pi, 100)\n"," x,y = np.sin(t), np.cos(t)\n"," coords = np.array([x,y]).T\n"," return coords*radius + centre\n","\n","#fetch atom vDW radii:\n","ptable = Chem.GetPeriodicTable()\n","radii = np.array([ptable.GetRvdw(at.GetAtomicNum()) for at in molH.GetAtoms()])\n","\n","#for each atom, get an array of points on the circumference. \n","#then remove any points that sit within the radius of any other atom\n","#the remaining points sit on the vdw surface. \n","vdw_pts = []\n","for p, r in zip(xyz[:,:2], radii):\n"," surfpts = circle(p, r)\n"," #remove points that fall within radius of another atom: \n"," surfpts = surfpts[(cdist(surfpts, xy)-radii).min(1)>-1e-3]\n"," canvas_coords = surfpts * np.array([gradx, grady]) + np.array([interceptx, intercepty])\n"," vdw_pts.extend([tuple(c) for c in canvas_coords])\n"," \n","\n","# requires installing concave_hull\n","# https://concave-hull.readthedocs.io/en/latest/\n","from concave_hull import concave_hull, concave_hull_indexes\n","# get concave hull points\n","conc_hull = concave_hull(\n"," vdw_pts)\n","conc_hull.append(conc_hull[0])\n","conc_hull = np.array(conc_hull)\n","ax.fill(conc_hull[:,0], conc_hull[:,1], c='whitesmoke', zorder=5, alpha=0.4)\n","ax.plot(conc_hull[:,0], conc_hull[:,1], c='k')\n","ax.set_xticks([])\n","ax.set_yticks([])"],"metadata":{}},{"id":"65331315-8a1d-43a0-88c1-cd63fcf11fc3","cell_type":"code","execution_count":null,"outputs":[],"source":[""],"metadata":{}}],"metadata":{},"nbformat":4,"nbformat_minor":5}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment