Simulating a polyanion–counterion–solvent electric double-layer capactor using the Gaussian core model with smeared electrostatics¶

This guide contains a complete example for simulating a salt-free polyelectrolyte supercapactor with perfectly conducting surfaces using the Gaussian core model with smeared electrostatics (GCMe) and evaluating the number and charge density profiles.

Parameters¶

The system is simulated at a temperature of \(300\,\mathrm{K}\) and contains a total of \(N=96,000\) particles at a number density of \(\rho=2.5d^{-3}\), where \(d\) is the particle size. With a polyanion ion fraction of \(x_\mathrm{p}\), there are \(x_\mathrm{p}N\) polyanion monomers, \(x_\mathrm{p}N\) counterions, and \((1-2x_\mathrm{p})N\) solvent particles. The polyanions have chain lengths of \(N_\mathrm{p}=60\) and have connectivity modeled by harmonic bonds with an equilibrium bond length of \(b=0.8d\) and a force constant of \(k=100k_\mathrm{B}T/d^2\), where \(k_\mathrm{B}\) is the Boltzmann constant.

Example Python script¶

Settings and variables¶

The following code cell sets up the logging format and defines constants, such as the fundamental quantities (length \(d\), mass \(m\), energy \(\epsilon\)), number of particles \(N\), number density \(\rho\), relative permittivity \(\varepsilon_\mathrm{r}\), etc., used in the simulation.

import logging

from openmm import unit

logging.basicConfig(
    format="{asctime} | {levelname:^8s} | {message}",
    style="{",
    level=logging.INFO
)

size: unit.Quantity = 0.275 * unit.nanometer
mass: unit.Quantity = 18.01528 * unit.amu
inverse_kappa_md: float = 15.9835
omega: float = 0.499

N: int = 96_000
temperature: unit.Quantity = 300.0 * unit.kelvin

N_m: int = 4
rho_md: float = 2.5
varepsilon_r: float = 78.0
N_p: int = 60
x_p: float = 0.2
b_md: float = 0.8
k_md: float = 100.0

Scale factors for coarse-grained simulations¶

To ensure consistency in scaling and units throughout the simulation, we first define the scale factors that will be used to convert between reduced simulation units and real units. This is accomplished by passing the fundamental quantities to the mdcraft.openmm.unit.get_scale_factors or mdcraft.openmm.unit.get_lj_scale_factors function. A dict containing the scale factors for molar energy \(N_\mathrm{A}\epsilon\), time \(\tau=\sqrt{md^2/\epsilon}\), force \(f=\epsilon/d\), pressure \(p=\epsilon/d^3\), etc. is returned.

from mdcraft.openmm.unit import get_lj_scale_factors

scales = get_lj_scale_factors(
    {
        "energy": (unit.BOLTZMANN_CONSTANT_kB * temperature)
                   .in_units_of(unit.kilojoule),
        "length": size * (N_m * rho_md) ** (1 / 3) if N_m > 1 else size,
        "mass": mass * N_m
    }
)
logging.info("Computed scaling factors for reducing physical quantities.\n"
             "  Fundamental quantities:\n"
             f"    Molar energy: {scales['molar_energy']}\n"
             f"    Length: {scales['length']}\n"
             f"    Mass: {scales['mass']}")

/home/docs/checkouts/readthedocs.org/user_builds/mdcraft/conda/stable/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
2025-01-15 22:46:29,781 |   INFO   | Computed scaling factors for reducing physical quantities.
  Fundamental quantities:
    Molar energy: 2.4943387854459718 kJ/mol
    Length: 0.5924695397587681 nm
    Mass: 72.06112 Da

Getting system dimensions¶

We use the number density \(\rho\) specified previously to determine the system dimensions.

We first get an estimate of the \(x\)- (\(L_x\)) and \(y\)-dimensions (\(L_y\)) by assuming that the \(z\)-dimension \(L_z\) will be approximately \(2.5\times\) the \(x\)- and \(y\)-dimensions:

\[\begin{equation*} L_x=L_y\approx [N/(c_L\rho)]^{1/3},\quad c_L=2.5 \end{equation*}\]

Then, we call the mdcraft.openmm.topology.create_atoms function to calculate system dimensions that will accomodate a infinitely repeating single-layer hexagonal close-packed (HCP) lattice in the \(xy\)-plane, which will serve as the surfaces in our simulation. The arguments required are, in order, the estimated \(L_x\) and \(L_y\) and \(L_z=0\) (to ensure that a single HCP lattice layer is generated), the type of lattice, the spacing between surface particles in the HCP lattice, and permission to adjust the estimated \(L_x\) and \(L_y\) to satisfy the system periodicity. Finally, we determine \(L_z\) using the number of particles \(N\), the number density \(\rho\), and the actual \(L_x\) and \(L_y\):

\[\begin{equation*} L_z=N/(\rho L_xL_y) \end{equation*}\]

from mdcraft.openmm.topology import create_atoms
import numpy as np

rho = rho_md / scales["length"] ** 3
L = ((N / (2.5 * rho)) ** (1 / 3))
positions_wall, dimensions = create_atoms(
    L * np.array((1, 1, 0)),
    lattice="hcp",
    length=scales["length"] / 2,
    flexible=True
)
dimensions[2] = N / (rho * dimensions[0] * dimensions[1])

Initializing the simulation system and topology¶

We create an OpenMM simulation system and topology using the system dimensions we previously calculated.

import openmm
from openmm import app

system = openmm.System()
system.setDefaultPeriodicBoxVectors(*(dimensions * np.diag(np.ones(3))))
topology = app.Topology()
topology.setUnitCellDimensions(dimensions)
logging.info("Created simulation system and topology with "
             f"dimensions {dimensions[0]} x {dimensions[1]} "
             f"x {dimensions[2]}.")

2025-01-15 22:46:29,792 |   INFO   | Created simulation system and topology with dimensions 14.811738493969203 nm x 14.879716499587463 nm x 36.23502900256258 nm.

Registering pair and bond potentials to the simulation system¶

We initialize and add the GCMe excluded volume and smeared electrostatic pair potentials and the harmonic bond potential to the simulation system.

Excluded volume pair potential¶

The excluded volume (Gaussian) potential has the form:

(1)¶\[\begin{equation} u_\mathrm{Gauss}(r_{ij})=\alpha_{ij}\exp{\left(-\beta_{ij}r_{ij}^2\right)},\quad\alpha_{ij}=A_{ij}\left(\frac{\beta_{ij}}{\pi}\right)^{3/2},\quad\beta_{ij}=\frac{3}{2\sigma_{ij}^2},\quad\sigma_{ij}^2=\sigma_i^2+\sigma_j^2, \end{equation}\]

where \(\sigma_{i}\) is the mass smearing radius of particle \(i\) and \(A_{ij}\) is the repulsion parameter between particles \(i\) and \(j\).

In the Gaussian potential, there are only two types of interactions that need to be accounted for: ion–ion and ion–surface. We use the same \(A_{ij}\) but different \(\sigma_{ij}\) for these interactions since the surface particles have no mass smearing radius (\(\sigma_{\mathrm{s},\,i}=0~\mathrm{nm}\)) such that the mixing rule will give an effective interaction range of \(\sigma_{ij}=\sigma_{i}\). This minimizes the surface thickness and accurately captures the image charge interactions.

The other possible interactions (ion–image charge, surface–surface, surface–image charge, and image charge–image charge) need to be disabled. This is done by setting \(\alpha_{ij}=0~\mathrm{kJ}\) for these pairs.

After evaluating \(\alpha_{ij}\) and \(\beta_{ij}\), we can determine a cutoff for the exponentially-decaying Gaussian potential by finding the \(r_{ij}\) value at which the strongest pair interaction energy drops below \(0.001\epsilon\).

Finally, we initialize the openmm.CustomNonbondedForce object for the Gaussian potential using the mdcraft.openmm.pair.gauss function by specifing tabulated functions containing the precalculated \(\alpha_{ij}\) and \(\beta_{ij}\).

from mdcraft.openmm.pair import gauss
from scipy import optimize

# excluded volume interaction potential
radius_nd = scales["length"].value_in_unit(unit.nanometer) / 2
sigmas_i_sq = (np.array((radius_nd, 0, radius_nd)) * unit.nanometer) ** 2
sigmas_ij_sq = sigmas_i_sq + sigmas_i_sq[:, None]
betas_ij = 3 / (2 * sigmas_ij_sq)
alphas_ij_coefs = 1 + np.array(
    (
        (0, 0, -1),     # pp, pw, pi;
        (0, -1, -1),    # wp, ww, wi;
        (-1, -1, -1)    # ip, iw, ii
    )
)
A_md = (N_m * inverse_kappa_md - 1) / (2 * omega * rho_md)
A = A_md * scales["molar_energy"] * scales["length"] ** 3
alphas_ij = alphas_ij_coefs * A * (betas_ij / np.pi) ** (3 / 2)
alphas_ij[np.isnan(alphas_ij)] = 0 * unit.kilojoule_per_mole
cutoff = optimize.fsolve(
    lambda r: np.max(alphas_ij)
              * np.exp(-np.min(betas_ij) * (r * unit.nanometer) ** 2)
              / scales["molar_energy"] - 1e-3,
    scales["length"].value_in_unit(unit.nanometer)
)[0] * unit.nanometer
pair_gauss = gauss(
    cutoff,
    mix="alpha12=alpha(type1,type2);beta12=beta(type1,type2);",
    per_params=("type",),
    tab_funcs={"alpha": alphas_ij, "beta": betas_ij}
)

/home/docs/checkouts/readthedocs.org/user_builds/mdcraft/conda/stable/lib/python3.12/site-packages/openmm/unit/quantity.py:443: RuntimeWarning: divide by zero encountered in reciprocal
  return Quantity(pow(self._value, exponent), pow(self.unit, exponent))
/home/docs/checkouts/readthedocs.org/user_builds/mdcraft/conda/stable/lib/python3.12/site-packages/openmm/unit/quantity.py:683: RuntimeWarning: invalid value encountered in multiply
  value = factor * self._value # works for number, numpy.array, or vec3, e.g.

Smeared electrostatic pair potential¶

The smeared electrostatic potential has the form:

\[\begin{equation*} u_\mathrm{elec}(r_{ij})=\frac{q_iq_j}{4\pi\varepsilon_0r_{ij}}\mathrm{erf}(\alpha_{ij}r_{ij}),\quad\alpha_{ij}=\sqrt{\frac{\pi}{2a_{ij}^2}},\quad a_{ij}^2=a_i^2+a_j^2, \end{equation*}\]

where \(\varepsilon_0\) is the vacuum permittivity and \(a_i\) is the charge smearing radius of particle \(i\).

Here, there are only two unique interaction types: charge–charge and charge–surface. “Charge” can refer to either an ion or an image charge. Like before, the surface particles has no charge smearing radius (\(a_{\mathrm{s},\,i}=0~\mathrm{nm}\)) so that potential differences, if necessary, can correctly be applied by giving surface particles charges.

Finally, we initialize the openmm.NonbondedForce and openmm.CustomNonbondedForce objects for the smeared electrostatic potential using the mdcraft.openmm.pair.coul_gauss function by specifing tabulated functions containing the precalculated \(a_{ij}\). The openmm.CustomNonbondedForce evaluates the electrostatic interactions before the cutoff in real space, and the openmm.NonbondedForce evaluates the electrostatic interactions after the cutoff in reciprocal space.

from mdcraft.openmm.pair import coul_gauss

# smeared electrostatic interaction potential
as_i_sq = (np.array((1, 0)) * scales["length"] / 2) ** 2
as_ij = (as_i_sq + as_i_sq[:, None]) ** (1 / 2) # pp, pw; wp, ww
pair_elec_real, pair_elec_fourier = coul_gauss(
    cutoff,
    mix="alpha12=alpha(type1,type2);",
    per_params=("type",),
    tab_funcs={"alpha": np.sqrt(np.pi / 2) / as_ij}
)

Harmonic bond potential¶

The harmonic bond potential has the form:

\[\begin{equation*} u_\mathrm{harm}(r_{ij})=\frac{1}{2}k_{ij}(r_{ij}-b_{ij})^2, \end{equation*}\]

where \(k_{ij}\) is the force constant and \(b_{ij}\) is the equilibrium bond length.

For now, we just initialize the openmm.HarmonicBondForce object for the harmonic bond potential.

# harmonic bond potential
if x_p > 0:
    b = b_md * scales["length"]
    k = k_md * scales["molar_energy"] / scales["length"] ** 2
    bond_harm = openmm.HarmonicBondForce()

system.addForce(pair_gauss)
system.addForce(pair_elec_real)
system.addForce(pair_elec_fourier)
if x_p > 0:
    system.addForce(bond_harm)
logging.info(f"Registered {system.getNumForces()} pair "
             "potential(s) to the simulation.")

2025-01-15 22:46:29,808 |   INFO   | Registered 4 pair potential(s) to the simulation.

Determining the identity and number of the different particles¶

For easy reference later, we first determine the number of particles for each particle type. We also assign arbitrary but logical particle identities to the polyanion monomers, counterions, and solvent particles so that visualization tools like OVITO can distinguish between them.

# Assign arbitrary particle identities
element_a = app.Element.getBySymbol("Cl")
element_c = app.Element.getBySymbol("Na")
element_s = app.Element.getBySymbol("Ar")
element_w = app.Element.getBySymbol("C")

# Determine the number of polyanions, counterions, and solvent particles
M = round(x_p * N / N_p)    # Number of polyanions
N_a = N_c = M * N_p         # Number of polyanion beads and/or counterions
N_s = N - N_a - N_c         # Number of solvent particles
if N_a != x_p * N:
    emsg = (f"The polyanion chain length {N_p=} is incompatible "
            f"with the total number of particles {N=} and the "
            f"polyanion number concentration {x_p=}.")
    raise RuntimeError(emsg)

Registering particles to the simulation system and force field¶

Now, we register the particles to the simulation system and the potentials we created previously. …

Note that the smeared electrostatic potential does not have a relative permittivity \(\varepsilon_\mathrm{r}\) term in the denominator. As such, we need to scale all charges by \(1/\sqrt{\varepsilon_\mathrm{r}}\) to get the correct electrostatic behavior.

from itertools import combinations

from mdcraft.openmm.system import register_particles

q_scaled = unit.elementary_charge / np.sqrt(varepsilon_r)

# Register polyanions to pair potentials
for _ in range(M):
    chain = topology.addChain()
    register_particles(
        system, topology, N_p, scales["mass"],
        chain=chain,
        element=element_a,
        name="PAN",
        nbforce=pair_elec_fourier,
        charge=-q_scaled,
        cnbforces={pair_elec_real: (-q_scaled, 0), pair_gauss: (0,)}
    )
logging.info(f"Registered {M:,} polyanion(s) with {N_p:,} monomer(s) "
                "to the force field.")

# Register polyanion bonds to bond potential and remove 1-2 interactions
if x_p > 0:
    atoms = list(topology.atoms())
    for m in range(M):
        for n in range(N_p - 1):
            i = m * N_p + n
            j = i + 1
            topology.addBond(atoms[i], atoms[j])
            bond_harm.addBond(i, j, b, k)
            pair_elec_real.addExclusion(i, j)
            pair_elec_fourier.addException(i, j, 0, 0, 0)
            pair_gauss.addExclusion(i, j)
    logging.info(f"Registered {topology.getNumBonds():,} bond(s) to "
                    "the force field.")

# Register counterions to pair potentials
register_particles(
    system, topology, N_c, scales["mass"],
    element=element_c,
    name="CAT",
    nbforce=pair_elec_fourier,
    charge=q_scaled,
    cnbforces={pair_elec_real: (q_scaled, 0), pair_gauss: (0,)}
)
logging.info(f"Registered {N_c:,} counterion(s) to the simulation.")

# Register solvent particles to pair potentials
register_particles(
    system, topology, N_s, scales["mass"],
    element=element_s,
    name="SOL",
    resname="SOL",
    nbforce=pair_elec_fourier,
    cnbforces={pair_elec_real: (0, 0), pair_gauss: (0,)}
)
logging.info(f"Registered {N_s:,} solvent particle(s) to the simulation.")

# Determine positions and number of wall particles
positions_wall = np.concatenate((
    positions_wall,
    positions_wall + np.array(
        (0, 0, dimensions[2].value_in_unit(unit.nanometer))
    ) * unit.nanometer
))
N_wall = positions_wall.shape[0]

# Register wall particles to pair potentials
for name in ("LWL", "RWL"):
    register_particles(
        system, topology, N_wall // 2, 0,
        element=element_w,
        name=name,
        nbforce=pair_elec_fourier,
        cnbforces={pair_elec_real: (0, 1), pair_gauss: (1,)}
    )
logging.info(f"Registered {N_wall:,} wall particles to the force field.")

# Remove wall–wall interactions
wall_indices = range(N, N + N_wall)
for i, j in combinations(wall_indices, 2):
    pair_elec_real.addExclusion(i, j)
    pair_elec_fourier.addException(i, j, 0, 0, 0)
    pair_gauss.addExclusion(i, j)
logging.info("Removed wall–wall interactions.")

2025-01-15 22:46:30,216 |   INFO   | Registered 320 polyanion(s) with 60 monomer(s) to the force field.
2025-01-15 22:46:30,390 |   INFO   | Registered 18,880 bond(s) to the force field.
2025-01-15 22:46:30,805 |   INFO   | Registered 19,200 counterion(s) to the simulation.
2025-01-15 22:46:31,340 |   INFO   | Registered 57,600 solvent particle(s) to the simulation.
2025-01-15 22:46:31,362 |   INFO   | Registered 5,800 wall particles to the force field.
2025-01-15 22:46:59,779 |   INFO   | Removed wall–wall interactions.

The rest of this guide is to be finished at a later date.