(264h) A Rejection-Free Markov Chain Monte Carlo to Simulate the Dynamics of Molecular Systems

Monte Carlo methods are today well established as numerical tools for unraveling the complex behavior of a large variety of systems.^[1-4] The broad class of MC algorithms can be divided into two main groups: methods built-up on the Metropolis scheme (MMC) and kinetic Monte Carlo methods (KMC).

The latter apply to systems evolving dynamically from state to state. As such, they can be used when the number of states is discrete and the rate constants for the state-to-state transitions are known.^[5] They are therefore not suitable for molecular systems evolving in a continuous space characterized by an infinite and unknown number of states. In such cases one typically has to revert to MMC.

In MMC the system evolves through randomly sampled configurations which are accepted or rejected according to specific rules which ensure that the distribution of sampled configurations satisfies the Boltzmann equiprobability principle. The possibility to let the system evolve through unphysical trajectories makes MMC algorithms particularly efficient for calculating equilibrium properties. On the other hand, the price to pay in MMC is the lack of a physically meaningful time scale. To date, in fact, with the exception of a few cases (e.g. Brownian and spin dynamics), MMC methods do not provide information about the true time evolution of molecular systems.^[3]

In this contribution a new MC algorithm for molecular systems is proposed that overcomes these limitations and allows simulating the correct dynamics of molecular systems too.

It is first shown that the equations of motion can be written equivalently as a set of equations describing a chemical reaction network, i.e., one can identify the current state of the system with the paths traveled by the atoms. This identification allows to use a KMC scheme, namely the stochastic simulation algorithm,^[6] where the state-to state transitions are determined on-the-fly and a direct link between the MC time step and the true physical time is established by using propensities proportional to the velocities. The method is applied to the classical problem of hard disks and validated against event driven Molecular Dynamics. The position probabilities as well as the relaxation time of the global orientational order parameter and of the velocities converge to the correct results for decreasing step size. Additionally, the Maxwell velocity distribution is recovered directly without any need of independent samplings, as instead in traditional MMC.^[7]

In perspective, the new algorithm can provide a unified numerical framework to simulate the dynamics of systems characterized by a competition between conformational relaxation and chemical reactions, as, for example, in organic light emitting diodes and enzyme catalysis.^[8,9]

References

[1] Potestio et al., Physical Review Letters, 111, 060601 (2013).

[2] Meimaroglou D. and Kiparissides C., Industrial & Engineering Chemistry Research, 53, 8963-8979       (2014).

[3] Landau D.P. and Binder K., Monte-Carlo Simulations in Statistical Physics, Cambdridge University       Press (2009).

[4] Zuniga et al., Physical Review Letters, 114, 155301 (2015).

[5] Voter A.F., Introduction to Kinetic Monte Carlo Method, in Radiation Effects in Solids, Springer            (2005).

[6] Gillespie D.T., Journal of Computational Physics, 22, 403-434 (1976).

[7] Costa L.I., manuscript in preparation.

[8] Daub J. et al., The Journal of Physical Chemistry A, 105, 5655-5665 (2001).

[9] Wolf-Watz M. et al., Nature Structural & Molecular Biology, 11, 945-949 (2004).