A Kinetic Monte Carlo Homework
Vacancy-assisted Ordering in 2-D Binary Alloy

Copyright D.D. Johnson, November 30, 1998	Homework Solutions are here.
updated May 30, 2001 for Summer School	No peeking before you do the HW!

GOAL: To simulate the temperature versus composition (T vs. c) phase diagram of a Binary Alloy (A_1-cB_c) on a Square Lattice.

METHOD: Use a 2D lattice Kinetic Monte Carlo code for Binary Alloy that uses a Diffusing Vacancy to assist in reaching Steady-State chemical arrangements in a c-T diagram, even under driving force. (We include the possibility of ballistic jumps which induce local disorder, such as possible under irradition of alloy).

N.B. In principle, we are looking at steady-state solutions which do not require dynamic or explicit time information. However, the vacancy is diffusing to assist in ordering the alloy, and, as such, we use KMC method to do this. (In the Homework, there is a question about this.)

SIMPLIFICATION: For this Lab Exercise, we will focus on c=50% alloys, which, with ordering-type interactions, yields CuAu ordering. That is, the thermal equilibrium ground state is 2D "checker board" arrangement of atoms A and B.

N.B. This code does work for c=(0,1), but then you must be careful that much of the simplicifcation, etc., done below is not correct. For example, off 50-50 the reduction to Ising model is in a field (chemical potential difference for alloy) and the transition temperature is not that given below.

Some Definitions:

The Model Hamiltonian: For a simple model of alloy, one can use pairwise interactions. For a binary alloy, this requires three types of interactions (which can still be long-ranged); namely, V^AA(r), V^BB(r), and V^AB(r). For a rigid lattice, the Hamiltonian is

H = -(1/4) _{i < j} (2V^AB_ij - V^AA_ij - V^BB_ij ) s_i s_j.
The Ising Hamiltonian (a simplification): We will consider only nearest-neighbor interactions. In addition, as explained below, we will choose V^AA = V^BB = 0. Hence, the Hamiltonian is that of a nearest-neighbor Ising model, or H = - (1/2) _{i < j} V^AB_ij s_i s_j, which favors ordering (i.e. A-B neighbors) if V^AB < 0, or favors clustering or phase segregation (i.e. A-A or B-B neighbors) if V^AB > 0. We want ordering (i.e., CuAu) so V^AB will be negative in what follows (see below). We may choose A atoms to be "spin up (+1 magnetization)" and B atoms "spin down (-1 magnetization)". Therefore, clustering interactions lead to Ferromagnetism (FM), whereas ordering interactions lead to Anti-Ferromagnetism (AFM).

The Vacancy: If a site is not occupied (i.e., it has a vacancy), then there is no atom or zero "spin" residing and the site has zero magnetization. In the KMC code, we use "+1", "-1", and "0" to identify the types of species occupying a lattice site., which are output in the *.res files. It is important to note that vancancy-assisted migration is important in alloys, especially to help order (or segregate) the atoms. Otherwise, there are large barriers to interstitial migration of atoms. Note that a vacancy helps assisted the alloy reach thermal equilibrium by finding the most favorable site for an atom to occupy.

The Short- and Long-Range Order: The degree of short-range order (SRO) in an alloy can be quantified by the so-called Warren-Cowley SRO parameters, _lmn(T), for a fixed (average) composition, C^B, which are the pair-correlations. For A-B pairs, _lmn^AB(T) = 1 - P_lmn^AB(T)/C^B , where P_lmn^AB(T) is the conditional probability that an atom A is at the origin and there is a B atom at the neighbor shell (site) lmn . Notice that if the sites are homogeneous random, then P(T) = C^B and the correlations are ZERO! However, even at high temperatures where the sites are random on average, there are local correlations, which are the SRO. The Long-range order (LRO) is that found below a phase transition The LRO = C_i^B - C^B, where the C_i^B is the local site probability (composition) of finding a B atom (found in the MC average below transition). Notice that this, in the Ising Model, is the magnetization, so the LRO is non-zero only below critical temperature. (See definitions in code.)

The Radiation Damage: If random radiative decay particles impinge on an alloy, then recoil damage creates disorder and lattice defects (like Schottky, Frenkel, interstitial defects). In our simple 2D Alloy model, we can do this by random swaps of atoms that are nearest neighbors. (See Question 4 below).

Getting Started

• Download the 2-D Kinetic Monte Carlo code (w/o or w/ ballistic jumps)
  • The source code without ballisitic jumps- ordering.f,      and Save As    ordering.f .
  • The source code with ballisitic jumps- ordering_ball.f,      and Save As    ordering_ball.f .
    Top lines of code give compiler options for most workstations.
    If you use ordering_ball.f exclusively, remember to make ballistic jump rate ZERO for Questions 1-3!
  • ordera.res, orderb.res, orderv.res, and order.res - The outputs are described in source.
    The code will not overwrite the output files, so, if you want to keep them, you must rename them.
• Download the input data
  • Example input DATA file - order.data,      and Save As    order.data .

NOTE: The total number of vacancy jumps and sub-blocks over which averages are performed are indicated in order.data. You can explore what these do, but notice that increasing vacancy jumps by order of magntiude affects the elapsed time between blocks. If you do anything, it is better to increase jumps. You MUST do this when nearing a phase transition, e.g., 80-100 Million, whereas much above or below this temperature 10-20 Million may be fine. Experiment a little.

NOTE: Do not forget that to get reasonable results one must always average over multiple MC runs. You cannot just run one temperature for one 80 Million Steps and stop.

Question 1 - Flowchart

• Describe how the time is incremented. Why is there no ln R as required for Poisson Process?
  (Hint: Thermal vs. Kinetic and Computational Savings.)
• How are the rates for the vacancy jumps are determined. Why?
• Describe what SRO and LRO is telling you.

Question 2 - Initialization with Random Configuration

NOTE: Run for c_B=0.50 A-B alloy, which is equivalent to square lattice Ising Model, or an 2-D FCC CuAu. The pairwise interactions interaction have been set to V_AA=V_BB=0 and _AB= -0.05 eV, so as to produce an ordering (AFM) type phase diagram. We know that 2D Ising Model has T_crit= (2/ln(1 + sqrt(2))/Z (with the number of nearest-neighbors Z=4 for a square lattice). Therefore, T_crit= 0.5673 eV = 658 K. Hence you know your answer here before you start!

NOTE: The nearer the T_crit the more the jump attempts must be increased, because fluctuations are large near first-order transitions.

NOTE: Do not forget that to get reasonable results in MC one must always average over multiple MC runs (to average over fluctuation). You cannot just run one temperature for one 80 Million Steps and stop. Required Vacancy hops can be very large near T_crit (80-100 MIllion hops) and much lower (1-10 Million hops) at high temperature. (Similar can be said in Question 3.)

NOTE: You may want to speed up lab by sharing MC runs amongst groups to do averaging at given temperature, especially near critical temperature.

Question 3 - Initialization with Ordered Configuration

NOTE: Recall that the nearer to T_crit the more the jump attempts must be increased.

Question 4 - Special Effects: Radiation Damage (Kinetic vs. Thermal Effects)

Typically, it is considered that radiation damage will disorder an alloy, such as the 2-dimensional one that you are studying. Why?

From what you have found from Questions 2 and 3, consider an additional rate in the KMC arising from radiation damage, a purely stochastic event.

• Explain physically what to expect to happen (due to both kinetic and thermal factors) to the low-T portion of the phase boundary, rather than the usual phase LOOP meeting at c=0.50 exactly at T=0 K.
• How would you incorporate this stochastic effect (the random radiation damage) into this code?
  HINT: Consider the jump portion of code where time is calculated and explain how to put in stochastic jumps. If you decide to try to first do this yourself, modify the MC portion of the code. The neighbor list update is really the tricky part, if you assume the ballistic jumps arise from a constant rate independent of environment.
• We have a 2-D Kinetic Monte Carlo code with Ballistic Jumps included (ordering_ball.f ). Importantly, due to scales set in code, ballistic jump rates are allowed to vary from 0 to 10 in the input data file, order.data )