1 of 16

3D Materials Lab: Tutorials

A collection of tutorials contributed by the members of the 3D Materials Lab

Overview

The 3D Materials Lab is led by Dr. Prashun Gorai. We utilize computational tools, including first-principles methods and materials informatics, to discover and design functional materials for thermoelectrics, microelectronics, energy storage, ferroelectrics, photovoltaics etc. We are based at the Colorado School of Mines and the National Renewable Energy Lab located in beautiful Golden, Colorado.

These tutorials are made possible through the wonderful contributions of the team members.

Defects

Intro to Defect Calculations

Computational workflow for performing point defect calculations

Written by Prashun Gorai

This tutorial describes the step-by-step process for calculating point defect formation energy, defect concentrations, and associated carrier (electron, hole) concentrations. The workflows use open-source Python packages, and . This tutorial is intended to be a practical guide. Best practices and common pitfalls are discussed. This tutorial does not cover the post-processing and interpretation of the results, nor does it cover the theory of defect calculations.

Pylada-defects is an open-source python package for automation of first-principles point defect calculations. For details, please read the paper:

The formation energy of point defects is calculated using a supercell approach, and is given by

The defect formation energies are obtained from a suite of calculations (e.g., supercell relaxations, phase stability etc.) and therefore, one must carefully strategize these calculations for efficiency and to avoid common pitfalls. A summary of step-by-step procedure:

Perform geometry optimization (relaxation) of the unit cell.
Calculate the electronic structure with a dense k-point grid, using relaxed unit cell as input.
Calculate electronic structure with GW, if band gap correction is needed.
Perform phase stability calculations to determine range of elemental chemical potentials.
Create defect supercells and perform geometry optimization (typically, ionic relaxation).
Identify likely sites for interstitials using Voronoi tessellation, and follow specific sequential procedure described below to calculate formation energy of interstitials.
Calculate electronic and ionic contributions to dielectric constant.

Steps 1 and 2 should be performed sequentially. Steps 3-7 can be performed simultaneously. When handling hypothetical starting structures, it is recommended to perform Step 4 after Steps 1 and 2 to ensure phase stability.

This tutorial focuses on defect calculation with density functional theory (DFT) as the electronic structure method. The choice of DFT functional should be carefully considered for the problem at hand. In this tutorial, it is assumed that DFT calculations are performed with a semi-local exchange correlation functional such GGA-PBE with the band gap "corrected" with quasiparticle energies obtained from GW.

The defect energetics are very sensitive to the band gap. The well-known band gap underestimation in semi-local DFT can lead to highly erroneous defect energetics.

Step 1: Structure Relaxation

Relax the primitive or conventional cell using the DFT functional of choice. Perform full geometry optimization i.e., volume, cell shape, and ionic positions. The relaxed structure will be used in all subsequent steps.

Step 2: Electronic structure on dense k-point grids

The calculation of the electronic band structure is an important, but often overlooked, step in the defect calculation workflow.

Under construction...

Interpreting Defect and Energy Level Diagrams

A beginner's guide to understanding defect and energy level diagrams obtained from first-principles calculations

Written by Prashun Gorai

The formation energy of point defects is calculated using the supercell approach described here. The results of the defect calculations are typically presented in the form of "defect diagrams" and "energy level diagrams" – both are discussed below. These diagrams are useful for qualitative assessment of defect and doping properties (e.g., dopability) as well as quantitative determination of defect and electronic charge carrier (electrons, holes) concentrations. The discussion here will focus on semiconductors and insulators i.e., materials with a band gap.

Defect Diagrams

In a defect diagram, the defect formation energy $(\Delta E_{D,q})$ is plotted against the Fermi energy $(E_\mathrm{F})$ , which is typically referenced to the valence band maximum i.e., $E_\mathrm{F(VBM)} = 0$ . A schematic defect diagram is shown below, where $E_g$ is the band gap.

The defect diagrams are plotted at a chosen elemental chemical potential ( $\mu_i$ ) condition. Here, $\mu_i$ for every element i should be within the range where the material is thermodynamically stable against decomposition into competing phases. These ranges are determined from a grand potential phase diagram.

How are these defect diagrams plotted?

Each set of connected lines (red solid, red dotted, blue solid, blue dotted) represents a defect D, which could be, for example, a cation vacancy or substitution dopant. The slope of the line is the charge state q of the defect D. Recall that $\Delta E_{D,q}$ has a linear dependence on $E_\mathrm{F}$ . Therefore, the equation for calculating $\Delta E_{D,q}$ represents an equation of a straight line with slope q.

\Delta E_{D,q} = E_{D,q} - E_{H} + \sum_i{n_i \mu_i} + E_{corr} + ~ qE_\mathrm{F}

For a defect D, conventionally only the charge state q with the lowest formation energy at a given $E_\mathrm{F}$ is plotted as shown below. This way of plotting the data is motivated by two reasons: (1) to improve readability of the defect diagrams, and (2) the charge state q with lowest $\Delta E_{D,q}$ has the highest concentration at that $E_\mathrm{F}$ . This is shown schematically for a defect.

Even though only the charge state q with lowest $\Delta E_{D,q}$ is plotted at a given $E_\mathrm{F}$ , a non-zero Boltzmann population of other charge states exists, in reality.

What does the charge state/slope of a defect signify and how is related to doping?

Charged defects create electronic carriers – electrons and holes. A neutral defect (e.g. isoelectronic doping, alloying) does not create electronic carriers. A positively-charged defect has a positive slope in a defect diagram and represent a donor-like defect. In other words, the defect ionizes to a positively-charged state by giving up (donating, hence donor) electron(s). Donors tend to make the material n-type. Similarly, negatively-charged defects are acceptor-like and they tend to make the p-type. Whether a material is actually doped n- or p-type depends on several other factors. We will return to the topic of doping once we have understood a few additional concepts.

What is the significance of a charge transition level in defect diagram?

The charge transition level (q/q') of a defect is the Fermi energy at which the formation energy of charge states q and q' are equal. On a defect diagram, charge transition levels (CTLs) are identified by points at which the slope of the defect lines charge. A defect can have multiple CTLs or they can be absent within the band gap i.e. $0 \leq E_\mathrm{F} \leq E_g$ . In the schematic below, defect D1 has no CTLs, D2 has one CTL, and D3 has multiple CTLs.

The charge transition levels represent the approximate energetic location of defect states (a.k.a. defect levels) inside the band gap. These defect states give rise to many interesting phenomena that have implications for the electronic, optical, magnetic properties and many more. It is beyond the scope of this tutorial to delve into these effects. At this point, the mere acknowledgement of the existence of these defect states is sufficient.

What are shallow and deep defects and how to find them in defect diagrams?

Charged defects need to be ionized to make the associated electronic carriers available for conduction. If this ionization energy is large, the electronic carriers remain "bound" to the defect site, but if it is small (typically within a few $k_\mathrm{B}T$ ), the bound electronic carriers are excited into conduction states. The energetic location of a CTL relative to a band edge is essentially the thermal energy needed to ionize the defect. For a donor defect (donates electrons), the CTL relative to the conduction band edge is relevant. Similarly, for an acceptor defect (accepts electron/donates holes), the CTL location from the valence band edge is relevant.

How far the CTL is from the relevant band edge determines if a defect is "shallow" or "deep" (and how deep). When a defect has no CTLs inside the band gap, it implies that the defect state(s) lies inside the bands. Such defects are "shallow" defects – they are readily ionized and the associated electronic carriers are available as conduction carriers (also, free carriers). Another case of a shallow defect is when the CTL is close to the corresponding band edge (donor CTL relative to CBM, acceptor CTL relative to VBM), typically within a few $k_\mathrm{B}T$ . The electrons(holes) in the defect states can be easily thermalized (remember thermal energy ~ $k_\mathrm{B}T$ ) into the conduction(valence) band. The schematic below illustrates these two cases for shallow donor defects – one with no CTL inside the band gap and another with CTL close to the conduction band edge.

By extension, it should be now clear that "deep" defects has CTLs far from the corresponding band edges, typically more than few $k_\mathrm{B}T$ . Electronic carriers are trapped at these deep defect sites; such carriers cannot contribute to the electrical conductivity, which is only due to free carriers. Deep defects states are also sometimes referred to as mid-gap states. Next time you see a defect diagram with CTLs far from the band edges, you can quickly conclude that there are deep defects present. However, it is also important to consider the concentration of such defects, which brings us to the next important topic – calculation of defect concentrations.

In rare cases, an extremely deep donor defect will have CTLs lie inside the valence band. Analogously, an extremely deep acceptor will have CTLs inside the conduction band.

What is the effect of elemental chemical potentials on the defect diagram and how is related to the synthesis/growth conditions?

The defect formation energy depends on the elemental chemical potentials ( $\mu_i$ ).

\Delta E_{D,q} = E_{D,q} - E_{H} + \sum_i{n_i \mu_i} + E_{corr} + ~ qE_\mathrm{F}

As discussed above, defect diagrams are typically plotted at a fixed set of $\mu_i$ s, where i are the elements in the chemical phase space. For example, the defect formation energy of ZnO will depend on the chemical potentials of Zn and O. Similarly, for KGaSb4, the defect energetics depend on the chemical potentials of K, Ga, and Sb. Mathematically, changing $\mu_i$ shifts the defect lines vertically in a defect diagram (changes y-axis intercept), but the dominant charge q at a given $E_\mathrm{F}$ and the CTLs remain unchanged. This is schematically shown in the figure below.

The elemental chemical potentials intimately depend on the chemical environment prevalent during growth/synthesis. Formation of defects involve exchange of atoms between the material and external elemental reservoirs. This exchange requires a specific amount of energy, which is determined by the chemical potentials of the elemental reservoirs ( $\mu_i$ ). Here, the chemical conditions during growth dictate $\mu_i$ .

Let us consider the example of ZnO grown under oxygen-rich (O-rich) and oxygen-deficient (O-poor) conditions. Under what conditions would oxygen vacancy formation in ZnO be more favorable? We can arrive at the answer by simply thinking about the process. Oxygen vacancy formation involves the removal of O atoms from ZnO and placing it in the external O reservoir. Under O-rich conditions, this external reservoir already has a high concentration of O, so placing additional O atoms would require a lot more energy. In contrast, under O-poor conditions, where the reservoir has a lower concentration of oxygen, removal of O from ZnO and placing it in this reservoir should be energetically cheaper. Therefore, O vacancy formation is more favorable under O-poor conditions. A detailed mathematical treatment on the effect of elemental chemical potentials and growth conditions on defect energetics will be added later.

How do we calculate the defect concentrations and associated electron and hole concentrations?

The defect concentration at a given temperature is calculated from the Boltzmann probability:

[D_q] = N_s.exp\left(- \frac{\Delta E_{D,q}}{k_\mathrm{b}T} \right)

where $[D_q]$ is the concentration of defect D in charge state q, $N_s$ is the number of atomic sites where D can form, $k_\mathrm{B}$ is the Boltzmann constant, and T is the temperature. The exponential term is the essentially the probability of forming the D, which when multiplied with the volumetric (or areal) site concentration, given the volumetric defect concentration. Typically, defect concentrations are expressed in numbers per unit volume ( $\mathrm{cm}^{-3}$ ).

The equation for calculating defect concentration using Boltzmann probability implicitly includes the configurational entropy contribution. However, it does not include the vibrational entropy contribution, which is generally negligible in most cases. The latter is an important consideration when defect formation involves elements that are gases under standard conditions. In such cases, chemical potential of the elemental (gaseous) reservoir has a strong T dependence.

Looking at the defect diagram, the defect formation energy depends on the Fermi energy ( $E_\mathrm{F}$ ), but it is not clear at which $E_\mathrm{F}$ should the defect concentration be calculated. For this, we need to consider charge neutrality.

The defect and electronic carrier (electron, hole) concentrations are calculated by imposing the condition of overall charge neutrality. This means, for the material to be charge neutral, the total number of positive and negative charges must be balanced. The positive charges include donor defects (+ve slope in defect diagrams) and holes, and negative charges are acceptors (-ve slope in defect diagrams) and electrons.

First, we need to choose a relevant temperature, typically, the synthesis or annealing temperature at which defects form and equilibrate. Next, we calculate the acceptor and donor defect concentrations as a function of $E_\mathrm{F}$ . Mathematically, $\Delta E_{D,q} = c + q.E_\mathrm{F}$ , where c is the y-axis intercept in a defect diagram. Neutral defects do not participate in charge balance. The electron and hole concentrations are given by the density of states g(E) and the Fermi-Dirac distribution f(E). Mathematically, $[e] = \int_{E_C}^{\infty}g(E)f_c(E)dE$ and $[h] = \int^{E_V}_{-\infty}g(E)f_v(E)dE$ , where $E_C$ and $E_V$ are the conduction and valence band extrema, respectively, and $f_c(E)$ and $f_v(E)$ are the Fermi-Dirac distribution functions for the conduction and valence bands, respectively. Here, $f_c(E)$ depends on $E_\mathrm{F}$ , $E_C$ , and $T$ . Similarly, $f_v(E)$ depends on $E_\mathrm{F}$ , $E_V$ , and $T$ . Please refer to standard semiconductor textbooks for equations, including the simplified formula to calculate electron and hole concentrations in the non-degenerate limit using band effective mass (instead of performing the integrals discussed above).

Note how the only unknown in the equations to calculate defect and electron and hole concentrations is $E_\mathrm{F}$ ! By imposing the condition of charge neutrality, we determine the equilibrium Fermi energy ( $E_\mathrm{F,eq}$ ).

\sum_{d}\left( q_d[D_{d}] \right) + [h] = \sum_{a}\left( q_a[D_{a}] \right) +[e]

where a and d denote acceptors and donors, respectively. $q_a$ and $q_d$ are the charge states of the corresponding acceptor and donor defects, respectively. Terms in [] denote concentrations. The left-hand side is the total number of positive charges and the right-hand side the negative charges. The charge neutrality equation is solved self-consistently (often, numerically) to obtain $E_\mathrm{F,eq}$ . Once $E_\mathrm{F,eq}$ is determined, we can calculate the defect and electron and hole concentrations.

There are several open-source Python scripts that are available for numerically solving $E_\mathrm{F,eq}$ . I recommend py-sc-fermi developed by the groups of David Scanlon and Benjamin Morgan. A paper describing the software is available here.

The equilibrium Fermi energy is also often referred to as pinned Fermi level. Strict practitioners of the field may argue that Fermi level pinning is not an appropriate term in this context.

This tutorial focuses on the equilibrium thermodynamics of defect formation. Defects can and do form under non-equilibrium condition, but determination of their concentrations is not straightforward without detailed experimental inputs. It is a specialized topic that is beyond the scope of this guide designed for beginners.

Can we guess the approximate location of the equilibrium Fermi energy? What can we learn?

We can guess, to a first order, the approximate location of the equilibrium $E_\mathrm{F}$ . The exact location must be determined mathematically by self consistently solving the charge neutrality condition, as described above. To locate the approximate equilibrium $E_\mathrm{F}$ , one must identify the "defect lines" (in a defect diagram) corresponding to the lowest-formation energy donor and acceptor defects. Remember that neutral defects do not participate in charge balance. The equilibrium $E_\mathrm{F}$ is, in most cases, in the vicinity of where the lowest-energy donor and acceptor defect lines intersect.

How do we qualitatively interpret the defect diagrams in the context of doping?

We learned above that donor-like defects tend to donate electrons while acceptor-like defects tend to accept electrons (create holes). Rarely does a material contain only donor (or acceptor) defects; often, both donor-like and acceptor-like defects are present. Whether a material is overall n- or p-type, depends on the relative concentrations of donor and acceptor defects. Intuitively, if the concentration of donors (+ve charged defects) is much higher than acceptors (-ve charged defects), one can imagine that the material will be n-type. Mathematically, this can be understood by charge balance discussed above. If the +ve charged defects are in higher concentration, the overall charge balance is restored by creating more electrons (-ve charges) than holes. As such, the material will have a high concentration of electrons, making it n-type. Here, we have discussed the scenario of "self-doping" by the native defects. However, in practical applications, materials are often doped with extrinsic dopants e.g., n-type Si doped with P, n-type ZnO doped with Al.

Materials cannot be extrinsically doped n- or p-type at will (otherwise life would have been easy!). A lot depends on the native defect chemistry. In other words, whether a material can be doped n- or p-type depends critically on the donors and acceptors natively present in material e.g., Si vacancy in Si or O vacancy in ZnO. Native defects, as the name suggests, are inherently present in the material and are not introduced through extrinsic doping. Therefore, one cannot bypass the existence of such native defects. The formation energetics of these native defects determine if a material is n- or p-type (or both, or neither) dopable. Here, dopable means the possibility of extrinsically doping a material. The figure below showcases four different scenarios for materials that are (1) p-type dopable, (2) n-type dopable, (3) both p- and n-type dopable ("ambipolar"), and (4) neither p- nor n-type dopable (typical insulators).

Scenario 1:

Transport

Boltzmann Transport Theory and Thermoelectric Performance

A guide to modeling charge transport properties and thermoelectric performance

Written by Michael Y. Toriyama (MichaelToriyama2024 [at] u [dot] northwestern [dot] edu)

Introduction

Thermoelectric properties in a material can often be described using Boltzmann transport theory under low applied field. The electrical and thermal current densities can be described in terms of the applied electric field and a thermal gradient. The relevant charge transport properties, namely the electrical conductivity, Seebeck coefficient, and the electronic contribution to the thermal conductivity, can be derived from electrical and thermal current densities.

Basic Theory

Thermoelectric Performance

The efficiency of a thermoelectric device, from a material engineering perspective, is parametrized by the figure-of-merit , which can be expressed as

where is the electrical conductivity, is the Seebeck coefficient, and are the electronic and lattice contributions to the thermal conductivity, and is the temperature.

Single Parabolic Band Model

The single parabolic band assumption forms the simplest, first-order description of a material's electronic structure. Often, this is a reasonable approximation in regards to charge transport properties.

We do not delve into the derivation of the transport parameters here. You can refer to the theory section of for more details. Rather, we list the final forms of the transport properties, namely the electrical conductivity (), Seebeck coefficient (), and Lorenz number ():

The electronic thermal conductivity can be expressed as

Using the Fermi-Dirac integral, the carrier concentration can also be expressed as

Two-Band Model

The two-band model described here involves the transport of both electrons and holes. In other words, information about the conduction and valence bands are necessary, as well as the band gap. It is convenient to represent the electron and hole contributions separately (using the single parabolic band model of transport as described above) and derive the total electrical conductivity, Seebeck coefficient, and electronic thermal conductivity from the two contributions.

The important aspect of the two-band model to keep in mind is that the model can be built from the individual contributions of electrons and holes. Therefore, a rational way to write a script for the two-band model would be to write functions relevant for the single-band model and to subsequently use those functions to build the two-band model.

Computational Details

Necessary Modules

Fundamental Constants

We need to define fundamental constants (in SI units) to describe transport properties.

Transport Properties for a Single Band

For example, if we set the following parameters for our model system:

we obtain the following plots:

Transport Properties in a Two-Band Model

When both holes and electrons contribute to thermoelectric transport, the transport coefficients need to be evaluated using the two-band model.

The electrical conductivity is simply the sum of the conductivities of the electrons and holes:

Note that Sigma_SPB is the electrical conductivity from a single parabolic band.

The Seebeck coefficient is a sum of the Seebeck coefficients of the electrons and holes, weighted by their respective conductivities:

Note that Seebeck_SPB is the Seebeck coefficient from a single parabolic band.

The electronic thermal conductivity can be codified as

Note that Lorenz_SPB is the Lorenz number from a single parabolic band.

If we set the following parameters in our model system:

then we obtain the following plots:

Polarization

Introduction to Polarization Calculations

This tutorial describes how to calculate spontaneous polarization using density functional theory and modern theory of polarization.

Spontaneous Polarization

Ferroelectrics play a significant role in various technological applications, including tunable capacitors, non-volatile random access memory devices, and electro-optical data storage. The phenomenon of ferroelectricity often arises due to a structural phase transition from a high-symmetry nonpolar phase to a low-symmetry polar phase as temperature decreases, leading to the spontaneous emergence of polarization.

Ferroelectrics exhibit a distinctive hysteresis loop in the polarization versus electric field curve.

The spontaneous polarization can be experimentally determined by taking half of the change in polarization at zero external fields.

The spontaneous polarization of a material is not directly observable. Instead, it is determined by measuring the change in polarization between two stable configurations of the material.

Following the modern theory of polarization, the polarization, P, of a crystal is defined as,

where e is the charge of the electron, ni is an integer, Ri is a primitive lattice vector, and Ω is the unit cell volume. The second term on the right-hand side of the equation is the quantum of polarization, , a consequence of translation symmetry.

Only differences in the computed polarization on the “same branch” are physically meaningful, so one needs a polar structure and a non-polar reference structure in order to perform the polarization calculations.

The following section describes how to compute spontaneous polarization values using density functional theory and the modern theory of polarization

Getting Polar and non-Polar structure

To compute the spontaneous polarization using the berry phase method implemented in VASP.

1: Relax the polar structure using VASP.

2: Create an anti-polar structure from the relaxed polar structure.

use nebmake.pl script to create 10 images between the polar and anti-polar endpoints with the following command line

4: Use these images to calculate the value of spontaneous polarization. Use the flags LCALPOL = .TRUE. and DIPOL = 0.5 0.5 0.5 .

5: After completing this step collect all the POSCAR and OUTCAR files in a directory with proper indexing. e.g POSCAR-00 .... POSCAR-10, OUTCAR-00....OUTCAR-10

If you are working with many Python packages, it is generally recommended you create a separate environment for each of your packages. For example:

Now using pymatgen utility you can extract the spontaneous polarization values from the POSCAR and OUTCAR files collected in step 5.

The spontaneous polarization can be calculated by using the formula

Dynamics in crystalline solids

First-principles Calculations of Atomic Diffusion in Crystalline Solids

This approach (nudged elastic band method) requires knowledge of diffusion mechanisms. A complementary approach using molecular dynamics can help explore diffusion pathways.

Written by Cheng-Wei Lee (clee2 [at] mines [dot] edu)

Atomic diffusion mechanisms

There are multiple textbook diffusion mechanisms and the vacancy-mediated diffusion is one of the most common diffusion mechanisms in ionic compounds. The associated diffusion coefficient can be written as:

D=\frac{1}{6}fZ_{f}Z_{m}l^{2}C_{D}\Gamma

where $\frac{1}{6}$ is for 3D system ( $\frac{1}{2n}$ and $n=3$ ) and $f$ , $Z_{f}$ , and $Z_{m}$ are geometry factors, which depend on the crystal structure. $l$ is the distance the vacancy travel for each jump, $C_{D}$ is the concentration of participating defects, e.g. vacancies. Lastly, $\Gamma$ is the successful jump frequency,

\Gamma = \nu^{*}e^{-\frac{\Delta E_{m}}{k\mathrm{_{B}}T}}

where $\nu^{*}$ is the attempted jump frequency, $\Delta E_{m}$ is the migration barrier of a diffusing atom (see figure below)

Therefore, the diffusion coefficient for vacancy-mediated diffusion in crystalline solids can be calculated from first principles:

The attempt frequency can be approximated by the local dynamic matrix of the moving atom and the force due to smaller change in moving atom can be calculated using DFT.
The defect concentration can be calculated using the supercell approach (see defect calculations)
The migration barrier can be calculated using nudged elastic band method

Tutorial on Using Nudged Elastic Band (NEB) method

This tutorial is based on VASP 5.4.4 and the VTST package developed and maintained by the Henkelman research group at UT Austin

Written by Cheng-Wei Lee (clee2 [at] mines [dot] edu)

Installing VASP with VTST package

The general information of how to install VASP with VTST can be found at the VTST website. Here is a demonstration of installation VASP 5.4.4 + VTST on NREL EAGLE HPC system:

Load modules needed for VASP

1) gcc/8.4.0  2) comp-intel/2020.1.217 3) intel-mpi/2020.1.217              
4) fftw/3.3.8/intel-impi 5) mkl/2020.1.217

Modify main.F file within the src/ folder of VASP 5.4.4:

Find and Replace the following lines

CALL CHAIN_FORCE(T_INFO%NIONS,DYN%POSION,TOTEN,TIFOR, &
    LATT_CUR%A,LATT_CUR%B,IO%IU6)

with

CALL CHAIN_FORCE(T_INFO%NIONS,DYN%POSION,TOTEN,TIFOR, &
    TSIF,LATT_CUR%A,LAT  T_CUR%B,IO%IU6)

Apply patches (Need to make sure the patch file is located at the correct folder by checking the first few lines of the patch file). The patch is included in the VTST package

patch -p0 < vasp-5.4.4-mpmd.patch

Copy source file from vtstcode5 folders. For example, the version 182 has the following source files:

bbm.F, bfgs.F, chain.F, dynamic.F, fire.F, lanczos.F, neb.F, qm.F, 
cg.F, dimer.F, dynmat.F, instanton.F, lbfgs.F, opt.F, sd.F.

Modify the compiling order within the hidden .objects file within the folder of src/. Note that the added .o file should be placed before chain.o

 ....
 hamil_rot.o \
 bfgs.o dynmat.o instanton.o lbfgs.o sd.o cg.o dimer.o bbm.o \
 fire.o lanczos.o neb.o qm.o opt.o \
 chain.o \
 dyna.o \ 
 ....

compile VASP as it is done typically (e.g. make std)

There are significant differences in the installing procedure for VASP 6, see details on the VTST website.

High-level idea of nudged elastic band method

Only very high-level picture is provided here, further details can be found at the original papers by Henkelman and Jónsson[1][2].

The general goal of NEB calculation is to find the minimal energy path (MEP) between two user-defined local minima on the potential energy surface (PES) of a given system (see figure below). The potential energy surface is a function of configuration coordinate and the relation can be described by semi-empirical force field or electronic-structure method like DFT under the Born-Oppenheimer approximation. Using the 2 dimensional PES as an example, one naive way to find the MEP is to sample a set of evenly distributed 2D grid points and get the PES for each point. In principle, we can find the MEP when the grid is dense enough.

To find the MEP more efficiently, nudged elastic band method (NEB) is commonly used. The method start with an initial discrete pathway constructed by linear interpolation between two end points (A and B sites in the figure). Such discrete linear pathway is generally a good initial guess but by no means has the MEP. In order to find MEP, ionic relaxation is needed but these discrete points will relax into either site A or site B, depending on its initial position on the PES. To make these points "hung" along the path between two local minima, springs are introduced to connect adjacent discrete points and thus the method has the name of "elastic band". The concept of "nudging" is to separate forces along and perpendicular to the tangence of the elastic band. The calculation only take into account spring forces along the tangence of the band and the true force (e.g. forces from DFT) perpendicular to the tangence of the band. As a result, the distance between discrete points will remain roughly the same while the path is nudged toward minimum energy path.

Building on top of the NEB method, the climbing image NEB (CI-NEB) method was introduced to find the saddle point with higher energy resolution. The saddle point is the maximum along the minimum energy path while the minimum along other orthogonal directions. The idea of climbing image is to make the discrete point or the intermediate image that has the higher energy climb to the maximum. In order to achieve this numerically, for this specific image, the sign of the forces along the elastic band is flipped such that the saddle point becomes a local minimum.

Furthermore, in addition to atomic degree of freedom (DOF), and unit cell DOF can be also included and such approach is referred to as the generalized solid-state NEB (SS-NEB). This is a powerful method to find solid-solid phase transformation, ferroelectric switching, etc. Generally speaking, much higher number of images may be needed for SS-NEB since the minimum energy path could be more complicated than phenomena that can already captured by NEB. (SS-NEB generally involves more DOF)

VASP input for NEB calculations

Here's an overview of the steps to perform NEB calculations and detail

Prepare end point/images. Using vacancy-mediated diffusion as an example. The end images are supercell with vacancy on two different but adjacent sites. The end images should be relaxed.
Prepare INCAR and initial intermediate images specific for NEB calculations.
Run NEB first to have a clear intermediate image with highest energy difference from end points. Then run CI-NEB to refine the energy barrier
Extract and analyze the data

If multiple local minima are found along the MEP, the MEP can be broken into smaller pieces and calculated separately.

Preparing intermediate images

The script (nebmake.pl) provided within the VTST package is handy in terms of making linear interpolation between end images.

nebmake.pl POSCAR_i POSCAR_f 5

where i and f indicates initial and final images and 5 is the user-defined number of intermediate images to create. Odd number is generally preferred for systems that are expected to have single peak along the minimum energy path. Once this is done, 5+2 folders (+2 are for end images) will be created with the associated linear interpolated images inside. Note: intermediate images are the configuration coordinate, i.e. positions of atoms, along the path and are embodied by the POSCAR files for VASP program.

To use the script properly, the order of atoms within the VASP POSCAR should be the same for end images. Otherwise, weird initial path will be created and lead to weird results.

VASP+VTST INCAR specific for NEB calculations

Only flags different from normal VASP relaxation are highlighted below

ISIF = 3  (need for SS-NEB; ISF=2 for NEB)
NELMIN = 5 (avoid ionic stpes with only two electronic scf steps)
ICHAIN =0  (for VTST package, this is the flag for NEB)
IMAGES = 5 (the number of intermediate images, end images not included)
SPRINGS = -5  (Spring constant for the elastic band. -5 generally works well)
LCLIMB = .FALSE. (FALSE is default value. Only turn on this flag when there is a clear maximum)
IOPT = 3  (3 for SS-NEB; 3→ 2 → 1 for CI-NEB)
IBRION  = 3 (required for VTST)
POTIM = 0 (required for VTST)
LNEBCELL = .TRUE. (True for SS-NEB; otherwise normal NEB)

Generally, higher the number of images, better sampling of the minimum energy path. Simple MEPs only require 3 images while complicated ones need 10+. Unable to reach force convergence criteria is usually an indication to increase the number of intermediate images.

Restart the simulation

The neb calculation generally requires multiple restart to converge the results. To restart the simulation, like other VASP relaxation calculations, we need to make CONTCAR the new POSCAR. Below is an example bash script to do so, assuming there are 13 intermediate images:

for i in `seq -w 1 13`; do mv $i/CONTCAR $i/POSCAR;done

The scripts provided by the VTST package are very useful! The following are few examples

Data extraction and analysis

nebbarrier.pl   (create neb.dat file which has the information on relative energy, distance between image and forces along the bands)
nebspline.pl (forced based cubic spline interpolation for the MEP)

There is a different set of scripts for solid-state NEB (e.g nebefs.pl and nebsplines.pl)

Calculating Polarization Switching Barrier using SS-NEB

This tutorial is based on VASP compiled with VTST package

Written by Cheng-Wei Lee (clee2 [at] mines [dot] edu)

Ferroelectric or not?

Pyroelectric materials by definition have spontaneous polarization ( $P_s$ ) due to lack of inversion symmetry but whether $P_s$ can be experimentally reversed by an external electric field before dielectric breakdown happens determines if a material is ferroelectric. The critical electric field to reverse $P_s$ is the coercive field ( $E_c$ ). While quantitative prediction of coercive field quires larger scale simulations of domain wall motions, polarization switching barrier ( $\omega_s$ ) at atomic scale can provide qualitative insights since $\omega_s$ describe the required energy density to change $P_s$ . $E_D$

A general workflow to calculate polarization switching barrier

Calculating polarization switching barrier ( $\omega_s$ ) can be broken into two general components:

Finding the polar (+ $P_s$ ) and antipolar (- $P_s$ ) structures of a given material (with polar point group)
Determined the minimum energy path (MEP) between them

Once the MEP is determined, the largest barrier along the path is the $\omega_s$ .

Finding antipolar structure

We can determine the antipolar structure once the nonpolar structure is determined. The nonpolar parent structure can be determined using Group-super group relations[1] and one examplar online tool to find the pseudosymmetry is PSEUDO. Here, a python script to find antipolar structure in wurtzite-type is provided. The script assumes the layered hexagonal structure as the nonpolar structure for wurtzite-type polar structures.

from pylada.crystal import read,write,neighbors
import numpy as np

def myneighbors(strc=None, howmany=12, pos=None, tol=0.2):
    """
    A function that replaces the original pylada
    neighbors function that seems to be broken.
    It still uses the original one, just makes sure
    the output is correct. Also, it returns only the 
    first coordination shell.

    The tolerance is in fraction of the distance to 
    the closest atom and is 20% by default (test first)
    """

    nghs=[]

    # Loop ove pylada neighbors (arbitrary small tolerance)
    for ng in neighbors(strc,howmany,pos,0.1):
        # If empty just add ng and the first one should have the shortest distance 
        if len(nghs)==0:
            nghs.append(ng)
            continue

        # Make sure the distance is within the tolerance
        elif nghs[0][-1]<=ng[-1]<=(1.+tol)*nghs[0][-1]:
            nghs.append(ng)

    return nghs

'''
This script assumes that it is the cations that are moving during polarization switching 
'''

# Anion list; make sure all the anions in the dataset are added 
anions = ['N','O','S','Se']

# read polar structure
strc = read.poscar('POSCAR_p')

# find polar direction

### first find all the cation-to-anion vectors
vectot = []
for atom in strc:
    if atom.type not in anions:
        nghs = myneighbors(strc,4,atom.pos,0.2)  # based on our knowledge for wurtzite-like structure
        for ngh in nghs:
            vectot.append(ngh[1])

### The average should reveal the polar direction for wurtzite-like structure       
tmp_dir = np.average(vectot,axis=0) 
### Get rid of noise
polar_dir = np.round(tmp_dir/max(tmp_dir,key=abs),0)
#print(polar_dir)

shifts=[]
# find distance from the hexagonal plane along polar direction
for atom in strc:
    if atom.type not in anions:  # only look at cations
        nghs = myneighbors(strc,4,atom.pos,0.2)
        vecs = []
        for ngh in nghs:
            if np.linalg.norm(ngh[1]*polar_dir/np.linalg.norm(ngh[1])) < 0.9: # i.e. not in polarization direction, only the neighbors in the polar direction can get large number
                vecs.append(ngh[1])

        assert len(vecs) == 3, "The code only works for wurtzite like structure"  # assert function to catch unexpected situtations
        shifts.append(np.average(vecs,axis=0))  # the average of distance vectors from anion to cation in non-polar direction gives the distance from cation to hexagonal plane
       
shift = np.average(shifts,axis=0)  # get the average value out of all the cations

for atom in strc:
    if atom.type not in anions:  # only move the cations
        atom.pos = atom.pos + shift*2.0   # factor of 2.0 is to get to antipolar structure; factor of 1.0 is roughly hexagonal structure

with open("POSCAR_a", "w") as file: write.poscar(strc, file, vasp5=True)// Some code

Calculating MEP using SS-NEB

Once the end images (polar and antipolar structures) are determined and created, here are the steps to perform the SS-NEB calculations

Fully relax (ionic position, cell size, and cell shape) the end images
Use linear interpolation (cell as well as ionic positions) between end images to get the initial intermediate images
Run SS-NEB calculations till convergence is reached.

nebmake.pl (provided within the VTST package) can be used to do the linear interpolation

Here's an exemplar VASP input file for the SS-NEB calculations:

ADDGRID = .TRUE.
ISPIN = 2 # depends on the system of interest
ISTART = 1
ICHARG = 0
ISYM = 0
LMAXMIX = 4
LVHAR = .FALSE.
LORBIT = 10
LWAVE = .TRUE.
LCHARG = .TRUE.
LVTOT = .FALSE.
LPEAD = .FALSE.
LPLANE = .TRUE.
NELMIN = 5
NPAR = 3 # depends on the total number of requested cores and number of images
NELECT = 288.0 # for a supercell with charged defects
ALGO = Normal
EDIFF = 1e-06
EDIFFG = -0.01
ENCUT = 340.0
PREC = Accurate
NSW = 10  # small number since we want to restart the calculations to have updated FFT grid
IBRION = 3
ISIF = 3
ISMEAR = 0
SIGMA = 0.01
POTIM = 0
ICHAIN = 0
IOPT = 3
SPRING = -5 #default value generally works fine
LCLIMB = .FALSE. 
LNEBCELL = .TRUE. # True is for SS-NEB; False is for NEB
IMAGES = 36 # number of intermediate images

A large number of intermediate images are usually needed to find the SS-NEB pathway between polar and antipolar structures of multinary wurtzite-type materials. Since the current VTST implementation is limited to 99 intermediate images and to speed up the process, segmentation of the pathway is generally needed (see below). We recommend to have 2-3x the number cations within the simulation cell as the initial number of intermediate images and run roughly 100-200 SS-NEB steps to capture the the unconverged MPE. We then separate the path into multiple segments based on the local minima and fully relaxed them to be the new end images. Further segmentation may be needed and we converge each segment independently and connect them at the end.

Calculating polarization switching barrier

Once the minimum energy path is determined by the SS-NEB method, we extract the energy density with respect to the polar structure. It is calculated by the energy difference from the polar structure and normalized by the volume of the polar structure. The polarization barrier is then determined by the largest energy density peak along the switching pathway (see the image below)

While the unit of eV/f.u. is commonly used in the literature, it is not a good choice when comparing materials with different number of cations[2] since the choice leads to spuriously lower barrier for materials with higher number of cations.

Installing Pylada on HPC

Instructions

This simplified tutorial use the Kestrel (HPC@NREL) as an example

Written by Cheng-Wei Lee (clee2 [at] mines [dot] edu)

The very first step is to log into Kestrel. Since there is no external login nodes available yet, you need to use login the SSH gateway first via

ssh -AY [user name]@hpcsh.nrel.gov

[user name] is your user account

-AY is for X11 forwarding and is optional

Once on the ssh gateway, you can log into Kestrel by

ssh [user name]@kl2.hpc.nrel.gov

there are multiple login nodes for Kestrel, e.g. kl1 and kl2

Once on the login nodes of Kestrel, we can load the following modules

source /nopt/nrel/apps/env.sh
module purge
module load anaconda3/2022.05
module load PrgEnv-intel/8.5.0
module swap cray-mpich cray-mpich-abi
module unload cray-libsci
module load intel-oneapi-compilers/2023.2.0
module load intel-oneapi-mpi/2021.10.0-intel 
module load intel-oneapi-mkl/2023.2.0-intel

We only need anaconda3 to set up pylada. The rest are the libraries used to compile VASP and thus needed to run VASP

We need to use the specific mpi, which can be different from the default one

Once anaconda is loaded, we can now use conda to create a python virtual environment

conda create -n [name] python=3.9

[name] is the user defined name for the virtual environment

python 3.9 has been working fine for pylada and pymatgen

Once the virtual environment is set up, we can activate it via

conda activate [name]

This requires the module of anaconda3 being loaded

With the environment set up, we also need git to install pylada

module load git

Now with everything set up, we can follow the instruction on the pylada's github website

(We use pip to install pylada)

pip install git+https://github.com/pylada/pylada-light

One common error message is likely due to lack of compiler

Once the pylada is installed, we need to create a file name "ipython_config.py" with the path of

~/.ipython/profile_default/ipython_config.py

And it needs to have the following content

c = get_config()
c.InteractiveShellApp.extensions = [ "pylada.ipython" ]

In order to submit jobs using pylada, we also need the following file at home directory

~/.pylada

with the following content


vasp_has_nlep = False

################## QDEL definition ################

mpirun_exe = "srun --mpi=pmi2 -n {n} {program}"   # Kestrel requires --mpi=pmi2 for vasp

qdel_exe = "scancel"

qsub_exe = "sbatch"

################### QSTAT definition ##############
def ipython_qstat(self, arg):

  """ Prints jobs of current user. """
  from subprocess import Popen, PIPE
  from IPython.utils.text import SList
  # get user jobs ids
  jobs   = Popen(['squeue', '-u', '[user name]','--format','%.18i %.9P %j %.8u %.2t %.10M %.6D %R'], stdout=PIPE, encoding='utf-8').communicate()[0].split('\n')

  names  = [lines.strip().split()[2] for lines in jobs[1:-1]]
                                                  
  mpps   = [lines.strip().split()[0] for lines in jobs[1:-1]]
                                                  
  states = [lines.strip().split()[4] for lines in jobs[1:-1]]
                                                  
  ids    = [lines.strip().split()[0] for lines in jobs[1:-1]]

  return SList([ "{0:>10} {1:>4} {2:>3} -- {3}".format(id, mpp, state, name)   \
                 for id, mpp, state, name in zip(ids, mpps, states, names)])

##  return SList([ "{0}".format(name)   \
##                 for id, mpp, state, name in zip(ids, mpps, states, names)])	

##################### PBSSCRIPT #####################

pbs_string =  '''#!/bin/bash -x
#SBATCH --account={account}   
#SBATCH --nodes={nnodes}
#SBATCH --ntasks-per-node={ppn}
#SBATCH --export=ALL
#SBATCH --time={walltime}
#SBATCH --job-name={name}
###SBATCH --mem=176G.  # uncomment to requrest specific memory (usually not needed)
#SBATCH --partition={queue}
####SBATCH --qos=high # uncommented if high qos is needed (usually not needed)
###SBATCH -o out
###SBATCH -e err

# Go to the directoy from which our job was launched
cd {directory}

# Make sure the library are loaded properly for vasp
source /nopt/nrel/apps/env.sh
module purge
module load anaconda3/2022.05
module load PrgEnv-intel/8.5.0
module swap cray-mpich cray-mpich-abi
module unload cray-libsci
module load intel-oneapi-compilers/2023.2.0
module load intel-oneapi-mkl/2023.2.0
module load intel-oneapi-mpi/2021.10.0-intel 

source activate [path to virtual environment] 
export OMP_NUM_THREADS=1 #turns off multithreading, needed for VASP

{header}
python {scriptcommand}
{footer}
'''

There are two places requiring modification

[user name] is your account
[path to virtual environment] is the path to your python environment. you can check the python by conda info -e

ALLOY MODELING

Create alloy models

This tutorial showcases how to use Alloy Theoretic Automated Toolkit (ATAT) to create random alloys based on the theory of special quasi-random structure (SQS)

Written by Cheng-Wei Lee (clee2 [at] mines [dot] edu)

Installing Alloy Theoretic Automated Toolkit

Download the source file from the official website (e.g. atat3_36.tar.gz)
upload the file to HPC
decompress the .tar.gz file

tar -zxvf atat3_36.tar.gz

"3_36" is the version number

change the directory into the folder

cd atat

check the makefile within the folder and change the path to binary file (BINDIR) if need

The default path way is ~/bin/ which should be added to $PATH (the default for Kestrel)
Make sure you have g++ , which is provided on Kestrel by default

Then it's time to compile the source code. With the provided makefile, we only need to

make

and if there are no error messages

make install

Sanity check of compiling

echo $?

If it returns 0, it indicates no error message. The binary like mcsqs can be found in ~/bin/

Preparing input for mcsqs

Here, we are using the Monte Carlo (mcsqs) approach to find the special quasi-random structures (SQS) that match the correlation of random distribution

To start with, we need a "lat.in" file which defines the parent structure. The one for wurtzite AlN is provided below as example

3.1284 0.0000 0.0000
-1.5642 2.7093 0.0000
0.0000 0.0000 5.0153
1 0 0
0 1 0
0 0 1
0.667 0.333 0.500 Al=0.08333333,Gd=0.91666667
0.333 0.667 0.000 Al=0.08333333,Gd=0.91666667
0.667 0.333 0.881 N
0.333 0.667 0.381 N

The first six lines define the lattice vectors. If the first three lines form a 3x3 matrix M1 and the next 3 lines for another 3x3 matrix M2. The final lattice vectors in matrix form will be M1xM2. Here's another example for rocksalt AlN

4.983 0.000000 0.000000
0.000000 4.983 0.000000
0.000000 0.000000 4.983
0.000000 0.500000 0.500000
0.500000 0.000000 0.500000
0.500000 0.500000 0.000000
0.500000 0.500000 0.500000 N
0.000000 0.000000 0.000000 Al=0.08333333,Gd=0.91666667

The remaining lines define the basis and the elements. "Al=0.08333333,Gd=0.91666667" is where we can define the Gd content for AlGdN alloys

we need sufficient decimal digit for the fractional numbers (e.g. 8 digit)
The chosen fraction is related to the number of atoms in SQS

Another required input is the files that describe the correlation between atoms. It can be generated using corrdump, which requires lat.in

corrdump -l=lat.in -ro -noe -nop -clus -2=4.5 -3=3.2; getclus

-2 and -3 defines the cutoff radius for what are included as two- and three-body clusters based on the given lat.in structure
The provided numbers for the given wurtizte AlN structure are to include the 1st and 2nd nearest neighbors (NN) for generating the two-body clusters and only the 1st NN for the three-body clusters

corrdump will create two files, sym.out and clusters.out, which will be needed for the monte carlo simulation.

With lat.in, sym.out, and clusters.out, we can run the mcsqs next

Running Monte Carlo simulation to find SQS

The typical prompt for mcsqs is shown below

mcsqs -l=lat.in -n=48

Make sure mcsqs is located in the directory whose path is added to $PATH
-n=48 defines the size of SQS. 48 is chosen since we use the increment of 1/12 and we want larger enough of SQS
We need to perform a convergence test with respect to SQS cell size to minimize the mismatch in correlation functions for those 2-body and 3-body clusters

To fully utilize all the cpus of a given node. Here's an example SLURM script for Kestrel, which has 104 cores/node

#!/bin/bash

#SBATCH --job-name="wz SQS Al=0.25,Sc=0.08333333,Gd=0.66666667"
#SBATCH --account=multiferro 
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=104
#SBATCH --time=04:00:00
#SBATCH -o stdout
#SBATCH -e stderr
#SBATCH -p short

# Go to the directoy from which the job was launched
cd $SLURM_SUBMIT_DIR

for i in {1..104}
do
    seed=$RANDOM
    mcsqs -l=lat.in -n=48 -ip=$i -sd=$seed & 
    echo $i, $seed >> seed_record
done

wait

The idea is to initial parallel MC simulations with different initial seeds.

Analyze the MC simulation results

The goal of the MC simulation is to find a structure that has a distribution of atoms that perfectly match the correlation function of selected 2- and 3-body clusters. During the mcsqs simulation, there are three outputs: mcsqs.log, bestsqs.out, and bestcorr.out

mcsqs.log records the evolution of the MC simulation but only the ones with lower objective functions will be added to it. The mismatches are shown in this file
bestsqs.out is the SQS that best match the correlation function of the clusters
bestcorr.out save the correlation function values for the best SQS

The goal is to find a SQS with not only low objective function but also small mismatches. It is more important to have lower mismatches if two SQS have similar objective function

The bestsqs.out can be converted to VASP POSCAR using the sqs2poscar code

Examples of mismatches in clusters for 48-atom, 72-atom, and 96-atom SQS created for wz-Al1Gd11N12 alloys

# 48-atom
Correlations_mismatch=  0.055556        -0.027778       -0.027778       -0.027778       -0.027778       -0.027778       0.004630        0.004630        -0.078704       -0.078704       0.087963        0.087963
# 72-atom
Correlations_mismatch=  0.027778        -0.027778       0.027778        -0.027778       -0.027778       0.027778        0.032407        0.032407        0.032407        0.032407        0.032407        0.032407
# 96-atom
Correlations_mismatch=  0.013889        -0.027778       0.013889        -0.027778       -0.027778       0.013889        0.004630        0.004630        0.004630        0.004630        0.004630        0.004630

48-atom has some clusters with large mismatch, i.e. number > 0.04
72-atom one is acceptable
96-atom SQS is the best among the three

Boltzmann Transport Theory and Thermoelectric Performance

A guide to modeling charge transport properties and thermoelectric performance

Written by Michael Y. Toriyama (MichaelToriyama2024 [at] u [dot] northwestern [dot] edu)

Introduction

Basic Theory

Thermoelectric Performance

The efficiency of a thermoelectric device, from a material engineering perspective, is parametrized by the figure-of-merit , which can be expressed as

where is the electrical conductivity, is the Seebeck coefficient, and are the electronic and lattice contributions to the thermal conductivity, and is the temperature.

Single Parabolic Band Model

\sigma(\eta) = \frac{8 \pi e}{3} \left(\frac{2 m_e k_B T}{h^2}\right)^{3/2} \mu_w \left(r + \frac{3}{2}\right) F_{r+\frac{1}{2}}(\eta)

S(\eta) = \frac{k_B}{e} \left[ -\eta + \left(\frac{r + \frac{5}{2}}{r + \frac{3}{2}} \right) \frac{F_{r+\frac{3}{2}}(\eta)}{F_{r+\frac{1}{2}}(\eta)} \right]

L(\eta) = \left(\frac{k_B}{e}\right)^2 \left[ \frac{r+\frac{7}{2}}{r + \frac{3}{2}}~ \frac{F_{r+\frac{5}{2}}(\eta)}{F_{r+\frac{1}{2}}(\eta)} - \left\{\frac{r+\frac{5}{2}}{r + \frac{3}{2}}\frac{F_{r+\frac{3}{2}}(\eta)}{F_{r+\frac{1}{2}}(\eta)} \right\}^2 \right]

The electronic thermal conductivity can be expressed as

\kappa_e(\eta) = L(\eta) \sigma(\eta) T

Fundamental constants are needed to describe transport; namely, is the electric charge, is the free electron mass, is the Boltzmann constant, and is Planck's constant.

The transport properties are all dependent on the reduced Fermi level () and the scattering parameter (). The reduced Fermi level is given as

\eta=\frac{E_F}{k_B T}

where is the Fermi level and is the temperature. is referenced to the band edge; in other words, at the band edge, when the Fermi level is in the band, and when the Fermi level is in the band gap. This description holds regardless of whether the majority carriers are electrons (conduction band) or holes (valence band).

is a parameter representing the scattering mechanism of charge carriers. Some common scattering mechanisms are:

Acoustic phonon scattering:
Polar optical phonon scattering:
Ionized impurity scattering:

The Fermi-Dirac integral is given as

F_i(\eta) = \int \frac{\varepsilon^i}{1 + \exp\left(\varepsilon - \eta\right)} d\varepsilon

Using the Fermi-Dirac integral, the carrier concentration can also be expressed as

n(\eta) = 4\pi \left( \frac{2 m_{\rm DOS}^* k_B T}{h^2} \right)^{3/2} F_{1/2}(\eta)

The electrical conductivity is also dependent on the weighted mobility , given as

\mu_w = \frac{e \tau_0}{m_b^*} N_V \left(\frac{m_b^*}{m_e}\right)^{3/2}

where is the effective band mass, and is the valley degeneracy.

Given that the transport properties can be written in closed forms, we can express in a more compartmentalized way:

zT = \frac{S(\eta)^2}{L(\eta) + \frac{(k_B/e)^2}{B \left( r+\frac{3}{2} \right)F_{r+1/2}(\eta)}}

where is the quality factor of the majority carriers defined as

B = \left(\frac{k_B}{e}\right)^2 \frac{8 \pi e}{3} \left(\frac{2 m_e k_B}{h^2}\right)^{3/2} \frac{\mu_{w} T^{5/2}}{\kappa_L}

Notice that , as expressed this way, is determined by two principal quantities of the material: the intrinsic property , and an extrinsic factor .

Two-Band Model

If we reference the Fermi level to the conduction band minimum (i.e. at the CBM), then the reduced Fermi levels of electrons () and holes () are given as

\eta_{n} = \frac{E_F}{k_B T}

\eta_{p} = -\eta_{n} - \frac{E_g}{k_B T}

where is the band gap. The form of can be justified as follows: if is in the conduction band (i.e. ), then from the perspective of holes, is far below the valence band edge (i.e. ). On the other hand, if is below the conduction band edge, then is at the valence band edge, and therefore .

In general, the reference energy for is arbitrary, meaning that can be defined as the CBM, VBM, mid-gap, etc. It is more convenient to define this reference point by the majority carrier type in the material. For example, can be defined as the CBM when the material is n-type, whereas it can be defined as the VBM when the material is p-type. The mathematical formalism for the transport properties can be generalized by using the "majority" and "minority" carrier labels; in other words, we can express the reduced Fermi levels as

\eta_{maj} = \frac{E_F}{k_B T}

\eta_{min} = -\eta_{maj} - \frac{E_g}{k_B T}

where at the edge of the majority carrier band. The thermoelectric transport properties can then be written as

\sigma=\sigma(\eta_{maj})+\sigma(\eta_{min})

S = \frac{\sigma(\eta_{maj})S(\eta_{maj}) + \sigma(\eta_{min})S(\eta_{min})}{\sigma}

L = \frac{ \left[ L_{maj} + S_{maj}^2 \right] \sigma_{maj} + \left[ L_{min} + S_{min}^2 \right] \sigma_{min} }{\sigma} - S^2

As before, the electronic thermal conductivity is given as .

Computational Details

Necessary Modules

Scientific computing can be handled solely by numpy. The Fermi-Dirac integral can be evaluated using the module.

import numpy as np
from fdint import *

From fdint, we use the function fdk(k=i, phi=eta) to evaluate Fermi integrals of the form .

Fundamental Constants

We need to define fundamental constants (in SI units) to describe transport properties.

k = 1.38E-23      # Boltzmann's constant
e = 1.602E-19     # Charge
me = 9.109E-31    # Free electron mass
h = 6.626E-34     # Planck's constant

Transport Properties for a Single Band

As suggested in the , we build the transport model from two single-band transport models (one for the valence band, and the other for the conduction band). The transport properties (, , and ) within a single-band model can be codified as follows:

# Electrical conductivity
def Sigma_SPB(mu_w, eta, T, r):
	return 8. * np.pi * e / 3 / h**3 * (2*me*k*T)**(3./2) * mu_w * (r+3./2) * fdk(k=r+1./2, phi=eta)

# Seebeck coefficient
def Seebeck_SPB(eta, r):
	return k/e * (-eta + (r+5./2)*fdk(k=r+3./2, phi=eta)/(r+3./2)/fdk(k=r+1./2, phi=eta))

# Lorenz number
def Lorenz_SPB(eta, r):
	return (k/e)**2 * ( (r+7./2)*fdk(k=r+5./2, phi=eta)/(r+3./2)/fdk(k=r+1./2, phi=eta) - ( (r+5./2)*fdk(k=r+3./2, phi=eta)/(r+3./2)/fdk(k=r+1./2, phi=eta) )**2 )

Here, eta is the reduced Fermi level andr is the scattering parameter. If needed, the electronic thermal conductivity () can be expressed straightforwardly as

# Electronic thermal conductivity
def Kappa_e_SPB(mu_w, eta, T, r):
        return Lorenz(eta, r) * sigma(mu_w, eta, T, r) * T

For example, if we set the following parameters for our model system:

temperature = 300 	# Temperature (K)
kappa_L = 1 		# Lattice thermal conductivity (W/mK)
band_mass = 0.5*me	# Effective electron mass (kg)
NV = 1 			# Valley degeneracy
r = -1/2 		# Scattering parameter (assumes acoustic phonon scattering)
tau0 = 1e-13		# Scattering prefactor (s)

we obtain the following plots:

Transport Properties in a Two-Band Model

When both holes and electrons contribute to thermoelectric transport, the transport coefficients need to be evaluated using the two-band model.

The electrical conductivity is simply the sum of the conductivities of the electrons and holes:

# Electrical conductivity within the two-band model
def Sigma_TwoBand(mu_w_maj, mu_w_min, eta_maj, eta_min, T, r):
	sigma_maj = Sigma_SPB(mu_w_maj, eta_maj, T, r)
	sigma_min = Sigma_SPB(mu_w_min, eta_min, T, r)
	sigma_total = sigma_maj + sigma_min
	return sigma_total

Note that Sigma_SPB is the electrical conductivity from a single parabolic band.

The Seebeck coefficient is a sum of the Seebeck coefficients of the electrons and holes, weighted by their respective conductivities:

# Seebeck coefficient within the two-band model
# NOTE: Assumes n-type material (electrons = majority)
def Seebeck_TwoBand(mu_w_maj, mu_w_min, eta_maj, eta_min, T, r):
	sigma_maj = Sigma_SPB(mu_w_maj, eta_maj, T, r)
	sigma_min = Sigma_SPB(mu_w_min, eta_min, T, r)

	Seebeck_maj = Seebeck_SPB(eta_maj, r)
	Seebeck_min = Seebeck_SPB(eta_min, r)

	sigma_total = sigma_maj + sigma_min
	Seebeck_total = (-sigma_maj*Seebeck_maj + sigma_min*Seebeck_min) / sigma_total

	return Seebeck_total

Note that Seebeck_SPB is the Seebeck coefficient from a single parabolic band.

The electronic thermal conductivity can be codified as

def Kappa_e_TwoBand(mu_w_maj, mu_w_min, eta_maj, eta_min, T, r):
	
	sigma_maj = Sigma_SPB(mu_w_maj, eta_maj, T, r)
	sigma_min = Sigma_SPB(mu_w_min, eta_min, T, r)

	Seebeck_maj = Seebeck_SPB(eta_maj, r)
	Seebeck_min = Seebeck_SPB(eta_min, r)

	Lorenz_maj = Lorenz_SPB(eta_maj, r)
	Lorenz_min = Lorenz_SPB(eta_min, r)	

	sigma_ratio = sigma_maj / sigma_min

	sigma_total = sigma_maj + sigma_min
	Seebeck_total = (-sigma_maj*Seebeck_maj + sigma_min*Seebeck_min) / sigma_total
	kappa_e_total = (Lorenz_maj*sigma_maj + Lorenz_min*sigma_min)*T + sigma_total*T*(Seebeck_maj**2/(1.+1./sigma_ratio) + Seebeck_min**2/(1.+sigma_ratio) - Seebeck_total**2)

	return kappa_e_total

Note that Lorenz_SPB is the Lorenz number from a single parabolic band.

If we set the following parameters in our model system:

temperature = 300 	# Temperature (K)
kappa_L = 1 		# Lattice thermal conductivity (W/mK)
r = -1/2 		# Scattering parameter (assumes acoustic phonon scattering)
Eg = 0.4 		# Band gap (eV)

band_mass_e = 0.5*me    # Effective electron mass (kg)
band_mass_h = 0.5*me	# Effective hole mass (kg)
tau0_e = 1e-13		# Scattering prefactor for electrons (s)
tau0_h = 1e-13		# Scattering prefactor for holes (s)
NV_e = 1 		# Valley degeneracy for electrons
NV_h = 1 		# Valley degeneracy for holes

then we obtain the following plots:

Tutorial on Using Nudged Elastic Band (NEB) method

This tutorial is based on VASP 5.4.4 and the VTST package developed and maintained by the Henkelman research group at UT Austin

Written by Cheng-Wei Lee (clee2 [at] mines [dot] edu)

Installing VASP with VTST package

The general information of how to install VASP with VTST can be found at the VTST website. Here is a demonstration of installation VASP 5.4.4 + VTST on NREL EAGLE HPC system:

Load modules needed for VASP

1) gcc/8.4.0  2) comp-intel/2020.1.217 3) intel-mpi/2020.1.217              
4) fftw/3.3.8/intel-impi 5) mkl/2020.1.217

Modify main.F file within the src/ folder of VASP 5.4.4:

Find and Replace the following lines

CALL CHAIN_FORCE(T_INFO%NIONS,DYN%POSION,TOTEN,TIFOR, &
    LATT_CUR%A,LATT_CUR%B,IO%IU6)

with

CALL CHAIN_FORCE(T_INFO%NIONS,DYN%POSION,TOTEN,TIFOR, &
    TSIF,LATT_CUR%A,LAT  T_CUR%B,IO%IU6)

Apply patches (Need to make sure the patch file is located at the correct folder by checking the first few lines of the patch file). The patch is included in the VTST package

patch -p0 < vasp-5.4.4-mpmd.patch

Copy source file from vtstcode5 folders. For example, the version 182 has the following source files:

bbm.F, bfgs.F, chain.F, dynamic.F, fire.F, lanczos.F, neb.F, qm.F, 
cg.F, dimer.F, dynmat.F, instanton.F, lbfgs.F, opt.F, sd.F.

Modify the compiling order within the hidden .objects file within the folder of src/. Note that the added .o file should be placed before chain.o

 ....
 hamil_rot.o \
 bfgs.o dynmat.o instanton.o lbfgs.o sd.o cg.o dimer.o bbm.o \
 fire.o lanczos.o neb.o qm.o opt.o \
 chain.o \
 dyna.o \ 
 ....

compile VASP as it is done typically (e.g. make std)

There are significant differences in the installing procedure for VASP 6, see details on the VTST website.

High-level idea of nudged elastic band method

Only very high-level picture is provided here, further details can be found at the original papers by Henkelman and Jónsson[1][2].

VASP input for NEB calculations

Here's an overview of the steps to perform NEB calculations and detail

Prepare end point/images. Using vacancy-mediated diffusion as an example. The end images are supercell with vacancy on two different but adjacent sites. The end images should be relaxed.
Prepare INCAR and initial intermediate images specific for NEB calculations.
Run NEB first to have a clear intermediate image with highest energy difference from end points. Then run CI-NEB to refine the energy barrier
Extract and analyze the data

If multiple local minima are found along the MEP, the MEP can be broken into smaller pieces and calculated separately.

Preparing intermediate images

The script (nebmake.pl) provided within the VTST package is handy in terms of making linear interpolation between end images.

nebmake.pl POSCAR_i POSCAR_f 5

To use the script properly, the order of atoms within the VASP POSCAR should be the same for end images. Otherwise, weird initial path will be created and lead to weird results.

VASP+VTST INCAR specific for NEB calculations

Only flags different from normal VASP relaxation are highlighted below

ISIF = 3  (need for SS-NEB; ISF=2 for NEB)
NELMIN = 5 (avoid ionic stpes with only two electronic scf steps)
ICHAIN =0  (for VTST package, this is the flag for NEB)
IMAGES = 5 (the number of intermediate images, end images not included)
SPRINGS = -5  (Spring constant for the elastic band. -5 generally works well)
LCLIMB = .FALSE. (FALSE is default value. Only turn on this flag when there is a clear maximum)
IOPT = 3  (3 for SS-NEB; 3→ 2 → 1 for CI-NEB)
IBRION  = 3 (required for VTST)
POTIM = 0 (required for VTST)
LNEBCELL = .TRUE. (True for SS-NEB; otherwise normal NEB)

Restart the simulation

for i in `seq -w 1 13`; do mv $i/CONTCAR $i/POSCAR;done

The scripts provided by the VTST package are very useful! The following are few examples

Data extraction and analysis

nebbarrier.pl   (create neb.dat file which has the information on relative energy, distance between image and forces along the bands)
nebspline.pl (forced based cubic spline interpolation for the MEP)

There is a different set of scripts for solid-state NEB (e.g nebefs.pl and nebsplines.pl)

Interpreting Defect and Energy Level Diagrams

A beginner's guide to understanding defect and energy level diagrams obtained from first-principles calculations

Written by Prashun Gorai

Defect Diagrams

How are these defect diagrams plotted?

\Delta E_{D,q} = E_{D,q} - E_{H} + \sum_i{n_i \mu_i} + E_{corr} + ~ qE_\mathrm{F}

Even though only the charge state q with lowest $\Delta E_{D,q}$ is plotted at a given $E_\mathrm{F}$ , a non-zero Boltzmann population of other charge states exists, in reality.

What does the charge state/slope of a defect signify and how is related to doping?

What is the significance of a charge transition level in defect diagram?

What are shallow and deep defects and how to find them in defect diagrams?

In rare cases, an extremely deep donor defect will have CTLs lie inside the valence band. Analogously, an extremely deep acceptor will have CTLs inside the conduction band.

What is the effect of elemental chemical potentials on the defect diagram and how is related to the synthesis/growth conditions?

The defect formation energy depends on the elemental chemical potentials ( $\mu_i$ ).

\Delta E_{D,q} = E_{D,q} - E_{H} + \sum_i{n_i \mu_i} + E_{corr} + ~ qE_\mathrm{F}

How do we calculate the defect concentrations and associated electron and hole concentrations?

The defect concentration at a given temperature is calculated from the Boltzmann probability:

[D_q] = N_s.exp\left(- \frac{\Delta E_{D,q}}{k_\mathrm{b}T} \right)

\sum_{d}\left( q_d[D_{d}] \right) + [h] = \sum_{a}\left( q_a[D_{a}] \right) +[e]

The equilibrium Fermi energy is also often referred to as pinned Fermi level. Strict practitioners of the field may argue that Fermi level pinning is not an appropriate term in this context.

Can we guess the approximate location of the equilibrium Fermi energy? What can we learn?

How do we qualitatively interpret the defect diagrams in the context of doping?

Scenario 1: