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=6
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("feasst < 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: feasst < file.txt
# For more information, use the command "feasst-menu"
# Exit with ctrl-c
FEASST version 0.25.19
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=6 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: 1781198285
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.66719 0.001785
2 3.98682 0.002033
3 3.59093 0.001868
4 3.31297 0.002114
5 3.10353 0.002382
energy average block_stdev
0 1.455599e-17 3.578654e-16
1 -6.057400e-04 7.702409e-10
2 -3.067295e-02 2.366927e-04
3 -9.009797e-02 3.669175e-04
4 -1.794613e-01 5.969186e-04
5 -2.971404e-01 1.503380e-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 | 4 | 0.666667 | 4 | -6.574620e-02 | -5.605440e-02 | -9.691840e-03 | 0 | 0 | 0 | 100000 | 0.924354 | 2.86267 | 0.041399 | 0.040860 |
| 1 | 512 | 3 | 0.666667 | 3 | -1.149280e-01 | -1.094770e-01 | -5.451660e-03 | 0 | 0 | 0 | 200000 | 0.938366 | 2.86267 | 0.423966 | 0.454913 |
| 2 | 512 | 4 | 0.666667 | 4 | -2.291060e-02 | -1.321870e-02 | -9.691840e-03 | 0 | 0 | 0 | 300000 | 0.944099 | 2.86267 | 0.552455 | 0.615041 |
| 3 | 512 | 1 | 0.666667 | 1 | -6.057400e-04 | 3.996800e-15 | -6.057400e-04 | 0 | 0 | 0 | 400000 | 0.948031 | 2.86267 | 0.618505 | 0.699765 |
| 4 | 512 | 4 | 0.666667 | 4 | -5.909450e-01 | -5.812530e-01 | -9.691840e-03 | 0 | 0 | 0 | 500000 | 0.950384 | 2.86267 | 0.657824 | 0.751987 |
| 5 | 512 | 1 | 0.666667 | 1 | -6.057400e-04 | 4.440890e-16 | -6.057400e-04 | 0 | 0 | 0 | 600000 | 0.951794 | 2.86267 | 0.683809 | 0.787278 |
| 6 | 512 | 1 | 0.666667 | 1 | -6.057400e-04 | 9.992010e-16 | -6.057400e-04 | 0 | 0 | 0 | 700000 | 0.952943 | 2.86267 | 0.701278 | 0.812808 |
| 7 | 512 | 3 | 0.666667 | 3 | -5.451660e-03 | 6.361230e-15 | -5.451660e-03 | 0 | 0 | 0 | 800000 | 0.953478 | 2.86267 | 0.714862 | 0.832068 |
| 8 | 512 | 4 | 0.666667 | 4 | -1.211070e-01 | -1.114150e-01 | -9.691840e-03 | 0 | 0 | 0 | 900000 | 0.954305 | 2.86267 | 0.725634 | 0.847293 |
| 9 | 512 | 5 | 0.666667 | 5 | -6.965650e-02 | -5.451300e-02 | -1.514350e-02 | 0 | 0 | 0 | 1000000 | 0.954824 | 2.86267 | 0.733594 | 0.859479 |
| 10 | 512 | 4 | 0.666667 | 4 | -6.309660e-02 | -5.340470e-02 | -9.691840e-03 | 0 | 0 | 0 | 1100000 | 0.955037 | 2.86267 | 0.740955 | 0.869446 |
| 11 | 512 | 0 | 0.666667 | 0 | -6.478110e-16 | 3.747000e-15 | -1.626300e-18 | 0 | 0 | 0 | 1200000 | 0.955470 | 2.86267 | 0.746228 | 0.877741 |
| 12 | 512 | 5 | 0.666667 | 5 | -1.514350e-02 | -1.901260e-15 | -1.514350e-02 | 0 | 0 | 0 | 1300000 | 0.955691 | 2.86267 | 0.751517 | 0.884949 |
| 13 | 512 | 2 | 0.666667 | 2 | -2.422960e-03 | -1.623700e-15 | -2.422960e-03 | 0 | 0 | 0 | 1400000 | 0.955936 | 2.86267 | 0.755277 | 0.891036 |
| 14 | 512 | 0 | 0.666667 | 0 | 2.439450e-17 | -1.651460e-15 | -1.626300e-18 | 0 | 0 | 0 | 1500000 | 0.956335 | 2.86267 | 0.759130 | 0.896458 |
| 15 | 512 | 5 | 0.666667 | 5 | 1.597370e-01 | 1.748800e-01 | -1.514350e-02 | 0 | 0 | 0 | 1600000 | 0.956408 | 2.86267 | 0.762546 | 0.901215 |
| 16 | 512 | 4 | 0.666667 | 4 | -3.710930e-01 | -3.614010e-01 | -9.691840e-03 | 0 | 0 | 0 | 1700000 | 0.956654 | 2.86267 | 0.765445 | 0.905320 |
| 17 | 512 | 3 | 0.666667 | 3 | -5.451660e-03 | 2.220450e-16 | -5.451660e-03 | 0 | 0 | 0 | 1800000 | 0.956853 | 2.86267 | 0.768070 | 0.909066 |
| 18 | 512 | 5 | 0.666667 | 5 | -2.790970e-01 | -2.639530e-01 | -1.514350e-02 | 0 | 0 | 0 | 1900000 | 0.956900 | 2.86267 | 0.770122 | 0.912469 |
| 19 | 512 | 4 | 0.666667 | 4 | -9.691840e-03 | -1.044300e-14 | -9.691840e-03 | 0 | 0 | 0 | 2000000 | 0.956982 | 2.86267 | 0.772176 | 0.915463 |
| 20 | 512 | 1 | 0.666667 | 1 | -6.057400e-04 | -4.718450e-15 | -6.057400e-04 | 0 | 0 | 0 | 2100000 | 0.956910 | 2.86267 | 0.773766 | 0.918128 |
| 21 | 512 | 4 | 0.666667 | 4 | -1.080630e+00 | -1.070930e+00 | -9.691840e-03 | 0 | 0 | 0 | 2200000 | 0.957016 | 2.86267 | 0.775535 | 0.920608 |
| 22 | 512 | 1 | 0.666667 | 1 | -6.057400e-04 | -2.102480e-15 | -6.057400e-04 | 0 | 0 | 0 | 2300000 | 0.957172 | 2.86267 | 0.777049 | 0.922857 |
| 23 | 512 | 2 | 0.666667 | 2 | -2.422960e-03 | 8.257280e-16 | -2.422960e-03 | 0 | 0 | 0 | 2400000 | 0.957370 | 2.86267 | 0.778519 | 0.924987 |
| 24 | 512 | 0 | 0.666667 | 0 | -2.982750e-15 | -1.373900e-15 | -2.493660e-18 | 0 | 0 | 0 | 2500000 | 0.957449 | 2.86267 | 0.780096 | 0.926912 |
| 25 | 512 | 5 | 0.666667 | 5 | -1.844300e-01 | -1.692870e-01 | -1.514350e-02 | 0 | 0 | 0 | 2600000 | 0.957587 | 2.86267 | 0.781245 | 0.928673 |
| 26 | 512 | 4 | 0.666667 | 4 | -1.107340e-01 | -1.010420e-01 | -9.691840e-03 | 0 | 0 | 0 | 2700000 | 0.957705 | 2.86267 | 0.782642 | 0.930332 |
| 27 | 512 | 5 | 0.666667 | 5 | -1.170030e-01 | -1.018600e-01 | -1.514350e-02 | 0 | 0 | 0 | 2800000 | 0.957745 | 2.86267 | 0.783507 | 0.931764 |
| 28 | 512 | 4 | 0.666667 | 4 | -7.942230e-02 | -6.973040e-02 | -9.691840e-03 | 0 | 0 | 0 | 2900000 | 0.957887 | 2.86267 | 0.784274 | 0.933213 |
| 29 | 512 | 5 | 0.666667 | 5 | -2.061260e-01 | -1.909830e-01 | -1.514350e-02 | 0 | 0 | 0 | 3000000 | 0.957983 | 2.86267 | 0.785203 | 0.934533 |
| 30 | 512 | 2 | 0.666667 | 2 | -2.775730e-02 | -2.533430e-02 | -2.422960e-03 | 0 | 0 | 0 | 3100000 | 0.957955 | 2.86267 | 0.785910 | 0.935723 |
| 31 | 512 | 5 | 0.666667 | 5 | -3.401790e-01 | -3.250350e-01 | -1.514350e-02 | 0 | 0 | 0 | 3200000 | 0.958009 | 2.86267 | 0.786847 | 0.936868 |
| 32 | 512 | 3 | 0.666667 | 3 | -4.438220e-02 | -3.893050e-02 | -5.451660e-03 | 0 | 0 | 0 | 3300000 | 0.957973 | 2.86267 | 0.787314 | 0.937935 |
| 33 | 512 | 1 | 0.666667 | 1 | -6.057400e-04 | 1.443290e-15 | -6.057400e-04 | 0 | 0 | 0 | 3400000 | 0.958029 | 2.86267 | 0.787905 | 0.938976 |
| 34 | 512 | 0 | 0.666667 | 0 | 9.958400e-16 | -1.221250e-15 | -1.626300e-18 | 0 | 0 | 0 | 3500000 | 0.958059 | 2.86267 | 0.788568 | 0.939989 |
| 35 | 512 | 3 | 0.666667 | 3 | -1.961830e-01 | -1.907310e-01 | -5.451660e-03 | 0 | 0 | 0 | 3600000 | 0.958147 | 2.86267 | 0.789102 | 0.940884 |
| 36 | 512 | 3 | 0.666667 | 3 | -1.011000e-01 | -9.564870e-02 | -5.451660e-03 | 0 | 0 | 0 | 3700000 | 0.958188 | 2.86267 | 0.789876 | 0.941777 |
| 37 | 512 | 1 | 0.666667 | 1 | -6.057400e-04 | 5.606630e-15 | -6.057400e-04 | 0 | 0 | 0 | 3800000 | 0.958289 | 2.86267 | 0.790623 | 0.942630 |
| 38 | 512 | 1 | 0.666667 | 1 | -6.057400e-04 | 2.949030e-15 | -6.057400e-04 | 0 | 0 | 0 | 3900000 | 0.958438 | 2.86267 | 0.791159 | 0.943417 |
| 39 | 512 | 3 | 0.666667 | 3 | -5.451660e-03 | -1.221250e-15 | -5.451660e-03 | 0 | 0 | 0 | 4000000 | 0.958509 | 2.86267 | 0.791675 | 0.944148 |
| 40 | 512 | 3 | 0.666667 | 3 | -5.451660e-03 | 4.135580e-15 | -5.451660e-03 | 0 | 0 | 0 | 4100000 | 0.958471 | 2.86267 | 0.792289 | 0.944843 |
| 41 | 512 | 0 | 0.666667 | 0 | -1.123130e-15 | -3.108620e-15 | -2.493660e-18 | 0 | 0 | 0 | 4200000 | 0.958443 | 2.86267 | 0.792664 | 0.945508 |
| 42 | 512 | 3 | 0.666667 | 3 | -5.451660e-03 | -4.954370e-15 | -5.451660e-03 | 0 | 0 | 0 | 4300000 | 0.958419 | 2.86267 | 0.792928 | 0.946161 |
| 43 | 512 | 5 | 0.666667 | 5 | -3.871480e-02 | -2.357130e-02 | -1.514350e-02 | 0 | 0 | 0 | 4400000 | 0.958501 | 2.86267 | 0.793267 | 0.946816 |
| 44 | 512 | 2 | 0.666667 | 2 | -2.422960e-03 | 1.262880e-15 | -2.422960e-03 | 0 | 0 | 0 | 4500000 | 0.958539 | 2.86267 | 0.793738 | 0.947404 |
| 45 | 512 | 4 | 0.666667 | 4 | -2.267950e-02 | -1.298760e-02 | -9.691840e-03 | 0 | 0 | 0 | 4600000 | 0.958512 | 2.86267 | 0.794078 | 0.947975 |
| 46 | 512 | 0 | 0.666667 | 0 | 5.547750e-15 | 2.331470e-15 | -2.385240e-18 | 0 | 0 | 0 | 4700000 | 0.958516 | 2.86267 | 0.794308 | 0.948526 |
| 47 | 512 | 5 | 0.666667 | 5 | -3.313640e-02 | -1.799290e-02 | -1.514350e-02 | 0 | 0 | 0 | 4800000 | 0.958538 | 2.86267 | 0.794764 | 0.949043 |
| 48 | 512 | 3 | 0.666667 | 3 | -5.451660e-03 | 1.276760e-15 | -5.451660e-03 | 0 | 0 | 0 | 4900000 | 0.958546 | 2.86267 | 0.795101 | 0.949521 |
| 49 | 512 | 5 | 0.666667 | 5 | -3.805840e-02 | -2.291490e-02 | -1.514350e-02 | 0 | 0 | 0 | 5000000 | 0.958578 | 2.86267 | 0.795487 | 0.949949 |
| 50 | 512 | 2 | 0.666667 | 2 | -2.422960e-03 | -2.831070e-15 | -2.422960e-03 | 0 | 0 | 0 | 5100000 | 0.958691 | 2.86267 | 0.795859 | 0.950385 |
| 51 | 512 | 2 | 0.666667 | 2 | -2.422960e-03 | 0.000000e+00 | -2.422960e-03 | 0 | 0 | 0 | 5100000 | 0.958691 | 2.86267 | 0.795859 | 0.950385 |
Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!