free_energy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate Free Energies\n",
    "Based on: https://github.com/freitas-rodrigo/FreeEnergyLAMMPS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyiron import Project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas\n",
    "from jinja2 import Template\n",
    "import matplotlib.pyplot as plt \n",
    "import scipy.constants as sc\n",
    "from scipy.integrate import cumtrapz\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pr = Project(\"free\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pr.remove_jobs_silently(recursive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "fe_potential = pandas.DataFrame({\n",
    "  'Name': ['Fe_Potential'],\n",
    "  'Filename': [[os.path.abspath(\"Fe.eam\")]],\n",
    "  'Model': ['Custom'],\n",
    "  'Species': [['Fe']],\n",
    "  'Config': [['pair_style eam\\n', \n",
    "              'pair_coeff * * Fe.eam\\n']]\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Frenkel-Ladd thermodynamic integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature_lst = [100, 400, 700, 1000, 1300, 1600]\n",
    "alat_lst = [2.8841, 2.9115, 2.9315, 2.9484, 2.9637, 2.9782]\n",
    "k_lst = [5.787, 4.866, 4.073, 3.373, 2.799, 2.443]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "lmp_frenkel_ladd_input_str = \"\"\"\\\n",
    "# This script runs a Frenkel-Ladd thermodynamic integration for the calculation of the free energy of bcc iron. See job.sh for the variables that need to be passed in the command line.\n",
    "\n",
    "\n",
    "#--------------------------- Simulation variables -----------------------------#\n",
    "  # Initalizes the random number generator.\n",
    "  variable         rnd equal round(random(0,999,{{ RANDOM }}))\n",
    "\n",
    "  # Simulation control parameters.\n",
    "  variable         t_eq equal 1000 # Equilibration time.\n",
    "  variable         t equal 5000 # Switching time.\n",
    "#------------------------------------------------------------------------------#\n",
    "\n",
    "\n",
    "#---------------------------- Atomic setup ------------------------------------#\n",
    "  units            metal\n",
    "\n",
    "  # Create atoms.\n",
    "  read_data        structure.inp\n",
    "\n",
    "  # Define interatomic potential.\n",
    "  include          potential.inp\n",
    "  neigh_modify     delay 0\n",
    "  timestep         0.002\n",
    "#------------------------------------------------------------------------------#\n",
    "\n",
    "\n",
    "#----------------------------- Run simulation ---------------------------------#\n",
    "  # Integrator, Frenkel-Ladd modifier, and thermostat.\n",
    "  fix              f1 all nve \n",
    "  fix              f2 all ti/spring {{ k }} ${t} ${t_eq} function 2\n",
    "  fix              f3 all langevin {{ T }} {{ T }} 0.1 ${rnd} zero yes\n",
    "\n",
    "  # Compute temperature using COM coordinates.\n",
    "  compute          c1 all temp/com\n",
    "  fix_modify       f3 temp c1\n",
    "\n",
    "  # Redefine simulation time to capture first and last steps on output files.\n",
    "  # Notice that this does not affect the switching time.\n",
    "  variable         t_eq equal ${t_eq}-1\n",
    "  variable         t equal ${t}+1\n",
    "\n",
    "  # Setup thermo output.\n",
    "  thermo_style     custom step pe\n",
    "  thermo           0\n",
    "  variable         dE equal pe-f_f2\n",
    "  variable         lambda equal f_f2[1]\n",
    "\n",
    "  # Initial temperature to accelerate equilibration.\n",
    "  variable         rnd equal round(random(1,666,0)) # Generates new rnd #.\n",
    "  variable         T0 equal 2*{{ T }}\n",
    "  velocity         all create ${T0} ${rnd} dist gaussian\n",
    "\n",
    "  # Forward integration.\n",
    "  run              ${t_eq}\n",
    "  fix              f4 all print 1 \"${dE} ${lambda}\" title \"# dE [eV] lambda\" &\n",
    "                   screen no file forward.dat\n",
    "  run              ${t}\n",
    "  unfix            f4\n",
    "\n",
    "  # Backward integration. \n",
    "  run              ${t_eq}\n",
    "  fix              f4 all print 1 \"${dE} ${lambda}\" title \"# dE [eV] lambda\" &\n",
    "                   screen no file backward.dat\n",
    "  run              ${t}\n",
    "#------------------------------------------------------------------------------#\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "tm_frenkel_ladd = Template(lmp_frenkel_ladd_input_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job lmp_frenkel_ladd_100 was saved and received the ID: 1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/pyiron/lammps/base.py:997: UserWarning: LAMMPS warning: No dump.out output file found.\n",
      "  warnings.warn(\"LAMMPS warning: No dump.out output file found.\")\n",
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/pyiron/lammps/base.py:583: UserWarning: LAMMPS warning: log.lammps does not contain the required pressure values.\n",
      "  warnings.warn(\"LAMMPS warning: log.lammps does not contain the required pressure values.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job lmp_frenkel_ladd_400 was saved and received the ID: 2\n",
      "The job lmp_frenkel_ladd_700 was saved and received the ID: 3\n",
      "The job lmp_frenkel_ladd_1000 was saved and received the ID: 4\n",
      "The job lmp_frenkel_ladd_1300 was saved and received the ID: 5\n",
      "The job lmp_frenkel_ladd_1600 was saved and received the ID: 6\n"
     ]
    }
   ],
   "source": [
    "for temp, alat, k in zip(temperature_lst, alat_lst, k_lst):\n",
    "    job_lmp = pr.create_job(pr.job_type.Lammps, \"lmp_frenkel_ladd_\" + str(temp))\n",
    "    fe_bcc = pr.create_ase_bulk(\"Fe\", a=alat ,cubic=True).repeat(5)\n",
    "    job_lmp.structure = fe_bcc\n",
    "    job_lmp.potential = fe_potential\n",
    "    job_lmp.input.control.load_string(\n",
    "        tm_frenkel_ladd.render(\n",
    "            RANDOM=np.random.randint(999), \n",
    "            T=temp, \n",
    "            k=k\n",
    "        )\n",
    "    )\n",
    "    job_lmp.run()\n",
    "    job_lmp.decompress()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = np.array(temperature_lst) # [K].\n",
    "a = np.array(alat_lst) # [A].\n",
    "k = np.array(k_lst) # [eV/A^2]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 55.845 # Iron mass [g/mol].\n",
    "natoms = len(fe_bcc) # Number of atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "250"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "natoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "kB = sc.value('Boltzmann constant in eV/K')\n",
    "eV = sc.value('electron volt')\n",
    "hbar = sc.value('Planck constant over 2 pi in eV s')\n",
    "mu = sc.value('atomic mass constant')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = np.zeros(len(T)) # Reversible work for each temperature.\n",
    "for i, job in enumerate(pr.iter_jobs(convert_to_object=False)):\n",
    "    if \"lmp_frenkel_ladd_\" in job.job_name:\n",
    "        dE, lamb = np.loadtxt(os.path.join(job.working_directory, \"forward.dat\"), unpack=True)\n",
    "        I_forw = np.trapz(dE,lamb)\n",
    "        dE, lamb = np.loadtxt(os.path.join(job.working_directory, \"backward.dat\"), unpack=True)\n",
    "        I_back = np.trapz(dE,lamb)\n",
    "        W[i] = (I_forw-I_back) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define harmonic reference system free energy [Eq.(15) in the paper].\n",
    "omega = np.sqrt(k*eV/(m*mu)) * 1.0e+10 # [1/s].\n",
    "F_harm = 3*natoms*kB*T * np.log(hbar*omega/(kB*T)) # [eV]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fixed center of mass correction [Eq.(24) in the paper].\n",
    "V = (a**3/2) * natoms # Total volume.\n",
    "F_CM = (kB*T)*np.log((natoms/V) * (2*np.pi*kB*T / (natoms*m*omega**2))**(3/2)) # [eV]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute absolute free energy per atom [Eq.(16) in the paper] and save data.\n",
    "F_fl = (F_harm + W + F_CM) / natoms # [eV/atom].\n",
    "T_fl = T.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = ['#E41A1C','#377EB8','#4DAF4A','#984EA3','#FF7F00','#FFFF33','#A65628','#F781BF','#999999']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Free energy [eV/atom]')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAETCAYAAAD+spv+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdvklEQVR4nO3de5gddZ3n8fdnGpAmGy6SBOhc7OijKAOK2LZBR+USWAxMMNGZB2fYwdUx7owXxGXTDXEU91mdBC9xHGcdA6ijqKjQUWRBJCqg7nLpBJokAiKCId1kaEYNTYxAmu/+UXXg0J5zurq76lzSn9fz9NNVv6pT9elA55tfXX4/RQRmZmZ5+JNGBzAzs72Hi4qZmeXGRcXMzHLjomJmZrlxUTEzs9zs0+gARZg1a1Z0dnY2OoaZ2V5p48aNj0bE7Erb9sqi0tnZSX9/f6NjmJntlST9uto2X/4yM7PcuKiYmVluXFTMzCw3LipmZpYbFxUzM8uNi4qZmeXGRSW1q289O7oXMThvATu6F7Grb32jI5mZtZy98j2VidrVt56dK3uI3bsBGB0cZOfKHgBmLF/WyGhmZi3FPRVgZPWaZwpKSezezcjqNQ1KZGbWmlxUgNGhoQm1m5lZZS4qQFtHx4TazcysMhcVYGZvD2pvf06b2tuZ2dvToERmZq3JN+p59mb8yOo1jA4N0dbRwczeHt+kNzObIBeV1Izly1xEzMymqKGXvySdLykkzaqwbX9Jt0kakLRV0kcbkdHMzLJrWE9F0nzgFGBblV2eAE6KiMcl7Qv8VNJ1EXFL3UKamdmENLKnshZYCUSljZF4PF3dN/2quK+ZmTWHhhQVSUuBwYgYGGe/Nkl3Ao8AN0TErfXIZ2Zmk1PY5S9JG4DDK2xaBVwInDreMSJiFDhW0sHAeklHR8SWKudbAawAWLBgwWRjm5nZFBRWVCJicaV2SccAC4EBSQDzgE2SuiNiR5Vj/U7SjcBpQMWiEhHrgHUAXV1dvkxmZtYAdb/8FRGbI2JORHRGRCewHThubEGRNDvtoSCpHVgM3FPvvGZmll1TvVEvqUPStenqEcCPJd0F3E5yT+WaxqUzM7PxNPzlx7S3UloeApaky3cBr2xQLDMzm4Sm6qmYmVlrc1Exz3ppZrlp+OUvayzPemlmeXJPZZrzrJdmlicXlWnOs16aWZ5cVKY5z3ppZnlyUZnmPOulmeXJN+qnOc96aWZ5clExz3ppZrnx5S8zM8uNi4qZmeXGRcXMzHLjomJmZrlxUTEzs9y4qJiZWW5cVMzMLDcuKmZmlhsXFTMzy01Di4qk8yWFpFk19mmTdIckz09vZtbkGlZUJM0HTgG2jbPrucDdxScyM7OpamRPZS2wEohqO0iaB5wOXFqvUGZmNnkNKSqSlgKDETEwzq6fISk8T2c45gpJ/ZL6h4eHc0hpZmYTVdgoxZI2AIdX2LQKuBA4dZzPnwE8EhEbJZ0w3vkiYh2wDqCrq6tq78fMzIpTWFGJiMWV2iUdAywEBiQBzAM2SeqOiB1lu74OWCppCbA/cKCkyyPi7KIym5nZ1NT98ldEbI6IORHRGRGdwHbguDEFhYi4ICLmpfucBfzIBcXMrLk11XsqkjokXdvoHGZmNjkNn/kx7YmUloeAJRX2uRG4sW6hzMxsUpqqp2JmZq3NRcXMzHLjomJmZrlxUTEzs9y4qJiZWW5cVMzMLDcuKmZmlhsXFTMzy42LipmZ5abqG/WSlmf4/B8iwsOqmJkZUHuYlkuA7wKqsc8bABcVMzMDaheV6yLiHbU+LOnynPOY1dWuvvWMrF7D6NAQbR0dzOztYcbyZY2OZdayqhaVLMPMeyh6a2W7+tazc2UPsXs3AKODg+xc2QPgwmI2SeOOUiypjWSe+M7y/SPi08XFMiveyOo1zxSUkti9m5HVa1xUzCYpy9D33wP+AGwmw1zxZq1idGhoQu1mNr4sRWVeRLy88CRmddbW0cHo4GDFdjObnCzvqVwn6dTCk5jV2czeHtTe/pw2tbczs7enQYnMWl+WnsotwHpJfwI8RfKIcUTEgYUmMytY6b6Jn/4yy0+WovIp4Hhgc0REnieXdD7wCWB2RDxaYfuDwAgwCuyJiK48z282Y/kyFxGzHGUpKvcBWwooKPOBU4Bt4+x6YqWCY2ZmzSdLUXkYuFHSdcATpcYcHileC6wkeWvfzMz2Allu1D8A/BDYD5hZ9jVpkpYCgxExMM6uAfxA0kZJK8Y55gpJ/ZL6h4eHpxLPzMwmadyeSkR8FEDSzGQ1Hs9yYEkbgMMrbFoFXAhkeaLsdRExJGkOcIOkeyLi5io51wHrALq6unK9VGdmZtlkeaP+aOCrwPPT9UeBv4mIrbU+FxGLqxzvGGAhMCAJYB6wSVJ3ROwYc4yh9PsjktYD3UDFomJmZo2X5fLXOuCDEfGCiHgB8N9JRjCelIjYHBFzIqIzIjqB7cBxYwuKpBlp7whJM0h6Nlsme14zMytelqIyIyJ+XFqJiBuBGUWEkdQhqTSU/mHATyUNALcB/ycivl/Eec3MLB9Znv76laR/ILkEBnA2yc37XKS9ldLyELAkXf4V8Iq8zmNmZsXL0lN5BzAb6Eu/ZgFvLzCTmZm1qCw9lcUR8f7yBkl/AXy7mEhmZtaqsvRULsjYZmZm01zVnoqkN5Hc35gr6bNlmw4E9hQdzMzMWk+ty19DQD+wFNhY1j4CnFdkKDMza0215qgfIHlB8esR8VQdM5mZWYvKcqO+U9I/AkcB+5caI+KFhaUyM7OWlOVG/ZeAz5PcRzkR+ArPvrNiZmb2jCxFpT0ifggoIn4dERcBJxUby8zMWlGWy19/SKcSvk/Se4FBYE6xsczMrBVl6al8ADgAeD/wKpJhWs4pMJOZmbWoWu+pXAB8PyJuT5seB/5rXVKZmVlLqnX56wHgXEmvAAaA64AfRMRv65LMzMxaTq33VK4ArgCQ9ErgNKBPUhuwgaQXc1tdUpqZWUuoek9F0r6l5Yi4IyL+MSJOBM4AtgJ/W4d8ZmbWQmrdqB+UdImkk5TO+wsQEY9FxFURsaIO+czMrIXUKiovIxn76x+AhyR9RtJr6hPLzMxaUdWiEhH/ERFfSC95dZPcuP+MpPslfSyPk0s6X1JImlVl+8GSrpR0j6S7JR2fx3nNzKwYWd5TKU3zexnJcC0j5HA/RdJ84BRgW43d/onkgYCXkkwtfPdUz2tmZsWpWVQk7S/pLyT1AfcDJ5NM0NWRw7nXAiuBqHLuA4E3kBQzIuLJiPhdDuc1M7OC1Hr58evAYuBm4OvAX0XEH/I4qaSlwGBEDJQ9AzDWC4Fh4EvpuzIbgXMjYlceGczMLH+1Xn68Hnh3RIxM5sCSNgCHV9i0CrgQODVDtuOA90XErZL+CegleXCg0vlWACsAFixYMJnIZmY2RYqoePXp2R2kw4CPA3Mj4jRJRwHHR8RlkzqhdAzwQ+D3adM8klkmuyNiR9l+hwO3RERnuv56oDciTh/vHF1dXdHf3z+ZeGZmNg5JGyOiq9K2LDfqv0zSazkiXf8FySCTkxIRmyNiTkR0pgVjO3BceUFJ99tB8ijzkWnTycDPJ3teMzMrXpaiMisivgU8DRARe4DRIsJI6pB0bVnT+4CvSboLOJakx2RmZk0qy3wquyQdSvqUlqRFwM68ApQub6XLQ8CSsvU7gYpdLDMzaz5ZisoHgauBF0n6GTAbeGuhqczMrCWNW1QiYpOkNwJHAgLujYinCk9mZmYtp9YoxceVliNiT0RsjYgt5QWlfB8zM7NaPZUvSTqBpHdSzWXAK/MMZGZmratWUTmI5C32WkVlON84ZlYvu/rWM7J6DaNDQ7R1dDCzt4cZy5c1Opa1uFozP3bWMYeZ1dGuvvXsXNlD7N4NwOjgIDtX9gC4sNiUZBql2Mz2LiOr1zxTUEpi925GVq9pUCLbW7iomE1Do0NDE2o3y8pFxWwaauuoPHtFtXazrMYtKpKuknS6JBcgs73EzN4e1N7+nDa1tzOzt6dBiWxvkaVQfB74K+A+SaslvbTgTGZWsBnLl3HQxWtomzsXJNrmzuWgi9f4Jr1N2bhD3z+zo3QQ8DaS+VAeAi4BLm/Gt+s99L2ZWXGmOvQ96YCSbyeZm/4OkrnjjwNuyCmjmZntBcYd+yudn/6lwFeBP4+Ih9NN35Tk7oCZmT0jyyjFn4uIH1XaUK37Y2Zm01OWonKwpOVj2nYCmyPikQIymZlZi8pSVN4JHA/8OF0/AbgFeImk/xkRXy0om5mZtZgsReVp4GUR8e8Akg4jecz4NcDNJPdazMzMMj391VkqKKlHgJdExG+AKT1OLOl8SSFpVoVtR0q6s+zrMUkfmMr5zMysWFl6Kj+RdA3w7XT9LcDNkmYAv5vsiSXNB04BtlXaHhH3Asem+7YBg8D6yZ7PzMyKl2U64b+X9Bbgz0jmVvkKcFUkb02eOIVzrwVWAt/NsO/JwP0R8espnM/MzApWs6ik433dFRFHA1fldVJJS4HBiBiQas0B9oyzgG+Mc8wVwAqABQsWTDmjmZlNXM2iEhFPSxqQtCAiKl6mqkbSBuDwCptWARcCp2Y8zn7AUuCCcbKuA9ZBMkzLRLKamVk+stxTOQLYKuk2YFepMSKW1vpQRCyu1C7pGGAhUOqlzAM2SeqOiB0VPvImYNOYhwXMzKwJZSkqH83zhBGxGZhTWpf0INAVEY9W+cjbGOfSl5mZNYdxHymOiJuAB4F90+XbgU1FhJHUIenasvUDSJ4Q6yvifGZmlq8sk3S9C7gS+ELaNBf4Tl4BIqKz1EuJiKGIWFK27fcRcWhE7MzrfGZmVpwsLz++B3gd8BhARNxH2eUrMzOzkixF5YmIeLK0ImkfwE9XmZnZH8lSVG6SdCHQLukUkjfrv1dsLDMza0VZikovMAxsBt4NXAt8qMhQZmbWmrIM0/I0yXz0lxQfx8zMWlmW6YRfB1wEvCDdX0BExAuLjWZmZq0my8uPlwHnARuB0WLjmJlZK8tSVHZGxHWFJzEzs5aXpaj8WNInSN5qf6LUGBGFvFVvZmatK0tReU36vausLYCT8o9jZmatLMvTX1OZiMvMzKaRLGN/HSbpMknXpetHSXpn8dHMzKzVZHn58cvA9UBHuv4L4AMF5TEzsxaWpajMiohvAU8DRMQe/GixmZlVkKWo7JJ0KOkgkpIWAR6K3szM/kiWp78+CFwNvEjSz4DZwFsLTWVmZi0py9NfmyS9ETiSZIiWeyPiqcKTmZlZy8nSUyndR9lacBYzM2txWe6pFEbS+ZJC0qwq28+TtFXSFknfkLR/vTOamVl2DSsqkuYDpwDbqmyfC7wf6IqIo4E24Kz6JTQzs4nK8vKjJJ0t6cPp+gJJ3Tmcey2wktpTE+9DMuPkPsABwFAO5zUzs4Jk6an8b+B44G3p+gjwL1M5qaSlwGBEDFTbJyIGgU+S9GQeJhkt+Qc1jrlCUr+k/uHh4anEMzOzScpSVF4TEe8B/gAQEb8F9hvvQ5I2pPdCxn6dCawCPjzO5w8BzgQWkrzNP0PS2dX2j4h1EdEVEV2zZ8/O8GOZmVnesjz99ZSkNp59+XE26dv1tUTE4krtko4hKRQDkgDmAZskdUfEjrJdFwMPRMRw+rk+4LXA5Rkym5lZA2QpKp8F1gOHSfoYyYuPH5rsCSNiMzCntC7pQZKb8Y+O2XUbsEjSAcBu4GSgf7LnNTOz4mV5+fFrkjaS/KUu4M0RcXcRYSR1AJdGxJKIuFXSlcAmYA9wB7CuiPOamVk+Mr38CMwCfh8RX5I0W9LCiHggjwAR0Vm2PAQsKVv/CPCRPM5jZmbFy/JI8UeAHuCCtGlffF/DzMwqyPL01zJgKbALnulNzCwylJmZtaYsReXJiAieffprRrGRzMysVWUpKt+S9AXgYEnvAjYAlxQby8ysOLv61rOjexGD8xawo3sRu/rWNzrSXqPmjXolL5J8E3gp8BjJ8Pcfjogb6pDNzCx3u/rWs3NlD7F7NwCjg4PsXNkDwIzlyxoZba+g5MpWjR2kjRHxqjrlyUVXV1f09/uVFjP7Yzu6FzE6OPhH7W1z53L4bbc0IFHrSetCV6VtWS5/3SLp1TlnMjNriNGhyuPSVmu3iclSVE4kKSz3S7pL0mZJdxUdzMysCG0dHRNqt4mpek9F0oKI2Aa8qY55zMwKNbO35zn3VADU3s7M3p4Gptp71LpR/x3guIj4taSrIuItdcpkZlaY0s34kdVrGB0aoq2jg5m9Pb5Jn5NaRUVlyy8sOoiZWb3MWL7MRaQgte6pRJVlMzOzimr1VF4h6TGSHkt7uky6HhFxYOHpzMyspVQtKhHRVs8gZmbW+rI8UmxmZpaJi4qZmeXGRcXMzHLT0KIi6XxJIWlWle3nStoiaaukD9Q5npmZTVDDioqk+cApwLYq248G3gV0A68AzpD04volNDOziWpkT2UtsJLq78C8DLglIn4fEXuAm0hmoTQzsybVkKIiaSkwGBEDNXbbArxB0qGSDgCWAPNrHHOFpH5J/cPDwzknNjOzLGpO0jUVkjYAh1fYtAq4EDi11ucj4m5Ja4AbgMeBAWBPjf3XAesgmU9lkrHNzGwKCisqEbG4UrukY4CFwEAysSTzgE2SuiNix5hjXAZcln7u48D2ovKamdnUFVZUqomIzcCc0rqkB4GuiHh07L6S5kTEI5IWAMuB4+sW1MzMJqyp3lOR1CHp2rKmqyT9HPge8J6I+G2DopmZWQZ176mMFRGdZctDJDfkS+uvb0QmMzObnKbqqZiZWWtzUTEzs9y4qJiZWW5cVMzMLDcuKmZmlhsXFTMzy42LipmZ5cZFxczMcuOiYmZmuXFRMTOz3LiomJlZblxUzMwsNy4qZmaWGxcVMzPLjYuKmZnlxkXFzMxy46JiZma5cVExM7PcNKSoSLpI0qCkO9OvJVX2O03SvZJ+Kam33jnNzGxiGjlH/dqI+GS1jZLagH8BTgG2A7dLujoifl6vgGZmNjHNfPmrG/hlRPwqIp4ErgDObHAmMzOroZFF5b2S7pL0RUmHVNg+F3iobH172laRpBWS+iX1Dw8P553VzMwyKKyoSNogaUuFrzOBzwMvAo4FHgY+VekQFdqi2vkiYl1EdEVE1+zZs/P4EczMbIIKu6cSEYuz7CfpEuCaCpu2A/PL1ucBQzlEMzOzgjTq6a8jylaXAVsq7HY78GJJCyXtB5wFXF2PfGZmNjmNevrrYknHklzOehB4N4CkDuDSiFgSEXskvRe4HmgDvhgRWxuU18zMMmhIUYmI/1KlfQhYUrZ+LXBtvXKZmdnUNPMjxWZm1mJcVMzMLDcuKmZmlhsXFTMzy42LipmZ5cZFxcxsmtjVt54d3YsYnLeAHd2L2NW3PvdzNHKUYjMzq5NdfevZubKH2L0bgNHBQXau7AFgxvJluZ3HPRUzs2lgZPWaZwpKSezezcjqNbmex0XFzGwaGB2qPHRitfbJclExM5sG2jo6JtQ+WS4qZmbTwMzeHtTe/pw2tbczs7cn1/P4Rr2Z2TRQuhk/snoNo0NDtHV0MLO3J9eb9OCiYmY2bcxYviz3IjKWL3+ZmVluXFTMzCw3LipmZpYbFxUzM8uNi4qZmeXGRcXMzHKjiGh0htxJGgZ+3egcZWYBjzY6RAbOmS/nzJdz5meqGV8QEbMrbdgri0qzkdQfEV2NzjEe58yXc+bLOfNTZEZf/jIzs9y4qJiZWW5cVOpjXaMDZOSc+XLOfDlnfgrL6HsqZmaWG/dUzMwsNy4qZmaWGxeVKZI0X9KPJd0taaukc9P250u6QdJ96fdDyj5zgaRfSrpX0n+uc942SXdIuqZZc0o6WNKVku5J/1yPb9Kc56X/zbdI+oak/Zshp6QvSnpE0paytgnnkvQqSZvTbZ+VpDrk/ET63/0uSeslHdyMOcu2nS8pJM1q1pyS3pdm2Srp4sJzRoS/pvAFHAEcly7PBH4BHAVcDPSm7b3AmnT5KGAAeB6wELgfaKtj3g8CXweuSdebLifwb8Dfpsv7AQc3W05gLvAA0J6ufwt4ezPkBN4AHAdsKWubcC7gNuB4QMB1wJvqkPNUYJ90eU2z5kzb5wPXk7xoPasZcwInAhuA56Xrc4rO6Z7KFEXEwxGxKV0eAe4m+QvnTJK/HEm/vzldPhO4IiKeiIgHgF8C3fXIKmkecDpwaVlzU+WUdCDJL8dlABHxZET8rtlypvYB2iXtAxwADDVDzoi4GfjNmOYJ5ZJ0BHBgRPy/SP6m+UrZZwrLGRE/iIg96eotwLxmzJlaC6wEyp92aracfwesjogn0n0eKTqni0qOJHUCrwRuBQ6LiIchKTzAnHS3ucBDZR/bnrbVw2dIfgmeLmtrtpwvBIaBL6WX6S6VNKPZckbEIPBJYBvwMLAzIn7QbDnLTDTX3HR5bHs9vYPkX8rQZDklLQUGI2JgzKamygm8BHi9pFsl3STp1UXndFHJiaT/BFwFfCAiHqu1a4W2wp/rlnQG8EhEbMz6kQpt9Xj+fB+SLvznI+KVwC6SyzXVNOrP8xCSf+0tBDqAGZLOrvWRCm3N8Dx/tVwNzStpFbAH+FqpqUqeuueUdACwCvhwpc1V8jTy9+kQYBHwP4BvpfdICsvpopIDSfuSFJSvRURf2vzvaVeS9Hup27md5FpsyTySyyZFex2wVNKDwBXASZIub8Kc24HtEXFrun4lSZFptpyLgQciYjgingL6gNc2Yc6SiebazrOXnsrbCyfpHOAM4K/TSzDNlvNFJP+YGEh/n+YBmyQd3mQ5Sc/bF4nbSK5SzCoyp4vKFKVV/zLg7oj4dNmmq4Fz0uVzgO+WtZ8l6XmSFgIvJrkxVqiIuCAi5kVEJ3AW8KOIOLsJc+4AHpJ0ZNp0MvDzZstJctlrkaQD0v8HTia5n9ZsOUsmlCu9RDYiaVH68/1N2WcKI+k0oAdYGhG/H5O/KXJGxOaImBMRnenv03aSh3V2NFPO1HeAkwAkvYTkwZdHC82Z59MH0/EL+DOS7uFdwJ3p1xLgUOCHwH3p9+eXfWYVydMW95LzEyAZM5/As09/NV1O4FigP/0z/Q5J970Zc34UuAfYAnyV5EmahucEvkFyn+cpkr/w3jmZXEBX+rPdD3yOdASOgnP+kuRaf+l36V+bMeeY7Q+SPv3VbDlJisjl6Xk3AScVndPDtJiZWW58+cvMzHLjomJmZrlxUTEzs9y4qJiZWW5cVMzMLDcuKrbXk3SopDvTrx2SBsvW92t0vnKSTpD02oKO3Slpt6Q7y9oeL1teomQU4wVKRmDeJulzRWSxvdc+jQ5gVrSI+A+Sd1+QdBHweER8slF5JO0Tzw6aONYJwOPA/53A8doiYjTj7vdHxLEVjnEy8M/AqRGxDVgr6bck7yyYZeaeik1L6ZwRN0naKOn6siFMbpS0VtLNSuZyebWkvvRf8P8r3adTyZwf/6Zk3o8r0/GgxjvuxyXdBJwr6c/TQf7ukLRB0mFKBiT9b8B5aS/q9ZK+LOmtZbkfT7+foGQen68Dm5XMk/MJSbenmd49gT+L1wOXAKdHxP25/AHbtOWiYtORSP5V/taIeBXwReBjZdufjIg3AP9KMkTFe4CjgbdLOjTd50hgXUS8HHgM+HslY8DVOu7BEfHGiPgU8FNgUSSDZl4BrIyIB9Nzro2IYyPiJ+P8HN3Aqog4iuTt6Z0R8Wrg1cC70uE3xvO89Gd8c0Tck2F/s5p8+cumo+eRFIkbkuGNaCMZ3qLk6vT7ZmBrpEPGS/oVySB8vwMeioifpftdDrwf+P44x/1m2fI84JtpT2Y/kgm/Juq2SObCgGRyq5eX9WoOIhnPabzjPkVyqe2dwLmTyGD2HC4qNh2JpFgcX2X7E+n3p8uWS+ul35mx4xuVhg2vddxdZcv/DHw6Iq6WdAJwUZXP7CG9opAO8Ff+YEH58QS8LyKur3Kcap4G/hLYIOnCiPj4BD9v9hy+/GXT0RPAbEnHQzJ1gaQ/neAxFpQ+D7yN5HLWvRM47kHAYLp8Tln7CMm01CUPAq9Kl88E9q1yvOuBv0svwSHpJUomNxtXJKMBnwH8taR3ZvmMWTUuKjYdPQ28FVgjaYBkNNyJPsZ7N3COpLuA55NMKvbkBI57EfBtST8hGYq85HvAstKNepIb6G+UdBvwGp7bOyl3KckUAZskbQG+wASuRETEb4DTgA9JOjPr58zG8ijFZhOUPqV1TUQc3egsEzHR3JLeDnRFxHuLzGV7F/dUzKaPUeCg8pcfq5F0HnAByZNtZpm5p2JmZrlxT8XMzHLjomJmZrlxUTEzs9y4qJiZWW5cVMzMLDf/Hyqa+Tr+UShzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Start figure.\n",
    "fig = plt.figure()\n",
    "ax = fig.add_axes([0.15,0.15,0.80,0.80])\n",
    "\n",
    "# Plot.\n",
    "ax.plot(T_fl, F_fl, 'o', c=c[0], lw=1)\n",
    " \n",
    "# Add details and save figure.\n",
    "ax.set_xlabel(r'Temperature [K]')\n",
    "ax.set_ylabel(r'Free energy [eV/atom]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reversible Scaling thermodynamic integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "reversible_scaling_input_str = \"\"\"\\\n",
    "# This script runs a Reversible Scaling thermodynamic integration for the calculation of the temperature dependence of the free energy of bcc iron. See job.sh for the variables that need to be passed in the command line.\n",
    "\n",
    "\n",
    "#--------------------------- Simulation variables -----------------------------#\n",
    "  # Initalizes the random number generator.\n",
    "  variable         rnd equal round(random(0,999,{{ RANDOM }}))\n",
    "\n",
    "  # Simulation control parameters.\n",
    "  variable         t_eq equal 1000 # Equilibration time.\n",
    "  variable         t equal 10000 # Switching time.\n",
    "  variable         T equal 100 # Simulation temperature.\n",
    "#------------------------------------------------------------------------------#\n",
    "\n",
    "\n",
    "#---------------------------- Atomic setup ------------------------------------#\n",
    "  units            metal\n",
    "\n",
    "  # Create atoms.\n",
    "  read_data        structure.inp\n",
    "\n",
    "  # Define Interatomic potential.\n",
    "  include          potential.inp\n",
    "  neigh_modify     delay 0\n",
    "  timestep         0.002\n",
    "#------------------------------------------------------------------------------#\n",
    "\n",
    "\n",
    "#----------------------------- Run simulation ---------------------------------#\n",
    "  # Setup thermo output.\n",
    "  thermo_style     custom step pe\n",
    "  thermo           0\n",
    "\n",
    "  # Integrator for barostat and thermostat.\n",
    "  fix              f1 all nph iso 0 0 1\n",
    "  fix              f2 all langevin ${T} ${T} 0.1 ${rnd} zero yes\n",
    "\n",
    "  # Compute temperature using center-of-mass coordinates.\n",
    "  compute          c1 all temp/com\n",
    "  fix_modify       f1 temp c1\n",
    "  fix_modify       f2 temp c1\n",
    "\n",
    "  # Initial temperature to accelerate equilibration.\n",
    "  variable         rnd equal round(random(0,666,0)) # Generates new rnd #.\n",
    "  variable         T0 equal 2*${T}\n",
    "  velocity         all create ${T0} ${rnd} dist gaussian\n",
    "\n",
    "  # Forward integration.\n",
    "  run              ${t_eq}\n",
    "  print            \"$(pe/atoms) 1\" file forward.dat\n",
    "  variable         lambda equal 1/(1+(elapsed/${t})*({{ t_max }}/$T-1))\n",
    "  fix              f3 all adapt 1 pair eam scale * * v_lambda\n",
    "  fix              f4 all print 1 \"$(pe/atoms) ${lambda}\" screen no &\n",
    "                   append forward.dat title \"# pe lambda\"\n",
    "  run              ${t}\n",
    "  unfix            f3\n",
    "  unfix            f4\n",
    "\n",
    "  # Backward integration. \n",
    "  run              ${t_eq}\n",
    "  print            \"$(pe/atoms) $(v_T/{{ t_max }})\" file backward.dat\n",
    "  variable         lambda equal 1/(1+(1-(elapsed/${t}))*({{ t_max }}/$T-1))\n",
    "  fix              f3 all adapt 1 pair eam scale * * v_lambda\n",
    "  fix              f4 all print 1 \"$(pe/atoms) ${lambda}\" screen no &\n",
    "                   append backward.dat title \"# pe lambda\"\n",
    "  run              ${t}\n",
    "#------------------------------------------------------------------------------#\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "alat = 2.8841"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "tm_reversible_scaling = Template(reversible_scaling_input_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The job lmp_reversible_scaling was saved and received the ID: 7\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/pyiron/lammps/base.py:997: UserWarning: LAMMPS warning: No dump.out output file found.\n",
      "  warnings.warn(\"LAMMPS warning: No dump.out output file found.\")\n",
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/pyiron/lammps/base.py:583: UserWarning: LAMMPS warning: log.lammps does not contain the required pressure values.\n",
      "  warnings.warn(\"LAMMPS warning: log.lammps does not contain the required pressure values.\")\n"
     ]
    }
   ],
   "source": [
    "job_reversible_scaling = pr.create_job(pr.job_type.Lammps, \"lmp_reversible_scaling\")\n",
    "fe_bcc = pr.create_ase_bulk(\"Fe\", a=alat ,cubic=True).repeat(5)\n",
    "job_reversible_scaling.structure = fe_bcc\n",
    "job_reversible_scaling.potential = fe_potential\n",
    "job_reversible_scaling.input.control.load_string(\n",
    "    tm_reversible_scaling.render(\n",
    "        RANDOM=np.random.randint(999),\n",
    "        t_max=temperature_lst[-1]\n",
    "    )\n",
    ")\n",
    "job_reversible_scaling.run()\n",
    "job_reversible_scaling.decompress()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "T0 = 100 # Reference temperature [K]\n",
    "kB = sc.value('Boltzmann constant in eV/K') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "F0 = F_fl[T==T0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "U_f, lamb_f = np.loadtxt(os.path.join(job_reversible_scaling.working_directory, \"forward.dat\"), unpack=True)\n",
    "U_b, lamb_b = np.loadtxt(os.path.join(job_reversible_scaling.working_directory, \"backward.dat\"), unpack=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fix adapt also scales the potential energy besides the forces, so we unscale.\n",
    "U_f /= lamb_f\n",
    "U_b /= lamb_b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute work done using cummulative integrals [Eq.(21) in the paper].\n",
    "I_f = cumtrapz(U_f,lamb_f,initial=0)\n",
    "I_b = cumtrapz(U_b[::-1],lamb_b[::-1],initial=0)\n",
    "W = (I_f+I_b) / (2*lamb_f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute free energy [Eq.(22) in the paper] and save results.\n",
    "T_rs = T0 / lamb_f\n",
    "F_rs = F0/lamb_f  + 1.5*kB*T_rs*np.log(lamb_f) + W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fdca4ec9810>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAETCAYAAAD+spv+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5OklEQVR4nO3de3zP9f//8dvDnEZMOaRZTEQfhpmFGaExZzJ8JJRD+EQ+PqmPUycdP6k++VB+fT5J8pHIuZKc5VTCnKWSLJn6OFRyLJvH74/3e/vO2uE9e7/3em97XC+X12Xv9+t439j7sefr8HyKqmKMMcZ4QxGnAxhjjCk4rKgYY4zxGisqxhhjvMaKijHGGK+xomKMMcZrijodwBcqVKigoaGhTscwxpgCKT4+/pSqVsxoWYEsKqGhoezYscPpGMYYUyCJyHeZLbPTX8YYY7zGiooxxhivsaJijDHGa6yoGGOM8RorKsYYY7zGiooxxhivsaJijDHGa6youJ1fvIQfGzclMaQqPzZuyvnFS5yOZIwx+Y4VFVwF5cyYsSQnJoIqyYmJnHpoNKtGP8yqVas4ePAgly5dcjqmMYVaQEAA4eHhhIWF0aVLF3755Zc8z3D8+HF69uwJwNtvv82DDz6Y4XrXXXddjvb71ltvUa9ePerXr09YWBjvv//+NeVLOW7anHmtQD5Rn1NnX5iEXrx41byiSUkEL32fUXt2c/ToUY4ePcqNN95IjRo1qFmzJmFhYdSvX5969epRoUIFh5IbU3gEBgaye/duAO677z6mTZvGo48+6pNjJSUlUbToHz8eg4ODWbhwoVePdezYMZ577jl27txJUFAQ586d4+TJk7napy9yespaKkDy8eMZzr/+999Zu3Ythw4d4vz582zYsIEJEybQsGFDvvrqK5544glq1KhBcHAw7du357HHHmPZsmWcOnUqj78DYwqXqKgoEhMTATh8+DDt27enUaNGtGjRgi+//JIzZ84QGhrKlStXALhw4QI333wzly9fznB9gAEDBjB69Ghat27N2LFj2bBhA+Hh4YSHh9OwYUPOnj1LQkICYWFhqTm+//572rdvT+3atXnqqacyzPrSSy9x++23U79+fZ588sk/LD9x4gRlypRJbWVcd911VK9eHYBvvvmGNm3a0KBBAyIiIjh8+DDnzp0jJiaGiIgI6tWrl2GrJm3Ot99+m7i4ONq3b8+tt97KmDFjUtebMWMGtWrVolWrVgwZMiTTlleOqGqBmxo1aqQ58cPtTfRYcMgfph9ub5LttleuXNGEhAR9//339bHHHtOYmBgtU6aM1qxZU/v376/Tp0/XhISEHOUxxvxR6dKlVVU1KSlJe/bsqR9//LGqqt5555369ddfq6rq1q1btXXr1qqq2rVrV123bp2qqs6bN08HDx6c5fr33XefdurUSZOSklRVtXPnzrp582ZVVT179qxevnxZjxw5onXr1lVV1ZkzZ2rlypX11KlTeuHCBa1bt65u3779qqwrV67UIUOG6JUrVzQ5OVk7deqkGzZsuOr7SkpK0tjYWL355pt1wIAB+sEHH6Qua9y4sS5evFhVVS9evKjnz5/Xy5cv65kzZ1RV9eTJk1qjRg29cuXKVcdNn7N69er6yy+/6MWLF7Vq1ap69OhRTUxM1GrVqunp06f1999/1+bNm+uIESM8+rcAdmgmn792+gsoM24sZ8aMveoUmAQGUmbc2Gy3FRGqVatGtWrV6Nq1KwDJyckcPHiQzz77jHXr1jFhwgTKlStH27Ztadu2La1btyYoKMhn348xviYiXt+n67MqcxcvXiQ8PJyEhAQaNWpE27ZtOXfuHJ9++im9evVKXe+3334DoHfv3rz33nu0bt2aefPmMXz48CzXB+jVqxcBAQEAREdHM3r0aPr27UtcXBwhISF/yNS2bVvKly8PQFxcHJs3byYyMjJ1+apVq1i1ahUNGzYE4Ny5cxw6dIg77rgjdZ2AgABWrFjB9u3bWbt2LQ899BDx8fE8/PDDJCYm0r17dwBKliwJwOXLl5kwYQIbN26kSJEiJCYm8r///Y/KlStn+rOLiYlJ/cypU6cO3333HadOnaJly5bccMMNqd/7119/neW/gSesqACl41z/aGdfmETy8eMEBAdTZtzY1Pk5FRAQQFhYGGFhYQwZMoQrV66wd+9eVq9ezf/7f/+P/v37Ex4eTrdu3ejatSu1atXy5rdjjM9lVwB8IeWaypkzZ+jcuTPTpk1jwIABlCtXLvVaS1pdu3Zl/Pjx/PTTT8THx3PnnXdy/vz5TNcHKF26dOrrcePG0alTJ5YvX07Tpk1Zs2ZN6gd7ivTFNf17VWX8+PEMGzYsy+9NRGjcuDGNGzembdu2DBw4kNGjR2e47pw5czh58iTx8fEUK1aM0NDQbG8kKlGiROrrgIAAkpKSfPZvaNdU3ErHdafytq1UOXaUytu2XnNByUiRIkUIDw/n73//O6tWreLEiRNMmDCBb775hlatWnHbbbcxduxYtmzZQnJysteOa0xBFBQUxNSpU3n55ZcJDAykevXqLFiwAHB9iO/ZswdwXZto3Lgxo0aNonPnzgQEBFC2bNlM10/v8OHD1KtXj7FjxxIZGZl67SWt1atX89NPP3Hx4kWWLl1KdHT0VcvbtWvHW2+9xblz5wBITEzkxIkTV61z/Phxdu7cmfp+9+7dVKtWjbJlyxISEsLSpUsBV4vqwoULnDlzhkqVKlGsWDHWr1/Pd99l2gt9lho3bsyGDRv4+eefSUpKYtGiRde0n/QcLSoi8oiIqIj84fYpESkpIttEZI+IHBCRjK+C5UOBgYF06NCBf//73xw7dozZs2dTvHhx/vKXv1C1alVGjx7N9u3bHflr0Jj8oGHDhjRo0IB58+YxZ84cZsyYQYMGDahbt+5VF6579+7NO++8Q+/evVPnZbV+Wv/6178ICwujQYMGqb+z6TVv3jz1zEOPHj2uOvUFEBsbyz333ENUVBT16tWjZ8+enD179qp1Ll++zCOPPMJtt91GeHg47733HlOmTAFg9uzZTJ06lfr169OsWTN+/PFH+vbty44dO4iMjGTOnDncdttt1/QzrFKlChMmTKBJkya0adOGOnXqeOW0vDj1wSUiNwNvArcBjVT1VLrlApRW1XMiUgzYDIxS1a3Z7TsyMlLz6yBdX3zxBfPmzWPu3LmoKnfffTd9+vShbt26TkczxhQw586d47rrriMpKYnu3bszaNCg1Gs4WRGReFWNzGiZky2VycAYIMOq5r7J4Jz7bTH3VOD/dK9Tpw5PP/00X3/9Ne+99x6XLl2iXbt21K9fn3/84x98//33Tkc0xhQQEydOTH2gtHr16tx111253qcjLRUR6QrEqOooEUkAItO3VNzrBQDxQE1gmqpmejuWiAwFhgJUrVq10bWeZ/RHV65cYcuWLbz77rvMnz+fJk2acP/999OlSxeKFSvmdDxjTCGTVUvFZ0VFRNYAGd3j9igwAYhV1TNZFZU0+yoHLAFGqur+7I6dn09/ZefixYssWrSI6dOn89VXX3Hvvfdy//332x1kxpg848jpL1Vto6ph6SfgW6A6sMddUEKAnSKS6U3WqvoL8AnQ3ld584vAwED69evHhg0b2LhxIyLCHXfcQcuWLZk9e7b1UWaMcVSeX1NR1X2qWklVQ1U1FDgGRKjqj2nXE5GK7hYKIhIItAH+eE9fIVarVi0mTZrE0aNHGTVqFHPmzKFq1ao8+uijHDt2zOl4xphCyK+eUxGRYBFZ7n57E7BeRPYC24HVqrrMuXT+q3jx4sTFxbFixQo2b97MuXPnqF+/Pr1792bLli12a7IxJs84XlTcLZZT7tfHVbWj+/VeVW2oqvXdp86edjZp/lCrVi2mTJlCQkICzZs3Z+DAgURGRjJr1iw7NWbytZSu71OmhISEXO+zVatWeHr9NX1nkmkNGDDgmnsFnjhxIi+//HKOjufPHC8qxjfKli3LyJEj+fLLL3n22WeZN28e1apV4/HHH+d4Jr0yG+Mtvhj0LqWblpQpNDQ0dZmqpvZIbJxlRaWAK1KkCB06dODjjz9m48aN/Pzzz6l9kqV0HmejXhpvymjQuzNjxnr9/1VCQgJ/+tOfGD58OBEREXz//fcZdjOfst6QIUOoW7cusbGxXEw3ftKVK1e47777eOyxx0hOTubvf/976n7+85//XFO+rLqof+6556hduzZt2rThq6++Sp0fHx9PgwYNiIqKYtq0add0XKdZUSlEateuzWuvvcahQ4cICQmhefPm/KNlK3565O8+/wAwhUdGg97pxYucfWFSrvab0ktxeHh46lPfKbfV79q1i6+++opDhw6xbds2du/eTXx8PBs3bgTg0KFDjBgxggMHDlCuXLmr+rlKSkqib9++1KpVi2effZYZM2YQFBTE9u3b2b59O9OnT+fIkSM5zluyZEmWLFnCzp07Wb9+PQ8//DCqSnx8PPPmzWPXrl0sXryY7du3p24zcOBApk6dymeffZarn5WTrKgUQuXLl+fJJ5/kyJEj9DhxkiJpuv4G73wAmMIrs0HvMpvvqbSnv5Yscf3RU61aNZo2bQpc3c18REQEX375JYcOHQKgevXqhIeHA9CoUaOrrscMGzaMsLCw1FEkV61axX//+1/Cw8Np0qQJp0+fTt1PTqgqEyZMoH79+rRp0ya1i/pNmzbRvXt3SpUqRdmyZVOHzDhz5gy//PILLVu2BKB///7X9HNymnV9X4iVLl2a0uk6t0uR2w8AU3gFBAe7Wr4ZzPe2tF3VZ9bNfEJCwh+6fk97+qtZs2apLYmSJUuiqrz66qu0a9fuD/tJMXDgQHbt2kVwcDDLly8nI1l1UZ/ReDSq6pNxavKatVQKucx+0U+IsHTpUrv4aXKszLixSGDgVfM8HfQuNzzpZj4jgwcPpmPHjvTq1YukpCTatWvH66+/zuXLlwH4+uuvOX/+/FXbzJw5k927d2daUIBMu6i/4447WLJkCRcvXuTs2bN8+OGHAJQrV46goCA2b94MuIpSfmRFpZDL6AOAwEDO3dufZ555hoYNG7JgwQIrLsZjpeO6E/TiJAKqVAERAqpUIejFSV4doygjnnQzn5nRo0cTERFB//79uf/++6lTpw4RERGEhYUxbNgwkpKSst3HsGHDCAkJISQkhKioqEy7qI+IiKB3796p3eW3aNEidR8zZ85kxIgRREVFEZj+9zKfcKzre18qyH1/+cL5xUsyHPVSVVm+fDnPPPMMv/76K4899hi9e/dOHW7VGFM4OdKhpJOsqHiXqrJ69WqefvppTpw4wcSJE7n77rspUsQausYURv46norJJ0SE2NhYNm3axOuvv87UqVNp2LAhH374oXUBY4y5ihUV4zERISYmhs8++4ynn36a8ePHEx0dzSeffOJ0NGOMn7CiYnJMROjWrRt79uxh+PDhDBo0iHbt2nnch5IxpuCyomKuWUBAAP369ePLL7/krrvuolu3bvTo0YMvvvjC6WjGGIdYUTG5Vrx4cR544AEOHTpEkyZNaNmyJQMHDiQxgwfgjDEFmxUV4zWlSpVizJgxfPPNN9x0003Ur1+fxx9/3ONnBYwx+Z8VFeN1QUFBPP/88+zatYuEhARq1arFG2+84dEDZMaY/M2KivGZqlWrMnv2bJYtW8a7775LeHg4H3/8sd2GbEwBZkXF+FyjRo1Yv349//jHP3jooYeIjY1lz549TscyxviAo0VFRB4RERWRClmsEyAiu0TExqfPx0SELl26sG/fPrp37067du144IEHOHXqlNPRjDFe5FhREZGbgbbA0WxWHQUc9H0ikxeKFSvG8OHDOXjwIMWLF6dOnTq89tprdr3FmALCyZbKZGAMkOkJdhEJAToBb+ZVKJM3rr/+eqZMmcK6detYsmQJDRs2ZP369U7HMsbkkiNFRUS6Aomqmt2J9X/hKjzZ9rsuIkNFZIeI7Dh58qQXUpq8EBYWxpo1a3jqqacYNGgQvXr1Sh13whiT//isqIjIGhHZn8HUDXgUeCKb7TsDJ1Q13pPjqeobqhqpqpEVK1b0wndg8oqIEBcXxxdffEH9+vWJiIjgySef5MKFC05HM8bkkM+Kiqq2UdWw9BPwLVAd2CMiCUAIsFNEKqfbRTTQ1b3OPOBOEXnHV3mN8wIDA3n88cfZtWsXX375JXXq1OH999+3W5CNyUccH0/FXTQiVTXT24BEpBXwiKp29mSfNp5KwbB27VoefPBBatSowdSpU7nlllucjmSMIR+NpyIiwSKS+aDPplCJiYlhz549NG/enMaNG/PMM89w6dIlp2MZY7LgeFFR1dCUVoqqHlfVjhms84mnrRRTsBQvXpxx48YRHx/Pzp07qVevHitXrnQ6ljEmE44XFWM8Ua1aNZYsWcKUKVMYPnw4vXr14tixY07HMsakY0XF5CsdO3Zk//791K1bl/DwcF566SUuX77sdCxjjJsVFZPvBAYGMnHiRLZu3cq6desIDw9n06ZNTscyxmBFxeRjNWvWZPny5Tz99NP06dOHYcOG8csvvzgdy5hCzYqKyddEhB49enDgwAGKFi1K3bp1WbBggT3bYoxDrKiYAiEoKIhp06axYMECnnrqKbp27cr333/vdCxjCh0rKqZAadasGTt37qRJkyZEREQwdepUkpOTnY5lTKFhRcUUOMWLF+exxx5j8+bNLF68mKioKBsUzJg8UjSzBSIS58H2l1TVnoA3fql27dqsW7eOmTNn0rZtWwYNGsSTTz5JYGCg09GMKbCyaqlMBzoDXbKYXvV1QGNyo0iRIgwePJi9e/eSkJBAvXr1WLNmjdOxjCmwMm2pAB+r6qCsNrZeg01+UblyZebNm8dHH33E4MGDadWqFa+88golN2zk7AuTSD5+nIDgYMqMG0vpuO5OxzUm38q0paKq/bLb2JN1jPEnnTp14sCBA9xwww08Uq8+p0c/THJiIqiSnJjImTFjOb94idMxjcm3su36XkQCcA3pG0qalo2qvuLTZLlgXd8bT3zXoCFFT/1xxIWAKlWovG2rA4mMyR+y6vo+q9NfKT4ELgH78GBYX2Pyi6KnT2c4P/n48TxOYkzB4UlRCVHV+j5PYkweCwgOdp36Su/GG/M+jDEFhCfPqXwsIrE+T2JMHiszbiyS7vbiywEBTEg8xuzZs62rF2OugSdFZSuwREQuisivInJWRH71dTBjfK10XHeCXpxEQJUqIEJAlSpU/NdkRi3/iJdffpkuXbpw3E6FGZMjnhSVfwJRQClVLauqZVS1rDcOLiKPiIiKSIVMlieIyD4R2S0iduXdeF3puO5U3raVKseOUnnbVkrHdSciIoLt27cTGRlJw4YNmTt3rrVajPGQJ0XlELBfvfxbJSI3A22Bo9ms2lpVwzO708AYXyhevDgTJ07ko48+4plnnqF3796cyuBOMWPM1TwpKj8An4jIeBEZnTJ54diTgTGA/Qlo/FZkZCTx8fHcfPPN1K9fn2XLljkdyRi/5klROQKsBYoDZdJM10xEugKJqppdL38KrBKReBEZms0+h4rIDhHZcfLkydzEM+YqgYGB/POf/2Tu3Ln89a9/ZfDgwfz6q11WNCYj2T78mLqiSBlAVfWch+uvASpnsOhRYAIQq6pnRCQBiFTVP5xbEJFgVT0uIpWA1cBIVd2Y3bHt4UfjK2fPnmX06NGsXr2amTNn0rp1a6cjGZPnsnr4MduWioiEicguYD9wwN1qqJvddqraRlXD0k/At0B1YI+7oIQAO0XkDwVIVY+7v54AlgCNszuuMb5UpkwZpk+fzrRp0+jXrx8PPfQQly5dcjqWMX7Dk9NfbwCjVbWaqlYDHsbVg/E1UdV9qlpJVUNVNRQ4BkSo6o9p1xOR0u7WESJSGojFVdiMcVynTp3Yu3cvx44dIzIykr179zodyRi/4ElRKa2q61PeqOonQGlfhBGRYBFJGZ/lRmCziOwBtgEfqeoKXxzXmGtRvnx55s+fzyOPPEJMTAyvvPIKV65YT0amcPOkQ8klwE5gtntWP1zXQO7ybbRrZ9dUTF779ttv6devH6VKlWLWrFlUqVLF6UjG+EyurqkAg4CKwGL3VAEY4LV0xhQAt9xyCxs3bqRly5ZERESwcOFCpyMZ4whPWiq9VHVBdvP8ibVUjJM+//xz+vXrR/PmzZkyZQply3qlAwpj/EZuWyrjPZxnjAGaNGnCrl27KFq0KA0bNuTTTz91OpIxeSbTru9FpAPQEagiIlPTLCoLJPk6mDH52XXXXcf06dNZunQpcXFxDB06lMcff5xixYo5Hc0Yn8qqpXIc2IFrgK74NNMHQDvfRzMm/7vrrrvYtWsX27Zto0WLFnzzzTdORzLGp7Iao36Pqs4CaqrqrDTTYlX9OQ8zGpOv3XTTTSxfvpy+ffsSFRXFm2++ab0emwLLk2sqoSKyUES+EJFvUyafJzOmAClSpAgjR47kk08+4bXXXiMuLs56PTYFkidFZSbwOq7rKK2B//J/z6wYY3Kgbt26fP7559SsWZMGDRqwevVqpyMZ41WeFJVAVV2L6/bj71R1InCnb2MZU3CVKFGCl156if/+978MGDCAsWPHcvnyZadjGeMVnhSVSyJSBDgkIg+KSHegko9zGVPgxcTEsHv3bvbv30/z5s359ls7q2zyP0+Kyt+AUsBfgUa4umm5z4eZjCk0KlasyLJly+jTpw9NmjRh3rx5TkcyJlcyfaJeRMYDK1R1V95Gyj17ot7kRzt37uTuu++mefPmvPrqq5Qu7ZN+W43JtWt9ov4IMEpEdonI2yLSW0Su901EY0xERATx8fFcuXKFRo0asXv3bqcjGZNjWT2nMk9VB6hqQ2AKcAuwWEQ2isgTImIDZhnjZWXKlOHtt9/m8ccfp23btkydOtWeaTH5SqZFRURS+5NQ1V2q+g9VbQ10Bg4A9+dBPmMKpb59+/LZZ58xe/ZsunXrZs+0mHwjq9NfiSIyXUTuFBFJmamqv6rqIlUdmgf5jCm0atasyZYtW6hduzbh4eGsX78++42McVhWReVPuPr+ehz4XkT+JSJN8iaWMQagePHivPTSS7z55pvcc889PP744yQlWX+uxn9ldU3ltKr+x33KqzGuC/f/EpHDIvKcNw4uIo+IiIpIhUyWl3N3EfOliBwUkShvHNeY/KZ9+/bs2rWLrVu3EhMTw/Hjx52OZEyGPHlOBVU9DszA1V3LWbxwPUVEbgbaAkezWG0KrtuabwMaAAdze1xj8qvKlSuzYsUKYmJiaNSoEWvWrHE6kjF/kGVREZGSItJLRBYDh4EYXAN0BXvh2JOBMUCGt7aISFngDlzFDFX9XVV/8cJxjcm3AgICeOKJJ3jnnXe49957mThxIsnJyU7HMiZVVnd/vYurFdEbeBeopqr3qerHqpqr/8Ui0hVIVNU9Wax2C3ASmOl+VuZNEcn0aTARGSoiO0Rkx8mTJ3MTzxi/FxMTQ3x8PBs2bKBdu3b873//czqSMUDWLZWVQA1V7amqC1X1Uk52LCJrRGR/BlM34FHgiWx2URSIAF53PytzHhiX2cqq+oaqRqpqZMWKFXMS1Zh86aabbmL16tVERUURERHBJ5984nQkY7K8UD9LVc+KyI0iMkNEVgCISB0RGZzdjlW1jaqGpZ+Ab4HqwB4RSQBCgJ0iUjndLo4Bx1T1c/f7hbiKjDHGrWjRojzzzDO89dZb3H333Tz//PNcuXLF6VimEPPkQv3buFotN7nff42rk8lroqr7VLWSqoaqaiiu4hGhqj+mW+9HXLcy13bPigG+uNbjGlOQtWvXjh07drB8+XI6depkD0sax3hSVCqo6nzgCoCqJgE+uTIoIsEisjzNrJHAHBHZC4QDz/viuMYUBCEhIaxfv5769esTERHB1q1bnY5kCiFPisp5ESmP+y4tEWkKnPFWAHeL5ZT79XFV7Zhm2W73dZL6qnqXqv7sreMaUxAVK1aMSZMm8eqrr9K1a1emTZtmfYeZPFXUg3VGAx8ANURkC1AR6OnTVMaYXOnWrRt169YlLi6OrVu38u9//9u60jd5ItuWiqruBFoCzYBhQF1V3evrYMaY3KlZsyZbt25FRIiKiuLQoUNORzKFQFbPqaTeaaWqSap6QFX3q+rljNYxxvifUqVKMWvWLB544AGio6N5//33nY5kCrisWiozReR6Ebkhswn30+7GGP8lIjzwwAN8+OGHjBw5kvHjx5OUlMT5xUv4sXFTEkOq8mPjppxfvMTpqKYAyOqaShAQD0gW69ij68bkE02aNCE+Pp4+ffowsWkUD/z8C1xyPdOcnJjImTFjASgd193BlCa/y3SM+vzMxqg3JnPJycl8/ae6lD1//g/LAqpUofI2uxXZZO1ax6g3xhRAAQEBlL1wIcNlydalvsklKyrGFEIBwRl3NJ7ZfGM8ZUXFmEKozLixSGDgVfMuAZeH5HqoJFPIZVtURGSRiHQSEStAxhQQpeO6E/TiJAKqVAERAqpUYWeH9jR+4nHr7djkSrYX6kWkDTAQaAosAN5W1S/zINs1swv1xlybNWvW0K9fP8aNG8eoUaMQyermT1NY5epCvaquUdW+uLqdTwBWi8inIjJQRIp5N6oxxklt2rRh69atzJo1i/79+3Mhkwv6xmTGo1Na7g4lB+Aam34XrrHjI4DVPktmjHFEaGgoW7ZsASA6OpqEhARnA5l8xZNrKouBTUApoIuqdlXV91R1JHCdrwMaY/JeqVKlmD17NgMGDKBp06asWbPG6Ugmn/CkpfKaqtZR1X+o6g9pF2R2Ts0Yk/+JCKNGjWLevHn079+fl19+2brRN9nypOv7ciISl27eGWCfqp7wQSZjjB9p1aoVn3/+OXFxcezYsYMZM2ZYN/omU560VAYDbwJ93dN0XGOsbBGR/j7MZozxE1WrVmXTpk2ULFmSqKgovv32W6cjGT/lSVG5AvxJVXuoag+gDvAb0AQY68twxhj/ERgYyMyZMxk2bBhRUVGsXLnS6UjGD3lSVEJV9X9p3p8AaqnqT8DlTLbxiIg8IiIqIhUyWFZbRHanmX4Vkb/l5njGmNwREUaMGMHChQsZOHAgL7zwgl1nMVfxpKhsEpFlInKfiNwHvA9sFJHSwC/XemARuRloCxzNaLmqfqWq4aoaDjQCLgA24IMxfqBFixZs27aNJUuW8Oc//5lz5845Hcn4CU8efhwOzATCgYbAf4ERqnpeVVvn4tiTgTGAJ3/mxACHVfW7XBzPGONFISEhbNy4kaCgIJo2bWrDFRsgm6IiIkVEZL+qLlLVh1T1b6q6UHPZ3hWRrkCiqu7xcJO7gbnZ7HOoiOwQkR0nT9rYYcbkhRIlSjB9+nRGjhxJdHQ0H330kdORjMOyLCqqegXYIyJVc7pjEVkjIvszmLoBjwJPeLif4kBXXP2OZZX1DVWNVNXIihUr5jSuMeYaiQjDhg1j6dKlDB06lGeffZYrV644Hcs4xJPnVG4CDojINiB1qDhV7ZrVRqraJqP5IlIPqI6rWAGEADtFpLGq/pjBJh2AneluFjDG+JlmzZqxfft2evbsSXx8PLNmzaJs2bJOxzJ5zJOi8pQ3D6iq+4BKKe9FJAGIVNVTmWzSh2xOfRlj/ENwcDCffPIJo0aNokmTJixdupTatWs7HcvkIU8u1G/A1TtxMffr7cBOX4QRkWARWZ7mfSlcd4gt9sXxjDHeV7x4cV5//XUefvhhWrRowQcffOB0JJOHsm2piMgQYChwA1ADqAL8G9cdWbmmqqFpXh8HOqZ5fwEo743jGGPy1v3330+9evXo0aMHBw4cYNy4cTY+SyHgyXMqI4Bo4FcAVT1EmtNXxhiTmSZNmvD555+zePFi+vfvz6VLl5yOZHzMk6Lym6r+nvJGRIri2bMlxhhDlSpV2LBhA0lJSbRs2ZIffvgh+41MvuVJUdkgIhOAQBFpi+vW3g99G8sYU5CUKlWKuXPn0qVLFxo3bkx8fLzTkYyPeFJUxgEngX3AMGA58JgvQxljCh4R4bHHHmPq1Km0b9+e9957z+lIxgeyvVDvfgByunsyxphc6d69O9WrV+euu+7iwIEDTJw4kSJFPBrZ3OQDngwnHC0iq0XkaxH5VkSOiIgNpmCMuWbh4eFs27aNdevW0atXL86fP5/9RiZf8OTPgxnAK0Bz4HYg0v3VGGOuWaVKlVi7di1BQUFER0dz9GiGHZabfMaTonJGVT9W1ROqejpl8nkyY0yBV6JECWbMmMG9995L06ZN+fTTT52OZHLJk6KyXkReEpEoEYlImXyezBhTKIgIo0eP5s033+Suu+7i7bffdjqSyQVP+v5q4v4amWaeAnd6P44xprDq2LEjn3zyCV27dmX//v1MmjSJgIAAp2OZHJKCOBRoZGSk7tixw+kYxphrcPr0af785z9TsmRJ3n33XYKCgpyOZNIRkXhVjcxomSd3f90oIjNE5GP3+zoiMtjbIY0xBqB8+fKsWLGC0NBQoqKi+Oabb5yOZHLAk2sqbwMrgWD3+6+Bv/kojzHGUKxYMaZNm8aDDz5IdHQ069evdzqS8ZAnRaWCqs4HrgCoahKQ7NNUxhgDDB8+nLlz59KnTx9ef/11p+MYD3hSVM6LSHncnUiKSFPgjE9TGWOM25133snmzZt59dVXGTFiBElJSU5HMlnwpKiMBj4AaojIFuC/wEifpjLGmDRq1qzJZ599xuHDh+ncuTNnztjftf7Kk5EfdwItgWa4OpSsq6p7fR3MGGPSCgoKYtmyZdSoUYPo6GgSEhKcjmQy4FEvbqqapKoHVHW/ql721sFF5BERURGpkMnyh0TkgIjsF5G5IlLSW8c2xuQ/RYsW5bXXXmPIkCE0a9aMzz//3OlIJh3HugYVkZtxjT+fYYc/IlIF+CsQqaphQABwd94lNMb4IxFh1KhR/Oc//6Fz587Mnz/f6UgmDSf7m54MjCHrUSSL4hocrChQCjieF8GMMf6vS5curF69mkceeYTnn3+egvggd37kycOPIiL9ROQJ9/uqItI4NwcVka5AoqruyWwdVU0EXsbVkvkBV8eWq7LY51AR2SEiO06ePJmbeMaYfCI8PJytW7eyaNEiBg4cyO+//579RsanPGmp/D8gCujjfn8WmJbdRiKyxn0tJP3UDXgUeCKb7a8HugHVcT14WVpE+mW2vqq+oaqRqhpZsWJFD74tY0xBEBwczMaNGzlz5gyxsbGcPm2dqDvJk6LSRFVHAJcAVPVnoHh2G6lqG1UNSz8B3+IqFHtEJAEIAXaKSOV0u2gDHFHVk+6bAxbjugPNGGOuUrp0aRYtWkTjxo2Jiori0KFDTkcqtDzppfiyiATwfw8/VsT9dP21UNV9QKWU9+7CEqmqp9KtehRoKiKlgItADGC9RBpjMlSkSBFefPFFbr31Vpo3b878+fNp2bKl07EKHU9aKlOBJcCNIvIcsBl43hdhRCRYRJYDqOrnwEJgJ7DPnfUNXxzXGFNwDBkyhDlz5tCrVy9mzZrldJxCx6Ou70XkNlwtBQHWqupBXwfLDev63hhz8OBBOnfuTJ8+fXj66acpUsTJm10Lllx1fe9WAbigqq8Bp0SkutfSGWOMD/zpT39i69atrFu3jn79+vHbb785HalQ8OSW4ieBscB496xiwDu+DGWMMd5QsWJF1q5dy++//05sbCw//fST05EKPE9aKt2BrsB5AFU9DpTxZShjjPGWwMBA5s+fz+23306zZs04cuSI05EKNE/u/vpdVVVEUu7+Ku3jTMYY41VFihTh5ZdfJjQ0lOjoaFaNfpgK8xeQfPw4AcHBlBk3ltJx3Z2OWSB4UlTmi8h/gHIiMgQYBEz3bSxjjPG+Bx98kAY//kjJVyaTLAJAcmIiZ8aMBbDC4gVZnv4SEQHew3Vr7yKgNvCEqr6aB9mMMcbrbl21hkB3QUmhFy9y9oVJDiUqWLJsqbhPey1V1UbA6jzKZIwxPpN8PON+aTObb3LGkwv1W0Xkdp8nMcaYPBAQHJzh/CI33ZTHSQomT4pKa1yF5bCI7BWRfSJiIz8aY/KlMuPGIoGBV837vUgRJv92CevhPPcyPf0lIlVV9SjQIQ/zGGOMT6VcjD/7wqTUu78qjh1D0fgdNGvWjOXLl3Prrbc6nDL/yuqaylIgQlW/E5FFqtojjzIZY4xPlY7r/oc7vZ7vEUdoaCgtWrRg8eLFNGtmnaJfi6xOf6W9PeIWXwcxxhinDR06lJkzZ9KtWzcWLlzodJx8KauWimby2hhjCqwOHTqwatUqunTpwnfffcfo0aORdLcgm8xlVVQaiMivuFosge7XuN+rqpb1eTpjjHFAw4YN+fTTT+nYsSNHjhxhypQpBAQEOB0rX8j09JeqBqhqWVUto6pF3a9T3ltBMcYUaFWrVmXLli0cPHiQnj17cvHiRacj5Qs2wIAxxmQiKCiI5cuXExgYSJs2bTh9+rTTkfyeFRVjjMlCiRIleOedd4iOjqZ58+Z89913Tkfya44WFRF5RERURCpksnyUiOwXkQMi8rc8jmeMMYCrl+MXX3yRv/zlL0RHR7N7926nI/ktx4qKiNwMtAWOZrI8DBgCNAYaAJ1FxJ5IMsY4ZtSoUUyePJnY2FjWrl3rdBy/5GRLZTIwhsxvV/4TsFVVL6hqErAB14BhxhjjmF69erFgwQLuuece5syZ43Qcv+NIURGRrkCiqu7JYrX9wB0iUl5ESgEdgZuz2OdQEdkhIjus/x5jjC+1bNmStWvXMn78eF588UVU7VG+FJ4M0nVNRGQNUDmDRY8CE4DYrLZX1YMiMglXl/vngD1AUhbrvwG8ARAZGWn/wsYYnwoLC+PTTz+lQ4cOHDt2jMmTJ9uzLPiwpaKqbVQ1LP0EfAtUB/aISAIQAuwUkT8UIFWdoaoRqnoH8BNwyFd5jTEmp0JCQti0aRP79u2jd+/eXLp0yelIjsvz01+quk9VK6lqqKqGAsdwdVz5Y/p1RaSS+2tVIA6Ym6dhjTEmG+XKlWPFihUEBAQQGxvLzz//7HQkR/nVcyoiEiwiy9PMWiQiXwAfAiNUtXD/axlj/FKJEiWYO3cukZGRNG/enO+//97pSI5xvKi4Wyyn3K+Pq2rHNMtaqGodVW2gqnb/njHGbxUpUoRXXnmFQYMGER0dzf79+52O5AifXag3xpjC6OGHH6Zy5crExMSwYMEC7rjjDqcj5SnHWyrGGFPQ9O3blzlz5tCzZ08WLVrkdJw8ZS0VY4zxgTZt2rBy5Uo6d+7Mjz/+yIgRI5yOlCesqBhjjI80bNiQTZs20b59e44fP86zzz5b4Af8stNfxhjjQ7fccgtbtmxh9erVDBo0iMuXLzsdyaesqBhjjI9VrFiR9evXc+LECbp168b58+edjuQzVlSMMSYPlC5dmqVLl3LjjTfSunVrCmofhVZUjDEmjxQrVoy33nqL2NhYoqOj+fbbb52O5HV2od4YY/KQiPDss88SHBxM8+bNWbZsGREREU7H8hprqRhjjAOGDx/Oq6++Svv27Vm9erXTcbzGiooxxjikR48eLFq0iH79+jF3bsHoL9dOfxljjINatGjB2rVrad++PadOnWLkyJFOR8oVKyrGGOOwsLAwNm/eTGxsLCdPnuSpp57Ktw9J2ukvY4zxA6GhoWzevJmPPvqI4cOHk5yc7HSka2JFxRhj/ESlSpVYv349X331FX369OG3335zOlKOWVExxhg/UrZsWZYvX05ycjKdO3fm7NmzTkfKESsqxhjjZ0qWLMn8+fMJDQ0lJiYmXz1970hREZGJIpIoIrvdU8dM1msvIl+JyDciMi6vcxpjjFMCAgJ44403aNu2LS1atODo0aNOR/KIk3d/TVbVlzNbKCIBwDSgLXAM2C4iH6jqF3kV0BhjnCQiPPfcc1SsWJHmzZuzYsUK6tSp43SsLPnzLcWNgW9U9VsAEZkHdAOsqBhjCpW//e1vlC9fnjvvvJOlS5fStGlTpyNlyslrKg+KyF4ReUtErs9geRXg+zTvj7nnGWNModO/f39mzJhBly5dWLlypdNxMuWzoiIia0RkfwZTN+B1oAYQDvwA/DOjXWQwT7M43lAR2SEiO/LTRS1jjPFUp06dWLp0Kffee6/fduvis9NfqtrGk/VEZDqwLINFx4Cb07wPAY5ncbw3gDcAIiMjMy0+xhiTn0VHR7NmzRo6dOjAzz//zPDhw52OdBVHrqmIyE2q+oP7bXdgfwarbQduFZHqQCJwN3BPHkU0xhi/Va9ePTZu3Ejbtm05c+YM48aN85tuXZy6UP+iiITjOp2VAAwDEJFg4E1V7aiqSSLyILASCADeUtUDDuU1xhi/csstt7Bp0yZiY2P5+eefmTRpkl8UFlEteGeKIiMjdceOHU7HMMYYnzt9+jQdO3akQYMGvP766wQEBPj8mCISr6qRGS2zJ+qNMSYfK1++PGvWrOHw4cPcc889/P77747msaJijDH5XJkyZfjoo4/47bff6NatGxcuXHAsixUVY4wpAEqWLMnChQupVKkSsbGx/PLLL47ksKJijDEFRNGiRZk5cyYRERG0bt2aEydO5HkGKyrGGFOAFClShClTptC1a1dHOqL0576/jDHGXAMR4amnnqJcuXK0aNGCVatWUbt27Tw5thUVY4wpoB566CGCgoJo1aoVy5cvp9aRBM6+MInk48cJCA6mzLixlI7r7tVjWlExxpgCbNCgQQQFBTG5fQeeK1WaIu5bjpMTEzkzZiyAVwuLXVMxxpgCrkePHjx1Y+XUgpJCL17k7AuTvHosKyrGGFMIFP/ppwznJx/PtJ/ea2JFxRhjCoGA4OAczb9WVlSMMaYQKDNuLBIYeNU8CQykzLixXj2OXag3xphCIOVivN39ZYwxxitKx3X3ehFJz05/GWOM8RorKsYYY7zGiooxxhivsaJijDHGa6yoGGOM8RorKsYYY7zGiooxxhivEVV1OoPXichJ4Dunc6RRATjldAgPWE7vspzeZTm9J7cZq6lqxYwWFMii4m9EZIeqRjqdIzuW07ssp3dZTu/xZUY7/WWMMcZrrKgYY4zxGisqeeMNpwN4yHJ6l+X0LsvpPT7LaNdUjDHGeI21VIwxxniNFRVjjDFeY0Ull0TkZhFZLyIHReSAiIxyz79BRFaLyCH31+vTbDNeRL4Rka9EpF0e5w0QkV0issxfc4pIORFZKCJfun+uUX6a8yH3v/l+EZkrIiX9IaeIvCUiJ0Rkf5p5Oc4lIo1EZJ972VQRkTzI+ZL7332viCwRkXL+mDPNskdEREWkgr/mFJGR7iwHRORFn+dUVZtyMQE3ARHu12WAr4E6wIvAOPf8ccAk9+s6wB6gBFAdOAwE5GHe0cC7wDL3e7/LCcwC7ne/Lg6U87ecQBXgCBDofj8fGOAPOYE7gAhgf5p5Oc4FbAOiAAE+BjrkQc5YoKj79SR/zemefzOwEteD1hX8MSfQGlgDlHC/r+TrnNZSySVV/UFVd7pfnwUO4vrA6YbrwxH317vcr7sB81T1N1U9AnwDNM6LrCISAnQC3kwz269yikhZXL8cMwBU9XdV/cXfcroVBQJFpChQCjjuDzlVdSPwU7rZOcolIjcBZVX1M3V90vw3zTY+y6mqq1Q1yf12KxDijzndJgNjgLR3O/lbzgeAF1T1N/c6J3yd04qKF4lIKNAQ+By4UVV/AFfhASq5V6sCfJ9ms2PueXnhX7h+Ca6kmedvOW8BTgIz3afp3hSR0v6WU1UTgZeBo8APwBlVXeVvOdPIaa4q7tfp5+elQbj+UgY/yykiXYFEVd2TbpFf5QRqAS1E5HMR2SAit/s6pxUVLxGR64BFwN9U9desVs1gns/v6xaRzsAJVY33dJMM5uXF/edFcTXhX1fVhsB5XKdrMuPUz/N6XH/tVQeCgdIi0i+rTTKY5w/382eWy9G8IvIokATMSZmVSZ48zykipYBHgScyWpxJHid/n64HmgJ/B+a7r5H4LKcVFS8QkWK4CsocVV3snv0/d1MS99eUZucxXOdiU4TgOm3ia9FAVxFJAOYBd4rIO36Y8xhwTFU/d79fiKvI+FvONsARVT2pqpeBxUAzP8yZIqe5jvF/p57Szvc5EbkP6Az0dZ+C8becNXD9MbHH/fsUAuwUkcp+lhP3cReryzZcZykq+DKnFZVcclf9GcBBVX0lzaIPgPvcr+8D3k8z/24RKSEi1YFbcV0Y8ylVHa+qIaoaCtwNrFPVfn6Y80fgexGp7Z4VA3zhbzlxnfZqKiKl3P8HYnBdT/O3nClylMt9iuysiDR1f3/3ptnGZ0SkPTAW6KqqF9Ll94ucqrpPVSupaqj79+kYrpt1fvSnnG5LgTsBRKQWrhtfTvk0pzfvPiiME9AcV/NwL7DbPXUEygNrgUPurzek2eZRXHdbfIWX7wDxMHMr/u/uL7/LCYQDO9w/06W4mu/+mPMp4EtgPzAb1500jucE5uK6znMZ1wfe4GvJBUS6v7fDwGu4e+Dwcc5vcJ3rT/ld+rc/5ky3PAH33V/+lhNXEXnHfdydwJ2+zmndtBhjjPEaO/1ljDHGa6yoGGOM8RorKsYYY7zGiooxxhivsaJijDHGa6yomAJPRMqLyG739KOIJKZ5X9zpfGmJSCsRaeajfYeKyEUR2Z1m3rk0rzuKqxfjquLqgfmoiLzmiyym4CrqdABjfE1VT+N69gURmQicU9WXncojIkX1/zpNTK8VcA74NAf7C1DVZA9XP6yq4RnsIwZ4FYhV1aPAZBH5GdczC8Z4zFoqplByjxmxQUTiRWRlmi5MPhGRySKyUVxjudwuIovdf8E/614nVFxjfswS17gfC939QWW33+dFZAMwSkS6uDv52yUia0TkRnF1SPoX4CF3K6qFiLwtIj3T5D7n/tpKXOP4vAvsE9c4OS+JyHZ3pmE5+Fm0AKYDnVT1sFd+wKbQsqJiCiPB9Vd5T1VtBLwFPJdm+e+qegfwb1xdVIwAwoABIlLevU5t4A1VrQ/8CgwXVx9wWe23nKq2VNV/ApuBpurqNHMeMEZVE9zHnKyq4aq6KZvvozHwqKrWwfX09BlVvR24HRji7n4jOyXc3+NdqvqlB+sbkyU7/WUKoxK4isRqV/dGBODq3iLFB+6v+4AD6u4yXkS+xdUJ3y/A96q6xb3eO8BfgRXZ7Pe9NK9DgPfcLZniuAb8yqlt6hoLA1yDW9VP06oJwtWfU3b7vYzrVNtgYNQ1ZDDmKlZUTGEkuIpFVCbLf3N/vZLmdcr7lN+Z9P0bpXQbntV+z6d5/Srwiqp+ICKtgImZbJOE+4yCu4O/tDcWpN2fACNVdWUm+8nMFeDPwBoRmaCqz+dwe2OuYqe/TGH0G1BRRKLANXSBiNTN4T6qpmwP9MF1OuurHOw3CEh0v74vzfyzuIalTpEANHK/7gYUy2R/K4EH3KfgEJFa4hrcLFvq6g24M9BXRAZ7so0xmbGiYgqjK0BPYJKI7MHVG25Ob+M9CNwnInuBG3ANKvZ7DvY7EVggIptwdUWe4kOge8qFelwX0FuKyDagCVe3TtJ6E9cQATtFZD/wH3JwJkJVfwLaA4+JSDdPtzMmPeul2Jgcct+ltUxVw5zOkhM5zS0iA4BIVX3Ql7lMwWItFWMKj2QgKO3Dj5kRkYeA8bjubDPGY9ZSMcYY4zXWUjHGGOM1VlSMMcZ4jRUVY4wxXmNFxRhjjNdYUTHGGOM1/x8OVnMc/uwKbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Start figure.\n",
    "fig = plt.figure()\n",
    "ax = fig.add_axes([0.15,0.15,0.80,0.80])\n",
    "\n",
    "# Plot.\n",
    "ax.plot(T_rs, F_rs, '-', c='k', lw=1, label='Reversible Scaling')\n",
    "ax.plot(T_fl, F_fl, 'o', c=c[0], label='Frenkel-Ladd')\n",
    " \n",
    "# Add details and save figure.\n",
    "ax.set_xlabel(r'Temperature [K]')\n",
    "ax.set_ylabel(r'Free energy [eV/atom]')\n",
    "ax.legend(loc='best', frameon=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}