Application of reaction path concept in intramolecular proton transfer

PII: S0097-8485(97)00057-0

Computers Chem. Vol. 22, No. 1, pp. 13±20, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0097-8485/98 $19.00 + 0.00

Application of Reaction Path Concept in Intramolecular Proton Transfer Andrzej Jaworski Institute of Physics, Warsaw University of Technology, Koszykowa 75, 00-662 Warsaw, Poland (Received 19 July 1997; Accepted 4 August 1997) AbstractÐIntramolecular proton transfer in the ground state and the ®rst singlet excited state are investigated on ab initio level using 2(2'-hydroxyphenyl)-imidazoline as a model compound. The reaction-path concept is applied and new descriptions of the proton positions are proposed. The full calculations show surprising di€erences between the method using mass-weighted internals and mass-weighted Cartesian coordinates. In the simpli®ed ``static'' approach the results are dependent on the choice of the reaction coordinate. If the choice is made with care, the results comparable with full intrinsic reaction coordinate search are obtained. The concerted reaction-path calculations are justi®ed only in the vicinity of the transition state. # 1998 Elsevier Science Ltd. All rights reserved

1. INTRODUCTION

path (CRP) concept (Miller et al., 1988), based on the straight-line intramolecular coordinate between optimized geometries of the reactant and the product. In the case of the IPT reaction one assumes that all geometrical parameters of the molecule change uniformly with the motion of the proton. The second approach (Fukui, 1981) is based on the classical equation of motion. By assumption, the nuclei start to move with an in®nitesimal velocity. This quasistatic nuclear displacement, ``intrinsic motion'', represents the vibrationless-rotationless motion path of the reacting system. Solution of the resulting equation in the ``mass-weighted'' Cartesian coordinates is called ``intrinsic reaction coordinate'' (IRC) (Fukui, 1981), ``intrinsic reaction path'' (IRP) (Mezey, 1980), or ``minimum energy path'' (MEP) (Gonzalez and Schlegel, 1989). In practice, we start from the transition state and go down to the reactants and to the products. The CRP approach has been frequently used in the most advanced theoretical calculations (Domcke et al., 1993; Sobolewski and Domcke, 1993). In such an approach the geometry of the molecule has to be optimized only twice for the considered electronic state. On the other hand, application of the full IRP/MEP approach, especially for investigations of excited states in larger systems, is frequently out of reach with present computational resources. Therefore, the simpli®ed version of the approach based on constrained optimization procedure is frequently used (Sobolewski and Adamowicz, 1996). In this ``static'' approach, one coordinate is chosen as the reaction coordinate and the remaining ones are optimized at each step of the reaction. In the case of the IPT, the reaction coordinate should be connected with the proton

In our previous paper (Jaworski and DegoÂrski, 1995) the AM1 Hamiltonian (Dewar et al., 1985) has been employed for the investigations of the intramolecular proton transfer (IPT) in the ground state and in the ®rst singlet excited states. In this paper we use the same model compound, 2(2'hydroxyphenyl)-imidazoline (HPI), for the investigation of the phenomenon on the ab initio level. We examine in detail the applicability of two kinds of reaction-path concept using Hartree-Fock (HF) (Roothan, 1951) and CI Singles (CIS) (Foresman et al., 1992) methods and make some calculations by means of density functional theory (DFT method with B3LYP functional; Becke, 1993). The IPT can be considered as a kind of chemical reaction. The position of a mobile proton determines the tautomeric form of the molecule. If two stable minima (in the ground or excited states) exist, they can be regarded as a substrate and a product; the motion of the proton is the reaction. The most interesting IPT processes take place in the excited states. However, any excited state calculations are much more dicult than corresponding ones in the ground state. In these type of calculations one should ®rstly construct the potential energy surface (PES) and then do the molecular dynamics computations. (We tacitly assume that the Born-Oppenheimer approximation is satis®ed.) Presently, such a scheme can be used only in the case of extremely small molecules. Therefore, we are forced to use the reaction-path idea. There are two kinds of reaction-path concept used in theoretical investigations of potential energy surfaces. The ®rst one is the concerted reaction13

14

A. Jaworski

methods or between the simpli®ed/full IRC/MEP approaches, thus we will mainly use parameters (1) or (2). The ab initio calculations were performed with the use of the Gaussian 94 program package (Frish et al., 1995). Scheme 1.

position. We verify such obtained results* by the full IRP/IRC approach following the method of Gonzalez and Schlegel (1989, 1990). 2. METHOD 2.1. Computational Details The structure of both enol (I) and keto (II) tautomers and employed numbering of our model compound (Jaworski and DegoÂrski, 1995) are the same as in the previous paper (Scheme 1). It should be stressed that the very promising class of laser dyes contains this compound as its main constituent (DegoÂrski et al., 1992). We consider our compound to be planar (see Concluding Remarks), therefore having 37 degrees of freedom. The tautomers di€er mainly in the position of the mobile hydrogen. We employ standard (6-31G(d) and 6-31G(d,p) basis (Krishnan et al., 1980) without additional functions at the mobile hydrogen. We have checked that appending polarization function to hydrogens has negligible e€ect on the ground state results. Since the number of basis functions is the crucial factor in estimating the computer time we decided to use 6-31G(d) basis set. For description of the ``reaction progress'' we employ three kinds of parameters: (1) Normalized projection of the moving proton position on the straight line connecting donor and acceptor atoms (Jaworski and DegoÂrski, 1995). (2) The normalized O120C10H13 angle (for the atom numbering see Scheme 1). (3) The length of the mobile proton path calculated as the integral of dl = r(f)df, where r(f) is the distance from C1 nucleus to the mobile proton at the angle f. Unfortunately, the last parameter, although the most intuitive one, cannot be used either for the comparison between the di€erent computational *The results can be dependent on the choice of the reaction coordinate, as noted by Sobolewski and Adamowicz (1996) for a compound without preexisting hydrogen bond (HB). We will see that even in the system with HB, the proper choice of the reaction coordinate can be extremely important and the most obvious one is frequently not the best. {Resulting from computations. {The identical results obtained for the CH3+H2 4CH4+H reaction with use of massweighted Cartesian and mass-weighted internal coordinates are regarded by Gonzalez and Schlegel (1990) as a proof for the proper implementation of their algorithm.

3. RESULTS AND DISCUSSION In all cases when the double minima with well de®ned transition states exist (HF method for the ground state and CIS one for the ®rst singlet excited state), the full IRC calculations were performed using the method of Gonzalez and Schlegel (1989, 1990). The full IRC approach cannot be used for DFT calculations because no double minimum (so no transitions state) appears in the ground state. 3.1. Geometrical Parameters In the ``static'' simpli®ed version of the MEP approach the use of the reaction coordinate of the type X0H . . . Y seems to be quite natural (Sobolewski and Adamowicz, 1996) (X and Y denote here donor and acceptor heteroatoms). It seems that one should use the O0H distance on the enol side of the barrier and the H0N span on the keto side and verify that the exact position of the switch point has no in¯uence on the results. It is clear that geometry of molecule at the minima and the transition states can be found without any constrains; the assumptions of the planarity can be checked by the frequency calculations. Rigorous calculations for the ground state of our molecule show that the real{ situation is much more complicated. Checking the consistency of both types of reaction coordinate in the vicinity of the transition state can lead to incorrect results as well. The relatively small basis set (6-31G*) used on our calculations enables the calculations in the ground state point by point without any interpolations, especially in the most interesting regions. Apart from this, our calculations indicate di€erences in the results of the full IRP/MEP approach, depending on the Cartesian/internal coordinates. 3.1.1. Ground state Figure 1 shows the dependence of the donor±acceptor distance on the proton position (PP), de®ned according to item (1) of Section 2.1. It presents three variants of the full IRP/MEP approach and the ``static'' results. The di€erence between the method using mass-weighted internals (+ symbols in Fig. 1) and the method using massweighted Cartesian coordinates (. symbols) is really surprising.{ The path calculated in the massweighted (MW) Cartesians is the intrinsic reaction coordinate of Fukui (1981) and is considered here as a standard. It was impossible to go from TS to the enol or keto minima using any of the mass-weighted algorithms. The calculations in Cartesians broke down after some number of steps and the use of the MW internals leads to smaller and smaller changes in the PP while approaching the minima (the reaction

Reaction path concept in IPT

15

Fig. 1. The dependence of the donor±acceptor (O120N4) distance on the proton position in the S0 state (HF method with 6±31G* basis). Proton position variable is de®ned according to item (1) of Section 2.1. All symbols are de®ned in text.

coordinate changes are due to the movement of the remaining nuclei). In calculations with the non mass-weighted (NMW) internal coordinates (x symbols) we get the complete proton path from the TS to both minima. Figure 1 shows that for all PP on the left-hand side* of the TS{ the donor±acceptor distance calculated with the use of NMW coordinates and the ``static'' calculations with N±H distance used as the parameter (r symbols) are indistinguishable. That is also almost true for the right-hand side (keto tautomeric form) of the barrier. For that form the points resulting from the ``static'' approach with O± H distance used as the parameter (q symbols) lie much nearer to the true (standard) mass-weighted Cartesian results than those calculated with massweighted internals. However, the O±H parameterized ``static'' calculations give incorrect results on the enol side of the TS. We think that the bias is the result of the rigid constrains (short O±H distance). When much larger N±H span was kept constant on the enol side of the barrier, the constrains were much softer and resulting bias smaller. This interpretation holds also for the O±H distance used as the parameter on the keto side. In the ``central'' region (proton positions between 0.4 and 0.6) the points calculated using massweighted Cartesians show parabolic shape similar to the non mass-weighted results. By contrast, the calculations done with use of the mass-weighted internals give almost constant donor±acceptor dis*Enol tautomeric form. {At TS all paths have to cross.

tance. Moreover, there is a region (PP in the range 0.37±0.41) where the results of ``static'' calculations parameterized by N±H distance are nearer to the true Cartesian results than those resulting from using the mass-weighted internals. It is evident that all paths cross at the points corresponding to the energy minima ($ symbols). All approaches used here predict pronounced changes in the molecular geometry during IPT. Although the results of the IRP/MEP calculations in the MW internal coordinates suggest that the distance between donor and acceptor remains almost constant in the central part of the moving proton path, the IRC (calculations in the MW Cartesians) predicts some changes even in that region. Presented results (Fig. 1) show that the CPR approach could be valid only in the central part (proton position 0.4±0.6) of the moving proton path. It means that using this approach one should optimize the geometry not in the region of the energy minima, but at the points much nearer to the TS. Therefore, one cannot avoid the problems connected with the validity of the constrained optimization using this method. Figure 2 presents the dependence of the perpendicular proton coordinate (PPC) on the position of the moving proton. PPC is de®ned as the distance from the proton to the proton-acceptor line and is measured in angstroms. We choose to use PPC instead of the curvature vector (see Gonzalez and Schlegel, 1989; Page and McIver, 1988) because the former can be used in connection with all methods. During the proton transfer our molecule remains planarÐif the calculated points lie on the straight

16

A. Jaworski

Fig. 2. The dependence of perpendicular proton coordinate (see text) on the proton position for the proton transfer in the S0 state (HF method with 6±31G* basis).

line on the plot (Fig. 2), they should also form a line of similar shape in the physical (three dimensional) space. One can see that both IRP/MEP methods employing the MW coordinates predict similar, straight-line like proton movement for PP from 0.44 to 0.54 (that is very close to TS). The visible di€erences appear for the values of PP smaller than 0.4. This is the side e€ect of the mentioned earlier di€erences in the values of the donor±acceptor distances.* It is interesting that on the keto side of the potential barrier the ``static'' approach based on the calculations with OH distance used as a parameter (q symbols) gives PPC values indistinguishable from the results of the IRP search in the massweighted internal coordinates (+ symbols). The ``static'' approach with the same (OH distance) parameter fails completely on the enol side of the barrier. In the case when the NH distance is used the results (r symbols) are close to the results of the NMW IRP/MEP search ( symbols). The linear change in the PPC is observed only for proton positions between 0.44 and 0.54, thus the CPR approach can be valid only in this region. 3.1.2. Excited state As was already mentioned ab initio calculations for the excited states need much more computational resources than calculations for the ground state. This is why for the ®rst singlet excited state *These values are used for normalization during the calculation of the proton position (Jaworski and DegoÂrski, 1995). {It is shifted towards the enol form in the S1 state.

we employed larger steps in the IRP/MEP calculations and we used fever points in the ``static'' ones. The IRP method without the mass-weighting, applied in the ground state mainly for comparison with the ``static'' approach, was not used here. The dependence of the donor±acceptor distance on the PP in the ®rst singlet excited state is presented in Fig. 3. All symbols have the same meaning as in the ground state. Because of the much larger spans between the points obtained with the ``static'' approaches, for the sake of clarity we join points resulting from these calculations by dashed lines. The main factor which in¯uences all ®gures is the positions of the TS.{ Due to this all paths cross at PP smaller than 0.5. Both MW MEP methods predict the minimal donor±acceptor distances at that value of PP. Similarly as in the ground state, the di€erences between the results obtained with the use of Cartesians or internal coordinates are remarkable. In the central part of the plot almost constant donor±acceptor distance is predicted by the internal coordinate MW method (+ symbols). By contrast, the use of Cartesians (. symbols) gives the parabolic shape. The donor±acceptor distances are now very similar for both minima ($ symbols), whereas in the ground state (compare with Fig. 1) the O12±N4 distance was shorter by over 0.1 AÊ in the keto case. As in the ground state, the ``static'' approach with the N±H distance used as the parameter (r symbols) works rather satisfactory on the enol side of the TS, whereas the O±H parameterization (q symbols) works fairly well on the keto side. The com-

Reaction path concept in IPT

17

Fig. 3. The dependence of the donor±acceptor (O12±N4) distance on the proton position in the S1 state (CIS method with 6±31G* basis).

mon feature of the plots (Figs 1 and 3) are large (about 0.3 AÊ) variations in the O12±N4 distances. Although the plots of the perpendicular proton coordinate (for the de®nition see Section 3.1.1.) are

in the ground state and S1 state similar (Figs 2 and 4), there are some di€erences. The moving proton in the enol minimum is nearer to the donor±acceptor line in the S1 state than in the ground state. For

Fig. 4. The dependence of perpendicular proton coordinate on the proton position for the proton transfer in the S1 state (CIS/6±31G*).

18

A. Jaworski

Fig. 5. Energy pro®le of the proton transfer in the S0 state (HF/6±31G*). The zero energy corresponds to ÿ529.2367 a.u.

the keto minimum, H13 is considerably further away in the S1 state. In the S1 state the ``static'' approach with N±H distance used as the parameter does not work on the keto side, but is much closer to the proper IRC results on the enol side of the barrier. 3.2. Energy Pro®les The observations made for the geometrical parameters apply also for the energy pro®les. For the ground state (Fig. 5) the results of both IRP/MEP methods which use the mass-weighted coordinates are virtually the same for the values of PP greater than 0.42, whereas in the S1 state (Fig. 6) those methods give the same results in the 0.42±0.56 range of PP. The di€erences between the values calculated with the use of the internal MW coordinates (+ symbols) and the Cartesian MW ones (. symbols) arise from the di€erent donor±acceptor distances calculated in this region by the considered methods. The ``static'' approach with the N±H distance used as a parameter (r symbols) is appropriate on the enol side of the barrier for PP greater than 0.44. For the values of PP smaller than 0.42 in the ground state the results of that approach are the same as arising from the use of the IRP/MEP search without mass-weighting (x symbols) and differ from the standard IRC (Cartesians with MW) results as much as those resulting from IRP/MEP calculations without MW. The same is true also for PP greater than 0.6 for the ``static'' approach with O±H distance used as the parameter (q symbols). The ``static'' approach works relatively well on the keto side of the barrier in the ground state even

when we use N±H distance as a parameter (r symbols). However, on that side of the barrier the O±H parameterization (q symbols) works a little bit better. The ``static'' method parameterized by N±H distance fails completely on the keto side in the S1 state. Irrespective of employed procedure, in the ground state the HF method gives a curve with the asymmetric double minima (the enol minimum situated at 14 kcal/mol below the keto one). The minima are separated by the 18 kcal/mol energy barrier, shifted towards the keto form. The total distance between minima, measured along the bottom of the potential valley, is equal to 0.9 AÊ [all distances are calculated according to item (3) of Section 2.1]. In contrast, the ``static'' approach with DFT method gives asymmetric single minimum curve. The minimum is clearly on the enol side, but situated only about 0.35 AÊ from the midpoint. The in¯ection on the keto side of the curve suggests the second minimum. The distance between enol minimum and the footprint of a keto one is below 0.6 AÊ. As the full IRC approach cannot be used in connection with DFT method, we carried out a series of single-point DFT/6±311G(d,p) calculations at geometries generated by the full IRC/6±31G(d) search and on the optimized geometries of the HF/ 6±31G(d) minima (Fig. 7, x symbols). Additionally, the HF calculations with the triple split-valence 6± 311G(d,p) basis were made for those geometries (+ symbols). For comparison, we present the IRC Cartesian results in the standard 6±31G(d) basis (. symbols) as well as the results of the standard calculations at enol and keto S0 energy minima ($ symbols). One can see (Fig. 7) that the use of a bet-

Reaction path concept in IPT

19

Fig. 6. Energy pro®le of the proton transfer in the S1 state (CIS/6±31G*). The zero energy corresponds to ÿ529.0541 a.u.

ter basis does not qualitatively change the results obtained with the smaller basis. In the (®rst singlet) excited state the CIS (Foresman et al., 1992) method with the same (6± 31G*) basis set gives curve with asymmetric double minimum. In contrast to the ground state the keto minimum is more stable by 5 kcal/mol. The energy barrier is situated almost exactly at the midpoint,

with its height 9.5 kcal/mol and 14.5 kcal/mol over the enol and keto minimum, respectively. The total distance between the minima (measured along the bottom of the potential valley) is almost equal to 1 AÊ. The results presented above are obtained with the use of the three di€erent numerical approaches. The use of the CIS method for excited states is compar-

Fig. 7. Energy pro®les of the proton transfer in the S0 state. The zero energy corresponds to: ÿ529.2367 a.u. (HF/6±31G* calculations, . and $ symbols), ÿ529.3622 a.u. (HF/6±311G** calculations, + symbols) and ÿ532.6471 a.u. (B3LYP/6±311G** calculations, x symbols).

20

A. Jaworski

able with the HF one for the ground state. In this case the trend is clearÐin the ®rst singlet excited state the keto tautomer is more preferred than the ground state and the barrier height over the preferred minimum is lower. Both methods do not take into account the main part of the correlation energy. In contrast, the DFT method (even with the B3LYP functional) exaggerate the correlation e€ects which gives too strong hydrogen bonds and reduced heights of the barriers.* 4. CONCLUDING REMARKS It is well known that the CIS method, especially without di€use function added to the basis set, cannot give good values of the spectroscopic parameters. However, the main goal of this paper is to ®nd the proper description of the IPT process. As we mentioned in Section 2.1, we consider our molecule to be planar. In contrast to semiempirical results (Jaworski and DegoÂrski, 1995), both quantum-mechanical methods used here for the ground state lead to the conclusion that for all the positions of the mobile hydrogen the molecule remains planar. A slightly di€erent situation is in the ®rst singlet excited state. When the proton is between the minima, the CIS method predicts perfect planar conformationsÐthe hydrogen bond stabilizes the planar arrangement of the phenol and imidazole rings. However, it can be shown that for the proton in the potential minimum or nearer to the heteroatom slightly non-planar conformations are preferred. The comparison of the total energies between the planar and the fully optimized conformations gives the magnitude of the e€ect; for the enol tautomer we have 1.3  10ÿ4 a.u. (0.08 kcal/mol) and 1.2  10ÿ4 a.u. for the keto one. Thus, the nonplanarity of the molecules can be neglected. On the other hand it agrees with the conclusion, derived from the semiempirical calculations (Catalan et al., 1990; Jaworski and DegoÂrski, 1995), that the molecule is extremely soft. The calculations for the HPI have shown that the CRP approach should be used with great care for IPT systems. We get the extreme values of some parameters when the proton is half-way between the donor and acceptor atoms. Those values can di€er much from ones at the minima, thus the optimization should be made nearer the TS.{ Acknowledgements ÐThe author is grateful to the anonymous referee for useful comments and sugges-

*The DFT calculations done for the same molecule in S1 state give the double minimum potential curve with a very low barrier (A. Jaworski, unpublished results). {On the other hand: the subject of a constrained optimization is a mine®eld (Bone, 1995).

tions. The main part of the computations was performed at the Interdisciplinary Centre for Mathematical and Computational Modelling of Warsaw University. Some of the calculations were done with the Cray 6400 computer of the Warsaw University of Technology. REFERENCES Becke, A. D., (1993) Journal of Chemical Physics 98, 5648. Bone, G. A. (1995) Computational Chemistry List. Discussion on the Internet. Catalan, J., Fabero, F., Soledad Guijarro, M., Claramunt, R. M., Santa Maria, M. D., de la ConceptioÂn FocesFoces, M., Cano, F. H., Elguero, J. and Saste, R., (1990) Journal of the American Chemical Society 112, 747. DegoÂrski, A., Jaworski, A., Szerewicz, E., Tykarski, L. and Wlostowski, M., (1992) Elektronika XXXIII (7), 18. Dewar, M. J. S., Zoebisch, E. G., Healy, E. F. and Stewart, J. J. P., (1985) Journal of the American Chemical Society 107, 3902. Domcke, W., Sobolewski, A. L. and Woywod, C., (1993) Chemical Physics Letters 203, 220. Foresman, J. B., Head-Gordon, M., Pople, J. A. and Frish, M. J., (1992) Journal of Chemistry and Physics 96, 135. Frish, M. J., Trucks, G. W., Schlegel, H. B., Gill, P. M. W., Johnson, B. G., Robb, M. A., Cheeseman, J. R., Keith, T., Petersson, G. A., Montgomery, J. A., Raghavachari, K., Al-Laham, M. A., Zakrzewski, V. G., Ortiz, J. V., Foresman, J. B., Cioslowski, J., Stefanov, B. B., Nanayakkara, A., Challacombe, M., Peng, C. Y., Ayala, P. Y., Chen, W., Wong, M. W., Andres, J. L., Replogle, E. S., Gomperts, R., Martin, R. L., Fox, D. J., Binkley, J. S., Defrees, D. J., Baker, J., Stewart, J. J. P., Head-Gordon, M., Gonzalez, C. and Pople, J. A. (1995) Gaussian 94, Revision D.3. Gaussian, Inc., Pittsburgh PA. Fukui, K., (1981) Accounts of Chemistry Research 14, 363. Gonzalez, C. and Schlegel, H. B., (1989) Journal of Physical Chemistry 90, 2154. Gonzalez, C. and Schlegel, H. B., (1990) Journal of Physical Chemistry 94, 5523. Jaworski, A. and DegoÂrski, A., (1995) Computers and Chemistry 19, 189. Krishnan, R., Binkley, J. S., Seeger, R. and Pople, J. A., (1980) Journal of Chemical Physics 72, 650. Mezey, P. G., (1980) Theor. Chimica Acta 54, 95. Miller, W. H., Ruf, B. A. and Chang, Y-T., (1988) Journal of Chemical Physics 89, 6298. Page, M. and McIver, J. W., (1988) Journal of Chemical Physics 88, 922. Roothan, C. C., (1951) Reviews of Modern Physics 23, 69. Sobolewski, A. L. and Adamowicz, L., (1996) Journal of Physical Chemistry 100, 3933. Sobolewski, A. L. and Domcke, W., (1993) Chemical Physics Letters 211, 82.