notebook/prob_nearly_ideal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7bbf254d",
   "metadata": {},
   "source": [
    "# Pair probabilities for a nearly ideal symmetric mixture\n",
    "\n",
    "In this document we compare results from our SAC equations with the ones obtained with classic models.\n",
    "\n",
    "We start by importing the necessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c848bfd",
   "metadata": {},
   "source": [
    "### Compound definitions\n",
    "\n",
    "For the symmetric case, we set both compounds as having an unitary area such that area fractions match molar fractions. Additionally, in order to compare the results with the cubic lattice system we need to set `Q_eff = 1/z`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ef2e3471",
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_1   = [1]\n",
    "Q_2   = [1]\n",
    "\n",
    "z = 6\n",
    "Q_eff = 1/z\n",
    "\n",
    "# The mixture area\n",
    "Q = [Q_1, Q_2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b473c56",
   "metadata": {},
   "source": [
    "## Interaction energies\n",
    "\n",
    "Here we study the case of a nearly ideal mixture, without much difference for the different pair contacts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "879b371a",
   "metadata": {},
   "outputs": [],
   "source": [
    "class MonomerStronglyNonIdeal(SAC):\n",
    "    def __init__(self):\n",
    "        super().__init__(Q_eff=Q_eff)\n",
    "\n",
    "    def calc_u(self, T, i, j, m, n):\n",
    "        if i == 0 and j == 0:\n",
    "            return -0.5 * RGAS * self.T\n",
    "        if i == 1 and j == 1:\n",
    "            return -0.4 * RGAS * self.T\n",
    "        \n",
    "        return -0.3 * RGAS * self.T\n",
    "\n",
    "# Create an instance of our\n",
    "sac = MonomerStronglyNonIdeal()\n",
    "\n",
    "# Set a temperature, actually the way we define the energies, any temperature\n",
    "# would produce the same results\n",
    "T = 300\n",
    "sac.set_compounds(Q)\n",
    "sac.set_temperature(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09927e6d",
   "metadata": {},
   "source": [
    "Now we compare the `Gamma_2/Gamma1` quotient for different models. For UNIQUAC this quotient is always 1 while for the NRTL model it is a constant different from 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6d008cdf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEMCAYAAADTfFGvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU1f3/8dcHCJuAAnHDoCAoVBYjBgguMEAWAogKFJVSBYpQpQsKVtRWrbi1X/yK/qy16NfiQjHKIgJJSAKyBigBg4CACCoEQQhCFJElyfn9kWRkJxOSmcnk/Xw8eJCTc2fmk/OYzDv3nnvPNeccIiIiJVUl0AWIiEjFouAQERGfKDhERMQnCg4REfGJgkNERHyi4BAREZ/4LTjM7E0z221m686wjcfMssxsvZkt9FdtIiJScuav6zjMrDNwAHjbOdf6FP0XABlAD+fcNjO7yDm32y/FiYhIifltj8M5twj47gybDASmO+e2FW2v0BARCULBNMdxNVDfzBaY2SozuzvQBYmIyMmqBbqAY1QDrge6A7WAZWa23Dn3+YkbmtlwYDjAeeedd33Lli39WqiISEW3atWqHOfchaV5bDAFRzaQ45z7EfjRzBYB1wInBYdzbiIwESAqKsplZmb6tVARkYrOzL4u7WOD6VDVTOBmM6tmZrWBjsCGANckIiIn8Nseh5lNATxAuJllA08AYQDOudeccxvMLAX4FCgA3nDOnfbUXRERCQy/BYdz7q4SbPM/wP/4oRwRESmlYJrjKDNHjx4lOzubQ4cOBbqUoFCzZk0iIiIICwsLdCkiEgJCMjiys7OpW7cuTZo0wcwCXU5AOefYu3cv2dnZNG3aNNDliEgICKbJ8TJz6NAhGjZsWOlDA8DMaNiwofa+RKTMhGRwAAqNY2gsRKQshWxwBINnnnmGVq1a0bZtWyIjI1mxYgUAe/bsISwsjH/961/HbX/gwAFGjBhBs2bNaNWqFZ07d/Y+RkQkWITkHEcwWLZsGbNnz2b16tXUqFGDnJwcjhw5AsAHH3xAdHQ0U6ZMYcSIEd7HDBs2jKZNm7J582aqVKnC1q1b2bBBl7KISHBRcJSTnTt3Eh4eTo0aNQAIDw/39k2ZMoUXXniBgQMHsmPHDi677DK2bNnCihUrmDx5MlWqFO4IXnnllVx55ZUBqV9E5HRCPjhGjRpFVlZWmT5nZGQkEyZMOOM2cXFxPPXUU1x99dXExMRwxx130KVLF7Zv386uXbvo0KEDAwYMIDExkQcffJD169cTGRlJ1apVy7RWEZGypjmOclKnTh1WrVrFxIkTufDCC7njjjuYNGkS7733HgMGDADgzjvvZMqUKQGuVETENyG/x3G2PYPyVLVqVTweDx6PhzZt2vDWW2+xY8cOvv32WyZPngzAN998w+bNm2nVqhVr1qyhoKDAe6hKRCQY6ROqnGzatInNmzd721lZWeTl5fHjjz+yY8cOvvrqK7766iseeeQR3nvvPZo1a0ZUVBRPPPEExXdl3Lx5MzNnzgzUjyAickoKjnJy4MAB7rnnHq655hratm3LZ599RrNmzbj99tuP265fv37ew1VvvPEGu3btonnz5rRp04Z7772XRo0aBaJ8EZHT8ts9x8vLqe7HsWHDBn7xi18EqKLgpDERkWOZ2SrnXFRpHqs9DhER8YmCQ0REfKLgEBERnyg4RETEJwoOERHxiYJDRER8ouAoJ2bG6NGjve3x48fz5JNPAvDkk09y2WWXERkZyTXXXOO9jmPkyJHe79WqVYvIyEgiIyOZOnUqgwcPZurUqYH4UUREjqPgKCc1atRg+vTp5OTknLL/gQceICsri5kzZzJixAiOHj3KP/7xD7KyskhKSqJZs2ZkZWWRlZVF//79/Vy9iMjpKTjKSbVq1Rg+fDgvvvjiGbe76qqrqF27Nvv27fNTZSIi58Zvixya2ZtAb2C3c671Kfo9wEzgy6JvTXfOPVUWr+3xeE763oABA7j//vs5ePAgPXv2PKl/8ODBDB48mJycnJP+4l+wYEGJXnfkyJG0bduWP/3pT6fdZvXq1Vx11VVcdNFFJXpOEZFA8+fquJOAV4C3z7DNYudcb/+UU/7q1avH3Xffzcsvv0ytWrWO63vxxRd5/fXX2bp1KykpKQGqUETEd34LDufcIjNr4q/XO9aZ9hBq1659xv7w8PAS72GcyqhRo2jXrh1Dhgw57vsPPPAAY8aMYfr06dx9991s2bKFmjVrlvp1RET8JdjmODqZ2RozSzazVoEupiw0aNCAAQMG8H//93+n7O/bty9RUVG89dZbfq5MRKR0gik4VgNXOOeuBf4f8OHpNjSz4WaWaWaZe/bs8VuBpTV69OjTnl0F8Pjjj/O///u/FBQUnPF5RowYQUREBBEREXTq1KmsyxQRKRG/LqtedKhq9qkmx0+x7VdAlHPu9J+4aFn1ktKYiMixQmJZdTO7xMys6OsOFNa2N7BViYjIifx5Ou4UwAOEm1k28AQQBuCcew3oD9xnZnnAT8CdrqLfZUpEJAT586yqu87S/wqFp+uKiEgQC5pDVSIiUjEoOERExCcKDhER8YmCoxw88MADTJgwwduOj49n2LBh3vbo0aN59tlnvWtgLViwgN69Q2alFREJcQqOcnDDDTeQkZEBQEFBATk5Oaxfv97bn5GRQffu3XV/DRGpkBQc5eDGG2/0Bsf69etp3bo1devWZd++fRw+fJgNGzZQv359Wrc++TrIhQsXem/gdN111/HDDz/gnOOhhx6idevWtGnThsTERKBwT8Xj8dC/f39atmzJr371K3QGs4iUt0oRHB6Ph0mTJgFw9OhRPB4P7777LgAHDx7E4/F4P4xzc3PxeDxMnz4dgJycHDweD7NmzQJg165dZ329Ro0aUa1aNbZt20ZGRgadOnWiY8eOLFu2jMzMTNq2bUv16tVP+djx48d7b+i0ePFiatWqxfTp08nKymLNmjWkp6fz0EMPsXPnTgA++eQTJkyYwGeffcbWrVtZunTpOY2ViMjZVIrgCITivY7i4OjUqZO3fcMNN5zxcQ8++CAvv/wy+/fvp1q1aixZsoS77rqLqlWrcvHFF9OlSxdWrlwJQIcOHYiIiKBKlSpERkby1Vdf+eknFJHKyp/34wiYY5dFDwsLO6594rLq559//nHtE5dVv+SSS0r0msXzHGvXrqV169Y0btyYF154gXr16jF06NDTPm7s2LH06tWLpKQkoqOjSU9PP+Phpxo1ani/rlq1Knl5eSWqT0SktLTHUU5uvPFGZs+eTYMGDahatSoNGjRg//79LFu27Iwr227ZsoU2bdrw8MMPExUVxcaNG+ncuTOJiYnk5+ezZ88eFi1aRIcOHfz404iI/KxS7HEEQps2bcjJyWHgwIHHfe/AgQOEh4dz4MCBUz5uwoQJfPzxx1StWpVrrrmGhIQEqlevzrJly7j22msxM/7+979zySWXsHHjRn/9OCIiXn5dVr08aFn1ktGYiMixQmJZdRERqRgUHCIi4hMFh4iI+CRkg6Oiz92UJY2FiJSlkAyOmjVrsnfvXn1gUhgae/fupWbNmoEuRURCREiejhsREUF2djZ79uwJdClBoWbNmkRERAS6DBEJESEZHGFhYTRt2jTQZYiIhKSQPFQlIiLlR8EhIiI+UXCIiIhP/BYcZvamme02s3Vn2a69meWbWX9/1SYiIiXnzz2OSUCPM21gZlWBvwFz/VGQiIj4zm/B4ZxbBHx3ls1+D0wDdpd/RSIiUhpBM8dhZpcBtwOvBboWERE5vaAJDmAC8LBzLv9sG5rZcDPLNLNMXeQnIuJfwXQBYBTwnpkBhAM9zSzPOffhiRs65yYCE6Hwfhx+rVJEpJILmuBwznkv9TazScDsU4WGiIgElt+Cw8ymAB4g3MyygSeAMADnnOY1REQqCL8Fh3PuLh+2HVyOpYiIyDkIpslxERGpABQcIiLiEwWHiIj4RMEhIiI+UXCIiIhPFBwiIuITBYeIiPhEwSEiIj5RcIiIiE8UHCIi4hMFh4iI+ETBISIiPlFwiIiITxQcIiLiEwWHiIj4RMEhIiI+UXCIiIhPFBwiIuITBYeIiPhEwSEiIj5RcIiIiE8UHCIi4hO/BYeZvWlmu81s3Wn6bzWzT80sy8wyzewmf9UmIiIl5889jklAjzP0zwOudc5FAkOBN/xRlIiI+MZvweGcWwR8d4b+A845V9Q8D3Cn21ZERAInqOY4zOx2M9sIzKFwr0NERIJMUAWHc26Gc64lcBsw7nTbmdnwonmQzD179vivQBERCa7gKFZ0WKuZmYWfpn+icy7KORd14YUX+rk6EZHKLWiCw8yam5kVfd0OqA7sDWxVIiJyomr+eiEzmwJ4gHAzywaeAMIAnHOvAf2Au83sKPATcMcxk+UiIhIk/BYczrm7ztL/N+BvfipHRERKKWgOVYmISMVwzsFhZg+XRSEiIlIx+HyoyszeP7YJRKJDTCIilUZp5ji+d84NK26Y2T/LsB4REQlypTlU9cwJ7cfKohAREakYzhocZjb82LZz7ssT2qddf0pEREJPSfY4flvuVYiISIVRkuCwcq9CREQqjJIEh67eFhERr5KcVdXWzL7n5z2P4iAxwDnn6pVLZSIiEpRKEhydgeXOuYLyLkZERIJfSQ5V3Q1kmtl7ZjbYzC4p76JERCR4nXWPwzn3WwAzawkkAJPM7HzgYyAFWOqcyy/XKkVEJGiU+AJA59xG59yLzrkeQDdgCfBLYEV5FSciIsGnVMuqO+d+ApKK/omISCXi85IjZhZrZq+bWWRRe/jZHiMiIqGjNHsc9wNDgD+bWQMKV8cVEZFKojSLHO5xzu13zo0B4oD2ZVyTiIgEMZ8XOQTmFH/hnBsLvF3WRYmISPn44osvePTRR8/pOUpyqOq3wMTihnNu5rGdzrn/d04ViIhIufrpp5+YPn06b7zxBgsWLKBKlXO7+asWORQRCVFr1qzh97//PY0aNWLQoEFs27aNZ599lu3bt5/T85Zkj0OLHIqIVBD79u1jypQpvPnmm6xatYrq1avTr18/hg0bhsfjOee9DShZcFxbtMjhiXxa5NDM3gR6A7udc61P0f8r4OGi5gHgPufcmpI8t4hIZZafn8+8efP497//zYwZMzh8+DBt27ZlwoQJDBo0iIYNG5bp65VkyZGqZfRak4BXOP1k+pdAF+fcPjNLoHBepWMZvbaISMj54osvePvtt5k0aRLbt2+nfv363HvvvQwZMoTrrrsOs/KZaSjVleOl4ZxbZGZNztCfcUxzORBR3jWJiFQ0ubm5vP/++7z11lssXboUMyMuLo7x48fTp08fatasWe41+C04fPQbIDnQRYiIBIP8/HzS0tJ46623+PDDDzl06BAtW7bkueeeY9CgQURE+Pfv7KALDjPrSmFw3HSGbYYDwwEuv/xyP1UmIuJfa9as4Z133uE///kPO3fupH79+gwdOpR77rmH9u3bl9uhqLMJquAws7bAG0CCc27v6bZzzk2k6NqSqKgonfUlIiFjx44d/Oc//+Gdd95h7dq1VKtWjYSEBO655x569+5NjRo1Al1i8ASHmV0OTAd+7Zz7PND1iIj4y/fff8+MGTN45513mD9/Ps45OnbsyCuvvMIdd9xBeHh4oEs8jt+Cw8ymAB4g3MyygSeAMADn3GvA40BD4NWi3a8851yUv+oTEfGnI0eOkJyczOTJk5k1axaHDh2iadOm/PnPf2bQoEFcffXVgS7xtPx5VtVdZ+kfBgzzUzkiIn5XUFDAkiVLmDx5Mh988AH79u0jPDyc3/zmNwwcOJBOnToFbN7CF0FzqEpEJBQ551i9ejVTpkwhMTGR7Oxsateuze23387AgQOJjY0lLCws0GX6RMEhIlIONm7cyJQpU5gyZQqbN28mLCyM+Ph4/va3v3Hrrbdy3nnnBbrEUlNwiIiUka1bt5KYmEhiYiJr1qzBzOjatSt/+tOf6Nu3Lw0aNAh0iWVCwSEicg62bdvG+++/T2JiIpmZmQBER0czYcIEBgwYwKWXXhrgCsuegkNExEc7duxg6tSpJCYmsmzZMgCuv/56/v73vzNgwACuuOKKAFdYvhQcIiIlsGPHDqZNm8b777/P0qVLAWjbti3PPPMMAwYMoHnz5gGu0H8UHCIip1EcFh988AFLliwBoE2bNowbN45f/vKXtGjRIsAVBoaCQ0TkGF9//TXTpk1j2rRpZGQULtqtsDiegkNEKr0tW7Ywbdo0pk6dysqVKwGIjIzk6aefpl+/frRs2TLAFQYXc65irxFYt25dd/311x/3vQEDBnD//fdz8OBBevbsedJjBg8ezODBg8nJyaF///4n9d93333ccccdbN++nV//+tcn9Y8ePZpbbrmFTZs2MWLEiJP6//znPxMTE0NWVhajRo06qf/ZZ5/lhhtuICMjg0cfffSk/gkTJhAZGUl6ejpPP/30Sf3/+te/aNGiBbNmzeKFF144qf+dd96hcePGJCYm8s9//vOk/qlTpxIeHs6kSZOYNGnSSf1JSUnUrl2bV199lffff/+k/gULFgAwfvx4Zs+efVxfrVq1SE4uXBF/3LhxzJs377j+hg0bMm3aNAAeeeQR78RisYiICN59910ARo0aRVZW1nH9V199NRMnTgRg+PDhfP758cuaRUZGMmHCBAAGDRpEdnb2cf2dOnXiueeeA6Bfv37s3Xv8Wprdu3fnL3/5CwAJCQn89NNPx/X37t2bMWPGAODxeDiR3nsV4703duxY5s2bR05ODnv27OHHH38EoH379vTv358NGzbw5ZdfHvf4UHvvLVy4cFVpl3XSHoeIVBqZmZlMnz6d119/ne+++w6AevXq0axZM9q0acOMGTMAThm68rMKv8cRFRXlis+dFhE5Vl5eHkuWLGHGjBl8+OGHbNu2japVq9K1a1f69u3LbbfdFpLXWZSEmWmPQ0QE4NChQ6SlpTFjxgw++ugj9u7dS40aNYiLi+PJJ5+kT58+NGzYMNBlVmgKDhGp8Pbt28ecOXOYOXMmycnJ/Pjjj9SrV4/evXtz++2306NHD+rUqRPoMkOGgkNEKqTt27czc+ZMPvzwQxYuXEheXh6XXHIJgwYN4vbbb6dr165Ur1490GWGJAWHiFQIzjk+/fRTZs6cycyZM1m9ejUALVu2ZMyYMdx22220b9+eKlWqBLjS0KfgEJGgdfToURYtWsTMmTP56KOP+PrrrzEzoqOjvcuT64I8/1NwiEhQ2b9/P8nJycyaNYukpCRyc3OpWbMmsbGx/OUvf6F3795cfPHFgS6zUlNwiEjAbdmyhVmzZvHRRx+xePFi8vLyCA8Pp2/fvvTp04fY2NgKfeOjUKPgEBG/y8/PZ/ny5cyaNYtZs2bx2WefAXDNNdcwZswYbrnlFjp27EjVqlUDXKmcioJDRPwiNzeXuXPnMnv2bJKSkti7dy/VqlXj5ptv5t577+WWW26hWbNmgS5TSkDBISLl5vPPP2fOnDnMnj2bRYsWkZeXR4MGDejZsye33HILcXFxXHDBBYEuU3yk4BCRMnPkyBEWL17sDYvNmzcDhYegihdojI6O1iGoCs5vwWFmbwK9gd3Oudan6G8J/BtoBzzmnBvvr9pEpPR27dpFcnIyc+bMITU1lR9++IEaNWrQtWtX/vjHP9KrVy+aNGkS6DKlDPlzj2MS8Arw9mn6vwP+ANzmr4JExHcFBQVkZmYyZ84ckpKSKF5ktFGjRtx111306tWL7t276yyoEOa34HDOLTKzJmfo3w3sNrNe/qpJREpm3759pKamkpSURHJyMnv27KFKlSpER0fzzDPP0LNnT6699lrMLNClih9UyDkOMxsODAe4/PLLA1yNSOgpXt4jKSmJpKQkli1bRn5+Pg0aNCA+Pp5evXrRo0cPrTJbSVXI4HDOTQQmQuH9OAJcjkhIyM3NJT09neTkZFJSUtixYwcA7dq145FHHqFnz5506NBBE9tSMYNDRM5d8V5FcnIyycnJZGRkkJeXx/nnn09MTIx3r6Ky3uhITk/BIVKJ7N+//7i9im+++QYovF/2Qw89REJCAtHR0YSFhQW4Uglm/jwddwrgAcLNLBt4AggDcM69ZmaXAJlAPaDAzEYB1zjnvvdXjSKhpqCggKysLG9QFM9VXHDBBcTExJCQkECPHj1o1KhRoEuVCsSfZ1XddZb+XUCEn8oRCVl79uwhLS2NlJQU5s6dy+7duwG4/vrrGTt2LAkJCXTs2JFq1XTAQUpH7xyRCi4vL48VK1aQkpJCSkoKq1atwjlHw4YNiY+P9/7TUuRSVhQcIhXQ119/zdy5c5k7dy7z5s0jNzeXKlWq0KlTJ5566ini4+Np166dzoCScqHgEKkADh48yMKFC71hsXHjRgAaN27ML3/5S+Lj44mJidGCgeIXCg6RIOScY+3atd6gWLx4MUeOHKFmzZp06dKFESNGEB8fT8uWLXW1tvidgkMkSOzevZv09HTmzp1Lamoqu3btAqBVq1b87ne/Iz4+nptvvplatWoFuFKp7BQcIgFy+PBhMjIyvEHxySefANCgQQNiY2OJj48nLi6Oyy67LMCVihxPwSHiJ845Nm7cSGpqKqmpqSxYsICDBw9SrVo1brjhBp5++mliY2O5/vrrNaktQU3BIVKOcnJymDdvnjcssrOzAWjevDlDhgwhLi4Oj8dDvXr1AlypSMkpOETKUPHhp7S0NNLS0rzXVBRfqR0bG0tsbCxNmzYNdKkipabgEDkHzjnWr1/vDYqFCxd6Dz9FR0fz17/+lbi4OKKionT4SUKGgkPERzt37iQ9PZ309HTS0tLYuXMnAC1atGDo0KHew09169YNcKUi5UPBIXIWBw8eZNGiRaSmppKWlsa6desAaNiw4XGHn3RTMaksFBwiJ8jPz2f16tXew08ZGRkcOXKEGjVqcNNNNzFo0CBiY2OJjIykSpUqgS5XxO8UHCLA1q1bvUExf/589u3bB8C1117LH/7wB2JiYrj55pupXbt2gCsVCTwFh1RKe/fuZf78+d55ii+//BKAiIgIbrvtNmJjY+nWrZtWlBU5BQWHVAqHDh1i6dKl3qBYvXo1zjnq1atH165defDBB4mJiaFFixZa+0nkLBQcEpKK73xXfPbT4sWLOXToENWqVaNTp048+eSTxMbG0r59e93QSMRH+o2RkPHll196g2LevHns3bsXgNatW/Pb3/6WmJgYOnfurNNkRc6RgkMqrGPnKdLT09m6dSsAjRo1olevXsTGxtK9e3cuvfTSAFcqEloUHFJh/PTTTyxZssS7R1E8T1G3bl26du3KAw88oHkKET9QcEjQKr6eoniPYunSpRw+fJiwsDDNU4gEkH7bJGg45/jiiy+8QTF//nz2798PQNu2bRk5cqT3eoo6deoEuFqRystvwWFmbwK9gd3Oudan6DfgJaAncBAY7Jxb7a/6JDC+/fbb4+Yptm3bBhTeS7tv3750796d7t2763oKkSDizz2OScArwNun6U8Arir61xH4Z9H/EkIOHDjA4sWLvUHx6aefAlC/fn26devG2LFjiYmJoXnz5pqnEAlSfltoxzm3CPjuDJvcCrztCi0HLjCzs54Os2nTJiZNmgTA0aNH8Xg8vPvuu0Dh4nQej4fExEQAcnNz8Xg8TJ8+HSi8yY7H42HWrFkA7Nq1C4/HQ0pKCgDbt2/H4/GQnp4OFC5L4fF4WLhwofe1PR4PGRkZAKxbtw6Px8PKlSsByMrKwuPxkJWVBcDKlSvxeDzeRfIyMjLweDxs2rQJgIULF+LxeLxnB6Wnp+PxeNi+fTsAKSkpeDwe772oZ82ahcfjIScnB4Dp06fj8XjIzc0FIDExEY/Hw8GDBwF499138Xg8HD16FIBJkybh8Xi8Y/n6668TExPjbb/66qskJCR42y+99BJ9+vTxtsePH0+/fv287eeff54777zT2x43bhwDBw4kIyODcePGccUVV1CvXj169uzJP/7xD3Jzc2nfvj0rV65kz549NGnShHXr1nHVVVdhZowaNYpRo0Z5n2/kyJGMGTPG2x4+fDiPPPKItz1kyBAef/xxb3vQoEGMGzfO277zzjt5/vnnve1+/foxfvx4b7tPnz689NJL3nZCQgKvvvqqtx0TE8Prr7/ubXs8Hr33gvi9N2jQIG/78ccfZ8iQId72I488wvDhw73tMWPGMHLkSG+7Mrz3zkUwzXFcBmw/pp1d9L2dJ25oZsOB4QA1atTwS3FScuvXr2fevHlMnjyZLVu2MGXKFMyMiy++mFatWvHiiy9y44038tRTT7F3716ioqICXbKI+MCcc/57MbMmwOzTzHHMAZ5zzi0pas8D/uScW3Wm54yKinKZmZnlUK2U1I4dO4678K74/hTNmzene/fuxMTE0LVrVxo2bBjgSkWkmJmtcs6V6q+2YNrjyAYaH9OOAL4JUC1yBrm5uSxcuNC77tPGjRsBCA8P9wZFTEwMTZo0CWyhIlIugik4PgJ+Z2bvUTgpnuucO+kwlfjfkSNHWL58uXev4r///S/5+fnUqlWLLl26MGzYMGJiYmjTpo3uTyFSCfjzdNwpgAcIN7Ns4AkgDMA59xqQROGpuF9QeDrukFM/k5S3goIC1q1b5w2K4vtoV6lShQ4dOjB27FhiY2OJjo7WHJNIJeS34HDO3XWWfgeMPNM2Un62bdt23DzF7t27AWjZsiVDhw4lJiaGLl26cMEFFwS4UhEJtGA6VCV+tH//fj7++GNvWHz++ecAXHLJJcTFxRETE0P37t2JiIgIcKUiEmwUHJXE4cOHWb58OWlpaaSnp7Ny5UoKCgo477zz6NKlC/fddx8xMTG0atVKF96JyBkpOEKUc45169Z5g6J4nqJq1ap06NCBxx57jNjYWDp27Ej16tUDXa6IVCAKjhBSfD1FcVh8++23ALRo0cI7T+HxeDj//PMDXKmIVGQKjgrswIEDLFiwgLS0NNLS0tiwYQMAF154ITExMcTGxhITE0Pjxo3P8kwiIiWn4KhA8vLyyMzM9AbFsmXLyMvLo2bNmnTu3JmhQ4cSGxur6ylEpFwpOILc1q1bSU1NJS0tzXt/CmmmON0AAAhoSURBVDPjuuuuY/To0cTGxnLjjTdSs2bNQJcqIpWEgiPI7N+/n/nz55OWlkZqaqp3tdLGjRvTr18/7320w8PDA1ypiFRWCo4Ay8vLY8WKFd6gWLFiBQUFBdSpU8d7H+3Y2FiuvvpqnSYrIkFBwREAW7ZsITU1ldTUVObPn8/3339PlSpVaN++vfc02ejoaMLCwgJdqojISRQcfpCbm8v8+fO9YVF8+KlJkybceeedxMXF0a1bN+rXrx/gSkVEzk7BUQ7y8/PJzMxk7ty5pKamsnz5cvLz86lTpw7dunXjwQcfJC4uTrdHFZEKScFRRrKzs5k7dy5z584lPT2dffv2YWa0a9eOhx9+mPj4eKKjo3WVtohUeAqOUjp06BCLFi1i7ty5pKSk8NlnnwFw6aWXcuuttxIfH09MTIzOfhKRkKPgKCHnHJs2bSIlJYW5c+eyYMECDh06RPXq1encuTNDhgwhPj6e1q1b6/CTiIQ0BccZfP/998yfP5/k5GRSUlLYtm0bULj20/Dhw+nRowddunShdu3aAa5URMR/FBzHcM6xZs0aUlJSSElJYenSpeTl5VG3bl26devGo48+Snx8vO6lLSKVWqUPjn379pGamuoNi127dgEQGRnJmDFj6NGjB506ddKktohIkUoXHAUFBXzyySckJyeTnJzM8uXLKSgooH79+sTFxdGjRw/i4+O59NJLA12qiEhQqhTBsX//flJTU71hUXyfiqioKB577DESEhJo37491apViuEQETknIflJ6Zxj7dq1JCUlkZSUREZGBvn5+dSvX5/4+HgSEhLo0aMHF110UaBLFRGpcEImOA4cOMC8efOYM2cOSUlJ7NixA4DrrruOsWPH0rNnTzp06KC9ChGRc+TXT1Ez6wG8BFQF3nDOPX9Cf33gTaAZcAgY6pxbd6bn3L17N3FxcSxcuJAjR45Qt25dYmNj6dmzJwkJCTRq1KicfhoRkcrJnHP+eSGzqsDnQCyQDawE7nLOfXbMNv8DHHDO/dXMWgL/cM51P8vzupYtW9KrVy969uzJTTfdpDOgRETOwsxWOeeiSvNYf+5xdAC+cM5tBTCz94Bbgc+O2eYa4DkA59xGM2tiZhc757493ZO2bt2atWvXlmPZIiJyLH8Gx2XA9mPa2UDHE7ZZA/QFlphZB+AKIAI4LjjMbDgwvKh52MzOeDirEgkHcgJdRJDQWPxMY/EzjcXPWpT2gf4MjlMt4HTicbLngZfMLAtYC3wC5J30IOcmAhMBzCyztLtboUZj8TONxc80Fj/TWPzMzDJL+1h/Bkc20PiYdgTwzbEbOOe+B4YAWOFKgV8W/RMRkSBRxY+vtRK4ysyamll14E7go2M3MLMLivoAhgGLisJERESChN/2OJxzeWb2O2AuhafjvumcW29mvy3qfw34BfC2meVTOGn+mxI89cTyqrkC0lj8TGPxM43FzzQWPyv1WPjtdFwREQkN/jxUJSIiIUDBISIiPqkwwWFmPcxsk5l9YWZjT9FvZvZyUf+nZtYuEHX6QwnG4ldFY/CpmWWY2bWBqNMfzjYWx2zX3szyzay/P+vzp5KMhZl5zCzLzNab2UJ/1+gPJfj9ON/MZpnZmqJxGBKIOv3BzN40s92nu9at1J+bzrmg/0fhZPoW4EqgOoUXCl5zwjY9gWQKrxeJBlYEuu4AjsUNQP2irxMq81gcs918IAnoH+i6A/i+uIDCk04uL2pfFOi6AzQOjwJ/K/r6QuA7oHqgay+n8egMtAPWnaa/VJ+bFWWPw7tciXPuCFC8XMmxbgXedoWWAxeYWSjejemsY+Gcy3DO7StqLqfwmplQVJL3BcDvgWnAbn8W52clGYuBwHTn3DYA51wojkdJxsEBdYuuFatDYXCcdKFxKHDOLaLw5zudUn1uVpTgONVyJZeVYptQ4OvP+RsK/6IIRWcdCzO7DLgdeM2PdQVCSd4XVwP1zWyBma0ys7v9Vp3/lGQcXqHw1P9vKFyh4o/OuQL/lBd0SvW5WVFuTlGS5UpKsk0oKPHPaWZdKQyOm8q1osApyVhMAB52zuUX/oEZskoyFtWA64HuQC1gmZktd859Xt7F+VFJxiEeyAK6UXgLhzQzW+wq58XGpfrcrCjBcdblSkq4TSgo0c9pZm2BN4AE59xeP9XmbyUZiyjgvaLQCAd6mlmec+5D/5ToNyX9Hclxzv0I/Ghmi4BrKbzdQagoyTgMAZ53hQf5vzCzL4GWwH/9U2JQKdXnZkU5VHXW5UqK2ncXnSUQDeQ653b6u1A/KMnSLZcD04Ffh9hfkyc661g455o655o455oAU4H7QzA0oGS/IzOBm82smpnVpnB16g1+rrO8lWQctlG414WZXUzhKrFb/Vpl8CjV52aF2ONwJVuuJInCMwS+AA5StFhiqCnhWDwONAReLfpLO8+F4IqgJRyLSqEkY+Gc22BmKcCnQAGFd+EMqVsSlPA9MQ6YZGZrKTxU87BzLiSXWjezKYAHCDezbOAJIAzO7XNTS46IiIhPKsqhKhERCRIKDhER8YmCQ0REfKLgEBERnyg4RETEJwoOERHxiYJDpIyY2cdmFlv09dNm9nKgaxIpDxXiAkCRCuIJ4Ckzuwi4DugT4HpEyoUuABQpQ0U3R6oDeJxzP5jZlcBjwPnOuZC9iZRULjpUJVJGzKwNcClw2Dn3A0DRfSF+E9jKRMqWgkOkDBTd/GYyhTfG+dHM4gNckki5UXCInKOilWanA6OdcxsoXETvyYAWJVKONMchUo7MrCHwDBBL4Wq0zwW4JJFzpuAQERGf6FCViIj4RMEhIiI+UXCIiIhPFBwiIuITBYeIiPhEwSEiIj5RcIiIiE8UHCIi4hMFh4iI+OT/A9ScbL6ixECzAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 51\n",
    "ue_pw = []\n",
    "ae_pw = []\n",
    "\n",
    "gamma2_gamma1 = []\n",
    "gamma2_gamma1_wilson = []\n",
    "gamma2_gamma1_nrtl = []\n",
    "x1 = []\n",
    "alpha = 0.3\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",
    "    gamma2_gamma1.append(sac.seg_gamma[1][0]/sac.seg_gamma[0][0])\n",
    "    gamma2_gamma1_wilson.append(1)\n",
    "    gamma2_gamma1_nrtl.append(math.pow(sac.psi[0][0][0][0]/sac.psi[0][1][0][0], 1-alpha))\n",
    "\n",
    "\n",
    "plt.plot(x1, gamma2_gamma1, 'k-', label='SAC')\n",
    "plt.plot(x1, gamma2_gamma1_nrtl, 'k--', label='NRTL')\n",
    "plt.plot(x1, gamma2_gamma1_wilson, 'k:', label='Wilson')\n",
    "\n",
    "plt.xlabel(r'$x_1$')\n",
    "plt.ylabel(r'$\\Gamma_2/\\Gamma_1$')\n",
    "plt.legend(loc='upper left')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0.9, 1.6])\n",
    "\n",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.plot()\n",
    "plt.savefig('pdf/prob-symmetric-nearly-ideal.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8baf968a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}