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_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("../../../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.25.2-4-g362d22fce9-user/qtask
# Usage: ./fst < file.txt
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/user/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: 1726668049
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.test_serial_5max)
Compare the free energies and potential energies with the NIST SRSW ... ok
----------------------------------------------------------------------
Ran 1 test in 0.004s
OK
delta_ln_prob delta_ln_prob_stdev
0 NaN 0.000000
1 4.67368 0.001822
2 3.98603 0.001865
3 3.58937 0.002211
4 3.31294 0.001840
5 3.10077 0.002147
energy average block_stdev
0 5.386599e-16 2.057114e-16
1 -6.057400e-04 7.702397e-10
2 -3.071674e-02 1.870853e-04
3 -8.982909e-02 4.264981e-04
4 -1.779654e-01 9.794333e-04
5 -2.966230e-01 2.148960e-03
[3]:
<unittest.main.TestProgram at 0x7fe2515f5f70>
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 | 512 | 5 | 0.666667 | 5 | -1.04842 | -1.03327 | -0.0151435 | 100000 | 0.921757 | 2.79198 | 0.046112 | 0.0461438 |
| 1 | 512 | 1 | 0.666667 | 1 | -0.00060574 | 4.45997e-15 | -0.00060574 | 200000 | 0.950677 | 2.79198 | 0.249782 | 0.250912 |
| 2 | 512 | 0 | 0.666667 | 0 | -7.30221e-15 | 1.27676e-15 | -1.6263e-18 | 300000 | 0.954167 | 2.79198 | 0.434975 | 0.46054 |
| 3 | 512 | 4 | 0.666667 | 4 | -0.198234 | -0.188542 | -0.00969184 | 400000 | 0.955229 | 2.79198 | 0.528943 | 0.577013 |
| 4 | 512 | 3 | 0.666667 | 3 | -0.782811 | -0.77736 | -0.00545166 | 500000 | 0.955989 | 2.79198 | 0.585786 | 0.650345 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 100 | 512 | 0 | 0.666667 | 0 | 2.51286e-15 | -1.76421e-15 | -1.51788e-18 | 4800000 | 0.958491 | 2.79046 | 0.791414 | 0.942253 |
| 101 | 512 | 0 | 0.666667 | 0 | -1.79869e-16 | -2.26555e-15 | -1.6263e-18 | 4900000 | 0.958516 | 2.79046 | 0.791782 | 0.942893 |
| 102 | 512 | 5 | 0.666667 | 5 | -0.382384 | -0.367241 | -0.0151435 | 5000000 | 0.958548 | 2.79046 | 0.792158 | 0.9435 |
| 103 | 512 | 1 | 0.666667 | 1 | -0.00060574 | -2.27509e-15 | -0.00060574 | 5100000 | 0.958609 | 2.79046 | 0.792561 | 0.944102 |
| 104 | 512 | 1 | 0.666667 | 1 | -0.00060574 | 0 | -0.00060574 | 5100000 | 0.958609 | 2.79046 | 0.792561 | 0.944102 |
105 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!