A Charged System: Counterion Condensation

Table of contents

• Introduction
• System setup
• First run and observable setup
• Production run and analysis
• Overcharging by added salt

Introduction

In this tutorial, we simulate a charged system consisting of a fixed charged rod with ions around it. This setup represents a simplified model for polyelectrolyte gels. We will investigate the condensation of ions onto the oppositely charged rod and compare the results to a meanfield analytical solution obtained from Poisson−Boltzmann (PB) theory. Finally we will go beyond the expected applicability of PB and add concentrated additional salt ions to observe an overcharging effect.

The tutorial follows "Deserno, Markus, Christian Holm, and Sylvio May. "Fraction of condensed counterions around a charged rod: Comparison of Poisson−Boltzmann theory and computer simulations. Macromolecules 33.1 (2000): 199-206, 10.1021/ma990897o". We refer to that publication for further reading.

import espressomd
import espressomd.electrostatics
import espressomd.observables
import espressomd.accumulators
import espressomd.math

espressomd.assert_features(['ELECTROSTATICS', 'P3M', 'WCA'])

import tqdm
import numpy as np
import scipy.optimize
%matplotlib inline
import matplotlib.pyplot as plt

np.random.seed(41)
plt.rcParams.update({'font.size': 18})

System setup

After importing the necessary ESPResSo features and external modules, we define a cubic system geometry and some physical parameters (which define our unit system).

# system parameters
ROD_LENGTH = 50
BJERRUM_LENGTH = 1.0

# we assume a unit system where the elementary charge and the thermal energy are both 1
system = espressomd.System(box_l=3 * [ROD_LENGTH])
KT = 1.
Q_E = 1.

system.time_step = 0.01
system.cell_system.skin = 0.4

We will build the charged rod from individual particles that are fixed in space. With this, we can use the particle-based electrostatics methods of ESPResSo. For analysis, we give the rod particles a different type than the counterions.

# interaction parameters
WCA_EPSILON = 1.0
ION_DIAMETER = 1.0
ROD_RADIUS = 1.0
# particle types
ROD_TYPE = 1
COUNTERION_TYPE = 2

Exercise:

• Setup the purely repulsive Weeks-Chandler-Anderson (WCA) interaction (Non-bonded Interactions) between the ions and between the ions and the rod particles. Use the parameters introduced in the cell above.

Hints:

• The WCA potential uses the same parameters as the Lennard-Jones potential, but the cutoff and shift are calculated automatically
• Use the Lorentz combining rule (arithmetic mean) to determine the sigma parameter of the interaction between the rod particles and the ions

# ion-ion interaction
system.non_bonded_inter[COUNTERION_TYPE, COUNTERION_TYPE].wca.set_params(
    epsilon=WCA_EPSILON, sigma=ION_DIAMETER)

# ion-rod interaction
system.non_bonded_inter[COUNTERION_TYPE, ROD_TYPE].wca.set_params(
    epsilon=WCA_EPSILON, sigma=ION_DIAMETER / 2. + ROD_RADIUS)

Now we need to place the particles in the box

Exercise:

• Implement a function to place the rod particles along the $x_3$ axis in the middle of the simulation box and the ions randomly distributed
• Use the signature setup_rod_and_counterions(system, ion_valency, counterion_type, rod_charge_dens, N_rod_beads, rod_type)
• Determine the number of counterions from the condition of neutrality for the whole system (the rod should be positive, the counterions negative)
• Assign the rod particles and counterions their correct type
• Give the counterions a charge q according to their ion_valency
• Give the rod particles a charge such that the rod_charge_dens is uniformly distributed along the N_rod_beads individual particles
• Fix the rod particles in space so they do not get moved if forces act upon them
• Return the newly created counterion particles

Hints:

• Look into espresso particle properties to find the keywords to set charges and to fix particles
• use np.random.random() to generate the counterion positions

def setup_rod_and_counterions(system, ion_valency, counterion_type,
                              rod_charge_dens, N_rod_beads, rod_type):

    # calculate charge of the single rod beads
    rod_length = system.box_l[2]
    total_rod_charge = rod_charge_dens * rod_length
    rod_charge_per_bead = total_rod_charge / N_rod_beads

    # number of counterions
    N_ions = int(total_rod_charge / ion_valency)

    rod_zs = np.linspace(0, rod_length, num=N_rod_beads, endpoint=False)
    rod_positions = np.column_stack(([system.box_l[0] / 2.] * N_rod_beads,
                                     [system.box_l[1] / 2.] * N_rod_beads,
                                     rod_zs))

    system.part.add(pos=rod_positions, type=[rod_type] * N_rod_beads,
                    q=[rod_charge_per_bead] * N_rod_beads,
                    fix=[3 * [True]] * N_rod_beads)

    ion_positions = np.random.random((N_ions, 3)) * system.box_l

    counter_ions = system.part.add(pos=ion_positions, type=[
                                   counterion_type] * N_ions, q=[-ion_valency] * N_ions)

    return counter_ions

COUNTERION_VALENCY = 1
ROD_CHARGE_DENS = 2

# number of beads that make up the rod
N_rod_beads = int(ROD_LENGTH / ROD_RADIUS)

setup_rod_and_counterions(system, COUNTERION_VALENCY, COUNTERION_TYPE,
                          ROD_CHARGE_DENS, N_rod_beads, ROD_TYPE)

# check that the particle setup was done correctly
assert abs(sum(system.part.all().q)) < 1e-10
assert np.all(system.part.select(type=ROD_TYPE).fix)

Now we set up the electrostatics method to calculate the forces and energies from the long-range Coulomb interaction. ESPResSo uses so-called actors for electrostatics, magnetostatics and hydrodynamics. This ensures that unphysical combinations of algorithms are avoided, for example simultaneous usage of two electrostatic interactions. Adding an actor to the system also activates the method and calls necessary initialization routines. Here, we define a P$^3$M object using the Bjerrum length and rms force error. This automatically starts a tuning function which tries to find optimal parameters for P$^3$M and prints them to the screen. For more details, see the user guide chapter on Electrostatics.

p3m_params = {'prefactor': KT * BJERRUM_LENGTH * Q_E**2,
              'accuracy': 1e-3}

For the accuracy, ESPResSo estimates the relative error in the force calculation introduced by the approximations of $P^3M$. We choose a relatively poor accuracy (large value) for this tutorial to make it run faster. For your own production simulations you should reduce that number.

Exercise:

• Set up a p3m instance and add it to the actors of the system

p3m = espressomd.electrostatics.P3M(**p3m_params)
system.actors.add(p3m)

WARNING: Statistics of tuning samples is very bad.

CoulombP3M tune parameters: Accuracy goal = 1.00000e-03 prefactor = 1.00000e+00
System: box_l = 5.00000e+01 # charged part = 150 Sum[q_i^2] = 3.00000e+02
mesh cao r_cut_iL    alpha_L     err       rs_err    ks_err    time [ms]
6    7   cao too large for this mesh
10   7   2.37352e-01 8.46940e+00 9.881e-04 7.071e-04 6.902e-04 0.27    
10   6   2.42156e-01 8.29105e+00 9.864e-04 7.071e-04 6.878e-04 0.24    
10   5   2.50805e-01 7.98771e+00 9.948e-04 7.071e-04 6.997e-04 0.23    
10   4   2.68102e-01 7.44127e+00 1.000e-03 7.071e-04 7.071e-04 0.23    
10   3   3.10383e-01 6.36820e+00 9.996e-04 7.071e-04 7.065e-04 0.26    
10   2   4.92000e-01 3.89721e+00 1.029e-03 7.071e-04 7.471e-04 accuracy not achieved
16   7   1.56753e-01 1.31492e+01 9.997e-04 7.071e-04 7.066e-04 0.36    
16   6   1.61651e-01 1.27276e+01 9.817e-04 7.071e-04 6.810e-04 0.32    
16   5   1.67530e-01 1.22551e+01 9.950e-04 7.071e-04 7.000e-04 0.30    
16   4   1.79286e-01 1.14053e+01 9.974e-04 7.071e-04 7.035e-04 0.29    
16   3   2.10637e-01 9.61380e+00 9.897e-04 7.071e-04 6.924e-04 0.30    
16   2   2.50805e-01 7.98771e+00 1.859e-03 7.071e-04 1.720e-03 accuracy not achieved
20   7   1.29562e-01 1.60861e+01 9.845e-04 7.071e-04 6.850e-04 0.56    
20   6   1.33064e-01 1.56387e+01 9.829e-04 7.071e-04 6.827e-04 0.50    
20   5   1.37966e-01 1.50515e+01 9.949e-04 7.071e-04 6.999e-04 0.48    
20   4   1.47771e-01 1.39968e+01 9.990e-04 7.071e-04 7.057e-04 0.46    
20   3   1.74384e-01 1.17454e+01 9.956e-04 7.071e-04 7.008e-04 0.46    
20   2   1.79286e-01 1.14053e+01 2.850e-03 7.071e-04 2.761e-03 accuracy not achieved
26   7   1.02747e-01 2.05532e+01 9.915e-04 7.071e-04 6.950e-04 1.29    
26   6   1.05633e-01 1.99606e+01 9.860e-04 7.071e-04 6.871e-04 1.25    
26   5   1.09674e-01 1.91846e+01 9.947e-04 7.071e-04 6.996e-04 1.22    
26   4   1.17755e-01 1.77961e+01 9.980e-04 7.071e-04 7.043e-04 1.19    
26   3   1.40267e-01 1.47904e+01 9.924e-04 7.071e-04 6.964e-04 1.18    
26   2   1.47771e-01 1.39968e+01 2.778e-03 7.071e-04 2.686e-03 accuracy not achieved
30   7   9.04065e-02 2.35257e+01 9.990e-04 7.071e-04 7.056e-04 1.82    
30   6   9.31461e-02 2.27961e+01 9.852e-04 7.071e-04 6.860e-04 1.77    
30   5   9.69815e-02 2.18455e+01 9.852e-04 7.071e-04 6.860e-04 1.77    
30   4   1.04104e-01 2.02703e+01 9.948e-04 7.071e-04 6.997e-04 1.71    
30   3   1.24377e-01 1.67960e+01 9.940e-04 7.071e-04 6.986e-04 1.67    
30   2   1.40267e-01 1.47904e+01 2.380e-03 7.071e-04 2.272e-03 accuracy not achieved
36   7   7.77359e-02 2.75877e+01 9.605e-04 7.071e-04 6.500e-04 3.42    
36   6   7.96793e-02 2.68787e+01 9.672e-04 7.071e-04 6.600e-04 3.17    
36   5   8.25944e-02 2.58792e+01 9.857e-04 7.071e-04 6.867e-04 3.13    
36   4   8.93963e-02 2.38064e+01 9.757e-04 7.071e-04 6.723e-04 3.01    
36   3   1.06887e-01 1.97134e+01 9.923e-04 7.071e-04 6.962e-04 2.95    
36   2   1.24377e-01 1.67960e+01 2.259e-03 7.071e-04 2.146e-03 accuracy not achieved

resulting parameters: mesh: (10, 10, 10), cao: 5, r_cut_iL: 2.5080e-01,
                      alpha_L: 7.9877e+00, accuracy: 9.9479e-04, time: 0.23

Before we can start the simulation, we need to remove the overlap between particles to avoid large forces which would crash the simulation. For this, we use the steepest descent integrator with a relative convergence criterion for forces and energies.

def remove_overlap(system, sd_params):
    # Removes overlap by steepest descent until forces or energies converge
    # Set up steepest descent integration
    system.integrator.set_steepest_descent(f_max=0,
                                           gamma=sd_params['damping'],
                                           max_displacement=sd_params['max_displacement'])

    # Initialize integrator to obtain initial forces
    system.integrator.run(0)
    maxforce = np.max(np.linalg.norm(system.part.all().f, axis=1))
    energy = system.analysis.energy()['total']

    i = 0
    while i < sd_params['max_steps'] // sd_params['emstep']:
        prev_maxforce = maxforce
        prev_energy = energy
        system.integrator.run(sd_params['emstep'])
        maxforce = np.max(np.linalg.norm(system.part.all().f, axis=1))
        relforce = np.abs((maxforce - prev_maxforce) / prev_maxforce)
        energy = system.analysis.energy()['total']
        relener = np.abs((energy - prev_energy) / prev_energy)
        if i > 1 and (i + 1) % 4 == 0:
            print(f"minimization step: {(i+1)*sd_params['emstep']:4.0f}"
                  f"    max. rel. force change:{relforce:+3.3e}"
                  f"    rel. energy change:{relener:+3.3e}")
        if relforce < sd_params['f_tol'] or relener < sd_params['e_tol']:
            break
        i += 1

    system.integrator.set_vv()

STEEPEST_DESCENT_PARAMS = {'f_tol': 1e-2,
                           'e_tol': 1e-5,
                           'damping': 30,
                           'max_steps': 10000,
                           'max_displacement': 0.01,
                           'emstep': 10}

remove_overlap(system, STEEPEST_DESCENT_PARAMS)

minimization step:   40    max. rel. force change:+4.122e-02    rel. energy change:+5.678e-03
minimization step:   80    max. rel. force change:+4.867e-02    rel. energy change:+5.896e-03
minimization step:  120    max. rel. force change:+1.576e-01    rel. energy change:+6.261e-03
minimization step:  160    max. rel. force change:+1.292e-01    rel. energy change:+6.096e-03
minimization step:  200    max. rel. force change:+3.474e-02    rel. energy change:+5.932e-03
minimization step:  240    max. rel. force change:+4.081e-02    rel. energy change:+6.393e-03
minimization step:  280    max. rel. force change:+4.900e-02    rel. energy change:+6.947e-03
minimization step:  320    max. rel. force change:+6.124e-02    rel. energy change:+7.692e-03
minimization step:  360    max. rel. force change:+6.884e-02    rel. energy change:+7.363e-03
minimization step:  400    max. rel. force change:+6.468e-02    rel. energy change:+7.311e-03
minimization step:  440    max. rel. force change:+8.435e-02    rel. energy change:+8.266e-03
minimization step:  480    max. rel. force change:+4.175e-02    rel. energy change:+7.374e-03
minimization step:  520    max. rel. force change:+4.906e-02    rel. energy change:+8.191e-03
minimization step:  560    max. rel. force change:+6.033e-02    rel. energy change:+9.252e-03
minimization step:  600    max. rel. force change:+5.565e-02    rel. energy change:+9.827e-03

After the overlap is removed, we activate a thermostat to simulate the system at a given temperature.

LANGEVIN_PARAMS = {'kT': KT,
                   'gamma': 0.5,
                   'seed': 42}
system.thermostat.set_langevin(**LANGEVIN_PARAMS)

First run and observable setup

Before running the simulations to obtain the histograms, we need to decide how long we need to equilibrate the system. For this we plot the total energy vs the time steps.

energies = []
STEPS_PER_SAMPLE_FIRST_RUN = 10
N_SAMPLES_FIRST_RUN = 1000
for i in range(N_SAMPLES_FIRST_RUN):
    system.integrator.run(STEPS_PER_SAMPLE_FIRST_RUN)
    energies.append(system.analysis.energy()['total'])

# plot time in time_steps so we can judge the number of warmup steps
ts = np.arange(0, N_SAMPLES_FIRST_RUN) * STEPS_PER_SAMPLE_FIRST_RUN
plt.figure(figsize=(10, 7))
plt.plot(ts, energies)
plt.xlabel('time steps')
plt.ylabel('system total energy')
plt.show()

WARMUP_STEPS = 5000
STEPS_PER_SAMPLE = 100

Now we are ready to implement the observable calculation. As we are interested in the condensation of counterions on the rod, the physical quantity of interest is the density of charges $\rho(r)$ around the rod, where $r$ is the distance from the rod. We need many samples to calculate the density from histograms.

From the last tutorial you should already be familiar with the concepts of observables and accumulators in ESPResSo. We will use the CylindricalDensityProfile observable and the MeanVarianceCalculator accumulator

Exercise:

• Write a function setup_profile_calculation(system, delta_N, ion_types, r_min, n_radial_bins) to create observables for $\rho(r)$
• delta_N is the number of integration steps between observable calculation
• ion_types is a list of types for which the radial distances should be calculated. For the moment we only have counterions, but later we will also add additional salt ions for which we would also like to calculate the density
• return a a dictionary of the accumulators radial_distances[counterion_type] = <accumulator> and the edges of the bins

Hints:

• Use system.part.select(type=...) to get only the particles of a specific type
• The azimuthal angle and the $x_3$ position are irrelevant, so you need only one big bin for these coordinates

def setup_profile_calculation(system, delta_N, ion_types, r_min, n_radial_bins):
    radial_profile_accumulators = {}
    ctp = espressomd.math.CylindricalTransformationParameters(center = np.array(system.box_l) / 2.,
                                                              axis = [0, 0, 1],
                                                              orientation = [1, 0, 0])
    for ion_type in ion_types:
        ion_ids = system.part.select(type=ion_type).id
        radial_profile_obs = espressomd.observables.CylindricalDensityProfile(
            ids=ion_ids,
            transform_params = ctp,
            n_r_bins=n_radial_bins,
            min_r=r_min,
            min_z=-system.box_l[2] / 2.,
            max_r=system.box_l[0] / 2.,
            max_z=system.box_l[2] / 2.)

        bin_edges = radial_profile_obs.bin_edges()

        radial_profile_acc = espressomd.accumulators.MeanVarianceCalculator(
            obs=radial_profile_obs, delta_N=delta_N)
        system.auto_update_accumulators.add(radial_profile_acc)

        radial_profile_accumulators[ion_type] = radial_profile_acc

    return radial_profile_accumulators, bin_edges

r_min = ROD_RADIUS + ION_DIAMETER / 2.
r_max = system.box_l[0] / 2.
N_RADIAL_BINS = 200
radial_profile_accs, bin_edges = setup_profile_calculation(
    system, STEPS_PER_SAMPLE, [COUNTERION_TYPE], r_min, N_RADIAL_BINS)
assert isinstance(
    radial_profile_accs[COUNTERION_TYPE], espressomd.accumulators.MeanVarianceCalculator)
assert len(bin_edges) == N_RADIAL_BINS + 1

To run the simulation with different parameters, we need a way to reset the system and return it to an empty state before setting it up again.

Exercise:

• Write a function clear_system(system) that
  • turns off the thermostat
  • removes all particles
  • removes all actors
  • removes all accumulators added to the auto-update-list
  • resets the system clock

Hints:

• The relevant parts of the documentation can be found here: Thermostats, ParticleList, Electrostatics, AutoUpdateAccumulators, System properties

def clear_system(system):
    system.thermostat.turn_off()
    system.part.clear()
    system.actors.clear()
    system.auto_update_accumulators.clear()
    system.time = 0.

clear_system(system)

Production run and analysis

Now we are finally ready to run the simulations and produce the data we can compare to the Poisson-Boltzmann predictions. First we define the parameters and then loop over them.

runs = [{'params': {'counterion_valency': 2, 'rod_charge_dens': 1},
         'histogram': None},
        {'params': {'counterion_valency': 1, 'rod_charge_dens': 2},
         'histogram': None}
        ]
N_SAMPLES = 1500

For longer simulation runs it will be convenient to have a progress bar

def integrate_system(system, n_steps):
    for i in tqdm.trange(100):
        system.integrator.run(n_steps // 100)

Exercise:

• Run the simulation for the parameters given above and save the histograms in the corresponding dictionary for analysis

Hints:

• Don't forget to clear the system before setting up the system with a new set of parameters
• Don't forget to add a new p3m instance after each change of parameters. If we reuse the p3m that was tuned before, likely the desired accuracy will not be achieved.
• Extract the radial density profile from the accumulator via .mean()

for run in runs:
    clear_system(system)
    setup_rod_and_counterions(
        system, run['params']['counterion_valency'], COUNTERION_TYPE,
        run['params']['rod_charge_dens'], N_rod_beads, ROD_TYPE)
    p3m = espressomd.electrostatics.P3M(**p3m_params)
    system.actors.add(p3m)
    remove_overlap(system, STEEPEST_DESCENT_PARAMS)
    system.thermostat.set_langevin(**LANGEVIN_PARAMS)
    print('', end='', flush=True)
    integrate_system(system, WARMUP_STEPS)
    radial_profile_accs, bin_edges = setup_profile_calculation(
        system, STEPS_PER_SAMPLE, [COUNTERION_TYPE], r_min, N_RADIAL_BINS)
    integrate_system(system, N_SAMPLES * STEPS_PER_SAMPLE)

    run['histogram'] = radial_profile_accs[COUNTERION_TYPE].mean()
    print(f'simulation for parameters {run["params"]} done\n')

CoulombP3M tune parameters: Accuracy goal = 1.00000e-03 prefactor = 1.00000e+00
System: box_l = 5.00000e+01 # charged part = 75 Sum[q_i^2] = 1.50000e+02
mesh cao r_cut_iL    alpha_L     err       rs_err    ks_err    time [ms]
4    7   cao too large for this mesh
10   7   2.17172e-01 8.90360e+00 9.870e-04 7.071e-04 6.886e-04 0.17    
10   6   2.21016e-01 8.73849e+00 9.866e-04 7.071e-04 6.881e-04 0.15    
10   5   2.28703e-01 8.42538e+00 9.868e-04 7.071e-04 6.883e-04 0.14    
10   4   2.43117e-01 7.89317e+00 9.929e-04 7.071e-04 6.970e-04 0.13    
10   3   2.76750e-01 6.87266e+00 9.955e-04 7.071e-04 7.007e-04 0.14    
10   2   4.20891e-01 4.38613e+00 9.980e-04 7.071e-04 7.043e-04 0.21    
10   1   4.92000e-01 3.70898e+00 5.274e-03 7.071e-04 5.226e-03 accuracy not achieved
14   7   1.61445e-01 1.22121e+01 9.942e-04 7.071e-04 6.989e-04 0.24    
14   6   1.65244e-01 1.19135e+01 9.958e-04 7.071e-04 7.012e-04 0.22    
14   5   1.71891e-01 1.14236e+01 9.888e-04 7.071e-04 6.912e-04 0.21    
14   4   1.83288e-01 1.06686e+01 9.853e-04 7.071e-04 6.862e-04 0.20    
14   3   2.08929e-01 9.27880e+00 9.991e-04 7.071e-04 7.059e-04 0.20    
14   2   2.43117e-01 7.89317e+00 1.736e-03 7.071e-04 1.585e-03 accuracy not achieved
20   7   1.19155e-01 1.68665e+01 9.781e-04 7.071e-04 6.758e-04 0.47    
20   6   1.21603e-01 1.65061e+01 9.955e-04 7.071e-04 7.007e-04 0.46    
20   5   1.26500e-01 1.58282e+01 9.889e-04 7.071e-04 6.913e-04 0.46    
20   4   1.34661e-01 1.48109e+01 9.938e-04 7.071e-04 6.983e-04 0.43    
20   3   1.55880e-01 1.26765e+01 9.895e-04 7.071e-04 6.922e-04 0.43    
20   2   2.08929e-01 9.27880e+00 1.320e-03 7.071e-04 1.115e-03 accuracy not achieved
24   7   1.01688e-01 1.99567e+01 9.790e-04 7.071e-04 6.770e-04 0.88    
24   6   1.04123e-01 1.94618e+01 9.853e-04 7.071e-04 6.862e-04 0.86    
24   5   1.07777e-01 1.87626e+01 9.963e-04 7.071e-04 7.019e-04 0.86    
24   4   1.15084e-01 1.75007e+01 9.945e-04 7.071e-04 6.993e-04 0.83    
24   3   1.33351e-01 1.49656e+01 9.990e-04 7.071e-04 7.056e-04 0.82    
24   2   1.55880e-01 1.26765e+01 1.871e-03 7.071e-04 1.733e-03 accuracy not achieved
30   7   8.33443e-02 2.46412e+01 9.943e-04 7.071e-04 6.991e-04 1.75    
30   6   8.54279e-02 2.40050e+01 9.973e-04 7.071e-04 7.033e-04 1.72    
30   5   8.90742e-02 2.29651e+01 9.841e-04 7.071e-04 6.844e-04 1.71    
30   4   9.48042e-02 2.14966e+01 9.987e-04 7.071e-04 7.053e-04 1.65    
30   3   1.10952e-01 1.81932e+01 9.937e-04 7.071e-04 6.981e-04 1.65    
30   2   1.33351e-01 1.49656e+01 1.803e-03 7.071e-04 1.659e-03 accuracy not achieved
34   7   7.54128e-02 2.73938e+01 9.606e-04 7.071e-04 6.502e-04 5.02    
34   6   7.71464e-02 2.67426e+01 9.706e-04 7.071e-04 6.649e-04 5.00    
34   5   7.97468e-02 2.58201e+01 9.888e-04 7.071e-04 6.911e-04 4.98    
34   4   8.58145e-02 2.38904e+01 9.726e-04 7.071e-04 6.679e-04 4.96    
34   3   9.96835e-02 2.03825e+01 9.997e-04 7.071e-04 7.067e-04 4.93    
34   2   1.10952e-01 1.81932e+01 2.250e-03 7.071e-04 2.136e-03 accuracy not achieved

resulting parameters: mesh: (10, 10, 10), cao: 4, r_cut_iL: 2.4312e-01,
                      alpha_L: 7.8932e+00, accuracy: 9.9287e-04, time: 0.13
minimization step:   40    max. rel. force change:+2.980e-02    rel. energy change:+6.903e-03
minimization step:   80    max. rel. force change:+3.472e-02    rel. energy change:+6.025e-03
minimization step:  120    max. rel. force change:+4.095e-02    rel. energy change:+5.899e-03
minimization step:  160    max. rel. force change:+4.868e-02    rel. energy change:+6.142e-03
minimization step:  200    max. rel. force change:+6.053e-02    rel. energy change:+6.782e-03

WARNING: Statistics of tuning samples is very bad.
100%|███████████████████████████████████████████████████████████████████| 100/100 [00:00<00:00, 128.00it/s]
100%|████████████████████████████████████████████████████████████████████| 100/100 [00:25<00:00,  3.99it/s]

simulation for parameters {'counterion_valency': 2, 'rod_charge_dens': 1} done

minimization step:   40    max. rel. force change:+7.683e-01    rel. energy change:+4.830e-01
minimization step:   80    max. rel. force change:+9.299e-01    rel. energy change:+7.917e-03
minimization step:  120    max. rel. force change:+5.163e-02    rel. energy change:+5.519e-03
minimization step:  160    max. rel. force change:+7.904e-02    rel. energy change:+5.604e-03
minimization step:  200    max. rel. force change:+5.227e-01    rel. energy change:+5.212e-03
minimization step:  240    max. rel. force change:+5.482e-02    rel. energy change:+5.425e-03
minimization step:  280    max. rel. force change:+1.176e-01    rel. energy change:+5.363e-03
minimization step:  320    max. rel. force change:+3.314e-02    rel. energy change:+5.596e-03
minimization step:  360    max. rel. force change:+3.778e-02    rel. energy change:+6.035e-03
minimization step:  400    max. rel. force change:+4.444e-02    rel. energy change:+6.598e-03
minimization step:  440    max. rel. force change:+5.373e-02    rel. energy change:+7.315e-03
minimization step:  480    max. rel. force change:+1.656e-02    rel. energy change:+7.980e-03
minimization step:  520    max. rel. force change:+8.291e-02    rel. energy change:+7.887e-03
minimization step:  560    max. rel. force change:+1.542e-01    rel. energy change:+7.834e-03
minimization step:  600    max. rel. force change:+5.499e-02    rel. energy change:+8.175e-03
minimization step:  640    max. rel. force change:+4.641e-02    rel. energy change:+8.432e-03
minimization step:  680    max. rel. force change:+5.584e-02    rel. energy change:+9.530e-03
minimization step:  720    max. rel. force change:+7.505e-02    rel. energy change:+9.429e-03
minimization step:  760    max. rel. force change:+1.952e-01    rel. energy change:+9.324e-03
minimization step:  800    max. rel. force change:+5.948e-02    rel. energy change:+9.910e-03
minimization step:  840    max. rel. force change:+5.323e-02    rel. energy change:+1.097e-02
minimization step:  880    max. rel. force change:+3.462e-02    rel. energy change:+1.276e-02
minimization step:  920    max. rel. force change:+1.828e-02    rel. energy change:+1.126e-02
minimization step:  960    max. rel. force change:+6.457e-02    rel. energy change:+1.061e-02

100%|████████████████████████████████████████████████████████████████████| 100/100 [00:01<00:00, 51.85it/s]
100%|████████████████████████████████████████████████████████████████████| 100/100 [00:55<00:00,  1.80it/s]

simulation for parameters {'counterion_valency': 1, 'rod_charge_dens': 2} done

Question

• Why does the second simulation take much longer than the first one?

The rod charge density is doubled, so the total charge of the counterions needs to be doubled, too. Since their valency is only half of the one in the first run, there will be four times more counterions in the second run.

We plot the density of counterions around the rod as the normalized integrated radial counterion charge distribution function $P(r)$, meaning the integrated probability to find an amount of charge within the radius $r$. We express the rod charge density $\lambda$ in terms of the dimensionless Manning parameter $\xi = \lambda l_B / e$ where $l_B$ is the Bjerrum length and $e$ the elementary charge

# With the notion of P(r) the probability to find the charge up to r,
# we only use the right side of the bin edges for plotting
rs = bin_edges[1:, 0, 0, 0]

fig, ax = plt.subplots(figsize=(10, 7))
for run in runs:
    hist = np.array(run['histogram'][:, 0, 0])
    # The CylindricalDensityProfile normalizes the bin values by the bin size.
    # We want the 'raw' distribution (number of ions within a radius)
    # so we need to multiply by the radii
    hist = hist * rs
    cum_hist = np.cumsum(hist)
    cum_hist /= cum_hist[-1]
    manning_xi = run['params']['rod_charge_dens'] * BJERRUM_LENGTH / Q_E
    ax.plot(rs, cum_hist, label=rf'$\xi ={manning_xi}, \nu = {run["params"]["counterion_valency"]}$')
ax.set_xscale('log')
ax.legend()
plt.xlabel('r')
plt.ylabel('P(r)')
plt.show()

In the semilogarithmic plot we see an inflection point of the cumulative charge distribution which is the indicator for ion condensation. To compare to the meanfield approach of PB, we calculate the solution of the analytical expressions given in 10.1021/ma990897o

def eq_to_solve_for_gamma(gamma, manning_parameter, rod_radius, max_radius):
    # eq 7 - eq 6 from 10.1021/ma990897o
    return gamma * np.log(max_radius / rod_radius) - np.arctan(1 / gamma) + np.arctan((1 - manning_parameter) / gamma)


def calc_manning_radius(gamma, max_radius):
    # eq 7 from 10.1021/ma990897o
    return max_radius * np.exp(-np.arctan(1. / gamma) / gamma)


def calc_PB_probability(r, manning_parameter, gamma, manning_radius):
    # eq 8 and 9 from 10.1021/ma990897o
    return 1. / manning_parameter + gamma / manning_parameter * np.tan(gamma * np.log(r / manning_radius))

For multivalent counterions, the manning parameter $\xi$ has to be multiplied by the valency $\nu$. The result depends only on the product of rod_charge_dens and ion_valency, so we only need one curve

rod_charge_density = runs[0]['params']['rod_charge_dens']
ion_valency = runs[0]['params']['counterion_valency']
manning_parameter_times_valency = BJERRUM_LENGTH * rod_charge_density * ion_valency

gamma = scipy.optimize.fsolve(eq_to_solve_for_gamma, 1, args=(
    manning_parameter_times_valency, r_min, r_max))
manning_radius = calc_manning_radius(gamma, r_max)

PB_probability = calc_PB_probability(
    rs, manning_parameter_times_valency, gamma, manning_radius)

ax.plot(rs, PB_probability, label=rf'PB $\xi \cdot \nu$ = {manning_parameter_times_valency}')
ax.legend()
ax.set_xscale('log')
fig

We see that overall agreement is quite good, but the deviations from the PB solution get stronger the more charged the ions are. Poisson Boltzmann makes two simplifying assumptions: Particles are points and there are no correlations between the particles. Both is not given in the simulation. Excluded volume effects can only lower the density, but we see in the figure that the simulated density is always larger that the calculated one. This means that correlation effects cause the discrepancy.

Overcharging by added salt

Above simulations were performed for a system where all ions come from dissociation off the polyelectrolyte. We can also investigate systems where there are additional salt ions present.

def add_salt(system, anion_params, cation_params):

    N_anions = anion_params['number']
    N_cations = cation_params['number']

    anion_positions = np.random.random((N_anions, 3)) * system.box_l
    cation_positions = np.random.random((N_cations, 3)) * system.box_l

    anions = system.part.add(pos=anion_positions, type=[anion_params['type']] * N_anions,
                             q=[-anion_params['valency']] * N_anions)
    cations = system.part.add(pos=cation_positions, type=[cation_params['type']] * N_cations,
                              q=[cation_params['valency']] * N_cations)

    return anions, cations

ANION_PARAMS = {'type': 3,
                'valency': 2,
                'number': 150}
CATION_PARAMS = {'type': 4,
                 'valency': 2,
                 'number': 150}
ROD_LENGTH = 10
N_rod_beads = int(ROD_LENGTH / ROD_RADIUS)
ROD_CHARGE_DENS = 1
COUNTERION_VALENCY = 1

STEPS_PER_SAMPLE_SALT = 20
N_SAMPLES_SALT = 1500
N_RADIAL_BINS = 100

all_ion_types = [COUNTERION_TYPE, ANION_PARAMS['type'], CATION_PARAMS['type']]

# set interactions of salt with the rod and all ions
for salt_type in [ANION_PARAMS['type'], CATION_PARAMS['type']]:
    system.non_bonded_inter[salt_type, ROD_TYPE].wca.set_params(
        epsilon=WCA_EPSILON, sigma=ION_DIAMETER / 2. + ROD_RADIUS)
    for ion_type in all_ion_types:
        system.non_bonded_inter[salt_type, ion_type].wca.set_params(
            epsilon=WCA_EPSILON, sigma=ION_DIAMETER)

clear_system(system)
system.box_l = 3 * [ROD_LENGTH]
counterions = setup_rod_and_counterions(
    system, COUNTERION_VALENCY, COUNTERION_TYPE,
    ROD_CHARGE_DENS, N_rod_beads, ROD_TYPE)
anions, cations = add_salt(system, ANION_PARAMS, CATION_PARAMS)
assert abs(sum(anions.q) + sum(cations.q)) < 1e-10

p3m = espressomd.electrostatics.P3M(**p3m_params)
system.actors.add(p3m)
remove_overlap(system, STEEPEST_DESCENT_PARAMS)
system.thermostat.set_langevin(**LANGEVIN_PARAMS)
print('', end='', flush=True)
integrate_system(system, WARMUP_STEPS)
radial_profile_accs, bin_edges = setup_profile_calculation(
    system, STEPS_PER_SAMPLE_SALT, all_ion_types, r_min, N_RADIAL_BINS)
integrate_system(system, N_SAMPLES_SALT * STEPS_PER_SAMPLE_SALT)

WARNING: Statistics of tuning samples is very bad.
WARNING: Statistics of tuning samples is very bad.

CoulombP3M tune parameters: Accuracy goal = 1.00000e-03 prefactor = 1.00000e+00
System: box_l = 5.00000e+01 # charged part = 150 Sum[q_i^2] = 3.00000e+02
mesh cao r_cut_iL    alpha_L     err       rs_err    ks_err    time [ms]
6    7   cao too large for this mesh
10   7   2.37352e-01 8.46940e+00 9.881e-04 7.071e-04 6.902e-04 0.28    
10   6   2.42156e-01 8.29105e+00 9.864e-04 7.071e-04 6.878e-04 0.24    
10   5   2.50805e-01 7.98771e+00 9.948e-04 7.071e-04 6.997e-04 0.23    
10   4   2.68102e-01 7.44127e+00 1.000e-03 7.071e-04 7.071e-04 0.23    
10   3   3.10383e-01 6.36820e+00 9.996e-04 7.071e-04 7.065e-04 0.26    
10   2   4.92000e-01 3.89721e+00 1.029e-03 7.071e-04 7.471e-04 accuracy not achieved
16   7   1.57091e-01 1.31192e+01 9.917e-04 7.071e-04 6.953e-04 0.37    
16   6   1.61280e-01 1.27587e+01 9.899e-04 7.071e-04 6.927e-04 0.33    
16   5   1.67563e-01 1.22525e+01 9.943e-04 7.071e-04 6.990e-04 0.30    
16   4   1.79607e-01 1.13837e+01 9.919e-04 7.071e-04 6.955e-04 0.29    
16   3   2.09978e-01 9.64580e+00 9.968e-04 7.071e-04 7.026e-04 0.29    
16   2   2.68102e-01 7.44127e+00 1.588e-03 7.071e-04 1.422e-03 accuracy not achieved
20   7   1.29093e-01 1.61480e+01 9.981e-04 7.071e-04 7.043e-04 0.54    
20   6   1.32601e-01 1.56965e+01 9.956e-04 7.071e-04 7.009e-04 0.52    
20   5   1.38213e-01 1.50230e+01 9.886e-04 7.071e-04 6.910e-04 0.49    
20   4   1.48036e-01 1.39703e+01 9.934e-04 7.071e-04 6.977e-04 0.47    
20   3   1.74696e-01 1.17231e+01 9.915e-04 7.071e-04 6.951e-04 0.48    
20   2   1.79607e-01 1.13837e+01 2.836e-03 7.071e-04 2.746e-03 accuracy not achieved
26   7   1.02931e-01 2.05145e+01 9.847e-04 7.071e-04 6.853e-04 1.35    
26   6   1.05244e-01 2.00386e+01 9.999e-04 7.071e-04 7.070e-04 1.28    
26   5   1.09870e-01 1.91484e+01 9.883e-04 7.071e-04 6.905e-04 1.24    
26   4   1.17966e-01 1.77625e+01 9.924e-04 7.071e-04 6.963e-04 1.22    
26   3   1.39940e-01 1.48270e+01 9.977e-04 7.071e-04 7.038e-04 1.23    
26   2   1.48036e-01 1.39703e+01 2.764e-03 7.071e-04 2.673e-03 accuracy not achieved
30   7   9.12392e-02 2.32992e+01 9.643e-04 7.071e-04 6.557e-04 2.00    
30   6   9.30824e-02 2.28126e+01 9.878e-04 7.071e-04 6.897e-04 1.93    
30   5   9.67688e-02 2.18962e+01 9.931e-04 7.071e-04 6.973e-04 1.84    
30   4   1.04142e-01 2.02627e+01 9.936e-04 7.071e-04 6.981e-04 1.83    
30   3   1.17966e-01 1.77625e+01 1.136e-03 7.071e-04 8.893e-04 accuracy not achieved
36   7   7.72927e-02 2.77545e+01 9.818e-04 7.071e-04 6.811e-04 3.46    
36   6   7.97335e-02 2.68594e+01 9.648e-04 7.071e-04 6.564e-04 3.32    
36   5   8.29879e-02 2.57498e+01 9.690e-04 7.071e-04 6.625e-04 3.66    
36   4   8.94968e-02 2.37781e+01 9.724e-04 7.071e-04 6.675e-04 3.37    
36   3   1.04142e-01 2.02627e+01 1.055e-03 7.071e-04 7.829e-04 accuracy not achieved

resulting parameters: mesh: (10, 10, 10), cao: 4, r_cut_iL: 2.6810e-01,
                      alpha_L: 7.4413e+00, accuracy: 9.9996e-04, time: 0.23
CoulombP3M tune parameters: Accuracy goal = 1.00000e-03 prefactor = 1.00000e+00
System: box_l = 1.00000e+01 # charged part = 320 Sum[q_i^2] = 1.22000e+03
mesh cao r_cut_iL    alpha_L     err       rs_err    ks_err    time [ms]
8    7   4.60000e-01 6.13067e+00 2.012e-02 7.071e-04 2.010e-02 accuracy not achieved
8    7   4.60000e-01 6.13067e+00 2.012e-02 7.071e-04 2.010e-02 accuracy not achieved
10   7   4.60000e-01 6.13067e+00 2.984e-03 7.071e-04 2.899e-03 accuracy not achieved
10   7   4.60000e-01 6.13067e+00 2.984e-03 7.071e-04 2.899e-03 accuracy not achieved
12   7   4.43828e-01 6.36120e+00 9.961e-04 7.071e-04 7.016e-04 2.22    
12   6   4.60000e-01 6.13067e+00 1.144e-03 7.071e-04 8.991e-04 accuracy not achieved
12   7   4.43828e-01 6.36120e+00 9.961e-04 7.071e-04 7.016e-04 2.20    
12   6   4.43828e-01 6.36120e+00 1.453e-03 7.071e-04 1.269e-03 accuracy not achieved
14   7   3.84882e-01 7.36815e+00 9.975e-04 7.071e-04 7.036e-04 1.74    
14   6   4.10021e-01 6.90279e+00 9.944e-04 7.071e-04 6.991e-04 2.46    
14   5   4.43828e-01 6.36120e+00 1.171e-03 7.071e-04 9.339e-04 accuracy not achieved
14   7   3.84882e-01 7.36815e+00 9.975e-04 7.071e-04 7.036e-04 1.70    
14   6   3.84882e-01 7.36815e+00 1.465e-03 7.071e-04 1.283e-03 accuracy not achieved
16   7   3.40531e-01 8.35944e+00 9.919e-04 7.071e-04 6.957e-04 1.40    
16   6   3.62331e-01 7.84143e+00 1.000e-03 7.071e-04 7.071e-04 1.51    
16   5   3.84882e-01 7.36815e+00 1.307e-03 7.071e-04 1.100e-03 accuracy not achieved
16   7   3.40531e-01 8.35944e+00 9.919e-04 7.071e-04 6.957e-04 1.40    
16   6   3.40531e-01 8.35944e+00 1.463e-03 7.071e-04 1.281e-03 accuracy not achieved
18   7   3.05280e-01 9.35607e+00 9.936e-04 7.071e-04 6.980e-04 1.41    
18   6   3.25233e-01 8.76502e+00 9.999e-04 7.071e-04 7.070e-04 1.44    
18   5   3.40531e-01 8.35944e+00 1.427e-03 7.071e-04 1.240e-03 accuracy not achieved
18   7   3.05280e-01 9.35607e+00 9.936e-04 7.071e-04 6.980e-04 1.37    
18   6   3.05280e-01 9.35607e+00 1.476e-03 7.071e-04 1.295e-03 accuracy not achieved
20   7   2.76660e-01 1.03550e+01 9.986e-04 7.071e-04 7.051e-04 1.18    
20   6   2.95740e-01 9.66727e+00 9.925e-04 7.071e-04 6.965e-04 1.26    
20   5   3.05280e-01 9.35607e+00 1.548e-03 7.071e-04 1.378e-03 accuracy not achieved
20   7   2.76660e-01 1.03550e+01 9.986e-04 7.071e-04 7.051e-04 1.17    
20   6   2.76660e-01 1.03550e+01 1.495e-03 7.071e-04 1.317e-03 accuracy not achieved
22   7   2.53425e-01 1.13346e+01 9.953e-04 7.071e-04 7.005e-04 1.90    
22   6   2.70716e-01 1.05894e+01 9.971e-04 7.071e-04 7.030e-04 1.41    
22   5   2.76660e-01 1.03550e+01 1.668e-03 7.071e-04 1.511e-03 accuracy not achieved
22   7   2.53268e-01 1.13418e+01 9.989e-04 7.071e-04 7.055e-04 1.42    
22   6   2.70716e-01 1.05894e+01 9.971e-04 7.071e-04 7.030e-04 1.40    
22   5   2.70716e-01 1.05894e+01 1.915e-03 7.071e-04 1.780e-03 accuracy not achieved
24   7   2.33704e-01 1.23212e+01 9.975e-04 7.071e-04 7.036e-04 1.40    
24   6   2.50095e-01 1.14901e+01 9.943e-04 7.071e-04 6.991e-04 1.40    
24   5   2.70716e-01 1.05894e+01 1.218e-03 7.071e-04 9.915e-04 accuracy not achieved
24   7   2.34464e-01 1.22800e+01 9.795e-04 7.071e-04 6.778e-04 1.41    
24   6   2.50095e-01 1.14901e+01 9.943e-04 7.071e-04 6.991e-04 1.40    
24   5   2.50095e-01 1.14901e+01 1.922e-03 7.071e-04 1.787e-03 accuracy not achieved
26   7   2.17857e-01 1.32454e+01 9.756e-04 7.071e-04 6.722e-04 1.81    
26   6   2.32511e-01 1.23864e+01 9.920e-04 7.071e-04 6.957e-04 1.72    
26   5   2.50095e-01 1.14901e+01 1.261e-03 7.071e-04 1.044e-03 accuracy not achieved
26   7   2.17070e-01 1.32948e+01 9.955e-04 7.071e-04 7.007e-04 1.80    
26   6   2.32511e-01 1.23864e+01 9.920e-04 7.071e-04 6.957e-04 1.75    
26   5   2.32511e-01 1.23864e+01 1.928e-03 7.071e-04 1.794e-03 accuracy not achieved
28   7   2.02539e-01 1.42783e+01 9.988e-04 7.071e-04 7.054e-04 2.01    
28   6   2.17070e-01 1.32948e+01 9.953e-04 7.071e-04 7.004e-04 1.84    
28   5   2.32511e-01 1.23864e+01 1.300e-03 7.071e-04 1.091e-03 accuracy not achieved
28   7   2.02656e-01 1.42698e+01 9.955e-04 7.071e-04 7.008e-04 2.05    
28   6   2.17070e-01 1.32948e+01 9.953e-04 7.071e-04 7.004e-04 2.11    
28   5   2.17070e-01 1.32948e+01 1.949e-03 7.071e-04 1.817e-03 accuracy not achieved
30   7   1.89990e-01 1.52505e+01 9.987e-04 7.071e-04 7.053e-04 2.50    
30   6   2.02656e-01 1.42698e+01 1.021e-03 7.071e-04 7.368e-04 accuracy not achieved
30   7   1.89990e-01 1.52505e+01 9.987e-04 7.071e-04 7.053e-04 2.43    
30   6   1.89990e-01 1.52505e+01 1.532e-03 7.071e-04 1.359e-03 accuracy not achieved
32   7   1.79600e-01 1.61598e+01 9.787e-04 7.071e-04 6.767e-04 2.39    
32   6   1.89990e-01 1.52505e+01 1.049e-03 7.071e-04 7.745e-04 accuracy not achieved
32   7   1.79600e-01 1.61598e+01 9.787e-04 7.071e-04 6.767e-04 2.44    
32   6   1.79600e-01 1.61598e+01 1.497e-03 7.071e-04 1.320e-03 accuracy not achieved
34   7   1.69778e-01 1.71231e+01 9.788e-04 7.071e-04 6.768e-04 5.99    
34   6   1.79600e-01 1.61598e+01 1.052e-03 7.071e-04 7.789e-04 accuracy not achieved

resulting parameters: mesh: (20, 20, 20), cao: 7, r_cut_iL: 2.7666e-01,
                      alpha_L: 1.0355e+01, accuracy: 9.9856e-04, time: 1.17
minimization step:   40    max. rel. force change:+9.562e-01    rel. energy change:+9.299e-01
minimization step:   80    max. rel. force change:+7.655e-01    rel. energy change:+1.979e-01
minimization step:  120    max. rel. force change:+6.491e-03    rel. energy change:+4.764e-03

100%|████████████████████████████████████████████████████████████████████| 100/100 [00:06<00:00, 15.65it/s]
100%|████████████████████████████████████████████████████████████████████| 100/100 [00:38<00:00,  2.62it/s]

rs = bin_edges[1:, 0, 0, 0]
cum_hists = {}
for ion_type in all_ion_types:
    hist = radial_profile_accs[ion_type].mean()
    hist = hist[:, 0, 0] * rs
    cum_hist = np.cumsum(hist)
    cum_hist /= cum_hist[-1]
    cum_hists[ion_type] = cum_hist

Exercise:

• Use the cumulative histograms from the cell above to create the cumulative charge histogram of the total ion charge

Hints

• You need to account for the fact that the cumulative histograms are all normalized, but the total charge of each ion type is different

counterion_charge = sum(counterions.q)
anion_charge = sum(anions.q)
cation_charge = sum(cations.q)
charge_hist = counterion_charge * cum_hists[COUNTERION_TYPE] + \
    anion_charge * cum_hists[ANION_PARAMS['type']] + \
    cation_charge * cum_hists[CATION_PARAMS['type']]

charge_hist /= charge_hist[-1]
fig2, ax2 = plt.subplots(figsize=(10, 7))
ax2.plot(rs, charge_hist)
ax2.set_xscale('linear')
plt.xlabel('r')
plt.ylabel('P(r)')
plt.show()

You should observe a strong overcharging effect, where ions accumulate close to the rod.