time step in kinetic montecarlo simulation

I have a doubt about how the time step is calculated in kinetic montecarlo simulations.

One state with index $i$ is connected to other $N$ states, indexed by $j=1...N$, by transitions that happen with rates $r_{ij}$ (form $i$ to $j$). At each iteration of the algorithm, one of the $N$ transitions is randomly chosen (with the right probability), and the time step is cleverly calculated using a random number $u$ in $(0,1]$ by

$\Delta t = \frac{ln(1/u)}{R}$

where $R$ is the sum of all the rates of the transitions leaving from state $i$ : $R=\Sigma_{j=1}^N r_{ij}$. In this way $\Delta t$ is exponentially distributed with mean $1/R$.

My doubt: I understand that the time between events (or transitions) $\Delta t$ should be exponentially distributed because the events are Poisson distributed. But why should the mean of the time between events be $1/R$? If the algorithm in a particular iteration picks up the particular transition $i$ to $j_o$, with the rate $r_{ij_o}$, I would say that the average residence time in state $i$ before leaving to state $j_o$ should be $1/r_{ij_o}$, and not $1/R$. So instead of the $\Delta t$ above, I would use in that particular iteration a time step

$\Delta t=\frac{ln(1/u)}{r_{ijo}}$.

Is this wrong?

• I think experts of this can be here, but I'm not 100% sure this fits PhysicsSE. In case can you suggest where to migrate? – scrx2 Jan 31 '16 at 13:55
• More people will probably see it here, but if you don't get an answer here in a few days, it might be worth migrating to Computational Science. – tpg2114 Jan 31 '16 at 14:04
• I'm voting to close this question as off-topic because of the reasons cited by @tpg2114 . – user36790 Jan 31 '16 at 14:54
• 1/R is not the same R as above but the sum resulting from the random transition that happens before. Then it makes a distribution on the cumulative rate from the beginning until the current random event ( and not only on the random event rate ) – user46925 Jan 31 '16 at 16:44
• read directly Monte Carlo and kinetic Monte Carlo methods - a tutorial by Peter Kratzer chapter 3 From MC to kMC: the N-fold way. – user46925 Jan 31 '16 at 19:25

From the wikipedia article you cited

"It is important to note that the timestep involved is a function of the probability that all events j, did not occur." (note they use i in wikipedia but we are using j so I changed the quote to match)

This probability is u which can be constructed from the multiple poisson distributions for each individual event $r_{i,j}$. $$P(r_{i,j}, k) = \lambda^k e^{-\lambda}/k!$$ If no events occur, k=0. $$P(r_{i,j}, 0) = e^{-\lambda}$$ From dimensional analysis, we can determine $\lambda = r_{i,j}\Delta t$

Finally, the probability of no event occurring is the product of all these probabilities (which we call u). $$u = P(r_{i,1})P(r_{i,2})...P(r_{i,N}) = \exp(-\Delta t (r_{i,1} + r_{i,2} + ... + r_{i,N})) = \exp(-\Delta t R)$$ You can easily solve for $\Delta t$ and arrive at the above expression.

So to understand the logic here, let's consider we start at $t=0$. After some time, you cause an event to occur at time $t_1$. Now, between time $t=0$ and $t=t_1$, no other event has occurred which is the exact probability that was just derived.

If you instead used the formula that you are suggestion, you are drawing from the probability that that single event has not occurred in a time $\Delta T$ without considering any other possible event.

In kinetic Monte Carlo, the idea is to describe a trajectory as a set of events, at which the system makes a transition from one state ($i$) to another ($j$). To generate such a trajectory, we need to randomly select both the states that are visited and the intervals between them.

The time interval $t$ between a pair of events is the time in which nothing happened—the system stayed in the current state because none of its possible transitions occurred. So the probability distribution for $t$ involves the total rate for anything happening, which is $R$. After we have determined the interval $t$, we can then decide which event it was that happened.

An example might help: Imagine you have a large number $n$ of radioactive atoms, with a decay rate $\gamma$. The typical time until any given atom decays is $\gamma^{-1}$, but the typical time until the first decay is $(\gamma n)^{-1}$. When you generate a trajectory, this is the mean time you wait until the first event. If $n$ is very large, you will not have to wait long to see a decay, even if $\gamma$ is small.