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_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
[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 ...
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
[3]:
<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:
[4]:
pd.read_csv("lj.txt", header=0).dropna(axis='columns')
[4]:
| 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!