{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Comparison of equations of state for model fluids\n", "\n", "In this notebook, we take a look at various model fluid equations of state and compare them.\n", "These EoS are:\n", "\n", "- The PeTS equation of state (from `feos.pets`),\n", "- the uv-theory equation of state (from `feos.uvtheory`), \n", "- the uv-B3-theory equation of state (from `feos.uvtheory` using `VirialOrder.Third`), and\n", "- the Thol equation of state (implemented in this notebook as Python class).\n", "\n", "The PeTS equation of state is for a Lennard-Jones fluid with cut off distance of $r_c = 2.5 \\sigma$ and a potential shift while the uv-theory and Thol equation of state model the full Lennard-Jones potential. Note that the uv-theory was developed to model Mie potentials and not primarily adjusted to the Lennard-Jones fluid. Thus, we don't expect similar results when comparing uv-theory and Thol to PeTS." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/bauer/opt/anaconda3/lib/python3.9/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.24.2\n", " warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n" ] } ], "source": [ "import numpy as np\n", "from feos.si import *\n", "from feos.uvtheory import UVParameters, Perturbation, VirialOrder\n", "from feos.pets import PetsParameters\n", "from feos.eos import State, PhaseEquilibrium, PhaseDiagram, EquationOfState, Contributions\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "import pandas as pd\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "sns.set_context('talk')\n", "sns.set_palette(\"hls\", 8)\n", "sns.set_style('ticks')\n", "colors = sns.palettes.color_palette('Dark2', 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implement the Thol equation of state as Python class" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Parameters for Thol EoS:\n", "A = np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 2.0,\n", " 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0])\n", "N = np.array([0.005208073, 2.186252000, -2.161016000, 1.452700000, -2.041792000, 0.186952860, -0.090988445, \n", " -0.497456100, 0.109014310, -0.800559220, -0.568839000, -0.620862500, -1.466717700, 1.891469000, \n", " -0.138370100, -0.386964500, 0.126570200, 0.605781000, 1.179189000, -0.477326790, -9.921857500, -0.574793200, 0.003772923])\n", "T = np.array([1.000, 0.320, 0.505, 0.672, 0.843, 0.898, 1.294, 2.590, 1.786, 2.770, 1.786,\n", " 1.205, 2.830, 2.548, 4.650, 1.385, 1.460, 1.351, 0.660, 1.496, 1.830, 1.616, 4.970])\n", "D = np.array([4.0, 1.0, 1.0, 2.0, 2.0, 3.0, 5.0, 2.0, 2.0, 3.0, 1.0,\n", " 1.0, 1.0, 1.0, 2.0, 3.0, 3.0, 2.0, 1.0, 2.0, 3.0, 1.0, 1.0])\n", "L = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 1.0, 2.0, 2.0,\n", " 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])\n", "ETA = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.067, 1.522,\n", " 8.82, 1.722, 0.679, 1.883, 3.925, 2.461, 28.2, 0.753, 0.82])\n", "BETA = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.625,\n", " 0.638, 3.91, 0.156, 0.157, 0.153, 1.16, 1.73, 383, 0.112, 0.119])\n", "GAMMA = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.71,\n", " 0.86, 1.94, 1.48, 1.49, 1.945, 3.02, 1.11, 1.17, 1.33, 0.24])\n", "EPSILON = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2053,\n", " 0.409, 0.6, 1.203, 1.829, 1.397, 1.39, 0.539, 0.934, 2.369, 2.43])\n", "\n", "class Thol:\n", " def __init__(self):\n", " self.sigma = 3.7039 * ANGSTROM\n", " self.eps_k = 150.03 * KELVIN\n", " self.tc = 1.32 * self.eps_k / KELVIN\n", " self.rhoc = 0.31 / self.sigma**3 * ANGSTROM**3\n", " \n", " def components(self): \n", " return 1\n", " \n", " def subset(self, components):\n", " return self\n", " \n", " def molar_weight(self):\n", " return np.array([1.0])\n", " \n", " def max_density(self, moles):\n", " return 0.04\n", " \n", " def helmholtz_energy(self, state):\n", " \"\"\"\n", " state (StateHD):\n", " temperature in Kelvin als Float, Dual oder HD oder HD3, HDD, HDD3,\n", " partial_density in # / Angstrom^3\n", " volume in Angstrom^3\n", " moles in mol\n", " \"\"\"\n", " tau = self.tc / state.temperature\n", " delta = np.sum(state.partial_density) / self.rhoc\n", " a = 0.0\n", " \n", " for i in range(6):\n", " a = a + N[i] * delta**D[i] * tau**T[i]\n", " for i in range(6, 12):\n", " a = a + N[i] * delta**D[i] * tau**T[i] * np.exp(-1.0 * delta**L[i])\n", " for i in range(12, 23):\n", " a = a + N[i] * delta**D[i] * tau**T[i] * np.exp(- 1.0 * ETA[i] * (delta - EPSILON[i])**2.0 - 1.0 * BETA[i] * (tau - 1.0 * GAMMA[i])**2.0)\n", " return a * np.sum(state.moles)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Instantiate all equations of state" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```rust\n", "let omega = if t_x.re() < 175.0 {\n", " (-t_x * CU_WCA[5] * (reduced_density - CU_WCA[6]).powi(2)).exp()\n", " * ((t_x * CU_WCA[7]).tanh().recip() - 1.0).powi(2)\n", " } else {\n", " (-t_x * CU_WCA[5] * (reduced_density - CU_WCA[6]).powi(2)).exp() * 0.0\n", " };\n", "\n", " -(-(reduced_density.powi(2) * ((rep_x + CU_WCA[4]).recip() * CU_WCA[3] + CU_WCA[2]) + omega))\n", " .exp()\n", " + 1.0\n", "```" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "sigma = 3.7039\n", "sigma_a = sigma * ANGSTROM\n", "eps_k = 150.03\n", "eps_k_k = eps_k * KELVIN\n", "rep = 12.0\n", "parameters = UVParameters.new_simple(rep, 6.0, sigma, eps_k)\n", "CU_WCA = [26.454, 1.8045, 1.7997, 161.96, 11.605, 12., 0.4, 2.0]\n", "\n", "uvtheory_b3 = EquationOfState.uvtheory(parameters, virial_order=VirialOrder.Third)\n", "\n", "def ufraction(state):\n", " t_x = state.temperature / eps_k\n", " reduced_density = state.density * sigma**3\n", " rep_x = rep\n", " omega = 0.0\n", " if t_x.value < 175:\n", " omega = (-t_x * CU_WCA[5] * (reduced_density - CU_WCA[6])**2).exp() \\\n", " * (1.0 / (t_x * CU_WCA[7]).tanh() - 1.0)**2\n", " ufrac = -(-(reduced_density**2 * (1.0 / (rep_x + CU_WCA[4]) * CU_WCA[3] + CU_WCA[2]) + omega)).exp() + 1.0\n", " return ufrac\n", "\n", "uvtheory_b3_uf = EquationOfState.uvtheory(parameters, virial_order=VirialOrder.Third, ufraction=ufraction)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare Helmholtz energies at single state point" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "u-fraction: 0.9996284573272382 + [-0.00011029124189623964]ε1 + [-0.000026521524493183303]ε1²\n", "u-fraction: 0.9988379797559929 + [-0.00027300538747495245]ε1 + [-0.00004990266174618351]ε1²\n", "u-fraction: 0.9966746965481094 + [-0.0006062062111358156]ε1 + [-0.00008146147063767258]ε1²\n", "u-fraction: 0.9912931419552962 + [-0.001203092119270068]ε1 + [-0.00011367301308275238]ε1²\n", "u-fraction: 0.9791404839344673 + [-0.0021238832321700614]ε1 + [-0.00013243123583561224]ε1²\n", "u-fraction: 0.9542743833113834 + [-0.0033138426175314743]ε1 + [-0.00012339313039826935]ε1²\n", "u-fraction: 0.9082877185095404 + [-0.004529474333716721]ε1 + [-0.00008324992172326038]ε1²\n", "u-fraction: 0.8316907324621025 + [-0.005354527151690742]ε1 + [-0.000026953387161395334]ε1²\n", "u-fraction: 0.7173811027783104 + [-0.005369248730264456]ε1 + [0.000019078372866783005]ε1²\n", "u-fraction: 0.5657828152983877 + [-0.004424758231726879]ε1 + [0.00003598597556005315]ε1²\n", "u-fraction: 0.38958466059040253 + [-0.0028312484896880603]ε1 + [0.000026773517218010987]ε1²\n", "u-fraction: 0.21484307935875258 + [-0.0012491170649092847]ε1 + [0.000010336886840306822]ε1²\n", "u-fraction: 0.07594022435560499 + [-0.00027430427708941533]ε1 + [0.000001465064860887849]ε1²\n", "u-fraction: 0.00492400656989278 + [-0.00000461539514104226]ε1 + [0.00000000648383857547817]ε1²\n", "u-fraction: 0.000012350254504078784 + [-0.0000000005797770223328869]ε1 + [0.00000000000004085836761109884]ε1²\n", "u-fraction: 0.000013087009891443735 + [-0.0000000006324654622769026]ε1 + [0.00000000000004588267534449147]ε1²\n", "u-fraction: 0.00001308704432334551 + [-0.0000000006324679601653612]ε1 + [0.000000000000045882916959367675]ε1²\n", "u-fraction: 0.000013044906825188107 + [-0.0000000006294135239265427]ε1 + [0.000000000000045587706390672894]ε1²\n", "u-fraction: 0.000013087044213766497 + [-0.0000000006324679522144209]ε1 + [0.000000000000045882916190291433]ε1²\n", "u-fraction: 0.00001308704432334551 + [-0.0000000006324679601653662]ε1 + [0.00000000000004588291695936817]ε1²\n", "u-fraction: 0.00001308704432334551 + [-0.0000000006324679601653662]ε\n", "u-fraction: 0.00001308704432334551 + [-0.0000000006324679601653665]ε\n" ] } ], "source": [ "s1 = State(uvtheory_b3, temperature=300*KELVIN, pressure=1*BAR)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "u-fraction: 0.9996284573272382 + [-0.00011029124189623964]ε1 + [-0.000026521524493183303]ε1²\n", "u-fraction: 0.9988379797559929 + [-0.00027300538747495245]ε1 + [-0.00004990266174618351]ε1²\n", "u-fraction: 0.9966746965481094 + [-0.0006062062111358156]ε1 + [-0.00008146147063767258]ε1²\n", "u-fraction: 0.9912931419552962 + [-0.001203092119270068]ε1 + [-0.00011367301308275238]ε1²\n", "u-fraction: 0.9791404839344673 + [-0.0021238832321700614]ε1 + [-0.00013243123583561224]ε1²\n", "u-fraction: 0.9542743833113834 + [-0.0033138426175314743]ε1 + [-0.00012339313039826935]ε1²\n", "u-fraction: 0.9082877185095404 + [-0.004529474333716721]ε1 + [-0.00008324992172326038]ε1²\n", "u-fraction: 0.8316907324621025 + [-0.005354527151690742]ε1 + [-0.000026953387161395334]ε1²\n", "u-fraction: 0.7173811027783104 + [-0.005369248730264456]ε1 + [0.000019078372866783005]ε1²\n", "u-fraction: 0.5657828152983877 + [-0.004424758231726879]ε1 + [0.00003598597556005315]ε1²\n", "u-fraction: 0.38958466059040253 + [-0.0028312484896880603]ε1 + [0.000026773517218010987]ε1²\n", "u-fraction: 0.21484307935875258 + [-0.0012491170649092847]ε1 + [0.000010336886840306822]ε1²\n", "u-fraction: 0.07594022435560499 + [-0.00027430427708941533]ε1 + [0.000001465064860887849]ε1²\n", "u-fraction: 0.00492400656989278 + [-0.00000461539514104226]ε1 + [0.00000000648383857547817]ε1²\n", "u-fraction: 0.000012350254504078784 + [-0.0000000005797770223328869]ε1 + [0.00000000000004085836761109884]ε1²\n", "u-fraction: 0.000013087009891443735 + [-0.0000000006324654622769026]ε1 + [0.00000000000004588267534449147]ε1²\n", "u-fraction: 0.00001308704432334551 + [-0.0000000006324679601653612]ε1 + [0.000000000000045882916959367675]ε1²\n", "u-fraction: 0.000013044906825188107 + [-0.0000000006294135239265427]ε1 + [0.000000000000045587706390672894]ε1²\n", "u-fraction: 0.000013087044213766497 + [-0.0000000006324679522144209]ε1 + [0.000000000000045882916190291433]ε1²\n", "u-fraction: 0.00001308704432334551 + [-0.0000000006324679601653662]ε1 + [0.00000000000004588291695936817]ε1²\n", "u-fraction: 0.00001308704432334551 + [-0.0000000006324679601653662]ε\n", "u-fraction: 0.00001308704432334551 + [-0.0000000006324679601653665]ε\n" ] } ], "source": [ "s2 = State(uvtheory_b3_uf, temperature=300*KELVIN, pressure=1*BAR)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "u-fraction: 1.9996000799840032 + [0.006665333599946677]ε\n" ] }, { "data": { "text/latex": [ "$96.684\\,\\mathrm{\\frac{ J}{molK}}$" ], "text/plain": [ "96.68386253040592 J/mol/K" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s2.molar_entropy()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[('Ideal gas (default)', -4.8170398860180577e-20 J),\n", " ('Hard Sphere', 9.709460360314697e-24 J),\n", " ('Reference Perturbation', 6.640075859812054e-25 J),\n", " ('Attractive Perturbation', -1.7058012153921695e-23 J)]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s1.helmholtz_energy_contributions()" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[('Ideal gas (default)', -4.8698028896541795e-20 J),\n", " ('Hard Sphere', 8.547391383441759e-24 J),\n", " ('Reference Perturbation', 5.844972740477965e-25 J),\n", " ('Attractive Perturbation', 5.475676104748454e-22 J)]" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s2.helmholtz_energy_contributions()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.0016213224956115704" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s1.molar_helmholtz_energy(Contributions.ResidualNvt) / (RGAS * 300 * KELVIN)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.0016138652202999679" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s2.molar_helmholtz_energy(Contributions.ResidualNvt) / (RGAS * 300 * KELVIN)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.0016144150962601586" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s3.molar_helmholtz_energy(Contributions.ResidualNvt) / (RGAS * 300 * KELVIN)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.000940704010881386" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s4.molar_helmholtz_energy(Contributions.ResidualNvt) / (RGAS * 300 * KELVIN)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare critical points\n", "- Reference data for the Lennard-Jones from Potoff and Panagiotopoulos (MC): $T_c^*=1.3120$, $\\rho_c^*=0.316$, $p_c^*=0.1279$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "uvtheory_b3:\n", "Tc*=1.31329, Error: 0.099%\n", "pc*=0.1286, Error: 0.548%\n", "rhoc*=0.31425, Error: 0.553%\n", "\n", "uvtheory_wca:\n", "Tc*=1.30977, Error: 0.17%\n", "pc*=0.12448, Error: 2.674%\n", "rhoc*=0.30233, Error: 4.327%\n", "\n", "uvtheory_bh:\n", "Tc*=1.32083, Error: 0.673%\n", "pc*=0.13486, Error: 5.442%\n", "rhoc*=0.31865, Error: 0.838%\n", "\n", "thol:\n", "Tc*=1.32, Error: 0.61%\n", "pc*=0.13006, Error: 1.689%\n", "rhoc*=0.3132, Error: 0.886%\n", "\n" ] } ], "source": [ "def evaluate_critical_points():\n", " names = ['uvtheory_b3', 'uvtheory_wca', 'uvtheory_bh', 'thol']\n", " print('')\n", " i = 0\n", " for eos in [uvtheory_b3, uvtheory_wca, uvtheory_bh, thol]:\n", " rep = 12\n", " att = 6.0\n", " sigma = 3.7039 # in Angstrom\n", " eps_k = 150.03 # eps / kB in K\n", " print('{}:'.format(names[i]))\n", " Tc = State.critical_point(eos).temperature / eps_k_k\n", " print('Tc*={}, Error: {}%'.format(np.round(Tc, 5), np.round(np.abs(Tc - 1.3120 ) / 1.3120 * 100 ,3)))\n", " pc = State.critical_point(eos).pressure() * (sigma * ANGSTROM)**3 / KB / eps_k_k\n", " print('pc*={}, Error: {}%'.format(np.round(pc, 5), np.round(np.abs(pc - 0.1279 ) / 0.1279 * 100 ,3)))\n", " rhoc = State.critical_point(eos).density * (sigma * ANGSTROM)**3 * NAV\n", " print('rhoc*={}, Error: {}%'.format(np.round(rhoc, 5), np.round(np.abs(rhoc - 0.316) / 0.316 * 100 ,3)))\n", " i += 1\n", " print('')\n", " return\n", "\n", "evaluate_critical_points()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compare Phase diagrams" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 20.3 ms, sys: 667 µs, total: 21 ms\n", "Wall time: 20.9 ms\n" ] } ], "source": [ "%%time\n", "vle_uv = PhaseDiagram.pure(uvtheory_wca, 150*KELVIN, 500)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 16.8 ms, sys: 0 ns, total: 16.8 ms\n", "Wall time: 16.7 ms\n" ] } ], "source": [ "%%time\n", "vle_uv_bh = PhaseDiagram.pure(uvtheory_bh, 150*KELVIN, 500)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 24.8 ms, sys: 0 ns, total: 24.8 ms\n", "Wall time: 24.8 ms\n" ] } ], "source": [ "%%time\n", "vle_uv_b3 = PhaseDiagram.pure(uvtheory_b3, 150*KELVIN, 500)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.79 s, sys: 40.7 ms, total: 1.83 s\n", "Wall time: 1.83 s\n" ] } ], "source": [ "%%time\n", "vle_thol = PhaseDiagram.pure(thol, 150*KELVIN, 500)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 7.01 ms, sys: 4.16 ms, total: 11.2 ms\n", "Wall time: 11.1 ms\n" ] } ], "source": [ "%%time\n", "vle_pets = PhaseDiagram.pure(pets, 150*KELVIN, 500)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
| \n", " | t | \n", "rhov | \n", "rhol | \n", "p | \n", "
|---|---|---|---|---|
| 0 | \n", "0.70 | \n", "0.002012 | \n", "0.843237 | \n", "0.001381 | \n", "
| 10 | \n", "0.71 | \n", "0.002284 | \n", "0.839106 | \n", "0.001586 | \n", "
| 20 | \n", "0.72 | \n", "0.002584 | \n", "0.834868 | \n", "0.001815 | \n", "
| 30 | \n", "0.73 | \n", "0.002912 | \n", "0.830543 | \n", "0.002069 | \n", "
| 40 | \n", "0.74 | \n", "0.003271 | \n", "0.826147 | \n", "0.002350 | \n", "