# 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())
"weight": "1.",
"tunable_param": "2.",
"num_steps": "2",       # attempt a number of configurational bias steps
"reference_index": "0", # using an (optimized) reference potential
})))

"steps_per" : str(steps_per),
"file_name": "movie",
"clear_file": "true"})))

# 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",
}))

# 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>