Modeling Vacancy Diffusion in Si0.5 Ge0.5 with AKMC¶
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
initial-structure.py (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
, 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( bravais_lattice=lattice, elements=elements, fractional_coordinates=fractional_coordinates ) # 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( bravais_lattice=lattice, elements=elements, cartesian_coordinates=coordinates ) return bulk_configuration # Setup the calculator. potentialSet = StillingerWeber_SiGe_1995() calculator = TremoloXCalculator(parameters=potentialSet) # Generate a random alloy. bulk_configuration = random_alloy( random_seed=1, na=4, nb=4, nc=4, num_vacancies=1, si_ge_ratio=0.5) # Set the calculator. bulk_configuration.setCalculator(calculator) # Optimize the geometry. bulk_configuration = OptimizeGeometry( bulk_configuration, max_forces=0.01*eV/Ang, max_stress=0.01*eV/Ang**3, max_steps=500, max_step_length=0.2*Ang, trajectory_filename='alloy.nc', optimizer_method=LBFGS(), ) # Save the result to alloy.nc. nlsave('alloy.nc', bulk_configuration, 'minimized')
Here you have generated a 4x4x4 supercell of SiGe, calculated the optimized
geometry and saved the structure to
alloy.nc. 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 . 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.
- Drag the minimized configuration from the LabFloor to the Scripter icon .
- In the Scripter, double click on .
- Choose the default output file to be
- Send the script to the Job Manager using the 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
initial-structure.py from the main VNL window to the
Editor icon . Then append the following code to the script and save it as
# -------------------------------------------------------------- # AKMC # -------------------------------------------------------------- # Reusing existing MarkovChain object if it already exists, otherwise a new one is created if os.path.isfile('akmc_markov_chain.nc'): markov_chain = nlread('akmc_markov_chain.nc') else: 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('akmc_kmc.nc'): kmc = nlread('akmc_kmc.nc') else: 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, kmc=kmc, calculator=calculator, kmc_temperature=300*Kelvin, md_temperature=1200*Kelvin, saddle_search_parameters=saddle_search_parameters, constraints=, confidence=0.99, write_searches=False, write_kmc=True, write_markov_chain=True, write_log=True) # Run 100 saddle searches. akmc.run(100)
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
akmc_markov_chain.nc file as well as a KineticMonteCarlo
object stored inside the
akmc_kmc.nc 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/akmc_kmc.nc # Title: akmc_kmc.nc - 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/akmc_log.nc # Title: akmc_log.nc - 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('akmc_markov_chain.nc') kmc = nlread('akmc_kmc.nc') 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.getStateConfiguration(state_number)) 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) delete_atoms = delete_atoms.union(set(diamond)) for configuration in configurations: configuration = configuration._fastClone() configuration._deleteAtoms(list(delete_atoms)) trajectory._addConfiguration(configuration) nlsave('trajectory.nc', trajectory)
The Trajectory output from this script is saved to
Opening it in the Viewer shows the state-to-state dynamics from our AKMC