Grand canonical ensemble Monte CarloΒΆ

In this example, a short grand canonical Monte Carlo simulation of Lennard Jones particles is conducted.

[1]:
import unittest
import feasst as fst

class TestMonteCarlo2LJGCMC(unittest.TestCase):
    """Test a grand canonical ensemble Lennard Jones Monte Carlo simulation"""
    def test(self):
        """Compute the average number of particles and assert that it is greater than 0"""
        monte_carlo = fst.MonteCarlo()
        monte_carlo.set(fst.lennard_jones())
        monte_carlo.set(fst.MakeThermoParams(fst.args({"beta": str(1./1.5), "chemical_potential": "-8.352321"})))
        monte_carlo.set(fst.MakeMetropolis())
        monte_carlo.add(fst.MakeTrialTranslate(fst.args({"weight": "1.", "tunable_param": "2."})))
        monte_carlo.add(fst.MakeTrialTransfer(fst.args({"weight": "1.", "particle_type": "0"})))
        steps_per = int(1e3)
        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)})))

        # Add an Analyze which computes the average number of particles.
        # Just before adding, store the number of existing Analyzers in order to remember the
        # index of the newly added Analyze.
        analyze_index = monte_carlo.num_analyzers()
        monte_carlo.add(fst.MakeNumParticles(fst.args(
            {"steps_per_write": str(steps_per), "file_name": "gcmc_num_particles.txt"})))

        # peform a short simulation
        monte_carlo.attempt(int(1e5))

        # assert that particles were added during the simulation
        self.assertTrue(monte_carlo.analyze(analyze_index).accumulator().average() > 0)

[2]:
%%time
unittest.main(argv=[''], verbosity=2, exit=False)
test (__main__.TestMonteCarlo2LJGCMC)
Compute the average number of particles and assert that it is greater than 0 ...
CPU times: user 188 ms, sys: 3.7 ms, total: 192 ms
Wall time: 190 ms
ok

----------------------------------------------------------------------
Ran 1 test in 0.187s

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

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