Protocols for transition state searches
Posted: Fri May 19, 2017 9:12 am
Dear all,
I would appreciate your thoughts on reliable protocols for TS optimisations. I am not so much after the actual methods that are available (NEB, dimer), but which steps you follow to set up TS searches and what to do when you encounter problems. I am currently trying to map out reaction profiles of organometallic species (viz. rhodium containing complexes) in the solid state, using periodic VASP calculations to do so. I somewhat despair at the moment about the fact that it seems to be awkwardly difficult to optimise to transition states. I have to do a large number of such calculations, and would need much more reliable protocols to efficiently advance this project.
From previous calculations in gas phase and other periodic codes I know exactly what the transition states should look like, they always involve transfer of hydrogen atoms from the metal center to carbon centres of a bound substrate, plus some rearrangements of the substrate. I use a number of protocols at the moment.
(i) One involves distorting the atoms to an approximate TS guess, starting from an optimised ground state structure. I then perform partial phonon calculations, only including relevant atoms that will be involved in the TS. In most cases this works pretty well, and I get very decent guess structures and one imaginary mode that corresponds to the reaction path. I then use the vector for this mode along with the same structure in a (advanced) dimer setup, hopefully converging to the desired TS. What I frequently see (in 90% of cases) is that the molecule explodes or selected atoms just float away, giving complete garbage structures. I don't understand why this happens, since the guess structures are pretty good (one imaginary mode). So no success with this in some cases (in a few cases this does work and I smoothly converge to the desired TS).
(ii) I have been told that the quasi-Newton method is an option to nail TSs, provided the structure is close enough. I have employed this approach numerous times, and in a few cases I had some success, but often the structure converges to the nearest minimum.
(iii) NEB seems to be an option, however, since I have a very good idea what the TS looks like, I think the dimer method should provide the best option. Just how can I make this more reliable?
I would appreciate your input on this problem, I am sure there are people with more expertise than I currently have that could give some advice. In particular I would be interested in
(i) Which parameters should I play around with? POTIM is certainly one of them, but so far changing this value in either direction has not changed much.
(ii) Do you normally preoptimize TS guess structures, relaxing only relevant atoms? Does this normally help to improve convergence?
(iii) Which part of the molecule do you include for the initial phonon calculation? I focus on 3-4 atoms, sometimes the whole substrate, depending on the problem. Normally this gives very nice imaginary modes that can be used as vector in the dimer method. Is it advisable to partially optimise these few atoms to the TS first, and then relax the full geometry to the TS? I want to add that in my case the remainder of the unit cell is normally fairly well converged to a minimum and doesn't do much else other than provide the crystal matrix.
Below is an INPUT file for the dimer method I use. Unfortunately I am unable to append any other file types such as structures so you could visualise what the problem is.
Your help is much appreciated.
Tobias
INCAR (dimer)
General:
SYSTEM = Rh-dcype-endo-nbd
! if enough K-points are present -5 should be used
ISMEAR = 0 ! 0: Gaussian -> semicond/insulator ; 1-N: MP; -5: Tetra+Blochl -> metal/very accurate energy/forces
SIGMA = 0.05
EDIFF = 1.0E-7
PREC = Normal ! to make ROPT=2.0E-4 for LREAL=Auto
LREAL = Auto ! cheaper for large cell's we are using
ENCUT = 600 !
ISYM = 2 # use symmetry as done for PAW PP
NELMIN = 8
LWAVE = .FALSE.
LCHARG = .FALSE.
! LVTOT = .TRUE. ! Write LOCPOT file for Toon
! LVHAR = .TRUE. ! but LOCPOT without exchange and correlation contribution for forcefields Toon
VOSKOWN = 1 ! important for GGA (PW91) for interpolation of XC...since we only use PBE=switch on (LDA: VOSKOWN=0)
LASPH = .TRUE. ! For VASP.5.X the aspherical contributions are properly accounted for in the Kohn-Sham potential
! as well. This is essential for accurate total energies and band structure calculations for f-elements
! (e.g. ceria), all 3d-elements (transition metal oxides), and magnetic atoms in the 2nd row (B-F atom),
! in particular if LDA+U or hybrid functionals or meta-GGAs are used, since these functionals often result
! in aspherical charge densities.
Van der Waals Interaction (vasp 5.3.3 patched verion):
IVDW = 12 ! switches between 0:off, 1: DFT-D3 and 2: TS-VDW (default=1)
dynamic:
IBRION = 44 ! -1: Fix atoms; 0: MD; 2: ConjGrad relax; 44: improved dimer method
NSW = 500 ! Number ionic steps
ISIF = 0 # relax ions, no change cell shape, no change cell volume (volume changes require ENCUT*1.3)
POTIM = 0.005
parallel:
LPLANE = .TRUE.
NPAR = 16
I would appreciate your thoughts on reliable protocols for TS optimisations. I am not so much after the actual methods that are available (NEB, dimer), but which steps you follow to set up TS searches and what to do when you encounter problems. I am currently trying to map out reaction profiles of organometallic species (viz. rhodium containing complexes) in the solid state, using periodic VASP calculations to do so. I somewhat despair at the moment about the fact that it seems to be awkwardly difficult to optimise to transition states. I have to do a large number of such calculations, and would need much more reliable protocols to efficiently advance this project.
From previous calculations in gas phase and other periodic codes I know exactly what the transition states should look like, they always involve transfer of hydrogen atoms from the metal center to carbon centres of a bound substrate, plus some rearrangements of the substrate. I use a number of protocols at the moment.
(i) One involves distorting the atoms to an approximate TS guess, starting from an optimised ground state structure. I then perform partial phonon calculations, only including relevant atoms that will be involved in the TS. In most cases this works pretty well, and I get very decent guess structures and one imaginary mode that corresponds to the reaction path. I then use the vector for this mode along with the same structure in a (advanced) dimer setup, hopefully converging to the desired TS. What I frequently see (in 90% of cases) is that the molecule explodes or selected atoms just float away, giving complete garbage structures. I don't understand why this happens, since the guess structures are pretty good (one imaginary mode). So no success with this in some cases (in a few cases this does work and I smoothly converge to the desired TS).
(ii) I have been told that the quasi-Newton method is an option to nail TSs, provided the structure is close enough. I have employed this approach numerous times, and in a few cases I had some success, but often the structure converges to the nearest minimum.
(iii) NEB seems to be an option, however, since I have a very good idea what the TS looks like, I think the dimer method should provide the best option. Just how can I make this more reliable?
I would appreciate your input on this problem, I am sure there are people with more expertise than I currently have that could give some advice. In particular I would be interested in
(i) Which parameters should I play around with? POTIM is certainly one of them, but so far changing this value in either direction has not changed much.
(ii) Do you normally preoptimize TS guess structures, relaxing only relevant atoms? Does this normally help to improve convergence?
(iii) Which part of the molecule do you include for the initial phonon calculation? I focus on 3-4 atoms, sometimes the whole substrate, depending on the problem. Normally this gives very nice imaginary modes that can be used as vector in the dimer method. Is it advisable to partially optimise these few atoms to the TS first, and then relax the full geometry to the TS? I want to add that in my case the remainder of the unit cell is normally fairly well converged to a minimum and doesn't do much else other than provide the crystal matrix.
Below is an INPUT file for the dimer method I use. Unfortunately I am unable to append any other file types such as structures so you could visualise what the problem is.
Your help is much appreciated.
Tobias
INCAR (dimer)
General:
SYSTEM = Rh-dcype-endo-nbd
! if enough K-points are present -5 should be used
ISMEAR = 0 ! 0: Gaussian -> semicond/insulator ; 1-N: MP; -5: Tetra+Blochl -> metal/very accurate energy/forces
SIGMA = 0.05
EDIFF = 1.0E-7
PREC = Normal ! to make ROPT=2.0E-4 for LREAL=Auto
LREAL = Auto ! cheaper for large cell's we are using
ENCUT = 600 !
ISYM = 2 # use symmetry as done for PAW PP
NELMIN = 8
LWAVE = .FALSE.
LCHARG = .FALSE.
! LVTOT = .TRUE. ! Write LOCPOT file for Toon
! LVHAR = .TRUE. ! but LOCPOT without exchange and correlation contribution for forcefields Toon
VOSKOWN = 1 ! important for GGA (PW91) for interpolation of XC...since we only use PBE=switch on (LDA: VOSKOWN=0)
LASPH = .TRUE. ! For VASP.5.X the aspherical contributions are properly accounted for in the Kohn-Sham potential
! as well. This is essential for accurate total energies and band structure calculations for f-elements
! (e.g. ceria), all 3d-elements (transition metal oxides), and magnetic atoms in the 2nd row (B-F atom),
! in particular if LDA+U or hybrid functionals or meta-GGAs are used, since these functionals often result
! in aspherical charge densities.
Van der Waals Interaction (vasp 5.3.3 patched verion):
IVDW = 12 ! switches between 0:off, 1: DFT-D3 and 2: TS-VDW (default=1)
dynamic:
IBRION = 44 ! -1: Fix atoms; 0: MD; 2: ConjGrad relax; 44: improved dimer method
NSW = 500 ! Number ionic steps
ISIF = 0 # relax ions, no change cell shape, no change cell volume (volume changes require ENCUT*1.3)
POTIM = 0.005
parallel:
LPLANE = .TRUE.
NPAR = 16