Skip to content

Physics API

Poisson-Boltzmann solver

edl_ml.physics.pb

Poisson-Boltzmann solver for the diffuse layer of a symmetric electrolyte.

The full nonlinear one-dimensional Poisson-Boltzmann equation for a :math:z:z symmetric electrolyte reads

.. math::

\frac{d^2 \psi}{dx^2} = \frac{2 z e n_0}{\epsilon_r \epsilon_0}
                         \sinh\!\left(\frac{z e \psi}{k_B T}\right),

with :math:\psi the electrostatic potential, :math:n_0 the bulk number density of each ion, :math:z the valence, :math:T the temperature and :math:\epsilon_r the relative permittivity of the solvent.

For a one-dimensional symmetric electrolyte in contact with a planar electrode there is a closed-form first integral

.. math::

\frac{d\psi}{dx} = -\,\mathrm{sgn}(\psi)\;\frac{2 k_B T \kappa}{z e}
                     \sinh\!\left(\frac{z e \psi}{2 k_B T}\right),

which lets us evaluate the diffuse-layer surface charge density analytically via Gauss's law while still integrating the full nonlinear ODE numerically to recover the spatial profile.

PBParameters dataclass

Input parameters to the Poisson-Boltzmann solver.

Parameters:

Name Type Description Default
concentration_mol_l float

Bulk electrolyte concentration in mol/L.

required
valence int

Ionic valence :math:z for the symmetric electrolyte.

1
temperature_k float

Absolute temperature in Kelvin.

ROOM_TEMPERATURE
relative_permittivity float

Relative permittivity of the solvent.

WATER_PERMITTIVITY
psi_diffuse_v float

Potential at the Stern / diffuse boundary, in volts.

0.05
domain_debye_lengths float

Length of the integration domain measured in Debye lengths.

30.0
n_points int

Number of points used to represent the solution on export.

400
Source code in src/edl_ml/physics/pb.py
@dataclass(frozen=True, slots=True)
class PBParameters:
    """Input parameters to the Poisson-Boltzmann solver.

    Parameters
    ----------
    concentration_mol_l
        Bulk electrolyte concentration in mol/L.
    valence
        Ionic valence :math:`z` for the symmetric electrolyte.
    temperature_k
        Absolute temperature in Kelvin.
    relative_permittivity
        Relative permittivity of the solvent.
    psi_diffuse_v
        Potential at the Stern / diffuse boundary, in volts.
    domain_debye_lengths
        Length of the integration domain measured in Debye lengths.
    n_points
        Number of points used to represent the solution on export.
    """

    concentration_mol_l: float
    valence: int = 1
    temperature_k: float = ROOM_TEMPERATURE
    relative_permittivity: float = WATER_PERMITTIVITY
    psi_diffuse_v: float = 0.05
    domain_debye_lengths: float = 30.0
    n_points: int = 400

    def __post_init__(self) -> None:
        if self.concentration_mol_l <= 0:
            raise ValueError("concentration must be positive")
        if self.valence <= 0:
            raise ValueError("valence must be a positive integer")
        if self.temperature_k <= 0:
            raise ValueError("temperature must be positive")
        if self.relative_permittivity <= 0:
            raise ValueError("relative_permittivity must be positive")
        if self.domain_debye_lengths <= 0:
            raise ValueError("domain_debye_lengths must be positive")
        if self.n_points < 10:
            raise ValueError("n_points must be at least 10")

PBResult dataclass

Output of the Poisson-Boltzmann solver.

Attributes:

Name Type Description
x_m NDArray[float64]

Position array in metres, measured from the Stern / diffuse boundary.

psi_v NDArray[float64]

Electrostatic potential in volts.

field_v_m NDArray[float64]

Electric field :math:-d\psi/dx in V/m.

cation_density_m3 NDArray[float64]

Cation number density as a function of position, in 1/m³.

anion_density_m3 NDArray[float64]

Anion number density as a function of position, in 1/m³.

surface_charge_c_m2 float

Diffuse-layer surface charge density, in C/m² (Grahame equation).

debye_length_m float

Debye screening length, in metres.

Source code in src/edl_ml/physics/pb.py
@dataclass(frozen=True, slots=True)
class PBResult:
    """Output of the Poisson-Boltzmann solver.

    Attributes
    ----------
    x_m
        Position array in metres, measured from the Stern / diffuse boundary.
    psi_v
        Electrostatic potential in volts.
    field_v_m
        Electric field :math:`-d\\psi/dx` in V/m.
    cation_density_m3
        Cation number density as a function of position, in 1/m³.
    anion_density_m3
        Anion number density as a function of position, in 1/m³.
    surface_charge_c_m2
        Diffuse-layer surface charge density, in C/m² (Grahame equation).
    debye_length_m
        Debye screening length, in metres.
    """

    x_m: NDArray[np.float64]
    psi_v: NDArray[np.float64]
    field_v_m: NDArray[np.float64]
    cation_density_m3: NDArray[np.float64]
    anion_density_m3: NDArray[np.float64]
    surface_charge_c_m2: float
    debye_length_m: float

debye_length

debye_length(concentration_mol_l: float, valence: int = 1, temperature_k: float = ROOM_TEMPERATURE, relative_permittivity: float = WATER_PERMITTIVITY) -> float

Return the Debye screening length for a symmetric electrolyte.

Parameters:

Name Type Description Default
concentration_mol_l float

Bulk concentration of each ion in mol/L.

required
valence int

Ionic valence.

1
temperature_k float

Temperature in Kelvin.

ROOM_TEMPERATURE
relative_permittivity float

Relative permittivity of the solvent.

WATER_PERMITTIVITY

Returns:

Type Description
float

Debye length in metres.

Notes

The Debye length in a symmetric electrolyte is

.. math::

\kappa^{-1} = \sqrt{\frac{\epsilon_r \epsilon_0 k_B T}
                             {2 N_A e^2 z^2 c}}

where :math:c is the bulk concentration in mol/m³.

Source code in src/edl_ml/physics/pb.py
def debye_length(
    concentration_mol_l: float,
    valence: int = 1,
    temperature_k: float = ROOM_TEMPERATURE,
    relative_permittivity: float = WATER_PERMITTIVITY,
) -> float:
    """Return the Debye screening length for a symmetric electrolyte.

    Parameters
    ----------
    concentration_mol_l
        Bulk concentration of each ion in mol/L.
    valence
        Ionic valence.
    temperature_k
        Temperature in Kelvin.
    relative_permittivity
        Relative permittivity of the solvent.

    Returns
    -------
    float
        Debye length in metres.

    Notes
    -----
    The Debye length in a symmetric electrolyte is

    .. math::

        \\kappa^{-1} = \\sqrt{\\frac{\\epsilon_r \\epsilon_0 k_B T}
                                     {2 N_A e^2 z^2 c}}

    where :math:`c` is the bulk concentration in mol/m³.
    """
    c_m3 = concentration_mol_l * 1e3
    denom = 2.0 * AVOGADRO * ELEMENTARY_CHARGE**2 * valence**2 * c_m3
    num = relative_permittivity * VACUUM_PERMITTIVITY * BOLTZMANN * temperature_k
    return float(np.sqrt(num / denom))

solve_poisson_boltzmann

solve_poisson_boltzmann(params: PBParameters) -> PBResult

Solve the nonlinear Poisson-Boltzmann equation for the diffuse layer.

The first-integral formulation of the Gouy-Chapman equation is integrated forward from x = 0 with psi(0) = psi_d. This formulation is equivalent to the full second-order PB equation under the boundary condition psi → 0 as x → infinity, and is numerically stable because the bulk fixed point is attracting on forward integration. The surface charge density is reported via the Grahame equation.

Parameters:

Name Type Description Default
params PBParameters

Solver parameters.

required

Returns:

Type Description
PBResult

Potential, field and ion density profiles with the derived surface charge density and Debye length.

Source code in src/edl_ml/physics/pb.py
def solve_poisson_boltzmann(params: PBParameters) -> PBResult:
    """Solve the nonlinear Poisson-Boltzmann equation for the diffuse layer.

    The first-integral formulation of the Gouy-Chapman equation is integrated
    forward from ``x = 0`` with ``psi(0) = psi_d``. This formulation is
    equivalent to the full second-order PB equation under the boundary
    condition ``psi → 0`` as ``x → infinity``, and is numerically stable
    because the bulk fixed point is attracting on forward integration. The
    surface charge density is reported via the Grahame equation.

    Parameters
    ----------
    params
        Solver parameters.

    Returns
    -------
    PBResult
        Potential, field and ion density profiles with the derived surface
        charge density and Debye length.
    """
    kd = debye_length(
        params.concentration_mol_l,
        params.valence,
        params.temperature_k,
        params.relative_permittivity,
    )
    thermal_scale = BOLTZMANN * params.temperature_k / (params.valence * ELEMENTARY_CHARGE)
    length = params.domain_debye_lengths * kd
    kappa = 1.0 / kd

    psi_d = params.psi_diffuse_v
    x_eval = np.linspace(0.0, length, params.n_points)
    sol = solve_ivp(
        _first_integral_rhs,
        (0.0, length),
        np.array([psi_d]),
        args=(kappa, thermal_scale),
        method="LSODA",
        rtol=1e-10,
        atol=1e-14,
        t_eval=x_eval,
    )
    if not sol.success:
        raise RuntimeError(f"PB solver failed: {sol.message}")

    psi = sol.y[0]
    # Derivative from the same first-integral relation.
    reduced_half = np.clip(psi / (2.0 * thermal_scale), -40.0, 40.0)
    dpsi = -np.sign(psi) * 2.0 * thermal_scale * kappa * np.sinh(np.abs(reduced_half))
    field = -dpsi

    c_m3 = params.concentration_mol_l * 1e3 * AVOGADRO
    reduced = np.clip(psi / thermal_scale, -80.0, 80.0)
    cation = c_m3 * np.exp(-reduced)
    anion = c_m3 * np.exp(reduced)

    sigma = _grahame_surface_charge(
        psi_d,
        params.concentration_mol_l,
        params.valence,
        params.temperature_k,
        params.relative_permittivity,
    )

    return PBResult(
        x_m=x_eval,
        psi_v=psi,
        field_v_m=field,
        cation_density_m3=cation,
        anion_density_m3=anion,
        surface_charge_c_m2=sigma,
        debye_length_m=kd,
    )

Gouy-Chapman-Stern model

edl_ml.physics.gcs

Gouy-Chapman-Stern composite double-layer model.

The Stern layer behaves as a molecular condenser with thickness equal to the closest approach of a hydrated ion to the electrode surface, and a dielectric permittivity significantly reduced compared with bulk water due to the orientation of water molecules near the interface.

The total differential capacitance is obtained by placing the Stern and diffuse-layer capacitances in series,

.. math::

\frac{1}{C_\text{dl}} = \frac{1}{C_H} + \frac{1}{C_d},

with :math:C_H = \epsilon_r^H \epsilon_0 / d_H the Helmholtz (Stern) capacitance and :math:C_d the diffuse-layer capacitance from the Gouy-Chapman model.

GCSParameters dataclass

Parameters defining a Gouy-Chapman-Stern electrode-electrolyte interface.

Parameters:

Name Type Description Default
concentration_mol_l float

Bulk electrolyte concentration, mol/L.

required
valence int

Ionic valence :math:z for a symmetric electrolyte.

1
temperature_k float

Temperature in Kelvin.

ROOM_TEMPERATURE
stern_thickness_m float

Thickness of the Stern (inner Helmholtz) layer, in metres. Typically 3–6 Å.

3e-10
stern_permittivity float

Relative permittivity of the Stern layer, dimensionless. For water the oriented interfacial value is ~6.

6.0
bulk_permittivity float

Relative permittivity of the bulk solvent, dimensionless.

WATER_PERMITTIVITY
Source code in src/edl_ml/physics/gcs.py
@dataclass(frozen=True, slots=True)
class GCSParameters:
    """Parameters defining a Gouy-Chapman-Stern electrode-electrolyte interface.

    Parameters
    ----------
    concentration_mol_l
        Bulk electrolyte concentration, mol/L.
    valence
        Ionic valence :math:`z` for a symmetric electrolyte.
    temperature_k
        Temperature in Kelvin.
    stern_thickness_m
        Thickness of the Stern (inner Helmholtz) layer, in metres. Typically
        3–6 Å.
    stern_permittivity
        Relative permittivity of the Stern layer, dimensionless. For water the
        oriented interfacial value is ~6.
    bulk_permittivity
        Relative permittivity of the bulk solvent, dimensionless.
    """

    concentration_mol_l: float
    valence: int = 1
    temperature_k: float = ROOM_TEMPERATURE
    stern_thickness_m: float = 3e-10
    stern_permittivity: float = 6.0
    bulk_permittivity: float = WATER_PERMITTIVITY

    def __post_init__(self) -> None:
        if self.concentration_mol_l <= 0:
            raise ValueError("concentration must be positive")
        if self.valence <= 0:
            raise ValueError("valence must be a positive integer")
        if self.temperature_k <= 0:
            raise ValueError("temperature must be positive")
        if self.stern_thickness_m <= 0:
            raise ValueError("stern_thickness must be positive")
        if self.stern_permittivity <= 0:
            raise ValueError("stern_permittivity must be positive")
        if self.bulk_permittivity <= 0:
            raise ValueError("bulk_permittivity must be positive")

stern_capacitance

stern_capacitance(params: GCSParameters) -> float

Return the Stern (Helmholtz) capacitance per unit area in F/m².

Parameters:

Name Type Description Default
params GCSParameters

GCS parameters.

required

Returns:

Type Description
float

:math:C_H = \epsilon_r^H \epsilon_0 / d_H.

Source code in src/edl_ml/physics/gcs.py
def stern_capacitance(params: GCSParameters) -> float:
    """Return the Stern (Helmholtz) capacitance per unit area in F/m².

    Parameters
    ----------
    params
        GCS parameters.

    Returns
    -------
    float
        :math:`C_H = \\epsilon_r^H \\epsilon_0 / d_H`.
    """
    return float(params.stern_permittivity * VACUUM_PERMITTIVITY / params.stern_thickness_m)

diffuse_capacitance

diffuse_capacitance(params: GCSParameters, psi_diffuse_v: float | NDArray[float64]) -> float | NDArray[np.float64]

Gouy-Chapman diffuse-layer differential capacitance.

Parameters:

Name Type Description Default
params GCSParameters

GCS parameters. Only the bulk-phase properties are used.

required
psi_diffuse_v float | NDArray[float64]

Potential at the Stern / diffuse boundary, V.

required

Returns:

Type Description
float or ndarray

:math:C_d = \epsilon_r \epsilon_0 \kappa \cosh(z e \psi_d / 2 k_B T), in F/m².

Notes

The closed form is obtained by differentiating the Grahame equation with respect to the diffuse-layer potential.

Source code in src/edl_ml/physics/gcs.py
def diffuse_capacitance(
    params: GCSParameters,
    psi_diffuse_v: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
    """Gouy-Chapman diffuse-layer differential capacitance.

    Parameters
    ----------
    params
        GCS parameters. Only the bulk-phase properties are used.
    psi_diffuse_v
        Potential at the Stern / diffuse boundary, V.

    Returns
    -------
    float or ndarray
        :math:`C_d = \\epsilon_r \\epsilon_0 \\kappa
        \\cosh(z e \\psi_d / 2 k_B T)`, in F/m².

    Notes
    -----
    The closed form is obtained by differentiating the Grahame equation with
    respect to the diffuse-layer potential.
    """
    from edl_ml.physics.pb import debye_length

    kd = debye_length(
        params.concentration_mol_l,
        params.valence,
        params.temperature_k,
        params.bulk_permittivity,
    )
    eps = params.bulk_permittivity * VACUUM_PERMITTIVITY
    kappa = 1.0 / kd
    thermal = BOLTZMANN * params.temperature_k
    argument = params.valence * ELEMENTARY_CHARGE * psi_diffuse_v / (2.0 * thermal)
    return eps * kappa * np.cosh(argument)

total_capacitance

total_capacitance(params: GCSParameters, psi_diffuse_v: float | NDArray[float64]) -> float | NDArray[np.float64]

Series combination of the Stern and diffuse-layer capacitances.

Parameters:

Name Type Description Default
params GCSParameters

GCS parameters.

required
psi_diffuse_v float | NDArray[float64]

Potential at the Stern / diffuse boundary, V.

required

Returns:

Type Description
float or ndarray

:math:C_\text{dl} in F/m².

Source code in src/edl_ml/physics/gcs.py
def total_capacitance(
    params: GCSParameters,
    psi_diffuse_v: float | NDArray[np.float64],
) -> float | NDArray[np.float64]:
    """Series combination of the Stern and diffuse-layer capacitances.

    Parameters
    ----------
    params
        GCS parameters.
    psi_diffuse_v
        Potential at the Stern / diffuse boundary, V.

    Returns
    -------
    float or ndarray
        :math:`C_\\text{dl}` in F/m².
    """
    c_h = stern_capacitance(params)
    c_d = diffuse_capacitance(params, psi_diffuse_v)
    return 1.0 / (1.0 / c_h + 1.0 / c_d)

gouy_chapman_stern

gouy_chapman_stern(params: GCSParameters, electrode_potentials_v: NDArray[float64], *, max_iter: int = 80, tol: float = 1e-10) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]

Solve the self-consistent Gouy-Chapman-Stern problem.

The electrode potential :math:E relative to the point of zero charge is split across the Stern and diffuse layers,

.. math::

E = \psi_H + \psi_d,\qquad \psi_H = \sigma / C_H,

with :math:\sigma given by the Grahame equation as a function of :math:\psi_d. We solve the resulting nonlinear equation in :math:\psi_d by bisection for each electrode potential.

Parameters:

Name Type Description Default
params GCSParameters

GCS parameters.

required
electrode_potentials_v NDArray[float64]

Grid of electrode potentials in volts, relative to the point of zero charge.

required
max_iter int

Maximum bisection iterations per potential.

80
tol float

Absolute tolerance on the residual E - (psi_H + psi_d) in volts.

1e-10

Returns:

Type Description
surface_charge

Array of diffuse-layer surface charge densities, C/m².

psi_diffuse

Array of diffuse-layer potentials, V.

differential_capacitance

Array of total differential capacitances at each electrode potential, F/m².

Source code in src/edl_ml/physics/gcs.py
def gouy_chapman_stern(
    params: GCSParameters,
    electrode_potentials_v: NDArray[np.float64],
    *,
    max_iter: int = 80,
    tol: float = 1e-10,
) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
    """Solve the self-consistent Gouy-Chapman-Stern problem.

    The electrode potential :math:`E` relative to the point of zero charge is
    split across the Stern and diffuse layers,

    .. math::

        E = \\psi_H + \\psi_d,\\qquad \\psi_H = \\sigma / C_H,

    with :math:`\\sigma` given by the Grahame equation as a function of
    :math:`\\psi_d`. We solve the resulting nonlinear equation in
    :math:`\\psi_d` by bisection for each electrode potential.

    Parameters
    ----------
    params
        GCS parameters.
    electrode_potentials_v
        Grid of electrode potentials in volts, relative to the point of zero
        charge.
    max_iter
        Maximum bisection iterations per potential.
    tol
        Absolute tolerance on the residual ``E - (psi_H + psi_d)`` in volts.

    Returns
    -------
    surface_charge
        Array of diffuse-layer surface charge densities, C/m².
    psi_diffuse
        Array of diffuse-layer potentials, V.
    differential_capacitance
        Array of total differential capacitances at each electrode potential,
        F/m².
    """
    c_h = stern_capacitance(params)
    thermal = BOLTZMANN * params.temperature_k
    n0 = params.concentration_mol_l * 1e3 * AVOGADRO
    eps = params.bulk_permittivity * VACUUM_PERMITTIVITY
    prefactor = np.sqrt(8.0 * eps * n0 * thermal)

    def sigma_of_psid(psi_d: float) -> float:
        """Grahame equation."""
        return float(
            np.sign(psi_d)
            * prefactor
            * np.sinh(params.valence * ELEMENTARY_CHARGE * abs(psi_d) / (2.0 * thermal))
        )

    def residual(psi_d: float, e: float) -> float:
        return e - psi_d - sigma_of_psid(psi_d) / c_h

    sigma = np.zeros_like(electrode_potentials_v)
    psi_d_arr = np.zeros_like(electrode_potentials_v)
    for idx, e in enumerate(electrode_potentials_v):
        lo, hi = -abs(e) - 1.0, abs(e) + 1.0
        f_lo = residual(lo, float(e))
        f_hi = residual(hi, float(e))
        if f_lo * f_hi > 0:
            # Expand bracket; the residual is monotone in psi_d so this
            # terminates quickly.
            for _ in range(20):
                lo *= 2.0
                hi *= 2.0
                f_lo = residual(lo, float(e))
                f_hi = residual(hi, float(e))
                if f_lo * f_hi <= 0:
                    break
        psi_d = 0.5 * (lo + hi)
        for _ in range(max_iter):
            psi_d = 0.5 * (lo + hi)
            f_mid = residual(psi_d, float(e))
            if abs(f_mid) < tol:
                break
            if f_lo * f_mid < 0:
                hi, f_hi = psi_d, f_mid
            else:
                lo, f_lo = psi_d, f_mid
        psi_d_arr[idx] = psi_d
        sigma[idx] = sigma_of_psid(psi_d)
    cap = total_capacitance(params, psi_d_arr)
    return sigma, psi_d_arr, np.asarray(cap, dtype=np.float64)

Constants

edl_ml.physics.constants

Physical constants used throughout the electric double layer model.

All constants are in SI units unless otherwise noted. Values follow CODATA 2018.

AVOGADRO module-attribute

AVOGADRO: Final[float] = 6.02214076e+23

Avogadro constant, 1/mol.

BOLTZMANN module-attribute

BOLTZMANN: Final[float] = 1.380649e-23

Boltzmann constant, J/K.

ELEMENTARY_CHARGE module-attribute

ELEMENTARY_CHARGE: Final[float] = 1.602176634e-19

Elementary charge, C.

FARADAY module-attribute

FARADAY: Final[float] = ELEMENTARY_CHARGE * AVOGADRO

Faraday constant, C/mol.

GAS_CONSTANT module-attribute

GAS_CONSTANT: Final[float] = BOLTZMANN * AVOGADRO

Ideal gas constant, J/(mol K).

ROOM_TEMPERATURE module-attribute

ROOM_TEMPERATURE: Final[float] = 298.15

Reference temperature, K.

VACUUM_PERMITTIVITY module-attribute

VACUUM_PERMITTIVITY: Final[float] = 8.8541878128e-12

Vacuum permittivity, F/m.

WATER_PERMITTIVITY module-attribute

WATER_PERMITTIVITY: Final[float] = 78.4

Relative permittivity of water at 298.15 K (dimensionless).