{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1afc4517",
   "metadata": {},
   "source": [
    "# Thermodynamic Derivatives\n",
    "\n",
    "## Helmholtz energy derivatives\n",
    "Thermodynamic derivatives are at the very heart of teqp.  All models are defined in the form $\\alpha^r(T, \\rho, z)$, where $\\rho$ is the molar density, and z are mole fractions.  There are exceptions for models for which the independent variables are in simulation units (Lennard-Jones and its ilk).\n",
    "\n",
    "Thereofore, to obtain the residual pressure, it is obtained as a derivative:\n",
    "\n",
    "$$ p^r = \\rho R T \\left( \\rho \\left(\\frac{\\partial \\alpha^r}{\\partial \\rho}\\right)_{T}\\right)$$\n",
    "\n",
    "and other residual thermodynamic properties are defined likewise.\n",
    "\n",
    "We can define the concise derivative\n",
    "\n",
    "$$ \\Lambda^r_{xy} = (1/T)^x(\\rho)^y\\left(\\frac{\\partial^{x+y}(\\alpha^r)}{\\partial (1/T)^x\\partial \\rho^y}\\right)$$\n",
    "\n",
    "so we can re-write the derivative above as \n",
    "\n",
    "$$ p^r = \\rho R T \\Lambda^r_{01}$$\n",
    "\n",
    "Similar definitions apply for all the other thermodynamic properties, with the tot superscript indicating it is the sum of the residual and ideal-gas (not included in teqp) contributions:\n",
    "\n",
    "$$\n",
    "\\frac{p}{\\rho R T}  = 1 + \\Lambda^r_{01}\n",
    "$$\n",
    "Internal energy ($u= a+Ts$):\n",
    "$$\n",
    "\t\\frac{u}{RT} = \\Lambda^{\\rm tot}_{10}\n",
    "$$\n",
    "Enthalpy ($h= u+p/\\rho$):\n",
    "$$\n",
    "\\frac{h}{RT} = 1+\\Lambda^r_{01} + \\Lambda^{\\rm tot}_{10}\n",
    "$$\n",
    "Entropy ($s\\equiv -(\\partial a/\\partial T)_v$):\n",
    "$$\n",
    "\\frac{s}{R} = \\Lambda^{\\rm tot}_{10}-\\Lambda^{\\rm tot}_{00}\n",
    "$$\n",
    "Gibbs energy ($g= h-Ts$):\n",
    "$$\n",
    "\t\\frac{g}{RT} = 1+\\Lambda^r_{01}+\\Lambda^{\\rm tot}_{00}\n",
    "$$\n",
    "Derivatives of pressure:\n",
    "$$\n",
    "\t\\left(\\frac{\\partial p}{\\partial \\rho}\\right)_{T} = RT\\left(1+2\\Lambda^r_{01}+\\Lambda^r_{02}\\right)\n",
    "$$\n",
    "$$\n",
    "\t\\left(\\frac{\\partial p}{\\partial T}\\right)_{\\rho} = R\\rho\\left(1+\\Lambda^r_{01}-\\Lambda^r_{11}\\right)\n",
    "$$\n",
    "Isochoric specific heat ($c_v\\equiv (\\partial u/\\partial T)_v$):\n",
    "$$\n",
    "\\frac{c_v}{R} = -\\Lambda^{\\rm tot}_{20}\n",
    "$$\n",
    "Isobaric specific heat ($c_p\\equiv (\\partial h/\\partial T)_p$; see Eq. 3.56 from Span for the derivation):\n",
    "$$\n",
    "\\frac{c_p}{R} = -\\Lambda^{\\rm tot}_{20}+\\frac{(1+\\Lambda^r_{01}-\\Lambda^r_{11})^2}{1+2\\Lambda^r_{01}+\\Lambda^r_{02}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "raw",
   "id": "4e621170",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "In teqp, these derivatives are obtained from methods like \n",
    "\n",
    "* :py:meth:`~teqp.teqp.AbstractModel.get_Arxy`\n",
    "* :py:meth:`~teqp.teqp.AbstractModel.get_Ar06n`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7622b9b7",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "where the ``A`` in this context indicates the variable $\\Lambda$ above. This naming is perhaps not ideal, since $A$ is sometimes the total Helmholtz energy, but it was a close visual mnemonic to the character $\\Lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2ced52fd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:06.372355Z",
     "iopub.status.busy": "2024-12-12T18:11:06.372201Z",
     "iopub.status.idle": "2024-12-12T18:11:06.388127Z",
     "shell.execute_reply": "2024-12-12T18:11:06.387651Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.22.0'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import teqp\n",
    "teqp.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ca35cc29",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:06.389785Z",
     "iopub.status.busy": "2024-12-12T18:11:06.389392Z",
     "iopub.status.idle": "2024-12-12T18:11:06.436438Z",
     "shell.execute_reply": "2024-12-12T18:11:06.435998Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8d254c40",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:06.438521Z",
     "iopub.status.busy": "2024-12-12T18:11:06.438292Z",
     "iopub.status.idle": "2024-12-12T18:11:06.442029Z",
     "shell.execute_reply": "2024-12-12T18:11:06.441590Z"
    }
   },
   "outputs": [],
   "source": [
    "Tc_K = [300]\n",
    "pc_Pa = [4e6]\n",
    "acentric = [0.01]\n",
    "model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "49ee3453",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:06.444231Z",
     "iopub.status.busy": "2024-12-12T18:11:06.443903Z",
     "iopub.status.idle": "2024-12-12T18:11:06.448964Z",
     "shell.execute_reply": "2024-12-12T18:11:06.448517Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.06836660379313926"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = np.array([1.0])\n",
    "model.get_Ar01(300,300,z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "498a2350",
   "metadata": {},
   "source": [
    "And there are additional methods to obtain all the derivatives up to a given order:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3ae755e1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:06.450838Z",
     "iopub.status.busy": "2024-12-12T18:11:06.450668Z",
     "iopub.status.idle": "2024-12-12T18:11:06.455945Z",
     "shell.execute_reply": "2024-12-12T18:11:06.455495Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-6.96613834e-02, -6.83666038e-02,  2.53578225e-03, -1.57011622e-04,\n",
       "        1.68186288e-05, -2.23059409e-06,  3.82592585e-07])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.get_Ar06n(300,300,z) # derivatives 00, 01, 02, ... 06"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd644a3e",
   "metadata": {},
   "source": [
    "But more derivatives are slower than fewer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4e0609ae",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:06.457956Z",
     "iopub.status.busy": "2024-12-12T18:11:06.457789Z",
     "iopub.status.idle": "2024-12-12T18:11:20.083575Z",
     "shell.execute_reply": "2024-12-12T18:11:20.083031Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "567 ns ± 0.747 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.11 μs ± 4.12 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit model.get_Ar01(300,300,z)\n",
    "%timeit model.get_Ar04n(300,300,z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d05e9070",
   "metadata": {},
   "source": [
    "Note: calling overhead is usually on the order of 1 microsecond"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "502e9830",
   "metadata": {},
   "source": [
    "## Virial coefficients\n",
    "\n",
    "Virial coefficients represent the thermodynamics of the interaction of two-, three-, ... bodies interacting with each other.  They can be obtained rigorously if the potential energy surface of interaction is fully known.  In general, such a surface can only be constructed for small rigid molecules.  Many simple thermodynamic models do a poor job of predicting the thermodynamics captured by the virial coefficients.  \n",
    "\n",
    "The i-th virial coefficient is defined by \n",
    "$$\n",
    "\tB_i = \\frac{(\\alpha^r)^{(i-1)}}{(i-2)!}\n",
    "$$\n",
    "with the concise derivative term\n",
    "$$\n",
    "\t(\\alpha^r)^{(i)} = \\lim_{\\rho \\to 0} \\left(\\frac{\\partial^{i}\\alpha^r}{\\partial \\rho^{i}}\\right)_{T,\\vec{x}}\n",
    "$$\n",
    "\n",
    "teqp supports the virial coefficient directly, there is the ``get_B2vir`` method for the second virial coefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "720e120c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:20.085342Z",
     "iopub.status.busy": "2024-12-12T18:11:20.085169Z",
     "iopub.status.idle": "2024-12-12T18:11:20.089129Z",
     "shell.execute_reply": "2024-12-12T18:11:20.088764Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.00023661263734465424"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.get_B2vir(300, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fba4369",
   "metadata": {},
   "source": [
    "And the ``get_Bnvir`` method that allows for the calculation of higher virial coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "92d01992",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:20.090673Z",
     "iopub.status.busy": "2024-12-12T18:11:20.090476Z",
     "iopub.status.idle": "2024-12-12T18:11:20.093822Z",
     "shell.execute_reply": "2024-12-12T18:11:20.093423Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{2: -0.0002366126373446542,\n",
       " 3: 3.001768410777936e-08,\n",
       " 4: -3.2409760373816364e-12,\n",
       " 5: 3.961781646633723e-16,\n",
       " 6: -4.5529239838367004e-20,\n",
       " 7: 5.375927851118494e-24}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.get_Bnvir(7, 300, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d50998b1",
   "metadata": {},
   "source": [
    "The ``get_Bnvir`` method was implemented because when doing automatic differentiation, all the intermediate derivatives are also retained.\n",
    "\n",
    "There is also a method to calculate temperature derivatives of a given virial coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fa4ff386",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:20.095282Z",
     "iopub.status.busy": "2024-12-12T18:11:20.095127Z",
     "iopub.status.idle": "2024-12-12T18:11:20.098604Z",
     "shell.execute_reply": "2024-12-12T18:11:20.098083Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0095625628421257e-10"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.get_dmBnvirdTm(2, 3, 300, z) # third temperature derivative of the second virial coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9a9823b",
   "metadata": {},
   "source": [
    "## Isochoric Thermodynamics Derivatives\n",
    "\n",
    "In the isochoric thermodynamics formalism, the EOS is expressed in the Helmholtz energy  density $\\Psi$ as a function of temperature and molar densities $\\vec\\rho$. This formalism is handy because it allows for a concise mathematical structure, well suited to implementation in teqp. For instance the pressure is obtained from (see https://doi.org/10.1002/aic.16074):\n",
    "$$\n",
    "p=-\\Psi + \\sum_{i=1}^N\\rho_i\\mu_i\n",
    "$$\n",
    "with the chemical potential $\\mu_i$ obtained from\n",
    "$$\n",
    "\\mu_i = \\left(\\frac{\\partial \\Psi}{\\partial \\rho_i}\\right)_{T,\\rho_{j\\neq i}}\n",
    "$$\n",
    "The molar densities $\\rho_i$ are related to the total density and the mole fractions:\n",
    "$$\n",
    "\\rho_i = x_i\\rho\n",
    "$$\n",
    "\n",
    "In teqp, the isochoric derivative functions like ``get_fugacity_coefficients``, ``get_partial_molar_volumes`` take as arguments the temperature ``T`` and the vector of molar concentrations ``rhovec``=$\\vec\\rho$, which are obtained by multiplying the mole fractions by the total density.\n",
    "\n",
    "Example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d56640d3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-12-12T18:11:20.100127Z",
     "iopub.status.busy": "2024-12-12T18:11:20.099872Z",
     "iopub.status.idle": "2024-12-12T18:11:20.144915Z",
     "shell.execute_reply": "2024-12-12T18:11:20.144518Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.97884567, 0.99866748])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "array([0.00470644, 0.00480351])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = teqp.build_multifluid_model([\"CO2\",\"Argon\"], teqp.get_datapath())\n",
    "T, rhovec = 300, np.array([0.3,0.4])*300 # K, mol/m^3\n",
    "display(model.get_fugacity_coefficients(T, rhovec))\n",
    "display(model.get_partial_molar_volumes(T, rhovec))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}