ESPResSo 4.0 – An Extensible Software Package for Simulating Soft Matter Systems

# ESPResSo 4.0 – An Extensible Software Package for Simulating Soft Matter Systems

Florian Weik Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany    Rudolf Weeber Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany    Kai Szuttor Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany    Konrad Breitsprecher Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany    Joost de Graaf Institute for Theoretical Physics, Center for Extreme Matter and Emergent Phenomena, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands    Michael Kuron Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany    Jonas Landsgesell Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany    Henri Menke Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany Department of Physics, University of Otago, P.O. Box 56, Dunedin 9054, New Zealand    David Sean Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany    Christian Holm Institut für Computerphysik, Universität Stuttgart, Allmandring 3, 70569 Stuttgart, Germany
July 1, 2019
###### Abstract

ESPResSo is an extensible simulation package for research on soft matter. This versatile molecular dynamics program was originally developed for coarse-grained simulations of charged systems [Limbach et al., Comput. Phys. Commun. 174, 704 (2006)]. The scope of the software has since broadened considerably: ESPResSo can now be used to simulate systems with length scales spanning from the molecular to the colloidal. Examples include, self-propelled particles in active matter, membranes in biological systems, and the aggregation of soot particles in process engineering. ESPResSo also includes solvers for hydrodynamic and electrokinetic problems, both on the continuum and on the explicit particle level. Since our last description of version 3.1 [Arnold et al., Meshfree Methods for Partial Differential Equations VI, Lect. Notes Comput. Sci. Eng. 89, 1 (2013)], the software has undergone considerable restructuring. The biggest change is the replacement of the Tcl scripting interface with a much more powerful Python interface. In addition, many new simulation methods have been implemented. In this article, we highlight the changes and improvements made to the interface and code, as well as the new simulation techniques that enable a user of ESPResSo 4.0 to simulate physics that is at the forefront of soft matter research.

\pdfstringdefDisableCommands

## I Introduction

Soft matter de Gennes (1992); Doi (2013); Barrat and Hansen (2003) is an established field of research that concerns itself with systems ranging from the nanometer scale up to tens of micrometers, existing on the interface between physics, chemistry, and biology. These systems are classified as “soft” because the relevant energy scale of interactions is typically comparable to the thermal energy, making them easy to deform through the application of external forces and fields. The broadness of the field is exemplified by the subjects covered therein: polymer science Doi and Edwards (1988); Rubinstein and Colby (2003), colloidal science Verwey and Overbeek (1948); Löwen (2001), liquid crystals Chandrasekhar (1992); Hamley (2003), biological physics Nelson (2004); Levental et al. (2007), food science Ubbink et al. (2008), material science Polarz and Antonietti (2002), to name but a few.

Statistical techniques may be employed to describe soft matter systems, as thermal fluctuations and thus Brownian motion play an important role and allow the system to readily explore phase space. The field has thus proven itself to be a fruitful playground for shaping the community’s understanding of statistical physics concepts Kroy and Frey (2005); Frenkel (2002), as well as emergent non-equilibrium phenomena Seifert (2012); Cates (2012); Fodor and Marchetti (2018). The delicate interplay of energy and entropy that is found in many colloidal suspensions leads to a plethora of complex structural properties that are difficult to grasp with pen and paper. Indeed, the soft matter community owes a great deal to computer simulations for its current level of understanding.

The use of computers in soft matter traces its root to first molecular-dynamics (MD) simulations of phase separation in a system of hard spheres Alder and Wainwright (1957, 1962). From these humble beginnings, computer simulation grew out to a new field that bridges experimental and theoretical efforts. These methods came to their own in soft matter, as the typical size and nature of the interactions found therein, lend themselves well to coarse graining. This is a process by which many of the microscopic freedoms that govern a system may be captured in terms of effective interactions on a larger scale. The use of effective interactions enables the simulation of emergent effects using simple models at a fraction of the computation effort that is required to account for the full microscopic details. This has played a decisive role in elucidating, e.g., the guiding principles for phase behavior of complex liquids Pusey et al. (1994); Bates and Frenkel (1998); van Roij et al. (1999); Leunissen et al. (2005), the structure of ferrofluids Camp et al. (2000); Wang et al. (2002); Klapp and Schoen (2004); Klinkigt et al. (2013), the structure of polyelectrolytes Stevens and Kremer (1993); Dobrynin et al. (1996); Micka et al. (1999); Limbach et al. (2002), and the swelling of hydrogels and ferrogels Schneider and Linse (2002); Yan and de Pablo (2003); Mann et al. (2004); Weeber et al. (2012, 2018a).

The MD simulation package ESPResSo is one of the most broad, efficient, and well-suited platforms for the study of coarse-grained models of soft matter. Only LAMMPS Plimpton (1995) and ESPResSo++ Guzman et al. (2018) approach ESPResSo in terms of the scope of their simulation methods and features. Here, we should emphasize that while ESPResSo is capable of simulating molecular fluids, its scope does not cover atomisticly resolved systems. There are a number of better-suited, open-source softwares available for atomistic simulations, such as GROMACSBerendsen et al. (1995); van der Spoel et al. (2005); Pronk et al. (2013), NAMDNelson et al. (1996); Phillips et al. (2005, 2014), and AMBERPearlman et al. (1995); Case et al. (2005), to name just a few of the most commonly used ones.

The aim of this paper is to present the improvements made in ESPResSo 4.0 compared to our previous release ESPResSo 3.1 Arnold et al. (2013a). We have also taken the opportunity to introduce readers to our software package and encourage them to consider ESPResSo, both for use in soft matter research, and as a platform for method and algorithm development. A strong community is key to maintaining and extending a strong and reliable software for research. We will provide examples of cutting-edge physics that can be simulated using ESPResSo 4.0, which will hopefully create interest in tyring out ESPResSo and potentially joining our thriving user and developer community.

In brief, this paper covers the following subjects. We start by giving a broad introduction to ESPResSo, discussing its history, goals, and community in Section II. This will set the background to the developments made for our ESPResSo 4.0 release. Next, we focus on the new interface (Section III), the interaction with external packages (Section IV), user interaction with the software in terms of visualization and output (Sections V and VI, respectively), and changes to the new code development procedures and practices (Section VII) that we adopted with ESPResSo 4.0. Finally, we introduce the standout physical methods and models that were recently introduced within ESPResSo: energy minimization routines (Section IX), cluster analysis (Section X), chemical reactions (Section XI), polarizable molecules (Section XII), and self-propelled particles (Section XIII). We end with a summary and outlook in Section XIV.

## Ii Extensible Simulation Package for Research on Soft Matter Systems

### ii.1 The Idea behind ESPResSo

ESPResSo development began in 2001. Its name is an acronym for Extensible Simulation Package for Research on Soft matter systems, which encapsulates the goal of the software project:

On the one hand, the original development team wished to create a highly parallel molecular dynamics program that could deal efficiently with charged bead-spring models of polymers. At that time, we had learned to master electrostatic interactions in periodic geometries Deserno and Holm (1998a, b) and had also developed other fast electrostatic solvers for partially periodic geometries Arnold and Holm (2002a, b); Arnold et al. (2002a, b); Arnold and Holm (2005a, b); Arnold et al. (2006). The resulting software saw its first use in the study of polyelectrolyte solutions and charged colloidal suspensions.

On the other hand, the ESPResSo framework was designed to be “extensible”. That is, its core architecture was created to facilitate the testing and integration of new algorithms. ESPResSo has seen many additions since 2001, living up to this extensible character Arnold et al. (2013a). This includes support for the simulation of hydrodynamic interactions via an efficient GPU implementation of a lattice Boltzmann solver, rigid-body dynamics, a scheme for the study of irreversible agglomeration, and electrostatic solvers that allow for dielectric inclusions or even locally varying dielectric permittivities Tyagi et al. (2007, 2008, 2010); Arnold et al. (2013a, b); Fahrenberger and Holm (2014).

From the start, ESPResSo has a two-fold architecture, which roughly separates the developer from the user and promotes ease of use: (i) Massively parallel simulation routines were implemented in the C programming language using the Message Passing Interface (MPI). These routines scale well up to hundreds of processors on large computing clusters. They also facilitate the addition of new algorithms, even by novice developers. (ii) The simulation core exposes a scripting interface to the user, such that they did not have to directly interact with the core code. This interface was based on the Tcl language—the most suitable choice at that time to provide easy-to-use scripting. This gave the user programmatic control over the simulation protocol, and enabled run-time analysis and visualization in a straightforward manner.

### ii.2 The Motivation for ESPResSo 4.0

Every simulation package has to cope with changing computer architectures and programming models, ensure continued ease of development, and cater to users’ expectations concerning interaction with the software. ESPResSo 4.0 is our answer to the changes in the field of scientific computing that have taken place over the past five years.

Today, the Python language has become the de-facto standard in scientific computing. There exists a vast ecosystem of third-party Python packages for numerical algorithms, visualization, and statistical analysis. This motivated us to port the scripting interface from Tcl to Python for the latest version of ESPResSo, as will be described in Section III. We have also updated our online visualization to make the most of ESPResSo 4.0’s new Python interface, see Section V. Post-processing and visualization of large data sets is facilitated by the introduction of the HDF5 output format and parallellized output routines, see Section VI. Finally, the ESPResSo core has been converted to the C++ programming language, which has had added benefits both in terms of software development practices and in terms of performance, see Section VII.

To stay at the forefront of the field, it is not only necessary to commit to the user and development experience, it is also critical to keep up with the ever shifting interests of the user base. ESPResSo has done so in the past Arnold et al. (2013a) and continues this effort with the ESPResSo 4.0 release. Our simulation package now covers an even wider range of physical systems of relevance to modern soft matter research. This includes non-equilibrium phenomena in active matter, which ESPResSo 4.0 captures using model self-propelled particles Fischer et al. (2015); de Graaf et al. (2015, 2016a, 2016b), see Section XIII. Going beyond soft matter, particle-based methods have attracted interest from the engineering community. Specifically, this community utilizes MD to model large-scale physical processes such as soot aggregation Inci et al. (2014, 2017) and air filtration Schober et al. (2016), with the aim to improve these processes. Section X describes the opportunities for these kinds of studies in our latest release. However, as will become clear shortly, ESPResSo 4.0 is capable of much more.

### ii.3 The User and Developer Community

The ESPResSo package was set up as an open software project, committed to open source development under the GNU General Public License (GPL), a practice that we continue to adhere to. Users of ESPResSo are supported through mailing lists and workshops, such as the annual ESPResSo summer schools, which have between 20 and 40 participants. These workshops are often run as CECAM tutorials, jointly funded through grants from the German Research Foundation through the SFB 716 and the cluster of excellence SimTech (EXC 310). The mailing lists and workshops also allow the developer team to announce new features implemented in the code base of ESPResSo to the global research community.

ESPResSo is used by scientists all over the world, which is documented by more than 431 citations on Web of Science since 2006 and 608 citations on Google Scholar, as of this writing. Several research groups have contributed algorithms to the core of ESPResSo, which include: (i) Continuum flow solvers (including ones on GPUs) and various methods to couple them to MD Röhm and Arnold (2012); Cimrák et al. (2012); Cimrak et al. (2014); Cimrák et al. (2013); Rempfer et al. (2016); Guckenberger et al. (2016); Bächer et al. (2017). (ii) Advanced algorithms for rare-event sampling and the reconstruction of free-energy landscapes Kratzer et al. (2013, 2014). (iii) Monte-Carlo methods Samin et al. (2013). (iv) An interface to Python analysis package MDAnalysis. (v) Extensions to the support for anisotropic particles. (vi) Methods to calculate dipolar interactions on the GPU. Significant effort has been made to accommodate the various sources of these contributions and ensure code quality of ESPResSo 4.0. This includes putting in place an extended test infrastructure and increasing the test coverage. Furthermore, all changes to the software, including those made by the core team, undergo peer review before they are merged, see Section VII.

The latest version of ESPResSo and its documentation can be found on the website http://espressomd.org. Ongoing development is organized through the social code hosting platform GitHub, which facilitates contribution handling through its collaborative features.

## Iii The Python interface

In recent years, Python has become increasingly popular in the scientific community as a scripting language, in areas such as machine learning (PyTorch, scikit-learn, Keras) Paszke et al. (2017); Pedregosa et al. (2011); Chollet (2017), symbolic mathematics (SymPy, SAGE) Meurer et al. (2017); Stein and Joyner (2005), data visualization (Matplotlib, Mayavi) Hunter (2007); Ramachandran and Varoquaux (2011), and numerical mathematics (NumPy, SciPy, pandas) Jones et al. (year); McKinney (2010). In ESPResSo 4.0, we have therefore decided to switch the interface scripting language from Tcl to Python. This allows for the quick generation of clear and concise scripts that can harness the power of third party numerical tools. The ESPResSo core functionality can be accessed like any other Python module by importing the spressomd module. The whole architecture is built around the ystem object which encapsulates the physical state of the high-performance simulation core.

ESPResSo is a particle-based molecular dynamics simulator, i.e., all fundamental operations in the simulation setup revolve around particles. Inspired by Python’s list datatype and the NumPy package Jones et al. (year), we introduced slicing operations on the list of particles. Many functions also accept lists and automatically broadcast themselves to the individual list elements. This often reduces the number of lines of code significantly, e.g., placing several particles at random positions in the simulation box can now be done in a single line:

In Python, functions can be passed to other functions as arguments. ESPResSo 4.0 makes use of this, e.g., for particle selection using a user-specified criterion:

harged_particles = system.part.select(lambda p: p.q != 0.0)

Because the NumPy package Jones et al. (year) is very popular in numerical computing, including in undergraduate courses, we designed the interface for maximal compatibility. This allows the user to directly pass on values returned from ESPResSo functions to a variety of NumPy routines, e.g., to compute the mean, standard deviation, or histograms.

We further make use of the object-oriented programming capabilities of Python. Instead of commands having an immediate effect, the user instantiates objects of a class and adds these objects to the system instance. For example, when the user creates an object for the electrostatic solver, it is not immediately active, but has to be added to the list of actors first. This new design aims to reduce dependencies between different parts of the user-provided simulation script and eliminate issues related to the order of commands.

## Iv Integration with External packages

### iv.1 Using ESPResSo with other Python packages

As we already mentioned earlier, the ESPResSo Python interface was designed with NumPy in mind and most functions return datatypes which are compatible with the routines offered by NumPy. These routines are well-maintained and supported by a large community with numerous online resources and introductory material. One example are various statistics routines. These include mean, variance and median calculations, but also umpy.histogram(). As an example, a probability distribution of inter-particle distances can be used to obtain effective potentials via Boltzmann inversion Reith et al. (2003). Furthermore, the confidence interval estimators from cipy.stats.t.interval() can be used to control the amount of sampling in a simulation. Random number generation from umpy.random and cipy.stats can be used, e.g., for setting up a randomized starting configuration. Linear algebra methods from umpy.linalg allow for easily calculating the principal axes of and moments of inertia of a set of particles. Lastly, new interpolated bonded and pair potentials in ESPResSo can be added by providing NumPy arrays for the energies and forces.

It is also possible to create Jupyter notebooks Kluyver et al. (2016) using ESPResSo. These interactive worksheets are popular for teaching and provide a great way to nicely present results and the generating simulation code alongside. For example the introductory tutorials for ESPResSo 4.0 are presented in this fashion.

ESPResSo 4.0 provides an interface to the MDAnalysis Python package, which provides a large number of analysis routines for particle-based data. This includes the calculation of radial distribution functions, density profiles, and the persistence length of polymer chains.

### iv.2 Algorithms for Coulomb and dipolar interactions from the ScaFaCoS library

The treatment of Coulomb and dipolar interactions has always been one of the strong points of ESPResSo. Electrostatic interactions are inherently important in many soft matter systems: e.g. colloidal suspensions are often stabilized by means of electrostatic repulsion, the structure of polyelectrolytes depends on the concentration of charged salt atoms, and electric fields are used for the analysis and separation of particles in electrophoresis applications.

Being long-ranged, electrostatic interactions have to be treated by specialized algorithms in systems with periodic and mixed boundary conditions. ESPResSo already contains implementations of the particle-particle-particle-mesh (PM) algorithm Hockney and Eastwood (1988) for fully periodic systems, as well as the electrostatic layer correction (ELC) Arnold et al. (2002a) and MMM2D Arnold and Holm (2002a) scheme for 2D-periodic slit-pore geometries. However, more algorithms have been developed, designed for particular applications. Many of these are available in the ScaFaCoS library Arnold et al. (2013c) (Scalable Fast Coulomb Solvers). ESPResSo 4.0 provides an interface to this library, exposing all provided electrostatic solvers. Moreover, an extension to the ScaFaCoS library adds solvers for dipolar interactions with an scaling for periodic, mixed, and open boundary conditions. With ESPResSo 4.0, this version of the library can be used as a solver for magnetostatic interactions Nestler (2016); Weeber et al. (2018b). The dipolar PNFFT solver for open boundaries is of particular relevance to research in magnetic soft matter: properties of such material can depend on the sample shape due to the demagnetization field Raikher and Stolbov (2003); Morozov et al. (2009).

## V Online Visualization

A graphical representation of the system is of great value when it comes to presenting results of a simulation study. Also, it is common practice to visualize intermediate results while setting up the simulation. This saves time, as possible pitfalls related to positions, directions or orientations can be rectified immediately when discovered via visualization. Usually, this is done by writing out trajectories or volumetric data and utilizing external tools like Paraview Ayachit (2015) or VMD Humphrey et al. (1996). ESPResSo with Python is well suited for this kind of workflow. Furthermore, ESPResSo features two options for visualizing simulations while they are running, i.e., live or on-line visualization. The first one uses Mayavi Ramachandran and Varoquaux (2011), a Python framework for “3D scientific data visualization and plotting” to visualize particles. Mayavi has a user-friendly graphical interface to manipulate the visual appearance and to interact with the resulting representation. Second, ESPResSo 4.0 comes with a 3D rendering engine that uses PyOpenGL Fletcher and Liebscher (2005); Shreiner (1999). Both visualizers are not primarily intended to produce print-quality renderings, but rather to give a quick impression of the system setup and the equilibration process. Two exemplary snapshots are shown in Figs. 1 and 2; the corresponding Python scripts are available in the distribution package of ESPResSo 4.0.

Our OpenGL visualizer has been crafted specifically for use with ESPResSo and is capable of not only visualizing particles, but also several specific features like constraints, detailed particle properties, the cell system, the domain decomposition across processors, or fluid flows computed with the lattice Boltzmann method Röhm (2011) (Fig. 1). It has a large number of options to adjust colors, materials, lighting, camera perspective, and many more. Vectorial particle properties can be visualized by 3D arrows on the particles to display forces, velocities, etc. The OpenGL visualizer can also be used for offline visualization to display completed simulations.

A unique feature is the ability to interact with the running simulation: Particles can be queried for their properties in real-time by simply clicking on them in the 3D window. Most of the system information like active interactions, actors and global properties can be viewed in the visualization. In the Python script, the user can assign callback functions to keyboard input, so that all exposed system and particle properties (e.g., the temperature of the Langevin thermostat or the external force on a particle group) can be manipulated in real-time for testing and demonstration purposes. This allows to use ESPResSo also as an educational tool for a broad range of skill levels from basic physics to advanced simulation methods.

## Vi Parallel Output

Writing trajectories of a large set of particles to disk can be a bottleneck of parallel MD simulations. For serial I/O the costs of writing the trajectory include the communication to gather the data on the node that actually performs the write. A second aspect of the performance cost is the amount of data that has to be written to disk. Thus, in ESPResSo 4.0 we added the possibility to write binary files in parallel using the H5MD specification de Buyl et al. (2014) which utilizes the HDF5 file format The HDF Group (). H5MD files are intended to be self-contained and therefore enhance the portability of simulation output. This more readily allows one to use a wide class of analysis programs. The implementation is based on a wrapper around the HDF5 library. In typical ESPResSo simulations, the I/O performance is improved by using parallel output for particle numbers greater than approximately . Writing in this binary format is significantly faster than more traditional ASCII text. The H5MD files written by ESPResSo 4.0 may conveniently be read into Python using the H5py module Collette (2013). Simulation trajectories may also be visualized in VMD Humphrey et al. (1996) by using its H5MD plugin. The H5MD output file contains a field that holds the source script that was used for its creation. Thus the origin of the raw simulation data can always be traced back to the ESPResSo Python script, which helps tremendously in reproducing data.

For convenience, we still support ESPResSo’s VTF format which outputs system and particle properties as ASCII text. However, the use of Python as a scripting language opens the door to a large number of alternatives. Many aspects of a simulation, as well as simulation parameters can, for instance, be stored in a structured format using Python’s pickle module.

## Vii Software Engineering and Testing

As outlined in the previous section, the scripting interface is powered by a high-performance simulation core. The core of previous versions of ESPResSo was implemented in the C programming language. It was only natural to also upgrade the core in a manner that mirrors the new object-oriented approach of the interface. Therefore, ESPResSo 4.0 has switched from C to C++. C++ is a modern object-oriented programming language that allows for a more compact, readable code style. Just like Python, C++ benefits from an active community and a rich ecosystem of third-party support libraries, e.g., the Boost libraries mis (year). Thus, functionality that was previously implemented in the simulation core can now be provided externally. This reduces the code complexity and consequently the maintenance effort.

ESPResSo is designed with high-performance computing in mind. On high-performance compute clusters, the environment is often times very heterogeneous. For example, external libraries provided by different operating systems can vary noticeably, and only very recent C++ compilers fully support modern standards like C++11. To ensure that ESPResSo can be built and run on a large variety of different setups, it turned out to be necessary to comprehensively test it on different combinations of operating systems and compiler versions. ESPResSo has always included a suite of test cases that could automatically run small simulations to formally verify the physical correctness of the implementation by comparing to established results. In ESPResSo 4.0, not only have these tests been extended to cover more code, but also to specifically test the correctness of certain core aspects, a method called unit testing.

### vii.1 Improvements in the simulation core

The ESPResSo code base has been developed continuously for more than 15 years by a host of contributors of various backgrounds and styles. Thanks to the substantial overlap of language features between C and C++, it was possible to switch the language with minimal effort. Even though C and C++ share a common heritage, the programming paradigms are radically different and we are gradually transforming the old C code into more modern C++ code.

New language and library features of C++, especially ones which became available with the C++11 standard, are employed for more concise and expressive design of our software. ESPResSo is continuously used in day-to-day research activities and mainly developed by interested domain researchers; hence a gradual modernization strategy was necessary. Among the many algorithms whose code has been refactored, one central method is the traversal of particle pairs used for the evaluation of short-range interactions. In ESPResSo 4.0, the traversal has been separated from the kernel that evaluates the interactions for each particle pair. Repetition of the complicated routine for the pair energy, virial, and force calculation in the three types of particle cell systems supported by ESPResSo is hence avoided. This has been achieved by replacing several loops over the particles by a single function that takes the kernels as template parameters. This allows us to test key routines in isolation, independent of any interaction potentials. Added benefits to this decoupling are, improved performance and the ease by which new methods may be introduced by (future) collaborators.

### vii.2 Testing

ESPResSo contains a test suite to verify the correctness of the code. Hitherto, these tests were run by the maintainers before accepting code contributions, but features without tests were not covered. A formal measurement of the code coverage was introduced because it proved hard to detect which features are untested. In preparation for this release, we systematically improved coverage of previously untested areas of code. It is our experience that adding tests to existing code not only ensures correctness, but also improves the user experience by streamlining interfaces and avoiding inconsistencies and confusion.

For example, ESPResSo’s GPU and CPU lattice Boltzmann (LB) implementations had few systematic tests. However, this is an important method for many users of ESPResSo, and central for the correctness of many soft matter simulations Dünweg and Ladd (2009). We therefore introduced multiple new tests in this area that compare the flow fields computed by the LB method, as well as its coupling to particles, against stationary and time-dependent analytical solutions and reference data. These tests include direct verification of momentum and mass conservation, as well as validation of the statistical properties of the fluctuating variant of the LB and particle coupling. This helped us to improve the existing code as well as its interface.

In addition to the improved integration tests, this release includes some unit tests. Unit tests help identify issues on the level of individual functions which make up the algorithms and give more precise information where an issue occurred. This further allows changing the building blocks of an algorithm or their reuse across different simulation methods.

### vii.3 Continuous Integration

As mentioned earlier, the development of ESPResSo is organized through GitHub. Contributors can submit patches to fix a bug or introduce new features by means of a pull request. Whenever such a pull request is filed, the aforementioned set of unit and integration tests is automatically run and the outcome is communicated to GitHub where it is displayed in an informative fashion. This method of automatically running the test suite for changes is called continuous integration. Many online services exist which are free to use for open source projects like ESPResSo. Unfortunately, their free service plans are too limited for the vast variety of setups we need to test. Therefore we mirror the ESPResSo repository to a private instance of GitLab999GitLab is a free, self-hosted platform similar to GitHub, see https://gitlab.com.. This instance manages a fleet of dedicated in-house computers and distributes jobs for each configuration to test. On top of the automated testing, other developers are tasked with peer review of the patches to ensure a consistent integration into the code base.

We also try to cover as many combinations of operating systems, compilers, Python versions, and with and without GPU processors as possible. Our use of the Clang compiler mis () is particularly noteworthy as it is able to perform static code analysis to detect a large number of common programming errors and it is able to generate code to automatically recognize undefined behavior during runtime, thus often detecting bugs that are hard to discover for a human software developer.

## Viii External Fields

ESPResSo 4.0 implements a new extensible framework for coupling particles to external fields, i.e., fields that are not caused by particle interactions. An “external field“ is composed of a field source and a coupling, which can be combined as needed. An important use case for this feature are particles advected in a prescribed flow field under the assumption that they do not back-couple to that flow field. An example are soot particles on the nanometer scale in combustion exhausts Inci et al. (2017); Smiljanic et al. (2018). Here the field source would be a flow field specified on a regular grid which is interpolated to the particle positions. The coupling term would be a viscous friction, combining the velocity difference between the particle and the flow field and a friction constant to a force. Further use cases for external fields are complex boundary shapes as they are found in super-capacitor electrodes as for example described by Breitsprecher et al. (2018). In this application, electrodes with non-trivial geometry are considered. The electric field caused by the electrodes is pre-computed by solving Poisson’s equation for the system, and the resulting field is then superimposed on the inter-particles forces in the MD simulations as an interpolated external electric field. A snapshot of particles and electric field lines in such a system is shown in Fig. 3. The feature can also be used for Lennard-Jones particles in front of convex surfaces. Such a particle can interact with an extended part of the surface. This cannot be described by ESPResSo’s constraint feature, which only accounts for an interaction along the surface normal vector.

Let us now look at the two components of external fields in detail. The field source is described as a global scalar or vectorial field. Either, user-provided data is interpolated from a regular grid onto the particle positions, a constant value throughout the full domain is used, or an affine map according to

 u(r)=Ar+b, (1)

with a user-provided matrix and shift is evaluated.

The coupling maps a particle property and the value of the field source at the particle position to a force or potential. In the latter case the force is calculated by invoking the coupling with the gradient of the field. For example a gravitational field can be defined as a constant scalar field source coupling to the particle mass.

Presently there are several possibilities for the particle coupling, including but not limited to a particle’s scalar mass or charge, or the particle’s velocity vector using a viscous coupling, e.g.,

 Fi=−γ(vi−u). (2)

As a simple example we consider a two-dimensional system comprising particles with excluded-volume interactions in a square simulation box of side length . We combine interpolated flow field data with a viscous coupling. The flow field is a discretized static Taylor-Green vortex, specified by

 ux=cos(x)sin(y),uy=−sin(x)cos(y), (3)

where is the flow velocity. It has the vorticity

 ω≡∇×u=2|cos(x)cos(y)|. (4)

In Fig. 4 we show a snapshot of a simulation where the particles were initially placed at random positions distributed equally over the simulation volume. They are driven out of the areas with high vorticity by inertial forces and accumulate elsewhere. A complete simulation script can be found in the supplementary information.

## Ix Energy Minimization

Since bulk molecular dynamics are typically set up with random particle positions, initial configurations often have a lot of overlap between neighboring particles. This leads to problems with numerical stability due to very high energies. To remove this overlap there are two common methods. Limiting the forces between particles to a maximum value and performing an energy minimization step after the initial setup. While it was possible to cap the the forces in previous version of ESPResSo, now also energy minimization by the steepest descent method is available. It allows relaxing the position as well as the orientation of the particles by running the following iterative scheme:

 Δx,y,zi =sgn(Fx,y,zi)min(γFx,y,zi,Δmax) (5) xx,y,zi =xx,y,zi+Δx,y,zi, (6)

while . Here the and are the components of the position and force of the th particle. That is the update rule and limits are applied by component. and are relaxation parameters provided by the user. The orientation of the particle is relaxed in a similar fashion. Steepest descent energy minimization is available as an alternative to the integrator in ESPResSo.

In soft matter simulations the solvent is often treated on a coarse-grained level using a lattice Boltzmann fluid, which interacts with the particles of the system via so-called point coupling Dünweg and Ladd (2009). This type of coupling limits the maximal coupling strength for one particle which can lead to unrealistic transport properties for larger objects like colloids. One way to overcome this problem is to use so-called raspberry models Lobaskin and Dünweg (2004), rigid bodies of multiple particles, which couple to the fluid in multiple points to reduce interpolation artifacts and can achieve higher friction with the fluid. This can help to obtain better transport coefficients Fischer et al. (2015); de Graaf et al. (2015). To setup such a raspberry particle, its volume is filled with particles in a homogeneous manner. A common way to do this is to put the particles into a spherical constraint, give them a repulsive interaction and relax the energy. This can now be easily implemented in ESPResSo. As an example, Fig. 5 show the initial and relaxed configuration of 1400 particles in a spherical constraint with a soft-core Weeks-Chandler-Andersen interaction. Initially the particles are placed randomly in a cube contained by the sphere, then the energy minimization is performed. For the simulation script and detailed parameters please refer to the supplementary information.

## X Cluster Analysis

An analysis of the spatial clustering of particles is an important tool in many fields of soft matter science. Some examples include the nucleation of crystals Radu and Schilling (2014); Kratzer and Arnold (2015), clustering of magnetic nanoparticles in magnetic soft matter systems Butter et al. (2003); Cerdà et al. (2008a); Weeber et al. (2013); Donaldson and Kantorovich (2015), and the irreversible agglomeration of soot particles in combustion processes Attili et al. (2014); Inci et al. (2014, 2017). Clusters are typically defined by means of a criterion for a given pair of particles which takes the form of an equivalence relation. That is, if particles and are “neighbors” and particles and are “neighbors”, all three particles belong to a cluster. The pair criterion is chosen based on the application. For crystallization studies, a local bond order parameter Lechner and Dellago (2008) above a certain threshold on both particles is a common choice. For magnetic systems, a combination of distance and either dipole configuration or pair energies are used Cerdà et al. (2008a); Weeber et al. (2013). For the agglomeration of soot particles, the bonds created by the collision algorithms are used Inci et al. (2017).

While it is possible to perform cluster analysis in the post-processing state of a simulation, this requires the storage of a significant number of snapshots of the simulation. Due to the potentially large amount of storage space this might take, as well as the extra time it takes to write the snapshot to disk, this is often not a good approach. ESPResSo 4.0 therefore introduces an online cluster analysis which can be run from within the simulation. The implementation is loosely based on the Hoshen-Kopelman scheme Weeks et al. (1971), but does not use a lattice. The procedure is as follows: All pairs of particles are examined. If they are “neighbors” as determined by a user-specified pair criterion, there are several possible cases.

• If neither particle belongs to a cluster, a new clusters is created and both particles are marked as members.

• If one particle is part of a cluster and the other is not, the free particle is added to the cluster.

• If the two particles belong to different clusters, a note is made that the two clusters are one and the same and need to be merged later.

The clusters are labeled by numeric IDs, assigned in ascending order. After all pairs of particles are examined, the clusters marked as identical in the previous step are merged. Cluster IDs are traversed in descending order and merged clusters receive the lower of the two cluster IDs. In this way, the merging can be done in a single pass. ESPResSo 4.0 currently contains pair criteria based on inter-particle distance, pairwise short-range energy, and the presence of a bonded interaction. Further criteria can be added easily. The cluster analysis is organized around a ClusterStructure object, which provided the methods to run analysis as well as access to the clusters found. Furthermore, some analysis routines are provided on a per-cluster level, e.g., for a cluster’s radius of gyration, fractal dimension, or inertia tensor. Lastly, direct access to the particles making up a cluster is provided.

As an example, let us consider the simulation of a ferrofluid monolayer, i.e, a two-dimensional suspension of soft spheres carrying a magnetic dipole at their center. The example is loosely based on Ref. Cerdà et al. (2008a). Magnetic particles tend to form chain and ring-like clusters due to the non-isotropic nature of the dipole-dipole potential. For simplicity, we define particles as “neighbors” if the distance between their centers is less than , where denotes the particle diameter in the purely repulsive Lennard-Jones potential Weeks et al. (1971). The sample system contains 1000 magnetic particles at an area fraction of . The relative strength of the dipolar interactions to the thermal energy is

 λ=μ0m24πσ3kBT, (7)

where denotes the vacuum permittivity and the particles’ dipole moment. The dipolar interactions are calculated by means of the dipolar PM method Cerdà et al. (2008b) and the dipolar layer correction Bródka (2004). While the particle positions are confined to a plane, the magnetic dipoles can rotate freely in three dimensions. This is called a quasi-2D system.

In Fig. 6, a snapshot from the simulation along with a plot of the cluster composition of the suspension averaged over 1000 snapshots are shown. The full simulation script is provided in the supplementary information. Please note that the script is meant as an example. For a scientific study, larger systems, higher accuracies for the dipolar interactions, and longer sampling times are needed.

## Xi Particle-based Reactions

Chemical reactions can take place in many systems, for example water can autodissociate in positively charged protons and negatively charged OH groups, or there are other groups that can dissociate only partially in solvents like water, like the carboxylic group COOH. These groups can be found, for example, in weak polyelectrolytes, charge stabilized colloids, or proteins Richtering (2006), and their charging state depends on the pH-value Berg (2015) of the solution. Therefore, the presence of weak groups in charged polymers (polyelectrolytes) promote phenomena like protonation-configuration coupling Castelnovo et al. (2000); Shi et al. (2012). Protonation-configuration coupling means that the configuration of a protein or polymer depends on the protonation state of the titratable groups within the polymer. However, the protonation state itself depends on the configuration. This tight coupling introduces the need to investigate this phenomenon through computer simulations. Reactions can also give rise to a net attraction of particles, which are on average charge neutral Lund and Jönsson (2005, 2013). Here, the net charge of the particles fluctuates around its mean, thus allowing a particle to be temporarily positive or negative. This fluctuation effect then induces a net attraction. Investigating physical phenomena like the above with the correct statistics, requires introducing particle based reaction schemes in ESPResSo 4.0 which are described in the following.

A typical acidic reaction is shown in Fig. 7: an acid particle () may release a proton (), or vice versa, the deprotonated form of the acid () may take up a proton. Such protonation reactions change the electrostatic charge of the particles taking part in the reaction. Therefore, ESPResSo with its various electrostatic solvers Limbach et al. (2006), is well-suited to investigate electrostatic effects involved in reactions.

The first scheme that allowed for acid-base reactions in thermodynamic equilibrium using Molecular Dynamics and Monte Carlo simulations appeared in the 1990s Reed and Reed (1992) in the form of the constant-pH ensemble. Soon after, the reaction ensemble Smith (1994); Johnson et al. (1994) was introduced which allows for the simulation of arbitrary chemical reactions. In both reaction schemes, particle properties (such as the charge) need to be changed and particles need to be created or deleted based on probabilistic criteria. The reaction ensemble and the constant-pH method only differ in the probabilistic criteria Landsgesell et al. (2017a) and both methods are implemented in ESPResSo 4.0.

In the reaction ensemble, arbitrary reactions can be implemented in the simulation:

 z∑i=1νisi=0, (8)

where chemical species of type with stoichiometric coefficients are reacting Atkins and de Paula (2010). The acceptance probability in the reaction ensemble for a reaction from state to is given by Heath Turner et al. (2008)

 accRE,ξ(l|r)=min{1,(VΓ)¯¯¯νξz∏i=1[Nri!(Nri+ξνi)!]exp(−βΔEpot,r→l)}, (9)

where is the number of particles prior to a reaction and the “extent” of the reaction which is selected randomly with and proportional to the inverse of the temperature. Further parameters are the concentration-dependent reaction constant where is the reference concentration for which the change in the free enthalpy is tabulated, the potential energy difference with , the volume of the system and the total change in the number of molecules due to the reaction. The corresponding protonation and deprotonation reactions are usually performed after a fixed number of MD simulation steps with constant particle numbers Heath Turner et al. (2008).

In addition to the above standard reaction ensemble algorithm, we implemented a rare event sampling method on top of the reaction ensemble Landsgesell et al. (2017b), called the Wang-Landau Wang and Landau (2001) reaction ensemble algorithm. The Wang-Landau scheme is a so-called flat histogramme technique that samples all states with equal probability. This modification of the algorithm might prove useful when investigating systems with metastable states and rare transitions between them.

In the other chemical equilibrium sampling method which is present in literature—the constant-pH ensemble—the acceptance probability for a reaction is given by:

 acccpH,ξ(l|r)=min(1,exp[−β(ΔEpot±(ln(10)/β)(pH−pKa))]), (10)

where is the pH of an implicitly imposed proton reservoir, and where is used in the case of a association/dissociation. Additionally the dissociation proposal probability is proportional to the number of dissociable groups , while the association proposal probability is proportional to the number of deprotonated groups .

As a showcase for our ability to simulate chemcial reactions in ESPResSo 4.0  we present in figure 8 the titration curve of a linear weak polyelectrolyte. A typical titration curve which is obtained from experiment shows the degree of dissociation as a function of . The degree of association measures how many particles are in the associated state or in the dissociated state : . At low pH, a weak acid is mostly associated () while at higher pH a weak acid becomes less and less associated ().

The reaction ensemble implementation in ESPResSo also allows for simulations in the grand-canonical ensemble. This is, because the grand canonical simulation scheme Frenkel and Smit (1996) can be represented as a reaction \ce ¡=¿ A with , where is the bulk concentration of the species and the excess chemical potential of the species in the bulk. The needed excess chemical potential can be obtained via the Widoms insertion method Frenkel and Smit (1996). We also implemented this method in ESPResSo 4.0. It resembles a reaction which is constantly rejected while observing potential energy changes due to the insertion of particles. The excess chemical potential in a homogeneous system is then given by Frenkel and Smit (1996)

 μex=−kBTln(⟨e−β(E% pot(N+1)−Epot(N))⟩). (11)

## Xii Including Explicit Dipolar Polarization with Drude Oscillators

A particle’s electron density is redistributed by the local electric field, leading to a nonadditive molecular interaction. Especially the study of ionic liquids (molten salts) polarization plays a large role . Often, this is effectively included in the force-field (e.g. via reduced charges) to match continuum properties Leontyev and Stuchebrukhov (2011); Schmidt et al. (2010); Dommert et al. (2012); Dommert and Holm (2013); Kohagen et al. (2016). In systems where polarization effects are supposed to be inhomogeneous, e.g. at interfaces, incorporating the dynamic nature of polarization can refine the physical picture. Also simulations of bulk systems can benefit from explicit polarization. For example, the experimental values for the static dielectric constant and the self-diffusion coefficient of water were accurately reproduced using a polarizable model Lamoureux et al. (2006). Thermalized cold Drude oscillators Lamoureux and Roux (2003), available in ESPResSo 4.0, can be used to model this dynamic particle polarization. The basic idea is to add a “charge-on-a-spring” (Drude charge) to a particle (Drude core) that mimics an electron cloud which can be elongated to create a dynamically inducible dipole (see Fig. 9)Bordin et al. (2016). The energetic minimum of the Drude charge can be obtained self-consistently, which requires several iterations of the system’s electrostatics solver and is usually considered computationally expensive Yu et al. (2003). However, with thermalized cold Drude oscillators, the distance between Drude charge and core is coupled to a thermostat, so that it fluctuates around the self-consistent field solution. This thermostat is kept at a lower temperature compared to the global temperature to minimize the heat flow into the system. A second thermostat is applied on the center of mass of the Drude charge and core system to maintain the global temperature. The downside of this approach is that usually a smaller time step has to be used to resolve the high frequency oscillations of the spring to get a stable system Mitchell and Fincham (1993). In ESPResSo, the basic ingredients to simulate such a system are split into three bonds; their combined usage creates polarizable compounds:

1. A harmonic bond between charge and core.

2. For the cold thermostat, a Langevin-type thermalized distance bond is used for the Drude-Core distance.

3. A subtract PM short-range bond to cancel unwanted electrostatic interaction.

The system-wide thermostat has to be applied to the center of mass and not to the core particle directly. Therefore, the particles have to be excluded from the global thermostat, which is possible in ESPResSo by setting the temperature and friction coefficient of the Drude complex to zero. It thus remains possible to use a global Langevin thermostat for non-polarizable particles. As the Drude charge should not alter the charge or mass of the Drude complex, both properties have to be subtracted from the core when adding the Drude particle. In the following convention, we assume that the Drude charge is always negative. It is calculated via the spring constant and polarizability (in units of inverse volume) with . For polarizable molecules (i.e., connected particles, coarse grained models etc.) with partial charges on the molecule sites, the Drude charges will have electrostatic interaction with other cores of the molecule. Often, this is unwanted, as it might be already part of the force-field (via partial charges or parametrization of the covalent bonds). Without any further measures, the elongation of the Drude particles will be greatly affected by the partial charges of the molecule that are in close proximity. To prevent this, one has to cancel the interaction of the Drude charge with the partial charges of the cores within the molecule. This can be done with the subtracts PM short-range bond, which is also used to cancel the electrostatics between Drude core and Drude charge . This ensures that only the dipolar interaction inside the molecule remains. The error of this approximation increases with the share of the long-range part of the electrostatic interaction. In most cases, this error is negligible comparing the distance of the charges and the real-space cutoff of the PM electrostatics solver. In ESPResSo, helper methods assist setting up this exclusion. In combination with particle polarizability, the Thole correction Thole (1981) is often used to correct for overestimation of induced dipoles at short distances. Ultimately, it alters the short-range electrostatics of PM to result in a damped Coulomb interaction potential

 V(r)=qiqjr[1−e−sr(1+sr2)]. (12)

The Thole scaling coefficient is related to the polarizabilities and Thole damping parameters of the interacting species via

 s=(ai+aj)/2(αiαj)1/6. (13)

Note that for the Drude oscillators, the Thole correction should be applied only for the dipole part added by the Drude charge and not on the total core charge, which can be different for polarizable ions. Also note that the Thole correction acts between all dipoles, intra- and intermolecular. Again, the accuracy is related to the PM accuracy and the split between short-range and long-range electrostatics interaction. ESPResSo assists with the bookkeeping of mixed scaling coefficients and has convenient methods to set up all necessary Thole interactions.

## Xiii Active Particles

In ESPResSo 4.0 it is possible to simulate active systems, wherein individual particles constantly transduce (internal) energy to perform work in the form of motion. The combination of this self-propulsion with particle interactions gives rise to unique out-of-equilibrium behaviors, of which living systems offer a wealth of examples: mammals, fish and birds Helbing et al. (2000); Ballerini et al. (2008); Katz et al. (2011); Zhang et al. (2013); Silverberg et al. (2013), as well as microorganisms Woolley (2003); Riedel et al. (2005); Sokolov et al. (2007); Polin et al. (2009); Geyer et al. (2013); Ma et al. (2014); Reufer et al. (2014); Schwarz-Linek et al. (2016). The last decade and a half have seen the rapid development of man-made counterparts to microorganisms, following the first experimental realizations of “chemical swimmers” by Paxton et al. (2004) and Howse et al. (2007), respectively. These artificial swimmers reduce the complexity of living matter, but still display signature features of being out of equilibrium, such as the ability to perform work Maggi et al. (2016) and motility-induced phase separation Theurkauff et al. (2012); Palacci et al. (2013).

Significant progress has been made both theoretically and computationally in understanding the individual and collective behavior of swimmers using simple models. Arguably the most famous of these is the Vicsek model Vicsek et al. (1995), followed closely by the active Brownian model Ebeling et al. (1999); Stenhammar et al. (2013); Zheng et al. (2013). In our developments, we have drawn inspiration from the Active Brownian Model (ABM) and created a variant thereof: The Active Langevin Model (ALM), which we will describe here. This ALM can also be coupled to the ESPResSo (GPU) lattice Boltzmann (LB) fluid dynamics solver Arnold et al. (2013a); Röhm and Arnold (2012) to account for the characteristic dipolar flow field that is associated with self-propulsion by microorganisms Drescher et al. (2010, 2011) as well as chemical swimmers Campbell et al. (2018). Our implementation de Graaf et al. (2016a) is similar to the sub-lattice approach introduced by Nash et al. Nash et al. (2008); Nash (2010). We refer to this extension as the hydrodynamic ALM (HALM) and we will briefly touch upon it here.

The ALM equation of motion for the translation of the th swimmer’s center of mass is given by:

 M∂2∂t2ri=−Γ––t∂∂tri+f^ui−∑j≠i∇V(rij,Qi,Qj)+ξi,t(t). (14)

Here, we assume that the swimmer is fully shape anisotropic and we have introduced the following quantities: the mass , the translational diffusion tensor , a (co-rotating) unit vector that gives the direction of self-propelled motion , the self-propulsion force , and swimmer-swimmer interaction potential . The latter depends on the separation and orientation of the swimmers, as specified by quaterions . Thermal noise is introduced via , which satisfies and . The equation of motion for quaternions is lengthy and well-described in Ref. Martys and Mountain (1999), it is therefore not reproduced here.

Figure 10a shows a representation of an ALM swimmer. Each particle is assigned a direction of self-propulsion , along which it experiences a constant force . This self-propulsion force balances against the friction exerted by the implicit solvent, , leading to persistent directed motion with swim speed . Individual swimmers can have different propulsion forces and/or friction matrices in ESPResSo 4.0, providing users with tremendous flexibility in creating mixtures with a range of mobilities.

A swimmer’s direction of self-propulsion may be changed by rotational Brownian motion or by torques, either from external fields or collisions. The former leads to enhanced diffusion Howse et al. (2007) as shown in Fig. 10b. ALM in ESPResSo 4.0 currently does not allow for the coupling of translational and rotational degrees of freedom, as required for the simulation of L-shaped Kümmel et al. (2013) and chiral swimmers Wensink et al. (2014); future releases will bring such functionality. Lastly, note that ALM is not fully damped (Fig. 10b), i.e., there is an inertial component to the dynamics, which is uncommon in simulating (large) colloidal particles. However, ALM approaches the fully damped ABM when the friction coefficient is chosen to be large. In this limit, the time step should be suitably small to ensure that the algorithm produces the correct physics.

HALM introduces coupling between ALM and an LB fluid to account for hydrodynamic interactions between swimmers, particles, and obstacles. In developing HALM we have drawn inspiration from the work of Nash et al. Nash et al. (2008); Nash (2010), who formulated a similar coupling to an LB fluid for the ABM. A directional self-propulsion force is applied to the particle, similar to the way this is done in ALM. In LB, a single point particle gains an effective radius Ahlrichs and Dünweg (1999) or rather an inherent friction by coupling to the grid. This friction counter-balances the propulsive force leading to a constant swim speed. However, this also leads to force onto the LB fluid, while self-propulsion in a fluid should be force free. Therefore, a counter-force is applied to the LB fluid to ensure that momentum is not transferred into it globally, see Fig. 11. This force/counter-force pair induces a flow field with the typical leading-order dipolar contribution that characterizes self-propulsion. Figure 11 illustrates two combinations of the force pair leading to the well-known puller and pusher dipolar flow fields, representing, e.g., algae Drescher et al. (2010); Harris (2013) and spermatozoa Kantsler et al. (2013), respectively. We refer to Ref. de Graaf et al. (2016a) for full details of the implementation.

The reason for introducing ALM and HALM here, rather than relying on Brownian dynamics, is related to the way particles in ESPResSo couple to the LB fluid. Using these algorithms in tandem, the effect of the hydrodynamic interactions between swimmers may be disentangled, in much the same way as was done for passive particles Röhm et al. (2014). Careful tuning of the relevant hydrodynamic parameters and the self-propulsion forces is key to reproducing low-Reynolds number hydrodynamic solutions, we refer to Ref. de Graaf and Stenhammar (2017a) for an in-depth discussion. These authors have verified that the Nash et al. implementation Nash et al. (2008); Nash (2010) and HALM have the same near-field flow characteristics, long-range hydrodynamic retardation effects, and particle dynamics de Graaf and Stenhammar (2017b). Thus, the inertial component to HALM does not significantly affect the dynamics of our model.

ALM and HALM can both be used in conjunction with the raspberry method of creating shape-anisotropic particles Lobaskin and Dünweg (2004); Chatterji and Horbach (2005); Fischer et al. (2015); de Graaf et al. (2015). This enables, for example, the study of the effect of polarization and roughness on motility-induced clustering Ilse et al. (2016). HALM combined with raspberry particles leads to hydrodynamic multipole moments beyond the leading dipole moment for the self-propulsion of shape-anisotropic particles in fluids de Graaf et al. (2016a). This property has been exploited to determine the effect of these hydrodynamic moments on the motion of these swimmers in confining geometries de Graaf et al. (2016b), as well as their interaction with passive particles in their surrounding de Graaf et al. (2016a); de Graaf and Stenhammar (2017a). The use of raspberry particles in combination with HALM removes most of the lattice artifacts that point-particle HALM suffers from de Graaf et al. (2016a). We therefore recommend employing this combination within ESPResSo to properly account for any rigid body dynamics in a fluid, including rotational coupling due to shape anisotropies.

## Xiv Conclusions

In this work, we outlined the major revisions that ESPResSo 4.0 has undergone, and the benefits that the user will experience from these changes. These include the reimplementation of the user interface in Python and an overhaul of the core architecture to modern software development paradigms. We have also reported new features and presented them accompanied by specific use cases. This enables the simulation of systems at the forefront of soft matter research, including active matter and catalytic reactions.

We foresee a bright future for ESPResSo as a software package and as a research project. Together with our collaborators, we will further improve the documentation and usability of this platform. In addition, algorithmic improvements are planned, which include: (i) Adaptive grid refinement to the electrokinetics code. (ii) A load-balancing scheme for efficient simulation of heterogeneous systems. (iii) Lees-Edwards boundary conditions for rheological measurements. We cordially invite new users to try out ESPResSo as a simulation tool for their research and participate in its continued development.

The current status of the package and the latest tutorials and documentation can be found on the project website http://espressomd.org, or on GitHub https://github.com/espressomd/.

## Acknowledgments

We would like to acknowledge more than 100 researchers who contributed over the last fifteen years to the ESPResSo software, and whose names can be found on our website http://espressomd.org or in the AUTHORS file distributed with ESPResSo. We are also grateful to our colleagues of the SFB 716 for numerous suggestions for extending the capabilities of ESPResSo and in helping us to improve our software. As part of a collaboration Smiljanic et al. (2018), Milena Smiljanic contributed code to the cluster analysis framework, particularly to the analysis routines on the single cluster level.

CH, KS, and FW gratefully acknowledge funding by the German Science Foundation (DFG) through the collaborative research center SFB 716 within TP C5, the SimTech cluster of excellence (EXC 310), and grants HO 1108/25-1, HO 1108/26-1, AR 593/7-1, HO 1108/28-1. MK, JdG, and CH thank the DFG for funding through the SPP 1726 Microswimmers–From Single Particle Motion to Collective Behavior. JdG further acknowledges funding from an NWO Rubicon Grant (#680501210) and a Marie Skłodowska-Curie Intra European Fellowship (G.A. No. 654916) within Horizon 2020. The simulations were partially performed on the bwUniCluster funded by the Ministry of Science, Research and Arts and the Universities of the State of Baden-Württemberg, Germany, within the framework program bwHPC.

## Author Contributions

Supervision: CH; Funding Acquisition, CH; Resources, CH; Conceptualization: CH, RW, and FW. Writing: All authors contributed to the preparation of the manuscript. Software and Validation: FW, RW, and KS are core developers of ESPResSo. The main code contributors of individual features discussed in this article are the following. Test infrastructure: KS, MK. Visualization: KB, MK. Particle-based reactions: JL. Drude oscillators: KB. External fields: FW. Steepest descent energy minimization: FW. Cluster analysis: RW. Active particles: JdG, HM. H5MD parallel output: KS. VTF output: DS. A full list of code contributions can be found at https://github.com/espressomd/espresso/graphs/contributors.

## References

• de Gennes (1992) P. G. de Gennes, Reviews of Modern Physics 64, 645 (1992).
• Doi (2013) M. Doi, Soft Matter Physics (Oxford University Press, 2013).
• Barrat and Hansen (2003) J.-L. Barrat and J.-P. Hansen, Basic Concepts fo Simple and Complex Liquids (Cambridge University Press, 2003) iSBN: 0-521-78953-2.
• Doi and Edwards (1988) M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Vol. 73 (Oxford University Press, 1988).
• Rubinstein and Colby (2003) M. Rubinstein and R. H. Colby, Polymer Physics (Oxford University Press, Oxford, UK, 2003).
• Verwey and Overbeek (1948) E. J. Verwey and J. T. G. Overbeek, Theory of the stability of Lyophobic Colloids (Elsevier, Amsterdam, 1948).
• Löwen (2001) H. Löwen, Journal of Physics: Condensed Matter 13, R415 (2001).
• Chandrasekhar (1992) S. Chandrasekhar, Liquid Crystals (Cambridge, 1992).
• Hamley (2003) I. V. Hamley, Introduction to Soft Matter (Wiley, 2003) iSBN: 0-471-89952-6.
• Nelson (2004) P. Nelson, Biological Physics - Energy, Information, Life (Freeman, 2004) iSBN: 0-7167-4372-8.
• Levental et al. (2007) I. Levental, P. C. Georges,  and P. A. Janmey, Soft Matter 3, 299 (2007).
• Ubbink et al. (2008) J. Ubbink, A. Burbidge,  and R. Mezzenga, Soft Matter 4, 1569 (2008).
• Polarz and Antonietti (2002) S. Polarz and M. Antonietti, Chemical Communications , 2593 (2002).
• Kroy and Frey (2005) K. Kroy and E. Frey, Annalen der Physik 14, 20 (2005).
• Frenkel (2002) D. Frenkel, Science 296, 65 (2002).
• Seifert (2012) U. Seifert, Reports on Progress in Physics 75, 126001 (2012).
• Cates (2012) M. E. Cates, Reports on Progress in Physics 75, 042601 (2012).
• Fodor and Marchetti (2018) É. Fodor and M. Marchetti, Physica A: Statistical Mechanics and its Applications 504, 106 (2018).
• Alder and Wainwright (1957) B. J. Alder and T. E. Wainwright, Journal of Chemical Physics 27, 1208 (1957).
• Alder and Wainwright (1962) B. Alder and T. Wainwright, Physical Review 127, 359 (1962).
• Pusey et al. (1994) P. Pusey, W. C. Poon, S. Ilett,  and P. Bartlett, Journal of Physics: Condensed Matter 6, A29 (1994).
• Bates and Frenkel (1998) M. A. Bates and D. Frenkel, The Journal of Chemical Physics 109, 6193 (1998).
• van Roij et al. (1999) R. van Roij, M. Dijkstra,  and J.-P. Hansen, Physical Review E 59, 2010 (1999).
• Leunissen et al. (2005) M. E. Leunissen, C. G. Christova, A. P. Hynninen, C. P. Royall, A. I. Campbell, A. Imhof, M. Dijkstra, R. van Roij,  and A. van Blaaderen, Nature 437, 235 (2005).
• Camp et al. (2000) P. J. Camp, J. C. Shelley,  and G. N. Patey, Physical Review Letters 84, 115 (2000).
• Wang et al. (2002) Z. Wang, C. Holm,  and H. W. Müller, Physical Review E 66, 021405 (2002).
• Klapp and Schoen (2004) S. H. L. Klapp and M. Schoen, Journal of Molecular Liquids 109, 55 (2004).
• Klinkigt et al. (2013) M. Klinkigt, R. Weeber, S. Kantorovich,  and C. Holm, Soft Matter 9, 3535 (2013).
• Stevens and Kremer (1993) M. J. Stevens and K. Kremer, Physical Review Letters 71, 2228 (1993).
• Dobrynin et al. (1996) A. V. Dobrynin, M. Rubinstein,  and S. P. Obukhov, Macromolecules 29, 2974 (1996).
• Micka et al. (1999) U. Micka, C. Holm,  and K. Kremer, Langmuir 15, 4033 (1999).
• Limbach et al. (2002) H. J. Limbach, C. Holm,  and K. Kremer, Europhysics Letters 60, 566 (2002).
• Schneider and Linse (2002) S. Schneider and P. Linse, The European Physical Journal E 8, 457 (2002).
• Yan and de Pablo (2003) Q. Yan and J. J. de Pablo, Physical Review Letters 91, 018301 (2003).
• Mann et al. (2004) B. A. Mann, R. Everaers, C. Holm,  and K. Kremer, Europhysics Letters 67, 786 (2004).
• Weeber et al. (2012) R. Weeber, S. Kantorovich,  and C. Holm, Soft Matter 8, 9923 (2012).
• Weeber et al. (2018a) R. Weeber, M. Hermes, A. M. Schmidt,  and C. Holm, Journal of Physics: Condensed Matter 30, 063002 (2018a).
• Plimpton (1995) S. J. Plimpton, Journal of Computational Physics 117, 1 (1995).
• Guzman et al. (2018) H. V. Guzman, H. Kobayashi, N. Tretyakov, A. C. Fogarty, K. Kreis, J. Krajniak, C. Junghans, K. Kremer,  and T. Stuehn, arXiv preprint arXiv:1806.10841  (2018).
• Berendsen et al. (1995) H. J. C. Berendsen, D. van der Spoel,  and R. van Drunen, Computer Physics Communications 91, 43 (1995).
• van der Spoel et al. (2005) D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark,  and H. J. C. Berendsen, Journal of Computational Chemistry 26, 1701 (2005).
• Pronk et al. (2013) S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess,  and E. Lindahl, Bioinformatics 29, 845 (2013).
• Nelson et al. (1996) M. Nelson, W. Humphrey, A. Gursoy, A. Dalke, L. Kale, R. D. Skeel,  and K. Schulten, International Journal of Supercomputer Applications 10, 251 (1996).
• Phillips et al. (2005) J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale,  and K. Schulten, Journal of Computational Chemistry 26, 1781 (2005).
• Phillips et al. (2014) J. C. Phillips, J. E. Stone, K. L. Vandivort, T. G. Armstrong, J. M. Wozniak, M. Wilde,  and K. Schulten, in Proceedings of the 1st First Workshop for High Performance Technical Computing in Dynamic Languages (IEEE Press, 2014) pp. 6–17.
• Pearlman et al. (1995) D. A. Pearlman, D. A. Case, J. W. Caldwell, W. R. Ross, I. T. E. Cheatham, S. DeBolt, D. Ferguson, G. Seibel,  and P. Kollman, Computer Physics Commications 91, 1 (1995).
• Case et al. (2005) D. A. Case, T. E. Cheatham, T. Darden, H. Gohlke, R. Luo, K. M. Merz, A. Onufriev, C. Simmerling, B. Wang,  and R. Woods, Journal of Computational Chemistry 26, 1668 (2005).
• Arnold et al. (2013a) A. Arnold, O. Lenz, S. Kesselheim, R. Weeber, F. Fahrenberger, D. Röhm, P. Košovan,  and C. Holm, in Meshfree Methods for Partial Differential Equations VI, Lecture Notes in Computational Science and Engineering, Vol. 89, edited by M. Griebel and M. A. Schweitzer (Springer Berlin Heidelberg, 2013) pp. 1–23.
• Deserno and Holm (1998a) M. Deserno and C. Holm, Journal of Chemical Physics 109, 7678 (1998a).
• Deserno and Holm (1998b) M. Deserno and C. Holm, Journal of Chemical Physics 109, 7694 (1998b).
• Arnold and Holm (2002a) A. Arnold and C. Holm, Computer Physics Communications 148, 327 (2002a).
• Arnold and Holm (2002b) A. Arnold and C. Holm, Chemical Physics Letters 354, 324 (2002b).
• Arnold et al. (2002a) A. Arnold, J. de Joannis,  and C. Holm, Journal of Chemical Physics 117, 2496 (2002a).
• Arnold et al. (2002b) A. Arnold, J. de Joannis,  and C. Holm, Journal of Chemical Physics 117, 2503 (2002b).
• Arnold and Holm (2005a) A. Arnold and C. Holm, in Advanced Computer Simulation Approaches for Soft Matter Sciences II, Advances in Polymer Sciences, Vol. II, edited by C. Holm and K. Kremer (Springer, Berlin, 2005) pp. 59–109.
• Arnold and Holm (2005b) A. Arnold and C. Holm, Journal of Chemical Physics 123, 144103 (2005b).
• Arnold et al. (2006) A. Arnold, B. A. Mann,  and C. Holm, in Computer Simulations in Condensed Matter: from Materials to Chemical Biology, Lecture Notes in Physics No. 1, edited by M. Ferrario, G. Ciccotti,  and K. Binder (Springer, Berlin, Germany, 2006) pp. 193–222.
• Tyagi et al. (2007) S. Tyagi, A. Arnold,  and C. Holm, Journal of Chemical Physics 127, 154723 (2007).
• Tyagi et al. (2008) S. Tyagi, A. Arnold,  and C. Holm, Journal of Chemical Physics 129, 204102 (2008).
• Tyagi et al. (2010) C. Tyagi, M. Süzen, M. Sega, M. Barbosa, S. S. Kantorovich,  and C. Holm, The Journal of Chemical Physics 132, 154112 (2010).
• Arnold et al. (2013b) A. Arnold, K. Breitsprecher, F. Fahrenberger, S. Kesselheim, O. Lenz,  and C. Holm, Entropy 15, 4569 (2013b).
• Fahrenberger and Holm (2014) F. Fahrenberger and C. Holm, Physical Review E 90, 063304 (2014).
• Fischer et al. (2015) L. P. Fischer, T. Peter, C. Holm,  and J. de Graaf, The Journal of Chemical Physics 143, 084107 (2015).
• de Graaf et al. (2015) J. de Graaf, T. Peter, L. P. Fischer,  and C. Holm, The Journal of Chemical Physics 143, 084108 (2015).
• de Graaf et al. (2016a) J. de Graaf, H. Menke, A. J. Mathijssen, M. Fabritius, C. Holm,  and T. N. Shendruk, The Journal of Chemical Physics 144, 134106 (2016a).
• de Graaf et al. (2016b) J. de Graaf, A. J. Mathijssen, M. Fabritius, H. Menke, C. Holm,  and T. N. Shendruk, Soft Matter 12, 4704 (2016b).
• Inci et al. (2014) G. Inci, A. Arnold, A. Kronenburg,  and R. Weeber, Aerosol Science and Technology 48, 842 (2014).
• Inci et al. (2017) G. Inci, A. Kronenburg, R. Weeber,  and D. Pflüger, Flow, Turbulence and Combustion , 1 (2017).
• Schober et al. (2016) C. Schober, D. Keerl, M. J. Lehmann,  and M. Mehl, in Proceedings of the VII International Conference on Coupled Problems in Science and Engineering, edited by M. Papadrakakis, E. Oñate,  and B. Schrefler (International Center for Numerical Methods in Engineering, 2016).
• Röhm and Arnold (2012) D. Röhm and A. Arnold, The European Physical Journal Special Topics 210, 89 (2012).
• Cimrák et al. (2012) I. Cimrák, M. Gusenbauer,  and T. Schrefl, Computers & Mathematics with Applications 64, 278 (2012).
• Cimrak et al. (2014) I. Cimrak, M. Gusenbauer,  and I. Jančigová, Computer Physics Communications 185, 900 (2014).
• Cimrák et al. (2013) I. Cimrák, I. Jancigová, K. Bachratá,  and H. Bachratỳ, in III International Conference on Particle-based Methods–Fundamentals and Applications PARTICLES, Vol. 2013 (2013) pp. 133–144.
• Rempfer et al. (2016) G. Rempfer, G. B. Davies, C. Holm,  and J. de Graaf, The Journal of Chemical Physics 145, 044901 (2016).
• Guckenberger et al. (2016) A. Guckenberger, M. P. Schraml, P. G. Chen, M. Leonetti,  and S. Gekle, Computer Physics Communications 207, 1 (2016).
• Bächer et al. (2017) C. Bächer, L. Schrack,  and S. Gekle, Physical Review Fluids 2, 013102 (2017).
• Kratzer et al. (2013) K. Kratzer, A. Arnold,  and R. J. Allen, Journal of Chemical Physics 138, 164112 (2013).
• Kratzer et al. (2014) K. Kratzer, J. T. Berryman, A. Taudt, J. Zeman,  and A. Arnold, Computer Physics Communications 185, 1875 (2014).
• Samin et al. (2013) S. Samin, Y. Tsori,  and C. Holm, Physical Review E 87 (2013), 10.1103/PhysRevE.87.052128.
• Paszke et al. (2017) A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga,  and A. Lerer, in Neural Information Processing Systems (2017).
• Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot,  and E. Duchesnay, Journal of Machine Learning Research 12, 2825 (2011).
• Chollet (2017) F. Chollet, Deep learning with Python (Manning Publications Co., 2017).
• Meurer et al. (2017) A. Meurer, C. P. Smith, M. Paprocki, O. Čertík, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh, T. Rathnayake, S. Vig, B. E. Granger, R. P. Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson, F. Pedregosa, M. J. Curry, A. R. Terrel, v. Roučka, A. Saboo, I. Fernando, S. Kulal, R. Cimrman,  and A. Scopatz, PeerJ Computer Science 3, e103 (2017).
• Stein and Joyner (2005) W. Stein and D. Joyner, ACM SIGSAM Bulletin 39, 61 (2005).
• Hunter (2007) J. D. Hunter, Computing In Science & Engineering 9, 90 (2007).
• Ramachandran and Varoquaux (2011) P. Ramachandran and G. Varoquaux, Computing in Science & Engineering 13, 40 (2011).
• Jones et al. (year) E. Jones, T. Oliphant, P. Peterson, et al.“SciPy: Open source scientific tools for Python,”  (2001–2019).
• McKinney (2010) W. McKinney, in Proceedings of the 9th Python in Science Conference, edited by S. van der Walt and J. Millman (2010) pp. 51 – 56.
• Reith et al. (2003) D. Reith, M. Pütz,  and F. Müller-Plathe, Journal of Computational Chemistry 24, 1624 (2003).
• Kluyver et al. (2016) T. Kluyver, B. Ragan-Kelley, F. Pérez, B. E. Granger, M. Bussonnier, J. Frederic, K. Kelley, J. B. Hamrick, J. Grout, S. Corlay, P. Ivanov, D. Avila, S. Abdalla,  and C. Willing, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, edited by F. Loiyides and B. Schmidt (IOS Press, 2016) pp. 87–90.
• Hockney and Eastwood (1988) R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles (IOP, London, 1988).
• Arnold et al. (2013c) A. Arnold, F. Fahrenberger, C. Holm, O. Lenz, M. Bolten, H. Dachsel, R. Halver, I. Kabadshow, F. Gähler, F. Heber, J. Iseringhausen, M. Hofmann, M. Pippig, D. Potts,  and G. Sutmann, Physical Review E 88, 063308 (2013c).
• Nestler (2016) F. Nestler, Applied Numerical Mathematics 105, 25 (2016).
• Weeber et al. (2018b) R. Weeber, F. Nestler, F. Weik, M. Pippig, D. Potts,  and C. Holm, Arxiv pre-print arXiv:1808.10341  (2018b).
• Raikher and Stolbov (2003) Y. L. Raikher and O. V. Stolbov, JMMM 258/259, 477 (2003).
• Morozov et al. (2009) K. Morozov, M. Shliomis,  and H. Yamaguchi, Physical Review E 79, 040801 (2009).
• Ayachit (2015) U. Ayachit, The ParaView Guide: A Parallel Visualization Application (Kitware, Inc., USA, 2015).
• Humphrey et al. (1996) W. Humphrey, A. Dalke,  and K. Schulten, Journal of Molecular Graphics 14, 33 (1996).
• Fletcher and Liebscher (2005) M. Fletcher and R. Liebscher, “PyOpenGL–the Python OpenGL binding,”  (2005).
• Shreiner (1999) D. Shreiner, OpenGL Reference Manual: The Official Reference Document to OpenGL, Version 1.2, 3rd ed. (Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1999).
• Röhm (2011) D. Röhm, Lattice Boltzmann Simulations on GPUs, Diplomarbeit, University of Stuttgart (2011).
• de Buyl et al. (2014) P. de Buyl, P. H. Colberg,  and F. Höfling, Computer Physics Communications 185, 1546 (2014).
• (103)
• Collette (2013) A. Collette, Python and HDF5 (O’Reilly, 2013).
• mis (year) “Boost C++ libraries,” https://www.boost.org (1998–2019).
• Dünweg and Ladd (2009) B. Dünweg and A. J. C. Ladd, in Advanced Computer Simulation Approaches for Soft Matter Sciences III, Advances in Polymer Science, Vol. 221 (Springer-Verlag Berlin, Berlin, Germany, 2009) pp. 89–166.
• (107)
• Smiljanic et al. (2018) M. Smiljanic, R. Weeber, D. Pflüger, C. Holm,  and A. Kronenburg, Submitted to European Physics Journal Special Topics  (2018).
• Breitsprecher et al. (2018) K. Breitsprecher, C. Holm,  and S. Kondrat, ACS Nano 12, 9733 (2018).
• Lobaskin and Dünweg (2004) V. Lobaskin and B. Dünweg, New Journal of Physics 6, 54 (2004).
• Radu and Schilling (2014) M. Radu and T. Schilling, Europhysics Letters 105, 26001 (2014).
• Kratzer and Arnold (2015) K. Kratzer and A. Arnold, Soft Matter 11, 2174 (2015).
• Butter et al. (2003) K. Butter, P. H. H. Bomans, P. M. Frederik, G. J. Vroege,  and A. P. Philipse, Nature Materials 2, 88 (2003).
• Cerdà et al. (2008a) J. J. Cerdà, S. Kantorovich,  and C. Holm, Journal of Physics: Condensed Matter 20, 204125 (2008a).
• Weeber et al. (2013) R. Weeber, M. Klinkigt, S. Kantorovich,  and C. Holm, Journal of Chemical Physics 139, 214901 (2013).
• Donaldson and Kantorovich (2015) J. G. Donaldson and S. S. Kantorovich, Nanoscale 7, 3217 (2015).
• Attili et al. (2014) A. Attili, F. Bisetti, M. E. Mueller,  and H. Pitsch, Combustion and Flame 161, 1849 (2014).
• Lechner and Dellago (2008) W. Lechner and C. Dellago, The Journal of Chemical Physics 129, 114707 (2008).
• Weeks et al. (1971) J. D. Weeks, D. Chandler,  and H. C. Andersen, The Journal of Chemical Physics 54, 5237 (1971).
• Cerdà et al. (2008b) J. J. Cerdà, V. Ballenegger, O. Lenz,  and C. Holm, Journal of Chemical Physics 129, 234104 (2008b).
• Bródka (2004) A. Bródka, Chemical Physics Letters 400, 62 (2004).
• Richtering (2006) W. Richtering, Smart Colloidal Materials 133, 9 (2006).
• Berg (2015) J. M. Berg, ed., Biochemistry, 8th ed. (Freeman, New York, NY (USA), 2015).
• Castelnovo et al. (2000) M. Castelnovo, P. Sens,  and J.-F. Joanny, European Physical Journal E: Soft Matter 1, 115 (2000).
• Shi et al. (2012) C. Shi, J. A. Wallace,  and J. K. Shen, Biophysical Journal 102, 1590 (2012).
• Lund and Jönsson (2005) M. Lund and B. Jönsson, Biochemistry 44, 5722 (2005).
• Lund and Jönsson (2013) M. Lund and B. Jönsson, Quarterly Reviews of Biophysics 46, 265 (2013).
• Limbach et al. (2006) H. J. Limbach, A. Arnold, B. A. Mann,  and C. Holm, Computer Physics Communications 174, 704 (2006).
• Reed and Reed (1992) C. E. Reed and W. F. Reed, Journal of Chemical Physics 96, 1609 (1992).
• Smith (1994) E. R. Smith, Journal of Statistical Physics 77, 449 (1994).
• Johnson et al. (1994) J. K. Johnson, A. Z. Panagiotopoulos,  and K. E. Gubbins, Molecular Physics 81, 717 (1994).
• Landsgesell et al. (2017a) J. Landsgesell, C. Holm,  and J. Smiatek, The European Physical Journal Special Topics 226, 725 (2017a).
• Atkins and de Paula (2010) P. W. Atkins and J. de Paula, Physical Chemistry (Oxford Univ. Press, Oxford (UK), 2010).
• Heath Turner et al. (2008) C. Heath Turner, J. K. Brennan, M. Lisal, W. R. Smith, J. Karl Johnson,  and K. E. Gubbins, Molecular Simulation 34, 119 (2008).
• Landsgesell et al. (2017b) J. Landsgesell, C. Holm,  and J. Smiatek, Journal of Chemical Theory and Computation 13, 852 (2017b).
• Wang and Landau (2001) F. Wang and D. P. Landau, Physical Review E 64, 056101 (2001).
• Frenkel and Smit (1996) D. Frenkel and B. Smit, Understanding Molecular Simulation, 1st ed. (Academic Press, San Diego, 1996).
• Leontyev and Stuchebrukhov (2011) I. Leontyev and A. Stuchebrukhov, Physical Chemistry Chemical Physics 13, 2613 (2011).
• Schmidt et al. (2010) J. Schmidt, C. Krekeler, F. Dommert, Y. Zhao, R. Berger, L. Delle Site,  and C. Holm, Journal of Physical Chemistry B 114, 6150 (2010).
• Dommert et al. (2012) F. Dommert, K. Wendler, R. Berger, L. D. Site,  and C. Holm, ChemPhysChem 13, 1625 (2012).
• Dommert and Holm (2013) F. Dommert and C. Holm, Physical Chemistry Chemical Physics 15, 2037 (2013).
• Kohagen et al. (2016) M. Kohagen, P. E. Mason,  and P. Jungwirth, The Journal of Physical Chemistry B 120, 1454 (2016).
• Lamoureux et al. (2006) G. Lamoureux, E. Harder, I. Vorobyov, B. Roux,  and A. MacKerell, Chemical Physics Letters 418, 245 (2006).
• Lamoureux and Roux (2003) G. Lamoureux and B. Roux, Journal of Chemical Physics 119, 3025 (2003).
• Bordin et al. (2016) J. R. Bordin, R. Podgornik,  and C. Holm, The European Physical Journal Special Topics 225, 1693 (2016).
• Yu et al. (2003) H. Yu, T. Hansson,  and W. F. van Gunsteren, The Journal of Chemical Physics 118, 221 (2003).
• Mitchell and Fincham (1993) P. Mitchell and D. Fincham, Journal of Physics: Condensed Matter 5, 1031 (1993).
• Thole (1981) B. T. Thole, Chemical Physics 59, 341 (1981).
• Helbing et al. (2000) D. Helbing, I. Farkas,  and T. Vicsek, Nature 407, 487 (2000).
• Ballerini et al. (2008) M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale,  and V. Zdravkovic, Proceedings of the National Academy of Sciences 105, 1232 (2008).
• Katz et al. (2011) Y. Katz, K. Tunstrøm, C. Ioannou, C. Huepe,  and I. Couzin, Proceedings of the National Academy of Sciences 108, 18720 (2011).
• Zhang et al. (2013) J. Zhang, W. Klingsch, A. Schadschneider,  and A. Seyfried, in Traffic and Granular Flow ’11, edited by V. V. Kozlov, A. P. Buslaev, A. S. Bugaev, M. V. Yashina, A. Schadschneider,  and M. Schreckenberg (Springer (Berlin/Heidelberg), 2013) p. 241.
• Silverberg et al. (2013) J. Silverberg, M. Bierbaum, J. Sethna,  and I. Cohen, Physical Review Letters 110, 228701 (2013).
• Woolley (2003) D. Woolley, Reproduction 126, 259 (2003).
• Riedel et al. (2005) I. Riedel, K. Kruse,  and J. Howard, Science 309, 300 (2005).
• Sokolov et al. (2007) A. Sokolov, I. S. Aranson, J. O. Kessler,  and R. E. Goldstein, Physical Review Letters 98, 158102 (2007).
• Polin et al. (2009) M. Polin, I. Tuval, K. Drescher, J. Gollub,  and R. Goldstein, Science 325, 487 (2009).
• Geyer et al. (2013) V. Geyer, F. Jülicher, J. Howard,  and B. Friedrich, Proceedings of the National Academy of Sciences 110, 18058 (2013).
• Ma et al. (2014) R. Ma, G. Klindt, I. Riedel-Kruse, F. Jülicher,  and B. Friedrich, Physical Review Letters 113, 048101 (2014).
• Reufer et al. (2014) M. Reufer, R. Besseling, J. Schwarz-Linek, V. Martinez, A. Morozov, J. Arlt, D. Trubitsyn, F. Ward,  and W. Poon, Biophysical Journal 106, 37 (2014).
• Schwarz-Linek et al. (2016) J. Schwarz-Linek, J. Arlt, A. Jepson, A. Dawson, T. Vissers, D. Miroli, T. Pilizota, V. A. Martinez,  and W. C. Poon, Colloids and Surfaces B: Biointerfaces 137, 2 (2016).
• Paxton et al. (2004) W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St. Angelo, Y. Cao, T. E. Mallouk, P. E. Lammert,  and V. H. Crespi, Journal of the American Chemical Society 126, 13424 (2004).
• Howse et al. (2007) J. R. Howse, R. A. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh,  and R. Golestanian, Physical Review Letters 99, 048102 (2007).
• Maggi et al. (2016) C. Maggi, J. Simmchen, F. Saglimbeni, J. Katuri, M. Dipalo, F. De Angelis, S. Sanchez,  and R. Di Leonardo, Small 12, 446 (2016).
• Theurkauff et al. (2012) I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert,  and L. Bocquet, Physical Review Letters 108, 268303 (2012).
• Palacci et al. (2013) J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine,  and P. M. Chaikin, Science 339, 936 (2013).
• Vicsek et al. (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen,  and O. Shochet, Physical Review Letters 75, 1226 (1995).
• Ebeling et al. (1999) W. Ebeling, F. Schweitzer,  and B. Tilch, BioSystems 49, 17 (1999).
• Stenhammar et al. (2013) J. Stenhammar, A. Tiribocchi, R. Allen, D. Marenduzzo,  and M. Cates, Physical Review Letters 111, 145702 (2013).
• Zheng et al. (2013) X. Zheng, B. ten Hagen, A. Kaiser, M. Wu, H. Cui, Z. Silber-Li,  and H. Löwen, Physical Review E 88, 032304 (2013).
• Drescher et al. (2010) K. Drescher, R. Goldstein, N. Michel, M. Polin,  and I. Tuval, Physical Review Letters 105, 168101 (2010).
• Drescher et al. (2011) K. Drescher, J. Dunkel, L. Cisneros, S. Ganguly,  and R. Goldstein, Proceedings of the National Academy of Sciences 108, 10940 (2011).
• Campbell et al. (2018) A. I. Campbell, S. J. Ebbens, P. Illien,  and R. Golestanian, arXiv preprint arXiv:1802.04600  (2018).
• Nash et al. (2008) R. Nash, R. Adhikari,  and M. Cates, Physical Review E 77, 026709 (2008).
• Nash (2010) R. W. Nash, Efficient lattice Boltzmann simulations of self-propelled particles with singular forces, Ph.D. thesis, The University of Edinburgh (2010).
• Martys and Mountain (1999) N. S. Martys and R. D. Mountain, Physical Review E 59, 3733 (1999).
• Kümmel et al. (2013) F. Kümmel, B. ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen,  and C. Bechinger, Physical Review Letters 110, 198302 (2013).
• Wensink et al. (2014) H. Wensink, V. Kantsler, R. Goldstein,  and J. Dunkel, Physical Review E 89, 010302 (2014).
• Ahlrichs and Dünweg (1999) P. Ahlrichs and B. Dünweg, Journal of Chemical Physics 111, 8225 (1999).
• Harris (2013) E. Harris, The Chlamydomonas Sourcebook: A Comprehensive Guide to Biology and Laboratory Use (Elsevier Science, 2013).
• Kantsler et al. (2013) V. Kantsler, J. Dunkel, M. Polin,  and R. E. Goldstein, Proceedings of the National Academy of Sciences 110, 1187 (2013).
• Röhm et al. (2014) D. Röhm, S. Kesselheim,  and A. Arnold, Soft Matter 10, 5503 (2014).
• de Graaf and Stenhammar (2017a) J. de Graaf and J. Stenhammar, Physical Review E 95, 023302 (2017a).
• de Graaf and Stenhammar (2017b) J. de Graaf and J. Stenhammar, Personal Communication (2017b).
• Chatterji and Horbach (2005) A. Chatterji and J. Horbach, Journal of Chemical Physics 122, 184903 (2005).
• Ilse et al. (2016) S. E. Ilse, C. Holm,  and J. de Graaf, The Journal of Chemical Physics 145, 134904 (2016).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters