Article pubs.acs.org/jcim
Low-Mode Conformational Search Method with Semiempirical Quantum Mechanical Calculations: Application to Enantioselective Organocatalysis Takashi Kamachi* and Kazunari Yoshizawa Institute for Materials Chemistry and Engineering, Kyushu University, Fukuoka 819-0395, Japan S Supporting Information *
ABSTRACT: A conformational search program for finding low-energy conformations of large noncovalent complexes has been developed. A quantitatively reliable semiempirical quantum mechanical PM6-DH+ method, which is able to accurately describe noncovalent interactions at a low computational cost, was employed in contrast to conventional conformational search programs in which molecular mechanical methods are usually adopted. Our approach is based on the low-mode method whereby an initial structure is perturbed along one of its low-mode eigenvectors to generate new conformations. This method was applied to determine the most stable conformation of transition state for enantioselective alkylation by the Maruoka and cinchona alkaloid catalysts and Hantzsch ester hydrogenation of imines by chiral phosphoric acid. Besides successfully reproducing the previously reported most stable DFT conformations, the conformational search with the semiempirical quantum mechanical calculations newly discovered a more stable conformation at a low computational cost.
1. INTRODUCTION Asymmetric organocatalysis is one of the most useful methods for a wide variety of enantioselective transformations because of operational simplicity, mild reaction conditions, and the environmentally benign nature of the reaction.1−3 Unlike traditional metal-based catalysts with well-defined metal−ligand coordination, noncovalent interactions, such as hydrogenbonding interactions, van der Waals interactions, and electrostatic interactions, play a central role in controlling the stereochemical outcome of the organocatalysts. Although these interactions are relatively weak compared with covalent bonds, they can collectively contribute to the well-ordered binding between the organocatalysts and substrates. The high enantioselective organocatalytic reactions have attracted considerable attention of theoretical chemists, and density functional theory (DFT) methods have been becoming a powerful tool to predict the structure and reactivity of organocatalysts.4−8 Although a great number of organocatalytic systems have been successfully investigated with DFT methods, it remains a daunting task to reveal the most stable structure of stereodeterming transition states due to their conformational flexibility. In a previous study, 9 we performed a systematic conformation analysis of binaphthyl-modifed chiral phasetransfer catalysts 1 (so called Maruoka catalyst10−13) by combining a conventional molecular mechanics (MM) based conformational search with semiempirical quantum mechanical (SQM) and DFT calculations (Figure 1). A total of 1436 initial conformations were obtained with the MacroModel program14 when all conformations within 20 kcal/mol above the global minimum were considered with the OPLS 200515 force field. The obtained structures were reoptimized by using the PM6© XXXX American Chemical Society
Figure 1. Asymmetric benzylation of glycine derivatives 2 with the Maruoka catalyst.
DH+ SQM method16 implemented in the MOPAC program.17 Then, we optimized 154 refined structures whose MOPAC energy was within 10 kcal/mol above the global minimum using the B97-D method.18 The C···C distance between the αcarbon of the glycine enolate ion 2 and the benzylic carbon of benzyl bromide was fixed to 2.7 Å in these calculations to generate a reasonable initial guess for subsequent transition state searches. Finally, the obtained top 20 conformations of the partially optimized structures were fully optimized to locate the transition states. The high enantioselectivity was found to arise from an energy difference of 5.3 kcal/mol between the most stable conformation of the transition state for R- and Sproduction. The computed transition structures well explain experimentally observed substituent effects on the enantioseReceived: November 5, 2015
A
DOI: 10.1021/acs.jcim.5b00671 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Information and Modeling
energy minimization. This process is repeated to find additional structures. The low-mode method has been successfully applied to protein loop optimization and flexible active site docking.22 This method is advantageous over torsional Monte Carlo methods,23 particularly for the organocatalytic system containing more than one molecule because torsional angle cannot be defined between the fragments. The procedure of the low-mode conformational search with SQM calculations is summarized in Figure 2. The conformation
lectivity. The computational results were reinforced by a recent artificial force induced reaction (AFIR) study.19 The conformational search provides a promising approach to reveal the origin of the enantioselctivity of organocatalytic reactions. However, we noticed that the MM part of the procedure sometimes causes a serious problem in the conformation search, leading to a failure to generate lowenergy conformations despite the extensive and timeconsuming computations. As summarized in Table 1, lowTable 1. Lowest Five DFT Conformations of Transition State for R-Production, and SQM and MM Ones Linked to the DFT Conformations in the Conformational Search Procedurea DFT (B97-D)b
SQM (PM6-DH+)b
MM (OPLS2005)
rank
energy
rank
energy
rank
energy
1 2 3 4 5
0.0 1.1 3.2 3.7 4.1
7 4 8 3 13
1.8 1.4 1.9 0.5 3.1
205 53 50 69 159
11.3 7.6 5.5 7.2 9.2
In the conformations, the key C···C distance was fixed to 2.7 Å. Energies (kcal/mol) were measured from the most stable conformation. bThese values are the gas-phase SCF energy without any thermal corrections.
a
energy conformations obtained in the final DFT part can be traced back20 to SQM and MM conformations; for example, the most stable conformation for the enantioselective benzylation by the Maruoka catalyst comes from conformation #7 in the SQM part and #205 in the MM part. All of the lowest five conformations were found in the low-energy top 13 SQM conformations, which indicates that the PM6-DH+ method is reliable and accurate for correctly describing weak interactions that are crucial in the organocatalytic reactions. In contrast, high-energy conformations up to #205 with an energy of 11.3 kcal/mol must be considered to obtain the lowest five conformations. This is because the MM force field is not parameterized to describe electronic structures of transition states and weak interactions such as face-to-face and T-shaped π−π interactions between the organocatalysts and substrates. We are aware that the 20 kcal/mol energy window of the MM part is not sufficient to find the most stable structure for some organocatalytic reactions due to the low accuracy of the MM calculations. In addition, such a high-energy conformation is rarely chosen as starting structures, and therefore, important conformations derived from the high MM energy conformations are overlooked with high probability. Inspired by the computational results, we developed a new program, ConFinder, for conformational search with reliable SQM calculations and demonstrated its effectiveness in searching the low-energy conformations of stereodeterming transition states for organcatalysis.
Figure 2. Flow diagram of the low-mode conformation search with SQM calculations.
search starts with energy minimization of an initial arbitrary structure with a SQM method using a convergence criterion of 0.01 kcal/(mol Å). The optimized structure is subjected to normal-mode analysis with the SQM method. All the SQM calculations were performed with the MOPAC program. The optimized structure is perturbed by moving along a lowfrequency vibrational mode in both directions. Lowest 10 lowfrequency modes were explored in this study. The eigenvector is scaled to move the fastest moving atom of the selected mode to a randomly chosen distance between a user-specified minimum and maximum value. We recommend a distance range of 6.0−9.0 Å for the first run to locate previously nonsampled regions and that of 3.0−6.0 Å for a restart run to find neighboring conformations of the previously identified structures. As shown in Figure S1, the long leap usually generates an unreasonable structure, leading to a SCF convergence failure in subsequent SQM calculations. To solve this problem, the perturbed structure is subjected to MM minimization with the TINKER program24 using a root mean
2. COMPUTATIONAL METHODS We employed the low-mode algorithm of Kolossváry and Guida21 to generate new structures. The low-mode method is based on an eigenvector following of low-frequency vibrational modes (typically
