Materials Science Applications

Quantum Espresso and VASP are two of the most widely used first‑principles codes for exploring the electronic structure and related properties of materials. Mastery of the terminology that underpins these packages is essential for anyone aiming to apply them to real‑world problems in materials science. The following exposition presents the core vocabulary, enriched with examples, practical uses, and common challenges that arise during simulations. The discussion is organized thematically, moving from the description of crystal geometry to the electronic‑structure formalism, then to the treatment of forces and dynamics, and finally to advanced topics such as defects, surfaces, and many‑body corrections.

The description of a crystal begins with the concept of a lattice. A lattice is an infinite array of points generated by translating a set of primitive vectors. The smallest repeating unit that captures the symmetry of the lattice is the unit cell. Unit cells are classified as primitive or conventional; the former contains the minimum number of atoms required to reproduce the crystal, while the latter often aligns with higher‑symmetry axes for easier visualisation. For example, the face‑centered cubic (fcc) conventional cell contains four atoms, whereas its primitive cell contains only one. Understanding the distinction is crucial when constructing input files, because the number of atoms entered into the code must correspond to the chosen cell definition.

Associated with the lattice is the concept of reciprocal space. Reciprocal vectors are defined as the Fourier‑conjugate of the real‑space lattice vectors and form the basis for describing electron wavefunctions in periodic systems. The Brillouin zone, the Wigner‑Seitz cell of the reciprocal lattice, contains all distinct k‑points that need to be sampled. In practice, a finite set of k‑points is selected using schemes such as the Monkhorst‑Pack grid or the Gamma‑centered mesh. The density of this grid directly influences the accuracy of the calculated total energy, forces, and derived properties. As a rule of thumb, metals typically require denser k‑point meshes than insulators because of the abrupt change in occupancy at the Fermi surface.

The electronic‑structure problem is solved within the framework of density functional theory (DFT). DFT rests on the Hohenberg‑Kohn theorems, which state that the ground‑state energy of a many‑electron system is a unique functional of the electron density. The practical implementation replaces the many‑body problem with a set of single‑particle equations, the Kohn‑Sham equations, that are solved self‑consistently. The key ingredient of DFT is the exchange‑correlation functional, which captures the many‑body effects beyond the classical Coulomb interaction. The most common approximations are the local density approximation (LDA) and the generalized gradient approximation (GGA). LDA assumes that the exchange‑correlation energy at each point depends only on the local density, while GGA includes the gradient of the density, providing better accuracy for systems with inhomogeneous electron distributions.

A central practical choice is the representation of the electron–ion interaction. Plane‑wave basis sets are employed in both Quantum Espresso and VASP because they naturally satisfy periodic boundary conditions and are systematically improvable by increasing the cutoff energy. The cutoff determines the maximum kinetic energy of plane waves retained in the expansion; higher cutoffs yield more accurate results but increase computational cost. To reduce the required cutoff, pseudopotentials are introduced, which replace the all‑electron potential with a smoother effective potential that reproduces the valence‑electron behavior. There are several families of pseudopotentials: norm‑conserving, ultrasoft, and the projector‑augmented wave (PAW) method. PAW, used exclusively in VASP, combines the accuracy of all‑electron calculations with the efficiency of pseudopotentials by reconstructing the full wavefunction from a smooth part and atom‑centered augmentation functions.

When setting up a calculation, the user must decide on the occupancy scheme that dictates how electronic states are filled. For insulators and semiconductors, a simple fixed‑occupation (or “cold‑smearing”) approach suffices, as the band gap ensures a clear separation between occupied and empty states. Metals, however, require a smearing technique to aid convergence. Common smearing methods include the Methfessel‑Paxton scheme, the Gaussian smearing, and the Fermi‑Dirac distribution. The smearing width must be chosen carefully: Too large a value artificially broadens the electronic states and can affect calculated properties, while too small a value may hinder SCF convergence.

The iterative process that solves the Kohn‑Sham equations is called the self‑consistent field (SCF) cycle. At each iteration, the code updates the electron density, constructs a new Kohn‑Sham Hamiltonian, and solves for the eigenvalues and eigenvectors (the Kohn‑Sham orbitals). Convergence is monitored by the change in total energy, charge density, or forces between successive iterations. When the SCF is converged, one may perform a non‑self‑consistent field (NSCF) calculation to generate a finer sampling of the band structure or density of states without recomputing the charge density. This two‑step approach saves time because the expensive SCF step is performed only once.

The output of a converged SCF calculation provides the band structure, which plots the eigenvalues as a function of k‑point along high‑symmetry directions in the Brillouin zone. The band structure reveals whether the material is metallic, semiconducting, or insulating by indicating the presence or absence of a band gap at the Fermi level. Complementary to the band structure is the density of states (DOS), which quantifies the number of electronic states per energy interval. The DOS can be projected onto atomic orbitals (projected DOS or PDOS) to identify the contributions of specific elements or orbitals to particular energy ranges, a useful tool for interpreting chemical bonding and magnetic behavior.

Magnetism is another essential aspect of many materials. In DFT calculations, spin polarization is introduced by treating the spin‑up and spin‑down electrons separately, resulting in two independent Kohn‑Sham equations. The resulting magnetic moment on each atom is obtained from the difference in spin densities. For transition‑metal oxides and rare‑earth compounds, the standard LDA or GGA often fails to capture strong electron correlation, leading to an underestimation of band gaps and magnetic moments. To address this, the DFT+U method adds a Hubbard‑type correction term that penalises fractional occupancies of localized d or f orbitals. The magnitude of the U parameter can be determined from linear‑response calculations or fitted to experimental data.

Beyond static electronic properties, both Quantum Espresso and VASP can compute forces and stresses, enabling structural relaxation and molecular dynamics. The force on each atom is derived from the Hellmann‑Feynman theorem, which states that the force equals the negative gradient of the total energy with respect to atomic positions, provided the wavefunctions are fully converged. The stress tensor describes the response of the system to changes in the simulation cell shape and volume. Geometry optimization proceeds by iteratively moving atoms along the direction of the forces while simultaneously adjusting the cell parameters to reduce the stress. Convergence criteria for forces (e.G., 0.01 EV/Å) and stresses (e.G., 0.5 Kbar) must be specified, and a careful balance between accuracy and computational cost is required.

A common practical challenge is the selection of an appropriate convergence threshold. Tight thresholds guarantee reliable results but increase the number of SCF iterations and the overall runtime. Conversely, loose thresholds may lead to spurious minima or inaccurate forces, especially when calculating vibrational properties. A systematic convergence test—varying the cutoff energy, k‑point density, and smearing width—should be performed for each new material system. The test typically involves plotting the total energy versus each parameter and identifying a plateau where further increases produce changes smaller than the desired accuracy (often a few meV per atom).

The next level of analysis concerns lattice dynamics. The harmonic approximation treats atomic vibrations as independent normal modes, each characterized by a frequency and polarization vector. Within this approximation, the dynamical matrix is constructed from the second derivatives of the total energy with respect to atomic displacements. In Quantum Espresso, the phonon module computes these derivatives using density‑functional perturbation theory (DFPT), which avoids the need for explicit finite‑difference supercells. VASP, on the other hand, typically employs finite‑difference approaches, where a supercell is generated and atoms are displaced by a small amount (often 0.01 Å) to evaluate the forces and build the force‑constants matrix. The resulting phonon dispersion curves reveal dynamical stability (absence of imaginary frequencies) and provide thermodynamic quantities such as free energy, entropy, and specific heat.

Phonon calculations also enable the study of electron‑phonon coupling, a key factor in conventional superconductivity. The electron‑phonon coupling constant λ can be extracted from the Eliashberg spectral function α²F(ω), which quantifies the interaction strength between electrons at the Fermi level and phonons of frequency ω. Accurate evaluation of λ requires dense k‑ and q‑point meshes, as well as fine smearing parameters, making these calculations computationally demanding.

Materials with defects, impurities, or surfaces demand a different modeling strategy. Defects are introduced by creating a supercell, a larger periodic cell that contains the defect while preserving the underlying lattice periodicity. For a vacancy, one simply removes an atom from the supercell; for an interstitial, an extra atom is added at a high‑symmetry site. The defect formation energy is then computed as the difference between the total energy of the defective supercell and that of the pristine supercell, corrected for the chemical potentials of added or removed species. Charge states are accommodated by adding or removing electrons and compensating with a uniform background charge, followed by corrections for spurious electrostatic interactions (e.G., The Makov‑Payne correction). These steps are essential for predicting defect concentrations and their impact on electronic properties.

Surface calculations employ a slab geometry, where a finite number of atomic layers are stacked along a chosen crystallographic direction, and a vacuum region separates periodic images. The slab thickness must be sufficient to converge surface energies and electronic states, while the vacuum thickness must eliminate interactions between opposite surfaces. Surface relaxations often differ significantly from bulk relaxations, leading to reconstruction patterns that can be captured by allowing all atoms in the slab to relax. Adsorption studies then place an adsorbate molecule on the surface and compute the adsorption energy as the difference between the total energy of the adsorbate‑surface system and the sum of the isolated surface and molecule energies.

A further extension involves molecular dynamics (MD) simulations, which trace the evolution of atomic positions over time according to Newton’s equations of motion. In the ab initio context, forces are obtained from DFT at each MD step, yielding so‑called Born‑Oppenheimer MD. Quantum Espresso provides a dedicated CP (Car‑Parrinello) module that propagates both electronic and ionic degrees of freedom simultaneously, allowing for larger time steps but requiring careful tuning of the fictitious electron mass and the thermostat. VASP implements Born‑Oppenheimer MD with efficient parallelisation, making it suitable for studying processes such as diffusion, phase transitions, and reaction pathways. Temperature control is achieved using thermostats (e.G., Nosé‑Hoover) or by scaling velocities, while pressure control utilizes barostats (e.G., Parrinello‑Rahman). The choice of ensemble (NVE, NVT, NPT) influences the thermodynamic quantities that can be extracted from the trajectory.

When dealing with strongly correlated materials, standard DFT functionals often prove insufficient. Hybrid functionals, such as HSE06 or PBE0, mix a fraction of exact Hartree‑Fock exchange with a GGA exchange‑correlation term, improving band‑gap predictions and magnetic properties. Both Quantum Espresso and VASP support hybrid functional calculations, though they are considerably more expensive due to the non‑local nature of the exchange operator. For even higher accuracy, many‑body perturbation methods like the GW approximation can be employed. In GW, the electron self‑energy Σ is approximated by the product of the Green’s function G and the screened Coulomb interaction W, leading to quasiparticle energies that correct the DFT band structure. VASP implements GW as a post‑processing step on top of a DFT calculation, while Quantum Espresso offers GW via the Yambo or BerkeleyGW packages. These methods require careful convergence with respect to the number of empty bands, the dielectric cutoff, and the k‑point sampling.

Time‑dependent DFT (TDDFT) extends the static DFT formalism to excited‑state phenomena. It is widely used for calculating optical absorption spectra, exciton binding energies, and response functions. In practice, TDDFT calculations are performed by linear‑response methods that perturb the ground‑state electron density and solve the resulting eigenvalue problem for excitation energies. Both Quantum Espresso (via the turboTDDFT module) and VASP (through the LOPTICS tag) can compute dielectric functions, enabling the analysis of materials for photovoltaic or photocatalytic applications.

The practical workflow for a typical materials‑science project using Quantum Espresso or VASP can be summarised in several stages. First, the crystal structure is obtained from experimental databases (e.G., ICSD) or generated using crystallographic software. Next, a suitable pseudopotential library is selected, ensuring that the chosen potentials include the required valence states and are compatible with the chosen functional (e.G., GGA‑PBE). The user then defines the plane‑wave cutoff, k‑point mesh, and smearing parameters, performing convergence tests to balance accuracy and cost. An SCF calculation yields the ground‑state electron density, from which the total energy, band structure, and DOS are extracted. Geometry optimisations are then carried out to relax atomic positions and cell parameters, using force and stress convergence criteria appropriate for the target property. Once the relaxed structure is obtained, further analyses such as phonon calculations, defect formation energies, surface studies, or MD simulations can be performed. For properties that demand higher‑level theory, hybrid functionals or GW post‑processing are applied, often on top of the converged DFT ground state.

During the course of these calculations, a number of challenges frequently arise. One common issue is SCF convergence failure, which may manifest as oscillating total energies or persistent residual forces. Strategies to resolve this include increasing the mixing parameter, using a different mixing scheme (e.G., Broyden), adjusting the smearing width, or providing a better initial guess for the charge density (e.G., From a previous calculation with a similar structure). Another difficulty is the appearance of spurious imaginary phonon modes in a nominally stable structure. These can result from insufficient supercell size, inadequate k‑point sampling, or too small a displacement amplitude in finite‑difference calculations. Careful convergence testing and symmetry analysis often eliminate the false instabilities. In defect calculations, the finite‑size error due to the artificial periodicity of the supercell can be significant, especially for charged defects. Corrections such as the Freysoldt, Makov‑Payne, or Kumagai‑Oba schemes must be applied to obtain reliable formation energies.

Parallelisation is a key consideration for large‑scale simulations. Both Quantum Espresso and VASP employ MPI (Message Passing Interface) to distribute k‑points, plane‑wave components, and band indices across multiple processors. Efficient parallel performance requires a balanced distribution of work; for example, assigning too many processors to a small k‑point set can lead to idle cores. In VASP, the KPAR and NPAR tags allow the user to control the division of k‑point and band parallelisation, respectively. In Quantum Espresso, the npools option splits the k‑point pool, while the nthreads environment variable enables OpenMP threading for intra‑pool parallelism. Optimising these parameters is essential for achieving good scaling on modern multi‑core clusters.

Another advanced topic is the treatment of van der Waals (vdW) interactions, which are poorly described by conventional LDA or GGA functionals. Several approaches exist to incorporate dispersion forces: The semi‑empirical DFT‑D3 correction adds pairwise C₆ terms with a damping function; the non‑local vdW‑DF functional introduces a fully non‑local correlation term; and the many‑body dispersion (MBD) method captures collective electronic fluctuations. VASP supports DFT‑D3 and vdW‑DF, while Quantum Espresso offers DFT‑D2/D3, vdW‑DF, and the Tkatchenko‑Scheffler scheme. Selecting the appropriate vdW method is critical for layered materials (e.G., Graphene, transition‑metal dichalcogenides) and for adsorption studies where weak interactions dominate.

In the realm of high‑throughput materials discovery, the vocabulary expands to include concepts such as workflow automation, database integration, and machine‑learning descriptors. Automated pipelines generate input files, launch calculations, monitor convergence, and parse outputs to store results in repositories like the Materials Project or the Open Quantum Materials Database (OQMD). Key descriptors extracted from DFT outputs—formation energy per atom, band gap, magnetic moment, elastic constants—serve as features for training predictive models. Understanding the terminology of these pipelines (e.G., “Job scheduler”, “error handler”, “post‑processing script”) enables researchers to scale their investigations from a single material to thousands.

Elastic properties are derived from the second derivative of the total energy with respect to strain. By applying a series of small deformations to the relaxed structure and computing the resulting stress tensors, one can extract the full elastic constant tensor Cᵢⱼ. The Voigt, Reuss, and Hill averages provide bulk and shear moduli, while the Poisson ratio and Young’s modulus describe the material’s response to mechanical loading. Both Quantum Espresso and VASP can automate this process through built‑in scripts (e.G., The elastic module in Quantum Espresso or the IBRION=6 tag in VASP). Accurate elastic constants require well‑converged k‑point meshes and high cutoff energies because stresses are sensitive to numerical noise.

Thermodynamic stability analysis often employs the concept of the convex hull in composition space. By calculating the formation energies of a series of compounds within a given chemical system, one constructs a hull that defines the lowest‑energy phases at each composition. Materials that lie on the hull are thermodynamically stable, while those above the hull are metastable and may decompose into a mixture of stable phases. Tools such as the pymatgen library can generate convex hulls from DFT data, facilitating the screening of novel compositions. Understanding terms such as “above‑hull energy”, “decomposition pathway”, and “phase diagram” is essential for interpreting these results.

Another practical aspect concerns the treatment of spin‑orbit coupling (SOC). SOC becomes significant for heavy elements where relativistic effects split degenerate bands, influencing band topology and magnetic anisotropy. In VASP, SOC is enabled with the LSORBIT tag; in Quantum Espresso, the noncolinear and lspinorb flags activate the relativistic Hamiltonian. Inclusion of SOC often modifies the band gap and may induce topological insulating behavior, as seen in bismuth selenide (Bi₂Se₃) or tungsten ditelluride (WTe₂). Accurate SOC calculations demand dense k‑point sampling and careful convergence testing because the spin‑splitting energies can be small.

The analysis of electronic topology introduces terminology such as Z₂ invariant, Chern number, and Weyl points. These quantities are derived from the Berry curvature of the occupied Bloch states and can be computed using post‑processing tools that read the wavefunction files generated by Quantum Espresso or VASP. For instance, the Z2Pack code can evaluate Z₂ invariants from VASP wavefunctions, while the Wannier90 package facilitates the construction of maximally localized Wannier functions, which serve as an efficient basis for evaluating topological indices and for interpolating band structures. Mastery of this vocabulary enables researchers to identify and design topological materials with robust surface states.

In experimental validation, the calculated quantities must be compared against measurable observables. For example, the band gap obtained from DFT is often compared with optical absorption edges measured by UV‑vis spectroscopy; however, standard GGA functionals typically underestimate the gap by 30–50 %. Hybrid functionals or GW corrections bring the theoretical values closer to experiment. The elastic constants can be compared with ultrasonic measurements or nanoindentation data. Phonon frequencies are validated against Raman or infrared spectroscopy; imaginary frequencies signal dynamical instability, prompting a re‑examination of the crystal structure or the computational parameters. Understanding the correspondence between computed and experimental descriptors is vital for assessing the reliability of the simulation.

When dealing with alloys or disordered systems, the special quasi‑random structure (SQS) approach is frequently employed. An SQS mimics the statistical distribution of atomic species in a random alloy within a finite supercell, allowing DFT to approximate the properties of the disordered phase. The construction of an SQS involves optimisation of pair, triplet, and higher‑order correlation functions, often using Monte‑Carlo algorithms. Once an SQS is generated, the same workflow of SCF, relaxation, and property calculation is applied. This technique is essential for studying composition‑dependent properties such as alloy phase stability, electronic conductivity, and lattice parameters.

The calculation of optical properties extends beyond the static dielectric constant. By evaluating the frequency‑dependent dielectric function ε(ω), one can derive the absorption coefficient, reflectivity, and loss function. These quantities are obtained from the linear‑response of the electron density to an external electric field, which in VASP is activated with the LOPTICS tag, while in Quantum Espresso the epsilon.X utility performs the same task. For materials with strong excitonic effects, the Bethe‑Salpeter equation (BSE) provides a more accurate description, albeit at a higher computational cost. BSE calculations are typically performed as a post‑processing step on top of GW quasiparticle energies.

The field of catalysis heavily relies on surface DFT calculations. Key descriptors such as the adsorption energy, reaction barrier, and activation energy are extracted from potential energy surfaces generated by nudged elastic band (NEB) calculations. The NEB method interpolates a series of images between initial and final states and optimises them to locate the minimum energy path. Both Quantum Espresso and VASP implement NEB, with VASP providing the IBRION=3 option and Quantum Espresso offering the neb.X driver. Accurate barrier heights require careful convergence of the forces and the inclusion of zero‑point energy corrections, especially for reactions involving light atoms such as hydrogen.

Thermal transport properties are accessed through the solution of the Boltzmann transport equation (BTE). The lattice contribution to thermal conductivity κₗ is computed from phonon lifetimes, group velocities, and specific heats, which are derived from the anharmonic force constants. The third‑order force constants are obtained by finite‑difference displacements of supercells, a procedure that can be automated with the thirdorder.Py script in the ShengBTE package. Electronic contributions to thermal conductivity, κₑ, are evaluated from the electronic band structure and scattering rates, often using the BoltzTraP code. Understanding terms such as “relaxation time”, “mean free path”, and “phonon‑phonon scattering” is essential for interpreting the computed transport coefficients.

In the context of battery materials, DFT is employed to calculate the intercalation voltage. The voltage is derived from the difference in total energies between the host material with and without the inserted ion (e.G., Li⁺), divided by the charge transferred. Accurate voltage predictions require consistent treatment of the reference state (often metallic Li) and careful convergence of the electronic and ionic degrees of freedom. Moreover, the migration barriers for ion diffusion are obtained from NEB calculations, providing insight into the rate capability of the electrode material.

Finally, the emerging area of machine‑learning potentials (MLPs) seeks to bridge the gap between quantum accuracy and classical efficiency. By training neural networks or Gaussian‑process models on DFT‑generated datasets of energies and forces, one can obtain interatomic potentials that reproduce DFT‑level accuracy at a fraction of the computational cost. Popular frameworks such as DeepMD, ANI, and Gaussian Approximation Potential (GAP) rely on descriptors like the smooth overlap of atomic positions (SOAP) or symmetry functions. The workflow typically involves generating a diverse set of DFT configurations (e.G., Via MD or random structure searches), extracting energies and forces, training the model, and validating it against a test set. Once validated, the MLP can be used for large‑scale simulations, enabling the study of phenomena such as crack propagation, phase transformations, and long‑time diffusion that are otherwise inaccessible to direct DFT.

Throughout these topics, the recurring theme is the careful balance between computational accuracy and resource expenditure. The practitioner must be fluent in the terminology that describes each approximation, method, and parameter, as well as the practical steps required to implement them in Quantum Espresso and VASP. Mastery of this vocabulary empowers researchers to design robust simulations, interpret results with confidence, and push the frontiers of materials discovery.