Modeling Vacancy Diffusion in Si0.5 Ge0.5 with AKMC

Version: 2015.1

In this tutorial, we will model vacancy diffusion in a Si0.5 Ge0.5 alloy using the Adaptive Kinetic Monte Carlo method (AKMC). AKMC is an algorithm for performing Kinetic Monte Carlo (KMC) simulations on-the-fly without having to have a predefined set of states and reaction mechanisms. The reaction mechanisms are determined by locating saddle points on the potential energy surface and then using Harmonic Transition State Theory (HTST) to calculate the reaction rate. You can learn more about the AKMC method in this tutorial: Adaptive Kinetic Monte Carlo Simulation of Pt on Pt(100).


For typical MD simulations, it is impractical to model more than a couple of hundred nanoseconds. However, by using AKMC, we will be able to model the dynamics of the system on a longer time scale. This is because the KMC simulation describes the dynamics of the system as transitions between states that – once found – can be computed very quickly.


Obtaining an Initial Structure

In order to generate the initial structure we will use the script (shown in the panel below). This script generates a silicon supercell and then randomly replaces Si atoms with Ge atoms in order to achieve a 50:50 ratio. The focus of this tutorial is on the AKMC algorithm, and we will therefore not go into detail with this part. You can download the script below the box containing the code. Run the script either by dragging it from the main VNL window and dropping it on the Job Manager icon job_manager_icon, or directly through the terminal.

from NanoLanguage import *

def random_alloy(random_seed, na, nb, nc, num_vacancies, si_ge_ratio=0.5):
    Generate a random SiGe alloy.

    @param random_seed   : The random seed used to generate this crystal.
    @type                : int

    @param na            : The number of times to repeat the primitive cell
                           along the a direction.
    @type                : int

    @param nb            : The number of times to repeat the primitive cell
                           along the b direction.
    @type                : int

    @param nc            : The number of times to repeat the primitive cell
                           along the c direction.
    @type                : int

    @param num_vacancies : The number of vacancies in the crystal.
    @type                : int

    @param si_ge_ratio   : The ratio of Si to Ge atoms.
    @type                : float

    @return BulkConfiguration

    rng = numpy.random.RandomState(random_seed)

    # Define a Silicon lattice.
    lattice = FaceCenteredCubic(5.4306*Angstrom)
    elements = [Silicon, Silicon]
    fractional_coordinates = [[ 0.  ,  0.  ,  0.  ],
                              [ 0.25,  0.25,  0.25]]

    bulk_configuration = BulkConfiguration(

    # Repeat the structure the requested number of times.
    bulk_configuration = bulk_configuration.repeat(na, nb, nc)

    n = len(bulk_configuration)
    keep_indices = rng.choice(range(n), size=n-num_vacancies, replace=False)

    elements = numpy.array(bulk_configuration.elements())[keep_indices]
    coordinates = bulk_configuration.cartesianCoordinates()[keep_indices]
    lattice = bulk_configuration.bravaisLattice()

    # Make the alloy.
    n = len(elements)
    num_germanium = int(n * (1.0 - si_ge_ratio))
    germanium_indices = rng.choice(range(n), size=num_germanium, replace=False)
    for i in germanium_indices:
        elements[i] = Germanium

    # Set up configuration
    bulk_configuration = BulkConfiguration(

    return bulk_configuration

# Setup the calculator.
potentialSet = StillingerWeber_SiGe_1995()
calculator = TremoloXCalculator(parameters=potentialSet)

# Generate a random alloy.
bulk_configuration = random_alloy(

# Set the calculator.

# Optimize the geometry.
bulk_configuration = OptimizeGeometry(

# Save the result to
nlsave('', bulk_configuration, 'minimized')

Here you have generated a 4x4x4 supercell of SiGe, calculated the optimized geometry and saved the structure to Note that this is a very small system, useful only for illustrative purposes. You will need a larger system for scientifically valid modelling. You can visualize the structure in VNL by selecting the item labeled minimized from the LabFloor and drag and drop it on the Viewer icon viewer_icon. The result is shown below.


There is a vacancy in the structure, but it is difficult to see. You will perform a LocalStructure calculation to see which atoms do not have the typical diamond-like packing.

  1. Drag the minimized configuration from the LabFloor to the Scripter icon script_generator_icon.
  2. In the Scripter, double click on Analysis ‣ LocalStructure.
  3. Choose the default output file to be
  4. Send the script to the Job Manager using the sendto_icon button and run the calculation. It should take no more than a minute.

Now you can get a better look at the structure around the vacancies. From the main window, open the LocalStructure item in the Viewer and expand the Local Structure panel in the sidebar. Select Diamond-like. Now all of the atoms that are not near the vacancies, and thus have the ideal crystal structure, have been selected. Open the Properties dialog and in the Atoms tab and change the opacity to 0.05 for the selected atoms. This will highlight the atoms near the vacancies, by making the others almost transparent. The result is shown in the figure below. Note the obvious vacancy in middle of the cluster on the left. The cluster on the right is also close to the vacancy, but shifted due to the periodic boundary conditions.


Running the AKMC Simulation

Now that you have an initial structure, you can run the AKMC simulation. Drag from the main VNL window to the Editor icon editor_icon. Then append the following code to the script and save it as

# --------------------------------------------------------------
# --------------------------------------------------------------

# Reusing existing MarkovChain object if it already exists, otherwise a new one is created
if os.path.isfile(''):
    markov_chain = nlread('')[0]
    markov_chain = MarkovChain(bulk_configuration, TotalEnergy(bulk_configuration).evaluate())

# Reusing existing KMC object if it already exists, otherwise a new one is created
if os.path.isfile(''):
    kmc = nlread('')[0]
    kmc = None

# Modify the default maximum number of NEB images to 5.
saddle_search_parameters = SaddleSearchParameters(max_neb_images=5)

# Setup the AKMC simulation.
akmc = AdaptiveKineticMonteCarlo(markov_chain,

# Run 100 saddle searches.

The first lines enable a restart of the simulation, by reusing the existing files if they exist. For a detailed description of the parameters for the AKMC object, please see the tutorial: Adaptive Kinetic Monte Carlo Simulation of Pt on Pt(100).

Now run the script. Running the script should take no more than a half an hour on a desktop machine. At the end you will have a MarkovChain object stored inside the file as well as a KineticMonteCarlo object stored inside the file. The MarkovChain object stores all of the states (configurations) that were found during the simulation as well as the connectivity between them. It can be used to calculate the full rate matrix for the whole system. The KineticMonteCarlo object stores the time and state-to-state trajectory and is also able to take steps along the Markov chain (i.e. KMC steps).

In the LabFloor, click on the KineticMonteCarlo object and select Text Representation from the panel on the right. The output should look like, but probably not be identical to, the one seen below:

# Item:  0
# File:  /home/$USER/your-projects/AKMC-empirical-SiGe/
# Title: - gID000
# Type:  KineticMonteCarlo
    Step #   Time (seconds)  State
         0   0.00000000e+00      0
         1   5.91307921e-08      3
         2   5.96128092e-08      7
         3   6.12494862e-08      3
         4   6.14761792e-08      7
         5   6.20331334e-08      3
         6   6.39890780e-08      7
         7   6.46268953e-08      3
         8   6.59538019e-08      7
         9   6.66472197e-08      3
        10   6.71560053e-08      7
        11   6.77515572e-08      3
        12   6.79111154e-08      0
        13   8.69634889e-07      3
        14   8.71340919e-07      7
        15   8.71394296e-07      3
       300   1.11970777e-06      7
       301   1.12135663e-06      3
       302   1.12355061e-06      7
       303   1.12434612e-06      3
       304   1.12444238e-06      7
       305   1.12448188e-06      3
       306   1.12470034e-06      7
       307   1.12537843e-06      3
       308   1.12544785e-06      6
       309   1.14161108e-06     13

Only the top and bottom lines have been shown here. The first column is the step number for the KMC algorithm, the second column represents the total elapsed time in the KMC simulation and the third column shows the current state number, which corresponds to a particular minimum energy configuration in the MarkovChain object.


Note how the final time is greater than a microsecond. This is a much longer timescale than we could probe using MD alone.

In this case there are more KMC steps than the 100 saddle searches we specified in the script. This shows that the KMC algorithm, to some degree, has been jumping back and forth between the same states. However, as this is a stochastic method it is also entirely possible to see the opposite, with the total number of KMC steps being smaller than the specified number of saddle searches. This behavior is also very system-dependent.

On the LabFloor you will also find an AKMCLog object. Open it by clicking on Text Representation. You should see a message equivalent to what is shown below. This shows the results of the 100 saddle searches. The first column is the id of the initial state used for the saddle search, the second column is the search number (out of the 100 saddle searches), the third column is the calculated confidence that the processes starting in the current basin are adequately sampled and the last column is the result of the search. In this case, the algorithm finds a total of 4 new states starting from state 0, 13 searches find a duplicate state and 14 searches find no new states. In the last line the initial state id has changed, showing that the KMC algorithm has taken a step to a new state, and the saddle searches then continue from there.

# Item:  0
# File:  /home/$USER/your-projects/AKMC-empirical-SiGe/
# Title: - gID000
# Type:  AKMCLog
state id search number confidence  message
       0             0   0.000000  No reaction occurred in 5 ps
       0             1   0.000000  Found new state
       0             2   0.000000  No reaction occurred in 5 ps
       0             3   0.000000  No reaction occurred in 5 ps
       0             4   0.000000  No reaction occurred in 5 ps
       0             5   0.000000  Found new state
       0             6   0.000000  Found new state
       0             7   0.333404  Found saddle connecting to state 3 again
       0             8   0.333404  No reaction occurred in 5 ps
       0             9   0.776412  Found saddle connecting to state 1 again
       0            10   0.809398  Found saddle connecting to state 3 again
       0            11   0.821533  Found saddle connecting to state 3 again
       0            12   0.821533  No reaction occurred in 5 ps
       0            13   0.821533  No reaction occurred in 5 ps
       0            14   0.821533  No reaction occurred in 5 ps
       0            15   0.821533  No reaction occurred in 5 ps
       0            16   0.819629  Found new state
       0            17   0.907678  Found saddle connecting to state 2 again
       0            18   0.907678  No reaction occurred in 5 ps
       0            19   0.912132  Found saddle connecting to state 3 again
       0            20   0.955860  Found saddle connecting to state 1 again
       0            21   0.964572  Found saddle connecting to state 2 again
       0            22   0.964572  No reaction occurred in 5 ps
       0            23   0.964572  No reaction occurred in 5 ps
       0            24   0.964572  No reaction occurred in 5 ps
       0            25   0.964572  No reaction occurred in 5 ps
       0            26   0.980659  Found saddle connecting to state 1 again
       0            27   0.983863  Found saddle connecting to state 2 again
       0            28   0.985042  Found saddle connecting to state 2 again
       0            29   0.985476  Found saddle connecting to state 2 again
       0            30   0.991394  Found saddle connecting to state 1 again
       3            31   0.000000  Found new state

Now, we can use the KineticMonteCarlo and MarkovChain objects to make a visualization of the simulation. Use the script (see below), which loads the two objects and uses the state-to-state trajectory to construct a Trajectory object, where all of the diamond-like atoms have been deleted. Drag the script onto the Job Manager to run it. The run time is highly dependent on the length of your KineticMonteCarlo object, but should take no more than a few seconds per entry, equivalent to no more than half an hour.

from NanoLanguage import *

markov_chain = nlread('')[0]
kmc = nlread('')[0]

trajectory = Trajectory()
kmc_trajectory = kmc.trajectory()

for i in xrange(len(kmc_trajectory)-1):
    print "%i/%i" % (i,len(kmc_trajectory))

    state_number = kmc_trajectory[i]
    next_state_number = kmc_trajectory[i+1]

    configurations = []

    configurations.append(markov_chain.getSaddleConfiguration(state_number, next_state_number))

    local_structures = [ LocalStructure(c).evaluate() for c in configurations ]
    delete_atoms = set()
    for local_structure in local_structures:
        diamond = numpy.where(local_structure == 5)[0]
        delete_atoms = delete_atoms.union(set(diamond))

    for configuration in configurations:
        configuration = configuration._fastClone()

nlsave('', trajectory)

The Trajectory output from this script is saved to Opening it in the Viewer shows the state-to-state dynamics from our AKMC simulation.