To work with this tutorial, you should be familiar with the following topics:
lennard_jones folder.Ionic liquid (IL) based capacitors have long been established as promising candidates in the area of efficient energy storage devices due to their extraordinary capacitance, which is why they are also termed super-capacitors. Typically, in these setups a fluid consisting of mobile charge carriers is placed between two electrodes and thus gets polarized upon application of an external potential formic an electric double layer at the interfaces [1]. Electric double-layer capacitors (EDLCs) can be constructed from electrodes of various geometries and materials where energy is stored by potential driven adsorption of counterions on the surface of the electrodes and forming the double layer. Thus, conducting, high surface area electrode materials can further maximize the energy per volume.
In this tutorial we are going to investigate the ionic layer formation between two conducting dielectric walls in presence of an applied voltage using ESPResSo's "Electrostatic layer correction with image charges" (ELC-IC) method [2]. We employ a primitive model of an aqueous salt solution, where the solvent is treated implicitly as a homogeneous dielectric background. Thus, within the limits of a mean-field treatment and for not too small slit widths, we can compare our findings with the analytical Gouy–Chapmann solution for a single charged plane since additivity is assumed to hold if the potential of both walls is screened sufficiently. While this mean-field formalism properly describes the behavior of Coulomb fluids composed of monovalent ions at low concentrations in the vicinity of weakly charged interfaces, for strongly charged systems, where correlation and finite size effects begin to dominate, Poisson–Boltzmann theory falls inadequate. Our goal in this tutorial is to demonstrate how coarse grained implicit solvent simulations can corroborate on some of the approximations and hint on extensions/deviations.
The inclusion of dielectric inhomogeneities appearing at the conducting interfaces demands to take into account image effects that involve the full solution of the Poisson equation. This is dealt with in a computationally cost-effective way using the ELC-IC method to treat the image charge effects in the presence of 2D dielectric bounding interfaces.
Charged surfaces in contact with a liquid containing free charges (ions) attract oppositely charged ions that form a diffusive electric double layer. The competition between electrostatic interactions and entropy of ions in solution determines the exact distribution of mobile ions close to charged membranes. Gouy [3] and Chapman [4] derived in the early 20th century the analytic solution for the case of a single planar wall within the mean-field approximation of the Poisson–Boltzmann (PB) equation. We will use it to describe our two-electrode system, which is justified if the electrodes are so far apart that one surface does not influence the ion distribution in front of the other surface.
In the case of a monovalent electrolyte, double integrating the PB equation and employing the corresponding boundary conditions for the charged plane located at $z=0$ yields the electrostatic potential: $$\phi(z) = -2\ln\left[ \frac{1-\tanh(\phi_\mathrm{s}/4)e^{-\kappa_\mathrm{D} z}} {1+\tanh(\phi_\mathrm{s}/4)e^{-\kappa_\mathrm{D}z}} \right].$$ Here, $\phi_\mathrm{s}=\phi(z=0)=$ const is the surface potential such that $\phi(z\rightarrow \infty)=0$. $\kappa_\mathrm{D} = \lambda_\mathrm{D}^{-1}$ is the inverse Debye screening length given by $$ \lambda_\mathrm{D} = \left(\frac{\varepsilon \, k_{\mathrm B} T}{\sum_{j = 1}^N n_j^0 \, q_j^2}\right)^{1/2}, $$ where $n_j^{(0)}$ and $q_j$ are the equilibrium number density and charge of the $j$-th ion species. For the monovalent salt this can conveniently be expressed in terms of the Bjerrum length $\ell_\mathrm{B}$ and the equilibrium salt concentration $\rho^{(0)}=\sum_j \rho_j^{(0)}$, $$ \lambda_{\mathrm D} = \left(4 \pi \, \ell_\mathrm{B} \, \rho^{(0)}/e\right)^{-1/2} . $$
The cationic and anionic density profiles then follow from the Boltzmann distribution as: $$n_{\pm}(z)=n_\pm^{(0)}\left(\frac{1\pm\gamma e^{-\kappa_\mathrm{D}z}} {1\mp\gamma e^{-\kappa_\mathrm{D}z}} \right)^2$$ Here, $\gamma$ is associated with the surface potential as $\phi_\mathrm{s}=-4\tanh^{-1}(\gamma)$. At large z, where the potential decays to zero, the ionic profiles tend to their bulk (reservoir) densities, $n_\pm(z\to\infty) = n_\pm^{(0)}$
The relation between the surface potential and the surface charge induced on the electrode is given by the Grahame Equation [5] : $$ \sigma = \sinh(\phi_\mathrm{s}/2) \sqrt{\frac{2 n_\mathrm{b}}{\pi \ell_\mathrm{B}}} $$ The latter expression thus yields the differential capacitance $C=\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}\phi_\mathrm{s}}$ within the mean-field solution for non-overlapping double layers.
In this tutorial we employ a parallel plate capacitor setup with constant potential boundary conditions which needs to be treated appropriately by the electrostatics solver. To simulate a two-dimensional partially periodic system, we combine the efficient scaling behavior of three-dimensional mesh-based solvers (typically $\mathcal{O}(N \log N)$ for P3M) with the Electrostatic Layer Correction (ELC) [1]. The latter removes the contributions from the periodic images in the non-periodic direction and its numerical cost grows linear with the number of point charges $N$, hence the performance overall depends on the underlying 3D Coulomb solver. Furthermore, ELC can be extended straightforwardly to metallic boundary conditions (or any other dielectric contrast) by using the method of image charges, which is referred to as the “Electrostatic Layer Correction with Image Charges” (ELC-IC) approach used in this tutorial.
A voltage difference can be applied between the electrodes by following considerations: The total potential drop $\Delta \phi$ across the simulated system is readily obtained from the ion distribution and integrating twice over the one-dimensional Poisson equation: $$-\varepsilon_{0}\varepsilon_{r}\phi_\mathrm{ind}(z)=\iint_{0}^{z}\rho(z^{\prime})(z-z^{\prime})dz^{\prime}$$ Here, the subscript 'ind' indicates that this is the potential due to the induced inhomogeneous charge distribution. In order to set up a constant potential difference $\Delta \phi$, a homogeneous electric field is superimposed such that $$ \Delta \phi = \Delta \phi_\mathrm{ind} + \Delta \phi_\mathrm{bat},$$ where $\Delta \phi_\mathrm{bat}$ corresponds to the potential of a (virtually) applied battery. In practice, the linear electric field in $E_z^\mathrm{(bat)}=-\phi_\mathrm{bat}/L_z$ in the $z$-direction normal to the surface that one needs to apply can be calculated straightforwardly, as the corresponding contribution from the induced charges is known: $$ E_z^\mathrm{(ind)} = \frac{1}{\varepsilon_0 \varepsilon_r L^2 L_z} P_z$$ Here, $L$ denotes the lateral system size, $L_z$ the distance between the plates and $P_z = \sum_i q_i z_i$ is the total dipole moment in $z$-direction due to the charges $q_i$. Then, to maintain $\Delta \phi$, a force $F_z^\mathrm{bat} = q E_z^\mathrm{(bat)}$ is applied on all charges. Since ELC already calculates $P_z$, the constant potential correction requires no additional computational effort.
Note: Apart from ELC-IC, ESPResSo also provides the ICC$^\star$ method [2] which employs an iterative numerical scheme with discretized surface particles to solve the boundary integrals at the dieletcric interface. The tutorial on Basic simulation of electrodes in ESPResSo part I addresses this in detail.
First we import all ESPResSo features and external modules
import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as constants
import scipy.stats
import tqdm
import espressomd
import espressomd.electrostatics
import espressomd.electrostatic_extensions
import espressomd.observables
import espressomd.accumulators
import espressomd.shapes
import espressomd.zn
espressomd.assert_features(['WCA', 'ELECTROSTATICS'])
rng = np.random.default_rng(42)
plt.rcParams.update({'font.size': 18})
We need to define system dimensions and some physical parameters related to length, time and energy scales of our system. As discussed in previous tutorials, all physical parameters are defined in terms of a length $\sigma$, mass $m$ and time $t$ and unit of charge $q$. Since we are not explicitly interested in the dynamics of the system, we set the mass to $m=1$ (particle mass). For convenience, we choose the elementary charge as fundamental unit ($q=1e$) and $\sigma = 1 \,\mathrm{nm}$.
With this we can now define the fundamental parameters of our system:
# water at room temperature
EPSILON_R = 78.4 # Relative dielectric constant of water
TEMPERATURE = 300.0 # Temperature in Kelvin
BJERRUM_LENGTH = constants.elementary_charge**2 / constants.nano / \
(4 * np.pi * constants.epsilon_0 * EPSILON_R * constants.Boltzmann * TEMPERATURE)
# BERRUM_LENGTH of water at room temperature is 0.71 nm; electrostatic prefactor passed to P3M KBT/e2
# Lennard-Jones parameters
LJ_SIGMA = 0.3 # Particle size in nanometers
LJ_EPSILON = 1.0
CONCENTRATION = 1e-2 # desired concentration in mol/l
DISTANCE = 10 # 10 Debye lengths
N_IONPAIRS = 500
POTENTIAL_DIFF = 5.0
# Elementary charge
q = 1.
types = {"Cation": 0, "Anion": 1, "Electrodes": 2}
charges = {"Cation": q, "Anion": -q}
We want to make use of the optimal performance of ESPResSo in this tutorial,
which is roughly at 1000 particles/core.
Thus, we fixed above the number of ion pairs to N_IONPAIRS = 500.
To be able to employ the analytical solution for the single plate also for the
double layer capacitor setup, the two electrodes need to be sufficiently far
away such that the additivity of the two surface potentials holds. In practice,
a separation of $d=10\lambda_\mathrm{D}$ is a good choice, represented by
DISTANCE = 10.
Our choice of $c=10\,\mathrm{mmol}$ is a compromise between a sufficiently low
concentration for the PB theory to hold and not too large distances $d$ such
that the equilibration/diffusion of the ions is sufficiently fast
(CONCENTRATION = 1e-2).
Note that in order to obtain results that we can interpret easily, we explicitly set a unit system using nanometers as length-scale above. The corresponding ion size of about 0.3 nm is a typical value for a simple salt; this, however, is in sharp contrast to the mean-field assumption of point-like ions. The latter are not easily studied within Molecular Dynamics simulations due to the required small time steps and are better suited for Monte Carlo type simulations. We instead focus here on analyzing deviations from PB theory due to the finite ion size.
Task
get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS)
that returns the lateral and normal box lengths box_l_xy and box_l_z (in
nanometers) for the given parameters.Hint: To account for the finite ion size and the wall interaction it is useful to define the effective separation $d^\prime = d-2\sigma$, such that the concentration is $\rho = N/(A \cdot d^\prime)$.
# SOLUTION CELL
def get_box_dimension(concentration, distance, n_ionpairs=N_IONPAIRS):
"""
For a given number of particles, determine the lateral area of the box
to match the desired concentration.
"""
# concentration is in mol/l, convert to 1/sigma**3
rho = concentration * (constants.Avogadro / constants.liter) * constants.nano**3
debye_length = (4. * np.pi * BJERRUM_LENGTH * rho * 2.)**(-1. / 2.) # desired Debye length in nm
l_z = distance * debye_length
box_volume = n_ionpairs / rho
area = box_volume / (l_z - 2. * LJ_SIGMA) # account for finite ion size in density calculation
l_xy = np.sqrt(area)
return l_xy, l_z
box_l_xy, box_l_z = get_box_dimension(CONCENTRATION, DISTANCE, N_IONPAIRS)
# useful quantities for the following calculations
DEBYE_LENGTH = box_l_z / DISTANCE # in units of nm
rho = N_IONPAIRS / (box_l_xy * box_l_xy * box_l_z) # in units of 1/nm^3
We now can create the ESPResSo system.
Note that for ELC to work properly, we need to add a gap of ELC_GAP in the
non-periodic direction.
The precise value highly affects the performance due to the tuning of the P3M
electrostatic solver.
For $d=10\lambda$ a gap value of $6 L_z$ is a good choice.
We also set the time-step $dt = 0.01 \tau$, which is limited by the choice of $\sigma$ and $\tau$ in the repulsive WCA interaction.
ELC_GAP = 6. * box_l_z
system = espressomd.System(box_l=[box_l_xy, box_l_xy, box_l_z + ELC_GAP])
system.time_step = 0.01
Task
types dictionary for the
particle type.# SOLUTION CELL
# Bottom wall, normal pointing in the +z direction
floor = espressomd.shapes.Wall(normal=[0, 0, 1])
c1 = system.constraints.add(
particle_type=types["Electrodes"], penetrable=False, shape=floor)
# Top wall, normal pointing in the -z direction
ceiling = espressomd.shapes.Wall(normal=[0, 0, -1], dist=-box_l_z)
c2 = system.constraints.add(
particle_type=types["Electrodes"], penetrable=False, shape=ceiling)
Task
Note, that unfavorable overlap can be avoided by placing the particles in the interval $[\sigma, d-\sigma]$ in the $z$-direction only.
# SOLUTION CELL
offset = LJ_SIGMA # avoid unfavorable overlap at close distance to the walls
init_part_btw_z1 = offset
init_part_btw_z2 = box_l_z - offset
ion_pos = np.empty((3,), dtype=float)
for i in range(N_IONPAIRS):
ion_pos[0] = rng.random(1)[0] * system.box_l[0]
ion_pos[1] = rng.random(1)[0] * system.box_l[1]
ion_pos[2] = rng.random(1)[0] * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1
system.part.add(pos=ion_pos, type=types["Cation"], q=charges["Cation"])
for i in range(N_IONPAIRS):
ion_pos[0] = rng.random(1)[0] * system.box_l[0]
ion_pos[1] = rng.random(1)[0] * system.box_l[1]
ion_pos[2] = rng.random(1)[0] * (init_part_btw_z2 - init_part_btw_z1) + init_part_btw_z1
system.part.add(pos=ion_pos, type=types["Anion"], q=charges["Anion"])
Task
Refer to the documentation to set up the WCA interaction under Non-bonded section.
# SOLUTION CELL
for t1 in types.values():
for t2 in types.values():
system.non_bonded_inter[t1, t2].wca.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA)
For the (2D+h) electrostatic with dielectrics we choose the ELC-IC with P3M.
Refer the documentation to set up ELCIC with P3M under the electrostatics section.
As later we will study different potential drops between the electrodes, write a
function that sets up the electrostatic solver for a given value
POTENTIAL_DIFF.
This function will take care of tuning the P3M and ELC parameters.
For our purposes, an accuracy of $10^{-3}$ is sufficient.
Task
setup_electrostatic_solver(potential_diff) that
returns the ELC instance.# SOLUTION CELL
def setup_electrostatic_solver(potential_diff):
delta_mid_top = -1. # (Fully metallic case both -1)
delta_mid_bot = -1.
p3m_accuracy = 1e-3
elc_accuracy = 1e-3
p3m = espressomd.electrostatics.P3M(prefactor=BJERRUM_LENGTH,
accuracy=p3m_accuracy,
mesh=[12, 12, 48], # pinned for tuning reproducibility
cao=4, # pinned for tuning reproducibility
tune=True,
verbose=False)
elc = espressomd.electrostatics.ELC(actor=p3m,
gap_size=ELC_GAP,
const_pot=True,
pot_diff=potential_diff,
maxPWerror=elc_accuracy,
delta_mid_bot=delta_mid_bot,
delta_mid_top=delta_mid_top)
return elc
# Set up steepest descent integration
MASS = 1.0
FMAX = 0.01 * LJ_SIGMA * MASS / system.time_step**2
system.integrator.set_steepest_descent(
f_max=FMAX,
gamma=30.,
max_displacement=0.01)
system.integrator.run(5000)
assert np.all(np.abs(system.part.all().f)<FMAX), "Overlap removal did not converge!"
We now switch to a Velocity Verlet integrator and set up a Langevin thermostat. Note, that we only analyze static properties, thus the damping and temperature chosen here only determine the simulation time towards the equilibrium distribution.
# Equilibration parameters
system.time_step = 0.003 # this time step limits particle movement to ~5% sigma per integration step
system.integrator.set_vv()
system.thermostat.set_langevin(kT=1., gamma=2., seed=42)
system.electrostatics.solver = setup_electrostatic_solver(POTENTIAL_DIFF)
times = []
e_kin = []
for i in tqdm.trange(100):
energy = system.analysis.energy()
e_kin.append(energy['kinetic'])
times.append(system.time)
system.integrator.run(10)
100%|██████████| 100/100 [00:07<00:00, 13.57it/s]
# Plot the convergence of the total energy
plt.figure(figsize=(10, 6))
plt.plot(times, len(e_kin) * [N_IONPAIRS * 2 * 3. / 2.], label='heat bath')
plt.plot(times, e_kin, label='kinetic energy')
plt.xlabel('Simulation time')
plt.ylabel('Energy')
plt.legend()
plt.show()
Convergence after $t \sim 5$ time units.
Now we let ions diffuse to the electrodes.
N_SAMPLES_EQUIL = 20
STEPS_PER_EQUIL = 200
system.thermostat.set_langevin(kT=1.0, gamma=0.1, seed=42)
for tm in tqdm.trange(N_SAMPLES_EQUIL):
system.integrator.run(STEPS_PER_EQUIL)
100%|██████████| 20/20 [00:23<00:00, 1.19s/it]
We now need to set up an espressomd.observables.DensityProfile observable to calculate the anion and cation density profiles.
The time average is obtained through a espressomd.accumulators.MeanVarianceCalculator.
Task
setup_densityprofile_accumulators(bin_width) that returns the
bin_centers and the accumulators for both ion species in the $z$-range $[0,d]$.
Since we are not estimating errors in this tutorial, the choice of delta_N is
rather arbitrary and does not affect the results. In practice, a typical value is
delta_N=20.# SOLUTION CELL
def setup_densityprofile_accumulators(bin_width):
cations = system.part.select(type=types["Cation"])
anions = system.part.select(type=types["Anion"])
n_z_bins = int(np.round((system.box_l[2] - ELC_GAP) / bin_width))
density_profile_cation = espressomd.observables.DensityProfile(
ids=cations.id, n_x_bins=1, n_y_bins=1, n_z_bins=n_z_bins, min_x=0, min_y=0, min_z=0,
max_x=system.box_l[0], max_y=system.box_l[1], max_z=system.box_l[2] - ELC_GAP)
density_accumulator_cation = espressomd.accumulators.MeanVarianceCalculator(
obs=density_profile_cation, delta_N=20)
density_profile_anion = espressomd.observables.DensityProfile(
ids=anions.id, n_x_bins=1, n_y_bins=1, n_z_bins=n_z_bins, min_x=0, min_y=0, min_z=0,
max_x=system.box_l[0], max_y=system.box_l[1], max_z=system.box_l[2] - ELC_GAP)
density_accumulator_anion = espressomd.accumulators.MeanVarianceCalculator(
obs=density_profile_anion, delta_N=20)
zs = density_profile_anion.bin_centers()[0, 0, :, 2]
return zs, density_accumulator_cation, density_accumulator_anion
zs, density_accumulator_cation, density_accumulator_anion = setup_densityprofile_accumulators(
bin_width=DEBYE_LENGTH / 10.)
Now we take some measurement sampling the density profiles.
N_SAMPLES_PROD = 20
STEPS_PER_SAMPLE = 200
# Add the accumulators
system.auto_update_accumulators.add(density_accumulator_cation)
system.auto_update_accumulators.add(density_accumulator_anion)
for tm in tqdm.trange(N_SAMPLES_PROD):
system.integrator.run(STEPS_PER_SAMPLE)
cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]
anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]
100%|██████████| 20/20 [00:23<00:00, 1.17s/it]
Since we assume additivity, the total ion density follows from $$ \rho (z) = \rho_+(z) - \rho_+ (d-z) + \rho_-(z) - \rho_-(d-z) .$$
def gouy_chapman_potential(x, debye_length, phi_0):
kappa = 1. / debye_length
return 2. * np.log((1. + np.tanh(1. / 4. * phi_0 * np.exp(-kappa * x))) \
/ (1. - np.tanh(1. / 4. * phi_0 * np.exp(-kappa * x))))
def gouy_chapman_density(x, c0, debye_length, phi_0):
phi = gouy_chapman_potential(x, debye_length, phi_0)
return c0 / 2. * np.exp(-phi)
fig, (ax1, ax2, ax3) = plt.subplots(figsize=(16, 4), nrows=1, ncols=3, sharey=True)
fig.subplots_adjust(wspace=0)
x = np.linspace(LJ_SIGMA, box_l_z-LJ_SIGMA, 100)
anion_profile_analytic = (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.) \
+ gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.))/2.
cation_profile_analytic = (gouy_chapman_density(x, CONCENTRATION, DEBYE_LENGTH,POTENTIAL_DIFF/2.) \
+ gouy_chapman_density(box_l_z-LJ_SIGMA-x, CONCENTRATION, DEBYE_LENGTH,-POTENTIAL_DIFF/2.))/2.
ax1.set_title('Cation')
ax2.set_title('Anion')
ax3.set_title('Total')
ax1.plot(x, cation_profile_analytic, label='analytic')
ax2.plot(x, anion_profile_analytic, label='analytic')
ax3.plot(x, cation_profile_analytic + anion_profile_analytic, label='analytic')
ax1.plot(zs[1:-1], cation_profile_mean[1:-1], 'o', mfc='none', label='simulation')
ax2.plot(zs[1:-1], anion_profile_mean[1:-1], 'o', mfc='none', label='simulation')
ax3.plot(zs[1:-1], cation_profile_mean[1:-1] + anion_profile_mean[1:-1], 'o', mfc='none', label='simulation')
ax1.legend(loc='upper center')
ax2.legend(loc='upper center')
ax3.legend(loc='upper center')
ax2.set_xlabel(r'$z\,\mathrm{[nm]}$')
ax1.set_ylabel(r'$\rho(z)\,\mathrm{[mol/l]}$')
plt.show()
We see good agreement between our simulation and the meanfield solution of Guy and Chapman. Low density and reasonably low potential make the assumptions of the analytical approach justified.
We now check how well the surface charge agrees with Grahame's equation. To this end we calculate $$\sigma = \int_0^{d/2} \rho(z) \,\mathrm{d}z .$$
sigma_left = np.sum((cation_profile_mean-anion_profile_mean)[:len(zs)//2]) * (zs[1] - zs[0])
sigma_right = np.sum((cation_profile_mean-anion_profile_mean)[len(zs)//2:]) * (zs[1] - zs[0])
def grahame_sigma(phi):
return np.sinh(phi / 4.) * np.sqrt(2. * rho / (np.pi * BJERRUM_LENGTH))
sigma_grahame = grahame_sigma(POTENTIAL_DIFF)
print(f'simulation: {sigma_right:.3f} e/nm^2')
print(f'grahame: {sigma_grahame:.3f} e/nm^2')
print(f'relative deviation: {abs(1. - sigma_right/sigma_grahame) * 100.:.0f}%')
simulation: 0.075 e/nm^2 grahame: 0.117 e/nm^2 relative deviation: 36%
The electric field is readily obtained from the integral $$E(z) = \int_0^{z} \frac{1}{\varepsilon_0 \varepsilon_r} \rho(z^\prime) \,\mathrm{d}z^\prime .$$
# plot the electric field
fig, ax = plt.subplots(figsize=(10, 6))
dz_SI = (zs[1] - zs[0]) * constants.nano
chargedensity = (cation_profile_mean - anion_profile_mean) * constants.elementary_charge / constants.nano**3
E_SI = 1. / (EPSILON_R * constants.epsilon_0) * np.cumsum(chargedensity * dz_SI)
# integration constant: zero field in the center
E_SI -= E_SI.min()
E = E_SI / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE) / constants.nano)
ax2 = plt.twinx()
ax.plot(zs, E_SI)
ax2.plot(zs, E)
ax.set_xlabel(r'$z\,\mathrm{[nm]}$')
ax.set_ylabel(r'$E_\mathrm{ind}\,\mathrm{[V/m]}$')
ax2.set_ylabel(r'$E_\mathrm{ind}\,\mathrm{[(k_\mathrm{B}T/e)/nm]}$')
plt.show()
We see that the electric field reduces to 0 in the middle of the channel, justifying the assumption that the two electrodes are far enough apart to not influence each other.
The electric potential can be calculated from $\phi(z) = \int_0^z -E(z^\prime)\,\mathrm{d}z^\prime$.
# plot the elecrostatic potential
fig, ax = plt.subplots(figsize=(10, 6))
ax2 = ax.twinx()
phi_SI = -np.cumsum(E_SI * dz_SI)
phi = phi_SI * (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE))
ax.plot(zs, phi_SI)
ax.set_xlabel(r'$z\,\mathrm{[nm]}$')
ax.set_ylabel(r'$\phi\,[V]$')
ax2.set_ylabel(r'$\phi\,[k_\mathrm{B}T/e]$')
ax2.axhline(-5, ls='--', color='k')
ax2.axhline(0, ls='--', color='k')
ax.set_xlim(0, 10. * DEBYE_LENGTH)
plt.show()
measured_potential_difference = -(phi[-1] - phi[0])
print(f'applied voltage: {POTENTIAL_DIFF:.2f} V')
print(f'measured voltage: {measured_potential_difference:.2f} V')
print(f'relative deviation: {abs(1. - measured_potential_difference / POTENTIAL_DIFF) * 100:.0f}%')
applied voltage: 5.00 V measured voltage: 4.80 V relative deviation: 4%
With the above knowledge, we can now assess the differential capacitance of the system, by changing the applied voltage difference and determining the corresponding surface charge density. We will use ZnDraw to visualize our system:
color = {
types["Cation"]: "#ff0000", #red
types["Anion"]: "#030ffc" #blue
}
radii = {
types["Cation"]: LJ_SIGMA,
types["Anion"]: LJ_SIGMA
}
vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)
#vis.draw_constraints([floor, ceiling])
# note: you may need to zoom out since the ELC gap region takes a significant portion of
# the box along the non-periodic direction. The particles are only in the smaller subsystem.
# note: The particles are shown bigger for visualization purpose.
Do the sampling from high to low potential:
sigma_vs_phi = []
MIN_PHI = 0.5
MAX_PHI = 10
N_PHI = 7
N_SAMPLES_EQUIL_CAP = 15
N_SAMPLES_CAP = 5
# sample from high to low potential to improve sampling
for potential_diff in tqdm.tqdm(np.linspace(MIN_PHI, MAX_PHI, N_PHI)[::-1]):
system.auto_update_accumulators.clear()
system.electrostatics.solver = setup_electrostatic_solver(potential_diff)
for i in range(N_SAMPLES_EQUIL_CAP * STEPS_PER_SAMPLE//50):
system.integrator.run(50)
vis.update()
sigmas = []
for tm in range(N_SAMPLES_CAP):
zs, density_accumulator_cation, density_accumulator_anion = \
setup_densityprofile_accumulators(bin_width=DEBYE_LENGTH / 10.)
system.auto_update_accumulators.clear()
system.auto_update_accumulators.add(density_accumulator_cation)
system.auto_update_accumulators.add(density_accumulator_anion)
for j in range(int(STEPS_PER_SAMPLE/50)):
system.integrator.run(50)
vis.update()
cation_profile_mean = density_accumulator_cation.mean()[0, 0, :]
anion_profile_mean = density_accumulator_anion.mean()[0, 0, :]
sigmas.append(np.sum((cation_profile_mean - anion_profile_mean)[:int(len(zs) / 2.)]) * (zs[1] - zs[0]))
sigma_vs_phi.append([potential_diff, np.mean(sigmas), scipy.stats.sem(sigmas)])
100%|██████████| 7/7 [02:41<00:00, 23.04s/it]
fig, ax = plt.subplots(figsize=(10, 6))
sigma_vs_phi = np.array(sigma_vs_phi)
x = np.linspace(0, sigma_vs_phi[:,0].max())
phi_SI = sigma_vs_phi[:,0] / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE))
plt.errorbar(-sigma_vs_phi[:,1] * constants.elementary_charge / constants.nano**2,
phi_SI, xerr=sigma_vs_phi[:,2] * scipy.constants.elementary_charge / scipy.constants.nano**2,
fmt='o',label='Simulation')
plt.plot(grahame_sigma(x) * constants.elementary_charge / constants.nano**2,
x / (constants.elementary_charge / (constants.Boltzmann * TEMPERATURE)), label='Grahame')
x = np.linspace(0, ax.get_ylim()[1])
plt.plot(EPSILON_R * constants.epsilon_0 * x / 2. / (DEBYE_LENGTH * constants.nano), x, label='linear PB',
ls='--')
plt.xlabel(r'$\sigma\,\mathrm{[C/m^2]}$')
plt.ylabel(r'$\phi_\mathrm{s}\,\mathrm{[V]}$')
plt.legend()
plt.show()
For small potential drops, one observes the expected Poisson–Boltzmann behavior. It also agrees with the linearized solution $\sigma(\phi_\mathrm{s}) = \varepsilon_r\varepsilon_0 \frac{\phi_\mathrm{s}}{2 \lambda_\mathrm{D}}$. However, we observe already for potentials $\sim 0.1\,\mathrm{V} = 4\,k_\mathrm{B}T / e$ a significant deviation which can be attributed to the fact that our ions are of finite size and thus the surface charge saturates.
[1] Conway, B. E. Electrochemical Supercapacitors; Springer US: Boston, MA, 1999. https://doi.org/10.1007/978-1-4757-3058-6.
[2] Tyagi, S.; Arnold, A.; Holm, C. Electrostatic Layer Correction with Image Charges: A Linear Scaling Method to Treat Slab 2D+h Systems with Dielectric Interfaces. J. Chem. Phys. 2008, 129 (20), 204102. https://doi.org/10.1063/1.3021064.
[3] Gouy, G. Constitution of the Electric Charge at the Surface of an Electrolyte. J. Phys. 1910, 9 (4), 457–467.
[4] Chapman, D. L. A Contribution to the Theory of Electrocapillarity. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 1913, 25 (148), 475. https://doi.org/10.1080/14786440408634187.
[5] Grahame, D. C. The Electrical Double Layer and the Theory of Electrocapillarity. Chem. Rev. 1947, 41 (3), 441–501. https://doi.org/10.1021/cr60130a002.
[6] Tyagi, S.; Süzen, M.; Sega, M.; Barbosa, M.; Kantorovich, S. S.; Holm, C. An Iterative, Fast, Linear-Scaling Method for Computing Induced Charges on Arbitrary Dielectric Boundaries. J. Chem. Phys. 2010, 132 (15), 154112. https://doi.org/10.1063/1.3376011.