IBRION

IBRION = -1 | 0 | 1 | 2 | 3 | 5 | 6 | 7 | 8 | 44

Default: IBRION	= -1	for NSW=−1 or 0
	= 0	else

Description: IBRION determines how the ions are updated and moved.

For IBRION=0, a molecular dynamics is performed, whereas all other algorithms are destined for relaxations into a local energy minimum. For difficult relaxation problems it is recommended to use the conjugate gradient algorithm (IBRION=2), which presently possesses the most reliable backup routines. Damped molecular dynamics (IBRION=3) are often useful when starting from very bad initial guesses. Close to the local minimum the RMM-DIIS (IBRION=1) is usually the best choice. IBRION=5 and IBRION=6 are using finite differences to determine the second derivatives (Hessian matrix and phonon frequencies), whereas IBRION=7 and IBRION=8 use density functional perturbation theory to calculate the derivatives.

IBRION=-1: no update.

The ions are not moved, but NSW outer loops are performed. In each outer loop the electronic degrees of freedom are re-optimized (for NSW>0 this obviously does not make much sense, except for test purposes). If no ionic update is required use NSW=0 instead.

IBRION=0: molecular dynamics.

Standard ab-initio molecular dynamics. A Verlet algorithm (or fourth-order predictor-corrector if VASP was linked with stepprecor.o) is used to integrate Newton's equations of motion. POTIM supplies the timestep in femto seconds. The parameter SMASS provides additional control.

Mind: At the moment only constant volume MD's are possible.

IBRION=1: ionic relaxation (RMM-DIIS).

For IBRION=1, a quasi-Newton (variable metric) algorithm is used to relax the ions into their instantaneous groundstate. The forces and the stress tensor are used to determine the search directions for finding the equilibrium positions (the total energy is not taken into account). This algorithm is very fast and efficient close to local minima, but fails badly if the initial positions are a bad guess (use IBRION=2 in that case). Since the algorithm builds up an approximation of the Hessian matrix it requires very accurate forces, otherwise it will fail to converge. An efficient way to achieve this is to set NELMIN to a value between 4 and 8 (for simple bulk materials 4 is usually adequate, whereas 8 might be required for complex surfaces where the charge density converges very slowly). This forces a minimum of 4 to 8 electronic steps between each ionic step, and guarantees that the forces are well converged at each step.

The implemented algorithm is called RMM-DIIS.^[1] It implicitly calculates an approximation of the inverse Hessian matrix by taking into account information from previous iterations. On startup, the initial Hessian matrix is diagonal and equal to POTIM. Information from old steps (which can lead to linear dependencies) is automatically removed from the iteration history, if required. The number of vectors kept in the iterations history (which corresponds to the rank of the Hessian matrix must not exceed the degrees of freedom. Naively the number of degrees of freedom is 3(NIONS-1). But symmetry arguments or constraints can reduce this number significantly.

There are two algorithms build in to remove information from the iteration history:

If NFREE is set in the INCAR file, only up to NFREE ionic steps are kept in the iteration history (the rank of the approximate Hessian matrix is not larger than NFREE).
If NFREE is not specified, the criterion whether information is removed from the iteration history is based on the eigenvalue spectrum of the inverse Hessian matrix: if one eigenvalue of the inverse Hessian matrix is larger than 8, information from previous steps is discarded.

For complex problems NFREE can usually be set to a rather large value (i.e. 10-20), however systems of low dimensionality require a carful setting of NFREE (or preferably an exact counting of the number of degrees of freedom). To increase NFREE beyond 20 rarely improves convergence. If NFREE is set to too large, the RMM-DIIS algorithm might diverge.

The choice of a reasonable POTIM is also important and can speed up calculations significantly, we recommend to find an optimal POTIM using IBRION=2 or performing a few test calculations (see below).

IBRION=2: ionic relaxation (conjugate gradient algorithm).

A conjugate-gradient algorithm (a simple discussion of this algorithm can be found for instance in Numerical Recipes, by Press et al.^[2]) is used to relax the ions into their instantaneous groundstate. In the first step ions (and cell shape) are changed along the direction of the steepest descent (i.e. the direction of the calculated forces and stress tensor). The conjugate gradient method requires a line minimization, which is performed in several steps:

First a trial step into the search direction (scaled gradients) is done, with the length of the trial step controlled by the POTIM tag. Then the energy and the forces are recalculated.
The approximate minimum of the total energy is calculated from a cubic (or quadratic) interpolation taking into account the change of the total energy and the change of the forces (3 pieces of information), then a corrector step to the approximate minimum is performed.
After the corrector step the forces and energy are recalculated and it is checked whether the forces contain a significant component parallel to the previous search direction. If this is the case, the line minimization is improved by further corrector steps using a variant of Brent's algorithm.^[2]

To summarize: In the first ionic step the forces are calculated for the initial configuration read from the POSCAR file, the second step is a trial (or predictor step), the third step is a corrector step. If the line minimization is sufficiently accurate in this step, the next trial step is performed.

NSTEP:

initial positions
trial step
corrector step, i.e. positions corresponding to anticipated minimum
trial step
corrector step

...

IBRION=3: ionic relaxation (damped molecular dynamics).

If a damping factor is supplied in the INCAR file by means of the SMASS tag, a damped second order equation of motion is used for the update of the ionic degrees of freedom:

{\displaystyle {\ddot {\vec x}} = -2 \alpha {\vec F} - \mu {\dot {\vec x}}, }

where SMASS supplies the damping factor μ, and POTIM controls α. A simple velocity Verlet algorithm is used to integrate the equation, the discretised equation reads:

{\displaystyle \begin{align} {\vec v_{N+1/2}} =& \Big((1-\mu/2) {\vec v_{N-1/2}} - 2\alpha {\vec F_N} \Big)/(1+\mu/2)\\ {\vec x_{N+1}} =& {\vec x_{N+1}} + {\vec v_{N+1/2}} \end{align} }

One may immediately recognize, that μ=2 is equivalent to a simple steepest descent algorithm (of course without line optimization). Hence, μ=2 corresponds to maximal damping, μ=0 corresponds to no damping. The optimal damping factor depends on the Hessian matrix (matrix of the second derivatives of the energy with respect to the atomic positions). A reasonable first guess for μ is usually 0.4.

Mind that our implementation is particular user-friendly, since changing μ usually does not require to re-adjust the time step POTIM. To choose an optimal time step and damping factor, we recommend the following two step procedure: First fix μ (for instance to 1) and adjust POTIM. POTIM should be chosen as large as possible without getting divergence in the total energy. Then decrease μ and keep POTIM fixed. If POTIM and SMASS are chosen correctly, the damped molecular dynamics mode usually outperforms the conjugate gradient method by a factor of two.

If SMASS is not set in the INCAR file (respectively SMASS<0), a velocity quench algorithm is used. In this case the ionic positions are updated according using the following algorithm: F are the current forces, and α equals POTIM. This equation implies that, if the forces are antiparallel to the velocities, the velocities are quenched to zero. Otherwise the velocities are made parallel to the present forces, and they are increased by an amount that is proportional to the forces.

Mind: For IBRION=3, a reasonable time step must be supplied by the POTIM parameter. Too large time steps will result in divergence, too small ones will slow down the convergence. The stable time step is usually twice the smallest line minimization step in the conjugate gradient algorithm.

IBRION=5 and 6: second derivatives, Hessian matrix and phonon frequencies (finite differences).

A how to for phonon calculations from finite differences is also found here: Phonons from finite differences.

IBRION=5, is available as of VASP.4.5, IBRION=6 starting from VASP.5.1. Both flags allow to determine the Hessian matrix (matrix of the second derivatives of the energy with respect to the atomic positions) and the vibrational frequencies of a system. Only zone centered (Γ-point) frequencies are calculated automatically and printed after

Eigenvectors and eigenvalues of the dynamical matrix

To calculate the Hessian matrix, finite differences are used, i.e. each ion is displaced in the direction of each Cartesian coordinate, and from the forces the Hessian matrix is determined. The two modes differ in the way symmetry is considered. For IBRION=5, all atoms are displaced in all three Cartesian directions, resulting in a significant computational effort even for moderately sized high symmetry systems. For IBRION=6, only symmetry inequivalent displacements are considered, and the remainder of the Hessian matrix is filled using symmetry considerations.

Selective dynamics are presently only supported for IBRION=5; in this case, only those components of the Hessian matrix are calculated for which the selective dynamics tags are set to .TRUE. in the POSCAR file. Contrary to the conventional behavior, the selective dynamics tags now refer to the Cartesian components of the Hessian matrix. For the following POSCAR file, for instance,

Cubic BN
   3.57
 0.0 0.5 0.5
 0.5 0.0 0.5
 0.5 0.5 0.0
   1 1
selective
Direct
 0.00 0.00 0.00  F F F
 0.25 0.25 0.25  T F F

atom 2 is displaced in the x-direction only, and only the x-component of the second atom of the Hessian matrix is calculated.

Three parameters influence the determination of the Hessian matrix: The parameter NFREE determines how many displacements are used for each direction and ion, and POTIM determines the step size. The step size is defaulted to 0.015 Å(starting from VASP.5.1), if too large values are supplied in the input file. Expertise shows that this is a very reasonable compromise.

NFREE=2 uses central differences, i.e., each ion is displaced by a small positive and negative displacement, ±POTIM, along each of the cartesian directions. For NFREE=4, four displacement along each of the cartesian directions are used: ±POTIM and ±2×POTIM.

For NFREE=1, only a single displacement is applied (it is strongly discouraged to use NFREE=1).

Finally, IBRION=6 and ISIF≥3 allows to calculate the elastic constants. The elastic tensor is determined by performing six finite distortions of the lattice and deriving the elastic constants from the strain-stress relationship.^[3] The elastic tensor is calculated both, for rigid ions, as well, as allowing for relaxation of the ions. The elastic moduli for rigid ions are written after the line

SYMMETRIZED ELASTIC MODULI (kBar)

The ionic contributions are determined by inverting the ionic Hessian matrix and multiplying with the internal strain tensor,^[4] and the corresponding contributions are written after the lines:

ELASTIC MODULI CONTR FROM IONIC RELAXATION (kBar)

The final elastic moduli including both, the contributions for distortions with rigid ions and the contributions from the ionic relaxations, are summarized at the very end:

TOTAL ELASTIC MODULI (kBar)

There are a few caveats to this approach: most notably the plane wave cutoff needs to be sufficiently large to converge the stress tensor. This is usually only achieved if the default cutoff is increased by roughly 30%, but it is strongly recommended to increase the cutoff systematically (e.g. in steps of 15%), until full convergence is achieved.

Born effective charges, piezoelectric constants, and the ionic contribution to the dielectric tensor can be calculated additionally by specifying LEPSILON=.TRUE. (linear response theory) or LCALCEPS=.TRUE. (finite external field).

Comments with respect to older releases (pre VASP.5.1): In some older versions, NSW (number of ionic steps) must be set to 1 in the INCAR file, since NSW=0 sets the IBRION=−1 regardless of the value supplied in the INCAR file.

Furthermore, although VASP.4.6 supports IBRION=5-6, VASP.4.6 does not change the set of k-points automatically (often the displaced configurations require a different k-point grid). Hence, the finite difference routine might yield incorrect results in VASP.4.6.

IBRION=7 and 8: second derivatives, Hessian matrix, and phonon frequencies (perturbation theory).

IBRION=7 and IBRION=8 are supported from VASP.5.1 and later versions. It determines the Hessian matrix (matrix of second derivatives with respect to the position of the ions) using density-functional-perturbation theory (DFPT). IBRION=7 does not apply symmetry, whereas IBRION=8 uses symmetry to reduce the number of displacements. The output is similar as for IBRION=5 and 6. Specifically, the second derivates with respect to ionic displacements (interatomic force constants) and the mixed second derivative with respect to the strain and the ionic displacement (internal strain tensor) are evaluated. Although the contributions from the ionic relaxations to the elastic tensor are calculated, the ion-clamped elastic tensor (rigid ion) is not determined. Born effective charges, piezoelectric constants, and the ionic contributions to the dielectric tensor are calculated if LEPSILON=.TRUE. is specified in the INCAR file.

In general, the DFPT routines in VASP are somewhat rudimentary and only support displacements commensurate with the supercell, i.e., so-called q=0 phonons. In other words, VASP can only determine phonon frequencies at the Gamma point of the supercell. Therefore, the code offers few advantages over the finite differences methods discussed above. In particular, the linear response is limited to LDA and GGA functionals and it does not determine the elastic tensors, since the linear response with respect to the strain tensor is not implemented. The only advantage of the linear response routines is that they eliminate the need to choose the magnitude of the finite displacement. Therefore, it might be helpful to first calculate phonon frequencies using linear response and then switch to finite differences and determine the largest displacement that will produce results compatible with the linear response routines.

A few technical comments are in order at this point. VASP solves the linear Sternheimer equation to determine the linear response of the orbitals. Hence, unoccupied orbitals are not required. Internally, the VASP routines for linear response rely on finite differences in two places. (i) The first place is the determination of the second derivative of the exchange-correlation functional: Since most functionals do not support an algebraic determination of second derivatives, VASP always resorts to finite differences to determine the second-order change of the exchange correlation-potential and the PAW one-center terms for each atomic displacement. (ii) Second, after VASP has determined the first-order change of the orbitals, it computes all second derivatives using finite displacements. To do this, VASP displaces the selected atom in the selected directions adds the calculated linear response to the orbitals, and finally determines the differences in the forces and the stress tensor for positive and negative displacements. It can be shown that this yields exactly the second-order force constants and the internal strain tensor, respectively. Furthermore, the Born effective charges are determined "analytically" by contracting the linear response of the orbitals over the "polarization" vector (30) in Ref. ^[4]. These should agree well with the Born effective charges that were previously determined when the linear response with respect to external fields was calculated (there are two different routes to calculate mixed derivatives). The final summary output towards the end of the OUTCAR file. writes the Born effective charges determined from the linear response with respect to external fields.

IBRION=44: the Improved Dimer Method.

This method is described in the Improved Dimer Method section.

Some general comments

For IBRION=1, 2, and 3, the flag ISIF determines whether the ions and/or the cell shape is changed. Update of the cell shape is supported for molecular dynamics (IBRION=0) only if the dynamics module of Tomas Bucko (precompiler flag -Dtbdyn) is used.

Within all relaxation algorithms (IBRION=1, 2, and 3) the parameter POTIM should be supplied in the INCAR file. For IBRION>0, the forces are scaled internally before calling the minimization routine. Therefore for relaxations, POTIM has no physical meaning and serves only as a scaling factor. For many systems, the optimal POTIM is around 0.5. Because the Quasi-Newton algorithm and the damped algorithms are sensitive to the choice of this parameter, use IBRION=2, if you are not sure how large the optimal POTIM is.

In this case, the OUTCAR file and stdout will contain a line indicating a reliable POTIM. For IBRION=2, the following lines will be written to stdout after each corrector step (usually each odd step):

 trial: gam=  .00000 g(F)=   .152E+01 g(S)=  .000E+00 ort = .000E+00
 (trialstep = .82)

The quantity gam is the conjugation parameter to the previous step, g(F) and g(S) are the norm of the force respectively the norm of the stress tensor. The quantity ort is an indicator whether this search direction is orthogonal to the last search direction (for an optimal step this quantity should be much smaller than (g(F) + g(S)). The quantity trialstep is the size of the current trialstep. This value is the average step size leading to a line minimization in the previous ionic step. An optimal POTIM can be determined, by multiplying the current POTIM with the quantity trialstep.

After at the end of a trial step, the following lines are written to stdout:

 trial-energy change:   -1.153185  1.order   -1.133   -1.527  -.739
 step:   1.7275(harm=  2.0557)  dis=  .12277  
                            next Energy= -1341.57 (dE= -.142E+01)

The quantity trial-energy change is the change of the energy in the trial step. The first value after 1.order is the expected energy change calculated from the forces: (F(start)+F(trial))/2×change of positions. The second and third value corresponds to F(start)×change of positions, and F(trial)×change of positions.

The first value in the second line is the size of the step leading to a line minimization along the current search direction. It is calculated from a third order interpolation formula using data from the start and trial step (forces and energy change). harm is the optimal step using a second order (or harmonic) interpolation. Only information on the forces is used for the harmonic interpolation. Close to the minimum both values should be similar. dis is the maximum distance moved by the ions in fractional (direct) coordinates. next Energy gives an indication how large the next energy should be (i.e. the energy at the minimum of the line minimization), dE is the estimated energy change.

The OUTCAR file will contain the following lines, at the end of each trial step:

 trial-energy change:   -1.148928  1.order   -1.126  -1.518  -.735
  (g-gl).g =  .152E+01      g.g   =  .152E+01  gl.gl    =  .000E+00
 g(Force)  =  .152E+01   g(Stress)=  .000E+00 ortho     =  .000E+00
 gamma     =    .00000
 opt step  =   1.72745  (harmonic =   2.05575) max dist = .12277085
 next E    = -1341.577507   (d E  =   1.42496)

The line trial-energy change was already discussed. g(Force) corresponds to g(F), g(Stress) to g(S), ortho to ort, gamma to gam. The values after gamma correspond to the second line (step: ...) previously described.

Related Tags and Sections

NSW, SMASS, POTIM, NFREE, ISIF, LEPSILON, LCALCEPS, Improved Dimer Method

Examples that use this tag

References

↑ P. Pulay, Chem. Phys. Lett. 73, 393 (1980).
↑ ^a ^b W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes (Cambridge University Press, New York, 1986).
↑ Y. Le Page and P. Saxe, Phys. Rev. B 65, 104104 (2002).
↑ ^a ^b X. Wu, D. Vanderbilt, and D. R. Hamann, Phys. Rev. B 72, 035105 (2005).

[pulay:cpl:80-1] P. Pulay, Chem. Phys. Lett. 73, 393 (1980).

[press-2] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes (Cambridge University Press, New York, 1986).

[lepage:prb:02-3] Y. Le Page and P. Saxe, Phys. Rev. B 65, 104104 (2002).

[wu:prb:05-4] X. Wu, D. Vanderbilt, and D. R. Hamann, Phys. Rev. B 72, 035105 (2005).

[1]

[2]

[3]

[4]