Electronic Structure Methods

Electronic structure methods form the foundation of modern computational materials science, providing a bridge between quantum mechanics and observable material properties. In the context of the Certificate in Quantum Espresso and VASP Theo…

Electronic Structure Methods

Electronic structure methods form the foundation of modern computational materials science, providing a bridge between quantum mechanics and observable material properties. In the context of the Certificate in Quantum Espresso and VASP Theory, a solid grasp of the terminology is essential for both practical simulations and the interpretation of results. The following exposition presents the most frequently encountered terms, organized thematically, and illustrated with examples and typical challenges faced by users of these two leading first‑principles codes.

Wavefunction and orbital concepts lie at the heart of any electronic‑structure calculation. In the Kohn‑Sham formulation of density functional theory (DFT), the many‑electron problem is mapped onto a set of non‑interacting particles moving in an effective potential. The resulting single‑particle entities are the Kohn‑Sham orbitals, denoted ψ_i(r). Their squared modulus yields the electronic density, ρ(r)=∑_i f_i|ψ_i(r)|², where f_i are the occupation numbers. The total energy is a functional of this density, and the self‑consistent solution of the Kohn‑Sham equations provides both the ground‑state density and the associated eigenvalues ε_i.

The exchange‑correlation functional (XC functional) encapsulates all many‑body effects beyond the classical Hartree term. Common approximations include the local density approximation (LDA), which assumes the XC energy at each point depends only on the local density, and the generalized gradient approximation (GGA), which incorporates the density gradient. Popular GGA forms are PBE (Perdew‑Burke‑Ernzerhof) and PBEsol, the latter being optimized for solids. Hybrid functionals, such as HSE06, mix a fraction of exact Hartree‑Fock exchange with a GGA functional, improving band‑gap predictions at increased computational cost.

In plane‑wave based codes like Quantum Espresso and VASP, the electronic wavefunctions are expanded in a set of plane waves, ψ_i(r)=∑_G c_i(G) e^{iG·r}, where G are reciprocal‑lattice vectors. The completeness of the basis is controlled by a single parameter, the plane‑wave cutoff (E_cut). A higher cutoff yields more accurate results but increases the number of plane waves dramatically, thus raising memory and CPU requirements. Users typically perform convergence tests, varying E_cut until total‑energy differences fall below a chosen tolerance (e.G., 1 MeV/atom).

Because the core electrons are tightly bound and require a very high cutoff for accurate representation, most calculations employ pseudopotentials. A pseudopotential replaces the all‑electron potential with a smoother effective potential that reproduces the scattering properties of the valence electrons while removing the need to explicitly treat core states. Several families exist: Norm‑conserving pseudopotentials enforce the same norm for pseudo‑ and all‑electron wavefunctions; ultrasoft pseudopotentials relax this constraint, allowing a lower cutoff; the projector‑augmented‑wave (PAW) method combines the advantages of both, reconstructing the all‑electron wavefunction from the pseudo‑wavefunction. VASP’s PAW datasets are widely used, while Quantum Espresso offers both norm‑conserving and ultrasoft pseudopotentials in its PSP library.

The Brillouin zone (BZ) is the primitive cell of reciprocal space. Since periodic systems have infinitely many k‑points, the electronic structure must be sampled at a finite set. The most common scheme is the Monkhorst‑Pack grid, defined by a set of integers (n₁,n₂,n₃) that specify how many k‑points lie along each reciprocal‑lattice direction. The k‑point density directly influences the accuracy of total‑energy and forces; metallic systems generally require denser sampling because of partially occupied states near the Fermi level. A typical convergence test might increase the grid from 4×4×4 to 8×8×8, monitoring changes in the total energy.

The self‑consistent field (SCF) cycle iteratively solves the Kohn‑Sham equations. Starting from an initial guess for the density (often a superposition of atomic densities), the Hamiltonian is constructed, eigenvalues and eigenvectors are obtained, a new density is generated, and the process repeats until convergence criteria are met. Convergence is usually assessed by the change in total energy (ΔE) and the root‑mean‑square (RMS) change in the density or potential. Typical thresholds are ΔE < 10⁻⁵ Ry and RMS < 10⁻⁶ e. In practice, achieving SCF convergence can be challenging for systems with near‑degenerate states, metallic occupations, or strong magnetic interactions. Techniques such as mixing schemes (Pulay, Broyden) and smearing methods (Methfessel‑Paxton, Gaussian, Fermi‑Dirac) are employed to stabilize the iteration.

Smearing introduces a finite electronic temperature to smooth the occupation function, alleviating convergence issues in metals. The smearing width (σ) is a tunable parameter; too large a σ can artificially broaden the density of states, while too small a σ may cause SCF oscillations. After convergence, the smearing contribution to the total energy is removed (extrapolated to zero temperature) to obtain the physical ground‑state energy.

The total energy itself is a central observable, from which many other properties derive. For structural relaxations, the forces on atoms, F_i = −∂E/∂R_i, must be computed accurately. In plane‑wave codes, forces are obtained analytically from the Hellmann‑Feynman theorem, supplemented by Pulay corrections when the basis set depends on atomic positions (as is the case for plane waves with a fixed cutoff). The stress tensor, σ_αβ, is likewise derived from the energy derivative with respect to strain, enabling cell‑shape optimization and the calculation of elastic constants.

Geometry optimization proceeds by iteratively updating atomic coordinates (and sometimes cell vectors) to minimize the forces and stresses. Algorithms such as the conjugate‑gradient method, BFGS (Broyden‑Fletcher‑Goldfarb‑Shanno), and the damped dynamics scheme are implemented in both Quantum Espresso’s ‘vc‑relax’ and VASP’s ionic relaxation. Convergence criteria for the ionic steps typically involve a maximum force threshold (e.G., 0.01 EV/Å) and a stress tolerance (e.G., 0.5 Kbar). For complex systems—surfaces, interfaces, or large supercells—relaxations can be computationally demanding, and careful pre‑conditioning (e.G., Fixing deep layers) often helps.

The band structure is a plot of eigenvalues ε_i(k) along high‑symmetry directions in the BZ. To generate a band structure, one first performs a self‑consistent calculation on a dense k‑mesh to obtain the converged charge density, then carries out a non‑self‑consistent calculation (often called ‘nscf’) on the desired path, using the same density. The resulting eigenvalues are plotted, revealing features such as band gaps, band crossings, and effective masses. In VASP, the ‘IBRION = -1’ and ‘LWAVE = .FALSE.’ Settings are common for static runs, while the ‘KPOINTS’ file defines the path; Quantum Espresso requires separate input files for the SCF and band‑structure steps, and the ‘bands.X’ utility extracts the data.

The density of states (DOS) provides a statistical distribution of electronic states as a function of energy. The total DOS, N(E), is computed by summing contributions from all k‑points and applying a smearing (often a Gaussian) to produce a continuous curve. The partial DOS (PDOS) decomposes N(E) into atomic or orbital contributions, enabling the identification of which atoms or orbitals dominate near the Fermi level. In VASP, the ‘DOSCAR’ file stores DOS data, while the ‘PROCAR’ file contains projected information. Quantum Espresso’s ‘dos.X’ and ‘projwfc.X’ utilities serve analogous functions.

The Fermi level (E_F) is the chemical potential at zero temperature, separating occupied from unoccupied states. For metals, E_F lies within a band, leading to a finite DOS at E_F; for insulators and semiconductors, a gap separates the valence‑band maximum (VBM) from the conduction‑band minimum (CBM). The size of this band gap is a key property; however, standard DFT (LDA/GGA) notoriously underestimates gaps due to the self‑interaction error and the lack of derivative discontinuity in the XC functional. Hybrid functionals, DFT+U, and many‑body perturbation techniques (GW) are employed to improve gap predictions.

Spin polarization is essential when dealing with magnetic materials. In collinear spin‑polarized calculations, the spin density is split into up‑ and down‑spin components, and the Kohn‑Sham equations are solved for each spin channel independently. The resulting magnetic moment on an atom is obtained by integrating the spin density over a chosen region (e.G., Bader volume). Non‑collinear magnetism, where the spin direction varies in space, is supported in both codes via additional tags (e.G., ‘LNONCOLLINEAR = .TRUE.’ In VASP). Spin‑orbit coupling (SOC) can be turned on to capture relativistic effects, crucial for heavy elements and topological insulators. Including SOC typically doubles the size of the Hamiltonian, increasing computational cost.

DFT+U is a corrective scheme that adds an on‑site Hubbard term to the DFT Hamiltonian, addressing the self‑interaction error for localized d or f electrons. The method introduces two parameters: The Hubbard U (Coulomb) and the exchange J, often combined as U_eff = U − J. In VASP, the ‘LDAU = .TRUE.’ Flag activates the correction, and the ‘LDAUU’, ‘LDAUJ’, and ‘LDAUL’ tags specify the values for each atomic species. Quantum Espresso implements DFT+U via the ‘Hubbard_U’ and ‘Hubbard_J’ variables in the input. Proper selection of U values is critical; they can be derived from constrained DFT calculations, linear response, or fitted to experimental data.

Hybrid functionals such as HSE06 mix a fraction (typically 25 %) of exact exchange with a GGA component. The screened nature of HSE reduces the long‑range exchange, making it more tractable for periodic systems. In VASP, the ‘LHFCALC = .TRUE.’ Flag together with ‘HFSCREEN = 0.2’ Defines the HSE range‑separation parameter. Quantum Espresso’s hybrid implementation (e.G., ‘EXX’ keyword) requires the use of the ‘EXX’ module and is computationally intensive due to the need for evaluating exact‑exchange integrals.

The GW approximation represents a higher level of many‑body theory, where the self‑energy Σ is approximated as the product of the Green’s function G and the screened Coulomb interaction W. The method corrects quasiparticle energies, yielding more accurate band gaps and effective masses. GW calculations are typically performed as a post‑processing step on top of a DFT ground state. VASP offers a GW implementation (e.G., ‘ALGO = GW0’), while Quantum Espresso provides the ‘Yambo’ interface for GW and Bethe‑Salpeter calculations. GW is computationally demanding, scaling as O(N⁴), and requires careful convergence with respect to empty bands, dielectric cutoff, and k‑point sampling.

Bethe‑Salpeter equation (BSE) extends GW by explicitly treating electron‑hole interactions, enabling the calculation of optical absorption spectra, exciton binding energies, and dielectric functions. The BSE Hamiltonian is constructed from GW quasiparticle energies and screened Coulomb matrix elements, then diagonalized to obtain excitation energies. Practical BSE calculations often involve a limited number of valence and conduction bands, and the use of a model dielectric function to reduce cost.

The van der Waals (vdW) corrections address the missing dispersion forces in standard DFT. Empirical schemes such as DFT‑D2/D3 (Grimme) add pairwise C₆ R⁻⁶ terms with damping functions, while non‑local functionals (e.G., VdW‑DF, optB88‑vdW) incorporate dispersion directly into the XC functional. In VASP, the ‘IVDW = 12’ flag activates the DFT‑D3 correction, while Quantum Espresso supports vdW‑DF via the ‘vdW_corr’ keyword. Selecting an appropriate vdW method is important for layered materials, molecular adsorption, and weakly bonded interfaces.

Supercell constructions are employed to model defects, surfaces, and interfaces. By replicating the primitive cell along one or more lattice vectors, a larger periodic cell is created that isolates the defect or surface from its periodic images. The size of the supercell influences spurious interactions; convergence tests with respect to cell size are essential. For surface calculations, a slab model with a vacuum region (typically >15 Å) is used to separate the periodic images along the surface normal. Dipole corrections may be applied to counteract artificial electric fields arising from asymmetrical slabs; both Quantum Espresso and VASP provide dipole‑correction options (e.G., ‘LDIPOL = .TRUE.’ In VASP).

The crystal symmetry of a system is exploited to reduce computational effort. Symmetry operations map k‑points onto each other, allowing the code to consider only the irreducible Brillouin zone (IBZ). Users can either let the program automatically detect symmetry (default behavior) or manually specify a lower symmetry to explore, for example, magnetic ordering that breaks the original space group. In VASP, the ‘ISYM = 0’ tag disables symmetry, while Quantum Espresso’s ‘nosym’ option serves the same purpose.

Reciprocal lattice vectors (b₁,b₂,b₃) define the periodicity in momentum space, and their magnitudes are inversely proportional to the real‑space lattice constants. The relationship between real and reciprocal lattices is fundamental for constructing k‑point meshes and interpreting band structures. The volume of the reciprocal cell appears in the expression for the density of states and in the conversion between real‑space and reciprocal‑space integrals.

The charge density is often visualized to gain insight into bonding characteristics. Tools such as VESTA, XCrySDen, and the VASP ‘CHGCAR’ file enable the generation of isosurfaces and planar cuts. In Quantum Espresso, the ‘charge density’ is stored in the ‘charge-density.Dat’ file, which can be processed by the ‘pp.X’ post‑processor. Visual analysis assists in identifying charge accumulation at interfaces, polarization in ferroelectrics, and the presence of lone‑pair electrons.

Electron localization function (ELF) is another scalar field derived from the Kohn‑Sham orbitals, indicating the degree of electron pairing. ELF values range from 0 (delocalized) to 1 (highly localized), and are useful for distinguishing covalent bonds, metallic regions, and lone pairs. Both VASP and Quantum Espresso can calculate ELF via post‑processing utilities.

Bader charge analysis partitions the electron density into atomic basins defined by zero‑flux surfaces in the gradient of ρ(r). The resulting Bader charges provide a quantitative measure of charge transfer, often more reliable than Mulliken populations for plane‑wave calculations. The external program ‘bader’ processes the charge density files produced by VASP (‘CHGCAR’) or Quantum Espresso (‘charge-density.Dat’). Challenges include ensuring sufficient grid density and handling noisy charge densities near vacuum regions.

Mulliken population analysis is traditionally used with localized basis sets (e.G., Gaussian‑type orbitals) and is less common in plane‑wave codes. However, Quantum Espresso’s ‘Atomic Projected Density of States’ (projwfc.X) can approximate population analysis by projecting onto atomic orbitals defined within the pseudopotential.

Effective mass extraction is a routine post‑processing task for semiconductor applications. By fitting a parabolic dispersion near the band extrema (VBM or CBM), the curvature yields the effective mass tensor. Accurate effective masses require dense k‑point sampling around the extrema, and often a non‑self‑consistent band‑structure calculation on a fine mesh. The resulting values are essential for transport modeling and carrier mobility predictions.

Phonons describe lattice vibrations and are central to thermal properties, electron‑phonon coupling, and superconductivity. In Quantum Espresso, the ‘ph.X’ module computes phonon frequencies via density‑functional perturbation theory (DFPT). The process involves a self‑consistent electronic calculation at a given q‑point, followed by the linear response of the electron density to atomic displacements. VASP implements phonon calculations through external tools such as Phonopy, which require a series of finite‑difference force calculations. Convergence with respect to supercell size and k‑point sampling is critical for accurate phonon spectra.

Electron‑phonon coupling (EPC) quantifies the interaction strength between electronic states and phonons, influencing phenomena such as superconductivity and carrier scattering. EPC constants (λ) can be extracted from DFPT calculations, and the Eliashberg spectral function α²F(ω) provides a frequency‑resolved picture of coupling. In Quantum Espresso, the ‘lambda.X’ utility processes DFPT data to compute λ and related quantities. VASP users typically rely on the EPW (Electron‑Phonon Wannier) package interfaced with Wannier90 to perform EPC calculations.

Wannier functions are maximally localized orbitals constructed from Bloch states, facilitating the interpolation of band structures, the calculation of transport properties, and the analysis of topological invariants. The Wannier90 code interfaces with both Quantum Espresso and VASP, requiring the generation of a set of seed projections and the extraction of overlap matrices (the ‘amn’ and ‘mmn’ files). Once obtained, Wannier90 can produce interpolated band structures on arbitrarily dense k‑meshes, enabling fine‑resolution Fermi‑surface studies.

Berry phase methods are employed to compute electronic polarization, orbital magnetization, and topological invariants such as Chern numbers. In VASP, the ‘LCALCPOL = .TRUE.’ Flag activates the modern theory of polarization, while Quantum Espresso’s ‘berry_phase.X’ utility performs similar calculations. The Berry phase is evaluated as a product of overlaps between wavefunctions at adjacent k‑points, and requires a well‑converged k‑mesh along the polarization direction.

Magnetic anisotropy energy (MAE) quantifies the energy difference between different spin‑orientation directions, arising from SOC and crystal field effects. Calculating MAE involves performing two SCF runs with the magnetization constrained along orthogonal axes and taking the energy difference. Since MAE values are often on the order of μeV per atom, highly converged settings (e.G., Dense k‑meshes, tight energy thresholds) are mandatory. Both VASP and Quantum Espresso support non‑collinear SOC calculations, but VASP’s implementation is generally more streamlined for MAE studies.

Defect formation energy is a key metric for assessing the thermodynamic stability of vacancies, interstitials, and substitutional impurities. The formation energy ΔE_f is computed as ΔE_f = E_defect − E_bulk + ∑_i n_i μ_i + q (E_F + E_VBM) + ΔV, where E_defect and E_bulk are the total energies of the defective and pristine supercells, n_i denotes the number of atoms added or removed, μ_i are the chemical potentials, q is the defect charge state, and ΔV accounts for potential alignment. Accurate defect calculations require large supercells, charge‑correction schemes (e.G., Makov‑Payne), and careful convergence of the electronic structure.

Charge‑state corrections address the spurious electrostatic interactions between periodic images of charged defects. Methods such as the Freysoldt–Neugebauer–Van de Walle (FNV) scheme involve aligning the electrostatic potential far from the defect and applying a correction term based on the dielectric constant. Both VASP and Quantum Espresso users can apply these corrections using external scripts (e.G., ‘Sxdefectalign’).

Elastic constants are derived from the stress‑strain relationship. In practice, a series of small deformations (e.G., ±1 %) Are applied to the cell, the resulting stresses are computed, and the elastic tensor C_ij is obtained by fitting. VASP’s ‘IBRION = 6’ flag automates this process, while Quantum Espresso’s ‘elastic.X’ utility performs similar calculations. The derived bulk modulus, shear modulus, and Poisson’s ratio provide insight into mechanical stability and hardness.

Thermodynamic integration and ab‑initio molecular dynamics (AIMD) extend static DFT to finite temperatures. In AIMD, the forces from DFT are used to propagate atoms via Newton’s equations of motion, generating trajectories that sample the canonical ensemble. VASP implements AIMD using the Nosé‑Hoover thermostat (e.G., ‘SMASS = -3’), while Quantum Espresso’s ‘cp.X’ module provides Car‑Parrinello and Born‑Oppenheimer dynamics. Convergence of AIMD requires adequate time steps (≈1 fs for light atoms) and sufficient simulation length to achieve statistical averaging.

Car‑Parrinello molecular dynamics (CPMD) treats the electronic wavefunctions as dynamical variables coupled to the ionic motion, enabling larger time steps but requiring careful control of the fictitious electron mass. The method is particularly useful for exploring reaction pathways and structural phase transitions. In Quantum Espresso, CPMD is activated via the ‘cp.X’ executable with appropriate input parameters (e.G., ‘Electron_mass’ and ‘time_step’). VASP does not natively support CPMD; users typically resort to standard Born‑Oppenheimer AIMD.

Transition state search methods such as the nudged elastic band (NEB) technique locate minimum‑energy paths between initial and final states, identifying activation barriers. VASP provides a built‑in NEB implementation (e.G., ‘IBRION = 3’), while Quantum Espresso users can employ external tools like ‘neb.X’ or interface with the ‘Atomic Simulation Environment’ (ASE). Convergence of NEB calculations hinges on the number of images, spring constants, and force criteria; inadequate settings may lead to spurious kinks or incomplete convergence.

Hybrid functional calculations often demand a reduced k‑point mesh because the exact‑exchange term scales poorly with the number of k‑points. A common strategy is to perform a coarse‑grid hybrid run to obtain band gaps, then refine the result with a GGA calculation for structural properties. Users must remain vigilant about the possible inconsistency between geometries optimized with different functionals, especially when lattice constants differ by several percent.

Convergence testing is an indispensable practice. Typical parameters to test include plane‑wave cutoff, k‑point density, smearing width, and the number of empty bands (relevant for GW and hybrid calculations). A systematic approach involves varying one parameter while keeping others fixed, plotting the target observable (energy, forces, or band gap) against the parameter, and identifying a plateau where changes fall below a predefined tolerance. Documenting these tests ensures reproducibility and builds confidence in the simulation protocol.

Input file structure differs between the two codes but follows a similar logical flow: Definition of the crystal lattice, specification of atomic positions, selection of pseudopotentials, and setting of calculation parameters. In VASP, the POSCAR file contains lattice vectors and atomic coordinates, POTCAR bundles the pseudopotentials, INCAR holds control flags (e.G., ‘ENCUT’, ‘EDIFF’, ‘ISMEAR’), and KPOINTS defines the k‑mesh. Quantum Espresso uses a single input file where sections such as ‘&control’, ‘&system’, ‘&electrons’, and ‘ATOMIC_SPECIES’ are delimited by ‘/’. Understanding the mapping between these sections and their VASP counterparts aids in transferring workflows between the two environments.

Parallelization strategies are crucial for scaling calculations to high‑performance computing resources. Both codes support MPI parallelism, with domain decomposition over k‑points or plane‑wave coefficients. VASP also offers OpenMP threading for shared‑memory parallelism, while Quantum Espresso can distribute plane‑wave FFTs across processor groups. Efficient parallel execution requires careful matching of the number of MPI ranks to the problem size; oversubscription may degrade performance due to communication overhead.

File formats for wavefunctions and charge densities are binary and code‑specific. VASP’s ‘WAVECAR’ and ‘CHGCAR’ store wavefunction coefficients and electron density, respectively; these files can be read by the VASP post‑processing suite. Quantum Espresso’s ‘wavefunction.Dat’ and ‘charge-density.Dat’ serve analogous purposes. Interoperability is facilitated through conversion tools (e.G., ‘Pw2bgw’, ‘pw2vasp’, or the ‘cif2cell’ utility) that translate structures between formats such as CIF, POSCAR, and Quantum Espresso’s input syntax.

Common pitfalls include neglecting to relax both atomic positions and cell shape when dealing with pressure‑dependent phenomena, using insufficient vacuum for slab models (leading to spurious interactions), and overlooking the need for dipole corrections in asymmetric slabs. Another frequent source of error is the mismatch between the pseudopotential type and the chosen XC functional; for example, using a GGA‑type PAW dataset with an LDA functional can introduce inconsistencies. Users should always verify that pseudopotentials are generated with the same functional as the calculation.

Best‑practice recommendations encompass: (1) Performing thorough convergence tests for each new material class; (2) employing symmetry analysis to reduce computational load while ensuring that magnetic ordering is correctly captured; (3) validating pseudopotential choices against all‑electron benchmark data when available; (4) documenting all input parameters, version numbers, and computational settings to guarantee reproducibility; (5) leveraging automated workflow tools (e.G., AiiDA, FireWorks, or Custodian) to manage job submission, error handling, and data provenance; and (6) staying current with code updates, as both Quantum Espresso and VASP frequently introduce performance improvements and new methodological features.

Practical example: Band‑gap calculation for silicon. A typical workflow begins with a self‑consistent SCF run using a GGA‑PBE functional, a plane‑wave cutoff of 500 eV (VASP) or 60 Ry (Quantum Espresso), and a 8×8×8 Monkhorst‑Pack k‑mesh. After convergence, a non‑self‑consistent calculation is performed along the high‑symmetry path Γ‑X‑W‑K‑Γ‑L, extracting eigenvalues for plotting. The resulting PBE gap (≈0.6 EV) underestimates the experimental value (≈1.1 EV). To improve accuracy, a hybrid HSE06 calculation is set up with a reduced k‑mesh (e.G., 4×4×4) And the same cutoff. The HSE band structure yields a gap of ≈1.2 EV, in good agreement with experiment. This example illustrates the trade‑off between computational cost and accuracy, and highlights the importance of selecting the appropriate functional for the property of interest.

Practical example: Defect calculation in TiO₂. A 2×2×3 supercell of rutile TiO₂ is constructed, and a neutral oxygen vacancy is introduced by removing an O atom. The defect supercell is relaxed with a 400 eV cutoff and a 2×2×2 k‑mesh. The total energy of the defect cell is compared to that of the pristine cell, and the formation energy is computed using the chemical potential of oxygen derived from O₂ gas (adjusted for temperature and pressure). Charge corrections are applied using the FNV scheme, and the resulting formation energy informs the likelihood of vacancy formation under various synthesis conditions. This workflow demonstrates the integration of supercell modeling, chemical potential considerations, and post‑processing corrections.

Practical example: Phonon dispersion of graphene. A 4×4 supercell of graphene is generated, and DFPT calculations are performed at each q‑point along the high‑symmetry directions Γ‑M‑K‑Γ. The resulting dynamical matrices are diagonalized to obtain phonon frequencies, revealing the characteristic linear acoustic branches near Γ and the out‑of‑plane flexural mode. Convergence with respect to supercell size and k‑point sampling is examined, ensuring that the acoustic sum rule is satisfied. The phonon spectrum can be further used to compute the thermal conductivity via the Boltzmann transport equation, illustrating the link between vibrational properties and macroscopic transport.

Practical example: Surface adsorption of CO on Pt(111). A slab model consisting of five Pt layers with a 15 Å vacuum gap is built; the bottom two layers are fixed while the top three are allowed to relax. A CO molecule is placed at a bridge site, and a geometry optimization is performed with a 400 eV cutoff and a 5×5×1 k‑mesh. After relaxation, the adsorption energy is calculated as E_ads = E_CO/Pt − (E_Pt + E_CO), where E_CO/Pt is the total energy of the slab with CO, E_Pt is the energy of the clean slab, and E_CO is the energy of an isolated CO molecule in a large cell. Van der Waals corrections (e.G., DFT‑D3) are included to capture dispersion interactions. The resulting adsorption energy (≈−1.5 EV) can be compared with experimental temperature‑programmed desorption data, providing validation of the computational approach.

Practical example: Magnetic anisotropy of FePt. A tetragonal FePt alloy is modeled, and non‑collinear SOC calculations are performed with the magnetization constrained along the [001] and [100] directions. The energy difference, ΔE = E_[100] − E_[001], yields the MAE, which is typically on the order of a few meV per formula unit. Because MAE is sensitive to k‑point sampling, a dense 12×12×12 mesh and a tight energy convergence (10⁻⁸ eV) are required. The computed MAE can be correlated with experimental coercivity measurements, providing insight into the material’s suitability for magnetic recording applications.

Practical example: GW correction of band structure for MoS₂. A monolayer of MoS₂ is first relaxed using GGA‑PBE with a 60 Ry cutoff and a 12×12×1 k‑mesh. A non‑self‑consistent calculation with 200 empty bands is then performed to obtain the dielectric matrix. The GW0 calculation is executed, updating the quasiparticle energies iteratively while keeping the screened interaction fixed. Convergence with respect to the number of empty bands, the dielectric cutoff, and the k‑mesh is carefully checked. The resulting quasiparticle gap (~2.8 EV) aligns closely with experimental optical measurements after accounting for excitonic binding via BSE. This case highlights the multi‑step nature of many‑body perturbation calculations and the necessity of systematic convergence studies.

Practical example: BSE optical spectrum of ZnO. After a GW calculation provides corrected quasiparticle energies, the BSE Hamiltonian is constructed using a limited set of valence and conduction bands (e.G., 4 Valence and 4 conduction). The excitonic eigenvalues are solved, yielding absorption peaks corresponding to bound excitons. The calculated exciton binding energy (~60 meV) can be compared with photoluminescence data. This workflow demonstrates how higher‑order electron–hole interactions refine the description of optical properties beyond the independent‑particle approximation.

Practical example: Charge analysis in LiFePO₄. A Bader analysis is performed on the charge density obtained from a GGA‑U (U_eff = 4.3 EV on Fe) calculation. The resulting Bader charges reveal a charge transfer from Li to Fe–O framework, consistent with the oxidation state Fe³⁺ in the fully lithiated phase. Comparing the Bader charges before and after lithium extraction provides quantitative insight into the redox mechanism, essential for battery material design.

Practical example: Elastic constants of cubic AlN. Small strains (±0.5 %) Are applied to the cubic cell, and the resulting stress tensors are extracted from VASP SCF runs. By fitting the stress–strain relationship, the three independent elastic constants C₁₁, C₁₂, and C₄₄ are obtained. These values are then used to compute the bulk modulus (B = (C₁₁ + 2C₁₂)/3) and shear modulus (G = (C₁₁ − C₁₂ + 3C₄₄)/5). The mechanical stability criteria (C₁₁ > C₁₂, C₄₄ > 0) are verified, confirming that the structure is mechanically stable.

Practical example: AIMD of liquid water. A 64‑molecule cubic cell is simulated at 300 K using VASP with the PBE functional, a 400 eV cutoff, and a 1 fs time step. The Nosé‑Hoover thermostat maintains the temperature, and a trajectory of 20 ps is generated. Structural properties such as the radial distribution function g_OO(r) are computed from the trajectory and compared with experimental neutron scattering data. The diffusion coefficient is extracted from the mean‑square displacement, providing a benchmark for the chosen functional’s ability to describe hydrogen bonding.

Practical example: NEB reaction pathway for H₂ dissociation on Cu(111).

Key takeaways

  • The following exposition presents the most frequently encountered terms, organized thematically, and illustrated with examples and typical challenges faced by users of these two leading first‑principles codes.
  • The total energy is a functional of this density, and the self‑consistent solution of the Kohn‑Sham equations provides both the ground‑state density and the associated eigenvalues ε_i.
  • Common approximations include the local density approximation (LDA), which assumes the XC energy at each point depends only on the local density, and the generalized gradient approximation (GGA), which incorporates the density gradient.
  • In plane‑wave based codes like Quantum Espresso and VASP, the electronic wavefunctions are expanded in a set of plane waves, ψ_i(r)=∑_G c_i(G) e^{iG·r}, where G are reciprocal‑lattice vectors.
  • A pseudopotential replaces the all‑electron potential with a smoother effective potential that reproduces the scattering properties of the valence electrons while removing the need to explicitly treat core states.
  • The k‑point density directly influences the accuracy of total‑energy and forces; metallic systems generally require denser sampling because of partially occupied states near the Fermi level.
  • Techniques such as mixing schemes (Pulay, Broyden) and smearing methods (Methfessel‑Paxton, Gaussian, Fermi‑Dirac) are employed to stabilize the iteration.
May 2026 intake · open enrolment
from £90 GBP
Enrol