Canonical ensemble Monte CarloΒΆ

In this example, a low density simulation is conducted with a constant number of particles. The average energy is compared to published results from the NIST SRSW. https://www.nist.gov/programs-projects/nist-standard-reference-simulation-website

The log file is designed to monitor simulations but may not be the best route to computing properties. For example, the energy is output in the log file, but this tutorial shows two alternate ways to compute average energy. First, an Analyze, Energy, to compute average energies. Second, the Monte Carlo simulation attempts were performed one step at a time to allow for accumulation of the ensemble average directly in the Python script.

[1]:
import sys
import unittest
import feasst as fst

class TestMonteCarlo1LJNVT(unittest.TestCase):
    """Test a canonical ensemble Lennard Jones Monte Carlo simulation"""
    def test_srsw(self, num_particles=500, density=0.001, steps_per=1e5):
        """Compare with the reported average energy from the NIST SRSW.
        https://mmlapps.nist.gov/srs/LJ_PURE/mc.htm
        https://www.nist.gov/programs-projects/nist-standard-reference-simulation-website

        num_particles -- number of LJ particles
        density -- number density
        steps_per -- steps between each Anaylze/Modify
        """

        monte_carlo = fst.MonteCarlo()
        # monte_carlo.set(fst.MakeRandomMT19937(fst.args({"seed": "1234"})))
        monte_carlo.set(fst.lennard_jones(fst.args({
            "cubic_box_length": str((num_particles/density)**(1./3.)),
            "dual_cut": "1"})))
        monte_carlo.set(fst.MakeThermoParams(fst.args({"beta": str(1./0.9), "chemical_potential": "1."})))
        monte_carlo.set(fst.MakeMetropolis())
        monte_carlo.add(fst.MakeTrialTranslate(fst.args({
            "weight": "1.",
            "tunable_param": "2.",
            "num_steps": "2",       # attempt a number of configurational bias steps
            "reference_index": "0", # using an (optimized) reference potential
        })))
        fst.SeekNumParticles(num_particles).with_trial_add().run(monte_carlo)

        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(1e7))

        # compute average energy using a stepper/analysis and output into file
        energy = fst.MakeEnergy(fst.args({
            "steps_per_update": "1",
            "steps_per_write": str(steps_per),
            "file_name": "lj_nvt_srsw_energy.txt",
        }))
        monte_carlo.add(energy)

        # compute average using this script
        energy_alt = fst.Accumulator()

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

        # test that the two methods to compute average energy are the same
        self.assertAlmostEqual(energy.energy().average(), energy_alt.average(), delta=1e-6)

        # test the average against the NIST SRSW
        stdev = (energy.energy().block_stdev()**2 + (1.89E-05)**2)**(1./2.)
        # print(energy.energy().average(), energy_alt.average())
        self.assertAlmostEqual(-9.9165E-03*num_particles, energy.energy().average(),
                               delta=2.576*stdev)

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

[2]:
%%time
unittest.main(argv=[''], verbosity=2, exit=False)
test_srsw (__main__.TestMonteCarlo1LJNVT)
Compare with the reported average energy from the NIST SRSW. ...
CPU times: user 21min 15s, sys: 976 ms, total: 21min 16s
Wall time: 21min 15s
ok

----------------------------------------------------------------------
Ran 1 test in 1275.158s

OK
[2]:
<unittest.main.TestProgram at 0x7f4d805104a8>

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