<jupyter_start><jupyter_text>PMF and virial coefficient calculation using VirtualTranslateLet molecule A (index 0) be centered in the the middle of a spherical simulation cell and allowed to _rotate only_.Place another molecule, B (index 1), in on the positive $z$-axis and let it translate (`dir=[0,0,1]`) and rotate in such a way that itcannot cross to the negative $z$-axis. Add the following to the analysis:~~~ yml- virtualtranslate: {file: vt.dat, dL: 0.1, molecule: B, nstep: 10, dir: [0,0,1]}- reactioncoordinate: {file: R.dat, nstep: 10, type: molecule, property: com_z, index: 1}~~~<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import numpy as np
from scipy.stats import binned_statistic

def getForces(xfile, yfile, bins):
    ''' get average force as a function of z-position of molecule B '''
    R = np.loadtxt(xfile, usecols=[1])
    steps, dL, du, force = np.loadtxt(yfile, unpack=True, skiprows=1)
    means, edges, bins = binned_statistic(R, np.exp(-du), 'mean', bins)
    return edges[1:], (np.log(means) / dL[0])
    #return (edges[:-1] + edges[1:]) / 2, np.log(means) / dL[0]<jupyter_output><empty_output><jupyter_text>Run Faunus from bash<jupyter_code>%%bash
if [[ -z "${FAUNUS_EXECUTABLE}" ]]; then
  yason.py virial.yml | faunus --nobar
else
  echo "Seems we're running CTest - use Faunus target from CMake"
  "${YASON_EXECUTABLE}" virial.yml | "${FAUNUS_EXECUTABLE}" --nobar
fi
# calculate mean force as a function of 
R, force = getForces('R.dat', 'vt.dat', 150)
def lj(r, sigma=5, eps=0.40090549):
    return 4*eps*( (sigma/r)**12 - (sigma/r)**6  )

# LJ contribution to B2 from zero to contact
R0 = np.linspace(1e-6, R.min(), 500)
B2hs = -2 * np.pi * np.trapz( np.expm1( -lj(R0) )*R0**2, R0 )

# integrate simulated force to get PMF and normalize
pmf = -np.cumsum(force)*np.diff(R)[:-2].mean()
pmf = pmf - pmf[-1]

# integrate to get second virial
print("B2_hs = ", B2hs)
print("B2_lj = ", B2hs -2 * np.pi * np.trapz( np.expm1(-lj(R))*R**2, R )) 
print("B2_vt = ", B2hs -2 * np.pi * np.trapz( np.expm1(-pmf)*R**2, R ))

plt.plot(R, np.expm1(-pmf)*R**2, label='from force calc')
plt.plot(R, np.expm1(-lj(R))*R**2, label='LJ')
plt.legend(frameon=False, loc=0)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16511318242343542, 'fraction_numerical': 0.032401242787394584, 'mean_word_length': 4.551724137931035, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '80396', 'license': 'MIT', 'n_tokens_mistral': 932, 'n_tokens_neox': 880, 'n_words': 254, 'path': 'examples/virial/virial.ipynb', 'repo_name': 'rc83/faunus', 'size': '2253'}