Skip to content

VASP

Jump to bottom
Lee Burton edited this page Sep 29, 2026 · 9 revisions

⚛️ VASP — fundamentals

This tutorial introduces the basics of VASP (the Vienna Ab initio Simulation Package): the four input files, the density-functional physics VASP solves, and the INCAR settings you will meet most often. It uses a simple FCC silicon crystal as a teaching example. 🧪

📘 Teaching page. The examples on this page are simplified for learning. They are not BMD Compute production settings. BMD Compute is authoritative for how BMD VASP calculations are generated and executed; to see the INCARs it actually generates, and why, read BMD Compute INCARs.

📂 Required Input Files

Before running VASP, you need four key files in your working directory. Each plays a unique role in the simulation process:

INCAR  POSCAR  KPOINTS  POTCAR
  • INCAR — what to calculate and how (the settings).
  • POSCAR — the crystal structure: lattice vectors and atomic positions.
  • KPOINTS — how to sample the Brillouin zone.
  • POTCAR — the pseudopotential (PAW dataset) for each element.

All four are required for VASP to run properly. Let's look at them one by one. 🔍

🔬 The Physics VASP Solves

VASP is a plane-wave implementation of Kohn–Sham density functional theory (DFT). Knowing what is being approximated, and where, is what lets you choose and check its settings.

Kohn–Sham DFT in one paragraph

Within the Born–Oppenheimer approximation the nuclei are fixed point charges, and the problem is the ground state of the interacting electrons in their potential. Hohenberg and Kohn showed that the ground-state energy is a functional of the electron density n(r) alone. Kohn and Sham made this practical by mapping the interacting system onto non-interacting electrons moving in an effective potential,

veff(r) = vext(r) + vHn + vxcn,

the external (ionic), Hartree and exchange–correlation potentials. Solving the single-particle Kohn–Sham equations gives orbitals ψi, and the density is built from the occupied orbitals. Because veff depends on n, which depends on the orbitals, the equations must be solved self-consistently.

Exchange–correlation: the approximation you choose

All the many-body physics beyond the classical Hartree term is in the exchange–correlation functional Exc[n], which is not known exactly. The choice of functional is a physical approximation, not a numerical setting:

  • LDA depends only on the local density.
  • GGA also depends on its gradient. PBE (GGA = PE in VASP) is the most widely used GGA in solids.
  • meta-GGAs (for example r²SCAN) add the kinetic-energy density.
  • Hybrid functionals replace a fraction of semilocal exchange with exact (Hartree–Fock, or Fock) exchange, which is built from the orbitals rather than the density. HSE06 uses a fraction of 0.25 (AEXX) and applies it only at short range: the Coulomb interaction is split with a screening parameter μ (HFSCREEN, in Å⁻¹), and only the short-range part is treated exactly.

Semilocal functionals such as PBE commonly underestimate semiconductor band gaps, in part because of self-interaction error. Hybrids often improve gaps, but the exact-exchange operator is non-local and couples pairs of orbitals at pairs of k-points, so it is far more expensive.

Plane waves and the cutoff energy

In a crystal, Bloch's theorem lets each orbital be written as a plane wave times a lattice-periodic function, which VASP expands in plane waves labelled by reciprocal-lattice vectors G:

ψnk(r) = ΣG cnk,G ei(k+G)·r.

The sum is truncated by keeping only plane waves with kinetic energy ħ²|k+G|²/2m below the cutoff ENCUT. This is a basis-set truncation: raising ENCUT enlarges the basis and, by the variational principle, lowers the total energy towards its complete-basis limit. The number of plane waves grows roughly as (cell volume) × ENCUT3/2, so cost rises quickly with both.

Near the nuclei, valence orbitals oscillate rapidly and would need an enormous number of plane waves. VASP avoids this with the projector augmented-wave (PAW) method: core electrons are frozen, and the valence orbitals are represented by smooth functions plus atom-centred corrections. The POTCAR supplies these PAW datasets, each with a recommended minimum cutoff, ENMAX.

Because the plane-wave basis is tied to the cell rather than to the atoms, moving atoms does not change the basis, but changing the cell does. That distinction matters when the cell is relaxed; see BMD Compute Workflows.

Reciprocal-space sampling

Densities and energies involve integrals over the first Brillouin zone. VASP replaces them by weighted sums over a finite k-point mesh (usually Monkhorst–Pack or Γ-centred), and uses crystal symmetry to reduce the mesh to the irreducible k-points it actually computes (listed in IBZKPT). Larger real-space cells have smaller Brillouin zones and need fewer k-points. Metals are the hard case: their occupations change abruptly at the Fermi surface, so the integrand is discontinuous and converges slowly with mesh density.

Occupations and smearing

At zero temperature each state is either full or empty. On a finite k-mesh, that step function makes metals converge badly, so VASP can broaden the occupations:

  • Gaussian (ISMEAR = 0) or Fermi–Dirac (ISMEAR = -1) smearing replaces the step by a smooth function of width SIGMA. The quantity then minimised is a free energy that includes a fictitious electronic entropy term; OUTCAR also reports an estimate extrapolated to SIGMA → 0.
  • Methfessel–Paxton (ISMEAR = 1, 2) is designed for metals; it can give unphysical occupations in gapped systems.
  • The tetrahedron method with Blöchl corrections (ISMEAR = -5) instead interpolates the bands linearly within tetrahedra of the k-mesh. It gives accurate total energies and densities of states without a smearing width (SIGMA is ignored), but it needs a uniform mesh and is not suited to relaxing metals.

The self-consistent field (SCF) cycle

For a fixed structure, VASP iterates:

  1. start from a trial density (for example a superposition of atomic densities, ICHARG = 2);
  2. build veff[n];
  3. partially diagonalise the Kohn–Sham Hamiltonian iteratively for the lowest bands (Davidson or RMM-DIIS, chosen with ALGO);
  4. form a new density from the occupied orbitals and mix it with previous densities (Pulay/Broyden mixing) for stability;
  5. repeat until the total-energy change between steps is below EDIFF, or until NELM steps have been used.

Each pass is an electronic step. A calculation that hits NELM has not converged, and its energy and forces should not be trusted.

Forces, stresses and ionic relaxation

Once the electrons are self-consistent, VASP computes the forces FI = −∂E/∂RI on each nucleus (via the Hellmann–Feynman theorem plus PAW corrections) and the stress, the derivative of the energy with respect to strain. Forces are sensitive to residual SCF error, which is why relaxations need a tight EDIFF.

A relaxation (geometry optimisation) uses these derivatives to move the nuclei, and optionally the cell, downhill on the Born–Oppenheimer energy surface. Each move starts a new SCF cycle; each move is an ionic step. The loop stops when the EDIFFG criterion is met or after NSW steps. A static (single-point) calculation has only the SCF cycle (NSW = 0).

You can follow both loops in OSZICAR: each electronic step is one line, and each completed ionic step is summarised on a line starting with its step number.

🧠 INCAR — The Brain of VASP

The INCAR file is the main input file. It controls what kind of calculation VASP performs and how it does it. 🧾 It contains various tags that define algorithms, accuracy settings, and convergence criteria.

Default values exist, but it's best to set parameters yourself. See the VASP INCAR wiki for full documentation.

Example: FCC Silicon INCAR

📘 Teaching example, not production policy. The FCC-Si input below is a simplified example for learning what each file and tag does. It is not BMD Compute production policy. For the INCARs BMD Compute actually generates, see BMD Compute INCARs.

System = fcc Si
ISTART = 0 ; ICHARG = 2
EDIFF = 1e-4
ENCUT = 240
ISMEAR = 0 ; SIGMA = 0.1
NCORE = 4
KPAR  = 1
NSW   = 0

This is a static calculation (NSW = 0) that starts from scratch (ISTART = 0 reads no old wavefunctions; ICHARG = 2 starts from a superposition of atomic charge densities). The loose EDIFF and low ENCUT keep it quick; the sections below explain what each choice means.

There are over 600 possible parameters! Explore them all in the VASP Manual.

⚡ Plane-Wave Cutoff: ENCUT and PREC

ENCUT (in eV) is the basis-set truncation described above.

  • Use at least the largest ENMAX among your POTCARs, and more when you need accurate stresses or relaxations of the cell.
  • Convergence test: compute the total energy at increasing ENCUT (for example in steps of 50 eV) and stop when the energy per atom changes by less than your target accuracy. Converge the quantity you actually need: energy differences and forces often converge faster than absolute energies.
  • Compare energies only between calculations that use the same ENCUT and POTCARs.
  • PREC sets the density of the real-space FFT grids used to represent the density and potentials; PREC = Accurate avoids aliasing errors and is the usual choice for reliable forces and stresses.

🧮 k-Point Sampling

The mesh in KPOINTS (below) is the reciprocal-space sampling described above.

  • Converge the mesh the same way as ENCUT: increase it until the energy per atom stops changing.
  • Scale the mesh inversely with the cell's lattice lengths, so that the k-point spacing, not the number of k-points, stays roughly constant.
  • Metals need much denser meshes than insulators.

🌡️ Smearing: ISMEAR and SIGMA

  • ISMEAR = 0 — Gaussian smearing. A safe general choice; for semiconductors and insulators use a small SIGMA (around 0.01–0.05).
  • ISMEAR = -5 — the tetrahedron method with Blöchl corrections. Accurate total energies and densities of states, but it needs a uniform k-point mesh, SIGMA is ignored, and it is not suitable for relaxing metals.
  • ISMEAR = 1 or 2 — Methfessel–Paxton smearing, for metals. Never use it for semiconductors or insulators.

Check the effect of smearing in OUTCAR: the difference between energy without entropy and free energy TOTEN should be small.

🎯 Convergence Criteria: EDIFF and EDIFFG

  • EDIFF — the SCF cycle stops when the total energy changes by less than this (in eV). 1e-4 is fine for learning; production and property calculations typically use 1e-5 to 1e-6.
  • EDIFFG — the ionic loop's stopping criterion:
    • a positive value is an energy criterion (stop when the energy change between ionic steps is below it);
    • a negative value is a force criterion (stop when all forces are below |EDIFFG| in eV/Å), for example EDIFFG = -0.01.
  • EDIFFG only matters when ions move (NSW > 0).
  • NELM is the maximum number of electronic steps; if it is reached, the SCF cycle has not converged.

🧗 Ionic Relaxation: IBRION, ISIF and NSW

  • NSW — maximum number of ionic steps (0 for a static calculation).
  • IBRION — how the atoms are moved: 2 is conjugate gradient (robust), 1 is quasi-Newton (efficient near the minimum), -1 means no movement.
  • ISIF — what may relax: 2 relaxes atomic positions only; 3 relaxes positions, cell shape and cell volume.

A relaxation that stops because it reached NSW has not converged; check the final forces in OUTCAR before using the structure. Relaxing the cell also raises a basis-set subtlety (Pulay stress) that is handled at the level of the whole workflow; see BMD Compute Workflows.

💾 Restarts and Output Files

  • ISTART and ICHARG control whether VASP starts from scratch or reads a previous WAVECAR (orbitals) or CHGCAR (charge density).
  • LWAVE and LCHARG control whether WAVECAR (large) and CHGCAR are written.
  • LORBIT and NEDOS control the projected and total density of states written to DOSCAR, PROCAR and vasprun.xml.

🖥️ Parallelisation: NCORE and KPAR

These settings change how the work is split across processors. They should change speed and memory use, not results.

  • NCORE: Number of MPI ranks that work together on a single band (one band group). The examples on this wiki use NCORE; see the note on NPAR below.
    • Choose a divisor of your ranks per k-point group (e.g. 2–6).
  • KPAR: Parallelization across k-points.
    • Must divide both the number of k-points and the total MPI ranks.
    • Safe to use in VASP 6; differences in energy are only roundoff-level.

📌 Note on NPAR: Older tutorials often set NPAR instead of NCORE. Both control the same part of VASP's parallel decomposition, the way the MPI ranks within one k-point group are split over bands, but they specify it differently:

  • NCORE sets how many ranks work together on each band (the size of a band group);
  • NPAR sets how many band groups there are.

For a given number of ranks per k-point group, choosing one therefore fixes the other (roughly, NPAR × NCORE = ranks per k-point group). Use NCORE, as in our examples, and do not set both in the same INCAR. Combine it with KPAR for efficient scaling.

🧱 POSCAR — The Structure File

POSCAR defines the atomic structure and lattice geometry. It’s often the first file you prepare for a new simulation. You can write it manually ✍️ or download it from databases like the Materials Project. 🌐

Example: FCC Silicon POSCAR 

fcc Si:
 3.9
 0.5 0.5 0.0
 0.0 0.5 0.5
 0.5 0.0 0.5
   1
cartesian
0 0 0

📌 Notes:

  • Lattice constant: 3.9 Å
  • Primitive FCC unit cell
  • 1 Si atom at origin
  • This is a simplified teaching structure. Real crystalline silicon has the diamond structure, with 2 atoms in its primitive cell.

🧮 KPOINTS — Brillouin Zone Sampling

The KPOINTS file tells VASP how to sample k-space, which is key to convergence in electronic structure calculations. ⚡

Example: Monkhorst-Pack Grid

k-points
 0
Monkhorst Pack
 11 11 11
 0 0 0

🧾 Format Breakdown

  1. Line 1: Comment
  2. Line 2: 0 = automatic generation
  3. Line 3: Grid type (e.g. Monkhorst Pack)
  4. Line 4: Mesh size in 3 directions
  5. Line 5 (optional): Mesh shift (usually 0 0 0)

🧠 Use finer meshes for better accuracy—at the cost of computation time. For hexagonal cells, use a Γ-centred grid (Gamma on line 3) rather than Monkhorst–Pack.

🧪 POTCAR — The Pseudopotentials

The POTCAR file contains element-specific pseudopotentials and exchange-correlation data.

  • Must match the order of elements in your POSCAR file ✅
  • Combine all needed elements into one file using cat 🐱

Example (using PAW PBE 64-bit potentials on the POWER cluster):

/bmd/shared/vasp/recommended-potpawPBE64

🔧 Generate POTCAR

cat /path/to/Si/POTCAR > POTCAR

🎉 Next Steps

Congrats! 🥳 You’ve now learned the basics of setting up a VASP calculation with:

  • 🧠 INCAR
  • 🧱 POSCAR
  • 🧮 KPOINTS
  • 🧪 POTCAR

Next, in order:

  1. Why calculations run as sequences of stages, and what passes between them: 👉 BMD Compute Workflows
  2. The INCARs BMD Compute actually generates for each stage, and why: 👉 BMD Compute INCARs
  3. Spin polarisation, DFT+U and other modifiers, including automatic DFT+U: 👉 BMD Compute Advanced Options

Ready to run it on a supercomputer yourself? 🚀 Head over to 👉 Working on the TAU Supercomputer

Clone this wiki locally