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
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-104-gc3a7e003a2-dirty-user/newtutorials
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 1e5
Log output_file ethane.txt trials_per 1e5
Movie output_file ethane.xyz trials_per 1e5
Tune
Run num_trials 1e5
# initializing random number generator with seed: 1653572948
RemoveModify name Tune
Run num_trials 1e6
MayerSampling
Run num_trials 1e7

 exit: 0

import math
import unittest
class TestMayerSamplingTraPPEethane(unittest.TestCase):
    def test(self):
        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'])
            self.assertAlmostEqual(-63, b2hs*b2['second_virial_ratio'],
                                   delta=20*b2hs*b2['second_virial_ratio_block_stdev'])
unittest.main(argv=[''], verbosity=2, exit=False)

test (__main__.TestMayerSamplingTraPPEethane) ...

{'second_virial_ratio': -0.908299, 'second_virial_ratio_block_stdev': 0.00480318}
b22(L/mol) -60.41324381361913
b22_block_stdev(L/mol) 0.319471544525205

ok

----------------------------------------------------------------------
Ran 1 test in 0.002s

OK

<unittest.main.TestProgram at 0x7fe7d825e670>

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