RMPCDMD: Simulations of colloids with coarse-grained hydrodynamics, chemical reactions and external fields

The RMPCDMD software package performs hybrid Molecular Dynamics simulations, coupling Multiparticle Collision Dynamics to model the solvent and Molecular Dynamics to model suspended colloids, including hydrodynamics, thermal fluctuations, and chemically active solvent particles and catalytic colloids. The main usage of RMPCDMD is the simulation of chemically powered nanomotors, but other setups are considered: colloids in the presence of a thermal gradients or forced flows. RMPCDMD is developed in Fortran 2008 with OpenMP for multithreaded operation and uses the HDF5-based H5MD file format for storing data. RMPCDMD comes with documentation and a tutorial for the simulation of chemically powered nanomotors.

RMPCDMD is implemented in Fortran, following the 2008 version of the standard, and uses OpenMP for multi-core execution. Trajectories are stored in the H5MD file format [9], and Python programs are provided to illustrate data analysis. The module fortran h5md [10] provides helper routines to write H5MD files from Fortran. H5MD files are HDF5 files (https://www.hdfgroup.org/HDF5/) with an internal organization tailored for molecular simulations [9].
In designing RMPCDMD, we tried to adopt practices for scientific computing that help to deliver reproducible results and facilitate the extension of the code and its maintenance. We also made special efforts on the usability of the software. The design choices are the results of the authors' experience on scientific programming, with further motivation by the article on "Best Practices for Scientific Computing" by Wilson et al [11]: • The algorithms are implemented in Fortran modules in the src/ directory and shared between the simulations found in programs/. Table I summarizes their content.
• Subroutines and functions are limited in size, so that understanding any single code unit remains reasonable.
• Data is encapsulated in derived types.
• No global variable is used, to clarify to the reader what data is passed to subroutines.
• The programming style strives for legibility, conciseness, and efficiency.
• All simulation parameters are given from a text configuration file, there is no hard-coded setting in the code that would require re-compilation between runs.
• The build procedure and partial testing are automatically executed on the Travis CI platform (https://travis-ci.org/) for continuous integration.

Module name Functionality
common Routines of general use: histogramming, parsing of command-line arguments, timing and minimumimage convention function.
hilbert Computation of the compact Hilbert index of a lattice cell and vice-versa.
cell system Cell data structure.
particle system Particles data structure and sorting routine.
particle system io Data structures to facilitate the storage of thermodynamic and particle data. interaction Computation of Lennard-Jones force and energy.
neighbor list Data structure and routines to update the neighbor list.
mpcd Algorithms for MPCD streaming and collision, and for bulk reactions.  Instead of building a generic program that could handle the logic of different geometries, boundary conditions, and colloid types, specific simulation programs have been implemented. Each program calls the algorithms from the modules and resembles a statically compiled "simulation script". The amount of administrative code (reading parameters, defining storage, etc.) is however larger than what could be achieved in a higher-level language such as Python.
CMake is used to build RMPCDMD and run automated tests in test/ that cover some of the base algorithms such as particle sorting or the histogram routine. The Fortran module fortran tester [13] provide convenience routines for testing (e.g. equality assertions). Automated testing using actual simulations is difficult and very much time intensive, and we rely on it only for a small part of the code. In simulations without a thermostat or wall interactions, the total momentum and energy of the system must be conserved, which is verified for a short run of the single dimer pbc simulation as part of the continuous integration testing. Further validation of the code is done via comparison of physical results: the velocity distribution of particles in equilibrium, the viscosity of the fluid, etc.

User perspective
For the end-user, the interface comes as a single command-line tool, rmpcdmd, that drives the simulations and provides access to utility programs. A typical command is rmpcdmd run single_dimer_pbc dimer.parameters dimer.h5 auto where single dimer pbc is the simulation program, dimer.parameters is the configuration file, dimer.h5 the output datafile, and auto is the setting for the random-number seed. While the individual programs can be run directly, the command-line tool rmpcdmd aims to provide a single file to place in the program "PATH" on the user's computer and to automate some actions. The commands available are: run: execute one of the simulation programs. The setting of the environment variable OMP NUM THREADS, and the start and end time of the run are shown. Bash's time provides the "user time" actually used by the program. seeder: print out a signed 64-bit integer value generated by the computer's /dev/urandom device. plot: display an observable from a simulation file. timers: print out (or display) the timing information from a simulation file.
A benefit of using separate simulation programs is that the configuration file for a simulation only contains directly relevant parameters. All simulation programs share the same command-line syntax, read parameters from a text file using the ParseText Fortran library [14], read the seed for the RNG from the command line and output the trajectory to a H5MD file. An example configuration file for single dimer pbc is listed in Appendix A.
A generic plotting utility is provided and can be invoked as rmpcdmd plot dimer.h5 --obs temperature to plot the temperature, for instance. This plotting program and the other analysis programs that come with RMPCDMD are written in Python and rely on NumPy for storing and processing numerical data, SciPy [15] for numerical algorithms, Matplotlib [16] to produce 2D figures, h5py [17] to read HDF5 files and Mayavi [18] to visualize 3D data.

C. Documentation
A documentation is found in the directory doc/. It is published on the web, http://lab.pdebuyl.be/rmpcdmd/ and also serves as the homepage of RMPCDMD. The documentation is built with the Sphinx documentation generator [19] and is written in reStructuredText. An extension to Sphinx is bundled with RMPCDMD to provide (i) automatic links from the documentation to the source code documentation that is annotated with Doxygen (http://www.stack.nl/~dimitri/doxygen/) and (ii) automatic inclusion of the Doxygen header of Fortran programs. Specific instructions for the installation and execution of RMPCDMD are provided, as well as general informations on the algorithms are presented. The description of the simulation programs is generated automatically from the beginning of the corresponding .f90 files and includes all parameters for the simulations.
A tutorial on the simulation of chemically powered nanomotors, based on RMPCMD, was given at the recent workshop "Modeling of chemically powered nanomotors" held at the KU Leuven in April 2016 (http://lab.pdebuyl.be/2016-nanomot This tutorial is reproduced in the documentation of RMPCDMD and represents a unique resource in the field of nanomotor simulation.

D. Algorithms
For completeness, we state here exactly what algorithms are implemented in RMPCDMD, with the corresponding reference to the literature.
RMPCDMD implements the original MPCD collision rule [2] and the Anderson thermostat [20]. Random shifting of the spatial grid is always performed [21] to ensure Galilean invariance. Flows are generated by applying a constant acceleration to solvent particles [22,23]. Stick boundary conditions at the wall use the bounce-back rule and virtual particles during the collision step [24].
Molecular Dynamics is performed with the velocity-Verlet algorithm [1,3,8]. The interaction potential, both for solvent-colloid and colloid-colloid interactions, are purely repulsive cut-off 12-6 Lennard-Jones potentials of the form V (r) = 4ǫ σ r 12 − σ r 6 + 1 4 for r < 2 1/6 σ where the energy parameter ǫ and the scale parameter σ depend on the species considered. Verlet and cell lists are used to compute the solvent-colloid forces and energies [1]. The sorting of solvent particles in cells is the same as for the MPCD collision step. The rigid body dynamics uses RATTLE [25] in its iterative implementation for the Janus particle and exactly for the dimer motor. Spatial sorting of the solvent particles is based on compact Hilbert indices [26]. Random numbers are generated with the Threefry 2x64 RNG [27], allowing each OpenMP thread to generate random numbers independently. The RNG is provided by the random module library [28], where it is implemented in C99. The choices are inspired by the nano-dimer [29] simulation code, which implements these algorithms in OpenCL C for large-scale parallel processors.
To obtain a consistent measure of the temperature T even in the presence of flows, we compute T by averaging the center-of-mass reference frame temperature over all MPCD cells as where N c is the number of cells, the variable ξ represents a cell, N ξ the number of particles in a cell, v ξ the centerof-mass velocity of the cell, and m i and v i the mass and velocity of particle.

E. Applications of RMPCDMD
The prototypical use of the MPCD algorithm to study chemically active particles is the dimer nanomotor, whose model was introduced by Rückner and Kapral [5]. It consists of two rigidly linked spheres evolving in a solvent. Chemical reactions are catalyzed by one of the spheres and create a local concentration gradient, effectively propelling the dimer along its symmetry axis. This simulation is implemented in the program single dimer pbc.
To facilitate usage of RMPCDMD by newcomers, ready-made simulation "experiments" are provided in experiments/. There, makefiles allow the user to simply type make simulation to start a simulation. The terminology that we have chosen reflects our opinion that computer simulations can be viewed as numerical experiments in which one prepares a setup and conducts a study to understand a certain physical phenomenon. The ability to modify a given experimental choice (i.e. thermostating or not, surface properties, etc.) to assess its effect is therefore critical. By providing several chemical models and collision rules, RMPCDMD provides such a framework.
Four programs are provided to analyze the data of this experiment. plot velocity.py displays the laboratory reference frame velocity or the directed velocity of the dimer (with the option --directed). plot histogram.py displays the cylindrical shell histogram of the reaction product concentration. plot msd.py displays the meansquared displacement of the dimer's center-of-mass position. This quantity is a practical interest in comparison to experimental studies [30]. view last simulation frame.py displays a 3D representation of the dimer motor using the Mayavi visualization library [18], either with its full trajectory as a line (with option --unwrap) or with the product solvent particles B (with option --show-B). The two first analysis programs are modeled after the results by Rückner and Kapral [5].
While the aims of single dimer pbc are to reproduce Ref. [5] and to provide a starting point for other developments, the other simulation programs in RMPCDMD are related to forward-looking research projects that will benefit from a robust codebase. These other setups are: poiseuille flow generates a Poiseuille flow between parallel plates, with no colloid. chemotactic cell generates a Poiseuille flow with a two-species chemical gradient perpendicular to the flow and an embedded spherical colloid or a dimer nanomotor. single janus pbc implements a composite Janus nanomotor [31]. single sphere thermo trap generates a thermal gradient between two plates, with an embedded spherical colloid.

II. AVAILABILITY
Operating System: RMPCDMD has been tested on Linux using the gcc/gfortran and icc/ifort compilers and on OS X (El Capitan) using gcc/gfortran (installed via MacPorts).
Programming Language: Fortran 2008 for the simulation code, C99 for the Random Number Generator, bash and make for the execution and Python (with NumPy, SciPy, matplotlib and h5py) for analysis. Besides nano-dimer [29] (by Peter Colberg, available under the MIT license), RMPCDMD is the only opensource software to provide support for chemically-powered nanomotor simulation. RMPCDMD has, in addition to nano-dimer, support for external fields or rigid-body dynamics for assemblies larger than the dimer. nano-dimer, on the other hand, is based on OpenCL C and on the Message Passing Interface for distributed parallel execution, and is able to handle many-motors simulations that are out of reach of RMPCDMD for performance reasons. Given the research activity on the topic of active colloids, we believe that providing a reference implementation of a chemically active MPCD fluid can serve as a starting point for further modeling, by us or by other groups. RMPCMD has not been used for published work yet. It is in use for ongoing research projects at the Instituut voor Theoretische Fysica of the KU Leuven.