Kinetic MC: Vacancy Migration in 2-D Binary Alloy
Answers and Solutions

Question 1 - Flowchart

Simple Flowchart
  1. Read Input
  2. Set up either RANDOM or ORDERED Configurations
  3. Random place Diffusing Vacancy
  4. Calculate initial SRO and LRO
  5. Determine jump rates for moving Vacancy to a nearest-neighbor site (i.e., dtime=( w(1)+w(2)+w(3)+w(4) )-1 )
  6. Choose jump stochastically from "n-fold-way" type method
  7. Update neighbor lists and SRO and LRO values
  8. Perform this over the number of steps or Blocks of data.
Note that time is not chosen as TIME= -ln(R)*( w(1)+w(2)+w(3)+w(4) )-1, where R is a random number, thus conforming to Possion rate. This is because there is no need in this simulation to continually calculate ln(R) because we are always performing just a given set of jumps that do not have different overall rates. We save some calculation time. However, for adsorption/desorption problem we did different types of moves with different rates depending on occupany of surface. The ln(R) must then be included.

SRO and LRO
Long-range order is given as values between 0 and 1, depending the state of order. It was counted like a checker board and any flaws reduced the value from 1. If P_BA is the probability of an A atom at a site given a B at the central site, and P_A is the probability of the occurance of an A atom (or the concentration of A), the usual Warren-Cowley SRO definition, 1 - P_BA/P_A, is:

 

Question 2 - Initialization with Random Configuration

INIT CONF parameter is set to 3 giving a random starting configuration. The 2D alloy is chosen as 50% A-B alloy, so it is like Ising Model. With the VAA=VAB=0 and VBB= -0.05 eV, we know that Tcrit= (2/ln(1 + sqrt(2))/4= 0.5673 eV = 658 K.

By trying various temperatures, such as T=360 K, 560 K, 660 K, 680 K, 760 K and 860 K, we may find the Tcrit as indicated by the SRO and LRO parameters. See kmc_random_out . Recall that this is not best way to find sensitive quantities like Tcrit. It is best to use fluctuations (or derivatives) in quantities like LRO, or other so-called cumulants. But, in any case you must average over several MC runs to obtain proper averages. Nontheless, you get the idea.

Importantly, if you would average over 10 or more configurations for these temperatures, you would obtain very reasonable representation for the SRO and LRO parameters versus Temperature with a curvature change around 650-680 K. See MC average:

N.B. After we ran this SRO vs. T, I corrected a definition of SRO within the code that was missing a 1/c (= 2 at 50%) factor. Therefore, at 50% you will have SRO that goes from 0 to -1 in the present case.

Notice also that more statistics are required above 700 K to get the LRO better.

How does temperature affect the KMC time?
The higher the temperature the faster the processes are and the shorter the time scale should be. Thus, time increments will be shorter as temperature is increased. This is mirrored in the fact that the simulation is modeling a Possion process. As the time, t, in a Possion processes goes like Rate-1 and the Rate scales as e-E/T, plotting ln(t) for several T's (keeping the number of jumps fixed and long enough) should give a linear behavior. For 100 Million vacancy jumps, the final times and various temperatures are given in the kmc_random_out file.

 

Question 3 - Initialization with Ordered Configuration

Starting with an ordered configuration confirms most of what you found from random configuration. However, by starting out ordered and increasing temperature, and by starting out disordered and decreasing temperature, you may sometimes see the transition and the hysterisis if a first-order transition. See kmc_ordered_out .

 

Question 4 - Special Effects: Radiation Damage

For driven systems like those subjected to neutron radiation damage, the atoms will continually be disordered, even if the systems chemically favors an ordered state. At high T, the vacancies also created or existing in the system can diffuse and "heal" the occassional disordering introduced via the neutron damage. As the T is lowered, however, the vacancy can no longer diffuse and there is a critical temperature (determined by the radiation dose rate and the chemical interactions) where an ordered phase can no longer be sustained. Thus, the ordered phase loop in a T vs. c diagram will "pinch off" and there will be order only at a finite temperature range.

I have altered the ordering.f to incorporate one model of athermal ballistic radiation damage; that is, there is a constant ballistic jump rate that is added to the other vacancy jump rates and it happens at random. See ordering_ball.f . Note that the additional input is the ballistic jump rate, w_ball, which according to my scales can be between 0 and 10, say. The output kmc_ballistic_out show that the LRO can be destroyed at low T and depends on w_ball.

You may see how this athermal stochastic effect (the random radiation damage) was incorporated into this code. I assumed that a site was swapped at random with a random neighbor, and either of these could be an A, B, or vacancy type site. The update of the neighbor tables is the nasty part. The TIME was take as ( w(1)+w(2)+w(3)+w(4) + w_ballistic )-1, where w_ballistic was an input. However, there are many scenarios one could explore. For example, just swapping two sites at random (which seems less plausible in some causes like neutron damage). Nonetheless, interesting things can happen: in the case of a phase separating alloy, rather than a ground state of all A on one side and B on the other, a new scale emerges such that the stationary state is domains of A and B interspersed. The scale is dependent on the swap distance of the atoms, apparently. Look for results from Enrique and Bellon (2000 and 2001) for such simulations.

Thanks to Raul Enrique for helping find a bug for this case.

Copyright D.D. Johnson, December 13, 1998