<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Equation of State using Generalized van der Waals Theory
This notebook illustrates the generalized van der Waals theory (gvdW) for the equation of state for interacting particles. Based on the lecture notes, Properties of Molecular Fluids in Equilibrium by Sture Nordholm.
Step4:  Pair potentials
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$.
Define your own potentials below
The name should start with potential_
The first parameter should be r
The docstring can be added for a nice displayname of the function (The raw python string like r"$ \mu $" is convenient when writing latex, because in normal strings the backslash acts as an escape character)
Step5:  Interaction parameter
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
Step6:  Interactive EOS plot
Step8:  With Data
<ASSISTANT_TASK:>
Python Code:
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})
Explanation: Equation of State using Generalized van der Waals Theory
This notebook illustrates the generalized van der Waals theory (gvdW) for the equation of state for interacting particles. Based on the lecture notes, Properties of Molecular Fluids in Equilibrium by Sture Nordholm.
End of explanation
# 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)
Explanation: Pair potentials
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$.
Define your own potentials below
The name should start with potential_
The first parameter should be r
The docstring can be added for a nice displayname of the function (The raw python string like r"$ \mu $" is convenient when writing latex, because in normal strings the backslash acts as an escape character)
End of explanation
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
Explanation: Interaction parameter
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$.
End of explanation
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()
Explanation: Interactive EOS plot
End of explanation
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)
df = pd.read_csv("datafile.csv", delimiter="(?:\s+|,)", engine="python")
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)'
Explanation: With Data
End of explanation
<END_TASK>