Subject: CCL: Re: Diffusion Equation from Kinetic Lattice Monte Carlo

From: Harold W Schranz <[email-masked]>
Date: Thu, 16 Sep 1999 11:03:39 +1000

Hi Andrew,

A brief answer (ignoring some details of notation and implementation):

(1) Time dependence of probability
 ---------------------------------
Andrew Horsfield wrote:
> I have a question about kinetic Monte Carlo simulations on a lattice, and
> how they relate to traditional diffusion and reaction equations.

The relationship between chemical kinetics, diffusion and Monte Carlo
is indeed a close one as can be observed from some of the treatments
in some physical chemistry textbooks (e.g. some editions of Atkins).

> I want to derive continuum equations from this algorithm. To do this I
> begin by defining the probability that my system is in configuration i to
> be P_i. I then ask how P_i varies with time. The answer (I think) is:
> 
> d P_i
> ----- = \sum_j L_{ij} P_j
>  dt
>
> where
> 
> L_{ij} = W_{ij}             i \= j
>        = -\sum_k W_{ki}     i  = j

The dP_i/dt equation is linear in the probabilities P_i unless the transition
rate matrix L_ij also depends on P_i. At the moment it is quite general
as you havent specified what the configurations  described by P_i are.
If L_ij is independent of P_i or P_j then you have a linear system.

To make this more concrete (there are other examples), this set of equations,
for the case that L_ij is independent of P_i or P_j,
is analogous to what is found in the description of unimolecular reactions
by a master equation. If the probability is conserved it time you could
be describing a closed isomerization system e.g. A->B, B->A or an energy
transfer system e.g. A(E) -> A(E') where the indices either tag different
species or energies of each molecule. In the former case, L_ij is related to
the unimolecular rate constant for isomerization from configuration j to
configuration i. In the latter case, L_ij is related to the transition
probabilty matrix for energy transfer from a reactant molecule initially
with energy E_j which ends up at a final energy E_i.

One could think of similar interpretation in terms of hopping rates on a
lattice, L_ij determining the rate of hopping from site j to site i.
 
(2) Time dependence of concentrations
-------------------------------------
> This all looks fine until one transforms this equation into a
> corresponding equation for the rate of change of concentration (C):
> 
> d C(r)_n
> -------- = \sum_{sm} R_{rn,sm} C(s)_m
>    dt
> 
> where r and s are indices for the type of atom (C, O, H or whatever), n
> and m are indices for spatial position (lattice site), and R_{rn,sm} is a
> concentration independent rate matrix that can be derived from L_{ij}. The
> feature that worries me about this is that it is *linear* in C (as it has
> to be, given the starting point).
> 
> I believe that this cannot be correct if there are chemical reactions
> taking place, since then we would expect the rate to depend on the
> *product* of two concentrations (at least).

Not true!!!

This misconception seems to be a consequence of people being taught
kinetics from the bimolecular viewpoint only (A+B->C+D) and most texts
disregard the huge literature on unimolecular and related reactions.
(Probably since the mathematical treatment is rather more complex
than the trivial level required for simple bimolecular reactions!)

Unimolecular reactions have a linear rate law:

d C(t) / dt = - k C(t)

Of course in most cases, the unimolecular rate constant k  depends on
a lot of factors, which may include concentrations of other species which
might be in excess.

> Can anyone tell me where I have gone wrong? Is there a standard reference
> for this kind of approximation?

You havent gone wrong except in interpretation. If you are indeed describing
bimolecular reactions on a lattice (or whatever) then your matrices
L_ij and R_rn,sm will have to have a dependence on concentration as well.

For example:
Basic chemical kinetics: if A + B -> C + D  is an elementary reaction
it has a rate law:

Rate = -dC_A/dt = -dC_B/dt = dC_C/dt = dC_D/dt = k C_A C_B

which of course be written in a linear form:

Rate = -d C_A/dt = k' C_A where k' = k C_B 

but that is likely to only be useful if C_A is constant or slowly varying
in time.

I hope the above deliberations help.

Regards,
Dr. Harold W. Schranz
Visiting Fellow, RSC, ANU
Fujitsu Area 3 Project, ANUSF

-- 
------------------------------------------------------------------------
Dr. Harold W. Schranz,            Office Ph.:         +61 (02) 6249 5988
Computational Chemist,            Dept Ph.:           +61 (02) 6249 3437
ANU Supercomputer Facility,       Fax:                +61 (02) 6279 8199
Australian National University,   Email:       [email-masked]
Canberra, ACT 0200, AUSTRALIA.    WWW:   http://anusf.anu.edu.au/~hws900
------------------------------------------------------------------------