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
[1]:
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
[2]:
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!