Grand canonical ensemble transition-matrix Monte Carlo

In this example, flat histogram methods are employed for a small macrostate range from 0 to 5 particles. To begin, the system is initialized with the minimum number of particles by setting Metropolis acceptance criteria with favorable conditions for adding particles. The Metropolis criteria are then replaced with the flat histogram criteria. At this point, typical analysis from the previous tutorials are added. In addition, we also add checkpoint files, criteria status, and average energy of a given macrostate. Finally, the simulation is run until the requested number of iterations of the flat histogram algorithm are complete.

A small macrostate range allows the simulation to run quickly with good sampling, and thus it is an ideal starting point to test the simulations. To begin, read the previous SRSW values from file for comparison.

[1]:

import pandas as pd
import numpy as np

srsw_file = "../test/data/stat150.csv"
fh=pd.read_csv(srsw_file)[:6]
fh['lnPI'] -= np.log(sum(np.exp(fh['lnPI'])))   # renormalize
fh['delta_ln_prob']=fh['lnPI'].diff()
fh

[1]:

	N	energy	lnPI	energystd	lnPIstd	delta_ln_prob
0	0	-2.312265e-10	-18.707570	6.689238e-10	0.037092	NaN
1	1	-6.057402e-04	-14.037373	6.709198e-10	0.036838	4.670197
2	2	-3.057422e-02	-10.050312	9.649147e-06	0.036964	3.987061
3	3	-8.992832e-02	-6.458921	1.387472e-04	0.037746	3.591391
4	4	-1.784571e-01	-3.145637	3.315245e-05	0.038097	3.313283
5	5	-2.961920e-01	-0.045677	1.348791e-05	0.038458	3.099960

Run a short flat-histogram simulation to compare with these values

[2]:

script="""
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /feasst/particle/lj.fstprt
Potential Model LennardJones
Potential VisitModel LongRangeCorrections
ThermoParams beta {beta} chemical_potential0 -2.352321
FlatHistogram Macrostate MacrostateNumParticles width 1 min 0 max 5 \
              Bias TransitionMatrix min_sweeps 50
TrialTranslate weight 0.25 tunable_param 1
TrialTransfer weight 1 particle_type 0
Log trials_per_write 1e5 output_file lj.txt
Movie trials_per_write 1e5 output_file lj.xyz
CheckEnergy trials_per_update 1e5 tolerance 1e-8
Tune
CriteriaUpdater trials_per_update 1e5
CriteriaWriter trials_per_write 1e5 output_file lj_fh.txt
Energy output_file lj_en.txt trials_per_update 1 trials_per_write 1e5 multistate true
Checkpoint checkpoint_file checkpoint.fst num_hours 0.0001
Run until_criteria_complete true
""".format(beta=1./1.5)

def run(script):
    with open('script.txt', 'w') as file: file.write(script)
    import subprocess
    syscode = subprocess.call("~/feasst/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.24.1-3-g9284623035-hwh/prefetch_txt2
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/hwh/feasst/particle/lj.fstprt
Potential Model LennardJones
Potential VisitModel LongRangeCorrections
ThermoParams beta 0.6666666666666666 chemical_potential0 -2.352321
FlatHistogram Bias TransitionMatrix Macrostate MacrostateNumParticles max 5 min 0 min_sweeps 50 width 1
TrialTranslate tunable_param 1 weight 0.25
TrialTransfer particle_type 0 weight 1
Log output_file lj.txt trials_per_write 1e5
Movie output_file lj.xyz trials_per_write 1e5
CheckEnergy tolerance 1e-8 trials_per_update 1e5
Tune
CriteriaUpdater trials_per_update 1e5
CriteriaWriter output_file lj_fh.txt trials_per_write 1e5
Energy multistate true output_file lj_en.txt trials_per_update 1 trials_per_write 1e5
Checkpoint checkpoint_file checkpoint.fst num_hours 0.0001
Run until_criteria_complete true
# initializing random number generator with seed: 1701178789

 exit: 0

[3]:

import pandas as pd
import unittest
class TestFlatHistogramLJ(unittest.TestCase):
    """Test flat histogram grand canonical ensemble Monte Carlo simulations"""
    def test_serial_5max(self):
        """Compare the free energies and potential energies with the NIST SRSW
        https://www.nist.gov/programs-projects/nist-standard-reference-simulation-website
        https://mmlapps.nist.gov/srs/LJ_PURE/eostmmc.htm
        """
        fh2=pd.read_csv('lj_fh.txt', comment='#')
        en=pd.read_csv('lj_en.txt')
        print(fh2[['delta_ln_prob', 'delta_ln_prob_stdev']])
        print('energy', en[['average', 'block_stdev']])
        for macro in range(1, len(fh2)):
            self.assertAlmostEqual(
                fh["delta_ln_prob"][macro],
                fh2["delta_ln_prob"][macro],
                delta=3*fh2["delta_ln_prob_stdev"][macro])
            self.assertAlmostEqual(
                fh['energy'][macro],
                en['average'][macro],
                delta=3*(fh["energystd"][macro]**2 + en['block_stdev'][macro]**2)**(1./2.))
unittest.main(argv=[''], verbosity=2, exit=False)

test_serial_5max (__main__.TestFlatHistogramLJ)
Compare the free energies and potential energies with the NIST SRSW ... ok

----------------------------------------------------------------------
Ran 1 test in 0.005s

OK

   delta_ln_prob  delta_ln_prob_stdev
0            NaN             0.000000
1        4.66725             0.002131
2        3.98548             0.001933
3        3.59168             0.002015
4        3.31427             0.001937
5        3.09959             0.002238
energy         average   block_stdev
0  1.016696e-16  3.522090e-16
1 -6.057400e-04  7.702390e-10
2 -3.043054e-02  1.836644e-04
3 -9.030030e-02  4.319611e-04
4 -1.779839e-01  7.443080e-04
5 -2.941329e-01  1.343163e-03

[3]:

<unittest.main.TestProgram at 0x7fd9caf09ae0>

A number of files should also have been created. If the flat histogram method is sampling perfectly, the simulation performs a random walk along the macrostate. For larger ranges of macrostates, or for more difficult sampling cases, monitoring the macrostate can help you determine what conditions are preventing convergence. For example, a plot of the macrostate as a function of the number of attempts may look like the following:

[4]:

pd.read_csv("lj.txt", header=0).dropna(axis='columns')

[4]:

	volume	num_particles_of_type0	beta	state	energy	LennardJones	LongRangeCorrections	trial	TrialTranslate	tunable	TrialAdd	TrialRemove
0	volume	num_particles_of_type0	beta	state	energy	LennardJones	LongRangeCorrections	trial	TrialTranslate	tunable	TrialAdd	TrialRemove
1	512	5	0.666667	5	-0.14133	-0.126186	-0.0151435	100000	0.925813	2.81401	0.0454977	0.0454591
2	512	2	0.666667	2	-0.00242296	9.95731e-16	-0.00242296	200000	0.939471	2.81401	0.360396	0.362235
3	512	0	0.666667	0	-3.77725e-15	1.60982e-15	-2.49366e-18	300000	0.945475	2.81401	0.505605	0.536263
4	512	2	0.666667	2	-0.0158973	-0.0134743	-0.00242296	400000	0.948501	2.81401	0.581509	0.636591
5	512	3	0.666667	3	-0.00545166	4.16334e-17	-0.00545166	500000	0.950833	2.81401	0.627296	0.699383
6	512	5	0.666667	5	-0.176458	-0.161314	-0.0151435	600000	0.952604	2.81401	0.65805	0.742501
7	512	5	0.666667	5	-0.098881	-0.0837375	-0.0151435	700000	0.953478	2.81401	0.679971	0.773808
8	512	4	0.666667	4	-0.0549583	-0.0452665	-0.00969184	800000	0.954208	2.81401	0.695987	0.797581
9	512	0	0.666667	0	1.60646e-15	-2.77556e-15	-1.6263e-18	900000	0.954738	2.81401	0.708714	0.816365
10	512	1	0.666667	1	-0.00060574	1.11022e-16	-0.00060574	1000000	0.955485	2.81401	0.719009	0.831751
11	512	1	0.666667	1	-0.00060574	2.91434e-15	-0.00060574	1100000	0.955818	2.81401	0.727015	0.84423
12	512	0	0.666667	0	-2.07115e-15	-4.996e-16	-2.49366e-18	1200000	0.95631	2.81401	0.734193	0.854884
13	512	5	0.666667	5	-0.0484292	-0.0332857	-0.0151435	1300000	0.956692	2.81401	0.740248	0.863686
14	512	3	0.666667	3	-0.00545166	1.80758e-15	-0.00545166	1400000	0.956989	2.81401	0.745378	0.871351
15	512	5	0.666667	5	-0.163555	-0.148411	-0.0151435	1500000	0.957176	2.81401	0.749456	0.877962
16	512	1	0.666667	1	-0.00060574	2.77556e-15	-0.00060574	1600000	0.957455	2.81401	0.753541	0.883987
17	512	4	0.666667	4	-0.00969184	1.16573e-15	-0.00969184	1700000	0.957672	2.81401	0.756736	0.889133
18	512	3	0.666667	3	-0.0332282	-0.0277765	-0.00545166	1800000	0.957975	2.81401	0.760206	0.893794
19	512	2	0.666667	2	-0.00242296	-2.498e-15	-0.00242296	1900000	0.95787	2.81401	0.762967	0.898006
20	512	1	0.666667	1	-0.00060574	-2.33147e-15	-0.00060574	2000000	0.957844	2.81401	0.765414	0.901732
21	512	3	0.666667	3	-0.00545166	2.9976e-15	-0.00545166	2100000	0.957939	2.81401	0.767779	0.905121
22	512	1	0.666667	1	-0.00060574	-1.44329e-15	-0.00060574	2200000	0.957777	2.81401	0.769498	0.908119
23	512	3	0.666667	3	-0.00545166	1.02696e-15	-0.00545166	2300000	0.957841	2.81401	0.771254	0.910924
24	512	0	0.666667	0	3.15828e-16	-1.60982e-15	-1.6263e-18	2400000	0.957928	2.81401	0.772858	0.913471
25	512	3	0.666667	3	-0.123957	-0.118505	-0.00545166	2500000	0.95803	2.81401	0.774332	0.915807
26	512	5	0.666667	5	-0.183254	-0.168111	-0.0151435	2600000	0.958095	2.81401	0.775958	0.918064
27	512	1	0.666667	1	-0.00060574	6.66134e-16	-0.00060574	2700000	0.958178	2.81401	0.777306	0.920045
28	512	5	0.666667	5	-0.227334	-0.21219	-0.0151435	2800000	0.958262	2.81401	0.778501	0.921877
29	512	5	0.666667	5	-0.217753	-0.20261	-0.0151435	2900000	0.958241	2.81401	0.779472	0.923595
30	512	5	0.666667	5	-0.920985	-0.905841	-0.0151435	3000000	0.958158	2.81401	0.780645	0.92526
31	512	2	0.666667	2	-0.195853	-0.19343	-0.00242296	3100000	0.958214	2.81401	0.781618	0.92677
32	512	3	0.666667	3	-0.0402955	-0.0348438	-0.00545166	3200000	0.95829	2.81401	0.782519	0.928227
33	512	1	0.666667	1	-0.00060574	3.37924e-15	-0.00060574	3300000	0.958397	2.81401	0.783489	0.929628
34	512	1	0.666667	1	-0.00060574	2.38351e-15	-0.00060574	3400000	0.958426	2.81401	0.784288	0.930922
35	512	4	0.666667	4	-0.0321393	-0.0224474	-0.00969184	3500000	0.958485	2.81401	0.78522	0.932092
36	512	2	0.666667	2	-0.00242296	-6.22506e-15	-0.00242296	3600000	0.9585	2.81401	0.785983	0.933251
37	512	4	0.666667	4	-0.0185971	-0.00890523	-0.00969184	3700000	0.958645	2.81401	0.786573	0.934358
38	512	2	0.666667	2	-0.00242296	-4.82253e-15	-0.00242296	3800000	0.958675	2.81401	0.787335	0.935387
39	512	5	0.666667	5	-0.0151435	2.22045e-15	-0.0151435	3900000	0.958693	2.81401	0.787925	0.936339
40	512	4	0.666667	4	-0.0233359	-0.0136441	-0.00969184	4000000	0.958785	2.81401	0.788628	0.937283
41	512	2	0.666667	2	-0.00242296	-3.7817e-15	-0.00242296	4100000	0.958769	2.81401	0.789218	0.93816
42	512	2	0.666667	2	-0.00242296	-5.57887e-15	-0.00242296	4200000	0.958809	2.81401	0.789883	0.938994
43	512	3	0.666667	3	-0.00545166	-1.88738e-15	-0.00545166	4300000	0.958903	2.81401	0.790377	0.939753
44	512	5	0.666667	5	-0.0446818	-0.0295383	-0.0151435	4400000	0.958927	2.81401	0.790762	0.940527
45	512	2	0.666667	2	-0.00838755	-0.00596459	-0.00242296	4500000	0.958949	2.81401	0.791285	0.94124
46	512	3	0.666667	3	-0.00545166	7.77156e-16	-0.00545166	4600000	0.959004	2.81401	0.791751	0.941895
47	512	2	0.666667	2	-0.00242296	6.09235e-15	-0.00242296	4700000	0.959022	2.81401	0.792205	0.942537
48	512	5	0.666667	5	-0.179433	-0.164289	-0.0151435	4800000	0.958979	2.81401	0.792588	0.943169
49	512	1	0.666667	1	-0.00060574	-5.33254e-15	-0.00060574	4900000	0.959012	2.81401	0.792982	0.943757
50	512	3	0.666667	3	-0.0200995	-0.0146478	-0.00545166	5000000	0.959074	2.81401	0.793388	0.944371
51	512	0	0.666667	0	4.38375e-15	1.80411e-15	-1.6263e-18	5100000	0.959078	2.81401	0.793606	0.944915
52	512	0	0.666667	0	-1.59486e-16	2.1155e-15	-2.49366e-18	5200000	0.959151	2.81401	0.793853	0.945424

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