ACSF-Dimer.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.2"
    },
    "colab": {
      "name": "ACSF-Dimer.ipynb",
      "provenance": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ABTjh0HkNVTp"
      },
      "source": [
        "# Atom Centered Symmetry Functions\n",
        "\n",
        "ACSFs are a convenient way of transforming atomic coortinates and types into a computer-friendly string of numbers. Each atom gets its own set of ACSFs, computed using itself as the center, and all other atomic coordinates, which encode its chemical environment.<br>\n",
        "<img src=\"https://github.com/reinimaurer1/ML-CSC-tutorial/blob/master/images/acsf-schema.png?raw=1\" width=\"400px\"><br>\n",
        "The two main type of ACSFs are two- an three-body. Each set of ACSFs becomes the input of a neural network that calculates the corresponding energy contribution. The only important quantity is the total energy of the system, given by the sum of all contributions.\n",
        "\n",
        "\n",
        "For more info see: Jörg Behler, <i>J. Chem. Phys.</i> <b>134</b>, 074106 (2011)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Tr3oh7naNVTt"
      },
      "source": [
        "# A pratical example\n",
        "\n",
        "We are going to see ACSFs in action for a simple dimer system.\n",
        "\n",
        "Since there is only one atomic species and only two atoms, we will not need the three-body terms. Due to the symmetry of the system we can just compute the ACSFs for one of the two atoms, feed them to a single NN to get the total energy directly.\n",
        "\n",
        "Here are some definitions we will need."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "emuP-HlFNVTv",
        "outputId": "8ce70446-83f7-4bf3-a57b-06e592f146f0",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "source": [
        "!pip install ase\n",
        "# --- INITIAL DEFINITIONS ---\n",
        "from sklearn.neural_network import MLPRegressor\n",
        "import numpy, math, random\n",
        "import matplotlib.pyplot as plt\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "from ase import Atoms\n",
        "from ase.visualize import view\n",
        "\n",
        "\n",
        "# cutoff for all ACSF\n",
        "Rcut = 5.0\n",
        "\n",
        "# cutoff function\n",
        "def fcut(r):\n",
        "    if r >= Rcut: return 0\n",
        "    return (math.cos(math.pi * r/Rcut)+1) * 0.5\n",
        "\n",
        "# G1 function definition\n",
        "def G1f(r, eta, Rs):\n",
        "     return math.exp(-eta*(r-Rs)*(r-Rs)) * fcut(r)"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Collecting ase\n",
            "  Downloading ase-3.22.0-py3-none-any.whl (2.2 MB)\n",
            "\u001b[K     |████████████████████████████████| 2.2 MB 5.2 MB/s \n",
            "\u001b[?25hRequirement already satisfied: scipy>=1.1.0 in /usr/local/lib/python3.7/dist-packages (from ase) (1.4.1)\n",
            "Requirement already satisfied: numpy>=1.15.0 in /usr/local/lib/python3.7/dist-packages (from ase) (1.19.5)\n",
            "Requirement already satisfied: matplotlib>=3.1.0 in /usr/local/lib/python3.7/dist-packages (from ase) (3.2.2)\n",
            "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.7/dist-packages (from matplotlib>=3.1.0->ase) (0.10.0)\n",
            "Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib>=3.1.0->ase) (2.8.1)\n",
            "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib>=3.1.0->ase) (1.3.1)\n",
            "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib>=3.1.0->ase) (2.4.7)\n",
            "Requirement already satisfied: six in /usr/local/lib/python3.7/dist-packages (from cycler>=0.10->matplotlib>=3.1.0->ase) (1.15.0)\n",
            "Installing collected packages: ase\n",
            "Successfully installed ase-3.22.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JnkxilxYNVTw"
      },
      "source": [
        "### Creating the training set\n",
        "\n",
        "In order to learn the relationship between ACSFs and the energy of the system, we need a database of ACSFs for several atomic configurations, and the corresponding energy.\n",
        "\n",
        "The sample configurations consist of the dimer, stretched and compressed. In reality the energy is calculated with quantum methods (DFT, CC, ...) but here we will use a simple Lennard-Jones function."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5uwxaFuaNVTx",
        "outputId": "ec1f4eb9-1646-4c10-bd7e-29bc99d55bc7",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        }
      },
      "source": [
        "# array of meaningful distances\n",
        "dists = numpy.arange(1.95, Rcut, Rcut/30)\n",
        "# LJ energy at those distances\n",
        "energy = numpy.power(dists/2,-12)-numpy.power(dists/2,-6) - 2\n",
        "\n",
        "plt.plot(dists, energy,'.' )\n",
        "plt.xlabel('Pair distance')\n",
        "plt.ylabel('Energy')\n",
        "plt.show()"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAVHUlEQVR4nO3dfZBdd33f8fdHkg3mIUZjuyG2QbLHPMRxDEFbV8RNCNjNGMXYExKKgSQ44LqAU3AhdUOdFkIZppRAOwFqj8ak0CIearCKARksU6c8FJnuGj8LGFdBiTEzFkKYBzPYa337xz0Ly/ru2avde/eeXb1fMzu6957fnvP96Urnc8/v3PM7qSokSZrPmnEXIEnqNoNCktTKoJAktTIoJEmtDApJUqt14y5gFI499tjauHHjuMuQpBVjamrqO1V1XL9lqzIoNm7cyOTk5LjLkKQVI8ne+ZY59CRJajWWoEjyoiR3JjmYZKKl3b9s2t2R5MNJHr2cdUqSxndEcQfwQuDz8zVIcgLwWmCiqk4D1gIXLE95kqQZYzlHUVW7AZIs1HQdcFSSh4DHAPeOuDRJ0hydPUdRVd8C/hL4O+DbwP1Vdf187ZNcnGQyyeS+ffuWq0xJWvVGFhRJbmjOLcz9OX/A318PnA+cBBwPPDbJH8zXvqq2VtVEVU0cd1zfb3hJkhZhZENPVXX2EldxNvC3VbUPIMk1wK8DH1xqbfOZ2nuAXXv2s/nkY9i0Yf2oNiNJK0qXr6P4O2BzkscAPwbOAkZ2ccTU3gO87KpdPDh9kCPXrWHbRZsNC0lifF+P/d0k9wDPBj6d5LPN68cn2QFQVTcBHwNuBm5vat06qpp27dnPg9MHOVjw0PRBdu3ZP6pNSdKKMq5vPW0Htvd5/V5gy6znbwLetBw1bT75GI5ct4aHpg9yxLo1bD75mOXYrCR1XpeHnpbVpg3r2XbRZs9RSNIcBsUsmzasNyAkaY7OXkchSeoGg0KS1MqgkCS1MigkSa0MCklSK4NCktTKoJAktTIoJEmtDApJUiuDQpLUyqCQJLUyKCRJrQwKSVIrg0KS1MqgkCS1MigkSa0MCklSK4NCktTKoJAktTIoJEmtDApJUiuDQpLUyqCQJLUyKCRJrQwKSVIrg0KS1MqgkCS1MigkSa0MCklSK4NCktTKoJAktTIoJEmtDApJUiuDQpLUaixBkeRFSe5McjDJREu71yW5o2l76XLWKEnqGdcRxR3AC4HPz9cgyWnAPwPOAJ4BnJvklOUpT5I0YyxBUVW7q+rrCzT7ZeCmqnqgqqaB/00vXCRJy6jL5yjuAH4jyTFJHgNsAZ40X+MkFyeZTDK5b9++ZStSkla7daNacZIbgCf2WXR5VX1iod+vqt1J3g5cD/wIuAV4uKX9VmArwMTERC2qaEnSI4wsKKrq7CGs433A+wCSvA24Z6nrlCQdmpEFxTAk+QdVdV+SJ9M7P7F53DVJ0uFmXF+P/d0k9wDPBj6d5LPN68cn2TGr6ceT3AV8Erikqr43hnIl6bA2liOKqtoObO/z+r30TlrPPP+N5axLkvRIXf7WkySpAwwKSVIrg0KS1MqgkCS1MigkSa0MCklSK4NCktTKoJAktTIoJEmtDApJUiuDQpLUyqCQJLUyKCRJrQwKSVIrg0KS1MqgkCS1MigkSa0MCklSK4NCktTKoJAktTIoJEmtDApJUiuDQpLUyqCQJLUyKCRJrQwKSVIrg0KS1MqgkCS1MigkSa0MCklSK4NCktTKoJAktTIoJEmtDApJUiuDQpLUyqCQJLUaS1AkeUeSryW5Lcn2JE+Yp905Sb6e5O4kf7bcdUqSBgyKJO9M8itD3O5O4LSqOh34BvDGPttcC7wXeD5wKvCSJKcOsQZJ0gAGPaLYDWxNclOSVyU5eikbrarrq2q6eboLOLFPszOAu6tqT1U9CHwEOH8p25UkHbqBgqKqrqqqM4E/AjYCtyX5UJLnDqGGVwDX9Xn9BODvZz2/p3mtryQXJ5lMMrlv374hlCVJK8fU3gO898a7mdp7YOjrXjdow2Yo6OnNz3eAW4HXJ/nnVXVBn/Y3AE/ss6rLq+oTTZvLgWlg2yJq/zlVtRXYCjAxMVFLXZ8kDWJq7wF27dnP5pOPYdOG9WNZx9TeA7zsql08OH2QI9etYdtFmxddSz8DBUWS/wS8APgc8Laq+kqz6O1Jvt7vd6rq7AXWeSFwLnBWVfXbsX8LeNKs5yc2r0nSUHRhBz2Mdezas58Hpw9ysOCh6YPs2rN/+YMCuA3486r6UZ9lZxzqRpOcA1wGPKeqHpin2f8FnpLkJHoBcQHw0kPdlqTVayk7+q7soIexjs0nH8OR69bw0PRBjli3hs0nH3NIv7+QQYPiVuBpSWa/dj+wt6ruX8R23wM8CtjZrHNXVb0qyfHAVVW1paqmk/wJ8FlgLfDXVXXnIrYlqYPG/Wm+KzvoYaxj04b1bLto85KHwOYzaFD8F+BZ9I4sApwG3AkcneTVVXX9oWy0qk6Z5/V7gS2znu8AdhzKuiV1Xxc+zXdlBz2snfymDeuHHhAzBg2Ke4FXznyib65neAu94aNrgEMKCkkr21KPBrrwab5LO+hR7uSHYdCgeOrsYZ+quivJ06tqz5zhKEmr3DCOBrr0ab7LO+iuGDQo7kpyBb2L3gBe3Lz2KOChkVQmaSS6cDTQpU/zWtigQfFy4DXApc3zLwF/Si8khnHRnaRl0JWjAXAnv5IsGBTNhXY7quq5wDv7NPnh0KuSNBJdOhrQyrFgUFTVw0kOJjl6kV+FlTQkSx028mhAizHo0NMPgduT7AR+etFdVb12JFVJeoRhDBt5NKDFGDQorml+JI3JsKZp8GhAh2qgoKiqDyQ5CnhyVfWd20nSaI16mgZpPoNOCvgC4C+BI4GTkjwTeEtVnTfK4qTVZKnnFxw20rgMOvT0ZnqT//0NQFXdkuTkEdUkrTrDmgbaYSONw6B3uHuozzeeDg67GGm16nd+QVopBg2KO5O8FFib5ClJ3g38nxHWJa0qM+cX1gbPL2jFGXTo6V8AlwM/AT5Mb+rvfz+qoqTVxvMLWsnS/+ZyK9vExERNTk6OuwytIsO43aXUZUmmqmqi37JBv/X0VHpzO22c/TtV9bxhFCh12ajvRyx13aBDT1cDVwJXAQ+Prhype0Z9P2Kp6wYNiumqumKklUgd5YVuOtwNGhSfTPIaYDu9E9oAVNV3R1KV1CGeiNbh7lDuRwHwr2a9VoAX3emw4IVuOpwNOtfTSaMuRJLUTa0X3CW5bNbjF81Z9rZRFSVJ6o6Frsy+YNbjN85Zds6Qa5FGYmrvAd57491M7T0w7lKkFWmhoafM87jfc6lzvAZCWrqFjihqnsf9nkud42R80tItdETxjCTfp3f0cFTzmOb5o0damTQEXgMhLV1rUFTV2uUqRBoFr4GQlm7Q6yikFctrIKSlGfR+FJKkw5RBIUlqZVBIkloZFJKkVgaFOs2rqqXx81tP6iyvqpa6wSMKdZZXVUvdYFCos2auql4bvKpaGqOxDD0leQfwAuBB4P8Bf1xV3+vT7q+Bc4H7quq05a1S4+ZV1VI3jOuIYidwWlWdDnyDR05hPuP9OJ35YW3ThvVc8txTDAlpjMYSFFV1fVVNN093ASfO0+7zgPfllqQx6sI5ilcA1y11JUkuTjKZZHLfvn1DKEuSBCM8R5HkBuCJfRZdXlWfaNpcDkwD25a6varaCmwFmJiY8F4ZkjQkIwuKqjq7bXmSC+mdqD6rqtyxS1JHjetbT+cAlwHPqaoHxlGDJGkw4zpH8R7g8cDOJLckuRIgyfFJdsw0SvJh4MvA05Lck+SV4ylXkg5fYzmiqKpT5nn9XmDLrOcvWbaiJEl9deFbT5KkDjMoNDLO/CqtDs4eq5Fw5ldp9fCIQiPhzK/S6mFQaCSc+VVaPRx60kg486u0ehgUGplNG9YbENIq4NCTJKmVQSFJamVQSJJaGRSSpFYGhSSplUEhSWplUEiSWhkUkqRWBoUkqZVBIUlqZVBIkloZFJKkVgaF+vLudJJmOHusHsG700mazSMKPYJ3p5M0m0GhR/DudJJmc+hJj+Dd6STNZlCoL+9OJ2mGQ0+SpFYGhSSplUEhSWplUEiSWhkUkqRWBoUkqZVBIUlqZVBIkloZFJKkVgaFJKmVQSFJajWWoEjyjiRfS3Jbku1JntCnzZOS3JjkriR3JnndOGqVpMPduI4odgKnVdXpwDeAN/ZpMw28oapOBTYDlyQ5dRlrlCQxpqCoquurarp5ugs4sU+bb1fVzc3jHwC7gROWr0pJEnTjHMUrgOvaGiTZCPwacFNLm4uTTCaZ3Ldv31ALlKTD2ciCIskNSe7o83P+rDaX0xti2taynscBHwcurarvz9euqrZW1URVTRx33HHD7Mohmdp7gPfeeDdTew+MrQZJGqaR3bioqs5uW57kQuBc4KyqqnnaHEEvJLZV1TVDL3LIpvYe4GVX7eLB6YMcuW4N2y7aPLab/0ztPeAd6iQNxVjucJfkHOAy4DlV9cA8bQK8D9hdVe9azvoWa9ee/Tw4fZCDBQ9NH2TXnv1j2Ul3KbAkrXzjOkfxHuDxwM4ktyS5EiDJ8Ul2NG3OBP4QeF7T5pYkW8ZU70A2n3wMR65bw9rAEevWsPnkY8ZSR7/AkqTFGssRRVWdMs/r9wJbmsdfBLKcdS3Vpg3r2XbR5rEP+cwE1kPTB8caWJJWh8xzemBFm5iYqMnJyXGXMVaeo5B0KJJMVdVEv2VjOaLQ6G3asN6AkDQUXbiOQpLUYQaFJKmVQSFJamVQSJJaGRSSpFYGhSSplUEhSWplUEiSWhkUkqRWBoUkqZVBIUlqZVBIkloZFB3k7VQldYmzx3aMd6eT1DUeUXSMd6eT1DUGRcd05XaqkjTDoaeO6crtVCVphkHRQd6dTlKXOPQkSWplUEiSWhkUkqRWBoUkqZVBIUlqZVBIklqlqsZdw9Al2QfsHeIqjwW+M8T1jdNq6ctq6QfYl65aLX0ZtB8bquq4fgtWZVAMW5LJqpoYdx3DsFr6slr6Afalq1ZLX4bRD4eeJEmtDApJUiuDYjBbx13AEK2WvqyWfoB96arV0pcl98NzFJKkVh5RSJJaGRSSpFYGRSPJk5LcmOSuJHcmeV2fNknyV0nuTnJbkmeNo9Y2A/bjt5Lcn+SW5uffjaPWhSR5dJKvJLm16ctf9GnzqCQfbd6Tm5JsXP5KFzZgXy5Msm/W+3LROGodRJK1Sb6a5FN9lq2I92TGAn1ZSe/JN5Pc3tQ52Wf5ovdf3o/iZ6aBN1TVzUkeD0wl2VlVd81q83zgKc3PPwKuaP7skkH6AfCFqjp3DPUdip8Az6uqHyY5AvhikuuqatesNq8EDlTVKUkuAN4OvHgcxS5gkL4AfLSq/mQM9R2q1wG7gV/os2ylvCcz2voCK+c9AXhuVc13cd2i918eUTSq6ttVdXPz+Af0/uGcMKfZ+cB/q55dwBOS/NIyl9pqwH6sCM3f8w+bp0c0P3O/fXE+8IHm8ceAs5JkmUoc2IB9WRGSnAj8DnDVPE1WxHsCA/VlNVn0/sug6KM5VP414KY5i04A/n7W83vo8E64pR8Az26GQa5L8ivLWtghaIYFbgHuA3ZW1bzvSVVNA/cDnbzR+AB9Afi9ZljgY0metMwlDuo/A5cBB+dZvmLeExbuC6yM9wR6HzyuTzKV5OI+yxe9/zIo5kjyOODjwKVV9f1x17NYC/TjZnrzujwDeDfwP5e7vkFV1cNV9UzgROCMJKeNu6bFGqAvnwQ2VtXpwE5+9qm8M5KcC9xXVVPjrmWpBuxL59+TWf5xVT2L3hDTJUl+c1grNihmacaOPw5sq6pr+jT5FjD7E8WJzWudslA/qur7M8MgVbUDOCLJsctc5iGpqu8BNwLnzFn00/ckyTrgaGD/8lZ3aObrS1Xtr6qfNE+vAjYtd20DOBM4L8k3gY8Az0vywTltVsp7smBfVsh7AkBVfav58z5gO3DGnCaL3n8ZFI1mDPV9wO6qetc8za4F/qj59sBm4P6q+vayFTmAQfqR5IkzY8ZJzqD376Bz/5GTHJfkCc3jo4B/AnxtTrNrgZc3j38f+F/VwatIB+nLnPHi8+idX+qUqnpjVZ1YVRuBC+j9ff/BnGYr4j0ZpC8r4T0BSPLY5ssrJHks8NvAHXOaLXr/5beefuZM4A+B25txZIB/AzwZoKquBHYAW4C7gQeAPx5DnQsZpB+/D7w6yTTwY+CCLv5HBn4J+ECStfTC7H9U1aeSvAWYrKpr6YXif09yN/Bdev/hu2iQvrw2yXn0vrn2XeDCsVV7iFboe9LXCn1PfhHY3nz+Wwd8qKo+k+RVsPT9l1N4SJJaOfQkSWplUEiSWhkUkqRWBoUkqZVBIUlqZVDosJPk4WaGzTuSXJ3kMS1tz0vyZ4vYxt8kmWge75i5hmKetpe21SCNm0Ghw9GPq+qZVXUa8CDwqvkaVtW1VfUf5r7eXHE8kKra0lyNPZ9LAYNCnWVQ6HD3BeCUJC9I794JX01yQ5JfhJ/ej+A9zeP3J7kyyU3Af5y9kiRHJflIkt1JtgNHzVr2zSTHNlfPfrqZjPGOJC9O8lrgeODGJDc27a9IMpk5961o1vMXSW5O774DT29ef1yS/9q8dluS32te/+0kX27aX93M/yUdMq/M1mGrOSp4PvAZ4IvA5qqq9G5Ocxnwhj6/diLw61X18JzXXw08UFW/nOR0ehMvznUOcG9V/U6z/aOr6v4kr+fn7yNweVV9t7mK+3NJTq+q25pl36mqZyV5DfCnwEXAv6U3HcOvNutd38zd9efA2VX1oyT/Gng98JZD/ovSYc+g0OHoqFnTm3yB3pQTTwM+2sztcyTwt/P87tV9QgLgN4G/Aqiq25Lc1qfN7cA7k7wd+FRVfWGebfzT9KaJXkdv6o9TgZn1zUzyOAW8sHl8NrOmyaiqA+nNjHoq8KVmWocjgS/Psz2plUGhw9GPm+m+fyrJu4F3VdW1SX4LePM8v/ujxW60qr6R3u0ntwBvTfK5qvq5T/hJTqJ3pPAPmx3++4FHz2oyM5Ppw7T//w29e168ZLH1SjM8RyH1HM3Pplx+eVvDeXweeClAeveZOH1ugyTH0xue+iDwDmDmnsU/AB7fPP4FemF0f3Oe5PkDbHsncMms7awHdgFnJjmlee2xSZ66iH5JBoXUeDNwdZIpYL57Dre5Anhckt30zgP0uxnOrwJfaYa93gS8tXl9K/CZJDdW1a3AV+lNQf4h4EsDbPutwPrmBPmt9M537KM30+mHm2GwLwNPX0S/JGePlSS184hCktTKoJAktTIoJEmtDApJUiuDQpLUyqCQJLUyKCRJrf4/owll997RWLYAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true,
        "id": "SSDyIgo3NVTy"
      },
      "source": [
        "Then we calculate the ACSFs for each dimer configuration. The results are formatted as a matrix: one row for each configuration, one column for each ACSF."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MTqnrp1ANVTz",
        "outputId": "6c35b186-1f6f-4700-d345-617b17ac05ec",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 317
        }
      },
      "source": [
        "# ACSFs G1 parameter pairs: this is a list of eta/Rs values\n",
        "params = [(0.4, 0.2),(0.4, 0.5)]\n",
        "\n",
        "# initialise a matrix that will store the ACSFs of the first atom in all dimer configurations\n",
        "nConfs = dists.shape[0]\n",
        "acsf = numpy.zeros((nConfs, 1+len(params)))\n",
        "\n",
        "print(\"Number of configurations: \" + str(nConfs))\n",
        "print(\"Number of ACSfs: \" + str(acsf.shape[1]))\n",
        "\n",
        "\n",
        "for k in range(nConfs): # for each configuration\n",
        "    \n",
        "    r = dists[k] # distance between atoms\n",
        "    # compute G0 - sum of cutoffs\n",
        "    acsf[k,0] = fcut(r)\n",
        "    \n",
        "    # compute all the G1\n",
        "    for p in range(len(params)):\n",
        "        # extract parameters\n",
        "        eta,rs = params[p]\n",
        "        # compute G1\n",
        "        acsf[k,1+p] = G1f(r, eta, rs)\n",
        "\n",
        "# plot the Gs as a function of distance\n",
        "for a in range(acsf.shape[1]):\n",
        "    plt.plot(dists, acsf[:,a])\n",
        "plt.xlabel('Pair distance')\n",
        "plt.ylabel('ACSFs')\n",
        "plt.show()"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Number of configurations: 19\n",
            "Number of ACSfs: 3\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVxU5f7A8c93hn0RUVDccAm13DVcsbqZllmpLZaWpVmZli33VvdX3W51bb/ddk0zNbPNytTMLM0yc1dwwR1xxwURFwSUbZ7fHzMYEiDCDAPM9/1qXpw55+GZ72lkvvOc8yxijEEppZTnsrg7AKWUUu6liUAppTycJgKllPJwmgiUUsrDaSJQSikPp4lAKaU8nEsTgYj0FZEdIpIoIk8XcfwdEdngeCSIyElXxqOUUuqvxFXjCETECiQAfYAkYC0wxBiztZjyjwAdjTEjXBKQUkqpIrmyRdAFSDTG7DbGZAMzgAEllB8CfOXCeJRSShXBy4V1NwAOFHieBHQtqqCINAaaAr9dqNKwsDDTpEkTZ8SnlFIeIy4u7pgxJryoY65MBBdjMDDTGJNX1EERGQmMBIiMjCQ2NrYiY1NKqSpPRPYVd8yVl4YOAo0KPG/o2FeUwZRwWcgYM8kYE22MiQ4PLzKhKaWUKiNXJoK1QHMRaSoiPtg/7OcWLiQilwKhwEoXxqKUUqoYLksExphcYAywANgGfGOM2SIiY0Wkf4Gig4EZRqdBVUopt3DpPQJjzHxgfqF9zxd6/qIrY1BKKVUyHVmslFIeThOBUkp5OE0ESinl4TwmEew9lsEbP28nz6b3pJVSqiCPSQQLtx5hwu+7GPV5HGeyixy3ppRSHsljEsHIKy/hP/1bs2hbMoM/XsWx9Cx3h6SUUpWCxyQCgGE9mvDR0MvZcSSNWz5cwa6UdHeHpJRSbudRiQDg2tYRfPVANzKycrl1wgrW7j3u7pCUUsqtPC4RAHSMDGXWQz0IDfDhrsmr+TH+sLtDUkopt/HIRADQuHYgs0b3oG2DEB7+ch0f/7EbneVCKeWJPDYRAIQG+vDF/V3p1zaCV+Zv48W5W7R7qVLK43h0IgDw87YybkgnHriiKZ+u3KfdS5VSHsfjEwGAxSL864ZWvHhTK+1eqpTyOJoIChge01S7lyqlPI4mgkIKdy+N1e6lSqlqThNBEQp2L71Tu5cqpao5TQTFaFw7kO+0e6lSygNoIihBLe1eqpTyAJoILkC7lyqlqjtNBKVQuHvp0CmrOZmZ7e6wlFLKKTQRXIThMU0Zf2cnNiWd4raJKzl48oy7Q1JKqXJzaSIQkb4iskNEEkXk6WLK3C4iW0Vki4h86cp4nKFf23p8OqILyWlnueXD5Ww/kubukJRSqlxclghExAqMB64HWgFDRKRVoTLNgWeAGGNMa+BxV8XjTN0vqc23o7oDMGjiSlbtTnVzREopVXaubBF0ARKNMbuNMdnADGBAoTIPAOONMScAjDFHXRiPU10aUYNZD8VQJ9iXe6au4adNOtZAKVU1uTIRNAAOFHie5NhXUAughYgsF5FVItLXhfE4XYOa/swc1YM29Wvw0JfrmL5yr7tDUkqpi+bum8VeQHPgb8AQ4GMRqVm4kIiMFJFYEYlNSUmp4BBLZp/KuhvXXFqX57/fwpsLtuvAM6VUleLKRHAQaFTgeUPHvoKSgLnGmBxjzB4gAXtiOI8xZpIxJtoYEx0eHu6ygMvK38fKxKGdGNKlEeMX7+KfM+PJybO5OyyllCoVVyaCtUBzEWkqIj7AYGBuoTJzsLcGEJEw7JeKdrswJpfxslp49ea2PHZNc76NS+KB6bFkZue6OyyllLoglyUCY0wuMAZYAGwDvjHGbBGRsSLS31FsAZAqIluBxcBTxpgq2wVHRPh7nxa8enNb/khIYcjHq0nVdQ2UUpWcVLXr2dHR0SY2NtbdYVzQwi1HeOSr9dSv6c/0EV1oVCvA3SEppTyYiMQZY6KLOubum8XV1rWtI/ji/q4cz8jmlgkr2HzwlLtDUkqpImkicKHoJrWYOao73hZh8KRVLE885u6QlFLqLzQRuFjzusF891APGtT0Z/gna5i78ZC7Q1JKqfNoIqgA9UL8+WZUdzpGhvLoV+uZsmyPu0NSSqlzNBFUkBB/b6aP6ML1bSJ4ad5WXv9JB54ppSoHTQQVyM/byrg7O3FX10gmLtnF/30XT64OPFNKuZmXuwPwNFaL8PLANoQF+fLerzs5npHDuDs74udtdXdoSikPpS0CN8gfePbSgNb8uj2Zu6es5lRmjrvDUkp5KE0EbnR39yZ8MKQjGw6c5PaPVpKcdtbdISmlPJAmAje7sV19pt3bhaQTmdzy4Qp2paS7OySllIfRRFAJxESFMWNkd87m5DFo4ko2Hjjp7pCUUh5EE0El0bZhCDNH9yDAx8qQj1exdGflWndBKVV9aSKoRJqGBTJrdA8iawUwYtpaHYWslKoQmggqmTo1/Pj6Qfso5MdmrGfach2FrJRyLU0ElVD+KOQ+l9XlxR+28tbCHToKWSnlMpoIKik/bysf3tWJO6Ib8cFviTw7ezN5Nk0GSinn05HFlZiX1cLrt7YlLNiH8Yt3cTwji/cG6yhkpZRzaYugkhMRnrruUl64qRULtiQzbOoa0s7qKGSllPNoIqgi7o1pynuDOxC37wR3fLSKozoKWSnlJJoIqpABHRowZXhn9qVmcOvEFew9luHukJRS1YBLE4GI9BWRHSKSKCJPF3F8uIikiMgGx+N+V8ZTHVzVIpwvH+hG+tlcbpu4UtdCVkqVm8sSgYhYgfHA9UArYIiItCqi6NfGmA6Ox2RXxVOddGhUk29H9cDHKgyZtIqVu1LdHZJSqgpzZYugC5BojNltjMkGZgADXPh6HiWqThDfPdSDuiF+DPtkDT9vPuLukJRSVZQrE0ED4ECB50mOfYXdKiLxIjJTRBq5MJ5qp16IP98+2J3W9Wvw0BdxzFiz390hKaWqIHffLP4BaGKMaQf8AnxaVCERGSkisSISm5Kik7EVFBrowxf3d+WK5uE8PWsT4xcn6ihkpdRFcWUiOAgU/Ibf0LHvHGNMqjEmy/F0MnB5URUZYyYZY6KNMdHh4eEuCbYqC/DxYvKwaAZ2qM+bC3bw0rxt2HQUslKqlFw5sngt0FxEmmJPAIOBOwsWEJF6xpjDjqf9gW0ujKda87ZaePv2DoQG+jB1+R6OZ2Tx5qD2eFvd3ehTSlV2LksExphcERkDLACswFRjzBYRGQvEGmPmAo+KSH8gFzgODHdVPJ7AYhGev7EVYUG+vLlgBycyc5gwtBMBPjqTiFKqeFLVridHR0eb2NhYd4dR6c1Ys59nZ2+iXcOafDK8M6GBPu4OSSnlRiISZ4yJLuqYXjeopgZ3ieTDuy5n6+E0Bn20kkMnz7g7JKVUJaWJoBrr2yaC6SO6kHzqLLdNWEHi0XR3h6SUqoQ0EVRz3ZrVZsaD3cjOMwyauIINB066OySlVCWjicADtK4fwnejuxPs582dH6/ijwQdi6GU+pMmAg/RuHYgM0d1p3HtQO77dC1zNx5yd0hKqUpCE4EHqVPDjxkju9ExMpTHZqzn0xV73R2SUqoS0ETgYUL8vZk+ogt9LqvLC3O38NbCHTolhVIeThOBB/LztvLhXZ0Y0qURH/yWyDOzNpGbZ3N3WEopN9Ehpx7Ky2rh1ZvbEh7ky/u/JZKakc0HQzri5211d2hKqQqmLQIPJiL849qW/Kd/axZtS+buKas5lZnj7rCUUhVME4FiWI8mfDCkIxsOnOT2j1Zy5NRZd4eklKpAmggUADe2q8+0e7tw8OQZbtVRyEp5FE0E6pyYqDBmjOxGVm6ejkJWyoNoIlDnadMghJmjehDs582QSav4fcdRd4eklHIxTQTqL5qEBTJzdHeahgVy/6exzF6f5O6QlFIupIlAFalOsB9fP9iNzk1q8fevNzJ56W53h6SUchFNBKpYwX7eTBvRmX5tI3j5x228Nn+bjkJWqhrSAWWqRL5eVj4Y0onagVv46I/dpKRn8cat7XQtZKWqEU0E6oKsFmHsgNaEB/vy9i8JnMjIZvxduhayUtWFfq1TpSIiPHpNc167pS1LElK4a/JqTmRkuzsspZQTaCJQF2WIYy3kLYfSuG3iCg7qWshKVXkuTQQi0ldEdohIoog8XUK5W0XEiEi0K+NRztG3TQSfjejC0dNZ3PLhcrYdTnN3SEqpcnBZIhARKzAeuB5oBQwRkVZFlAsGHgNWuyoW5Xxdm9Xm21HdEYTbJ65kxa5j7g5JKVVGrmwRdAESjTG7jTHZwAxgQBHlXgLeAHSmsyrm0ogazHqoB/Vq+jFs6hpd/lKpKsqViaABcKDA8yTHvnNEpBPQyBjzowvjUC5Uv6Y/3z7Yg46RoTz61XodeKZUFeS2m8UiYgHeBp4oRdmRIhIrIrEpKSmuD05dlJAA+/KXN7Stx8s/bmPsD1ux2XTgmVJVhSsTwUGgUYHnDR378gUDbYDfRWQv0A2YW9QNY2PMJGNMtDEmOjw83IUhq7Ly87bywZCO3BvThKnL9/DIjPWczclzd1hKqVJw5YigtUBzEWmKPQEMBu7MP2iMOQWE5T8Xkd+BJ40xsS6MSbmQxSI8f2Mr6of488r8bRw7ncWke6IJ8fd2d2hKqRK4rEVgjMkFxgALgG3AN8aYLSIyVkT6u+p1lXuJCA9c2Yz3Bndg3f4TDJq4gkM61kCpSk2q2iRi0dHRJjZWGw1VwYrEYzz4WRyBvl5MG9GZSyNquDskpTyWiMQZY4ocq6Uji5XL9IgK45tR3TEYBk1cycpdqe4OSSlVhAsmAhEJdPTwQURaiEh/EdGLvqpULqtXg1kPxVC3hn2swbx4HWugVGVTmhbBH4CfiDQAFgJ3A9NcGZSqXhrU9GfmqO60bxTCmC91rIFSlU1pEoEYYzKBW4APjTGDgNauDUtVNzUDfPjsvq5c38a+yM3L83SsgVKVRakSgYh0B+4C8kcAW10Xkqqu/LytjLuzE8N7NGHysj08OmM9Wbk61kApdyvNOILHgGeA2Y7un82Axa4NS1VXVovwwk2tqBfix2s/bedYehYf3a1jDZRyp2K7j4qIl2MsQKWi3UerjznrD/LUzI00qR3IJ/d2pmFogLtDUqraKmv30TUFKvjA6VEpjzewYwM+vbcLR9LOMnD8CuKTTro7JKU8UkmJQApsx7g6EOWZekSFMWt0D/y8Ldz+0UoWbjni7pCU8jglJQLt0qEqRPO6wcx+KIaWdYN58PM4pi7b4+6QlPIoJSWCS0UkXkQ2FdiOF5FNIhJfUQEqzxAe7MuMkd25tlVdxs7byotzt5Cn3UuVqhAl9Rq6rMKiUArw97Hy4V2X89r8bUxetocDxzN5f0hHAn1dOUmuUqrYFoExZl/BB5AOdALCHM+VcjqrRXjuxla8NKA1i3cc5Y5JK0lO01VMlXKlYhOBiMwTkTaO7XrAZmAE8JmIPF5B8SkPdXf3JkweFs3ulAxuHr+c7UfS3B2SUtVWSfcImhpjNju27wV+McbcBHTFnhCUcqlel9blmwe7k2cMt01YyR8JukypUq5QUiLIKbB9DTAfwBhzGrC5Miil8rVpEMKch2NoGOrPvdPW8tWa/e4OSalqp6REcEBEHhGRm7HfG/gZQET8AZ0PQFWYeiH+fDuqOz2jwnhm1iZe/2m7TlinlBOVlAjuwz7L6HDgDmNM/rDPbsAnLo5LqfME+3kzZVg0d3aNZOKSXTwyYz1nc3TCOqWcoaR+eWnA88aYo4X2bwZWui4kpYrmZbXwysA2NK4VwGs/befwyTN8fE80tYN83R2aUlVaSS2C94GeRezvCbzjmnCUKpmI8OBVl/DhXZ3YciiNmz9cwa6UdHeHpVSVVlIiuNwYM6vwTmPMbODK0lQuIn1FZIeIJIrI00UcH+UYqbxBRJaJSKvSh648Wb+29fhqZDcysnK55cMVuh6yUuVQUiIoaU7g0qx1bAXGA9cDrYAhRXzQf2mMaWuM6QD8F3j7QvUqla9TZCizH4ohLMiHu6es5svV2qNIqbIo6QP9qIh0KbzTsa80Hbq7AInGmN3GmGxgBjCgYAFjTMFRQoHoRHfqIkXWDmD2wzHERIXx7OxNvPD9ZnLztHezUhejpJvFTwHfiMg0IM6xLxq4BxhcirobAAcKPE/CPhjtPCLyMPAPwAfoVYp6lTpPDT9vpg7vfG6Ool0pGYy/sxMhAdrLWanSKGmuoTXYP7gFexfSYY5Dw7AnA6cwxow3xlwC/B/wXFFlRGSkiMSKSGxKio4uVX+VP0fRf29rx+o9qQz8cDmJR/UmslKlUeK1fmNMsjHmBeAVYA/2JPAfYFsp6j4INCrwvKFjX3FmAAOLiWOSMSbaGBMdHh5eipdWnur26EZ8+UA30s7kcPOHy/l9R+Hez0qpwkqadK6FiLwgItuxdyXdj32N46uNMeNKUfdaoLmINBURH+yXk+YWeo3mBZ7eAOy86DMoJWMMu0/tdlX1qhLp3KQW34+JoWFoACOmrWXKsj0Utza3UqrkFsF27NfsbzTG9DTGfACUeiinY+H7McAC7C2Ib4wxW0RkrIj0dxQbIyJbRGQD9vsEw4qprtwmbpzI4HmDSTiR4KqXUJVIw9AAZo7qTp9WdXlp3lae/m4T2bl6E1mpokhx35REZCD2b/Ex2OcZmgFMNsY0rbjw/io6OtrExsZe9O8dzTzK4HmD8bX6MuPGGYT4hrggOlXZ2GyGdxYl8MFviXRuEsqEoZcTpiORlQcSkThjTHRRx0q6WTzHGDMYuBRYDDwO1BGRCSJyrWtCdZ06AXV45+p3SM5M5qklT5Fry3V3SKoCWCzCE9e25P0hHYlPOsWAccvZdljXNlCqoAsODDPGZBhjvnSsRdAQWI+9h0+V0z68Pc91e46Vh1fy3rr33B2OqkD929fnmwe7k2uzceuEFSzccsTdISlVaVwwERRkjDnh6MFzjasCcrVbmt/C4JaDmbZlGj/u/tHd4agK1L5RTeaO6UnzOkGM/CyO8YsT9SayUlxkIqgu/tnln1xe93JeWPECW1O3ujscVYHq1vDj6we70799fd5csIPHv96g01krj+eRicDb4s1bV71FqF8ojy1+jNQzOmGZJ/HztvLe4A48dV1Lvt9wiDs+Wkly2ll3h6WU23hkIgCo7V+bd69+lxNnT/DEkifIseVc+JdUtSEiPHx1FJPuvpydR9PpP24Z6/efcHdYSrmFxyYCgNa1W/NijxeJS47jzbVvujsc5QbXto7gu9E98LZauP2jlXy+ap/eN1Aex6MTAcCNzW5kWKthfLX9K2bvnO3ucJQbXFavBvMe6UlMVBjPzdnMk9/G630D5VE8PhEAPH7543Sv152XVr3ExpSN7g5HuUHNAB+mDuvMY9c057t1Sdzy4Qr2p2a6OyylKoQmAsDL4sWbV71J3YC6/H3x3zmaqROVeSKLRfh7nxZMHR5N0olMbhq3jMU6aZ3yAJoIHEJ8Q3iv13uk56Tz99//TnZetrtDUm7S69K6zHvkCurX9GfEtLW8uygBm03vG6jqSxNBAS1CW/BKz1eIT4nnldWv6E1DDxZZO4BZo3twc8cGvLtoJ/d9upaTmfrlQFVPmggK6dO4Dw+0fYBZO2fx9Y6v3R2OciN/HytvDWrPSwPbsCzxGDeNW8bmg6fcHZZSTqeJoAhjOo7hqoZX8caaN4g9cvEznarqQ0S4u1tjvn6wOzm5hlsnrGBmXJK7w1LKqTQRFMEiFl674jUaBjfkiSVPcDj9sLtDUm7WKTKUeY/2pFNkKE9+u5Hn5mwiK1e7mKrqQRNBMYJ9gnm/1/tk52Xz2OLHOJurUxB4urAgXz67rwsPXtWMz1ft546PVnH41Bl3h6VUuWkiKEHTkKa8fsXrbD++nRdXvqg3jxVeVgvPXH8ZE+7qxM7k09z4/jJWJB5zd1hKlYsmggu4qtFVjOk4hh93/8j0rdPdHY6qJK5vW4/vx/QkNNCHoVNWM3HJLv2ioKosTQSl8EDbB+jTuA9vx73NikMr3B2OqiSi6gQx5+EYrm9Tj9d/2s7oz9eRdlYnL1RVjyaCUhARXo55mWYhzXjy9yd1DQN1TpCvF+Pu7MhzN1zGL9uS6ffeUtbpLKaqitFEUEoB3gGMv2Y8wT7BjPxlJDuO73B3SKqSEBHuv6IZ3zzYHWNg0MSVjF+cqKORVZXh0kQgIn1FZIeIJIrI00Uc/4eIbBWReBH5VUQauzKe8qofVJ/J103Gz+rHAwsfIPFEortDUpXI5Y1Dmf/YFfRtE8GbC3Zw99TVuuCNqhJclghExAqMB64HWgFDRKRVoWLrgWhjTDtgJvBfV8XjLI2CGzHluil4Wby4f+H97Dm1x90hqUokxN+bcUM68satbVm37yTXv7eUxdt14jpVubmyRdAFSDTG7DbGZAMzgAEFCxhjFhtj8uf6XQU0dGE8TtO4RmMmXzcZg+H+BfezP22/u0NSlYiIcEfnSH54JIY6wb7cO20tY3/YqgPQVKXlykTQADhQ4HmSY19x7gN+cmE8kJXutKqahTRjyrVTyLHlMGLBCJJO67QD6nxRdYKZ83AMw3s0YeryPdzy4Qp2pzjv36BSzlIpbhaLyFAgGihyvUgRGSkisSISm5KSUrYXif0EPuwOqbvKHmghUaFRfHztx5zNO8t9C+7TqSjUX/h5W3mxf2sm3xPNoZNnuPGDZXwbe0DHHKhKxZWJ4CDQqMDzho595xGR3sC/gP7GmKyiKjLGTDLGRBtjosPDw8sWTf2OkJMBn1wPyc7r/tmyVksm9ZnE6ZzTjFgwgiMZR5xWt6o+ereqy0+PXUm7hiE8NTOex2Zs4LSOOVCVhCsTwVqguYg0FREfYDAwt2ABEekIfIQ9Cbj2jlr9DjB8PogFpvWDg3FOq7pV7VZ81PsjTmad5P6F95OSWcZWi6rWIkL8+OL+bjx5bQt+3HSYfu8vZcOBk+4OSynXJQJjTC4wBlgAbAO+McZsEZGxItLfUexNIAj4VkQ2iMjcYqpzjjqXwr0/gW8N+HQA7F3utKrbhrdlQu8JpGSmcN/C+zh2RuefUX9ltQhjejXnmwe7YbPBbRNWMHHJLh1zoNxKqtq1yujoaBMbW841AtIOwfSBcHIf3PEFNO/tnOCA2COxPPTrQzQIasDU66YS6hfqtLpV9XLqTA7PztrEj5sO0zMqjLdvb0+dGn7uDktVUyISZ4yJLupYpbhZXOFq1Id750NYC/hqMGyZ47SqoyOi+aDXBxw4fYAHFj7AqSxd0UoVLcTfm3F3duT1W9oSu++4fczBDh1zoCqeZyYCgMAwGPYDNLgcZt4L679wWtVd63Xl/avfZ8+pPYz8ZSRp2WlOq1tVLyLC4C6RzHukJ+HBvtz7yVqenb2J9Kxcd4emPIjnJgIA/5pw9yxoeiV8/xCsnuS0qns06ME7V79DwokERv8ymvRs7T+uipc/5mDklc2YsWY/173zB8t26n0mVTE8OxEA+ATCkK+h5Q3w01Ow9C2nVX1lwyt566q32Jq6ldGLRpOZk3nhX1Iey8/byrP9LuPbUT3w9bYwdMpqnp29SbuZKpfTRADg7Qe3fwptb4dfx8KiF8FJN9F7Rfbiv1f9l03HNvHQrw9pMlAXdHnjUOY/esW51kHfd5dq60C5lCaCfFZvuPkjiB4By96B+U+CzeaUqvs07sNrV7zG+qPrefS3R3X9Y3VBRbUOnpmlrQPlGpoICrJY4Ia3ocejsHay/b5BnnNu2l3f9HpejnmZNUfWMObXMdqbSJVKwdbB12u1daBcQxNBYSLQZyxc/Rxs/ApmDofcIme+uGg3XXITr/R8hbijcdz5453sOum8eY9U9aWtA+VqmgiKIgJXPQV9X4dtP8BXQyDbOdf2b7rkJqZeN5WMnAzumn8Xvx/43Sn1qupPWwfKVTQRlKTbaOg/DnYvhs9vhbPOuZzTsU5HZtw4gyY1mvDob4/y0caPdDZKVSraOlCuoIngQjrdDbdOgaQ18Gl/yEh1SrURgRFM6zuNG5rdwLgN43hiyRPao0iVWlGtg6U7dbJDVTaaCEqjzS0w+EtI2Q6f9IWj25xSrZ+XH6/2fJUno5/k1/2/MvSnobrAjSq1wq2Du6es0daBKhNNBKXV4joYOgvOnIBJV8O66U4ZayAiDGs9jAnXTOBIxhGG/DiE1YdXOyFg5SkKtw6ue+cPftp0WC83qlLTRHAxmsTAqOXQqAvMfQRmPQBZp51SdY8GPZhxwwxq+9XmwV8e5IttX+gfsiq1/NbBzNE9qOHvzegv1nHP1DXs0qUxVSl45jTU5WXLg2Vvw+JXIbQJ3PaJfeEbJ8jIyeCZpc+w+MBiBkYN5N/d/o2P1ccpdSvPkJtn47NV+3h7YQJnc/O4r2czHukVRaCvl7tDU25U0jTUmgjKY98KmHkfZB6Da1+GLiPtXU/LyWZsTNg4gYkbJ9IurB3vXP0OdQLqOCFg5UlSTmfxxs/bmRmXRL0QP/51w2Xc0LYe4oR/o6rq0UTgShmp9hHICT/DpTfCgHHg75zFaH7Z9wv/WvYvgryDePfqd2kX3s4p9SrPErfvOP+es4Wth9OIiarNf/q3JqpOsLvDUhVMF6ZxpcDaMGQGXPcqJCyAiVfCgTVOqbpP4z583u9zfKw+DP95OHMSnbeAjvIclzeuxQ+P9OSlAa3ZlHSKvu8u5dX523TNA3WOtgicKSnOvsjNqSS45nn7nEWW8ufak2dP8uSSJ1l9ZDVDLxvKE9FP4GXR673q4qWmZ/Hfn3fwdewB6tbw5dl+l9G/fX29XOQB9NJQRTp7CuY+ClvnwCXX2Gc0DQovd7W5tlzein2Lz7d9Ttd6Xfnflf+jpl9NJwSsPNH6/Sd4/vstbDp4iq5NazF2QBtaRujlourMbZeGRKSviOwQkUQRebqI41eKyDoRyRWR21wZS4XxC4FB0+DGd2DvMpjYE/YsLXe1XhYv/q/L//FSzEusS17HrXNvZeHehdrFVJVJx8hQ5jwcwys3t2FH8mn6vb+UsT9sJU0Ho3kklyUCEbEC44HrgVbAEBFpVajYfmA48KWr4nALEV8ndeAAABbkSURBVPu6Bg/8Br7BML0/LH7N3u20nAZGDeSzfp9Ry78WTyx5god/fVhHI6sysVqEu7o2ZvETf+OOzo34ZMUeev1vCbPWJekXDA/jyhZBFyDRGLPbGJMNzAAGFCxgjNlrjIkHnLMCTGUT0QZG/g7tBsOS12H6AEg7XO5qW9duzVc3fMVT0U8RmxzLzd/fzORNk8nJ029z6uKFBvrw6s1t+f7hGBqE+vOPbzZy28SVrNlz3N2hqQriykTQADhQ4HmSY59n8Q2CmyfAwIlwcB1MjIEdP5W7Wi+LF/e0voe5A+fSs0FP3lv3HrfPu511yeucELTyRO0a1mT26B68cWtbDhzP5PaPVjJs6ho2JekiStVdleg+KiIjRSRWRGJTUqroDIsdhthbB8H14KvB9mmtk7eWu9qIwAjeufodxvUaR0ZOBsN+HsYLK17g5NmT5a5beR6LRbijcyRLnrqaZ/tdysakk9w0bhmjP49jZ7JzplNRlY/Leg2JSHfgRWPMdY7nzwAYY14rouw0YJ4xZuaF6q30vYYuJDcL1nwMf/zXPk9Rx7vh6n9BcN1yV52Zk8nE+IlM3zKdGj41eCL6Cfpf0l+7BqoyO302h8lL9zBl2R4ys3MZ2LEBj1/TgsjaAe4OTV0kt3QfFREvIAG4BjgIrAXuNMZsKaLsNDwlEeTLPA5L/gtrPwYvP4h5HLo/DD7l/wNLOJHA2JVj2Ziykei60fy7279pVrOZE4JWnup4RjYfLdnFtBV7ybMZBndpxCO9mlO3hp+7Q1Ol5LZxBCLSD3gXsAJTjTGviMhYINYYM1dEOgOzgVDgLHDEGNO6pDqrTSLIl7oLfnkets+DGg2g17+h3R3lHohmMzZm7ZzF23Fvcyb3DPe2vpeR7Ubi56V/uKrsktPO8sFvO5mx5gBWizCsRxNGXXUJtQJ1YsTKTgeUVQV7l8PCf8Gh9VCvPVz7CjS9otzVpp5J5a3Yt/hh9w80Cm7Ec12fo0eDHk4IWHmy/amZvPtrAnPWHyTAx4sRPZty/xVNqeHn7e7QVDE0EVQVNhtsngmL/gNpSdCyH/QZC2HNy1316sOreXnVy+xN20vfJn35Z+d/Eh5Q/hHPyrPtTD7N278k8NPmI9QM8GbUVZcwrHsT/H2s7g5NFaKJoKrJOQOrPoSl70DuGfvgtKuetk9wVw7ZedlM2TyFyfGT8bH6MLLdSAa1GESQT5CTAleealPSKf63cAdLElIID/blkV5RDO4ciY9XleiY6BE0EVRV6Snw+6sQNw18guHKJ6DrKPDyLVe1+9L28drq11h+aDlB3kEMajGIOy+7k4jACOfErTzWmj3HeXPBdtbuPUG9ED+G9WjCkM6RhAToJSN300RQ1R3dDr/8G3YuhJqR0PtFaH1LuRfB2XJsC9O2TGPhvoVYsNCvWT/uaXUPLWu1dErYyjMZY1iSkMLEJbtYtfs4/t5Wbru8IffGNKFZuLY+3UUTQXWxazEsfA6SN0PdttD5Pmg7yD56uRySTifx2dbPmJ04mzO5Z4ipH8PwNsPpGtFVxyCoctly6BRTl+3lh42HyLHZ6NWyDvf1bEr3S2rrv60KpomgOrHlwcYZsHI8HN1iv2TU7nb7fYSINuWq+lTWKb7e8TVfbvuS1LOpXFbrMoa1Hsa1Ta7F26JNe1V2R0+f5fNV+/li1T5SM7K5NCKYET2bMqBDfXy99MZyRdBEUB0ZA0lrIXYqbJ4FeVnQsIs9IbQeCN7+Za46Ky+Lebvm8enWT9lzag/1Ausx9LKh3NriVgK9A514EsrTnM3J4/sNB5m6bC87kk8TFuTD0G6NGdqtMWFB5bv3pUqmiaC6yzwOG7+yJ4XURPCrCR3uguh7y9X11GZs/JH0B9O2TCMuOY5g72AGtRzEXZfdRZ2AOk48AeVpjDEsT0xlyrLdLN6Rgo+XhYEd6jOiZ1Mujajh7vCqJU0EnsIY2LvUnhC2/QC2XGh6pb2V0PIG8Cr76M/4lHimbZnGr/t/xSIWbmx2I8NaDSMqNMqJJ6A8UeLRdKat2MPMuCTO5tiIiarNfT2b8rcWdbBY9D6Cs2gi8ESnk2HD5xA7DU7th8Bw+wR3lw+D0CZlrvZA2gGmb53OnMQ5nM07S6c6nejduDe9I3tTL6ie08JXnudkZjZfrtnP9BX7OJJ2lmZhgQzt1pib2tcnPFgvG5WXJgJPZsuDXb/ZWwkJP9tbDVG97a2E5teC1atM1Z48e5JvE77l570/k3AiAbAvmJOfFJqENHHiSShPkpNnY/6mw0xdvpeNB05iEYiJCmNghwZc27ouwTqNRZloIlB2p5Jg3XSI+xTSj4B/KFzSC6L6QNQ1EFS26/770vbx6/5fWbRvEZuObQIgqmbUuaTQIrSFdhVUZbIz+TTfbzjEnA0HSTpxBl8vC71b1WVghwZc1SJcRy5fBE0E6nx5OZCwwD7jaeIiyHAs9lOvvT0pNO8DDaLL1Fo4knHkXFJYd3QdNmMjMjiSaxpfQ+/I3rQJa4NF9I9XXRxjDOv2n2DO+kP8uOkwxzOyCfH3pl/begzsUJ/OTWrp/YQL0ESgimezwZF4SPwFdi6CpDVgbOAXUqC10LtMC+eknkll8YHFLNq3iNWHV5NrcqkbUJdrIq+hd+PedKrTCatF+5Cri5OTZ2PZzmPM2XCQhVuSOZOTR/0QP27qUJ+BHRpwaUSwtkCLoIlAld6ZE/YRzImL7I/0ZPv+iLZ/thYadrno1sKprFP8kfQHi/YtYvmh5WTlZVHLrxZXN7qabvW70T6sPRGBEfoHrC5KRlYui7YlM2f9Qf7YeYw8m6FF3SAGdGhA//b1aVRLV1LLp4lAlY0xcGTTn62FA6vB5IFvCFzyN3tiaNwDQpte1EI6mTmZLDu4jEX7FrEkaQmZuZkAhPuH0y68nf0R1o7WYa3x9yr7wDjlWVLTs5i/6TBzNhwibt8JAKIbh9K3TQQ9m4fRsq5ntxQ0ESjnOHMSdv/+Z2vh9GH7fp9ge4uhXjuIaGf/GX4pWC/cuyPHlkPCiQTiU+LPPfaf3g+AVay0CG1Bu/B2tA9vT7vwdkQGR3r0H7MqnQPHM5m78RDfbzhIQnI6AGFBPvS4JIyeUWH0iKpNw1DPai1oIlDOZwwc3QoH4+BwPBzeaJ8ML8f+7R6rD9RpVSA5tIe6rcHnwlNUHD97nE0pm9iYspH4Y/FsPraZjJwMAGr61jzXYmgX3o62YW11PQVVooMnz7A88Zjjkcqx9CwAmtQOICYqjJioMLo3q01oNV9uUxOBqhi2PPsazEccieHwRvv2GXszHbFA7ebntxwi2kFArRKrzbPlsfvUbuJT4u3JISWeXad22atEaBbSjEtqXkLjGo2JrBFJ4xqNaRTciNp+OsOlOp8xhoTk9HOJYfWe46Rn5SICbeqH0COqNj2jwujcpBZ+3tWrI4MmAuU+xtjHL5xLDvH27bSDf5YJqG1fZ+Hco/Gf2yGNipxmOy07jc3HNhOfEs+mY5vYe2ovB9MPkmfyzpUJ8g6iUXCjc4mhcY3G55JFqG+oJglFTp6N+KSTLNuZyvJdx1i//wQ5eQYfLwuXR4bSs7m9xdC6fg28rVW727MmAlX5ZBxzJId4OLEXTu7/85GXdX7ZUiaKHFsOh9MPsy9tH/tP7z/3c3/afg6lHzovSQR7B9OoRiMaB9sTQ2SNSML9w6ntX5vafrWp6VtTu7Z6oIysXNbsPc6KxGMsS0xl2+E0ALytQrOwIFpEBNOybhAt6gbTMiKYRqEBVWb8gtsSgYj0Bd4DrMBkY8zrhY77AtOBy4FU4A5jzN6S6tREUM3ZbPYBbif3w8l95yeI4hKFfy0IqmtPGAG1HD/Pf+T41eAQuezLPc3+M0fZl7afA6cPsC9tH4cyDmEztvOqtIiFUN/Qc4nhvJ+F9oX6heJlKdtUHapyS03PYuXuVLYcSiPhyGl2JJ8m6cSZc8f9vC00rxPsSAx/JoiIGn6VrsXplkQgIlYgAegDJAFrgSHGmK0FyjwEtDPGjBKRwcDNxpg7SqpXE4GHKy5RZB6zT8edmep4HLd3dS2K1adAkqhFjn8oh3z9OWb14pgFUrGRSi6ptmxSbWc5nptJak46qdlpnLVl/6U6QajpW5NafrUI9gkm0CeQQK9AgnyCCPQOJMjb/rPgdv6x/H0B3gG6+E8VkZ6Vy87k0+xMTmdH8mkSkk+z48hpjp7+8wtKsK8XLSIcCcLRgogI8aNWoA81/Lzd0opwVyLoDrxojLnO8fwZAGPMawXKLHCUWSkiXsARINyUEJQmAlUqNhtknfprcji3XcTz7Iy/tjYKMECmCKlWK6lWC6nefhzz8SPVx49ULy9OWC2ki5AhkI4hA0M6eWRiK7bOgvzEiwCrL74Wb3wt3vjk/7R642vxwcfqg6/Vx7Hti6/VBx+rH75evvhaffHx8sPX6oePly9WixdeFm+sFm+8LF54We3b555bvO1lrN5YxerY53Vu2ypWRASLWP7cxoJFLIhIkfs83cnMbBLyk4Oj9ZCQfJqTmTnnlbMI1AzwITTAm9AAH0IDHduBPoQG+FCriH0h/t5Yy5k8SkoErmzPNgAOFHieBHQtrowxJldETgG1gWMujEt5AovFPqmefyjUvqT0v5eXY08I5x7p57YlO53A7AwCszOILHTs3HZeNuRk2RNKbjbkZZGXm8sZWxbpeTlkmBzSTS4ZFrEnDYuFDIuFdIuQIRYyLEKWCNly/s90x8/8fX8eh6yLGMznShZjsMC5h+Q/zJ/b9vYTf3lwgX2U8JPijhvHKxbz+Vnsx2oRX0PL9BHsDZENIRIwGBz/2V/CsZEJZGYZkrKAExeu8vrQvjx62//KEk2JqsSFTREZCYwEiIyMdHM0qlqzeoN/TfvDWVUCQY7HOXm5jmSRZU8eBX+aPPuiQjab42f+I6/AsT/3mbxccvKyyHI8svOyyDM2ck0eubZc8kweebY8ckz+tuOYySPP2MgzeeTa/nye/9NgyDMGg8FmbNgwGGPII3+fse/DkFfguA0cxx0/HR9/9pL2D8H8bfL3Ye/a+Zd95/ZQ4Oefn9Tnf7CaQuUK7imdokpfbB2lrbfIcgb7/1dj7G9//rYBm81Qt0ZEuWMpiisTwUGgUYHnDR37iiqT5Lg0FIL9pvF5jDGTgElgvzTkkmiVqkhWL/ujFAPsLkQAH8cjuNy1KU/kyjblWqC5iDQVER9gMDC3UJm5wDDH9m3AbyXdH1BKKeV8LmsROK75jwEWYG8dTzXGbBGRsUCsMWYuMAX4TEQSgePYk4VSSqkK5NJ7BMaY+cD8QvueL7B9FhjkyhiUUkqVrHJ0N1BKKeU2mgiUUsrDaSJQSikPp4lAKaU8nCYCpZTycFVuGmoRSQH2ObHKMKrPlBZ6LpVTdTmX6nIe4Jnn0tgYE17UgSqXCJxNRGKLm4ipqtFzqZyqy7lUl/MAPZfC9NKQUkp5OE0ESinl4TQROCazqyb0XCqn6nIu1eU8QM/lPB5/j0AppTydtgiUUsrDeUQiEJFGIrJYRLaKyBYReayIMiIi74tIoojEi0gnd8R6IaU8l7+JyCkR2eB4PF9UXe4mIn4iskZENjrO5T9FlPEVka8d78tqEWlS8ZGWrJTnMVxEUgq8J/e7I9bSEhGriKwXkXlFHKv070lBFziXKvO+iMheEdnkiPMv6/WW5zOsSqxQ5gS5wBPGmHUiEgzEicgvxpitBcpcDzR3PLoCE/jr0pqVQWnOBWCpMeZGN8R3MbKAXsaYdBHxBpaJyE/GmFUFytwHnDDGRInIYOAN4A53BFuC0pwHwNfGmDFuiK8sHgO2ATWKOFYV3pOCSjoXqFrvy9XGmOLGDJT5M8wjWgTGmMPGmHWO7dPY/1E0KFRsADDd2K0CaopIvQoO9YJKeS5VguP/dbrjqbfjUfim1QDgU8f2TOAaqWQrpZfyPKoMEWkI3ABMLqZIpX9P8pXiXKqTMn+GeUQiKMjRjO0IrC50qAFwoMDzJCr5B2wJ5wLQ3XGp4icRaV2hgV0ER7N9A3AU+MUYU+z7YozJBU4BtSs2ygsrxXkA3Oposs8UkUZFHK8s3gX+CdiKOV4l3hOHC50LVJ33xQALRSTOsY57YWX+DPOoRCAiQcB3wOPGmDR3x1MeFziXddiHk7cHPgDmVHR8pWWMyTPGdMC+pnUXEWnj7pjKohTn8QPQxBjTDviFP79RVyoiciNw1BgT5+5YyquU51Il3heHnsaYTtgvAT0sIlc6q2KPSQSOa7ffAV8YY2YVUeQgUPDbQEPHvkrnQudijEnLv1ThWCXOW0TCKjjMi2KMOQksBvoWOnTufRERLyAESK3Y6EqvuPMwxqQaY7IcTycDl1d0bKUUA/QXkb3ADKCXiHxeqExVeU8ueC5V6H3BGHPQ8fMoMBvoUqhImT/DPCIROK5fTgG2GWPeLqbYXOAex533bsApY8zhCguylEpzLiISkX/NVkS6YH+fK90fqoiEi0hNx7Y/0AfYXqjYXGCYY/s24DdTyQa/lOY8Cl2r7Y/93k6lY4x5xhjT0BjTBPsa4r8ZY4YWKlbp3xMo3blUlfdFRAIdnUMQkUDgWmBzoWJl/gzzlF5DMcDdwCbHdVyAZ4FIAGPMROxrK/cDEoFM4F43xFkapTmX24DRIpILnAEGV8Y/VKAe8KmIWLEnq2+MMfNEZCwQa4yZiz3pfSYiicBx7H/QlU1pzuNREemPvdfXcWC426Itgyr4nhSrir4vdYHZju93XsCXxpifRWQUlP8zTEcWK6WUh/OIS0NKKaWKp4lAKaU8nCYCpZTycJoIlFLKw2kiUEopD6eJQFUrIpLnmJ1xs4h8KyIBJZTtLyJPl+E1fheRaMf2/PwxBMWUfbykGJSqDDQRqOrmjDGmgzGmDZANjCquoDFmrjHm9cL7HaNlS8UY088xmrg4jwOaCFSlpolAVWdLgSgRuUns8+avF5FFIlIXzs1FP86xPU1EJorIauC/BSsREX8RmSEi20RkNuBf4NheEQlzjPz80THR32YRuUNEHgXqA4tFZLGj/AQRiZVC6xY46vmPiKwT+5zzlzr2B4nIJ4598SJyq2P/tSKy0lH+W8fcU0qViaeMLFYexvGt/nrgZ2AZ0M0YY8S+8Mg/gSeK+LWGQA9jTF6h/aOBTGPMZSLSDvukfoX1BQ4ZY25wvH6IMeaUiPyD8+eQ/5cx5rhjFPKvItLOGBPvOHbMGNNJRB4CngTuB/6NfaqAto56Qx3zRj0H9DbGZIjI/wH/AMZe9P8opdBEoKof/wJTbyzFPh1CS+Brx7wyPsCeYn732yKSAMCVwPsAxph4EYkvoswm4C0ReQOYZ4xZWsxr3C72KYS9sE9N0QrIry9/AsE44BbHdm8KTOFgjDkh9lk1WwHLHVMO+AAri3k9pS5IE4Gqbs44poM+R0Q+AN42xswVkb8BLxbzuxllfVFjTILYlwbsB7wsIr8aY877hi4iTbF/0+/s+ECfBvgVKJI/C2YeJf9tCvY1D4aUNV6lCtJ7BMoThPDndLzDSipYjD+AOwHEvs5Au8IFRKQ+9stHnwNvAvnrxZ4Ggh3bNbAnm1OO+xTXl+K1fwEeLvA6ocAqIEZEohz7AkWkRRnOSylAE4HyDC8C34pIHFDceq8lmQAEicg27Nfhi1ropC2wxnFZ6gXgZcf+ScDPIrLYGLMRWI99iuovgeWleO2XgVDHDeiN2O83pGCfJfMrx2WqlcClZTgvpQCdfVQppTyetgiUUsrDaSJQSikPp4lAKaU8nCYCpZTycJoIlFLKw2kiUEopD6eJQCmlPJwmAqWU8nD/D7bifndcp3bXAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qZmz1cEbNVT0"
      },
      "source": [
        "#### OPTIONAL TRICK\n",
        "\n",
        "We can center the ACSFs around their mean and rescale them so that their standard deviation is 1. This is a common trick in ML with neural networks, to make the learning easier."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zE1NQn8TNVT1",
        "outputId": "95c9701c-26fa-4ab3-decf-76d59c101db4",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        }
      },
      "source": [
        "acsf_mean = numpy.mean(acsf, axis=0)\n",
        "for a in range(acsf.shape[1]):\n",
        "    acsf[:,a] -= acsf_mean[a]\n",
        "acsf_std = numpy.std(acsf, axis=0)\n",
        "for a in range(acsf.shape[1]):\n",
        "    acsf[:,a] /= acsf_std[a]\n",
        "\n",
        "# plot the Gs as a function of distance\n",
        "for a in range(acsf.shape[1]):\n",
        "    plt.plot(dists, acsf[:,a])\n",
        "plt.xlabel('Pair distance')\n",
        "plt.ylabel('ACSFs - scaled and shifted')\n",
        "plt.show()"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xUVdrA8d+T3kMIndC7NOlNAyIWimAXdO1dfAVZdXUXXXXVXV07FuxiBVHsUkTpAkqv0ntLQgkJqTPzvH/MwEZMwpBkMknm+b7v/cyde+/c81xnmSfnnnPPEVXFGGNM4AnydwDGGGP8wxKAMcYEKEsAxhgToCwBGGNMgLIEYIwxASrE3wGcjho1amjjxo39HYYxxlQqS5cuTVPVmidvr1QJoHHjxixZssTfYRhjTKUiIjsK2263gIwxJkBZAjDGmABlCcAYYwKUJQBjjAlQlgCMMSZA+S0BiEiEiPwqIitFZK2IPOavWIwxJhD5sxtoLtBfVTNFJBSYLyJTVXWRH2MyxpiA4bcagLplet6GehbfjE29+SeY97xPTm2MMZWVX9sARCRYRFYAKcCPqrq4kGNuE5ElIrIkNTW1ZAVtnQ2znoRjaaWK1xhjqhK/JgBVdarqmUAS0F1E2hVyzJuq2lVVu9as+acnmb3TcTi4HLBmSukCNsaYKqRC9AJS1SPALOBCnxRQuy3UbgerJvrk9MYYUxn5sxdQTRGp5lmPBM4DfvdZgR2ugj1LIW2Tz4owxpjKxJ81gLrALBFZBfyGuw3gO18Vlt7qApAgWPWZr4owxphKxW/dQFV1FdCpPMp6avFTzN09lx8aJxO0ahKc83cQKY+ijTGmwqoQbQC+1rlWZ/Zk7mF+025wZAfstEcNjDEmIBLAuQ3PJTEikUk5uyA0yhqDjTGGAEkAocGhXN7ycubtXcjulgNg7ZeQn+PvsIwxxq8CIgEAXN7ycoIkiM8SqkNOOmya7u+QjDHGrwImAdSJrsM5Dc7hy7Rl5MbUtt5AxpiAFzAJAOCq1ldxJPcI05v3hI3TIeuQv0Myxhi/CagE0KNODxrHNWaiHgVXPqy1oSGMMYEroBKAiDC89XBWH93C2jqtYOUkf4dkjDF+E1AJAGBos6FEhkQysVYS7P4VDm7xd0jGGOMXAZcAYsNiGdx0MFOP7SA9KNgag40xASvgEgDA8FbDyXXl8VWjdrBqEqhv5qExxpiKLCATQKvqrehUqxOTwhTX4W2w+zd/h2SMMeUuIBMAuGsBu/KO8EtMHKy0oSGMMYEnYBPAgEYDqB5RnUl1Gru7gzry/B2SMcaUq4BNAGHBYVzW4jLmONPZk38UNs3wd0jGGFOuAjYBAFzR8gpEhMnVa9kIocaYgBPQCaBuTF36JfVjSkw0eRunQ/Zhf4dkjDHlJqATALjHBzqseUyPCIG1X/k7HGOMKTcBnwB61u1J47hGTEqs5X4mwBhjAkTAJ4AgCeLKVlexMtjF+v1L4PB2f4dkjDHlIuATAMCw5sOICA5nUlysDQ1hjAkYlgCAuLA4BjcdwvcxsaSv/NSGhjDGBARLAB7DWw8nR5RvHKmwZ6m/wzHGGJ+zBODRunprOia2Y1JcLC4bGsIYEwAsARQw/Ixr2BEawqJNX4Ez39/hGGOMT1kCKOD8RudTPSSaiWEKm2f6OxxjjPEpSwAFhAWHcWmrK5kTFcm+5R/4OxxjjPEpSwAnuaL1cBBhcsoiyD7i73CMMcZnLAGcpF5MPZJrdOKLmEjy1nzh73CMMcZnQoraISLjgCI7xKvqPT6JqAIY3vE2Zv+0nB/XfMjgbjf7OxxjjPGJ4moAS4ClQATQGdjkWc4Ewnwfmv/0qt+bhiGxTMrbA0d2+jscY4zxiSITgKpOUNUJQAegn6qOU9VxwLm4k0CV5R4f6EqWR0Sw4bfx/g7HGGN8wps2gAQgrsD7GM+2Ku3i9jcSocLE7T/Y0BDGmCrJmwTwH2C5iLwvIhOAZcBTvg3L/+LD4xmY0JbvQ/I5umO+v8Mxxpgyd8oEoKrvAT2AL4EpQC/PraEqb3i3UWQHBfHt0lf8HYoxxpS5UyYAERFgANBRVb8GwkSku88jqwDOqNeTDhLFxCNrUEeev8Mxxpgy5c0toNeAXsAIz/sM4FWfRVTBDG8yhO0hQSxe/qa/QzHGmDLlTQLooaojgRwAVT1MFe8GWtD53UeT4FImbrDpIo0xVYs3CSBfRILxPBQmIjUBV2kLFpEGIjJLRNaJyFoRGVXac/pCeHgsl0Q3YZbjMPsPbfZ3OMYYU2a8SQAv424AriUiTwLzgX+XQdkO4K+qegbQExgpImeUwXnL3BVn3okCkxc/4+9QjDGmzHjTC+hj4AHcP/r7gItVtdQT56rqPlVd5lnPANYD9Ut7Xl9IajGQvo4gJh9YTFZ+lr/DMcaYMuFNL6APVfV3VX1VVV9R1fUi8mFZBiEijYFOwOJC9t0mIktEZElqampZFus9EW5pchGHxcUnvz7rnxiMMaaMeXMLqG3BN572gC5lFYCIxABfAKNV9ejJ+1X1TVXtqqpda9asWVbFnraOfR7grJx83t88hWP5x/wWhzHGlJUiE4CIPCQiGUAHETnqWTKAFODrsihcREJx//h/rKpTyuKcPhMRz11J55GOk0+W2YNhxpjKr7gawGZVjQWmqGqcZ4lV1URVfai0BXseMHsHWK+qz5f2fOWhffI/SM7O4f0NE8nMy/R3OMYYUyrFJYDjP/LNfVR2H+BaoL+IrPAsg3xUVtmIrcNdtXpzVB18vNIeDDPGVG5FTggDHBSRGUATEfnm5J2qOrQ0BavqfEBKcw5vOZwuRITgoNIX1zZ5LP0+G8iE9R9xdcdbiQ2LLYMIjTGm/BVXAxgMPAKkAc8VslQab8zdyuXjf2HjgYzSn6xGC+6s1p4Mzeej1e+W/nzGGOMnxU0Ik6eqi4Deqjrn5KUcYyy1BtWj2J52jMEvz+P5GRvIyXeW6nxnJP+Dc45l8eG6Dzia96eOS8YYUykU1wvoRc/quyLyzclLOcVXJoZ2rMfMMX0Z0qEeL/+8mUEvz+PXbYdKfsKkrtwV2YQMVx4frQmIkbGNMVWQaBGzXYlIF1VdKiJ9C9vvj1pA165ddcmSJaU6x5yNqfzjy9XsPpzN1T0a8uDA1sRFhJ7+iTbOYPTMO1kcm8C0K38iPjy+VHEZY4yviMhSVe168vbibgEt9bz+6fZPZbsFVFDfljWZcW8yt5zVhIm/7mTAc3OYtmbf6Z+oxXncKYlkunL5YO0HZR+oMcb4mDdDQfQRkR9FZKOIbBWRbSKytTyC85WosBDGDjmDr0b2oUZMOHd8tIzbP1zCgaM53p9EhFa9x3DesSw+XjeB9Nx03wVsjDE+4M1QEO8AzwNnAd2Arp7XSq9DUjW+vrsPf7uwNbM3pDLguTl8tGgHLpeXk8C3u4w78yPIcuQyYa21BRhjKhdvEkC6qk5V1RRVPXh88Xlk5SQ0OIg7+zVj+uhk2ifFM/arNVz15kI2p3jRZTQ4lBbdR3L+sWN8vO4DDucc9n3AxhhTRorrBdRZRDoDs0TkvyLS6/g2z/YqpXGNaD6+pQf/vbwDGw9kMuil+bw0cxN5jlPMfdP5Ou7MVrKdVgswxlQuxfUCmlXM51RV+/smpKKVRS8gb6Rm5PL4d+v4duVeWtSK4T+XtadLo+pFf+DnJ3hg/bvMjk9k2uXTqR5RzLHGGFPOiuoFVGQCqIjKKwEc9/PvBxj75Rr2Hc3h2p6NuP+CVsQW1mU0M5Wtr3Tk4rqJ3NDuRsZ0GVNuMRpjzKmcdjfQAh8cJSJx4va2iCwTkfN9E2bF0r91bWaM6cv1vRrz4aIdnP/CXGb9nvLnA2Nq0rT9CAYey2bi+k84mF1lmkiMMVWYN43AN3kmajkfSMQ9gud/fBpVBRITHsKjQ9sy5c7exEaEcOP7vzFm0gqOZOX98cBed3PH4XRynbm8t+Y9/wRrjDGnwZsEcHwIzUHAB6q6lnIaxbMi6dQwgW//7yzu6d+cb1buZcDzc5m6usADZNWb0KTlEAZn5TFpw0TSstP8F6wxxnjBmwSw1DMs9CBguojEAqfoGlM1hYcEM+b8Vnxz91nUjgvnzo+XcdfHS0nNyHUf0GcUtx9MI9+Zx7trbKRQY0zF5k0CuBl4EOimqllAGHCjT6Oq4M6oF8dXI/tw/wWtmLkuhfNemMOXy3ejdTvSqOFZDM5x8NmGz0jN8tMk9sYY44VTJgBVdanqMlU94nl/UFVX+T60ii00OIiR5zTnh1Fn0bRGNPdOWsnNE5Zw8My7uCP1AA6rBRhjKjhvagCmGM1rxTL5jt48MuQMFm45SL8vlOioFlyUJ3y24TNSsgrpNWSMMRWAJYAyEBwk3HRWE6aPTqZd/Wo8kjaA2w7swuly8M7qd/wdnjHGFKq4oSCqF7eUZ5CVRcPEKD65tQd9LroJ8hM56yhM+v0z9maUYLhpY4zxseJqAEuBJZ7XVGAjsMmzvtT3oVVOIsKIXk2J7X8vD6XvwqlOrpz0JJtTMv0dmjHG/EFxE8I0UdWmwEzgIlWtoaqJwBBgRnkFWFlV630j9cISGOyIJj1kPoNe+5bXZm/G4QzIHrTGmArImzaAnqr6w/E3qjoV6O27kKqI0Eikx+2M2reBkCBo1HQhz0zbwCWv/cLv+20ieWOM/3mTAPaKyFgRaexZ/gHs9XVgVUK3W6gbFM6lITVIk3k8dXl99qVnc9G4+Yz7aRP5VhswxviRNwlgBFAT+NKz1PJsM6cSVR06X8+t21YiwKb8r5lxb18GtqvLcz9u5OJXF7B+n9UGjDH+4c2DYIdUdZSqdvIso1T1UHkEVyX0Gkkdp4tLw+vz1aavyNZUXh7RifF/6cKBo7lcNG4+L87ceOqJZ4wxpox5Mxx0SxF5U0RmiMjPx5fyCK5KqNYA2l3OLVtXICK8vuJ1AC5sV4cf701mSIe6vDhzE8NeXcDavTaxvDGm/HhzC2gysBwYC9xfYDHe6nMPdXIyuCamBV9v+Zol+92T2iREh/Hi8E68eW0X0jJzGfbKAp7/0WoDxpjyccoZwTwzyXQpp3iKVd4zgpWpj68ga88yLm3ajNDgcD4f+jnhweEndh/JyuPxb9cxZfkeWteJ5dkrOtKufrwfAzbGVBUlnhEM+FZE7hKRuvYkcCn0GUVUVhpjq3dj+9HtvL367T/srhYVxvNXncnb13Xl0LE8hr26gOdmbCDX4fRTwMaYqs6bGsC2Qjar5yGxclWpawCq8M75cGQnf+symBm7ZzF5yGSaJzT/06HpWfn86/t1fL50N61qx/LfKzrQIamaH4I2xlQFJa4BeJ4IPnkp9x//Sk8EBv4HMg/wQJYSHRrNYwsfw6V/vt8fHxXKs1d05L0bunEkO49LXvuF/07/3WoDxpgy5dVooCLSTkSuFJHrji++DqxKqt8Fut1M4pL3ub/5VaxIXcHkDZOLPPyc1rWYcW9fLutcn1dnbWHIy/NZuetIOQZsjKnKvOkG+k9gnGc5B3gGGOrjuKqu/g9DVA2GLp9Cjzo9eGHZCxw4dqDIw+MjQ3nm8o68f2M3MnMdXPLaAv4z9Xdy8q02YIwpHW9qAJcD5wL7VfVGoCNg3VNKKrIaXPAksmcZj8S0weFy8NTip075sX6tajH93mSu6NKA8XO2cNG4+ayw2oAxphS8SQDZquoCHCISB6QADXwbVhXX/gpokkzDeS9xZ5tr+XnXz8zcMfOUH4uLCOXpyzsw4abuZOY6uNRqA8aYUvAmASwRkWrAW7jnAVgGLPRpVFWdCAx6DvKyuG77aloltOKpxU+RkZfh1cf7tqzJ9HuTubKruzYwxGoDxpgS8KYX0F2qekRVxwPnAdd7bgWZ0qjZEs4aTejqyTzWeBgHcw7y4tIXvf54XEQo/7nMXRs4ZrUBY0wJnNacwKq6XVVXlVXhIvKuiKSIyJqyOmelcvZfIaExbee+zDWtRvDZxs9YdmDZaZ2isNrA8p2HfRSwMaYq8fek8O8DF/o5Bv8JjYRBz0LaRu7OVupF1+PRhY+S58w7rdMcrw18cFN3snIdXPb6L/x76nqrDRhjiuXXBKCqc4HAHlq6xXnQZihR819kbNtb2Za+jXdWv1OiUyW3rMk0T23gjTlbrTZgjClWkQmg4Lg/hS3lFaCI3CYiS0RkSWpqankVW74u/A8EhXD2skkMbDyQt1a/xdYjW0t0KqsNGGO8VVwNYCmwxPOaCmwENnnWl/o+NDdVfVNVu6pq15o1a5ZXseUrvj6c83fYNIO/VTuTyJBIHl34aKHDRHjreG3gqm7u2sDgl+dZbcAY8wdFJoACY/7MBC5S1RqqmggMAWaUV4ABo/vtULs9iT89wX1n/h/LU5bz+cbPS3XKuIhQ/n2puzaQnee02oAx5g+8aQPoqao/HH+jqlOB3r4LKUAFh8CQ5+HoHi7esco9TMTSF0jJSin1qZM9PYUK1gaWWW3AmIDnTQLYKyJjRaSxZ/kHsLcsCheRT3E/VNZKRHaLyM1lcd5Kq0F36Hw9sng8jzQfTr4rn38v/neZnDr2pNrA5a//wpPfr7PagDEBzJsEMAKoCXwJTPGsjyiLwlV1hKrWVdVQVU1S1ZJ1f6lKBjwKkdVoOOtp7uhwOzN3zuSnnT+V2emP1waGd2/IW/O2MfClefy2PbA7YhkTqLx5EviQqo4CzlLVzqo6WlXtF8NXoqrD+U/A7l+5Pi+YlgkteWqR98NEeCM2IpSnLmnPx7f0IN/p4so3FvLYt2vJynOUWRnGmIrPm+Gge4vIOmC9531HEXnN55EFso4joFEfQn96nEc73UtqdiovLXupzIvp07wG00cnc13PRry3YDsXvjiPRVsPlnk5xpiKyZtbQC8AFwAHAVR1JZDsy6ACnggMfg5yM2i/9FOuaXMNn234jBUpK8q8qOjwEB4b1o6Jt/VEBIa/uYiHv1rDsVyrDRhT1Xn1JLCq7jppk7Uc+lqtNtD7/2DFR/xfYjfqRNfh0V9Of5gIb/Vsmsi0Ucnc1KcJHy3ewfkvzGXB5jSflGWMqRi8SQC7RKQ3oCISKiL34bkdZHws+QGIb0jUtL8zttuDbEnfwjtrfNdOHhkWzCMXncHk23sRHhLENW8v5qEpq8nIyfdZmcYY//EmAdwBjATqA3uAMz3vja+FRcGgZyB1Pcm7VnFh4wt5a9VbbE0v2TAR3urauDo/jDqb25ObMum3nVzwwlzmbKyiw3AYE8C86QWUpqrXqGptVa2lqn9RVWspLC+tBkKrwTDnaf7W6i9EhkTy4NwHycrP8mmxEaHBPDSoDV/c2Zvo8BCuf/dX7p+8kvRsqw0YU1WIqha+Q2QcUPhOQFXv8VVQRenatasuWbKkvIv1vyO74NXu0LQfc8++k7t/upsBjQbwbN9nCRLfD+ia63Dy8k+bGD9nKzViwnjqkvac26a2z8s1xpQNEVmqql1P3l7cr8fxgeCKWkx5qdYA+j0IG34gOTOTv3b9Kz/u+JHXVpRPb9zwkGDuv6A1X93Vh2qRYdw8YQn3TlrBkSzfNEgbY8pHkTWAiihgawAAznx4IxlyM9C7FvHIkqf5avNXPJP8DAObDCy3MPIcLl6ZtZnXZm2mWlQojw9rx8B2dRCRcovBGHN6SlIDOP7BmiLyrIj8ICI/H198E6YpUnAoDH4e0nchP/+Lh3s+TOdanXl4wcOsTl1dbmGEhQQx5ryWfHP3WdSNj+Suj5dxx0dLSTmaU24xGGPKhjc3kD/G3e2zCfAYsB34zYcxmaI06gU97oTF4wlb8SkvnPMCNSJrMGrWKA4cO1CuoZxRL44v7+rNQwNbM3tDKgOen8Nnv+2iMtUojQl03iSARM8gbfmqOkdVbwL6+zguU5Tzn4Bm58L3Y6i+bw3j+o/jWP4x7pl1D9mO7HINJSQ4iNv7NmPqqLNpXTeOB75YxbXv/MrOg77toWSMKRveJIDj/f72ichgEekElNuUkOYkwSFwxXuQ2BwmXUsLp/BM8jOsP7iesfPHlmoWsZJqWjOGibf25ImL27Fi1xEueHEu78zfhtNltQFjKjJvEsATIhIP/BW4D3gbuNenUZniRcTD1ZMgKBg+uYK+ie0Z02UMM3bMYPzK8X4JKShI+EvPRsy4N5lezRL513fruOz1X9h4oOxGMTXGlC1vHgT7TlXTVXWNqp6jql1U9ZvyCM4UI6ExDP8E0nfDZ9dxfasRXNz8Yl5f+TrTtk/zW1j1qkXyzvVdeWn4mew4eIzBL8/jpZmbyHOUf83EGFM8b3oBTRCRagXeJ4jIu74Ny3ilYU8Y9ipsn4f8MIaHe4ylc63OjJ0/ljVpa/wWlogw7Mz6zBzTl4Ht6vLCzI1cNG4+K3cd8VtMxpg/8+YWUAdVPfEvV1UPA518F5I5LR2udA8at/wjwhaP5/l+z5MYkcg9P99T7j2DTpYYE87LIzrx9nVdSc/O55LXFvDk9+vIzrPBZI2pCLxJAEEiknD8jYhUB0J8F5I5bf0egraXwI//JHH7Qsad67+eQYUZcEZtZoz53zSUF7w4l1+22FDTxvibNwngOWChiPxLRJ4AfgGe8W1Y5rQEBcHFr0P9zjDlVlrmZPN08tOsP7iehxc8XCH65sd5pqH89Fb3xDNXv7WYh6asssHljPEjbxqBPwAuBQ4A+4FLVfVDXwdmTlNoJAz/FCKrwyfD6RffitFdRjN9+3S/9QwqTK9m7olnbktuyqTfdjHg+Tl8v2pfhUhSxgQabxqBmwFbVPUVYA0woGCjsKlAYmvD1RMh9yh8OpwbW1zJ0GZDeW3la0zfPt3f0Z0QGRbM3we14euRZ1ErNpyRnyzj5glL2H3YHiAzpjx5cwvoC8ApIs2BN4AGwCc+jcqUXJ32cNk7sG8l8tUd/LPHw5xZ80zGzh/L2oNr/R3dH7RPiufrkX0YO7gNC7cc5Lzn5/L2vK04nNZl1Jjy4E0CcKmqA/dtoFdU9X6grm/DMqXS6kK44ElY/y1hc57mxXNeJCEigXt+uoeUrBR/R/cHIcFB3HJ2U34c436A7Inv13PxawtYvTvd36EZU+V5NRSEiIwArgO+82wL9V1Ipkz0vAu63ADznydxw3TG9R9HRn4Go34eRY6j4o3cmZQQxTvXd+XVqztz4Gguw16dz+PfruNYrsPfoRlTZXmTAG4EegFPquo2EWkCWCNwRScCg56FJsnwzT20OprG02c/zdqDa3lkwSMVstFVRBjcoS4zx/RlRPeGvLtgG+c9P4eZ6/z7PIMxVZVNCFPVZR+Gt8+DrDS45Sfe2TubF5e9yLVnXMt9Xe8rlyklS2rpjkM8NGU1Gw9kMrBdHR4d2pbacRH+DsuYSqfEE8KYSi4ywT1wHMAnV3FTs0u4uvXVfLjuQx6c9yB5zoo7rWOXRtX57v/O5v4LWvHT7ykMeG4OHy7agctGGTWmTFgCCASJzeCqj+DwduTzG3mwy18Z3Xk0U7dN5Y6Zd3A076i/IyxSWEgQI89pzvTRybRPiufhr9Zw+fhf2LDfRhk1prROKwGISB1fBWJ8rPFZcNGLsHU2MvUBbm53E/8++98sT1nO9VOvZ/+x/f6OsFhNakTz8S09eO6KjmxLc48y+sy038nJt3GFjCmp060B/OCTKEz56PQX6DMalr4Hs55iSJPBjB8wnv3H9nPN99ew4dAGf0dYLBHhsi5J/PTXfgw7sz6vzd7CBS/OZfaGitW11ZjK4nQTgPgkClN+zv0ndLwa5j4DX9xMjxodmDBwAghcP+16Fu1b5O8IT6l6dBjPXdmRT27pQbAIN7z3G3d8uJQ9R/w/8J0xlcnpJoC3fBKFKT9BQXDxa+5EsGYKvDeQlsGxfDzoY+pG1+XOH+/k2y3f+jtKr/RuXoOpo92NxLM3pnDuc7N5ddZmm3zGGC9ZN9BA9vv38MWtEBEHwz8mo2ZLRs8aza/7f2VU51Hc3O5mRCpHpW/34Sz+9d06pq89QNOa0Tw+tB1ntajh77CMqRCsG6j5s9aD4eYZEBQK7w0iduMMxg8Yz+Cmg3lp2Us8segJHK7K8SRuUkIUb1zblfdu7IbTpfzlncWM/HgZ+9LttpAxRbEagIHMVPjsWti5EJIfwNX3b7y84hXeWfMO/ZL68XTy00SFRvk7Sq/l5Dt5a+5WXpm1meAgYdS5LbixTxPCQuzvHROYiqoBnDIBiEg0kK2qLhFpCbQGpqpquc/kYQnAhxy58N0YWPERtBkKl4xn4tZv+fev/6ZtYlvG9R9HYmSiv6M8LbsOZfHYt+uYuf4AzWvF8PiwtvRuZreFTOApzS2guUCEiNQHZgDXAu+XbXjG70LCYdgrcP6T8Pt38O4FDK97Fi/0e4FNhzdx7dRr2XF0h7+jPC0Nqkfx9vVdeef6ruQ6nFz91mLu+XQ5B45WvMHwjPEHbxKAqGoW7uGgX1PVK4C2ZVG4iFwoIhtEZLOIPFgW5zSlIAK974YRk+DQdnjzHPpLLG9f8DaZeZlc+8O1rExd6e8oT9u5bWrz4719GXVuC6at3U//Z2fz9ryt5Nu8AybAeZUARKQXcA3wvWdbcGkLFpFg4FVgIHAGMEJEzijteU0ZaHk+3DITwqLg/cF03LueDwd9SExYDLdMv4Wfd/7s7whPW0RoMPee15If702me5PqPPH9eoa8PJ/FWw/6OzRj/MabBDAKeAj4UlXXikhTYFYZlN0d2KyqW1U1D5gIDCuD85qyUKs13DoLGnSHL2+n0eJ3+fCCCTSv1px7Z9/LxN8n+jvCEmmUGM27N3TjzWu7kJnr4Ko3F3HvpBV2W8gEpCIbgUUkxDMTmG8KFrkcuFBVb/G8vxbooap3n3TcbcBtAA0bNuyyY0flug9d6Tnz4Yf73cNHtBpE1tCXeGDR48zZPYdhzYZxX9f7qBZROaeIzs5z8trszbwxZyvBQcJd/Zpxa3JTIkJLXcE1pkIpSSPwrwU+PM4nUXlBVd9U1a6q2rVmzZr+CiNwBYfCkBdg4H9h43SiJlzMi2eO4db2t/L91u8Z+tVQvtnyTYWcYOZUIsOC+ev5rZg5pi/9WtXkuR830v/Z2Xyzcm+lvB5jTldxCaDgI6B9fFD2HkKM+5kAABiYSURBVNwTzB+X5NlmKhoR6HEb/OVzOLqbkHfO457Ebky6aBIN4xryj/n/4NYZt7I9fbu/Iy2RholRvP6XLky8rScJ0WHc8+lyLnv9F5bvPOzv0IzxqeISgK//BPoNaCEiTUQkDBgOfOPjMk1pNOsPt/wMkdVgwkW0XP01Hwx4g4d7Psy6g+u47JvLGL9yfIWeZKY4PZsm8s3dZ/HMZR3YeSibS177hXsnrbCniU2VVVwbQBawGXdNoJlnHc97VdUOpS5cZBDwIu5eRe+q6pPFHW8PglUQ2Yfhm3tg/TcQlwTnPkJa83N4ZsmzTN0+lSbxTXik5yN0rfOnW46VRmaug9dnb+atedsIErijbzNuT25GZJi1D5jK57SfBBaRRsWdUFXLvTXWEkAFs30+TP877FsJ9TrBBU8xP8TFE4ueYE/mHi5pfgljuoyptI3E4H6a+D/Tfuf7VfuoGx/B3y5szdCO9QgKqhyD5BkDpRgKosAJEoFkYKeqLi3j+LxiCaACcrlg9Wcw8zHI2AttLiL7nL/zxq4ZTFg7gdiwWO7vdj9Dmg6pNCOLFua37Yd4/Nt1rN6TzpkNqvHwkDPo0ijB32EZ45WS1AC+Ax5U1TUiUhdYBizBfTvoTVV90ZcBF8YSQAWWlwULX4X5L4AzD7rfysaOl/L48pdYmbqSHnV6MLbnWBrHN/Z3pCXmcilTlu/hmWm/k5KRy9CO9fjbwNbUrxbp79CMKVZJEsBaVW3rWf870FpVrxORWGBBWbQBnC5LAJVAxn6Y9SQs/wjC43AlP8DnCdV5cfk4cp253NrhVm5qdxNhwWH+jrTEjuU6GD9nC2/O3QrA7clNub1vM6LDQ/wcmTGFK0kCWKGqZ3rWfwLeUtWJJ+8rT5YAKpH9a2DGWNg6C6o3Ja3v/TxzZDlTt0+rEo3E4J6E5ulpG/h25V5qxYYzekBLruiaRGiwDTttKpaSJIBvcY/+uRt4F2iiqkdEJBJYcrx2UJ4sAVQyqrB5Jkz/B6RtgEZ9mN/5Sp7YPJE9mXu4uPnFjDxzJHWi6/g70lJZuuMQT/3wO0t3HKZJjWjGnNeSwe3rWkOxqTBKkgBqAY8DdYFXVXWGZ/s5QBdVfdaH8RbKEkAl5XTAsgkw6ynISiO7/eW8Wa8J72/6AkXp37A/I1qPoGvtrpW2oVhV+Wl9Cv+dvoENBzJoWy+OBy5sTXKLGpX2mkzVUZIEEAHEqWrKSdtrAhmqWu6jZ1kCqORyjsL852HhayDC3m43MDEujik7ppKem06LhBZc3fpqBjUZVKlmICvI6VK+XrGH53/cyO7D2fRsWp0HLmxN54bWY8j4T0kSwJvANFWdctL2S4DzVfVOn0RaDEsAVcSRne5uo2s+Bwkmp1l/fkhqxSfp69lwZBOxYbFc2vxSrmp9FQ1iG5z6fBVQrsPJp4t3Mu7nzRw8lsf5Z9Tmvgta0bJ2rL9DMwGoJAlgqap2KWLfWmsDMKWWsh5WfgqrJkPGXjQ8juUt+/FJZDAz01bgUhd9k/oyos0IetXtVSlvpRzLdfDu/G28MXcrWXkOLu2cxOgBLUhKqJw1HFM5lSQBrFfVNqe7z5csAVRRLidsnwcrJ8G6ryH/GAeqNeSzhm34PHcvh/LSaRzXmBGtRzCs+TCiQ6P9HfFpO3Qsj9dnb2bCwh2gcE3Phow8pzk1YsL9HZoJACVJAHOA+1X115O2dwOeU9Vkn0RaDEsAASDvGPz+PaycCFtnkacupiedwSexUazJSSE6NJphzYYxovWISvlQ2d4j2bw0cxOTl+4iMjSYW85uyi1nNyE2ItTfoZkqrCQJoDvwGe4J4I8P/dAVuA4YrqqLfRNq0SwBBJiM/bB6srtmcGA1qyMi+aRec6aRiUOd9KnXhxGtR9C7Xm9CgyvXD+jmlEyem7GBqWv2Uz06jJHnNOeaHg1tMhrjEyUaC8jTFXQk0M6zaS3wysk9g8qLJYAAdmCtu1awejJpxw7weUINJsfHkaJ5RIZE0qNOD3rX702fen1oGNfQ39F6beWuI/x3+gbmb06jXnwEd/RrxpVdG1giMGWq1IPBFTjRWcAIVR1ZVsF5yxKAweWEbXNg5STy13/LghAX86NjWBATy25xz2CaFF2PPkln06deH7rX7V4p2gzmb0rjhZkbWbrjMDVjw7k9uSlX92hIVJgNL2FKr1QJQEQ6ASOAK4FtwBRVLfdpIi0BmD/IOwabfoSdi2DXYnamrWVBeCi/REawOCqKbIEQCeLMhNb0aTiA3kl9aF29NUFSMYdqUFUWbj3IKz9v5pctB0mICuWWs5tyba9GxFkbgSmFkrQBtMT9oz8CSAMmAfeparHzBPiSJQBTrLws2Lscdi0ib+ciVqQsY0FQPr9ERvJ7uHvwuepBEfSu0YHeTQfRu2E/EiMT/Rx04ZbuOMQrP29m1oZU4iJCuKFPE27q05hqUZV3ED3jPyVJAC5gHnCzqm72bNuqqk19GmkxLAGY06IKBzfDrsWk7ZjLL/uXsMB5hIWRERwOdt9jbxMSR6/E9rSp1ZEW9XrSsFY7QoMqzl/bq3en88qsTUxfe4DosGD+0qsRt57d1LqPmtNSkgRwMe55evsA04CJwNuq2sSXgRbHEoAptaxDuHb9yvqt01lwYAkL8lJZGRaC0/OQWagqTTSEFqHxtIiuT4vqrWhRuzN16nZGYutCkH9uH23Yn8Grszbz3aq9hIUEMaJ7Q25Pbkad+Ai/xGMqlxK3AYhINDAM962g/sAHwJfHB4crT5YATJlz5pO7fzXb9v3KptQ1bDq6nU05KWxyZXOgwG99jMtF83wHLYKiaBFegxbxTWhRox3xNdtAQhOo1gBCfP9X+dbUTF6bvYUvl+8hWIQruiZxR99mNKhuTxabopVJLyARSQCuAK5S1XPLMD6vWAIw5Sk96yBb9ixk0/6lbDq8gU2Ze9iUn04GzhPH1HI4aJGXT9P8fGoTSu2QKGqFxVMrIpFaUbUJi64BUYmFL5HVIKhk3T13Hcri9Tlb+HzJblyqXNypPiPPaU6TGhW/x5Mpf2XWDdSfLAEYf1NVDmQdYPPhTWw6sJxNqavZlL6N7blp5KjzT8dXc7qo5XBQy+mkttNJLYeTWk6H59VFrZAYEiKrI1E1IDIBwqIgNArCoj2vURAafdLr//an5AbzwdJUPlmWRrozlIHt63PzWU3oZKOPmgIsARjjQ6rK0byjpGSlnFgOZB1wr2fuI+XYflKyUzmUdxTlj//mQhFqaRCJLoh1uYhzOohxOoh15BHrdBLjchFbYIlx6Yn1KFUKDpGXq6HkEYIGhRIWFkZ4RAQSHAbBYRAUCsGh7vXgk9cL7A8KcddMJNjzGuRe/rAt2N0eIkEnbStwvAgg7lcJ+t86UsR+CtmP+33BfSfWvdhX5ACCRWwvdHMZDEJYFgMZ1ukI0SXrtVZUArCnTIwpAyJCfHg88eHxtEhoUeRx+a580rLSTiSH1OzUE+tp2Wkcyctkd34mGXkZZORlkO/KL7bcIISYoDBig0KJkRAiCcKZ5yI/x4k4XESKUCM8lMTwYKKAcFXCNJsI1zHC8lxEuJyEuVyEuxyEO91LmMtBsMtFiDrdry4nwep5dTkJURfBqoQAIaoEAxXzyYoq5povoMWAMj2lJQBjylFoUCh1Y+pSN6buKY9VVXKduWQWSAiZeZkczT9KZl6mez3v6In9mXmZHHMcw+XMRR25HMrO4nB2FrnOXCQ/m6DgfPR4+0UQRfxqh1CSnwVBCAkKIUSCCZZggoOCCCIIESEIIUiCEIQgcS/i3utZ/+P+YM+6eM4LnkrCiZLE814KbKPw7fLHGN3/YSl8e8HrOX6AaqH7izy+UN6d41TujIqhrMfgtwRgTAUlIkSERBAREkGNyBolPs/GAxm8t2A7U5btJtfhoFfzOEZ0r0u3pnE4XPnkOHPIc+aR48wh15lLvjMfhzpwupw41YnD5cDhcuBUJ06XE4ee9N7l+NPxTnXiUheqigvPq7rc2yiw7tl//P3xbU51nrhVpiju/1dUleP/d2Jbge0u/d9nCt7eLnjb7Q+34P6w6vms6okdKvrHgwpx8i29P+0vo9vs+aFl3+XX2gCMCRCHj+Xx6W87+eCXHew/mkPjxChu6N2Yy7s2ICbc/hasyqwR2BgDQL7TxdQ1+3lvwTaW7zxCbHgIV3ZrwPW9GtMw0Z4nqIosARhj/mT5zsO8t2A7P6zeh1OVAW1qc0PvxvRqmkhQUOWbgtMUzhKAMaZI+9Nz+HDRdj5ZvJPDWfkkJURyeZckLuucZE8ZVwGWAIwxp5ST72T62v1MXrKbBVvSUIXezRK5smsDLmhbh8gwm6imMrIEYIw5LbsPZ/HF0j18vmwXuw5lExsewpCO9biiaxKdGlQ7RddHU5FYAjDGlIjLpSzedojJS3bxw5p95OS7aF4rhiu6JHFJ5/rUirURSSs6SwDGmFLLyMnn+1X7+GzJLpbtPEJwkHBOq5pc3qUB/VvXIizEngmuiCwBGGPK1OaUTD5fupspy3aTkpFL9egwLulUnyu6JtG6Tpy/wzMFWAIwxviEw+li7qZUJi/Zzcz1B8h3Km3qxnFh2zoMbF+HFrVirL3AzywBGGN87tCxPL5avocfVu9j6c7DqELTGtFc0K4OF7atQ4ekeEsGfmAJwBhTrlKO5jBj3QGmrdnPwq0HcbqUevERnN+2DgPb1aFr4+oE28Nm5cISgDHGb45k5TFzfQrT1uxn7qZU8hwuEqPDOL9tbS5oW4fezWpYA7IPVagEICJXAI8CbYDuqurVr7olAGMqv2O5DmZvSGXa2v38vP4Ax/KcxEaEcG7rWlzYri59W9a0B87KWEWbEGYNcCnwhp/KN8b4SXR4CIM71GVwh7rk5DtZsDmNaWv28+P6A3y1Yi8RoUH0a1mL/q1r0bt5IkkJNhSFr/glAajqejjVJArGmKouIjSYc9vU5tw2tXE4Xfy67RDT1u5n+tr9TFu7H4CG1aPo3SyRXp7FHjwrO35tAxCR2cB9xd0CEpHbgNsAGjZs2GXHjh3lFJ0xxl9UlU0pmSzYnMYvWw6yaOtBMnIcALSoFUPvZon0bl6Dnk0SiY8K9XO0FV+5twGIyEygTiG7/qGqX3uOmc0pEkBB1gZgTGByupS1e9P5ZctBftlykN+2HSI734kItKsXf6KG0K1xdaJtcps/qVCNwCcKtwRgjCmBPIeLlbuPnKghLN95mHynEhIknNmg2okaQsekatagjCUAY0wVlp3nZMmOQydqCKt3H8GlECTQrGYMbevF0bZe/InXQLttVKESgIhcAowDagJHgBWqesGpPmcJwBjjjaM5+fy69RCr9qSzdk86a/ceZf/RnBP7kxIiaVsvjnb14mlb350UasWGV9mOKRUqAZSUJQBjTEmlZeaydu9R1u51J4R1e4+yLe3Yif01YsL+UEtoVz+OhtWjqkRSqGjPARhjTLmqERNO35Y16duy5oltGTn5rN+Xwdq96azZ404OCzan4XC5/zCODQ+hcY1okhIiPUvUidf6CZHEVPIG58odvTHGlEJsRCjdm1Sne5PqJ7bl5DvZdCCTNXvTWbf3KDsPZbHxQAY//55CrsP1h88nRIUWSAqVL0FU7OiMMaacRYQG0z4pnvZJ8X/YrqqkZeax+3AWuw9nexb3elEJolpUKPWrRZIQFUZMeAixESHERIQQG+55jQglxrMeFxFCTHgoMREh7m3hIT4fLM8SgDHGeEFEqBkbTs3YcDo1TPjT/qISxJ4j2RzNziclI4eMHAeZOQ4y8xx40/waHRZ8IlE8eXE7ejRNLNNrsgRgjDFl4FQJoiCXS8nKd5KRk09mjoOMXHdiyMhxkJmb73l1/C9h5DqIjSj7rquWAIwxppwFBcmJ2zzEn/p4n8Xhv6KNMcb4kyUAY4wJUJYAjDEmQFkCMMaYAGUJwBhjApQlAGOMCVCWAIwxJkBZAjDGmABVqYaDFpFUoCwnBa4BpJXh+fzJrqXiqSrXAXYtFZW319JIVWuevLFSJYCyJiJLChsjuzKya6l4qsp1gF1LRVXaa7FbQMYYE6AsARhjTIAK9ATwpr8DKEN2LRVPVbkOsGupqEp1LQHdBmCMMYEs0GsAxhgTsCwBGGNMgKryCUBEGojILBFZJyJrRWRUIceIiLwsIptFZJWIdPZHrKfi5bX0E5F0EVnhWR7xR6ynIiIRIvKriKz0XMtjhRwTLiKTPN/LYhFpXP6RFs/L67hBRFILfCe3+CNWb4lIsIgsF5HvCtlX4b+T405xHZXmOxGR7SKy2hPnkkL2l/j3KxBmBHMAf1XVZSISCywVkR9VdV2BYwYCLTxLD+B1z2tF4821AMxT1SF+iO905AL9VTVTREKB+SIyVVUXFTjmZuCwqjYXkeHA08BV/gi2GN5cB8AkVb3bD/GVxChgPRBXyL7K8J0cV9x1QOX6Ts5R1aIe+Crx71eVrwGo6j5VXeZZz8D9P4j6Jx02DPhA3RYB1USkbjmHekpeXkul4Plvnel5G+pZTu6RMAyY4Fn/HDhXRKScQvSKl9dRaYhIEjAYeLuIQyr8dwJeXUdVUuLfryqfAAryVFc7AYtP2lUf2FXg/W4q+A9rMdcC0MtzS2KqiLQt18BOg6eKvgJIAX5U1SK/F1V1AOlAYvlGeWpeXAfAZZ7q+eci0qCcQzwdLwIPAK4i9leK74RTXwdUnu9EgRkislREbitkf4l/vwImAYhIDPAFMFpVj/o7ntI4xbUswz3uR0dgHPBVecfnLVV1quqZQBLQXUTa+TumkvDiOr4FGqtqB+BH/vcXdIUiIkOAFFVd6u9YSsPL66gU34nHWaraGfetnpEiklxWJw6IBOC5N/sF8LGqTinkkD1Awb8AkjzbKpxTXYuqHj1+S0JVfwBCRaRGOYd5WlT1CDALuPCkXSe+FxEJAeKBg+UbnfeKug5VPaiquZ63bwNdyjs2L/UBhorIdmAi0F9EPjrpmMrwnZzyOirRd4Kq7vG8pgBfAt1POqTEv19VPgF47k++A6xX1eeLOOwb4DpPa3pPIF1V95VbkF7y5lpEpM7xe7Ii0h33d1zR/oEiIjVFpJpnPRI4D/j9pMO+Aa73rF8O/KwV7MlFb67jpPuxQ3G33VQ4qvqQqiapamNgOO7/3n856bAK/514cx2V5TsRkWhPhw9EJBo4H1hz0mEl/v0KhF5AfYBrgdWe+7QAfwcaAqjqeOAHYBCwGcgCbvRDnN7w5louB+4UEQeQDQyvaP9APeoCE0QkGHeS+kxVvxORx4ElqvoN7mT3oYhsBg7h/sdc0XhzHfeIyFDcvbgOATf4LdoSqITfSaEq6XdSG/jS8zddCPCJqk4TkTug9L9fNhSEMcYEqCp/C8gYY0zhLAEYY0yAsgRgjDEByhKAMcYEKEsAxhgToCwBmCpDRJyeERPXiMhkEYkq5tihIvJgCcqYLSJdPes/HH8GoIhjRxcXgzH+ZgnAVCXZqnqmqrYD8oA7ijpQVb9R1f+cvN3zdKtXVHWQ5+nfoowGLAGYCssSgKmq5gHNReQicY9bv1xEZopIbTgxHvwrnvX3RWS8iCwGnil4EhGJFJGJIrJeRL4EIgvs2y4iNTxPa37vGYBvjYhcJSL3APWAWSIyy3P86yKyRE6aN8BznsdEZJm4x31v7dkeIyLvebatEpHLPNvPF5GFnuMne8aGMua0BcKTwCbAeP6KHwhMA+YDPVVVxT3pxwPAXwv5WBLQW1WdJ22/E8hS1TYi0gH3YHsnuxDYq6qDPeXHq2q6iIzhj+O4/0NVD3meGv5JRDqo6irPvjRV7SwidwH3AbcAD+N+rL+957wJnnGdxgIDVPWYiPwNGAM8ftr/oUzAswRgqpLIAkNkzMM9bEErYJJn7JcwYFsRn51cyI8/QDLwMoCqrhKRVYUcsxp4TkSeBr5T1XlFlHGluIfzDcE9hMQZwPHzHR/YbylwqWd9AAWGWlDVw+Ie6fIMYIFneIAwYGER5RlTLEsApirJ9gzLfIKIjAOeV9VvRKQf8GgRnz1W0kJVdaO4p+EbBDwhIj+p6h/+IheRJrj/su/m+SF/H4gocMjxkSmdFP/vUnDPOTCipPEac5y1AZiqLp7/DY17fXEHFmEucDWAuMf573DyASJSD/dtoo+A/wLH52TNAGI963G4k0y6px1ioBdl/wiMLFBOArAI6CMizT3bokWkZQmuyxhLAKbKexSYLCJLgaLmVC3O60CMiKzHfZ+9sElG2gO/em4//RN4wrP9TWCaiMxS1ZXActxDRX8CLPCi7CeABE/D8krc7QmpuEeu/NRzO2oh0LoE12WMjQZqjDGBymoAxhgToCwBGGNMgLIEYIwxAcoSgDHGBChLAMYYE6AsARhjTICyBGCMMQHq/wHek3MYxH7AkAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3rdaI0VKNVT3"
      },
      "source": [
        "### Training\n",
        "\n",
        "We create a neural network and train it on the ACSF database we just constructed."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1BZlNPamNVT4"
      },
      "source": [
        "# setup the neural network\n",
        "# the network uses tanh function on all hidden neurons\n",
        "\n",
        "nn = MLPRegressor(hidden_layer_sizes=(5,), activation='tanh')"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "clZbfbDkNVT4"
      },
      "source": [
        "The fitting may not be trivial since our database is small... the next instruction can be executed multiple times let the NN train more and hopefully improve."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2KcIPVQqNVT4",
        "outputId": "2cc4ee1b-6fc1-4597-f2a8-8c55db169f32",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "source": [
        "# change some training parameters\n",
        "nn.set_params(solver='lbfgs', alpha=0.001, tol=1.0e-10, learning_rate='constant', learning_rate_init=0.01)\n",
        "# do some training steps\n",
        "nn.fit(acsf, energy);"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/sklearn/neural_network/_multilayer_perceptron.py:470: ConvergenceWarning: lbfgs failed to converge (status=1):\n",
            "STOP: TOTAL NO. of ITERATIONS REACHED LIMIT.\n",
            "\n",
            "Increase the number of iterations (max_iter) or scale the data as shown in:\n",
            "    https://scikit-learn.org/stable/modules/preprocessing.html\n",
            "  self.n_iter_ = _check_optimize_result(\"lbfgs\", opt_res, self.max_iter)\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ujgh8CENNVT5",
        "outputId": "408d2e6b-568f-4537-91dc-e16059ce954f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 558
        }
      },
      "source": [
        "# evaluate the training error\n",
        "energyML = nn.predict(acsf)\n",
        "\n",
        "print (\"Mean Abs Error (training)  : \", (numpy.abs(energyML-energy)).mean())\n",
        "\n",
        "# energy curve\n",
        "plt.plot(dists, energy,'-.' )\n",
        "plt.plot(dists, energyML,'o' )\n",
        "plt.xlabel('Pair distance')\n",
        "plt.ylabel('Energy')\n",
        "plt.show()\n",
        "\n",
        "# regression plot\n",
        "plt.plot(energy,energyML,'o')\n",
        "plt.plot([-2.3,-1.7],[-2.3,-1.7]) # perfect fit line\n",
        "plt.xlabel('correct energy')\n",
        "plt.ylabel('NN energy')\n",
        "plt.show()"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Mean Abs Error (training)  :  0.0007966465027036246\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXzU1b3/8dcnC5NAJmGbIBAQFMEFcYvoda3LrUuL3nqtS1ervV6rLWqrtlZvq9y293a1tb3VHw9bbatVsRWXuou2tdYtICAI7iKbJGwhCSHr5/fHTCDEmckkmcx3knk/H488MjPfM9/v58uQ72fOOd9zjrk7IiIiieQFHYCIiGQ3JQoREUlKiUJERJJSohARkaSUKEREJKmCoAPoD6NHj/ZJkyYFHYaIyICxcOHCje4eibdtUCaKSZMmUVVVFXQYIiIDhpmtSrRNTU8iIpJUIInCzD5tZsvNrN3MKpOUuzJWbpmZ3W1mRZmMU0REgqtRLAPOAv6eqICZjQdmA5XuPh3IB87LTHgiItIhkD4Kd18BYGbdFS0Ais2sBRgKrOvn0EREpIus7aNw97XAT4APgPVArbs/mai8mV1sZlVmVlVTU5OpMEVEBr1+SxRm9nSsb6Hrz5kpvn8EcCYwGRgHDDOzzyUq7+5z3b3S3Ssjkbh3eImISC/0W9OTu5/cx12cDLzn7jUAZnY/cBRwZ19ji2vpPFgwB2rXQFkFnPQdmHFOvxxKRGQgyeZxFB8AR5rZUKAROAnon8ERS+fBw7OhpTH6vHZ19DkoWYhIzgvq9thPmdka4F+AR8zsidjr48zsUQB3fwn4E7AIeC0W69x+CWjBnF1JokNLY/R1EZEcF9RdT/OB+XFeXwec3un5d4Hv9ntAtWt69rqISA7J2rueMqqsomevi4jkECUKiHZcFxbv/lphcfR1EZEcp0QB0Q7rWTezdcgetLtB2QSYdbM6skVEyO67njJrxjncs+Uw/vexlSy/9BSGhfRPIyICqlHsJlISAqCmringSEREsocSRSeRcCxR1CtRiIh0UKLopLw0miiqtylRiIh0UKLopDxcxL57hCnI73ZWWxGRnKEe205GDhvC41ccF3QYIiJZRTUKERFJSomii2v+tISr71sSdBgiIllDTU9dlIe1LLeISGdKFF1cdcq0oEMQEckqanoSEZGklCi6mP/qGg74zuNU1+0IOhQRkaygRNFFcWE+Dc1tmsZDRCRGiaKLSKwzu1qJQkQEUKL4iPKwJgYUEelMiaKLiBKFiMhulCi6KCrMJ1xUoEQhIhKjRBFHJBxSohARiVGiiKM8HNLtsSIiMUoUcUTCRapRiIjEaAqPOI6ZMorRJUOCDkNEJCsoUcRx7uETgw5BRCRrqOkpgfZ2p63dgw5DRCRwShRxvPDOJva5/jEWrtoSdCgiIoFToohjz1FDueT4vXaO0hYRyWXqo4hj3PBirj5l36DDEBHJCqpRJLB1ezMb63WLrIiIEkUCp/78OX742MqgwxARCZwSRQLlpSFqVKMQEVGiSCRSEqJ6mxKFiIgSRQKRsGoUIiKgRJFQeTjEpvomDboTkZynRJFAJByi3WFTg2oVIpLbAkkUZvZpM1tuZu1mVpmk3OVmtixW9opMxqiV7kREooKqUSwDzgL+nqiAmU0H/gOYCRwEfNLMpmQmvOhU4wDVShQikuMCSRTuvsLd3+im2H7AS+6+3d1bgb8RTS4ZUa4ahYgIkN19FMuAY81slJkNBU4HJiQqbGYXm1mVmVXV1NT0+eCRcIirT5nG9HFlfd6XiMhA1m9zPZnZ08AecTZd5+4Pdvd+d19hZj8EngQagMVAW5Lyc4G5AJWVlX2+VamoMJ/LTshYS5eISNbqt0Th7ienYR+/AX4DYGY/ANb0dZ89Ub1tB40tbew5algmDysiklWyevZYMyt392ozm0i0f+LITB5/9j2v0tbu3HfJUZk8rIhIVgkkUZjZp4BfAhHgETNb7O6nmNk44DZ3Pz1W9M9mNgpoAS5z962ZjHP2iftk8nAiIlnJ3AffyOPKykqvqqoKOgwRkQHDzBa6e9xxbdl811Pgqut28MzKDTS1JuxDFxEZ9JQoknjuzY1ceEcV67fuCDoUEZHAKFEkUV4aHXSn0dkiksuUKJLQfE8iIkoUSZXH5nuqqVPTk4jkLiWKJIYXF1KQZ2p6EpGcpkSRRF6eMbokpKYnEclpShTdKC8NqUYhIjlNiaIbEdUoRCTHKVF0IxIOUVOvRCEiuSurJwXMBl88ahJnHDwu6DBERAKjRNGN/caWBh2CiEig1PTUjc0NzTyydL36KUQkZylRdOO9jQ1c9sdFLFtXG3QoIiKBUNNTN/YfW8rjVxzLniO1yp2I5CYlim4UD8ln3z3UTyEiuUtNTymYV7WaZ1ZuCDoMEZFAKFGk4P/97R3+vHBt0GGIiARCiSIFkXCIas0gKyI5SokiBZFwkW6PFZGcpUSRgvKwJgYUkdylRJGCSDjE9uY2Gppagw5FRCTjlChSECnRkqgikruUKFJQXhpNFGp+EpFcpESRgkhYNQoRyV1KFCkoDxcBUKNbZEUkB2kKjxQMLy7k71efsLMJSkQklyhRpCAvz5g4amjQYYiIBEJNTyl64NW13P3yB0GHISKScUoUKXp4yTrufWV10GGIiGScmp5SdMvnDqMw34IOQ0Qk45QoUjSkQJUvEclNuvqlaOGqLVx572I2NzQHHYqISEYpUaSopq6J+a+uZd3WxqBDERHJKCWKFO0cnV2v0dkikluUKFJUrmk8RCRHBZIozOzHZrbSzJaa2XwzG56g3Klm9oaZvW1m38p0nJ1pvicRyVUpJQoz+6mZHZDG4z4FTHf3GcCbwLVxjpkP/B9wGrA/cL6Z7Z/GGHqkqDCfcFGBEoWI5JxUaxQrgLlm9pKZXWJmZX05qLs/6e4dqwC9CFTEKTYTeNvd33X3ZuAe4My+HLevyrV2toik29J5cNN0uGF49PfSecHsI4mUEoW73+buRwNfACYBS83sj2Z2QhpiuBB4LM7r44HOQ6HXxF6Ly8wuNrMqM6uqqalJQ1gfFQmHVKMQkV36eoFeOg8eng21qwGP/n54ds/2k459dCPlAXexpqB9Yz8bgSXA183sP939vDjlnwb2iLOr69z9wViZ64BW4K5exL4bd58LzAWorKz0vu4vnki4iNfWbO2PXYtITy2dBwvmQO0aKKuAk74DM87J7Psfng0tsVvma1fjD81my/Zmtu79b+xoaaeptW3n72P3iZCfZyxctYW3q+s49/CJ0eO3dLnlvqWRTQ9dz7WL96K5rZ2WtnYK8vL43YUzAZjz8Ou8VV3HHy46AoCND17H6LaP7oMFc3p2PkmklCjM7CZgFrAA+IG7vxzb9EMzeyPee9z95G72eQHwSeAkd493YV8LTOj0vCL2WmDGhEO82t4vOUhEeiLORZqHZ0cfd3Nx3NHSRv0rdzPq2auwTu9vffBr/PPNGpaPPpXG5la2N7exvaWNxuY2tje3csMZBzC2rJg/L1zD/z37Nk/n3Uhel4u8tTay/dHvcmLziI+GfMPHKS0q5InlH/L7F96PJoraNXFjHNFazQebtzOkII/C/DxKQrsafyaOLCav02xCo9oStKAk2HdvpFqjWApc7+4NcbbN7OlBzexU4BrgeHffnqDYK8A+ZjaZaII4D/hMT4+VTtd9Yj+u/2Rg/ekig0cfvs1vb26l4IkbGBLnm/iWh69nzoqp1O1opb6phfqmVup3tFLf1MrN5x3CUVNG88zKamY8/l0sb/f3F7TtYK+lP+MLzZMwg+LCfIYOyad4SD7Fhfk0NrcBMHLYEPYfV4q9Gf976/i8Tfz83IMpKswjVJBPKPa7uDAfgK+eOIWLj9srWrisItZktLu8sgoev+K4uPu/4OjJuz23BPugLF7Xb++kmiiWANPMdpsUrxZY5e61vTjur4AQ8FRsny+6+yVmNg64zd1Pd/dWM/sq8ASQD/zW3Zf34lhp0+X8RXJTGpps/KHZWOtHv80vGv6vbN3eQm1jS7QJp9Pja0/bl3MPn8g71Q0cUL8W4vw5lrVUU7VqMyWhQsKhAsrDRew1uoBwUQEjS4YAMKOijPF5m+KGNj5vEyvmnEpRYV7Cv/cT9i3nhH3L4ab4F2grq+DfDknYnUppUeGuJyd9Z/eaEUBhcfT1VKVjH91INVH8GjiUaM3CgOnAcqDMzL7i7k/25KDuPiXB6+uA0zs9fxR4tCf77k9vV9fz0yff4KsnTuGAcX268UtkYErS5NM+/dNs3t5M9bYmauqbqKlrorpuBzV10cfH7jOacw+fSNvTN5LfmvjbfLiogOFDCxlePIThQwuZMHIoI4YWMnl0CQCTI8NoGjaO4u3rPhJeXlkFz115YtJTqBgxNOE3eSuroHhIfmr/Fum4QHck2L4k3nTsoxupJop1wEUd3+hj4xnmEG0+uh/oUaIYuJy3q+vZ1tjafVGRQaSlrZ36Ha2MSND5uv7+aznmnhLa4vThlYQKiIRD7D+uFIC8bYmbbN76/mkU5ie/GbMkVACn3ti3i3S2XOQ79tPXi3o69pFEqoliaudmH3d/3cz2dfd3c6k5Zkp5mKe+fnzQYYj0XoJmo7Z258NtO1i9eTtrtjTu/L1mS/T3+tpG/mXvUdyVoIN0DzbxleP3JhIOEQmHKI/9Hl0SYlho98tMojZ1K6voNkns1NeLdDZd5AeAVBPF62Z2C9FBbwDnxl4LAS39EpmIpFXb4nvJf+TyuM1GZ/x1LMvXbdtZ1gz2KC2iYkQxMyePZMKIYvYbWwp1iS/yV50yLbVA0tWm3teLdI5c5NMh1UTxReBS4IrY8+eBq4gmiXQMuhswLr1rIZNGDeOaU/cNOhTJNSl2Ire2tfP+pgbe3FDPGx/W8VZ1HW98WMcf6r7NOOLfb3/R8Y/T1NpOxYhiJowYytjhRYQK4rTVexY12UjGdJsoYgPtHnX3E4CfxilSn/aostiaLY00NLUFHYbkmniDu2K1AZtxDs+s3MADr67jzQ11vFvTQHNbe3SbwZ4jhzJ1TJixdfHv9KF2DWcdmuKtlGqyyUndJgp3bzOzdjMr6+WtsINKpCTE+lrN9ySZ5Qtu3DU4LMZaGml96kYKZpzDuzUNLFy1haljSjh+WoRpY8JMHRNm70jJrrt4EtzO2eP77XWRzzmpNj3VA6+Z2VPAzkF37j67X6LKYuWlIZasyfl8Kf2strGFRR9sYeH7W6hatZk/bl1LvPtG8uuidxBddMxkvnzsXsl3moH77WVwSjVR3B/7yXmRkhCbG5poa3fy83Lnji/po276F7Y3t/LE8g+pen8LC1dt4Y0NdbhDfp6x/9hStoXKGd684SO7tVhtIKW7D9U3IL2UUqJw99+ZWTEw0d3jzu2UKyLhEO0OmxqaKA8XBR2ODAQJ+heee7OGhmlncdqBY2ltd74+bwnDhhRwyMThnH7gWCr3HMFBE4ZHby9d+r3suFNIclKqkwLOAn4CDAEmm9nBwBx3P6M/g8tGkVhyqN6mRCGp8QVz4vYvTFt+Ez9oPYrTDhxLaVEhT115PJNHD4tfU1VtQAKUatPTDUQn//srgLsvNrNuGkQHp51LotZrXQpJrL3dWbJmK08s38A1tWviTUtEuW/kF+cdsvP5lPKS5DtVbUACkmqiaHH32i7toO39EE/WK9fa2ZJAS1s7L727mSeWf8iTr3/Ihm1NFOQZXy6OMLqt+iPlLY2ze4r0p1QTxXIz+wyQb2b7ALOBf/ZfWNkrEg4xo6KMoalOHCYDX5KO6MbmNsyia6rf+eIqbnz4dYoK8zh+aoRTp+/BidPGUPb293W3kQxoFn/NoC6FzIYC1wEfJzp77BPAf7t7Vg4oqKys9KqqqqDDkMGga0c04IXF2KybeW/cJzjtF3/nf846kE8dUsGGbTtYvHorx+0T+egMpH2dmlukn5nZQnevjLstlUQx0ChRSNrcND3BILUJtF/+Gj98YiVnHDRO087LgJcsUaR619NUonM7Ter8HndPPvH7IPXt+a+xub6ZWz9/WNChSD/zBB3R1K4hL8+49rT9Mh2SSMal2kdxH3ArcBuQ8xMdTRw5lOHFhd0XlAGprd15esUGbn/+PX7SPoqKvI0fLaSOaMkhqSaKVne/pV8jGUAuOX7voEOQfvLPtzfyzfuXsnpzI+OHF/PG9CsY/+b3dy3bCeqIlpyTaqJ42MwuBeYDO+8LdffN/RLVANDRt5NLCzcNVqs2NdDW7uwVKaG8tIgx4SKuPW0/Pr7/GAryT4SlY9QRLTkt1bue3ovzsrt7Vg666+/O7EdfW88V9yzmiSuPY/LoYf12HOl/za3tHPGDpzlq79H832cPDTockcD0uTPb3SenN6SBLVxUQHNbOzV1TUoU2a7LbanNH7ueB1qPZsHKDdzy2cMYUpDHL847hGl7hIOOVCRrJV2g1syu6fT40122/aC/gsp2HXM8aXR2lusYA1G7GnCoXU3bg1/jH/N/zapN23dOw3Lc1AhjSjVvl0gi3a1kfl6nx9d22XZqmmMZMDrme6quy8rxhtJhwZzdR0MDxTTzoxEP8Njlxyo5iKSou0RhCR7He54zhhcXUpBnqlFksfZ2x2vXxN1W1LBeNyGI9EB3icITPI73PGfk5RmRcEiJIot99e5FrG0fFX+jxkCI9Eh3ieIgM9tmZnXAjNjjjucHZiC+rBUJh6hWosgqtY0tNLdGJzX+t4PHs/awq/HC4t0LaQyESI8lTRTunu/upe4edveC2OOO5zk9NDlSohpFNqnetoMTf/JXbn8+eif3xw/YgyPOvASbdTOUTQAs+nvWzRoDIdJDqQ64ky7KS0MsWVMbdBg5b3NDMyOHDaG8tIjzZ07k6Cmjdy+gxX5E+kyJopeO2ns0oYJ83F0dowHYtqOFnz35Jve+sprHLj+WSaOHcdUp04IOS2RQUqLopVkHjWPWQeOCDiPnuDsPLF7L9x9ZyaaGJj57xESGD83pVlCRfqdE0QfNre2YQWF+d/cESK90GVW9tvJqrlw+lZff38xBE4Zz+wWHc2CF1oEQ6W+6wvXS4tVbmXr9Y/zjrThTUEvfxRlVPeLpq9j7w0f537MOZP5XjlKSEMkQJYpemjCimCtPnsqEkUODDmVwijOqeqg1873S+zlv5kTy8tQvJJIpanrqpVElIS4/eZ+gwxi8Eoyqzt+2NsOBiIhqFH1QXbeDD2s131O67Whpo75oj/gbNapaJOMCSRRm9mMzW2lmS81svpkNT1Dut2ZWbWbLMh1jKs6f+yJz/rI86DAGnYeXrOPb2z5FW75GVYtkg6BqFE8B0919BvAmH52ZtsMdZPEstZrvKb12tESXYz/7sAou/Mo3yT9To6pFskEgfRTu/mSnpy8CZyco93czm5SJmHqjPFzEkjVbgw5jUHhw8Vq+/8gK7rn4SPaKlHDwhOEwQaOqRbJBNvRRXAg8FnQQvaEaRd+1tzs/enwll9+zmEmjh1FWrMFzItmm32oUZvY0EK9H8jp3fzBW5jqgFbgrDce7GLgYYOLEiX3dXUoi4RDbm9toaGplWEg3kPVU3Y4Wrrx3MU+vqOb8mRO58YwDGFKQDd9dRKSzfru6ufvJybab2QXAJ4GT3L3Pa1u4+1xgLkBlZWVG1soo37nSXROTlSh6ZNWmBv7j91W8U9PAnDMP4PNH7qk5s0SyVCBXNzM7FbgGON7dtwcRQzp0LIlaU9fE5NHDAo5m4Pjn2xu59I+LAPjDhTM5quuMryKSVYKq5/8KCANPmdliM7sVwMzGmdmjHYXM7G7gBWCama0xs4uCCTe+8nB0zWX1U6TuDy+u4vO/fZlISYgHLztaSUJkAAjqrqcpCV5fB5ze6fn5GQuqFyI7m5406C5V7s7Hpkb4+XkHEy5Sx7XIQKCG9T4YXlzIt0/fl8MnjQw6lOzTaebXttLxvDfjG0w5+UI+f+SefO6IPTVXk8gAokTRB3l5xsXH7R10GNmnY+bX2KR++dvWMP4f36Rp5FBCh56H+qxFBhbdi9hH67Y28uaGuqDDyC5xZn4tppnQ374XUEAi0hdKFH30Xw8s4/J7FgcdRnZJMPNrwtdFJKup6amPLj1hCi1t7UGHkVVaSsZRWB9nOnDN/CoyIKlG0UeH7TmCI/caFXQYWWPNlu38d+OnaSS0+wbN/CoyYClR9FF13Q4eX/YhDU2tQYcSuI31TXzhNy/zQNtRbD7xx5r5VWSQUNNTHy1atZVL7lzIX752DNPH5+4aznU7Wrjg9pdZV9vInRcdwfhJI+G4LwYdloikgWoUfdR5Go9ctaOljf/4fRUr19dxy2cPo1LjSkQGFdUo+qg8xxNFa1s7s+9+lRff3cwvzjuYE/YtDzokEUkz1Sj6KNen8WhobmPDth3cMGt/zjx4fNDhiEg/UI2ij4oK8yktKsjJGkVrWztlxYXcd8lRWkdCZBDTX3caRMIhqnMsUdz6t3f40h2vsKOlTUlCZJDTX3ga5OKSqCOHDSFSEmJIvv4LiQx2anpKg/JwEYtXbw06jIyobWyhrLiQcyon8OnDKrQqnUgO0NfBNOioUaRhRdes9vzbGznmh8/wwjubAJQkRHKEahRp8JkjJnLq9D2CDqNfLVm9lYt/X0XFiKHsNzYcdDgikkFKFGmwd6SEvSNBR9F/3q6u54LbX2bEsCH8/qKZDB86JOiQRCSDlCjSYOv2Zv76Rg0zJ49k3PDioMPpu06r07WGx/OH7WeTn3cMd150BGNKi4KOTkQyTH0UabBhWxNX3LuYRR9sCTqUvutYna52NeAU1K3hW62/Zv6xa5k0eljQ0YlIAJQo0mDy6GE8843jOXm/MUGH0nfxVqezZiYs+klAAYlI0NT0lAZDCvLYK1ISdBjpodXpRKQL1SjSZN4rq3lk6fqgw+i7RKvQaXU6kZylRJEmd760inlVq4MOo8+aP3Z9dDW6zrQ6nUhOU6JIk/JBMN9Tdd0OTniinKoZN2p1OhHZSX0UaRIJh1i8ujboMPpkSH4eMyrKGFZ5Psz6z6DDEZEsoUSRJpGSEJsammhta6dgAE6U197uDB86hFs+d1jQoYhIlhl4V7QsFSktwh02NzQHHUqPPb5sPf9+6z/ZWD+wm85EpH8oUaRJpKRjpbuBdbF9t6aeq+5bSrtDuEgVTBH5KCWKNCkvHXhrZ29vbuUrdy6iMN/49WcPJVSQH3RIIpKF9BUyTXbVKAbG2tnuznXzl/FmdR2/+9JMxg+GOapEpF+oRpEmkfDAqlHc+dIHzH91LVeePJXjpg7iqW9FpM9Uo0iTosJ8Xrj2REYNCwUdSrde/WALcx5ezgnTInz1hClBhyMiWU6JIo3GlmV/882m+iYuvWsRY0qLuOncg8nL0yp1IpKcmp7S6IFX13Lbc+8GHUZSf3zpAzY1NHPLZw/TAkQikhLVKNLomZXVvFNTz5eP3SvoUBK67IQpnLTfGPYfVxp0KCIyQASSKMzsx8AsoBl4B/iSu2/tUmYC8HtgDODAXHf/RaZj7Ymbzj2Y/Cxtynnx3U1UjCimYsRQJQkR6ZGgmp6eAqa7+wzgTeDaOGVagW+4+/7AkcBlZrZ/BmPssfxl98FN0+GG4dHfS+cFHRIALW3tXP2nJXzzz0uDDkVEBqBAahTu/mSnpy8CZ8cpsx5YH3tcZ2YrgPHA6xkJsqeWzqP9wdnktcVWh6tdHV1SFAKfebUwP48/XHgEeZadtR0RyW7Z0Jl9IfBYsgJmNgk4BHgpSZmLzazKzKpqamrSGmBKFszZlSQ6tDRGlxYN0N/frMHdmTR6GBNHDQ00FhEZmPotUZjZ02a2LM7PmZ3KXEe0iemuJPspAf4MXOHu2xKVc/e57l7p7pWRSAADyLJpCdGl8+Cm6fgNw9nrziNZ+Je5mY9BRAaNfmt6cveTk203swuATwInubsnKFNINEnc5e73pz3IdCqriDY3xXs9k5bOizZ5tTRiQEXeRsYv/S7sOSLwJjARGZgCaXoys1OBa4Az3H17gjIG/AZY4e4/y2R8vXLSd7JjCdEFc6JNXp1YFjSBicjAFVQfxa+AMPCUmS02s1sBzGycmT0aK3M08HngxFiZxWZ2ekDxdm/GOTDrZrysAsdY66NZfth/Z/xbvGdTE5iIDApB3fUUd4Ihd18HnB57/A9gYN2mM+McbMY5NDS18p9zX+Cd5xu454CtHDRheEYOX9/Uyva8COXt1R/dmOkmMBEZNLLhrqdBZ1iogN9ecDijw0O48I5XWL05butaWjU2t3HRHa/wg+ZP05pftPvGIJrARGTQUKLoJ+XhIn73pZmcMn2PnYsa9afvPrSMl9/fzAlnX0bBmb+EsgmARX/Pulkd2SLSa5bghqMBrbKy0quqqoIOYzdbtzczpCCPoUP6p7Vv7dZGFq3awqyDxvXL/kVkcDOzhe5eGW+bahQZ0Nzazrn/70W+fu+StO63rd255+UPaGt3xg8vVpIQkX6h2WMzYEhBHhcdO5mJI9M7Mvqp1zfwrftfIxIOcdJ+Y9K6bxGRDmp6CsCK9dvYb2x6ZnB9+b3NzJw8Mi37EpHcpaanLPLsympO+8Vz3PPyB716v7vzi6ffYvm6WgAlCRHpd2p6yrBj9hnN8VMjXPfAMspLQ5y4b8+ajG56+i1uXvAWjS1tHDCurJ+iFBHZRTWKDCvMz+PXnz2U/ceWctldr7Jk9dbu3xTz67++zc0L3uKcygquOWVaP0YpIrKLEkUAug7Ie39jQ7fvuf359/jR429wxkHj+J+zZpCXpSvpicjgo0QRkEg4xO++NJN2d754+8tsrG9KWPbulz/gxodf55QDxvDTcw7K2uVWRWRwUqII0F6REn5zweFs2LaDi+54he3NrR8pM//VNXx7/mt8bFqEm88/hMJ8fWQiklm66gTs0Ikj+OX5h/La2lqu/tPua1o/+tp6vjFvCUdOHsWtnzuMUEF+QFGKSC7TXU9Z4F/3H8OPzj6IqWNKdnv9g83bOWTiCG77YiVFhUoSIhIMJYoscfZhsWnAl86j6YkbCDWs45KyCr58wn9REDoq2OBEJKcpUWSTpfNoffBrhNp2RJ/XrqbgkZT6Gp4AAAeQSURBVMshzzT7q4gERn0U2WTBHAo6kkQHLWMqIgFTosgmWsZURLKQEkU2SbRcqZYxFZEAKVFkk5O+E122tDMtYyoiAVOiyCYzzokuW6plTEUki+iup2wz4xwlBhHJKqpRiIhIUkoUIiKSlBKFiIgkpUQhIiJJKVGIiEhS5u5Bx5B2ZlYDrErjLkcDG9O4vyANlnMZLOcBOpdsNVjOJdXz2NPdI/E2DMpEkW5mVuXulUHHkQ6D5VwGy3mAziVbDZZzScd5qOlJRESSUqIQEZGklChSMzfoANJosJzLYDkP0Llkq8FyLn0+D/VRiIhIUqpRiIhIUkoUIiKSlBJFjJlNMLNnzex1M1tuZpfHKWNmdrOZvW1mS83s0CBiTSbF8/iYmdWa2eLYT1YueGFmRWb2spktiZ3LjXHKhMzs3thn8pKZTcp8pN1L8VwuMLOaTp/Ll4OINRVmlm9mr5rZX+JsGxCfSYduzmUgfSbvm9lrsTir4mzv9fVL04zv0gp8w90XmVkYWGhmT7n7653KnAbsE/s5Argl9jubpHIeAM+5+ycDiK8nmoAT3b3ezAqBf5jZY+7+YqcyFwFb3H2KmZ0H/BA4N4hgu5HKuQDc6+5fDSC+nrocWAGUxtk2UD6TDsnOBQbOZwJwgrsnGlzX6+uXahQx7r7e3RfFHtcR/Y8zvkuxM4Hfe9SLwHAzG5vhUJNK8TwGhNi/c33saWHsp+vdF2cCv4s9/hNwkplZhkJMWYrnMiCYWQXwCeC2BEUGxGcCKZ3LYNLr65cSRRyxqvIhwEtdNo0HVnd6voYsvggnOQ+Af4k1gzxmZgdkNLAeiDULLAaqgafcPeFn4u6tQC0wKrNRpiaFcwH491izwJ/MbEKGQ0zVz4FrgPYE2wfMZ0L35wID4zOB6BePJ81soZldHGd7r69fShRdmFkJ8GfgCnffFnQ8vdXNeSwiOq/LQcAvgQcyHV+q3L3N3Q8GKoCZZjY96Jh6K4VzeRiY5O4zgKfY9a08a5jZJ4Fqd18YdCx9leK5ZP1n0skx7n4o0Samy8zsuHTtWImik1jb8Z+Bu9z9/jhF1gKdv1FUxF7LKt2dh7tv62gGcfdHgUIzG53hMHvE3bcCzwKndtm08zMxswKgDNiU2eh6JtG5uPsmd2+KPb0NOCzTsaXgaOAMM3sfuAc40czu7FJmoHwm3Z7LAPlMAHD3tbHf1cB8YGaXIr2+filRxMTaUH8DrHD3nyUo9hDwhdjdA0cCte6+PmNBpiCV8zCzPTrajM1sJtH/B1n3h2xmETMbHntcDPwrsLJLsYeAL8Yenw0841k4ijSVc+nSXnwG0f6lrOLu17p7hbtPAs4j+u/9uS7FBsRnksq5DITPBMDMhsVuXsHMhgEfB5Z1Kdbr65fuetrlaODzwGuxdmSAbwMTAdz9VuBR4HTgbWA78KUA4uxOKudxNvAVM2sFGoHzsvEPGRgL/M7M8okms3nu/hczmwNUuftDRJPiH8zsbWAz0T/4bJTKucw2szOI3rm2GbggsGh7aIB+JnEN0M9kDDA/9v2vAPijuz9uZpdA369fmsJDRESSUtOTiIgkpUQhIiJJKVGIiEhSShQiIpKUEoWIiCSlRCE5x8zaYjNsLjOz+8xsaJKyZ5jZt3pxjL+aWWXs8aMdYygSlL0iWQwiQVOikFzU6O4Hu/t0oBm4JFFBd3/I3f+36+uxEccpcffTY6OxE7kCUKKQrKVEIbnuOWCKmc2y6NoJr5rZ02Y2BnauR/Cr2OM7zOxWM3sJ+FHnnZhZsZndY2YrzGw+UNxp2/tmNjo2evaR2GSMy8zsXDObDYwDnjWzZ2PlbzGzKuuybkVsPzea2SKLrjuwb+z1EjO7PfbaUjP799jrHzezF2Ll74vN/yXSYxqZLTkrVis4DXgc+AdwpLu7RRenuQb4Rpy3VQBHuXtbl9e/Amx39/3MbAbRiRe7OhVY5+6fiB2/zN1rzezr7L6OwHXuvjk2inuBmc1w96WxbRvd/VAzuxS4Cvgy8F9Ep2M4MLbfEbG5u64HTnb3BjP7JvB1YE6P/6Ek5ylRSC4q7jS9yXNEp5yYBtwbm9tnCPBegvfeFydJABwH3Azg7kvNbGmcMq8BPzWzHwJ/cffnEhzjHItOE11AdOqP/YGO/XVM8rgQOCv2+GQ6TZPh7lssOjPq/sDzsWkdhgAvJDieSFJKFJKLGmPTfe9kZr8EfubuD5nZx4AbEry3obcHdfc3Lbr85OnA98xsgbvv9g3fzCYTrSkcHrvg3wEUdSrSMZNpG8n/fo3omhfn9zZekQ7qoxCJKmPXlMtfTFYwgb8DnwGw6DoTM7oWMLNxRJun7gR+DHSsWVwHhGOPS4kmo9pYP8lpKRz7KeCyTscZAbwIHG1mU2KvDTOzqb04LxElCpGYG4D7zGwhkGjN4WRuAUrMbAXRfoB4i+EcCLwca/b6LvC92OtzgcfN7Fl3XwK8SnQK8j8Cz6dw7O8BI2Id5EuI9nfUEJ3p9O5YM9gLwL69OC8RzR4rIiLJqUYhIiJJKVGIiEhSShQiIpKUEoWIiCSlRCEiIkkpUYiISFJKFCIiktT/B4CD9yXdwF8yAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5gUVdbH8e+ZBAMoCJhAESNmRWcVcwAUSSKK4q4BE++ua2DNiHGNiDmti2LGLAgoklQWRVFJEgVBDIAiguQw6bx/VKEDzvQ0THdXz8zv8zz92NV1p+pch5kzdavuuebuiIiIlCUj6gBERCS9KVGIiEhMShQiIhKTEoWIiMSkRCEiIjFlRR1AMjRs2NCbNm0adRgiIpXGhAkTfnX3bUvbVyUTRdOmTRk/fnzUYYiIVBpm9n1Z+zT0JCIiMSlRiIhITEoUIiISkxKFiIjEpEQhIiIxRZIozKyLmU03s2IzyyujTTMzm1zitcLMeqQ6VhGR6i6qx2OnAZ2B/5bVwN1nAQcDmFkmsAAYmJLoRETkd5FcUbj7zDARxKslMNfdy3zOV0SkWvv+M/jk4aQcurLco+gKvBqrgZl1N7PxZjZ+8eLFKQpLRCRi61fCe9fAc21gwnOQvzrhp0ja0JOZjQJ2KGVXL3cftBnHyQE6Aj1jtXP3vkBfgLy8PK3GJCJV35xRMKQHLJ8Ph/8DTrwJcmon/DRJSxTu3ipBhzoFmOjuixJ0PBGRym3NUhh+I3z1KjRsBheNgJ0PS9rpKkOtp7MpZ9hJRKRacIcZg2DoNbD2Nzj22uCVVSOpp43q8djTzGw+cATwnpkNDz9vZGZDS7SrDbQGBkQRp4hI2lj5M7x+Drx5PmzdGLqPDoaakpwkIKIrCncfSCmPurr7QqBtie3VQIMUhiYikl7cYXL/YKipcD20uh2OuAwyU/fruzIMPYmIVE+/fQdDroRvR0OTI6HjY9Bwj5SHoUQhIpJuiovgi77wwb/BMqHdA3DohZARzYwGJQoRkXTyy9cw+HKY/wXs0Ro6PAx1d4o0JCUKEZF0UFQQzKwecx/k1IHOT8MBXcAs6siUKEREIrdwEgy6DBZNg/06wyn3QZ1Sl6+OhBKFiEhUCtbC6Hvg08eg9nbQ9RXYu13UUf2JEoWISBS+Gxvci1g6Fw45H1r/G3LrRR1VqZQoRERSad0KGHUbjO8H2zSF8wbDbsdFHVVMShQiIqkyewS82wNW/hRMmjvhxqQU8Us0JQoRkWRbvQSG3QBT34Bt94YzX4SdSl3cMy0pUYiIJIs7TB8AQ6+DdcvguBvgmKtSUp8pkZQoRESSYcVP8N5VMGsoNGoOpw6G7feLOqotokQhIpJI7jDxRRhxMxTlw0l3BosKpbCIX6JV3shFRNLN0m+DIn7zxkDTY6DDI9Bg96ijqjAlChGRiiougnH/gQ/vhMxsaP9wMDcioiJ+iaZEISJSEYtmwODLYMEE2KsNtHsQ6jaOOqqEUqIQEdkShfnwyYMw5n6ouTWc3g/2Pz0tivglmhKFiMjmWjAhKOL3y4ygwmube6F2w6ijSholChGReOWvgY/ugnFPQp0d4OzXoVmbqKNKOiUKEZF4zBsTFPH77Ts49AJofTvUrBt1VCmhRCEiEsu65TDyFpjwPGyzK5z/Lux6TNRRpZQShYhIWWa9D+/+C1YtgiMvh+NvhJxaUUeVckoUIiKbWv0rvH89THsLttsPuvaHxodGHVVklChERDZwh6lvwfvXwfqVcEIvOKoHZOVEHVmkIpk2aGZdzGy6mRWbWZm1ds3sX2G7aWb2qpnVTGWcIlKNLF8Ar3aFARdD/d3g7x/DcddV+yQBESUKYBrQGRhTVgMzawxcAeS5+/5AJtA1NeGJSLVRXAzjn4UnDg+ebDr5HrhoBGy3T9SRpY1Ihp7cfSaAlT+DMQvINbMCoBawMMmhiUh1smQuDL4Cvv8Edj0uKOJXf9eoo0o7aXuPwt0XmNn9wA/AWmCEu48oq72ZdQe6AzRp0iQ1QYpI5VRUGEya++guyKwBHR+D5udWyfIbiZC0RGFmo4AdStnVy90HxfH12wCnArsCy4A3zewcd3+5tPbu3hfoC5CXl+dbHLiIVEnvTFpAn+Gz2Hr5LB6s+Qz7+Bxo1g7aPQBb7xh1eGktaYnC3VtV8BCtgHnuvhjAzAYARwKlJgoRkbK8M2kBtwyYyEX+NpfmDGZ5cW3+5T04rtkldFKSKFfaDj0RDDm1MLNaBENPLYHx0YYkIpXR0PcH85Y9xl6ZC3i76GjuKDiXZWzFFyNm0+mQnaIOL+1F9XjsaWY2HzgCeM/MhoefNzKzoQDu/jnwFjARmBrG2jeKeEWkkspfDcN68tT6ntS2dXTLv46rCy5lGVsBsHDZ2ogDrByieuppIDCwlM8XAm1LbN8K3JrC0ESkqvh2dPBE07LvGZjVhltXn8EqNi6/0ahebjSxVTLpPPQkIrL51i6DETfBpJeg/u7QbSiZvzWlaMBUKCj6vVludibXntwswkArDyUKEak6vn4P3r0KVi8OSm8cfwNk59KpabC7z/BZLFy2lkb1crn25GZ0al61lixNFiUKEan8Vv0S1GeaPhC2PwD++ho0ar5Rk07NGysxbCElChGpvNxhyusw7IbgxvWJN8NRV0JmdtSRVSlKFCJSOS37MVgrYs5I2OkwOPVx2Fb3HJJBiUJEKpfiYhjfD0bdFlxRnHIf/OViyMiMOrIqS4lCRCqPX+cE61b/8CnsdkJQxG+bXaKOqspTohCR9FdUCJ89Bh/dA9k14dQn4eC/qohfiihRiEh6+2kKDL4MfvoK9ukAbe+HrUqrNyrJokQhIumpYB2MuQ8+eRhqNYAzX4R9T406qmpJiUJE0s8PnwdXEb/OhoP+CiffBbXqRx1VtaVEISLpY/0q+ODf8EVfqLsTnPM27FHRFQukopQoRCQ9zPkAhvSA5T/CYZdAy1ugxlZRRyUoUYhI1NYsDYr4Te4PDfaEC4dBkxZRRyUlKFGISHRmDIL3roE1S+CYq+HY64LHXyWtKFGISOqtXARDr4GZg2GHA4N7ETseGHVUUgYlChFJHXeY/AoMvxEK1kLLW+HIy1XEL80pUYhIavz2PbzbA+Z+CE2OgI6PQcM9o45K4qBEISLJVVwMXz4No24PSm60vR/yLoKMjKgjkzgpUYhI8iyeHRTx+3FcMB+i/UNQr0nUUclmUqIQkcQrKoCxj8D/ekNObTjtv3DgWSriV0kpUYhIYi2cHJTf+Hkq7NsJ2vaBOttFHZVUgBKFiCRGwdrgCmLso1C7IZz1clDtVSo9JQoRqbjvPwuuIpbMgebnwEl3Qu42UUclCRLJYwdm1sXMpptZsZnlxWh3pZlNC9v2SGWMIhKH9SuDmdXPtYGifDj3HTj1CSWJKiaqK4ppQGfgv2U1MLP9gUuAw4B8YJiZvevuc1IToojE9M3IoIjfigXQ4lI4oRfUqBN1VJIEkSQKd58JYLGfgNgH+Nzd14Rt/0eQXO5LeoAiUrY1S2FYT5jyGjRsBheNgJ0PizoqSaJ0vkcxDbjLzBoAa4G2wPiyGptZd6A7QJMmek5bJOHcYcY7MPRaWPtbUMDv2Gsgq0bUkUmSJS1RmNkooLSFbXu5+6Dyvt7dZ5pZb2AEsBqYDBTFaN8X6AuQl5fnWxS0iJRu5c/w3tXw9buw48Fw7kDY4YCoo5IUSVqicPcKL0vl7v2AfgBmdjcwv6LHFJHN4A6TXobhvaBoPbT+N7T4J2Sm82CEJFpaf7fNbDt3/8XMmhDcn9BqJiKpsnReUMTv29Gwy1HQ4VFouEfUUUkEIkkUZnYa8BiwLfCemU1295PNrBHwjLu3DZu+Hd6jKAD+6e7LoohXpFopLoLP/wsf3gGWCe0ehEMvUBG/aiyqp54GAgNL+XwhwU3rDdvHpDIukWrvl6+DiXPzv4Q9TwqK+NXdKeqoJGJpPfQkIilSmA9jH4YxfSCnDnR+Gg7ooiJ+AihRiMiCiUEp8EXTYP/ToU1vqLNt1FFJGik3UZjZA8Cz7j49BfGISKrkr4HR98Bnj0Od7aHrq7B32/K/TqqdeK4oZgJ9zSwLeA541d2XJzcsEUmq7z4JriKWfguHnB889ppbL+qoJE2V+xiDuz/j7kcB5wFNgSlm9oqZnZDs4EQkwdatgHf/Bc+3Ay+G8wZDx0eVJCSmuO5RmFkmsHf4+hX4CrjKzP7P3bsmMT4RSZTZw4MksfInOOKyoIhfTq2oo5JKIJ57FA8BHYAPgLvd/YtwV28zm5XM4EQkAVYvgWE3wNQ3YNt94MwXYacyq/uL/Ek8VxRTgJvcfXUp+1QyUiRducO0t+H964Ihp+NugGOuhqycqCOTSiaeRPEV0GyTkuDLge91U1skTa1YGBTxmzUUGh0Cpz4O2+8XdVRSScWTKJ4EDiG4sjBgf2A6UNfM/uHuI5IYn4hsDneY+AKMuBmKCuCku6DFPyAjM+rIpBKLp3jLQqC5u+e5+6FAc+BboDVaREgkfSz9Fl7oAEOuhB0Pgks/hSMvU5KQCovnimKvkpPt3H2Gme3t7t+Ws0KdiKRCcRGM+w98eCdkZkOHR4K5Efr5lASJJ1HMMLP/AK+F22eFn9UgqOoqIlFZNCMo4rdgAux1CrR/ELZuFHVUUsXEkyjOBy4FeoTbY4FrCJKEJt2JRKEwHz5+IHjV3BpO7xfUadJVhCRBzEQRTrQb6u4nAA+U0mRVUqISkbLNnwCD/gmLZwYVXtv0htoNoo5KqrCYicLdi8ys2Mzq6lFYkYjlr4GP7oJxT0KdHeDs16FZm6ijkmognqGnVcBUMxsJ/D7pzt2vSFpUIrKxeWOCIn6/fQd5F0Kr24MhJ5EUiCdRDAhfIpJq65YHcyImvgD1d4Nu70HTo6OOSqqZchOFu79gZrlAE3dXbSeRVJn1flDEb9UiOPIKOL6nivhJJMqdcGdmHYDJwLBw+2AzG5zswESqrVWL4a0L4dWukFsfLv4ATrpDSUIiE8/Q020Exf9GA7j7ZDPbLYkxiVRP7jD1TXj/eli/MigDflQPFfGTyMWTKArcffkms7CLkxSPSPW0fD68exV8Mxx2+gt0fAy22yfqqESA+BLFdDP7K5BpZnsCVwCfJjcskWqiuBgmPAcjbwUvgjb3wmHdVZ9J0ko8ieJyoBewHngVGA7ckcygRKqFJXNh8BXw/Sew63FBjab6u0YdlcifxPPU0xqCRNErUSc1sz4Eq+blA3OBC9x9WSnt2gCPAJnAM+5+b6JiEIlMUSGMewI+uhsya0DHx6H5OSq/IWkrnqVQ9yKo7dS0ZHt3P7EC5x0J9HT3QjPrDfQErt/kvJnAEwTlzOcDX5rZYHefUYHzikTr56kw6DL4aTI0awftHoCtd4w6KpGY4hl6ehN4CngGKErESTdZ7GgccEYpzQ4D5rj7twBm9hpwKqBEIZXCTe9Mpf+4H3AghwKuqjGI7hmDyai1DXR5HvbtpKsIqRTiSRSF7v6fJMZwIfB6KZ83Bn4ssT0fOLysg5hZd6A7QJMmTRIZn8hm+9vTnzF27lIADrHZ9M5+mj1tAQMKj6HGsb1pt5+WJZXKI55EMcTMLgUGEtzQBsDdl8b6IjMbBexQyq5e7j4obNMLKAT6xx1xGdy9L9AXIC8vzyt6PJHN9c6kBfQZPosFy9YCkMs6rs16g26Zw/mJ+nTLv47RxQfT+KNFtDtciUIqj3jXowC4tsRnDsScdOfurWLtN7NuQHugpbuX9ot9AbBzie2dws9E0s47kxbQc8BU1hYEo7NHZUzl3qxn2DljMS8Utua+wq6sJheAhWEiEaks4nnqKeHP64VPM10HHBc+VVWaL4E9zWxXggTRFfhromMRSYTbh0xnbUERW7OKXlmvcFbWaL4t3oEu62/hS997o7aN6uVGFKXIlonnqadawFUERQG7h5Pumrn7uxU47+NADWBkOON7nLv/3cwaETwG2zZ8IuoygnkbmcCzJdfuFonapkNNJ2d8yR3Zz1GfFTxZ2JFHCjuzno3Lb2RmGNee3CyKcEW2WDxDT88BE4Ajw+0FBE9CbXGicPc9yvh8IdC2xPZQYOiWnkckWUoONTVkObdlP0/7zM+ZXrwLFxRcy/RSLsRrZGXQ+/QD6dS8ccrjFamIeBLF7u5+lpmdDcEEPDM90yfV08ZXEU7njI+5JfslclnPfQVn0reoPYWb/FgZ8LcWTbiz0wGRxCxSUfEkivxwPQoHMLPdKfH0k0h1UfIqohG/cnd2P47P/IrxxXtxfcElzPU/Xyk0rpfLtSc301WEVGrxJIpbCdai2NnM+gNHAd2SGZRIOuozfBbrCgo4N3MU12e9huHcWnA+Lxa1xjdZ2iU3O5N7Oh+gBCFVQjxPPY00s4lAC4Kr6Cvd/dekRyYSsQ3DTAuXraVRvVxqLJ/L6zlPc1jGLMYUHcCNhRcz37f9vb0RXHbrKkKqmniuKHD3JcB7SY5FJG2UHGbKopCOK1+nR87brCObawr+j7eKjiVIDQElB6nK4koUItVNn+GzWFtQxH72Hb2z+7J/xncMLTqMWwu6sZh6v7fTEJNUB0oUIqX4ddlyrskayN8zh/AbW/H3/B4MKz4MCK4eNgxH6SpCqgMlCqn2Nr0XcU/eakbk9mIXX8CbhcdyZ+E5LKcOECSJsTdUpMK+SOVTZqIws5WEj8Tyx2Csh1+T4+5KMlLplbwXUYt1XLzqBY7+eATLc7bnonU38kHh/r+3zc3O1KxqqZbK/GXv7luV3DazOsA/gf8jqCQrUultuBdxbMZX3J3dj0Ys4YWik+ifcT6XdT6Yr0tcaWiYSaqreGo91QN6AOcBrwB/CZ+CEqn0Vi9bzP3ZL3NG5hjmFDeiS8EtTPBm2HLo1LyxEoMIsYeeGgJXA2cBzwLN3X15qgITSboZg/ig5nXU9RU8VtiJxws7/V7ETxVeRf4Q64rie2AxQVHANcBFJUs8ufuDyQ1NJDE2vVl903H1OeWH+2HmEDLr7kOXpT2ZVPjH0ie6FyGysViJog9/3MzeKkY7kbS18YJCzhErh3HksJcpyiwgs9Vt1Dvics6fsohfdC9CpEyxbmbflsI4RJJiw83qnWwxd2c9w7GZU/miuBkPZV3Oq0efC+hehEh5Yt2juCXG17m735GEeEQS6udlq+mWOYJrs17HMW4quID+RS0hP6P8LxYRIPbQ0+pSPqsNXAQ0AJQoJL0tnsXA3Ds50L9mdNFB3FhwEQtpCAQT50QkPrGGnh7Y8N7MtgKuBC4AXgMeKOvrRCJXVABjH4b/3cfeWblcv/6fvF5wJBvmjepmtcjmiTmPwszqE6yX/TfgBeAQd/8tFYGJbJGFk2HQZbBoKux3Gjmn3McR3xTwiW5Wi2yxWPco+gCdgb7AAe6+KmVRiWyugrUw+l749DGo3RDO6g/7tAegU3OUGEQqINYVxdUES57eBPQqMYfCCG5mb53k2ERKtem8iHsPXckxM2+HpXOh+blw0h2Qu03UYYpUGbHuUeixEEk7JedF1GEN/7fqOY4ZO5LVtRpT+9x3YPcTog5RpMpRBVipVDbMizg+YzJ3ZfdjR5bSr/AUXik6jw+UJESSQolCKpU1y37hgeyXOD3zE2YXN+b0gtuY5Htiy738LxaRLRJJoghvlHcA8oG5wAXuvqyUds8C7YFf3H3/TfdLNeIO0wfyYc1rqeOreaTwNJ4o7EQ+2YCK+IkkU1T3IUYC+7v7gcBsoGcZ7Z4H2qQqKElTK36C1/4Gb10AdXfmjOJ7eKiwy+9JQvMiRJIrkkTh7iPcvTDcHAfsVEa7McDSlAUm6cUdJr4ITxwOcz+A1v9mmyvGcEHn9jSul4sRzLC+p/MBevxVJInS4R7FhcDrFT2ImXUHugM0adKkooeTqC2dB0OugHljYJejoONj0GB3QEX8RFItaYnCzEYBO5Syq5e7Dwrb9AIKgf4VPZ+79yWYHEheXp7ubFZWxUXw+X/hwzvAMqH9Q3BIN8jQ09oiUUlaonD3VrH2m1k3ghvVLd1dv9gFfpkZlN9YMB72PDlIEnV15SAStaieemoDXAcc5+5roohBolVydnWTulk8tevH7DP7KaixFXR+Bg44A0qsqCgi0Ynqev5xglXzRprZZDN7CsDMGpnZ0A2NzOxV4DOgmZnNN7OLoglXEmnD7OoFy9ZygM3lqbVXs8/XjzF/x9Zw2ZdwYBclCZE0EskVhbvvUcbnC4G2JbbPTllQkjJ9hs/CC9bQM+stLs4cymLqcXH+1cxccjRjazeMOjwR2UQ6PPUk1czOKybwcs7T7JqxiFcKT+Sewr+yklrYsrVRhyYipVCikNRZtxxG3sprOc/xXfH2nJ3fi8+K9/t9t2ZXi6QnJQpJjdnDYUgPWPUz3+zejTNnn8hvxX/889PsapH0pYfTJblW/wpvXwyvnAm59eCiUex57iPc2jlPs6tFKgldUUhyuMO0t+H962DdCji+Jxx9FWTlAJpdLVKZKFFI4i1fAO9dBbOHQeNDoePjsP2+UUclIltIiUISp7gYJr4AI2+BogI46S5o8Q/IyIw6MhGpACUKSYwlc2HIlfDdx9D0GOj4KNTfLeqoRCQBlCikYoqLYNyT8OFdkJkNHR6FQ87TzGqRKkSJQrbcoulBEb+FE2GvU6D9g7B1o6ijEpEEU6KQzVe4Hj5+IHjVrAdnPAv7ddZVhEgVpUQh5SpZ6bXV1j/yQI2n2XrlHDjgTGhzL9RuEHWIIpJEShQS04ZKrxSsplfWm1y4fhiL1tdneov/cMQpf406PBFJASUKianP8FkcXDSFe3OeZpeMX3ipsBW9C7tS96sGjD0l6uhEJBWUKKRsa5dx+apH6ZrzEfOKt+es9Tfzue8DwGpVehWpNpQopHRfvwfvXkWXrF94qrADDxWeznpyft+tSq8i1YcShWxs1eKgPtP0AbDdfnx8yKM88pGznqLfm6jSq0j1okQhAXeY8gYMux7yV8MJN8FRV3J8Vg73bPPHU0+N6uVy7cnNVNBPpBpRohBYPh/e/Rd8MwJ2+ktQxG+7vX/frUqvItWbEkV1VlwME56FkbeBFwVzIg7rriJ+IrIRJYrq6tc5MOQK+H4s7HY8dHgEtmkacVAiko6UKKqbokL47HEYfQ9k1giGmZqfo/IbIlImJYrq5OepMOif8NNXsHd7aHs/bL1j1FGJSJpToqgOCtfDmD7wyUOQuw10eQH2PVVXESISl0gShZn1AToA+cBc4AJ3X7ZJm52BF4HtAQf6uvsjqY61sihZuG+jR1h//CIoBf7rLDjobDj5bqhVP+pwRaQSMXdP/UnNTgI+dPdCM+sN4O7Xb9JmR2BHd59oZlsBE4BO7j6jvOPn5eX5+PHjkxF6WtpQuG9twR+T4hpkF/DGnqPY/duXoe5O0P5h2LNVhFGKSDozswnunlfavkiuKNx9RInNccAZpbT5CfgpfL/SzGYCjYFyE0V102f4rI2SxNEZU7nHnmHnbxfDXy6BVrdCja0ijFBEKrOMqAMALgTej9XAzJoCzYHPY7TpbmbjzWz84sWLExpgOntn0gIWhAX6tmYV92X9l5dz7iGfLM5cfwu0u19JQkQqJGlXFGY2CtihlF293H1Q2KYXUAj0j3GcOsDbQA93X1FWO3fvC/SFYOipAqFXGr+vFQGcnPEld2Q/R31W8GRhRx4p7EzDenUjjlBEqoKkJQp3jzkgbmbdgPZASy/jRomZZRMkif7uPiDhQVZyfYbPok7BEu7Pfp52mV8wvXgXLii4lum+qwr3iUjCRPXUUxvgOuA4d19TRhsD+gEz3f3BVMZXKbjTYsVwbq7xErnkc1/BWfQtakdh+C29p/MBqs8kIgkR1T2Kx4GtgJFmNtnMngIws0ZmNjRscxRwLnBi2GaymbWNKN70suwHePl0Hsh5ijnemLb5d/Nk0am/J4nG9XKVJEQkYaJ66mmPMj5fCLQN338CaEZYScXF8OUzMOo2AL46oBfnTd6fNSVG7jTkJCKJppnZlcWv3wQT534cB7ufCO0f5qBtduHu3bRWhIgklxJFuisqgE8fhdG9ITsXOv0nmGEdlt/QWhEikmxKFOnsp6+Cq4ifpwS1mU7pA1ttH3VUIlLNKFGko4J18L/eMPYRqNUAznwJ9u0YdVQiUk0pUaSb7z+DwZfDkm/g4HPg5DuDiq8iIhFRokgX61fCqNvhy6ehbhM4ZwDs0TLqqERElCjSwpxRMKQHLJ8Ph/8dTrwZatSJOioREUCJIlprlsLwG+GrV6HhXnDhMGjSIuqoREQ2okSRYhsWGDpoxWjuzHmeeraajGOugWOvheyaUYcnIvInShQp9M6kBTw4YAw30482OV8ytbgpF/mNnF+/A52UJEQkTSlRpIo70997kiEZz1KTAu4t6MrTRe0oIpNfhs/SpDkRSVtKFKnw23cw5Ep6FY7mc9+bGwouYZ7v+PvuheHCQyIi6UiJIpmKi+CLp+GD28EyuD+rO0+sOhbfpGhvo3q5EQUoIlI+JYpkWTwrKL8x/wvYozW0f4g95mVQc8DUjda3VrVXEUl3ShSJVlQAYx+G/90HObXhtL5w4JlgRqfmQRNVexWRykSJIpEWTgquIhZNg/1OC4r41dl2oyaq9ioilY0SRSIUrIXR98Cnj0Ht7eCs/rBP+6ijEhFJCCWKivpubFDEb+lcaH4unHQn5NaLOioRkYRRothS61YES5KO7wf1doHzBsFux0cclIhI4ilRbInZI+Ddf8GKBdDin3Bir+DGtYhIFaREsTlWL4HhPWHK67Dt3nDRSNj5L1FHJSKSVEoU8XCH6QNg6HWwbhkcdz0cczVk1Yg6MhGRpFOiKM+Kn+C9q2DWUGjUHDoOgh32jzoqEZGUUaIoiztMfBFG3AxF66H1HdDiUsjU/zIRqV4i+a1nZn2ADkA+MBe4wN2XbdKmJjAGqEEQ51vufmtKAlw6D4ZcAfPGwC5HQ8dHocHuKTm1iEi6ySi/SVKMBPZ39wOB2UDPUtqsB05094OAg4E2Zpbc5d+Ki+CzJ+DJI2DBJGj/EJw/RElCRKq1SK4o3H1Eic1xwBmltHFgVbiZHb48aUGt/Q1ePgMWjIc9Tw6SRCBLE10AAAiLSURBVF2V2hARSYcB9wuB10vbYWaZwARgD+AJd/+8rIOYWXegO0CTJk02P4qa9aD+rnD43+GAM8Bs848hIlIFWfCHexIObDYK2KGUXb3cfVDYpheQB3T2GIGYWT1gIHC5u08r79x5eXk+fvz4LQtcRKQaMrMJ7p5X2r6kXVG4e6tY+82sG9AeaBkrSYTHWmZmHwFtgHIThYiIJE4kN7PNrA1wHdDR3deU0Wbb8EoCM8sFWgNfpy5KERGB6J56ehzYChhpZpPN7CkAM2tkZkPDNjsCH5nZFOBLYKS7vxtNuCIi1VdUTz3tUcbnC4G24fspQPNUxiUiIn8W1RWFiIhUEkoUIiISkxKFiIjEpEQhIiIxJW3CXZTMbDHw/RZ+eUPg1wSGE6Wq0peq0g9QX9JRVekHVKwvu7j7tqXtqJKJoiLMbHxZsxMrm6rSl6rSD1Bf0lFV6Qckry8aehIRkZiUKEREJCYlij/rG3UACVRV+lJV+gHqSzqqKv2AJPVF9yhERCQmXVGIiEhMShQiIhJTtU8UZtbHzL42sylmNnBDafNN2tQ0sy/M7Cszm25mt0cRayxx9mNnM/vIzGaE/bgyiljLE09fwnbPmtkvZpa2a5RsRl/amNksM5tjZjekOs54mFmX8N9NsZmV+QimmV1pZtPCtj1SGWM8NqMf/wrbTTOzV82sZirjjEc8fTGzZmGV7g2vFZv9fXH3av0CTgKywve9gd6ltDGgTvg+G/gcaBF17FvQjx2BQ8L3WwGzgX2jjn1L+hLuOxY4BJgWdcwV/L5kAnOB3YAc4Ks0/b7sAzQDRgN5ZbTZn2BxsVoE1alHAXtEHfsW9KMxMA/IDbffALpFHfuW9KWUf2s/E0yui/s81f6Kwt1HuHthuDkO2KmUNu7uq8LN7PCVVk8BxNmPn9x9Yvh+JTCT4AcircTTl7DdGGBpygLbAnH25TBgjrt/6+75wGvAqamKMV7uPtPdZ5XTbB/gc3dfE/b7f0Dn5EcXvzj7AUGiyzWzLILEtzC5kW2+zejLBi2Bue6+WZUrqn2i2MSFwPul7TCzTDObDPxCsIjS5ymNbPOU2Y8NzKwpwXof6dwPiKMvlUhZfWkM/Fhiez5pmMDjNA04xswamFktgvVldo44ps3m7guA+4EfgJ+A5e4+ItqoEqIr8OrmflEkCxelmpmNAnYoZVcvdx8UtukFFAL9SzuGuxcBB4djzAPNbH93T+nYeCL6EbapA7wN9HD3FcmItTyJ6ks6qG59icXdZ5pZb2AEsBqYDBQlNsryVbQfZrYNwVXdrsAy4E0zO8fdX05spOWraF9KHCcH6Aj03NwYqkWicPdWsfabWTegPdDSw4G8GMdaZmYfAW0I/npKmUT0w8yyCZJEf3cfkPAg45TI70nUEtCXBWz8V/dO4WcpV15f4jxGP6AfgJndTXCFlFIJ6EcrYJ67LwYwswHAkUDKE0UiviehU4CJ7r5oc7+w2g89mVkb4Dqgo7uvKaPNthueVjGzXKA18HXqoixfnP0wgh/gme7+YCrj2xzx9KWyiLMvXwJ7mtmu4V99XYHBqYox0cxsu/C/TQjuT7wSbURb5AeghZnVCn9uWhLc06vMzmYLhp0APfUEzCEYH54cvp4KP28EDA3fHwhMAqYQXEXcEnXcW9iPowluwk8p0a5t1LFvSV/C7VcJxo8LCP5qvSjq2CvQl7YET6HNJRhSiDz2UvpyWvj/eT2wCBheRl8+BmYQPL3VMuq4K9CP2wn+IJwGvATUiDr2CvSlNrAEqLsl51EJDxERianaDz2JiEhsShQiIhKTEoWIiMSkRCEiIjEpUYiISExKFCJJYmb1zOzSqOMQqSglCpEyhMXgytyOQz0gpYliC2IUKZf+UUmVZ2bnAdcQTjZ093PDoojPAg2BxcAF7v6DmT0PrCMomDjWzOpvsv0E8ASwLbAGuMTdvzaz7YGnCEqFA/wDuALYPSwmOdLdr90krnPCNjkExRkvdfciM1sFPEJQ9mMtcKq7LzKzbcNzNAkP0cPdx5rZbcDu4bl/MLMrCGZDNwI+I6gkcGh4rqXu/nB4/ruAX9z9kQr9D5aqL+qZhXrplcwXsB/BjOeG4Xb98L9DgPPD9xcC74TvnwfeBTLL2P4A2DN8fzjwYfj+dYJf3BDU/K8LNKWMtTIIynEPAbLD7SeB88L3DnQI398H3BS+fwU4OnzfhKAUC8BtwAT+WDvhcaBn+L5NeLyGYTwTw88zCGaBN4j6e6RX+r90RSFV3YnAm+7+K4C7b1i/4gj+WCfhJYJfyBu86UG14I22w6q7RxJUEt2wr0aJ85wXnqMIWB5WIC1LS4K/8r8Mj5VLUMIeIJ8gOUGQAFqH71sB+5Y499ZhTACD3X1t+P5ogtIOuPswM/stfP+dmS0xs+bA9sAkd18SI0YRQENPIqVZXcZ2BrDM3Q9OwDkMeMHdSyv5XODuG2rrFPHHz2kGwcqK6zY6UJA4No25LM8A3QjKVj+7mTFLNaWb2VLVfQh0MbMGAOE9B4BPCaq0AvyNoJBdTB6s3THPzLqExzIzOyjc/QHBfYkNi1zVBVYSLDlbmg+AM0pUWq1vZruUE8II4PING2ZWVsIaC5wZtjkJKHllM5BgOOovwPByzicCKFFIFefu04G7gP+Z2VfAhvLqlwMXmNkU4FzgyjgP+TfgovBY0/ljydIrgRPMbCrBcNG+4bDOWDObZmZ9NolrBnATMCKMYSTBmuaxXAHkmdkUM5sB/L2MdrcDJ5nZNKALwRrJK8Pz5gMfAW9sMrwmUiZVjxWpYsysBlDk7oVmdgTwnw3DZWaWAUwEurj7N1HGKZWH7lGIVD1NgDfCpJAPXAJgZvsS3CQfqCQhm0NXFCIiEpPuUYiISExKFCIiEpMShYiIxKREISIiMSlRiIhITP8P0PQHdiaii4UAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "w5_JsMXcNVT6"
      },
      "source": [
        "### Remarks\n",
        "\n",
        "Do not be fooled! Real systems are much more difficult to model, requiring more ACSFs, larger NNs, and much larger datasets for training."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oL2OXRhONVT6"
      },
      "source": [
        "## Exercises"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true,
        "id": "nmYIMOxxNVT6"
      },
      "source": [
        "### 1. Create a vaidation set and test the NN performance\n",
        "\n",
        "For simplicity we just checked the error on training data, but it is better to check performance on a validation set not included in the training.\n",
        "\n",
        "Create different dimer configurations and test NN performance on those."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IYN-YwzLNVT6"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "k_CRAqm2NVT7"
      },
      "source": [
        "### 2. Craft you own energy\n",
        "\n",
        "Make the dimer energy expression more complex and attempt to machine-learn it."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bcBRiYVJNVT7"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0vA2YeuDNVT7"
      },
      "source": [
        "### 3. Add/edit the ACSFs parameters\n",
        "\n",
        "Try to change the ACSFs parameters to get better model performance."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k-0vXEPcNVT7"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dU18snMKNVT8"
      },
      "source": [
        "### 4. A real molecule\n",
        "\n",
        "Here is a real organic molecule... try to compute the ACSFs for its atoms using the DScribe package.\n",
        "Documentation can be found here: https://singroup.github.io/dscribe/tutorials/acsf.html"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "scrolled": true,
        "id": "hhFtxKvLNVT8"
      },
      "source": [
        "# atomic positions as matrix\n",
        "molxyz = numpy.load(\"./data/molecule.coords.npy\")\n",
        "# atom types\n",
        "moltyp = numpy.load(\"./data/molecule.types.npy\")\n",
        "\n",
        "atoms_sys = Atoms(positions=molxyz, numbers=moltyp)\n",
        "view(atoms_sys)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Bn57x8MvNVT8"
      },
      "source": [
        "from dscribe.descriptors import ACSF\n",
        "\n",
        "# Setting up the ACSF descriptor\n",
        "acsf = ACSF(\n",
        "    species=[\"H\", \"C\", \"N\", \"O\"],\n",
        "    rcut=6.0,\n",
        "    # configure parameters for desired ACSFs\n",
        "    g2_params=[[1, 1], [1, 2], [1, 3]],\n",
        "    g4_params=[[1, 1, 1], [1, 2, 1], [1, 1, -1], [1, 2, -1]],\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0saXFEmCNVT8"
      },
      "source": [
        "# calculate the descriptor"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}