.. _tutorial: Tutorial ======== This tutorial introduces the reader to particle-based simulations of nanomotors. The main simulation method consists of a coupled scheme of Multiparticle Collision Dynamics (MPCD), for the fluid, and Molecular Dynamics (MD), for the motor. Chemical activity is introduced in the fluid via a surface-induced catalytic effect and bulk kinetics. To start working with RMPCMD, visit the :ref:`install` section first. Introduction ------------ Nano- to micro-meter scale devices that propel themselves in solution are built and studied since about a decade. They represent a promise of future applications at scales where usual control strategies reach their limits and, ideally, autonomous action replaces manual control. It is possible to compute, via phoretic theory, the stationary regime of operation of nanomotors in simple geometries. Still, testing geometrical effects or including fluctuating behaviour is best done using numerical simulations. A successful modeling strategy was started by Rückner and Kapral in 2007 :cite:`ruckner_kapral_prl_2007`. It builds on a particle-based fluid, explicitly the Multiparticle Collision Dynamics (MPCD) algorithm. The flexibility of particle-based simulations allowed for numerous extensions of their work to Janus particles, polymer nanomotors, various chemical kinetics, thermally active motors, among others. The principle of the hybrid scheme is very close to full Molecular Dynamics (MD), with the major difference that solvent-solvent interactions are not explicitly computed and are replaced by cell-wise collisions at fixed time intervals. This saves computational time and renders otherwise untractable problems feasible. In Ref. :cite:`ruckner_kapral_prl_2007`, the authors introduce a computational model for a dimer nanomotor that is convenient thanks to its simple geometry. There are two spheres making up the dimer, linked by a rigid bond, one of which being chemically active and the other not. Solvent particles in contact with the chemically active sphere are converted from product to fuel. The active sphere thus acts as a sink for reagent particles (the "fuel") and a source for product particles. Preliminary remarks ------------------- This tutorial does not intend to cover all *possible* manners to conduct nanomotor simulations. Rather, it aims at presenting one strategy for modeling chemically powered nanomotors, that is the combination of a chemically active MPCD fluid coupled to possibly catalytic colloidal beads. While this tutorial relies on the RMPCDMD software, it is worth mentioning that Peter Colberg has also made his software for dimer nanomotor simulations available openly :cite:`colberg_nanodimer_web`. ``nano-dimer`` is based on OpenCL and can benefit from GPU acceleration. Literature ---------- The MPCD algorithm was introduced in :cite:`malevanets_kapral_mpcd_1999` and :cite:`malevanets_kapral_mpcd_2000`. General reviews on the MPCD simulation method are available in the literature. - Raymond Kapral, *Multiparticle collision dynamics: simulation of complex systems on mesoscales*, Adv. Chem. Phys. **140**, 89 (2008). :cite:`kapral_adv_chem_phys_2008` - G. Gompper, T. Ihle, D. M. Kroll and R. G. Winkler, *Multi-Particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids*, Adv. Polymer Sci. **221**, 1 (2008). :cite:`gompper_et_al_adv_polym_sci_2008` An overview of chemically powered synthetic nanomotors has been published by Kapral :cite:`kapral_perspective_jcp_2013`. There is no literature on the practical conduct of nanomotor simulations, however. The MPCD fluid and Molecular Dynamics ------------------------------------- A MPCD fluid consists of point particles with a mass (set to unity here for convenience), a position :math:`x` and a velocity :math:`v`. The particles evolve in two step: (i) indepedent streaming of the particles for a duration :math:`\tau` and (ii) cell-wise collision of the particles’ velocities. For particle :math:`i` this results in the following equations: .. math:: :label: stream x_i' = x_i + v_i \tau and .. math:: :label: collide v_i' = v_\xi + \omega_\xi ( v_i - v_\xi ) where the prime denotes the quantities after the corresponding step, :math:`\xi` is a cell, :math:`\omega_\xi` is a rotation operator and :math:`v_\xi` is the center-of-mass velocity in the cell. The cell consists in a regular lattice of cubic cells in space. Equations :eq:`stream` and :eq:`collide` conserve mass, energy and linear momentum. The viscosity for a MPCD fluid can be computed from its microscopic properties: .. math:: \eta = \frac{k_BT\tau\rho}{2m} \left( \frac{ 5\gamma - (\gamma-1+e^{-\gamma})(2-\cos\alpha-\cos 2\alpha) }{(\gamma - 1 + e^{-\gamma})(2-\cos\alpha-\cos 2\alpha)} \right) + \frac{m}{18 a \tau} (\gamma -1 + e^{-\gamma})(1-\cos\alpha) One can embed a body in a MPCD fluid by using a explicit potential energy. Then, the streaming step is replaced by the velocity-Verlet integration scheme. Collision involve only fluid particles and not the colloid. The dimer nanomotor ------------------- Physical setup ^^^^^^^^^^^^^^ In this section, we review the propulsion of the dimer nanomotor presented by Rückner and Kapral. The geometry of the motor and the chemical kinetics are presented in :numref:`simpledimer`. The solvent consists of particles of types A and B, initially all particles are set to A (the fuel). Fuel particles that enter the interaction range of the catalytic sphere are flagged for reaction but the actual change of A to B only occurs when the solvent particle is outside of any interaction range. Else, the change would generate a discontinuous jump the in the potential energy and disrupt the trajectory. This chemical activity generates an excess of product particles "B" around the catalytic sphere and a gradient of solvent concentration is established. .. _simpledimer: .. figure:: simple_dimer.png Geometry and chemistry for the dimer nanomotor. The graph sketched below represents the local excess of “B” particles that is asymmetric for the “N” sphere. Many more “A” and “B” particles not shown. In this type of simulation, the total energy is conserved but the system is maintained in nonequilibrium by *refueling*, that is by changing B particles to species A when they are far enough from the colloid. The solvent and colloids interact via a purely repulsive Lennard-Jones potential of the form .. math:: V(r) = 4 \epsilon \left( \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} - 1 \right) where :math:`\epsilon` and :math:`\sigma` can be different depending on the combination of solvent and colloid species. Simulation setup ^^^^^^^^^^^^^^^^ Within RMPCDMD, the simulation program for the dimer is called ``single_dimer_pbc``. This program requires a configuration file that contains the physical parameters, an example of which is given in the listing below. :: # physical parameters T = .16666666 L = 32 32 32 rho = 9 tau = 1.0 probability = 1.0 bulk_rmpcd = F bulk_rate = 0.001 # simulation parameters N_MD = 200 N_loop = 50 equilibration_loops = 50 colloid_sampling = 50 # interaction parameters sigma_N = 4.0 sigma_C = 2.0 d = 6.8 epsilon_N = 1.0 0.1 epsilon_C = 1.0 1.0 epsilon_N_N = 1.0 epsilon_N_C = 1.0 epsilon_C_C = 1.0 The configuration allows one to set the size of both spheres in the dimer as well as the interaction parameters. The setting ``epsilon_N`` contain the prefactor to the Lennard-Jones potential for the "N" sphere and all solvent species on a single line. In this example, all the interaction parameters are set to 1 except for the interaction between the "N" sphere and "B" solvent particles, as was done in :cite:`ruckner_kapral_prl_2007`. Running the simulations ^^^^^^^^^^^^^^^^^^^^^^^ .. note:: Make sure that you have built the code properly (see :ref:`install`) and that the command-line tool ``rmpcdmd`` is available at your command-line prompt. You will also need a working scientific Python environment (see :ref:`install_python`). An example simulation setup is provided in the directory ``experiments`` of RMPCDMD. There, the sub-directory ``01-single-dimer`` contains a parameter file. Review the parameters in the file ``dimer.parameters`` then execute the code .. code:: bash make dimer.h5 The actual commands that are executed will be shown in the terminal. Analyzing the data ^^^^^^^^^^^^^^^^^^ The output of the simulation is stored in the file ``dimer.h5``, that follows the H5MD convention for storing molecular data :cite:`h5md_cpc_2014`. H5MD files are regular HDF5 files and can be inspected using the programs distributed by the HDF Group. Issue the following command and observe the output: .. code:: bash h5ls dimer.h5 HDF5 files have an internal directory-like structure. In ``dimer.h5`` you should find :: fields Group h5md Group observables Group parameters Group particles Group timers Group The elements are called "groups" in HDF5 terminology. Here, there is data about the particles (positions, velocities, etc), observables (e.g. temperature) and fields (here, the histogram of "B" particles). The ``h5md`` group contains metadata (simulation creator, H5MD version, etc.), the ``timers`` group contains timing data that is collected during the simulation and ``parameters`` contains all the parameters with which the simulation was run. The command .. code:: bash h5ls -r dimer.h5 will visit all groups recursively. The output is then rather large. Let us focus first on the velocity of the dimer, it is located at ``/particles/dimer/velocity``, where it is stored in ``value`` and the time step information of the dataset is stored in ``step`` and ``time``. In the present case, the velocity is sampled at regular time interval of 100 timesteps or equivalently 1 in units of :math:`\tau`. All the data analysis in this tutorial is done using the Python language and a set of libraries: NumPy for storing and computing with array data, h5py for reading HDF5 files, matplotlib for plotting and SciPy for some numerical routines. For installation, see appendix [install-py]. Some generic programs are provided with as an introduction to reading the files, such as ``h5md_plot.py``. Its usage is .. code:: bash rmpcdmd plot dimer.h5 --obs temperature (the ``obs`` option is preceded by two dashes) to display the temperature in the course of time. This program can also display the trajectory of the dimer .. code:: bash rmpcdmd plot dimer.h5 --traj dimer/position More specific information on the dimer nanomotor can be obtained via Python programs located in the ``experiments/01-single-dimer`` directory. The directed velocity, that is the velocity in the direction of the motor's propulsion axis, can be obtained via .. code:: bash python plot_velocity.py dimer.h5 --directed A further option ``--histogram`` show an histogram instead of the time-dependent value. Finally, an important quantity to assess both for theory and experiments is the effective diffusion that results from both thermal fluctuations and the combination of self-propelled motion and random reorientation. .. code:: bash python plot_msd.py dimer.h5 This latter program can take several simulation files as input to obtain better statistics. It is also important to use a much simulation time (``N_loop``) than the default one to produce meaningful results. The Janus nanomotor ------------------- A Janus motor comprises two hemispheres with one chemically active surface and one inactive surface as shown in :numref:`janusfigure` (a). Because chemical reactions happen asymmetrically on the Janus motor surface, a concentration gradient of product particles is generated giving rise to self-propulsion. In 2013, de Buyl and Kapral introduced a composite model for Janus motor :cite:`de_buyl_kapral_nanoscale_2013`, see :numref:`janusfigure` (b). The active (blue, C) and inactive (red, N) parts are composed of spheres linked by rigid bonds. These spheres have the same radius, and interact with the surrounding solvent particles through repulsive Lennard-Jones potentials :math:`V_{\alpha C}` and :math:`V_{\alpha N}`, where :math:`\alpha = A, B` is the type of solvent species. An example simulation setup for self-proplusive Janus motor is provided in the directory ``experiments/03-single-janus``. The parameter file is show below: :: # physical parameters T = .333333333 L = 32 32 32 rho = 9 tau = 1.0 probability = 1 # simulation parameters N_MD = 100 N_loop = 512 colloid_sampling = 50 equilibration_loops = 20 data_filename = janus_structure.h5 link_treshold = 2.5 do_read_links = F polar_r_max = 10 bulk_rate = 0.1 # interaction parameters sigma_colloid = 2 epsilon_colloid = 2 do_lennard_jones = T do_elastic = F do_rattle = T rattle_pos_tolerance = 1d-8 rattle_vel_tolerance = 1d-8 sigma = 3 epsilon_N = 1.0 5.0 epsilon_C = 1.0 1.0 To run the simulation, use ``make simulation``, and check the propulsion speed :math:`V_z` with .. code:: bash python plot_velocity.py janus.h5 --directed .. _janusfigure: .. figure:: JP.png (a) The sketch of a Janus particle with active (C) and inactive (N) hemispherical surface moving along particle axis with velocity :math:`V_z`. (b) The composite model for Janus particle. Chemical reactions, :math:`A \to B`, take place at the active surface. Controls of motor speed ^^^^^^^^^^^^^^^^^^^^^^^ For Janus particles, one can obtain from phoretic theory that the propulsion speed is determined by factors such as the system temperature, fluid properties (viscosity and density) but also the chemical kinetics and the specific solvent-colloid interactions :cite:`kapral_perspective_jcp_2013`. In this section, we explore the effects of the interaction parameters and of the reaction rate :math:`k_2` on the direction and strength of the propulsion. By modifying only the lines below in the parameter file for the janus simulation, it is possible to observe a reversal of propulsion and the effect of the reaction rate :math:`k_2`. Example 1, forward moving Janus motor. :: epsilon_N = 1.0 0.5 epsilon_C = 1.0 0.5 bulk_rate = 0.001 Example 2, backward moving Janus motor. :: epsilon_N = 0.5 1.0 epsilon_C = 0.5 1.0 bulk_rate = 0.001 Example 3, changing the bulk reaction rate. :: epsilon_N = 1.0 0.5 epsilon_C = 1.0 0.5 bulk_rate = 0.0001