Improved dimer method

The dimer method is a technique for the optimization of transition states. In VASP, the method improved by Heyden et al.~\cite{Heyden:05} (IDM) is implemented, detailed presentation of the method can be found in Ref. ~\cite{Heyden:05}. Algorithm for IDM consists of the following cyclically repeated steps: \begin{itemize} \item curvature along the dimer axis is computed using finite diferences. The initial dimer direction must be provided by user (see below). \item dimer is rotated such as its axis is parallel with the direction of the maximal negative curvature \item optimization step is taken, potential energy is maximized along the unstable direction (i.e. dimer axis) while it is minimized in all other directions \end{itemize}

The method is invoked by setting {\tt IBRION}=44 (see Sec.~\ref{sec:IBRION}).\index{INCAR!I!IBRION|textit} %(Martijn, please add a note about the dimer method in the section describing IBRION parameter).

Furthermore, user must specify direction of the unstable mode. Corresponding 3N dimensional vector is defined in the POSCAR file after the lines with atomic coordinates and a separating blank line. Note that the dimer direction is automatically normalized, i.e. the norm of the dimer axis defined by user is irrelevant. Example of POSCAR file for simulation with dimer method:

\begin{verbatim} ammonia flipping 1. 6. 0. 0. 0. 7. 0. 0. 0. 8. H N 3 1 cart

      -0.872954        0.000000       -0.504000        ! coordinates for atom 1
       0.000000        0.000000        1.008000
       0.872954        0.000000       -0.504000
       0.000000        0.000000        0.000000        ! coordinates for atom N
      ! here we define trial unstable direction:
       0.000001    0.522103   -0.000009        ! components for atom 1
      -0.000006    0.530068    0.000000
      -0.000005    0.522067   -0.000007
       0.000001   -0.111442    0.000001        ! components for atom N

\end{verbatim}

As in the other structural optimization algorithms in VASP, convergence is controlled through the {\tt EDIFFG} tag.

Experienced users can affect the performance of the dimer method by modifying the numerical values of the following parameters:

\begin{itemize} \item {\tt FINDIFF}= 1 $|$ 2 \\ Use a forward (1) or central (2) difference formula for the numerical differentiation to compute the curvature along the dimer direction\\

Default: {\tt FINDIFF}=1

\item {\tt DIMER\_DIST}=[real] the step size for a numerical differentiation ({\AA})\\

Default: {\tt DIMER\_DIST}=0.01

\item {\tt MINROT}=[real] dimer is rotated only if the predicted rotation angle is greater than MINROT (rad.)\\

Default: {\tt MINROT}=0.01

\item {\tt STEP\_SIZE}=[real] trial step size for optimization step ({\AA})\\

Default: {\tt STEP\_SIZE}=0.01

\item {\tt STEP\_MAX}=[real] trust radius (upper limit) the optimization step ({\AA})\\

Default: {\tt STEP\_MAX}=0.1 \end{itemize}

Important information about the progress of optimization is written in the file OUTCAR after the expression 'DIMER METHOD'.

In particular, it is useful to check the curvature along the dimer direction, which should be a negative number (long sequence of positive numbers usually indicates that the algorithm fails to converge to the correct transition state).\\

IMPORTANT NOTE: The current implementation does not support lattice optimizations ({\tt ISIF}$>$2)\index{INCAR!I!ISIF|textit} and can be used only for the relaxation of atomic positions.

\subsubsection{Initial dimer axis}\label{hint} The direction of unstable vibrational mode can be obtained by performing vibrational analysis ({\tt IBRION}=5, see Sec.~\ref{sec:IBRION=5and6}) and taking the x-, y-, and z- components of the imaginary vibrational mode (after division by SQRT(mass)!) parallel with the reaction coordinate. Note that in order to plot "Eigenvectors after division by SQRT(mass)", {\tt NWRITE}=3\index{INCAR!N!NWRITE|textit} should be used.

\subsubsection{Practical example}\label{example} In this example, transition state for the ammonia flipping is computed. All calculations discussed here were performed using the PBE functional, Brillouin zone sampling was restricted to the gamma point. This practical example can be completed in a few seconds on a standard desktop PC. The starting structure for IDM simulation should be a reasonable guess for the transition state. POSCAR with the initial guess for the ammonia flipping:

\begin{verbatim} ammonia flipping 1. 6. 0. 0. 0. 7. 0. 0. 0. 8. H N 3 1 cart

      -0.872954        0.000000       -0.504000
       0.000000        0.000000        1.008000
       0.872954        0.000000       -0.504000
       0.000000        0.000000        0.000000

\end{verbatim}

As an input for the dimer method, direction of unstable mode (dimer axis) is needed. This can be obtained by performing vibrational analysis. The INCAR file should contain the following lines:

\begin{verbatim} NSW = 1 Prec=Normal IBRION=5  ! perform vibrational analysis NFREE=2  ! select central differences algorithm POTIM=0.02  ! step for the numerical differenciation NWRITE=3  ! write down eigenvectors of dynamical matrix after division by SQRT(mass) \end{verbatim}

After completing the vibrational analysis, we look up the hardest imaginary mode (Eigenvectors after division by SQRT(mass)!) in the OUTCAR file:

\begin{verbatim}

12 f/i=   23.224372 THz   145.923033 2PiTHz  774.681641 cm-1    96.048317 meV
            X         Y         Z           dx          dy          dz
     5.127046  0.000000  7.496000     0.000001    0.522103   -0.000009
     0.000000  0.000000  1.008000    -0.000006    0.530068    0.000000
     0.872954  0.000000  7.496000    -0.000== References ==

---

Contents 005 0.522067 -0.000007

     0.000000  0.000000  0.000000     0.000001   -0.111442    0.000001

\end{verbatim}

and use the last three columns to define the dimer axis in POSCAR:

\begin{verbatim}== References ==

---

Contents

ammonia flipping 1. 6. 0. 0. 0. 7. 0.== References ==

---

Contents

0. 0. 8. H N 3 1 cart

      -0.872954        0.000000       -0.504000        ! coordinates for atom 1
       0.000000        0.000000        1.008000
       0.872954        0.000000       -0.504== References ==

---

Contents 000

       0.000000        0.000000        0.000000        ! coordinates for atom N
      ! here we define trial unstable direction:
       0.000001    0.522103   -0.000009      == References ==

---

Contents== References ==

---

Contents

 ! components for atom 1
      -0.000006    0.530068    0.000000
      -0.000005    0.522067   -0.000007== References ==

---

Contents== References ==

---

Contents

References

---

Contents

References

---

Contents

References

---

Contents

       0.000001   -0.111442    0.000001        ! components for atom N

\end{verbatim}

In order to perform IDM calculation, INCAR should contain the following lines:

\begin{verbatim} NSW = 100 Prec=Normal IBRION=44  ! use the dimer method as optimization engine EDIFFG=-0.03 \end{verbatim}

With this setting, algorithm converges in just a few relaxation steps. Further vibrational analysis can be performed to prove that the relaxed structure is indeed a first order saddle point (one imaginary frequency).

References

---

Contents