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 MASS
#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:

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.

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.

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.

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 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.

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:

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).

Figure 1: VMD Snapshot of the Salt System

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:

Additionally, we append all particle configurations in the core with system.analysis.append() for a very convenient analysis later on.

5 Analysis

Now, we want to calculate the averaged radial distribution functions $g_{++}(r)$ and $g_{+-}(r)$ with the rdf() command from system.analysis:

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 of the stored configurations, denoted by the chevrons in rdf_type='$<\mathrm{rdf}>$'. Using rdf_type='rdf' would simply calculate the rdf of the current particle configuration. The results are two NumPy arrays containing the $r$ and $g(r)$ values. We can then write the data into a file with standard python output routines.

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

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 parameter taken from [1]:

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).

Figure 2: Snapshot and RDF of the parameterized NaCl crystal.


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