IDM-IRC
go back to Main Page, Computational Resources, Chemistry & More, Computational Codes, VASP
To see newer comments about the Improved Dimer Method, please see IDM, written by Rodrigo García-Muelas.
DESCRIPTION[edit]
Effective algorithm for transition state search (→IDM) and subsequent path sampling (→IRC) was implemented into VASP by Dr. Tomáš Bučko (1/2) in 2010.
Using the corresponding VASP binary ("module load vasp/5.2.11" or later version), the improved dimer method (IDM) (IBRION=44) and intrinsic reaction coordinate (IRC) search (IBRION=40) become available.
Manual page for the IDM method has already appeared in the official VASP manual: http://cms.mpi.univie.ac.at/vasp/vasp/Improved_dimer_method.html in March 2011 and contains a nice example added. The IRC manual is, however, still missing there. Here, in the following, the manual based on older text for both methods, supplied by the developer, is presented.
Note: The implementation of Dimer Method by Univ.Texas within their VTST tools is different, it is not the "improved" version and uses different setup.
Improved Dimer Method (IDM)[edit]
The current implementation does not support lattice optimisation and can be used only for atomic relaxations (ISIF=2).
| IBRION=44 | # invokes relaxation with the dimer method |
| FINDIFF=1|2 | # forward (1) of central (2) differences formula for numerical differentiation |
| DIMER DIST=0.01 | # step for numerical differentiation (Å) |
| STEP SIZE=0.01 | # trial step size for energy minimisation(Å) |
| STEP MAX=0.1 | # maximal step size for energy minimisation (Å) |
| MINROT=0.01 | # minimal rotation of dimer (rad.) |
The method is described in detail in Ref.(2). The only input which has to be defined by user is a 3N-dimensional vector defining direction of negative curvature on potential energy surface. It must be specified in POSCAR on place of ionic velocities (see documentation for the POSCAR file). Note that this vector is automatically normalised so it is the only direction that really matters. Output including maximal gradient in the current step, curvature along the dimer direction and the angle through which the dimer is rotated is written in OUTCAR after key words "DIMER METHOD".
(1) G. Henkelman and H. Jónsson. A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives. J. Chem. Phys., 111(15):7010–7022, 1999. [1]
(2) A. Heyden, A.T. Bell, and F.J. Keil. Efficient methods for finding transition states in chemical reactions: Comparison of improved dimer method and partitioned rational function optimization method. J. Chem. Phys., 123(22):224101, 2005. [2]
Intrinsic Reaction Cooordinate (IRC)[edit]
The current implementation does not support lattice optimization and can be used only for atomic relaxations (ISIF=2).
| IBRION=40 | # perform IRC search using the DVV method |
| DVVEHISTORY=5 | # number of subsequent increasing energy steps taken before simulation is terminated |
| DVVDELTA0=1.5e-3 | # the error tolerance Δ0 |
| DVVVNORM0=0.01 | # magnitude of the velocity vector v0 |
| DVVMINPOTIM=0.025 | # lower limit for the time step Δt |
| DVVMAXPOTIM=3.0 | # upper limit for the time step Δt |
| DVVMINUS=.FALSE.|.TRUE. | negative direction of the initial velocity vector |
The method is described in detail in Ref.(3). Velocity verlet algorithm is used to drive the system from a converged saddle point into a nearby minimum along intrinsic reaction coordinate (IRC):
- xi = xi-1 + vi-1Δti + ½ai-1Δti2
- vi = vi-1 + ½(ai-1+ai)Δti .
Direction of atomic velocities (i.e. vector of negative curvature) must be defined in POSCAR. Note that if DVVMINUS is set to =.FALSE., opposite direction of vector defined in POSCAR is used to initiate atomic velocities. In order to keep intermediate structures close to IRC, atomic velocities are damped such as the norm of velocity vector is constant:
- |vi| = v0 .
The size of integration step is varied during the simulation such as to ensure optimal performance:
- Δti+1 = Δti|Δ0/Δi |1/3 ,
Δti is constrained to interval DVVMINPOTIM < Δti < DVVMAXPOTIM. Output (in particular value of intrinsic reaction coordinate (IRC) and corresponding potential energy) is written into OUTCAR file after key words "DAMPED VELOCITY VERLET ALGORITHM:"
(3) H.P. Hratchian and H.B. Schlegel. Following reaction pathways using a damped classical trajectory algorithm. J. Phys. Chem. A, 106(1):165–169, 2002. [3]
Brief Hands-on Guide for IDM+IRC (by David)[edit]
Finding a transition state (TS) using the improved dimer method (IDM)
- 1. Suggest the geometry of the TS.
Do not care much about the structure of initial state (IS) and final state (FS) of the reaction as it is for (CI-)NEB calculations, simply stretch out bond(s) to be broken or shorten bond(s) to be formed, by hand in a visualizer like P4V. Of course, you can create a reasonable TS guess using tools like nebmake.pl (from VTST tools) from IS and FS, or re-use a TS result structure from (CI-)NEB calculation.
- 2. Run a frequency job for the TS guess,
i.e. ensure IBRION=5, NFREE=2, very small POTIM=0.02 and NWRITE=3 (important!) in your INCAR file. Use low k-point grid in KPOINTS, because there is no need for high accuracy at this step. Hopefully you have a "hard" imaginary vibration mode (no matter that there are some very soft ones). Optionally, you can check the modes using nfreq.sh script in MOLDEN. Now, extract the displacement vector of the imaginary mode, preferably using get-dimer-dir.sh script.
- 3. Submit the IDM job.
The most important keyword switching on the IDM method is IBRION=44 tag. Don't worry to setup the NSW-tag in the INCAR file to a very high number such as 1000+ or so, the calculation may need a lot of ionic steps. The input files are the same as for usual relaxation run (INCAR+KPOINTS+POSCAR+POTCAR), the POSCAR containing the eigenvector. Use VASP 5.2.11 (or later) binary. Sometimes, it is needed to refine the TS structure via resubmitting for classical optimization with IBRION=1 and small POTIM=0.1.
- 4. Check the TS you've obtained,
i.e. run a frequency job (IBRION=5) and check that you have 1 "hard" imaginary vibration mode. Well done! Afterwards, you may want to find the corresponding IS and TS (incl. full energy profile) using the IRC method (vide infra).
Finding the initial state (IS) and final state (FS) using the improved dimer method (IRC)
- 1. Slide down the energy saddle from the TS.
Use the CONTCAR from a preceding IDM calculation as POSCAR (contains the dimer vector inside) for the two IRC runs. The most important INCAR keyword switching on the IRC method is IBRION=40 tag. Furthermore, you need 2 directories, in which the INCARs differ by the setting for sliding down the forward (DVVMINUS=.FALSE.) and reverse (DVVMINUS=.TRUE.) direction from the saddle point (the TS).
- 2. Explore the full energy path.
Put the energies of the forward and reverse runs together and enjoy your full energy trajectory by pasting the output of the script get-irc-path.sh to your favourite graph-plotting utility.
Detailed Step-by-Step Example of IDM+IRC (by David)[edit]
Let's see how it works. Supposed that you have relaxed structures of the initial (IS) and final state (FS), you want to find the transition state (TS). The example is a surface reaction (simplified) from David's PhD thesis.
Finding a transition state (TS) using the improved dimer method (IDM)
- 1. Suggest the geometry of the TS.
We have a methanethiol molecule adsorbed on Pt(111) surface and the reaction occuring is the breaking of S-H bond. So, for our inital guess of the TS, let's simply elongate the S-H bond manually in p4vasp (from 1.36 Å to 1.50 Å). No matter how the IS and FS look like, we suppose that in the transition state the S atom is bound over bridge site with the H atom fcc-likepointing out towards a surface atom.
- 2. Run a frequency job for the TS guess
Let's fix the adsorbate atoms and perform vibrational analysis, create a "PV" (preliminary vibrations) directory and setup the calculation.
| INCAR | KPOINTS | POSCAR | POTCAR |
| general: ISTART = 0 ICHARG = 2 SYSTEM = IDM PV ISMEAR = 1 SIGMA = 0.2 ALGO = V ENCUT = 400 EDIFFG = -0.01 EDIFF = 1E-5 LREAL = Auto NWRITE = 3 dynamic: NSW = 1 POTIM = 0.02 NFREE = 2 IBRION = 5 |
K-points 0 Monkhorst Pack 3 3 1 0 0 0 |
Methane thiol on Pt(111) - Pt15CH4S 1.0 |
concatenate: Pt C H S |
Check the "quality" of your initial guess:
grep 'cm-1' OUTCAR
We want 1 hard imaginary mode, 0 or few soft ones. We may visualize them (e.g. in MOLDEN using nfreq.sh script) to be sure that we got what we wanted. When finished, use the get-dimer-dir.sh script in the "PV" directory to extract the hardest imaginary mode:
get-dimer-dir.sh
A new up-directory "TS" with the desired POSCAR for IDM run will be created.
- 3. Submit the IDM job.
Allow some surface atoms to relax if needed. Use the full k-point grid as you do for IS and FS.
| INCAR | KPOINTS | POSCAR | POTCAR |
| general: ISTART = 0 ICHARG = 2 SYSTEM = IDM PV ISMEAR = 1 SIGMA = 0.2 ALGO = V ENCUT = 400 EDIFFG = -0.01 EDIFF = 1E-5 LREAL = Auto dynamic: NSW = 1000 POTIM = 0.2 IBRION = 44 |
K-points 0 Monkhorst Pack 6 6 1 0 0 0 |
generated by script enable surface atoms nr. 13-15 to be relaxed!!! |
no changes |
- 4. Check the TS you've obtained.
Run vibrational analysis with the same settings like in step "2.". Check the modes. If your TS looks reasonable you're lucky and you might refine it by re-optimizing with IBRION=1 and POTIM=0.1. If the TS is a non-sense, than go back to step "1.".
Finding the initial state (IS) and final state (FS) using the improved dimer method (IRC)
- 1. Slide down the energy saddle from the TS.
Use the output structure of the TS to be the input for the two IRC runs, one "forward" (DVVMINUS=.FALSE.) and "backward" (DVVMINUS=.TRUE.):
mkdir 1 2 for i in 1 2 ; do cp TS/CONTCAR TS/KPOINTS TS/POTCAR TS/INCAR TS/runscript $i/ ; done for i in 1 2 ; do sed -i 's/IBRION = 44/IBRION = 40' $i/INCAR ; done echo 'DVVMINUS=.FALSE.' 1/INCAR echo 'DVVMINUS=.TRUE.' 2/INCAR
- 2. Explore the full energy path.
Use the get-irc-path.sh script:
get-irc-path.sh
to get some simple statistics and the IRC trajectory. Visualise it e.g. with:
xmgrace trajectory.dat
If you are suspicious, that you IS and FS are still not yet relaxed to the local minima, you may resubmit them as normal relaxations with IBRION=1 and POTIM=0.1.
Dkarhanek 19:10, 3 June 2011 (CEST)