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, pylada and pylada-defects. 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.

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:

1. Perform geometry optimization (relaxation) of the unit cell.

2. Calculate the electronic structure with a dense k-point grid, using relaxed unit cell as input.

3. Calculate electronic structure with GW, if band gap correction is needed.

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.

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...

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:

where is for 3D system ( and ) and , , and are geometry factors, which depend on the crystal structure. is the distance the vacancy travel for each jump, is the concentration of participating defects, e.g. vacancies. Lastly, is the successful jump frequency,

where is the attempted jump frequency, 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:

1. 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.

2. The defect concentration can be calculated using the supercell approach (see defect calculations)

3. The migration barrier can be calculated using

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.

3: One can use to create several images in between these two endpoints, for example for AlN wurtzite structure, the interpolation path is such that the middle image looks like hexagonal boron nitride.

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

to get the polarization branch use

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

The relative energy per formula shows the switching barrier from polar to anti-polar structures assuming a smooth homogeneous switching from polar to anti-polar structure via non-polar h-BN type intermediate. In order to get a proper switching barrier one needs to do ss-NEB calculations. One can find this useful on how to perform ss-NEB calculations

Overview

The 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 and the located in beautiful Golden, Colorado.

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

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

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

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

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

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

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:

where is the quality factor of the majority carriers defined as

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

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.

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

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

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

As before, the electronic thermal conductivity is given as .

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

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

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.

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:

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

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:

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 is plotted against the Fermi energy , which is typically referenced to the valence band maximum i.e., . A schematic defect diagram is shown below, where is the band gap.

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 has a linear dependence on . Therefore, the equation for calculating represents an equation of a straight line with slope q.

For a defect D, conventionally only the charge state q with the lowest formation energy at a given 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 has the highest concentration at that . This is shown schematically for a defect.

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. . 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 ), 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 . The electrons(holes) in the defect states can be easily thermalized (remember thermal energy ~ ) 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 . 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.

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 ().

As discussed above, defect diagrams are typically plotted at a fixed set of 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 shifts the defect lines vertically in a defect diagram (changes y-axis intercept), but the dominant charge q at a given 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 (). Here, the chemical conditions during growth dictate .

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:

where is the concentration of defect D in charge state q, is the number of atomic sites where D can form, 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 ().

Looking at the defect diagram, the defect formation energy depends on the Fermi energy (), but it is not clear at which 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 . Mathematically, , 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, and , where and are the conduction and valence band extrema, respectively, and and are the Fermi-Dirac distribution functions for the conduction and valence bands, respectively. Here, depends on , , and . Similarly, depends on , , and . 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 ! By imposing the condition of charge neutrality, we determine the equilibrium Fermi energy ().

where a and d denote acceptors and donors, respectively. and 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 . Once 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 . I recommend developed by the groups of David Scanlon and Benjamin Morgan. A paper describing the software is .

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 . The exact location must be determined mathematically by self consistently solving the charge neutrality condition, as described above. To locate the approximate equilibrium , 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 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:

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:

1. Load modules needed for VASP

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

Find and Replace the following lines

with

1. 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

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

1. 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

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

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.

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 () 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 . 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

1. 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.

2. Prepare INCAR and initial intermediate images specific for NEB calculations.

3. 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

Preparing intermediate images

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

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.

VASP+VTST INCAR specific for NEB calculations

Only flags different from normal VASP relaxation are highlighted below

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:

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

Data extraction and analysis

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 .

A general workflow to calculate polarization switching barrier

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

1. Finding the polar (+) and antipolar (-) structures of a given material (with polar point group)

2. Determined the minimum energy path (MEP) between them

Once the MEP is determined, the largest barrier along the path is the .

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 and one examplar online tool to find the pseudosymmetry is . 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.

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

1. Fully relax (ionic position, cell size, and cell shape) the end images

2. Use linear interpolation (cell as well as ionic positions) between end images to get the initial intermediate images

3. Run SS-NEB calculations till convergence is reached.

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

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)

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

Installing Alloy Theoretic Automated Toolkit

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

This tutorial assumes that users already have the relaxed endpoints for (SS-)NEB calculations, which can be also done in a high-throughput fashion using Pylada scripts like HT_relax.py or any other Python workflow that can automate the VASP calculations. The tutorials cover three topics: (1) modification to existing Pylada package, (2) customized Pylada workflows, and (3) user defined inputs.

Modification to exiting Pylada

Since nudged elastic band (NEB) calculations are performed differently from the typical VASP calculations, i.e., there are multiple working subfolders within a working folder, there are two required modifications to the published Pylada package.

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

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