01 / FOUNDATIONS
Quantum Chemistry & DFT
From the wavefunction of an electron to the calculable and testable energy and properties of crystals. Understand the model first, then decide the parameters.
The input examples in this manual have been checked against the data and have not yet been verified by VASP actual calculations case by case. When using, it must be combined with software version, pseudopotential and specific system verification.
Why Describe Electrons with a Wavefunction?
WAVEFUNCTION & UNITS
The starting point of computational chemistry is to find out how the energy changes with the structure given the atomic species, the number of electrons and the geometric arrangement. It does not need to first specify a spring for each chemical bond, but instead tries to obtain bonding, repulsion, magnetism and other results from the quantum behavior of electrons. First answer three questions: what function are we looking for, what units do the numbers we find have, and what can be compared with experiments.
A wavefunction is not an electron's path around a nucleus
For a non-relativistic electron, the space state can be written as ψ(r,t). r represents position, t represents time; ψ is generally a complex number, and the directly drawn positive and negative or color represents amplitude and phase, which cannot be understood as positive and negative charges. The probability of finding an electron in a small volume d³r is |ψ(r,t)|²d³r, so the normalization condition is ∫|ψ|²d³r = 1. The dimension of the three-dimensional normalized wavefunction is length −3/2 ; probability density has dimensions of length −3 ; Probability itself is dimensionless.
With multiple electrons, the real wavefunction is Ψ(x₁,x₂,…,x N ), x contains space coordinates and spin variables, and N is the number of electrons. It describes the joint state of all electrons. You cannot think of each electron as a small ball carrying a name and moving along a fixed route. The common "orbital" φᵢ in actual calculations is a single-electron function that constructs an approximate multi-electron state or auxiliary system. The orbital diagram has explanatory value, but it is not a photographed electron trajectory.
This is the stationary Schrödinger equation: Ĥ is the energy operator, E is the eigenenergy, ℏ = h/(2π) is the reduced Planck constant, m e is the electron mass, ∇² is the second-order differential operator with respect to position, and V is the potential energy. If using SI units, use m for position and J for V and E; you cannot plug Å directly into an expression written in terms of SI constants only. The stationary state also does not mean that the electron is at rest: the time component is exp(−iEt/ℏ), where i² = −1; it causes the overall phase to evolve, while the stationary state probability density does not change with time.
Allowed energies in a one-dimensional box
Let the electron locality be 0 < x < The ideal region of L: V = 0 inside and infinite external potential energy at both ends. The boundary conditions require ψ(0) = ψ(L) = 0, so only an integer number of half-waves can be put down. The nth stationary state is:
x and L must use the same unit of length, and n is a dimensionless positive integer. The dimension of one dimension ψ is length −1/2 . n = 0 results in a function that is zero everywhere and cannot be normalized, so the lowest energy is not zero either. As the box becomes narrower, the allowed space changes faster and the kinetic energy increases; this can help understand quantum confinement. However, real quantum dots also have finite barriers, effective masses, Coulomb interactions and material structures, and their absorption spectrum cannot be predicted directly by applying this formula.
| quantity | Commonly used units and conversions | Be careful when reading output |
|---|---|---|
| length | 1 Å = 10−10 m;1 bohr ≈ 0.529177 Å | VASP structures and forces are usually expressed around Å |
| Energy | 1 Hartree ≈ 27.211386 eV; 1 eV/particle ≈ 96.485 kJ/mol | eV/unit cell, eV/atom and eV/reaction are not the same normalization |
| force | eV/Å | It is the negative derivative of energy with respect to position and cannot be mixed with the eV threshold |
| thermal energy | kB ≈ 8.61733×10−5 eV/K | 300K of k B T ≈ 0.0259 eV, which is the energy scale but not the upper limit of all potential barriers |
Exercise: When the box width shrinks from 1 nm to 0.5 nm, how do the ground state energy and wavefunction nodes change?
E₁ and L −2 Proportional, so from about 0.376 eV to about 1.504 eV. The ground state has no nodes inside the box, and the two ends are still nodes. Size reduction changes the energy and amplitude scales without changing the ground state into a second excited state. Can be further used ∫₀ L (2/L)sin²(πx/L)dx = 1 Verify normalization.
Extended learning: MIT 5.61: Free particles and particles in one-dimensional boxes .
From Electronic Energy to Potential Energy Surfaces
BORN–OPPENHEIMER & PES
Finding the primary electron energy of a water molecule will not automatically obtain its equilibrium bond length. The O—H distance and bond angle you gave are the external conditions of this electronic problem; changing these coordinates will lead to another electronic problem. After connecting many results, the potential energy surface of "energy changes with structure" appears.
Why can we temporarily hold nuclei fixed?
Atomic nuclei are usually much heavier than electrons, and electrons and nuclei move on different time scales. The Born–Oppenheimer approximation treats the two separately: first solve the electron problem under a set of fixed nuclear coordinates R, and then let the nucleus move on the potential energy surface generated by the corresponding electronic state. "Fixed" here is the condition when solving for an electron step, and does not deny the vibration of the molecule. Common ground state DFT structure optimization is performed within this framework.
R stands for all nuclear coordinates, R I is the position vector of the I-th core; E elec is the electron energy under fixed nuclear conditions, E NN is the internuclear repulsion energy, E BO is the total potential energy, the unit can be eV; force F I The unit of is determined by the coordinate unit. If R is in Å, the force is eV/Å. The "total energy" of actual software usually already includes inter-core interactions and cannot be added manually. This article breaks these terms apart for the purpose of understanding the physics.
Geometry optimization has an inner and an outer loop
- The outer layer is given the atomic coordinates and unit cell, and the inner layer is self-consistent to find the electronic state.
- Calculate the energy, atomic forces, and stresses required for the task from the converged electronic states.
- The optimizer proposes the next set of coordinates or unit cells based on these quantities.
- Repeat until the force, stress, or energy change for the selected degree of freedom meets the criteria.
When the inner layer is not converged, the force received by the outer layer may cause obvious noise in the band structure. Blindly increasing the number of geometric steps at this time is like finding the bottom of the valley along an ever-changing map. In contrast, even if the electronic solution is accurate at each step, the optimizer usually only finds a local minimum near the initial structure; starting from one adsorption site is not guaranteed to find another site through the high potential barrier by itself.
Curvature, vibrations, and transition states
Near the equilibrium position, the energy can be approximately expanded into a quadratic function. For an effective vibration coordinate q, we can write E(q) ≈ E₀ + ½κq²: κ is the force constant, if q is expressed in Å, and κ is expressed in eV/Ų. The vibration angular frequency satisfies ω = √(κ/μ), and μ is the effective mass of the mode; when using this formula, energy, length, and mass must be converted into a consistent set of units. Real molecular vibrations arise from the eigenproblem of the mass-weighted Hessian, which is a matrix of second derivatives of energy with respect to coordinates.
Negative curvature along a direction gives negative ω², which is represented as an imaginary frequency. A molecular ground-state minimum should have no genuinely unstable modes, apart from zero modes associated with translation and rotation. A first-order transition state is expected to have one unstable mode along the reaction direction. Small imaginary frequencies can also arise from geometry, force precision, supercell size, or violations of acoustic sum rules. Inspect the displacement pattern and convergence trends before interpreting a negative plotted frequency as material instability.
The single-potential-energy-surface Born–Oppenheimer picture requires extension when electronic states approach or cross, or when strong nonadiabatic transitions occur. Ground-state optimization cannot replace excited-state dynamics in photochemistry. Tunneling and zero-point motion of light nuclei are not automatically included in classical geometry optimization. A harmonic oscillator, for example, has zero-point energy ½ℏω. Calling the minimum electronic energy an experimental energy that includes all thermal effects would therefore be misleading.
Exercise: After the optimization is completed, the maximum force is very small, but there is still an obvious imaginary frequency in a certain direction. How to judge?
First confirm that the electronic and geometric thresholds are strict enough and that there are no unexpected constraints on the structure before looking at the mode shapes. If the displacement corresponds to a clear structural distortion, it can be perturbed in the positive and negative directions and then optimized separately to see if it can be reduced to a new structure; if it is an overall translation close to zero or a long-wave mode of layered materials, check the supercell, sampling and numerical errors first. "The program exits normally" alone is not enough to determine the nature of the stagnation point.
Extended learning: MIT: Born–Oppenheimer approximation and calculus of variations . For actual structural optimization input, see Getting Started with VASP , vibration and path issues see Advanced tasks .
Variational Methods, Basis Sets, and Electron Correlation
VARIATION, BASIS & CORRELATION
The ideal model of the hydrogen atom can be solved analytically, but the multi-electron system quickly becomes complicated due to the interaction between electrons. Practical algorithms often reduce the problem in two steps: first specifying what multi-electron states are allowed, and then using a finite number of mathematical functions to represent these states. The former is a physical or wavefunction approximation, the latter is a basis set truncation; the two cannot be confused.
The variational principle provides an optimization target
Ĥ is the Hermitian Hamiltonian of the chosen quantum model, Ψ is a non-zero trial wavefunction that satisfies boundary and particle exchange requirements, and E₀ is the exact ground state energy of the same Hamiltonian. The brackets represent the inner product of all electron space and spin variables, and the denominator is normalized. When Ψ is expanded into exact energy eigenstates, the average energy is the result of weighting the eigenvalues with non-negative probabilities, so it cannot go below the lowest eigenvalue. This also explains why "increasing expressiveness and minimizing energy" can gradually improve many wavefunction methods.
But please pay attention to the two conditions of "same Hamiltonian" and "true minimization". Changing pseudopotential, functional, element or electron number, it is no longer the above simple upper bound comparison. There is also no general guarantee that the total energy given by the approximate DFT functional must be higher than the real energy. Even if the calculation of the same method stops at different local electronic solutions, a larger basis set cannot be used to ensure that the number obtained each time becomes monotonically smaller.
A basis set converts functions into expansion coefficients
χμ is the basis function, C μi is the expansion coefficient and φᵢ is the i-th orbital. H μν = ⟨χμ|Ĥeff|χν ⟩ is the effective one-electron Hamiltonian matrix, and S μν = ⟨χμ|χν is the overlap matrix, and ε is the orbital energy diagonal matrix. If the basis function is normalized, S is dimensionless, and H and ε use energy units; S will appear for non-orthogonal basis, and S will be the identity matrix for orthogonal basis.
The local atomic basis set places the function near the atoms and is suitable for understanding bonds in the language of molecular orbitals. Double ζ and triple ζ represent more flexible radial expressions, the polarization function allows the electron cloud to change shape, and the dispersion function helps describe the loose distribution of electrons in space. Plane waves expand the periodic function into waves throughout the unit cell, and rely on the cutoff energy system to increase resolution. Both are just tools for expressing functions. We cannot equate "Gaussian-type functions" with a specific electron theory, nor can we assume that period calculations must only use plane waves.
Hartree–Fock includes exchange, but not all correlation
Electrons are fermions: exchanging the coordinates of two electrons must reverse the sign of the many-electron wavefunction. Hartree–Fock (HF) enforces this antisymmetry with a Slater determinant of spin orbitals, then optimizes the orbitals variationally. Each electron moves in the mean field of the others. The determinant produces exchange terms that constrain the joint distribution of same-spin electrons. Exchange is not a mechanical process of two small balls swapping positions, and HF does not ignore electron–electron repulsion.
In the context of the same non-relativistic Hamiltonian, complete basis set and ground state, E is usually corr = Eexact − EHF It is called correlation energy and it is not positive. It is also possible to define the correlation contribution relative to the HF in a finite basis set, but the range must be stated. Dynamic correlation describes the detailed linkage of electrons in order to avoid each other; nearly degenerate or static correlation requires multiple electron configurations to be described together. For example, when H₂ is pulled very far, a single restricted shell determinant is difficult to maintain the correct limit of two neutral atoms at the same time.
Increasing basis functions reduces basis set errors but does not automatically turn single determinant HF into an accurate multi-electron theory. MP2, configuration interaction, coupled-cluster methods, multi-reference methods, etc. add correlation from different directions, with different applicable systems, costs and failure scenarios. Understanding this level explains why sometimes "calculating longer" simply solves a model more accurately that is still inappropriate.
Exercise: The HF energy of the same water molecule decreases significantly as the basis set increases. Does this mean that the electronic correlation has been completed?
No. You mainly improve the representation of orbitals, still constrained by a single determinant during convergence to the HF limit. To evaluate the correlation, it is necessary to compare appropriate correlation methods with a clear basis set and Hamiltonian, or use a known high-precision benchmark; it is also necessary to avoid mixing changes in different structures into method comparisons.
Extended learning: MIT Variation principle and matrix form , basis set , Electronic correlation .
DFT: A Many-electron Theory Based on Density
DENSITY FUNCTIONAL THEORY
An N-electron space wavefunction depends on 3N space variables, while the electron density n(r) depends only on three-dimensional position. Density somewhere represents the average number of electrons per unit volume; integrating over the entire space for an isolated system, or integrating over a unit cell for a periodic system, the corresponding number of electrons is obtained. The common unit is the number of electrons/ų, and the electron charge density must be multiplied by the negative basic charge. Density does not contain an identity tag for each electron.
A functional maps a function to a number
An ordinary function f(x) maps one number to another; an energy functional E[n] maps an entire spatial density distribution to an energy. Ground-state DFT states that, under appropriate conditions, the ground-state electron density determines the external potential and the system's ground-state properties. The exact ground-state energy follows from a variational search over admissible densities. This existence theorem does not supply a practical exact functional, nor does it make every arbitrarily drawn density physically admissible.
Each term is energy: T s is the kinetic energy of the auxiliary non-interacting electron system that generates this density; v ext It is the potential energy given to an electron by external sources such as the nucleus; E H is the classical Coulomb energy of density; E xc is the exchange–correlation energy, including the difference between real and auxiliary kinetic energy and the electron interaction beyond the classical Coulomb; E NN for internuclear repulsion. The exchange–correlation term is not an optional “little correction after the fact”, it takes on key multi-electron physics.
Kohn–Sham theory introduces auxiliary orbitals
This formula uses atomic units: ℏ, m e The combination of related constants such as basic charge and basic charge is 1, the unit of length is bohr, and the unit of energy is Hartree. v H with v xc They are the Hartree potential and the exchange–correlation potential; φᵢ is the normalized auxiliary orbital, εᵢ is the eigenvalue, and fᵢ is the occupation number. If written in terms of space orbitals and regardless of spin, fᵢ can be between 0 and 2; if i already includes spin, the maximum is 1. The periodic system also needs to include k points and their weights in the summation.
The density determines the effective potential; the potential determines the orbitals; the orbitals generate a new density. These equations therefore require a self-consistent solution. Hybrid functionals often use a generalized Kohn–Sham form with nonlocal operators, so the local-potential expression above is not the only possible implementation. The total energy is not simply the sum of occupied orbital eigenvalues εᵢ: that would double-count some interactions and omit the corresponding corrections and internuclear energy.
Choosing a functional means choosing an approximation
| approximate type | The information used or the main idea | Questions to ask |
|---|---|---|
| LDA | Organize exchange–correlation with local density | Are the bonding and structures studied suitable for this approximation? |
| GGA, like PBE | Add density gradient | Are there known deviations in structure, energy differences, magnetism, and band gaps? |
| meta-GGA | Introduce additional local information such as kinetic energy density | Is the mesh, PAW data compatible and converged with the version used? |
| hybrid functional | Mix in a certain form and range of HF exchanges | Are the exchange ratios, shielding range, costs, and nature of the target appropriate? |
| DFT+U | Add correction to selected local subspace | Are the definitions, values, orbital subspaces and references of U clear? |
| Dispersion correction or non-local vdW functional | Added description of long-range related effects | Is the method repeatedly counted in a certain item, and is it suitable for adsorption and interlayer interaction? |
This table is not an accuracy ranking. More expensive methods may also perform poorly in certain types of reactions, metals, or strongly correlated systems; the errors of different approximations sometimes cancel out. You should perform benchmarks or cross-checks on your own target quantities. You cannot adjust U at will in order to fit a certain experimental number, and then treat the parameter adjustment results as independent predictions.
Exercise: The band gap calculated by PBE is 1.0 eV, and the experimental absorption edge is 1.8 eV; can all conduction bands be moved up by 0.8 eV as a whole?
The empirical scissors correction can be explicitly employed for specific demonstration purposes, but it does not resolve the electronic structure, nor does it prove that the effective mass, defect energy levels, optical intensity, and interface band structure arrangement are all correct. First verify the experimental material phase, temperature, defects and measured quantities; then select the appropriate hybrid functional, quasiparticle or exciton method and benchmark according to the problem, and do not regard fitting as verification.
Theoretical starting point: Kohn and Sham (1965) ; Functional type see VASP exchange–correlation functional description ;For the conceptual boundary of the energy gap, see Perdew and Levy (1983) .
Self-consistency, Occupations, and Spin
SCF, OCCUPATIONS & SPIN
Solving the Kohn–Sham equation looks like finding the matrix eigenvalues, but the matrix itself depends on the electron density which has not yet been found. This is a feedback problem: first guess the density, construct the effective potential, find the orbital and energy, and then obtain the new density based on the orbital and occupation; when the input and output are consistent enough, it is called self-consistent.
Self-consistency is more than repeating the same equation
- Initialization: For example superimposing atomic densities, or reading charges from a compatible previous calculation with wavefunction .
- Solve the current effective problem: Optimize the orbital step by step to obtain the required number of eigenstates.
- Assign occupations: Meet the total number of electrons and use occupation or smearing methods appropriate to the task.
- Construct new density: Check the energy and residuals, mix old and new information, and continue iterating.
This is a simple linear mixture formula to help understand, j is the number of iterations, α is the dimensionless mixing ratio, usually 0 < α ≤ 1;nin with n out are the input and output densities of the round respectively. The actual mixing and preprocessing of VASP is more complicated than this equation, but the purpose is the same: to avoid over-correction of the old and new charge distributions back and forth. Long unit cells, metallic surfaces, and difficult-to-treat magnetic states are prone to "charge oscillations," and not all problems can be solved by increasing the maximum number of electron steps.
You can map x with toys out = 1 − xin Understanding oscillations: x is a dimensionless scalar, direct substitution will jump between two values, α = ½ then reaches the fixed point ½ in one step. There are many spatial patterns of real density, each requiring a different numerical treatment; this example illustrates why blending is useful and is not a formula for predicting the rate of convergence of real materials.
Occupation smearing and electronic temperature require care
There are both occupied states and empty states near the Fermi energy of metals. Small energy level changes will change the occupation, making k-space integration and self-consistent difficult. Smearing allows the occupation to change smoothly around the Fermi energy. The Fermi–Dirac form is:
Here, counting is based on a single spin state, f is between 0 and 1; ε is the energy of the state, μ is the electron chemical potential, T e is the electron temperature, k BTe Uses the same energy units as ε. Only when this statistical form is adopted, the smearing parameter can be directly interpreted as the corresponding temperature; Gaussian and Methfessel–Paxton smearing are mainly numerical integration tools. The electron temperature is not automatically equal to the ion temperature in molecular dynamics.
In VASP, ISMEEAR selects the mode, and SIGMA gives the eV width. Larger SIGMA is sometimes easier to converge, but may change the target energy, magnetic state, and force; too small may require a dense k-grid. The unknown system can be started by smearing a small-width Gaussian and testing the target quantity. Bandgap systems should not mechanically copy the positive value ISMEAR commonly used for metals; tetrahedral integration requires a suitable mesh and is not a universal option when there is only a Γ point. Free energy, deentropy energy and extrapolated energy with smearing have different meanings and the same and appropriate definitions should be adhered to when comparing.
Spin is part of the solution space
Collinear spin DFT treats n↑ and n↓ separately: total density n = n↑ + n↓, spin density difference m = n↑ − n↓. m is calculated here in terms of electrons/volume, and after integration and multiplication by the corresponding magnetic moment units, the total spin magnetic moment can be discussed. The absence of a net magnetic moment does not mean that there is no magnetism everywhere: different sites in an antiferromagnetic system can have opposite local magnetic moments, which eventually cancel each other out.
ISPIN = 2 allows collinear spin polarization; MAGMOM gives the initial local magnetic moment and affects the preserved symmetry, it does not permanently lock the final magnetic moment. Multiple ferromagnetic, antiferromagnetic or different initial amplitudes should be tried according to structural and chemical possibilities, and the energy and local state after convergence should be compared. If you study non-collinear magnetism or spin–orbit coupling, you need a more complete setup, which cannot be completed by just adding a line of MAGMOM.
Exercise: The total magnetic moment is 0, can the spin polarization calculation be turned off?
You cannot decide based on this number alone. First check the local magnetic moment and spin density: the net magnetic moment of two antiferromagnets with equal opposite magnetic moments is also 0. Changing it to not allow spin-polarized calculations would change the physical model, possibly changing the energy gaps, bond lengths, and total energy. Also make sure that the supercell can accommodate the magnetic sequence you want to study.
For setting definition see Self-consistent cycle , smearing method , MAGMOM with EDIFF . When encountering non-convergence problems, first press Parameter and diagnostic manual Check the structure and methods, and then adjust the algorithm.
Crystals, Reciprocal Space, and k Points
BLOCH WAVES & BRILLOUIN ZONE
Instead of suspending a small piece of crystal in a vacuum, periodic calculations stipulate that it repeats infinitely along the direction of the crystal lattice. Atoms, defects and adsorbates in the unit cell are also repeated. Therefore "putting only one defect" usually means that each supercell has one defect; increasing the supercell reduces the concentration and interaction of repeated defects, rather than changing the box boundaries of the same finite cluster.
Distinguish the lattice from the atomic basis
a₁, a₂, a₃ are three lattice vectors, with units such as Å; nᵢ is an integer, and R is any translation vector. fᵢ is the dimensionless fractional coordinate of an atom and r is its Cartesian position. A primitive cell is the smallest translational repeating unit. A regular unit cell may be more convenient for showing symmetry, while a supercell is composed of multiple primitive cells. With an oblique unit cell, a one-dimensional change in the fractional coordinates is not equivalent to a separate movement along Cartesian x, y, or z; looking at POSCAR you must also look at the lattice vectors.
Reciprocal space describes spatial variations of waves
bⱼ is the inverse lattice vector, the unit is Å −1 ; δᵢⱼ is 1 when i = j, otherwise it is 0; hᵢ is an integer, and G is a reciprocal vector combination. Physics conventions that include 2π are used here and should not be confused with the conventions of some crystallographic tools that do not include 2π. For the direction with orthogonal length a, the reciprocal space scale is 2π/a, so the real space repetition length becomes larger, and the corresponding reciprocal space repetition scale becomes smaller.
The one-electron solution in the periodic potential can take the Bloch form. It does not mean that the wavefunction is exactly the same in adjacent unit cells, but that adjacent unit cells can differ by one phase:
n is the band structure number, k is the Bloch wave vector (Å −1 ), r and R are represented by Å, and the dot product k·r is dimensionless; u is a function with lattice periodicity. Independent values of k can be restricted to the first Brillouin zone. The Γ point is k = 0. The crystal momentum label ℏk is useful for understanding translational symmetry, but it cannot be directly regarded as the ordinary mechanical momentum of electrons at each moment in the crystal.
Total energies require Brillouin-zone integration
A(k) represents a certain k-related quantity to be integrated, Ω BZ is the Brillouin zone volume (Å −3),wk is a dimensionless weight. A regular k-grid approximates this integral, and crystal symmetries merge equivalent points into irreducible regions, with corresponding weight increases. The band structure diagram is sampled along the Γ—X contour symmetry lines, shows only a few directions, and cannot replace the volume integration grid for total energy and density.
The k-grid requirement depends on the rate of change of the reciprocal space. Small primitive cells usually require more points; metals change quickly near the Fermi surface and are often more difficult to converge than insulators of similar size. For orthorhombic unit cells, generally fewer points can be taken along the longer direction of real space. Non-orthogonal lattice must also take into account the reciprocal lattice vector, grid center and crystal symmetry. The "reciprocal length ratio" cannot be regarded as a strict rule for all crystal systems.
It is also important to understand supercell folding: when the same ideal crystal is expanded twice in each direction, the Brillouin zone shrinks, and part of the original k information will be folded into more band structures. Equivalent physical results should be obtained under compatible sampling and the same numerical conditions, but "any large supercell only needs Γ point" is not true, especially for metals, shorter directions and small energy differences.
Exercise: The cubic primitive cell uses an 8×8×8 grid. After constructing the 2×2×2 supercell, what grid can you start testing from?
4×4×4 is a reasonable starting point for maintaining similar reciprocal space point spacing, since each real space direction becomes twice as large. It should also be confirmed that the grid centers and folds are compatible, that the structure is still the original crystal, and that the target properties converge. If defects are added or the magnetic sequence is changed, re-verification is required, and the same accuracy cannot be declared based on point scaling alone.
Reference implementation: VASP KPOINTS and mesh symmetry , Brillouin zone points . band structure and density of states workflow see electronic structure tasks .
Plane Waves, PAW, and the Energy Cutoff
PLANE WAVES & PAW
After Bloch's theorem, the periodic part u can be expressed by Fourier series. Similar to synthesizing complex sounds using musical notes, we synthesize tracks using plane waves of different spatial frequencies. High-frequency components can express finer spatial changes, but retaining more components makes the calculation more expensive.
Ω is the unit cell volume, G is the inverse lattice vector, n and k still represent the band structure and Bloch wave vector, and C is the expansion coefficient; the tilde on the wavefunction reminds us that this can refer to the smooth auxiliary orbital optimized by the PAW calculation. If Ω is in ų, the unit of the factor before normalization is Å −3/2 . An infinite number of G is not feasible on a computer, so it is truncated by the kinetic energy of the corresponding plane wave:
ENCUT of VASP is E cut , unit eV. If the wave vector is expressed in Å −1,ℏ²/(2me ) ≈ 3.80998 eV·Å². Increasing ENCUT expands the range of allowed wavevectors, improving spatial resolution; it does not remove actual electrons with energies above this value, nor is it an upper limit on how large a bandgap is allowed to be calculated. The number of unoccupied bands is controlled by another type of setting, not the same thing as the number of basis functions.
Why not resolve every rapid variation near a nucleus?
When close to the nucleus, the valence electron orbital is affected by the orthogonal condition between the nucleus and the core state, and the space changes rapidly. Direct representation by plane waves requires a large number of high-frequency components. PAW (Projected Additive Wave) uses smooth auxiliary functions to undertake large-scale representation, and then uses local information near the atoms to reconstruct the required full electronic features. It includes processing such as projection, local partial waves, and compensation. It cannot be understood as just flattening the area near the nucleus and ignoring errors.
The VASP common PAW data set also uses a frozen core approximation: some inner electrons remain in the core state description when the data was generated, and the valence electrons are explicitly solved for. Different data sets may include semicore electrons in the valence shell, or use different hardnesses and reference configurations. Therefore two POTCARs of "the same element" do not necessarily have the same number of electrons, energy basis, and accuracy requirements. The specific data set identification, version, and valence electron configuration should be recorded, and data for related species should be consistent when comparing energy differences.
How the cutoff, mesh, and cell determine computational cost
Approximating the available wave vectors as lattice points in a sphere, under a three-dimensional fixed truncation, the number of plane waves roughly increases with ΩE cut3/2 growth. This is a scalable estimate of the basis set size, not a strict formula for total runtime; the true cost also depends on the number of electrons, number of unoccupied bands, k-points, algorithm, and parallelism. Although adding a vacuum to an isolated molecule does not add atoms, it does increase the volume of the unit cell and the space that needs to be represented, so it will significantly increase the amount of memory and calculation.
The real space FFT grid is responsible for converting between real space and reciprocal space; settings such as PREC affect the relevant numerical grid, and it cannot be considered that just increasing ENCUT will cover all grid errors. For stress and unit cell optimization, the finite basis set will also produce Pulay stress as the unit cell changes; it will always appear to be stable, but there is no guarantee that the pressure or volume will also achieve the required accuracy. Checking the final structure with uniform and fully converged parameters after optimization can help identify such problems.
How to perform a meaningful ENCUT test
Fixed structure, k-grid, electronic methods, and other conditions, select an incremental set of ENCUTs spanning the POTCAR recommended scales, and record the target energy difference, force, or stress. When using a mix of elements, keep an eye out for the most demanding species among them. The final threshold is determined by property requirements and is not always 400 eV or 520 eV. If you want to compare the formation energies or adsorption energies of different systems, you should use explicit cutoff settings that have been validated for the entire set of comparison objects to avoid using different implicit default values for each.
Exercise: Why would it still be slower if we increase the vacuum layer from 15 Å to 25 Å, keeping the same number of atoms?
As the volume increases, the plane wave and FFT mesh sizes available within a given cutoff generally increase; the orbitals need to be represented in a larger space. At the same time, increasing vacuum may also change the long-wavelength charge pattern and self-consistent difficulty. The required size should be determined based on the convergence of the target quantity to the vacuum, rather than increasing it infinitely. The k-point in the vacuum direction can often be maintained at one, but this does not eliminate the cost of volume.
Further reading: PAW principle , ENCUT definition , Pulay stress .
From a Finished Run to a Reliable Conclusion
RELIABILITY & SCIENTIFIC JUDGMENT
The first calculation report should not only have a band structure diagram and a minimum energy. What really needs to be answered is: what is compared, what physical model is used, how big the numerical error is, and how strong the current evidence can support the conclusion. Writing these questions clearly is closer to reliable research than adjusting all parameters to the maximum.
Distinguish at least four types of error
| level | Typical questions | Verification method |
|---|---|---|
| physical model | Wrong crystal phase, magnetic order, surface termination, defect concentration; missing solvent or temperature | Compare experimental conditions and construct reasonable alternative models |
| electronic approximation | functional, frozen core, strong correlation or insufficient excited state description | Benchmarks, cross-checking of different methods, known limits |
| Numerical discreteness | ENCUT, k-point, vacuum, supercell, FFT grid insufficient | Increase the accuracy item by item to see if the target quantity is stable. |
| Solving and Sampling | SCF does not converge, local minimum, MD is too short, initial state is single | Check reasons for stopping, multiple starting points, time and size tests |
EDIFF from 10 −5 Change to 10 −8 eV, mainly tightens the stop criteria for electronic solution, and will not fix the wrong crystal phase or unsuitable functional. When EDIFFG is a negative value, the geometric stop is determined according to the atomic force norm, and the corresponding dimension is eV/Å; small energy change, small force, and target property convergence are three different checks. The actual task also requires confirming that the processing of fixed atoms, unit cells, and constrained degrees of freedom conforms to the settings.
Converge the target property
First clarify the allowed variation. For example, if you want to distinguish adsorption configurations that differ by about 0.1 eV, you can set the numerical accuracy target to be significantly smaller than this difference before conducting the test. Below is a set Artificially constructed teaching data, not real VASP calculation results , used to demonstrate the recording method:
| Test to fix remaining conditions | Adsorption energy/eV | Change from previous gear / eV |
|---|---|---|
| ENCUT = 400 eV | −0.620 | — |
| ENCUT = 450 eV | −0.648 | −0.028 |
| ENCUT = 500 eV | −0.653 | −0.005 |
| ENCUT = 550 eV | −0.654 | −0.001 |
This only shows that the toy data is stable against truncation, and does not prove that the k-point, coverage, or functional errors are also small. The next step should be to fix sufficient cutoffs, test the k-grid, and review the combination parameters. If the difference between two configurations is only 0.003 eV, but the convergence test still drifts by 0.01 eV, there is not enough basis to give a robust ranking of advantages and disadvantages. Avoid mistaking an output with many decimal places for a result with the same number of significant digits.
Here E slab+A is the total energy of electrons in the adsorption system, E slab is the corresponding cleaning surface energy, E A is the adsorbate energy in eV using a well-defined reference state; according to this sign convention, negative values indicate a more favorable electron energy relative to the selected reference. It must be stated whether each structure, atomic number, magnetic state, and calculation setup are comparable. Electron adsorption energy does not automatically equal experimental adsorption free energy; zero-point energy, entropy, coverage, solvent and electrode conditions may all affect comparisons.
AI-assisted workflows still need physical validation
A common direction is machine learning potential: training models with energy, force, and stress calculated from high-quality electronic structures, allowing the model to quickly estimate the potential energy surface and derivatives of new structures, thereby expanding the time and size that molecular dynamics can achieve. Active learning prioritizes supplementing first-principles data for uncertain or novel structures. In this way, quantum computing no longer needs to be re-executed at every step, but it still bears the training label and verification basis.
Two levels of error need to be asked simultaneously: whether the model approximates the DFT used for training, and whether this DFT is appropriate for the real research problem. Being well trained on a perfect crystal at a certain temperature does not guarantee that the model can describe melting, bond breaking, defect migration, or new elements. If the training set and test set contain a large number of adjacent time frames, optimistic errors may be obtained due to too similar structures. It is more meaningful to conduct verification based on independent trajectories, different structures and target processes.
AI can also screen candidate materials, suggest reaction pathways, or generate initial structures. The generated result is first a candidate; verification via DFT is computational evidence; further experiments answer whether it actually exists or works under specified conditions. When writing a report, use clear wordings such as "predicted candidate", "optimized by a certain method", "simulated under certain conditions", and "experimental observation", etc., which is more reviewable than writing the conclusions of each layer as "discovered high-performance materials".
A checklist for your first VASP task
- Write a scientific question to be answered and identify a specific comparable target quantity.
- List the structural origin, composition, number of electrons, period boundaries, magnetic states and environmental assumptions.
- Document software version, PAW data, functional, all key inputs and post-processing definitions.
- Keep convergence data, exception logs, failed attempts, and final output instead of just pretty graphs.
- Check the most critical conclusions against reasonable references, alternative configurations, or another method.
Comprehensive exercise: The machine learning potential reduces the test energy error per atom to 2 meV. Can you directly publish "the diffusion barrier is accurate to 2 meV"?
No. The average error per atom is not the error of the local migration barrier, and the training and test distributions may not cover the transition state; the errors at both ends and on the path may be related or amplified. The relative energies and forces of the actual migration paths need to be verified and errors in the original electron theory, supercell and path solution need to be checked. Also report how the test set composition is separated from training, rather than just extracting a minimum error metric.
Reference: Geometric stopping conditions , VASP machine learning potential theory , Machine learning potential verification suggestions . Next step to enter Getting Started with VASP ;You can also return to Learning Path Select a topic.
Primary Sources & Further Reading
READ THE PRIMARY SOURCES
The main text is an independently organized study instruction. For label meanings and version differences, refer to the official documentation corresponding to your program version.
- MIT 5.61 Lecture 5: Free Particle and Particle in a 1D Box
- MIT 5.61: Born-Oppenheimer approximation and variational principle
- MIT 5.61 Lecture 24: Matrix Mechanics and Variational Principle
- MIT 5.61 Lecture 28: Basis Sets
- MIT 5.61 Lecture 29: Electron Correlation
- Kohn & Sham (1965): Self-Consistent Equations Including Exchange and Correlation Effects
- VASP Wiki: Exchange-correlation functionals
- Perdew & Levy (1983): Physical Content of the Exact Kohn-Sham Orbital Energies
- VASP Wiki: Self-consistency cycle
- VASP Wiki: Smearing technique
- VASP Wiki: MAGMOM
- VASP Wiki: EDIFF
- VASP Wiki: KPOINTS
- VASP Wiki: Integrating over all orbitals
- VASP Wiki: Projector-augmented-wave formalism
- VASP Wiki: ENCUT
- VASP Wiki: Pulay stress
- VASP Wiki: EDIFFG
- VASP Wiki: Machine learning force field theory
- VASP Wiki: Best practices for machine-learned force fields