# 1. Introduction¶

ESPResSo is a simulation package designed to perform Molecular Dynamics (MD) and Monte Carlo (MC) simulations. It is meant to be a universal tool for simulations of a variety of soft matter systems. It features a broad range of interaction potentials which opens up possibilities for performing simulations using models with different levels of coarse-graining. It also includes modern and efficient algorithms for treatment of Electrostatics (P3M, MMM-type algorithms, constant potential simulations, dielectric interfaces, …), hydrodynamic interactions (DPD, Lattice-Boltzmann), and magnetic interactions, only to name a few. It is designed to exploit the capabilities of parallel computational environments. The program is being continuously extended to keep the pace with current developments both in the algorithms and software.

The kernel of ESPResSo is written in C++ with computational efficiency in mind. Interaction between the user and the simulation engine is provided via a Python scripting interface. This enables setup of arbitrarily complex systems which users might want to simulate in future, as well as modifying simulation parameters during runtime.

## 1.1. Guiding principles¶

ESPResSo is a tool for performing computer simulation and this user guide describes how to use this tool. However, it should be borne in mind that being able to operate a tool is not sufficient to obtain physically meaningful results. It is always the responsibility of the user to understand the principles behind the model, simulation and analysis methods he or she is using.

It is expected that the users of ESPResSo and readers of this user guide have a thorough understanding of simulation methods and algorithms they are planning to use. They should have passed a basic course on molecular simulations or read one of the renown textbooks, [FS02]. It is not necessary to understand everything that is contained in ESPResSo, but it is inevitable to understand all methods that you want to use. Using the program as a black box without proper understanding of the background will most probably result in wasted user and computer time with no useful output.

To enable future extensions, the functionality of the program is kept as general as possible. It is modularized, so that extensions to some parts of the program (implementing a new potential) can be done by modifying or adding only few files, leaving most of the code untouched.

To facilitate the understanding and the extensibility, much emphasis is put on readability of the code. Hard-coded assembler loops are generally avoided in hope that the overhead in computer time will be more than compensated for by saving much of the user time while trying to understand what the code is supposed to do.

Hand-in-hand with the extensibility and readability of the code comes the flexibility of the whole program. On the one hand, it is provided by the generalized functionality of its parts, avoiding highly specialized functions. An example can be the implementation of the Generic Lennard-Jones potential described in section Generic Lennard-Jones interaction where the user can change all available parameters. Where possible, default values are avoided, providing the user with the possibility of choice. ESPResSo cannot be aware whether your particles are representing atoms or billiard balls, so it cannot check if the chosen parameters make sense and it is the user’s responsibility to make sure they do. In fact, ESPResSo can be used to play billiard (see the sample script in samples/billiard.py)!

On the other hand, flexibility of ESPResSo stems from the employment of a scripting language at the steering level. Apart from the ability to modify the simulation and system parameters at runtime, many simple tasks which are not computationally critical can be implemented at this level, without even touching the C++-kernel. For example, simple problem-specific analysis routines can be implemented in this way and made interact with the simulation core. Another example of the program’s flexibility is the possibility to integrate system setup, simulation and analysis in one single control script. ESPResSo provides commands to create particles and set up interactions between them. Capping of forces helps prevent system blow-up when initially some particles are placed on top of each other. Using the Python interface, one can simulate the randomly set-up system with capped forces, interactively check whether it is safe to remove the cap and switch on the full interactions and then perform the actual productive simulation.

## 1.2. Basic program structure¶

As already mentioned, ESPResSo consists of two components. The simulation engine is written in C and C++ for the sake of computational efficiency. The steering or control level is interfaced to the kernel via an interpreter of Python scripting languages.

The kernel performs all computationally demanding tasks. Before all, integration of Newton’s equations of motion, including calculation of energies and forces. It also takes care of internal organization of data, storing the data about particles, communication between different processors or cells of the cell-system. The kernel is modularized so that basic functions are accessed via a set of well-defined lean interfaces, hiding the details of the complex numerical algorithms.

The scripting interface (Python) is used to setup the system (particles, boundary conditions, interactions, …), control the simulation, run analysis, and store and load results. The user has at hand the full readability and functionality of the scripting language. For instance, it is possible to use the SciPy package for analysis and PyPlot for plotting. With a certain overhead in efficiency, it can also be used to reject/accept new configurations in combined MD/MC schemes. In principle, any parameter which is accessible from the scripting level can be changed at any moment of runtime. In this way methods like thermodynamic integration become readily accessible.

The focus of the user guide is documenting the scripting interface, its behavior and use in the simulation. It only describes certain technical details of implementation which are necessary for understanding how the script interface works. Technical documentation of the code and program structure is contained in the online wiki.

## 1.3. Basic python simulation script¶

In this section, a brief overview is given over the most important components of the Python interface. Their usage is illustrated by short examples, which can be put together to a demo script.

Imports

As usual, the Python script starts by importing the necessary modules. The ESPResSo interface is contained in the espressomd Python module, which needs to be imported, before anything related can be done.

import espressomd


This should be followed by further necessary imports of the example at hand:

from espressomd.interactions import HarmonicBond
from espressomd.electrostatics import P3M


espressomd.System

Access to the simulation system is provided via the System class. As a first step, an instance of this class needs to be created.

system = espressomd.System(box_l=[10, 10, 10])


Note that only one instance of the System class can be created due to limitations in the simulation core. Properties of the System class are used to access the parameters concerning the simulation system such as box geometry, time step or cell-system:

print("The box dimensions are {}".format(system.box_l))
system.time_step = 0.01
system.cell_system.skin = 0.4


Particles

The particles in the simulation are accessed via system.part, an instance of the ParticleList class. Use the add method to create new particles:

system.part.add(id=0, pos=[1.0, 1.0, 1.0], type=0)


Individual particles can be retrieved by their numerical id using angular brackets:

system.part[1].pos = [1.0, 1.0, 2.0]


It is also possible to loop over all particles:

for p in system.part:
print("Particle id {}, type {}".format(p.id, p.type))


An individual particle is represented by an instance of ParticleHandle. The properties of the particle are implemented as Python properties.

particle = system.part[0]
particle.type = 0
print("Position of particle 0: {}".format(particle.pos))


Properties of several particles can be accessed by using Python slices:

positions = system.part[:].pos


Interactions

In ESPResSo, interactions between particles usually fall in three categories:

• Non-bonded interactions are short-ranged interactions between all pairs of particles of specified types. An example is the Lennard-Jones interaction mimicking overlap repulsion and van-der-Waals attraction.

• Bonded interactions act only between two specific particles. An example is the harmonic bond between adjacent particles in a polymer chain.

• Long-range interactions act between all particles with specific properties in the entire system. An example is the Coulomb interaction.

Non-bonded interaction

Non-bonded interactions are represented as subclasses of NonBondedInteraction, e.g. LennardJonesInteraction. Instances of these classes for a given pair of particle types are accessed via the non_bonded_inter attribute of the System class. This sets up a Lennard-Jones interaction between all particles of type 0 with the given parameters:

system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=1, sigma=1, cutoff=5.0, shift="auto")


Bonded interaction

Next, we add another pair of particles with a different type to later add a harmonic bond between them:

system.part.add(id=2, pos=[7.0, 7.0, 7.0], type=1)


To set up a bonded interaction, first an instance of the appropriate class is created with the desired parameters:

harmonic = HarmonicBond(k=1.0, r_0=0.5)


Then, the bonded interaction is registered in the simulation core by adding the instance to bonded_inter:

system.bonded_inter.add(harmonic)


Finally, the bond can be added to particles using the add_bond() method of ParticleHandle with the instance of the bond class and the id of the bond partner particle:

system.part[2].add_bond((harmonic, 3))


Charges

Now we want to setup a pair of charged particles treated by the P3M electrostatics solver. We start by adding the particles:

system.part.add(id=4, pos=[4.0, 1.0, 1.0], type=2, q=1.0)
system.part.add(id=5, pos=[6.0, 1.0, 1.0], type=2, q=-1.0)


Long-range interactions and other methods that might be mutually exclusive are treated as so-called actors. They are used by first creating an instance of the desired actor:

p3m = P3M(accuracy=1e-3, prefactor=1.0)


and then adding it to the system:

print("Tuning p3m...")


Integration

So far we just added particles and interactions, but did not propagate the system. This is done by the integrator. It uses by default the velocity Verlet algorithm and is already created by the system class. To perform an integration step, just execute:

system.integrator.run(1)


Usually, the system is propagated for a number of steps in a loop alongside with some analysis. In this last snippet, the different energy contributions of the system are printed:

num_configs = 10
num_steps = 1000

for i in range(num_configs):
system.integrator.run(num_steps)
energy = system.analysis.energy()
print("System time: {}".format(system.time))
print("Energy of the LJ interaction: {}".format(energy["non_bonded"]))
print("Energy of the harmonic bond: {}".format(energy["bonded"]))
print("Energy of the Coulomb interaction: {}".format(energy["coulomb"]))


## 1.4. Tutorials¶

There are a number of tutorials that introduce the use of ESPResSo for different physical systems. You can also find the tutorials and related scripts in the directory /doc/tutorials or online on GitHub. Currently, the following tutorials are available:

• 01-lennard_jones: Modelling of a single-component and a two-component Lennard-Jones liquid.

• 02-charged_system: Modelling of charged systems such as ionic crystals.

• 04-lattice_boltzmann: Simulations including hydrodynamic interactions using the lattice-Boltzmann method.

• 05-raspberry_electrophoresis: Extended objects in a lattice-Boltzmann fluid, raspberry particles.

• 06-active_matter: Modelling of self-propelling particles.

• 07-electrokinetics: Modelling electrokinetics together with hydrodynamic interactions.

• 08-visualization: Using the online visualizers of ESPResSo.

• 11-ferrofluid: Modelling a colloidal suspension of magnetic particles.

• 12-constant_pH: Modelling the titration of a weak acid using the constant pH method

## 1.5. Sample scripts¶

Several scripts that can serve as usage examples can be found in the directory /samples, or in the git repository.

• billiard.py

ESPResSo 8Ball billiards game.

• chamber_game.py

Game based on Maxwell’s demon, a thought experiment used to teach statistical thermodynamics. The user has to scoop particles from a chamber and guide them to another chamber through a channel with the help of a snake controlled by a gamepad or the keyboard. The particle imbalance between chambers creates a pressure gradient that makes it harder to move particles to the chamber with an excess of particles.

• constraints.py

Confine a polymer between two slabs and check that it cannot escape them during the entire simulation.

• diffusion_coefficient.py

Compare the diffusion coefficient of a single thermalized particle obtained from the particle’s mean square displacement and the auto correlation function of its velocity to the expected value. Uses the Observables/Correlators framework.

• dpd.py

Set up a DPD fluid and calculate pressure as a function of the varying density. The fluid is thermalized using a DPD thermostat.

• drude_bmimpf6.py

Particle polarization with cold Drude oscillators on a coarse-grained simulation of the ionic liquid BMIM PF6.

• ekboundaries.py

Set up an electrokinetics (LB) fluid confined between charged walls.

• electrophoresis.py

Simulate electrophoresis of a linear polymer using the P3M electrostatics solver.

• espresso_logo.py

Build the ESPResSo logo with particles.

• grand_canonical.py

Perform a grand canonical simulation of a system in contact with a salt reservoir while maintaining a constant chemical potential.

• h5md.py

Write ESPResSo trajectories in the H5MD format. See Writing H5MD-files.

• immersed_boundary/sampleImmersedBoundary.py

Simulate the motion of a spherical red blood cell-like particle advected in a planar Poiseuille flow, with or without volume conservation. For more details, see Immersed Boundary Method for soft elastic objects.

• lb_profile.py

Simulate the flow of a lattice-Boltzmann fluid past a cylinder, obtain the velocity profile in polar coordinates and compare it to the analytical solution.

• lbf.py

Set up a lattice-Boltzmann fluid and apply an external force density on it.

• lj-demo.py

Simulate a Lennard-Jones fluid in different thermodynamic ensembles (NVT, NpT). Sliders from a MIDI controller can change system variables such as temperature and volume. Some thermodynamic observables are analyzed and plotted live.

• lj_liquid.py

Simulate a Lennard-Jones fluid maintained at a fixed temperature by a Langevin thermostat. Shows the basic features of how to:

• set up system parameters, particles and interactions.

• warm up and integrate.

• write parameters, configurations and observables to files.

The particles in the system are of two types: type 0 and type 1. Type 0 particles interact with each other via a repulsive WCA interaction. Type 1 particles neither interact with themselves nor with type 0 particles.

• lj_liquid_distribution.py

Set up a Lennard-Jones fluid maintained at a fixed temperature by a Langevin thermostat. The particles in the system are of two types: type 0 and type 1. Type 0 particles interact with each other via a repulsive WCA interaction. Type 1 particles neither interact with themselves nor with type 0 particles. The distribution of minimum distances between particles of type 0 and type 1 is recorded with distribution(). See Particle distribution.

• lj_liquid_structurefactor.py

Set up a Lennard-Jones fluid maintained at a fixed temperature by a Langevin thermostat. The particles in the system are of two types: type 0 and type 1. Type 0 particles interact with each other via a repulsive WCA interaction. Type 1 particles neither interact with themselves nor with type 0 particles. The spherically averaged structure factor of particles of type 0 and type 1 is calculated with structure_factor(). See Structure factor.

• load_checkpoint.py

Basic usage of the checkpointing feature. Show how to load the state of:

• custom user variables.

• non-bonded interactions.

• particles.

• P3M parameters.

• thermostat.

• MDAnalysisIntegration.py

Show how to expose configuration to MDAnalysis at run time. The functions of MDAnalysis can be used to perform some analysis or convert the frame to other formats (CHARMM, GROMACS, …). For more details, see Writing various formats using MDAnalysis.

• minimal-charged-particles.py

Simulate an equal number of positively and negatively charged particles using the P3M solver. The system is maintained at a constant temperature using a Langevin thermostat.

• minimal-diamond.py

Set up a diamond-structured polymer network.

• minimal-polymer.py

Set up a linear polymer.

• object_in_fluid/motivation.py

Simulate the motion of flexible red blood cells in a lattice-Boltzmann fluid with solid obstacles. For more details, see Object-in-fluid.

• observables_correlators.py

Measure mean square displacements using the Observables/Correlators framework.

• p3m.py

Simulate a Lennard-Jones liquid with charges. The P3M method is used to calculate electrostatic interactions. The ELC method can be optionally added to subtract the electrostatic contribution from the z-direction. For more details, see Electrostatic Layer Correction (ELC).

• reaction_ensemble.py

Guide for the reaction ensemble and the constant pH ensemble.

• rigid_body.py

Demonstrates the construction of a rigid object by means of the VIRTUAL_SITES_RELATIVE feature.

• save_checkpoint.py

Basic usage of the checkpointing feature. Show how to write the state of:

• custom user variables.

• non-bonded interactions.

• particles.

• P3M parameters.

• thermostat.

• slice_input.py

Illustrate how particles of interest can be accessed via slicing.

• visualization_bonded.py

Visualize the simulation of a linear polymer.

• visualization_cellsystem.py

Visualize the system cells and MPI domains. Run ESPResSo in parallel to color particles by node.

• visualization_charged.py

Visualize a simulation with a pool of particles with various charges, LJ parameters and masses.

• visualization_constraints.py

Visualize shape-based constraints interacting with a Lennard-Jones gas.

• visualization_interactive.py

Visualize a simulation box where the particles can be repositioned via the mouse and timed callbacks, and the temperature of the thermostat can be changed via the keyboard.

• visualization_lbboundaries.py

Visualize lattice-Boltzmann boundary nodes.

• visualization_ljliquid.py

Visualize a Lennard-Jones liquid with live plotting via matplotlib.

• visualization_mmm2d.py

Visualize charged particles confined between two plates of a capacitor with an applied potential difference.

• visualization_npt.py

Visualize particle dumbbells in the NpT ensemble (constant temperature, constant pressure, variable volume).

• visualization_poiseuille.py

Visualize the Poiseuille flow in a lattice-Boltzmann fluid with an external force applied.

• wang_landau_reaction_ensemble.py

Simulate two reacting monomers which are bonded via a harmonic potential. The script aborts as soon as the abortion criterion in the Wang-Landau algorithm is met: the Wang-Landau simulation runs until the Wang-Landau potential has converged and then raises a Warning that it has converged, effectively aborting the simulation.

With the setup of the Wang-Landau algorithm in this script, you sample the density of states of a three-dimensional reacting harmonic oscillator as a function of the two collective variables 1) degree of association and 2) potential energy.

The recorded Wang-Landau potential (which is updated during the simulation) is written to the file WL_potential_out.dat.

In this simulation setup the Wang-Landau potential is the density of states. You can view the converged Wang-Landau potential e.g. via plotting with gnuplot: splot "WL_potential_out.dat". As expected the three-dimensional harmonic oscillator has a density of states which goes like $$\sqrt{E_{\text{pot}}}$$.

For a scientific description and different ways to use the algorithm please consult https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00791

• widom_insertion.py

Measure the excess chemical potential of a charged WCA fluid via Widom’s insertion method.

## 1.6. On units¶

What is probably one of the most confusing subjects for beginners of ESPResSo is, that ESPResSo does not predefine any units. While most MD programs specify a set of units, like, for example, that all lengths are measured in Ångström or nanometers, times are measured in nano- or picoseconds and energies are measured in $$\mathrm{kJ/mol}$$, ESPResSo does not do so.

Instead, the length-, time- and energy scales can be freely chosen by the user. Once these three scales are fixed, all remaining units are derived from these three basic choices.

The probably most important choice is the length scale. A length of $$1.0$$ can mean a nanometer, an Ångström, or a kilometer - depending on the physical system, that the user has in mind when he writes his ESPResSo-script. When creating particles that are intended to represent a specific type of atoms, one will probably use a length scale of Ångström. This would mean, that the parameter $$\sigma$$ of the Lennard-Jones interaction between two atoms would be set to twice the van-der-Waals radius of the atom in Ångström. Alternatively, one could set $$\sigma$$ to $$2.0$$ and measure all lengths in multiples of the van-der-Waals radius. When simulation colloidal particles, which are usually of micrometer size, one will choose their diameter (or radius) as basic length scale, which is much larger than the Ångström scale used in atomistic simulations.

The second choice to be made is the energy scale. One can for example choose to set the Lennard-Jones parameter $$\epsilon$$ to the energy in $$\mathrm{kJ/mol}$$. Then all energies will be measured in that unit. Alternatively, one can choose to set it to $$1.0$$ and measure everything in multiples of the van-der-Waals binding energy of the respective particles.

The final choice is the time (or mass) scale. By default, ESPResSo uses a reduced mass of 1 for all particles, so that the mass unit is simply the mass of one particle. Combined with the energy and length scale, this is sufficient to derive the resulting time scale:

$[\mathrm{time}] = [\mathrm{length}]\sqrt{\frac{[\mathrm{mass}]}{[\mathrm{energy}]}}$

This means, that if you measure lengths in Ångström, energies in $$k_B T$$ at 300K and masses in 39.95u, then your time scale is $$\mathring{A} \sqrt{39.95u / k_B T} = 0.40\,\mathrm{ps}$$.

On the other hand, if you want a particular time scale, then the mass scale can be derived from the time, energy and length scales as

$[\mathrm{mass}] = [\mathrm{energy}]\frac{[\mathrm{time}]^2}{[\mathrm{length}]^2}.$

By activating the feature MASS, you can specify particle masses in the chosen unit system.

A special note is due regarding the temperature, which is coupled to the energy scale by Boltzmann’s constant. However, since ESPResSo does not enforce a particular unit system, we also don’t know the numerical value of the Boltzmann constant in the current unit system. Therefore, when specifying the temperature of a thermostat, you actually do not define the temperature, but the value of the thermal energy $$k_B T$$ in the current unit system. For example, if you measure energy in units of $$\mathrm{kJ/mol}$$ and your real temperature should be 300K, then you need to set the thermostat’s effective temperature to $$k_B 300\, K \mathrm{mol / kJ} = 2.494$$.

As long as one remains within the same unit system throughout the whole ESPResSo-script, there should be no problems.

## 1.7. Available simulation methods¶

ESPResSo provides a number of useful methods. The following table shows the various methods as well as their status. The table distinguishes between the state of the development of a certain feature and the state of its use. We distinguish between five levels:

Core

means that the method is part of the core of ESPResSo, and that it is extensively developed and used by many people.

Good

means that the method is developed and used by independent people from different groups.

Group

means that the method is developed and used in one group.

Single

means that the method is developed and used by one person only.

None

means that the method is developed and used by nobody.

Experimental

means that the method might have side effects.

In the “Tested” column, we note whether there is an integration test for the method.

If you believe that the status of a certain method is wrong, please report so to the developers.

Feature

Development Status

Usage Status

Tested

Integrators, Thermostats, Barostats

Velocity-Verlet Integrator

Core

Core

Yes

Langevin Thermostat

Core

Core

Yes

Isotropic NPT

Experimental

None

Yes

Quaternion Integrator

Core

Good

Yes

Interactions

Short-range Interactions

Core

Core

Yes

Constraints

Core

Core

Yes

Relative Virtual Sites

Good

Good

Yes

RATTLE Rigid Bonds

Single

Group

Yes

Gay-Berne Interaction

Experimental

Experimental

No

Coulomb Interaction

P3M

Core

Core

Yes

P3M on GPU

Single

Single

Yes

Dipolar P3M

Group

Good

Yes

MMM1D

Single

Good

No

MMM2D

Group

Good

Yes

MMM1D on GPU

Single

Single

No

ELC

Good

Good

Yes

ICC*

Group

Group

Yes

Hydrodynamic Interaction

Lattice-Boltzmann

Core

Core

Yes

Lattice-Boltzmann on GPU

Group

Core

Yes

Input/Output

VTF output

Core

Core

Yes

VTK output

Group

Group

No

Checkpointing

Experimental

Experimental

Yes

Visualization

Online visualisation (Mayavi)

Good

Good

No

Online visualisation (OpenGL)

Good

Good

No

Miscellaneous

Electrokinetics

Group

Group

Yes

Collision Detection

Group

Group

Yes

Swimmer Reactions

Single

Single

Yes

Reaction Ensemble

Group

Group

Yes

Constant pH Method

Group

Group

Yes

Object-in-fluid

Group

Group

Yes

Immersed boundary method

Group

Group

Yes

DPD

Single

Good

Yes

## 1.8. Intended interface compatibility between ESPResSo versions¶

We use the following versioning scheme: major.minor.patch_level

With regards to the stability of the Python interface, we plan the following guidelines:

• patch_level: The Python interface will not change, if only the patch_level number is different. Example: 4.0.0 -> 4.0.1.

• minor: There will be no silent interface changes between two versions with different minor numbers. I.e., a simulation script will not silently produce different results with the new version. The interface can, however, be extended. In important cases, the interface can change in such a way that using the old interface produces a clear error message and the simulation is terminated. Example: 4.0.1 -> 4.1.0

• major: No guarantees are made for a transition between major versions. Example 4.1.2 -> 5.0.

• No guarantees are made with regards to the development branch on GitHub.

• No guarantees are made with respect to the C++ bindings in the simulation core.