# 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['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('script1.txt', 'w') as file: file.write(script)
import subprocess
syscode = subprocess.call("~/feasst/build/bin/fst < script1.txt > script1.log", shell=True, executable='/bin/bash')
with open('script1.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
"""
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