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

params={"beta": 1./423.2} # 1/K

script="""
MonteCarlo
Configuration cubic_side_length 1e100 particle_type0 ethane_with_ref.fstprt \
    add_particles_of_type0 2 cutoff 1e4 \
    group0 first first_particle_index 0
Potential Model LennardJones
RefPotential Model HardSphere sigma0 4 sigma1 0
ThermoParams beta {beta}
MayerSampling
TrialTranslate new_only true reference_index 0 tunable_param 0.1 group first
TrialRotate new_only true reference_index 0 tunable_param 40
set_variable trials_per 1e5
CriteriaWriter trials_per_write trials_per output_file ethane_b2.txt
Log trials_per_write trials_per output_file ethane.txt
Movie trials_per_write trials_per output_file ethane.xyz

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

# Usage: ./fst < file.txt
FEASST version 0.25.6
MonteCarlo
Configuration add_particles_of_type0 2 cubic_side_length 1e100 cutoff 1e4 first_particle_index 0 group0 first particle_type0 ethane_with_ref.fstprt
Potential Model LennardJones
RefPotential Model HardSphere sigma0 4 sigma1 0
ThermoParams beta 0.0023629489603024575
MayerSampling
TrialTranslate group first new_only true reference_index 0 tunable_param 0.1
TrialRotate new_only true reference_index 0 tunable_param 40
CriteriaWriter output_file ethane_b2.txt trials_per_write 1e5
Log output_file ethane.txt trials_per_write 1e5
Movie output_file ethane.xyz trials_per_write 1e5
Tune
Run num_trials 1e5
# initializing random number generator with seed: 1734453604
Remove name Tune
Run num_trials 1e6
MayerSampling
Run num_trials 1e7

 exit: 0

import math
import numpy as np
with open("ethane_b2.txt") as f:
    firstline = f.readline().rstrip()
    #print(firstline)
    b2=eval(firstline)
    print(b2)
    b2hs = 2./3.*math.pi*3.75**3  # A^3
    b2hs *= 1e-30*1e3*1e3*6.02214076E+23 # L/mol
    print('b22(L/mol)', b2hs*b2['second_virial_ratio'])
    print('b22_block_stdev(L/mol)', b2hs*b2['second_virial_ratio_block_stdev'])
    assert np.abs(-63 - b2hs*b2['second_virial_ratio']) < 20*b2hs*b2['second_virial_ratio_block_stdev']

{'second_virial_ratio': -0.92192, 'second_virial_ratio_block_stdev': 0.00626382}
b22(L/mol) -61.31921067473569
b22_block_stdev(L/mol) 0.4166223731002939

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