The Monte Carlo method was originally developed by von Neumann, Ulam, and Metropolis to study the diffusion of neutrons in fissionable material on the Manhattan Project9,10 and was first applied to simulate radiation damage of metals more than 40 years ago by Besco,11 Doran,12 and later Heinisch and coworkers.13,14

Monte Carlo utilizes random numbers to select from probability distributions and generate atomic configurations in a stochastic process,15 rather than the deterministic manner of MD simulations. While different Monte Carlo applications are used in com­putational materials science, we shall focus our atten­tion on KMC simulation as applied to the study of radiation damage.

The KMC methods used in radiation damage studies represent a subset of Monte Carlo (MC) methods that can be classified as rejection-free, in contrast with the more classical MC methods based on the Metropolis algorithm.9,10 They provide a solu­tion to the Master Equation which describes a physi­cal system whose evolution is governed by a known set of transition rates between possible states.16 The solution proceeds by choosing randomly among the various possible transitions and accepting them on the basis of probabilities determined from the corresponding transition rates. These probabilities are calculated for physical transition mechanisms as Boltzmann factor frequencies, and the events take place according to their probabilities leading to an evolution of the microstructure. The main ingredients of such models are thus a set of objects (which can resolve to the atomic scale as atoms or point defects) and a set of reactions or (rules) that describe the manner in which these objects undergo diffusion, emission, and reaction, and their rates of occurrence.

Many of the KMC techniques are based on the residence time algorithm (RTA) derived 50 years ago by Young and Elcock17 to model vacancy diffusion in ordered alloys. Its basic recipe involves the following: for a system in a given state, instead of making a number of unsuccessful attempts to perform a transi­tion to reach another state, as in the case of the Metropolis algorithm,9, 0 the average time during which the system remains in its state is calculated. A transition to a different state is then performed on the basis ofthe relative weights determined among all possible transitions, which also determine the time increment associated with the selected transition. According to standard transition state theory (see for instance Eyring1 ) the frequency Гx of a ther­mally activated event x, such as a vacancy jump in an alloy or the jump of a void can be expressed as:

Гх =exp(-йг) [11

where nx is the attempt frequency, kB is Boltzmann’s constant, Г is the absolute temperature, and Ea is the activation energy of the jump.

During the course of a KMC simulation, the probabilities of all possible transitions are calculated and one event is chosen at each time-step by extract­ing a random number and comparing it to the relative probability. The associated time-step length dt and average time-step length At are given by:

= and At = 1 [2]

Гх Г x И

n n

where r is a random number between 0 and 1. The RTA is also known as the BKL (Bortz, Kalos, Lebowitz) algorithm.19 Other techniques are possible, as described by Chatterjee and Vlachos.20 The basic steps in a KMC simulation can be summarized thus:

1. Calculate the probability (rate) for a given event to occur.

2. Sum the probabilities of all events to obtain a cumulative distribution function.

3. Generate a random number to select an event from all possible events.

4. Increase the simulation time on the basis of the inverse sum of the rates of all possible events

At — , where — is a random deviate that

N, R,

assures a Poisson distribution in time-steps and N and R are the number and rate of each event i.

5. Perform the selected event and all spontaneous events as a result of the event performed.

6. Repeat Steps 1-4 until the desired simulation condition is reached.