Editing IDM-IRC (section)

== Detailed Step-by-Step Example of IDM+IRC ''(by David)'' ==

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 [http://othes.univie.ac.at/10117/ 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 &#197; to 1.50 &#197;). No matter how the IS and FS look like, we suppose that in the transition state the S atom is bound over <i>bridge</i> site with the H atom <i>fcc</i>-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.
{| border="1" cellpadding="4"
 |+ INCAR, KPOINTS, POSCAR, POTCAR
 |-
 | INCAR || KPOINTS || POSCAR || POTCAR
 |- valign="top" |
 | <TT>general:<br />   ISTART = 0<br />   ICHARG = 2<br />   SYSTEM = IDM PV<br />   ISMEAR = 1<br />   SIGMA = 0.2<br />   ALGO   = V<br />   ENCUT  = 400<br />   EDIFFG = -0.01<br />   EDIFF  = 1E-5<br />   LREAL  = Auto<br />   NWRITE = 3<br />dynamic:<br />   NSW    = 1<br />   POTIM  = 0.02<br />   NFREE  = 2<br />   IBRION = 5</TT>
 | valign="top" | <TT>K-points<br />0<br />Monkhorst Pack<br />3  3  1<br />0  0  0</TT>
 | valign="top" | <TT>Methane thiol on Pt(111) - Pt15CH4S<br />
1.0<br />
+4.8847729580  +0.0000000000  +0.0000000000 <br />
+2.4423844848  +4.2303363440  +0.0000000000 <br />
+0.0000000000  +0.0000000000 +30.2160342448 <br />
Pt C H S<br />
15 1 4 1<br />
Selective<br />
Cartesian<br />
+5.6988999391  +1.4101110007  +0.0000000000  F F F<br />
+3.2565110124  +2.8202209861  +0.0000000000  F F F<br />
+3.2565120000  +0.0000000000  +0.0000000000  F F F<br />
+0.0000000000  +0.0000000000  +2.3026968902  F F F<br />
+2.4423830055  +1.4101110007  +2.3026968902  F F F<br />
+4.8847669880  +2.8202209861  +2.3026968902  F F F<br />
+6.5130230124  +2.8202209861  +4.6054279245  F F F<br />
+1.6282560244  +0.0000000000  +4.6054279245  F F F<br />
+4.0706389811  +1.4101110007  +4.6054279245  F F F<br />
+0.8175384611  +1.3965016466  +6.9003379459  F F F<br />
+3.2445703658  +2.8010128820  +6.8723385942  F F F<br />
+5.6719750456  +4.1939255595  +6.9125466994  F F F<br />
+4.8043915705  +0.0768282988  +9.2088498811  F F F<br />
+2.3123278652  +1.4423788632  +9.2437528525  F F F<br />
+4.8502102632  +2.9085629024  +9.2590920347  F F F<br />
+4.0047175525  +1.6615844468 +12.7666445850  T T T<br />
+5.0834210571  +1.7888803172 +12.9041598137  T T T<br />
+3.7347823663  +0.5994492394 +12.7699390555  T T T<br />
+3.4353761210  +2.2061660991 +13.5293434876  T T T<br />
+2.8327281158  +3.6631835395 +11.1212059991  T T T<br />
+3.5215200590  +2.3281542653 +11.0836576860  T T T</TT>
 | valign="top" style="font-style:italic; color:green;" | <i>concatenate:<br /> Pt<br /> C<br /> H<br /> S</i>
 |}

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 <i>k</i>-point grid as you do for IS and FS.
{| border="1" cellpadding="4"
 |+ INCAR, KPOINTS, POSCAR, POTCAR
 |-
 | INCAR || KPOINTS || POSCAR || POTCAR
 |- valign="top" |
 | <TT>general:<br />   ISTART = 0<br />   ICHARG = 2<br />   SYSTEM = IDM PV<br />   ISMEAR = 1<br />   SIGMA = 0.2<br />   ALGO   = V<br />   ENCUT  = 400<br />   EDIFFG = -0.01<br />   EDIFF  = 1E-5<br />   LREAL  = Auto<br />dynamic:<br />   NSW    = 1000<br />   POTIM  = 0.2<br />   IBRION = 44</TT>
 | valign="top" | <TT>K-points<br />0<br />Monkhorst Pack<br />6  6  1<br />0  0  0</TT>
 | valign="top" style="font-style:italic; color:green;" | <i>generated<br />by script<br /><br />enable surface atoms<br /> nr. 13-15<br />to be relaxed!!!</i>
 | valign="top" style="font-style:italic; color:green;" | <i>no changes</i>
 |}

: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 <TT>IBRION=1</TT> and <TT>POTIM=0.1</TT>. 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" (<TT>DVVMINUS=.FALSE.</TT>) and "backward" (<TT>DVVMINUS=.TRUE.</TT>):
 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 <TT>IBRION=1</TT> and <TT>POTIM=0.1</TT>.

[[User:Dkarhanek|Dkarhanek]] 19:10, 3 June 2011 (CEST)