{
 "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: /home/user/feasst/build/bin/fst < file.txt\n",
      "FEASST version 0.25.13\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: 1748436045\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(\"../../../build/bin/fst < 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.28280704118529504 0.28269896272478096 0.00010807846051408276\n",
      "0.3125716966244494 0.31244307289505413 0.0001286237293952608\n",
      "0.34547034653255504 0.3453177336655509 0.00015261286700413157\n",
      "0.38183303177462946 0.38165257153600674 0.0001804602386227261\n",
      "0.4220244921106349 0.42181192333104445 0.00021256877959047182\n",
      "0.4664477843615707 0.4661984840716928 0.0002493002898779073\n",
      "0.5155482686199676 0.5152573347805578 0.00029093383940981266\n",
      "0.5698179984963251 0.5694803887943721 0.0003376097019529567\n",
      "0.6298005554495112 0.6294112989521635 0.00038925649734766754\n",
      "0.6960963724179439 0.695650872279282 0.00044550013866184734\n",
      "0.7693685988173394 0.7688630436346048 0.0005055551827346516\n",
      "0.8503495683146425 0.8497814654371211 0.0005681028775214036\n",
      "0.9398479437373484 0.9392167773410379 0.0006311663963104452\n",
      "1.0387566314485 1.0380646280144998 0.0006920034340001102\n",
      "1.1480615821327125 1.1473145315787034 0.0007470505540090322\n",
      "1.2688516277131583 1.2680596545244074 0.0007919731887509052\n",
      "1.4023295457556206 1.4015076459041516 0.0008218998514690767\n",
      "1.5498245917662752 1.548992645140384 0.0008319466258912556\n",
      "1.7128067904295579 1.7119886283904553 0.0008181620391025657\n",
      "1.8929033155988746 1.8921242856350553 0.0007790299638192888\n",
      "2.091917290667462 2.0911996541790443 0.0007176364884178277\n",
      "2.311849265784693 2.3112047645079863 0.0006445012767066771\n",
      "2.554921421878632 2.55434057111087 0.0005808507677618202\n",
      "2.8236041555272093 2.8230424278659183 0.0005617276612910338\n",
      "3.1206440810345684 3.1200063039286365 0.0006377771059318249\n",
      "3.4490916966378316 3.448217801082639 0.000873895555192572\n",
      "3.8123262159536124 3.8109838190477983 0.0013423969058141694\n",
      "4.214074826040011 4.211966446490894 0.002108379549116357\n",
      "4.658424613796119 4.655218417077476 0.003206196718642751\n",
      "5.1498283637039 5.14521942319203 0.00460894051187033\n",
      "5.693110632789214 5.686912943496283 0.006197689292930519\n",
      "6.293486864983308 6.285744198960884 0.00774266602242335\n",
      "6.956612734546251 6.947701383505162 0.008911351041088977\n",
      "7.68867900465214 7.679364009113267 0.009314995538872495\n",
      "8.496556312780378 8.48796326233004 0.008593050450338424\n",
      "9.387978117404057 9.381458908048847 0.006519209355209554\n",
      "10.371739360968139 10.368635469168693 0.0031038917994461457\n",
      "11.457891247681914 11.459218204317192 -0.0013269566352782647\n",
      "12.657916833631345 12.664007768756589 -0.006090935125243391\n",
      "13.984847176088419 13.995030016061225 -0.010182839972806335\n",
      "15.453222546857663 15.4656906518093 -0.01246810495163686\n",
      "17.078805328181264 17.090913936011155 -0.012108607829890872\n",
      "18.878140293987375 18.887242884257727 -0.009102590270352096\n",
      "20.868393862736035 20.872904678409167 -0.0045108156731323845\n",
      "23.0679586056754 23.067896032937664 6.257273773613292e-05\n",
      "25.49744441372899 25.494158143956025 0.00328626977296409\n",
      "28.178750282503444 28.17578450802109 0.0029657744823552434\n",
      "31.126281609320195 31.138819772364624 -0.012538163044428785\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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
