One-dimensional Poisson-Nernst-Planck systems

The matscipy.electrochemistry.PoissonNernstPlanckSystem solver provides an interface for solving classical, one-dimensional electrochemical systems with an arbitrary number of \(N\) ionic species by means of a controlled-volumes method.

The inert electrode

An archetypical system is the inert electrode at the open half-space governed by the Poisson equation and \(N\) Nernst-Planck equations with Dirichlet and zero flux’ boundary conditions as shown in Figure 1 below.

../_images/0c4c36a4ddf1f5bc4eb5cb59eafb5da547efaa5641b0678b582fa863441ed569.png

The electrochemical cell

Another archetypical system is the one-dimensional electrochemical cell governed by flux and Robin boundary conditions as shown below in Figure 2.

../_images/94cc97439dd4d31331fdce2df4e4b6c4b66b3347be84e0684cbc7b9eb2f4c3a5.png

From continuous double layer theory to discrete coordinate sets

Poisson-Nernst-Planck-System with Stern layer boundary conditions

In this example We want to fill a gap of 5 nm between gold electrodes with 0.2 wt % NaCl aqueous solution, apply a small potential difference and generate an initial configuration for LAMMPS within a cubic box.

box_Ang = np.array([50.,50.,50.]) # Angstrom
box_m = box_Ang*sc.angstrom # meters

vol_AngCube = box_Ang.prod() # Angstrom^3
vol_mCube = vol_AngCube*sc.angstrom**3

With a concentration of 0.2 wt %, we are close to NaCl’s solubility limit in water. We estimate molar concentrations and atom numbers in our box.

weight_concentration_NaCl = 0.2 # wt %
# calculate saline mass density g/cm³
saline_mass_density_kg_per_L  = 1 + weight_concentration_NaCl * 0.15 / 0.20 # g / cm^3, kg / L
# see e.g. https://www.engineeringtoolbox.com/density-aqueous-solution-inorganic-sodium-salt-concentration-d_1957.html

saline_mass_density_g_per_L = saline_mass_density_kg_per_L*sc.kilo

molar_mass_H2O = 18.015 # g / mol
molar_mass_NaCl  = 58.44 # g / mol

cNaCl_M = weight_concentration_NaCl*saline_mass_density_g_per_L/molar_mass_NaCl # mol L^-1

cNaCl_mM = np.round(cNaCl_M/sc.milli) # mM

n_NaCl = np.round(cNaCl_mM*vol_mCube*sc.value('Avogadro constant'))

Nearly \(4\,\mathrm{M}\) of saline solution correspond to about 300 ions of each species in our cubic \(5^3\,\mathrm{nm}^3\) box.

We initialize the system assuming a Stern layer of 5 Å width, close to the order of magnitude of ion size.

Solving with Robin and zero flux boundary conditions

c = [cNaCl_mM,cNaCl_mM]
z = [1,-1]
L = box_m[2]
lambda_S = 5.0e-10 # Stern layer width in m
delta_u  = 0.5

Other than above, we modify a few solver parameters as well.

pnp = PoissonNernstPlanckSystem(c,z,L, lambda_S=lambda_S, delta_u=delta_u, N=200, maxit=20, e=1e-6)

The following applies zero flux and, particularly, Robin boundary conditions as shown in Figure 2.

pnp.use_stern_layer_cell_bc()

pnp.output = True
xij = pnp.solve()

Visualization

x = np.linspace(0, L, 100)

deb = debye(c, z)

fig, (ax1,ax4) = plt.subplots(nrows=2,ncols=1,figsize=[16,10])
ax1.plot(pnp.grid/sc.nano, pnp.potential, marker='', color='tab:red', label='potential, PNP', linewidth=1, linestyle='-')

ax2 = ax1.twinx()
ax2.plot(x/sc.nano, np.ones(x.shape)*c[0], label='average concentration', color='grey', linestyle=':')
ax2.text(1, c[0]*1.5, "average concentration")

ax2.plot(pnp.grid/sc.nano, pnp.concentration[0], marker='', color='tab:orange', label='Na+, PNP', linewidth=2, linestyle='-')

ax2.plot(pnp.grid/sc.nano, pnp.concentration[1], marker='', color='tab:blue', label='Cl-, PNP', linewidth=2, linestyle='-')
ax1.axvline(x=deb/sc.nano, label='Debye length', color='grey', linestyle=':')
ax1.text(deb/sc.nano*1.2, 0.01, "Debye length", rotation=90)

ax3 = ax1.twinx()
ax3.spines["right"].set_position(("axes", 1.1))
make_patch_spines_invisible(ax3)

ax3.spines["right"].set_visible(True)
ax3.plot(pnp.grid/sc.nano, pnp.charge_density, label='charge density, PNP', color='grey', linewidth=1, linestyle='-')

ax4.semilogy(x/sc.nano, np.ones(x.shape)*c[0], label='average concentration', color='grey', linestyle=':')
ax4.text(1, c[0]*1.5, "average concentration")
ax4.semilogy(pnp.grid/sc.nano, pnp.concentration[0], marker='', color='tab:orange', label='Na+, PNP', linewidth=2, linestyle='-')
ax4.semilogy(pnp.grid/sc.nano, pnp.concentration[1], marker='', color='tab:blue', label='Cl-, PNP', linewidth=2, linestyle='-')

ax1.set_ylabel('potential $\phi$ (V)')
ax2.set_ylabel('concentration $c$ (mM)')
ax3.set_ylabel(r'charge density $\rho \> (\mathrm{C}\> \mathrm{m}^{-3})$')
ax4.set_ylabel('concentration $c$ (mM)')
ax4.set_xlabel('distance $x$ (nm)')

ax1.legend(loc='upper right',  bbox_to_anchor=(-0.1,1.02), fontsize=15)
ax2.legend(loc='center right', bbox_to_anchor=(-0.1,0.5),  fontsize=15)
ax3.legend(loc='lower right',  bbox_to_anchor=(-0.1,-0.02), fontsize=15)

fig.tight_layout()
plt.show()

../_images/b778d1f2c21598af1a33dcb30a2fb3efe78769b4b6016eb93dee00af54849cd6.png

Notice the nearly linear behavior of potential close to the electrodes arising from the Robin boundary conditions.

Sampling from continuous distributions

We want to generate three-dimensional atomistic configurations from the one-dimensional concentration distributions retrieved above. matscipy.electrochemistry provides the sampling interface continuous2discrete for this purpose.

from matscipy.electrochemistry import continuous2discrete

continuous2discrete expects a callable distribution function. Hence, we convert the physical concentration distributions into callable density functions.

from scipy import interpolate
distributions = [interpolate.interp1d(pnp.grid,pnp.concentration[i,:]) for i in range(pnp.concentration.shape[0])]

Normalization is not necessary here. Now we can sample the distribution of our \(Na^+\) ions in z-direction.

na_coordinate_sample = continuous2discrete(distribution=distributions[0], box=box_m, count=n_NaCl)

matscipy.electrochemistry provides the utility function get_histogram to count coordinate sets within bins along all three spatial dimensions.

from matscipy.electrochemistry import get_histogram
histx, histy, histz = get_histogram(na_coordinate_sample, box=box_m, n_bins=51)

For visualization purposes, we define two little helper functions.

# helper functions
def get_centers(bins):
    """Return the center of the provided bins.

    Example:
    >>> get_centers(bins=np.array([0.0, 1.0, 2.0]))
    array([ 0.5,  1.5])
    """
    bins = bins.astype(float)
    return (bins[:-1] + bins[1:]) / 2

def plot_dist(histogram, name, reference_distribution=None):
    """Plot histogram with an optional reference distribution."""
    hist, bins = histogram
    width = 1 * (bins[1] - bins[0])
    centers = get_centers(bins)

    fi, ax = plt.subplots()
    ax.bar( centers, hist, align='center', width=width, label='Empirical distribution',
            edgecolor="none")

    if reference_distribution is not None:
        ref = reference_distribution(centers)
        ref /= sum(ref)
        ax.plot(centers, ref, color='red', label='Target distribution')

    ax.set_title(name)
    ax.legend()
    ax.set_xlabel('Distance ' + name)

Eventually, we show the distribution of samples sodium ions along the non-uniform z-direction.

plot_dist(histz, 'Distribution of Na+ ions in z-direction', reference_distribution=distributions[0])

../_images/c09f77e1a29da5305fd6bb2014bdddbbb9b01a98c883aad6c62ec0fbeff97263.png

Similarly, we sample chlorine ions…

cl_coordinate_sample = continuous2discrete(distributions[1], box=box_m, count=n_NaCl)

… and show their distribution …

histx, histy, histz = get_histogram(cl_coordinate_sample, box=box_m, n_bins=51)
plot_dist(histx, 'Distribution of Cl- ions in x-direction', reference_distribution=lambda x: np.ones(x.shape)*1/box_m[0])
plot_dist(histy, 'Distribution of Cl- ions in y-direction', reference_distribution=lambda x: np.ones(x.shape)*1/box_m[1])
plot_dist(histz, 'Distribution of Cl- ions in z-direction', reference_distribution=distributions[1])

../_images/89d67ec2e4cd5f1d923d350a6de7f9499fe5cd7097ac1918681ac332ebcbbe86.png

../_images/a93d653e29a5412f0b1784c3728c6de95c8f472087dd4634747d6ced5f412d82.png

../_images/e2a9c9574bd39953506a7547b534ae7a60cc644f010c0083d578639fc785e016.png

Notice the uniform distribution along x- and y-direction.

Coordinates export

Eventually, we use ASE to export our sample to some standard format, i.e. .xyz or LAMMPS data file.

import ase
import ase.io

ASE speaks Ångström per default, thus we convert our SI units:

sample_size = int(n_NaCl)

na_atoms = ase.Atoms(
    symbols='Na'*sample_size,
    charges=[1]*sample_size,
    positions=na_coordinate_sample/sc.angstrom,
    cell=box_Ang,
    pbc=[1,1,0])

cl_atoms = ase.Atoms(
    symbols='Cl'*sample_size,
    charges=[-1]*sample_size,
    positions=cl_coordinate_sample/sc.angstrom,
    cell=box_Ang,
    pbc=[1,1,0])

system = na_atoms + cl_atoms

ase.io.write('NaCl_c_4_M_u_0.5_V_box_5x5x5nm_lambda_S_5_Ang.xyz', system, format='xyz')

# LAMMPS data format, units 'real', atom style 'full'
# before ASE 3.19.0b1, ASE had issues with exporting atom style 'full' in LAMMPS data file format, so do not expect this line to work for older ASE versions
ase.io.write('NaCl_c_4_M_u_0.5_V_box_5x5x5nm_lambda_S_5_Ang.lammps', system,
             format='lammps-data', units="real", atom_style='full')

Looking at the configuration, most ions are found within a Debye length next to the electrodes.

import matplotlib.pyplot as plt
from ase.visualize.plot import plot_atoms
fig, ax = plt.subplots()
plot_atoms(system, ax, radii=0.3, rotation=('80x,10y,10z'))
ax.set_axis_off()

../_images/b7626d5bfd59c2fcf230824a4d8982e49c9d0f24121231b8c8708b6a70cd75b8.png

References

[1] M. Z. Bazant, K. T. Chu, and B. J. Bayly, “Current-Voltage Relations for Electrochemical Thin Films,” SIAM Journal on Applied Mathematics, Jul. 2006, doi: 10.1137/040609938.

[2] O. Stern, “Zur Theorie Der Elektrolytischen Doppelschicht,” Zeitschrift für Elektrochemie und angewandte physikalische Chemie, vol. 30, no. 21–22, pp. 508–516, 1924, doi: 10.1002/bbpc.192400182.

[3] J. N. Israelachvili, Intermolecular and surface forces. Academic Press, London, 1991.

One-dimensional Poisson-Nernst-Planck systems

The inert electrode

Preparations

Solving with Dirichlet and zero flux boundary conditions

Validation: Analytical half-space solution vs. numerical finite-interval PNP system

The electrochemical cell

Solving with Dirichlet and zero flux boundary conditions

Visualization

From continuous double layer theory to discrete coordinate sets

Poisson-Nernst-Planck-System with Stern layer boundary conditions

Solving with Robin and zero flux boundary conditions

Visualization

Sampling from continuous distributions

Coordinates export

References