Author: Naveen Agrawal
Introduction: Nudged Elastic Band calculation is one of the methods to search for a transition state which considers a chain of images attached through elastic bands between reactant and product.[1] This chain should ideally represent a set of states along the reaction path after the NEB( Nudged Elastic Band) calculation has converged to the minimum energy path. However,we have to guess these set of images to initialize the search. One way to generate these first set of images is a linear interpolation of atomic coordinates between the reactant and product. This scheme provides a faster way to create intermediate states when the reaction itself does not involve a complex rearrangement of atoms, I am presenting a case of isomerization where using linear interpolation creates chemically unreasonably intermediate states and initializing the NEB calculation with these images would not lead to desired reaction trajectory. Specifically for the isomerization of HCN, the linear interpolation leads to unreasonable geometric configurations as described in figure 2 which would be unreasonably difficult or rather impossible to optimize to desired reaction path as described in figure 4.
Methods: DFT (Density Functional Theory) calculations were performed within the Vienna Ab initio simulation package (VASP, version 5.4.4), using the periodic supercell approach [2]. The projector augmented wave (PAW) method was used for electron-ion interactions [3,4]. The Perdew-Burke-Ernzerhof functional described the electron-electron exchange and correlation energies [5]. A plane wave basis set is used with an energy cutoff of 450 eV. For geometry optimization, the convergence criteria of the forces acting on atoms were 0.02 eV Å-1, while the energy threshold-defining self-consistency of the electron density was 10-5 eV. A cubic unit cell with 20A dimension was used for periodic supercell calculations. Considering that it is a unimolecular rearrangement and the electronic states should be discreet and localized therefore only gamma point was included for the Brillouin zone sampling. Zero Point Vibrational energy (ZPVE) corrections were included in the electronic energy calculated through VASP. Transition state search was performed using the Climbing Image Nudge Elastic band Method (CI-NEB). CI-NEB is an improved version of NEB method which improves the definition of reaction coordinate near the highest energy image.[6]
Geometry optimization of Isomers and linear interpolation: HCN and CNH molecules are optimized using periodic supercell approach as mentioned in the previous section. The optimized structures are linear and the distinction between the two structures is the position of hydrogen atom as attached to Carbon or Nitrogen. The optimized structures as obtained are shown in Figure 1.
Figure 1: Reactant and Product state for HCN to CNH isomerization
A linear interpolation of geometric co-ordinates between the reactant and product state resulted in a series of chemically unreasonable intermediate states as shown in Figure 2. Linear interpolation prompts the transfer of hydrogen from carbon to the nitrogen through the carbon-nitrogen bond itself which would basically optimize to dissociated atoms due to huge repulsion during the NEB calculation.
Figure 2: Reaction trajectory through linear interpolation from 00 to 05
Transition state search and confirmation: A more reasonable set of intermediate images are considered as presented in Figure 3 as an initial guess to the Climbing Image Nudge Elastic Band (CI-NEB) method. These images were created based on nature of the geometry for the transition state presented in the literature [7] allowing the reasonable transition of hydrogen through C-H bending.
Figure 3: Initial guess of intermediate images for NEB calculation
This sequence of intermediate states resulted into a converged NEB calculation with image 03 as the transition state as indicated by the change of sign of tangent forces along the reaction path and the converged normal forces perpendicular to the path for all the images to be less than 0.02 eV/A as can be seen in the attached snapshot of summary of VASP output attached in table below. Ideally the tangent forces should change sign across the transition state only once, however in the calculation summary the initial images show very small positive sign. This is possibly due to the close proximity of the image 01 and 02 to the reactant state, and the gradient of energy at such small energy differences (0.02 eV) is probably not accurate. However, it is still highly likely that indeed 03 is the transition state as it is the maxima along a reasonable reaction trajectory as indicated by figure 4 and table below. The relative energy of the highest energy image is 1.99 eV (activation barrier) without ZPVE corrections. I also increased the number of intermediate images to seven while gradually bending the C-H bond to allow better sampling of reaction path and it was possible to see a slightly better transition to the steep barrier of 1.99 eV as can be seen in the attached animation in appendix.
| Image number | Energy relative to initial state | Tangent forces | Normal forces |
|---|---|---|---|
| 0 | 0 | ||
| 1 | 0.022862 | 0.0004 | 0.008946 |
| 2 | 0.022567 | 0.00515 | 0.007189 |
| 3 | 1.992472 | -0.006065 | 0.011774 |
| 4 | 0.772263 | 1.023981 | 0.019798 |
| 5 | 0.616651 |
Figure 4: Converged NEB reaction trajectory
A transition state should only have one imaginary normal mode of vibration and should lead to the product or reactant state considered.To confirm that 03 indeed is a valid transition state, a vibrational analysis was performed using IBRION =5 in VASP (finite difference energy derivative with harmonic approximation).Vibrational analysis revealed a high frequency unstable normal mode (negative curvature of energy surface) for the 03 image as shown in the first animation below. There were two additional low frequency mode obtained from vibrational analysis which represented rotation of C-N bond and are very weakly associated with reaction coordinate of H shuttle through rotation of C-H bonds. Therefore, the activation barrier and the transition state geometry are precise within the convergence bounds of calculation mentioned.
Conclusion: HCN to CNH isomerization involves shuttling of hydrogen atom through bending of C-H bond in the transition state as can be seen in Figure 5 and Figure 4. Linear interpolation between the reactant and the product state provides very unreasonable guess for the reaction trajectory because of the linear geometry of both product and reactant state. Allowing the C-H bending in initial guess images provides a well converged NEB calculation with 1.84 eV of kinetic barrier with ZPVE corrections included. Previous theoretical studies are in agreement with the barrier calculated. [8]The presence of first order saddle point was confirmed with vibrational analysis.
References:
1. Henkelman, G., Jóhannesson, G., & Jónsson, H. (2002) (pp. 269-302). Springer, Dordrecht.
2. Kresse, G. and Furthmüller, J.Physical review B, 54(16), p.11169.
3. P. E. Blochl, Phys Rev B, 1994, 50, 17953-1797
4. G. Kresse.et.al Phys Rev B, 1999, 59, 1758-1775
5. J. P. Perdew. et.al , Phys Rev B, 1992, 46, 6671-6687.
6. Henkelman.et.al , The Journal of chemical physics, 2000,, 113(22), 9901-9904.
7. Koch.et.al., The Journal of Physical Chemistry C, 2007,111(41), 15026-15033.
8.Gong.et.al, The Journal of chemical physics,2005, 122(14), 144311.
Appendix:
