Grand canonical ensemble Monte Carlo

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

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.add(fst.MakeConfiguration(fst.args({"cubic_box_length": "8",
                                       "particle_type0": fst.install_dir() + "/forcefield/lj.fstprt"})))
        monte_carlo.add(fst.MakePotential(fst.MakeLennardJones()))
        monte_carlo.add(fst.MakePotential(fst.MakeLongRangeCorrections()))
        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)

%time  # Note: any line starting with % is only to be used with ipynb
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 3 µs, sys: 0 ns, total: 3 µs
Wall time: 7.15 µs

ok

----------------------------------------------------------------------
Ran 1 test in 0.165s

OK

<unittest.main.TestProgram at 0x7f0e6c8bc040>

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