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.

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.

4. Perform phase stability calculations to determine range of elemental chemical potentials.

5. Create defect supercells and perform geometry optimization (typically, ionic relaxation).

6. Identify likely sites for interstitials using Voronoi tessellation, and follow specific sequential procedure described below to calculate formation energy of interstitials.

7. Calculate electronic and ionic contributions to dielectric constant.

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

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, $eR_{i}/\Omega$, 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 vtst scripts 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

nebmake.pl POSCAR_polar POSCAR_anti-polar 10 

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

ISYM = 0
NSW = 0
ISTART = 0
SYSTEM = AlN
PREC = Accurate
LCALCPOL = .TRUE.
DIPOL = 0.5 0.5 0.5
LORBIT = 10
ADDGRID = .TRUE.
IBRION = -1
ISMEAR = 0
LWAVE = .TRUE.
SIGMA = 0.01
ALGO = Normal

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 pymatgen utility

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

conda create --name my_pymatgen python
source activate my_pymatgen  # OSX or Linux
activate my_pymatgen  # Windows

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

from pymatgen.core.structure import Structure
from pymatgen.analysis.ferroelectricity.polarization import Polarization, calc_ionic, zval_dict_from_potcar, get_total_ionic_dipole
from pymatgen.io.vasp.inputs import Poscar, Potcar
from pymatgen.io.vasp.outputs import Outcar

num_structs = 10
poscars = [Structure.from_file("POSCAR-" + str(i)) for i in range(num_structs)]
outcars = [Outcar("OUTCAR-" + str(i)) for i in range(num_structs)]
potcar=Potcar.from_file("POTCAR")
for i in range(num_structs):
    outcars[i].zval_dict = zval_dict_from_potcar(potcar)


pol_from_out_struct_method = Polarization.from_outcars_and_structures(outcars,poscars,
    calc_ionic_from_zval = True)
print(pol_from_out_struct_method.get_same_branch_polarization_data(convert_to_muC_per_cm2=True))
~                                                                                                    

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 tutorial on how to perform ss-NEB calculations

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

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 $zT$, which can be expressed as

where $\sigma$ is the electrical conductivity, $S$ is the Seebeck coefficient, $\kappa_e$ and $\kappa_L$ are the electronic and lattice contributions to the thermal conductivity, and $T$ 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 this paper for more details. Rather, we list the final forms of the transport properties, namely the electrical conductivity ($\sigma$), Seebeck coefficient ($S$), and Lorenz number ($L$):

The electronic thermal conductivity can be expressed as

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

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

where $E_F$ is the Fermi level and $T$ is the temperature. $E_F$is referenced to the band edge; in other words, $E_F=0$ at the band edge, $E_F>0$ when the Fermi level is in the band, and $E_F<0$ 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).

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

• Acoustic phonon scattering: $r = -1/2$

• Polar optical phonon scattering: $r=1/2$

• Ionized impurity scattering: $r=3/2$

The Fermi-Dirac integral $F_i(\eta)$ is given as

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

The electrical conductivity $\sigma$ is also dependent on the weighted mobility $\mu_w$, given as

where $m_b^*$ is the effective band mass, and $N_V$ is the valley degeneracy.

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

where $B$ is the quality factor of the majority carriers defined as

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

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 $E_F$ to the conduction band minimum (i.e. $E_F = 0$ at the CBM), then the reduced Fermi levels of electrons ($\eta_n$) and holes ($\eta_p$) are given as

where $E_g$ is the band gap. The form of $\eta_p$ can be justified as follows: if $E_F$ is in the conduction band (i.e. $\eta_n>0$), then from the perspective of holes, $E_F$ is far below the valence band edge (i.e. $\eta_p << 0$). On the other hand, if $E_F$ is $E_g$ below the conduction band edge, then $E_F$ is at the valence band edge, and therefore $\eta_p = 0$.

In general, the reference energy for $E_F$ is arbitrary, meaning that $E_F = 0$ 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, $E_F = 0$ 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 $\eta$ as

where $E_F = 0$ 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 $\kappa_e=L \sigma T$.

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 fdint 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 $F_i(\eta)$.

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 two-band model section, 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 ($\sigma$, $S$, and $L$) 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 ($\kappa_e$) 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:

Written by Prashun Gorai

The formation energy of point defects is calculated using the supercell approach described . 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)

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

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

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

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

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

conda create -n [name] python=3.9

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

conda activate [name]

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

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}
'''

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

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

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 VTST package

patch -p0 < vasp-5.4.4-mpmd.patch

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

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

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

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

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

4. Extract and analyze the data

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.

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

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)

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

Installing Alloy Theoretic Automated Toolkit

1. Download the source file from the official website (e.g. atat3_36.tar.gz)

2. upload the file to HPC

3. decompress the .tar.gz file

tar -zxvf atat3_36.tar.gz

1. change the directory into the folder

cd atat

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

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. The first three lines form a 3x3 matrix M1, which consists of the lattice parameters of the unit cell you want to run SQS with. The following 3 lines form another 3x3 matrix M2. M2 is usually the identity matrix. The final lattice vectors in matrix form will be M1xM2. Here's another example of rocksalt AlN with different M2.

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

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

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

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.

Generate SQS with a specific supercell

mcsqs can generate supercell with very slanted shapes and one way to address this is to specify the supercell via a sqscell.out file but with the caveat that it could require larger supercell to reach desirable mismatches.

Here's the modified flag to run mcsqs:

mcsqs -l=lat.in -rc

and here's an exemplar sqscell.out file for the 3x3x2 supercell

1

3 0 0
0 3 0
0 0 2

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

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

2. bestsqs.out is the SQS that best match the correlation function of the clusters

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

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

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:

1. Finding the polar (+$P_s$) and antipolar (-$P_s$) 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 $\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

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.

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

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)

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.

modification 1 applies to "pylada/vasp/functional.py"

At line 865-868 (see blow) of the python script, there is a IF condition block to ensure that the typical VASP runs are executed properly. Since the file structure are different for NEB calculations from typical VASP calculations, we need to comment these lines out as a temporary solution.

modification 2 applies to "pylada/vasp/keywords.py"

For the two blocks that determine ICHARG and ISTART for continuing VASP calculations, the original Pylada package only consider the main working folder and we need to add one of the working subfolders to the checklist:

original ones at line 702 and 796

modified one for both

Customized Pylada workflow

Here's the customized Pylada workflow that takes the number of NEB iterations and the number of intermediate NEB images as input.

User-defined inputs

Here's an exemplar Pylada submission script to utilize the customized Pylada workflow above. Pylada will prepare the INCAR, KPOINTS, POTCAR for each provided pair of endpoints and interact with HPC system to submit VASP jobs. User still need to create the initial intermediate images.

Here an exmeplar script based on Pylada and vtstscript to generate initial neb intermediate images