From Structure to Properties: Advanced VASP Tasks

Bulk phase calculations allow the crystal to repeat infinitely in three directions; the real surface terminates. The commonly used slab model puts several atomic layers and vacuum into the same periodic box to obtain a repeating system of "sheet-vacuum-sheet". It is not a semi-infinite crystal, and the number of layers, vacuum, and lateral repeating dimensions all contribute to model errors. It is recommended to complete it first Silicon self-consistent and optimization exercises , then enter this section.

Build a consistent family of models

1. Determine the material and crystal planes. First use the same functional to optimize the bulk phase, and then cut out the required (hkl) crystal plane. Record surface termination, surface stoichiometry, lateral supercell, number of layers, and vacuum thickness. Polar surfaces may remodel, change termination, or require compensation that cannot be solved by increasing vacuum alone.
2. Align the normal to the third lattice vector. For starter models, let c be normal to the surface to facilitate the use of 1D vacuum, k-grid, and dipole corrections. Do not indiscriminately mix Cartesian z and tertiary lattice directions in tilted cells.
3. Fixation and relaxation serve the problem. For example, fix several layers at the bottom to simulate bulk support, and relax the top and adsorbents. The fixed number of layers is a model choice that needs to be tested. It cannot be assumed that the "bottom two layers" are true for all surfaces; a symmetrical slab can be fixed with a central layer and terminated identically on both sides.
4. Check periodic adjacent atoms. Adsorbed molecules may extend from boundaries to neighboring cells; it is not sufficient to see that individual boxes do not overlap. Repeat the structure along the in-plane direction to check shortest distances, molecular orientations, and coverage.

POSCAR's Selective dynamics Row enables per-atom degrees of freedom: F F F fixed, T T T Relax. In usual ion optimization, these flags correspond to direct lattice coordinate directions, even if the positions are written in Cartesian. If you want to "fix only the vertical height" in an oblique unit cell, you must first figure out the lattice direction. Official POSCAR definition .

A starting point for neutral, nonmagnetic slab relaxation

The following templates do not contain actual structures, k-grids or POTCAR; ENCUT and center coordinates must be replaced. If you study magnetic transition metals, DFT+U or weak molecular adsorption, the corresponding settings must be verified separately. ISMEAR=0 Is the initial occupation scheme, the metal surface must be tested for smearing and energy difference convergence.

# TEMPLATE: neutral, nonmagnetic slab; normal parallel to c
# Supply POSCAR, converged KPOINTS and licensed POTCAR.
# Replace ENCUT; add validated magnetic/U/vdW settings when needed.
SYSTEM = neutral-slab-relax
GGA = PE
ENCUT = REPLACE_WITH_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-6
NELM = 120
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.05
IBRION = 2
NSW = 150
ISIF = 2
EDIFFG = -0.02
POTIM = 0.3
ISYM = 0
LDIPOL = .TRUE.
IDIPOL = 3
DIPOL = REPLACE_WITH_SLAB_CENTER_IN_FRACTIONAL_COORDINATES
LWAVE = .FALSE.
LCHARG = .FALSE.

ISIF=2 Fixing the entire unit cell avoids treating the vacuum as a volume that should be compressed; the atoms still move according to selective constraints. EDIFFG=-0.02 is an exemplary force threshold of 0.02 eV/Å, which should be tightened further when ultimately comparing small adsorption energy differences. It is often helpful for difficult systems to pre-converge with self-consistent results without dipole potential before enabling correction; this is not a reason to skip the final correction.

For a neutral asymmetric slab with a net dipole moment, IDIPOL=3 specifies the energy correction along the third lattice direction, LDIPOL=.TRUE. Modify potential and force simultaneously; DIPOL Use fractional coordinates to give the reference center close to the slab centroid. Don’t fill in mechanically 0.5 0.5 0.5 , the sheet should be placed in the corresponding position first. Symmetric models should also be checked to see if there really is no net dipole. See IDIPOL, LDIPOL with DIPOL.

When the vacuum is along the c direction and the different lamellae are sufficiently decoupled, the Γ center is often used Nₓ × Nᵧ × 1 The slab k grid; the in-plane Nₓ and Nᵧ need to re-converge with the length of the supercell and the reciprocal lattice vector. Larger in-plane supercells generally allow for a sparser grid, but the same integer grid does not represent the same sampling density. When checking the vacuum platform, use LVHAR=.TRUE. The Hartree potential is added to the output ions and then averaged in the normal plane; this potential does not include the exchange–correlation term. LVHAR.

Define the references behind an adsorption energy

The convention in this tutorial: relative to the selected reference state, Negative values represent exothermic adsorption at the electron energy level . If the literature defines "binding energy" with opposite signs, the definition must be unified before comparison. For H adsorption, ½E(H₂) can be selected as the reference for each H; selecting isolated H atoms will result in different values and cannot be mixed for the same adsorption energy. If n X's are adsorbed in the box, whether to report the total adsorption energy or the average value for each X must also be stated.

1. Separately optimize clean surfaces, multiple adsorption sites/orientations, and independent reference molecules. Gas-phase molecules use a sufficiently large box and Γ point, and confirm the spin state; open-shell molecules such as O₂ cannot be treated in the default non-magnetic state.
2. Do final static calculations with consistent functional, POTCAR data set, ENCUT, in-plane k-point density and energy extraction definitions. Slab maintains the same lateral unit cell and fixed layer strategies as slab+X; the gas phase reference should be box-converged separately.
3. Don't just test the total energy; increase the number of layers, vacuum, lateral area, and k-grid one by one to see if the adsorption energy and key bond lengths meet the research requirements. Vacuum can be tested starting from about 15 Å, which should not be used as a general fit line.
4. For example, compare experimental temperature and pressure, electrode potential, or adsorption free energy, plus consistent zero-point energy, thermal contribution, entropy, solvent, and electrochemical references. Only ΔE is not equal to ΔG, and it does not directly give the reaction rate.

Vacancy, interstitial, and substitution defects alter local bonding, charge trapping, and ion migration. A periodic supercell containing vacancies represents a material with an array of periodic defects and does not automatically equal an infinitely diluted isolated defect. First practice with the simplest neutral vacancies and then discuss charged defects.

Start with a neutral vacancy in silicon

1. Construct supercells of multiple sizes from fully optimized bulk phases. Calculate the complete matrix reference for each size and do not directly splice primitive cell energies of different numerical precisions.
2. Delete a Si atom and update the element count and coordinate row number; the total number of electrons also changes with the number of atoms by default. Give several groups of slight disturbances around the vacancies to prevent the highly symmetrical initial state from locking the structure at a local stationary point.
3. Usually the lattice of the optimized matrix is fixed and only the atoms are relaxed, i.e. ISIF=2 ; This corresponds to the dilute defects being in a bulk lattice environment. If macroscopic volume changes caused by limited concentration or defects are studied, a separate lattice relaxation scheme needs to be defined.
4. Get E with consistent static settings _defect with E _host , while recording local geometry, spin, and whether new occupied states appear, and cannot only save a total energy.

μ_Si is the chemical potential of each Si atom when returned to the reference library; the bulk per atom energy commonly calculated consistently for pure silicon in the zero-temperature approximation. The "+μ" here comes from taking an atom from the supercell. Replacing the energy of isolated atoms changes the physical reference and is not a numerical fix.

Terms in the charged-defect formation energy

here n _i> 0 means that the host has added atoms of species i; removal corresponds to a negative value. q > 0 means that the neutral defect supercell has had electrons removed .E_F Calculated from the parent price band top, E _VBM It needs to be at an electrical potential reference. This equation unifies the potential alignment and finite size processing of the selected method into E _corr ;If another way of writing explicitly adds qΔV, you must check whether it is consistent with E _corr Repeat count.

Different charge states should be optimized separately, and the formation energy changes with E _F The slope of is q, and the intersection point of the lowest envelope gives the thermodynamic charge transition position. Optical transitions are vertical processes discussed under different structural constraints, and the thermodynamic intersection cannot be directly used as a spectral peak. The standard formation energy framework and the dilution limit problem can be obtained from A review of defects by Freysoldt et al. Keep learning.

# TEMPLATE: charged bulk defect at fixed host lattice
# Not a charged-slab template. Supply the actual supercell and references.
# Replace ENCUT and NELECT. Set a justified spin state / MAGMOM if required.
SYSTEM = bulk-defect-relax
GGA = PE
ENCUT = REPLACE_WITH_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-7
NELM = 150
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.05
NELECT = REPLACE_WITH_NEUTRAL_DEFECT_ELECTRONS_MINUS_Q
IBRION = 2
NSW = 200
ISIF = 2
EDIFFG = -0.01
ISYM = 0
LWAVE = .FALSE.
LCHARG = .TRUE.
# No universal correction tag: specify and validate correction externally.

NELECT should be current defective structure and the number of valence electrons of POTCAR to find the neutral N₀, and then take N₀−q ;The number of parent electrons of the undeleted atom cannot be used to determine q. Neutral calculations usually remove NELECT OK, let the program determine by element count. The non-neutral period calculation introduces a compensation background, and the fact that the program can run does not mean that the total charged energy can be directly compared. See NELECT official description .

After optimization, copy the CONTCAR of this charge state to a separate static directory to IBRION=-1; NSW=0 Get ultimate energy. If the correction process requires a potential, it can be set consistently in the Defect and Parent Static tasks. LVHAR=.TRUE. and save LOCPOT; also ensuring that the grid, potential definition, and reference area used for comparison match. Saving CHGCAR alone does not replace the potential data required by the selected correction tool. Official potential output definition .

Three independent sources of error

• chemical potential range: The "A-rich/A-poor" condition of a compound is constrained by the stability of the parent phase and competing phases; all elements cannot be set to the highest chemical potential at the same time. List the conditions of mass conservation and non-precipitation of competing phases first.
• Limited sizes: Increase the supercell, check whether defect charges are localized, and apply corrections appropriate to the unit cell shape and dielectric tensor. FNV and its extensions each have assumptions, and size trends still need to be compared after the software outputs the corrected numbers. See Anisotropy-Relaxation-Related Charge Correction Study .
• Electronic structure: The bandgap and localization errors of PBE affect charge state ordering. After changing the HSE or U value, the relevant references must be recalculated consistently, and it should be explained that the change is the approximation method, not just "drawing the band structure better".

The optimization results in a structure with "little force", which only shows that the first derivative is close to zero. It could be a valley, or it could be a saddle point. Resonance phonon studies the second derivative near the equilibrium position and writes the collective vibration of atoms as the wave vector q and the frequency on the mode branch. It can discuss dynamic stability, vibrational free energy, and heat capacity, but harmonic phonon alone cannot give true thermal conductivity, strong anharmonic phase transitions, or all finite temperature effects.

The finite displacement method finds the force constant from "how the force on each atom changes after being pushed." If the displacement is too small, the force noise will be amplified, and if the displacement is too large, it will deviate from the harmonic range; therefore, electronic convergence, structural residual force, displacement amplitude and supercell size all need to be checked.

Route A: Generate displacements within VASP

IBRION=5 without using symmetry to reduce the displacements, IBRION=6 Take advantage of symmetry; NFREE=2 using central differences with positive and negative displacements, POTIM=0.015 Here it is Å, never the time step of MD. Just ask for the following Enter unit cell Γ point of vibration. If the input is a supercell, the result is the Γ point of the supercell, not the dispersion of the entire Brillouin zone of the primitive cell.

# TEMPLATE: internal finite-displacement Gamma modes of the INPUT cell
# Start from a tightly optimized, nonmagnetic semiconductor/insulator.
# Supply POSCAR / KPOINTS / POTCAR and replace ENCUT.
SYSTEM = gamma-phonon-internal
GGA = PE
ENCUT = REPLACE_WITH_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-8
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.01
IBRION = 6
NFREE = 2
POTIM = 0.015
NSW = 1
ISIF = 2
LWAVE = .FALSE.
LCHARG = .FALSE.

Check eigenmodes and frequencies at OUTCAR, preserving atomic masses and eigenvectors. Do not arbitrarily fix atoms in a crystal that should be fully vibrating; the partial Hessian of the local adsorbate vibration is another approximation, and the built-in finite displacement selective constraint support is different from the direction definition and ordinary optimization. The current official document states that only IBRION=5 This constraint is supported, acting on the Cartesian Hessian component. finite displacement official description , NFREE.

Route B: Phonopy supercells and VASP forces

This route is easy to parallelize: each displacement structure is an independent static task. Phonopy reconstructs multiple sets of forces into force constants, and then performs Fourier transform to obtain the dynamical matrix at any q. The key is to retain the one-to-one correspondence between the displacement file and the output. You cannot reorder atoms, renormalize the unit cell, or optimize the displaced structure again.

# TEMPLATE: one externally displaced Phonopy supercell, SINGLE POINT ONLY
# Do not relax POSCAR-xxx. Use identical settings for every displacement.
SYSTEM = phonopy-force
GGA = PE
ENCUT = REPLACE_WITH_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-8
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.01
IBRION = -1
NSW = 0
ISYM = 0
LWAVE = .FALSE.
LCHARG = .FALSE.

# TEMPLATE: Phonopy 4.5-style workflow. Check your installed --version/--help.
# Need optimized POSCAR-unitcell; this does not run VASP or supply POTCAR.
phonopy-init -d --dim 2 2 2 --pa auto -c POSCAR-unitcell
# Copy each generated POSCAR-xxx to its own disp-xxx/POSCAR.
# Run STATIC VASP forces with the supplied INCAR-force.template there.
# Replace this sample list with ALL files, in phonopy_disp.yaml order:
phonopy-init --sp -f disp-001/vasprun.xml disp-002/vasprun.xml
# --band auto requires seekpath in the SAME active Python environment.
# Check: python -c "import seekpath"
# If missing, follow the official install guide; e.g. conda install -c conda-forge seekpath
# When phonopy_params.yaml has been created successfully:
phonopy --band auto -p phonopy_params.yaml
# Example post-processing q mesh; independently converge it:
phonopy --mesh 20 20 20 -t phonopy_params.yaml
# Older releases: see the manual; setup may use phonopy -d / phonopy -f,
# and postprocessing may use phonopy-load. Do not mix CLI generations.

Phonopy 4.5 documents verified as of this tutorial, used in the preparation and collection stages phonopy-init , post-processing using phonopy ;Old versions often put the preparation phase in phonopy , post-processing may use phonopy-load . First record the installation version and --help , do not splice the two generations of commands. Command migration instructions , VASP interface , Post-processing commands .

--band auto Also required to be available in the same Python environment seekpath ; Installing Phonopy does not mean that every installation method includes this optional dependency. Use first python -c "import seekpath" Check, press if missing Official installation instructions Complete and then generate an automatic path; the self-specified path must also be consistent with the reciprocal basis of the actual primitive cell.

Long-range dipole interactions in polar insulators also influence LO–TO splitting near Γ. The Born effective charge and electron dielectric tensors need to be calculated by the same method, and the BORN data must be provided according to Phonopy's non-analytical term correction requirements; the values cannot be borrowed from another functional or unrelated literature. Phonopy input files and BORN .

How to investigate an imaginary frequency

1. First check the plotting convention: a negative plotted frequency often represents ω² < 0, indicating an imaginary frequency; it does not mean a vibration in a negative direction.
2. Check the structural residual force, stress and magnetic state, increase the ENCUT/k point, tighten the EDIFF, and test the displacement amplitude and supercell. Small acoustic offsets near Γ may contain numerical or rotational invariance errors, but they cannot be deleted altogether.
3. Observe the corresponding eigenvectors, make small perturbations in the positive and negative directions of the soft mode and then optimize the structure to look for possible low-energy distortions; when the wave vector is not Γ, a supercell must be used to accommodate the distortion.
4. If a significant imaginary frequency still exists after convergence, "disharmonic dynamic stability under current approximation and structural conditions" should be reported. If it is claimed that temperature stabilizes the structure, finite temperature/anharmonic evidence is also required.

The out-of-plane acoustic branch of two-dimensional materials is particularly sensitive near Γ, and conclusions cannot be drawn from just a coarse-grid phonon diagram. Official imaginary frequency troubleshooting and soft model operation .

The elastic constant describes the crystal's resistance to small uniform deformations; mathematically, it is the second derivative of energy with respect to strain. It answers "how difficult it is to stretch in a certain direction" and is not the same physical quantity as local vibration of atoms, plastic deformation, and fracture strength. Before calculation, first clarify whether the object is a bulk phase without external stress, a crystal under a specified pressure, or a two-dimensional thin layer containing vacuum.

Two calculation approaches

One method applies several independent, positive and negative small strains to the unit cell, fixes the deformed lattice, calculates the stress or energy and fits it; the other uses VASP's built-in finite difference response. The former is convenient for seeing the linear interval, and the latter is convenient for obtaining the complete tensor. No matter which route, clamped ions and allow internal atomic positions to respond relaxation ions Both elastic constants must be distinguished: the latter are usually closer to low-frequency elastic measurements.

1. Using ENCUT and k-grids appropriate for the stress accuracy, the reference structure was fully optimized to confirm that residual stresses were small relative to the target accuracy. Just because the energy has converged does not mean that the stress must have converged.
2. Built-in routes start with IBRION=6 cooperate ISIF=3 Compute finite difference and internal strain responses. This is to find the derivation near the optimized structure, not to do another ordinary CG structure optimization. IBRION official description .
3. The explicit strain route can select several groups of positive and negative small strain tests, such as 0.3%, 0.5%, and 1%, which are just the starting points for testing. When finding the relaxed ion response in each deformed unit cell, use ISIF=2 Optimize the internal coordinates; you cannot also relax the unit cell back to its original shape.
4. Compare whether the fitting slope is stable under different strain ranges and calculation accuracy; too small strains are affected by noise, and too large strains are mixed with nonlinearity. Don't just take two points without giving uncertainty.

# TEMPLATE: internally relaxed elastic tensor of a 3D BULK solid
# Input must already be tightly relaxed at the target pressure.
# Not suitable for directly reporting a vacuum-containing slab in GPa.
SYSTEM = bulk-elastic-finite-difference
GGA = PE
ENCUT = REPLACE_WITH_STRESS_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-8
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.05
IBRION = 6
ISIF = 3
NFREE = 2
POTIM = 0.015
NSW = 1
LWAVE = .FALSE.
LCHARG = .FALSE.

Three common pitfalls when reading OUTCAR

• kBar is not GPa. TOTAL ELASTIC MODULI (kBar) Multiply by 0.1 to get GPa. Do not label the output number directly as GPa.
• Look at the cutting order of the header. VASP regular list XX YY ZZ XY YZ ZX , with common Voigt xx yy zz yz xz xy Different. When organizing the matrix, both rows and columns need to be mapped synchronously, and then the corresponding contingency conventions are used; you cannot just rename the last three labels.
• There is a convention for stress signs. The diagonal normal stress of VASP corresponds to the compressive tendency. When defining the fit yourself with a "stretched to positive" continuum, you should convert it first, otherwise you may get an overall negative C. ISIF stress definition .

Built-in output may list rigid ion contributions, internal relaxation contributions, and total elastic modulus; indicate which table is selected when reporting. Official output examples are available Discussion on elastic matrix column order . For two-dimensional materials containing vacuum, GPa will change with the height of the artificial box; if the in-plane two-dimensional stiffness is reported, a clear thickness convention should be used or converted to N/m, and the effective height should be stated.

Mechanical and thermodynamic stability are different

These conditions only apply to corresponding symmetry and no external stress conditions. More generally, it should be checked whether the quadratic form of the elastic energy consistent with the adopted strain convention is positive definite; in the presence of finite external pressure, correspondingly corrected stability conditions need to be used. The cubic formula cannot be applied to all crystal systems, nor can we only look at the diagonal element as positive. Mouhat and Coudert's stability criterion .

The energy difference between two stable configurations describes the reaction energy, but not the barrier separating them. NEB places intermediate configurations between known initial and final states, maintains their spacing along the path, and seeks lower energy in directions perpendicular to it. Each image is a complete atomic structure, not an image file. NEB searches for a minimum-energy path near the initial guess, so the starting path and possible alternative mechanisms matter.

Match each directory to its structure

neb-native/
  INCAR       # rename a completed INCAR.template
  KPOINTS     # same converged sampling for all images
  POTCAR      # supplied locally by a licensed user
  00/POSCAR   # relaxed initial state
  01/POSCAR   # intermediate image 1
  02/POSCAR
  03/POSCAR
  04/POSCAR
  05/POSCAR   # intermediate image 5
  06/POSCAR   # relaxed final state

IMAGES=5 Represents five intermediate configurations, so there are seven structure directories in total. Endpoints should be independently fully optimized beforehand and save energy with consistent accuracy; VASP is launched in the parent directory. Different images must have the same crystal lattice, atomic number, element order, atomic identity order, and constraints. It is not enough that the total number of elements is the same: if the number of atom A is exchanged with another equivalent atom B during interpolation, the path may pass through an extremely high-energy overlapping structure.

1. Select the initial and final states to confirm force convergence and local stability respectively. For diffusion, consider the periodic shortest path first; direct interpolation of the 0.99 and 0.01 fractional coordinates might go around the box.
2. Construct a small number of intermediate structures and watch the animation frame by frame, eliminating atomic collisions, overall translations, and unnecessary rotations. Complex reactions could use better geometric interpolation or segmentation mechanisms; adding image is not a way to fix incorrect atomic correspondences.
3. Pre-converge common paths first, then improve accuracy or increase critical interval sampling. When submitting MPI tasks, follow the cluster rules to make the image parallel division reasonable. The total number of processes usually needs to be grouped by the number of images; do not copy the total number of cores of others.
4. Check electron convergence and NEB projection force image by image; reaching the NSW upper limit does not mean the path has converged. You can continue running, but save the results of the previous round and make each intermediate CONTCAR its own new POSCAR.

For catalog conventions and fixed unit cell restrictions, see Native NEB with IMAGES.SPRING=-5 The negative sign of is a switching convention that enables the NEB form and does not represent physically negative stiffness; see SPRING.

Route A: Standard NEB in unmodified VASP

# TEMPLATE: NATIVE VASP NEB, fixed common cell; five intermediate images
# Supply 00/POSCAR through 06/POSCAR, parent KPOINTS and licensed POTCAR.
# Endpoints must be separately relaxed and calculated with matching settings.
SYSTEM = native-neb
GGA = PE
ENCUT = REPLACE_WITH_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-7
NELM = 150
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.05
IMAGES = 5
SPRING = -5
IBRION = 1
NFREE = 2
POTIM = 0.3
ISIF = 2
NSW = 250
EDIFFG = -0.03
ISYM = 0
LWAVE = .FALSE.
LCHARG = .FALSE.
# Native template: NO LCLIMB, NO IOPT, and NOT POTIM=0.

here IBRION=1 with shorter history NFREE=2 is the pre-convergence starting point, POTIM=0.3 It is a finite ion step scale and needs to be adjusted according to path stability. NEB's projection constraints may make conventional IBRION=2 Convergence is difficult and general structural optimization settings cannot be copied by default. You can also research suitable native ones IBRION=3 configuration, but it also requires a non-zero effective step size.

Route B: CI-NEB with a verified VTST installation

CI-NEB drives the highest-energy image uphill along the path toward the saddle point while minimizing its energy in perpendicular directions. It is usually more robust to preconverge an ordinary NEB path and activate the climbing image after the path and peak position are reasonable. Record the VASP and VTST versions and build configuration. Only an extension snippet is provided here; the tutorial does not assume that VTST is installed at your computing center.

# TEMPLATE: VTST-ENABLED VASP ONLY; not a standalone INCAR
# Replace native ionic-optimizer lines with these (do not leave duplicates).
# Retain validated electronic settings and fixed-cell NEB directory structure.
# Confirm the executable was built with a compatible VTST code version.
IBRION = 3
POTIM = 0
IOPT = 7
ICHAIN = 0
IMAGES = 5
SPRING = -5
ISIF = 2
NSW = 250
EDIFFG = -0.03
LCLIMB = .FALSE.
# First preconverge the path. For the separate climbing-image continuation:
# change the existing LCLIMB line to .TRUE. after inspecting the path.

IOPT=7 Refers to VTST's FIRE optimizer. The fragment is explicitly closed first LCLIMB , and then turn it on in the separate continuation stage after confirming the path. Source: VTST NEB Documentation , Optimizer documentation directory , Developer's note on POTIM=0 .

Validating the barrier

Use a clearly defined reaction coordinate, such as cumulative configurational displacement, on the horizontal axis; image number is not time. The highest discrete point in ordinary NEB only approximates the barrier. Refine the saddle point and verify one imaginary mode connecting the initial and final states within the relevant degrees of freedom. Two imaginary modes may indicate that a first-order saddle point has not been reached. A constrained Hessian supports conclusions only within the subspace examined. Finite-temperature free-energy barriers, tunneling, and kinetic prefactors require additional theory; a zero-temperature NEB curve alone does not establish an experimental rate.

Geometric optimization searches for stationary points along the energy decreasing direction; first-principles molecular dynamics substitutes the electronic structure force at each moment into the nuclear motion equation. It can observe local rearrangements, bond breaking and diffusion, but the absence of an event in a short trajectory does not prove that the event will not occur. First clarify the temperature, volume, composition, initial phase, simulation duration and physical quantities to be measured.

Getting started with fixed-cell Langevin NVT

The following examples are based on Exactly two POTCAR atom types Demonstration of 300 K fixed volume trajectories. POTIM=1.0 The unit in MD is fs; high-frequency systems containing H, strong collisions or high-temperature systems may require smaller step sizes and must be verified in practice. The example 2000 steps only corresponds to 2 ps, which is suitable for short debugging and cannot be regarded as a universal production sampling length.

# TEMPLATE: fixed-cell Langevin NVT; EXAMPLE for exactly TWO species
# Supply relaxed POSCAR, tested KPOINTS/POTCAR and electronic settings.
# Replace ENCUT and adapt one LANGEVIN_GAMMA value per POTCAR species.
SYSTEM = nvt-langevin-test
GGA = PE
ENCUT = REPLACE_WITH_FORCE_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-6
NELM = 120
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.05
IBRION = 0
MDALGO = 3
ISIF = 2
ISYM = 0
POTIM = 1.0
NSW = 2000
TEBEG = 300
TEEND = 300
LANGEVIN_GAMMA = 10.0 10.0
LWAVE = .FALSE.
LCHARG = .FALSE.
# No SMASS: it does not set the Langevin thermostat strength.
# 2000 * 1 fs = 2 ps; this is a short test, not a production guarantee.

LANGEVIN_GAMMA The order is consistent with the order of POTCAR types; if the same elements are divided into different types, they should also be listed separately by type. SMASS Used with other thermostats or some ion algorithms, it cannot be used as a replacement for the Langevin friction parameter. MD supports dependencies -Dtbdyn The compile option, currently official templates usually include it, should still be checked on your own executable. See MDALGO, Friction coefficient with SMASS.

Separate equilibration from production sampling

1. Short test: Run a small number of steps first to check for overlapping atoms, temperature anomalies, and depleting NELM of electron cycles at each step. When the input structure, force accuracy, or time steps are substandard, increasing the total time only magnifies the problem.
2. Equilibration: Observe whether the temperature, potential energy, volume (if variable), coordination and structural quantities leave the initial state under the target conditions. It is normal for the temperature to fluctuate around the target value, and a single step equal to exactly 300 K is not an equilibrium criterion.
3. Production section: Discard the balanced segments that have reason to be selected before counting. Block averaging, different time starting points, and independent initial velocities are used to estimate correlation and error, and each adjacent MD frame cannot be treated as an independent sample.
4. Verification step size: Compare critical quantities with smaller time steps and tighter electron convergence; check for energy drift under suitable NVE conditions if necessary. Thermostats exchange energy, and NVT physics cannot always be required to be strictly invariant like an ideal NVE.

Record the correspondence between OUTCAR/OSZICAR, REPORT, XDATCAR and the final state CONTCAR; if you want a continuous trajectory when restarting, you must save the speed and the restart state required by the program, instead of just copying the geometry and re-randomly allocating the speed. Official molecular dynamics workflow , Initial speed and temperature .

Atomic motion alone does not establish a diffusion coefficient

Before calculating the mean square displacement, the trajectory must be unpacked across period boundaries and the overall drift removed as needed; the molecular system should also distinguish between center-of-mass diffusion and internal vibration of the molecule. Only after entering the approximately linear long-term diffusion region and experiencing enough migration events does the slope become statistically significant. Use d=3 for three dimensions, and d=2 for strict two-dimensional diffusion; you cannot give an anisotropic system just a number without specifying the direction.

There are two categories of computational parameters: one determines numerical accuracy, such as ENCUT, k-point, and EDIFF; the other defines approximate physical models, such as exchange–correlation functional, U, dispersion, and spin–orbit coupling. Improving the former category does not automatically correct systematic biases in the latter category. When doing research, you should first write down "what physics you want to improve" and then decide on the method.

Dispersion in molecules, layered materials, and weak adsorption

Semi-local functional often cannot adequately describe long-range dispersion interactions. Issues such as interlayer distance, molecular orientation, weak adsorption energies, etc. can be compared with a consistent choice of dispersion correction or non-local vdW functional. D3(BJ) is a geometrically dependent dispersion correction. This tutorial is supported by PBE. IVDW=12 demonstration; demonstration; IVDW=11 It corresponds to another damping form, and the values cannot be interchanged as "accuracy levels" at will.

# TEMPLATE OVERLAY ONLY: PBE with D3(BJ) dispersion
# Add to ONE otherwise complete, validated PBE INCAR.
# Do not stack onto another dispersion-corrected functional uncritically.
GGA = PE
IVDW = 12
# Apply the same method to all comparable structures / references.

There are two things to check when using it: whether the actual output recognizes the expected dispersion scheme, and whether the functional has a supported damping parameter. Don't adjust a set of empirical coefficients individually to be close to the experiment but still call it standard PBE-D3(BJ). If you change to another functional, or use a new solution that requires an external library, you must recheck the version and parameter support. D3 official method and parameter description .

1. Compare PBE to selected vdW methods for the same system on structures such as layer spacing or molecule-to-surface height; then look at the change in the target energy difference rather than just comparing two total energy values.
2. If the vdW method is used in the final report, the structure corresponding to the problem should be obtained under this method or the "single point correction on the PBE geometry" should be explicitly stated. Geometric response and energy correction are not the same thing.
3. When calculating adsorption or phase change energy, the same method must be used for each term; you cannot just add dispersion to the adsorbed state and use different definitions for the gas phase/clean surface.
4. Do not unconditionally superimpose D3 on functionals that already contain non-local dispersion, otherwise they may be counted repeatedly; the applicability of U, magnetic or surface models must also be independently verified.

HSE06: An Orbital-dependent Screened Hybrid Functional

HSE06 replaces a fraction of semilocal exchange with short-range exact exchange. Its standard definition uses PBE with AEXX=0.25 with HFSCREEN=0.2 Å⁻¹; write only LHFCALC=.TRUE. It does not mean that HSE06 has been selected. Hybrid functional is often used to examine semiconductor band gaps, local electrons and certain defects, but "more expensive" does not mean more accurate for every energy, metal system or experimental condition. HFSCREEN.

# TEMPLATE: HSE06 static self-consistent calculation from scratch
# For a nonmagnetic semiconductor; regular k mesh REQUIRED.
# Supply POSCAR, KPOINTS and suitable licensed PAW-PBE POTCAR.
SYSTEM = hse06-static
GGA = PE
LHFCALC = .TRUE.
AEXX = 0.25
HFSCREEN = 0.2
ENCUT = REPLACE_WITH_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-7
NELM = 150
ALGO = Damped
TIME = 0.4
LREAL = .FALSE.
LASPH = .TRUE.
ISMEAR = 0
SIGMA = 0.05
ISTART = 0
ICHARG = 2
IBRION = -1
NSW = 0
LWAVE = .TRUE.
LCHARG = .TRUE.
# PBE WAVECAR restart is possible with compatible settings and ISTART=1;
# review ICHARG/restart settings then. NEVER set ICHARG=11 for this task.

Template explicitly enabled LASPH=.TRUE. , containing the non-spherical component of the density gradient within the PAW sphere; this is especially important for mixed functional systems containing d/f electrons, transition metal oxides, etc. PBE calculations used for method comparisons should also consistently enable this setting to avoid setting differences being mixed into functional differences. LASPH official description .

ALGO=Damped with TIME Controls electron orbital optimization; TIME is not equal to POTIM for ion optimization or MD. A more conservative TIME can be tested when oscillation occurs. It is officially not recommended to use it for mixed functional ALGO=Fast ; Specific optimizer and version performance should be tested with small tasks. Electronic optimization suggestions for mixed functional .

Why PBE CHGCAR + ICHARG=11 is insufficient for HSE bands

A semilocal Hamiltonian can be determined from the density. A hybrid functional also needs occupied orbitals on a regular k-point mesh to construct nonlocal exchange. Fixing the PBE charge density and solving only along high-symmetry lines omits this information. HSE must therefore be solved self-consistently as an orbital-dependent method. Retain a compatible WAVECAR for restarts. PBE wavefunctions can provide an initial guess, but they are not the final HSE wavefunctions.

The band structure can be used as "weighted regular grid + zero-weighted path points", or using the supported version of KPOINTS_OPT to have the regular grid defined separately from the paths. Neither of them can delete the regular grid, nor can they set it ICHARG=11 . The convergence check and memory cost of path points are also different from ordinary SCF, read first Official HSE band structure process , then use Electronic Structure Chapter drawing principles.

For the HSE band structure of the bandgap system, the official recommendation is to study HFRCUT=-1 Exchange Coulomb truncation to avoid some path discontinuities. The orbital of the explicit zero-weight point does not directly determine the SCF total energy criterion and may stop prematurely when restarting from a converged HSE wavefunction; iterative convergence of the path band structure should be independently checked and the effect of increasing NELMIN verified if necessary. This setting belongs to the special process of band structure, and "always converged" cannot be used to replace the path track check.

The dielectric function describes the response of a material to an electric field perturbation and is typically a frequency- and direction-dependent complex tensor. At very high frequencies, the atomic nuclei have no time to move; at lower frequencies, lattice vibrations can also contribute to polarization. Therefore, the "dielectric constant" must at least indicate the electron contribution or the electron plus ion contribution, static or dynamic, which crystal orientation, and whether it includes local field and exciton effects.

A workflow for independent-particle optical spectra

1. Establish electronic ground states with explicit geometry and functionality. Optical integration uses dense regular k grids, not the highly symmetric paths used in band structure diagrams.
2. Add unoccupied bands and see the range of transition energies covered by the highest band. NBANDS Including occupied and empty states, the default empty state margin of the program is often not an optically convergent value. Increase it step by step until the spectrum shape of the target energy region is stable. Don't just look at whether the total energy changes.
3. Encrypt k-points, increase frequency sampling, change smearing and compare respectively. A smoother spectrum may just mean greater smearing, but does not mean that the microscopic calculations are more accurate.
4. Extracts real and imaginary tensors in OUTCAR or a supported output, indicating components, energy units, and smearing. Each curve must be traceable to the same set of inputs.

# TEMPLATE: independent-particle PBE optics of a nonmagnetic 3D insulator
# Supply optimized POSCAR, DENSE regular KPOINTS, POTCAR; replace two fields.
# This is a self-consistent static calculation, not a band-line calculation.
SYSTEM = pbe-independent-particle-optics
GGA = PE
ENCUT = REPLACE_WITH_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-8
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.01
IBRION = -1
NSW = 0
LOPTICS = .TRUE.
NBANDS = REPLACE_WITH_CONVERGED_NUMBER_OF_BANDS
NEDOS = 2000
CSHIFT = 0.10
LWAVE = .TRUE.
LCHARG = .TRUE.

LOPTICS The frequency response is calculated after the electronic ground state, and the standard independent particle approximately ignores the corresponding local field changes; the spectral peak cannot be directly called the experimental exciton peak. It can be combined with a hybrid functional, but still requires the hybrid function to be self-consistent and the applicable optimizer. The electronic band gap, optical absorption onset, and experimental thresholds containing phonons/excitons may differ. Approximate and applicable conditions for LOPTICS .

CSHIFT=0.10 The unit is eV, used for Lorentz smearing; at the same time, there is also SIGMA effect when using Gaussian occupation, so the template makes SIGMA significantly smaller than CSHIFT, and the spectrum tail and frequency range still need to be checked. NEDOS=2000 Just an example of the output frequency samples, not the k-point count. See CSHIFT, NEDOS with NBANDS.

Static electronic response and Born effective charges

# TEMPLATE: electronic (ion-clamped) static dielectric tensor / Born charges
# For a nonmetallic system with supported semilocal DFT, NOT HSE.
# Supply optimized POSCAR, converged KPOINTS/POTCAR, replace ENCUT.
SYSTEM = pbe-dielectric-dfpt
GGA = PE
ENCUT = REPLACE_WITH_CONVERGED_ENCUT
PREC = Accurate
EDIFF = 1E-8
ALGO = Normal
LREAL = .FALSE.
ISMEAR = 0
SIGMA = 0.01
IBRION = -1
NSW = 0
LEPSILON = .TRUE.
LRPA = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.
# Not the full ionic+electronic static permittivity of a relaxed crystal.
# Do not append LOPTICS or excessive NBANDS to this separate DFPT task.

The Born effective charge is the response tensor between atomic displacements and macroscopic polarization changes and is not equal to the Bader net charge, nor is it necessarily equal to the formal oxidation state. LEPSILON The DFPT route does not require the massive addition of unoccupied bands like LOPTICS does, LRPA=.FALSE. Correspondence preserves the settings for proper exchange–correlation responsiveness; it is clearer to do the two tasks separately. LEPSILON.

When a complete static dielectric response is required, the applicable vibrational and ionic response calculations are further combined and the electronic ε∞ and ionic contributions are reported separately. Track-dependent methods such as HSE can study finite electric field routes, but the field strength needs to be tested in the linear region and the insulation must be maintained. The official current description is that LEPSILON is not suitable for mixed functional; refer to LCALCEPS.

The most common failure in advanced computing is not a program crash, but the completion of a set of calculations that do not match the scientific problem. Write a goal before starting each task, such as "Compare the electron adsorption energies of two adsorption sites under the same coverage", and then decide how strict the energy difference, force, size and temperature sampling standards are needed based on this.

Reduce a failure to a diagnosable problem

When debugging, save the original input and logs, and only change a set of settings at a time in the new directory that makes sense. Errors can be localized with cheaper setups, but the final physical results must be recovered and verified to the required accuracy. Do not pile INCAR with a dozen tags from unknown sources at the same time, making it impossible to attribute the problem. Official Electronic Convergence Troubleshooting .

Benchmark a small calculation before estimating total cost

• phonon: The total cost depends on the number of unequal displacements multiplied by the static stress cost of each supercell; each displacement can be submitted independently, but file ordering errors must be avoided.
• NEB: The number of intermediate structures, the number of most difficult electron cycles and ion steps per frame will affect the waiting time. Parallel images do not make the total computational resource consumption disappear.
• AIMD: Use test runs to measure average cost per step, multiply by the number of production steps, and reserve resources for thermal balancing, repeat trajectories, and failed restarts.
• HSE/Optics: More orbitals, k-points, and swap calculations may significantly increase memory; doubling the number of cores does not guarantee that the time will be halved, and may even be slower due to communication.

Use this machine model, this program version and approximately the same system to test in parallel. NCORE, KPAR , image grouping and GPU/OpenMP limitations are related to each other, and the best parameters on a certain CPU node cannot be used as a universal template. Check dry-run first or ALGO=None Yes, check actual bands, k-points, memory and short timing. Official parallel optimization process .

Where can AI accelerate the work?

Machine learning potential learns potential energy surfaces from DFT energies, forces, and stresses to expand dynamic time and space scales within proven configurations; active learning prioritizes computation of local environments unfamiliar to the model. VASP itself provides an MLFF route based on online learning. The specific tags depend on the version and mode. You should first learn to evaluate data quality when getting started. Official MLFF Overview .

A reasonable validation set should cover independent temperatures, structural perturbations, defects, or reaction environments, rather than randomly dividing adjacent frames of the same trajectory into training and testing. Low mean force errors do not guarantee that migration barriers, rare events, and new phases are also correct; key candidate structures and saddle points should still be returned for the same baseline DFT review. The model inherits the bias of the training function and cannot automatically turn it into the experimental true value based on the amount of data.

Large language models are suitable for interpreting output, helping generate inspection scripts and organizing documentation. It may also make up non-existent parameters or treat VTST as a native function, so before turning automatic suggestions into calculations, you should check the official definitions, units, compatibility, and check the output that actually works. The research reports are respectively marked "AI prediction", "DFT review", "limited temperature simulation" and "experimental verification".

What to include in a calculation record

1. Questions and models: Structural origin, unit cell, composition, charge, spin, surface termination, pinned layers and reference states.
2. Software and inputs: VASP/external tool version, build or extension information, POTCAR identification and verification information, functional, ENCUT, k-point, smearing, optimization and response parameters.
3. Evidence: Electron/ion convergence, supercell and numerical accuracy tests, output files used, post-processing scripts and unit conversions.
4. Scope of the conclusions: Which are definite numerical trends, which are predictions of current models, and which still require larger models, more sampling, or experimental testing.

Don’t just keep the screenshots and final form. Preserve original input, software output, structure, and analysis procedures; do not distribute license-restricted POTCAR files when making reproducible material public. return to Course Map Choose an exercise, or Parameters and Checklist Check item by item.

Read the file description and corresponding chapters first, and then complete the model, POTCAR and pre-calculation. Files with the suffix .template contain parameters to be changed and cannot be submitted as they are.

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.

1. Official POSCAR definition
2. IDIPOL
3. LDIPOL
4. DIPOL
5. LVHAR
6. A review of defects by Freysoldt et al.
7. NELECT official description
8. Anisotropy-Relaxation-Related Charge Correction Study
9. Official charging system discussion
10. finite displacement official description
11. NFREE
12. Command migration instructions
13. VASP interface
14. Post-processing commands
15. Phonopy input files and BORN
16. Official imaginary frequency troubleshooting and soft model operation
17. IBRION official description
18. ISIF stress definition
19. Discussion on elastic matrix column order
20. Mouhat and Coudert's stability criterion
21. Native NEB
22. IMAGES
23. SPRING
24. Official troubleshooting instructions
25. VTST NEB Documentation
26. Optimizer documentation directory
27. Developer's note on POTIM=0
28. MDALGO
29. Friction coefficient
30. SMASS
31. Official molecular dynamics workflow
32. Initial speed and temperature
33. D3 official method and parameter description
34. HFSCREEN
35. Electronic optimization suggestions for mixed functional
36. Official HSE band structure process
37. Approximate and applicable conditions for LOPTICS
38. CSHIFT
39. NEDOS
40. NBANDS
41. LEPSILON
42. LCALCEPS
43. Official Electronic Convergence Troubleshooting
44. Official parallel optimization process
45. Official MLFF Overview
46. Phonopy installation and seekpath dependencies
47. LASPH: PAW non-spherical component