A Visualization of the Intrinsic Reaction Coordinate and Global

ABSTRACT: A classical multidimensional scaling (CMDS) method is employed to visualize an intrinsic reaction coordinate (IRC) and a global reaction rou...
0 downloads 0 Views 2MB Size
Download PDF
Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Visualization of the Intrinsic Reaction Coordinate and Global Reaction Route Map by Classical Multidimensional Scaling Takuro Tsutsumi,† Yuriko Ono,† Zin Arai,‡ and Tetsuya Taketsugu*,† †

Department of Chemistry, Faculty of Science, Hokkaido University, Sapporo 060-0810, Japan Chubu University Academy of Emerging Sciences, Kasugai, Aichi 487-8501, Japan

Downloaded via KAOHSIUNG MEDICAL UNIV on July 31, 2018 at 11:48:06 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: A classical multidimensional scaling (CMDS) method is employed to visualize an intrinsic reaction coordinate (IRC) and a global reaction route map consisting of the equilibrium minima and transition state structures connected by the IRC network. As demonstrations, the method was applied to the IRCs of the intramolecular proton transfer in malonaldehyde and the SN2 reaction of OH− + CH3F → CH3OH + F−, which are both well described by two principal coordinates. Next, the method was applied to the global reaction route map of the Au5 cluster; the resulting map shows appropriate positions of five minima and 14 transition states in a reduced 2- or 3-dimensional coordinate space successfully.

1. INTRODUCTION In a quantum chemical study, the intrinsic reaction coordinate (IRC) has been used as a reference reaction pathway to describe an elementary reaction process.1,2 The IRC is defined as a minimum energy pathway on the potential energy surface, which connects reactant and product minima via a transition state (TS) geometry (corresponding to a first-order saddle point). Although the IRC itself is determined only based on a geometrical feature of the potential energy surface, it can be used as a starting point for investigating the reaction dynamics.3 Recently, Maeda, Ohno, and Morokuma developed an automated reaction-path search method, which leads to a concept of a global reaction route map consisting of the IRC network for a given molecular system.4 The global reaction route map gives all the information on the IRCs and the connectivity of the minima and TSs, and thus, one can discuss the unknown reaction mechanism for a given molecular system under a given reaction condition, considering possible complicated reaction routes. Once a global reaction map is obtained, it is possible to carry out kinetic simulations on the assumption of a transition state theory5 and to examine dynamical reaction routes on the global reaction route map by a trajectory on-the-fly molecular dynamics simulation.6 The global reaction route map contains vast of information on atomic coordinates of a lot of minima and TSs, the connectivity of stationary points, and the energetics, and thus, it is a hard task to grasp an entire picture of the connectivity of the reaction routes even after the global reaction route map is obtained. In previous studies, the minima and TSs in the global reaction route map were schematically illustrated to represent the connectivity, and there has been no attempt to make the reduced-dimensional reaction route map from a mathematical approach. © XXXX American Chemical Society

Nowadays, the reduction of multidimensional data to the smaller dimension is employed as a useful approach to save a huge amount of data and to visualize complex data, in a variety of research fields of natural and social sciences, e.g., analyses of big data, pattern recognition, genome analysis, psychology, market research, etc.; as standard methods, principal component analysis,7−9 multidimensional scaling (MDS),7−10 and selforganizing maps11,12 are well-known. In the MDS approach, multidimensional data are placed in a reduced dimensional space so that similar data are located close to each other, while the nonsimilar data are located far from each other.13,14 The MDS method that employs a linear distance (Euclidean distance) to measure the similarity of data is called a classical MDS (CMDS), or a principal coordinate analysis.7−10 By employing the CMDS approach for a given data set, the principal coordinates to represent the largest dispersion in geometrical structures can be chosen from all the degrees of freedom, based on the distances between each data. Very recently, the CMDS approach is applied to the classification of protein conformers and the analysis of molecular dynamics simulation.15,16 In this study, we propose a methodology to reduce the dimension of structural data in a global reaction route map using the CMDS approach. The method is first applied to the IRC pathway of the intramolecular proton transfer in malonaldehyde (C3O2H4) and the SN2 reaction of OH− + CH3F → CH3OH + F−, and then, it is applied to a set of structures of a small gold cluster, Au5, to visualize the global reaction route map in a reduced dimension. Received: February 17, 2018 Published: July 12, 2018 A

DOI: 10.1021/acs.jctc.8b00176 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

2. CLASSICAL MULTIDIMENSIONAL SCALING THEORY For a given N atomic molecular system, geometrical structures can be expressed by 3N Cartesian coordinates denoted as x = (x1, y1, z1, ..., zN) in which the center of mass is set as the origin of the x−y−z coordinate system. Here, we employ the massweighted Cartesian coordinates defined as ξ = (ξ1, ..., ξ3N) = ( m1 x1, ..., mN zN ) where mN is the Nth atomic mass; note that the Newton’s equation of motion becomes atomic-mass free in this coordinate system. For a pair of two geometrical structures, ξ(i) and ξ(j), a linear distance (Euclid distance) in mass-weighted coordinates, dij, can be expressed as

which can be used to measure the accuracy of the p-dimensional representation quantitatively. There is another index to measure the validity of a dimensional reduction, referred to as a strain function.9 The strain function approaches to zero when the reduced dimension covers the relative positions of n structures in a full-dimensional coordinate space.

3. RESULTS AND DISCUSSION 3.1. Applications to the IRC Pathway. 3.1.1. Intramolecular Proton Transfer in Malonaldehyde. The malonal-

3N

dij =

∑ (ξk(i) − ξk(j))2 k=1

(1)

The orientation of the x−y−z coordinate axes are determined to minimize dij for a pair of ξ(i) and ξ(j), by the Kabsch algorithm.17 When there are multiple structures in the coordinate space, one can evaluate dij for each pair and can make a distance matrix D with dij as the ijth element. D is of course a real symmetric matrix. Hence, we describe the way to reduce a dimension of n structures of an N atomic system by the CMDS algorithm.7−10 First, generate a distance matrix, D, for n structures and make a squared proximity matrix D(2) whose ijth element is dij2. Then, by employing the Young−Householder transformation,18 D(2) is transformed to the inner-product matrix, B, 1 B = − JD(2)Jt 2

Figure 1. Energy profile along the IRC for the intramolecular proton transfer in malonaldehyde.

(2)

where the centering matrix J is defined as J = E − n−11

(3)

Here, E is a unit matrix of n dimension, and 1 is a matrix of n dimension whose elements are all 1. Through a diagonalization of B in eq 2, one can obtain eigenvalues, {λ12, λ22, ..., λn2}, and the corresponding eigenvectors. The eigenvalues are arranged in descending order as λ12 ≥ ... ≥ λn2. The eigenvector with the largest eigenvalue, λ12, corresponds to a direction in the coordinate space in which a given set of n geometrical structures are most widely varied. By choosing the first p eigenvectors with eigenvalues, {λ12, λ22, ..., λp2}, as the coordinate axes, the positions of n structures can be expressed most effectively. The corresponding coordinates, {X1, ..., Xp}, are referred to as principal coordinates. When the CMDS analysis is applied to n structures of an N atomic system, the maximum number of positive eigenvalues is 3N − 3 that involves 3N − 6 internal degrees of freedom and three rotational degrees of the entire system, although the eigenvalues for three rotational modes should be almost zero. The number of relatively large eigenvalues determines the degree of freedom to represent the relative positions of n structures. The significance of the kth principal coordinate, Xk, can be measured by the proportion of variance defined as8 λk n ∑i λi

Figure 2. IRC profile of the intramolecular proton transfer in malonaldehyde in the 2-dimensional principal coordinate space. Relative energies are distinguished by a color change.

Figure 3. Energy profile along the IRC for the SN2 reaction of OH− + CH3F → CH3OH + F−.

(4)

The cumulated proportion can be defined as8 p

dehyde (C3O2H4) has an intramolecular hydrogen bonding which shows an intramolecular proton transfer via tunneling mechanism and has been used to test various semiclassical

∑k λk n ∑i λi

(5) B

DOI: 10.1021/acs.jctc.8b00176 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 6. Proportion of variance of 12 eigenvectors for 19 structures on the global reaction route map for Au5, obtained from the CMDS analysis.

Figure 4. IRC profile of the SN2 reaction of OH− + CH3F → CH3OH + F− in the 2-dimensional principal coordinate space. The additional points, Ai (i = 1−5), denote the structures where the C−F interatomic distance r is changed from 2.28 Å at the structure A to 2.50, 2.75, 3.00, 3.25, and 3.50 Å while the other geometrical parameters are fixed. Relative energies are distinguished by a color change.

other, indicating that the IRC pathway is sharply curved at two points before and after the transition state. Therefore, the onedimensional tunneling path model breaks down, and the multidimensional theory is required to reproduce the tunneling splitting value measured experimentally.19 The malonaldehyde has a planar structure with Cs symmetry which is conserved along the IRC for the intramolecular proton transfer, while only the TS structure has a more symmetric structure of C2v symmetry. We first optimized the TS geometry under a C2v symmetry restriction, followed by a normal-mode analysis to verify that the optimized TS structure has only one imaginary frequency mode. Then, an IRC calculation was carried out, starting from the TS. Figure 1 shows an energy profile along the IRC, with geometrical structures of the reactant, TS, and product. The energy variation shows a symmetric profile with respect to the TS, since the geometries of

tunneling theories.19 First, we calculated an IRC pathway for this intramolecular proton transfer process in malonaldehyde and employed the CMDS approach to extract significant coordinates to represent positions of several reference points along the IRC. To determine the IRC pathway, quantum chemical calculations were performed at the B3LYP/6-31G** level, using the GAMESS program package.20 The IRC pathway for this proton transfer involves three processes, i.e., (1) an approach of OH and O, (2) an H transfer between O and O, and (3) a departure of O and HO; also, the single and double bonds in bond alternation are exchanged as OH and O approach to (or depart from) each

Figure 5. All minima and TS structures for Au5, calculated by PBE/LanL2DZ.27 The point group, relative energy (kcal/mol), and principal coordinates (X1, X2, X3) (bohr amu1/2) are also given for each structure. C

DOI: 10.1021/acs.jctc.8b00176 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 7. Distributions of 19 structures of Au5 consisting of 5 MINs and 14 TSs in the global reaction route map, along the axis of the principal coordinates: (a) X1, (b) X2, and (c) X3.

reactant and product minima are symmetric with each other. As shown in Figure 1, 29 structures were located along the IRC, involving two minima and one TS. The malonaldehyde involves 9 atoms, with 21 internal degrees of freedom. Since the molecular structure conserves Cs symmetry along the IRC, 29 structures on the IRC can be expressed only by 15 totally symmetric coordinates of Cs symmetry. We applied the CMDS analysis to the 29 structures located along the IRC related to the intramolecular proton transfer in malonaldehyde. Through a diagonalization of the inner-product matrix B in eq 2, we obtained 15 positive eigenvalues, corresponding to the configuration space spanned by 15 totally symmetric coordinates of Cs symmetry in a planar geometry of malonaldehyde. Among 15 eigenvalues, the proportions of variance from the largest and second largest eigenvalues are 0.657 and 0.340, respectively, and thus, the cumulated proportion is 0.997 (almost one) by these two dimensions. The strain function for the reduction to one dimension (the largest eigenvalue) was calculated to be 0.0478, while this value is decreased to 0.0001 in a reduction to two dimensions. These results indicate that the IRC for the intramolecular proton transfer in malonaldehyde can be expressed by only two principal coordinates, X1 and X2. Figure 2 shows positions of the 29 points along the IRC in (X1, X2). The profile of the IRC pathway represents very well the feature of the reaction pathway for the intramolecular proton transfer where X1 corresponds to the atomic movement of proton transfer and exchange of bond alternation between CC single and double bonds, while X2

corresponds to the atomic movement of approach and departure of O and O atoms. The present result is consistent with the intuitive picture that the two degrees of freedom (H transfer and O···O approaching) contribute largely to the reaction process. A choice of the two reaction coordinates is of course possible by another approach, e.g., using the normal modes related to H transfer and O···O approaching, but such an approach requires a previous knowledge for the target reaction. It is noted that the present analysis requires no previous knowledge about the reaction pathway and automatically finds the significant principal coordinates based on the reference structures along the IRC pathway. Also, the evaluation index such as the cumulated proportion and strain function can give the information on the completeness of descriptions based on the reduced coordinate space. In the next section, we applied the present CMDS analysis to the more complicated reaction. 3.1.2. SN2 Reaction of OH− + CH3F → CH3OH + F−. The second application is an IRC pathway of the SN2 nucleophilic substitution reaction, OH− + CH3F → CH3OH + F− (involving 7 atoms with 15 degrees of freedom). According to the IRC calculation at the MP2/6-31+G* level,21 the product minimum of this SN2 reaction is a CH3OH···F− hydrogen-bonded complex which is 28.6 kcal/mol more stable than the dissociative products of CH3OH + F−. Hase and co-workers21 carried out ab initio direct dynamics simulations for this reaction, starting from the transition state toward the product minimum and found that 90% of trajectories lead to the dissociative products of CH3OH D

DOI: 10.1021/acs.jctc.8b00176 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

decreased to 0.0014 in a reduction to two dimensions. These results indicate that the IRC for the target SN2 reaction can be mostly expressed by two principal coordinates, X1 and X2. We also added five structures of Ai (i = 1−5), where the C−F interatomic distance r is changed from 2.28 Å at the structure A to 2.50, 2.75, 3.00, 3.25, and 3.50 Å while the other geometrical parameters are fixed those of A, and carried our CMDS analysis to the resultant 114 structures. Then, the proportions of variance from the largest and second largest eigenvalues become 0.801 and 0.146, respectively, and the cumulated proportion is 0.947; the strain function for the reduction to 1-dimension was calculated to be 0.0126, while this value is decreased to 0.0022 in a reduction to two dimensions. Figure 4 shows positions of the 114 reference points consisting of 109 IRC points and five points on the direct dissociation pathway (A1−A5) in the two principal coordinates (X1, X2). The IRC pathway is obtained as a continuous curve in this reduced dimensional space. It was verified that the IRC profile in Figure 4 is almost the same as that generated from only the 109 reference points on the IRC. The TS is located at the early stage of the entire IRC pathway, and in the region from R to the middle point of A and B, X2 mainly varies from negative to positive values, with a small change in X1; then, the IRC pathway changes its direction suddenly so that both X1 and X2 decrease linearly. This feature indicates that X1 is correlated with the O− C−F bond angle, while X2 is correlated with the C−F interatomic distance. Around the sharply curved region of the IRC pathway, the direct dissociation pathway to CH3OH and F− represented by A1−A5 appears as a branch, and the dynamics effect such as an inertial force (a centrifugal force) pushes the molecule off the IRC pathway toward this direct dissociation pathway effectively. 3.2. Application to the Global Reaction Route Map of Au5 Cluster. Gold is an inert material as bulk, but the gold nanoparticle with a small radius shows a catalytic activity; the activity increases very sharply as the radius of nanoparticle becomes smaller.25 The number of theoretical studies on the catalytic activity of the gold cluster increases more and more.26 Recently, we investigated a global reaction route map for the structural transformation of Au5 cluster and discussed features of the reaction-path bifurcation by surveying valley-ridge transition points along the IRC pathways.27 A methodology to examine the traces of dynamical trajectories based on a static reaction route network of IRCs was also proposed, taking the global reaction route map of Au5 as an example, and it was demonstrated that the dynamical trajectories do not always follow the IRC pathway and possibly jump to the other IRCs.6 Such dynamical effects beyond the static reaction pathway have been discussed by the Born−Oppenheimer MD approach.21 In this study, we applied the CMDS analysis to the stationary points on the global reaction route map of Au5, to reduce the dimensionality. It was previously shown that the Au5 cluster has five minima (referred to as MINi; i = 1−5) and 14 transition states (referred to as TSi-j which connects MINi and MINj) in the global reaction route map.27 Among 14 TSs, seven TSs are categorized as TSi−i that connect equivalent but distinct minima of MINi. Figure 5 shows the geometrical structures of these minima and TSs, point groups, and relative energies, as well as the principal coordinates, X1, X2, and X3, determined by the CMDS analysis. The global minimum, MIN1, is directly linked to nine TSs, consisting of four TSi−i (TS1−1a, TS1−1b, TS1−1c, and TS1−1d) and five TSi−j (TS1−2, TS1−3a, TS1− 3b, TS1−4, and TS1−5).

Figure 8. Five MINs (denoted by circle) and 14 TSs (denoted by square) in the global reaction route map for Au5, reduced to (a) two and (b) three dimensions, by the CMDS analysis. The relative energies are distinguished by a color change.

+ F− (direct dissociation) while only 10% of trajectories lead to the CH3OH···F− complex, indicating the significance of the nonIRC dynamical pathway in the actual reaction process.22,23 Such dynamics effects can be understood by considering the curvature of the IRC pathway; when the molecule proceeds along the curved pathway, the centrifugal force pushes the molecule off the pathway in a negative direction of the curvature vector.24 We carried out an IRC calculation for this SN2 reaction at the MP2/6-31+G* level using GAMESS.20 Figure 3 shows an energy profile along the IRC pathway, which consists of 109 geometries including the reactant minimum (R), transition state (TS), and product minimum (P). As shown in Figure 3, a geometrical change from R to A via TS corresponds to an intuitive SN2 reaction where OH− attacks the carbon of CH3F− from the backside of F−, and then, OH is bound to the CH3 via an inversion of the CH3 umbrella, with a departure of F− atom. In this reaction, however, the IRC continues; F− loosely bound to C atom moves to the position loosely bound to H atom of the CH3 side (structure B in Figure 3), then comes around behind CH3OH, and finally binds to the H atom of the OH. We applied the CMDS analysis to the 109 structures which are located on the IRC. The proportions of variance from the largest and second largest eigenvalues are 0.815 and 0.137, respectively, and thus, the cumulated proportion is 0.952. The strain function for the reduction to one dimension (the largest eigenvalue) was calculated to be 0.0033, while this value is E

DOI: 10.1021/acs.jctc.8b00176 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 9. Two-dimensional reduced global reaction route map for Au5. Five MINs (in red square) and 14 TSs (in blue circle) are connected by the IRC pathways. TSi−i that connects equivalent but distinct minima of MINi. The relative energy (kcal/mol) are given for each structure.

counterclockwise rotation as MIN1 → TS1−3a → TS3−3. The characteristic atomic movements from MIN2 to TS3−3 are denoted by a blue arrow in Figure 7b. The feature of geometrical change in X3 direction can be understood from the structures of MIN1, TS1−1b, TS1−3b, TS1−5, TS1−1c, and TS1−1d as shown in Figure 7c. These structures have a common tendency that Au5 is composed of a nearly equilateral triangle Au3 (denoted by a blue triangle) and two other Au atoms. As X3 increases, the two Au atoms change to a linear form and is oriented symmetrically to the Au3 triangle. Figure 8 shows (a) 2- and (b) 3-dimensional reduced global reaction route maps for Au5 which involve 19 structures. In the 2-dimensional map in Figure 8a, the global minimum, MIN1, is located almost at the origin, and the other four minima are separated from each other. Each minimum makes its own region with several TSs. It is clearly shown that TSi−i is located near MINi, while TSi−j is located between MINi and MINj. The location of the respective MINs and TSs is very informative, which coincides with an intuitive picture for the reaction route network. The origin and the region with negative X2 values near the X2 axis correspond to the low energy region. In the 3dimensional map, one can note that TS1−1d, TS5−5, and TS1−1c are separated from other structures which are located closely in the 2-dimensional map. Especially, TS1−1d and TS1− 1b are located very closely in Figure 8a, but they are separated far away from each other in Figure 8b. This feature can be understood by considering the meaning of X3 coordinate discussed above; in TS1−1b, one Au atom is bound to a vertex of the diamond shape of Au4, showing a large negative X3 value. Finally, we show a 2-dimensional reduced global reaction route map for Au5 in Figure 9 where the structures of five minima and 14 TSs are connected by the IRC pathways. In previous studies for Au5, the global reaction route map was illustrated by chemical intuition in a 2-dimensional figure, in which the similarity of geometrical structures was considered

The CMDS analysis was employed to determine principal coordinates for 19 structures (five MINs and 14 TSs) of Au5. The distance matrix D was generated for the set of these structures, and the inner-product matrix B in eq 2 was diagonalized. Figure 6 shows a proportion of variance of each eigenvector from the CMDS analysis. The internal degrees of freedom for Au5 is 3N − 6 = 9, and thus, 12 positive eigenvalues (including three rotational modes) were obtained. The largest proportion of variance is 0.466 for λ1, with the second 0.200 and the third 0.143, and then, the cumulated proportion for the first two eigenvalues is 0.666, while the one for the first three is 0.809. The strain functions for reduction to two and three dimensions are calculated to be 0.041 and 0.016, respectively. Thus, the representation of MINs and TSs placed in the reduced 2- or 3dimensional configurational space should work well to provide an intuitive picture of the positional relations of the respective structures. To give an insight into the meaning of each principal coordinate Xi, distributions of the reference structures are shown along the principal coordinates, (a) X1, (b) X2, and (c) X3, in Figure 7 where the structures of several MINs and TSs with relatively small values for |Xj| (j ≠ i) are given. As X1 increases from negative to positive values, geometrical structures change from quasi-linear (TS3−5) to 3-dimensional structures (MIN4, TS1−4, TS4−4) via 2-dimensional planar structures, as shown in Figure 7a. In other words, as X1 increases, the number of Au− Au bonds increases (4 for TS3−5, 7 for MIN1, and 9 for TS4− 4). It is noted that a small gold cluster has tendency to take a 2dimensional planar structure due to the relativistic effects, as shown in Figure 5. As for a change of geometries along the X2 axis, MIN2, TS1−2, MIN1, TS1−3a, and TS3−3 are shown in the order from a negative to positive value in Figure 7b. According to the structural transformation for these geometries, the center Au atom in MIN2 goes up to make a triangle in MIN1, and then, the upper right and left Au atoms show a F

DOI: 10.1021/acs.jctc.8b00176 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation only for the respective local regions.6,27 On the other hand, the reduced global reaction route map in Figure 9 is generated based on the distance functions for all the pair of nodes (structures) by the CMDS method, and thus, the validity of this map is assured by a mathematical foundation. The location of each structure and the connectivity of the respective structures can provide a wealth of information to the structural transformations of Au5. Very recently we proposed a new approach to analyze dynamical trajectories from on-the-fly molecular dynamics simulations based on the global reaction route map and applied it to Au5.6 In total, 200 trajectories were run starting from TS1−1d, and 198/ 200 trajectories first come across TS1−1b that is close to the valley-ridge transition points on the IRC between TS1−1d and MIN1. After passing TS1−1b, 86/198 trajectories directly reach MIN1, while 62/198 trajectories approach to TS1−3a before reaching MIN1.6 This dynamical feature is now understandable from the reduced global reaction route map shown in Figure 9. The combination of a dynamical trajectory analysis with a reduced global reaction route map definitely give the deeper insight and understanding to chemical reaction dynamics.

structural transformation of the Au5 cluster from the positions of MINs and TSs and their connectivity. The proposed analysis provides a general approach to visualize multiple multidimensional structures in a reduced dimensional space based on the distance matrix and will be used to understand a geometrical feature of the IRC and an entire picture of the global reaction route map.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Tetsuya Taketsugu: 0000-0002-1337-6694 Funding

T.Tsutsumi was supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT, Japan) through Program for Leading Graduate Schools (Hokkaido University “Ambitious Leader’s Program”). This work was supported by JSPS KAKENHI with Grant Number 16KT0047 and partly supported by MEXT as “Priority Issue on Post-K computer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use). A part of calculations was performed using the Research Center for Computational Science, Okazaki, Japan.

4. CONCLUSION To visualize the IRC pathway and global reaction route map, we employed the CMDS analysis to reduce the dimensionality of the multidimensional reference structures on the IRC and those in the global reaction route map which are important to understand the elementary reaction processes. The CMDS analysis can be used to choose a set of principal coordinates that show the largest (and second largest, etc.) distribution for a given set of geometrical structures, and the significance of each principal coordinate can be estimated by the proportion of variance. As an evaluation index to judge the validity of the representation in a reduced dimension, the strain function is also used. In the present CMDS analysis, the distance functions for each pair of structures are used based on the mass-weighted Cartesian coordinates derived from the Kabsh method. As the first application, the method was employed for reference structures on the IRC pathway for the intramolecular proton transfer in malonaldehyde. This reaction is categorized as a heavy−light−heavy mass-combination system in which the IRC shows a sharply curved region before and after the TS. The CMDS analysis verified this picture, and the IRC was wellrepresented by only two principal coordinates, i.e., the Htransfer and the approach of two O atoms. The second application is the analysis of the IRC pathway of the SN2 reaction, OH− + CH3F → CH3OH + F−, for which the previous dynamics studies suggested that non-IRC direct dissociation pathway becomes important by dynamics effects. The IRC profile in the two-dimensional reduced principal coordinate space clearly shows that the IRC pathway is sharply curved in the branching region of the IRC pathway and direct dissociation pathway in the product side, providing a very informative picture on this reaction. In the third application, the CMDS method was applied to 19 structures of Au5 in the global reaction route map. The cumulated proportion by the first three principal coordinates was calculated as 0.809, indicating that the relative positions of the target structures can be appropriately described by three dimensions. The analysis of the structural change along each principal coordinate provides an intuitive picture. By plotting 19 structures in a reduced 2- or 3-dimensional coordinate space, the global reaction route map was visualized in an appropriate manner, and we got an intuitive insight into the

Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Fukui, K. Formulation of the Reaction Coordinate. J. Phys. Chem. 1970, 74, 4161−4163. (2) Maeda, S.; Harabuchi, Y.; Ono, Y.; Taketsugu, T.; Morokuma, K. Intrinsic Reaction Coordinate: Calculation, Bifurcation, and Automated Search. Int. J. Quantum Chem. 2015, 115, 258−269. (3) Truhlar, D. G.; Gordon, M. S. From Force Fields to Dynamics: Classical and Quantal Paths. Science 1990, 249, 491−498. (4) Maeda, S.; Ohno, K.; Morokuma, K. Systematic Exploration of the Mechanism of Chemical Reactions: The Global Reaction Route Mapping (GRRM) Strategy Using the ADDF and AFIR Methods. Phys. Chem. Chem. Phys. 2013, 15, 3683−3701. (5) Sumiya, Y.; Nagahata, Y.; Komatsuzaki, T.; Taketsugu, T.; Maeda, S. Kinetic Analysis for the Multistep Profiles of Organic Reactions: Significance of the Conformational Entropy on the Rate Constants of the Claisen Rearrangement. J. Phys. Chem. A 2015, 119, 11641−11649. (6) Tsutsumi, T.; Harabuchi, Y.; Ono, Y.; Maeda, S.; Taketsugu, T. Analyses of Trajectory On-the-Fly Based on the Global Reaction Route Map. Phys. Chem. Chem. Phys. 2018, 20, 1364−1372. (7) Ingwer, B.; Groenen, P. J. F. Modern Multidimensional Scaling; Springer Series in Statistics; Springer New York: New York, NY, 2005. (8) Härdle, W. K.; Simar, L. Applied Multivariate Statistical Analysis; Springer Berlin Heidelberg: Berlin, Heidelberg, 2012. (9) Koch, I. Analysis of Multivariate and High-Dimensional Data; Cambridge University Press: Cambridge, 2013. (10) Torgerson, W. S. Multidimensional Scaling: I. Theory and Method. Psychometrica 1952, 17, 401−419. (11) Kohonen, T. Self-Organized Formation of Topologically Correct Feature Maps. Biol. Cybern. 1982, 43, 59−69. (12) Horton, D. E.; Johnson, N. C.; Singh, D.; Swain, D. L.; Rajaratnam, B.; Diffenbaugh, N. S. Contribution of Changes in Atmospheric Circulation Patterns to Extreme Temperature Trends. Nature 2015, 522, 465−469. (13) Agrafiotis, D. K.; Rassokhin, D. N.; Lobanov, V. S. Multidimensional Scaling and Visualization of Large Molecular Similarity Tables. J. Comput. Chem. 2001, 22, 488−500.

G

DOI: 10.1021/acs.jctc.8b00176 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation (14) Feher, M.; Schmidt, J. M. Metric and Multidimensional Scaling: Efficient Tools for Clustering Molecular Conformations. J. Chem. Inf. Comput. Sci. 2001, 41, 346−353. (15) Pisani, P.; Caporuscio, F.; Carlino, L.; Rastelli, G. Molecular Dynamics Simulations and Classical Multidimensional Scaling Unveil New Metastable States in the Conformational Landscape of CDK2. PLoS One 2016, 11, 1−23. (16) Li, X.; Xie, Y.; Hu, D.; Lan, Z. Analysis of the Geometrical Evolution in On-the-Fly Surface-Hopping Nonadiabatic Dynamics with Machine Learning Dimensionality Reduction Approaches: Classical Multidimensional Scaling and Isometric Feature Mapping. J. Chem. Theory Comput. 2017, 13, 4611−4623. (17) Kabsch, W. A Solution for the Best Rotation to Relate Two Sets of Vectors. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1976, 32, 922−923. (18) Young, G.; Householder, A. S. Discussion of a Set of Points in Terms of Their Mutual Distances. Psychometrika 1938, 3, 19−22. (19) Yagi, K.; Taketsugu, T.; Hirao, K. Generation of FullDimensional Potential Energy Surface of Intramolecular Hydrogen Atom Transfer in Malonaldehyde and Tunneling Dynamics. J. Chem. Phys. 2001, 115, 10647−10655. (20) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (21) Sun, L.; Song, K.; Hase, W. L. A SN2 Reaction That Avoids Its Deep Potential Energy Minimum. Science 2002, 296, 875−878. (22) Ma, X.; Hase, W. L. Perspective: Chemical Dynamics Simulations of Non-Statistical Reaction Dynamics. Philos. Trans. R. Soc., A 2017, 375, 20160204. (23) Pratihar, S.; Ma, X.; Homayoon, Z.; Barnes, G. L.; Hase, W. L. Direct Chemical Dynamics Simulations. J. Am. Chem. Soc. 2017, 139, 3570−3590. (24) Taketsugu, T.; Gordon, M. S. Dynamic Reaction Path Analysis Based on an Intrinsic Reaction Coordinate. J. Chem. Phys. 1995, 103, 10042−10049. (25) Haruta, M. Gold Rush. Nature 2005, 437, 1098−1099. (26) Gao, M.; Lyalin, A.; Takagi, M.; Maeda, S.; Taketsugu, T. Reactivity of Gold Clusters in the Regime of Structural Fluxionality. J. Phys. Chem. C 2015, 119, 11120−11130. (27) Harabuchi, Y.; Ono, Y.; Maeda, S.; Taketsugu, T. Analyses of Bifurcation of Reaction Pathways on a Global Reaction Route Map: A Case Study of Gold Cluster Au5. J. Chem. Phys. 2015, 143, 014301.

H

DOI: 10.1021/acs.jctc.8b00176 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX