Second virial coefficient calculation of a Trimer using Mayer-Sampling

Here, we reproduce the virial coefficient of a trimer at the Boyle temperature as reported Table 1 of the following manuscript for \(L=0.4\sigma\)

https://dx.doi.org/10.1063/1.4918557

params={"beta": 1./0.815, "rwca": 2.**(1./6.)}

script="""
MonteCarlo
Configuration cubic_side_length=1e100 particle_type=trimer:/feasst/particle/trimer_0.4L.txt \
    add_num_trimer_particles=2 cutoff0_1={rwca} cutoff1_1={rwca} \
    group=first first_particle_index=0
Potential Model=LennardJonesForceShift
RefPotential ref=hs Model=HardSphere sigma0=1 sigma1=0
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=1e4 output_file=trimer
Log [write].txt
Movie [write].xyz
CriteriaWriter [write]_b2.csv
Run num_trials 1e6
""".format(**params)

with open('script1.txt', 'w') as file: file.write(script)
import subprocess
syscode = subprocess.call("../../../build/bin/fst < script1.txt > script1.log", shell=True, executable='/bin/bash')
with open('script1.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_trimer_particles=2 cubic_side_length=1e100 cutoff0_1=1.122462048309373 cutoff1_1=1.122462048309373 first_particle_index=0 group=first particle_type=trimer:/feasst/particle/trimer_0.4L.txt
Potential Model=LennardJonesForceShift
RefPotential Model=HardSphere ref=hs sigma0=1 sigma1=0
ThermoParams beta=1.2269938650306749
MayerSampling
TrialTranslate group=first new_only=true ref=hs tunable_param=1
TrialRotate new_only=true ref=hs tunable_param=40
Log output_file=trimer.txt trials_per_write=1e4
Movie output_file=trimer.xyz trials_per_write=1e4
CriteriaWriter output_file=trimer_b2.csv trials_per_write=1e4
Run num_trials=1e6
# Initializing random number generator with seed: 1749647129

 exit: 0

import math
import numpy as np
with open("trimer_b2.csv") as f:
    firstline = f.readline().rstrip()
    #print(firstline)
    b2=eval(firstline)
    #b2={"second_virial_ratio": -0.0166853, "second_virial_ratio_block_stdev": 0.0159632}
    print(b2)
# print(mayer.mayer().str())
    assert np.abs(b2['second_virial_ratio']) < 4*b2['second_virial_ratio_block_stdev']

{'second_virial_ratio': 0.0129171, 'second_virial_ratio_block_stdev': 0.0152898}

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