Analysis of a hard sphere simulation

# compare to EOS in SRSW: https://www.nist.gov/mml/csd/chemical-informatics-research-group/hard-sphere-thermodynamic-and-transport-properties
import numpy as np
import pandas as pd
import feasst as fst
import pyfeasst
import subprocess

proc = subprocess.Popen(['python', 'tutorial_10_hs.py', '--num_procs', '4', '--max_particles', '100', '--min_sweeps', '20'], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
proc.wait()
df = pd.read_csv('ln_prob.txt', header=None)
ln_prob = fst.LnProbability(fst.DoubleVector(df[0]))
gce = fst.GrandCanonicalEnsemble(
    fst.Histogram(fst.args({"width": "1", "max": str(len(df)-1)})),
    ln_prob,
    -2.352321) # original conjugate, beta_mu = lnz

volume = 8**3
df = pd.read_csv(fst.install_dir() + '/plugin/flat_histogram/test/data/stat_hs.csv')
df = df[:6] # truncate high density because 100 particles in this example isn't enough. Instead, try 512 or more.

def objective_fn(target_density, dbetamu):
    gce.reweight(dbetamu)
    return (target_density - gce.average_macrostate()/volume)**2

from scipy.optimize import minimize
pressure=list()
for target_density in df['dens']:
    res = minimize(lambda beta_mu_rw: objective_fn(target_density, beta_mu_rw[0]), 1., tol=1e-8)
    gce.reweight(res.x[0])
    pressure.append(gce.betaPV()/volume)
df['P_FST'] = pressure
df.to_csv('tutorial_10_hs.csv')
assert(abs(df['P_MC'][4] - df['P_FST'][4]) < 1e-3)
print(df[['dens', 'P_MC', 'P_FST', '+/-']])

# plt.plot(df['dens'], df['P_FST'], label='fst')
# plt.plot(df['dens'], df['P_MC'], linestyle='dashed', label='srsw')
# plt.xlabel('density')
# plt.ylabel('pressure')
# plt.legend()

    dens      P_MC     P_FST       +/-
0  0.025  0.026352  0.026376  0.000000
1  0.050  0.055593  0.055662  0.000000
2  0.075  0.088023  0.088048  0.000002
3  0.100  0.123969  0.123898  0.000001
4  0.125  0.163795  0.163772  0.000002
5  0.150  0.207877  0.207817  0.000003

# Use chemical potential from Carnahan-Starling to compare expected average density
# http://www.sklogwiki.org/SklogWiki/index.php/Carnahan-Starling_equation_of_state
import math
import subprocess
import pandas as pd
import feasst as fst
rho = 0.1
cubic_box_length=8
eta = fst.PI/6*rho
betamu_ex = (8*eta-9*eta**2+3*eta**3)/(1-eta)**3
betamu = betamu_ex + math.log(rho)
print('betamu', betamu)
proc = subprocess.Popen(['python', 'tutorial_10_hs.py', '--num_procs', '4', '--max_particles', '115',
                         '--cubic_box_length', str(cubic_box_length), '--min_sweeps', '20',
                         '--mu', str(betamu)], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
proc.wait()
df = pd.read_csv('ln_prob.txt', header=None)
ln_prob = fst.LnProbability(fst.DoubleVector(df[0]))
gce = fst.GrandCanonicalEnsemble(
    fst.Histogram(fst.args({"width": "1", "max": str(len(df)-1)})),
    ln_prob,
    betamu)
print('target density', rho)
print('density', gce.average_macrostate()/cubic_box_length**3)
assert(abs(rho - gce.average_macrostate()/cubic_box_length**3) < 1e-4)

betamu -1.838854233754734
target density 0.1
density 0.09997942215978713

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