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

params={"beta": 0.1, "trials_per": 1e4} # 1/K

script="""
MonteCarlo
Configuration cubic_side_length 100 particle_type0 /feasst/plugin/patch/particle/janus.fstprt \
    add_particles_of_type0 2 \
    group0 first first_particle_index 0 \
    group1 centers centers_site_type0 0
Potential Model HardSphere group centers
Potential Model SquareWell VisitModel VisitModel VisitModelInner VisitModelInnerPatch group centers
RefPotential Model HardSphere group centers
ThermoParams beta {beta}
MayerSampling
TrialTranslate new_only true reference_index 0 tunable_param 1 group first
TrialRotate new_only true reference_index 0 tunable_param 40
CriteriaWriter trials_per_write {trials_per} file_name patch_b2.txt
Log trials_per_write {trials_per} file_name patch.txt
MoviePatch trials_per_write {trials_per} file_name patch.xyz

# tune trial parameters
Tune
Run num_trials 1e5
RemoveModify name Tune

# equilibrate
Run num_trials 1e6

# production
MayerSampling
Run num_trials 1e7
""".format(**params)

with open('script.txt', 'w') as file: file.write(script)
import subprocess
syscode = subprocess.call("../../../build/bin/fst < script.txt > script.log", shell=True, executable='/bin/bash')
with open('script.log', 'r') as file: print(file.read(), '\n', 'exit:', syscode)

# FEASST version: v0.19.0-107-gda43557a2e-dirty-hwh/newtutorials
Configuration add_particles_of_type0 2 centers_site_type0 0 cubic_side_length 100 first_particle_index 0 group0 first group1 centers particle_type0 /home/hwh/feasst/plugin/patch/particle/janus.fstprt
Potential Model HardSphere group centers
Potential Model SquareWell VisitModel VisitModel VisitModelInner VisitModelInnerPatch group centers
RefPotential Model HardSphere group centers
ThermoParams beta 0.1
MayerSampling
TrialTranslate group first new_only true reference_index 0 tunable_param 1
TrialRotate new_only true reference_index 0 tunable_param 40
CriteriaWriter file_name patch_b2.txt trials_per 10000.0
Log file_name patch.txt trials_per 10000.0
MoviePatch file_name patch.xyz trials_per 10000.0
Tune
Run num_trials 1e5
# initializing random number generator with seed: 1653578111
RemoveModify name Tune
Run num_trials 1e6
MayerSampling
Run num_trials 1e7

 exit: 0

import math
import unittest

class TestMayerSamplingTrimer(unittest.TestCase):
    def test(self):
        b2_reduced = 1-0.5**2*(1.5**3-1)*(math.exp(params['beta'])-1)
        with open("patch_b2.txt") as f:
            firstline = f.readline().rstrip()
            b2=eval(firstline)
            print(b2)
        print(b2['second_virial_ratio'])
        self.assertAlmostEqual(b2_reduced, b2['second_virial_ratio'],
                               delta=10*b2['second_virial_ratio_block_stdev'])
unittest.main(argv=[''], verbosity=2, exit=False)

test (__main__.TestMayerSamplingTrimer) ...

{'second_virial_ratio': 0.937239, 'second_virial_ratio_block_stdev': 0.00101948}
0.937239

ok

----------------------------------------------------------------------
Ran 1 test in 0.001s

OK

<unittest.main.TestProgram at 0x7f9134066310>

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