Ideal gas equation of state using grand canonical ensemble transition-matrix Monte Carlo

In this example, the ideal gas equation of state is obtained as a test of the flat histogram method.

params={"cubic_side_length": 8, "beta": 1./1.2, "mu": -3}
script="""
MonteCarlo
Configuration cubic_side_length {cubic_side_length} particle_type /feasst/particle/atom.fstprt
Potential Model IdealGas
ThermoParams beta {beta} chemical_potential {mu}
FlatHistogram Macrostate MacrostateNumParticles width 1 min 0 max 50 \
              Bias TransitionMatrix min_sweeps 100
TrialTransfer particle_type 0
CriteriaUpdater trials_per_update 1e5
CriteriaWriter trials_per_write 1e5 output_file id_fh.txt
Run until_criteria_complete true
""".format(**params)

def run(script):
    with open('script.txt', 'w') as file: file.write(script)
    import subprocess
    syscode = subprocess.call("../../../build/bin/fst < script.txt > script.log", shell=True, executable='/bin/bash')
    with open('script.log', 'r') as file: print(file.read(), '\n', 'exit:', syscode)
run(script)

# FEASST version: v0.20.0-4-g08d5a7de4e-user/pyfeasst
Configuration cubic_side_length 8 particle_type /home/user/feasst/particle/atom.fstprt
Potential Model IdealGas
ThermoParams beta 0.8333333333333334 chemical_potential -3
FlatHistogram Bias TransitionMatrix Macrostate MacrostateNumParticles max 50 min 0 min_sweeps 100 width 1
TrialTransfer particle_type 0
CriteriaUpdater trials_per 1e5
CriteriaWriter file_name id_fh.txt trials_per 1e5
Run until_criteria_complete true
# initializing random number generator with seed: 1655833271

 exit: 0

Check the ideal gas relationship, \(\beta P = \rho\)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

fh=pd.read_csv('id_fh.txt', comment="#")
#print(fh)
print('N betaPV difference')
for delta_conjugate in np.arange(-5, 5, 0.1):
    fh['ln_prob_rw'] = fh['ln_prob'] + fh['state']*delta_conjugate - fh['ln_prob'].min()  # avoid negative log
    fh['ln_prob_rw'] -= np.log(sum(np.exp(fh['ln_prob_rw'])))   # renormalize
    if fh['ln_prob_rw'].values[-1] < -6:
        #plt.plot(fh['ln_prob_rw'])
        N = (np.exp(fh["ln_prob_rw"]) * fh["state"]).sum()
        betaPV = -fh["ln_prob_rw"][0] - np.log(np.exp(fh["ln_prob_rw"]).sum())
        print(N, betaPV, N-betaPV)
        assert np.abs(1 - betaPV/N) < 1e-2

N betaPV difference
0.28260651938960046 0.28243688472004025 0.00016963466956021467
0.3123640140080462 0.312160569899379 0.00020344410866718343
0.3452576004050282 0.3450141889097207 0.0002434114953074884
0.3816180498224627 0.3813276134222697 0.0002904364001930504
0.4218109655460796 0.421465505072794 0.00034546047328559704
0.46624040381102455 0.4658309768551357 0.0004094269558888586
0.5153528518422282 0.5148696342203101 0.00048321762191816386
0.5696415928809494 0.5690740330727739 0.0005675598081754663
0.6296514900217144 0.6289885949347819 0.0006628950869325134
0.6959842234510876 0.6952150228614556 0.0007692005896319243
0.7693040199728768 0.7684182653461628 0.0008857546267140304
0.8503439206649607 0.8493330796640493 0.0010108410009114221
0.9399126438055039 0.9387712512061481 0.0011413925993557505
1.0389021180881615 1.0376295318984148 0.0012725861897466828
1.1482957884091773 1.1468983695744732 0.0013974188347041228
1.269177836159279 1.2676715122664641 0.0015063238928147982
1.4027435104269264 1.4011565881655705 0.0015869222613558165
1.5503108362055416 1.5486867850098522 0.0016240511956893133
1.7133340465602604 1.7117337832504989 0.001600263309761596
1.8934191660330588 1.8919221361398764 0.0014970298931824555
2.09234222935907 2.091045335779736 0.0012968935793340464
2.312070616138792 2.3110838530945963 0.0009867630441955377
2.554787869609652 2.554225483266434 0.0005623863432178133
2.822922094598291 2.8228883528533832 3.374174490788917e-05
3.1191775633209327 3.119746932642518 -0.0005693693215853557
3.4465685173845846 3.447761332294398 -0.0011928149098134888
3.808453449588775 3.810210016221623 -0.0017565666328476226
4.208567601348703 4.210725878602589 -0.002158277253886176
4.651051332028505 4.653335379732325 -0.002284047703819958
5.140472721569988 5.142500239164774 -0.0020275175947865876
5.681844447734084 5.683161092907086 -0.0013166451730013407
6.280637572452482 6.2807826492018295 -0.00014507674934716164
6.94279797735872 6.941400294607094 0.0013976827516257728
7.6747739799675525 7.671668823014407 0.0031051569531452827
8.4835648776486 8.478914893035899 0.0046499846127012745
9.37679881587365 9.371195747576984 0.005603068296666791
10.362846015682903 10.357367440547035 0.0054785751358679136
11.45097620934702 11.447166482505988 0.0038097268410322016
12.651582969631738 12.651310305232823 0.0002726643989152677
13.976513162600328 13.981625078653195 -0.005111916052866761
15.43952872436892 15.45121302699289 -0.011684302623971021
17.056871909304128 17.07467169000677 -0.01779978070264221
18.8477860252465 18.86836956370471 -0.020583538458211592
20.834521436787426 20.850754288226288 -0.016232851438861218
23.040776499430144 23.042593765009727 -0.0018172655795822834
25.487818908344565 25.466943711260193 0.020875197084372132
28.190314554370506 28.148681958503055 0.04163259586745127
31.153222301399527 31.113712284208187 0.03951001719133984

Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!