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=atom:/feasst/particle/atom.txt
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=atom
CriteriaUpdater trials_per_update=1e5
CriteriaWriter trials_per_write=1e5 output_file=id_fh.txt
Run until=complete
""".format(**params)

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

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration cubic_side_length=8 particle_type=atom:/feasst/particle/atom.txt
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=atom
CriteriaUpdater trials_per_update=1e5
CriteriaWriter output_file=id_fh.txt trials_per_write=1e5
Run until=complete
# Initializing random number generator with seed: 1748436045

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

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