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_type=ethane:ethane_with_ref.txt \
    add_num_ethane_particles=2 cutoff=1e4 \
    group=first first_particle_index=0
Potential Model=LennardJones
RefPotential ref=hs Model=HardSphere sigmaC1=4 sigmaC2=0
ThermoParams beta={beta}
MayerSampling
TrialTranslate new_only=true ref=hs tunable_param=0.1 group=first
TrialRotate new_only=true ref=hs tunable_param=40
Let [write]=trials_per_write=1e5 output_file=ethane
Log [write].csv
Movie [write].xyz
CriteriaWriter [write]_b2.csv

# 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: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration add_num_ethane_particles=2 cubic_side_length=1e100 cutoff=1e4 first_particle_index=0 group=first particle_type=ethane:ethane_with_ref.txt
Potential Model=LennardJones
RefPotential Model=HardSphere ref=hs sigmaC1=4 sigmaC2=0
ThermoParams beta=0.0023629489603024575
MayerSampling
TrialTranslate group=first new_only=true ref=hs tunable_param=0.1
TrialRotate new_only=true ref=hs tunable_param=40
Log output_file=ethane.csv trials_per_write=1e5
Movie output_file=ethane.xyz trials_per_write=1e5
CriteriaWriter output_file=ethane_b2.csv trials_per_write=1e5
Tune
Run num_trials=1e5
# Initializing random number generator with seed: 1749647178
Remove name=Tune
Run num_trials=1e6
MayerSampling
Run num_trials=1e7

 exit: 0

import math
import numpy as np
with open("ethane_b2.csv") 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.910745, 'second_virial_ratio_block_stdev': 0.00584095}
b22(L/mol) -60.575933406328275
b22_block_stdev(L/mol) 0.3884962291636991

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