{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of a two-phase (phase separated) simulation\n",
    "\n",
    "In this [previous tutorial](tutorial_02_reweight_supercritical.rst), we analyzed a simulation conducted above the critical temperature.\n",
    "In this tutorial, we now consider an [LennardJones](../../system/doc/LennardJones_arguments.rst) simulation with two stable phases at $T^*=1.2$.\n",
    "\n",
    "This tutorial uses data from the [NIST SRSW](https://www.nist.gov/programs-projects/nist-standard-reference-simulation-website). For an example of running and analyzing a phase separated system, see [tutorial 10](tutorial_10_wltm_lowt.rst) in the [flat_histogram plugin](../README.rst)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from feasst import macrostate_distribution\n",
    "\n",
    "srsw = macrostate_distribution.MacrostateDistribution(\n",
    "    file_name=\"../test/data/stat120.csv\",\n",
    "    macrostate_header=\"N\",\n",
    "    ln_prob_header=\"lnPI\")\n",
    "\n",
    "# Set show_plot to True to see the plots. Plots take up too much space in the git repo.\n",
    "show_plot = False\n",
    "if show_plot:\n",
    "    %matplotlib inline\n",
    "    import matplotlib.pyplot as plt\n",
    "    plt.plot(srsw.ln_prob())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this macrostate distribution may contain a second maximum depending on the value of the chemical potential.\n",
    "But at the chemical potential shown above, the low density phase is unstable and there is no minimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(srsw.minimums()) == 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, a minima appears after reweighting to a lower chemical potential. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "srsw_rw = srsw.reweight(delta_beta_mu=-0.1)\n",
    "minima=srsw_rw.minimums()\n",
    "assert len(minima) == 1 and minima.values[0] == 152\n",
    "\n",
    "if show_plot:\n",
    "    plt.plot(srsw_rw.ln_prob())\n",
    "    plt.gca().axvline(minima.values[0], color='black')\n",
    "    plt.text(20, -40, 'low density phase')\n",
    "    plt.text(200, -40, 'high density phase\\n(more probable)')\n",
    "    plt.xlabel('N')\n",
    "    plt.ylabel(r'$\\ln\\Pi$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the chemical potential at phase equilibrium subject to the constraint that the probability of observing the two phases are equal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_beta_mu_equilibrium = srsw.equilibrium()\n",
    "mins=srsw.minimums()\n",
    "assert len(mins)==1\n",
    "mn=mins.values[0]\n",
    "\n",
    "if show_plot:\n",
    "    plt.plot(srsw.ln_prob())\n",
    "    plt.gca().axvline(mn, color='black')\n",
    "    plt.text(20, -40, 'saturated vapor')\n",
    "    plt.text(225, -40, 'saturated liquid')\n",
    "    plt.xlabel('N')\n",
    "    plt.ylabel(r'$\\ln\\Pi$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare chemical potential and equilibrium properties with [published results](https://www.nist.gov/mml/csd/chemical-informatics-research-group/sat-tmmc-liquid-vapor-coexistence-properties-long-range)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pressure: 0.07722557238264459\n",
      "vapor_density: 0.10035100340092687\n",
      "liquid_density: 0.563186732420252\n"
     ]
    }
   ],
   "source": [
    "assert len(mins) == 1 and mins.values[0] == 166\n",
    "assert abs(-0.128071 - delta_beta_mu_equilibrium) < 0.01\n",
    "\n",
    "volume = 8**3\n",
    "vapor, liquid = srsw.split()\n",
    "betaPV = -vapor.ln_prob()[0]\n",
    "beta=1/1.2\n",
    "pressure=betaPV/beta/volume\n",
    "print('pressure:', pressure)\n",
    "assert abs(0.07722557238264459 - pressure) < 1e-8\n",
    "\n",
    "vapor_density = vapor.average_macrostate()/volume\n",
    "print('vapor_density:', vapor_density)\n",
    "assert abs(0.10035100340092683 - vapor_density) < 1e-8\n",
    "\n",
    "liquid_density = liquid.average_macrostate()/volume\n",
    "print('liquid_density:', liquid_density)\n",
    "assert abs(0.5631867324202517 - liquid_density) < 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please [contact](../../../CONTACT.rst) us. Any feedback is appreciated!"
   ]
  }
 ],
 "metadata": {
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
