# Tutorial 2: A Simple Charged System, Part 1¶

## 1 Introduction¶

This tutorial introduces some of the basic features of ESPResSo for charged systems by constructing a simulation script for a simple salt crystal. In the subsequent task, we use a more realistic force-field for a NaCl crystal. Finally, we introduce constraints and 2D-Electrostatics to simulate a molten salt in a parallel plate capacitor. We assume that the reader is familiar with the basic concepts of Python and MD simulations. Compile ESPResSo with the following features in your myconfig.hpp to be set throughout the whole tutorial:

#define EXTERNAL_FORCES
#define MASS
#define ELECTROSTATICS
#define WCA


## 2 Basic Set Up¶

The script for the tutorial can be found in your build directory at /doc/tutorials/02-charged_system/scripts/nacl.py. We start by importing numpy, pyplot, espressomd and setting up the simulation parameters:

In :
from espressomd import System, electrostatics
import espressomd
import numpy
import matplotlib.pyplot as plt

In :
plt.ion()

# Print enabled features
required_features = ["EXTERNAL_FORCES", "MASS", "ELECTROSTATICS", "WCA"]
espressomd.assert_features(required_features)
print(espressomd.features())

# System Parameters
n_part = 200
n_ionpairs = n_part / 2
density = 0.5
time_step = 0.01
temp = 1.0
gamma = 1.0
l_bjerrum = 7.0

num_steps_equilibration = 200
num_configs = 100
integ_steps_per_config = 250

# Particle Parameters
types       = {"Anion":          0, "Cation": 1}
numbers     = {"Anion": n_ionpairs, "Cation": n_ionpairs}
charges     = {"Anion":       -1.0, "Cation": 1.0}
wca_sigmas   = {"Anion":        1.0, "Cation": 1.0}
wca_epsilons = {"Anion":        1.0, "Cation": 1.0}

['ADDITIONAL_CHECKS', 'BLAS', 'BMHTF_NACL', 'BOND_CONSTRAINT', 'BROWNIAN_PER_PARTICLE', 'BUCKINGHAM', 'COLLISION_DETECTION', 'CUDA', 'DIPOLAR_BARNES_HUT', 'DIPOLAR_DIRECT_SUM', 'DIPOLES', 'DP3M', 'DPD', 'EK_BOUNDARIES', 'ELECTROKINETICS', 'ELECTROSTATICS', 'ENGINE', 'EXCLUSIONS', 'EXPERIMENTAL_FEATURES', 'EXTERNAL_FORCES', 'FFTW', 'GAUSSIAN', 'GAY_BERNE', 'GSL', 'H5MD', 'HAT', 'HERTZIAN', 'LANGEVIN_PER_PARTICLE', 'LAPACK', 'LB_BOUNDARIES', 'LB_BOUNDARIES_GPU', 'LB_ELECTROHYDRODYNAMICS', 'LENNARD_JONES', 'LENNARD_JONES_GENERIC', 'LJCOS', 'LJCOS2', 'LJGEN_SOFTCORE', 'MASS', 'MMM1D_GPU', 'MORSE', 'NPT', 'P3M', 'PARTICLE_ANISOTROPY', 'ROTATION', 'ROTATIONAL_INERTIA', 'SCAFACOS', 'SCAFACOS_DIPOLES', 'SMOOTH_STEP', 'SOFT_SPHERE', 'STOKESIAN_DYNAMICS', 'TABULATED', 'THOLE', 'VIRTUAL_SITES', 'VIRTUAL_SITES_INERTIALESS_TRACERS', 'VIRTUAL_SITES_RELATIVE', 'WCA']


These variables do not change anything in the simulation engine, but are just standard Python variables. They are used to increase the readability and flexibility of the script. The box length is not a parameter of this simulation, it is calculated from the number of particles and the system density. This allows to change the parameters later easily, e.g. to simulate a bigger system. We use dictionaries for all particle related parameters, which is less error-prone and readable as we will see later when we actually need the values. The parameters here define a purely repulsive, equally sized, monovalent salt.

The simulation engine itself is modified by changing the espressomd.System() properties. We create an instance system and set the box length, periodicity and time step. The skin depth skin is a parameter for the link--cell system which tunes its performance, but shall not be discussed here.

In :
# Setup System
box_l = (n_part / density)**(1. / 3.)
system = System(box_l=[box_l, box_l, box_l])
system.periodicity = [True, True, True]
system.time_step = time_step
system.cell_system.skin = 0.3


We now fill this simulation box with particles at random positions, using type and charge from our dictionaries. Using the length of the particle list system.part for the id, we make sure that our particles are numbered consecutively. The particle type is used to link non-bonded interactions to a certain group of particles.

In :
for i in range(int(n_ionpairs)):
pos=numpy.random.random(3) * box_l, q=charges["Anion"])
for i in range(int(n_ionpairs)):
pos=numpy.random.random(3) * box_l, q=charges["Cation"])


Before we can really start the simulation, we have to specify the interactions between our particles. We already defined the WCA parameters at the beginning, what is left is to specify the combination rule and to iterate over particle type pairs. For simplicity, we implement only the Lorentz-Berthelot rules. We pass our interaction pair to system.non_bonded_inter[\*,\*] and set the pre-calculated WCA parameters epsilon and sigma.

In :
def combination_rule_epsilon(rule, eps1, eps2):
if rule == "Lorentz":
return (eps1 * eps2)**0.5
else:
return ValueError("No combination rule defined")

def combination_rule_sigma(rule, sig1, sig2):
if rule == "Berthelot":
return (sig1 + sig2) * 0.5
else:
return ValueError("No combination rule defined")

# WCA interactions parameters
for s in [["Anion", "Cation"], ["Anion", "Anion"], ["Cation", "Cation"]]:
wca_sig = combination_rule_sigma("Berthelot", wca_sigmas[s], wca_sigmas[s])
wca_eps = combination_rule_epsilon("Lorentz", wca_epsilons[s], wca_epsilons[s])

system.non_bonded_inter[types[s], types[s]].wca.set_params(
epsilon=wca_eps, sigma=wca_sig)


## 3 Equilibration¶

With randomly positioned particles, we most likely have huge overlap and the strong repulsion will cause the simulation to crash. The next step in our script therefore is a suitable WCA equilibration. This is known to be a tricky part of a simulation and several approaches exist to reduce the particle overlap. Here, we use the steepest descent algorithm and cap the maximal particle displacement per integration step to 1% of the smallest sigma. We use system.analysis.min_dist() to get the minimal distance between all particles pairs. This value is used to stop the minimization when particles are far away enough from each other. At the end, we activate the Langevin thermostat for our NVT ensemble with temperature temp and friction coefficient gamma.

In :
# WCA Equilibration
max_sigma = max(wca_sigmas.values())
min_sigma = min(wca_sigmas.values())
min_dist = 0.0

system.integrator.set_steepest_descent(f_max=0, gamma=10.0,
max_displacement=min_sigma * 0.01)
while min_dist < max_sigma:
min_dist = system.analysis.min_dist([types["Anion"], types["Cation"]],
[types["Anion"], types["Cation"]])
system.integrator.run(10)
system.integrator.set_vv()

# Set thermostat
system.thermostat.set_langevin(kT=temp, gamma=gamma, seed=42)


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 with parameters 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:

In :
p3m = electrostatics.P3M(prefactor=l_bjerrum * temp,
accuracy=1e-3)

P3M tune parameters: Accuracy goal = 1.00000e-03 prefactor = 7.00000e+00
System: box_l = 7.36806e+00 # charged part = 200 Sum[q_i^2] = 2.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   4.59284e-01 6.50911e+00 8.77128e-02 7.071e-04 8.771e-02 accuracy not achieved
8    7   4.59284e-01 6.50911e+00 8.77128e-02 7.071e-04 8.771e-02 accuracy not achieved
8    7   4.59284e-01 6.50911e+00 8.77128e-02 7.071e-04 8.771e-02 accuracy not achieved
10   7   4.59284e-01 6.50911e+00 1.39625e-02 7.071e-04 1.394e-02 accuracy not achieved
10   7   4.59284e-01 6.50911e+00 1.39625e-02 7.071e-04 1.394e-02 accuracy not achieved
10   7   4.59284e-01 6.50911e+00 1.39625e-02 7.071e-04 1.394e-02 accuracy not achieved
12   7   4.59284e-01 6.50911e+00 2.51196e-03 7.071e-04 2.410e-03 accuracy not achieved
12   7   4.59284e-01 6.50911e+00 2.51196e-03 7.071e-04 2.410e-03 accuracy not achieved
12   7   4.59284e-01 6.50911e+00 2.51196e-03 7.071e-04 2.410e-03 accuracy not achieved
14   7   4.44034e-01 6.73901e+00 9.96813e-04 7.071e-04 7.026e-04 1.12
14   6   4.59284e-01 6.50911e+00 1.30343e-03 7.071e-04 1.095e-03 accuracy not achieved
14   7   4.44034e-01 6.73901e+00 9.96813e-04 7.071e-04 7.026e-04 1.18
14   6   4.44034e-01 6.73901e+00 1.65807e-03 7.071e-04 1.500e-03 accuracy not achieved
16   7   3.92866e-01 7.64271e+00 9.88469e-04 7.071e-04 6.907e-04 1.08
16   6   4.26689e-01 7.02075e+00 9.95097e-04 7.071e-04 7.002e-04 1.14
16   5   4.44034e-01 6.73901e+00 1.71326e-03 7.071e-04 1.561e-03 accuracy not achieved
16   7   3.92099e-01 7.65808e+00 9.99269e-04 7.071e-04 7.061e-04 1.09
16   6   3.92866e-01 7.64271e+00 1.65362e-03 7.071e-04 1.495e-03 accuracy not achieved
16   7   3.92099e-01 7.65808e+00 9.99269e-04 7.071e-04 7.061e-04 1.08
16   6   3.92099e-01 7.65808e+00 1.67861e-03 7.071e-04 1.522e-03 accuracy not achieved
18   7   3.51510e-01 8.56819e+00 9.98800e-04 7.071e-04 7.054e-04 1.25
18   6   3.82909e-01 7.84703e+00 9.95996e-04 7.071e-04 7.014e-04 1.29
18   5   3.92099e-01 7.65808e+00 1.91148e-03 7.071e-04 1.776e-03 accuracy not achieved
18   7   3.51510e-01 8.56819e+00 9.98800e-04 7.071e-04 7.054e-04 1.25
18   6   3.51510e-01 8.56819e+00 1.69109e-03 7.071e-04 1.536e-03 accuracy not achieved
18   7   3.51510e-01 8.56819e+00 9.98800e-04 7.071e-04 7.054e-04 1.25
18   6   3.51510e-01 8.56819e+00 1.69109e-03 7.071e-04 1.536e-03 accuracy not achieved
20   7   3.19243e-01 9.45923e+00 9.90212e-04 7.071e-04 6.932e-04 1.28
20   6   3.47391e-01 8.67261e+00 9.98911e-04 7.071e-04 7.056e-04 1.30
20   5   3.51510e-01 8.56819e+00 2.08854e-03 7.071e-04 1.965e-03 accuracy not achieved
20   7   3.19243e-01 9.45923e+00 9.90212e-04 7.071e-04 6.932e-04 1.29
20   6   3.19243e-01 9.45923e+00 1.68314e-03 7.071e-04 1.527e-03 accuracy not achieved
22   7   2.91808e-01 1.03740e+01 9.97604e-04 7.071e-04 7.037e-04 1.69
22   6   3.18619e-01 9.47825e+00 9.94623e-04 7.071e-04 6.995e-04 1.72
22   5   3.19243e-01 9.45923e+00 2.22942e-03 7.071e-04 2.114e-03 accuracy not achieved
22   7   2.91808e-01 1.03740e+01 9.97604e-04 7.071e-04 7.037e-04 1.71
22   6   2.91808e-01 1.03740e+01 1.71200e-03 7.071e-04 1.559e-03 accuracy not achieved
22   7   2.91808e-01 1.03740e+01 9.97604e-04 7.071e-04 7.037e-04 1.70
22   6   2.91808e-01 1.03740e+01 1.71200e-03 7.071e-04 1.559e-03 accuracy not achieved
24   7   2.69580e-01 1.12536e+01 9.88918e-04 7.071e-04 6.913e-04 1.83
24   6   2.91808e-01 1.03740e+01 1.03216e-03 7.071e-04 7.519e-04 accuracy not achieved
24   7   2.69054e-01 1.12762e+01 9.99666e-04 7.071e-04 7.066e-04 1.83
24   6   2.69580e-01 1.12536e+01 1.70164e-03 7.071e-04 1.548e-03 accuracy not achieved
24   7   2.69054e-01 1.12762e+01 9.99666e-04 7.071e-04 7.066e-04 1.84
24   6   2.69054e-01 1.12762e+01 1.72746e-03 7.071e-04 1.576e-03 accuracy not achieved
26   7   2.50136e-01 1.21531e+01 9.91724e-04 7.071e-04 6.954e-04 2.50
26   6   2.69054e-01 1.12762e+01 1.07209e-03 7.071e-04 8.058e-04 accuracy not achieved
26   7   2.50136e-01 1.21531e+01 9.91724e-04 7.071e-04 6.954e-04 2.50
26   6   2.50136e-01 1.21531e+01 1.71807e-03 7.071e-04 1.566e-03 accuracy not achieved
26   7   2.50136e-01 1.21531e+01 9.91724e-04 7.071e-04 6.954e-04 2.50
26   6   2.50136e-01 1.21531e+01 1.71807e-03 7.071e-04 1.566e-03 accuracy not achieved
28   7   2.33525e-01 1.30417e+01 9.91061e-04 7.071e-04 6.944e-04 2.72
28   6   2.50136e-01 1.21531e+01 1.09949e-03 7.071e-04 8.419e-04 accuracy not achieved
28   7   2.33525e-01 1.30417e+01 9.91061e-04 7.071e-04 6.944e-04 2.72
28   6   2.33525e-01 1.30417e+01 1.72548e-03 7.071e-04 1.574e-03 accuracy not achieved
30   7   2.18930e-01 1.39353e+01 9.93493e-04 7.071e-04 6.979e-04 3.41
30   6   2.33525e-01 1.30417e+01 1.13156e-03 7.071e-04 8.834e-04 accuracy not achieved

resulting parameters: mesh: (16 16 16), cao: 7, r_cut_iL: 3.9287e-01,
alpha_L: 7.6427e+00, accuracy: 9.8847e-04, time: 1.08



Before the production part of the simulation, we do a quick temperature equilibration. For the output, we gather all energies with system.analysis.energy(), calculate the "current" temperature from the ideal part and print it to the screen along with the total and Coulomb energies. Note that for the ideal gas the temperature is given via $1/2 m \sqrt{\langle v^2 \rangle}=3/2 k_BT$, where $\langle \cdot \rangle$ denotes the ensemble average. Calculating some kind of "current temperature" via $T_\text{cur}=\frac{m}{3 k_B} \sqrt{ v^2 }$ won't produce the temperature in the system. Only when averaging the squared velocities first one would obtain the temperature for the ideal gas. $T$ is a fixed quantity and does not fluctuate in the canonical ensemble.

We integrate for a certain amount of steps with system.integrator.run(100).

In :
# Temperature Equilibration
system.time = 0.0
for i in range(int(num_steps_equilibration / 50)):
energy = system.analysis.energy()
temp_measured = energy['kinetic'] / ((3.0 / 2.0) * n_part)
print("t={0:.1f}, E_total={1:.2f}, E_coulomb={2:.2f},T={3:.4f}"
.format(system.time, energy['total'], energy['coulomb'], temp_measured), end='\r')
system.integrator.run(200)

print()

t=6.0, E_total=-471.27, E_coulomb=-849.76,T=1.0584


## 4 Running the Simulation¶

Now we can integrate the particle trajectories for a couple of time steps. Our integration loop basically looks like the equilibration, but we also want to calculate the averaged radial distribution functions $g_{++}(r)$ and $g_{+-}(r)$ using the RDF observable:

In :
from espressomd.observables import RDF
from espressomd.accumulators import MeanVarianceCalculator
import numpy as np

# Calculate the averaged rdfs
rdf_bins = 100
r_min = 0.0
r_max = system.box_l / 2.0
pids_anion = system.part.select(type=types["Anion"]).id
pids_cation = system.part.select(type=types["Cation"]).id
rdf_00_obs = RDF(ids1=pids_anion, ids2=pids_anion, min_r=r_min, max_r=r_max, n_r_bins=rdf_bins)
rdf_01_obs = RDF(ids1=pids_anion, ids2=pids_cation, min_r=r_min, max_r=r_max, n_r_bins=rdf_bins)
rdf_00_acc = MeanVarianceCalculator(obs=rdf_00_obs, delta_N=integ_steps_per_config)
rdf_01_acc = MeanVarianceCalculator(obs=rdf_01_obs, delta_N=integ_steps_per_config)

# Integration
system.time = 0.0
for i in range(num_configs):
energy = system.analysis.energy()
temp_measured = energy['kinetic'] / ((3.0 / 2.0) * n_part)
print("progress={:.0f}%, t={:.1f}, E_total={:.2f}, E_coulomb={:.2f}, T={:.4f}"
.format((i + 1) * 100. / num_configs, system.time, energy['total'],
energy['coulomb'], temp_measured), end='\r')
system.integrator.run(integ_steps_per_config)

print()

progress=100%, t=247.5, E_total=-480.94, E_coulomb=-856.97, T=0.9963


The shown RDF commands return the radial distribution functions for equally and oppositely charged particles for specified radii and number of bins. In this case, we calculate the averaged rdf using the MeanVarianceCalculator accumulator.

## 5 Analysis¶

The results are two NumPy arrays containing the $g(r)$ values.

In :
rdf_00 = rdf_00_acc.get_mean()
rdf_01 = rdf_01_acc.get_mean()
r = rdf_00_obs.bin_centers()


We can now write the $g_{++}(r)$ and $g_{+-}(r)$ data into a file with standard python output routines.

In :
with open('rdf.data', 'w') as rdf_fp:
for i in range(rdf_bins):
rdf_fp.write("%1.5e %1.5e %1.5e\n" % (r[i], rdf_00[i], rdf_01[i]))


Finally we can plot the two radial distribution functions using pyplot.

In :
# Plot the distribution functions
plt.figure(figsize=(10, 6), dpi=80)
plt.plot(r[:], rdf_00[:], label='Na$-$Na')
plt.plot(r[:], rdf_01[:], label='Na$-$Cl')
plt.xlabel('$r$', fontsize=20)
plt.ylabel('$g(r)$', fontsize=20)
plt.legend(fontsize=20)
plt.show() ## 6 Task - Real Units¶

So far, the system has arbitrary units and is not connected to any real physical system. Simulate a proper NaCl crystal with the force field parameters taken from :

Ion $q/\mathrm{e}$ $\sigma/\mathrm{\mathring{A}}$ $(\epsilon/\mathrm{k_B})/\mathrm{K}$ $m/\mathrm{u}$
Na +1 2.52 17.44 22.99
Cl -1 3.85 192.45 35.453

Use the following system parameters:

Parameter Value
Temperature $298\ \mathrm{K}$
Fiction Coeff. $10\ \mathrm{ps}^{-1}$
Density $1.5736\ \mathrm{u \mathring{A}}^{-3}$
Bjerrum Length (298 K) $439.2\ \mathrm{\mathring{A}}$
Time Step $2\ \mathrm{fs}$

To make your life easier, don't try to equilibrate randomly positioned particles, but set them up in a crystal structure close to equilibrium. If you do it right, you don't even need the WCA equilibration. To speed things up, don't go further than 1000 particles and use a P$^3$M accuracy of $10^{-2}$. Your RDF should look like the plot in figure 2. If you get stuck, you can look at the solution script /doc/tutorials/02-charged_system/scripts/nacl_units.py (or nacl_units_vis.py with visualization).

 R. Fuentes-Azcatl and M. Barbosa, Sodium Chloride, NaCl/$\epsilon$ : New Force Field, J. Phys. Chem. B, 2016, 120(9), pp 2460-2470.