IDM: Difference between revisions

From Wiki
Jump to navigation Jump to search
Dkarhanek (talk | contribs)
mNo edit summary
Rgarcia (talk | contribs)
 
(41 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 =


'''VASP_IDM''' is our internal abbreviation for a set of extra routines for VASP, developed by Dr. Tomáš Bučko ([http://www.cms.tuwien.ac.at/cms/person/33/ Univ.Vienna]), which is useful for setting up calculations on transition state search.
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 loaded with the module ''5.2_IDM'', 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.


In the following, the (currently not published elsewhere) manual 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:


== Improved Dimer Method (IDM) ==
frequencies:
  IBRION =   5    #   
  EDIFF  = 1E-6  #
  POTIM  =   0.010 #                 
  NFREE  =   2    #


The current implementation does not support lattice optimisation and can be used only for atomic relaxations (<TT>ISIF=2</TT>).
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:


<p align="center" style="width:100%">
  getdimer ../freq/POSCAR ../freq/OUTCAR POSCAR
{| border="1"
  |+ '''Flags and parameters'''
|-
| <TT>IBRION=44</TT>
| invokes relaxation with the dimer method
|-
| <TT>FINDIFF=1&#124;2</TT>
| forward (1) of central (2) differences formula for numerical differentiation
|-
| <TT>DIMER DIST=0.01</TT>
| step for numerical differentiation (&#197;)
|-
| <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>".
4. Specify the following in your INCAR file:


(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]
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.


(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]
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


== Intrinsic Reaction Cooordinate (IRC) ==
= Script and Executable =


The current implementation does not support lattice optimization and can be used only for atomic relaxations (<TT>ISIF=2</TT>).
''' Instructions '''


<p align="center" style="width:100%">
1. Download the program [[Image:dimerdiag.tgz]]. Extract it and put it in your ~/bin/exe/ folder
{| border="1"
|+ '''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):
2. Copy the script in the ~/bin folder. Recommended name: getdimer.  


: ''  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>  ''
3. Edit the script; change "rgarcia" by your username in this part:
: ''  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:
# PART 3 - DIAGONALIZE
/home/rgarcia/bin/exe/dimerdiag


: ''  &#124;v<sub>i</sub>&#124; = v<sub>0</sub>  .''
4. Now you can go to any folder and generate your dimer:


The size of integration step is varied during the simulation such as to ensure optimal performance:
$ getdimer <input:POSCAR> <input:OUTCAR> <output:POSCAR>


: ''  &#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> ,''
''''' 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


&#916;t<sub>i</sub> is constrained to interval <TT>DVVMINPOTIM</TT> < &#916;t<sub>i</sub> < <TT>DVVMAXPOTIM</TT>.
= Alternative (old) versions =
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]
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'')


== Hands-on Guide for IDM+IRC ==
--[[User:Rgarcia|Rgarcia]] 12:32, 18 November 2014 (CET)


'''Finding a transition state (TS) using the improved dimer method (IDM)'''
go back to [[Main Page]], [[Computational Resources]], [[Chemistry & More]], [[Computational Codes]], [[VASP]], [[Scripts for VASP]]
 
: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.
:Of course, you can create a reasonable guess using tools like <TT>nebmake.pl</TT> (VTST tools) from IS and FS, or resubmit a TS resulted 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>). Hopefully you have a "hard" imaginary vibration mode (no matter that there are some soft ones). Extract the displacement vector of the imaginary mode into a file DIMERDIRECTION, e.g. using [[get-dimer-dir.sh]] script.
 
:3. ''Submit the IDM job''. In comparison to a usual relaxation run (INCAR+KPOINTS+POSCAR+POTCAR), an additional input file for the module <TT>VASP/5.2_IDM</TT> is the previously created DIMERDIRECTION file. In the INCAR file, set the <TT>NSW</TT>-tag to a very high number such as 4000 or so, the calculation needs a lot of ionic steps. The most important keyword switching on the IDM method is <TT>IBRION=44</TT> tag.
 
: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 hill from the TS''. Use the CONTCAR from a preceding IDM calculation as POSCAR for IRC (contains the dimer vector inside). 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.
 
('''Warning:''' ''Under Construction: Not tested everything yet!'')
 
--[[User:Dkarhanek|Dkarhanek]] 18:29, 29 September 2010 (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