Grand canonical ensemble transition-matrix Monte Carlo with RPM

This example is similar to the SPC/E tutorial, except this time we simulate the RPM model.

import unittest
import feasst as fst

class TestFlatHistogramRPM(unittest.TestCase):
    """Test flat histogram grand canonical ensemble Monte Carlo simulations"""
    def test_serial_5max(self):
        """Compare the free energies and potential energies with the previously
        published values: https://doi.org/10.1063/1.5123683
        """
        monte_carlo = fst.MonteCarlo()
        monte_carlo.set(fst.rpm(fst.args({
            "cubic_box_length": "12",
            "cutoff": "4.891304347826090",
            "alphaL": "6.87098396396261"})))
        monte_carlo.add_to_reference(fst.MakePotential(fst.MakeDontVisitModel()))
        temperature = 0.047899460618081;
        beta_mu = -13.94;

        monte_carlo.set(fst.MakeThermoParams(fst.args({"beta": str(1/temperature),
                      "chemical_potential0": str(beta_mu*temperature),
                      "chemical_potential1": str(beta_mu*temperature)})))
        criteria = fst.MakeFlatHistogram(
            fst.MakeMacrostateNumParticles(
                fst.Histogram(fst.args({"width": "1", "max": "2", "min": "0"})),
                fst.args({"particle_type": "0"})),
            fst.MakeTransitionMatrix(fst.args({"min_sweeps": "100"})))
        monte_carlo.set(criteria)
        monte_carlo.add(fst.MakeTrialTranslate(fst.args({"weight": "0.25", "tunable_param": "0.1"})))
        monte_carlo.add(fst.MakeTrialTransferMultiple(fst.args({
            "weight": "4",
            "particle_type0": "0",
            "particle_type1": "1",
            "reference_index": "0"})))
        steps_per = str(int(1e5))
        monte_carlo.add(fst.MakeCriteriaUpdater(fst.args({"steps_per": steps_per})))
        monte_carlo.add(fst.MakeCriteriaWriter(fst.args({"steps_per": steps_per, "file_name": "crit.txt"})))
        monte_carlo.add(fst.MakeLogAndMovie(fst.args({"steps_per": steps_per, "file_name": "rpm"})))
        monte_carlo.add(fst.MakeCheckEnergyAndTune(fst.args({"steps_per": steps_per})))
        energy = fst.MakeEnergy(fst.args({"file_name": "rpm_fh_energy",
                                          "steps_per_update": "1",
                                          "steps_per_write": steps_per,
                                          "multistate": "true"}))
        monte_carlo.add(energy)
        monte_carlo.run_until_complete()

        lnpi_previous = [
            [-1.2994315780357, 0.05],
            [-1.08646312498868, 0.05],
            [-0.941850889679828, 0.05]]
        energy_previous = [
            [0, 1e-14],
            [-0.939408, 0.02],
            [-2.02625, 0.04]]
        for macro in range(criteria.num_states()):
            self.assertAlmostEqual(
                lnpi_previous[macro][0],
                criteria.bias().ln_prob().value(macro),
                delta=5.*lnpi_previous[macro][1]
            )
            energy_analyzer = monte_carlo.analyze(monte_carlo.num_analyzers() - 1)
            energy_accumulator = energy_analyzer.analyze(macro).accumulator()
            stdev = (energy_previous[macro][1]**2 + energy_accumulator.block_stdev()**2)**(1./2.)
            self.assertAlmostEqual(
                energy_previous[macro][0],
                energy_accumulator.average(),
                #criteria.bias().ln_prob().value(macro),
                delta=10*stdev
            )

%time
unittest.main(argv=[''], verbosity=2, exit=False)

test_serial_5max (__main__.TestFlatHistogramRPM)
Compare the free energies and potential energies with the previously ...

CPU times: user 3min 30s, sys: 253 ms, total: 3min 30s
Wall time: 3min 30s

ok

----------------------------------------------------------------------
Ran 1 test in 210.298s

OK

<unittest.main.TestProgram at 0x7fda8c209cf8>

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