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
import espressomd.zn
espressomd.assert_features(['ELECTROSTATICS', 'P3M', 'WCA', 'EXTERNAL_FORCES'])
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})
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
MASS=1.0
# particle types
ROD_TYPE = 1
COUNTERION_TYPE = 2
Exercise:
Hints:
sigma parameter of the interaction between the rod particles and the ions# SOLUTION CELL
# 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:
setup_rod_and_counterions(system, ion_valency, counterion_type, rod_charge_dens, N_rod_beads, rod_type)typeq according to their ion_valencyrod_charge_dens is uniformly distributed along the N_rod_beads individual particlesHints:
# SOLUTION CELL
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 solvers for electrostatics, magnetostatics and hydrodynamics. They are attached to system slots in such a way that unphysical combinations of algorithms are avoided, for example simultaneous usage of two electrostatics solvers. 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:
p3m instance and add it to the electrostatics slot of the system# SOLUTION CELL
p3m = espressomd.electrostatics.P3M(**p3m_params)
system.electrostatics.solver = p3m
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 8 7 2.89242e-01 6.86443e+00 9.976e-04 7.071e-04 7.038e-04 1.01 8 6 2.93086e-01 6.76873e+00 9.965e-04 7.071e-04 7.021e-04 0.65 8 5 3.02695e-01 6.54040e+00 9.956e-04 7.071e-04 7.008e-04 0.43 8 4 3.23836e-01 6.08705e+00 9.989e-04 7.071e-04 7.056e-04 0.43 8 3 3.73805e-01 5.22445e+00 9.996e-04 7.071e-04 7.066e-04 0.45 8 2 4.92000e-01 3.89721e+00 1.399e-03 7.071e-04 1.207e-03 accuracy not achieved 10 7 2.37072e-01 8.48000e+00 9.924e-04 7.071e-04 6.964e-04 0.35 10 6 2.41802e-01 8.30397e+00 9.916e-04 7.071e-04 6.951e-04 0.30 10 5 2.50670e-01 7.99229e+00 9.966e-04 7.071e-04 7.023e-04 0.33 10 4 2.68406e-01 7.43231e+00 9.965e-04 7.071e-04 7.022e-04 0.33 10 3 3.02695e-01 6.54040e+00 1.064e-03 7.071e-04 7.951e-04 accuracy not achieved 12 7 2.02131e-01 1.00435e+01 9.876e-04 7.071e-04 6.895e-04 0.35 12 6 2.06854e-01 9.80044e+00 9.884e-04 7.071e-04 6.906e-04 0.31 12 5 2.14410e-01 9.43439e+00 9.995e-04 7.071e-04 7.064e-04 0.29 12 4 2.29523e-01 8.77644e+00 9.988e-04 7.071e-04 7.054e-04 0.27 12 3 2.41802e-01 8.30397e+00 1.324e-03 7.071e-04 1.120e-03 accuracy not achieved 14 7 1.76625e-01 1.15875e+01 9.884e-04 7.071e-04 6.906e-04 0.38 14 6 1.81108e-01 1.12838e+01 9.885e-04 7.071e-04 6.907e-04 0.34 14 5 1.88280e-01 1.08287e+01 9.907e-04 7.071e-04 6.939e-04 0.32 14 4 2.01729e-01 1.00647e+01 9.886e-04 7.071e-04 6.908e-04 0.30 14 3 2.29523e-01 8.77644e+00 1.053e-03 7.071e-04 7.801e-04 accuracy not achieved 16 7 1.56813e-01 1.31439e+01 9.982e-04 7.071e-04 7.046e-04 0.41 16 6 1.61541e-01 1.27369e+01 9.841e-04 7.071e-04 6.845e-04 0.36 16 5 1.67845e-01 1.22307e+01 9.885e-04 7.071e-04 6.908e-04 0.36 16 4 1.79665e-01 1.13799e+01 9.909e-04 7.071e-04 6.941e-04 0.33 16 3 2.01729e-01 1.00647e+01 1.103e-03 7.071e-04 8.467e-04 accuracy not achieved 18 7 1.41767e-01 1.46249e+01 9.917e-04 7.071e-04 6.953e-04 0.57 18 6 1.45276e-01 1.42513e+01 9.978e-04 7.071e-04 7.039e-04 0.53 18 5 1.51592e-01 1.36236e+01 9.867e-04 7.071e-04 6.882e-04 0.51 18 4 1.62119e-01 1.26887e+01 9.941e-04 7.071e-04 6.988e-04 0.48 18 3 1.79665e-01 1.13799e+01 1.159e-03 7.071e-04 9.183e-04 accuracy not achieved 20 7 1.29189e-01 1.61353e+01 9.952e-04 7.071e-04 7.003e-04 0.52 20 6 1.32989e-01 1.56481e+01 9.849e-04 7.071e-04 6.856e-04 0.47 20 5 1.38055e-01 1.50413e+01 9.927e-04 7.071e-04 6.967e-04 0.46 20 4 1.48187e-01 1.39552e+01 9.902e-04 7.071e-04 6.932e-04 0.44 20 3 1.62119e-01 1.26887e+01 1.211e-03 7.071e-04 9.831e-04 accuracy not achieved 22 7 1.19244e-01 1.75613e+01 9.829e-04 7.071e-04 6.828e-04 0.83 22 6 1.22139e-01 1.71217e+01 9.913e-04 7.071e-04 6.947e-04 0.77 22 5 1.26770e-01 1.64611e+01 9.999e-04 7.071e-04 7.070e-04 0.76 22 4 1.36031e-01 1.52781e+01 9.998e-04 7.071e-04 7.069e-04 0.72 22 3 1.48187e-01 1.39552e+01 1.250e-03 7.071e-04 1.030e-03 accuracy not achieved 24 7 1.09994e-01 1.91256e+01 9.998e-04 7.071e-04 7.068e-04 0.88 24 6 1.13182e-01 1.85568e+01 9.908e-04 7.071e-04 6.940e-04 0.81 24 5 1.17965e-01 1.77627e+01 9.851e-04 7.071e-04 6.858e-04 0.81 24 4 1.26467e-01 1.65028e+01 9.919e-04 7.071e-04 6.956e-04 0.76 24 3 1.36031e-01 1.52781e+01 1.302e-03 7.071e-04 1.093e-03 accuracy not achieved 26 7 1.02754e-01 2.05518e+01 9.912e-04 7.071e-04 6.946e-04 1.13 26 6 1.05718e-01 1.99437e+01 9.830e-04 7.071e-04 6.829e-04 1.07 26 5 1.09670e-01 1.91853e+01 9.949e-04 7.071e-04 6.998e-04 1.05 26 4 1.18562e-01 1.76681e+01 9.770e-04 7.071e-04 6.741e-04 1.00 26 3 1.26467e-01 1.65028e+01 1.327e-03 7.071e-04 1.123e-03 accuracy not achieved 28 7 9.63319e-02 2.20011e+01 9.886e-04 7.071e-04 6.908e-04 1.34 28 6 9.91107e-02 2.13503e+01 9.805e-04 7.071e-04 6.792e-04 1.27 28 5 1.02816e-01 2.05387e+01 9.931e-04 7.071e-04 6.973e-04 1.25 28 4 1.11152e-01 1.89151e+01 9.776e-04 7.071e-04 6.750e-04 1.18 28 3 1.18562e-01 1.76681e+01 1.338e-03 7.071e-04 1.136e-03 accuracy not achieved 30 7 9.11796e-02 2.33153e+01 9.667e-04 7.071e-04 6.592e-04 1.54 30 6 9.29163e-02 2.28556e+01 9.946e-04 7.071e-04 6.994e-04 1.43 30 5 9.72582e-02 2.17799e+01 9.752e-04 7.071e-04 6.716e-04 1.43 30 4 1.04205e-01 2.02496e+01 9.917e-04 7.071e-04 6.953e-04 1.33 30 3 1.11152e-01 1.89151e+01 1.368e-03 7.071e-04 1.172e-03 accuracy not achieved 32 7 8.54808e-02 2.49583e+01 9.953e-04 7.071e-04 7.004e-04 1.40 32 6 8.79231e-02 2.42274e+01 9.881e-04 7.071e-04 6.902e-04 1.33 32 5 9.19937e-02 2.30976e+01 9.714e-04 7
WARNING: Statistics of tuning samples is very bad. in function void check_statistics(Utils::Statistics::RunningAverage<double>&) (/builds/espressomd/espresso/src/core/tuning.cpp:62) on node 0 WARNING: Statistics of tuning samples is very bad. in function void check_statistics(Utils::Statistics::RunningAverage<double>&) (/builds/espressomd/espresso/src/core/tuning.cpp:62) on node 0 WARNING: Statistics of tuning samples is very bad. in function void check_statistics(Utils::Statistics::RunningAverage<double>&) (/builds/espressomd/espresso/src/core/tuning.cpp:62) on node 0 WARNING: Statistics of tuning samples is very bad. in function void check_statistics(Utils::Statistics::RunningAverage<double>&) (/builds/espressomd/espresso/src/core/tuning.cpp:62) on node 0 WARNING: Statistics of tuning samples is very bad. in function void check_statistics(Utils::Statistics::RunningAverage<double>&) (/builds/espressomd/espresso/src/core/tuning.cpp:62) on node 0
.071e-04 6.660e-04 1.30
32 4 9.85065e-02 2.14886e+01 9.922e-04 7.071e-04 6.960e-04 1.23
32 3 1.04205e-01 2.02496e+01 1.418e-03 7.071e-04 1.230e-03 accuracy not achieved
34 7 8.15757e-02 2.62201e+01 9.692e-04 7.071e-04 6.628e-04 3.14
34 6 8.31148e-02 2.57083e+01 9.986e-04 7.071e-04 7.051e-04 3.02
34 5 8.69628e-02 2.45098e+01 9.814e-04 7.071e-04 6.806e-04 3.00
34 4 9.38890e-02 2.26058e+01 9.778e-04 7.071e-04 6.753e-04 2.95
34 3 9.85065e-02 2.14886e+01 1.445e-03 7.071e-04 1.261e-03 accuracy not achieved
resulting parameters: mesh: (12, 12, 12), cao: 4, r_cut_iL: 2.2952e-01,
alpha_L: 8.7764e+00, accuracy: 9.9881e-04, time: 0.27
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):
FMAX = 0.01 * ION_DIAMETER * MASS / system.time_step**2
system.integrator.set_steepest_descent(
f_max=FMAX,
gamma=10,
max_displacement=0.01)
system.integrator.run(5000)
assert np.all(np.abs(system.part.all().f) < FMAX), "Overlap removal did not converge!"
system.integrator.set_vv()
remove_overlap(system)
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)
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:
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 calculationion_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 densityradial_distances[counterion_type] = <accumulator> and the edges of the binsHints:
system.part.select(type=...) to get only the particles of a specific type# SOLUTION CELL
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:
clear_system(system) thatHints:
# SOLUTION CELL
def clear_system(system):
system.thermostat.turn_off()
system.part.clear()
system.electrostatics.clear()
system.auto_update_accumulators.clear()
system.time = 0.
clear_system(system)
# Visualizer Parameters
color = {ROD_TYPE: "#0911e8", COUNTERION_TYPE: "#f70519"} # Particle color by type
radii = {ROD_TYPE: ROD_RADIUS, COUNTERION_TYPE: ION_DIAMETER/2.0} # Particle size by type
print("System visualization for case (a) counterion valency: 2, rod_charge_density: 1")
# Initializing Visualizer
vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)
System visualization for case (a) counterion valency: 2, rod_charge_density: 1
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
Exercise:
Hints:
p3m instance after each change of parameters. If we reuse the p3m that was tuned before, likely the desired accuracy will not be achieved. .mean()# SOLUTION CELL
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.electrostatics.solver = p3m
remove_overlap(system)
system.thermostat.set_langevin(**LANGEVIN_PARAMS)
print('', end='', flush=True)
# For longer simulation runs it will be convenient to have a progress bar
for i in tqdm.trange(100):
system.integrator.run(WARMUP_STEPS // 100)
radial_profile_accs, bin_edges = setup_profile_calculation(
system, STEPS_PER_SAMPLE, [COUNTERION_TYPE], r_min, N_RADIAL_BINS)
for i in tqdm.trange(100):
system.integrator.run(N_SAMPLES * STEPS_PER_SAMPLE // 100)
# Updating frames for the visualizer
if(run['params']['counterion_valency']==2):
vis.update()
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 6 7 cao too large for this mesh 8 7 2.65219e-01 7.19254e+00 9.878e-04 7.071e-04 6.897e-04 0.20 8 6 2.68102e-01 7.10991e+00 9.881e-04 7.071e-04 6.902e-04 0.18 8 5 2.74828e-01 6.92406e+00 9.971e-04 7.071e-04 7.031e-04 0.17 8 4 2.92125e-01 6.48657e+00 9.991e-04 7.071e-04 7.058e-04 0.17 8 3 3.32484e-01 5.64760e+00 9.989e-04 7.071e-04 7.055e-04 0.17 8 2 4.92000e-01 3.70898e+00 1.035e-03 7.071e-04 7.561e-04 accuracy not achieved 10 7 2.16894e-01 8.91577e+00 9.914e-04 7.071e-04 6.948e-04 0.21 10 6 2.20141e-01 8.77554e+00 9.998e-04 7.071e-04 7.069e-04 0.18 10 5 2.27934e-01 8.45574e+00 9.973e-04 7.071e-04 7.033e-04 0.17 10 4 2.42869e-01 7.90177e+00 9.958e-04 7.071e-04 7.012e-04 0.17 10 3 2.76637e-01 6.87566e+00 9.965e-04 7.071e-04 7.021e-04 0.17 10 2 3.32484e-01 5.64760e+00 1.530e-03 7.071e-04 1.357e-03 accuracy not achieved 12 7 1.84306e-01 1.06058e+01 9.997e-04 7.071e-04 7.067e-04 0.22 12 6 1.88758e-01 1.03395e+01 9.917e-04 7.071e-04 6.953e-04 0.20 12 5 1.95880e-01 9.93926e+00 9.903e-04 7.071e-04 6.933e-04 0.20 12 4 2.08346e-01 9.30650e+00 9.941e-04 7.071e-04 6.987e-04 0.17 12 3 2.27934e-01 8.45574e+00 1.112e-03 7.071e-04 8.584e-04 accuracy not achieved 14 7 1.61956e-01 1.21711e+01 9.835e-04 7.071e-04 6.836e-04 0.27 14 6 1.65212e-01 1.19159e+01 9.965e-04 7.071e-04 7.021e-04 0.24 14 5 1.71722e-01 1.14355e+01 9.920e-04 7.071e-04 6.957e-04 0.24 14 4 1.83116e-01 1.06792e+01 9.881e-04 7.071e-04 6.901e-04 0.21 14 3 2.08346e-01 9.30650e+00 1.006e-03 7.071e-04 7.153e-04 accuracy not achieved 16 7 1.43775e-01 1.38148e+01 9.956e-04 7.071e-04 7.009e-04 0.32 16 6 1.47351e-01 1.34585e+01 9.963e-04 7.071e-04 7.019e-04 0.29 16 5 1.53074e-01 1.29240e+01 9.941e-04 7.071e-04 6.988e-04 0.27 16 4 1.63088e-01 1.20812e+01 9.933e-04 7.071e-04 6.976e-04 0.26 16 3 1.83116e-01 1.06792e+01 1.052e-03 7.071e-04 7.787e-04 accuracy not achieved 18 7 1.29961e-01 1.53808e+01 9.930e-04 7.071e-04 6.972e-04 0.50 18 6 1.33146e-01 1.49900e+01 9.965e-04 7.071e-04 7.021e-04 0.46 18 5 1.38242e-01 1.44033e+01 9.961e-04 7.071e-04 7.016e-04 0.44 18 4 1.47161e-01 1.34770e+01 9.990e-04 7.071e-04 7.057e-04 0.42 18 3 1.63088e-01 1.20812e+01 1.103e-03 7.071e-04 8.468e-04 accuracy not achieved 20 7 1.18419e-01 1.69778e+01 9.995e-04 7.071e-04 7.065e-04 0.46 20 6 1.21868e-01 1.64680e+01 9.880e-04 7.071e-04 6.900e-04 0.42 20 5 1.26467e-01 1.58326e+01 9.897e-04 7.071e-04 6.925e-04 0.41 20 4 1.34515e-01 1.48280e+01 9.972e-04 7.071e-04 7.031e-04 0.38 20 3 1.47161e-01 1.34770e+01 1.151e-03 7.071e-04 9.086e-04 accuracy not achieved 22 7 1.09293e-01 1.84864e+01 9.910e-04 7.071e-04 6.943e-04 0.76 22 6 1.11920e-01 1.80262e+01 9.970e-04 7.071e-04 7.029e-04 0.71 22 5 1.16124e-01 1.73343e+01 9.993e-04 7.071e-04 7.061e-04 0.71 22 4 1.24006e-01 1.61666e+01 9.953e-04 7.071e-04 7.004e-04 0.67 22 3 1.34515e-01 1.48280e+01 1.188e-03 7.071e-04 9.540e-04 accuracy not achieved 24 7 1.01723e-01 1.99492e+01 9.778e-04 7.071e-04 6.753e-04 0.82 24 6 1.04630e-01 1.93619e+01 9.689e-04 7.071e-04 6.624e-04 0.76 24 5 1.08505e-01 1.86290e+01 9.739e-04 7.071e-04 6.696e-04 0.75 24 4 1.15287e-01 1.74680e+01 9.891e-04 7.071e-04 6.916e-04 0.72 24 3 1.24006e-01 1.61666e+01 1.220e-03 7.071e-04 9.939e-04 accuracy not achieved 26 7 9.45710e-02 2.15527e+01 9.890e-04 7.071e-04 6.914e-04 1.06 26 6 9.72730e-02 2.09185e+01 9.797e-04 7.071e-04 6.781e-04 0.99 26 5 1.00876e-01 2.01271e+01 9.848e-04 7.071e-04 6.854e-04 0.98 26 4 1.08081e-01 1.87065e+01 9.758e-04 7.071e-04 6.724e-04 0.94 26 3
100%|██████████| 100/100 [00:00<00:00, 141.27it/s] 100%|██████████| 100/100 [00:22<00:00, 4.50it/s] WARNING: Statistics of tuning samples is very bad. in function void check_statistics(Utils::Statistics::RunningAverage<double>&) (/builds/espressomd/espresso/src/core/tuning.cpp:62) on node 0
simulation for parameters {'counterion_valency': 2, 'rod_charge_dens': 1} done
1.15287e-01 1.74680e+01 1.243e-03 7.071e-04 1.023e-03 accuracy not achieved
28 7 8.86603e-02 2.30787e+01 9.887e-04 7.071e-04 6.910e-04 1.25
28 6 9.11935e-02 2.23999e+01 9.794e-04 7.071e-04 6.776e-04 1.21
28 5 9.45710e-02 2.15527e+01 9.849e-04 7.071e-04 6.856e-04 1.19
28 4 1.01326e-01 2.00322e+01 9.775e-04 7.071e-04 6.749e-04 1.12
28 3 1.08081e-01 1.87065e+01 1.253e-03 7.071e-04 1.035e-03 accuracy not achieved
resulting parameters: mesh: (8, 8, 8), cao: 3, r_cut_iL: 3.3248e-01,
alpha_L: 5.6476e+00, accuracy: 9.9890e-04, time: 0.17
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
8 7 2.89242e-01 6.86443e+00 9.976e-04 7.071e-04 7.038e-04 0.39
8 6 2.93086e-01 6.76873e+00 9.965e-04 7.071e-04 7.021e-04 0.36
8 5 3.02695e-01 6.54040e+00 9.956e-04 7.071e-04 7.008e-04 0.35
8 4 3.23836e-01 6.08705e+00 9.989e-04 7.071e-04 7.056e-04 0.35
8 3 3.73805e-01 5.22445e+00 9.996e-04 7.071e-04 7.066e-04 0.40
8 2 4.92000e-01 3.89721e+00 1.399e-03 7.071e-04 1.207e-03 accuracy not achieved
10 7 2.37072e-01 8.48000e+00 9.924e-04 7.071e-04 6.964e-04 0.33
10 6 2.41802e-01 8.30397e+00 9.916e-04 7.071e-04 6.951e-04 0.29
10 5 2.50670e-01 7.99229e+00 9.966e-04 7.071e-04 7.023e-04 0.33
10 4 2.68406e-01 7.43231e+00 9.965e-04 7.071e-04 7.022e-04 0.34
10 3 3.02695e-01 6.54040e+00 1.064e-03 7.071e-04 7.951e-04 accuracy not achieved
12 7 2.02131e-01 1.00435e+01 9.876e-04 7.071e-04 6.895e-04 0.34
12 6 2.06854e-01 9.80044e+00 9.884e-04 7.071e-04 6.906e-04 0.30
12 5 2.14410e-01 9.43439e+00 9.995e-04 7.071e-04 7.064e-04 0.29
12 4 2.29523e-01 8.77644e+00 9.988e-04 7.071e-04 7.054e-04 0.27
12 3 2.41802e-01 8.30397e+00 1.324e-03 7.071e-04 1.120e-03 accuracy not achieved
14 7 1.76625e-01 1.15875e+01 9.884e-04 7.071e-04 6.906e-04 0.38
14 6 1.81108e-01 1.12838e+01 9.885e-04 7.071e-04 6.907e-04 0.33
14 5 1.88280e-01 1.08287e+01 9.907e-04 7.071e-04 6.939e-04 0.31
14 4 2.01729e-01 1.00647e+01 9.886e-04 7.071e-04 6.908e-04 0.30
14 3 2.29523e-01 8.77644e+00 1.053e-03 7.071e-04 7.801e-04 accuracy not achieved
16 7 1.56813e-01 1.31439e+01 9.982e-04 7.071e-04 7.046e-04 0.41
16 6 1.61541e-01 1.27369e+01 9.841e-04 7.071e-04 6.845e-04 0.37
16 5 1.67845e-01 1.22307e+01 9.885e-04 7.071e-04 6.908e-04 0.35
16 4 1.79665e-01 1.13799e+01 9.909e-04 7.071e-04 6.941e-04 0.32
16 3 2.01729e-01 1.00647e+01 1.103e-03 7.071e-04 8.467e-04 accuracy not achieved
18 7 1.41767e-01 1.46249e+01 9.917e-04 7.071e-04 6.953e-04 0.58
18 6 1.45276e-01 1.42513e+01 9.978e-04 7.071e-04 7.039e-04 0.53
18 5 1.51592e-01 1.36236e+01 9.867e-04 7.071e-04 6.882e-04 0.51
18 4 1.62119e-01 1.26887e+01 9.941e-04 7.071e-04 6.988e-04 0.47
18 3 1.79665e-01 1.13799e+01 1.159e-03 7.071e-04 9.183e-04 accuracy not achieved
20 7 1.29189e-01 1.61353e+01 9.952e-04 7.071e-04 7.003e-04 0.55
20 6 1.32989e-01 1.56481e+01 9.849e-04 7.071e-04 6.856e-04 0.49
20 5 1.38055e-01 1.50413e+01 9.927e-04 7.071e-04 6.967e-04 0.47
20 4 1.48187e-01 1.39552e+01 9.902e-04 7.071e-04 6.932e-04 0.44
20 3 1.62119e-01 1.26887e+01 1.211e-03 7.071e-04 9.831e-04 accuracy not achieved
22 7 1.19244e-01 1.75613e+01 9.829e-04 7.071e-04 6.828e-04 0.83
22 6 1.22139e-01 1.71217e+01 9.913e-04 7.071e-04 6.947e-04 0.79
22 5 1.26770e-01 1.64611e+01 9.999e-04 7.071e-04 7.070e-04 0.77
22 4 1.36031e-01 1.52781e+01 9.998e-04 7.071e-04 7.069e-04 0.73
22 3 1.48187e-01 1.39552e+01 1.250e-03 7.071e-04 1.030e-03 accuracy not achieved
24 7 1.09994e-01 1.91256e+01 9.998e-04 7.071e-04 7.068e-04 0.88
24 6 1.13182e-01 1.85568e+01 9.908e-04 7.071e-04 6.940e-04 0.83
24 5 1.17965e-01 1.77627e+01 9.851e-04 7.071e-04 6.858e-04 0.81
24 4 1.26467e-01 1.65028e+01 9.919e-04 7.071e-04 6.956e-04 0.79
24 3 1.36031e-01 1.52781e+01 1.302e-03 7.071e-04 1.093e-03 accuracy not achieved
26 7 1.02754e-01 2.05518e+01 9.912e-04 7.071e-04 6.946e-04 1.14
26 6 1.05718e-01 1.99437e+01 9.830e-04 7.071e-04 6.829e-04 1.07
26 5 1.09670e-01 1.91853e+01 9.949e-04 7.071e-04 6.998e-04 1.08
26 4 1.18562e-01 1.76681e+01 9.770e-04 7.071e-04 6.741e-04 1.00
26 3 1.26467e-01 1.65028e+01 1.327e-03 7.071e-04 1.123e-03 accuracy not achieved
28 7 9.63319e-02 2.20011e+01 9.886e-04 7.071e-04 6.908e-04 1.34
28 6 9.91107e-02 2.13503e+01 9.805e-04 7.071e-04 6.792e-04 1.27
28 5 1.02816e-01 2.05387e+01 9.931e-04 7.071e-04 6.973e-04 1.25
28 4 1.11152e-01 1.89151e+01 9.776e-04 7.071e-04 6.750e-04 1.19
28 3 1.18562e-01 1.76681e+01 1.338e-03 7.071e-04 1.136e-03 accuracy not achieved
30 7 9.11796e-02 2.33153e+01 9.667e-04 7.071e-04 6.592e-04 1.54
30 6 9.29163e-02 2.28556e+01 9.946e-04 7.071e-04 6.994e-04 1.46
30 5 9.72582e-02 2.17799e+01 9.752e-04 7.071e-04 6.716e-04 1.42
30 4 1.04205e-01 2.02496e+01 9.917e-04 7.071e-04 6.953e-04 1.34
30 3 1.11152e-01 1.89151e+01 1.368e-03 7.071e-04 1.172e-03 accuracy not achieved
32 7 8.54808e-02 2.49583e+01 9.953e-04 7.071e-04 7.004e-04 1.41
32 6 8.79231e-02 2.42274e+01 9.881e-04 7.071e-04 6.902e-04 1.36
32 5 9.19937e-02 2.30976e+01 9.714e-04 7.071e-04 6.660e-04 1.31
32 4 9.85065e-02 2.14886e+01 9.922e-04 7.071e-04 6.960e-04 1.20
32 3 1.04205e-01 2.02496e+01 1.418e-03 7.071e-04 1.230e-03 accuracy not achieved
34 7 8.15757e-02 2.62201e+01 9.692e-04 7.071e-04 6.628e-04 3.22
34 6 8.31148e-02 2.57083e+01 9.986e-04 7.071e-04 7.051e-04 3.08
34 5 8.69628e-02 2.45098e+01 9.814e-04 7.071e-04 6.806e-04 3.05
34 4 9.38890e-02 2.26058e+01 9.778e-04 7.071e-04 6.753e-04 2.98
34 3 9.85065e-02 2.14886e+01 1.445e-03 7.071e-04 1.261e-03 accuracy not achieved
resulting parameters: mesh: (12, 12, 12), cao: 4, r_cut_iL: 2.2952e-01,
alpha_L: 8.7764e+00, accuracy: 9.9881e-04, time: 0.27
100%|██████████| 100/100 [00:01<00:00, 83.72it/s] 100%|██████████| 100/100 [00:43<00:00, 2.32it/s]
simulation for parameters {'counterion_valency': 1, 'rod_charge_dens': 2} done
Question
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.
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)
# Resetting system with salt ions
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.electrostatics.solver = p3m
remove_overlap(system)
system.thermostat.set_langevin(**LANGEVIN_PARAMS)
print('', end='', flush=True)
for i in tqdm.trange(100):
system.integrator.run(WARMUP_STEPS // 100)
radial_profile_accs, bin_edges = setup_profile_calculation(
system, STEPS_PER_SAMPLE_SALT, all_ion_types, r_min, N_RADIAL_BINS)
for i in tqdm.trange(100):
system.integrator.run(N_SAMPLES_SALT * STEPS_PER_SAMPLE_SALT // 100)
100%|██████████| 100/100 [00:05<00:00, 18.22it/s] 100%|██████████| 100/100 [00:33<00:00, 3.00it/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:
Hints
# SOLUTION CELL
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.