Flat Histogram simulation of a Kern-Frenkel patch

Here, we reproduce the phase behavior of a Kern-Frenkel patchy model using conditions in Figure 2 of https://doi.org/10.1063/1.1569473

import math
import unittest
import feasst as fst

class TestFlatHistogramKFPatch(unittest.TestCase):
    def test(self):
        chi = 0.7 # patch coverage
        mx = 5
        patch_angle = 2*math.asin(math.sqrt(chi/2))*180/math.pi
        print('patch_angle', patch_angle)
        mc = fst.MonteCarlo()
        config = fst.MakeConfiguration(fst.args({"cubic_box_length": "8",
            "particle_type0": fst.install_dir() + "/plugin/patch/forcefield/two_patch_linear.fstprt",
            "patch_angle1": str(patch_angle)}))
        config.add(fst.MakeGroup(fst.args({"site_type0": "0"})))
        mc.add(config)
        mc.add(fst.MakePotential(fst.MakeHardSphere(), fst.args({"group_index": "1"})))
        mc.add(fst.MakePotential(fst.MakeSquareWell(),
            fst.MakeVisitModel(fst.MakeVisitModelInnerPatch()),
            fst.args({"group_index": "1"})))
        mc.set(fst.MakeThermoParams(fst.args({"beta": str(1/0.7), "chemical_potential": "-1.5"})))
        crit = fst.MakeFlatHistogram(
            fst.MakeMacrostateNumParticles(fst.Histogram(fst.args({"width": "1", "max": str(mx), "min": "0"}))),
            fst.MakeTransitionMatrix(fst.args({"min_sweeps": "1000"})))
        mc.set(crit)
        mc.add(fst.MakeTrialTranslate())
        mc.add(fst.MakeTrialRotate())
        mc.add(fst.MakeTrialTransfer(fst.args({"particle_type": "0", "weight": "4"})))
        steps_per = "1e4"
        mc.add(fst.MakeCheckEnergyAndTune(fst.args({"steps_per": steps_per})))
        mc.add(fst.MakeLog(fst.args({"steps_per": steps_per, "file_name": "patch.txt"})))
        mc.add(fst.MakeMoviePatch(fst.args({"steps_per": steps_per, "file_name": "patch.xyz"})))
        mc.add(fst.MakeCriteriaUpdater(fst.args({"steps_per": steps_per})))
        mc.add(fst.MakeCriteriaWriter(fst.args({"steps_per": steps_per, "file_name": "crit.txt"})))
        mc.add(fst.MakeEnergy(fst.args({"steps_per_write": steps_per, "file_name": "en.txt", "multistate": "True"})))
        mc.run_until_complete()

        z_factor = 10
        en = mc.analyze(mc.num_analyzers()-1)
        self.assertAlmostEqual(en.analyze(0).accumulator().average(), 0, delta=fst.NEAR_ZERO)
        self.assertAlmostEqual(en.analyze(1).accumulator().average(), 0, delta=fst.NEAR_ZERO)
        self.assertAlmostEqual(en.analyze(2).accumulator().average(), -0.038758392176564, delta=z_factor*en.analyze(2).accumulator().block_stdev())
        self.assertAlmostEqual(en.analyze(3).accumulator().average(), -0.116517384264731, delta=z_factor*en.analyze(3).accumulator().block_stdev())
        self.assertAlmostEqual(en.analyze(4).accumulator().average(), -0.232665619265520, delta=z_factor*en.analyze(4).accumulator().block_stdev())
        self.assertAlmostEqual(en.analyze(5).accumulator().average(), -0.387804181572135, delta=z_factor*en.analyze(5).accumulator().block_stdev())
        self.assertAlmostEqual(crit.ln_prob().value(0), -15.9976474469475, delta=0.035)
        self.assertAlmostEqual(crit.ln_prob().value(1), -11.9104563420586, delta=0.03)
        self.assertAlmostEqual(crit.ln_prob().value(2), -8.48324267323538, delta=0.025)
        self.assertAlmostEqual(crit.ln_prob().value(3), -5.42988602574393, delta=0.02)
        self.assertAlmostEqual(crit.ln_prob().value(4), -2.64984051640555, delta=0.015)
        self.assertAlmostEqual(crit.ln_prob().value(5), -0.07824246342703, delta=0.01)

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

%time  # Note: any line starting with % is only to be used with ipynb
unittest.main(argv=[''], verbosity=2, exit=False)

test (__main__.TestFlatHistogramKFPatch) ...

CPU times: user 3 µs, sys: 0 ns, total: 3 µs
Wall time: 5.01 µs
patch_angle 72.5423968762779

ok

----------------------------------------------------------------------
Ran 1 test in 37.814s

OK

<unittest.main.TestProgram at 0x7f96540a6fa0>

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