In recent years the lattice-Boltzmann method (LBM) has proven itself to be a viable way to introduce hydrodynamic interactions into coarse-grained MD simulations with moderate computational cost. The success of the GPU LBM implementation in ESPResSo and similar developments in other software packages created demand for further developments in this area. ESPResSo features two such algorithms, namely ELECTROHYDRODYNAMICS, and ELECTROKINETICS (EK). Both of these make use of the LBM and extend it to coarse-grain not only the solvent molecules but also ionic solutes. ELECTROHYDRODYNAMICS does so using a slip layer coupling for charged particles valid in the thin Debye layer (large salt concentration) limit [1], while EK explicitly treats the ionic solutes in a continuum fashion and is valid for a wide range of salt concentrations [2-4].

To make our first steps using ELECTROKINETICS we will work on one of the few systems for which analytic solutions for the electrokinetic equations exist: the slip pore geometry with a counterion-only electrolyte. The same slit pore system is also treated in the LBM tutorial, but there, the ionic species were modeled as explicit particles. For this system, the two approaches lead to exactly the same results [5]. Differences became significant for multivalent ions, very high salt concentrations, and very high surface charge, since then the mean-field approach the EK employs, is basically solving the Poisson-Nernst-Planck formalism plus the Navier-Stokes equation on a lattice. This leads to significantly different results from explicit ion approaches [6-8]. This tutorial is mainly divided into two sections.

**Theoretical Background**introduces the electrokinetic equations and the analytical solution for the slit pore system.**Simulation using ESPResSo**deals exclusively with the simulation.

If you already know about simple diffusion-migration-advection equations, continuum electrostatics, and Navier-Stokes, then you can skip the first section.

In the following, we will derive the equations modeling the time evolution of the concentrations of dissolved species as well as the solvent in the standard electrokinetic model. We do so, neglecting the salt ions' contribution to the overall mass density, which allows us to treat the dynamics of the ions and the fluid separately [8]. The solvent fluid will be modeled using the Navier-Stokes equations while we use a set of diffusion-migration-advection equations for the ionic species.

The description starts with the ionic species' concentrations $c_{k}(\vec{r}, t)$ (number density) and the associated flux densities $\vec{j}_{k}(\vec{r}, t)$, for which mass conservation holds

\begin{equation} \label{eq:model_continuity} \partial_{t} c_{k} = -\nabla \cdot\vec{j}_{k} . \end{equation}Here $\vec{r}$ denotes the spatial coordinate and $t$ the time, while $k$ enumerates the ionic species. The fluxes are caused by diffusion (due to density variations and external forces) and advection.

The advective contribution to the flux is given by

\begin{equation} \label{eq:adv_flux} \vec{j}_{k}^{\mathrm{adv.}} = c_{k} \vec{u} , \end{equation}where $\vec{u}(\vec{r}, t)$ denotes the fluid velocity (advective velocity). This equation models advection as a simple co-movement of the dissolved ions with the surrounding fluid. All inertial effects of the ions are neglected.

The diffusive behavior of the ions is best described in a reference frame co-moving with the local advective velocity $\vec{u}$. We assume that the species' relative fluxes instantaneously relax to a local equilibrium. This assumption allows us to derive the diffusive fluxes from a local free-energy density, which we define as

\begin{equation} \label{eq:model_free_energy} f \big( c_{k}(\vec{r}) \big) = \sum_{k} \underbrace{k_{\mathrm{B}}T c_{k}(\vec{r}) \left[ \log \left\lbrace \Lambda_{k}^{3} c_{k}(\vec{r}) \right\rbrace - 1 \right] }_{\mathrm{ideal~gas~contribution}} + \underbrace{z_{k} e c_{k}(\vec{r}) \Phi(\vec{r})}_{\mathrm{electrostatic~contribution}} , \end{equation}with the $\Lambda_{k}$ the species' thermal de Broglie wavelengths, $z_{k}$ their valencies, and $\Phi(\vec{r})$ the electrostatic potential. This free-energy density consists of only an ideal-gas and an electrostatic contribution. The same assumptions form the basis of Poisson-Boltzmann (PB) theory. Hence, the limitations of this model are the same as those of PB. That means this model applies to monovalent ions at low to intermediate densities and surface charges [6,7,11,12].

The species' chemical potentials $\mu_{k}$ implied by the free-energy density read

\begin{equation} \label{eq:chempot} \mu_{k}(\vec{r}) = \delta_{c_k} f(c_{k}\big( \vec r ) \big) = k_{\mathrm{B}}T \log(\Lambda_{k}^{3} c_{k}(\vec{r})) + z_{k} e \Phi(\vec{r}) . \end{equation}This in turn allows us to formulate the first-order approximation to the thermodynamic driving force as the gradient of the chemical potential, which we use to define an expression for the diffusive flux

\begin{equation} \begin{aligned} \label{eq:model_jdiff} \vec{j}_{k}^{\mathrm{diff}} &= \xi_{k} \left( -c_{k} \nabla \mu_{k} \right) = -k_{\mathrm{B}}T \xi_{k} \nabla c_{k} - \xi_{k} z_{k} e c_{k} \nabla\Phi \\ & = -D_{k} \nabla c_{k} - \xi_{k} z_{k} e c_{k} \nabla \Phi . \end{aligned} \end{equation}Here, $\xi_{k}$ and $D_{k}$ denote the mobility and the diffusion coefficient of species $k$, which are related by the Einstein-Smoluchowski relation $D_{k} / \xi_{k} = k_{\mathrm{B}}T$ [12,13].

Finally, the total number density flux combining effects of diffusion and advection reads

\begin{equation} \label{eq:model_fluxes} \vec{j}_{k} = \vec{j}_{k}^{\mathrm{diff}} + \vec{j}_{k}^{\mathrm{adv.}} = -D_{k} \nabla c_{k} - \xi_{k} z_{k} e c_{k} \nabla \Phi + c_{k} \vec{u} , \end{equation}where the first term represents Fick's law of diffusion in the absence of an external potential, the second term gives the additional flux due to an external (in this case electrostatic) potential, and the last term introduces the influence of the motion of the underlying solvent.

The dynamics of the charged species in a typical micro- or nanofluidic system are slow compared to the relaxation of the electromagnetic fields. This allows us to use stationary equations to model electromagnetic effects. We further assume that the modeled species do not carry permanent magnetic dipoles and that electric currents in the system are small. Under these conditions, the full set of Maxwell's equations reduces to the Poisson equation

\begin{equation} \label{eq:model_poisson} \nabla^2 \Phi = - \frac{1}{\varepsilon} \sum_{k} z_{k} e c_{k} = -4 \pi l_\mathrm{B} k_{\mathrm{B}}T \sum_{k} z_{k} c_{k} . \end{equation}Here $\varepsilon = \varepsilon_{0} \varepsilon_r$ denotes the product of the vacuum permittivity $\varepsilon_{0}$ with the relative permittivity of the solvent $\varepsilon_r$. We have also used the Bjerrum-length

\begin{equation} \label{eq:lbjerr} l_\mathrm{B} = \frac{e^{2}}{4 \pi \varepsilon k_{\mathrm{B}}T}. \end{equation}Finally, we have assumed that the permittivity is spatially homogeneous, since this will allow us to use efficient spectral methods to solve this equation.

As said before, since the ionic species' contribute at most a few percent to the overall mass, we can safely approximate the overall fluid's mass by the mass of the solvent (typically water) and model the solvents velocity field $\vec{u}(\vec{r}, t)$ using the Navier-Stokes equations for an isotropic, incompressible Newtonian fluid

\begin{align} \label{eq:NS} \rho \big( \partial_t \vec{u} + \left(\vec{u} \cdot \nabla \right) \vec{u} \big) &= -\nabla p_H + \eta \nabla^{2} \vec{u} + \vec{f} ,\\ \nabla \cdot \vec u &= 0 . \end{align}where $p_H$ denotes the hydrostatic pressure, $\eta$ the shear viscosity, $\rho$ the density of the fluid, and $\vec{f}$ an external body force density. For the assumption of incompressibility to hold, the Mach number needs to be small – a condition that is fulfilled for micro- and nanofluidic systems with flow velocities on the order of μm/s.

Earlier we assumed that the ions' velocity relative to the fluid instantaneously relaxes to a stationary state and that this stationary velocity is given by the product of their mobility and the force exerted on them. For this state to be stationary, all the momentum transferred into the ions by the external force needs to be dissipated into the fluid immediately. From this we can conclude that the force density acting on the fluid must read

\begin{equation} \label{eq:forcedens} \vec{f} = \sum_{k} \vec{j}^\mathrm{diff}_k / \xi_{k} = -\sum_{k} (k_\mathrm{B}T \nabla c_{k} + z_{k} e c_{k} \nabla \Phi) . \end{equation}Summarizing, the set of electrokinetic equations we solve is given by

\begin{align} \vec{j}_{k} ={}& -D_{k} \nabla c_{k} - \xi_{k} z_{k} e c_{k} \nabla \Phi + c_{k} \vec{u} ,\label{eq:flux}\\ \partial_{t} c_{k} ={}& -\nabla \cdot\vec{j}_{k} ,\label{eq:conti}\\ \nabla^2 \Phi ={}& -4 \pi l_b k_\mathrm{B}T \textstyle\sum_{k} z_{k} c_{k} ,\label{eq:poisson}\\ \rho \big( \partial_t \vec{u} + (\vec{u} \cdot \nabla ) \vec{u} \big) ={}& -\nabla p_H + \eta \nabla^{2} \vec{u} - \textstyle\sum_{k} (k_\mathrm{B}T \nabla c_{k} + z_{k} e c_{k} \nabla \Phi) ,\label{eq:ns}\\ \nabla \cdot \vec{u} ={}& 0 .\label{eq:incomp} \end{align}The slit pore system depicted in Fig. 1 consists of two like charged parallel plates of infinite extent, confining a solution of water and the plates' counterions.

Due to the net neutrality of the system and the translational symmetry in directions parallel to the plates, the potential outside the two plates must be constant. This means that using periodic or non-periodic boundary conditions makes no difference. As the system is in equilibrium in the normal direction, the electrokinetic equations for this dimension reduce to the Poisson-Boltzmann equation for the electrostatic potential, which reads \begin{align} \partial_x^2 \Phi(x) = -4 \pi \, k_\mathrm{B}T \, l_\mathrm{B} \, ze \, c_0 \cdot \exp{\left(-\frac{ze\Phi(x)}{k_\mathrm{B}T}\right)} \; , \end{align} where $x$ denotes the direction normal to the plates. The constant $c_0$ has to be chosen such that charge neutrality is fulfilled. Multiplying by $2 \partial_x \Phi(x)$ and applying the inverse chain rule further reduces the equation to first order. Subsequent separation of variables yields the solution \begin{align} \Phi(x) = -\frac{k_\mathrm{B}T}{ze} \cdot \log \left[ \frac{C^2}{8 \pi \, k_\mathrm{B}T \, l_\mathrm{B}} \cdot \cos^{-2}\left( \frac{zeC}{2 k_\mathrm{B}T} \cdot x\right) \right], \quad \left| \frac{zeC}{2 k_\mathrm{B}T} \cdot x \right| < \frac \pi 2\; . \label{eq:validation_pb_counterions} \end{align} Refer to [5] for details on this calculation. Knowing that the counterion density $c$ resembles a Boltzmann distribution in the potential $ze \Phi$ leads to the expression \begin{align} c(x) = \frac{C^2}{8 \pi \, k_\mathrm{B}T \, l_\mathrm{B}} \cdot \cos^{-2} \left( \frac{zeC}{2 k_\mathrm{B}T} \cdot x \right) \; . \label{eq:validation_pb_density} \end{align} The constant $C$ is determined by fixing the number of counterions or requiring the E-field to vanish outside the volume contained by the plates. Both yields \begin{align} C \cdot \tan \left( \frac{zed}{4 k_\mathrm{B}T} \cdot C \right) = -4 \pi \, k_\mathrm{B}T \, l_\mathrm{B} \sigma \; , \end{align} where $d$ denotes the distance between the plates and $\sigma$ their (constant) surface charge density.

Applying an electric field along one of the directions parallel to the plates does not influence the charge distribution in the normal direction, which allows us to write down the hydrodynamic equations for the parallel direction. After eliminating all terms from the Navier-Stokes Equations, which vanish due to symmetry, we are left with \begin{align} \frac{\partial_x^2 v_y(x)}{\partial x^2} = -\frac{q E C^2}{8 \, k_\mathrm{B}T \, l_\mathrm{B} \, \eta} \cdot \cos^{-2} \left( \frac{qC}{2 k_\mathrm{B}T} \cdot x \right) \; , \end{align} which yields, by means of simple integration and the application of no-slip boundary conditions \begin{align} v_y(x) = \frac{E}{2 \pi \, l_\mathrm{B} \, \eta \, ze} \cdot \log \left[ \frac{\cos \left( \frac{zeC}{2 k_\mathrm{B}T} \cdot x \right)}{\cos \left( \frac{zeC}{2 k_\mathrm{B}T} \cdot \frac d 2 \right)} \right] \; . \end{align}

With this tutorial comes a Python script `eof_analytical.py`, which evaluates all these expressions on the same grid as is used in the upcoming simulation.

To use the electrokinetics solver in ESPResSo enable the features `ELECTROKINETICS` and `EK_BOUNDARIES` during the build process.

ESPResSo does not predefine any unit system. This makes it more flexible but also requires us to spend some thought on the conversion from SI units to simulation units and back. Since most first time users have trouble with this, we will go through that process in detail here.

Important to note is that ESPResSo's unit system is nothing more than a rescaled variant of the SI system. All physical formulas you are used to in the SI system remain valid and you can use them to find relations between your units. Let's start by choosing a unit of length. Since we are going to deal with Debye layers with extensions of nanometers, a sensible choice is

\begin{align*} [x]=1\mathrm{nm}. \end{align*}The involved energies are of the magnitude of $k_{\mathrm{B}}T$. We will simulate our system at room temperature ($300\mathrm{K}$), hence we use as unit of energy \begin{align*} [E]=k_\mathrm{B}\cdot 300\mathrm{K}\approx 4.14 \cdot 10^{-21}\mathrm{J}. \end{align*}

By default ESPResSo has no concept for particle masses (but the feature can be activated). That means all particle masses are assumed to be $1\,[\mathrm{m}]$, which forces us to use the particle mass as mass unit. For this simulation we use the mass of sodium ions, which is \begin{align*} [m]=23\mathrm{u}\approx 3.82\cdot 10^{-26}\mathrm{kg}. \end{align*}

For the relation \begin{align*} E=\frac 1 2 mv^2 \end{align*}

to hold, the unit of time has to be defined so that \begin{align*} [E]=[m]\cdot\frac{[x]^2}{[t]^2}. \end{align*}

From that we get the missing unit of time \begin{align*} [t]=[x]\cdot\sqrt{\frac{[m]}{[E]}}=1\mathrm{nm}\cdot\sqrt{\frac{23\mathrm{u}}{k_B\cdot 300\mathrm{K}}}\approx 3.03760648\cdot 10^{-12}\mathrm{s}\approx 3.04\mathrm{ps}. \end{align*}

The last unit we need is the one of charge. We choose it to be the elementary charge \begin{align*} [q]=e\approx 1.60\cdot 10^{-19}\mathrm{C}. \end{align*}

We now have all the units necessary to convert our simulation parameters.

parameter | value (SI units) | value (simulation units) |
---|---|---|

channel width $d$ | $50\mathrm{nm}$ | $50\mathrm{[x]}$ |

thermal energy $k_B T$ | $k_B\cdot 300\mathrm{K}$ | $1\mathrm{[E]}$ |

Bjerrum length $l_B$ | $0.7095\mathrm{nm}$ | $0.7095\mathrm{[x]}$ |

counterion charge $q$ | $1e$ | $1\mathrm{[q]}$ |

counterion diffusion coefficient $D$ | $2.0\cdot 10^{-9}\mathrm{m^2/s}$ | $0.006075\mathrm{[x]^2/[t]}$ |

solvent density $\rho$ | $1.0\cdot 10^{3}\mathrm{kg/m^3}$ | $26.18\mathrm{[m]/[x]^3}$ |

solvent dynamic viscosity $\eta$ | $1.0\cdot 10^{-3}\mathrm{Pa}\mathrm{s}$ | $79.53\mathrm{[m]/([x][t])}$ |

external electric field $E$ | $2.585\cdot 10^{6}\mathrm{V/m}$ | $0.1\mathrm{[E]/([q][x])}$ |

ESPResSo determines the strength of the electrostatic interactions via the Bjerrum-length $l_\mathrm{B}$. That is the length for which the electrostatic interaction energy of two elementary charges equals the thermal energy

\begin{align*} k_\mathrm{B} T=\frac{e^2}{4\pi\varepsilon_0\varepsilon_r}\cdot\frac 1 {l_\mathrm{B}}. \end{align*}This yields for water at $300K$ with $\varepsilon_r = 78.54$, a Bjerrum length of $l_\mathrm{B}\approx 0.7095\mathrm{nm}$.

The script for this simulation comes with this tutorial and is called `eof_electrokinetics.py`.
All used commands are documented in the User's Guide in the section called **Electrokinetics**.

We first set up a periodic simulation box of the desired dimensions. Note that the dimensions are, of course, given in simulation units.

In [1]:

```
# Initializing espresso modules and the numpy package
import numpy as np
import espressomd
espressomd.assert_features(['CUDA', 'ELECTROKINETICS'])
from espressomd import System, electrokinetics, shapes, ekboundaries
# Set the slit pore geometry where the width is the non-periodic part of the geometry
# the padding is used to ensure that there is no field outside the slit since the
# electrostatics is used with a 3D periodic FFT solver.
box_y = 6
box_z = 6
width = 50
padding = 1
box_x = width + 2 * padding
system = System(box_l=[box_x, box_y, box_z])
```

We then store all the parameters we calculated earlier. At this point, these parameters only reside in Python variables. They will only be used by ESPResSo once they are being passed to the respective initialization functions.

In [2]:

```
# Set the electrokinetic parameters
agrid = 1.0
dt = 0.5
kT = 1.0
bjerrum_length = 0.7095
D = 0.006075
valency = 1.0
viscosity_dynamic = 79.53
density_water = 26.15
sigma = -0.05
ext_force_density = [0.0, 0.1, 0.0]
```

Before we initialize the actual electrokinetics algorithm, we need to set the time step and some other parameters that are not actually used, but would otherwise lead to error messages.

In [3]:

```
# Set the simulation parameters
system.time_step = dt
system.cell_system.skin = 0.2
system.thermostat.turn_off()
integration_length = 600
```

We can now set up the electrokinetics algorithm.
All functionality pertaining to this algorithm is available through the `electrokinetics` submodule of `espressomd`.
Please note that the fluid viscosity is specified as a kinematic viscosity, which is the dynamic viscosity divided by the fluid density.
The kinematic viscosity is also required if you initialize the pure lattice-Boltzmann method.

In [4]:

```
# Set up the (LB) electrokinetics fluid
viscosity_kinematic = viscosity_dynamic / density_water
ek = electrokinetics.Electrokinetics(agrid=agrid,
lb_density=density_water,
viscosity=viscosity_kinematic,
friction=1.0,
T=kT,
prefactor=bjerrum_length)
```

The value of the friction parameter in the previous setup command is irrelevant, since we don't include any explicit particles in our simulation, but it's needed to pass the sanity check of the LB.

Next, we set up the individual ionic species. In this case, we only set up one species of positively charged counterions. The charge density is chosen in such a way, that it will cancel out the charges of the walls which are being inserted in the step afterwards. After setting up the species, we have to add it to the electrokinetics instance.

In [5]:

```
# Set up the charged and neutral species
density_counterions = -2.0 * sigma / width
counterions = electrokinetics.Species(density=density_counterions,
D=D,
valency=valency,
ext_force_density=ext_force_density)
ek.add_species(counterions)
```

The `EKBoundary` command takes the keyword `charge_density` and the numerical charge density in simulation units as arguments.
The `shape` keyword takes an instance of a shape, which is provided by the `shapes` submodule and is the same as for the `LBBoundary` command.
Here we initialize two charged `Wall` boundaries.
To initialize the boundaries, we have to add them to the `ekboundaries` instance of the system class.
Finally, we initialize the electrokinetics algorithm with our setup by adding the electrokinetics instance as an actor to the system.

In [6]:

```
# Set up the walls confining the fluid
ek_wall_left = ekboundaries.EKBoundary(charge_density=sigma / agrid,
shape=shapes.Wall(normal=[1, 0, 0], dist=padding))
ek_wall_right = ekboundaries.EKBoundary(charge_density=sigma / agrid,
shape=shapes.Wall(normal=[-1, 0, 0], dist=-(padding + width)))
system.ekboundaries.add(ek_wall_left)
system.ekboundaries.add(ek_wall_right)
system.actors.add(ek)
```

After setting up the system, we integrate a sufficient number of time steps to relax the system into the stationary state and output the counterion density profile, the velocity profile, and the shear stress.
Since this system has translational symmetry in the x- and y-direction, we iterate over a line in the z direction and use the `species[node].quantity` command, to output local quantities.
You can instead also use the `electrokinetics.write_vtk_quantity` command to output the whole field at once in a ParaView-compatible format.

Density and velocity are not the only fields available for output. Please refer to the User's Guide for all available options.

In [7]:

```
# Integrate the system
for i in range(100):
system.integrator.run(integration_length)
print("integration step: %03i" % i, end="\r")
print("Integration finished.")
# Output
position_list = []
density_list = []
velocity_list = []
pressure_xy_list = []
for i in range(int(box_x / agrid)):
if (i * agrid >= padding) and (i * agrid < box_x - padding):
position = i * agrid - padding - width / 2.0 + agrid / 2.0
position_list.append(position)
node_idxs = (i, int(box_y / (2 * agrid)), int(box_z / (2 * agrid)))
# density
density_list.append(counterions[node_idxs].density)
# velocity
velocity_list.append(ek[node_idxs].velocity[1])
# xz component pressure tensor
pressure_xy_list.append(ek[node_idxs].pressure_tensor[0, 1])
np.savetxt("eof_electrokinetics.dat",
np.column_stack((position_list,
density_list,
velocity_list,
pressure_xy_list)),
header="#position calculated_density calculated_velocity calculated_pressure_xy")
```

With this tutorial also came a Python matplotlib script `plot.py`.
If everything went well, running this script with Python from a folder containing the output files `eof_analytical.dat` and `eof_electrokinetics.dat` should produce the result shown in Figure 2.

[1] O. A. Hickey, C. Holm, J. L. Harden and G. W. Slater *Implicit Method for Simulating Electrohydrodynamics of Polyelectrolytes* Physical Review Letters, 2010

[2] F. Capuani, I. Pagonabarraga and D. Frenkel *Discrete solution of the electrokinetic equations* The Journal of Chemical Physics, 2004

[3] G. Rempfer *A Lattice based Model for Electrokinetics* Master's thesis, University of Stuttgart, 2013

[4] G. Rempfer, G. B. Davies, C. Holm and J. de Graaf *Reducing spurious flow in simulations of electrokinetic phenomena* The Journal of Chemical Physics, 2016

[5] G. Rempfer *Lattice-Boltzmann simulations in complex geometries* Bachelor's thesis, University of Stuttgart, Institute for Computational Physics, 2010

[6] M. Deserno and C. Holm and S. May, *Fraction of Condensed Counterions around a Charged Rod: Comparison of Poisson-Boltzmann Theory and Computer Simulations* Macromolecules, 2000

[7] C. Holm, P. Kékicheff and R. Podgornik *Electrostatic Effects in Soft Matter and Biophysics* Kluwer Academic Publishers, 2001

[8] M. Deserno and C. Holm *Cell-model and Poisson-Boltzmann-theory: A brief introduction* Electrostatic Effects in Soft Matter and Biophysics, Kluwer Academic Publishers, 2001

[9] J de Graaf., G. Rempfer and C. Holm *Diffusiophoretic Self-Propulsion for Partially Catalytic Spherical Colloids* IEEE T. Nanobiosci., 2014

[10] M. Deserno *Counterion condensation for rigid linear polyelectrolytes* Universität Mainz, 2000

[11] J. de Graaf, N Boon, M Dijkstra and R. van Roij *Electrostatic interactions between Janus particles* The Journal of Chemical Physics, 2012

[12] A. Einstein *Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen* Annalen der Physik, 1905

[13] M. von Smoluchowski *Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen* Annalen der Physik, 1906

Solution from scripts/eof_analytical.py

In [8]:

```
# Set the slit pore geometry the width is the non-periodic part of the geometry
# the padding is used to ensure that there is no field outside the slit
import numpy as np
box_x = 6
box_y = 6
width = 50
padding = 6
box_z = width + 2 * padding
# Set the electrokinetic parameters
agrid = 1.0
temperature = 1.0
bjerrum_length = 0.7095
valency = 1.0
viscosity_dynamic = 79.53
density_water = 26.15
sigma = -0.05
force = 0.1
viscosity_kinematic = viscosity_dynamic / density_water
density_counterions = -2.0 * sigma / width
# Calculate the inverse length xi, which is a combination of various
# constants (xi = zeC/2kBT), with C a constant that needs to be
# solved for, or equivalently, xi needs to be solved for
# root finding function
def solve(xi=None, d=None, bjerrum_length=None, sigma=None, valency=None):
el_char = 1.0
return xi * np.tan(xi * d / 2.0) + 2.0 * np.pi * \
bjerrum_length * sigma / (valency * el_char)
size = np.pi / (2.0 * width)
pnt0 = 0.0
pntm = pnt0 + size
pnt1 = pnt0 + 1.9 * size
# the bisection scheme
tol = 1e-8
while size > tol:
val0 = solve(xi=pnt0, d=width,
bjerrum_length=bjerrum_length, sigma=sigma, valency=valency)
val1 = solve(xi=pnt1, d=width,
bjerrum_length=bjerrum_length, sigma=sigma, valency=valency)
valm = solve(xi=pntm, d=width,
bjerrum_length=bjerrum_length, sigma=sigma, valency=valency)
if (val0 < 0.0) and (val1 > 0.0):
if valm < 0.0:
pnt0 = pntm
size = size / 2.0
pntm = pnt0 + size
else:
pnt1 = pntm
size = size / 2.0
pntm = pnt1 - size
elif (val0 > 0.0) and (val1 < 0.0):
if valm < 0.0:
pnt1 = pntm
size = size / 2.0
pntm = pnt0 - size
else:
pnt0 = pntm
size = size / 2.0
pntm = pnt0 - size
else:
raise Exception(
"Bisection method fails:\nTuning of domain boundaries may be required.")
# obtain the desired xi value
xi = pntm
# function to calculate the density
def density(x=None, xi=None, bjerrum_length=None):
return (xi**2) / (2.0 * np.pi * bjerrum_length * np.cos(xi * x)**2)
# function to calculate the velocity
def velocity(x=None, xi=None, d=None, bjerrum_length=None, force=None,
viscosity_kinematic=None, density_water=None):
return force * np.log(np.cos(xi * x) / np.cos(xi * d / 2.0)) / \
(2.0 * np.pi * bjerrum_length * viscosity_kinematic * density_water)
# function to calculate the nonzero component of the pressure tensor
def pressure_tensor_offdiagonal(x=None, xi=None, bjerrum_length=None,
force=None):
return force * xi * np.tan(xi * x) / (2.0 * np.pi * bjerrum_length)
position_list = []
density_list = []
velocity_list = []
pressure_xy_list = []
for i in range(int(box_z / agrid)):
if (i * agrid >= padding) and (i * agrid < box_z - padding):
position = i * agrid - padding - width / 2.0 + agrid / 2.0
position_list.append(position)
# density
density_list.append(density(x=position, xi=xi,
bjerrum_length=bjerrum_length))
# velocity
velocity_list.append(velocity(x=position, xi=xi,
d=width, bjerrum_length=bjerrum_length,
force=force, viscosity_kinematic=viscosity_kinematic,
density_water=density_water))
# xz component pressure tensor
pressure_xy_list.append(pressure_tensor_offdiagonal(x=position, xi=xi,
bjerrum_length=bjerrum_length, force=force))
np.savetxt(
"eof_analytical.dat", np.column_stack(
(position_list, density_list, velocity_list, pressure_xy_list)),
header="#position calculated_density calculated_velocity calculated_pressure_xy")
```

Solution from scripts/plot.py

In [9]:

```
import numpy as np
import matplotlib.pyplot as plt
```

In [10]:

```
data_anal = np.loadtxt("eof_analytical.dat")
data_ek = np.loadtxt("eof_electrokinetics.dat")
fig1 = plt.figure(figsize=(16, 4))
ax = fig1.add_subplot(131)
ax.plot(data_ek[:, 0], data_ek[:, 1], 'o',
mfc='none', color="r", label="electrokinetics")
ax.plot(data_anal[:, 0], data_anal[:, 1], color="b", label="analytical")
ax.set_xlabel("x-position")
ax.set_ylabel("density")
ax.legend(loc="best")
ax = fig1.add_subplot(132)
ax.plot(data_ek[:, 0], data_ek[:, 2], 'o',
mfc='none', color="r", label="electrokinetics")
ax.plot(data_anal[:, 0], data_anal[:, 2], color="b", label="analytical")
ax.set_xlabel("x-position")
ax.set_ylabel("velocity")
ax.legend(loc="best")
ax = fig1.add_subplot(133)
ax.plot(data_ek[:, 0], data_ek[:, 3], 'o',
mfc='none', color="r", label="electrokinetics")
ax.plot(data_anal[:, 0], data_anal[:, 3], color="b", label="analytical")
ax.set_xlabel("x-position")
ax.set_ylabel("shear stress xz")
ax.legend(loc="best")
plt.show()
```