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("feasst < 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: feasst < file.txt
# For more information, use the command "feasst-menu"
# Exit with ctrl-c
FEASST version 0.25.19
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: 1781198252

 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.2829214450028523 0.28276142935893817 0.00016001564391415757
0.3127110344609224 0.3125181844159339 0.00019285004498853509
0.3456402287398188 0.34540825464818914 0.00023197409162967197
0.3820402993452214 0.3817618872672399 0.0002784120779814603
0.4222774542216826 0.421944174504675 0.0003332797170076396
0.46675648379027573 0.4663587266639026 0.00039775712637313276
0.5159247687089498 0.5154517276198478 0.0004730410891019554
0.570276678959833 0.5697164104095647 0.0005602685502683613
0.6303583950694145 0.6296979935758966 0.0006604014935178704
0.696773183838372 0.695999122080295 0.0007740617580769271
0.770187163519557 0.7692858599548528 0.0009013035647041923
0.8513355980068428 0.8502942855654231 0.001041312441419695
0.9410297679934123 0.9398377446931344 0.0011920233002779002
1.0401644817643028 1.0388148221072249 0.0013496596570778951
1.1497263127984105 1.14821809969524 0.0015082131031705082
1.2708026900299596 1.2691437797302167 0.001658910299742855
1.4045920241185472 1.4028022671395293 0.0017897569790179535
1.552415132895656 1.5505298267730436 0.0018853061226122847
1.715728330676701 1.7138014629021896 0.0019268677745114537
1.8961386595276704 1.8942452102998013 0.0018934492278690485
2.0954218391723356 2.093658079286071 0.0017637598862645376
2.315543542873657 2.314023957202639 0.0015195856710179534
2.558684481248569 2.5575338250073445 0.001150656241224457
2.8272693721483746 2.826608678578439 0.0006606935699355887
3.123999060529628 3.1239255153079615 7.354522166647115e-05
3.4518837585846183 3.452446612157545 -0.0005628535729265138
3.8142737410538046 3.8154520379532544 -0.0011782968994498155
4.214882414380421 4.216574899646911 -0.0016924852664903511
4.6577966791061085 4.659838295506624 -0.0020416164005157
5.147472750233575 5.149692566843759 -0.002219816610184111
5.688723927088134 5.691051625222572 -0.0023276981344384495
6.286720429404237 6.289328426302422 -0.0026079968981846946
6.947036426917338 6.950472444007084 -0.0034360170897462794
7.675785778802292 7.681015960238509 -0.005230181436216341
8.47987119024142 8.48813955101325 -0.008268360771829464
9.367320311579844 9.379767360418134 -0.012447048838289732
10.347603415096874 10.364696286167293 -0.017092871070419235
11.431761006992256 11.452748968406713 -0.020987961414457246
12.63219441581632 12.654923396269563 -0.022728980453242897
13.962166478423628 13.98350540546874 -0.021338927045112754
15.435377036087809 15.452130769149123 -0.016753733061314335
17.066095896279855 17.075829672218084 -0.009733775938229883
18.869827844735347 18.871115140125454 -0.0012872953901066353
20.86363652307769 20.8561301260069 0.007506397070791593
23.065498198040224 23.05077730812369 0.01472088991653564
25.493907587206714 25.47677929014812 0.0171282970585942
28.16900640379734 28.15778058598642 0.011225817810920802
31.109107602401853 31.11944202979696 -0.010334427395108037

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