IDM: Difference between revisions

From Wiki
Jump to navigation Jump to search
Dkarhanek (talk | contribs)
Adapting to VASP 5.2.11+
Rgarcia (talk | contribs)
 
(40 intermediate revisions by 2 users not shown)
Line 1: Line 1:
go back to [[Main Page]], [[Computational Resources]], [[Chemistry & More]], [[Computational Codes]], [[VASP]]
go back to [[Main Page]], [[Computational Resources]], [[Chemistry & More]], [[Computational Codes]], [[VASP]], [[Scripts for VASP]]


== DESCRIPTION ==
= Instructions =


Effective algorithm for ''transition state search (→IDM)'' and subsequent ''path sampling (→IRC)'' was implemented into VASP by Dr. Tomáš Bučko ([http://www.cms.tuwien.ac.at/cms/person/33/ 1]/[http://www.fns.uniba.sk/index.php?id=2391 2]) in 2010.
Before start, please read about the Improved Dimer Method by Henkelman & Jónsson [https://aip.scitation.org/doi/10.1063/1.480097] and implemented in VASP by Heyden[http://dx.doi.org/10.1063/1.2104507]. Then kindly follow the procedures detailed here:


Using the corresponding VASP binary ("<TT>module load vasp/5.2.11</TT>" or later version), the '''improved dimer method (IDM)''' (<TT>IBRION=44</TT>) and '''intrinsic reaction coordinate (IRC)''' search (<TT>IBRION=40</TT>) become available.
1. Generate your initial guess for the structure (POSCAR). To this end, you can either pre-converge with a NEB calculation, or do it by hand.


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.
2. Be sure the metal or solid phase is frozen and the molecule is relaxed. Do a frequency calculation based on your initial guess. Specify the following in your INCAR file:


'''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.
frequencies:
  IBRION =  5    #   
  EDIFF  =  1E-6  #
  POTIM  =  0.010 #                 
  NFREE =  2    #


3. Based on your POSCAR and OUTCAR of the frequencies calculation, generate in a new folder the input of the IDM with the getdimer script. For example:


== Improved Dimer Method (IDM) ==
getdimer ../freq/POSCAR ../freq/OUTCAR POSCAR


The current implementation does not support lattice optimisation and can be used only for atomic relaxations (<TT>ISIF=2</TT>).
4. Specify the following in your INCAR file:


<p align="center" style="width:100%">
improved dimer method:
{| border="1"
  IBRION 44    # IDM   
  |+ '''Flags and parameters'''
  EDIFF =  1E-6 #  
  |-  
  EDIFFG = -0.030 # use force criteria
  | <TT>IBRION=44</TT>
  POTIM =   0.050 # Between 0.010 and 0.100.                 
| # invokes relaxation with the dimer method
  NFREE =  2    # Between 2&5
  |-  
  FINDIFF=  2    #
  | <TT>FINDIFF=1&#124;2</TT>
  DIMER_DIST=0.010 #              
  | # forward (1) of central (2) differences formula for numerical differentiation
  MINROT =   0.010 #  
|-
  STEP_SIZE= 0.010 #              
| <TT>DIMER DIST=0.01</TT>
  STEP_MAX= 0.100 #               
| # step for numerical differentiation (&#197;)
  NSW    = 100    # Never use more than 200. Check periodically for divergence.
|-
| <TT>STEP SIZE=0.01</TT>
| # trial step size for energy minimisation(&#197;)
|-
| <TT>STEP MAX=0.1</TT>
| # maximal step size for energy minimisation (&#197;)
  |-
| <TT>MINROT=0.01</TT>
| # minimal rotation of dimer (rad.)
|}
</p>


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 hypersurface. It must be specified in <TT>POSCAR</TT> on place of ionic velocities (see documentation for the <TT>POSCAR</TT> 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 <TT>OUTCAR</TT> after key words "<TT>DIMER METHOD</TT>".
5. To check that it is converging well, (1) check with p4vasp that all forces are small, and (2) check that the curvature is always negative. Normal values are between -1.0 and -30.0:
grep curv OUTCAR


(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. [http://dx.doi.org/10.1063/1.480097]
= Script and Executable =


(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. [http://dx.doi.org/10.1063/1.2104507]
''' Instructions '''


1. Download the program [[Image:dimerdiag.tgz]]. Extract it and put it in your ~/bin/exe/ folder


== Intrinsic Reaction Cooordinate (IRC) ==
2. Copy the script in the ~/bin folder. Recommended name: getdimer.


The current implementation does not support lattice optimization and can be used only for atomic relaxations (<TT>ISIF=2</TT>).
3. Edit the script; change "rgarcia" by your username in this part:


<p align="center" style="width:100%">
  # PART 3 - DIAGONALIZE
{| border="1"
  /home/rgarcia/bin/exe/dimerdiag
|+ '''Flags and parameters'''
|-
| <TT>IBRION=40</TT>
  | # perform IRC search using the DVV method
|-
| <TT>DVVEHISTORY=5</TT>
| # number of subsequent increasing energy steps taken before simulation is terminated
|-
| <TT>DVVDELTA0=1.5e-3</TT>
| # the error tolerance &#916;<sub>0</sub>
|-  
  | <TT>DVVVNORM0=0.01</TT>
| # magnitude of the velocity vector v<sub>0</sub>
|-
| <TT>DVVMINPOTIM=0.025</TT>
| # lower limit for the time step &#916;''t''
|-
| <TT>DVVMAXPOTIM=3.0</TT>
| # upper limit for the time step &#916;''t''
|-
| <TT>DVVMINUS=.FALSE.&#124;.TRUE.</TT>
|  # take positive|negative direction of the initial velocity vector
|}
</p>


The method is described in detail in Ref.(3). Velocity verlet algorithm is used to drive system from (well relaxed!!!) saddle point to nearby minimum along intrinsic reaction coordinate (IRC):
4. Now you can go to any folder and generate your dimer:


: '' x<sub>i</sub> = x<sub>i-1</sub> + v<sub>i-1</sub>&#916;t<sub>i</sub> + &#189;a<sub>i-1</sub>&#916;t<sub>i</sub><sup>2</sup>  ''
  $ getdimer <input:POSCAR> <input:OUTCAR> <output:POSCAR>
: ''  v<sub>i</sub> = v<sub>i-1</sub> + &#189;(a<sub>i-1</sub>+a<sub>i</sub>)&#916;t<sub>i</sub> .''


Direction of atomic velocities (i.e. vector of negative curvature) must be defined in POSCAR. Note that if DVVMINUS is set to =.FALSE., oposite 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:
''''' SCRIPT '''''
#!/bin/bash
#
# Generates the POSCAR of a dimer calculation from the non mass-weighted Hessian matrix
# Rodrigo García-Muelas
# November 14th, 2014 - Created
# March    9th, 2015 - Debugged 1 (Thanks to Marçal & Luca)
#
# Use with VASP 5.x or higher
#
# INPUT
# $1 Path of the POSCAR file containing the coordinates
# $2 Path of the OUTCAR file containing the frequencies of the structure
# $3 <OUTPUT> POSCAR file
#
# PART 0 - SECURITY CHECKS
if [ -e $3 ] ; then
echo "Warning! $3 already exist. Overwrite? (y/Y for yes)"
read  overwrite
case $overwrite in
y|Y|Yes|yes|YES) echo "$3 will be overwriten"          ;;
*)              echo "No actions taken "    ; exit 1 ;;
esac
fi
IBRION=`grep IBRION $2 | awk '{print $3}'`
if <tt><nowiki> [[ </nowiki></tt> $IBRION -gt "8" || $IBRION -lt "5" ]] ; then
echo "$2 is not an OUTCAR of a frequency job."
echo "No actions taken" ; exit 1 ; fi
# PART 1 - READ POSCAR AND PREPARE INPUTS
cp $1 diagonalizer_poscar.tmp
# PART 2 - READ OUTCAR
grep "Degrees of freedom DOF " $2 > temp01.tmp
freedom=`awk '{print $6}' temp01.tmp`
echo "$freedom" > diagonalizer_matrix.tmp
grep -A $(($freedom+2)) 'SECOND DERIVATIVES' $2 > temp02.tmp
cut -b 5- temp02.tmp > temp03.tmp
tail -n $freedom temp03.tmp >> diagonalizer_matrix.tmp
# PART 3 - DIAGONALIZE
/home/rgarcia/bin/exe/dimerdiag
# PART 4 - CLEAN
head -n 9 $1 > $3
cat    diagonalizer_taildm.tmp >> $3
mv -f  diagonalizer_output.tmp eigenvectors.dat
rm -fv *.tmp


: ''  &#124;v<sub>i</sub>&#124; = v<sub>0</sub>  .''
= Alternative (old) versions =


The size of integration step is varied during the simulation such as to ensure optimal performance:
WARNING: Some of the topics discussed by David Karhánek are outdated. For instance, it is not necessary to set NWRITE=3 for the frequency calculation: [[IDM-IRC]] [[get-dimer-dir.sh]] [[get-irc-path.sh]] (''by David Karhanek'')


: ''  &#916;t<sub>i+1</sub> = &#916;t<sub>i</sub>&#124;&#916;<sub>0</sub>/&#916;<sub>i</sub> &#124;<sup>1/3</sup> ,''
--[[User:Rgarcia|Rgarcia]] 12:32, 18 November 2014 (CET)


&#916;t<sub>i</sub> is constrained to interval <TT>DVVMINPOTIM</TT> < &#916;t<sub>i</sub> < <TT>DVVMAXPOTIM</TT>.
go back to [[Main Page]], [[Computational Resources]], [[Chemistry & More]], [[Computational Codes]], [[VASP]], [[Scripts for VASP]]
Output (in particular value of intrinsic reaction coordinate (IRC) and corresponding potential energy) is written into <TT>OUTCAR</TT> file after key words "<TT>DAMPED VELOCITY VERLET ALGORITHM:</TT>"
 
(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. [http://dx.doi.org/10.1021/jp012125b]
 
 
== Hands-on Guide for IDM+IRC ''(by David)'' ==
 
'''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 <TT>nebmake.pl</TT> (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. (<TT>IBRION=5</TT>, <TT>NFREE=2</TT>, very small <TT>POTIM</TT>). Use low ''k''-point grid in <TT>KPOINTS</TT>, 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. 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 <TT>IBRION=44</TT> tag. Don't worry to setup the <TT>NSW</TT>-tag in the INCAR file to a very high number such as 4000 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.
 
:4. ''Check the TS you've obtained'',  
i.e. run a frequency job (<TT>IBRION=5</TT>) 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 <TT>IBRION=40</TT> tag. Furthermore, you need 2 directories, in which the INCARs differ by the setting for sliding down the forward (<TT>DVVMINUS=.FALSE.</TT>) and reverse (<TT>DVVMINUS=.TRUE.</TT>) direction from the saddle point (the TS).
 
:2. ''Explore the full energy path''.
Put the energies of the forward and reverse VASP run together and enjoy your full energy trajectory by pasting the output of the script [[get-irc-path.sh]] to your favourite plotting utility.
 
[[User:Dkarhanek|Dkarhanek]] 17:58, 1 June 2011 (CEST)

Latest revision as of 18:20, 5 November 2019

go back to Main Page, Computational Resources, Chemistry & More, Computational Codes, VASP, Scripts for VASP

Instructions[edit]

Before start, please read about the Improved Dimer Method by Henkelman & Jónsson [1] and implemented in VASP by Heyden[2]. Then kindly follow the procedures detailed here:

1. Generate your initial guess for the structure (POSCAR). To this end, you can either pre-converge with a NEB calculation, or do it by hand.

2. Be sure the metal or solid phase is frozen and the molecule is relaxed. Do a frequency calculation based on your initial guess. Specify the following in your INCAR file:

frequencies:
  IBRION =   5     #     
  EDIFF  =  1E-6   # 
  POTIM  =   0.010 #                  
  NFREE  =   2     # 

3. Based on your POSCAR and OUTCAR of the frequencies calculation, generate in a new folder the input of the IDM with the getdimer script. For example:

getdimer ../freq/POSCAR ../freq/OUTCAR POSCAR

4. Specify the following in your INCAR file:

improved dimer method:
  IBRION =  44     # IDM    
  EDIFF  =   1E-6  # 
  EDIFFG =  -0.030 # use force criteria
  POTIM  =   0.050 # Between 0.010 and 0.100.                  
  NFREE  =   2     # Between 2&5
  FINDIFF=   2     # 
  DIMER_DIST=0.010 #               
  MINROT =   0.010 # 
  STEP_SIZE= 0.010 #                
  STEP_MAX=  0.100 #                
  NSW    = 100     # Never use more than 200. Check periodically for divergence.

5. To check that it is converging well, (1) check with p4vasp that all forces are small, and (2) check that the curvature is always negative. Normal values are between -1.0 and -30.0:

grep curv OUTCAR

Script and Executable[edit]

Instructions

1. Download the program File:Dimerdiag.tgz. Extract it and put it in your ~/bin/exe/ folder

2. Copy the script in the ~/bin folder. Recommended name: getdimer.

3. Edit the script; change "rgarcia" by your username in this part:

# PART 3 - DIAGONALIZE 
/home/rgarcia/bin/exe/dimerdiag 

4. Now you can go to any folder and generate your dimer:

$ getdimer <input:POSCAR> <input:OUTCAR> <output:POSCAR>

SCRIPT

#!/bin/bash
# 
# Generates the POSCAR of a dimer calculation from the non mass-weighted Hessian matrix
# Rodrigo García-Muelas
# November 14th, 2014 - Created
# March     9th, 2015 - Debugged 1 (Thanks to Marçal & Luca)
# 
# Use with VASP 5.x or higher
# 
# INPUT
# $1 Path of the POSCAR file containing the coordinates
# $2 Path of the OUTCAR file containing the frequencies of the structure
# $3 <OUTPUT> POSCAR file
# 

# PART 0 - SECURITY CHECKS

if [ -e $3 ] ; then
echo "Warning! $3 already exist. Overwrite? (y/Y for yes)"
read  overwrite
case $overwrite in
y|Y|Yes|yes|YES) echo "$3 will be overwriten"          ;;
*)               echo "No actions taken "     ; exit 1 ;;
esac
fi

IBRION=`grep IBRION $2 | awk '{print $3}'`
if  [[  $IBRION -gt "8" || $IBRION -lt "5" ]] ; then
echo "$2 is not an OUTCAR of a frequency job."
echo "No actions taken" ; exit 1 ; fi

# PART 1 - READ POSCAR AND PREPARE INPUTS 

cp $1 diagonalizer_poscar.tmp 

# PART 2 - READ OUTCAR 

grep "Degrees of freedom DOF " $2 > temp01.tmp
freedom=`awk '{print $6}' temp01.tmp` 
echo "$freedom" > diagonalizer_matrix.tmp

grep -A $(($freedom+2)) 'SECOND DERIVATIVES' $2 > temp02.tmp
cut -b 5- temp02.tmp > temp03.tmp
tail -n $freedom temp03.tmp >> diagonalizer_matrix.tmp

# PART 3 - DIAGONALIZE 
/home/rgarcia/bin/exe/dimerdiag 

# PART 4 - CLEAN
head -n 9 $1 > $3
cat    diagonalizer_taildm.tmp >> $3
mv -f  diagonalizer_output.tmp eigenvectors.dat
rm -fv *.tmp

Alternative (old) versions[edit]

WARNING: Some of the topics discussed by David Karhánek are outdated. For instance, it is not necessary to set NWRITE=3 for the frequency calculation: IDM-IRC get-dimer-dir.sh get-irc-path.sh (by David Karhanek)

--Rgarcia 12:32, 18 November 2014 (CET)

go back to Main Page, Computational Resources, Chemistry & More, Computational Codes, VASP, Scripts for VASP