Molecular Dynamics (MD) is a computational technique to study condensed matter systems by numerically solving Newton’s equations of motion for a large number of particles . Often, due to the number of constituents in a system, so-called coarse-grained or multiscale strategies are used to reduce the computational cost of a simulation.
One simulation strategy for modelling colloids embedded in a solvent is to couple the Multiparticle Collision Dynamics (MPCD) algorithm  for the solvent with MD for the colloids . It differs from traditional MD by the use of a coarse-grained representation of the solvent’s structure and dynamics and by the possible occurrence of chemical species transformations in the course of the simulation. In the latter case, the extension Reactive MPCD (RMPCD) provides bulk chemical kinetics that reduce to reaction–diffusion equations . The name RMPCDMD results from the concatenation of the different algorithms’ acronyms: Reactive Multiparticle Collision Dynamics and Molecular Dynamics.
The focus of RMPCDMD lies on the modeling of chemically active colloids (e.g. nanomotors ) or of colloids subjected to external gradients (e.g. flows  or temperature gradients ). A RMPCDMD simulation consists in Nsolvent solvent particles and Ncolloids colloidal particles with a mass, position, and velocity, placed in a cuboid domain. Colloid–solvent and colloid–colloid interactions follow Lennard–Jones type potentials, but there is no interaction potential acting among the solvent particles. As Nsolvent >> Ncolloids, this considerably speeds up a simulation in comparison to a full MD simulation. The time evolution of colloidal and solvent particles follows Newton’s equations of motion, which are solved iteratively using the velocity-Verlet algorithm [1, 8]. The output from a simulation is a trajectory representing the coordinates of all particles in the course of time. Typically, solvent data is preaveraged during the simulation as its storage requires significant amounts of space.
Implementation and architecture
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 , and Python programs are provided to illustrate data analysis. The module fortran_h5md  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 .
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 :
- The algorithms are implemented in Fortran modules in the src/ directory and shared between the simulations found in programs/. Table 1 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 code is tracked with the git revision control software  and published on GitHub (https://github.com/pdebuyl-lab/RMPCDMD).
- The build procedure and partial testing are automatically executed on the Travis CI platform (https://travis-ci.org/pdebuyl-lab/RMPCDMD) for continuous integration.
|common||Routines of general use: histogramming, parsing of command-line arguments, timing and minimum-image 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.|
|md||Algorithms for velocity-Verlet and rigid body MD.|
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 test programs in test/ that cover some of the base algorithms such as particle sorting or the histogram routine. For spatial sorting, for instance, the compact Hilbert index of the particles is re-computed after sorting to check the consistency of the sorted list. The Fortran module fortran_tester  provide convenience routines for testing (e.g. equality assertions). All programs in test/ whose filename starts with test_ are executed automatically by Travis CI when the code repository is updated. Users can also execute the make test command on their machine after having compiled RMPCDMD.
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. This is verified for a short run of the single_dimer_pbc simulation as part of the continuous integration testing: A Python program test/dimer_test.py reads the simulation output and validates, to a set accuracy, these conservation laws. This also insures that the common foundation of RMPCDMD, such as the force computation, neighbor listing, particle sorting, etc., fulfil their role.
Further validation of the code is done via comparison of physical results. For pure solvent simulations in equilibrium, one can verify that the velocity distribution of the solvent particles is a Boltzmann distribution at the set temperature. The Poiseuille flow simulations allows one to verify the viscosity, via the velocity profile. For colloids, one can check the diffusion coefficient against the literature results . These validations must be performed by the user. The presence of analysis programs in Python in the experiments/ directory gives a good starting point for this work.
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 , read the seed for the RNG from the command line and output the trajectory to a H5MD file.
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  for numerical algorithms, Matplotlib  to produce 2D figures, h5py  to read HDF5 files and Mayavi  to visualize 3D data.
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  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-nanomotor/). This tutorial is reproduced in the documentation of RMPCDMD and represents a unique resource in the field of nanomotor simulation.
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  and the Anderson thermostat . Random shifting of the spatial grid is always performed  to ensure Galilean invariance. Flows are generated by applying a constant acceleration to solvent particles [23, 24]. Stick boundary conditions at the wall use the bounce-back rule and virtual particles during the collision step .
Molecular Dynamics is performed with the velocity-Verlet algorithm [1, 3, 8]. The interaction potentials, both for solvent–colloid and colloid–colloid interactions, are purely repulsive cut-off 12–6 Lennard-Jones potentials of the form
where the energy parameter ε and the length scale parameter σ depend on the species considered. Verlet and cell lists are used to compute the solvent–colloid forces and energies . The sorting of solvent particles in cells is the same as for the MPCD collision step. The rigid body dynamics uses RATTLE  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 . Random numbers are generated with the Threefry 2x64 RNG , allowing each OpenMP thread to generate random numbers independently. The RNG is provided by the random_module library , where it is implemented in C99. The choices are inspired by the nano-dimer  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 Nc is the number of cells, the variable ξ represents a cell, Nξ the number of particles in a cell, υξ the center-of-mass velocity of the cell, and mi and υi the mass and velocity of particle.
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 . 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
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 the first experiment 01-single-dimer. 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 mean-squared displacement of the dimer’s center-of-mass position. This quantity is a practical interest in comparison to experimental studies . view_last_simulation_frame.py displays a 3D representation of the dimer motor using the Mayavi visualization library , 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 . Their output is shown in Figures 1 and 2 for a simulation of the “experiment” 01-single-dimer.1 The tree-dimensional visualization of view_last_simulation_frame.py is shown in Figure 3. The assessment of the correctness for a simulation, discussed in section quality control, can be started by using the provided plotting tool to observe the conservation of the center of mass velocity (rmpcdmd plot dimer.h5 --obs center_of_mass_velocity) or of the internal energy (rmpcdmd plot dimer.h5 --obs internal_energy).
While the aims of single_dimer_pbc are to reproduce Ref.  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 . The rigid-body RATTLE or elastic network version can be chosen via an option in the configuration file. single_sphere_thermo_trap generates a thermal gradient between two plates, with an embedded spherical colloid.
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).
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.
Additional system requirements
The largest RAM usage in RMPCDMD comes from the solvent coordinates and requires about Nsolvent×250 bytes of storage. The example simulations in the experiments/ directory require between 100MB and 200MB of RAM.
The CMake and Make tools, and the HDF5 library. For data analysis, NumPy, SciPy, matplotlib and h5py.
List of contributors
Pierre de Buyl: Project creator, code, documentation, build system, tutorial.
Laurens Deprez: Code: RATTLE for the dimer, angle-dependent MPCD, boundary conditions and forcing for flows, Lennard-Jones 9-3 walls, chemotactic cell setup.
Mu-Jie Huang: Tutorial: Janus particle.
Peter Colberg: Initial OpenMP parallelization, continuous integration testing and build system.
Persistent identifier: DOI: 10.5281/zenodo.198599
Licence: BSD 3-clause
Publisher: Pierre de Buyl
Version published: 1.0
Date published: 9 December 2016
Persistent identifier: https://github.com/pdebuyl-lab/RMPCDMD
Licence: BSD 3-clause
Date published: 9 December 2016
(3) Reuse Potential
RMPCDMD’s most valuable content consists in the collection of algorithms in src/, many of which are not available widely, and the simulation protocols implemented in programs/. All steps of a simulation appear in the body of the corresponding program, including calls to the force computation, RATTLE, neighbor listing, etc, so that following the full execution of the code is explicit.
There are two directions to modify RMPCDMD. The first one is to add an application, or a virtual experiment, via an additional program in programs/, using the closest existing program as a basis. The second one is to implement a new algorithm or feature in an existing or new module under src/, ideally with dedicated tests in test/. Ideas for further extensions include alternative chemical schemes and geometries, performance tuning for many-motor simulations, more efficient rigid-body implementation, etc.
Given the scarcity of open-source nanomotor simulation code, we briefly compare RMPCDMD to nano-dimer  (by Peter Colberg, available under the MIT license), these two being the only open-source softwares 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 both as a tool to reproduce earlier results from the literature and 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.