{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# freud.diffraction.StaticStructureFactorDirect and freud.diffraction.StaticStructureFactorDebye\n", "\n", "The `freud.diffraction` module provides two methods for calculating a one-dimensional [static structure factor](https://en.wikipedia.org/wiki/Structure_factor) $S(k)$ which can be used to characterize structure of crystals, liquids or amorphous phases.\n", "\n", "The `freud.diffraction.StaticStructureFactorDirect` class implements a \"direct\" $S(k)$ method. First, the following expression is computed over a $k$-space (reciprocal space) grid:\n", "\n", "$$S(\\vec{k}) = {\\frac{1}{N}}\\sum_{i=1}^{N}\\sum_{j=1}^{N}\\mathrm{e}^{-i\\vec{k}\\cdot(\\vec{r}_{i} - \\vec{r}_{j})}$$\n", "\n", "Then, the angular dependence is integrated out, resulting in $S(|\\vec{k}|)$, otherwise denoted $S(k)$. For an excellent introduction to the theory of scattering and $S(k)$, please refer to the documentation of the [dynasor package](https://dynasor.materialsmodeling.org/), which performs a number of calculations related to scattering. The **freud** library implements the core method of static structure factor calculation based on the dynasor package, with some additional performance optimizations in parallelized C++ code, as well as an interface to compute $S(k)$ that aligns with the APIs and conventions of the **freud** library.\n", "\n", "The `freud.diffraction.StaticStructureFactorDebye` class computes static structure factor based on the Debye scattering equation:\n", "\n", "$$ S(k) = {\\frac{1}{N}} \\sum_{i=1}^{N}\\sum_{j=1}^{N}{\\frac{\\sin(kr_{ij})}{kr_{ij}}} $$\n", "\n", "which is obtained by integrating out the angular dependence from the original formula. This implementation provides a much faster algorithm, but gives worse results than the \"direct\" method at low $k$ values.\n", "\n", "Note that freud employs the usual physics convention, as opposed to the crystallographic convention, with the following expression linking the two: $k = 2\\pi q$. The static structure factor is related to the radial distribution function, $g(r)$, by a Fourier Transform:\n", "\n", "$$S(k) = 1 + \\rho \\int_{V}\\mathrm{d}\\vec{r}e^{-i\\vec{k}\\cdot\\vec{r}}g(r). $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lennard-Jones Liquid Example\n", "\n", "One of the use cases for $S(k)$ is to characterize structure of liquids. The example shown here uses data generated by a HOOMD-blue simulation of a 1000 particle system subject to the Lennard-Jones potential. See the HOOMD-blue [documentation](https://hoomd-blue.readthedocs.io/en/latest/) and [examples](https://github.com/glotzerlab/hoomd-examples) for more information." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEYCAYAAACgDKohAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA80UlEQVR4nO3deXxU5dn/8c81k8m+k4QthF3ZQWRRsRYRUamKVVtx19YHcW2ttWr112rrVp9qW8S61X3FituDuOEOLiyKIKDskBBIQsieSTLL/fvjTELALBMya3K9X695ZZYz51xnBvLNfd/n3EeMMSillFL+sIW7AKWUUtFDQ0MppZTfNDSUUkr5TUNDKaWU3zQ0lFJK+U1DQymllN80NFSXIyLni8h74a5Dqa5IQ0OFnIgcKyKfi0iFiOwTkWUiMtH32iUisrQD6xogIkZEYhqfM8Y8b4yZcQh15YrIQhHZ66ttrYhc0tp2AklEpopIQTDW3cY2t4uIU0Sqm936HOK6Ql6/Co+g/AdQqjUikgosAq4AXgZigZ8A9eGsy+dZ4FugP1Y9o4Fe/r5ZRGKMMe4g1RasbZ9mjFkS8II6KJyfneogY4ze9BayGzABKG/lteFAHeABqhuXA34GfANUAvnAbc3esxMwvuWrgaOBS4ClzZYZCbwP7AOKgD+2sv1qYFwrr7W2nWXAP3zrvgO4DXiu2fsG+N4X43ucCTwJFAJlwOtAEuAEvM3W3wd4Crij2bqmAgXNHm8HbgTWYIVcDHAU8DlQjhWAU9v4LrYD0w96LgMr1Et89S0Ccpu93pH644B/+pYt9N2Pa74vvvr3AM+G+9+m3vy7afeUCrWNgEdEnhaRU0Qko/EFY8wGYC7whTEm2RiT7nupBrgISMcKkCtE5Azfa8f5fqb73vNF842JSAqwBHgH6xfZEOCDVmr7EnhQRGaLSN5Br7W2ncnAViAHuNOP/X8WSMQKshzgH8aYGuAUoNC37mRjTKEf6wI4F+szSQd6Am9hhVcm8HtgoYhk+7kusLqsn8RqbeVhhcH8Q6z/FqwQGweMBSYBtzZbVy9fnf2BOR2oUYWRhoYKKWNMJXAs1l/fjwElIvKmiPRs4z0fG2PWGmO8xpg1wIvAT/3c5KnAHmPMfcaYOmNMlTHmq1aW/QXwGfD/gG0isrpxrKUNhcaYB4wxbmOMs60FRaQ31i/XucaYMmOMyxjziZ/70Zp5xph837YvABYbYxb7Pqv3gZXAzDbe/7qIlPturxtjSo0xC40xtcaYKqwg/Okh1n8+8BdjTLExpgS4Hbiw2ete4M/GmPr2PjsVOTQ0VMgZYzYYYy4xxuQCo7BaAP9sbXkRmSwiH4lIiYhUYLVGsvzcXD9gi591lRljbjLGjMT6q3011i9VaeNt+X7W0VjLPmNMWQfe057m2+8P/KJZCJRjBXTvNt5/hjEm3Xc7Q0QSReQREdkhIpXAp0C6iNgPof4+wI5mj3f4nmtUYoyp83NdKkJoaKiwMsZ8j9V3P6rxqRYWewF4E+hnjEkDHgakjeWbywcGH0Jde4G/Y/2Sy2xjOwc/X4PVfdOo+UB6PpApIul+rKe9dbX0vnyssYH0ZrckY8w9rdTekuuBw4HJxphU9nfLySHUX4gVZI3yfM+19R4V4TQ0VEiJyDARuV5Ecn2P+2H1y3/pW6QIyBWR2GZvS8H6C7dORCYB5zV7rQSrm2NQK5tcBPQSkd+KSJyIpIjI5FZq+5uIjBKRGN9YyBXAZmNMqR/babQaOE5E8kQkDbi58QVjzG7gbeDfIpIhIg4RafylXAT08L2n+bpmikimiPQCftvOtp8DThORk0TELiLxvkNhc9t5X3MpWOMY5SKSCfy5E/W/CNwqItkikgX8yVejimIaGirUqrAGj78SkRqssPgO6y9cgA+BdcAeEdnre+5K4C8iUoX1i+flxpUZY2qx+t2X+bpkjmq+MV+//InAaVhH6WwCjm+ltkTgNawjj7Zi/ZV8uj/baba994EFWEc0rcIKreYuBFzA90AxviDwtbheBLb61t+H/YcAbwfe8623VcaYfGAW8EeskMsHbqBj/8//CSQAe7G+m3c6Uf8dWGMqa4C1wNe+51QUE2O0haiUUso/2tJQSinlNw0NpZRSftPQUEop5TcNDaWUUn7r8hMWZmVlmQEDBoS7DKWUihqrVq3aa4xpcfqZLh8aAwYMYOXKleEuQymlooaI7GjtNe2eUkop5TcNDaWUUn7T0FBKKeW3Lj+moZRSzblcLgoKCqir0wl24+Pjyc3NxeFw+P0eDQ2lVLdSUFBASkoKAwYMoO1Z77s2YwylpaUUFBQwcOBAv9+n3VNKqW6lrq6OHj16dOvAABARevTo0eEWl4aGUqrb6e6B0ehQPgcNjWi2+QMoWhfuKpRS3YiGRjR7/Ur476Xg9Ya7EqVUB9jtdsaNG8fIkSMZO3Ys999/P952/h9//PHHnHrqqSGqsHU6EB6t3A1Qvce6rX8dRp0Z7oqUUn5KSEhg9erVABQXF3PeeedRUVHB7bffHt7C/KAtjWhV1exSy5/+r7Y2lIpSOTk5PProo8yfPx9jDB6PhxtuuIGJEycyZswYHnnkkaZlKysr+fnPf86IESOYO3cuXq+Xxx9/nOuuu65pmccee4zf/e53ADz33HNMmjSJcePGcfnll+PxeDpdb0S0NHzXiX4G6IV1HeZHjTH/OmgZAf4FzARqgUuMMV+HutaIUbHL+jn2PPj2Bfh+EYw4Pbw1KRVlbv+/dawvrAzoOkf0SeXPp43s0HsGDRqE1+uluLiYN954g7S0NFasWEF9fT1TpkxhxowZACxfvpz169fTv39/Tj75ZF599VVmz57NmDFjuPfee3E4HDz55JM88sgjbNiwgQULFrBs2TIcDgdXXnklzz//PBdddFGn9i8iQgNwA9cbY74WkRRglYi8b4xZ32yZU4Chvttk4CHfz+6p0tfSOOYaKFgOn9wLw08DPSpEqajUeOnt9957jzVr1vDKK68AUFFRwaZNm4iNjWXSpEkMGjQIgHPPPZelS5dy9tlnM23aNBYtWsTw4cNxuVyMHj2a+fPns2rVKiZOnAiA0+kkJyen03VGRGgYY3YDu333q0RkA9AXaB4as4BnjPXJfiki6SLS2/fe7qeywPqZ3g9+8nt4fS788DYMmxneupSKIh1tEQTL1q1bsdvt5OTkYIzhgQce4KSTTjpgmY8//vhHh8g2Pr7sssu46667GDZsGJdeeilghdDFF1/M3XffHdBaI25MQ0QGAEcAXx30Ul8gv9njAt9z3VNlIcSlQVwKjP4FpPWDr58Od1VKqQ4qKSlh7ty5XH311YgIJ510Eg899BAulwuAjRs3UlNTA1jdU9u2bcPr9bJgwQKOPfZYACZPnkx+fj4vvPAC5557LgAnnHACr7zyCsXFxQDs27ePHTtanfHcbxHR0mgkIsnAQuC3xpiDOxpb6ncxraxnDjAHIC8vL6A1RoyKXZDmy0x7DPSfAls/DmtJSin/OJ1Oxo0bh8vlIiYmhgsvvLBp8Pqyyy5j+/btjB8/HmMM2dnZvP766wAcffTR3HTTTaxdu5bjjjuOn//8503r/OUvf8nq1avJyMgAYMSIEdxxxx3MmDEDr9eLw+HgwQcfpH///p2qXRr70cJNRBzAIuBdY8z9Lbz+CPCxMeZF3+MfgKntdU9NmDDBdMmLMD1yHCRlwwULrcdf/BvevRmu/wFSeoW3NqUi2IYNGxg+fHi4ywi4U089leuuu44TTjihQ+9r6fMQkVXGmAktLR8R3VO+I6MeBza0FBg+bwIXieUooKLbjmeA1T2V2qx3rvdY6+fuNeGpRykVFuXl5Rx22GEkJCR0ODAORaR0T00BLgTWishq33N/BPIAjDEPA4uxDrfdjHXI7aWhLzNCuOuhpgTScvc/12u09XP3t3DYjPDUpZQKufT0dDZu3Biy7UVEaBhjltLymEXzZQxwVWgqinCVvnM0Uvvg9Rr+smg9g3OSuTBzMOxeHdbSlFJdW0SEhuqgxnM0Uvvy4oqdPPX5dgCO6D2Akbu/bTt9lVKqEyJiTEN1kO9s8CLJ4u7F3zNlSA9+fexA/q84B6nIx1leEuYClVJdlYZGNPKd2Penj/bh8Rru/vkY/t+pIzhi8k8BeOn/FoWzOqVUF6ahEY0qC2lwpPLupmp+f9Lh5PVIBODk6dYZpN7C1WEsTinVnramRl+5ciXXXnttQLZz1113BWQ9zWloRCFveQHbXRkckZfOJccM2P9CYiblsb3pVfMDbo/OeqtUpGqcGn3dunW8//77LF68uGla9AkTJjBv3rwfvcftdnd4OxoaCoD60p0UeDL41ZSB2G0HDnvXZI5kONvYXlobpuqUUh1x8NTozS+2dNtttzFnzhxmzJjBRRddRElJCWeddRYTJ05k4sSJLFu2DIDq6mouvfRSRo8ezZgxY1i4cCE33XRT05nn559/fsDq1aOnolFVIbvNkUwfmPmjl2Jyj6DvniW8X7CbITlDw1CcUlHk7Ztgz9rArrPXaDjlng69pfnU6AdbtWoVS5cuJSEhgfPOO4/rrruOY489lp07d3LSSSexYcMG/vrXv5KWlsbatda+lJWVcdZZZzF//vymiz0FioZGtHE5SXCVU5vQi56p8T96OWPwBFgJ5dtWwXgNDaWiRWtTOp1++ukkJCQAsGTJEtav3z/5d2VlJVVVVSxZsoSXXnqp6fnG+aeCQUMjypiKXQiQnN3yRIyx/cZbd3Z/C8wOWV1KRaUOtgiCpfnU6Bs2bDjgtaSkpKb7Xq+XL774oilEGhljfjRterDomEaUKczfDECvfoNbXiA5hzJ7D9LL17f8ulIqohw8NXpbZsyYwfz585seN3Y9Hfx8WVkZAA6Ho2mK9UDR0IgyO7dZoTF4yOGtLlORMoSchnzqXJ2/HrBSKvAaB6hHjhzJ9OnTmTFjBn/+85/bfd+8efNYuXIlY8aMYcSIETz88MMA3HrrrZSVlTFq1CjGjh3LRx99BMCcOXMYM2ZMQAfCI2Zq9GDpalOjv/Xg9fys5D+YPxYisUktLrP96TmkbX2L3ZdvYESf1BBXqFRk66pTox+qqJwaXfmvYV8+1bbUVgMDIKnnEDKkmm0Fu0JYmVKqO9DQiCKF5U5SG4qpT+rd5nLpuVbXVXF+6KZLVkp1DxoaUWTF9n30ln04MnLbXM7RYyAAzj2bQ1GWUlGnq3fL++tQPgcNjSjy1TYrNJKz27nGb4bv9bLtQa9JqWgTHx9PaWlptw8OYwylpaXEx//4fK+26HkaUeTrrcVkSBWk9Gx7wfg0nDFppNXtoqbeTVKcfs1KNcrNzaWgoICSEr2EQHx8PLm5bfdcHEx/m0SJqjoXZSWFEA8k57S7fENqf/rVF7OpuJpx/dKDXp9S0cLhcDBw4MBwlxG1tHsqSuyuqCNLKqwHSe2Hhr3HQPKkmI17qoJcmVKqO9HQiBJFlXVkN4aGHy2NxJ5D6Ct72binLMiVKaW6Ew2NKLGnoo5sKbce+BEatswBOMTDvsKtwS1MKdWtaGhEiaLKOrKotB740T1Fhq/PVo+gUkoFkIZGlCiqrKevowpikyE2sf03ZAwAIKE6H4+3ex9aqJQKHA2NKLGnso5cRxUkZfv3htQ+eCSGXIooqqwLbnFKqW5DQyNKFFXWkWOrhOR2ztFoZLNTn5xLPykmf59e+lUpFRgRExoi8oSIFIvId628PlVEKkRkte/2p1DXGE5FlXVkUg7JfrY0ADKsw27zy5xBq0sp1b1ETGgATwEnt7PMZ8aYcb7bX0JQU0Rwe7yUVNWT5inzbxDcJy57EP2liJ3a0lBKBUjEhIYx5lNgX7jriER7qxuwGTcJ7gr/u6ewTvBLk1pKS4qCWJ1SqjuJmNDw09Ei8q2IvC0iI8NdTKgUVdbRo/Fw2w51Tw0AwF2q52oopQIjmuae+hrob4ypFpGZwOvA0JYWFJE5wByAvLy8kBUYLHsqOzaFSBPfuRr2ip1BqEop1R1FTUvDGFNpjKn23V8MOEQkq5VlHzXGTDDGTMjO7sBf5hHKmkKk3Hrgx9ngTXxTpKfVFVDv1uuFK6U6L2pCQ0R6iYj47k/Cqr00vFWFRtPhttCx0IhLoS42k34UsUuPoFJKBUDEdE+JyIvAVCBLRAqAPwMOAGPMw8DZwBUi4gacwGzTTa6isqeinsFxNeChY91TgDs1j/5O67DbQdnJwSlQKdVtRExoGGPObef1+cD8EJUTUYoq6/hpbDV4/JxCpBlb9lAGlyxhiR52q5QKgKjpnurOiirr6G2v6FjXlE987xH0kjKKS4qDUJlSqrvR0IgCeyrr6EFFh7umAGw5hwPgKf4h0GUppbohDY0IV9vgpqrOTbq3vGPnaDTKHgaAo2xTYAtTSnVLGhoRbk+FNUNtkqu0Q2eDN0nvj1scpFfrCX5Kqc7T0IhwRZX1xOAmznVo3VPYY6hIzCPXk091vTvwBSqluhUNjQh3yFOINFOXPpQhUqhTpCulOk1DI8Id8hQizdhzhpEnxewq0fkglVKdo6ER4Yoq6+jnqLIeHMqYBpDUdwQ2MVTt+j6AlSmluiMNjQhXVFnHoMQa68Ehdk8l544AwFOkoaGU6hwNjQi3p6KOfg5faBxi95T0GIoHG7HlmwNYmVKqO9LQiHBFlfX0iamE2JQOTyHSxBHPXkdvUqu2BLY4pVS3o6ERwbxeQ3FVHdm2ikPummrkTB1Mb9dO9tU0BKg6pVR3pKERwfbVNuDyGDK85YfcNdUotvcIBspuVm8vCUxxSqluSUMjghVX1gOQ7C7rdEsja+BoYsXDlo3rAlGaUqqb0tCIYCXVVmgkNOztfEuj53AAyneu7XRdSqnuS0MjghVX1pFCLTH15ZDeyWudZ1mXU7ft3YjL4+18cRGq3u1h6aa9FFfVhbsUpbqkiLkIk/qxkup6Bspu64Hvl/4hi0/FmdCLAdUFrC+sZGy/9E7XF2m+3lnGTf9dzdllj5FPLVtTJ5EwbBqzjh7FYL1qoVIBoaERwUqq6hnuKLIe9BjS6fXZeo5gSu03fLhxM2P7Tej0+iJFncvDve/8wJOfb+POxJc4L+YtXLYEHLUf4V11L7tWZbEtvgcZOX1JH38WHHF+uEtWKmpp91QEK6mqZ2RcEYgdMgZ2en1x028lXWqYsvwKqK8OQIXhV9vg5ldPreCJZduYP3g553nehEmX47ilAH79Ps4pN1CVM4HC+jj27VgHb1xJxfv3hrtspaKWhkYEK66qZ6h9D2T0h5jYzq8w90ie7nsbfes3Y16+ENzRfc5Gdb2bSx5fztatm1l41FZ+VvAvOPxncPLdYI+BfpNImnELI656iXE3f8g7U99gkXcKacvu5KunbqTB3XXHdpQKFg2NCLa3qp5+3kLo0cnxjGbiR/6Mm12XIVs+hLf/ELD1hlpRwWa23jeNp/fM4su4qzhy9a3Qdzyc9R+w2X+0fFJcDFdOG8YRv32ZL5JnMHn7w3xy33mU7dIrGirVERoaEWxvlZOeroKAjGc0OrJ/Bv/1TGVH/7Ph2xfBXR+wdQebs8HDyyvyuW3+49gfO56B9RvZM/RcmPl3OH8hXPJWu1Ot9M1M5ujfLWDLoAs5vvZdUh+bRNXTs2HrJ+D1hGhPlIpeOhAeoZwNHpLri3HE10NW4EJjWK8UEmPtfG6fSH/3K1CwEgZMCdj6g2VLcRW3P7uYIfs+4xbHi1Qn9KbyrOcYOHRcx1dmszH4ovms23AlK/57Lz/f+h5se9uaen7EGTD+Iug1KtC7oFSXoKERoUqq6hlk8x1uG8DuqRi7jXH90lm4N4/ZCLJ9aeSFhjGw4U0o/Aaq9rC3cBs9ir/jGakGB5hBx5PxiyfJSMjo1GZGDh9BxrUP8etnlpGz5xOusK1m1KqnkJVPwKwHYew5Adqh8Cst2Mie5QvZV11PccJgShIG0ddTyDTPUpK2vAX9JsPs5w98k7MMGmohrW94ilYRSUMjQpVU1wXuHI2DnDqmD398rZSKnsNI3/4ZcGNA199pWz6Aly8CWwzVsdnsrE1ib+IUJh1zAulDJiG9xoItMD2rfdITeOHKafxzSS6zPjmKEemX8ULaQ6S+NgfKtsFPbwSRgGwrpDwutn37GTtXvkWfPR8y1LuVHi0s5jSxbI3tz6DvF1G27n0yRp7Y9H6eOg1KNsCEX1ufQ1JLa1DdTcSEhog8AZwKFBtjftQ3ICIC/AuYCdQClxhjvg5tlaFTUlXPINmNx5GM/RCv2Neacyb248XlO3mrbAjnVb+HuOrAER/QbRwyY+DDOyCtHzvO/5QZ877i2CFZPHTBkcTGBGcILjbGxh9OHsbxw3K4bsFqJuy4gv/2yWTsx3dD8QY46U5Iyw3KtgPCWQ6rn4fSLTRU7WVvcSFpZd8xECf9jbAlbjhf9LuOtPFnMji3J7GlG5Di7ykxKbxUPoKF3+zheXMVJQv+wFkp9zO8bxpH73meCyrX8okZz7HLH6Nh5XN8O+RKhs26gfSkuHDvcWBUFcHal+G7VyF7mHXUXUL6/te9Huvfoz0GYwzV9W5sIiTFRcyvzbAQY0y4awBARI4DqoFnWgmNmcA1WKExGfiXMWZye+udMGGCWblyZaDLDbpnv9hO/8UXcHQfG44rPg34+lfnlzP/4Xn8x3GfNYA84NiAb+OQbFgEC87HnP4Al6wexsrt+/jg+qn0SgtNqFXXu7nn7Q089+UObkl9h1+7F2ATgclz4NjfQWJmSOrwS00pfPUw5quHkPoqauxpFLsT2WeSKU4cSuLh0xj7k9NIz+rV5mq8XkP+h4/Sf+kfeCjnT3xYmcuzddewKelI/jv0Xur3bGBW0b85xnzD+xxF6fR/8ItjhmO3RWALrHA1LL4BKgvB+H7pH3UF3mN+w9a91Xy3qxLjcjLpu9vpk/8WYjzUZ40ktvR7vCl92DP9ASoT84hZ9Th9Nr0AXhef2ibzesNE1jT0JdtexVG9YMhhIxkzdjyDs5Nx2K0/Zowx5O/chnvNQvrkLyLWU4Pt/Jchc1B4P5NDICKrjDEtngEcMaEBICIDgEWthMYjwMfGmBd9j38Aphpjdre1zmgNjfve+4Fzls2kz+jjsZ39n6Bs4/b/LuPW735G2YTfknXabUHZRod4vfDwFPA08NZPXuOql9byp1NH8KtjO39iY0ct27yXP7yyhrjqfJ4d/CF9d7wBqX3hys8hPi34Bez4AuJTIXOw1QosWg/rXqXuu/9Dqouwu6qJMS4A3vZM4gH3GWyNGcQZ4/py4dH9GdmngzV6PfCQ9dnTYzBsXwZXfQXp/azXjaHo3b+T/eVdbPb25r70Wzh9+jRmjOzZ9Esz6OqrYc1LMOxUSDkoCL1e+Hye1UpNysY98KeU1nqpL9lKXsUK5pnZ3F9/OonU8R/H3znKtoEnPCfzomcaW0xfxslm5jkeoI+U4sFOnLj4wDseE5fGFM9yErw1B2yuwdi5x30ezzGTQdnJDIyr5NSihzjZLMMuhrXeAfSVvXhtDt454hFOnjaVrOQ2Wmgrn4TP7ofhp8L4iyFnGOzbCuvfhD1rYdIcyGv3b+SA6SqhsQi4xxiz1Pf4A+BGY0ybiRCtoXHrf1fwl3UnYpt6M0wNzphDWU0Du/93MsaRxIg/foaEu+9+7Suw8NfUnv4oU9/uQU5qHK9fOYWYUP1SOkhpdT2XPbOS1fnlPDjFycyVl8FRV1jdGMH0/Vvw0nnWfbFhErOQmmI82PjSM5ytpje1kogkpLEr56f0GDiOUX1TOTIvk7REx6Fv94e34cXZ1v0Zd8IxV/9oEbP1YxpeuoS4hjIKTBbrbYeT3mcwuTHlJNcVEeOIpX7iFcQePoPEuJhD+zfl9ULRd1Z4xSZZz+3+Fl75NZRugsQsKmf+m1UxR7CusAIKVnB8wUOMbFjDB3IUt3r+h90NCQDEiJdHUx5nWsNHrB86l4GVK4kv/obCqfdT2P90quvcVNa5qKpz43BXM2HrQzjEQ90Rl5F3+DjiHXbrsPQtH0FVISRmQUIGtZ89QOLWd9mUPoVVjOT0imdx4GHzwPOpG3ku26QvNfnfcdq3V+DxuPiV5xYGjT6aw3qmkJMSR0aSg/JaF8VV9Qze/BQn5s+jLHEAqc4C7MZNfVIf4moKAfA4krG5atg15FxWDL4aitYzuPBNBpR9SVH8QDYmT2BzwhhS63aTW7uePOcGFibNZkvqZB6/ZOIh/VPoKqHxFnD3QaHxB2PMqhaWnQPMAcjLyztyx44dQa07GP706AL+UjgHzn4CRp0VtO1sePpaBm19gQ/OWMHMI0L/F30TjxsenAQx8dyc829eWrWL16+cEvaJFZ0NHn7z0je8t76IBX1fZtK+/0Mu/zR4h+S6G+DfR4EtBs9PbmD92pXkb93AF/UDWZd2PKdNGcvxh+fQNyMh8H/hGwPPnWn9RX/p29ZZ9S2p3I33u1cp2fAZjsKVpLpL2UMmhaYHfaSUXNnLN94hPOOZQUxiOrkZiWRm9yau/yR6pSfQOy2e9MRY4j2VxO/5mpjBP0V8Y2rG5cTz6lxiNryO1x5PUa/jKHbkMnLHs9TY03gx5RJmlC9ggLeApz0zOEwKmGJfR5Uk80bOXNblzCIpLobk+BhG901jwoBM0uJs8Npca/zC5oCzH4cRszr/WS1/FN671WqdDZkOp9xrBV1zpVtwPXkatuo9LGcEi11H8rl3JJUmkXocXGx/j+sdr/CWZzK/cV1FKrWcZf+UybYNfOEdwbveSewzKVwf818usb+DFxsO8VBt4vmC0QyWQgaxq2lzDcSwNWYwryWfR0H2cTx4/vhD2r2uEhrdqnvqnvvu4aaqu+HyT6H32KBtx/P9O9hfOodrYv/C/95wtfWXVTh8txBe+RWfHnE/F33RiyumDubGk4eFp5aDeLyGe97ewMufreWT+N9j73kYKXOXBOeoqi8fhnduZOWUR7hxTS+2lNQwPi+da08YynFDs7EFexzB47Z+thYYLSiprGNPZT3lzgbKq2vI2fIqIzY9Qkr9ngOW+8Izgrvc57HODOAc+0fcELOATKlmu+nJ3d6LWGsbzjz+xgT5gfnuWaRQyyn2FeRIOR94j+S+hGuJT8smL1W4vOZRhu9+DW9yT2zHXAtHXgJxbcxk7PVY3T+5E2Dw8YfwwbSieIM1fjJ4Wuv/HioL4auHrRZk6eYfvz72PMzp83B6hEqnm9KaevZWN7C3yjrxNjHWTnysnd4135O79WViBx5D7OhZ+1th5flQsBzSB1h/zMR0/kCFrhIaPwOuZv9A+DxjzKT21hmtofHwX69grucF+GPh/n8cwVBXgfnbAOa5zsAx/RaunBq4Ewn9Zgw8OpW6mnJG772TY4bk8MQlEyNuoPXzzXv5+KX7+aP7QZ7vfTPesecyqk8qw3undjpsG9xevt+2kyEv/YS13gGc47yRQdnJ/OGkYZw0smf4uw47yt0Axet9g9Hgzl+BfHovdmcptQl9SHQWsjt9POtyTueI/Kfp4dxGnT2JGG8DHw//C2WDTiMrOY6eyQ76yF7Seg9GDj7Mumjd/jGfaFGyEXatAlctuOsgIRPGnBOwQ8gDpa3QiJhjx0TkRWAqkCUiBcCfAQeAMeZhYDFWYGzGOuT20vBUGnxeryGnIZ/KhBxSgxkYAPFpSO+xnFK6iTM/2sIvjuxHdkqID6ncsQx2r+Z+2xz6ZiQxb/YRERcYAMcMyWLM729j57yPOXX3PK7YIfw/7yhiY2xMGpDJcYdlccLwnn5du6OyzsUnP5Tw2aYSvttVyabiKv4gzzLKXsX7edfw8OQJnDA8J3SDzIEWEwt9xu1/mHskHHEefD6PxG2fwuQ76D3yTHqLgOdqWP4Y8Wtfhhl3Mv1HJ5u2chJnz5FBKz9osg+zblEsoloawRCNLY19NQ3s+NvR9MzqQZ9r3wv+Bj+6C/PJvUxvuI9JEyZx95mjg7/NZmqePIuGnSs4wTOfBVcdz9CeKSHdfoft24p5YTbs3cjmEdewIOEcPt1cysYia7r5n43qxc2H7yJ39xLr6KPsYdBzFM7kPBav3c2r3xTw1dZSzpSPuM7xKok2Dza7gyT3PtyjZhN71r/DvIOqu4uKlobar6SyjkFSSFlaiC6UNOHXyNJ/cE+fzzhnRS+umTaEPukJIdn0F18s4+gdS3hWzuHhS4+N/MAAyByE/M+HsOg6hq6dx615K+GYMyjpcSTvrN3N4NW/JXfzGpySQIJxNr3tG8bwasPPqM4YwaKeTzCs7CNMv6OQniOsM7Bjk4k97oYw7phS7dPQiECVuzdyuNSyLydEze+UnjB2Nkd++zKZzODllfn8dnrwm9CPfLKFtPfvZXxMLLPm/InevSPoxLn2xCXDmY9C/2Pgs/vg7T+QDVwIeBMyWNLzOu7dO4Uap5Ps+h0ca1vH/8S9x/PcjamLQ5xeOPEvyNHXRFx/tlJt0dCIQJL/FQAxA48O3UaPvgbb189wS/Yy/ndFT66ZNjSo4wovfLWT5e8+z8Nxy5AjLqR37wiepqM1IjDhUutWvhN2fA41JdiOuIDpCRlM9y1mjMHjNdbJeGteRrZ/Zp3v0eeIsJav1KHQ0IhAiUUrqTSJpPcfE7qNZh8Gh8/kZ9sWcXPlCXy6sYTjh+UEZVOfrlpL2qLreTz2K0zWMOSnvw/KdkIqPc+6tUBEiLELEAfjL7RuSkUpbRdHoKyy1azmMJLjA3CJ14445lpiG8q5JGEpLy7fGZRNrPtuNWPePIkT7V/TcNwfkcs/06m3lYoiGhqd8M3OMlZu3xfYlTrLyXJuY2PsiNAfm593FOROZK7jLT7/fifFlXUBXX1ZTQObFv6VBHFRc+lHxE67MTDXPldKhYyGRif87Z3vuXPxhsCutGAFNgz5yaE97BWw+uin30ZaQxE32Z7jv6sKArr6+xd+zEzvR9SMmE1G/zDsn1Kq0zQ0OqGsxkVFrSuwK935JR5slKWHcDyjuQHHIsdcwwUxH7Dzy4V4vYE5j+fttbvJ3fgUMQKZJ14fkHUqpUJPQ6MTyp0NVNYFODTyv+J7BpCenh7Y9XbEtFupSBvG7+vm8/byNW0uesDJofnL4ZHjrPl4mimtrudvr33JhTEfwqgzIWNAEIpWSoWChkYnlNe6qHS6CdhZ9R4XpmAly91DyW5r7v1gi4kj5dwnSBMnQ969GM/Ll8KL58Kb11pXifN5Y/UuJtyxhHkfbLJaJB/daU1hveACqKsArFC55bXvmOV6i0Sc2H5yXZh2SikVCB0ODRFJEpEwTYUaOepcHurdXho8Xurd3sCsdM8axO1kpffw0M//dBBbr5HsPOZOEj3VVG5bZc2kufoFeHImnvJd3L14A795aTUOu43739/I7f9ZAFs/tqac3rcNXrsC4/Vw51sbWLVuA3Pj34ehJ0XnfEFKqSbtnqchIjZgNnA+MBGoB+JEpARrEsFHjTGbglplBCpvNpZR4XQFZkrxndZJfau8QzkrNfzXYR4y43Iu230kX24t5ZOrptKj6HO8Cy6g/IGf8mHtDVxw1DH86dSRLFixk+S3r6bWFs/nQ27h+NyjsL93M98+fjWH7yjgi4QvsHsM6BQZSkU9f1oaHwGDgZuBXsaYfsaYHOAnwJfAPSJyQRBrjEjlzoam+5XOAI1r5H9JdXxv9tCD4b1TA7POTrrplGE4XR7+umg9N3+bxayaW/C4XbyV9FfuGF9LbIyNC0c4OMP+BYtjpnPZy5s57pPD+TptOuN2vcAsx1fYj7wYuWo59Du0q4gppSKHP2eETzfG/Oi3ojFmH7AQWCginbjGZHRq3tIIyGC4MbDzKzbGjiQ7JY5eqZFxjYAhOcmcNymPZ7/cQazdxuxJx8ERJxD7+rnwzCw45znY/imClzPm3kHyngSe+nwb5269gOv7TeTSiy5DUrLCvRtKqQBpNzQaA0NE7jDG3Nr8NRGxG2M8LYVKV3dAaDjdnV9h2Xao3sNncbMYm5sWURfdueHkwxmQlcTM0b3oneab/fZX71iXBn1xNthjYfjpxGQN5OQsOHlULwrKaumVOits1/dWSgVHR/5H9xWRcxsfiEgOsCTwJUWHiubdU4FoaWxfCsBbVYMZk5ve+fUFUGq8g18fO3B/YAAk58Alb0HuRHDVwDHXHPCe3IxEDQyluqCOTFh4OfCuiGwBDPAkcGNQqooCB7Y0AhMarrhMNtb1ZXRuWufXFwrxaXDR61C6BXqOCHc1SqkQ8OfoqWeAr4FvgKuAFwA3cIYxpoWrpHcPFU4XNgGvgcq6TnZPGQM7lrEzdTxUCGMjrKXRppg4DQyluhF/+g+e9i33K6zAGACUAReIyNnBKy2ylTtdZCbFERdj63xLo3wHVOSzkhHkZiSQmaST+CmlIpM/A+EfAB80PhaRGGAEMBY4CnglaNVFsIpaF+mJDkSsVken+MYzFlUNZmy/9M4Xp5RSQdJuS0MOOozHGOM2xqwxxjxrjPl9S8t0B+XOBtITHKTGx3R+IHz7UrwJPfisPCt6xjOUUt2SXyf3icg1InLAZclEJFZEponI08DFwSkvcpXXukhLcJCW4Oj8Ibfbl1GaNQEQxmhoKKUimD+hcTLgAV4UkUIRWS8i24BNwLnAP4wxTwWxxohUXusiLdFBaoLD75bGtr013PfeD9S7PfufLNsBFTtZHzsGERjdV0NDKRW5/BnTqAP+Dfzbd+Z3FuA0xpQHubaIVuF0kZ4Qi9tj2L63xq/3vPp1AQ98uJmCMif3/3KsdQKfbzzjg7rDGZSVREp8tzu5XikVRfwZ0/hAREZC09nhE4GrRWRSsIuLVC6Pl+p6N+mJDlITYvw+5HZzcTV2m/DaN7v4xxLfHI87lkFiD94tSou4k/qUUupg/nRP5Rpj1gGIyDHAs0Ae8JSI/DxQhYjIySLyg4hsFpGbWnh9qohUiMhq3+1Pgdp2RzUeLZWe6CA13kGl0+XXNTU2F1czbVgOvzgyl3kfbOL5r3bg2vIpe3tMoKjapeMZSqmI588Z4ZXN7l8EPGyMudE3jcibwGudLcJ3fY4HgROBAmCFiLxpjFl/0KKfGWNO7ez2OqvxbPC0BAe1DR7cXoPT5SExtvWP0+Xxsm1vDSeO6Ml1Jx5GYYWTt994kfNj8/nnvukATOifGZL6lVLqUPkTGpt9J/F9CpwBnAlgjCkWkUBd9GESsNkYsxVARF4CZgEHh0ZE2N/SiKWm3tP0XFuhsaO0FrfXMCQnGYfdxqO/PBweuoQq+0Cmn/I7ZqenMUoHwZVSEc6f7qnrsOad2gV8bYz5HMA3KJ4coDr6AvnNHhf4njvY0SLyrYi83TjOEg6NkxU2HnIL7c90u7m4CoChOSkAJH36V5Kce0j55SNMHZmngaGUigr+HD21BzhRRGzGmObXNT0e6wJNgdDSyYEHDxJ8DfQ3xlSLyEzgdWBoiysTmQPMAcjLy2tpkU5p7J5KT3BQ5Tvctr3DbjcXVwMwOCcJtn0KKx+Ho66CvMkBr08ppYLF77mrDwoMjDHvGWPmBKiOAqBfs8e5QOFB26s0xlT77i8GHCLS4tV9jDGPGmMmGGMmZGdnB6jE/ZpCwzcQDu3PdLupuJq+6QkkihveuBoyB8G0W9t8j1JKRZpIueDBCmCoiAwUkVisa5K/2XwBEenVOF2J73BfG1Aa8kqxJisUgZR46+Q+8K+lMSQnGbZ9Yk1QeNJdEJsYinKVUipgOnI9jaAxxrhF5GrgXcAOPGGMWScic32vPwycDVwhIm7ACcw2/hznGgQVtQ2kxjuw24TUeOsjbGtMw+s1bCmp5qhBPWDji+BIgsHTQlWuUkoFTESEBjR1OS0+6LmHm92fD8wPdV0tKXdaM9wC+1sabXRP7Sp3UufyMjQ7CZa9C4OPt65DoZRSUSZSuqeiSnmti3RfWDjsNhJj7W1Oj944CD7akQ+Vu+Cwk0NSp1JKBZqGxiEod7pIS9x/oaTU+LYnLdzkO9x2UJk1zxRDZwS1PqWUChYNjUNQUdvQdH4GYM0/1caYxubiarKS40jYtgT6jIeUnqEoUymlAk5D4xCUO/d3T0H7LY3NxdWM7+GCgpXaNaWUimoaGh3k9Roqmw2EA21eU8MYw6biak6KWwsYOOykEFWqlFKBp6HRQVX1bryGA7un4lvvniqpqqeqzs2R9cshpTf0HhuqUpVSKuA0NDqoonb/ZIWN2mppbCquxoGb3H1fWAPg3e9y6kqpLkRDo4PKfZMVNh/TsK4T3vI1NTYXV3O87RtiXNVw+MyQ1amUUsGgodFBzeedapQa78BroLr+x11UBWW1XOD4EJPSG4ZMD1mdSikVDBoaHVTu3H8BpkapCb6pRFq47Kt33w6OlTXIEReCPWJOwFdKqUOiodFBFbW+a2kc1NKAlqcSGVvyBgaB8ReFpkCllAoiDY0Oan6p10atzj/lcXFs1TusTZgE6f1QSqlop6HRQeVOF4mxduJi7E3PNbU0Du6e+mExmaaMb7LPCGGFSikVPNrJ3kHNJytslNZKS8OsfJLdpgd7ex8XsvqUUiqYtKXRQRUHTVYIzQfCm4XG+jeQrR/xkvt4MlP0YktKqa5BQ6ODKpwNP2ppJMc1uxBTQw28eQ28fBF1OWN5xnMiWcmxLa1KKaWijoZGB1U63aTEH9irF2O30T+umqE7XoRHjoOvn4Vjr2PNiQsoJ4UeSXrBJaVU16BjGh1U0+BuallgDGxeAsv+xUeyFFu+gezhcPGbMPA49q7dDUAPbWkopboIDY0Oqm3wkBBrh/zlsOQ22LEM0vN4IW42m7Knc/tlZzctW1pdD0BWsrY0lFJdg4ZGR7icnOj6kLlblsK3ayApB2b+HcZfzJv/WYXNe+Die6sbEIGMREfL61NKqSijodGWje/B6ufAWQ515Zh92/ibrZIydx7MuAOOvBTikgHrBL9d5c4D3l5aU09GYiwxdh06Ukp1DRoabaktheINEJ8OyT1x9xzLBV/154Tjz2DOMUMOWDQ7JY6vd5Yd8NzeqgZ6JOl4hlKq69DQaMu4c62bT1llHV99+QGnxf24u6l/j0T21TRQVecixXeGeGlNvQ6CK6W6FO036YCaBg8AibH2H72Wl2mdwLdzX23Tc6XVDfTQQXClVBeiodEBtQ3W3FKJsT9uoDWFRun+0NhbXU+2hoZSqguJmNAQkZNF5AcR2SwiN7XwuojIPN/ra0RkfKhrrPW1NJLiftzS6N/DCo0dvpZGg9tLZZ1bxzSUUl1KRISGiNiBB4FTgBHAuSIy4qDFTgGG+m5zgIdCWiRQU9/Y0vhxaKTEO8hMimWHr6VRWmOdo6HdU0qpriQiQgOYBGw2xmw1xjQALwGzDlpmFvCMsXwJpItI71AW6Wwa02j5+IG8zER27qsBrPEM0LPBlVJdS6SERl8gv9njAt9zHV0mqBoHwpPaDA2rpbG36WxwDQ2lVNcRKaEhLTxnDmEZa0GROSKyUkRWlpSUdLq4Rk7fQHhCC91TYI1rFJbX4fJ497c0dLJCpVQXEimhUQA0vx5qLlB4CMsAYIx51BgzwRgzITs7O2BF1rQxEA5WS8PjNewqczaNaWSlaGgopbqOSAmNFcBQERkoIrHAbODNg5Z5E7jIdxTVUUCFMWZ3sAraWFTF7ooDpwWp9Q2Ex8e01tJIAqwjqPZWNxAXYyOplVaJUkpFo4gIDWOMG7gaeBfYALxsjFknInNFZK5vscXAVmAz8BhwZTBrOu2BpTy1bPsBz9U2eEiMtWOztdRT1vxcjRr2VteTlRyHSMvLKqVUNIqYaUSMMYuxgqH5cw83u2+Aq0JVT0q8g8o69wHP1TR4Wj1yCiAnJY64GBs799X6zgbXQXClVNcSES2NSJSaEHPgNb+xBsJbOkejkc0m5GUmsqO01pp3Sk/sU0p1MRoarUiJd1DVYkuj7TGK/j0Sm7U0dBBcKdW1aGi0IjU+hqqDWhq1DW6S4tru0cvLTLJaGtUNesU+pVSXo6HRipT4GCqdB4dG+y2NvMwEnC4PDR6vntinlOpyNDRakdpC91RtvT/dU0lN93UgXCnV1WhotCIlPqaFMQ13m0dPAeT5ZrsFPRtcKdX1aGi0IiXegdPlweXxNj3n9KN7KjcjgcZTM7SloZTqajQ0WpEab7Uomrc2avwYCI+LsdM7NR5AB8KVUl2OhkYrGq/z3XgElcdrqHN5SXC0Py1IYxdVpp6noZTqYjQ0WpHia2lUOq2WhtPV9mSFzQ3NSSEnJQ6HXT9epVTXEjHTiESa1IQDWxq19a1fH/xgvzvxMC46un/wilNKqTDR0GhFU0vDN6ZR23TVvvZbGhlJsWRo15RSqgvS/pNWpPrGNBrnn6pp8L+loZRSXZWGRitSmwbCfWMaHWhpKKVUV6Wh0YrkpkNuG1sa/g+EK6VUV6Wh0Qq7TUiKtTcdPdWRgXCllOqqNDTakJrg2H/0lHZPKaWUhkZbms8/VasD4UoppaHRFuuSrzqmoZRSjTQ02pB6QEvDCo34GA0NpVT3paHRBuuSr/vPCE+MtWOzSZirUkqp8NHQaENKfMz+M8Jd7U+LrpRSXZ2GRhsaj54yxvhaGjoIrpTq3jQ02pASH4PLY6h3e6nx4wJMSinV1WlotKHxmhqVTpdfV+1TSqmuLuz9LSKSCSwABgDbgV8aY8paWG47UAV4ALcxZkKwa0ttNtNtTYOb5Hau2qeUUl1dJLQ0bgI+MMYMBT7wPW7N8caYcaEIDGg+aaG2NJRSCiIjNGYBT/vuPw2cEb5SDpRyUEtDB8KVUt1dJIRGT2PMbgDfz5xWljPAeyKySkTmhKKw5tcJr63XloZSSoXkT2cRWQL0auGlWzqwminGmEIRyQHeF5HvjTGftrK9OcAcgLy8vA7X2yg1Yf91wmsbPCTpmIZSqpsLyW9BY8z01l4TkSIR6W2M2S0ivYHiVtZR6PtZLCKvAZOAFkPDGPMo8CjAhAkTzKHW3djSqHC6cLo8JDi0paGU6t4ioXvqTeBi3/2LgTcOXkBEkkQkpfE+MAP4LtiFJcXasQkUVdZZj3WyQqVUNxcJoXEPcKKIbAJO9D1GRPqIyGLfMj2BpSLyLbAceMsY806wCxMRUuIdTaGRoAPhSqluLuy/BY0xpcAJLTxfCMz03d8KjA1xaYB1BFVTS0MHwpVS3VwktDQimtXSqAf0AkxKKaWh0Y7UZi0NPeRWKdXdaWi0IyXegdtrHYClA+FKqe5OQ6MdjfNPgXZPKaWUhkY7UhMcTfe1e0op1d1paLQjRVsaSinVREOjHQeGhrY0lFLdm4ZGOxqnRwd0GhGlVLenodGOxvmnEmPt2GwS5mqUUiq8NDTa0dg9pV1TSimlodGuxqOndBBcKaU0NNqlLQ2llNpPQ6MdGhpKKbWfhkY7Go+e0qv2KaWUhka74mJsxNpteritUkqhodEu60JMMdrSUEopIuAiTNHghpMOZ3BOcrjLUEqpsNPQ8MPsSXnhLkEppSKCdk8ppZTym4aGUkopv2loKKWU8puGhlJKKb9paCillPKbhoZSSim/aWgopZTym4aGUkopv4kxJtw1BJWIlAA7wl3HQbKAveEuIsi6+j7q/kW/rr6Pndm//saY7JZe6PKhEYlEZKUxZkK46wimrr6Pun/Rr6vvY7D2T7unlFJK+U1DQymllN80NMLj0XAXEAJdfR91/6JfV9/HoOyfjmkopZTym7Y0lFJK+U1DQymllN80NEJMRLaLyFoRWS0iK8NdT2eJyBMiUiwi3zV7LlNE3heRTb6fGeGssbNa2cfbRGSX73tcLSIzw1ljZ4hIPxH5SEQ2iMg6EfmN7/ku8T22sX9d6TuMF5HlIvKtbx9v9z0f8O9QxzRCTES2AxOMMV3ipCIROQ6oBp4xxozyPXcvsM8Yc4+I3ARkGGNuDGedndHKPt4GVBtj/h7O2gJBRHoDvY0xX4tICrAKOAO4hC7wPbaxf7+k63yHAiQZY6pFxAEsBX4DnEmAv0NtaahOMcZ8Cuw76OlZwNO++09j/QeNWq3sY5dhjNltjPnad78K2AD0pYt8j23sX5dhLNW+hw7fzRCE71BDI/QM8J6IrBKROeEuJkh6GmN2g/UfFsgJcz3BcrWIrPF1X0Vl183BRGQAcATwFV3wezxo/6ALfYciYheR1UAx8L4xJijfoYZG6E0xxowHTgGu8nV9qOjzEDAYGAfsBu4LazUBICLJwELgt8aYynDXE2gt7F+X+g6NMR5jzDggF5gkIqOCsR0NjRAzxhT6fhYDrwGTwltRUBT5+pEb+5OLw1xPwBljinz/Sb3AY0T59+jrB18IPG+MedX3dJf5Hlvav672HTYyxpQDHwMnE4TvUEMjhEQkyTcQh4gkATOA79p+V1R6E7jYd/9i4I0w1hIUjf8RfX5OFH+PvkHUx4ENxpj7m73UJb7H1vavi32H2SKS7rufAEwHvicI36EePRVCIjIIq3UBEAO8YIy5M4wldZqIvAhMxZqGuQj4M/A68DKQB+wEfmGMidqB5Fb2cSpWt4YBtgOXN/YdRxsRORb4DFgLeH1P/xGr3z/qv8c29u9cus53OAZroNuO1Rh42RjzFxHpQYC/Qw0NpZRSftPuKaWUUn7T0FBKKeU3DQ2llFJ+09BQSinlNw0NpZRSftPQUEop5TcNDaWUUn7T0FAqxERkuog8G+46lDoUGhpKhd5Y4JtwF6HUodDQUCr0xgLfiEiciDwlInf55kdSKuLFhLsApbqhsVizjb4L/McY81yY61HKbzr3lFIh5Juiey+wA2uCvC/CXJJSHaLdU0qF1ghgBeAGPGGuRakO09BQKrTGAp8Ds4EnRaRnmOtRqkM0NJQKrbHAd8aYjcCNwMu+LiulooKOaSillPKbtjSUUkr5TUNDKaWU3zQ0lFJK+U1DQymllN80NJRSSvlNQ0MppZTfNDSUUkr57f8DVM6cLhhT3OoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import freud\n", "import gsd.hoomd\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "bins = 100\n", "k_max = 30\n", "k_min = 3\n", "sfDirect = freud.diffraction.StaticStructureFactorDirect(\n", " bins=bins, k_max=k_max, k_min=k_min\n", ")\n", "sfDebye = freud.diffraction.StaticStructureFactorDebye(\n", " num_k_values=bins, k_max=k_max, k_min=k_min\n", ")\n", "\n", "with gsd.hoomd.open(\"data/LJsampletraj.gsd\", \"rb\") as traj:\n", " for frame in traj:\n", " sfDebye.compute(frame, reset=False)\n", " sfDirect.compute(frame, reset=False)\n", "\n", "plt.plot(sfDebye.k_values, sfDebye.S_k, label=\"Debye\")\n", "plt.plot(sfDirect.bin_centers, sfDirect.S_k, label=\"Direct\")\n", "plt.title(\"Static Structure Factor\")\n", "plt.xlabel(\"$k$\")\n", "plt.ylabel(\"$S(k)$\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Crystal Comparison Example\n", "\n", "The static structure factor $S(k)$ can also be used to characterize and compare crystal structures. In the below example we compare the computed static structure factors $S(k)$ of a face-centered cubic (fcc) crystal and simple cubic (sc) crystal." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEYCAYAAAC6MEqvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsJUlEQVR4nO3de5wddX3/8dfnXPaeze4mu7mwhCRAACEGMAFBpSgCIlWo1VarFvuzpVrro+2vtkprf7VqW36/1mqtt6K2Tb0VHqUqRaBCREEFTLglIEhCyGWzIbvZ7Gbve26f3x8zJzkJ2WSze267834+HudxzplzZuY7e5J5z/cyM+buiIhINMUqXQAREakchYCISIQpBEREIkwhICISYQoBEZEIUwiIiESYQkCqnpm908y+X+lyiMxFCgGZMTN7tZn91MwOmtkBM/uJma0LP3uPmf34JJa13MzczBL5ae7+DXe/ahrl6jSz281sf1i2LWb2nsnWU0xmdrmZdZVi2cdZ5w4zGzOz4YLH0mkuq+zll8ooyX8AiQ4zawbuBN4P3AbUAK8BJipZrtDXgCeB0wjKsxpYPNWZzSzh7pkSla1U636Tu99X9AKdpEr+7eQkubseekz7AawFBib57BxgHMgCw/nvAdcCjwODwG7gYwXz7AI8/P4wcAnwHuDHBd85F7gXOADsA/5skvUPA+dP8tlk6/kJ8Olw2Z8EPgZ8vWC+5eF8ifB9G/CvQDfQD3wHaATGgFzB8pcC/wZ8smBZlwNdBe93AB8GNhOEVgJ4JfBTYIAg0C4/zm+xA3j9UdNaCUK6NyzfnUBnwecnU/5a4DPhd7vD17WF2xKW/0Xga5X+t6nH1B5qDpKZeg7Imtl6M7vGzFrzH7j7M8D7gIfcvcndW8KPRoDfBFoIAuH9ZnZ9+Nll4XNLOM9DhSszs3nAfcA9BDumM4ANk5TtYeDzZvZ2M1t21GeTrediYDvQAfz1FLb/a0ADQTB1AJ929xHgGqA7XHaTu3dPYVkA7yD4m7QAi4DvEYRRG/Ah4HYza5/isiBo8v1XgtrQMoKd++emWf4/Jwil84E1wEXARwuWtTgs52nAjSdRRqkghYDMiLsPAq8mODr+MtBrZneY2aLjzPNDd9/i7jl33wx8C/ilKa7yl4EX3f1T7j7u7kPu/sgk330b8CDwF8ALZvZEvq/iOLrd/Z/cPePuY8f7opktIdhZvs/d+9097e4/muJ2TOaz7r47XPe7gLvc/a7wb3UvsAl443Hm/46ZDYSP77h7n7vf7u6j7j5EEGy/NM3yvxP4uLv3uHsv8FfAuws+zwF/6e4TJ/rbSfVQCMiMufsz7v4ed+8EziM4Qv/MZN83s4vN7H4z6zWzgwS1hYVTXN2pwPNTLFe/u3/E3c8lOKp+gmAnaceZbfcUy5EvywF37z+JeU6kcP2nAW8r2KkPEATukuPMf727t4SP682swcz+2cx2mtkg8ADQYmbxaZR/KbCz4P3OcFper7uPT3FZUiUUAlJU7v4sQdv3eflJx/jaN4E7gFPdfT7wJcCO8/1Cu4HTp1Gu/cDfE+y02o6znqOnjxA0l+QVdizvBtrMrGUKyznRso41326CtvWWgkeju988SdmP5Y+Bs4CL3b2Zw81gNo3ydxMEU96ycNrx5pEqpxCQGTGzs83sj82sM3x/KkG79sPhV/YBnWZWUzDbPIIj0HEzuwj4jYLPegmaFVZOsso7gcVm9odmVmtm88zs4knK9n/N7DwzS4R9Ce8Htrl73xTWk/cEcJmZLTOz+cBN+Q/cfS9wN/AFM2s1s6SZ5Xey+4AF4TyFy3qjmbWZ2WLgD0+w7q8DbzKzq80sbmZ14dDNzhPMV2geQT/AgJm1AX85g/J/C/iombWb2ULg/4RllFlMISAzNUTQmfqImY0Q7PyfIjgCBfgB8DTwopntD6f9HvBxMxsi2JHcll+Yu48StFv/JGwCeWXhysJ27SuBNxGMQtkKvHaSsjUA3yYYWbOd4Cj2zVNZT8H67gVuJRix8yhBCBV6N5AGngV6CHfsYY3oW8D2cPlLOTxkdQfw/XC5k3L33cB1wJ8RhNZu4E84uf+3nwHqgf0Ev809Myj/Jwn6JDYDW4DHwmkyi5m7anAiIlGlmoCISIQpBEREIkwhICISYQoBEZEIm3UXkFu4cKEvX7680sUQEZlVHn300f3u/pJLjsy6EFi+fDmbNm2qdDFERGYVM9t5rOlqDhIRiTCFgIhIhCkEREQibNb1CYiIlEI6naarq4vx8dl9IdS6ujo6OztJJpNT+r5CQEQE6OrqYt68eSxfvpzjX228erk7fX19dHV1sWLFiinNo+YgERFgfHycBQsWzNoAADAzFixYcFK1GYWAiEhoNgdA3slug0IgirJpePzrkMtVuiQiUmEKgSh64QH47gdg7xOVLomIFPjsZz/LOeecwzvf+c6yrVMdw1GUTQfPuUxlyyEiR/jCF77A3XffPeVO3WJQTSCKPHfks4hU3Pve9z62b9/Om9/8Zj7xiU/wW7/1W6xevZqXv/zl3H777QDcc889XHjhhaxZs4YrrriiKOtVTSCKFAIix/VX//00P+8eLOoyX7a0mb9807mTfv6lL32Je+65h/vvv5+/+7u/Y/78+WzZsgWA/v5+ent7+Z3f+R0eeOABVqxYwYEDB4pSLtUEokghIFLV7rvvPj7wgQ8cet/a2srDDz/MZZdddqipqK2trSjrUk0gihQCIsd1vCP2cnD3lwz1PNa0YlBNIIoUAiJV7aqrruJzn/vcoff9/f1ccskl/OhHP+KFF14AUHOQzIBCQKSqffSjH6W/v5/zzjuPNWvWcP/999Pe3s4tt9zCW97yFtasWcOv//qvF2VdZWsOMrMdwBCQBTLuvtbM2oBbgeXADuDX3L2/XGWKLPfwWSEgUk127Nhx6PX69etf8vk111zDNddcU9R1lrsm8Fp3P9/d14bvPwJscPczgQ3heym1QzUBr2w5RKTiKt0cdB2Qj7v1wPWVK0qEqDlIRELlDAEHvm9mj5rZjeG0Re6+FyB87ihjeaJLISAioXIOEX2Vu3ebWQdwr5k9O9UZw9C4EWDZsmWlKl90KAREJFS2moC7d4fPPcC3gYuAfWa2BCB87plk3lvcfa27r21vby9XkecuhYCIhMoSAmbWaGbz8q+Bq4CngDuAG8Kv3QB8txzliTyFgIiEytUctAj4dni2WwL4prvfY2YbgdvM7L3ALuBtZSpPtCkERCRUlhBw9+3AmmNM7wOKcyk8mTqFgIiEKj1EVCrh0MliOk9ApJqMjIxw7bXXsmbNGs477zxuvfVWNm7cyKWXXsqaNWu46KKLGBoaKuo6dQG5KFJNQOT47v4IvLiluMtcvBquufm4X7nnnntYunQp3/ve9wA4ePAgF1xwAbfeeivr1q1jcHCQ+vr6ohZLNYEoUgiIVKXVq1dz33338eEPf5gHH3yQXbt2sWTJEtatWwdAc3MziURxj91VE4gihYDI8Z3giL1UVq1axaOPPspdd93FTTfdxFVXXVWSy0cXUk0gihQCIlWpu7ubhoYG3vWud/GhD32Ihx9+mO7ubjZu3AjA0NAQmUxx7w2umkAUKQREqtKWLVv4kz/5E2KxGMlkki9+8Yu4Ox/84AcZGxujvr6e++67j6ampqKtUyEQRQoBkap09dVXc/XVV79k+sMPP1yydao5KIoUAiISUghEkW4qIyIhhUAU6aYyIsfkc+D/xMlug0IgitQcJPISdXV19PX1zeogcHf6+vqoq6ub8jzqGI4ihYDIS3R2dtLV1UVvb2+lizIjdXV1dHZ2Tvn7CoEoUgiIvEQymWTFihWVLkbZqTkoihQCIhJSCESRQkBEQgqBKFIIiEhIIRBFCgERCSkEokgni4lISCEQRTpZTERCCoEoUnOQiIQUAlGkEBCRkEIgihQCIhJSCESRQkBEQgqBKFIIiEhIIRBFCgERCSkEokjnCYhISCEQRTpPQERCCoEoUnOQiIQUAlGkEBCRUFlDwMziZva4md0Zvm8zs3vNbGv43FrO8kSWQkBEQuWuCfwB8EzB+48AG9z9TGBD+F5KTSEgIqGyhYCZdQLXAl8pmHwdsD58vR64vlzliTSFgIiEylkT+Azwp0DhnmeRu+8FCJ87jjWjmd1oZpvMbNNsvwl0VVAIiEioLCFgZr8M9Lj7o9OZ391vcfe17r62vb29yKWLIA0RFZFQokzreRXwZjN7I1AHNJvZ14F9ZrbE3fea2RKgp0zliTadLCYiobLUBNz9JnfvdPflwNuBH7j7u4A7gBvCr90AfLcc5Yk81QREJFTp8wRuBq40s63AleF7KTX1CYhIqFzNQYe4+w+BH4av+4Aryl2GyFMIiEio0jUBqQSFgIiEFAJRpBAQkZBCIIoUAiISUghEkUJAREIKgSjSeQIiElIIRJFqAiISUghEkU4WE5GQQiCS1BwkIgGFQBSpOUhEQgqBKFIIiEhIIRBFCgERCSkEokghICIhhUAUKQREJKQQiCKdLCYiIYVAFOk8AREJKQSiSM1BIhJSCESRQkBEQgqBKFIIiEhIIRBFCgERCSkEokghICIhhUAUKQREJKQQiCKdJyAiIYVAFOk8AREJKQSiSM1BIhJSCESRQkBEQgqBKFIIiEhIIRBFCgERCSkEokghICKhsoSAmdWZ2c/M7Ekze9rM/iqc3mZm95rZ1vC5tRzliTyFgIiEylUTmABe5+5rgPOBN5jZK4GPABvc/UxgQ/heSk3nCYhIqCwh4IHh8G0yfDhwHbA+nL4euL4c5Yk8nScgIqGy9QmYWdzMngB6gHvd/RFgkbvvBQifOyaZ90Yz22Rmm3p7e8tV5LlLzUEiEipbCLh71t3PBzqBi8zsvJOY9xZ3X+vua9vb20tWxshQCIhIqOyjg9x9APgh8AZgn5ktAQife8pdnkhSCIhIqFyjg9rNrCV8XQ+8HngWuAO4IfzaDcB3y1GeyFMIiEgocbIzmFkjMO7u2ZOYbQmw3sziBMFzm7vfaWYPAbeZ2XuBXcDbTrY8Mg0KAREJnTAEzCwGvB14J7COYLhnrZn1AncBt7j71uMtw903AxccY3ofcMU0yi0zoRAQkdBUmoPuB04HbgIWu/up7t4BvAZ4GLjZzN5VwjJKMRUOC1UIiETeVJqDXu/u6aMnuvsB4HbgdjNLFr1kUhqFO36dJyASeSesCeQDwMw+efRnYRs/xwoJqVJHhIBqAiJRdzKjg04xs3fk35hZB3Bf8YskJaUQEJECJzM66HeB/zGz5wku+fCvwIdLUiopHYWAiBSYyuigfwceAx4HPgB8E8gA17v7ttIWT4ouv+O3mEJARKbUHLQ+/N7/IgiA5UA/8C4ze2vpiiYlkd/xxxIKARE5cU3A3TcQXOYZADNLAC8D1gCvBP6zZKWT4isMgcxEZcsiIhU3leYgcz88ltDdM8Dm8PG1Y31HqtgRNYGxypZFRCpuSieLmdkHzWxZ4UQzqzGz15nZeg5f/0eqXT6rY3HAda6ASMRNZXTQGwj6A75lZiuAAaCeIEC+D3za3Z8oVQGlyAprAhCEgFnlyiMiFTWVPoFx4AvAF8IzgxcCY+EloWW2eUkI5KjAFcVFpEqc8H+/mW0ws3Ph0JnB64DfN7OLSl04KYFjhoCIRNVUDgE73f1pADO7lKAzeBnwb2b2K6UsnJTAoRCIH/leRCJpKiEwWPD6N4EvufuNwOXojOHZRzUBESkwlRDYZmZvDa8VdD3h3b/cvQeoLWHZpBQUAiJSYCoh8EcE1w3aAzzm7j8FCDuJm0pYNikFhYCIFJjK6KAXgSvNLOZ+xB7jtQQ3nJHZRH0CIlJgylcRPSoAcPfvE5wnILPJoZPFVBMQEQ0Qj55jnSwmIpGlEIga9QmISAGFQNSoT0BECigEoqZaagK7fwa5bGXWLSKHKASiphpC4MB2+OqVsE23qBapNIVA1FRDCEwMBc/jg8f/noiUnEIgaqqhTyCbCZ5z6fKvW0SOoBCImmqoCeR3/lmFgEilKQSiphpOFsvv/FUTEKm4soSAmZ1qZveb2TNm9rSZ/UE4vc3M7jWzreFzaznKE2nVcLLYoZpApvzrFpEjlKsmkAH+2N3PAV4JfMDMXgZ8BNjg7mcCG8L3UkrqE5ixj//3z/nknT+vdDFEimLK1w6aCXffC+wNXw+Z2TPAKcB1BPclAFgP/BDdo6C0qqJPIAyBWdon8MTufuIx3ZdZ5oay9wmY2XLgAuARYFEYEPmg6JhknhvNbJOZbert7S1bWeekqgiBfJ/A7DxZLJNz0lldc0nmhrKGgJk1AbcDf+juUx4k7u63uPtad1/b3t5eugJGQTWEwCzvGE5nnXRWl9uQuaFsIRDehOZ24Bvu/l/h5H1mtiT8fAnQU67yRFY1hMAsbw5KZ3NkVBOQOaJco4MM+CrwjLv/Q8FHdwA3hK9vILx1pZRQVXQMz+6aQCabI51TTUDmhrJ0DAOvAt4NbDGzJ8JpfwbcDNxmZu8FdgFvK1N5oqsqagKze4hoOuvEdB8GmSPKNTrox8BkwymuKEcZJPSSk8UqsDOb5UNE09kccdfoIJkbylUTkGpRFTWB2d0nkMk5OVUEZI5QCERNNfQJzPIhoulsjpipJiBzg0IgaqqhJjDLO4bT2RxxhYDMEQqBqKmGEJjtzUFZJ6dLL8ocoX/KUVMNIVDGmsCPt+4v6old7h6eMawhojI3KASippr6BEo8RHTH/hHe9dVHuP/Z4p2DmL9chDtk1Tssc4BCIGqqoiZQniGiwxOZI56LIVNwkphqAzIXKASiphpuKlOmPoGJTLBtqUzxtrHwwnEKAZkLFAJRU003lSnxENH8TrqYO+vCZen6QTIXKASipiqag8rTMZzfYaeKuLPOqCYgc4xCIGqqomO4PM1BqZI0BxX0CahjWOYAhUDURLAmUMwj9kzBjj+jmoDMAQqBqKmGECjTENF8x3Cp+gTUHCRzgUIgaqohBMo0RDQ/kqdkzUHqGJY5QCEQNVHsEyhmc1DWj/laZLZSCERNNdQEDg0RLW1z0KHRQSWqCRQzXEQqRSEQNUefLEYlbipT3hAobp+AOoZlblEIRE1VnCxW3jOGi9l2f+RlI9QcJLOfQiBqqqFPoMw1gdKdJ6CagMx+CoGoqaY+gVnYMZxWx7DMMQqBqKmKEAivGVSuy0YUsSagy0bIXKMQiJpqCIHC5qAS9kmkdLKYyAkpBKKmGkKgsAZQwn6B/IXjdBVRkckpBKKmKjqGC3b8JQyBkjQH5dQcJHOLQiBqqq0mUMLO4cMdw8U7YtdVRGWuUQhETTXcWSxbnuYgnSwmcmIKgaipipPF0hBLBq/LURMo6uggdQzL3KIQiJp8CFjsyPfllMtCsiF8XcIQKPH9BHTGsMwFCoGo8VwQAJUMgWwaknWHX5dIKTqGC5el0UEyF5QlBMzsX8ysx8yeKpjWZmb3mtnW8Lm1HGWJvGoIgVwaEmEIlHKIaCkuJZ3LEY8ZMVNzkMwN5aoJ/BvwhqOmfQTY4O5nAhvC91JqlQ6BXC5Y56HmoFJ2DBf/PIFM1knEjGQ8pmsHyZxQlhBw9weAA0dNvg5YH75eD1xfjrJEXsVDIGz+SdYHz7OsYziVzVETjwUhkFFzkMx+lewTWOTuewHC547JvmhmN5rZJjPb1NvbW7YCzkmVDoH8Tr8MHcP5GkDOIVukMf2ZrJOIG4m4HXFZaZHZalZ0DLv7Le6+1t3Xtre3V7o4s5t7ldUEStccNFFQAyhWbSCTy5HI1wTUMSxzQCVDYJ+ZLQEIn3sqWJboqHRNIH8F0XwIlKEmAMXrHE5lPGgOipk6hmVOqGQI3AHcEL6+AfhuBcsSHZ4Ds4IQKPPRbLZ8fQLpbA6zw6+LIagJGIl4TGcMy5xQriGi3wIeAs4ysy4zey9wM3ClmW0FrgzfS6kdqgnY4ffllD/yL9MQ0aaaxKHXxXB4dJDp2kEyJyRO/JWZc/d3TPLRFeVYvxQ4IgSsCjqGSztEtLUxwdBEpmg1gVQ2RzIeHDulizjqSKRSZkXHsBRRPgQgeC57TSDc6Ze4OcjdSWVzNNYGxzlFaw4KQyAYHaSagMx+ZakJSBWpdAiUaYhofuROUxgCE0UbHRQMEQV1DMvcoBCImkqHQJmGiOZ30I01+ZpAcY7aU5kcyVjsiHWIzGZqDoqaiodAeYaI5juC881BxTtPwEkmwpPFdJ6AzAEKgajJnywGFW4OKm2fwKGaQG38iPczlcnmSMTyJ4upJiCzn0IgavLnCUAYAmU+mi3TENGJo2sCRdphp7NOMh4OEZ1pTSCXg2GdIymVpRCImko3Bx3dMVzimkBTkZuD0uHooGQ8NvNrB/38O/CZ1TB69LUVRcpHIRA1R4RABc4TOHqIaIlqAvkj/4aaIjcH5ZyVqedYkX5u5jWB/hcgMw7D+4pSNpHp0OigqKm2mkCphohmjhwiWsyawK/2fh6AO/j4zBY21n/ks0gFKASiptIhcKgmkL+9ZGlrAsU+WSydzdFMPzliZJhhTWBUISCVpxCImmoJgXhtsP5yDREt0nDOTNZpZIAccdLM8G+X3/mrT0AqSH0CUVPpEMg3B8UTEEuWtGO4lUFef+8bONd2FK05KJdNU58dpj47RHamtRg1B0kVUAhETaXPE8gf+ccSwaNUHcOZHGfaHhqGd3J+bFvRmoMaskMAxMhRnxue2cLGwhqAQkAqSCEQNe5HnSdQoZpALBnUBkpYE1hoBwFot4GiXfGzOXfw0Oum7ODMFqaagFQBhUDUvKQ5qNwni+X7BJJBEJRwiGj7oRA4WJSTxdydebnDO/55uUF8un8/d4WAVAWFQNRU+jyBbEFzUDxZ0o7hfE1gUZFCIJNzWm3o0PtWG5r+DexTI5BNBa/H1DEslaMQiJpKdwwfXRMo2VVEnXYGgLAmUITmoEzWaTsiBIanf8JY4dG/agJSQQqBqKl4CBzVJ1CymkD2yD6BItQE0rkcrRSEAEOkp3vpiPyOv2YejA3MuGwi06UQiJpKh0D+yD9W6iGifqhPYCH9pNNFCIFMjjYbIh2vJ2uJoCYw3RpGvglowUrVBKSiFAJRU+kQyKWD9cZipR0iWjA6qIYM8fQMR/JwuE9goqaVVE0LrQxN/xaT+R1/20pIDUMmNb3lDOyGFx6Y3rwztLlrgINjpQlxKR+FQNRUOgSy6aAGACUdIppKZ2lnAJ/fCUDdxP4ZLzOdzdHGEKl8CNjw9JuZDoXA6cHz+MD0lvODT8A3fq1kf8fJjExkeOsXH+Lz928r63ql+BQCUVPxk8UyQacwlHSIaDw1QI1lsUWrAWhMzTwEMtmgJpCqbSVd00KrDc28Y7htZfA83UtHdG2EzBj0/Hx680/T5q5gxNWmHRrZNNspBKKm0jeVyaaDZiAo6RDR5Fhf8GJxEAIN6b4ZLzNfE8jUtpKubQ2ag6ZbExg9AIl6aF4SvJ9Ov8DoATiwPXi957HplWOaHt8dlPepPYNMZLJlXbcUl0Igaip9nsARNYFEyYaI1k30Bi/CEGhKz/yINZ11Wm2YTG0bmdrWGQ4RHYCGNqhvDd9PIwT2PHr4dXeZQ2DXABD0vfy8+yT6W9Lj8Mg/B+dJnKxyH7BEhEIgairdJ5Ar7BMoXU2gLhUe+S9cRYok8zIzrwlk02PMszEydW1k61ppYZj0dI+Cx/qDAJhxCBh0XgR7Hp9eOabB3Xl8Zz+far+LC2zroUCYko1fgbv/FH52y8mt9PGvwz+ugaEXT24+OSGFQNRUOgSymaBDGMKaQGlCoH4i3Ok3dXAw3kZzZubDMD1st8/VByGQtCy5iYMnmGsSYwdmHgJdm6D9bFhxWdAnkBqdXllO0p6BMS4c+wm/OvR1vlT7WZ7ZsWdqM6ZG4SefCV7/7MtT/+3TY7Dh4zCwE+79P1Oa5c7N3Vx58/e488kplu1kjU+x9uM+vVrPSxZTulqQ7icQNZUOgVxBn0AJO4Yb0gdIkyBZ38pgopX52eM0BxUG07E+zjlD42lyw0Gw5OrbsEwNAPt79rK3YxGLm+t46Pk+tvcOc3HLAGeceQ6WqIVcFmJxRiYy/MfG3dQlY7z8lBbOHj7ASPMKntw5wWUWx8YOQP8O0j3beOSZ53lyVz+vv+AMVp23jmfHW3lwyzZy2Qyvu/AcFtQ61reVlq5NcPa17KhZxUrPctM//QuxumZ+p/lnLDjtZfS0X8L2n91NsrGFVWsupT7dR66uhbF5KzmwZxvj42NMxJsYHBygNf0iyxMHsKYO9gzm6Nn9C84efYz5fpDR5pX01J5GX/PLWLL6Cn78fB9/lPhP0g0dtI/28sbtn2D0nkup8QmGs0kOxBewj4UctCbOb51gcVsrzFtM7rF/JzbSy4YFv8EVfd/kkX/+PTp9Hwdp4MXEKQw3LWdi/krmxydY0fcg46kUPbFFzGeItcP72D7/ElZuvpVtvaPU50ZIxetJN3Qw3nQaz04spLa+gbNbYcu2HbRu/2/ujT/OjtsXseVH55OsrWdvspN0TQtnNo7yYrqJ3elmkjU11NbUMHKwj7EdGzklu4d5iQzjp11OvKmd3MFu5o/tonVsFw3pAbqbzqV1fBedw1sYaFhOf+tqxpPzOdiwgp7BEWr2bWaodjFjzSupSziX9XyDxePbGantYG/LK9hZu4qa0RdpS0xQn4yxayTJ4HiGhKdY0ZIg1tjG9sxCRvv30ZJI0dpUz70vNpJIDXBN7Wbs2k9xzuq1Rf2/ohCImgqFgLuzs2+U07Jp7KghosMTGVKZHI21cWoTccjlgnvv5h8TQzC0l1xDOyPNK6mLQ3p8hKGxFAdpIrvzYVpe+B4dp6wgPq+D0YkUpw5vZsBaaDdjOLGA01I/J/3sPezr6aV/73b6R1PU9m+lc+xZlqR389T8X+Jzjb/PZTxOR80ENXUNpEaHeHKgjk3747wp9lPWxreCgde3QaYWgMw9f8HWe8bpT4xxMNPK+dbLmbEd9Fob+xvP5Izhx+hPLGRj7izGU83EGWev9bEqtp27e5Zy07MbeayukfkPfpr4g58iCbw6fLAheMRynfyW7SVGjs0Pnc5S66LJxgH40MYmfpRNs7EO/nboz2EIUj1xap7PMg84Pf8DPHXk79F5nN8q7KqmyxeyzReysu8+XhGec9H7s2Ze7W2cE9tN5uovs/3JB3jt9q+ReehnQbMbKVrMWXmM5caAu7Pr+Gj/9XzX7ufintvo8oXMB15mdx/x3QlPkCHOapsA4OHcOby393f5bmInC7t/wF5vo5Fxllo/tZbh5QXzngWM1jQzccH7GX56I/P3P0a9pTjbDtfaVgCXHFW+NEn66ztJZbJ0bv/0oek93sIOX8wQbVw49iMGaOKLues4d2g7K0ceooMhGsNyDsWaaRweIjYcHLnvYjFf8LeydLSb14z/mDPsbsap4aAHt1e92EYxM1LUMNYdp5UhzrKgiTFFgtj+HBda8H/0hdRp5MaLPxrLSlnNmFIBzN4A/CMQB77i7jcf7/tr1671TZs2laVs1c7d6R2aYEFTLfGYMZHO8ND2A+wZGOPa1Utoaah5yTypz78Kb+6k9t23wlevJjU+zPbT3sbizpU0Ln8Fu9PzsYlB6saCm5/v3rOb1MEeOhqM5LwOepjPxp1DdO/rIZnq56KOHBMjBxkYHGR5c4xceoyBkQmSnedzZsMI8/ueYE/NSra8OErLyPO8Jv5zsnWt3Lb2Vl7/zJ+zqu8HbMmdxrbcKcRixtWJx464UueJjHgtjTbBhCeotSNrFZtrX8HLb/oBX/vKP/DLuz9Fqx15/f8DzOe5+Jlsn5jH2xM/xDHik9wtLB2r45n4WTw93salv/dlFrOf2n9+Jel4AwcaVtCTruXUeD/1DY08u/AqbPsPaZ3oZkvdK1iQ28+ZmW20+gCebGSkbhGDNR10nft+DnasI/s/H6VlvItfNK7FF57NeatO52VLmvn6jzbTsf8RLrWnaDr9YjyWYOLZ++ibdxb7FlzMUMMynmUFDXVJbmh7mvqhneRq5vHj2tcw0vUU7UPPcMYr38TQYD9dzz3BSE07dak+WkZ3Em8/ndr6JpLpQRrntTBc0862iTbiY/tprzdWrTqLsbpF9AynSGVyLIiNENv5IONb7qAxtZ+GBadQ+9ZbwJ3nnt3ChhdryVqSxc01nJYcZEmsj5qJAR7ZX0N37wESQ91MLDyH5WdfyJXnLiHZ8xTZvhfoXvw62uc3UEcK73ueVM9zjKWyHFjyGtpa25g/upPUU3dgZ19LzeKzGU9lGE5lGUtlScSNsfEUmYN7WG49DAyP8osBOG/lqbQuPQOS9cFvl80xls7SlB1iePAAm/trWF4/xinJ4MZA4xMp4jUN1HeuhkQQ7j27foFn07QvXkasvrnwPx2Y4e4cHEszmsoSN0gM7aapLkntguVB88/BLkiPwKLVeDxJ30iKpDnNDGMNCzgwmqZnaJwzO+YRjwWj9XbsHyGTTrG8dojEvHYysVq6+oY4zfZhNQ0w/3jRfWJm9qi7v6QaUdEQMLM48BxwJdAFbATe4e6TDnouaQikxwE/9I8HoHdogrGJDE21MZqSMJaL0TOUom6km/p4lpqOM2hIxkmkh/DhffSOZNg9bIxk48wf76Yx1ctoNsZINkGKBLV1DWQsyfDIMNmux2lMH6AhliGemyDtMcasgUGvJ5VoIlHbSCw1SH3mIM0+SHaoj3R6gjGrZ3n2BZoyA2zPtFOTiHNqcoCOid087ct5ILeatfFtnBofwHC6E0up8xSLM10szPVxn6/la8v+mv+954+40J8+4k+w19tYRD8xm96/iwmSxIAkQXtvt7exOFzeUO1ihrJJ7sys42/G38prEk/z23X3s7whxaKxbZBL85PYOp4eX8hoLkFNXSONjY3UNTYzVruQpfF+TsnuZYIknqijLhGjLdVNrnUFW095C7c+/Dzx1EEuPWMRq09t4dxVZ9BQ30BX/yg/2LKLeS8+REt7J8tWvZxTW+qoaWgGM0ZTGeqe/z6xp2+HC99NduE5jIwMMW9eM3Zwd/AfeuXleG0zB8fSh8N1pA/qWyAWn9ofp/BeDiJlVq0hcAnwMXe/Onx/E4C7/+1k8xQjBB76979gxfZvYubECB6Nf/o0T911C2u3fIwJkmSJEyNH3LMk7fAIkJTHGaSRhRYcrfb5POpIHaoOTkfa40yQJEGWOjt2Z9mEJzkYa8ZiCep9lO7EqQzXtHMq+xiccF7MNuGtK7kwtYmGoRd4se509sSXgjuL0l1MWC37apYxv62Dhxteyzf3tPPLzc9zybweas56Pdt37KTlwJOsym5ltHkl/Q3L8ZyzcNESmhcspXsoR3a4h+bsAVa01VFT3wyNCxlJtFDTOJ9kbSOZWA2xWJyYZxnc+SS9mVpy809jUe0EzTWxYEhkKJPNEY8Zlt8pugePWIxczknnckHTkIgUxWQhUOk+gVOA3QXvu4CLj/6Smd0I3AiwbNmyGa+0rn0luw+sI+tGDiPnxsXxOPXL13Hv/vcRmxgAzxKLJ2htaqCmpobxrDGRhfrcKG0M8uT8sxnPxmnp38yINdKfWMBAbAEdTUmWNOSotzTDtYsYrltMY8JpiGdJepqJ8VFiuRSNdbW0rlxLrHUZRozG/A4xk4KJweCRGoW6ZmhYQG2ygY6Co8hVBdvTTkHbby4HqWEW1zWz+Kjtzn/nXOC9AFx+6LMLL1gHvHXSv9nRy8prLHidKHjVvPIVHK5Ez3vJfIn4UQPTzA4dJcdiRu1Uj65FZEYqXRN4G3C1u/92+P7dwEXu/sHJ5lGfgIjIyZusJlDp8wS6gFML3ncC3RUqi4hI5FQ6BDYCZ5rZCjOrAd4O3FHhMomIREZF+wTcPWNmvw/8D8EQ0X9xP2q4ioiIlEylO4Zx97uAuypdDhGRKKp0c5CIiFSQQkBEJMIUAiIiEaYQEBGJsIpfQO5kmVkvsBNYCMz8xrGzU5S3HbT92v7obv9Mtv00d28/euKsC4E8M9t0rLPfoiDK2w7afm1/dLe/FNuu5iARkQhTCIiIRNhsDoGTvFP1nBLlbQdtv7Y/uoq+7bO2T0BERGZuNtcERERkhhQCIiIRNutCwMzeYGa/MLNtZvaRSpen3Mxsh5ltMbMnzGzO313HzP7FzHrM7KmCaW1mdq+ZbQ2fWytZxlKaZPs/ZmZ7wn8DT5jZGytZxlIxs1PN7H4ze8bMnjazPwinR+L3P872F/X3n1V9AtO5Mf1cY2Y7gLXuHomTZczsMmAY+Hd3Py+c9v+AA+5+c3gg0OruH65kOUtlku3/GDDs7n9fybKVmpktAZa4+2NmNg94FLgeeA8R+P2Ps/2/RhF//9lWE7gI2Obu2909BfwHcF2FyyQl5O4PAAeOmnwdsD58vZ7gP8acNMn2R4K773X3x8LXQ8AzBPclj8Tvf5ztL6rZFgLHujF90f8oVc6B75vZo2Z2Y6ULUyGL3H0vBP9RgI4Kl6cSft/MNofNRXOyOaSQmS0HLgAeIYK//1HbD0X8/WdbCNgxps2e9qzieJW7XwhcA3wgbC6QaPkicDpwPrAX+FRFS1NiZtYE3A78obsPVro85XaM7S/q7z/bQiDyN6Z39+7wuQf4NkETWdTsC9tL8+2mPRUuT1m5+z53z7p7Dvgyc/jfgJklCXaA33D3/wonR+b3P9b2F/v3n20hEOkb05tZY9hBhJk1AlcBTx1/rjnpDuCG8PUNwHcrWJayy+8AQ7/CHP03YGYGfBV4xt3/oeCjSPz+k21/sX//WTU6CCAcDvUZDt+Y/q8rW6LyMbOVBEf/ENwf+ptzffvN7FvA5QSX0N0H/CXwHeA2YBmwC3ibu8/JztNJtv9ygqYAB3YAv5tvI59LzOzVwIPAFiAXTv4zgnbxOf/7H2f730ERf/9ZFwIiIlI8s605SEREikghICISYQoBEZEIUwiIiESYQkBEJMIUAiIiEaYQEBGJMIWAyAyZ2evN7GuVLofIdCgERGZuDfB4pQshMh0KAZGZWwM8bma1ZvZvZvY34XVfRKpeotIFEJkD1hBcyfJ/gK+4+9crXB6RKdO1g0RmILzU735gJ8GFvB6qcJFEToqag0Rm5mUElzjPANkKl0XkpCkERGZmDfBTgntb/KuZLapweUROikJAZGbWAE+5+3PAh4HbwiYikVlBfQIiIhGmmoCISIQpBEREIkwhICISYQoBEZEIUwiIiESYQkBEJMIUAiIiEfb/AZ02mm+YyQakAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sf = freud.diffraction.StaticStructureFactorDirect(bins=200, k_max=25, k_min=1)\n", "fcc_system = freud.data.UnitCell.fcc().generate_system(10, sigma_noise=0.10)\n", "sf.compute(fcc_system)\n", "plt.plot(sf.bin_centers, sf.S_k, label=\"fcc\")\n", "\n", "sc_system = freud.data.UnitCell.sc().generate_system(10, sigma_noise=0.10)\n", "sf.compute(sc_system)\n", "plt.plot(sf.bin_centers, sf.S_k, label=\"sc\")\n", "\n", "plt.title(\"Static Structure Factor\")\n", "plt.xlabel(\"$k$\")\n", "plt.ylabel(\"$S(k)$\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Calculation of Partial Structure Factors\n", "\n", "Both methods support calculation of partial structure factors according to [Faber-Ziman decomposition](https://freud.readthedocs.io/en/latest/modules/diffraction.html#freud.diffraction.StaticStructureFactorDirect). In the conventions adopted in **freud**, the summation of partials reproduces the total scattering. In this example we load a simulation trajectory of $\\text{GeS}_2$ and calculate the Ge-Ge partial, the S-S partial and the mixed Ge-S partial (which is the same as S-Ge partial). The calculation of the partials requires the usage of `query_points` and `N_total` parameters for the compute method." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEYCAYAAACgDKohAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABiPklEQVR4nO3dd3iUVdrA4d+Zkpn0npBGQu89NAWliSKK6Oq66qqAirjqt+ja1rLr6tpd115RUbGxNhBRUESKgtTQe00hlfRkkinn++OdhAApk2QmM4FzX9dc094582SSzPOeLqSUKIqiKIordN4OQFEURWk/VNJQFEVRXKaShqIoiuIylTQURVEUl6mkoSiKorhMJQ1FURTFZSppKGccIcR1Qoil3o5DUc5EKmkobU4IMUoI8ZsQolgIcVwI8asQYqjzuWlCiNXNKCtFCCGFEIaax6SUH0spJ7YgrkQhxJdCiHxnbNuEENMaeh93EkKMEUJkeKLsRt7zsBCiUghRVucS38Ky2jx+xTs88g+gKA0RQoQAi4DbgPmAHzAaqPJmXE4fAVuAZLR4+gEdXH2xEMIgpbR5KDZPvfelUsqf3B5QM3nzs1OaSUqpLurSZhcgFShq4LlegAWwA2U1xwGTgc1ACZAOPFrnNUcB6Ty+DBgJTANW1zmmD/AjcBzIAR5s4P3LgIENPNfQ+/wK/NdZ9r+BR4F5dV6X4nydwXk/AngfyAIKgW+AQKAScNQpPx6YC/y7TlljgIw69w8D9wNb0ZKcARgB/AYUoSXAMY38Lg4DE055LBwtqec541sEJNZ5vjnxm4AXncdmOW+b6v4szvizgY+8/bepLq5dVPOU0tb2AnYhxAdCiElCiPCaJ6SUu4BZwBopZZCUMsz5VDlwAxCGlkBuE0JMdT53nvM6zPmaNXXfTAgRDPwE/ID2RdYVWNZAbGuB14QQfxJCdDzluYbeZzhwEIgBnnDh5/8ICEBLZDHAf6WU5cAkIMtZdpCUMsuFsgCuQftMwoBY4Du05BUB3AN8KYSIdrEs0Jqs30erbXVESwavtjD+h9CS2EBgADAMeLhOWR2ccSYDM5sRo+JFKmkobUpKWQKMQjv7fgfIE0IsFELENvKaX6SU26SUDinlVuBT4HwX3/ISIFtK+R8ppUVKWSql/L2BY68CVgGPAIeEEGk1fS2NyJJSviKltEkpKxs7UAgRh/blOktKWSiltEopV7j4czTkZSlluvO9/wwsllIudn5WPwIbgIsbef03Qogi5+UbKWWBlPJLKWWFlLIULRGe38L4rwMek1LmSinzgH8B19d53gH8U0pZ1dRnp/gOlTSUNiel3CWlnCalTAT6otUAXmzoeCHEcCHEciFEnhCiGK02EuXi2yUBB1yMq1BK+YCUsg/aWXsa2peqaORl6S7GURPLcSllYTNe05S6758MXFUnCRShJei4Rl4/VUoZ5rxMFUIECCHeEkIcEUKUACuBMCGEvgXxxwNH6tw/4nysRp6U0uJiWYqPUElD8Sop5W60tvu+NQ/Vc9gnwEIgSUoZCrwJiEaOrysd6NKCuPKB59G+5CIaeZ9THy9Ha76pUbcjPR2IEEKEuVBOU2XV97p0tL6BsDqXQCnl0w3EXp+/AT2A4VLKEE40y4kWxJ+FlshqdHQ+1thrFB+nkobSpoQQPYUQfxNCJDrvJ6G1y691HpIDJAoh/Oq8LBjtDNcihBgGXFvnuTy0Zo7ODbzlIqCDEGK2EMIkhAgWQgxvILZnhBB9hRAGZ1/IbcB+KWWBC+9TIw04TwjRUQgRCvy95gkp5THge+B1IUS4EMIohKj5Us4BIp2vqVvWxUKICCFEB2B2E+89D7hUCHGhEEIvhDA7h8ImNvG6uoLR+jGKhBARwD9bEf+nwMNCiGghRBTwD2eMSjumkobS1krROo9/F0KUoyWL7WhnuAA/AzuAbCFEvvOxvwCPCSFK0b545tcUJqWsQGt3/9XZJDOi7ps52+UvAC5FG6WzDxjbQGwBwNdoI48Oop0lT3Hlfeq834/A52gjmjaiJa26rgeswG4gF2cicNa4PgUOOsuP58QQ4MPAUme5DZJSpgOXAQ+iJbl04F6a93/+IuAP5KP9bn5oRfz/RutT2QpsAzY5H1PaMSGlqiEqiqIorlE1DUVRFMVlKmkoiqIoLlNJQ1EURXGZShqKoiiKy874BQujoqJkSkqKt8NQFEVpNzZu3Jgvpax3+ZkzPmmkpKSwYcMGb4ehKIrSbgghjjT0nGqeUhRFUVymkoaiKIriMpU0FEVRFJed8X0aiqJ4n9VqJSMjA4tFLWrrS8xmM4mJiRiNRpdfo5KGoigel5GRQXBwMCkpKTS+0rzSVqSUFBQUkJGRQadOnVx+nWqeUhTF4ywWC5GRkSph+BAhBJGRkc2u/amkoShKm1AJw/e05HeimqfOUkVff4PObMJ/yBCMMTHeDkdRlHZCJY2zUNX+/Rz7e+3eQBiTOxKQmkrAkFQChqZiTExUZ4XKGeeJJ57gk08+Qa/Xo9PpeOuttxg+/OT9uBwOB7Nnz+bnn39GCIHZbGb+/PnNavN3RVpaGllZWVx8sbZ9+8KFC9m5cycPPPBAg6+ZO3cuGzZs4NVXX3VrLM2lksZZqGzFCgAS33yD6oOHqNiwgbKfllH85VcARMyYQex993ozREVxqzVr1rBo0SI2bdqEyWQiPz+f6urq0477/PPPycrKYuvWreh0OjIyMggMDHRrLDabjbS0NDZs2FCbNKZMmcKUKVPc+j6eopLGWahs+S+YevQgeMwYGDOGyBnTkQ4HVfv3k/3PRyn7+WeVNJQzyrFjx4iKisJkMgEQFRXV4HFxcXHodFp3b2Ji/TvlpqSkcPXVV7N8+XIAPvnkE7p27cq3337Lv//9b6qrq4mMjOTjjz8mNjaWRx99lKysLA4fPkxUVBSrV6+msrKS1atX8/e//53KysraWkRDZfgKlTTOMvbiYio2bybypptOelzodJi7dydozBjy/vtfbIWFGMLDvRSlcib717c72JlV4tYye8eH8M9L+zT4/MSJE3nsscfo3r07EyZM4Oqrr+b8888/7bg//vGPjBo1ilWrVjF+/Hj+/Oc/M2jQoHrLDAkJYd26dXz44YfMnj2bRYsWMWrUKNauXYsQgjlz5vDss8/yn//8B4CNGzeyevVq/P39T2tqmjt3bm25jZXhC9ToqbNM+a+/gt1O0Jgx9T7vP3AgAJatW9suKEXxsKCgIDZu3Mjbb79NdHQ0V1999Ulf1DUSExPZs2cPTz31FDqdjvHjx7Ns2bJ6y7zmmmtqr9esWQNo81EuvPBC+vXrx3PPPceOHTtqj58yZQr+/v5NxtpYGb7AZ2oaQoj3gEuAXCll33qevw6433m3DLhNSrmlDUM8I5StWIE+LAz/Af3rfd6/X1/Q66lISyOonjMxRWmtxmoEnqTX6xkzZgxjxoyhX79+fPDBB/Tq1Ytbb70VgMcee4wpU6ZgMpmYNGkSkyZNIjY2lm+++Ybx48efVl7dwSI1t++8807uvvtupkyZwi+//MKjjz5ae4yrfSONleELfKmmMRe4qJHnDwHnSyn7A48Db7dFUGcSabdTtnIVgaNHI/T6eo/RBQRg6tGdyrS0tg1OUTxoz5497Nu3r/Z+WloaycnJDB8+nLS0NNLS0pgyZQqbNm0iKysL0EZSbd26leTk5HrL/Pzzz2uvR44cCUBxcTEJCQkAfPDBBw3GExwcTGlpab3PuVqGt/hMTUNKuVIIkdLI87/VubsWqL+HSmmQZds27IWFTdYgAgYOpPibBUi7vcHkoijtSVlZGXfeeSdFRUUYDAa6du3K22+fft6Zm5vLLbfcQlVVFQDDhg3jjjvuqLfMqqoqhg8fjsPh4NNPPwXg0Ucf5aqrriIhIYERI0Zw6NChel87duxYnn76aQYOHMjf6wx/b04Z3iKklN6OoZYzaSyqr3nqlOPuAXpKKW9uqszU1FSpNmHS5L74IgVvv0P3335FHxbW4HHFCxaQdf8DdFrwDeYePdouQOWMtWvXLnr16uXtMNymZnO3hkZhtSf1/W6EEBullKn1He9LzVMuEUKMBW7iRP9GfcfMFEJsEEJsyMvLa7vgfFzZipX4Dx7UaMKAE53hlZvTPB6ToijtS7tKGkKI/sAc4DIpZUFDx0kp35ZSpkopU6Oj693mts1IKckrrfJqDADWnByqdu1yqXPb2LEj+vBw1a+hKA2omW9xNmo3SUMI0RH4CrheSrnX2/G46rcDBQx/8id+P9hgjmsTNbPAXUkaQgj8Bw6kcosanKYoysl8JmkIIT4F1gA9hBAZQoibhBCzhBCznIf8A4gEXhdCpAkh2kVHRWZhJQ4Jzy7Zgzf7j8pWrMQQH4epWzeXjvcfOJDqQ4ewFxV5NjBFUdoVXxo9dU0Tz98MNNnx7WtKLFYANh4pZPmeXMb1bPvlABxVVZT/9huhUy9zeSHC2n6NLVvUfA1FUWr5TE3jTFVSaUUISI4M4Lkle3E42r62UbFuPbKysllf/nUn+SmKotRQScPDSiw2gk0G7r6gO7uOlbBo27E2j6FsxQqE2UzgiBEuv0ZN8lPONDk5OVx77bV07tyZIUOGMHLkSL7++utmlbFnzx7GjBnDwIED6dWrFzNnzvRIrHPnzq2dZAhw8803s3PnzkZfM2bMGNpieoFKGh5WUmklxN/Ipf3j6dkhmBeW7sFqd7TZ+0spKVuxgsDhw9GZzc16bcDAgVi2bEXa7R6KrvUczdyqUjk7SSmZOnUq5513HgcPHmTjxo189tlnZGRkNKuc//u//+Ouu+4iLS2NXbt2ceedd7o9VrvdflrSmDNnDr1793b7e7WEShoeVmKxEmI2otMJ7pnYg8MFFXyxsXl/qK1RfegQ1vR0gsY0v1/Cf+BAHBUVVO3f74HIWq9kyVL2DB5C5n33UX3kiLfDUXzYzz//jJ+fH7Nmzap9LDk5ufZL3263c++99zJ06FD69+/PW2+9VW85x44dO2m59H79+p12zC+//MJ5553H5ZdfTu/evZk1axYOh3aieNttt5GamkqfPn345z//WfualJQUHnvsMUaNGsWnn37Khg0buO666xg4cCCVlZUn1SIaKqOt+ExH+JmqpNJGsFn7mMf3imFwxzBe+mkflw9KwGz0/BIdZb+4PtT2VP4DBgDaJD9fnBle+PHH6IKDKV36IyXfLSb08qlE33YbRue6PYqP+v4ByN7m3jI79INJTzf49I4dOxg8eHCDz7/77ruEhoayfv16qqqqOPfcc5k4ceJpO/bdddddjBs3jnPOOYeJEycyffp0wuqZLLtu3Tp27txJcnIyF110EV999RVXXnklTzzxBBEREdjtdsaPH8/WrVvp319bPNRsNrN69WpAq1k8//zzpKaePim7sTLagqppeFiJRWueAm3+w70X9iS7xMK8tW1zZlz2yy+YunXDGB/f7Nf68iS/6vR0KtatI3L6NLr+uJTwa6+lZMFC9l80iezHHseak+vtEBUfdvvttzNgwACGDh0KwNKlS/nwww8ZOHAgw4cPp6Cg4KQFDmtMnz6dXbt2cdVVV/HLL78wYsSI2nWq6ho2bBidO3dGr9dzzTXX1CaD+fPnM3jwYAYNGsSOHTtO6qe4+uqrXYq9sTLagqppeFipxUaI2Vh7f2SXSEZ3i+K15fu5emgSwXWeczd7aSkVmzYROX16i17vy5P8ir/+GoQg9LLLMERH0+GhB4mcMZ38N9+icP58ir78ksBzzyXo/PMJOm80xrg4b4es1GikRuApffr04csvv6y9/9prr5Gfn197Ji+l5JVXXuHCCy886XUPPfQQ3333HaCtjAsQHx/PjBkzmDFjBn379mX79u0MGTLkpNedOrRdCMGhQ4d4/vnnWb9+PeHh4UybNg1LnT45V5ZOb6qMtqBqGh5WUlFNXMXJs8HvvbAHhRVW3l3t2dUry3/9DWw2gs4/r8Vl+OIkP2m3U/T1NwSee+5JycAYF0fcvx6ly/eLCbvySqp27yb7n/9k/9hxHJxyGbn/+Q8V69cjHW03EEHxDePGjcNisfDGG2/UPlZRUVF7+8ILL+SNN97AatXmVe3du5fy8nKeeOKJ2qXTAX744YfaY7KzsykoKKhdxryudevWcejQIRwOB59//jmjRo2ipKSEwMBAQkNDycnJ4fvvv28w3oaWTm9OGZ6iahoeZHdIuqfvYPL8OVQO/LR2wlz/xDAm9e3AnFWHuGFkChGBfh55/7KVK9GFhNS+b0v44iS/8rVrsR07Ruy999T7vF9SEh0eeRj58ENUHzxI2YqVlK1YQcH7cyl4Zw6hU6cS//RTbRy14k1CCL755hvuuusunn32WaKjowkMDOSZZ54BtCGthw8fZvDgwUgpiY6O5ptvvjmtnKVLl/LXv/4Vs3Mk4nPPPUeHDh1OO27kyJE88MADbNu2rbZTXKfTMWjQIPr06UPnzp0599xzG4x32rRpzJo1C39//9pdAQEGDBjgchkeI6U8oy9DhgyR3lJUXi3/78r75c4ePWXm3x886bmdWcUy+f5F8uO1Rzzy3g6HQ+4ZNUqmz57dqnLs5eVyZ+8+MufFF90UWetl3P03uXvYcGm3WJr1Oltpqcx++hm5s0dPWfTtIg9Fp9Rn586d3g6hzSxfvlxOnjzZ22G4rL7fDbBBNvCdqpqnPKjEYiWysli7/f332MvKap/rHhuMXic4Vlzpkfeu2rULe14+QaNb3jQFvjfJz15cTOmPPxI6eTI6k6lZr9UHBRHzt7vxHzSI7H/9C2tmpoeiVJQzl0oaHlRcaSXKUow0GJGVlZQsXlz7nF4niA4ycazYM51YZStXARA0elSry/KlSX4lixcjq6sJ/cMVLXq9MBiIf/YZcDjIvP9+n/iZlDPLmDFjWLRokbfD8BiVNDyoxGIlqrIYe4/emLp1o6jO6A2A2FAzOSUeShqrVmHu3RuDG/YT8aVJfkVffoWpRw/MrZgd65eUROwjD1O5YSMF78xxY3SKcuZTScODSiptRFqK0cfEEHblH7Bs2Yplz4mtQOJCzGR7oKZhLy6mcvNmAs8b7ZbyfGUnP8uevVi2byfsD1e4vFpvQ0Ivu4yQiyeR9+qrVG5z80QzL5EOB+Xr1lH++zqqjxxRS6woHqFGT3lQaWU13SqL8YvrQMiUKeQ+/x+KvvyCDg8+CECHUDO/7s93+/uW//YbOBwEneee0U7GpCRtkt+WLYT/ybUJSJ5Q/NVXYDQScumlrS5LCEGHf/6Tis1pZN1zL52++hKdC+PkfZGjupqShQspeO99qg8ePOk5fVgYhg4dMMbG4j9wACGXTsEvUc2YV1pOJQ0PqigoxOSw4R/XAUN4OEETxlOyYCEx99yDzs+P2BAzpVU2yqtsBJrc96soW7ESXWgo/gPcs7SAEAL/QYOo3LTJLeW1hKyupvjbbwkeOxZDeLhbytSHhhL/9NMcnTaNnKefJu7xx91Sbluxl5RQ+PnnFH74Eba8PEy9exH/3LMYoqKwZudgy8nGmp2NLTsH67FjlL20gryXXiZg6FBCp04l+MIL0Qe1z0SpeI9qnvKg6uwcAAITtAloYVdeib24mLKffgKgQ6g2+ifbjf0a0uGgbPVqgs49F6F339pWAUOHUn3kCNacHLeV2RylK1ZgP36c0Csud2u5gcOHEXnzzRT97wtKfvzRrWV7iu34cXKefY79Y8eR958XMHXrRsf33qXTl18SeumlBI4cSdjlU4maNYu4Rx8l6c036LzgG7ou+4no2X/FlpvLsYceYt+oUWTeex/la9d6dVfJtuKOpdEB9u3bxyWXXEKXLl0YMmQIY8eOZeXKlR6IuGHeXDrdZ5KGEOI9IUSuEGJ7A88LIcTLQoj9QoitQoiGVx/zEY487QvWzzn5J3DkSIzx8RR9oXWIdwjxB3Brv4Zl5y7s+flu68+oEThiOAAVa9e6tVxXFX/1NYboaIJGtX402Kmi77wDc58+HLv/AZ9cMqWGo6KC/Dfe4MAFEzk+dy5BY8bQ6asv6fjeuwSec06T/TzGhASiZs2i8w/fk/LZp4ROvYyyFSs4Om06h6+8ipLvvz/jRpNJKZE2G3aLhalTpjBq2DD2btzI70uWMO/1N0g/fLhZ5VksFiZPnszMmTM5cOAAGzdu5JVXXuHgKc2CnuTtpdN9JmkAc4GLGnl+EtDNeZkJvNHIsT5B5Gv9FcYO2havQqcj9A9XUP7bb1RnZNAhVJtV6s6kUb5KO+Nx95erqUcP9GFhlK/93a3lusKWl0fZypXadrUG97eoCj8/El9/HX1kJEdn3nrSYAVfIG02Cj+fz4ELLyLvpZcJPGcknRctIuE/z7doFFnNmmJxjz5Kt1Ur6fD4YzjKy8m8624OTLqYws8+85lOdCklDosF2/HjVGdkYNm3j6p9+6g6eJCqI0eoTk+nOisLa3Y21mPHqE5Pp+rQIar27ceyezeWHTux7N7NknnzMNjtTJ8wAWtmJtbsbOKNBm654AKq9u3DkpXFPbPvanBp9JrkM2/ePEaOHMmUKVNqn+vbty/Tpk0DoLy8nBkzZjB06FAGDRrEggULTvuZ3L10enlJCWPOP5/1a9ci7XZmzZrl0aXTfaZPQ0q5UgiR0sghlwEfOmcrrhVChAkh4qSUbb8Vnov0x/NwIE4a9hp2+eXkv/oaxV99RYdZtwPubZ4qW7kKc9++GKKi3FYmaAkvYPjw2qaM1o5ecpWUkoJ33wO7ndDLWzY3wxXG2Bg6vv8eR669jqM33UTKvI/wS0nx2Pu5QkpJ2bJl5L7wX6oPHsR/0CASXnqRgEaW+G4unclE+FVXEXbFFZQuW0bBnHfJfvRf5L3yKuHXXUvYZZe5fan5Z9Y9w+7juxs+wOHQajx2u7ZOWE3TmRAInQ6E0JrTah6Xku7+yfwteRoYDAi9AeFnRBj8tdsGPXvy8xk8fDimzp1Br0cYDEi7HUdJKfbSEua89TaBdhur5s3DajAy5vKpnN+rF53i4pF2G9JmA2DrqlX0S0yk6tAhhJ8fwmhE+PmhM5kQJhNPPPEE48aN47333qOoqIhhw4YxYcKE0xYjbMnS6SaTiZVLl+KorOSdfv148t57Gdy9Oxw9qg2JP3IES1AQj1x/PRF33IFDCCbffrvbl073maThggQgvc79DOdjpyUNIcRMtNoIHTt2bJPg6mMqLKA8IBhhPLGSrTE+nsBRoyj66muibr+dUH+j22oa9qIiKrdsIWrWrW4p71SBI4ZTumQJ1qNH8UtO9sh71CXtdnKeeprCefMIvfxyTJ07Nf2iVvBLTKTje+9y5M/Xc2TGDFI+/rjNV8e1ZmVRvmYN5WvWUr52Lfb8fPw6dSLx1VcIGj/eY8la6PWETJxI8AUXULFuPQXvziH/5VfIf/kV/AcOJGTyZEIuutAt837qU5sk7HaoWVBSp9P65fR6LVnoTjSMnPopGCIiG6116QMC0Pn5oQsIALSl0VevXo2fnx/r16/nl+3b2LplC9/8/DPY7RSXlbH/4CE6J6doycdgQOj16Pz9tf9nhwNHaSl/vP12Dhw9StfkZD578UWWLPyWBV9+yXNPP43Q6bBUVnJw+3Z6de6MtNmQVivVmVmk9utHgsOB7cgRrrzgAlYsXsyU4cP5ZN483v3kE2w2G9m5uWxZvpweAQFIm42pgwZT5WwGk3Y7OqMRQ0yMFpufH4bISIwdOvDN998z56OPsNntZOfmsnPnzrM2adT331Jv752U8m3gbYDU1FSv9fAFlBynIiTitMfD/vAHMmfPpvzXX+kQYnZbTaPs11+dQ21bt3RIQwKce4yXr1nr8aThqKoi6977KF26lIjp04lpYHFCdzN17UrSnDkcnTaNo9NnkPzxPAyRkR57PyklFWvXUrJ0KRW/randgVAfGUngyJEEnX8eIZMmeaRZrj5CCAKHDyNw+DCq09MpWfw9JYsXk/PEE+Q89RQBw4cRfMEF6Pz8sJeU4igr1a5LS7CXluGXlETA0FT8hwxpcJTbfUPvQ1ZV4aisxFFejqO0VEsWQqALDEQfEoIuOBid0X3bBjS5NDrwymuvNbk0er+hQ1m5ciWmLl0AWPDjj6xfs4Z7778fQ3Q0UsCnL71Etzq7+wFYc3K0xGcwgk6g0+sRRqNWk3JIhJQc3LuX/775Jqu/+IKI8HBuvv9+LNXVWm1GCEISE/BLTET4+6MLCMAYF4cxJgYAYTRiCAsjvbSU/771lkeXTm9PSSMDSKpzPxHIauBYnxBUWoglJva0x4PHjUUfHk7RF18SO/A6t80KL1+5En1YGOZ6tqB0B7+UFAyxsZT/vtaj8zXsRUWk33EHlRs3Efv3B4i48UaPvVd9/Pv2IemtNzl6080cvfkWkj+Yiz4kxK3v4aispHjhtxTO+4iqffsRAQEEDh1K+LXXEDBiJKbu3dqsCbAhfklJRN06k6hbZ1K1fz8lixdT/N135DxWZ2iyEOiCg9EHB6MLCKD81185/sEHAJi6dcU/NZWA1FQccXFYs7NxVFYiKytrl6cXev2J1wcFuXXEX13jxo3jwQcf5I033uC2224D6l8afdy4cRiNRvbu3UtCQgJPPPEETzzxRO1x1157LU899RQLFy5kypQpWm3CbkcYDBhjY7no0kt5e9EiXn7xRaiuZtPmzQxOTa2tqQD4ZWWyfssWMh0OkpOT+Wr5z8ycOZOqqCiCwsKIHTKEvLw8lq5ezfhLL9VO0PR6jDExtX+HzVk6fcyYMW79LNtT0lgI3CGE+AwYDhT7cn8GQGh5EbkRfU97XPj5EXrZZRyfN49Og69kcfHpO381l3Q4KFu1mkA3D7WtSwhB4IgRlK1ciXQ4tDMnN7NmZXH0lplYjx4l4YX/EDJpktvfwxUBQ4aQ+MorpP/lLxy66irMvXpjiAhHHxGJITICfUQkxvg4zH37NuvL3ZqdTeHHn1A0fz724mJMvXoR9+SThEy+uNkLMLYlU9euRP/f/xF1551Y09MRBgO6kBB0AQEn/R04qquxbN9OxfoNVGzYQMnCbyn69DPsr72KzWRCZzajDwvTzpb9/REmU5skR3ctje7v78+iRYu4++67mT17NrGxsQQHB/Pwww8D8MgjjzB79mwGDBqElJKUlJR616Fq10unN7T8bVtfgE/R+iesaLWKm4BZwCzn8wJ4DTgAbANSXSnXW0ujW8sr5M4ePeU3dz9e7/OWffvkzh495Zf3PyNTHlgkq232Vr1fxdat2pLfCxa0qpymFH71tdzZo6es3L3b7WVX7t4t944+T+5OHSrLfv/d7eW3RMnPP8vDf75e7p90sdwzbLjc2aPnSZfD06ZJy/79TZZjOXBQZvztHrmzdx+5s1dvmX7HnbJ83TrpcDja4KfwHofVKiu2b5c70tKkw966v/Ezha8tnd7cpdF9pqYhpbymieclcHsbhdNqJZlaJUg00HFo6toVc//+pGz4BTmkN3mlVcSH+bf4/cpWrgQhCPTAPIa6auZrlK9Zg7lHD7eVK6Uk6557QAiSP56HuXt3t5XdGsFjxxI8dmztfWm1YissxH78OBXr1pH3yqscvGwqkdNuJOq2205bisSamUnea69T/M03CLOZiBtuIPy6686apTyEwYB/nz6IXbs8UjNV2p76LXpISbrW3WKop0+jRujUy/DPOESX4qxWL5FevnIV5n79MESc3vHuTsa4OPySk6lw83yNqr17qdq3n6jbbvOZhFEfYTRijInB3LMnETfcQJcfvid0yhQK5rzLgYsna0u3S4k1N5fsxx5n/0WTKFm0iIjr/0zXH5cSe/99Z03CUOrX3pdOV0nDQyqcNQ2/+NO3gqwRevHFSKORCUfXt6oz3FZYSOXWrQSNdu8s8IYEjBih7bXtHLvuDiXfLQa9nuALJ7qtzLZgiIwk/sknSP70E/SREWTe/TcOX3kVByZeSOH8+YRdfjldli4h9u9/9+goLEVpKyppeIglOxuAgPiGx/nrw8IwnzeGsRmbySk4fSSEq8pXrQIpCXLz0iENCRw5Akd5OZbt9a740mxSSkoWLybwnHPcthhhWwsYNIhO//sfsf94BHtpKcETL6DL4u+Ie+xfGOvZQ1pR2iuVNDzElp1DhcFESGRYo8dFX3kFodXlsG5No8c1RNrtFLwzB2PHjpj7nj5SyxMChg0DoNxN61BZtm3DmpHhtZFS7iL0eiKuvZauS5eQ8Oyz+HlxYqmieIpKGh7iyMulwBxCiLnxCUpBo0dR7B9CzG8/teh9ir9ZQNW+fcTcfZfHhtqeyhARgalnT7etQ1Xy3WKE0UjwhPFuKU9RFM9RScNDdAX5FJhDCfFvfICaMBjY0WsEHfduxlZY2Kz3cFRWkvfSS5gH9Cf4lJmsnhY4fDiVmzbhqGrdHBPpcFDyww8Ennee2yfQKUqNgoICBg4cyMCBA+nQoQMJCQm196urq0869sUXXzxp4l9D3LXUeHujkoaHGI7nk+8fSpALmyulDxuHwWGnZNF3zXqP4x98iC03l9h7723z2cMBI0cgq6up3Ly5VeVUbtqELSeHkIvbd9OU4tsiIyNJS0sjLS2NWbNmcdddd9Xe9/PzO+lYV5PG2UolDQ+QDgem4uOUBIZj0Df9Eft1786BsASKmrEhjO34cQreeYeg8eMJcK6f05YCUlNBr291v0bJ4sUIs5lgNy91oChNWbZsGYMGDaJfv37MmDGDqqoqXn75ZbKyshg7dixjnfNzGlqu/GzlM5P7ziT2ggJ0DjsVoa6NBIoNMfNjUipdti3Asmcv5h5Nz1PIf+11HBYLMX+7u7Xhtog+KAj/fv2oWLMWZresDGmzUfLDEoLGjmm3+3MrzZf95JNU7WpkafQWMPXqSYcHH3T5eIvFwrRp01i2bBndu3fnhhtu4I033mD27Nm88MILLF++nCjn9gKNLVd+NlI1DQ+w5uQCUBXq2p4WcaFmfkkcBHo9xfWsd3Oq6sOHKfz8c8KuulLbH8BLAkYMp3L7duxlZS16ffnvv2M/fpyQiy92c2SK0ji73U6nTp3o7pxIeuONNza4Zev8+fMZPHgwgwYNYseOHU1uq3qmUzUND7Dlatu8WiNdSxqxIWaKTUFYhoyk+Ntvifnb3Y0uhZ373xcRfn5E3+7dVVUCR4yg4M23qFi//qSlNlxVsngxusBAjy3lrvim5tQIPOXUTZEacujQIZ5//nmPLjXe3qiahgdYnRP7ZKRrG9bUbPt6bOR47Pn5lP/6a4PHVmzeTOmSJUTeNMNjG+K4yn/QIISfX4uWFJHV1ZT++BPBE8b79OquypnJYrFw+PBh9u/fD8BHH33E+eefD5y87Hh9S42f7VRNwwNsObnYhc7lZSNigk0IAXtT+tE1LIyir78hyPkHXJeUktznnkcfHUWkc09ib9KZTPgPHtyizvCyX3/FUVKimqYUrzCbzbz//vtcddVV2Gw2hg4dyqxZswCYOXMmkyZNIi4ujuXLl3t+qfF2RiUND7Dl5FDkH0JwgF/TBwNGvY6oIBPHyu2EXHIJRZ9/jr24GH1oKAD2khKsmZmU/7aGyk2b6PCvf/lMx3HgOeeQ98ILFH76KWF/+pPLQ39LFn+PPjSUwJEjPRyhopzs0Ucfrb29uZ4h43feeSd33nln7f25c+fWW84vv/zi5sjaB5U0PMCak0OeKYQQf9e3q6zZ9jX08qkUzpvH0ZtvQVqtWDMzcdTZocvUuxdhf7jCE2G3SPi111KxcQPZ/3qMyh076PCPf6DzazxZOiorKVu2jJDJkxFNHKsoim9RScMDqnNyyPcPbXIJkbpiQ8xkFFZg7j2UoAnjsR45gjE+gYDBgzEmJNReTN27tdl+0a7QBwWS9Prr5L36KgVvvEnVvn0kvvwyxtiGl4QvW7ESR0WFmtCnKO2Qz3z7CCEuAl4C9MAcKeXTpzwfCswDOqLF/byU8v02D9QFtuwc8jsMopPZ9Y83LtTM+sPHEUKQ9OqrHozO/YReT8xf/4q5Vy+OPfB3Dv3hShJf/G+9kw6l3U7xom/RR0XVLnyonB2klF7f91w5mba3XfP4RNIQQujRtnK9AG2r1/VCiIVSyroDom8HdkopLxVCRAN7hBAfSymr6ynSa+xl5VBRToF/KAOa0zwVaqa40orFasdsbJuFB90tZOJETJ07k3H7HRyZNp2Yu+/G2CGWqgMHqTp4gOoDB6k+dAhptRJ+w/VttsCi4n1ms5mCggIiIyNV4vARUkoKCgowm83Nep1PJA1gGLBfSnkQQAjxGXAZUDdpSCBYaH9xQcBxwH27ALlJzRyNAnPzm6cAsostpET5Rid3S5i6diXlf/PJuu9+cp99VntQCIxJSZg6dyZw9ChMXboScslk7waqtKnExEQyMjLIy8vzdihKHWazmcTExGa9xleSRgKQXud+BjD8lGNeBRYCWUAwcLWU0tE24bnOlqMljXz/ple4rSvOOVcju8S7SSOvtAqTUdeshHcqfUgIia+/RsX6DejDQvFLSVFzMc5yRqORTp06eTsMxQ18JWnUV189tbHtQiANGAd0AX4UQqySUpacVpgQM4GZAB3beCMca7YzabSipuEtBWVVTHppJXqd4PXrBjMkueX7jQudjsDhqs9CUc40vjIjPANIqnM/Ea1GUdd04Cup2Q8cAnrWV5iU8m0pZaqUMjW6jWdN19Q0CvxDmzfktk5NwxuklDz8zXaKK634GXT86e21fLjmcIs6yhRFOXP5StJYD3QTQnQSQvgBf0JriqrrKDAeQAgRC/QADrZplC6w5eZgDQiiWm8kuBmjp4JMBoJNBq/VNBZtPcb327OZPaE7i+4YzXndovnHgh38bf4WKqvtXolJURTf4xNJQ0ppA+4AlgC7gPlSyh1CiFlCiFnOwx4HzhFCbAOWAfdLKfO9E3HDrDm5VIREEOCnx+jCXhp1xYaavZI0ckstPLJgOwOSwrj1vM6EBhh554ZU7prQna/TMrnijd84WqA2pVEUxXf6NJBSLgYWn/LYm3VuZwET2zqu5rLl5FAaHNGijuSaWeFtSUrJQ19vp6Lazn+u6l+7aZROJ/jrhG70Twpl9mdpXPLKKl65djDnd/fuIomKoniXT9Q0ziTWnGyKA0Ob1TRVo0OomZw2ThoL0rL4cWcO90zsTteY4NOeH9sjhm/vGEV8mD93fLJJNVUpyllOJQ03klYr9vyCZneC1+gQYia3tAq7o206n3NKLPxz4Q4GdwzjplENb+bUMTKAf1zam1KLjSU7stskNkVRfJNKGm5ky88HKckzhRLSgppGbKgZu0OSX1blgehOJqXkwa+2YbHaef6qAeh1jc/SHdEpkqQIfz5fn97ocYqinNlU0nCjmuG2x/yCW1TTiGvDuRpfbspk2e5c7ruoJ52jg5o8XqcT/HFIEmsOFqhOcUU5i6mk4UY1E/sy9EEt6wiv2cHPw0mjrMrGv77dwbCUCKafk+Ly6/4wJBEh4IuNqrahKGcrlTTcqGbdqaO6wGYtIVKjZla4pzvDf9qZQ6nFxr0X9UDXRLNUXfFh/pzXLZr/bcxos36XxhwrrmTJjmyqbKpzXlHais8MuT0TWHNywM+PQkNAi2oakYF+GPXC48NuF209RocQM0M6hjf7tX9MTeL2Tzaxen++14bfbs8sZs6qgyzaegybQ5IcGcDDk3szoVeMWkFVUTxM1TTcyJaTiy4qGoRoUZ+GTieICfbsBL8Si5WVe/O4uF9cs2oZNSb0jiE8wMj8Nu4Ql1KyfE8u176zlkteWc2PO3O4YWQKr107GKNexy0fbuCG99axP7e06cIURWkxVdNwI1t2No4o7ey7pavEdvDwrPCfduZQbXcwuX9ci15vMuiZOiiBeWuPcLy8mohAz27XKqXk263HeGXZPvblltEhxMwDk3pyzbCOhDoT88Q+sXy05gj//WkvF764ihtGJjN7Qvfa5xVFcR+VNNzImpuLrVN3gBZN7gNtrsauY6ct3Os2i7YeIyHMn8Edw1pcxtVDk3j/18N8szmTGaM8t9z15qOFPL5oJ5uOFtEjNpgX/jiAS/rH42c4uYJs1OuYMaoTlw2M5/mle5n722EWpGVxy+jOXDeiY6uWeXcHu0NyIK+MLelFlFhsCEAnQAiBENoSz/Fh/ozqFoXJ0P42piqrsvH15kwWbM6ke4dg7hjblfgwf2+HpXiIShpuIqXElpODZdAIgBY1T4FW01i+J9cjW2MWV1hZtS+PaeektKrsnh1C6J8YyvwN6Uw/t3Vl1SezqJJnf9jNgrQsooNNPPuH/vxhSGKTc0kig0w8dUU/rhvekWd+2M0zP+zmteX7uW54R2aM6lQ70KA1pJQcKaggLb2I/bllBJj0hPobCfU3EubvR6i/kQCTnn05ZWzJKCLtaBHbMospq2p6v7AQs4FJfeO4bGA8wztHnvbzSinJLKpka0YxB3LLiAwykRThT2J4APFh5kYTjpSSY8UW9ueWaZe8strbdodkaEoEI7tEMqJzBL06hLjUdLk7u4R5a4/w9aZMyqvtdIsJ4n8b0vliQwbXDEvi9rFdiWnkM88utrAnp5TiSivFlVZKnNfFFVYqrXYiAv2IDjYRHWQiKtiP6CAzkUF+SMBitVNldVBls1Nlc1Blc6AXApNRh9mgx2TUYTLoMBn0BJr0BJkMzfo7lVJSUmkjo6iCzMJKMgoryS6x4KfX1f6+Q/yNhAVotw06gdUusTkcWO0O7bZdYm9glWidgAA/A0EmA4EmPYF+BgJNBvwMOu27xCG1cmySartWpl4nMOp1+Bl0GPUCo07Xoibm1lJJw00cxcXIqioqQiKglBZN7gOtplFRbafEYnN788rSndlY7ZLJ/eNbXdYfU5N4+JvtbMsspn9iWOuDA8qrbLy54gBvr9QWL75jbFdmjelCkKl5n2XfhFA+umk42zOLeWvlQd5ZdZD3fj3E5YMSmHleF7rGND0vpUZ+WRXbMovZkl5EWnoRW9KLKKywAiAENLZyvFEv6BUXwuWDEhiYFMaApFCig8xIJFKCQ0ok2vWOrBK+Tcti0dYsPt+QTkywicn940hNjmBPdglbM4vZllFMQXn9uxsLATHBJhLC/BFCUF5lo9Jqp6LaToXzdt0Bb2EBRrpGBzGxdyxSwtpDBfy0Sxv9F+pvZHinCAZ1DK9deNOgF/g5r8urbHyxMYP1hwsxGXRcOiCe60ckMyApjMyiSl79eR8f/36Uz9anc/2IZGaN6UJUkInsYgtrDxbUXg7XM9/HbNS+lM1GPcfLqym1uGdzTqNeEBbgR3iAkfAAP8ID/Ag0Gai2O6iy1iQe7bq8ysaxIgulpyR6k0GHzSE9OnLQoBPYpWz07+rU40P8jUQHmYgJ0RJstPM6LtS/xc3Qjb6n20s8S1lzcgEoCXYmjRZ+4ceGnhh26+6ksWjrMRLD/RmQGNrqsi4dEM/ji3Yyf0O6W5LGkYJybnxvHYcLKpgyIJ77J/UkoZVNHH0TQnnlmkHcO7EHc1Yf5PP16czfkEGX6EA6RQWSHBlISlQgKZEBpEQGYnNIdmaVsPNYsfO6hJwSbXa+ENAtJogLescyMCmcgUlhdI8NwuaQFFVYa8+WiyqqKauykRIVSO+4EJf3e4/pYWZsjxgsVjvLduWycEsmH689yvu/HkYnoHtsMON7xdAvMYz+CaF0jw3meEU1GccryHCeCacXVpBVVIlOCCIC/Qjw0zsvBgL89MSGmOkaE0TXmCAiA/1OO/POKqqs/UJfc7CApTtzGow3JTKAhy7uxZVDEgmv06+VEObPU1f057bzu/Lyz/t479dDfPz7UWJCTBxxJolgs4HhnSL584hk+ieGER5w4sz91M/LYrWTX1ZFXql2KSivrq1RmE6qUehwyBM1EIvtxHWZxUZhhfa7OV5eTVGFlQN5ZVRU2zEZtDN3k1GPyaAjyGQgOsjEOV2iSAjzJzHcn4RwrTYXHqD9P5ZV2Zy/6xO1I5tDamf/eh0GvQ6jTmDQ69DW/zy9NuCQkvIqG+VVdu262kZ5lY2KajsGnXDWJpwXg1aeXUqsNq0WU213UG1zUG13UFxpJa+0itzSKg7mlZNXWkW13VF74uFuKmm4iS1HW5OpyF/7Qm5pn0bttq/FFrrHnr6AYEsVllfz6/58bhrdyS3NSaH+Ri7uF8eCtCwentzb5S/H+mzLKGb63HXYHJJPbxnByC6RrY6vro6RATx2WV/+Or4bn647yrbMYo4UVLB6fz4W6+k7Bht0gq4xQZzbJYre8SH0jg+hX0IowfX0jRj00CFUXzsxs7XMRj2T+8cxuX8cJRYrh/LK6R4bjL/f6Z9vgp8/CWH+p+2L3FLxYf5cMTiRKwZre0aXWqxU2xzYHLL22mbXPq8u0UGNNo10jAzg+asG8JcxXXj9lwMUVVi5fkQyIzpH0isupMmmxhpmo57E8AASwwNa/wO6SbDZSLDZSGLzR6y3CSll7UmMJ6ik4SbWOnuDm40VLe7Q7OChpUSW7szG5pBc0q/1TVM1rkpN5OvNmfywPZupgxJaVMaKvXncNm8j4QF+fDZjWLOajporMsjEHeO61d53OCS5pVUcyi/ncEE5ep2gd1wI3WKDfKJDOsRsZEBSmNfev74k2Vydo4N4/qoBbohGcZUQWlNcWIBnRjaqpOEm1ixtd9o8YzAh5pZn+JgQE+D+bV8XbT1GcmQAfRNC3FbmiE6RdIwI4LP1R7lsYHyzazBfbcrgvi+20i02mLnTh7qlo7o5dDpBh1AzHULNbq/dKMqZSk3ucwNHZSXFX3yJf+oQim0t788AbR5EVJAfR4+7b1HAgrIqfjtQwOR+cW4d6aTTCf40LIm1B48z/oUVvLZ8P1lFlU2+TkrJG78c4O75WxjWKYLPbx3R5glDUZSW8ZmahhDiIuAlQA/MkVI+Xc8xY4AXASOQL6U8vw1DbFDhJ59gy8sj4b8vULLF2uKRUzUGdwxnzYECtw27XbIjB7tDeqRTbObozkQG+vHlxkyeW7KH55fu4ZwukfxhcCIX9e2AXifILdE66fJKLeSWVrHxSCEL0rK4dEA8z1/V3yeaghRFcY1PJA0hhB54DbgAyADWCyEWSil31jkmDHgduEhKeVQIEeOpeKSUFH32Gf5DhmDu3r3RY+2lpRS8/Q6Bo0cTkJpKydrVJ40maYnR3aNZujOHwwUVdIoKbFVZAN9ty6KzczSPuxn0Oq4e2pGrh3bkaEEFX23O4KtNmdw9fwv3fbEVWz3DEw06wa3ndeb+i3p6ZZy5oigt5xNJAxgG7JdSHgQQQnwGXAbsrHPMtcBXUsqjAFLKXE8F4ygpIe+119GHhdLpf/9D59/w0M+C997DXlxMzF2zASix2OgY2bov+tFdowBYtS+v1Ukjv6yKNQcKuH1sV48v5tcxMoDZE7rz1/HdWH+4kJ935xLopycmxERMsLn2OiLQz+XRM4qi+BZfSRoJQN0V8DLgtJGE3QGjEOIXIBh4SUr5YX2FCSFmAjMBOnbs2Oxg9KGhJDz7DEdvupmcJ58i7vHH6j3OVlDA8Q8+JHjSRZh79wagpLL1zVPJkQEkhvuzal8+N4xMaVVZ32/PxiHxSNNUQ4QQDOsUwbBOEW32noqitA1f6Qiv77Tz1HYNAzAEmAxcCDwihKi37UhK+baUMlVKmRod3bLluwPPOYfIW26h6H//o2Tx4nqPyX/rLWRVFdF3/l/N+1JisbaqIxy0L93R3aJZe6AAq/30eQTN8d3WLLpEB9LDjXM+FEU5e/lK0sgAkurcTwSy6jnmBylluZQyH1gJeHQAePSdd+A/cCDH/vFPqtNPXgrcmpVF0aefEXr5VEydtUX7LFZttqY7Fsgb3S2K0iobW9KLWlxGbqmF3w8d55L+zR8OqyiKUh9fSRrrgW5CiE5CCD/gT8DCU45ZAIwWQhiEEAFozVe7PBmUMBpJ+M/zoNORefffkNUn1v3Je+01AKL/8pfax0os2vyMluzad6pzukSiE7BqX36Ly5i35ghSwpSB7pvQpyjK2c0nkoaU0gbcASxBSwTzpZQ7hBCzhBCznMfsAn4AtgLr0Iblbvd0bMaEBOL+/TiWbdvIffElAKoOHqL4628Iv/YajPEnvpBLnNP23VHTCAvwo19iGKv25bXo9SUWK+//dpiL+nSgS7TnZlkrinJ2afYpsRAiELBIKd26MbOUcjGw+JTH3jzl/nPAc+58X1eETJxIxbXXcPy99wgcMZyir75GZzYTeeutJx1X4lyRs7V9GjVGd43ijRUHtH6SZiaij9YcodRi445xXd0Si6IoCrhQ0xBC6IQQ1wohvhNC5AK7gWNCiB1CiOeEEN2aKuNMEHP//Zh69CDznnsp/eEHIqbdiCHi5NFBtc1TrRw9VWN0tyjsDsmaAwXNel1FtY13Vx9iTI9o+ia0fkVbRVGUGq40Ty0HugB/BzpIKZOklDHAaGAt8LQQ4s8ejNEn6EwmEv77AtJqRR8aSsT06acdU9s85aaaRs1+Bs1tovrk96McL6/mTlXLUBTFzVw5JZ4gpTxtBT4p5XHgS+BLIcRZsRmzqXNnkue+D4A++PQhrDXNUy1dFv1UfgYdIztHsroZneEWq513Vh1kZOdIhiSreRKKorhXkzWNmoQhhPj3qc85l/+gvqRypvIfMAD/AfWP9HVnR3iNUd2iOFxQQbqLCxh+sTGDnJIq1ZehKIpHNGf0VIIQ4pqaO861n35yf0jtV4nFip9B16oNiU41ulvNkiJN1zasdgdvrjjAoI5hnKOW+lYUxQOakzRuBWYKIYYJIYYCPwPPeyas9qmk0ubWWgZoO6TFhZpZvb/pfo0FaVlkFFZyRxusM6UoytmpycZ3IcSHwCZgM3A78AlgA6ZKKfd7Nrz2RVtCxL3LeQkhGNU1iqU7teXNG1roz+6QvP7LfnrFhTCup8cWAFYU5SznSk3jA+dxM9ASRgpQCPxZCHGl50Jrf7TFCt0/JmB092iKK61syyxu8Jjvtx/jYF65qmUoiuJRTZ4WSymXActq7gshDEBvtHWfRgBfeCy6dqbUYnPbcNu6znX2T6zel8fAevaMllLy6s/76RIdyEV9O7j9/RVFUWq4MrnvpNNWKaVNSrlVSvmRlPKe+o45W2kzt92/2nxkkIk+8SGsrKcz3O6QvLPqILuzS7l9bFe1T4WiKB7l0uQ+IcSdQoiTNqYQQvgJIcYJIT4AbvRMeO1LSaVnahoAo7tFs/loIWVVttrHVuzN4+KXVvHk4t2c2zWSKQPUwoSKoniWK6fFF6H1Z3wqhOgEFAH+aAlnKfBfKWWapwJsT0osVrdN7DvV6G5RvLniAL8fLCAxPIAnFu9i5d48OkYE8Pp1g5nUt4Pqy1AUxeNc6dOwoO3N/bpz5ncUUCmlLPJwbO2KxWqn2uZwuSNcSsny9OWMiBtBgDGgyeOHJIdjNur417c7ySisIMhk4OHJvbh+ZDImg/vmhSiKojTGlT6NZUKIPlA783socIcQYping2tPTuyl4VrSWHRwEX9d/lfe3vq2S8ebjXpGdY0mq6iSG89JYcW9Y7l5dGeVMBRFaVOutKUkSil3AAghzgE+Aj4H5gohHpJSfu3JANuLkkrnsuguNE8VWYp4br22wvuX+75k1oBZmA3mJl/3nz8OoMpqJyak6WMVRVE8wZWO8JI6t28A3pRSzgTGAPd7Iqj2qDk1jec3PE9pdSkPDHuAoqoivj/0vUvvEepvVAlDURSvciVp7BdCXOlca2oq2rarSClzAZO7AhFCXCSE2COE2C+EeKCR44YKIey+NrHQ1cUK1x1bx4IDC7ixz41c2/NauoZ15ZPdnyClbIswFUVRWsWVpHEX2rpTmcAmKeVvAM5OcbfsI+pcLfc1YBLaxMFrhBC9GzjuGbRtYX1KqXNZ9NBGlhGpslfx2NrHSApOYtaAWQghuLbXtew+vpvNuZvbKlRFUZQWc2Vp9Gwp5QWASUp5cZ2nxqJt0OQOw4D9UsqDUspq4DPgsnqOuxNtD49cN72v25zYta/hmsbbW9/mSMkRHhnxSG0fxuROkwnxC+HjXR+3SZyKoiit4fIqt1JKxyn3lzr7NtwhAUivcz/D+VgtIUQCcDlw0r7hvqK2I7yBPo0DRQd4b/t7XNr5UkbGj6x9PMAYwBXdrmDZ0WVkl2e3SayKoigt1Zyl0T2pvllppzbyvwjcL6W0N1mYEDOFEBuEEBvy8pq3VWpLHS+vwk+vw2Q4/SN1SAf/WvMvgoxB3DP0ntOev7rH1Tikg/l75rdFqIqiKC3mK0kjA0iqcz8RyDrlmFTgMyHEYeBKtMmGU+srTEr5tpQyVUqZGh0d7YFwT7c9s4SeccH1zsr+Yu8XbM7dzD2p9xBhPn0L1sTgRM5POp8v9n5Blb2qLcJVFEVpEV9JGuuBbkKITkIIP+BPwMK6B0gpO0kpU6SUKWgr6/5FSvlNm0daD7tDsjWjqN4VaPMq8nhx44sM7zCcKV2mNFjGdb2uo7CqkB8O/eDBSBVFUVrHJ5KGlNIG3IE2KmoXMF9KuUMIMUsIMcu70TVtX24p5dV2BnUMO+lxKSWP/PYI1Y5qHhn5SKNrQw3vMJwuoV3U8FtFUXyaTyQNACnlYilldyllFynlE87H3pRSntbxLaWcJqX0mX08Nh8tAmBQUvhJj3+y+xN+zfyVe1LvITkkudEyhBBc0/MadhbsZEveFk+FqiiK0io+kzTas7SjRYQHGEmOPLHw4P7C/byw4QXOSzyPq3tc7VI5l3a5lGBjMJ/s+sRToTaqtLqUHw79wLKjyyipLmn6BYqinHU8s473WWZzeiEDk8Jqm5+q7dXcv+p+gvyCeOycx1xesjzAGMDUblP5dNen5FbkEhPg+b2+iyxFLE9fzo9HfmTNsTXYHNrQYZ3Q0TeqLyPjRjIyfiT9o/tj1DU8B0VKSbWjGovNol3sFuwOO0khSY2+rkZpdSkbsjcAMChmEGHmMLf8fHaHnazyLI6WHOVwyWHSS9MJMgbRObQzncM6kxKS4tK6X23B7rCTW5FLVnkWWWXapdpRzYi4EQyKGYRB1z7/Xa0OKw7pwKR32wISihe1z79CH1JisbIvt4xL+p/YAOnlTS+zt3Avr41/jUj/yGaVd02Pa5i3cx7Prn+Wx899HH+Dv7tDJr8yn2VHlvHj0R/ZkL0Bu7QTHxjPtT2v5YLkC7BLO2uy1rAmaw3vbHuHt7a+RYAhgM6hnbE6rFQ7qqm2Oy/O2xabBXnaKGkw6830j+7P4NjBDI4ZzIDoAQQYA7A5bOwo2MFvWb+xJmsNW/O2Yq8zmrpzaOfa1wyOHUx8YDxCCKSU2KUdq8OKzWHDYrOQX5lPXmUeeRV55Fbmkl+RT25FLkdLj5Jemo7VYa0t19/gT5W9Codz2pFAkBCUQOewzsQFxhHiF0KwX/BJF73Qk1uRS15lHrkVubWX4qpiYgJiiA+K1y6B2nVcYBx+ej8c0oGUEgcOHFK7FFcVk1ORQ25FLjnlOeRUaJfs8mxyynOwSdtJn59e6Hl769sE+wUzKmEU5yeez6iEUYSaQgEotBSyv2g/+wr3sa9oHweLDlJlr6r9XdTtHzPoDPjp/TDpTbUXP70fRp0RiURKedLrhBD46fxqX+OnP3Hb3+BPgCFAuzYGEGAIwGwwk1+Zz5GSI7WXo6VHySzNxCZt+On8CDGd+HxD/ELwN/hjdVixO+zYHDbs0l574hJmCiPKP4pI/0jt2hxJhH8EdoedKntV7cVis1Blr8Iu7bWfc93PvdpeTYW1ggpbBeXWcipsFVRYK7DYLdpx0lH789fc1gu9dtFp1wadofa+QWfAqDPWPl5z36Q3YdQbT/rMBAKrw3ryxa4l0QBjAIHGQAIM2nWgMRCzwYzFZqHMWka5tZyyaue1taz2b/ZUNf8PVvvJ7xNkDOLFsS829+uhSeJM73RNTU2VGzZs8Fj5q/fl8+d3f+ejm4Yxuls0a7LWMPPHmVzd42oeHvFwi8p8Y8sbvJH2BskhyTx93tP0iezT6jizy7NZdnQZPx75kU05m5BIUkJSuCD5AsYnj6d3RO96a0TFVcWsz17Pmqw1ZJZn1v5D1FwbdUb89H6YDWb8Df6Y9WZMBhNmvXb2vqNgB5tyNrGncA8O6UAv9HQN60pWeRal1aUIBH0i+zAyXqvR6IWeTbmb2JSzibTcNEqtpYCWfGzSVvuF0piaL5uOwR1JDk0mJSSFjsEdSQlNIdIcSbWjmiMlRzhYfJCDRQe16+KD5FbkUlpd2uA/J0CQMYjogGhiAmII8QshryKPrLIsciubv0iBUWckNiCWmIAYYgNjSQhKOC351CTwFRkrWJmxkuOW4+iFnu7h3cmtyKXAUlBbXrBfMF3DuhJk1Fb3qfl9Cuc0KJu0UWWrotpeXfuFW22vxuqwascI7dia10kptZME5/F1k29T/A3+2ucfkkxySDL+Bn9Kq0spqS6htLq09rbFZqn94tXr9BiEdlsiKbQUkl+ZT1FVUbM/2/riCTAE1Ca4QGMgJr0JnU6HQKATOnToan92h3Rgl3bsDnttIrNJW21yq0lwNScvNZ9TzYlUfQSi9n8HoMJWcdKJUkPMejMBxgAM4vRzfIlEr9Nj1BlPvuiNhJvCeWncSy36vIQQG6WUqfU+p5JG/ZYdXcaA6AFE+Uc1etyrP+/j+aV72fLPiUhRzh8W/oFAv0A+v+TzVtUSfj/2Ow+ufpDjlce5fdDtTO8zHb3O9b0ziixF7Cvax478Hfx09KfazvWuYV2ZmDyRCckT6BrWtc12+yurLmNL3hY25W5iW9424oLiGBk/khEdRjTYFGV32NlftJ/NuZvJKM046ayu5rZJbyLKP4rogGii/aOJ8o/CT+/X4jillFTaKk/6crNLO9H+WqJoaMOsKnsV2eXZZJZl1tYYdOjQCe2LSCe0L6dQU2htoggzhTXr83dIB9vzt7MiYwVbcrfQIbAD3cK70TWsK93CuxHtH+3R32fNWXuVvYpKWyUVtgrt2lpRez/SHElySLJbY7E6rByvPE6BpYBCSyF6nR6z3qydrDhPUkx6E0adUfusT/ncDcLQrP+d1pJSYnPYtBotDvx02snVqTHUNOmWW8u1GpDzc/Q3+BPkF0SQMYgAY4BLzbvuppJGM5NGkaWISV9NIjogmncnvkt0QMMTBG+au54jxyv48a7z+NuKv7E8fTmfXPwJvSJ7tTZ0iquKeXzt4yw5vITBMYN5avRTxAedaAZzSAe5Fbmkl6ZzpOQIB4oOsK9oHweKDpBfmV97XM+InlyQfAETkifQObRzq+NSFOXMppJGC2oaG3M28pef/kJ0QDRzJs6hQ2CH046RUjLk3z8xrkcU3Xuu55XNr3DXkLuY0XeGO0KvfY9vD37Lk78/iUBwcaeLyanIIb00ncyyzJNmkPsb/OkS2oWu4V3pGta19gy0LTrUFUU5c6ik0cI+jbTcNG776TbCTGG8d+F7xAXFnfT80YIKzn/hW/oP+p6D5Ru5MOVCnhn9jEeqwhmlGfzjt3+wI38HicGJJAUn1V4SgxPpGNyR+KB4dEKNolYUpXVU0mhFR/i2vG3c+tOthPiFMGfiHBKDE2ufe2HVd7y75wlMpiruG3ovV/e4us36CBRFUTylsaShTkub0C+6H3MmzqG0upTpS6ZztOQoNoeNlze9zPsH/46QZj6aNI8/9fyTShiKopzx1DwNF/SO7M17F77HzUtvZvoP04kPiictL43g6nNIlNfQJ6r1nd6KoijtgappuKhHRA/evfBdbNLGvqJ9PH7Ok+QdvozBSad3kCuKopypVE2jGbqHd+fLKV8CcCRXj9X+22kr2yqKopzJVNJopprJft8cPQjAoHr20FAURTlTqeapFtqcXkRCmD8xIb6x2J2iKEpbUEmjhdKOFjFQNU0pinKWUUmjBXJLLGQWVaqmKUVRzjo+kzSEEBcJIfYIIfYLIR6o5/nrhBBbnZffhBADvBEnaE1TAIM6hjd+oKIoyhnGJ5KGEEIPvAZMAnoD1wghep9y2CHgfCllf+Bx4O22jfKEtPQijHpBn/gQb4WgKIriFT6RNIBhwH4p5UEpZTXwGXBZ3QOklL9JKQudd9cCiXjJ5qOF9I4LwWxsu+WWFUVRfIGvJI0EIL3O/QznYw25Cfi+oSeFEDOFEBuEEBvy8vLcFKLGZnewNaNYNU0pinJW8pWkUd+iTfWupCiEGIuWNO5vqDAp5dtSylQpZWp0dMN7YbTE3pwyKqrtDFSd4IqinIV8ZXJfBpBU534ikHXqQUKI/sAcYJKUsuDU59tCWm0neJg33l5RFMWrfKWmsR7oJoToJITwA/4ELKx7gBCiI/AVcL2Ucq8XYgRg09FCIgL96BhR/7afiqIoZzKfqGlIKW1CiDuAJYAeeE9KuUMIMcv5/JvAP4BI4HXnEuS2htZ795S9OaUs3JLF5H5xahl0RVHOSj6RNACklIuBxac89mad2zcDN7d1XDWqbQ5mf5ZGsMnAQ5PVUuiKopydfCZp+LoXf9rLzmMlvH39EKKCTN4OR1EUxSt8pU/Dp204fJw3Vxzg6tQkJvZR+2coinL2UkmjCWVVNu6ev4WEcH8eufTUSeqKoihnF9U81YR/L9pJemEF828dSZBJfVyKopzdVE2jET/tzOGz9enMOr8LQ1MivB2OoiiK16lT5wbkl1XxwFdb6RUXwl0Tuns7nJZx2OHwKijPh6rSky82Cwy9CTr083aUiqK0Iypp1ENKyRcfvoqpMo4Xb74UP0M7rJBlbYZvZ8OxtFOeEGAKAYcVdnwN0xapxKEoistU0qhH6fFcrs99nmuDwwkxDgeCvR2S66pK4ecnYN1bEBgNl78F8YO0RGEKBmMA6HRQeBjenwwfXgY3LoJY1cmvKErT2uEptOeFRMZinLGQYJ0F3p0IWWneDsk1uxbBq8Pg9zchdQbcvg4G/Amie0BIHJiCtIQBEJ4CNy4EnVFLHHleW5lFUZR2RCWNBvh1HIqYsUQ7M597CRxc4e2QGlZ4BD69Bj6/DgIi4OafYPJ/wD+s8ddFdoEbv9Vuf3ApFBzweKiKorRvKmk0Jqob3LQEwpLg4yu1PgBfYimBnx6FV4fCwV/ggsdh5i+Q2IwluaK7ww0LwF4NH0zRmq0URVEaoJJGU0LiYfpiSBgC/5sO697xdkRgt8GG9+GVwbD6v9DncrhjA5z7f6A3Nr+82N5a4qgu02ocRelNv0ZRlLOSShqu8A+H67+G7hfB4nvg15e8F8uBn+Gt0bBoNkR2g1uWwxVvQWhjGx26IK6/9jNWFsEXM0DWuweWoihnOTV6ylVGf7h6HnwxHZY9Bj0u1pqv2tJPj2o1i/AU+OOH0GsKuHOJ9oTBcOGTsPAO2PY/6P9H95XtCTsXwtKHtJqXzgA6vVbT0hlA6EE6QNq1+SrSDg6blgw7j4GxD2mDAxRFaRYhz/AzytTUVLlhwwb3FViWC68MgaRhcN0X7v3SbvR98+C/faDnZLj8TTB4aKVdhwPmjIPSbK3JyxTkmfdprQM/w8d/hOieED9ASwwOW52LHYROu9QkFKEHexXs/k57bOQdWpOeqR0NqVaUNiCE2NjQfkWqptFcQTEw5gFY8iDsXQI9Lmqb993wrvaFN/ZBzyUM0IbkTnoW3r1Aq9WMf8Rz79VS6evgs+u0hDFtUdOjxE51/JBWW1z5LGycq32mg64Hvfp3UJSmqJpGS9it8MY52hntX9Z69kscwGrRahkJg+G6/3n2vWp8eQvsXAC3/w4RndrmPV2RswPev1gbWjxjiZbEWypjAyx9GI6ugageMPHf0H2i+2L1NCnh+EHI3ARZm7Tr7G1aLSogwnmJPHGJ7Aopo7Wh1mrnSaURjdU0fCZpCCEuAl5C2+51jpTy6VOeF87nLwYqgGlSyk1NleuRpAGwfxnMuwIm/AtGzXZ/+XVt+kjrZ7hhgdYe3xZKsuCVVOg6TuvL8QXHD8F7F2pNTjOWQHhy68uUUmuu+vEfcPwAdJ8EFz0JEZ1bX7YnFKVr8e5bApkbwVKsPW4wQ4f+ED9Q+3wqCupcjmvrj9kqtWODOkDKKOflLEsiDru2akJ1GSC00ZFny8/eDD6fNIQQemAvcAGQAawHrpFS7qxzzMXAnWhJYzjwkpRyeFNleyxpgDah7tBKuHMjBHtocyYp4Y1ztT/sWavb9g985fPw8+Nww0LofH7bvW99SrO1hGEphuk/QExP95Zvq9Zm0q94RqtJnvtXGHUX+AW4931aIm8v7FoIuxdpa4qB1jTXcQTED9aWiYnp1fhwaymhYL+2gOXhX+HwaijL1p4zhWpfnqddEiGyM4R29F7TXVUZHNuiJcjMjZCzXVvFwBx6yiXkREI46VKiJYiqUq2smsRZwxSifXYxvSG2j3Yd3UNLwkJ3oi9M6E6spuBtDgeU50FJBhRnav8TNos218pW5by2aD/DmAda9BbtIWmMBB6VUl7ovP93ACnlU3WOeQv4RUr5qfP+HmCMlPJYY2V7NGkcPwivDYe+f9A6pz3hwHL4aCpc9hoM+rNn3qMhVgu8Ngz8AuHWVd774ijL1ZY6KTyizWBPHOK59yo5ptU6ts2H0CRtNFmvS9s2WZcXwNHftC/3Az9D/h7t8YRU6HUJ9LwUorq27j2k1FYAOLxK+yIuOQalWVoNsywXqPO9oDNqI/Yiu2q1ksgu2u2o7hAU2/rPprpCe9+SzBOX44e1Jre83dooOICwjhA3QLttKT79IvRa8jAFOy8h4BekDeYwBTtvB5+4ba/Wys/ZCbk7TtTaGqMzaslZbzz5tikUAiO19d4CoiDQeTH4g7UCrJVgLdeuqyu05GWr0r7c617bq7Wf49T30Bm0WmNxuvZZ2aubjjMkDmZva9GvpD10hCcAdWeUZaDVJpo6JgE4LWkIIWYCMwE6duzo1kBPEtFZG4Gz+gVIvQmShrr/Pda+rv0h9r3S/WU3xWiGC5+Az/8MG9+HYbe0zfs67NpZ5f5lsP8n7ctDZ9D6czyZMED7R/vDO5A6HRbfC/Ovh8ShWjNO/EDtrD40yX1JxOGA4qPaz3vEmSjydmnPGcyQNByG3qyNmmvtXJy6hNAST33Jx27VanbF6VpiOX5Aq6UUHISDy7UvuBp+wVoZkd20RCJ0UHXqF3qJ9praoc9253Boh1YLqCw8PYbAGO3z7jVFm1ibMFj7Em6IlK37nUipfRnn7tR+Vru1TqzyxG2HTfvCrrm2W7VrSwlU5GsrKpQXQHVp/e8j9NpJmMGs/X8ZzKA3af2iBrOWzKRdG0ZurdRWo7bbtGv/CO2z6H2ZVgsMTYCQBK3vqrYMk3bbg7UiX6lpXAVcKKW82Xn/emCYlPLOOsd8BzwlpVztvL8MuE9KubGxsj1a0wCtyvtqqtY8dfPP7v1l5e3RzvTHPAhj7ndfuc0hJXw4BY5thf/brP2Besr+n2DzPK12ZSnSvoAShkCX8do/SluvxGu3acly0weQu0v7ogDtnzd+0MmXptrGHXYoOqr9TvN217ns1c5AQfvCSBoOyedo/Q3xg8Hg5/mfszkcDq0mULAP8vc7r/dqt0sytGOMAdpZft2mI4P5RFOPztncI/Ta/KeQeAhNdDaJOb8IjWbv/pytZbVoScRq0Zo4jc6L3tgu+lDaQ00jA0iqcz8RyGrBMW3PFAQXPAZf3QJpH8Pg691X9to3tLOG1BnuK7O5hICLnoE3R2kjjS592f3NVFLCmle18oNioecl0HW81unvySTVFL1Bq10Nu0X758/ZAcc2a/0KWWnakGRp144NjDmRQDr01WbWF+zTztQL9mtNmXWbFILjtH6JwTdobegd+mtNL74+7Fen09ZiC0uCLuNOfs5q0ZKBryU6bzCatUR4BvKVv9D1QDchRCcgE/gTcO0pxywE7hBCfIbWdFXcVH9Gm+l3Fax/V+s0HnitdibVWhXHYctn2qzsoOjWl9casb3hnDvh1xe16vuUV9y3cZPDDj88AOveht5Ttf0/fPEs02jWmsbqNo9ZKyF7uzOJOC/7llLbH6D305owI7tqS9BEdtGG9kb3aP7ckvbAF39vitv5RNKQUtqEEHcAS9CG3L4npdwhhJjlfP5NYDHayKn9aENup3sr3tMIodUGvp6pnY3G9W99mRve0zrLRvyl9WW5w4RHtTPh7++Dt87XRhedf5/WvNBS1eXw5c2wZ7GWlCY85jsjVFxh9Nf6ser2ZVWVaU1QARFax607TiAUxYf4RNIAkFIuRksMdR97s85tCdze1nG5LOVc7frIr61PGrZqbTXdLuN8Z0c9IaDvFVqT0dJHtM7/nQvg0peg0+jml1eWC59crW1He/HzbdfJ7mmmIM931iuKF7Wj0zofF5oIYcna+PfW2vGVNoZ+hA/myIAImPoaXP+N1p7/wSXw1UwtgZTmuFZG/j6YM0HrXL764zMnYSjKWcBnahpnhJRRsOd7bYRJa5pZfn9La/vuOt59sblbl7Fw2xr45SmtP2Lr59rj4Z2g40joOBziBmozkQv217kc0IZyBkbB9O+00VGKorQbKmm4U/K52giqvN0tb1aqLNTmJYx92PeH5vkFwMTHYdzD2qzdo2u1y74lsOWTk481hWgdwh1HQOSftQEDYUn1l6sois9SScOd6vZrtDRpZDinnSQNc09MbcFg0uJNGqYtNV6zZEX2Nm3+SmRXbYKirydBRVGapJKGO4UlazM1D69ueTt9xjrnpLbB7o2tLQmhbVDV1ptUKYricaoj3J2E0GobR35t+XapGeu1RdPUxkCKovgglTTcLflcbQXK/H3Nf63DoTVPJXpgDStFURQ3UEnD3VJGaddHWjD0Nn+PttibShqKovgolTTcLaKztslNS+ZrZKzXrttTJ7iiKGcVlTTcraZf43AL+jXS14F/uDbaSFEUxQeppOEJyedqM7qPH2ze6zLWa01Tamiqoig+SiUNT6jp12hOE1VlkTYpUPVnKIriw1TS8ISo7tpktiO/uv6aTOekPpU0FEXxYSppeIIQ2u5rzenXyFgPCLUWk6IoPk0lDU9JHqVtf1l0xLXjayb1mUM8G5eiKEorqKThKTXrUB12oYnK4XB2gte7Ja+iKIrP8HrSEEJECCF+FELsc16H13NMkhBiuRBilxBihxDir96ItVmie2nDZ13p1yjYB5ZiNT9DURSf5/WkATwALJNSdgOWOe+fygb8TUrZCxgB3C6E8JEt7Rqg02lDb10ZQVUzqS9RJQ1FUXybLySNy4APnLc/AKaeeoCU8piUcpPzdimwC0hoqwBbLPlcrU+jOKPx49LXgTlUTepTFMXn+ULSiJVSHgMtOQAxjR0shEgBBgG/N3LMTCHEBiHEhry8PHfG2jyu9mvUTOprzW5/iqIobaBNvqWEED8JIbbXc7msmeUEAV8Cs6WUJQ0dJ6V8W0qZKqVMjY6Obm34LRfbV6tBNLZ4oaVY2ytbNU0pitIOtMkmTFLKCQ09J4TIEULESSmPCSHigNwGjjOiJYyPpZRfeShU99LpoeM5jdc0MjcBEpLUpD5FUXyfL7SHLARudN6+EVhw6gFCCAG8C+ySUr7QhrG1Xsq5cPyAVpuoj5rUpyhKO+ILSeNp4AIhxD7gAud9hBDxQojFzmPOBa4Hxgkh0pyXi70TbjP1uQICIuGTq6GsnkpU+jqI7qk1YymKovg4r+8RLqUsAMbX83gWcLHz9mqgfS79GpoA186HuZfAx1fBtO/AFKQ9J6VW0+g9xbsxKoqiuMgXahpnvsRUuGouZG+F/00Du1V7vGA/WIrUIoWKorQbKmm0lR4XweQXYP+PsGi2VstIX6c9p0ZOKYrSTni9eeqskjodSjJh5XMQkghlOVpfRlR3b0emKIriEpU02trYh6AkC1Y8DX5BkDRcTepTFKXdUN9WbU0IuPQl6DIeqstUf4aiKO2KShreoDfCHz+AEbfDgD95OxpFURSXqeYpbzEFw0VPejsKRVGUZlE1DUVRFMVlKmkoiqIoLlNJQ1EURXGZShqKoiiKy1TSUBRFUVymkoaiKIriMpU0FEVRFJeppKEoiqK4TEgpvR2DRwkh8oAj3o4DiALyvR1EM7S3eKH9xazi9bz2FrOvxJsspYyu74kzPmn4CiHEBillqrfjcFV7ixfaX8wqXs9rbzG3h3hV85SiKIriMpU0FEVRFJeppNF23vZ2AM3U3uKF9hezitfz2lvMPh+v6tNQFEVRXKZqGoqiKIrLVNJQFEVRXKaShgcJIZKEEMuFELuEEDuEEH/1dkyuEELohRCbhRCLvB2LK4QQYUKIL4QQu52f9Uhvx9QYIcRdzr+H7UKIT4UQZm/HdCohxHtCiFwhxPY6j0UIIX4UQuxzXod7M8a6Goj3OeffxFYhxNdCiDAvhnia+mKu89w9QggphIjyRmyNUUnDs2zA36SUvYARwO1CiN5ejskVfwV2eTuIZngJ+EFK2RMYgA/HLoRIAP4PSJVS9gX0gC/u+TsXuOiUxx4AlkkpuwHLnPd9xVxOj/dHoK+Usj+wF/h7WwfVhLmcHjNCiCTgAuBoWwfkCpU0PEhKeUxKucl5uxTtyyzBu1E1TgiRCEwG5ng7FlcIIUKA84B3AaSU1VLKIq8G1TQD4C+EMAABQJaX4zmNlHIlcPyUhy8DPnDe/gCY2pYxNaa+eKWUS6WUNufdtUBimwfWiAY+Y4D/AvcBPjlKSSWNNiKESAEGAb97OZSmvIj2B+vwchyu6gzkAe87m9TmCCECvR1UQ6SUmcDzaGeRx4BiKeVS70blslgp5THQToiAGC/H0xwzgO+9HURThBBTgEwp5RZvx9IQlTTagBAiCPgSmC2lLPF2PA0RQlwC5EopN3o7lmYwAIOBN6SUg4ByfKvZ5CTOfoDLgE5APBAohPizd6M6swkhHkJrKv7Y27E0RggRADwE/MPbsTRGJQ0PE0IY0RLGx1LKr7wdTxPOBaYIIQ4DnwHjhBDzvBtSkzKADCllTQ3uC7Qk4qsmAIeklHlSSivwFXCOl2NyVY4QIg7AeZ3r5XiaJIS4EbgEuE76/qS0LmgnE1uc/4OJwCYhRAevRnUKlTQ8SAgh0Nrad0kpX/B2PE2RUv5dSpkopUxB65z9WUrp02fBUspsIF0I0cP50HhgpxdDaspRYIQQIsD59zEeH+64P8VC4Ebn7RuBBV6MpUlCiIuA+4EpUsoKb8fTFCnlNilljJQyxfk/mAEMdv6N+wyVNDzrXOB6tDP2NOflYm8HdQa6E/hYCLEVGAg86d1wGuasEX0BbAK2of0P+tzSEUKIT4E1QA8hRIYQ4ibgaeACIcQ+tNE9T3szxroaiPdVIBj40fm/96ZXgzxFAzH7PLWMiKIoiuIyVdNQFEVRXKaShqIoiuIylTQURVEUl6mkoSiKorhMJQ1FURTFZSppKIqiKC5TSUNRFEVxmUoaitLGhBAThBAfeTsORWkJlTQUpe0NADZ7OwhFaQmVNBSl7Q0ANgshTEKIuUKIJ53rUCmKzzN4OwBFOQsNQFshdgkwR0rp6ysJK0ottfaUorQh51L5+cAR4FYp5Rovh6QozaKapxSlbfUG1qNtCmT3ciyK0mwqaShK2xoA/Ia2X8n7QohYL8ejKM2ikoaitK0BwHYp5V60DYLmO5usFKVdUH0aiqIoistUTUNRFEVxmUoaiqIoistU0lAURVFcppKGoiiK4jKVNBRFURSXqaShKIqiuEwlDUVRFMVl/w88Gig26KwsUAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import freud\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# read in xyz file\n", "# read number of particles\n", "N_particles = int(np.genfromtxt(\"data/ges2.xyz\", max_rows=1, dtype=int))\n", "cleaned_data = []\n", "# remove lines that don't contain particle data\n", "with open(\"data/ges2.xyz\") as f:\n", " for line in f:\n", " if line[0] != \"\\n\" and line[0] != \" \":\n", " cleaned_data.append(line)\n", "\n", "positions = np.genfromtxt(cleaned_data)[:, 1:4].reshape(-1, N_particles, 3)\n", "particle_types = np.genfromtxt(cleaned_data, dtype=str)[:, 0].reshape(-1, N_particles)\n", "\n", "box = freud.Box.cube(19.21)\n", "\n", "# max_k_points is the number of k-points used in the calculation,\n", "# higher values give better S(k) but takes longer\n", "k_max = 15\n", "k_min = 1\n", "bins = 50\n", "sfGe_Ge = freud.diffraction.StaticStructureFactorDirect(\n", " bins=bins, k_max=k_max, k_min=k_min\n", ")\n", "sfGe_S = freud.diffraction.StaticStructureFactorDirect(\n", " bins=bins, k_max=k_max, k_min=k_min\n", ")\n", "sfS_S = freud.diffraction.StaticStructureFactorDirect(\n", " bins=bins, k_max=k_max, k_min=k_min\n", ")\n", "sfTotal = freud.diffraction.StaticStructureFactorDirect(\n", " bins=bins, k_max=k_max, k_min=k_min\n", ")\n", "\n", "for frame_positions, frame_types in zip(positions, particle_types):\n", " Ge_positions = frame_positions[frame_types == \"Ge\"]\n", " S_positions = frame_positions[frame_types == \"S\"]\n", " sfGe_Ge.compute(\n", " (box, Ge_positions),\n", " query_points=Ge_positions,\n", " N_total=N_particles,\n", " reset=False,\n", " )\n", " sfGe_S.compute(\n", " (box, S_positions),\n", " query_points=Ge_positions,\n", " N_total=N_particles,\n", " reset=False,\n", " )\n", " sfS_S.compute(\n", " (box, S_positions),\n", " query_points=S_positions,\n", " N_total=N_particles,\n", " reset=False,\n", " )\n", " sfTotal.compute(\n", " (box, frame_positions),\n", " reset=False,\n", " )\n", "\n", "plt.plot(sfS_S.bin_centers, sfS_S.S_k, label=\"S-S partial\")\n", "plt.plot(sfGe_S.bin_centers, sfGe_S.S_k, label=\"Ge-S partial\")\n", "plt.plot(sfGe_Ge.bin_centers, sfGe_Ge.S_k, label=\"Ge-Ge partial\")\n", "\n", "# Note that the Ge-S partial must be included twice\n", "S_tot = sfS_S.S_k + 2 * sfGe_S.S_k + sfGe_Ge.S_k\n", "assert np.allclose(S_tot, sfTotal.S_k, atol=1e-5, rtol=1e-5)\n", "\n", "plt.plot(sfGe_Ge.bin_centers, S_tot, label=\"Total\")\n", "plt.title(\"Static Structure Factor\")\n", "plt.xlabel(\"$k$\")\n", "plt.ylabel(\"$S(k)$\")\n", "plt.legend(loc=\"upper right\")\n", "plt.show()" ] } ], "metadata": { "interpreter": { "hash": "1c5ea5ec1f260001eb9a4e66a3ab43a0f997dbe5fe2d97f8292a6d25178f0968" }, "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.10.4" } }, "nbformat": 4, "nbformat_minor": 4 }