Using the Sentaurus Materials Workbench for studying point defects

Version: P-2019.03

In this tutorial, we will use the comprehensive, and highly-automated, framework provided by the QuantumATK module Sentaurus Materials Workbench (SMW) to calculate the formation energies and transition levels for a variety of charged defects in GaN. SMW is designed to make complex workflows of ab initio calculations more accessible, and is focused on specific applications of interest. In this tutorial, we will use the defect list functionality to easily study several point defects in GaN, with 3 different structure types and 2 charge states, for a total of 14 distinct defects.

introbar

Background

The central feature of this tutorial is called Defect Lists and is explained in more detail here: Operating with defect lists. It can be used to generate a list of specific defects of a particular type, in a number of charge states, and calculate the formation energies of these defects and the transition (or trap) levels between each of these charge states. After a type of defect is chosen, e.g. vacancies, QuantumATK requires no input from the user to determine which defects exist in a particular host material, as it will automatically determine all of them from symmetry analysis. However, SubstitutionalList and InterstitialList do require the user to specify the element that will be added to the host material as a substitutional or interstitial defect. QuantumATK will use the ChargedPointDefect study object for the individual defect calculations, which is explained in detail in this tutorial: Formation energies and transition levels of charged defects.

Getting started

Open the QuantumATK GUI, NanoLab, and select Create New… in the Project Dialog. Create a new project in an empty folder with an appropriate name, and make sure to select Sentaurus Materials Workbench as the Type.

../../_images/project-dialog.png

In QuantumATK, projects are a way to group related calculations together, and is closely related to the folder/directory structure of the OS. The chosen type of project helps QuantumATK determine what kind of work you want to do, and modify the GUI accordingly. Specifically, choosing the SMW project type will hide many of the, for this purpose, less-relevant plugins, allowing you to focus on these advanced features.

Material specifications

Before we can do the defect calculations, we first need what is called a Material Specification. This is the definition of the material we wish to study and the computational model that will be used. QuantumATK includes a database with example settings for some relevant materials, but if you wish to study a material that is not covered, you can provide your own specification as part of the input. In this case, we will use a material specification that has been generated particularly for use in this tutorial, and tested to ensure good results for the specific purpose we will look at here.

Tip

Creating a good MaterialSpecifications requires ab initio expertise to ensure that the computational model is sufficiently accurate. However, the philosophy behind the design of SMW is that it should not require ab initio expertise to use it, simplifying collaboration between ab-initio and higher-level modeling.

Setting up the calculation

To set up the calculation, open the Defect/Diffusion Scripter in the panel on the right.

../../_images/defect-scripter-plugin.png

Click the Database button to open the QuantumATK database of material specifications for automated defect calculations. It contains material specifications for a total of 29 materials, with parameters chosen to provide general applicability and high precision for studies of defects in this material. Choose the GaN_300_DFT_SCAN entry by double-clicking or selecting it and then clicking the OK button.

../../_images/database-widget.png

Now add three successive defects by clicking Add and then selecting first Substitutional, then Interstitial, and finally Vacancy. We wish to add Magnesium atoms and investigate the 0 and -1 charge states, so select the substitutional defect, click on Phosphorous to open a window with the periodic table, and select Mg there. After this window has closed, click on the downwards-pointing arrow in the left Charge states box to change the value to -1. The window should now look like this:

../../_images/defect-scripter-widget.png

Now click on Interstitial and do the same, then Vacancy and change the charge state here as well. With this, the script is ready to run, but let us first take a look at it. Click the sendto_icon and choose the Editor. You should now see the following script:

from SMW import *

# Material
material = MaterialSpecificationsDatabase.GaN_300_DFT_SCAN

# Defects
# Defect list 0
charge_states = [-1, 0]
defect_list_0 = SubstitutionalList(
    material_specifications=material,
    substitutional_element=Magnesium,
    charge_states=charge_states
)
nlsave('GaN_300_DFT_SCAN.hdf5', defect_list_0)
defect_list_0.update()

# Defect list 1
charge_states = [-1, 0]
defect_list_1 = InterstitialList(
    material_specifications=material,
    interstitial_element=Magnesium,
    charge_states=charge_states
)
nlsave('GaN_300_DFT_SCAN.hdf5', defect_list_1)
defect_list_1.update()

# Defect list 2
charge_states = [-1, 0]
defect_list_2 = VacancyList(
    material_specifications=material,
    charge_states=charge_states
)
nlsave('GaN_300_DFT_SCAN.hdf5', defect_list_2)
defect_list_2.update()

The script is rather simple, as the SMW module handles all details of the workflow without needing further input from the user. First, it extracts a material specification from the QuantumATK database, and then defines the three defect lists that we set up previously. Note that they are run sequentially, so it first calculates the substitutions, then the interstitials and then the vacancies.

Replacing the material specification

The material specifications in the database are designed to be generally applicable and give high-precision results for a wide range of defects. However, this also makes them a bit time-consuming to use, and we will therefore use a different material specification for GaN in this tutorial. It is specifically designed and tested for the magnesium defects we are investigating here, and should be used with care for other defects. In the original script, remove the line material = MaterialSpecificationsDatabase.GaN_300_DFT_SCAN with the material specification from the database, and replace it with the contents of this file: mat-spec-zero-entropies.py. The resulting script should be similar to this: SMW_zero-entropies.py.

Contents of the material specification

We will now go through the main points of this material specification. If you are not interested in these details, you may skip ahead to the next section: Run the calculation. We will look only at the material specification itself, which looks like this:

material = MaterialSpecifications(
    pristine_configuration=gan_primitive,
    supercell_repetitions=(3, 3, 3),
    formation_energy_calculator=calculator,
    atomic_chemical_potentials=[ga_mu, mg_mu],
    dielectric_constant=gan_dielectric_constant,
    bulk_modulus=173 * GPa,
    assumed_formation_entropy=0.0 * Joule/Kelvin,
)

Some of the parameters are defined further up in the script, and you may consult the full material specification file for those details.

Note

For more information, see the manual page: MaterialSpecifications

The supercell parameters

material = MaterialSpecifications(
    pristine_configuration=gan_primitive,
    supercell_repetitions=(3, 3, 3),
    formation_energy_calculator=calculator,

The first few lines specify the host material and the calculator, or computational model, to be used for it. Specifically, they define the unit cell of the pristine host material, the number of times it is repeated to create a supercell containing a defect, and the main calculator to use when calculating the formation energy of the defects.

Atomic chemical potentials

    atomic_chemical_potentials=[ga_mu, mg_mu],

This line specifies the atomic chemical potentials that will be used when calculating the formation energies for the defects. This can be given either as an energy and/or entropy per atom, a specific atomic configuration for the element, or left out entirely. In all cases, QuantumATK will use the calculator specified above to calculate any missing properties.

Host material properties

    dielectric_constant=gan_dielectric_constant,
    bulk_modulus=173 * GPa,

These lines specify some properties of the host material which are used to correct finite-size effects in the results. For that reason, it is more important that they are consistent with the computational model specified above than with experimental results.

Formation entropy

    assumed_formation_entropy=0.0 * Joule/Kelvin,
)

Finally, it contains the definition of of the assumed formation entropy. In this case, we do the calculation without entropic contributions and therefore define the formation entropy as zero. If we had not defined it, it would be calculated, which is a rather time-consuming calculation. If we knew the value from somewhere else we could have defined it as that value instead of 0. We will revisit these contributions later.

Run the calculation

The script is now ready to be run. Submit it to your favorite cluster - in our case, it ran in about 6.5 hours on a 24-core node with 256 GB of memory. It requires about 1.7 GB of disk space to store all the results.

Tip

You can use the QuantumATK Job Manager for submitting jobs to clusters. See this tutorial for more details: Job Manager for remote execution of QuantumATK scripts.

The SMW framework will now run a series of individual defect calculations, and set up a folder structure to keep track of them. If it does not finish, either due to external circumstances, such as a power failure, or simply running out of walltime, you can easily restart it by simply running the same script again. QuantumATK will check the status of the calculations and start from the interrupted task.

In the end, you should see a log-file containing blocks like these:

+------------------------------------------------------------------------------+
| Defect List                                                                  |
+------------------------------------------------------------------------------+
| 2 defect(s) will be calculated.                                              |
+------------------------------------------------------------------------------+
+------------------------------------------------------------------------------+
| Calculating defect 1 / 2:                                                    |
|   Results and log will be saved to the following folder:                     |
|   substitutional/substitutional_Mg^Ga_0/000/3848cac9fcfea4adc57af73c58700111/|
+------------------------------------------------------------------------------+
+------------------------------------------------------------------------------+
| Calculation 1 / 2 finished.                                                  |
|   Convergence status: Converged. All results will be available.              |
+------------------------------------------------------------------------------+

As could be seen in the script, each type of defect is defined as a separate list. We can see that the first defect list contains substitutions, and that the first of two is a Mg substitution at a Ga site. The location of this particular sub-calculation is also specified, so we know where to go if we want to take a closer look at it.

Analyzing the results

After the calculations have finished, we can now look at the results. Make sure to select the SMW_zero-entropies.hdf5 file in the panel on the left, so it appears on the LabFloor in the center. Click that line, and you should see something like this:

../../_images/labfloor-screenshot.png

Note

If the project type is QuantumATK instead of Sentaurus Materials Workbench, you will see more plugins in the panel on the right than what is shown above.

Select all three defect list objects and click on the Defect List Analyzer. The following window will now appear.

../../_images/defect-list-analyzer-screenshot.png

This window is composed of three main parts. On the left, there is a table showing the different defects and their charge states, and the calculated formation energy for each. Below this table, you can choose the Fermi level to be used when calculating the formation energy; the default is the midpoint of the band gap. You also find filters that can be applied to the list of defects, in cases where the defect lists are longer than those we are studying here. On the right, there is a plot showing the trap levels for each defect type, with the valence and conduction bands indicated as blue areas and the Fermi level as the dotted line in the center.

If we look at the formation energies, we quickly see that two values are close to 0 eV, while all the others are much higher and no smaller than 3.5 eV. We will therefore focus on the first two, which are the neutral and negative charge states for a Mg substitution at a Ga site. In the plot on the right, we see that the transition between them, which signifies a trap level, occurs slightly above the top of the valence band. If you hover over it, you will see that it is 249 meV above the top of the valence band. This is in excellent agreement with literature, where it has been found to be 260 meV above the top of the valence band [1]. You may also click SubstitutionalList in the table, in which case the plot will show only the substitutional defects, or select a subset of specific defects (hold down ctrl to add entries) to see only those in the plot. Click in the empty area at the bottom of the table area to reset the plot, or press Esc.

We see that the transition between the same two charge states occur only a little higher for the Ga vacancy (hover over it to identify the vacancy site). However, we can see from the table that the formation energy for the Ga vacancy is much higher, so those defects will be much rarer under the modeled conditions.

The other defects also have transitions between the two charge states, but they are not shown, as they are outside the band gap. To find them, click the Edit button Plot Editor Tool. Edit the plot layout, axes, titles, lines etc. to change the properties. Find the y-axis and change the limits to -0.1 and 3.5, for example, and the x-axis limits to -0.5 and 7.5 to produce this plot:

../../_images/defect-list-analyzer-plot-zoom.png

We now see that the trap levels for the other defects are (high) in the conduction band.

Tip

You can also use the mouse scroll to quickly zoom out and get an overview.

Expanding the data set

After taking a first look at your results, you might decide that you would like to expand the data set to include some more defects or charge states. This is quite easy to do, as there is intelligent restart functionality in the SMW defect list and its constituent features. Open the script (SMW_zero-entropies.py) and edit lines 79, 89 and 99 to also include the +1 charge state, so they look like this:

charge_states = [-1, 0, +1]

Now re-run the script in the same folder. QuantumATK will automatically detect that some of these calculations have already been done, and only do the necessary new calculations. In this case, the new calculations took about the same time as the original set, but it depends on the system.

Now look at the labfloor again - you will see that three more objects have appeared in the file. These are the new defect lists, containing both the original results, and the new results for the +1 charge states. Select all three, and open them in the Defect List Analyzer, as before. It should now look like this.

../../_images/defect-list-analyzer-screenshot-add-p1.png

Note that the original trap levels are still there, while 3 new ones appear, close to the band gap. In this way, we could easily expand the data set in the study, based on preliminary analysis of the original data.

Discussion and summary

In this section, we have shown how to use the Sentaurus Materials Workbench (SMW) feature of Defect Lists to easily study several different point defects in GaN, and even consider them in different charge states. Specifically, we studied Mg substitutions and interstitials, as well as vacancies in neutral and negatively charged states. We found that the Mg substitution at a Ga site has a low formation energy and a trap level 249 meV above the valence band maximum, consistent with literature.

Note that all of these calculations have been done without taking vibrational entropy and zero-point energy into account. These terms require the computation of a DynamicalMatrix and a PhononDensityOfStates for each sub-system, the first of which is very demanding. The computational cost normally increases by more than an order of magnitude when including these terms, and they are often small. If possible, it is therefore a good idea to evaluate whether their inclusion is relevant or not for the system at hand. In the next section, we will include these terms and see how it changes the results.

Including vibrational contributions

It is quite straightforward to change the original script to also include, and calculate, the vibrational terms. Simply delete the lines:

    assumed_formation_entropy=0.0 * Joule/Kelvin,
    assumed_transition_prefactor=10**13*Second**(-1),

There are now two ways to run the script. One is to parallelize over the defects in each defect list, and the other is to not do so. In this case, we choose to parallelize over defects, which is slightly more complicated, but provides a faster time-to-result at full computational efficiency. To achieve this, we need to ensure that the choice of compute nodes is compatible with the number of defects for each defect list. From before, we know that there are 2 substitutional defects, 2 vacancies and 3 interstitials. It will therefore be most efficient if we first run the substitutional defect list on 2 nodes, then the vacancy list on 2 nodes and finally the interstitial list on 3 nodes. This means that we will be running these three scripts in succession. In principle, they can also be run concurrently, but this would give a small risk of corrupting the central .hdf5 file if more than one job wrote to it at the same time. We ran these three scripts in succession in a new and empty directory, using the hardware and computational time specified in the table. This amounts to about a factor of 35 more CPU-time required for this calculation, compared to the previous one without vibrational terms.

Note

If you have already run the calculations without vibrations and wish to add these terms, you may save a little time by running it in the same directory. QuantumATK will detect that some calculations have already been done, and not re-calculate them.

Table 2 Phonon calculations

Script

CPU’s

Walltime (hours)

Total CPU-time (hours)

smw-gan-mg-sub.py

2x24

42.0

2014

smw-gan-mg-vac.py

2x24

42.0

2017

smw-gan-mg-in.py

3x24

22.1

1593

After the calculations have run, the main results will be saved in the file SMW_with-entropies.hdf5. In addition to the sub-folders and files created before, the phonon calculations will also produce several thousand log-files, representing the many different displacements that are calculated along the way. They can be studied if needed for troubleshooting or verification.

Tip

For more details on how vibrational terms are included, see the manual page of HarmonicChargedPointDefect.

Analyzing the results

Go to the LabFloor and make sure the file SMW_with-entropies.hdf5 is selected. As before, select the three items in it, and click on the Defect List Analyzer on the right. You should now see the following window:

../../_images/defect-list-analyzer-phonons-screenshot.png

This looks very familiar, as the inclusion of vibrational quantities changes very little for the defects investigated here. It turns out that the 0/-1 trap level moves 2 meV, for Mg substituted at a Ga site, while the underlying formation energies change about 5 and 2 meV respectively. Most other formation energies and trap levels also change very little, except for the formation energy for neutral Mg substitution at a N site. This particular formation energy changes by more than 1 eV with the inclusion of vibrational terms.

When seeing surprising results, it is good practice to investigate further, and in this case we can search through the log-files for warnings. It turns out that the log-file for elemental nitrogen (cpd_atomic_chemical_potential_Nitrogen.log) contains this message:

################################################################################
#                                                                              #
# WARNING                                                                      #
#                                                                              #
# ============================================================================ #
#                                                                              #
# The sum rules are not obeyed in this PhononDensityOfStates calculation.      #
# There are  147 non-positive eigenmodes. Be careful if you want to calculate  #
# thermodynamic properties.                                                    #
#                                                                              #
################################################################################

This means that the vibrational spectrum of nitrogen is, most likely, poorly described using the computational parameters in our current material specification. It follows that formation energies which contain vibrational terms for nitrogen, such as the Mg substitution at a nitrogen site, are a bit uncertain and should be interpreted with care.

Tip

For more information on phonon calculations, see this tutorial: Calculating and using Dynamical Matrix.

Conclusions

We found that the Mg substitution at a Ga site has a low formation energy and a trap level 249 meV above the valence band maximum, when not considering the vibrational terms. Including the vibrational terms takes approximately 35 times longer and changes the trap level position to 246 meV above the valence band maximum. This means that inclusion of the vibrational terms was not important in this case, and certainly not important enough to justify the large computational cost.

References