{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ideal gas equation of state using grand canonical ensemble transition-matrix Monte Carlo\n",
    "\n",
    "In this example, the ideal gas equation of state is obtained as a test of the flat histogram method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Usage: feasst < file.txt\n",
      "# For more information, use the command \"feasst-menu\"\n",
      "# Exit with ctrl-c\n",
      "FEASST version 0.25.19\n",
      "MonteCarlo\n",
      "Configuration cubic_side_length=8 particle_type=atom:/feasst/particle/atom.txt\n",
      "Potential Model=IdealGas\n",
      "ThermoParams beta=0.8333333333333334 chemical_potential=-3\n",
      "FlatHistogram Bias=TransitionMatrix Macrostate=MacrostateNumParticles max=50 min=0 min_sweeps=100 width=1\n",
      "TrialTransfer particle_type=atom\n",
      "CriteriaUpdater trials_per_update=1e5\n",
      "CriteriaWriter output_file=id_fh.txt trials_per_write=1e5\n",
      "Run until=complete\n",
      "# Initializing random number generator with seed: 1781198252\n",
      " \n",
      " exit: 0\n"
     ]
    }
   ],
   "source": [
    "params={\"cubic_side_length\": 8, \"beta\": 1./1.2, \"mu\": -3}\n",
    "script=\"\"\"\n",
    "MonteCarlo\n",
    "Configuration cubic_side_length={cubic_side_length} particle_type=atom:/feasst/particle/atom.txt\n",
    "Potential Model=IdealGas\n",
    "ThermoParams beta={beta} chemical_potential={mu}\n",
    "FlatHistogram Macrostate=MacrostateNumParticles width=1 min=0 max=50 \\\n",
    "              Bias=TransitionMatrix min_sweeps=100\n",
    "TrialTransfer particle_type=atom\n",
    "CriteriaUpdater trials_per_update=1e5\n",
    "CriteriaWriter trials_per_write=1e5 output_file=id_fh.txt\n",
    "Run until=complete\n",
    "\"\"\".format(**params)\n",
    "\n",
    "def run(script):\n",
    "    with open('script0.txt', 'w') as file: file.write(script)\n",
    "    import subprocess\n",
    "    syscode = subprocess.call(\"feasst < script0.txt > script0.log\", shell=True, executable='/bin/bash')\n",
    "    with open('script0.log', 'r') as file: print(file.read(), '\\n', 'exit:', syscode)\n",
    "run(script)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the ideal gas relationship, $\\beta P = \\rho$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N betaPV difference\n",
      "0.2829214450028523 0.28276142935893817 0.00016001564391415757\n",
      "0.3127110344609224 0.3125181844159339 0.00019285004498853509\n",
      "0.3456402287398188 0.34540825464818914 0.00023197409162967197\n",
      "0.3820402993452214 0.3817618872672399 0.0002784120779814603\n",
      "0.4222774542216826 0.421944174504675 0.0003332797170076396\n",
      "0.46675648379027573 0.4663587266639026 0.00039775712637313276\n",
      "0.5159247687089498 0.5154517276198478 0.0004730410891019554\n",
      "0.570276678959833 0.5697164104095647 0.0005602685502683613\n",
      "0.6303583950694145 0.6296979935758966 0.0006604014935178704\n",
      "0.696773183838372 0.695999122080295 0.0007740617580769271\n",
      "0.770187163519557 0.7692858599548528 0.0009013035647041923\n",
      "0.8513355980068428 0.8502942855654231 0.001041312441419695\n",
      "0.9410297679934123 0.9398377446931344 0.0011920233002779002\n",
      "1.0401644817643028 1.0388148221072249 0.0013496596570778951\n",
      "1.1497263127984105 1.14821809969524 0.0015082131031705082\n",
      "1.2708026900299596 1.2691437797302167 0.001658910299742855\n",
      "1.4045920241185472 1.4028022671395293 0.0017897569790179535\n",
      "1.552415132895656 1.5505298267730436 0.0018853061226122847\n",
      "1.715728330676701 1.7138014629021896 0.0019268677745114537\n",
      "1.8961386595276704 1.8942452102998013 0.0018934492278690485\n",
      "2.0954218391723356 2.093658079286071 0.0017637598862645376\n",
      "2.315543542873657 2.314023957202639 0.0015195856710179534\n",
      "2.558684481248569 2.5575338250073445 0.001150656241224457\n",
      "2.8272693721483746 2.826608678578439 0.0006606935699355887\n",
      "3.123999060529628 3.1239255153079615 7.354522166647115e-05\n",
      "3.4518837585846183 3.452446612157545 -0.0005628535729265138\n",
      "3.8142737410538046 3.8154520379532544 -0.0011782968994498155\n",
      "4.214882414380421 4.216574899646911 -0.0016924852664903511\n",
      "4.6577966791061085 4.659838295506624 -0.0020416164005157\n",
      "5.147472750233575 5.149692566843759 -0.002219816610184111\n",
      "5.688723927088134 5.691051625222572 -0.0023276981344384495\n",
      "6.286720429404237 6.289328426302422 -0.0026079968981846946\n",
      "6.947036426917338 6.950472444007084 -0.0034360170897462794\n",
      "7.675785778802292 7.681015960238509 -0.005230181436216341\n",
      "8.47987119024142 8.48813955101325 -0.008268360771829464\n",
      "9.367320311579844 9.379767360418134 -0.012447048838289732\n",
      "10.347603415096874 10.364696286167293 -0.017092871070419235\n",
      "11.431761006992256 11.452748968406713 -0.020987961414457246\n",
      "12.63219441581632 12.654923396269563 -0.022728980453242897\n",
      "13.962166478423628 13.98350540546874 -0.021338927045112754\n",
      "15.435377036087809 15.452130769149123 -0.016753733061314335\n",
      "17.066095896279855 17.075829672218084 -0.009733775938229883\n",
      "18.869827844735347 18.871115140125454 -0.0012872953901066353\n",
      "20.86363652307769 20.8561301260069 0.007506397070791593\n",
      "23.065498198040224 23.05077730812369 0.01472088991653564\n",
      "25.493907587206714 25.47677929014812 0.0171282970585942\n",
      "28.16900640379734 28.15778058598642 0.011225817810920802\n",
      "31.109107602401853 31.11944202979696 -0.010334427395108037\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fh=pd.read_csv('id_fh.txt', comment=\"#\")\n",
    "#print(fh)\n",
    "print('N betaPV difference')\n",
    "for delta_conjugate in np.arange(-5, 5, 0.1):\n",
    "    fh['ln_prob_rw'] = fh['ln_prob'] + fh['state']*delta_conjugate - fh['ln_prob'].min()  # avoid negative log\n",
    "    fh['ln_prob_rw'] -= np.log(sum(np.exp(fh['ln_prob_rw'])))   # renormalize\n",
    "    if fh['ln_prob_rw'].values[-1] < -6:\n",
    "        #plt.plot(fh['ln_prob_rw'])\n",
    "        N = (np.exp(fh[\"ln_prob_rw\"]) * fh[\"state\"]).sum()\n",
    "        betaPV = -fh[\"ln_prob_rw\"][0] - np.log(np.exp(fh[\"ln_prob_rw\"]).sum())\n",
    "        print(N, betaPV, N-betaPV)\n",
    "        assert np.abs(1 - betaPV/N) < 1e-2"
   ]
  },
  {
   "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
}
