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.

[1]:
import feasst as fst

monte_carlo = fst.MonteCarlo()
monte_carlo.add(fst.Configuration(fst.MakeDomain(fst.args({"cubic_box_length": "8"})),
                                  fst.args({"particle_type": fst.install_dir()+'/forcefield/atom.fstprt'})))
monte_carlo.add(fst.MakePotential(fst.MakeDontVisitModel()))
monte_carlo.set(fst.MakeThermoParams(fst.args({"beta": str(1./1.2), "chemical_potential": "-3"})))
criteria = fst.MakeFlatHistogram(
    fst.MakeMacrostateNumParticles(fst.Histogram(fst.args({"width": "1", "min": "0", "max": "50"}))),
    fst.MakeTransitionMatrix(fst.args({"min_sweeps": "100"})))
monte_carlo.set(criteria)
monte_carlo.add(fst.MakeTrialTransfer(fst.args({"particle_type": "0"})))
monte_carlo.add(fst.MakeCriteriaUpdater(fst.args({"steps_per": str(1e5)})))
monte_carlo.add(fst.MakeCriteriaWriter(fst.args({"steps_per": str(1e5), "file_name": "id_fh.txt"})))
[2]:
%time
monte_carlo.run_until_complete()
CPU times: user 13.7 s, sys: 32.6 ms, total: 13.7 s
Wall time: 13.7 s

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

[3]:
import numpy as np

volume = monte_carlo.configuration().domain().volume()
gce = fst.GrandCanonicalEnsemble(criteria, monte_carlo.thermo_params().beta_mu())
print('rho beta*P difference')
for delta_conjugate in np.arange(-5, 5, 0.1):
        lnpirw = gce.reweight(delta_conjugate)
        if lnpirw.value(lnpirw.size()-1) < -6:
            betaPV = gce.betaPV()
            N = gce.average_macrostate()
            print(N, betaPV, np.abs(1 - betaPV/N))
            assert(np.abs(1 - betaPV/N) < 1e-2)
rho beta*P difference
0.28352371179938896 0.2834109250174366 0.00039780370127273645
0.3133659767306003 0.3132305130399782 0.0004322858915171146
0.3463511041451104 0.3461888474487959 0.00046847460387045015
0.38281022254492336 0.3826164949224828 0.0005060670040436532
0.4231092809374098 0.4228788472982002 0.0005446196753213783
0.46765266421644996 0.46737978153695275 0.0005835157166365068
0.5168871656966174 0.516565698838592 0.0006219284968937977
0.5713063473665461 0.570929980132204 0.0006587835687051902
0.631455320931917 0.6310178983430001 0.000692721360351034
0.6979359862846308 0.697432031310502 0.0007220647509688094
0.771412769454396 0.7708382231478559 0.0007447975056810163
0.8526189105229772 0.85197214642801 0.0007585617524838995
0.9423633649624151 0.9416465232340795 0.0007606850552431821
1.0415384013082887 1.0407590703673615 0.0007482498388424741
1.151128006056329 1.1503012436171116 0.0007182193768786593
1.2722172448799705 1.2713678689055983 0.0006676343822491493
1.4060027781173632 1.4051677653898988 0.00059389123582132
1.553804785548966 1.553035488203192 0.0004951055325153053
1.7170806133044854 1.7164443469548176 0.00037055123955032077
1.8974404989316966 1.8970208897643575 0.0002211448356749246
2.096665732999686 2.0965610788082754 4.991458092895584e-05
2.3167295381427033 2.317048416181871 0.00013764146134342958
2.5598207392917223 2.560674297916247 0.00033344468674045835
2.8283699111229637 2.829860863660274 0.0005271419878449812
3.1250770935657326 3.1272865500817195 0.0007070086432543476
3.4529394082330316 3.4559144280807605 0.0008615905163686932
3.8152761774003796 3.8190231990755215 0.0009821102066835952
4.215748861013312 4.220240465294717 0.0010654345003666688
4.658373975902694 4.663577650034231 0.0011170580461026258
5.147530006393303 5.153465882359845 0.001153150337962039
5.687964826353012 5.69479251384809 0.0012003744227540025
6.284818132906159 6.29293897008084 0.001292135588166321
6.9436819457861985 6.95382252632704 0.0014604039499528376
7.670727114890588 7.68394720157967 0.0017234463553552626
8.472918308514009 8.490471605695918 0.0020716943729139903
9.358315789513808 9.381302804120729 0.0024563196117701747
10.336412753898973 10.365222839615774 0.0027872421896013577
11.41838675833331 11.452046007393868 0.0029478112602896456
12.617086537971835 12.652789611353448 0.002829739914532725
13.946605819631271 13.97982345527495 0.002381772029214524
15.421495767263217 15.446956881952453 0.0016510146015333849
17.056018442469547 17.06944424040934 0.0007871589717773908
18.86409972506837 18.863945156385697 8.19380118466384e-06
20.86039343798047 20.84853512853993 0.0005684604883314792
23.062002864519698 23.042863474283944 0.0008299101490963912
25.489623715743093 25.468465733436904 0.0008300625596572297
28.166798365996392 28.14910842417289 0.0006280423352927889
31.112973356044062 31.110824232384275 6.907484010587606e-05

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