Skip to content

3.5 Neutral and charged defect formation energies

These are unexecuted teaching inputs and starting settings to test. Original figures are schematics, not computed results. Use licensed VASP and PAW data, replace every placeholder, record the executable version and validate convergence.

3.5.1 Model, units and provenance

Use eV for energy, Å for length and eV/Å for force; 1 kbar = 0.1 GPa. State normalization per atom, molecule, primitive cell or simulation cell. Record PAW identifiers, release, ZVAL, ENMAX and permitted hashes; never redistribute POTCAR. SCF convergence addresses the chosen electronic problem; convergence of the target property requires separate tests.

Neutral and charged defect formation energies — original conceptual schematic

Original schematic. Curves explain concepts; blank data areas await verified learner results. No calculation is claimed.

\[ E_f(D^q)=E_D^q-E_{\mathrm{bulk}}-\sum_i n_i\mu_i+q(E_F+E_{\mathrm{VBM}})+E_{\mathrm{corr}} \]

3.5.2 Worked case: procedure, interpretation and checks

The energy-accounting problem

Removing an atom does not simply subtract two comparable total energies: the removed atom must go somewhere. A defect formation energy is the cost of exchanging atoms and electrons with specified reservoirs. Define ni = Ndefect,i−Nbulk,i, so an added atom has positive ni and a vacancy has ni=−1 for that species. Define q>0 as electrons removed from the neutral defect cell. One explicit convention is

Ef(Dq) = Edef(Dq)−Ebulk−Σi ni μi + q(EVBM+EF) + Ecorr.

Here EF is measured upward from the host valence-band maximum, EVBM is on the chosen aligned electrostatic reference, μi are atomic chemical potentials, and Ecorr includes the chosen finite-size/alignment correction convention. Some implementations separate qΔV from image-charge correction; others include alignment within Ecorr. Do not add it twice. All energies must share composition normalization, functional, PAW datasets, and a compatible energy reference.

Start neutral, then charge

First learn with a neutral vacancy in a single-element crystal. Its formula is Ef(V⁰)=E(N−1)−E(N)+μhost, with μhost obtained per atom from the appropriate bulk reservoir. This is not E(N−1)−(N−1)E(N)/N unless that reservoir choice is exactly intended and all reference conditions match.

  1. Relax the host lattice with the chosen physical model. Construct increasingly large, reasonably isotropic supercells and compute pristine references.
  2. Remove one identified atom, preserving metadata that records its original site. Keep the host supercell lattice fixed for a dilute-defect formation energy and relax internal coordinates. Full cell relaxation corresponds to a different finite-concentration mechanical condition unless carefully treated.
  3. Try plausible spin states and symmetry-breaking local distortions. Defects can have several metastable reconstructions; one SCF/relaxation run is not a global search.
  4. Compute converged static energies for pristine and defect cells using compatible mesh density and electronic settings. Inspect local structure, spin density, and defect-state localization.
  5. Evaluate the formation-energy bookkeeping, then repeat across supercell sizes. Even neutral defects have elastic and wavefunction-overlap finite-size errors.

Charged-defect extension

Determine Nneutral_defect from the actual defective POSCAR and its POTCAR valence counts. Do not subtract q from the pristine supercell's electron count after deleting an atom. VASP uses a compensating background for a charged periodic cell; that does not make the raw energy an isolated-defect formation energy.

For a three-dimensional bulk insulator, choose a suitable correction framework, such as a Freysoldt-type or anisotropic Kumagai–Oba-type approach implemented in a maintained external package. Record dielectric input, defect position, potential sampling, alignment region, correction components, and the resulting size dependence. Use the dielectric response appropriate to the process: relaxed thermodynamic charge states generally involve ionic screening, whereas vertical electronic processes need different treatment. Check the method's assumptions rather than inserting a generic dielectric constant.

Obtain the host VBM and band gap with the same stated electronic-structure method and align electrostatic references carefully. For semilocal DFT, underestimated gaps can seriously affect accessible charge states and transition levels. A charged calculation can place a carrier in an extended band instead of a localized defect state; then a localized-charge correction may not apply. Inspect charge density/localization and band filling.

Chemical potentials, plots, and limits

For a compound AB, reservoir choices satisfy μA+μB=μAB at host equilibrium and inequalities that prevent precipitation of elemental or competing phases. “A-rich” and “B-rich” are defined limiting conditions, not arbitrary offsets. Gas reservoirs introduce temperature/pressure thermochemistry beyond one DFT molecular energy; molecules may also need appropriate spin and systematic reference corrections.

Plot formation energy versus EF over the declared host gap. Each charge-state line has slope q. The lowest envelope identifies the thermodynamically preferred charge state for the chosen reservoirs. A transition ε(q/q′) occurs where two corrected formation-energy lines cross; it is not generally equal to a Kohn–Sham defect eigenvalue. Formation energy controls equilibrium propensity; migration barriers in 8.1 NEB diffusion barriers control kinetics. Concentrations additionally need degeneracy, entropy, charge neutrality, and competing species.

Do not transfer standard 3D bulk correction formulas uncritically to charged slabs, surfaces, 2D materials, or delocalized metallic defects. Vacuum-size dependence and electrostatic boundary conditions become fundamental. A dipole correction by itself is not a universal charged-defect cure.

Exercise

3.5.3 Unexecuted inputs and analysis scaffolds

These are unexecuted teaching inputs and starting settings to test. Original figures are schematics, not computed results. Use licensed VASP and PAW data, replace every placeholder, record the executable version and validate convergence.

3.5.3.1 Input block 1

# Defect relaxation delta; fixed host supercell
IBRION = 2
NSW    = 200
ISIF   = 2
EDIFFG = -0.01
ISYM   = 0
# Set ISPIN/MAGMOM when justified and test competing spin solutions.
# Do not set NELECT for the neutral case unless you have a specific reason.

# Charged case only:
# NELECT = N_neutral_defect - q
# Substitute a verified numeric electron count; this expression is not INCAR syntax.

3.5.5 Technical sources