Density Functional Theory

Density functional theory (DFT) is the most widely used framework for calculating the electronic structure of solids, molecules and nanostructures. The method rests on the idea that the ground‑state properties of a many‑electron system can be expressed uniquely in terms of its electron density ρ(r). In the context of the Certificate in Quantum Espresso and VASP Theory, a solid grasp of the terminology is essential for both setting up calculations and interpreting results. The following exposition defines the principal concepts, explains their practical use, and highlights common difficulties that arise in real‑world simulations. Throughout the text short phrases are emphasized with bold or italic tags; each tag encloses no more than a few words, as required.

---

Kohn‑Sham equations lie at the heart of DFT. They replace the interacting many‑electron problem by a set of single‑particle equations that generate the same ground‑state density. The Kohn‑Sham orbitals ψ_i(r) obey

[−½∇² + V_eff(r)] ψ_i(r) = ε_i ψ_i(r)

Where V_eff(r) is the effective potential comprising the external, Hartree and exchange‑correlation parts. The eigenvalues ε_i are often interpreted as band energies, though only the total energy and density are rigorously defined. In practical calculations the Kohn‑Sham equations are solved iteratively until self‑consistency is achieved.

---

Exchange‑correlation functional (XC functional) encodes all many‑body effects beyond the classical electrostatic interaction. The exact form is unknown; therefore approximations are employed. The most common families are

* Local density approximation (LDA) – the XC energy at a point depends only on the local density. * Generalized gradient approximation (GGA) – includes the density gradient ∇ρ(r). Popular GGA parametrizations are PBE and PW91. * Meta‑GGA – adds the kinetic‑energy density as an additional variable (e.G., SCAN). * Hybrid functionals – mix a fraction of exact Hartree–Fock exchange with a GGA component (e.G., PBE0, HSE06).

Each functional has strengths and weaknesses. LDA typically underestimates lattice constants but gives good bulk moduli; GGA improves structural parameters but may over‑delocalize electrons; hybrids correct band gaps at higher computational cost.

---

Pseudopotential replaces the all‑electron potential of a nucleus and its tightly bound core electrons by a smoother effective potential that acts only on valence electrons. The main advantages are reduced plane‑wave cutoff energies and faster convergence. Pseudopotentials are classified as

* Norm‑conserving – preserve the norm of the all‑electron wavefunction within a chosen radius. * Ultrasoft – relax the norm‑conservation condition, allowing much lower cutoffs. * Projector‑augmented wave (PAW) – reconstructs the all‑electron wavefunction from a pseudo‑wavefunction using augmentation functions.

Quantum Espresso (QE) and VASP both support these types. When selecting a pseudopotential, the user must ensure that the valence configuration matches the chemical environment of interest (e.G., Including semicore states for transition metals).

---

Plane‑wave basis is the default representation for Kohn‑Sham orbitals in periodic codes. A plane wave is defined by a reciprocal‑lattice vector G, and the expansion reads

Ψ_i(r) = Σ_G c_i(G) e^{iG·r}

The completeness of the basis is controlled by an energy cutoff E_cut; all G with ½|G|² ≤ E_cut are retained. Larger E_cut yields more accurate results but increases memory and CPU usage. Convergence tests with respect to E_cut are mandatory for reliable total energies, forces and stress tensors.

---

Brillouin zone is the primitive cell of the reciprocal lattice. Integration over the Brillouin zone is required to obtain quantities such as the total energy, charge density and forces. Because the Brillouin zone is continuous, it is sampled by a finite set of k‑points. The quality of the sampling determines the accuracy of the calculated properties.

---

K‑point sampling defines the discrete grid of reciprocal‑space points used for Brillouin‑zone integration. Common schemes include

* Monkhorst‑Pack grids – uniformly spaced points generated by three integers (n_x, n_y, n_z). * Gamma‑centered grids – include the Γ point at the origin, useful for systems with low symmetry.

In QE the grid is specified by the K_POINTS card, while in VASP it is defined in the KPOINTS file. Convergence with respect to the k‑point density must be checked separately for metals (where the Fermi surface demands fine sampling) and insulators (where a coarser grid may suffice).

---

Self‑consistent field (SCF) cycle is the iterative process that updates the electron density until the input and output densities agree within a prescribed tolerance. The steps are

1. Start from an initial guess for ρ(r) (e.G., Superposition of atomic densities). 2. Build V_eff(r) from the current density. 3. Solve the Kohn‑Sham equations to obtain ψ_i(r) and eigenvalues ε_i. 4. Construct a new density ρ_new(r) from the occupied orbitals. 5. Mix ρ_new with the previous density using a mixing scheme (simple linear mixing, Broyden, Pulay, or Kerker).

The mixing parameters are crucial: Too aggressive mixing may cause divergence, while overly conservative mixing leads to slow convergence. QE allows the user to set the mixing beta and the number of previous steps; VASP offers the ALGO tag for different algorithms (e.G., Normal, Fast, All).

---

Total energy in DFT is the sum of kinetic, external, Hartree, and exchange‑correlation contributions, plus a correction for the pseudopotential. The expression is

E_tot = Σ_i f_i ε_i – ½ ∫∫ ρ(r) ρ(r′)/|r−r′| dr dr′ + E_xc[ρ] + E_ion

Where f_i are occupation numbers. The total energy is the quantity minimized during the SCF cycle; it is also used to compare different crystal structures, adsorption configurations or defect states.

---

Force calculation derives from the Hellmann‑Feynman theorem, which states that the force on an atom α is

F_α = −∂E_tot/∂R_α

In plane‑wave DFT the forces consist of a Pulay term (arising from the incompleteness of the basis) and a contribution from the non‑local part of the pseudopotential. Accurate forces are essential for geometry optimization, molecular dynamics and phonon calculations. Both QE and VASP compute forces automatically; the user must ensure that the energy cutoff and k‑point mesh are converged to obtain forces within a few meV/Å.

---

Geometry optimization (also called structural relaxation) seeks the atomic positions and possibly the cell parameters that minimize the total energy. The algorithm iteratively updates the atomic coordinates according to the forces, often using a quasi‑Newton method (BFGS) or a conjugate‑gradient approach. In QE the ATOMIC_POSITIONS and CELL_PARAMETERS cards define the initial geometry, while the calculation mode “relax” triggers the optimizer. In VASP the INCAR tags IBRION, ISIF and NSW control the optimizer type, stress handling and the number of ionic steps.

---

Band structure represents the eigenvalues ε_i(k) along high‑symmetry paths in the Brillouin zone. It provides insight into electronic transport, optical properties and the nature of the band gap. To compute a band structure, one first performs an SCF calculation to obtain a converged charge density, then a non‑self‑consistent field (NSCF) run on a dense set of k‑points following the desired path. In QE the calculation mode “bands” generates the needed data; VASP uses the “ICHARG = 11” setting for a frozen‑density NSCF run.

---

Density of states (DOS) is the distribution of electronic states per energy interval. The total DOS integrates contributions from all k‑points, while the projected DOS (PDOS) decomposes the DOS onto atomic orbitals or angular momentum channels. DOS is useful for identifying the character of the valence band maximum, the conduction band minimum, and any mid‑gap states introduced by defects. In QE the dos.X utility processes the eigenvalues, while VASP provides the DOSCAR file for post‑processing.

---

Projector‑augmented wave (PAW) method combines the efficiency of pseudopotentials with the accuracy of all‑electron calculations. It defines a smooth pseudo‑wavefunction and augments it with atom‑centered partial waves that restore the true nodal structure near the nucleus. PAW datasets contain the projector functions, augmentation charges and reference all‑electron data. VASP’s default potentials are of PAW type; QE can use PAW via the “upf” format, though the implementation is less mature than the norm‑conserving and ultrasoft options.

---

Hybrid functional calculations mix a fraction α of exact exchange from Hartree–Fock with a GGA exchange term. The most common hybrid is HSE06, which screens the long‑range exchange to reduce computational cost. In VASP the tag “HFSCREEN” controls the screening length, and “AEXX” sets α. In QE, hybrid functionals are accessed through the “EXX” module and require the “EXX_ALPHA” and “EXX_SCREENING_PARAMETER” variables. Hybrid calculations are substantially more expensive because the exchange operator is non‑local and scales as O(N³) with the number of electrons.

---

Van der Waals corrections address the deficiency of standard LDA/GGA functionals in describing dispersion forces. Several schemes are implemented:

* DFT‑D2/D3 – empirical pairwise corrections (Grimme). * Tkatchenko‑Scheffler (TS) – environment‑dependent C6 coefficients. * Non‑local functionals (vdW‑DF, optB88‑vdW) – add a kernel that depends on the electron density.

In QE, the “vdW_corr” variable selects the correction type; VASP uses the “IVDW” tag. Including van der Waals interactions is essential for layered materials, molecular adsorption and weakly bound complexes.

---

Spin polarization allows the treatment of magnetic systems by assigning separate densities for spin‑up and spin‑down electrons, ρ↑(r) and ρ↓(r). The total magnetization M = ∫[ρ↑(r) − ρ↓(r)] dr. In spin‑polarized calculations the Kohn‑Sham equations are solved for each spin channel, and the XC functional must be spin‑dependent (e.G., LSDA, spin‑GGA). In QE the “nspin = 2” flag activates spin polarization; VASP uses “ISPIN = 2”. For antiferromagnetic ordering one defines a magnetic unit cell with appropriate initial magnetic moments via the “starting_magnetization” tag (QE) or “MAGMOM” (VASP).

---

Non‑collinear magnetism extends spin polarization by allowing the spin direction to vary continuously in space. The magnetization is described by a three‑component vector field m(r). This is required for spin‑orbit coupling (SOC) studies and complex magnetic textures such as skyrmions. In QE, the “noncolin = .True.” Flag enables the formalism, and “lspinorb = .True.” Activates SOC. VASP uses “LSORBIT = .TRUE.” Together with “SAXIS” to define the spin quantization axis.

---

Spin‑orbit coupling (SOC) couples the electron’s spin to its orbital motion, introducing relativistic effects that split degenerate bands. SOC is critical for heavy elements, topological insulators and materials with strong anisotropy. The inclusion of SOC increases computational cost because the Hamiltonian becomes a 2×2 matrix for each Kohn‑Sham state. In VASP the “LORBIT” tag controls the projection of orbital character with SOC; QE treats SOC via the relativistic pseudopotential and the aforementioned “lspinorb” flag.

---

Born‑Oppenheimer molecular dynamics (BOMD) propagates ionic positions on the potential energy surface defined by the electronic ground state. At each MD step a full SCF cycle is performed to obtain forces. The time step must be small enough to resolve the fastest vibrational modes (≈ 1 fs). In QE the “md” module implements BOMD; VASP uses the “IBRION = 0” and “SMASS” settings for velocity‑Verlet dynamics. BOMD provides insight into finite‑temperature properties, diffusion mechanisms and phase transitions, but the high cost limits the accessible simulation timescales.

---

Car‑Parrinello molecular dynamics (CPMD) treats electronic and ionic degrees of freedom on an equal footing by introducing a fictitious electronic mass. The equations of motion are integrated simultaneously, avoiding the need for a full SCF at each step. CPMD is efficient for long trajectories but requires careful tuning of the fictitious mass and the time step to maintain adiabatic separation. QE includes a CPMD implementation; VASP does not provide a native CPMD mode.

---

Phonon calculations determine the vibrational spectrum of a crystal. The most common approaches are

* Finite‑difference method – displace atoms in a supercell and compute forces to build the dynamical matrix. * Density‑functional perturbation theory (DFPT) – directly evaluates the linear response of the electron density to atomic displacements.

DFPT is implemented in QE through the “ph.X” program and yields phonon frequencies, eigenvectors and electron‑phonon coupling constants without constructing large supercells. VASP can perform finite‑difference phonons via the “IBRION = 8” setting or interface with external tools such as Phonopy. Phonon data are used to assess dynamical stability, calculate thermal conductivity and evaluate zero‑point energy corrections.

---

Electron‑phonon coupling quantifies the interaction between electrons and lattice vibrations. It is essential for understanding superconductivity, carrier mobility and temperature‑dependent band structures. In DFPT the electron‑phonon matrix elements g_{mn,ν}(k,q) are computed for electronic states m,n and phonon mode ν. The Eliashberg spectral function α²F(ω) can be derived from these matrix elements and integrated to obtain the superconducting transition temperature via the McMillan‑Allen‑Dynes formula. QE’s “epw.X” module provides a workflow for electron‑phonon calculations; VASP users typically rely on external packages such as EPW or Wannier90.

---

Wannier functions are localized orbitals obtained by a unitary transformation of Bloch states. Maximally localized Wannier functions (MLWFs) enable accurate interpolation of band structures, calculation of Berry phases, and construction of tight‑binding models. The workflow usually involves:

1. Perform a dense‑grid NSCF calculation to generate Bloch states. 2. Use Wannier90 to select initial projections (e.G., D‑orbitals on transition metals). 3. Optimize the spread functional to obtain MLWFs.

Both QE and VASP can export the necessary wavefunction data for Wannier90. Wannier interpolation reduces the cost of dense‑grid electronic‑structure calculations, making it feasible to evaluate transport coefficients and topological invariants.

---

Berry phase is a geometric phase accumulated by a wavefunction upon adiabatic transport around a closed loop in parameter space. In solids it underlies polarization, orbital magnetization and topological invariants such as the Chern number. The modern theory of polarization expresses the macroscopic polarization as a Berry phase of the occupied Bloch states. QE implements Berry‑phase polarization via the “electric field” module; VASP uses the “LCALCPOL” tag. Calculations of the Z₂ topological invariant for time‑reversal‑invariant insulators also rely on Berry‑phase techniques.

---

Charge density analysis provides insight into bonding, oxidation states and electronic localization. Common tools include

* Bader charge analysis – partitions space into atomic basins defined by zero‑flux surfaces of ∇ρ(r). * Mulliken population analysis – projects the density onto atomic orbitals (more natural for localized basis sets). * Löwdin charges – orthogonalize the basis before population analysis.

In plane‑wave calculations Bader analysis is typically performed on a fine real‑space grid extracted from the CHGCAR (VASP) or charge density file (QE). The resulting charges help interpret catalytic activity, defect charge states and ionic character.

---

Defect formation energy quantifies the thermodynamic cost of creating a vacancy, interstitial or substitutional impurity.

E_f = E_defect – E_bulk + Σ_i n_i μ_i + q(E_F + ε_VBM) + ΔV

Where E_defect and E_bulk are the total energies of the defective and pristine supercells, n_i are the numbers of atoms added or removed, μ_i their chemical potentials, q the charge state, E_F the Fermi level, ε_VBM the valence‑band maximum, and ΔV the potential alignment correction. Accurate defect calculations require large supercells, careful k‑point convergence, and correction schemes for finite‑size electrostatic interactions (e.G., Makov–Payne, Freysoldt). Both QE and VASP can be used to compute E_defect; VASP’s “NELECT” and “LVTOT” tags facilitate charged‑defect setups.

---

Surface slab models are employed to study adsorption, catalysis and surface reconstructions. A slab consists of a finite number of atomic layers separated by a vacuum region thick enough to prevent spurious interactions between periodic images. Convergence tests involve varying the slab thickness and vacuum size, as well as the k‑point sampling perpendicular to the surface. The dipole correction (“dipole = .True.” In QE, “LDIPOL = .TRUE.” In VASP) removes artificial electric fields caused by asymmetric terminations. Adsorption energies are computed as

E_ads = E_slab+adsorbate – E_slab – E_adsorbate

With appropriate reference states for the adsorbate (e.G., Gas‑phase molecule). Spin polarization and van der Waals corrections are often essential for realistic surface chemistry.

---

Band gap problem refers to the systematic underestimation of electronic band gaps by semilocal LDA/GGA functionals. The origin lies in the derivative discontinuity of the exact XC functional, which semilocal approximations miss. Remedies include

* Hybrid functionals – improve gaps but increase cost. * GW approximation – a many‑body perturbation method that corrects quasiparticle energies. * DFT+U – adds an on‑site Hubbard term to localized d or f states, partially opening the gap.

In VASP the “ALGO = GW0” or “ALGO = GW” settings invoke a GW calculation; QE requires external GW codes (e.G., Yambo) interfaced with its output. DFT+U is accessed via the “LDAU = .TRUE.” Flag (VASP) or “lda_plus_u = .True.” (QE). Understanding the limitations of each approach is crucial for interpreting electronic spectra.

---

Strong correlation arises when electron–electron interactions dominate over kinetic energy, leading to phenomena such as Mott insulators, charge ordering and magnetism that standard DFT fails to capture. Methods to treat strong correlation include

* DFT+U – adds a Hubbard‑like penalty to selected orbitals. * Dynamical mean‑field theory (DMFT) – solves a local impurity problem self‑consistently with DFT. * Hybrid functionals with a larger fraction of exact exchange.

The DFT+U formalism requires specifying the on‑site Coulomb (U) and exchange (J) parameters, often obtained from constrained DFT or experimental fits. In VASP the “LDAUVAL” tag sets U, while QE’s “Hubbard_U” variable does the same. DMFT implementations exist as external packages (e.G., TRIQS) that read charge densities from QE or VASP and return a self‑energy for further iterations.

---

Convergence criteria are thresholds that determine when an SCF or geometry optimization is considered complete. Typical criteria include

* Energy tolerance – change in total energy below 10⁻⁵ Ry (QE) or 10⁻⁶ eV (VASP). * Force tolerance – maximum force component below 0.01 EV/Å. * Stress tolerance – for cell relaxations, stresses below 0.5 KBar.

Users must balance strict criteria (ensuring accurate results) against computational expense. For high‑throughput screening, looser tolerances may be acceptable for initial screening, followed by tighter settings for promising candidates.

---

Parallelization strategies differ between QE and VASP. QE uses MPI and OpenMP; the “-np” flag distributes k‑points and plane‑wave workload across processors. VASP employs a hybrid MPI/OMP scheme with “NPAR” controlling the division of the plane‑wave basis. Efficient parallel scaling requires matching the number of processors to the size of the k‑point mesh and the plane‑wave cutoff. For large supercells, a “band parallel” approach (splitting over Kohn‑Sham states) can improve performance.

---

Input file structure varies between the two codes. QE uses a single text file with namelists (e.G., &Control, &system, &electrons) and cards (ATOMIC_SPECIES, ATOMIC_POSITIONS, K_POINTS). VASP relies on a set of fixed‑name files:

* INCAR – controls the calculation parameters (e.G., ENCUT, ISMEAR, IBRION). * POSCAR – contains lattice vectors and atomic positions. * POTCAR – concatenated pseudopotential data for each element. * KPOINTS – defines the k‑point grid.

Understanding the mapping between QE namelists and VASP tags is helpful when translating a workflow from one code to the other. For example, QE’s “ecutwfc” corresponds to VASP’s “ENCUT”, while “occupations = smearing” maps to “ISMEAR” and “SIGMA”.

---

Smearing methods are employed to improve convergence for metallic systems by smoothing the occupation of electronic states near the Fermi level. Common schemes are

* Methfessel‑Paxton – order‑1 smearing, useful for total‑energy calculations. * Gaussian – simple broadening, suitable for density of states. * Fermi‑Dirac – physically motivated, used when temperature effects are of interest.

In QE the “smearing” variable selects the method and “degauss” sets the width (in Ry). VASP uses “ISMEAR” (−5 for tetrahedron, 0 for Gaussian, 1 for Methfessel‑Paxton) and “SIGMA” (in eV). The smearing width must be chosen carefully: Too large a value can artificially raise the total energy, while too small a value may hinder SCF convergence.

---

Electric field and polarization calculations simulate the response of a material to an external bias. In QE the “electric field” module adds a sawtooth potential and computes the Berry‑phase polarization. VASP implements a finite‑field method via the “EFIELD” tag and the “LCALCPOL” flag for Berry‑phase analysis. These capabilities enable the study of ferroelectric switching, dielectric constants and piezoelectric coefficients.

---

Time‑dependent DFT (TDDFT) extends DFT to excited‑state properties by solving the linear‑response equations for the density‑density response function. TDDFT provides optical absorption spectra, exciton binding energies and dynamic polarizabilities. QE’s “turboTDDFT” module implements a Sternheimer approach, while VASP offers the “LOPTICS” tag to compute the frequency‑dependent dielectric function using the independent‑particle approximation. For strongly bound excitons, many‑body methods such as the Bethe‑Salpeter equation (BSE) are required, often interfaced with GW calculations.

---

Charge‑state corrections are necessary when modeling charged defects in periodic supercells. The spurious interaction between the defect charge and its periodic images leads to an error that scales inversely with the supercell size. Correction schemes include

* Makov‑Payne – analytic correction based on the supercell’s dielectric constant. * Freysoldt–Neugebauer–Van de Walle (FNV) – aligns the electrostatic potential far from the defect. * Kumagai–Oba – extends FNV to anisotropic materials.

Both QE and VASP users must apply these corrections manually, typically by post‑processing the electrostatic potential (e.G., LOCPOT in VASP) and extracting the alignment term.

---

Finite‑size effects affect any calculation performed in a supercell, whether for defects, surfaces or molecular adsorbates. The key sources of error are

* Interaction between periodic images (electrostatic and elastic). * Brillouin‑zone sampling that becomes coarse as the cell grows.

Mitigation strategies involve increasing the supercell dimensions, using dipole corrections for slab geometries, and employing extrapolation techniques (e.G., Plotting formation energy versus 1/L where L is the supercell length). For charged defects, the aforementioned charge‑state corrections are essential.

---

High‑throughput screening leverages automated DFT workflows to evaluate large material libraries for properties such as stability, band gap, and catalytic activity. Platforms such as Aflow, Materials Project and OQMD rely on standardized QE or VASP input files, robust convergence criteria, and databases of pre‑computed pseudopotentials. Automation scripts typically generate a series of INCAR or &control settings, submit jobs to a scheduler, parse the OUTCAR or output files, and store results in a relational database. Key challenges include handling failures gracefully, ensuring consistent convergence across diverse chemistries, and managing the computational cost of hybrid or GW calculations.

---

Machine‑learning potentials are emerging as a bridge between DFT accuracy and molecular‑dynamics efficiency. Training data are generated from DFT calculations (often using QE or VASP) that provide energies, forces and stresses for a representative set of atomic configurations. The resulting models (e.G., SNAP, GAP, DeepMD) can reproduce DFT‑level accuracy at a fraction of the cost, enabling simulations of thousands of atoms over nanosecond timescales. Careful selection of training structures, inclusion of diverse chemical environments, and validation against out‑of‑sample DFT data are critical for reliable predictions.

---

Visualization tools help interpret the output of DFT calculations. Common formats include

* XSF – extended structure format used by XCrySDen. * POSCAR/CONTCAR – VASP’s position files, readable by VESTA. * Cube files – volumetric data (charge density, electrostatic potential) for visualization in programs such as VMD or Chimera.

Quantum Espresso provides the “pp.X” utility to generate charge density, potential and wavefunction cube files. VASP’s “CHGCAR” and “ELFCAR” serve the same purpose. Visual inspection of charge accumulation at defect sites, bonding orbitals, or surface states often reveals insights not captured by numeric tables alone.

---

Benchmarking and validation are essential steps before trusting DFT predictions. Standard benchmark sets (e.G., The G2 set for molecules, the SSSP library for solids) provide reference data for total energies, lattice constants and bulk moduli. Comparing calculated values against experimental measurements or higher‑level theories (Quantum Monte Carlo, coupled‑cluster) quantifies the accuracy of a chosen functional and pseudopotential combination. Consistent under‑ or over‑estimation trends can be corrected empirically, but awareness of the underlying approximations remains crucial.

---

Common pitfalls encountered by new users include

1. **Insufficient k‑point density** – leads to noisy total energies and inaccurate forces. 2. **Too low plane‑wave cutoff** – causes basis‑set incompleteness, especially for hard pseudopotentials. 3. **Incorrect pseudopotential choice** – mismatched valence configurations produce erroneous chemistry. 4. **Neglecting spin polarization** – for magnetic materials, results may be qualitatively wrong. 5. **Improper smearing width** – can mask metallic behavior or introduce artificial metallicity. 6. **Ignoring dipole corrections** – for asymmetric slabs, spurious electric fields distort adsorption energies. 7. **Overlooking convergence of forces** – geometry optimizations may stop prematurely, leaving residual stresses.

Addressing these issues typically involves systematic convergence tests, careful inspection of intermediate outputs (e.G., Charge density plots, band structures), and consultation of the extensive documentation provided with QE and VASP.

---

Practical example: Bulk silicon

1. **Select pseudopotential** – use a PAW dataset for Si that includes 3s and 3p valence. 2. **Set cutoff** – start with 30 Ry (≈ 400 eV) and increase in 5 Ry steps until total energy changes < 1 meV/atom. 3. **Define k‑point mesh** – a 8×8×8 Monkhorst‑Pack grid is typical; verify convergence by testing 10×10×10. 4. **Choose functional** – PBE is a common GGA choice; for more accurate lattice constants, consider SCAN. 5. **Run SCF** – monitor the energy and charge residuals; adjust mixing if convergence stalls. 6. **Optimize geometry** – enable “relax” mode; converge forces below 0.01 EV/Å. 7. **Compute band structure** – perform an NSCF run along Γ‑X‑W‑K‑Γ‑L‑U‑W. 8. **Extract DOS** – use the DOS utility to obtain total and projected contributions, confirming the indirect gap.

Repeating the same workflow in VASP involves setting ENCUT to the same value, ISMEAR = −5 (tetrahedron), and IBRION = 2 for relaxation. The resulting lattice constant differs by less than 0.5 % Between the two codes, illustrating the consistency of well‑converged calculations.

---

Practical example: Adsorption of CO on Pt(111)

1. **Build slab** – create a 4‑layer Pt(111) slab with a 15 Å vacuum. Use a (2×2) surface supercell. 2. **Select pseudopotentials** – PAW for Pt and C, O. Include semicore states for Pt if high accuracy is needed. 3. **Apply dipole correction** – enable “dipole = .True.” (QE) or “LDIPOL = .TRUE.” (VASP) because the slab is asymmetric after CO adsorption. 4. **Choose functional** – PBE‑D3 to capture dispersion between CO and the metal surface. 5. **K‑point sampling** – 6×6×1 for the (2×2) cell; verify convergence by testing 8×8×1. 6. **SCF convergence** – tighten energy tolerance to 10⁻⁶ Ry, use a mixing beta of 0.4. 7. **Geometry optimization** – relax the top two Pt layers and the CO molecule; fix the bottom layers to mimic bulk. 8. **Compute adsorption energy** – subtract the energies of the clean slab and isolated CO (computed in a large box).

Typical challenges include the need for a dense k‑point mesh to resolve the metal’s Fermi surface, and the sensitivity of the adsorption energy to the chosen vdW correction. Benchmarking against experimental heats of adsorption helps validate the chosen setup.

---

Practical example: Vacancy formation in MgO

1. **Supercell construction** – generate a 2×2×2 cubic supercell (64 atoms). 2. **Introduce vacancy** – remove one O atom, creating a neutral oxygen vacancy. 3. **Charge corrections** – for the neutral vacancy, no correction is needed; for charged states (e.G., V_O^{2+}), apply the FNV scheme. 4. **Functional choice** – use PBE‑U with U = 5 eV on O‑2p to improve the band gap. 5. **K‑points** – a 2×2×2 grid suffices for this large cell; verify convergence by comparing with a 3×3×3 grid. 6. **SCF and relaxation** – relax all atomic positions; monitor forces below 0.02 EV/Å. 7. **Formation energy calculation** – use the chemical potential of O derived from O₂ gas (including the DFT‑GGA overbinding correction).

The resulting formation energy can be compared with experimental defect formation enthalpies, providing insight into the reliability of the DFT+U approach for oxides.

---

Advanced topic: GW quasiparticle corrections

1. **Initial DFT run** – perform a converged PBE calculation with a high plane‑wave cutoff (e.G., 80 Ry) and dense k‑point sampling. 2. **Prepare wavefunctions** – write the Kohn‑Sham orbitals to disk (VASP: “LWAVE = .TRUE.”). 3. **Run GW** – in VASP, set “ALGO = GW0” for a one‑shot GW calculation; specify the number of empty bands (NBANDS) to ensure sufficient unoccupied states. 4. **Self‑consistency** – optionally iterate the GW step (GW₀ → GW) to improve quasiparticle energies. 5. **Band gap extraction** – compare the GW corrected band gap with the DFT value; note the typical increase of 0.5–1.5 EV for semiconductors.

The GW method is computationally demanding, scaling as O(N³) with the number of bands, and requires careful convergence with respect to the number of empty states and the dielectric cutoff. Nevertheless, it provides the most reliable first‑principles band gaps for a wide range of materials.

---

Advanced topic: DFT+DMFT workflow

1.