On Benchmarking of Automated Methods for Performing Exhaustive

(4) Their ultimate goal was to cover reaction paths comprehensively by ..... Quality of ab initio and density functional theory (DFT) calculations sho...
1 downloads 0 Views 767KB Size
Download PDF
Letter pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

On Benchmarking of Automated Methods for Performing Exhaustive Reaction Path Search Satoshi Maeda*,†,‡,§ and Yu Harabuchi†,∥ †

Department of Chemistry, Faculty of Science, Hokkaido University, Kita 10, Nishi 8, Kita-ku, Sapporo 060-0810, Japan Institute for Chemical Reaction Design and Discovery (WPI-ICReDD), Hokkaido University, Sapporo 001-0021, Japan § Research and Services Division of Materials Data and Integrated System (MaDIS), National Institute for Materials Science (NIMS), Tsukuba 305-0044, Japan ∥ JST, PRESTO, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan J. Chem. Theory Comput. Downloaded from pubs.acs.org by MURDOCH UNIV on 03/18/19. For personal use only.

‡

S Supporting Information *

ABSTRACT: In recent years, the importance of computational chemistry approaches has grown rapidly because of recent advances in computational software and hardware. Automated reaction path search is one of the promising techniques which would allow exploration of unknown chemistry using computers. Several methods have been developed so far, and comparison of the performance of existing methods would be a subject of interest to many chemists. In this paper, we present the performance of our single-component artificial force induced reaction method (SC-AFIR) to rectify the result shown in a recent report [Grambow et al. J. Am. Chem. Soc. 2018, 140, 1035−1048] on the performance of the SC-AFIR method compared to four other methods. We discuss the flaws in their benchmark procedure and our thoughts on benchmarking of automated reaction path search methods.

I

KHP, and the adopted computational level is B3LYP/631+G(d).4 In these calculations, all geometrical displacements (including TS structure optimization) were treated by the GRRM17 program based on energy and gradient computed by the Gaussian 16 program.11 When trying to execute an exhaustive search of the paths as aimed by Grambow et al., it is required to set the model collision energy parameter γ extremely high as has been shown in our previous reports.10,12 In our past research, we have set γ = 1000 kJ/mol or higher to increase the comprehensiveness. Grambow et al., however, set γ as low as 400 kJ/mol. In addition, when the target is the first step of the unimolecular thermal isomerization and decomposition reactions of KHP, the NoBondRearrange option should be used. In the study of Grambow et al., however, the number of paths they tested was only 74 paths, indicating that they used a different option called FirstOnly. In other words, they performed calculations using the combination of γ = 400 kJ/mol and FirstOnly options, which was not suitable for this purpose, and using these options SC-AFIR only found TSs for 7 distinct product channels, significantly fewer than the number Grambow et al. found for the same system using other methods.4 Grambow et al. result for this particular case with these particular options might cause some readers to mistakenly think that the SC-AFIR method is incapable of finding more TSs.

n recent years, several groups have been working on the development of a computational code for automatically searching reaction paths while performing quantum chemical calculations.1,2 Automated reaction path search has attracted attention as a computational tool that accelerates prediction and design of new reactions as well as a conventional mechanism study.3 Recently, Grambow et al. studied performances of five different methods by applying them to a unimolecular reaction of the simplest γ-ketohydroperoxide (KHP), 3-hydroperoxypropanal.4 Their ultimate goal was to cover reaction paths comprehensively by integrating paths obtained through these five methods, i.e., freezing string method (FSM),5 growing string method (GSM),6 single-ended growing string method (SSM),7 single-component artificial force induced reaction method (SC-AFIR),8 and KinBot.9 For four out of five methods they used, the developers of those methods worked as co-authors or gave advice as stated in the acknowledgment, supporting to provide appropriate ways of their application. On the other hand, as for the SC-AFIR method we developed, they reported the result based on inappropriate usage of the method, thus leading to its incorrect interpretation. Below, for the same reaction system they considered, we demonstrate the correct performance of the SC-AFIR method. In addition, we would like to take this opportunity to write our thoughts on how to evaluate automated reaction path search methods. First, we show the results that can be obtained by the SCAFIR method implemented in the GRRM17 program when it is properly used.10 The target is the first step of the thermal unimolecular isomerization and decomposition reactions of © XXXX American Chemical Society

Received: November 22, 2018 Published: March 12, 2019 A

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

Letter

Journal of Chemical Theory and Computation

Table 1. 30 Lowest Barrier Channels and Relative Energy Values in kcal/mol of the Best (Energetically Most Preferable) TSsc

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

[H]OC1([H])OOC([H])([H])C1([H])[H] [H]C(O)C([H])([H])[H].[H]OC([H])O [H]OOC([H])([H])C([H])C([H])O[H] [H]OC([H])C([H])[H].[H]C([H])OO [H]C([H])C = O.[H]C([H])O.[H]O[H] [H]C(O)C([H])([H])C([H])O.[H]O[H] [H]C(O)C([H])([H])C([H])([H])O ([H])O [H]C(O)C([H])C([H])[H].[H]OO[H] [H]OOC([H])O.[H]C([H])C([H])[H] [H]C(O)C([H])C([H])[H].[H]O([H])O [H]OC1([H])O(O)C([H])([H])C1([H])[H] [H]C1([H])OC(O)C1([H])[H].[H]O[H] [H]OC([H])C([H])C([H])O.[H]O[H] [H]OC([H])(O[H])C([H])([H])C([H])O [H]OOC([H])C([H])[H].[H]C([H])O [H]C1([H])OC1([H])[H].[H]OC([H])O [H]C1([H])OOC(O)C1([H])[H].[H][H] [H]C(O)C([H])([H])[H].[H]C([H])OO [H]OC([H])([H])OC([H])([H])C([H])O [H]C(O)C([H])([H])C([H])OO.[H][H] [H]C(O)C([H])([H])[H].[H]C = O.[H]O [H]OC([H])C([H])[H].[H]O[H].CO [H]OOC([H])([H])OC([H])C([H])[H] [H]OC([H])([H])C([H])([H])C([H])OO [H]C([H])C([H])[H].[H]OO[H].CO [H]OOC([H])([H])[H].[H]C([H])CO [H]OOC1([H])OC([H])([H])C1([H])[H] [H]C(O)C([H])([H])[H].[H]COO[H] [H]OC([H])C([H])[H].[H]OC([H])O [H]OC([H])C([H])[H].[H]C = O.[H]O

SCAFIR1

SCAFIR1+

SCAFIR2

34.2 37.7 43.8

34.2 37.7 43.8 48.1 49.8 50.1 50.7

34.2 37.7 43.8 48.8 49.8 50.1 50.2

53.1

53.1 54.8 55.4 55.4 57.6 59.6 64.7 62.7 63.4 64.3 64.4 65.6 66.5 67.2 67.3 67.6 68.6 69.2 69.6 70.1 70.2 71.7 73.3

49.8 50.1 50.2 53.1 54.8 55.4

55.4

57.6 59.6 64.7 62.7 63.4

57.6 59.6 62.5 62.7

64.4 65.6 66.5 67.2 67.3 67.6 68.6 69.2

65.9 65.6 66.8 67.2 68.4 68.5

70.1 76.4 71.7

FSM

GSM

SSM

KinBot1

KinBot2

attributea

39.1 38.8 68.8

35.3 38.8 45.8

39.1 38.8 45.8

50.5 51.0 51.7

50.5 53.1 53.0

50.5 51.0 51.4

34.9 48.8 70.6 48.7 50.5 52.6 52.7

34.9 48.8 70.6 48.7 50.5 52.6 52.7

R37, R38 R69 R47, R53 R17 R10,b R28 R1−R5 R18, R19, R26

55.9 56.6

56.5 58.5

56.8 56.6

55.4 56.6

55.4 56.6

R74 R20,b R71, R72

76.0 61.6 65.8 69.9

58.1

58.1 61.6 63.6 66.4

58.3 60.3

58.3 60.3

63.4

63.4

65.1

65.9 65.1

69.5 68.6

69.5 68.6

R35, R36 R63 R54, R57 R67, R68, R73 R43 R39 R16 R52 R15 R14b

78.8

R56

69.9 80.4

R30 R65, R66 R44

66.8 65.6 66.8 69.9 71.0

69.7

63.6 72.8 65.5 65.4 65.6 67.2 67.4 75.2

66.4 67.1 66.8 70.0 69.4

78.2

68.7

69.7 89.6

70.3 71.6

69.9 80.4

a

Corresponding labels are defined in ref 4. bThe IRC calculation reached a product minimum different from the label in ref 4 after it completely terminated. cThe complete list of discovered channels is shown in the Supporting Information, Table SI1.

surface with weak artificial force of γ = 10/(N − 1) kJ/mol added between all pairs of atoms. Hence, a search with this option was also performed (the calculation result obtained with added weak force is referred to as SC-AFIR1+). In SCAFIR1+, 121 unique TSs directly connected between KHP and non-KHP structures and 50 TSs of conformational change of KHP were found after 346692 gradient calls. All the coordinates of minimum and TS geometries obtained by the SC-AFIR1+ search are shown in the Supporting Information, Table SI10 and Table SI11, respectively. Although many other TSs for second or later steps were also found during the search, we do not discuss them. Now let us compare the above calculation results with the data of ref 4. Because so many TSs were obtained, we organize and discuss each reaction channel. That is, SMILES was used to convert the molecular structure into a character string, and the structures having the same character strings were regarded as identical. Table 1 shows the character strings representing the products of the 30 lowest barrier channels. The full list of all 85 channels is shown in the Supporting Information, Table SI1. For 75 TSs tested by Grambow et al., the corresponding labels (R1−R75) are indicated in the table. In their report, there were several cases where the IRC calculation was terminated before reaching the local minimum structure completely when the system decomposed (R10, R12, R13, R14, R20, R31, and R62). This study continued the IRC

To prove the capability of SC-AFIR, we show the result obtained by combining γ = 1000 kJ/mol and NoBondRearrange. This search obtained 209 unique TSs (TS at the first step of reaction) directly connected between KHP and nonKHP structures. At the same time, 72 TSs of conformational change of KHP were obtained. All the coordinates of minimum and TS geometries obtained by the search are shown in the Supporting Information, Table SI6 and Table SI7, respectively. This calculation required 510960 gradient calls (the number of gradient calls required in the IRC calculations was subtracted to match the report by Grambow et al.). In this and the other calculations shown below, a Hessian matrix was never computed except during the IRC calculations. If a user intends to search more comprehensively, the SC-AFIRn (n > 1) calculation should be performed, which is also available in the GRRM17 program.10 The above standard SC-AFIR is denoted as SC-AFIR1 hereafter. The result of SC-AFIR2 calculation with γ = 1000 kJ/mol and NoBondRearrange is also shown: 401 unique TSs directly connected between KHP and nonKHP structures and 99 TSs of conformational change of KHP were found after 2033336 gradient calls. All the coordinates of minimum and TS geometries obtained by the SC-AFIR2 search are shown in the Supporting Information, Table SI8 and Table SI9, respectively. We know that, depending on the purposes as discussed below, the performance of the SC-AFIR method improves when searching on the potential energy B

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

Letter

Journal of Chemical Theory and Computation

stable TS, additional calculations are required to predict correct reactivity and selectivity. Grambow et al. reported that many search trials were not successful because of SCF convergence failures or geometry optimization failures, or because the same TS was discovered multiple times.4 However, in the SC-AFIR search, SCF convergence hardly failed; SCF convergence failed 4 times in SC-AFIR1, 18 in SC-AFIR2, and twice in SC-AFIR1+. As the ratio of all gradient calls, those failures were 0.0008% for SCAFIR1, 0.0009% for SC-AFIR2, and 0.0006% for SC-AFIR1+. The reason why SCF convergence rarely fails in SC-AFIR is extrapolated that starting from a single input structure, SCAFIR performs a structure search continuously with a relatively small step size (the norm of geometrical displacement in the Cartesian coordinate is limited to less than 0.5 Å). Moreover, because the error-handling is fully automated in GRRM 17 when SCF convergence fails, the user did not need to bother restarting the calculation. Meanwhile, the optimization of the TS structure was performed 1460 times, 5675 times, and 686 times in SC-AFIR1, SC-AFIR2, and SC-AFIR1+, respectively, of which, convergence did not occur 166 times, 828 times, and 80 times. By ratio, the optimization of the TS structure did not converge 11% in SC-AFIR1, 15% in SC-AFIR2, and 12% in SC-AFIR1+. These ratios can be significantly reduced by specifying an option to calculate exact Hessian in the first step of the optimization of the TS structure. On the other hand, the SC-AFIR also found the same TS multiple times. This happens when a reactant molecule has many conformers. Let us discuss a case where a reactant has two conformers A and B. The SCAFIR tries to find a path to a product X from both A and B. In case there are direct paths to X from both A and B, the two trials find two unique TSs leading to X; these two TSs have different conformation and energy. While, in the case where a direct path exists only from A, the trial from B finds a two-step path through two TSs: one is a TS for conformational change from B to A, and the other is a TS from A to X. In the latter case, SC-AFIR finds the same TS multiple times. We emphasize that such trials are necessary, because it is required to identify the lowest energy TS for each chemical transformation as discussed in the previous paragraph. Quality of ab initio and density functional theory (DFT) calculations should also be discussed. Some paths may appear or disappear depending on the ab initio method or DFT functional and basis set adopted. One of the serious examples is the roaming channel.14 At TSs for roaming channels, molecules take a shape in which two radical species are bound weakly to each other. Therefore, TSs for roaming channels are overlooked when a spin-restricted closed-shell single-reference theory is adopted. Hence, the very first paper that discussed roaming channels even declared that there is no TS in the roaming channel. Later studies have proved that the corresponding TS actually exists. In unimolecular decomposition of KHP, roaming channels may also play roles when the excess energy is sufficiently high. However, TSs for most of the roaming channels are missed in the above data set composed of hundreds of TSs optimized at the B3LYP/631+G(d) level. It was suggested that TSs for roaming channels can be explored through automated reaction path searches based on spin-unrestricted open-shell DFT calculations and subsequent reoptimizations based on multireference calculations.15 It is interesting to show roaming channels in the unimolecular decomposition of KHP, but it is beyond the scope of this paper. Since many methods use relative energy of

calculation until the local minimum structure was reached; therefore, the SMILES notation of the product is different from the one reported by Grambow et al. for these channels. Also, R32 is indicated in the column of entry 38 following the previous attribution4 because the R32 geometry in the Supporting Information4 could not be reproduced in this study. The energy of the TS with the lowest energy obtained by each method is also listed in Table 1. Because the energy reference is different in each channel and the KHP conformer used as the reference is not reported in ref 4, we copied values taken from the Supporting Information of ref 4 to Table 1. On the other hand, we consistently used the highest energy conformer as the energy reference in Table 1. There are 22 conformers for KHP at the present computational level, and the gap between the highest and lowest energy conformers is 2.1 kcal/mol. Therefore, when the TS energy of the other methods (i.e., FSM, GSM, SSM, KinBot1, and KinBot2) is higher than ours by more than 2.1 kcal/mol, they missed the lowest TS. On the other hand, when the other methods’ TS energy is lower than ours, ours missed the lowest TS. When the TS energy of the other methods falls between ours and ours plus 2.1 kcal/mol, the same TS was assumed to be found. As shown in Table SI1, our best result, SC-AFIR2, detected 80 reaction channels. It also shows that the number of channels obtained decreases in the order of SC-AFIR2 (80) > SC-AFIR1 (41) > SSM (37) = KinBot2 (37) > GSM (36) > SC-AFIR1+ (33) > FSM (30) > KinBot1 (25), where the number of detected channels is indicated in the parentheses. Based on this method of benchmarking, it can be concluded that our method can find many more than 7 channels if it is used properly, as described above. All three calculations described above were performed using the GRRM17 program; therefore, reproduction is possible by using this program with the combination of γ = 1000 kJ/mol and NoBondRearrange options. All input files used are listed in the Supporting Information, Tables SI2−SI5. Next, we would like to share our opinion on the way to compare automated reaction path search methods. We emphasize that the more versatile the automated reaction path search method is, the more options are available. The AFIR method is a highly versatile method. So far, it has been applied to analyze various chemical reactions, including thermal reaction of small molecule, organic reaction, organometallic catalysis, enzyme catalysis with QM/MM method, crystal phase transition and surface reaction using periodic boundary conditions, photoreaction taking multiple electronically excited states into consideration, etc. However, versatility comes at a cost: there are many options, and if wrong options are selected, a false conclusion can be drawn. It should also be noted that comments on the results may be biased according to the interests of those who make a comparison. Grambow et al. focused on the comprehensiveness of the reaction channels detected.4 Although this was not their focus, finding the TSs with the lowest energy is especially important in many practical applications. In these cases, it is necessary to identify the one having the lowest energy among multiple TSs that generate the same product.13 For example, there are multiple TSs with different conformations in the reaction channels of entries 1 and 2. Actually, we found five TSs for entry 1 and eight TSs for entry 2. They differ in conformations and energy of TSs, although the products generated are identical. If a method fails to identify the most C

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

Letter

Journal of Chemical Theory and Computation

showed the results side by side. We would like to point out that in using FSM, GSM, SSM, and KinBot, Grambow et al. made efforts to select the combination of the best options upon their main developers’ advice; however, in using SCAFIR, they selected a combination of options that seems obviously odd to its main developers. This paper confirmed that the SC-AFIR search results improved significantly depending on the options selected. Through this result, we want to show the readers that the choice of options greatly impacts the cost and comprehensiveness of the automated reaction path search methods. In the meantime, we would like to emphasize that it is not our intent to discuss the superiority of any of these methods in this paper. Moreover, further discussion would be necessary on how to evaluate the performance of the automated reaction path search methods. This method is considered to become important more than ever. We sincerely hope that this paper will help in research and development of the method and in its fair evaluation in the future.

obtained local minima and TSs during automated reaction path searches in deciding PES areas to be explored, errors in those values may also cause a serious problem. To solve this problem, a systematic procedure for estimating and controlling the errors has been proposed.16 Since automated reaction path search in general is computationally expensive, there are cases where the initial exhaustive search is made with a low-level computational method such as semiempirical methods. However, one needs to keep in mind that this could cause overlooking of important TSs. In this regard, the result by KinBot could have been affected by this problem, because the search by KinBot included a prescreening step in which AM1 or HF/STO-3G was employed to reduce the computational cost.4 Calculation cost is also an important criterion. Calculation volume is large in the order of SC-AFIR2 (2033336) > GSM (756227)4 > SSM (649589)4 > SC-AFIR1 (510960) > SCAFIR1+ (346692) > FSM (65420)4 > KinBot2 (23458/ 25333)4 > KinBot1 (8046/10093)4, where the number of gradient calls is indicated in the parentheses. The second number in the parentheses of KinBot2 and KinBot1 indicates the number of gradient calls required in AM1 and HF/STO3G constrained minimization steps, respectively. We note here that the CPU time is roughly proportional to the number of gradient calls. At a glance, some people might feel that methods with low cost are better. However, this is not always true. Among five methods, GSM, SSM, and SC-AFIR provide approximate paths that are obtained by optimizing the entire path. On the other hand, FSM avoids optimizing the entire path accurately and KinBot tries to optimize TSs directly. For example, in the case of SC-AFIR, approximate kinetics can be discussed even with the network of approximate paths. In the SC-AFIR1+ calculation, the approximate path network was constructed through 139196 gradient calls. This path network itself is very useful in rough evaluations of each channel since the network includes all the channels obtained by SC-AFIR1+, which are listed in Table 1. Then, bond rearranging paths with barriers less than 40 kcal/mol (refer to the lowest energy reactant in the network made by SC-AFIR1+) were selected from the network and optimized to actual TSs through additional 1496 gradient calls. With this procedure, the lowest TSs for entry 1 and entry 2 were obtained. In other words, the total number of 140692 gradient calls were performed in this calculation, and the calculation gave the approximate path network that comprehensively includes reaction channels and the actual TSs for the lowest lying paths. For many purposes, such a result is good enough to discuss or predict unknown reaction mechanisms. In other words, users should choose methods and options based on not only cost performance but also various perspectives. Many other options to reduce computational costs for specific purposes should be available in each automated reaction path search method; see for example AFIR.17 We again emphasize that the more versatile the method is, the more useful options are available, and proper usage of options will greatly enhance usefulness of automated reaction path search methods. We have utilized our AFIR method in various practical applications making the best use of such options.10,12,18 As described above, performance comparison of automated reaction path search methods needs various considerations. Grambow et al. performed the search using five different methods with a single combination of options respectively and

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b01182. Complete list of discovered channels and relative energy values of best (energetically most preferable) TSs (PDF) Input files for SC-AFIR1, SC-AFIR2, and SC-AFIR1+ calculations (PDF) Cartesian coordinates of minima and TSs obtained by SC-AFIR1, SC-AFIR2, and SC-AFIR1+ searches (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Satoshi Maeda: 0000-0001-8822-1147 Funding

S.M. is supported by JST, CREST with grant number JPMJCR14L5. Y.H. is supported by JST, PRESTO with grant number JPMJPR16N8. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We thank Ms. Takako Homma for editing a draft of this manuscript. REFERENCES

(1) 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. Dewyer, A. L.; Argüelles, A. J.; Zimmerman, P. M. Methods for exploring reaction space in molecular systems. WIREs Comput. Mol. Sci. 2018, 8, e1354. Simm, G. N.; Vaucher, A. C.; Reiher, M. J. Phys. Chem. A 2019, 123, 385−399. (2) Maeda, S.; Ohno, K. Global Mapping of Equilibrium and Transition Structures on Potential Energy Surfaces by the Scaled Hypersphere Search Method: Applications to Ab Initio Surfaces of Formaldehyde and Propyne Molecules. J. Phys. Chem. A 2005, 109, 5742−5753. Maeda, S.; Morokuma, K. Communication: A Systematic Method for Locating Transition Structures of A + B → X Type Reactions. J. Chem. Phys. 2010, 132, 241102. Rappoport, D.; Galvin, D

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

Letter

Journal of Chemical Theory and Computation C. J.; Zubarev, D. Y.; Aspuru-Guzik, A. Complex Chemical Reaction Networks from Heuristics-Aided Quantum Chemistry. J. Chem. Theory Comput. 2014, 10, 897−907. Kim, Y.; Choi, S.; Kim, W. Y. Efficient Basin-Hopping Sampling of Reaction Intermediates through Molecular Fragmentation and Graph Theory. J. Chem. Theory Comput. 2014, 10, 2419−2426. Habershon, S. Sampling reactive pathways with random walks in chemical space: Applications to molecular dissociation and catalysis. J. Chem. Phys. 2015, 143, 094106. Martínez-Núñez, E. An Automated Method to Find Transition States Using Chemical Dynamics Simulations. J. Comput. Chem. 2015, 36, 222−234. Bergeler, M.; Simm, G. N.; Proppe, J.; Reiher, M. Heuristics-Guided Exploration of Reaction Mechanisms. J. Chem. Theory Comput. 2015, 11, 5712−5722. Zhang, X.-J.; Liu, Z.-P. Reaction sampling and reactivity prediction using the stochastic surface walking method. Phys. Chem. Chem. Phys. 2015, 17, 2757− 2769. Wang, L.-P.; McGibbon, R. T.; Pande, V. S.; Martinez, T. J. Automated Discovery and Refinement of Reactive Molecular Dynamics Pathways. J. Chem. Theory Comput. 2016, 12, 638−649. Yang, M.; Zou, J.; Wang, G.; Li, S. Automatic Reaction Pathway Search via Combined Molecular Dynamics and Coordinate Driving Method. J. Phys. Chem. A 2017, 121, 1351−1361. (3) Aspuru-Guzik, A.; Lindh, R.; Reiher, M. The Matter Simulation (R)evolution. ACS Cent. Sci. 2018, 4, 144−152. (4) Grambow, C. A.; Jamal, A.; Li, Y. P.; Green, W. H.; Zador, J.; Suleimanov, Y. V. Unimolecular Reaction Pathways of a gammaKetohydroperoxide from Combined Application of Automated Reaction Discovery Methods. J. Am. Chem. Soc. 2018, 140, 1035− 1048. (5) Behn, A.; Zimmerman, P. M.; Bell, A. T.; Head-Gordon, M. Efficient exploration of reaction paths via a freezing string method. J. Chem. Phys. 2011, 135, 224108. Suleimanov, Y. V.; Green, W. H. Automated Discovery of Elementary Chemical Reaction Steps Using Freezing String and Berny Optimization Methods. J. Chem. Theory Comput. 2015, 11, 4248−4259. (6) Mallikarjun Sharada, S.; Zimmerman, P. M.; Bell, A. T.; HeadGordon, M. Automated Transition State Searches without Evaluating the Hessian. J. Chem. Theory Comput. 2012, 8, 5166−5174. Zimmerman, P. M. Automated discovery of chemically reasonable elementary reaction steps. J. Comput. Chem. 2013, 34, 1385−1392. (7) Zimmerman, P. M. Single-Ended Transition State Finding with the Growing String Method. J. Comput. Chem. 2015, 36, 601−611. (8) Maeda, S.; Taketsugu, T.; Morokuma, K. Exploring transition state structures for intramolecular pathways by the artificial force induced reaction method. J. Comput. Chem. 2014, 35, 166−173. (9) Zádor, J.; Najm, H. N. In KinBot 1.0: a code for automatic PES exploration; Report No. SAND2012-8095; Sandia National Laboratories: 2013. (10) Maeda, S.; Harabuchi, Y.; Takagi, M.; Saita, K.; Suzuki, K.; Ichino, T.; Sumiya, Y.; Sugiyama, K.; Ono, Y. Implementation and performance of the artificial force induced reaction method in the GRRM17 program. J. Comput. Chem. 2018, 39, 233−251. Maeda, S.; Harabuchi, Y.; Sumiya, Y.; Takagi, M.; Suzuki, K.; Hatanaka, M.; Osada, Y.; Taketsugu, T.; Morokuma, K.; Ohno, K. GRRM17. See: http://iqce.jp/GRRM/index_e.shtml (accessed Jan 24, 2019). (11) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A. V.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J. J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Keith, T. A.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.;

Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 16, Revision B.01; Gaussian, Inc.: Wallingford, CT, 2016. (12) Maeda, S.; Morokuma, K. Finding Reaction Pathways of Type A + B → X: Toward Systematic Prediction of Reaction Mechanisms. J. Chem. Theory Comput. 2011, 7, 2335−2345. Maeda, S.; Harabuchi, Y.; Takagi, M.; Taketsugu, T.; Morokuma, K. Artificial Force Induced Reaction (AFIR) Method for Exploring Quantum Chemical Potential Energy Surfaces. Chem. Rec. 2016, 16, 2232−2248. (13) Maeda, S.; Ohno, K. Lowest transition state for the chiralitydetermining step in Ru((R)-BINAP)-catalyzed asymmetric hydrogenation of methyl-3-oxobutanoate. J. Am. Chem. Soc. 2008, 130, 17228−17229. Donoghue, P. J.; Helquist, P.; Norrby, P. O.; Wiest, O. Prediction of enantioselectivity in rhodium catalyzed hydrogenations. J. Am. Chem. Soc. 2009, 131, 410−411. Seal, P.; Papajak, E.; Truhlar, D. G. Kinetics of the Hydrogen Abstraction from Carbon-3 of 1Butanol by Hydroperoxyl Radical: Multi-Structural Variational Transition-State Calculations of a Reaction with 262 Conformations of the Transition State. J. Phys. Chem. Lett. 2012, 3, 264−271. Zheng, J.; Meana-Pañeda, R.; Truhlar, D. G. MSTor version 2013: A new version of the computer code for the multi-structural torsional anharmonicity, now with a coupled torsional potential. Comput. Phys. Commun. 2013, 184, 2032−2033. Hatanaka, M.; Maeda, S.; Morokuma, K. Sampling of Transition States for Predicting Diastereoselectivity Using Automated Search Method-Aqueous Lanthanide-Catalyzed Mukaiyama Aldol Reaction. J. Chem. Theory Comput. 2013, 9, 2882−2886. Lime, E.; Lundholm, M. D.; Forbes, A.; Wiest, O.; Helquist, P.; Norrby, P. O. Stereoselectivity in Asymmetric Catalysis: The Case of Ruthenium-Catalyzed Ketone Hydrogenation. J. Chem. Theory Comput. 2014, 10, 2427−2435. (14) Townsend, D.; Lahankar, S. A.; Lee, S. K.; Chambreau, S. D.; Suits, A. G.; Zhang, X.; Rheinecker, J.; Harding, L. B.; Bowman, J. M. The Roaming Atom: Straying from the Reaction Path in Formaldehyde Decomposition. Science 2004, 306, 1158−1161. Herath, N.; Suits, A. G. Roaming Radical Reactions. J. Phys. Chem. Lett. 2011, 2, 642−647. Bowman, J. M.; Houston, P. L. Theories and Simulations of Roaming. Chem. Soc. Rev. 2017, 46, 7615−7624. (15) Maeda, S.; Taketsugu, T.; Ohno, K.; Morokuma, K. (Perspective) From Roaming Atoms to Hopping Surfaces: Mapping Out Global Reaction Routes in Photochemistry. J. Am. Chem. Soc. 2015, 137, 3433−3445. (16) Proppe, J.; Husch, T.; Simm, G. N.; Reiher, M. Uncertainty Quantification for Quantum Chemical Models of Complex Reaction Networks. Faraday Discuss. 2016, 195, 497−520. Simm, G. N.; Reiher, M. Systematic Error Estimation for Chemical Reaction Energies. J. Chem. Theory Comput. 2016, 12, 2762−2773. Simm, G. N.; Reiher, M. Error-Controlled Exploration of Chemical Reaction Networks with Gaussian Processes. J. Chem. Theory Comput. 2018, 14, 5238−5248. (17) AFIR-web. 2018. See: https://afir.sci.hokudai.ac.jp/ (accessed Oct 31, 2018). (18) Sameera, W. M.; Maeda, S.; Morokuma, K. Computational Catalysis Using the Artificial Force Induced Reaction Method. Acc. Chem. Res. 2016, 49, 763−773.

E

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