Second virial coefficient calculation of a Kern-Frenkel patch using Mayer-Sampling

Here, we reproduce the virial coefficient of a Kern-Frenkel patchy particle and compare to Eq. 8 of https://doi.org/10.1063/1.1569473

import math
import unittest
import feasst as fst

class TestMayerSamplingTrimer(unittest.TestCase):
    def test(self):
        mc = fst.MonteCarlo()
        config = fst.MakeConfiguration(fst.args({"cubic_box_length": str(fst.NEAR_INFINITY)}))
        config.add_particle_type(fst.install_dir() + "/plugin/patch/forcefield/janus.fstprt")
        config.add_particle_type(fst.install_dir() + "/plugin/patch/forcefield/janus.fstprt", "2")
        config.add(fst.MakeGroup(fst.args({"site_type0": "0", "site_type1": "2"})))
        config.add_particle_of_type(0)
        config.add_particle_of_type(1)
        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.add_to_reference(fst.MakePotential(fst.MakeHardSphere(), fst.args({"group_index": "1"})))
        beta = 0.1
        mc.set(fst.MakeThermoParams(fst.args({"beta": str(beta)})))
        mayer = fst.MakeMayerSampling()
        mc.set(mayer)
        mc.add(fst.MakeTrialTranslate(fst.args({"new_only": "true", "reference_index": "0",
            "tunable_param": "1", "particle_type": "1"})))
        mc.add(fst.MakeTrialRotate(fst.args({"new_only": "true", "reference_index": "0",
            "tunable_param": "40"})))
        steps_per = "1e4"
        mc.add(fst.MakeLogAndMovie(fst.args({"steps_per": steps_per, "file_name": "patch"})))
        mc.attempt(int(1e6))
        b2_reduced = 1-0.5**2*(1.5**3-1)*(math.exp(beta)-1)
        print(mayer.second_virial_ratio())
        self.assertAlmostEqual(b2_reduced, mayer.second_virial_ratio(),
                               delta=10*mayer.second_virial_ratio_block_stdev())

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__.TestMayerSamplingTrimer) ...

CPU times: user 3 µs, sys: 0 ns, total: 3 µs
Wall time: 6.68 µs
0.93850099356307

ok

----------------------------------------------------------------------
Ran 1 test in 2.336s

OK

<unittest.main.TestProgram at 0x7f16052668e0>

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