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.

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

	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

script="""
MonteCarlo
Configuration cubic_box_length 8 particle_type0 /feasst/forcefield/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 1e5 file_name lj.txt
Movie trials_per 1e5 file_name lj.xyz
CheckEnergy trials_per 1e5 tolerance 1e-8
Tune
CriteriaUpdater trials_per 1e5
CriteriaWriter trials_per 1e5 file_name lj_fh.txt
Energy file_name lj_en.txt trials_per_update 1 trials_per_write 1e5 multistate true
Checkpoint file_name checkpoint_file_name 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.19.0-59-gb2fd5186f4-dirty-hwh/newtutorials
Configuration cubic_box_length 8 particle_type0 /home/hwh/feasst/forcefield/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 file_name lj.txt trials_per 1e5
Movie file_name lj.xyz trials_per 1e5
CheckEnergy tolerance 1e-8 trials_per 1e5
Tune
CriteriaUpdater trials_per 1e5
CriteriaWriter file_name lj_fh.txt trials_per 1e5
Energy file_name lj_en.txt multistate true trials_per_update 1 trials_per_write 1e5
Checkpoint file_name checkpoint_file_name num_hours 0.0001
Run until_criteria_complete true
# initializing random number generator with seed: 1652468091

 exit: 0

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 ...

   delta_ln_prob  delta_ln_prob_stdev
0            NaN             0.000000
1        4.67154             0.001700
2        3.98743             0.001968
3        3.59350             0.002018
4        3.31255             0.001676
5        3.10120             0.002028
energy         average   block_stdev
0 -1.539438e-16  2.872237e-16
1 -6.057400e-04  7.109906e-10
2 -3.064926e-02  2.102718e-04
3 -8.984079e-02  4.283351e-04
4 -1.785064e-01  7.480973e-04
5 -2.973343e-01  1.535943e-03

ok

----------------------------------------------------------------------
Ran 1 test in 0.012s

OK

<unittest.main.TestProgram at 0x7f72f16c2310>

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:

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

	volume	p0	beta	state	energy	LennardJones	LongRangeCorrections	attempt	TrialTranslate	tunable	TrialAdd	TrialRemove
0	volume	p0	beta	state	energy	LennardJones	LongRangeCorrections	attempt	TrialTranslate	tunable	TrialAdd	TrialRemove
1	512	5	0.666667	5	-0.923994	-0.90885	-0.0151435	100000	0.925429	2.85092	0.0431069	0.0428732
2	512	2	0.666667	2	-0.00242296	-1.22125e-15	-0.00242296	200000	0.943075	2.85092	0.377961	0.416915
3	512	5	0.666667	5	-0.299839	-0.284696	-0.0151435	300000	0.948685	2.85092	0.52155	0.591715
4	512	3	0.666667	3	-0.307744	-0.302292	-0.00545166	400000	0.952078	2.85092	0.593087	0.683166
...	...	...	...	...	...	...	...	...	...	...	...	...
151	512	2	0.666667	2	-0.0161133	-0.0136903	-0.00242296	4700000	0.959029	2.83079	0.78844	0.939091
152	512	2	0.666667	2	-0.00242296	1.08247e-15	-0.00242296	4800000	0.959033	2.83079	0.788803	0.939807
153	512	3	0.666667	3	-0.00545166	-2.56739e-15	-0.00545166	4900000	0.95903	2.83079	0.789328	0.940546
154	512	3	0.666667	3	-0.00545166	-4.38538e-15	-0.00545166	5000000	0.959004	2.83079	0.789747	0.941208
155	512	4	0.666667	4	-0.0780117	-0.0683199	-0.00969184	5100000	0.959087	2.83079	0.790239	0.941852

156 rows × 12 columns

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