notebook/egner_polar_inert.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7bbf254d",
   "metadata": {},
   "source": [
    "# Excess properties for a polar-inert mixture\n",
    "\n",
    "In this document we compare excess properties calculated by our `SAC` implementation for a polar-inert system. The pair interaction energies are the same as given by [Egner et al. (1997)](https://doi.org/10.1002/bbpc.19971010208). These energies are set to make the system behave similarly to the acetone/-n-heptane mixture.\n",
    "\n",
    "We start by importing the necessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3f57911f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math, os, sys\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "parent_directory = os.path.dirname(os.path.abspath('.'))\n",
    "sys.path.append(parent_directory)\n",
    "\n",
    "from pysac import SAC, RGAS, RGAS_SI, DEFAULT_Q_EFF"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2c848bfd",
   "metadata": {},
   "source": [
    "## Compound definitions\n",
    "\n",
    "Ther original article considered a simple-cubic lattice (z = 6). Then, in order to compare the results with the cubic lattice system we need to set the areas with 6 segments each.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ef2e3471",
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_eff = DEFAULT_Q_EFF;\n",
    "Q_1   = [3*Q_eff, 3*Q_eff]\n",
    "Q_2   = [6*Q_eff]\n",
    "\n",
    "# The mixture area\n",
    "Q = [Q_1, Q_2]\n",
    "nu = 6"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8b473c56",
   "metadata": {},
   "source": [
    "## Interaction energies\n",
    "\n",
    "It is assumed that the polar compound has two types of segments: interaction between segment types 1 and 2 is attractive; between 1 and 1 or 2 and 2 is repulsive. The interaction energies with the inert compound segments are always zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "879b371a",
   "metadata": {},
   "outputs": [],
   "source": [
    "T0 = 300\n",
    "\n",
    "class EgnerPolarInert(SAC):\n",
    "    def __init__(self):\n",
    "        super().__init__(Q_eff=Q_eff, tol=1e-10)\n",
    "\n",
    "    # Interaction energies according to Table 3 of the original article\n",
    "    def calc_u(self, T, i, j, m, n):\n",
    "        if i == 1 or j == 1:\n",
    "            # no interaction with inert compound\n",
    "            return 0\n",
    "        if (m == n):\n",
    "            # e11 and e22\n",
    "            return 340 * RGAS\n",
    "        \n",
    "        # e12 and e21\n",
    "        return -340 * RGAS\n",
    "\n",
    "# Create an instance of our SAC class\n",
    "sac = EgnerPolarInert()\n",
    "\n",
    "# Set a temperature, we use the reference `T0` so we get the same\n",
    "# results as in the original work.\n",
    "T = T0\n",
    "sac.set_compounds(Q)\n",
    "sac.set_temperature(T)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b1d62130",
   "metadata": {},
   "source": [
    "## Excess properties\n",
    "\n",
    "Here we calculate the excess properties on multiple compositions and compare our results with Monte-Carlo simulations from the original article."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6d008cdf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAELCAYAAAAP/iu7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABRFklEQVR4nO2deZzN1f/Hn2cWs1jHNlPCZN/KNioUERKK+qFC+SoqyVKpiEgpvi3yrVQIlSihBW1UoiLLYGTPkrVhGPsw6/v3x7nDzJg7c2fmc+/n3pnzfDzO4977+ZzP+bw+d/m87znnfd5vJSIYDAaDweBO/OwWYDAYDIbCjzE2BoPBYHA7xtgYDAaDwe0YY2MwGAwGt2OMjcFgMBjcToDdAnyF8uXLS2RkpN0yDAaDwWNER0cfF5EKVrRljI2LREZGsn79ertlGAwGg8dQSu23qi0zjGYwGAwGt2OMjcFgMBjcjjE2BoPBYHA7xtgYDAaDwe0YY2MwGAwGt2OMjcFgMBjcjnF9NhgMPkFiYiLx8fGcPXuW1NRUu+X4NMWKFaN8+fKULl3aY+c0xsZgMHg9iYmJHDhwgLCwMCIjIwkMDEQpZbcsn0REuHDhAocOHSIoKIjg4GCPnNcMoxkMBq8nPj6esLAwypcvT7FixYyhKQBKKUJDQylfvjxxcXEeO68xNgaDwes5e/YspUqVsltGoaJkyZJcvHjRY+czxsZgMHg9qampBAYG2i2jUBEQEEBKSorHzmeMjcFg8AnM0Jm1ePr9NMbGYDAYDG7HGBuDwWAwuB2nrs9Kqbb5bHONiJzP57EGg8FgKITktM7mJ0AAVwb20usJ0AzYUHBpBoPBYABYu3Ytt956K7Vq1bq0LTQ0lD/++MNn5rJyW9T5BLAtD20tLZgcg8FgMGQlOjqazp07M3/+fLul5JvcjE20iKx1pSGllD+u9YIMBoPBkAeio6Np1qyZ3TIKRE7GpjJwzNWGRCRVKVUZOFpgVQaDwWC4RHR0NCtWrGDu3LmXti1evJjKlSvbqCpvODU2InI4r43l5xiDwWAwOCcxMZFt27YRFxfn01EUjOuzwWAw2EhSUhJjxoyhatWqhISE0K5dO/bs2UNQUBCTJk1i8+bNRERE+LShgZxdn/ehvctcQUSkujWSDAaDoWggIvTs2ZPly5czatQomjRpwnfffUf79u1JSkqiSZMmREdHU79+fbulFpic5mxW4LqxMRgMBkMemT59Ot9++y2rV68mKioKgHbt2rFy5Ur++ecfGjduzJw5c1izZg2NGjW6dFzHjh2ZOHGiTarzR05zNv9xxwmVUtcAzwFRQEMgBLhWRP7JUCcS2OekiTAROZWhbjDwMtAHKANsAp4TkZVZzuvnOO+jQASwE3hJRBYW/KoMBouJiICj2fjahIdDbKzn9Xgpw4YNY9OmTbZqaNSoEZMnT87Xsa+99ho9e/a8ZGjSqVGjBqdPn6Z06dJMnz6d6dOnW6DUXuyYs6kB9AROAr/lUncC0DxLOZulzgxgADAG6AL8C/yolGqUpd7LwIvAu8AdwJ/AfKVUp3xeh8HgPrIzNDltN/gce/fuZc+ePXTu3PmKfQcPHqRp06Y2qHIfLmfqVEpdB4wFWgNhQDzwK/CyiPyVh3OuFJFwR5v9gQ451N0rIn/moKkh0At4SERmObatALYCLwF3ObZVBIYDE0XkDcfhy5VSNYCJwHd50G8wGLyE/PYovIEjR44AEBERkWn7iRMniImJoWvXrnbIchsu9WyUUs2ANUAbYAnwOvAt0Bb4UynlsgkWkbR86HTGXUAyMC9D+ynA58DtSqkgx+bbgWLAp1mO/xS4Til1rYWaDAaDIVcqVKgAwM6dOzNtf/nllzl//nyR7dlMALYAt4nIpWEspVRJdAy1CeTcQ8kvE5RSHwDn0Q4Lo7L0ouoD+0QkIctxW9HGpYbjeX0gEdidTT2AejifIzIYDAbLqVmzJvXr12fcuHEUL16cSpUqMW/ePL788kuUUjRp0sRuiZbi6pzNTcCEjIYGwPH6v+i5FCtJBKaiJ/PboIfArgNWKaXqZqhXFj33k5X4DPvTH0+JSFbvuqz1DAaDwSP4+fmxYMECateuzcCBA3nooYcIDw/nvvvuo2HDhoSFhdkt0VJc7dnk5gJtqYu0iPwLPJZh029KqR/QPZFRaM8zuBxpOitZY7S5Wi/zTqUeAR4BqFKlSu7CDYacSEmB/ft1OXAg8/PDh+HcOUhI0MUZSkHFilC8OJQuDVWqQNWquqQ/j4zUdXwkGnBRpk6dOqxYseLS63379hEVFcXLL79soyr34KqxWQM8r5T6KcswWnG0O7HTSXyrEJGDSqnf0SkM0okHsrMCYRn2pz+GKaVUlt5N1npZzzkNmAYQFRVl1hwVJQrqepySAtu3Q3Q0rF+vH2Ni4MKFzPWuukobiAYNoGRJbURCQ3UpXhxCQiA19bIROn/+8vMTJ+Cff2DFCjhz5kqdTZvqEhWlH6++2nUDZFyv3c7q1av5/vvvuemmmwgKCmLDhg1MmjSJhg0b8sgjj9gtz3JcNTbPoz3P9iullqDdiyOAzuh1Mre6Q1w2ZO2hbAXuVkqFZpm3qQckcXmOZisQBFQn87xNPcejq2kUDEWFvLoep6Vpg/LDD7B0qX6eblhKlIDGjeHRR+G663TPo0oVqFwZgoKyby+vnDp1ube0dy9s3KiN3A8/aG2gDcXNN0PHjrpcc43z9ozrtdtJSEjgm2++YdKkSSQnJ1O9enUGDRrEM888Q0CAy47CPoNLVyQia5VSN6HXstyOnuOIB34h767P+UIpVQVoCXyVYfMiYBzQA/jYUS8AuBdYKiKJjno/oI1Pb0f9dPoAW0TEOAcY8s7Ro9qwpBuY48f19qZNtWFJ71nUqgX+/u7VUqaMLtdfn3n7+fOwadPlHtYvv8BCxzrmBg0uG56bb7bO8Blc4rbbbiMmJsZuGR7DZfMpIpuB7lacVCmV3k66b98dSqk4IE5EViil3kQ7L6wG4oDawEggDXg1g6ZNSql5wGSlVCDao2wgcC3asKTXO6aUegsYqZQ6i84kei/adbtwObMb3M+MGTB7th6+Aj0/cscdcPvt0L69fu0tFC8OLVvqAiACW7dqA/nDD/D22/DGG7r39X//Bw88ALfeaqtkQ+FEXemg5YGTKuXspCtE5Fal1ENoo1EDKAkcR/eixolIJqd0pVQI8Ap6cWcZIAYdrubXLPX80QZrAJnD1SxwRXNUVJSsX7/epeszFAJym9uoVQt694YuXaBRI/Dz0QDq58/Dr7/Cl1/CggV67qdSJe2w4Awb7hnbt2+nbt26uVc05Inc3lelVLSIRDmtkAdcNjZKqY7o4arKQHCW3SIira0Q5K0YY1PEyMnYrFkDzZoVPm+vCxdg0SL49FNYssR5PWNsCg2eNDauRhB4Fh3SpQtQHEjNUqyMCmAw2Mvq1c7nL8LD4YYbCp+hAe35du+9sHgxOFa3X0H58p7VZCg0uNr3fwK9yPJqEWkpIm2yFjdqNBjcT0oKzJ8PzZtDixb6xvvcc3DwoP4nn16KitvvsWP6etPS9NxU167awJ45A/36webNdis0+BiuGptSwHwRSXWnGIPB46SkwKxZULMm9OwJcXHwzjvayEycmLN7cFFAKWjVCr7+GnbuhAED4IsvoGFD7RQRHW23QoOP4Kqx+REdssZgKByIaBfg666Dhx7Sw0NffaVvqE88ob2zDJmpWRPefVcb4ldfhbVr9YLRnj31+2Yw5EBehtE6KaVGKqWaKqWqZS3uFGkwWMpPP+l5l+7d9T/3hQv1jbNbN/evhykMlC0LI0fqxaMvvADffQf160P//toQGQzZ4KqxEXTSsleAtcDf2RSDwbvZuBFuu02vhTl2DGbO1HMP99xTOCf83U3p0vDSS9roPPGEXntUsyY8/TScPm23OoOX4aqx+QhoAbyFDpD5UDbFYPBOzpyBoUP1kM/mzTB5MuzapSe6C2FYEI9TseLl97RXL/28dm2YO9cWN2mDd+LqL+1W4AkR+ch9UgwGixGBefPgqae0F9nAgfDKKzqsi8F6qlbVvcVBg/R73bu3jrYwZQrUqWO3OoPNuNqzOQ6YCHwG32HXLujQAe6/X0c7XrNG3/SMoXE/TZvqtUrvvae91a6/HkaPvjLitSFPrF27ltDQUBo1anSptGjRAjuiwOQHV43N28DjSikfjclhKDIkJ8O4cdrLbN067T2VvuLf4Dn8/XXvZudOuO8+3aOsX187ZxjyRXR0NJ07d2bTpk2XyqpVq1A+Mt/oqvEIAxoA25RS7yilXspSxuXWgMHgdnbu1AsyX3xRe5rt2KGHdIyHmX2Eh8Mnn+j4a0FB2jlj6FDTy8nAtm3bUEqxbNmyHOtFR0fTzIf/NLk6ZzMqw/Na2ewXYGzB5RgM+UBED9k884xOOrZggY5gbPAeWreGDRtgxAgdaXrZMh2DrUkTu5XZzoYNGwCIiso5BFl0dDQrVqxg7ty5l7YtXryYypUru1WfVbiaz8YMnxm8kyNH9KLMH3/UeVlmztTZLw3eR0gI/O9/OlL2f/4DN96oXaeffbZI9z6jo6OpXr06YWFhTuskJiaybds24uLiKFWqlAfVWYdTI6KUul8p5ZtXZSgapEcAWLlST/5/950xNL5A+/bw1196fdPzz+tez969l/dHROh1TxnL/v06rXYhJH14bPbs2TRp0oSQkBDq1avH8uXLL9XZvHkzERERPmtoIOc5m/eBY0qpH5VSA5VSV3tKlMGQ7Q1HKb09KUkvIuzeHapV04s1H3/cLMz0JcqWhc8/10NpW7botNnffKP3OUs9nZzsOX0eQkQuTfTPmTOH0aNH88UXX5CWlsaDDz54qV50dDT169e3UWnByWkYrQLQDp3JcjTwrlJqPfA18LWIbHe/PEORxdkN5+hRnUly9Wq9fmbiRAgM9Kg0g0Uopdfi3HyznmPr1k33dIoQu3bt4uzZs7Rv356F6em6gYMHDzJo0CAuXLhASEgI0dHRrFmzhkaNGl2q07FjRyZOnGiD6vzh1NiISDLwvaM8ppRqjjY8DwLjlVK70YbnKxH50wNaDQbN5s16sWbPnnYrMVhB1arw+++6t/rqq7nXz8iwYbBpkztUuU6jRjpqQj6IdkTNfjXLdR8/fpxSpUoREhICwPTp05k+fXpBVNqOyxP/IrJaREaISF2gPjqETWvgD6XUv0qpD9yk0WDIzNq1xtAUNoKD4cMPwcdvqHllw4YNREZGUrt27UzbN27cyPXXX2+TKveQr8BQIrIDmABMUEpdBXQD7rJQl8HgnHr17FZgcBf9++ucOa6Szx6FtxAdHU2TbNy/N27cSNeuXW1Q5D4K7NIsIv+KyPsicocVggwGQxEnPDz77YVsbk5E2LhxI40bN860/eTJk+zfv/+K7b6O056NUuolF9sQIBHYB/woIvFWCDMUccqVgxMnrtzu7EZkKDykp95OTYUnn9SZU0NDtZt7IWLPnj2cPn36ip7Nxo0bAbLt8fgyOQ2jjc5He2eVUreJyPr8CjIY+OUX7eZ61VV67UwGDxxDEcLfXy8CvfZaSEjQwVWrVy80PZx054DsjE1QUBD1CtlwsbIiYqhSKhgdxuYdR5utCtyolxEVFSXr1xsb6nZmz4aHH4ZatbShqVLFbkUGL2B7dDR1AYoV0wnagoPtllQo2L59O3Xr1nW6XykVLSI5x9FxEUvC0IjIRRHZjHYaaGRFm4Yihgi8/DI8+CC0agV//GEMjeEyoaE6IVtqqg6weu6c3YoMecTqmGfrgMEWt2ko7KSkaA+kMWOgb1/doyld2m5VBm+jRAmdhM3fX0f4PnXKbkWGPJBTbLRUpdQNrjaklPIHjgF/WSHMUERITr6c0fGFF2DWLD1UYjBkR3CwNjihobB7N8QbfyRfIScHAQVEKKVcHcsIcBxjMLhGYiLce6+OifXGG/D003YrMvgCgYF6Tm/3bh3AMy0Nype3W5UhF3Jb1PlVHtvzjfykBvtJSNBRf3/8UUdsfvxxuxUZfAl/f6hRA/bsgX/+0QanYkW7VRlyICdj0y+fbe7L53GGosK5c3DnnbBihR4+e+ghuxUZfJGMBufAAe1kYtZheS05BeL82JNCDEWE06fhjjt0fLM5c+D+++1WZPBl/Pz02pt9++DgQd3DMTmNvBKTgdPgOeLj4bbbYP16+OILY2gM1uDnp/MalSsHhw/rYsH6QYO15CsQp8GQZ86e1T2av/6Cr7+GTp3sVmQoTCgFkZH68d9/tQEyPRyvwhgbg/u5cAHuuguio+HLL42hMbgHpXRuHBHdu/HzM3M4XoQxNgb3kpSk0zevWKHnaO4ymSgMbiS9h5Oaqudw/P2NW7SXYOZsDO4jNRX69NERAaZONXM0Bs+glJ7DKVVKu0WbhZ9egUvGRik1QClV3N1iDIWItDQdgmb+fHjzzbwlxDIYCkq6l1qJEtpT7fRpuxW5lbVr1xIaGkqjRo0ulRYtWmBFoGWrcHUY7QPgDaXUp8BUR9BNgyF7RHQeklmzYOxYeOopuxUZiiLp63B27dLRBmrVgpIl7VblFqKjo+ncuTPz58+3W4pTXB1Gqw68B9wDbFRKrVZK9XWkFjAYMjNuHLz9tjY4Y8farcZQlAkI0CkJgoLg77/h/Hm7FV1ix44dKKVyLZ07d861rejoaJo1a+YB1fknT/lslFIBwN3AY8CtwGngY2CaiGx3h0BvweSzcZFZs3REgH79dHQAZcLlGQpObnlXciUpSacmENGBPIOCrBOXT06dOsWOHTsuvV6/fj2DBw/mlVdeoW3btpe2R0REEBkZmWNbjRs35ty5cxQvfnm2Y/HixVSuXDnH4zyZzyZP3mgikgLMB+YrpWoAHwJDgCFKqd+B10TkWyuEGXyQn36CRx6B9u21Q4AxNAZvIT3p2o4duodTp47u9dhImTJluOmmmy69XrduHQDdunXLU5bOxMREtm3bRlxcHKVKlbJcp1Xk2RtNKVVSKfU4sBBoBWwERqEN1yKl1EvWSjT4BFu2wP/9H9Stq50CCknqXkMhIiREOw0kJup4amlpdivKxKZNmwgODqZ27dpX7FuzZg133XUXV111FcHBwVStWpW+ffsCsHnzZiIiIrza0EAejI1SKkopNR04ArwBbAKai0iUiEwUkZbAi8Agdwg1eDFHjuiFmsWLw7ffmsRnBu+lVCm9DufsWdi/36vC2sTExNCgQQP8/f0zbV+3bh233HIL5cqV48MPP+S7777j+eefv2RcoqOjqV+/vh2S84RL/UilVDQ63fM+4CVghohk57y+DBhnmTqD95MewTk+Hn77DXIZIzYYrGTYMNi0Ka9HlYOkEpCYBEFJUKxg8zeNGsHkyQVqgtTUVLZu3Urv3r2v2Pfpp58SGRnJrFmzLm3LOKcTHR3NmjVraNSo0aVtHTt2ZOLEiQUTZTGuDloeAUYDP0jOHgUbgGsLrMrgG6SkwH336V/74sXQuLHdigwG1yhWTA+jJSaB8rN92HfHjh1cvHgxk8FIp2LFiuzevZvhw4fTt29frrvuukz7p0+fzvTp0z2kNP+4amxeBzZkZ2iUUiWAJiKyUkSSgP1WCjR4McOG6WGz99838c4MtpD/HoWCtCD4+4DundesqYfYbGKTo3vWsGHDK/YNHz6cpKQkZs+ezZtvvkmNGjUYOnQoTzzxhIdVFgxX52yWA87cI2o79huKElOn6gybw4fDY4/ZrcZgyDvpUQaCgnR66cREvT0mRqfByFpiYtwmJcbRdnbGJigoiHHjxrF37162bt1Kw4YNGTx4MKtWrXKbHnfgqrHJyYc1CEi1QIvBV/jjDxg8WKcM8LJxYYMhTwQE6CgDIjrKQGoqJCdnX9fZdgvYtGkT1157ba4eZfXq1WPYsGEApKSkuE2PO3A6jKaUigSqZdgU5Rgyy0gI8BBwwHppBq/kyBEdxblKFR3FOYvnjMHgcwQH68Cdf/+tPdRsICYmhpYtW16xfdCgQVy4cIF27dpRqVIl9uzZwyuvvMINN9yQbX1vJqc5m77AWEAc5R0y93DE8TqFPLg7K6WuAZ4DooCGaIN1rYj8k6VeGHquqJujzmrgSRH5K0u9YOBloA9QBu2S/ZyIrMxSz89x3keBCGAn8JKILHRVe5EiIgKOHr1ye2CgnlxdtgzCwjyvy2BwB6VLQ6VKOg+Oh4mNjeXYsWPZDqHVqVOHzz//nEWLFnH+/HmqVq1K7969eeaZZ65wkfZ2cjI2HwG/og3KL2iDsi1LnURglxM3aGfUAHoC0cBvQIesFZRSCliE9mwbDJwERgLLlVKNRORQhuozgM7AM8Beh84flVLNRWRThnovA8PRC1CjgfvQkRC6iMh3edBfNMjO0IAeSpg7Fxo08Kweg8HdRERAQgKcPOnh00Y4jc48ePBgBg8e7FE97sKpsRGR/cB+pVQgOgDnnqy9inyyUkTCAZRS/cnG2AB3ATcDbUVkuaPuavQ6n2fRIXJQSjUEegEPicgsx7YVwFb0eqC7HNsqog3NRBF5w3GO5Y6QOxMBY2zyQvfudiswGKwnPfGah41NUSFXBwERSQbmAZakuxMRV2JE3AUcSTc0juNOA4uBrlnqpetLr5cCfA7crpRKX611O1AM+DTLeT4FrlNKmbVBBoNBz0E6i5lmQjAVCFfX2ewFKrpTSBbqA1uy2b4VeFApVUJEzjnq7RORhGzqFUMP2W111EsEdmdTD7Rb9z6LtBsMBl+mUSOdbO3vv6FsWbj2WhNU1gJcdX1+DRillKrgTjEZKIuep8lK+txQmIv1ymZ4PJXNotSs9TKhlHpEKbVeKbU+Li7OJeEGg6EQkO4wEB8P5rdvCa72bNqib8j7lFJ/Av+ivdHSERHpa6EulaX9jNvdWS8TIjINmAY6n01OdQsd4eHZOwmEh3teSxEkIiKCo9m8/+Hh4cTGxtqgqAgSEaEDdh48qNNLh4barcincbVnczN6biQOnbXzZuCWLMVK4sm+t5HeoznpYr34DI9hDi+3nOoZ0pk3T6+wfuABHUNKRBdzo/MI2RmanLYXBfKS6NESlNJDaAEBOsJAauFau+7p99Olno2IeHoCfSvZe6nVAw445mvS692tlArNMm9TD0ji8hzNVnSkg+pknrdJD8GT1aW7aBMXB7166ZXV771nxqsNtlOsWDEuXLhAqKd7F4GB2uDs2gUHDujnhYQLFy4Q6EGnhzwnT/MQi4BKSqnW6RuUUqWAOx37MtYLBHpkqBcA3AssFRFHsCN+QBufrPG7+wBbRMQ4B6STlgb/+Q+cOKF7NyWyBo0wGDxP+fLlOXToEPHx8SQnJ3v2X3mpUnDVVfo3ceKE587rJkSEhIQEDh8+TMWKnvP7cjkvqlKqOPAwOjtnOeAREflbKXUfsElEduTYQOa20hdqNHU83qGUigPiRGQF2oisBj5VSj3D5UWdCu2sAICIbFJKzQMmO9YD7QMGoheD9s5Q75hS6i1gpFLqLDoVwr3ouaiMrtSGt96C776Dd9/VXjkGgxdQunRpgoKCiIuL48SJE56PCyai52/Wr9eGx8fdoAMDAwkPD/dodk9Xk6dVRkcTuAbYATQASjp2twHaAf3zcN75WV6/53hcAdwqImlKqS7ojKDvAcFo49NGRA5mObYf8AowHh2uJgboKCIbstQbBZwDhnI5XE1PEVmcB92Fm7VrYcQIuPtuePxxu9XYiicm6JOSkjhx4gRxcXHEx8dz8eJFEhMTSUyPPuyExYsXExQURPHixSlfvjwVKlSgTJky+PlZN1DhjQ4KwcHBVLYzOV+pUvoPWOXKsHq1jqlmcBnlSndUKfUF2sB0RCdSSwKiRGSDUqoX8KKI1HKrUpuJioqS9evX2y3DfZw6BU2a6GG0jRuLfNyzK31JLuPKbyYlJYWDBw+yd+9e9uzZc+nx0KFDHD9+nLi4OE6fPm2ZXn9/f8qVK0eFChWoUKECkZGRVK9enWrVql16LF++fI7XlZGCXn+hZckSnZn2iSfgnXfsVuN2lFLRIhJlRVuuDqO1Rw+bHVBKZY3+dhioZIUYg02I6Jw0Bw/q1M5F3NDkBREhNjaWTZs2sXHjRjZt2kRMTAx79+7NNNQTGBjItddeS+XKlYmMjLxkFNJ7JmXLliUkJISgoKBMpVixYqSlpV3q8SQmJpKUlERiYiJnz5691DOKi4u7ZMSOHj3K0qVLOXLkSCatJUuWpG7dujRq1OhSuf766ylevLin3zbfpUsXeOopmDQJ2rWDrmYU3lVcNTbFgLNO9pVGu0UbfJXPP9fOAK+8AjfdZLcar2f58uWsXLmSVatWsWnTJo4dO3ZpX7Vq1WjUqBHdu3fP1KuoVKmSx6P0JiQk8M8//1zqVe3Zs4etW7cyf/58pk2bBugeTK1atWjSpAm33HILrVq1om7duh7V6XNMmADLl8OAAdC8OXhwkt2XcXUYbTXaa2uAo2eTzOVhtPeBWiJym5u12kqhHUY7fFhHcK5bF1audB4XqoiR23CTUorrr7+eJk2a0Lhx40u9hNKlS3tIYf4REQ4ePMimTZsulTVr1lzqCZUrV44TOXhdFelhtHS2boWmTXUCwS+/LLTLA+wYRnsdWOD4Ac51bKunlOqK9lC7ywoxBg8jAv37Q1ISfPyxMTQO9u7dm+P+b7/9lhYtWlCmTBnPCLIYpRRVqlShSpUq3HWX/umKCPv27WPlypWsWLGCjz76yF6R3k79+jB+PDzzDMyeDQ8+aLcir8elng2AUuoxdDj+klwO83IWeMYR1qVQUyh7NlOn6rmaKVOKtPdZWloaq1atYvHixSxZsoRt25yv8S0q4WKceaMBtG/fni5dutC1a1eqVq3qYWVeRGoqtGkDMTGwZYv2UitkWNmzcdnYOE5cHGgBVABOAKtExNlcTqGi0BmbPXugYUNo0QJ+/LHQDgM4Q0SIiYlh7ty5fPbZZxw6dIjAwEBat25Nly5d6Ny5MzVq1LBbpteQkpLCqlWrWLJkCUuWLGH79u0AtGzZkl69etGjRw8qVPBUnF4vYu9euP56Pde5dKkO8VSIsNLYICKmuFCaNm0qhYaUFJGWLUVKlxY5eNBuNR5lz549Mn78eKlbt64AEhAQIF26dJE5c+bI6dOn7ZbnM+zevVsmTJggDRo0EED8/f3ljjvukNmzZ8uZM2fsludZpk7VkQPfecduJZYDrBeL7qGuV9QLJscBS9GxxpYCLwJlrBLjzaVQGZv//ld/9LNn263EI1y8eFE+++wzad26taCjf0urVq3kgw8+kOPHj9stz+fZvHmzjBgxQqpUqSKAhIaGysMPPyzr1q2zW5pnSEsTueMOkZAQkZ077VZjKR43NkBDdMTnJGAlOgLASsfrY8B1Vgny1lJojM3mzSLFioncc4/+kRRi9uzZI88995xUqFBBAKlWrZpMmDBB9u/fb7e0Qklqaqr8/vvv0r9/fwkNDRVAmjZtKh9++KGcO3fObnnu5fBhkbAwkRtvFElOtluNZdhhbFagw9RUzbI90rH9V6sEeWspFMYmMVGkUSORihVFjh2zW41bSE1NlUWLFsntt99+aXinW7du8uOPP0pqaqrd8ooMp06dknfffVfq168vgJQqVUqeeOIJ2bFjh93S3Mfnn+tb6vjxdiuxDDuMTQLQw8m+e4EEqwR5aykUxmb8eP2Rf/213Uos5+LFizJjxgypU6eOAFKpUiV58cUX5WARm5PyNtLS0uS3336T3r17S7FixUQpJXfffbesXr3abmnuoWdPPXKwbZvdSizBDmNzELjLyb6uwEGrBHlr8Xljs3OnSFCQSI8editxifDw8EvzKxlLeHh4pnqnT5+W119/Xa6++moBpGHDhjJnzhxJSkqySbnBGbGxsTJ69GgJCwu7NG+2ZMkSScsynOvqZ++VxMbq4bSbbxYpBD1pO4zN88BvQHCW7SGO7SOsEuStxaeNTWqqSOvWImXKiPz7r91qXCK7m016ERGJi4uTESNGSOnSpQWQtm3byg8//HDFjcvgfZw9e1beeustqVy5sgDSoEEDmT17tqSkpIhI7p+91zNzpr61fvCB3UoKjB3GZjxwwOEM8BHwX8fjMWA/8DLwkqOMs0qcNxWfNjbTp+uPevp0u5W4TE43nBdeeEFKlCghSinp0aOHrF271m65hnyQlJQkn3zyySX36Xr16smCBQt839ikpYm0bStSqpR2HPBhrDQ2rsZGS8u10mVERDwbcdAD+Oyizn//hXr19ALO5ct9ZvFmbrHJunfvzksvvWSCRhYC0tLSWLhwIWPGjGHHjpxzMLpyv/IKdu+G666DTp1g4UK71eQbKxd1urTcVUT88lAKnaHxaYYMgQsXYNo0nzE0ubFhwwbmz59vDE0hwc/Pjx49erBlyxY+/vhju+VYQ40a8OKLOkjnV1/ZrcYryFO4mqKMT/ZsFi3S+TZeeQWef95uNXnCJO8quhSazz45GZo1g2PHYPt28IGI4FnxeM8mw4m7KKVeV0rNUEq9ppTqbIUIgxs4c0YH17zuOh2Z1ofYvHkzgU5yvIeHh3tYjcHTOPuMQ0NDOXPmjIfVFIDAQPjwQzh6VKdbL+K4ZGyUUiWVUiuARcBQoBMwDFiklPpVKVXCfRIN+eL55+HIEZg+XX/pfYC4uDgGDhxI48aNKVmyJFOmTCE5OTnTJGNRiLhc1ImNjc30mR85coR+/fqRkJBAzZo1mTFjBqmpqXbLdI2oKBg6FD74AH7/3W419uKKFwHwDjqdQG/A37HN3/H6DPC2VR4L3lp8yhtt1SoRpUSGDLFbiUskJyfL5MmTpXTp0uLv7y9DhgyREydO2C3L4GWsW7dOWrRoIYA0btxY/vjjD7slucbZsyJVq4rUqSNy8aLdavIENrg+HwGGOtk3FDhslSBvLT5jbFJSRBo2FLnmGhEfiL67YcMGadq0qQDSoUMH2bp1q92SDF5MWlqafPbZZ3LNNdcIIAMHDpRTp07ZLSt3vvtO324nTLBbSZ6ww9gkAu2d7GsPXLRKkLcWnzE2U6boj/WLL+xWkiPnzp2T4cOHi7+/v4SHh8u8efPMgkyDy5w9e1aefPJJ8fPzk6uuukoWLlzo/d+fbt1EQkN9Kq2HHcZmB/A/J/veArZbJchbi08Ym2PHdJSAtm29OqLzDz/8IJGRkQLIgAEDJD4+3m5JBh9l3bp10qhRIwGka9eu3h0Lb+9ekeBgkXvvtVuJy1hpbFz1RpsKDHZ4obVVStVVSrVRSk0FhgAf5GvCyGAto0bBuXPwzjteuabm+PHj9O7dm44dOxIcHMzKlSuZNm0aYWFhdksz+ChRUVGsXbuW1157jaVLl1KvXj2mTJlCWlpe1qF7iGuv1V5p8+bpBdZFDVetEvAqcAFIdZQ0x+tXrLJ83ly8vmezdq12CnjqKbuVZMuSJUskPDxcAgMDZezYsXLRxyZKDd7P3r17pUOHDgLIbbfdJgcOHLBb0pUkJIhERorUry/iA8FisWEYrTQQDIQBd6C90O4AwqwS4u3Fq41NaqrIDTeIhIeLeFlq47Nnz8qAAQMEkOuuu05iYmLslmQoxKSlpcm0adOkePHiUrp0aZk9e7b3zeV8/bW+9b71lt1KcsWjxgYIAFKAO606qS8WrzY2M2boj/KTT+xWkonff/9dqlWrJkopefbZZ01vxuAxdu/efclNunv37t6V/jstTaRjRx2o08ujsFtpbHKdsxGRFOCoY+jM4G2cPKnHgVu2hD597FYDQFJSEiNHjqRVq1akpaWxYsUK/vvf/xIUFGS3NEMRoXr16qxcuZIJEybwzTff0KBBA77//nu7ZWmUgv/9T8csLEKRBVx1EPgU6O9OIYZ8MnYsnDgB777rFU4Be/bsoXnz5kycOJF+/fqxefNmbrnlFrtlGYog/v7+jBgxgrVr11KuXDk6derE0KFDSUxMtFsa1KoFTz8NH38Mq1bZrcYjuJpiYCA6gVos8A3wLzq/xCVEZKY7BHoLXhmIc/NmaNwYHnsMpkyxWw3z58+nf//++Pn5MWvWLLp162a3JIMBgIsXL/Lcc8/x9ttv06xZM+bNm8e1115rr6hz56BOHahYEdatA3/vC5hvZSBO1yZ2tOdZTiXVqnE9by1eN2eTlibSqpVIuXIiNod2uXjxogwaNEgAueGGG2Tfvn226jEYnLFw4UIpXbq0lC5dWr788ku75Yh8/rl4c1ZPbPBGq5pbsUqQtxavMzZffaU/vvfft1XG7t27pUmTJgLIU089JYmJibbqMRhyY+/evRIVFSWADB061N7vbPqfxooVvTK8lMeNjSleZmySkkRq1RKpW1ckOdk2GQsWLJBSpUpJmTJl5Ouvv7ZNh8GQVy5evChDhw4VQJo1ayZ79+61T8zatfpWPHq0fRqcYKWxyWs+mxpKqV5KqWccjzUKNIZnyB/TpsGuXfDaaxAQ4PHTp6amMmLECLp3706dOnXYuHEjXbt29bgOgyG/BAUFMXnyZBYuXMiuXbuIiori559/BiAiIgKl1BUlIiLCPWKaNYP774c334RDh9xzDm/AFYuEXtA5E0gm81xNMvAhEGSV9fPW4jU9m1OnRMqXF2nTxq3xz8LDwwXtBJKpVKhQQe644w4B5NFHHzXDZgafZ/fu3VK/fn3x9/eXyZMnZ/u9Ty9uY98+kWLFRP7zH/edIx9gw5zNu+jIz88D1YDijsdRju0mn42nGDFCf2zr17v1NDn94AICAuQDL53QNBjyw5kzZ6Rbt245fu/damxERIYP1yGnNm1y73nygJXGxlXX5+PAJBF5NZt9o4AnRaR8rg35MF7h+nzgANSuDd27w+zZbj1VTnngV65cadbOGAodaWlpvPzyy7z44otO67hyv8w3J09CjRrQpAksXeoV6+asdH12dc4mCFjrZN8aoJgVYgy5MHo0iMArr9gqwxgaQ2HEz8+PsWPH2icgLAzGjIGffgI/P21sMhZ3zRl5CFeNzU9AByf7OgC/WCPHkImIiMxfttmzITERbrjBbmUGg8EdDBzofN/Ro57T4QZcdWWaBMxWShUH5qNjpYUDPYFOQB+lVLX0yiKy12qhRRJnXy4f/9IZDN5MeHg4R7P5jVWsWNH9Jy9WeAeJXDU2KxyPA4HHMmxXWfan431xFwwu888//+Dv709q6pWxV8PDw21QZDB4jtjY2EvPExMT6devH5999hndunUjJSWFABuWGxQGXH3X+rlVhcFrWL9+PV26dKFkyZJ8/fXXtG7d2m5JBoNtBAUF8emnnxIZGcmECRM4ePAg8+bNo2TJknZL8zlcMjYi8rG7hRjsZ/Hixdx3331UrFiR5cuXU7duXbslGQy24+fnx6uvvkpkZCSPP/44rVu3ZsmSJVx99dV2S/Mp8hRBwFB4ef/99+nWrRv16tXjzz//NIbGZrL6hhQShySf5pFHHmHx4sX8/fff3HTTTWzdutU9J3I2VO3jQ9jG2HgzzrrqFn7pRITx48fz+OOP07lzZ3799VczL+MFGN8Q7+SOO+5g5cqVpKSk0KpVK9atW2f9SWJj9RIHEejbF4KDdRibDHNJvohLizoNNizqTEiA6tWhZk1YscItC7xEhOHDhzNp0iQefPBBZsyYYSY/3Uxqqk7QmJAA58/rx4QEvS0lRZfUVOjUyXkby5bp1CfFikFo6JUlOLjgX5eIiOwNW3i4z9/zLGHv3r20a9eOuLg4Fi1aRJs2bdx1Ir2Qe8AAeO8995wjB6xc1GmMjYt43Ni88QY884w2NK1aWd58SkoKjz76KDNnzmTw4MFMnjwZPz/T0U3HlZttairExel6sbGXy9GjOnnqyZNw6pR+TC8JCe7XrhSULq3XCKaXMmX0Y4UK+trCw/VjeilVKrOByslYmVuG5siRI9x+++38/fffzJs3z33BaAcOhBkzYOdO8HDCtyJhbJRStwLLs9l1WkTKZKgXBrwOdANCgNXo8Dl/ZWkvGHgZ6AOUATYBz4nISlf0eNTYnD2rv1RNm8KPP1refGJiIr169eLLL79k7NixjB07NsfwNEWRnN6Oli1h/344cgTS0q7cX7w4lCuX/c2+VCm9P2NPpHhx3RsJDNQ9Fn9/aN7c+fl/+033gJKSLveSMpZz57SRy2roTp6E48f1sVkJCYEqVaBqVf344YfOz++ltwxbiI+Pp1OnTqxfv56ZM2fy4IMPWn+SQ4d0GJtevWCmZxMiW2ls8j1mopQqC1wLbBERdyb1HgJkHBi99FNR+g65yKFjMHASGAksV0o1EpGM8bpnAJ2BZ4C9wCDgR6VUcxHZ5Eb9eWfyZP3XePx4y5s+f/48d999N8uWLWPy5MkMHTrU8nP4EiK6N7Jli/7j+PffOntDTgQGQtu2ULkyXH115l5CeDiUKOFezTffnP9j09IgPv7K3tiRIzr03v79EBOTcxujR+vR3Tp1oH5991+vN1O2bFl++uknunXrRt++fTl9+jSDBw+29iTXXAOPPw7/+x8895weVvNFXInWCYwGJmR43Qo4A6QCB4CaVkUGzXCOW9GRVtvlUKero06bDNtKA/FkiEQNNHTU65dhWwCwE1jkih6PRX0+cUKkVCmRrl0tb/rkyZPSvHlz8fPzk1mzZlnevrdz7pzI6tUi06aJDB4scuutOqv25dlYkeLFRRo3zrwta/EE4eHZnzs83DPnz+n6/f0vP1dKpFo1/XUdPVpk3jyR7dtFUlM9o9NbuHDhwqWo0S+99JL1Jzh6VCQ0VOS++6xvOwewIcXADmBAhtd/Ar8Bd6J7HZ9bJSjDOVwxNjOAw9ls/xjYn+H1C0ASEJql3jh0ioRc8/F4zNiMHKl/wTExljZ78uRJadasmQQGBnpH7nU3c/GiToD43nsi/fqJXHediJ9fZqNy440i/fuL/O9/Ij//LHLkyOUUQXYbG7vJ6foTE0V27BD58kuRceNEuncXqV078/tbooRI69Y6av7nn4vs2ePW9EteQXJysjz44IMCyNixY60/wciR+s21+N6QE3YYm7PArY7nFRw9mvTX/5fdDb/Awi4bm6OO850A5gJVMtT5E/gxm2OfdRxbwvH6c2BnNvV6OurVz02PR4yNm/69ZDQ0ixYtsrRtd5HXf/axsSILF4o8+aRIs2YigYGXjylfXqRjR5EXXhD56it948vtn3dRNzb56VklJIhER4vMnCny+OMiN9yg84GlH1u2rEinTiKvviqycqXIhQvWndtbSE1NlX79+rnH4MTHi5Qu7ZZRD2fYYWxOAR3ksnE5DwQ6XrcCLlglKMM5GwNvOHpPrYFhwDHgMFDRUWdXdr0qoL/DiFR2vF4K/JlNvXaOerfkpscjxmbYMP33cMcOy5rMaGgWL15sWbvuJreb/b59IjNm6F5LzZqX9wUH63/Uzz4rMn++yD//5O8ftS/f8LyJxERtgKZO1b3IevUuv5fFiom0aKE/q8WLRc6c0cf4uqHPaHDGjBkjaVZ26V56Sb8Ra9ZY12YO2GFs/gA+AUoA3wLfZdjXO+OQlTsL0ATtIDDe8fpv4LNs6g3IYmyWAauzqdc+J2MDPAKsB9ZXqVKlQB9arhw8KBIUpO+eFhEfHy9RUVFSrFgxnzI0IjnfcKpXv/y8XDn9R+/11/V8jMlS7f0cPy6yaJE2Mi1aXO6FBgSItGzp+8ZGRBuchx56yHqDc+aM/tJ36GBNe7lgh7G5HbjoGM66CLTOsG8O8JVVglzQsi196AyduM2VYbR5Xj+MNnCg/tXt22dJc/Hx8dK0aVMpVqyYLFmyxJI2PUlON5w779TzLFu2FP55gKJAQoLIL7+IPP+8HnorDMZGRBuchx9+WAB54YUXrDM4r7+u34zffrOmvRzwuLHR5+RaxxBa9SzbHwVuskqQCzq2Az84ns8EDmVT5yMyOwiMIXsHgRfxBgeBw4f1mMKAAZY0l9HQfPvtt5a06QlOn9ZDXw8+WHhuOIa8k9NnP3KkyB9/iKSk2K3SNVJTU6V///4CyOjRo60xOOfPi1SoIHL77QVvKxdsMTbeUIAoR+9qnON1N0fPpHWGOqXQzgTvZNjWyFGvb4ZtAQ7DtdiVc7vV2Dz1lPYn3bOnwE2dOXNGbrzxRp8xNEePinzwgUi7dnoYBUTCwoyxKcrk9Nmnu12XL69HnL/7zvuHTlNTU2XAgAECyCuvvGJNoxMm6Ddi3Tpr2nOCHcNoXcm8RqUqeqX+WWBB+nCVlcUxPDceuAdoCzwNHEev6ynvqOMHrAIOAvc5hvt+Ra+zqZylvc/Riz77A7c5dF8Emriix23G5tgx7YH2wAMFburChQvSpk0b8ff3l6+//toCce4hNla7JLdpc9ldtmZNkWee0V5Kyclmgr4ok9NnHx8v8tlnIr166eVoIFKmjEjfviJLlniv4UlNTZUHHnhAAHn77bcL3uDp0/rCu3UreFs5YIexWQc8m+H1QscN/k20h9gbVgnKcI6RwGbgNJDsON804Kos9co6htPigQTgZ6BhNu2FoNNbxzqMzBoc7tuuFLcZm+ef1+tqtm0rUDNJSUly5513ilJKPv30U4vEWcfJkyLTp+uFlOkGpnZtvRAwJsbMvRjyzsWL2outb1/tEQz6sW9fkWXLvG+oLTk5+dLCz48++qjgDY4dqy968+aCt+UEO4xNPNBRLt+0LwA9HK/7A3usEuStxS3G5uRJ/fesR48CNZOSkiL333+/APL+++9bo80CkpL0zaBnT+1oByJ16ujfyF9/GQNjsI7ERD2k1q/fZcNTqZL2eNuyxW51l7l48aK0a9dO/Pz8ZOHChQVr7MQJvXrWjVEF7DA2CenzIo4hqGSgtOP1LbhhnY23FbcYm5df1h/Bpk35biItLU0eeeQRAWTixIkWiss/GzeKDBmi5zDTx9cHD9bDy8bAGNzNhQsiX3wh0qXL5XnAxo1F3npLJC7ObnUi586dk+bNm0tgYKD88MMPBWvsuef0UMHOndaIy4IdxmYbMMrx/G0yLJBEe6jFWiXIW4vlxubsWb2k+s47891EWlqaDB8+XAAZOXKkheLyzrlzIh9+qFfvpy/Y69FDr6dISrJVmqEIc/SodpOPirr8vbz/fpFff3X+x8cT84UnT56Uhg0bSkhIiPxWEBfmo0dFQkIsXZ+XETuMzVBHb2Y9elFlxjhpbwC/WCXIW4vlxibdV/7PP/PdxPjx4wWQQYMGWbtKOQ/ExOjQJOmTtfXq6R93fLwtcgwGp/z1l+5hpw+z1akjMmmSHo3KiKc8IY8ePSq1atWSUqVKyYYNG/Lf0JAhugv3zz/WiXPgcWOjz0lv4B3gwSzbp2bdVhiLpcYmIUEkIkL7++aTGTNmCCB9+vSRVDeE2M3p311ysg6u2Ly53hYUJNKnj15jZobJDN7O+fMis2aJ3HTT5e9v376XR7M96XZ/4MABqVKlioSHh8u+/C7oPnhQLwh//HFLtYnYZGyKerHU2Lz7rn7rf/01X4d/99134u/vLx06dJAkN41R5fSDq1xZP9aoof8ZHj/uFgkGg9vZtEkH7wgN1d/pNm08a2xERLZu3SplypSR2rVry/H8/pgeeURbzcOHLdVmV89GAXc5hs1mAVUd21sDV1slyFuLZcYmMVHfrW++OV/dgHXr1knx4sWlcePGciY9cqEbyOkHd+utIt98432upQZDfomPF/nvf0WuucbzxkZEZOXKlRIUFCQtWrSQihUrCnoReqYSntOk0Z49esXrU09ZqstKY+NSWmhH6uXvgBvRSdNKAs1EZINS6lMgXkSG5NqQD2NZWugZM6B/f/jhB7j99hyrRkREcPTo0Su2+/n5cejQIa666qqC63GCyUFvKIokJ0OxYs73u/O7v2DBAnr27ElO9+Qc79d9+8KCBTrdavnylmiyMi20n4v1XgcqAy2B8uheTjo/od2hDbmRmgoTJ0JUFHTokGv17AwNQFpamlsNzerVbmvaYPBqAgN1au/sCA3V6bTdRffu3Zk0aVL+Gxg5Ei5c0GnlvRBXjU1XtOvzanSXLiMH0IbIkBvffAO7d8OIETl3HWxABH75Bdq2hRYt7FZjMNhHbGzmwbPNm+H+++HiRYiMhCFD4OBB95x72LBh+T+4Th245x547z04f94yTVbhqrEpgU5alh3BZO7pGJzxxhtQvTp062a3kkysXAm33AK33QY7dsCbb0LFitnXdfavz2AorFx3Hcydq38bvXrB++/rn/Hjj8O//9qtLgtPPw0nT0KJEvoPbcYSEWGrNFeNzU7A2bhPa+Ava+QUYlat0uNTTz4J/v52qwFgwwbo2BFat4Z9+2DKFNi7F556Sg8XZDdFGhtrt2qDwR5q1tRTrrt3w8MPw/Tp2uiMGAHx8Xarc9C8ufN97hwDdAFXjc0UYJhSahRQxbGtjFKqH/CEY78hKxERl/9VtGyptz3xhO3/MHbsgB49oGlTWLcOXn9d/4AefxyCg22VZjB4PVWr6t7Njh3wf/8Hr70G1arBK6/AuXMFbz/cyfCBs+0+g6tua8BEdPSAVCDN8ZgMvGKVa5w3l3y5PhfAh3L58uXZuj+SmwtkDhw+LPLQQzqUUokSImPGiJw6la+mDAaDg82bdWpyEKlYUeSdd6wL0bRjxw4pXbq0NGjQwPWlDhb6bmPXok50Hpv+wPPoDJ3VrBLi7cWTxmb37t1StmxZqVOnjpw8eTLv581CQoLI+PEixYvr2FBPPqnT6BgMBuv488/Li0Lr1hX5/ntr2l26dKn4+/tLly5dJMWVxW2FwdgU5eIpY3Pq1CmpU6eOlC1bVv7++++8nzMDaWki8+aJVK2qT3nPPZYkAzUYDE5IS9MLnmvU0L+5Tp1EduwoeLtTpkwRQIYPH557ZV82NkA/4EUn+14kQ7rlwlo8YWySk5Pl9ttvl4CAAFm+fHnez5eB6GiRW27Rp2rYUKSAzRkMhjyQmCjyxhs6QG1AgMiwYQUPTjto0CABZMaMGTlXtDBstR3GZhMw2Mm+x4GNVgny1uIJYzNkyBABZPr06Xk/l4P4eB0mSSmdT2baNBNWxmCwi6NHRR59VM+TlisnMnNm/oPVJicnS7t27SQwMFBWrFiR+wGxsTpe2mOP5e+EYo+xOQe0dbKvDXDGKkHeWvJlbEqWzN7QZPMPY+rUqQLIk08+mffziP4Cz5mjJyjTQySZyX+DwTuIidHhEEGkdWuR7dvz1058fLzUrl1bypUrJ3v37s39gP79RYKD8z1Ja6WxcdX1OQUdpiY7KrjYRtEiLQ0qVYImTfTzjOYmy2KVP//8kyeeeIKOHTvy+uuv5/lUu3frMGu9e+sVzuvX64WZpUtbdC0Gg6FAXH89rFgBH36oIxJcfz2MGaOjEuSFsLAwFi9eTGpqKvfccw8JCQk5H/DUU/ok77+ff/FW4YpFApbiJEEa8Avws1XWz1tLnns2ixdr0zJnTo7VYmNj5eqrr5Zq1arJiaxZnHIhMVF7mQUF6bHhd981Q2YGg7dz9KjO/wTakWDZsry38e2334pSSh544IHcEyd27qzH1BMS8nwebBhGa41eU7MGGAB0cjyucWy/1SpB3lrybGxat9apBHJwuE9KSpJWrVpJSEiIbErP3OQi69eL1K+vP8EePSxPY2EwGNzMsmWXvdb+85+8D3uPGzdOAHn33Xdzrrh8uT7J1Kl51uhxY6PPSVdgH3pBZ3rZC9xllRhvLnkyNuvW6bf2zTdzrDZs2DAB5NNPP3W56cREvRjT31/k6qtFlixxXZbBYPAuLlwQef55/Xu+5hqRH390/djU1FTp0qWLBAQEyO+//+68YlqaSNOmIrVqieQxq68txubSAVAbaAHUskqEL5Q8GZsHH9RL9HP4qzJnzhwBZOjQoU7rOPNgBJEHHii4K6XBYPAO1qwRqVNH/7ZDQrL/zWfnuXzy5EmpUaOGREREyJEjR5yf4NNPdSM//JAnXXYMo1XKZX9zqwR5a3HZ2Bw7ppfpDxrktEpMTIyEhITILbfckmNaZwvXZhkMBi8nIUHk6afz/rv/66+/JDQ0VFq2bCmJiYnZV0pM1NaqS5c8abLS2LjqjRajlOqWdaPSjAFWuNhO4efDDyEpCQYNynb3yZMnufvuuwkLC+OLL74gMDDQwwINBoM3EhKis5DklQYNGjBz5kz++OMPnn766ewrFSsGjzwC336rQ7vbgKvG5kdgoVLqPaVUMIBSqjLwKzAaGOceeT5GSop2MbztNqhb94rdaWlp9OnTh4MHD7JgwQIibI7+bDAYCgf33nsvTz/9NO+++y6zZ8/OvtKjj4Kfn31u0K52gYD/AGeALcAQIB7tIHCTVd0sby4uDaMtXKj7ul99le3uCRMmCCBTpkzJsZlz50QeftgMoxkMRZH8/u6Tk5OldevWAjlEi+/RQyQsTOT8eRe1WDeMpnR7rqGUagj8CRQDotFRBSzI4OD9REVFyfr163Ou1LYt7Nmju6n+/kRERHA0m4RFFStWzHY7wJYt0LOnzpWR00eTh4/NYDD4EDlljH/oIXj7bShePPv9R44coVKlSk6Pl5UroVUrnfmtf38XtKhoEYnKtaILuDqMhlLqemAuel3NMiAKeFUpVcwKIT7P1q2wfLnOQObIxOnMoBw7duyKbSIwbRo0a6az/i1d6jwFs6/nUDIYDM5x9vsuXhxmzdL3iL+c5Ea++uqrc2785pt1+IJ33/X4P1aXjI1SajB6Aed5oLGIdETntekHrFNKXTlBUdSYMgWCgnS+2DySkAB9+ugh1VatICYG2rXTUW2y60yb1MwGQ+HF2e/+3DlYtgxOnoQbboCPP85H40rpbMExMfDHH5ZrzwlXezaTgbeBFiKyB0BEZgJN0T2dXMaXCjmnT8Mnn8D990N5ZyHksmffPmjRAj77TKeV/f5703MxGAzZc9tt2k60aAH/+Q8MGQLJyXlspFcvKFNG9248iKvG5nYReU5EUjJuFJFdQHPgPcuV+RIffQTnz8PgwXk67OefISoK9u/XHonPP6+dRQwGg8EZFSvCjz/Ck0/CO+9A+/aQzci8c4oX1yMwCxfCkSNu05kVl25tIvJTDvuSReQZ6yT5GGlpegiteXMd4dlBfHx8jodNmgQdOsBVV8G6dXDHHe4WajAYCgsBAfoeMns2rFmj/7RGR+t94U6GRjJtHzgQUlNh6lQPqNXkxUGgsVLqS6XUcaVUilKqiWP7q0qpju6T6OUsWwZ//63HQR2ICP369XNyQAjBwQt4+mno1g1Wr4YaNTyi1GAwFDL69NFTL0pBy5Z6ND82NjaTy/Fzzz0HwP/+97/LB1avDp06aWOTlOQZsa74RwM3AxeBbei5mzSgiWPfeOBrq3yxvbU4XWfTpYsOA5EhTMTbb78tgJQsOdmpv/z48fnP2GcwGAwZOXZM5NZb9b0lNDTr/SZJoIUoVTJzwrXvv9cV5s512i42hKuZiI4iUB94Ksu+DUCTK44oCuzdqydbHnlEh4MAtm7dyjPPPEPnzp05e3aI00NHjcrZn95gMBhcpUIFvVxi0CDt3ZqZQGAuIn706dOHlBTH1HuHDlCzpsccBVw1Nk2A9x2WLqtz9nGKarbO99/Xa2oeewyAxMREevXqRalSpZgxYwZgrInBYPAMgYHaYSB7qgLvsWrVKiZMmKA3+flp67RqFWzY4HZ9rhqbi0Cok31XAaetkeNDJCbqFVZdu4JjIdWoUaPYvHkzM2fOdDpJZzAYDO4i59GSXvTq1Ytx48axZs0avalvXx0B1AOOAq4am9+BYUop/wzb0ns4D6NTQxctvvoKTpzQQ2jAzz//zJtvvsnAgQPp0qWLCSdjMBi8jilTplCpUiX69OnDuXPn9Hqbnj1h7ly9atSNuGpsXkAPpcU4ngvQVym1HLiJohj1efp0iIyEdu2Ij4+nb9++1KlThzfeeIPkZB3DyGAwGLyJMmXKMHv2bPbs2cOwYcP0xgEDtKGZN8+t53Z1nU0M0Ao4CoxCT0ak+/q2FpGd7pHnpezZA7/8Ag8/jCjFo48+yrFjx5gzZw5+fqH83//pdZ7OguWZETaDweAucrq/zJgBrVq1YsSIEcyYMYOvvvpKhyOoW1f/gXYneXVfA4KBq4FQq1zifKFkcn0eMULEz0/k0CGZNWuWADJx4kQ5e1bktttElBJ5/32n3oQGg8HgUc6fF+nYUXs6T54skpiYKE2bNpWyZcvK4cOHRSZN0js3b850HJ5OC21KBmOTlKTX1dx5p+zevVtKlCghrVq1kuPHU6R5c22DPvnEySduMBgMNnHxosg998ildX7btm2XkJAQad++vaQeParT2Q8enOkYK42NicSVV5YsgaNHSX34Yfr27Yu/vz+TJ8+hfXt/1q+H+fPhgQfsFmkwGAyZCQrS0zIPPACjR8PHH9fhzTcnsWzZMqbMmwf/9386/s2FC245vzE2eWX6dKhUibd37eKPP/5g3Ljp9O59Ddu3w6JFcM89dgs0GAyG7AkI0PPJjz0G//0vbNnyKB07dmLEiBEc7tQJTp3SATrdQJ4ydRZloqKiZP2XX0JkJMcHDqTyzJncfHMf9u6dxrFjiiVLoHVru1UaDAZD7ojAs8/CG29A9+7nWbq0Cg2vr8eK2FjUVVfBypWATZk6DcCMGQgwYPVqgoJqsH37+8THK376yRgag8HgOygFr70G48bBggXFadBgDb/9/gd/1KkDv/2m89JbjDE2eWHmTP6pWZOvN8YSHLyKc+cC+OknuPFGu4UZDAZD3lAKxoyBl1+GVatqULnyd/Re9hMSEAAffmj5+YqUsVFKVVZKLVBKnVZKnXGkTKji0sGnT8OhQzy55wwlS67h3LkSfP89NG3qZtEGg8HgRkaNgpEj4eDBjsTKW/xaqjTy8cc6JJeFFBljo5QKRYfVqQP0BR4AagLLlVJOll9mIC6OPf4V+JalpKRcw7ffKpo3d6tkg8FgcDtK6ZT0w4ZBUtJjPBw/HI4fh2++sfQ8AZa25t0MAKoBtUVkN4BSajPwN/AoMCmng1NPn+NWVkFAPb7+Wpk5GoPBUGhQSmf+vHBBmDp1BMNJ4KWMydYsoCgZm7uAP9MNDYCI7FNK/QF0JRdjs5saJKgmfLXQjw4d3KzUYDAYPIxS8N57ipMnE5j0xUtErLoArLKs/SIzjIZO/LYlm+1bgXq5HXyOEkyddp677jI5agwGQ+HEzw/mzg2lZZNdPMvrlrZdlHo2ZYGT2WyPB8KyO0Ap9QjwiONl4oABZbYMGOAmdd5NeXSSvKKKuf6ie/1F+doBalvVUFEyNnBlllHIIZ2miEwDpgEopdZbtbjJ1yjK1w7m+ovy9Rflawd9/Va1VZSG0U6iezdZCSP7Ho/BYDAYLKIoGZut6HmbrNQDtnlYi8FgMBQpipKxWQTcpJSqlr5BKRUJtHTsy41pbtLlCxTlawdz/UX5+ovytYOF119kAnE6Fm7GABeA0ej5m5eBksD1IuLeBNwGg8FQhCkyPRsROQ+0BXYBs4E5wD6grTE0BoPB4F6KTM/GYDAYDPZRZHo22VGQwJxKqWCl1OtKqX+VUheUUquVUq3crdlK8nv9SqkopdQ0pdQOpVSCUuqAUmqOUupaT+i2ggIFZc3czkillCilfneHTndR0OtXStVVSs1XSh13fP93KqWGulOzlRTwt19FKfWx43ufoJTapZQa71KMRS9AKXWNUuodxz0rwfH9jXTx2Hzf94qssSlwYE6YgY63NgboAvwL/KiUauQWwRZTwOu/D+3Z9zZwBzACaAKsV0pVdptoi7Dgs09vpxowCjjmDp3uoqDXr5SKAtYAQUB/oBPwJuDvLs1WUpDrd+z/CWgFvAB0Bj4EngZmulG2ldQAeqKXfPyWx2Pzf98TkSJZgKFAKlAjw7ZrgRTgqVyObYh2MOiXYVsAsBNYZPe1eeD6K2SzrSqQBrxk97W589qztPMjMBX4Ffjd7uvy0Gfvh15G8JXd12HT9Xdw/PY7ZNk+0XF8qN3X58L1+2V43t9xPZEuHFeg+16R7dngJDAnkB6YM7djk4F5GY5NAT4HbldKBVkv13Lyff0iEpfNtv1AHFDJYp3uoCCfPQBKqV7o3txItyh0LwW5/lvRa9NyDFzr5RTk+os5Hs9k2X4KbYi9PniiiKTl89AC3feKsrEpSGDO+sA+EUnI5thi6G6qt1OgwKRZUUrVBSoC2wuoyxMU6NqVUmHAW8CzIhJvsTZPUJDrv9nxGKyU+lMplayUOqaUelspFWKpSvdRkOv/CZ2W5L9KqXpKqRJKqbbo3tIHor1eCysFuu8VZWOT58CcLh6bvt/bKcj1Z0IpFQB8gO7ZzCi4NLdT0Gt/He1C/5GFmjxJQa7/asfjPGAp0B54DT0cM9cqgW4m39cvIhfRBjd9OPEs8DOwBHjCWpleR4Hue0UtEGdW8hSYM0ud/B7rTVh1De8CLYDOIuIrcebyde1KqVuAB4Em4hi09lHy+9mn/0H9VETGOJ7/qpTyByYqpeqJiC+Ef8rv5x+MNrQV0Y4FB4Ab0BPmKcBACzV6GwW67xVlY1OQwJzxQHZukmEZ9ns7lgQmVUpNQKdh6CsiSy3S5m4Kcu1T0b23Q0qpMo5tAYC/4/UFEbE2ebv1FOT6Tzgel2XZvhQ9Sd4I7481WJDrfxg9b1VDRPY4tq1USp0GpimlPhCRGMuUehcFuu8V5WG0ggTm3Apc63ChzHpsErD7ykO8jgIHJlVKjUK7PQ8VkdkWanM3Bbn2usBj6JtSemkJ3OR47gv/bAv63Ycr/+Gm/7vN7+SzJynI9V8HnMxgaNJZ63isW0Bt3kyB7ntF2dgUJDDnIiAQ6JHh2ADgXmCpD/yzhQIGJlVKDQHGA6NE5B13iXQTBbn2NtmUGPSEcxtggRv0Wk1Brv97IBHomGX77Y5Hy/KfuJGCXH8sEKaUyjoZfqPj8bBVIr2Qgt337Pb5ttHXvDjaEv+Fdne8C33T2AuUyFCvKnosdkyW4z9H/5PtD9yGvslcRI/l23597rx+9KLONPSN56YspZ7d1+buzz6b9n7Ft9bZFPS7P9ax/VWgHbp3ewH4yO5rc/f1A5Fot+dd6AWhbYBnHNvWk2ENizcXoLujvI/upQ50vG6dy2ef7/ue7Rdt8xteBVjo+KKcBb4my+Imx5dLgBezbA9BrzWIdbzZa4Bb7b4mT1w/2gtLnJRf7b4ud3/22bTlU8amoNePHjJ7ynHDTgL2Ay8BgXZfl4euvx7wBXAQbWR3AW8AYXZfVx6uP8ffrzvueyYQp8FgMBjcTlGeszEYDAaDhzDGxmAwGAxuxxgbg8FgMLgdY2wMBoPB4HaMsTEYDAaD2zHGxmAwGAxuxxgbg8FgMLgdY2wMBoPB4HaMsTEYvAClVA1HIrJxWba/r5Q6q5SKskubwWAFxtgYDF6A6BTFHwJPKqXKAyilxgAPAXeLiC8EuDQYnGLC1RgMXoJSKgLYA7wH7ACmAfeLyBe2CjMYLMD0bAwGL0FEYoHJwGB0krYhGQ2NUup5pdROpVSaUqqbPSoNhvxhjI3B4F38DQQBq0VkSpZ9PwOdgJUeV2UwFBBjbAwGL0Ep1Rbdo1kNtFRKNcy4X0TWyJUZIg0Gn8AYG4PBC1BKNUHnVPkQneP+ADo5mcFQKDDGxmCwGUeK4e+BpcBgEUkCxgGdlFKtbBVnMFiEMTYGg404PNCWAtuB3iKS5tj1CdojbaJd2gwGKwmwW4DBUJRxeKBVy2Z7KlDX84oMBvdg1tkYDD6CUmo08BhQATiLzgEf5TBYBoNXY4yNwWAwGNyOmbMxGAwGg9sxxsZgMBgMbscYG4PBYDC4HWNsDAaDweB2jLExGAwGg9sxxsZgMBgMbscYG4PBYDC4HWNsDAaDweB2/h97L+Ve2PUHYQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 51\n",
    "ue = []\n",
    "ae = []\n",
    "tse = []\n",
    "\n",
    "# Monte-Carlo data from Fig. 6 (image recognition) of the original work\n",
    "x1_exp = [0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95]\n",
    "ge_exp = [189, 362.4, 625.9, 819.7, 923.6, 965.4, 930.6, 819.8, 639.8, 369.4, 195.9]\n",
    "he_exp = [320.7, 598.1, 1021.1, 1298, 1462, 1519.7, 1462, 1307.9, 1038.5, 606, 336.6]\n",
    "tse_exp = [133.78, 230.73, 390.3, 487.2, 528.5, 556.3, 535.4, 487.2, 403.7, 230.7, 133.8]\n",
    "\n",
    "x1 = []\n",
    "\n",
    "for l in range(N):\n",
    "    x1_ = l/(N-1)\n",
    "    x1.append(x1_)\n",
    "\n",
    "    x = [x1_, 1-x1_]\n",
    "    \n",
    "    sac.set_composition(x)\n",
    "    sac.calc_ln_gamma()\n",
    "\n",
    "    ae_ = sac.get_helmholtz()\n",
    "    #print(ae_ - sac.get_helmholtz2())\n",
    "    ue_ = sac.get_energy()\n",
    "    se_ = sac.get_entropy()\n",
    "\n",
    "    # subtract the pure compound energies to get the excess\n",
    "    for i in range(len(x)):\n",
    "        ue_ -= x[i]*sac.get_energy_pure(i)\n",
    "        ae_ -= x[i]*sac.get_helmholtz_pure(i)\n",
    "        se_ -= x[i]*sac.get_entropy_pure(i)\n",
    "    \n",
    "    ae.append(ae_ * RGAS_SI * T)\n",
    "    ue.append(ue_ * RGAS_SI * T)\n",
    "    #tse.append((ue_ - ae_) * RGAS_SI * T)\n",
    "    tse.append(se_ * RGAS_SI * T)\n",
    "\n",
    "    #print(se_ - (ue_ - ae_))\n",
    "\n",
    "plt.plot(x1_exp, ge_exp, 'ks')\n",
    "plt.plot(x1_exp, he_exp, 'rs')\n",
    "plt.plot(x1_exp, tse_exp, 'bs')\n",
    "plt.plot(x1, ae, 'k-', label=r'$g^E$')\n",
    "plt.plot(x1, ue, 'r-', label=r'$h^E$')\n",
    "plt.plot(x1, tse, 'b-', label=r'$Ts^E$')\n",
    "\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel(r'$x_1$')\n",
    "plt.ylabel(r'excess property [J/mol]')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1750])\n",
    "\n",
    "plt.rcParams.update({'font.size': 16})\n",
    "plt.tight_layout()\n",
    "plt.plot()\n",
    "plt.savefig('pdf/egner_polar_inert.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b79db323",
   "metadata": {},
   "source": [
    "## Explaining the positive excess entropy\n",
    "\n",
    "Preferential contacts generate negative contributions to the residual entropy.\n",
    "\n",
    "For the present case, the polar compound has a higher order (nonrandomness) when pure than when in mixture. Compound 2 has no preferential interactions when pure, so zero residual entropy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d188ce81",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5 mol of pure compound 1, T*s [J]:  -1427.8685181170763\n",
      "0.5 mol of pure compound 2, T*s [J]:  0.0\n",
      "Equimolar mixture, T*s [J/mol]:  -876.8056302184567\n",
      "Equimolar mixture excess, T*s [J/mol]:  551.0628878986196\n"
     ]
    }
   ],
   "source": [
    "x = [0.5, 0.5]\n",
    "sac.set_composition(x)\n",
    "sac.calc_ln_gamma()\n",
    "\n",
    "print(\"0.5 mol of pure compound 1, T*s [J]: \", x[0]*sac.get_entropy_pure(0) * RGAS_SI*T)\n",
    "print(\"0.5 mol of pure compound 2, T*s [J]: \", x[1]*sac.get_entropy_pure(1) * RGAS_SI*T)\n",
    "print(\"Equimolar mixture, T*s [J/mol]: \", sac.get_entropy() * RGAS_SI*T)\n",
    "print(\"Equimolar mixture excess, T*s [J/mol]: \", (sac.get_entropy()-x[0]*sac.get_entropy_pure(0)-x[1]*sac.get_entropy_pure(1)) * RGAS_SI*T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c99fceb6",
   "metadata": {},
   "source": [
    "### Pair contact entropy contributions for pure compound 1\n",
    "\n",
    "We can investigate the *order* using the individual contributions to the residual entropy of every pair contact in mixture.\n",
    "\n",
    "When pure, compound 1 has a preferential (attractive) contact between segments 1-2. This generates a strong negative residual entropy contribution for the pairs (1-2 and 2-1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a10ea49d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entropy terms (times T) for 0.5 mols of pure compound 1 [J]:\n",
      "segments 1-1: 293.8032823224299\n",
      "segments 1-2: -1007.737532539546\n",
      "segments 2-1: -1007.737532539546\n",
      "segments 2-2: 293.8032646395858\n"
     ]
    }
   ],
   "source": [
    "print(\"Entropy terms (times T) for 0.5 mols of pure compound 1 [J]:\")\n",
    "i = j = 0\n",
    "alpha = sac.get_nonrandom_pure(i)\n",
    "theta_i = sac.theta_pure[i]\n",
    "seg_gamma_i = sac.seg_gamma_pure[i]\n",
    "theta_j = sac.theta_pure[j]\n",
    "seg_gamma_j = sac.seg_gamma_pure[j]\n",
    "psi_ij = sac.psi[i][j]\n",
    "for m in range(len(theta_i)):\n",
    "    for n in range(len(theta_j)):\n",
    "        alpha_mn = seg_gamma_i[m] * seg_gamma_j[n] * psi_ij[m][n]\n",
    "        si = -nu/2 * theta_i[m]*theta_j[n] * alpha_mn * (math.log(alpha_mn)) * RGAS_SI*T\n",
    "\n",
    "        print(f\"segments {m+1}-{n+1}: {x[0]*si}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f04e749f",
   "metadata": {},
   "source": [
    "### Pair contact entropy contributions for an equimolar mixture\n",
    "\n",
    "When in an equimolar mixture with an inert compound, the preferential compound 1 contacts (1-2 and 2-1)\n",
    "become less likely and their negative residual entropy contribution is reduced. This, along with other smaller positive contributions for contacts with compound 2, lead to a positive excess entropy.\n",
    "\n",
    "Note the preferential contact between neutral segments of compound 2. This produces a negative contribution to the miture residual entropy, but still not enough to make the excess entropy negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "bf192fd6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Entropy terms (times T) for 0.5 mols of each compound [J]:\n",
      "Compounds: 1-1\n",
      "segments 1-1: 154.05355397190132\n",
      "segments 1-2: -692.0927456159067\n",
      "segments 2-1: -692.0927456159067\n",
      "segments 2-2: 154.05355483558998\n",
      "Compounds: 1-2\n",
      "segments 1-1: 116.45616731233804\n",
      "segments 2-1: 116.45616180960346\n",
      "Compounds: 2-1\n",
      "segments 1-1: 116.45616731233804\n",
      "segments 1-2: 116.45616180960346\n",
      "Compounds: 2-2\n",
      "segments 1-1: -266.55190603801753\n"
     ]
    }
   ],
   "source": [
    "print(\"\\nEntropy terms (times T) for 0.5 mols of each compound [J]:\")\n",
    "for i in range(sac.ncomps):\n",
    "    theta_i = sac.theta[i]\n",
    "    seg_gamma_i = sac.seg_gamma[i]\n",
    "\n",
    "    for j in range(sac.ncomps):\n",
    "        theta_j = sac.theta[j]\n",
    "        seg_gamma_j = sac.seg_gamma[j]\n",
    "\n",
    "        psi_ij = sac.psi[i][j]\n",
    "\n",
    "        print(f\"Compounds: {i+1}-{j+1}\")\n",
    "\n",
    "        for m in range(len(theta_i)):\n",
    "            for n in range(len(theta_j)):\n",
    "                alpha = seg_gamma_i[m] * seg_gamma_j[n] * psi_ij[m][n]\n",
    "                si = -nu/2 * theta_i[m]*theta_j[n] * alpha * (math.log(alpha)) * RGAS_SI*T\n",
    "\n",
    "                print(f\"segments {m+1}-{n+1}: {si}\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}