Second virial coefficient calculation of TraPPE ethane using Mayer-Sampling

Here, we reproduce the virial coefficient of TraPPE ethane

Table 6 of https://doi.org/10.1021/jp972543+ reports −73 +/- 1 ml/mol at 423.2K

[7]:
import unittest
import feasst as fst

class TestMayerSamplingTraPPEethane(unittest.TestCase):
    def test(self):
        mc = fst.MonteCarlo()
        config = fst.MakeConfiguration(fst.args({"cubic_box_length": str(fst.NEAR_INFINITY)}))
        config.add_particle_type("ethane_with_ref.fstprt")
        config.add_particle_type("ethane_with_ref.fstprt", "2")
        config.add_particle_of_type(0)
        config.add_particle_of_type(1)
        for site in range(config.num_site_types()):
            config.set_model_param("cutoff", site, 1e4)
        mc.add(config)
        mc.add(fst.MakePotential(fst.MakeLennardJones()))

        # reference
        ref = fst.MakePotential(fst.MakeHardSphere())
        ref_params = mc.configuration().model_params().deep_copy()
        for itype in range(mc.configuration().num_site_types()):
            ref_params.set("sigma", itype, 0)
        for itype in [0, 2]:
            for jtype in [0, 2]:
                ref_params.set("sigma", itype, jtype, 4)
        ref.set(ref_params)
        mc.add_to_reference(ref)

        mc.set(fst.MakeThermoParams(fst.args({"beta": str(1./423.2)}))) # 1/K
        mayer = fst.MakeMayerSampling()
        mc.set(mayer)
        mc.add(fst.MakeTrialTranslate(fst.args({"new_only": "true", "reference_index": "0",
            "tunable_param": "0.1", "particle_type": "1"})))
        mc.add(fst.MakeTrialRotate(fst.args({"new_only": "true", "reference_index": "0",
            "tunable_param": "40"})))
        steps_per = "1e5"
        mc.add(fst.MakeLogAndMovie(fst.args({"steps_per": steps_per, "file_name": "ethane"})))
        mc.add(fst.MakeTune(fst.args({"steps_per": steps_per})))
        mc.attempt(int(1e5))
        mc.run(fst.MakeRemoveModify(fst.args({"name": "Tune"})))
        mc.attempt(int(1e6))
        mayer = fst.MakeMayerSampling()
        mc.set(mayer)
        mc.attempt(int(1e7))
        b2hs = 2./3.*fst.PI*mc.configuration().model_params().select("sigma").value(0)**3  # A^3
        b2hs *= 1e-30*1e3*1e3*mc.configuration().physical_constants().avogadro_constant() # L/mol
        print('mayer', mayer.mayer().str())
        print('mayer_ref', mayer.mayer_ref().str())
        print('b22(L/mol)', b2hs*mayer.second_virial_ratio())
        print('b22_block_stdev(L/mol)', b2hs*mayer.second_virial_ratio_block_stdev())
        self.assertAlmostEqual(-63, b2hs*mayer.second_virial_ratio(),
                               delta=3*b2hs*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.

[8]:
%time  # Note: any line starting with % is only to be used with ipynb
unittest.main(argv=[''], verbosity=2, exit=False)
test (__main__.TestMayerSamplingTraPPEethane) ...
CPU times: user 5 µs, sys: 0 ns, total: 5 µs
Wall time: 10 µs
mayer average,stdev,block_stdev,moment0,moment1,moment2,moment3,moment4,
0.28111720000000001,0.95967349237200528,0.023024244809358819,10000000,2811172,10000000,2811172,10000000,
mayer_ref average,stdev,block_stdev,moment0,moment1,moment2,moment3,moment4,
-0.32790070317005948,0.73209108309130144,0.012242208407583844,10000000,-3279007.031700595,6434761.7148547894,-335718717.02519226,89271297422.134458,
b22(L/mol) -57.02278043180012
b22_block_stdev(L/mol) 5.1326695646922404
ok

----------------------------------------------------------------------
Ran 1 test in 25.811s

OK
[8]:
<unittest.main.TestProgram at 0x7f4b2853d940>

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