notebook/fsac_u_check.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7bbf254d",
   "metadata": {},
   "source": [
    "# Energy calculations with F-SAC\n",
    "\n",
    "In this notebook, we use our SAC equations as the basis for F-SAC calculations.\n",
    "\n",
    "The purpose of this notebook is to check the energy calculations by Eq. (35) for\n",
    "the case of temperature dependent interaction energies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "3f57911f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys, math\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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c848bfd",
   "metadata": {},
   "source": [
    "### Compound definitions\n",
    "\n",
    "In F-SAC, compounds are defined with just a few surface areas and their respective charge densities.\n",
    "For a mixture of *n*-hexane and acetone it would be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "ef2e3471",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_1 = [0.0]\n",
    "Q_1     = [191.93]\n",
    "\n",
    "sigma_2 = [0,      0.0132, -0.00363]\n",
    "Q_2     = [20.167, 21.97,  80.23]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b473c56",
   "metadata": {},
   "source": [
    "## F-SAC implementation\n",
    "\n",
    "Now we implement the F-SAC model by deriving the basic `SAC` class and implementing the method calculating the interaction energies, `calc_u`, and its temperature derivative `calc_du_dT`.\n",
    "\n",
    "Here we use **temperature dependent** interaction energies by the following multiplying factor:\n",
    "\\begin{equation}\n",
    "f(T) = \\exp \\left (-\\beta \\frac{T-T_0}{1000 \\cdot \\text{K}} \\right )\n",
    "\\end{equation}\n",
    "\n",
    "This function returns 1 when $T=T_0$, being an exponential decay with the temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "879b371a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# F-SAC parameters\n",
    "alpha = 8544.6\n",
    "Q_eff = 3.597\n",
    "\n",
    "# The mixture area and sigma arrays\n",
    "Q = [Q_1, Q_2]\n",
    "sigma = [sigma_1, sigma_2]\n",
    "\n",
    "T0 = 308.15\n",
    "beta = 1\n",
    "\n",
    "class FSACNonHB(SAC):\n",
    "    def __init__(self):\n",
    "        super().__init__(Q_eff=Q_eff)\n",
    "\n",
    "    def calc_u(self, T, i, j, m, n):\n",
    "        '''\n",
    "        Interaction energies, only the electrostatic contribution (no hydrogen bonding)\n",
    "        and temperature dependent.\n",
    "        '''\n",
    "        f = math.exp(-beta * (T - T0)/1000)\n",
    "        return f*(alpha/2)*(sigma[i][m]+sigma[j][n])**2\n",
    "    \n",
    "    def calc_du_dT(self, T, i, j, m, n):\n",
    "        '''\n",
    "        Returns the derivative of the interaction energy.\n",
    "        '''\n",
    "        f = math.exp(-beta * (T - T0)/1000)\n",
    "        u = self.calc_u(T, i, j, m, n)\n",
    "        \n",
    "        return -beta/1000 * u * f\n",
    "\n",
    "\n",
    "# Create an instance of our F-SAC class\n",
    "sac = FSACNonHB()\n",
    "\n",
    "# Set a temperature and the areas of our compounds\n",
    "T = T0\n",
    "sac.set_compounds(Q)\n",
    "sac.set_temperature(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20b105e6",
   "metadata": {},
   "source": [
    "### Check the excess energy\n",
    "\n",
    "We can get the total internal and Helmholtz energies and subtract the pure values to get excess quantities. Additionally we can use standard thermodynamic relations to get the entropy.\n",
    "\n",
    "We first do this comparison with $\\beta = 0$ (interaction energies are not temperature dependent)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "accaa9bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [0.4, 0.6]\n",
    "beta = 0\n",
    "sac.set_composition(x)\n",
    "sac.calc_ln_gamma()\n",
    "\n",
    "ue = sac.get_energy()\n",
    "ae = sac.get_helmholtz()\n",
    "for i in range(len(Q)):\n",
    "    ue -= x[i] * sac.get_energy_pure(i)\n",
    "    ae -= x[i] * sac.get_helmholtz_pure(i)\n",
    "\n",
    "se_R = ue - ae"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4243509",
   "metadata": {},
   "source": [
    "We can also calculate the entropy by deriving the Helmholtz energy and then calculate back the energy. The values should match, within the accuracy of our finite differences derivative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "6d008cdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ae/RT: 0.4110998824992649\n",
      "ue/RT: 0.5121568520653366\n",
      "se/R: 0.10105696956607169\n",
      "ue/RT (check): 0.511953072515641\n",
      "se/R (check): 0.10085319001637616\n",
      "difference on ue/RT: -0.00020377954969552636\n",
      "difference on se/R: -0.00020377954969552636\n"
     ]
    }
   ],
   "source": [
    "dT = 0.1;\n",
    "T2 = T + dT\n",
    "sac.set_temperature(T2)\n",
    "sac.calc_ln_gamma()\n",
    "\n",
    "ae_T2 = sac.get_helmholtz()\n",
    "for i in range(len(Q)):\n",
    "    ae_T2 -= x[i] * sac.get_helmholtz_pure(i)\n",
    "\n",
    "se_R_check = -T*(ae_T2 - ae)/dT - ae;\n",
    "ue_check = ae + se_R_check\n",
    "\n",
    "print(f\"ae/RT: {ae}\")\n",
    "print(f\"ue/RT: {ue}\")\n",
    "print(f\"se/R: {se_R}\")\n",
    "print(f\"ue/RT (check): {ue_check}\")\n",
    "print(f\"se/R (check): {se_R_check}\")\n",
    "\n",
    "print(f\"difference on ue/RT: {ue_check-ue}\")\n",
    "print(f\"difference on se/R: {se_R_check-se_R}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4984cee8",
   "metadata": {},
   "source": [
    "Now we can do the same with $\\beta = 1$.\n",
    "The Hemlholtz excess energy should not change, since we are evaluating at the reference temperature $T_0$.\n",
    "The excess internal energy and entropy should change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "99db1705",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = 1\n",
    "sac.set_temperature(T)\n",
    "sac.calc_ln_gamma()\n",
    "\n",
    "ue = sac.get_energy()\n",
    "ae = sac.get_helmholtz()\n",
    "for i in range(len(Q)):\n",
    "    ue -= x[i] * sac.get_energy_pure(i)\n",
    "    ae -= x[i] * sac.get_helmholtz_pure(i)\n",
    "\n",
    "se_R = ue - ae"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51e7c1b3",
   "metadata": {},
   "source": [
    "The values calculated by finite differences should match again our new values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "e9274387",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ae/RT: 0.4110999207325188\n",
      "ue/RT: 0.6699780710468444\n",
      "se/R: 0.2588781503143256\n",
      "ue/RT (check): 0.6698595087236736\n",
      "se/R (check): 0.25875958799115484\n",
      "difference on ue/RT: -0.00011856232317075666\n",
      "difference on se/R: -0.00011856232317075666\n"
     ]
    }
   ],
   "source": [
    "dT = 0.1;\n",
    "T2 = T + dT\n",
    "sac.set_temperature(T2)\n",
    "sac.calc_ln_gamma()\n",
    "\n",
    "ae_T2 = sac.get_helmholtz()\n",
    "for i in range(len(Q)):\n",
    "    ae_T2 -= x[i] * sac.get_helmholtz_pure(i)\n",
    "\n",
    "se_R_check = -T*(ae_T2 - ae)/dT - ae;\n",
    "ue_check = ae + se_R_check\n",
    "\n",
    "print(f\"ae/RT: {ae}\")\n",
    "print(f\"ue/RT: {ue}\")\n",
    "print(f\"se/R: {se_R}\")\n",
    "print(f\"ue/RT (check): {ue_check}\")\n",
    "print(f\"se/R (check): {se_R_check}\")\n",
    "\n",
    "print(f\"difference on ue/RT: {ue_check-ue}\")\n",
    "print(f\"difference on se/R: {se_R_check-se_R}\")"
   ]
  }
 ],
 "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
}