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_type=janus:/feasst/plugin/patch/particle/janus.txt add_num_janus_particles=2 \
    group=first,centers first_particle_index=0 centers_site_type=0
Potential Model=HardSphere group=centers
Potential Model=SquareWell VisitModel=VisitModel VisitModelInner=VisitModelInnerPatch group=centers
RefPotential ref=hs Model=HardSphere group=centers
ThermoParams beta={beta}
MayerSampling
TrialTranslate new_only=true ref=hs tunable_param=1 group=first
TrialRotate new_only=true ref=hs tunable_param=40
Let [write]=trials_per_write={trials_per} output_file=patch
Log [write].csv
MoviePatch [write].xyz
CriteriaWriter [write]_b2.txt

# tune trial parameters
Tune
Run num_trials=1e5
Remove name=Tune

# equilibrate
Run num_trials=1e6

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

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

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration add_num_janus_particles=2 centers_site_type=0 cubic_side_length=100 first_particle_index=0 group=first,centers particle_type=janus:/feasst/plugin/patch/particle/janus.txt
Potential Model=HardSphere group=centers
Potential Model=SquareWell VisitModel=VisitModel VisitModelInner=VisitModelInnerPatch group=centers
RefPotential Model=HardSphere group=centers ref=hs
ThermoParams beta=0.1
MayerSampling
TrialTranslate group=first new_only=true ref=hs tunable_param=1
TrialRotate new_only=true ref=hs tunable_param=40
Log output_file=patch.csv trials_per_write=10000.0
MoviePatch output_file=patch.xyz trials_per_write=10000.0
CriteriaWriter output_file=patch_b2.txt trials_per_write=10000.0
Tune
Run num_trials=1e5
# Initializing random number generator with seed: 1749647244
Remove name=Tune
Run num_trials=1e6
MayerSampling
Run num_trials=1e7

 exit: 0

import math
import numpy as np
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'])
assert np.abs(b2_reduced - b2['second_virial_ratio']) < 10*b2['second_virial_ratio_block_stdev']

{'second_virial_ratio': 0.93762, 'second_virial_ratio_block_stdev': 0.000762112}
0.93762

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