Grand canonical ensemble transition-matrix Monte Carlo

In this example, FlatHistogram 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 FlatHistogram. Typical analysis from the previous tutorials are added, and we add Checkpoint files, CriteriaUpdater, CriteriaWriter, and average Energy of a given Macrostate. Finally, the simulation is run until the requested convergence criteria for the FlatHistogram algorithm are complete.

A small Macrostate range allows the simulation to run quickly with ample sampling, and thus it is an ideal starting point to test the simulations. To begin, read the previous NIST SRSW values for LJ 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 FlatHistogram simulation to compare with these NIST SRSW values.

[2]:
script="""
MonteCarlo
Configuration cubic_side_length=8 particle_type=lj:/feasst/particle/lj_new.txt
Potential Model=LennardJones
Potential VisitModel=LongRangeCorrections
ThermoParams beta={beta} chemical_potential=-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=lj
Let [write]=trials_per_write=1e5 output_file=lj
Log [write].txt
Movie [write].xyz
CheckEnergy trials_per_update=1e5 decimal_places=8
Tune
CriteriaUpdater trials_per_update=1e5
CriteriaWriter [write]_fh.txt
Energy [write]_en.txt multistate=true
Checkpoint checkpoint_file=checkpoint.fst num_hours=0.0001
Run until=complete
""".format(beta=1./1.5)

def run(script):
    with open('script1.txt', 'w') as file: file.write(script)
    import subprocess
    syscode = subprocess.call("../../../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)
# Usage: /Users/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration cubic_side_length=8 particle_type=lj:/feasst/particle/lj_new.txt
Potential Model=LennardJones
Potential VisitModel=LongRangeCorrections
ThermoParams beta=0.6666666666666666 chemical_potential=-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=lj weight=1
Log output_file=lj.txt trials_per_write=1e5
Movie output_file=lj.xyz trials_per_write=1e5
CheckEnergy decimal_places=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_write=1e5
Checkpoint checkpoint_file=checkpoint.fst num_hours=0.0001
Run until=complete
# Initializing random number generator with seed: 1749408693

 exit: 0
[3]:
import pandas as pd

"""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)):
    assert np.abs(fh["delta_ln_prob"][macro] - fh2["delta_ln_prob"][macro]) < 4*fh2["delta_ln_prob_stdev"][macro]
    assert np.abs(fh['energy'][macro] - en['average'][macro]) < 4*(fh["energystd"][macro]**2 + en['block_stdev'][macro]**2)**(1./2.)
   delta_ln_prob  delta_ln_prob_stdev
0            NaN             0.000000
1        4.66914             0.001645
2        3.98803             0.002257
3        3.59332             0.001912
4        3.31309             0.001980
5        3.10125             0.001972
energy         average   block_stdev
0  1.317466e-16  1.953789e-16
1 -6.057400e-04  7.702399e-10
2 -3.047790e-02  2.260309e-04
3 -9.005908e-02  4.629765e-04
4 -1.788966e-01  7.342915e-04
5 -2.965805e-01  1.231712e-03

A number of files should also have been created. If the FlatHistogram 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, 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_lj beta state energy LennardJones LongRangeCorrections BondTwoBody BondThreeBody BondFourBody trial TrialTranslate tunable TrialAdd TrialRemove
0 512 5 0.666667 5 -5.615610e-02 -4.101260e-02 -1.514350e-02 0 0 0 100000 0.927751 2.83079 0.047373 0.046741
1 512 0 0.666667 0 -3.669700e-15 -3.330670e-16 -1.517880e-18 0 0 0 200000 0.944879 2.83079 0.256928 0.346506
2 512 2 0.666667 2 -2.422960e-03 1.682680e-16 -2.422960e-03 0 0 0 300000 0.950235 2.83079 0.441354 0.570366
3 512 4 0.666667 4 1.532020e-01 1.628940e-01 -9.691840e-03 0 0 0 400000 0.953083 2.83079 0.534310 0.676512
4 512 4 0.666667 4 -9.691840e-03 3.108620e-15 -9.691840e-03 0 0 0 500000 0.954784 2.83079 0.589223 0.737878
5 512 1 0.666667 1 -6.057400e-04 -3.747000e-15 -6.057400e-04 0 0 0 600000 0.956285 2.83079 0.625951 0.778506
6 512 1 0.666667 1 -6.057400e-04 1.304510e-15 -6.057400e-04 0 0 0 700000 0.957065 2.83079 0.652960 0.807427
7 512 0 0.666667 0 -1.987020e-15 -1.935950e-15 -1.517880e-18 0 0 0 800000 0.957227 2.83079 0.672873 0.828548
8 512 5 0.666667 5 -3.277240e-01 -3.125800e-01 -1.514350e-02 0 0 0 900000 0.957637 2.83079 0.688387 0.844791
9 512 1 0.666667 1 -6.057400e-04 5.630910e-15 -6.057400e-04 0 0 0 1000000 0.957736 2.83079 0.700375 0.857707
10 512 4 0.666667 4 -4.035140e-02 -3.065960e-02 -9.691840e-03 0 0 0 1100000 0.957840 2.83079 0.709995 0.868286
11 512 2 0.666667 2 -1.205740e-01 -1.181510e-01 -2.422960e-03 0 0 0 1200000 0.958064 2.83079 0.718449 0.877129
12 512 3 0.666667 3 -3.439830e-02 -2.894660e-02 -5.451660e-03 0 0 0 1300000 0.958033 2.83079 0.725688 0.884632
13 512 4 0.666667 4 -2.340630e-02 -1.371450e-02 -9.691840e-03 0 0 0 1400000 0.958037 2.83079 0.731605 0.891108
14 512 0 0.666667 0 -4.331390e-16 -2.373100e-15 -2.493660e-18 0 0 0 1500000 0.958290 2.83079 0.736525 0.896649
15 512 3 0.666667 3 -2.461910e-02 -1.916740e-02 -5.451660e-03 0 0 0 1600000 0.958445 2.83079 0.741085 0.901540
16 512 5 0.666667 5 -1.807250e-01 -1.655820e-01 -1.514350e-02 0 0 0 1700000 0.958455 2.83079 0.745267 0.905789
17 512 1 0.666667 1 -6.057400e-04 3.070460e-15 -6.057400e-04 0 0 0 1800000 0.958639 2.83079 0.748933 0.909629
18 512 0 0.666667 0 1.944730e-15 1.887380e-15 -1.626300e-18 0 0 0 1900000 0.958730 2.83079 0.752052 0.912985
19 512 1 0.666667 1 -6.057400e-04 -4.413140e-15 -6.057400e-04 0 0 0 2000000 0.958743 2.83079 0.755071 0.916069
20 512 1 0.666667 1 -6.057400e-04 2.248200e-15 -6.057400e-04 0 0 0 2100000 0.958821 2.83079 0.757694 0.918892
21 512 5 0.666667 5 -1.044350e-01 -8.929170e-02 -1.514350e-02 0 0 0 2200000 0.958917 2.83079 0.760386 0.921492
22 512 3 0.666667 3 -3.782400e-02 -3.237240e-02 -5.451660e-03 0 0 0 2300000 0.959054 2.83079 0.762617 0.923845
23 512 5 0.666667 5 -1.086300e-01 -9.348650e-02 -1.514350e-02 0 0 0 2400000 0.958905 2.83079 0.764488 0.925978
24 512 1 0.666667 1 -6.057400e-04 -1.942890e-16 -6.057400e-04 0 0 0 2500000 0.958914 2.83079 0.766521 0.927944
25 512 2 0.666667 2 -4.778790e-02 -4.536500e-02 -2.422960e-03 0 0 0 2600000 0.958939 2.83079 0.768050 0.929698
26 512 2 0.666667 2 -9.900980e-01 -9.876750e-01 -2.422960e-03 0 0 0 2700000 0.959031 2.83079 0.769619 0.931341
27 512 4 0.666667 4 -1.657710e-02 -6.885270e-03 -9.691840e-03 0 0 0 2800000 0.959084 2.83079 0.770725 0.932810
28 512 1 0.666667 1 -6.057400e-04 -6.938890e-18 -6.057400e-04 0 0 0 2900000 0.959179 2.83079 0.772262 0.934336
29 512 3 0.666667 3 -5.451660e-03 -3.025360e-15 -5.451660e-03 0 0 0 3000000 0.959204 2.83079 0.773755 0.935689
30 512 2 0.666667 2 -2.422960e-03 -5.773160e-15 -2.422960e-03 0 0 0 3100000 0.959197 2.83079 0.774697 0.936865
31 512 1 0.666667 1 -6.057400e-04 -3.400060e-16 -6.057400e-04 0 0 0 3200000 0.959207 2.83079 0.775970 0.938091
32 512 0 0.666667 0 7.018800e-15 3.663740e-15 -2.385240e-18 0 0 0 3300000 0.959206 2.83079 0.777194 0.939220
33 512 4 0.666667 4 -8.814790e-01 -8.717870e-01 -9.691840e-03 0 0 0 3400000 0.959216 2.83079 0.778111 0.940227
34 512 2 0.666667 2 -2.422960e-03 -1.734720e-18 -2.422960e-03 0 0 0 3500000 0.959218 2.83079 0.778960 0.941170
35 512 0 0.666667 0 -1.876860e-15 3.663740e-15 -1.517880e-18 0 0 0 3600000 0.959229 2.83079 0.779697 0.942059
36 512 1 0.666667 1 -6.057400e-04 -1.075530e-15 -6.057400e-04 0 0 0 3700000 0.959234 2.83079 0.780548 0.942917
37 512 4 0.666667 4 -9.224960e-02 -8.255780e-02 -9.691840e-03 0 0 0 3800000 0.959233 2.83079 0.781348 0.943726
38 512 0 0.666667 0 -1.333570e-17 -2.550040e-15 -2.493660e-18 0 0 0 3900000 0.959321 2.83079 0.782267 0.944475
39 512 0 0.666667 0 4.490010e-15 -2.248200e-15 -1.517880e-18 0 0 0 4000000 0.959337 2.83079 0.782988 0.945215
40 512 4 0.666667 4 -1.042640e+00 -1.032940e+00 -9.691840e-03 0 0 0 4100000 0.959398 2.83079 0.783727 0.945905
41 512 0 0.666667 0 5.485300e-15 3.486790e-15 -2.493660e-18 0 0 0 4200000 0.959446 2.83079 0.784386 0.946561
42 512 1 0.666667 1 -6.057400e-04 -2.026160e-15 -6.057400e-04 0 0 0 4300000 0.959459 2.83079 0.784971 0.947184
43 512 3 0.666667 3 -5.451660e-03 7.261550e-15 -5.451660e-03 0 0 0 4400000 0.959442 2.83079 0.785598 0.947789
44 512 2 0.666667 2 -2.422960e-03 -8.049120e-16 -2.422960e-03 0 0 0 4500000 0.959441 2.83079 0.786153 0.948361
45 512 0 0.666667 0 -3.111660e-17 -2.220450e-15 -2.493660e-18 0 0 0 4600000 0.959466 2.83079 0.786648 0.948910
46 512 1 0.666667 1 -6.057400e-04 -3.441690e-15 -6.057400e-04 0 0 0 4700000 0.959506 2.83079 0.787197 0.949461
47 512 1 0.666667 1 -6.057400e-04 9.364040e-15 -6.057400e-04 0 0 0 4800000 0.959587 2.83079 0.787692 0.949976
48 512 2 0.666667 2 -2.422960e-03 -5.273560e-16 -2.422960e-03 0 0 0 4900000 0.959619 2.83079 0.788031 0.950492
49 512 0 0.666667 0 -1.609280e-15 4.440890e-16 -1.517880e-18 0 0 0 5000000 0.959615 2.83079 0.788432 0.950945
50 512 2 0.666667 2 -2.422960e-03 9.020560e-16 -2.422960e-03 0 0 0 5100000 0.959605 2.83079 0.788997 0.951434
51 512 2 0.666667 2 -2.422960e-03 0.000000e+00 -2.422960e-03 0 0 0 5100000 0.959605 2.83079 0.788997 0.951434

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