Key Considerations for Point Defect Calculations in VASP

Key points to consider during defect calculations.

Posted Apr 13, 2025

By Shibu Meher

10 min read

Point defects – vacancies, interstitials, and substitutional atoms – are the elementary excitations of crystalline materials. They might seem small, but they profoundly influence a material’s electronic, optical, magnetic, and transport properties. Understanding them is crucial for applications ranging from semiconductors and catalysts to batteries and structural alloys.

Density Functional Theory (DFT), particularly as implemented in the Vienna Ab initio Simulation Package (VASP), is a powerful tool for investigating these defects. However, getting reliable and meaningful results requires careful consideration of the setup and analysis. It’s easy to get a number, but much harder to get the right number.

Here are some critical things to keep in mind when embarking on point defect calculations with VASP:

1. The Supercell Approach: Size Matters, A Lot!

Since VASP uses periodic boundary conditions (PBCs), you can’t simulate an isolated defect in an infinite crystal. Instead, you create a supercell – a periodically repeated unit containing the defect.

Problem: PBCs mean the defect interacts with its periodic images in adjacent supercells. This interaction is artificial and needs to be minimized.
Solution: Use a sufficiently large supercell. How large? Large enough that the defect formation energy and relevant properties (like charge density or atomic relaxations) converge with respect to supercell size.
Action: Systematically test convergence by calculating the formation energy (see below) for increasing supercell sizes (e.g., 2x2x2, 3x3x3, 4x4x4 multiples of the primitive cell). This can be computationally expensive but is non-negotiable for accuracy.

2. Convergence is King (Still!)

Just like any DFT calculation, convergence with respect to computational parameters is essential.

ENCUT (Plane-wave Energy Cutoff): Test this carefully for your bulk system first. Use the same, converged ENCUT for your defect calculations. Don’t skimp here.
K-points: Defects break the symmetry of the primitive cell, often requiring a denser k-point mesh for the supercell than what might seem sufficient based on scaling the primitive cell’s mesh. Test k-point convergence for the defect supercell. A Gamma-point-only calculation (KPOINTS file with 1x1x1 mesh) is rarely sufficient, even for large supercells, unless the cell is truly enormous.
Forces: Ensure atomic positions are fully relaxed. Use a strict force convergence criterion (e.g., EDIFFG = -0.01 or smaller for forces in eV/Å).

3. Creating the Defect: Modifying POSCAR

Vacancy: Remove the corresponding atom’s line from the POSCAR file and adjust the total atom count in the line above the coordinates.
Interstitial: Add the atom’s coordinates to the POSCAR file (choosing a sensible starting position, e.g., a high-symmetry interstitial site) and update the atom count.
Substitutional: Replace the coordinates and type of the original atom with the substituting atom in the POSCAR. Update atom types and counts accordingly.
Initial Relaxation: Give the atoms near the defect a small “nudge” or displacement from their ideal lattice sites before starting the relaxation. This helps break symmetry and can prevent the calculation from getting stuck in a local energy minimum representing an unstable, high-symmetry configuration.

4. Handling Charged Defects: The NELECT Tag and Corrections

Defects can trap or donate electrons, becoming charged. Simulating these requires special care.

NELECT Tag: To simulate a defect with charge q, you need to adjust the total number of electrons in the system using the NELECT tag in the INCAR file. NELECT = (electrons in perfect supercell) - q. Remember q is positive for electron removal (positive charge) and negative for electron addition (negative charge).
Charge Neutrality: VASP automatically adds a neutralizing background charge (jellium) when NELECT doesn’t match the total ionic charge. This is artificial.
Finite-Size Corrections: The interaction of the defect charge with the background charge and its own periodic images leads to significant finite-size errors, especially for smaller supercells. These errors scale roughly as 1/L (where L is the supercell dimension). You must apply corrections.
- Common Schemes: Makov-Payne, Freysoldt-Neugebauer-Van de Walle (FNV), Kumagai-Oba (KO). The FNV and KO schemes are generally preferred as they account for both electrostatic potential alignment and charge interactions.
- VASP Implementation: Check VASP documentation for tags like IDIPOL, LDIPOL, LMONO, EFIELD_PEAD. Dedicated post-processing tools (like pycdt or implementing the correction formulas yourself using VASP outputs like the LOCPOT file) are often necessary. VASPsol can also be relevant for implicit solvation environments.

5. Relaxation Strategy: ISIF Tag

When relaxing the defect structure, you usually want the atoms to move, but the supercell lattice vectors to remain fixed (corresponding to the defect embedded in a bulk crystal).

Recommendation: Use ISIF = 2. This allows ions to relax (NSW > 0) but keeps the cell shape and volume constant. Relaxing the cell shape/volume (ISIF = 3 or higher) is generally not appropriate for isolated point defect calculations, as it allows the artificial periodic strain field to relax the simulation cell, which doesn’t happen for a single defect in a macroscopic crystal.

6. Calculating the Formation Energy: The Central Quantity

The key thermodynamic quantity is the defect formation energy (E_form). For a defect X in charge state q, it’s calculated as:

\[E_{form}[X^q] = (E_{tot}[X^q] - E_{tot}[perfect]) - Σ n_i * \mu_i + q * (E_F + E_{VBM} + \Delta V) + E_{corr}\]

Where:

\(E_{tot}[X^q]\): Total energy of the supercell containing the defect X with charge q.
\(E_{tot}[perfect]\): Total energy of the equivalent perfect supercell (calculated with identical k-points and ENCUT).
\(n_i\): Number of atoms of type i added (\(n_i > 0\)) or removed (\(n_i < 0\)) to create the defect.
\(\mu_i\): Chemical potential of species i. This is crucial and non-trivial! It represents the energy reservoir from which atoms are added or to which they are removed. It depends on experimental conditions (e.g., O-rich vs. O-poor). Often determined from calculations of elemental phases or stable compounds under thermodynamic equilibrium constraints.
\(q\): Charge state of the defect.
\(E_F\): Fermi level, treated as a variable, referenced to the Valence Band Maximum (VBM) of the perfect bulk material.
\(E_{VBM}\): Energy of the VBM in the perfect bulk calculation.
\(\Delta V\): Potential alignment correction between the defect supercell and the perfect supercell (essential for charged defects). This term aligns the electrostatic potential (e.g., deep core levels or average potential far from the defect) in the charged supercell calculation with that of the perfect bulk calculation. Often obtained as part of the finite-size correction scheme (FNV, KO).
\(E_{corr}\): The finite-size charge correction term mentioned in point 4.

7. Choosing Chemical Potentials (\(\mu_i\)): The Thermodynamic Environment

The choice of \(\mu_i\) dramatically affects the formation energy. It reflects the abundance of the atomic species.

Example: For an oxygen vacancy (\(V_O\)) in ZnO, you remove one O atom (\(n_O = -1\)). \(μ_O\) can range from the energy of an O atom in an O₂ molecule (O-rich limit) down to a value determined by equilibrium with Zn metal (Zn-rich limit, where \(μ_{Zn} + μ_O = E_{tot}[ZnO~bulk]\)).
Action: Calculate the total energies of relevant elemental phases (e.g., bulk Zn, O₂ molecule) and compounds (e.g., bulk ZnO) to determine the thermodynamically allowed range for \(μ_i\).

8. Functional Choice: Beyond Standard GGA

Standard DFT functionals like PBE (a GGA functional) often severely underestimate band gaps. This impacts:

Defect Level Positions: The energy levels introduced by defects within the band gap might be incorrectly positioned relative to the VBM and CBM.
Formation Energies: Errors in band edge positions directly affect the \(q * (E_F + E_{VBM})\) term, especially for defects whose levels are near the band edges.
Action: Consider using hybrid functionals (like HSE06) or DFT+U for systems where standard GGA fails significantly (e.g., many oxides and semiconductors). These are more computationally expensive but often yield more accurate band gaps and defect levels. Remember to recalculate all energies (perfect bulk, defect cells, reference phases for chemical potentials) with the chosen functional.

9. Symmetry Considerations (ISYM)

Defects often lower the symmetry of the crystal. While VASP can detect symmetry (ISYM tag), sometimes manually setting ISYM = 0 or ISYM = -1 (no symmetry) for defect calculations is safer to ensure that the relaxation explores the true energy minimum without being constrained by potentially incorrect symmetry assumptions. Be mindful this increases computational cost.

Conclusion

Calculating point defects with VASP is a powerful technique, but it’s far from a black box. Careful consideration of supercell size, convergence parameters, charge corrections, relaxation strategies, chemical potentials, and the choice of DFT functional is paramount for obtaining physically meaningful and quantitatively reliable results. Always question your setup, test convergence rigorously, and be critical of your results. Happy calculating!

Key Considerations for Point Defect Calculations in VASP

Key Considerations for Point Defect Calculations in VASP

Further Reading & References:

Trending Tags