Block Averaging Example

In this example, the block averaging approach to computing error bars is considered. First, please read the original method: https://doi.org/10.1063/1.457480. Then, also consider the discussion on block averages in Appendix D of Frenkel and Smit. Then, we’ll reproduce Frenkel and Smith’s case study 4 (page 98-100 and Figure 4.4).

import sys
import unittest
import feasst as fst

class TestMonteCarlo4Block(unittest.TestCase):
    def test_srsw(self, num_particles=108, density=0.8442, temperature=1.5, steps_per=1e5):
        monte_carlo = fst.MonteCarlo()
        # monte_carlo.set(fst.MakeRandomMT19937(fst.args({"seed": "1234"})))
        monte_carlo.add(fst.MakeConfiguration(fst.args({
            "cubic_box_length": str((num_particles/density)**(1./3.)),
            "particle_type0": fst.install_dir() + "/forcefield/lj.fstprt",
            "cutoff0": "2.5",
        })))
        monte_carlo.add(fst.MakePotential(fst.MakeLennardJones()))
        monte_carlo.add(fst.MakePotential(fst.MakeLongRangeCorrections()))
        monte_carlo.set(fst.MakeThermoParams(fst.args({"beta": str(1./temperature),
                                                       "chemical_potential": "1."})))
        monte_carlo.set(fst.MakeMetropolis())
        monte_carlo.add(fst.MakeTrialTranslate(fst.args({"weight": "1.", "tunable_param": "2."})))
        monte_carlo.add(fst.MakeTrialAdd(fst.args({"particle_type": "0"})))
        monte_carlo.run(fst.MakeRun(fst.args({"until_num_particles": str(num_particles)})))
        monte_carlo.run(fst.MakeRemoveTrial(fst.args({"name": "TrialAdd"})))
        monte_carlo.add(fst.MakeLogAndMovie(fst.args({
            "steps_per" : str(steps_per),
            "file_name": "movie",
            "clear_file": "true"})))
        monte_carlo.add(fst.MakeCheckEnergyAndTune(fst.args({"steps_per" : str(steps_per)})))

        # equilibrate
        monte_carlo.attempt(int(1e6))

        # compute average using this script
        energy = fst.MakeAccumulator()

        # production
        for _ in range(int(1e8)):
            monte_carlo.attempt()
            energy.accumulate(monte_carlo.criteria().current_energy())

        print('average energy', energy.average())
        for op in range(energy.max_block_operations()):
            print(op, energy.block_stdev(op), energy.block_std_of_std(op))

        #self.assertAlmostEqual(energy.block_stdev(13)/num_particles, 0.0012, delta=0.0004)
        #self.assertAlmostEqual(energy.average()/num_particles, -4.4190, delta=3*energy.block_stdev(13)/num_particles)

If the test passes, the energy is within the tolerance of the SRSW value and the two ensemble average methods agreed.

%time  # Note: any line starting with % is only to be used with ipynb
unittest.main(argv=[''], verbosity=2, exit=False)

test_srsw (__main__.TestMonteCarlo4Block) ...

CPU times: user 3 µs, sys: 0 ns, total: 3 µs
Wall time: 7.15 µs
average energy -576.4376128069971
0 0.07510636040234885 0.04322275467550779
1 0.07359011506962478 0.04749864777552218
2 0.07344944697207455 0.0589605210359109
3 0.07239626945616733 0.06608570396503526
4 0.05376806537982441 0.0
5 0.0 0.0

ok

----------------------------------------------------------------------
Ran 1 test in 1757.834s

OK

<unittest.main.TestProgram at 0x7f1b641185e0>

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