In principle, the implementation of chemical reactions is not a problem in a particle-based method; there exist several well-established stochastic implementation schemes, such as for example the Gillespie-algorithm and its descendents.
Unimolecular reactions are handled as simple stochastic events, in which a particle changes from one type to another. The reaction probability per unit time step is proportional to the macroscopic unimolecular rate coefficient. For bimolecular reactions, we proceed in two stages. Firstly, we describe the situation for spherically symmetric particles. For each reactive particle type, a reactive radius (less than the cutoff radius) is defined. Reactions between two particles take place with a certain probability per unit time step (derived from the rate coefficient) provided that their distance is less than or equal to the sum of their reactive radii.
With the multipolar particles introduced below, orientational constraints on reactions can also be described effectively. This extension also allows association reactions to be modeled directly using mpRDPD, without introducing new particle types for each complex. This is vital for dealing with the combinatorial complexity of supramolecular complexes. The key trick is to use the multipolar potential to correspond to sterically hindered saturation of the number of binding partners (down to just one partner for an appropriate dipole-dipole interaction). Thus, multipolar particles are important not only for dealing with slower, non-diffusion controlled reactions, but also for allowing a correct treatment of saturable association-dissociation reactions.