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.



In [1]:

    
import inspect
from math import sqrt, pi

import numpy as np
import pandas as pd
from bokeh.io import show, output_notebook, push_notebook
from bokeh.plotting import figure
from ipywidgets import interact
from scipy.integrate import quad



In [2]:

    
output_notebook()









    





    
        
        Loading BokehJS ...

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)



In [3]:

    
# 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)

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 Equation of State for an ideal system (van 't Hoff), $\beta \mu_{ideal}= \ln (1/y)$, where $\beta = 1/k_BT$.



In [29]:

    
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 mu_ideal(n):
    return np.log(n)

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

Interactive plot



In [30]:

    
n = np.linspace(1e-4, 1e-3, 100);
def update(eps=1.0, sigma=4.0, z=1.0, Cs=0.3, potential=potential_Combined):
    D = 3.04/sqrt(Cs)
    gvdw_line.data_source.data["y"] = mu_gvdw(n, z, D, eps, sigma, potential)
    push_notebook(handle=potfig_handle)



In [31]:

    
potfig = figure(
    title="Chemical potential",
    plot_height=300, 
    plot_width=600,
    x_axis_label="Number density n", 
    y_axis_label="Potential (k T)")
ideal_line = potfig.line(n, mu_ideal(n), legend="Ideal")
gvdw_line = potfig.line(n, mu_gvdw(n, z=0, D=3.04/sqrt(0.3), eps=1, sigma=4, potential=potential_Combined), color="green", legend="GvdW")
potfig_handle = show(potfig, notebook_handle=True)



In [19]:

    
_potentials = {fname[10:]: func for fname, func in globals().items() if fname.startswith("potential_")}
interact(update, eps=(0.0, 10.0, 0.1), sigma=(0, 5, 0.1), z=(0.0, 3, 1.0), Cs=(1e-3, 1.0, 0.1), potential=_potentials);

With Data



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")