IBRION: Difference between revisions

Latest revision as of 12:53, 18 October 2024

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

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

Description: determines how the crystal structure changes during the calculation:

no update
- IBRION=-1 (Avoid setting IBRION=-1 with NSW>0 to prevent recomputing the same structure NSW times).

Molecular dynamics
- IBRION=0

Structure optimization
- IBRION=1 RMM-DIIS
- IBRION=2 conjugate gradient
- IBRION=3 damped molecular dynamics

Computing phonon modes
- IBRION=5 finite differences without symmetry
- IBRION=6 finite differences with symmetry
- IBRION=7 perturbation theory without symmetry
- IBRION=8 perturbation theory with symmetry

Analyzing transition states
- IBRION=40 intrinsic-reaction-coordinate calculations
- IBRION=44 improved dimer method

User-supplied interactive changes
- IBRION=11 from standard input
- IBRION=12 from Python plugin

Molecular dynamics

In molecular-dynamics (MD) simulations the positions of the ions are updated using a classical equation of motion for the ions. There are several algorithms for the time propagation in MD controlled by selecting MDALGO and the choice of the thermostats. The MD run performs NSW timesteps of length POTIM.

Frequently, performing an ab-initio calculations in every step of an MD simulation is too expensive so that machine-learned force fields are needed.

Tip: In order to limit the output of the MD simulation, control the verbosity by setting NWRITE=0,1, or reduce the frequency of output using ML_OUTBLOCK, NBLOCK, or KBLOCK.

Structure optimization

VASP optimizes the structure based on the degrees of freedom selected with the ISIF tag and (if used) the selective dynamics POSCAR file. Generally, the larger the number of degrees of freedom, the harder it is to find the optimal solution. To find the solution, VASP provides multiple algorithms:

RMM-DIIS (IBRION=1) reduces the forces by linear combination of previous positions. It is the method of choice for larger systems (>20 degrees of freedom) that are reasonably close to the ground-state structure.
Conjugate gradient (IBRION=2) finds the optimal step size along a search direction. It is a robust default choice but may need more iterations than RMM-DIIS.
Damped molecular dynamics (IBRION=3) runs a MD simulation with decreasing velocity of the ions. Use this for large systems far away from the minimum to get to a better starting point for the other algorithms.

Consult the structure optimization page for advise on how to choose the optimization algorithm.

Computing the phonon modes

The second-order derivatives of the total energy $E$ with respect to ionic positions $R_{\alpha i}$ of ion $\alpha$ in the direction $i$ , is computed using a first-order derivative of the forces $F_{\beta j}$ . Then, the dynamical matrix $D_{\alpha i \beta j}$ is constructed, diagonalized, and the phonon modes and frequencies of the system are reported in the OUTCAR file and vaspout.h5. Also see theory on phonons.

Tip: It may be necessary to set EDIFF because the default (EDIFF = 1E-4) often results in unacceptably large errors.

VASP implements two different methods to compute the phonon modes and can use symmetry to reduce the number of computed displacements:

IBRION = 5 finite differences without symmetry
IBRION = 6 finite differences with symmetry
IBRION = 7 density-functional-perturbation theory without symmetry
IBRION = 8 density-functional-perturbation theory with symmetry

For finite differences, the elastic tensors and internal strain tensors is computed for ISIF>=3. Compute Born-effective charges, piezoelectric constants, and the ionic contribution to the dielectric tensor by specifying LEPSILON = .TRUE. (linear response theory) or LCALCEPS = .TRUE. (finite external field).

Also see computing the phonon dispersion and DOS.

Analyzing transition states

To study the kinetics of chemical reactions, one may want to construct transition states or follow the reaction path. For the analysis of transition states the following methods are available:

Setting IBRION = 40, you can start from a transition state and monitor the energy along an intrinsic-reaction coordinate (IRC). The IRC calculations section describes this method.
With the improved dimer method (IBRION = 44), you can search for a the transition state starting from an arbitrary structure in the investigated phase space.
The nudged elastic bands method finds an approximate reaction path based on the initial and final structure, i.e., reactant and product.

Interactively supplied positions and lattice vectors

Occasionally, you may want to run VASP for related structures where the overhead of restarting VASP is significant. In these scenarios, VASP provides the following alternatives

With IBRION = 11, you can provide new structures via the standard input. For ISIF>=3, a complete POSCAR file is read, otherwise just the positions in fractional coordinates.

@@ Line 1: / Line 1: @@
-{{TAGDEF|IBRION|-1 {{!}} 0 {{!}} 1 {{!}} 2 {{!}} 3 {{!}} 5 {{!}} 6 {{!}} 7 {{!}} 8 {{!}} 44}}
+{{TAGDEF|IBRION|-1 {{!}} 0 {{!}} 1 {{!}} 2 {{!}} 3 {{!}} 5 {{!}} 6 {{!}} 7 {{!}} 8 {{!}} 11 {{!}} 12 {{!}} 40 {{!}} 44}}
 {{DEF|IBRION|-1|for {{TAG|NSW}}{{=}}&minus;1 or 0|0|else}}
-Description: {{TAG|IBRION}} determines how the ions are updated and moved.
+Description: determines how the crystal structure changes during the calculation:
-----
+:::* no update
-For {{TAG|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 ({{TAG|IBRION}}=2), which presently possesses the most reliable backup routines. Damped molecular dynamics ({{TAG|IBRION}}=3) are often useful when starting from very bad initial guesses. Close to the local minimum the RMM-DIIS ({{TAG|IBRION}}=1) is usually the best choice. {{TAG|IBRION}}=5 and {{TAG|IBRION}}=6 are using finite differences to determine the second derivatives (Hessian matrix and phonon frequencies), whereas {{TAG|IBRION}}=7 and {{TAG|IBRION}}=8 use density functional perturbation theory to calculate the derivatives.
+:::** {{TAG|IBRION}}=-1 (Avoid setting {{TAG|IBRION}}=-1 with {{TAG|NSW}}>0 to prevent recomputing the same structure {{TAG|NSW}} times).
+:::* [[#Molecular dynamics|Molecular dynamics]]
+:::** {{TAG|IBRION}}=0
-== {{TAG|IBRION}}=-1: no update.==
+:::* [[#Structure optimization|Structure optimization]]
-:The ions are not moved, but {{TAG|NSW}} outer loops are performed. In each outer loop the electronic degrees of freedom are re-optimized (for {{TAG|NSW}}>0 this obviously does not make much sense, except for test purposes). If no ionic update is required use {{TAG|NSW}}=0 instead.
+:::** {{TAG|IBRION}}=1 RMM-DIIS
+:::** {{TAG|IBRION}}=2 conjugate gradient
+:::** {{TAG|IBRION}}=3 damped molecular dynamics
-== {{TAG|IBRION}}=0: molecular dynamics.==
+:::* [[#Computing the phonon modes|Computing phonon modes]]
-: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. {{TAG|POTIM}} supplies the timestep in femto seconds. The parameter {{TAG|SMASS}} provides additional control.
+:::** {{TAG|IBRION}}=5 finite differences without symmetry
+:::** {{TAG|IBRION}}=6 finite differences with symmetry
+:::** {{TAG|IBRION}}=7 perturbation theory without symmetry
+:::** {{TAG|IBRION}}=8 perturbation theory with symmetry
-:'''Mind''': At the moment only constant volume MD's are possible.
+:::* [[#Analyzing transition states|Analyzing transition states]]
+:::** {{TAG|IBRION}}=40 [[IRC calculations|intrinsic-reaction-coordinate calculations]]
+:::** {{TAG|IBRION}}=44 [[improved dimer method]]
-== {{TAG|IBRION}}=1: ionic relaxation (RMM-DIIS).==
+:::* [[#User-supplied interactive changes|User-supplied interactive changes]]
-:For {{TAG|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 {{TAG|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 {{TAG|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.
+:::** {{TAG|IBRION}}=11 from standard input
+:::** {{TAG|IBRION}}=12 from Python plugin
-:The implemented algorithm is called RMM-DIIS.<ref name="pulay:cpl:80"/> 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 {{TAG|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:
+== Molecular dynamics ==
-:#If {{TAG|NFREE}} is set in the {{FILE|INCAR}} file, only up to {{TAG|NFREE}} ionic steps are kept in the iteration history (the rank of the approximate Hessian matrix is not larger than {{TAG|NFREE}}).
-:#If {{TAG|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 {{TAG|NFREE}} can usually be set to a rather large value (i.e. 10-20), however systems of low dimensionality require a carful setting of {{TAG|NFREE}} (or preferably an exact counting of the number of degrees of freedom). To increase {{TAG|NFREE}} beyond 20 rarely improves convergence. If {{TAG|NFREE}} is set to too large, the RMM-DIIS algorithm might diverge.
+In [[:Category:Molecular dynamics|molecular-dynamics (MD) simulations]] the positions of the ions are updated using a classical equation of motion for the ions. There are several algorithms for the [[time-propagation algorithms in molecular dynamics|time propagation in MD]] controlled by selecting {{TAG|MDALGO}} and the choice of the [[thermostats]]. The MD run performs {{TAG|NSW}} timesteps of length {{TAG|POTIM}}.
-:The choice of a reasonable {{TAG|POTIM}} is also important and can speed up calculations significantly, we recommend to find an optimal {{TAG|POTIM}} using {{TAG|IBRION}}=2 or performing a few test calculations (see below).
+Frequently, performing an [[electronic minimization|ab-initio calculations]] in every step of an MD simulation is too expensive so that [[:Category:Machine-learned_force_fields|machine-learned force fields]] are needed.
+{{NB|tip|In order to limit the output of the MD simulation, control the verbosity by setting {{TAG|NWRITE}}{{=}}0,1, or reduce the frequency of output using {{TAG|ML_OUTBLOCK}}, {{TAG|NBLOCK}}, or {{TAG|KBLOCK}}.}}
-== {{TAG|IBRION}}=2: ionic relaxation (conjugate gradient algorithm).==
+== Structure optimization ==
-:A conjugate-gradient algorithm (a simple discussion of this algorithm can be found for instance in ''Numerical Recipes'', by Press ''et al.''<ref name="press"/>) 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 {{TAG|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.<ref name="press"/>
-:To summarize: In the first ionic step the forces are calculated for the initial configuration read from the {{FILE|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.
+VASP optimizes the structure based on the degrees of freedom selected with the {{TAG|ISIF}} tag and (if used) the selective dynamics {{FILE|POSCAR}} file.
+Generally, the larger the number of degrees of freedom, the harder it is to find the optimal solution.
+To find the solution, VASP provides multiple algorithms:
-:NSTEP:
+* RMM-DIIS ({{TAG|IBRION}}=1) reduces the forces by linear combination of previous positions. It is the method of choice for larger systems (>20 degrees of freedom) that are reasonably close to the ground-state structure.
-:#	initial positions
+* Conjugate gradient ({{TAG|IBRION}}=2) finds the optimal step size along a search direction. It is a robust default choice but may need more iterations than RMM-DIIS.
-:#	trial step
+* Damped molecular dynamics ({{TAG|IBRION}}=3) runs a MD simulation with decreasing velocity of the ions. Use this for large systems far away from the minimum to get to a better starting point for the other algorithms.
-:#	corrector step, i.e. positions corresponding to anticipated minimum
-:#	trial step
-:#	corrector step
-::...
-== {{TAG|IBRION}}=3: ionic relaxation (damped molecular dynamics).==
-== {{TAG|IBRION}}=5 and 6: second derivatives, Hessian matrix and phonon frequencies (finite differences).==
-:{{TAG|IBRION}}=5, is available as of VASP.4.5, {{TAG|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 (&Gamma;-point) frequencies are calculated automatically and printed after
- Eigenvectors and eigenvalues of the dynamical matrix
+Consult the [[structure optimization]] page for advise on how to choose the optimization algorithm.
-: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 {{TAG|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 {{TAG|IBRION}}=6, only symmetry inequivalent displacements are considered, and the remainder of the Hessian matrix is filled using symmetry considerations.
+== Computing the phonon modes ==
-:Selective dynamics are presently only supported for {{TAG|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 {{FILE|POSCAR}} file. Contrary to the conventional behavior, the selective dynamics tags now refer to the Cartesian components of the Hessian matrix. For the following {{FILE|POSCAR}} file, for instance,
+The second-order derivatives of the total energy <math>E</math> with respect to ionic positions <math>R_{\alpha i}</math> of ion <math>\alpha</math> in the direction <math>i</math>, is computed using a first-order derivative of the [[forces]] <math>F_{\beta j}</math>. Then, the dynamical matrix <math>D_{\alpha i \beta j}</math> is constructed, diagonalized, and the phonon modes and frequencies of the system are reported in the {{FILE|OUTCAR}} file and {{FILE|vaspout.h5}}. Also see [[Phonons: Theory|theory on phonons]].
+{{NB|tip|It may be necessary to set {{TAGO|EDIFF|<= 1E-6}} because the default ({{TAGO|EDIFF|1E-4}}) often results in unacceptably large errors.}}
+VASP implements two different methods to compute the phonon modes and can use symmetry to reduce the number of computed displacements:
- Cubic BN
+* {{TAGO|IBRION|5}} [[Phonons from finite differences|finite differences]] '''without''' symmetry
-.57
+* {{TAGO|IBRION|6}} [[Phonons from finite differences|finite differences]] '''with''' symmetry
-.0 0.5 0.5
+* {{TAGO|IBRION|7}} [[Phonons_from_density-functional-perturbation_theory|density-functional-perturbation theory]] '''without''' symmetry
-.5 0.0 0.5
+* {{TAGO|IBRION|8}} [[Phonons_from_density-functional-perturbation_theory|density-functional-perturbation theory]] '''with''' symmetry
-.5 0.5 0.0
-1
- selective
- Direct
-.00 0.00 0.00  F F F
-.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.
+For finite differences, the elastic tensors and internal strain tensors is computed for {{TAG|ISIF}}>=3.
+Compute Born-effective charges, piezoelectric constants, and the ionic contribution to the dielectric tensor by specifying {{TAGO|LEPSILON|.TRUE.}} ([[Linear response|linear response theory]]) or {{TAGO|LCALCEPS|.TRUE.}} (finite external field).
-:Three parameters influence the determination of the Hessian matrix: The parameter {{TAG|NFREE}} determines how many displacements are used for each direction and ion, and {{TAG|POTIM}} determines the step size. The step size is defaulted to 0.015 &Aring;(starting from VASP.5.1), if too large values are supplied in the input file. Expertise shows that this is a very reasonable compromise.
+Also see [[computing the phonon dispersion and DOS]].
-:{{TAG|NFREE}}=2 uses central differences, ''i.e.'', each ion is displaced by a small positive and negative displacement, &plusmn;{{TAG|POTIM}}, along each of the cartesian directions. For {{TAG|NFREE}}=4, four displacement along each of the cartesian directions are used: &plusmn;{{TAG|POTIM}} and &plusmn;2&times;{{TAG|POTIM}}.
-:For {{TAG|NFREE}}=1, only a single displacement is applied (it is strongly discouraged to use {{TAG|NFREE}}=1).
+== Analyzing transition states ==
-:Finally, {{TAG|IBRION}}=6 and {{TAG|ISIF}}&ge;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 [4]. 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)
+To study the kinetics of chemical reactions, one may want to construct [[transition states]] or follow the reaction path.
+For the analysis of transition states the following methods are available:
-:The ionic contributions are determined by inverting the ionic Hessian matrix and multiplying with the internal strain tensor [5], and the corresponding contributions are written after the lines:
+* Setting {{TAGO|IBRION|40}}, you can start from a transition state and monitor the energy along an intrinsic-reaction coordinate (IRC). The [[IRC calculations]] section describes this method.
+* With the [[improved dimer method]] ({{TAGO|IBRION|44}}), you can search for a the transition state starting from an arbitrary structure in the investigated phase space.
+* The [[nudged elastic bands]] method finds an approximate reaction path based on the initial and final structure, i.e., reactant and product.
- ELASTIC MODULI CONTR FROM IONIC RELAXATION (kBar)
+== Interactively supplied positions and lattice vectors ==
-: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:
+Occasionally, you may want to run VASP for related structures where the overhead of restarting VASP is significant.
+In these scenarios, VASP provides the following alternatives
- TOTAL ELASTIC MODULI (kBar)
+* With {{TAGO|IBRION|11}}, you can provide new structures via the standard input. For {{TAG|ISIF}}>=3, a complete {{FILE|POSCAR}} file is read, otherwise just the positions in fractional coordinates.
-: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.
+<!--
+* If you [[Makefile.include#Plugins_(optional)|linked VASP with Python]], you can [[Writing a Python plugin|write a Python plugin]] to modify the structure. Set {{TAGO|IBRION|12}} or {{TAGO|PLUGINS/STRUCTURE|T}} to activate it.
+-->
-:'''Mind''': In some older versions, {{TAG|NSW}} (number of ionic steps) must be set to 1 in the {{FILE|INCAR}} file, since {{TAG|NSW}}=0 sets the {{TAG|IBRION}}=&minus;1 regardless of the value supplied in the {{FILE|INCAR}} file.
+== Related tags and articles ==
-:A final problem concerns the symmetry treatment in VASP.4.6: it determines the symmetry for the displaced configurations correctly, but unfortunately VASP.4.6 does not change the set of '''k'''-points automatically (often the lower symmetry of configurations with displaced ions would require one to use more '''k'''-points). Hence, for accurate calculations, the symmetry must be switched off, or a '''k'''-point set which has not been reduced using symmetry considerations must be applied. VASP.5.X is able to change the '''k'''-point set on the fly and the previous restriction does not apply.
+Related tags: {{TAG|NSW}},
+{{TAG|POTIM}},
+{{TAG|MDALGO}},
+{{TAG|SMASS}},
+{{TAG|NFREE}},
+{{TAG|ISIF}},
+{{TAG|LEPSILON}},
+{{TAG|LCALCEPS}}
-== {{TAG|IBRION}}=7 and 8: second derivatives, Hessian matrix and phonon frequencies (perturbation theory).==
+Related files: {{FILE|POSCAR}}, {{FILE|CONTCAR}}, {{FILE|XDATCAR}}, {{FILE|vaspout.h5}}
-{{TAG|IBRION}}=7 and {{TAG|IBRION}}=8 are only supported as of VASP.5.1. It determines the Hessian matrix (matrix of second derivatives) using density functional perturbation theory. As for {{TAG|IBRION}}=5, {{TAG|IBRION}}=7 does not apply symmetry, whereas {{TAG|IBRION}}=8 uses symmetry to reduce the number of displacements. The output is similar as for [[#|{{TAG|IBRION}}=5 and 6]]. The only exception is that the ionic relaxation contributions to the elastic moduli are presently not determined. Born effective charges and piezoelectric constants can be calculated by specifying {{TAG|LEPSILON}}=.TRUE.
-== {{TAG|IBRION}}=44 ==
+Related topics and how-to pages:  [[Time-propagation algorithms in molecular dynamics]],
+[[Structure optimization]],
+[[POSCAR#Full_format_specification|Selective dynamics]],
+[[Computing the phonon dispersion and DOS]],
+[[Transition states]],
+[[IRC calculations]],
+[[Improved Dimer Method]],
+[[Writing a Python plugin]]
-== Related Tags and Sections ==
+{{sc|IBRION|Examples|Examples that use this tag}}
-{{TAG|NSW}},
-{{TAG|SMASS}},
-{{TAG|POTIM}},
-{{TAG|NFREE}}
-== References ==
-<references>
-<ref name="pulay:cpl:80">[http://dx.doi.org/10.1016/0009-2614(80)80396-4 P. Pulay, Chem. Phys. Lett. 73, 393 (1980).]</ref>
-<ref name="press">W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, ''Numerical Recipes'' (Cambridge University Press, New York, 1986).</ref>
-</references>
-----
-[[The_VASP_Manual|Contents]]
-[[Category:INCAR]][[Category:Dynamics]]
+[[Category:INCAR tag]][[Category:Ionic minimization]][[Category:Molecular dynamics]][[Category:Phonons]][[Category:Transition states]]