delocalization of muonium in the glycine–K+ complex

Chemical Physics 440 (2014) 135–141

Contents lists available at ScienceDirect

Chemical Physics journal homepage: www.elsevier.com/locate/chemphys

Quantum localization/delocalization of muonium in the glycine–K+ complex Takehiro Yoshikawa, Tomohiro Honda, Toshiyuki Takayanagi ⇑ Department of Chemistry, Saitama University, 255 Shimo-Okubo, Sakura-ku, Saitama City, Saitama 338-8570, Japan

a r t i c l e

i n f o

Article history: Received 1 April 2014 In final form 18 June 2014 Available online 24 June 2014 Keywords: Nuclear quantum effect Intramolecular proton-transfer Muonium Path-integral simulation Isotope effect

a b s t r a c t Previous electronic structure studies have revealed that the glycine–K+ complex has a low-barrier intramolecular proton-transfer pathway between zwitterionic and neutral forms. We have theoretically calculated quantum molecular structures of this complex including the proton-transfer process using a path-integral molecular dynamics technique on an interpolated potential energy surface developed at the B3LYP level of theory. When the transferring proton is substituted by muon, it was found that the muonium atom showed a broad distribution around the proton(muon)-transfer transition state region between the neutral and zwitterionic structures due to extreme nuclear quantum effects of a very light particle although the distribution peak is slightly deviated from the transition state. The present study demonstrates that Mu can be employed to probe transition-state regions of potential energy surfaces of proton-transfer chemical reactions. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Muonium (Mu) is an ultra-light isotope of hydrogen, consisting of a positive muon l+ and an electron e, since the mass of l+ is about one ninth of the proton mass. Although the lifetime of Mu is relatively short (2.2 ls), this timescale is long enough for studying the muonium chemistry experimentally [1,2]. In fact, Mu has long been used to understand nuclear quantum effects in chemical reactions including tunneling and quantized vibrations and there have been reported many studies showing very large isotope effects so far [3,4]. In particular, due to its extremely light mass and very large vibrational amplitudes, Mu-containing molecules frequently show a significantly delocalized character in molecular structures compared to hydrogen- and deuteriumcontaining molecules [5–15]. Alternatively, it has been reported that Mu can be localized in a very different region of the potential energy surface compared to hydrogen and deuterium cases due to multi-dimensional effects. In fact, it has been reported that Mu is sometimes localized around the saddle point region of the potential energy surface due to anomalously large zero-point energies at both reactant and product regions [16–19]. In this work, from a quantum mechanical point of view, we discuss localization/ delocalization of Mu in the simplest glycine–metal cation complex having an intramolecular proton-transfer pathway between amino ⇑ Corresponding author. Fax: +81 48 858 3700. E-mail address: [email protected] (T. Takayanagi). http://dx.doi.org/10.1016/j.chemphys.2014.06.013 0301-0104/Ó 2014 Elsevier B.V. All rights reserved.

group and carboxylic group. In particular, we are interested in whether Mu is localized around the proton-transfer transition state region between the neutral and zwitterionic forms due to its quantum mechanical nature. There have been extensive studies on possible conformers of amino acid molecules in the past. It has been generally accepted that amino acids exist in their neutral form (with NH2 and COOH groups) in the gas phase while they are found as the zwitterionic form (with NH+3 and COO groups) in aqueous solution [20–28]. This fact indicates that polar water environment significantly affects the stability of these different conformers. It is also known that metal-cation binding to amino acid molecules significantly affects the relative stability of various amino acid conformers [29–42]. In particular, there have been many theoretical studies on structures and stability of the complexes between metal-cation and glycine [29–42], the simplest amino acid molecule. Among various metal ions, we here focus on the glycine–K+ complex since previous electronic structure studies have shown that energy levels of the neutral and zwitterionic forms are comparable (nearly thermoneutral) and that the proton-transfer barrier height is known to be relatively low of a few kcal/mol. In addition, it should be pointed out that the interaction between potassium cation and protein plays an especially important role in the selective transport phenomenon across cell membrane. Although many structures of the glycine–K+ complex have been extensively studied in the past, we have also characterized the potential energy surface profile of this system using the GRRM (Global Reaction Route Mapping)

136

T. Yoshikawa et al. / Chemical Physics 440 (2014) 135–141

computational algorithm [43–46] (which will be described below in detail), that can automatically find all possible local minima and transition states on a given multi-dimensional potential energy surface, to confirm the previous results [29–42]. Fig. 1(a) displays the schematic potential energy profile of the glycine–K+ complex showing eleven low-lying minimum structures and transition states between those minima obtained at the B3LYP/ 6-31+G(d) level of theory. It is seen that four structures (EQ0, EQ1, EQ2, EQ3) are more stable than other structures and have comparable relative energies at this B3LYP/6-31+G(d) level without zero-point correction. This result is in good agreement with

previous theoretical studies [29–42]. The EQ2 structure is stabilized by the presence of the five-membered ring and by attractive electrostatic interactions between K+ and negative charges on carboxylic oxygen and the amine nitrogen. On the other hand, EQ0, EQ1, and EQ3 are stabilized with attractive interactions between K+ and carboxylic oxygen. Notice that the EQ3 structure takes the zwitterion form. It is interesting to notice that there exist large barriers between EQ2 and the latter group (EQ0, EQ1, and EQ3). Fig. 1(b) shows the schematic energy diagram with zero-point correction where harmonic zero-point vibrational energies were calculated with all the mass of five hydrogen atoms being 1.008 amu. It is

Fig. 1. Schematic potential energy profile of the glycine–K+ complex at the B3LYP/6-31+G(d) level of theory: (a) bare potential without zero-point energy correction, (b) with zero-point energy correction using harmonic frequency analyses, and (c) with zero-point energy correction but for the muoniated glycine–K+ complex where one of the H atoms is substituted by Mu. Notice that all five hydrogen atom sites were examined and the most stable result is presented.

137

T. Yoshikawa et al. / Chemical Physics 440 (2014) 135–141

seen that the zero-point corrected diagram is not largely different from the diagram without zero-point correction (bare potential energy diagram shown in Fig. 1(a). On the other hand, the schematic energy diagram shown in Fig. 1(c) is obtained when one of the five hydrogen atoms is substituted by Mu whose mass is 0.114 amu (the most stable site is shown). It is interesting to notice that the obtained energy diagram shows a very different behavior; the proton(muon)-transfer transition state between EQ1 and EQ3 is a global minimum on the free-energy surface. This is simply because Mu-containing molecules generally have anomalously large quantized vibrational energies due to very small reduced mass. In addition to this, the vibrational zero-point energy is significantly reduced at the transition state due to the lack of the vibrational energy along the proton(muon)-transfer reaction coordinate mode (because of the imaginary frequency). Thus, it is suggested that Mu substitution can largely affect the free-energy profile of the glycine–K+ system. This phenomenon is also interesting since Mu atom may act as a probe of proton-transfer transition state structure. Needless to say, proton-transfer reactions are known to be ubiquitous in biochemical systems. Motivated by this simple prediction, we here discuss quantum localization/ delocalization of Mu in the glycine–K+ complex and whether Mu can be used to probe the transition-state region of proton-transfer reaction. We employ path-integral molecular dynamics (PIMD) simulations to obtain thermal equilibrium structures including nuclear quantum effects. 2. Computational method 2.1. GRRM calculations As was already mentioned in Introduction, a global feature of the potential energy surface for the glycine–K+ complex has been firstly explored using the GRRM computer code developed by Maeda and Ohno [43–46]. The GRRM method is based on the scaled hypersphere search algorithm that can automatically find all local minimum structures from a given initial equilibrium structure using a uphill-walking technique. With this technique many low-lying structures on the glycine–K+ potential energy surface have been automatically explored. It should be mentioned that the GRRM code can simultaneously optimize all the transitionstate structures between minima. We have chosen to use the B3LYP/6-31+G(d) level of theory from compromise between accuracy and computational time. All the B3LYP-level calculations were performed using the Gaussian 03 program [47]. The GRRM code finally found 162 minimum structures and 216 transition state structures but only eleven low-lying equilibrium structures (EQ0–EQ10 denoted with an energetical order) and the transition states between these minimum structures, that are important in the present study, are displayed in Fig. 1. The transition state structure is denoted as TSij for the EQi M EQj interconversion.

and z0 is a set of internal coordinate values optimized at a given value of si with respect to energy V0, V1, and V2 are the minimum potential energy value, first derivative vector and second derivative matrix at si. Notice that our approach is essentially similar to the reaction path Hamiltonian idea [48]. We found that the above quadratic functional form works well for the most of the internal coordinate; however, for the z5 and z13 coordinates, which are corresponding to the O2–K stretch coordinate and \C–O2–K angle, respectively, we have used different functional forms (only for z5 and z13 diagonal terms) as

" i

V ðz5 Þ ¼

ci1



ci2 1þ 0 z5  z5 ðsi Þ þ ci2

12



ci2 2 0 z5  z5 ðsi Þ þ ci2

6 #

and 2

V i ðz13 Þ ¼ ci1 f1  expðci2 ðz13  z013 ðsi ÞÞÞg þci3 ðz13



z013 ðsi ÞÞ expfci4 ðz13

 z013 ðsi ÞÞg

where cij is a fitting parameter. The failure of the quadratic representation is simply due to very weak non-covalent interaction between K+ and carboxylic group The accuracy of these functional forms are demonstrated in Fig. 3, where fitted potential energy curves are compared to the B3LYP results at the neutral, transition-state and zwitterionic stationary points. A total of 31 points were sampled in the OAH distance (s) range of 1.0–1.8 Å and then we have optimized other internal coordinates with a fixed value of the OAH distance at the B3LYP/631+G(d) level of theory, where this level was chosen from the compromise of accuracy and computational cost. The quality of the present B3LYP results is compared to results obtained by other electronic structure calculations in Table 1. The first and second derivatives were then obtained at an each point. The optimized internal coordinate values, minimum potential energy values, first and second derivative values and fitting parameters for z5 and z13 were finally interpolated using the standard spline-interpolation as a function of s. In order to check the accuracy of our simple scheme, two-dimensional contour plots as a function of the OAH and NAO distances are compared in Fig. 4. It is seen that our approach can reasonably reproduce the original B3LYP-level potential energy surface at least for the proton-transfer pathway between the neutral and zwitterionic forms of the glycine–K+ complex.

2.2. Potential energy function of proton-transfer process In order to perform fully converged PIMD simulations we have constructed a full-dimensional potential energy surface that can describe only a limited region around the intramolecular protontransfer reaction pathway. In our approach, a local potential energy surface at a certain reaction coordinate point si is expanded up to the second-order polynomial as

1 T V i ðzÞ ¼ V 0 ðsi Þ þ V1 ðsi Þðz  z0 ðsi ÞÞ þ ðz  z0 ðsi ÞÞ V2 ðsi Þðz  z0 ðsi ÞÞ; 2 where s is the reaction coordinate and taken to be the OAH distance in this work. z is a set of the internal coordinates defined in Fig. 2

Fig. 2. Definition of internal coordinates of the glycine–K+ complex used in potential energy surface construction.

138

T. Yoshikawa et al. / Chemical Physics 440 (2014) 135–141

Fig. 3. Potential energy curves along the internal coordinates, z5 and z13, (see Fig. 2 for their definition) at the neutral, proton-transfer transition-state, and zwitterionic structures. Symbols and solid lines correspond to the B3LYP/6-31+G(d)-level results and fitted curves (see text), respectively.

Table 1 Comparison of barrier height for forward and reverse proton transfer processes (neutral form M zwitterion form).

a b

Level of theory

Forward barrier

Reverse barrier

B3LYP/6-31+G(d)a MP2/6-311+G(2d, 2p)a MP2/6-311++G(2df, 2p)a CCSD(T)/6-311++G(2df, 2p)b CCSD(T)/6-311++G(3df, 3pd)b xB97XD/6-311++G(3df, 3pd)b

3.95 3.93 3.53 5.09 4.77 4.79

3.02 2.10 1.71 2.13 1.44 1.61

Geometry was optimized at the B3LYP/6-31+G(d) level. Geometry was optimized at the CCSD/6-311+G(d, p) level.

2.3. PIMD simulations Standard imaginary-time PIMD simulations have been carried out to obtain quantum and thermal equilibrium structures of the glycine–K+ complex and its Mu-substituted variant. The massive Nóse–Hoover chain technique combined with the velocity Verlet algorithm was employed to control the system temperature. The time increment for solving the equation of motion was set to Dt = 5 atomic unit (0.12 fs). The PIMD simulations were performed with 64–256 ring polymer beads up to t = 500–1000 ps (4  106– 8  106 steps) depending on the temperature. Although the potential energy surface is constructed with internal coordinates, PIMD simulations are performed using the standard Cartesian coordinate system. Numerical convergence was carefully checked and we had

Fig. 4. Plots of two-dimensional potential energy surfaces of the intramolecular proton-transfer process between neutral and zwitterionic forms of the glycine–K+ complex as a function of the OAH and NAO distances: (a) B3LYP/6-31+G(d)-level calculations and (b) the present interpolated potential energy surface. Notice that all other internal coordinates are fully optimized with respect to energy. Zero energy is defined as the energy level of the neutral minimum. Contour increment is set to 1 kcal/mol.

to employ 256 beads for obtaining fully converged results at T = 100 K for the Mu-substituted complex. The details of our computational code used in this study are also described in Refs. [49–51]. 3. Results Fig. 5(a) displays three-dimensional perspective plots of the nuclear density distributions for the glycine–K+ complex obtained from the PIMD simulation at T = 300 K. Notice that density distributions of hydrogen atoms and K+ ion are displayed using 3D surface plots but the positions of other atoms are displayed using ball models with their positions being at their averaged positions. It is seen that the nuclear distribution of hydrogen atoms indeed shows that the molecular structure mostly takes the neutral form although a small contribution of the proton-transferred structure (zwitterionic form) is seen. The obtained result is reasonable since the neutral form is energetically more stable than the zwitterionic form (see Figs. 1 and 4). Fig. 5(b) shows a similar plot to Fig. 5(a) but for the Mu-substituted glycine–K+ complex, where the carboxylic hydrogen atom is substituted by Mu. It is interesting to notice that Mu is broadly distributed around the midpoint of the N atom of amino group and the O atom of carboxylic group. Thus, H and Mu show very different nuclear distributions. It is interesting to

T. Yoshikawa et al. / Chemical Physics 440 (2014) 135–141

Fig. 5. Three-dimensional perspective plots of quantum nuclear density distributions in the all-hydrogen glycine–K+ complex (a) and its Mu-substituted variant (b) obtained from PIMD simulations at T = 300 K.

note that K+-ion distributions also show a broad feature in both allhydrogen and Mu-substituted glycine–K+ complexes, but the obtained distribution is slightly broader in the former complex. This is simply because, in the neutral form, K+ ion is more loosely bound to carboxylic group compared to the zwitterionic form. This trend can be seen in Fig. 3(b); the angular potential energy curve for the neutral structure is very flat around the potential minimum

139

compared to that for the zwitterionic form. In the case of the zwitterionic form, the negative charge is highly localized in the –CO2 moiety and K+-ion binding becomes therefore tight. The harmonic vibrational frequencies of C–O2–K bending for the neutral and zwitterionic forms were calculated to be 56 and 106 cm1, respectively, at the B3LYP/6-31+G(d) level, thus showing a similar trend to the PIMD simulation results. In order to further understand the obtained quantum nuclear density distributions at a more quantitative level, the nuclear distributions obtained from the PIMD simulations have been projected on the two-dimensional surface consisting of the OAX (X = H or Mu) and NAO internuclear distances. The distributions are integrated over all other coordinates. Fig. 6 shows the twodimensional contour plots of H and Mu distribution functions at three temperature values of T = 100, 200 and 300 K. Also plotted in Fig 6 are the two-dimensional potential energy surface and minimum energy path. The most interesting point seen in the obtained distributions is that Mu is significantly delocalized and that the densities around the transition state region are always large. In addition, one cannot see the Mu densities around the neutral minimum nor zwitterionic minimum at all three temperatures. The Mu density peak for the T = 100 K result can be seen at R(OAMu)  1.18 Å and R(NAO)  2.42 Å close to the transition state but does not exactly match the saddle point. It is found that the Mu density distribution becomes broader as temperature increases. In addition, the density peak position is slightly shifted toward to the neutral minimum. In contrast to the Mu results, the H densities are mostly localized around the neutral minimum (R(OAH)  1.00 Å and R(NAO)  2.53 Å) although the peak positions are slightly deviated from the potential minimum presumably due to vibrational anharmonicity. The nuclear density distributions projected on the two-dimensional surface as a function of the OAX distance and OAXAN angle (X = H or Mu) are presented in Supplementary materials.

Fig. 6. Two-dimensional nuclear density distributions as a function of the OAH(Mu) and NAO distances obtained from the PIMD simulations at T = 100, 200, and 300 K for the glycine–K+ complex and its Mu-substituted variant. The two-dimensional potential energy surface is also shown.

140

T. Yoshikawa et al. / Chemical Physics 440 (2014) 135–141

Fig. 7. Two-dimensional free-energy surfaces as a function of the OAH(Mu) and NAO distances obtained from the PIMD simulations at T = 100, 200, and 300 K for the glycine–K+ complex and its Mu-substituted variant. Contour increment is 0.5 kcal/mol. The two-dimensional potential energy surface is also shown.

It should be emphasized that, due to the distribution difference between Mu and H, the averaged OAN distance is also different between the H and Mu complexes. The averaged OAN distance at T = 100 K in the Mu-substituted glycine–K+ complex is calculated at 2.43 Å, which is about 0.1 Å shorter than that for the neutral or zwitterionic structure. This suggests that Mu acts like an electron in chemical bonding although the corresponding bonding energy in the present system is rather small (an order of 2 kcal/ mol, see below). This is very interesting since the bonding energy is positive despite that there is no potential energy minimum. This bonding phenomenon comes from the fact that large quantized vibrational zero-point energies at the neutral and zwitterion structures and the simultaneous reduction of zero-point energy at the transition state, as was already discussed in Fig. 1. We also notice that this type of bonding is similar to ‘‘vibrational bonding’’ phenomenon proposed in previous theoretical studies on reactive scattering calculations of heavy-light-heavy atom–diatom chemical reactions [52–54]. It should be more useful to present a free-energy surface, FT(q) along appropriate reaction coordinates q at a certain temperature of T for more quantitative discussion. The free-energy surface FT(q) corresponds to the reversible work to move the centerof-masses of quantum particle to q, and can be defined as FT(q) = kBT ln qT(q), where qT(q) is the nuclear distribution probability to be observed at a given set of centroid coordinates q obtained from the PIMD simulation at T. Fig. 7 shows twodimensional free-energy surfaces as a function of the same coordinates as Fig. 6. The contour value shown is the relative free-energy value in unit of kcal/mol, which is measured from the free-energy minimum point. In the case of the Mu-substituted glycine–K+ complex, the T = 100 K free-energy value at its minimum is more stable than that at the neutral structure position (R(OAMu) = 1.01 Å and R(NAO) = 2.56 Å) by 1.9 kcal/mol. It is found that this free-energy value decreases only slightly with an increase of temperature;

1.4 and 1.2 kcal/mol at T = 200 and 300 K, respectively. As mentioned above, these free-energy values roughly correspond to the bonding energy due to quantum localization of Mu. Although the free-energy minimum position is somewhat shifted from the proton-transfer transition-state position, the disagreement suggests the importance of vibrational anharmonicity. On the other hand, in the case of the all-hydrogen glycine–K+ complex, the obtained free-energy surface feature is close to that of the bare potential energy surface as shown in Fig. 7. 4. Discussion and conclusions In this work, using a quantum path-integral simulation technique, we have demonstrated that Mu can locate in a very different region from the minimum of the potential energy surface due to an extreme nuclear quantum effect. We have focused on the Mu distribution in the glycine–K+ complex which has been known to show low-energy intramolecular proton-transfer between the neutral and zwitterionic forms. Preliminary prediction based on the 0 K free-energy profile using harmonic vibrational frequencies indicated that, when the carboxylic hydrogen atom is substituted by Mu, the proton(muon)-transfer transition state becomes a global minimum of the free-energy surface due to the significant reduction of the quantized vibrational energy at the transition state. The present quantum PIMD simulations have revealed that this simple prediction is qualitatively correct but quantitatively incorrect. The Mu distribution peak position is slightly deviated from the transition state region, but it should be emphasized that Mu is distributed in a very different region from the potential energy minimum regions and this behavior is thus completely different from all-hydrogen cases. As a result of this anomalous distribution, the average distance between N in amino group and O in carboxylic group becomes somewhat short in the Mu-substituted complex. The present study also suggests that Mu can be used as

T. Yoshikawa et al. / Chemical Physics 440 (2014) 135–141

a useful probe for different potential energy surface regions (close to the transition-state region) far from potential energy minima. It should be, however, mentioned that quantum Mu localization is highly sensitive to the proton-transfer barrier height of the reaction system. In the case of the present glycine–K+ complex, the intramolecular proton-transfer barrier between the neutral and zwitterionic forms is relatively low (4 kcal/mol, see Fig. 1 and Table 1). If the proton-transfer barrier would be much higher, it is expected that Mu may be located around the original neutral or zwitterionic minimum although vibrational amplitudes should be still large. On the other hand, when the barrier height is very low, quantum localization can be observed even for hydrogen/proton. For example, Kühn and his coworkers [55] have shown that hydrogen/proton is significantly localized around the transitionstate region of the potential energy surface for the protonated ammonia cluster N2 Hþ 7 using reduced dimensionality quantum calculations, where the corresponding proton-transfer barrier has been predicted at 0.8–1.0 kcal/mol. It should be mentioned that reaction exothermicity (endothermicity) is also playing an essential role in quantum localization phenomena. The relative energies of the neutral and zwitterionic forms in the glycine–K+ complex are close to each other, as mentioned in Introduction. If one of these structures are much more stable than the other structure, the present quantum localization phenomenon cannot be seen, of course. Finally, it should be emphasized that the present theoretical calculations have been done within a so-called Born–Oppenheimer approximation. The accuracy of this approximation for Mu-containing molecules should be carefully checked using a more sophisticated theory such as multi-component theory [56–59] in the near future. Nevertheless, we would like to emphasize that Mu localization phenomenon around transition-state regions of proton-transfer reactions due to purely quantum mechanical nature of light nuclei is very interesting both from theoretical and experimental points of view. The present result suggests that Mu can be used to probe transition-state regions of proton-transfer reactions as another candidate of transition-state spectroscopy [60–62]. Recently, many Mu-containing species (reaction products) have been experimentally observed with high detection sensitivity due to the advances in muon beam intensity as well as in spectroscopic techniques [63–65] although detailed magnetic hyperfine coupling information is available only for muoniated radical species and thus cannot be obtained for the present diamagnetic glycine–K+ chemical system. Nevertheless, we hope that the present theoretical study would stimulate experimental studies. In addition, we will extend the current study to more general chemical systems including intermolecular proton(hydrogen)transfer. Conflict of interest There is no conflict of interest among authors. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.chemphys.2014. 06.013. References [1] D.C. Walker, Muon and Muonium Chemistry, Cambridge University Press, Cambridge, 2009.

141

[2] E. Roduner, in: A. Kohen, H.-H. Limbach (Eds.), Isotope effects in Chemistry and Biology, Taylor and Francis, CRC Press, Nov. 1, 2005. [3] D.C. Walker, J. Chem. Soc., Faraday Trans. 94 (1998) 1. [4] C.J. Rhodes, J. Chem. Soc. Perkin Trans. 2 (2002) 1379. [5] R.M. Macrae, I. Carmichael, J. Phys. Chem. A 105 (2001) 3641. [6] I. McKenzie, J.-C. Brodovitch, P.W. Percival, J. Am. Chem. Soc. 125 (2003) 11565. [7] I. McKenzie, J.-C. Brodovitch, K. Ghandi, B.M. McCollum, P.W. Percival, J. Phys. Chem. A 111 (2007) 10625. [8] B.M. McCollum, T. Abe, J.-C. Brodovitch, J.A.C. Clyburne, T. Iwamoto, M. Kira, P.W. Percival, R. West, Angew. Chem. Int. Ed. 47 (2008) 9772. [9] B.M. McCollum, J.-C. Brodovitch, J.A.C. Clyburne, A. Mitra, P.W. Percival, A. Tomasik, R. West, Chem. Eur. J. 15 (2009) 8409. [10] R. West, P.W. Percival, Dalton Trans. 39 (2010) 9209. [11] D.G. Fleming, M.D. Bridges, D. Arseneau, J. Phys. Chem. A 115 (2011) 2778. [12] M.C. Böhm, R. Ramírez, J. Schultec, Mol. Phys. 103 (2005) 2407. [13] R. Ramírez, J. Schultec, M.C. Böhm, Chem. Phys. Lett. 402 (2005) 346. [14] R.M. Valladares, A.J. Fisher, W. Hayes, Chem. Phys. Lett. 242 (1995) 1. [15] B.S. Hudson, S.K. Chafetz, Molecules 18 (2013) 4906. [16] T. Tiyake, T. Ogitsu, S. Tsuneyuki, Phys. Rev. Lett. 81 (1998) 1873. [17] S. Tsuneyuki, Curr. Opin. Solid State Mater. Sci. 6 (2002) 147. [18] T. Takayanagi, Chem. Phys. 334 (2007) 109. [19] T. Yoshikawa, S. Sugawara, T. Takayanagi, M. Shiga, M. Tachikawa, Chem. Phys. 394 (2012) 46. [20] J.H. Jensen, M.S. Gordon, J. Am. Chem. Soc. 117 (1995) 8159. [21] Y. Ding, K. Krogh-Jespersen, J. Comput. Chem. 17 (1996) 338. [22] L. Gontrani, B. Mennucci, J. Tomasi, J. Mol. Struct. (THEOCHEM) 500 (2000) 113. [23] E. Kassab, J. Langlet, E. Evleth, Y. Akacem, J. Mol. Struct. (THEOCHEM) 531 (2000) 267. [24] A. Fernández-Ramos, Z. Smedarchina, W. Siebrand, M.Z. Zgierski, J. Chem. Phys. 113 (2000) 9714. [25] W. Wang, X. Pu, W. Zheng, N.-B. Wong, A. Tian, J. Mol. Struct. (THEOCHEM) 626 (2003) 127. [26] B. Balta, V. Aviyente, J. Comput. Chem. 25 (2004) 690. [27] C.M. Aikens, M.S. Gordon, J. Am. Chem. Soc. 128 (2006) 12835. [28] S.M. Bachrach, J. Phys. Chem. A 112 (2008) 3722. [29] S. Hoyau, G. Ohanessian, J. Am. Chem. Soc. 119 (1997) 2016. [30] S. Hoyau, G. Ohanessian, Chem. Eur. J. 8 (1998) 1561. [31] J. Bertrán, L. Rodíguez-Santiago, M. Sodupe, J. Phys. Chem. B 103 (1999) 2310. [32] F. Rogalewicz, G. Ohanessian, N. Gresh, J. Comput. Chem. 11 (2000) 963. [33] T. Wyttenbach, M. Witt, M.T. Bower, J. Am. Chem. Soc. 122 (2000) 3458. [34] T. Shoeib, C.F. Rodriquez, K.W. Michael Siu, A.C. Hopkinson, Phys. Chem. Chem. Phys. 3 (2001) 853. [35] C.H. Wong, F.M. Siu, N.L. Ma, C.W. Tsang, J. Mol. Struct. (THEOCHEM) 588 (2002) 9. [36] R.M. Moision, P.B. Armentrout, J. Phys. Chem. A 106 (2002) 10350. [37] E. Constantino, L. Rodíguez-Santiago, M. Sodupe, J. Tortajada, J. Phys. Chem. A 109 (2005) 224. [38] I. Corral, O. Mó, M. Yáñez, J.-Y. Salpin, J. Tortajada, D. Moran, L. Radom, Chem. Eur. J. 12 (2006) 6787. [39] M.H. Khodabandeh, M.D. Davari, M. Zahedi, G. Ohanessian, Int. J. Mass Spectrom. 291 (2010) 73. [40] T. Marino, N. Russo, M. Toscano, J. Inorg. Biochem. 79 (2000) 179. [41] T. Marino, N. Russo, M. Toscano, Inorg. Chem. 40 (2001) 6439. [42] N. Vyas, A.K. Ojha, Int. J. Quantum Chem. 112 (2012) 1526. [43] K. Ohno, S. Maeda, Chem. Phys. Lett. 384 (2004) 277. [44] S. Maeda, K. Ohno, J. Phys. Chem. A 109 (2005) 5742. [45] K. Ohno, S. Maeda, J. Phys. Chem. A 110 (2006) 8933. [46] S. Maeda, K. Ohno, J. Phys. Chem. A 111 (2007) 4527. [47] M.J. Frisch et al., Gaussian 03 (Revision D.02), Gaussian Inc., Pittsburgh, PA, 2003. [48] W.H. Miller, N.C. Handy, J.E. Adams, J. Chem. Phys. 72 (1980) 99. [49] M. Shiga, M. Tachikawa, S. Miura, Chem. Phys. Lett. 332 (2000) 396. [50] M. Shiga, M. Tachikawa, S. Miura, J. Chem. Phys. 115 (2001) 9149. [51] M. Shiga, M. Tachikawa, Chem. Phys. Lett. 374 (2003) 229. [52] J. Manz, Chem. Phys. Lett. 96 (1983) 607. [53] D.C. Clary, J.N.L. Connor, J. Phys. Chem. 88 (1984) 2758. [54] J. Manz, R. Meyer, H.H.R. Schor, J. Chem. Phys. 80 (1984) 1562. [55] Y. Yang, O. Kühn, G. Santambrogio, D.J. Goebbert, K.R. Asmis, J. Chem. Phys. 129 (2008) 224302. [56] M. Tachikawa, K. Mori, H. Nakai, K. Iguchi, Chem. Phys. Lett. 290 (1998) 437. [57] M. Tachikawa, Chem. Phys. Lett. 360 (2002) 494. [58] M.F. Shibl, M. Tachikawa, O. Kühn, Phys. Chem. Chem. Phys. 7 (2005) 1368. [59] A.D. Bochevarov, E.F. Valeev, C.D. SheRrill, Mol. Phys. 102 (2004) 111. [60] D.M. Neumark, Annu. Rev. Phys. Chem. 43 (1992) 153. [61] K. Prozument, R.G. Shaver, M.A. Ciuba, J.S. Muenter, G.B. Park, J.F. Stanton, H. Guo, B.M. Wong, D.S. Perry, R.W. Field, Faraday Disc. 163 (2013) 33. [62] R.M. Bowman, M. Dantus, A.H. Zewail, Chem. Phys. Lett. 589 (2013) 42. [63] C.J. Rhodes, Sci. Prog. 95 (2012) 101. [64] I. McKenzie, Annu. Rep. Prog. Chem. Sect. C 109 (2013) 65. [65] N.J. Clayden, Phys. Scr. 88 (2013) 068507.