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 cycles 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 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)
# 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 complete
# 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!