In [7]:
import inspect
from math import sqrt, pi
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from ipywidgets import interact
from ipywidgets.widgets import FloatSlider, IntSlider
from scipy.integrate import quad
plt.rcParams.update({'font.size': 16})
The particles are here assumed to interact via a Lennard-Jones and a screened Coulomb potential,
$$ \beta w(r) = \frac{\lambda_B z^2}{r} e^{-r/\lambda_D} + 4\beta \varepsilon_{LJ} \left ( \left ( \frac{\sigma}{r}\right )^{12} - \left ( \frac{\sigma}{r}\right )^{6}\right ) $$where $\lambda_B$ and $\lambda_D$ are the Bjerrum and Debye lengths, respectively. Any potential may in principle be given and must return the energy in units of $k_BT$.
potential_rr"$ \mu $" is convenient when writing latex, because in normal strings the backslash acts as an escape character)
In [8]:
# Debye-Huckel
def potential_Debye_Huckel(r, z, D):
r"""$\frac{\lambda_B z^2}{r} e^{-r/\lambda_D}$"""
lB = 7.0 # Bjerrum length, angstroms
return lB * z**2 * np.exp(-r/D) / r
# Lennard-Jones
def potential_Lennard_Jones(r, eps, sigma):
r"""$4\beta \varepsilon_{LJ} \left ( \left ( \frac{\sigma}{r}\right )^{12} - \left ( \frac{\sigma}{r}\right )^{6}\right )$"""
return 4 * eps * ( (sigma/r)**12 - (sigma/r)**6 )
# Total potential
def potential_Combined(r, z, D, eps, sigma):
r"""$\frac{\lambda_B z^2}{r} e^{-r/\lambda_D} + 4\beta \varepsilon_{LJ} \left ( \left ( \frac{\sigma}{r}\right )^{12} - \left ( \frac{\sigma}{r}\right )^{6}\right )$"""
return potential_Debye_Huckel(r, z, D) + potential_Lennard_Jones(r, eps, sigma)
Here we integrate the above pair potential to get the average interaction energy per particle, assuming that the pair correlation function, $g(r)$, can be described by a simple step function, zero when $r<\sigma$, unity otherwise: $$ \hat{a} = -\frac{1}{2} \int_{\sigma}^{\infty} 4\pi w(r) r^2 dr $$
In this Notebook we simply do the integration numerically so that we can use arbitrary pair potentials. From this we calculate the potential, $μ$, versus density, $n$, using,
$$ \beta \, \mu_{gvdW} = \ln \left( \frac{1}{y - y_0} \right) + \frac{y}{y-y_0} - 2 \frac{\hat{a}}{y} $$From this we calculate the pressure, $p$, versus density, $n$, using,
where $y=1/n$ and $y_0=\pi \sigma^2/2$ is the particle area.
For reference we'll also plot the chemical potential for an ideal system (van 't Hoff), $\beta \mu_{ideal}= \ln (1/y)$, where $\beta = 1/k_BT$.
In [20]:
def ahat(potential, **parameters):
sigma = parameters['sigma']
# extract the relevant parameters for the potential
parameters = {k:v for k,v in parameters.items() if k in inspect.signature(potential).parameters}
def integrand(r):
return potential(r, **parameters) * r**2
integral, error = quad(integrand, sigma, np.infty, limit=50)
return -2 * pi * integral
def ahatexact(z, D, eps, sigma):
return -2 * pi * (-8/9 * eps * sigma**3 + 7 * np.exp(-sigma/D) * z**2 * (D + sigma))
def mu_ideal(n):
return np.log(n)
def mu_gvdw_backup(n, z, D, eps, sigma, potential=potential_Combined):
y0 = pi*sigma**2 / 2
y = 1 / n
a = ahat(potential, z=z, D=D, eps=eps, sigma=sigma)
return np.log(n) + np.log(y / (y-y0)) + y / (y-y0) - 2 * a/y
def mu_gvdw(n, z, D, eps, sigma, potential=potential_Combined):
y0 = pi*sigma**2 / 2 # excluded volume
y = 1 / n
a = ahat(potential, z=z, D=D, eps=eps, sigma=sigma)
return - np.log(y-y0) - y0 / (y - y0) - 2 * a / y
In [21]:
n = np.linspace(1e-4, 1e-3, 100)
_potentials = {fname[10:]: func for fname, func in globals().items() if fname.startswith("potential_")}
sliders = {
"eps": FloatSlider(min=0, max=10, step=0.1, value=1, description=r'LJ , $\varepsilon_{LJ}$ ($\beta$)'),
"sigma": FloatSlider(min=0, max=10, step=0.1, value=4, description=r'Size, $\sigma_{LJ}$ (Å)'),
"z": IntSlider(min=0, max=3, step=1, value=1, description=r'Charge, $z$ ($e$)'),
"Cs": FloatSlider(min=1e-3, max=1, step=0.1, value=0.3, description="Salt, $c_s$ (M)"),
"potential": _potentials
}
@interact(**sliders)
def plot_EOS( eps=1.0, sigma=4.0, z=1.0, Cs=0.3, potential=potential_Combined):
D = 3.04/sqrt(Cs)
plt.figure(figsize=(10, 10/1.618), )
plt.plot(n, 10 + mu_ideal(n), 'k--', label='ideal', lw=2)
plt.plot(n, 10 + mu_gvdw(n, z, D, eps, sigma, potential=potential), 'r-', label=potential.__doc__ or potential.__name__, lw=2)
plt.title('Equation of State')
plt.xlabel(r'Number density, $n$')
plt.ylabel(r'Potential, $\beta \mu$')
plt.ylim([-2,10])
plt.legend(loc=4, frameon=False)
plt.show()
In [10]:
with open("datafile.csv", "wt") as stream:
stream.write("""\
length potential proteins density
1.414000000000000000e+03 1.429585460731450097e-01 2.000000000000000000e+01 1.000302091231552026e+03
9.990000000000000000e+02 2.990882091428900269e-01 2.000000000000000000e+01 2.004006008010011783e+03
8.160000000000000000e+02 4.684432472751309806e-01 2.000000000000000000e+01 3.003652441368703876e+03
6.320000000000000000e+02 8.629727385734929923e-01 2.000000000000000000e+01 5.007210382951449901e+03
5.340000000000000000e+02 1.353602607621670062e+00 2.000000000000000000e+01 7.013704779138435697e+03
4.710000000000000000e+02 1.970895704549270100e+00 2.000000000000000000e+01 9.015466031977857710e+03
4.260000000000000000e+02 2.788653065634310035e+00 2.000000000000000000e+01 1.102074103462716812e+04
3.920000000000000000e+02 3.842403663548089821e+00 2.000000000000000000e+01 1.301541024573094546e+04""")
In [11]:
df = pd.read_csv("datafile.csv", delimiter="(?:\s+|,)", engine="python")
In [12]:
def plot_EOS( eps=1.0, sigma=4.0, z=0.0, Cs=0.3, potential=potential_Combined):
D = 3.04/sqrt(Cs)
# plt.title(potential.__doc__)
plt.plot(n, 10+mu_ideal(n), 'k--', label='ideal', lw=2)
plt.plot(n, 10+mu_gvdw(n, z, D, eps, sigma, potential=potential), 'r-', label=potential.__doc__ or potential.__name__, lw=2)
plt.plot(df.density/10**8, df.potential + mu_ideal((df.density/10**8)) + 10, label="data")
plt.title('Equation of State')
plt.xlabel(r'Number density, $n$')
plt.ylabel(r'Potential, $\beta \mu$')
plt.legend(loc=0, frameon=False)
plt.show()
data_density = df.density/10**8
n = np.linspace(data_density[0], list(data_density)[-1], 100)
_potentials = {fname: func for fname, func in globals().items() if fname.startswith("potential_")}
i = interact(plot_EOS,
eps=(0.0, 10.0, 0.1),
sigma=(0, 10, 0.1),
z=(0.0, 3, 1.0),
Cs=(1e-3, 1.0, 0.1),
potential = _potentials )
i.widget.children[0].description=r'$\beta\varepsilon_{LJ}$'
i.widget.children[1].description=r'$\sigma_{LJ}$'
i.widget.children[2].description=r'$z$'
i.widget.children[3].description=r'$c_s$ (M)'