CHEMPHYSCHEM ARTICLES DOI: 10.1002/cphc.201402105
A Density-Based Adaptive Quantum Mechanical/Molecular Mechanical Method Mark P. Waller,* Sadhana Kumbhar, and Jack Yang[a] We present a density-based adaptive quantum mechanical/molecular mechanical (DBA-QM/MM) method, whereby molecules can switch layers from the QM to the MM region and vice versa. The adaptive partitioning of the molecular system ensures that the layer assignment can change during the optimization procedure, that is, on the fly. The switch from a QM molecule to a MM molecule is determined if there is an absence of noncovalent interactions to any atom of the QM core
region. The presence/absence of noncovalent interactions is determined by analysis of the reduced density gradient. Therefore, the location of the QM/MM boundary is based on physical arguments, and this neatly removes some empiricism inherent in previous adaptive QM/MM partitioning schemes. The DBA-QM/MM method is validated by using a water-in-water setup and an explicitly solvated l-alanyl–l-alanine dipeptide.
1. Introduction Modeling complex chemical systems in the gas or liquid phase, in the solid state, or on surfaces typically requires advanced mutliscale techniques.[1–3] Hybrid quantum mechanical/ molecular mechanical (QM/MM) methods partition a system into an active region that is treated with a high-level electronic structure method, and the surrounding environment can be inexpensively treated by a lower-level force field.[4, 5] One can therefore achieve good accuracy, while retaining efficiency, for studying complex chemical systems. For this reason, QM/MM methods have emerged as the method of choice for studying reaction mechanisms in large biomolecular systems. Furthermore, if certain parts of a system are either missing from the experimental data or if significant parts of the system have been modified, that is, by using chemical biology tools, then the whole system can be refined/remodeled with QM/MM methods.[6, 7] Multiscale methods for complex chemical systems are traditionally based on the rigid partitioning of the active region and the environment. To overcome this shortcoming and to investigate nonlocalized active regions more accurately, on-the-fly switching (i.e. adaptive QM/MM) is of great interest to the modeling community.[8–22] In adaptive schemes, the QM and MM assignments of molecules are not rigid, and molecules can therefore freely change their layers by repartitioning the system. Accurate modeling of complex chemical systems is largely reliant upon a solid experimental foundation, that is, X-ray [a] Dr. M. P. Waller, S. Kumbhar, Dr. J. Yang Theoretische Organische Chemie, Organisch-Chemisches Institut and Center for Multiscale Theory and Computation Westfalische Wilhelms-Universitt Mnster Corrensstrasse 40, 48149 Mnster (Germany) Fax: (+ 49) 0251 83 33202 E-mail: [email protected] Supporting Information for this article is available on the WWW under http://dx.doi.org/10.1002/cphc.201402105.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
crystallography or elegant NMR spectroscopy experiments to provide a starting structure. Traditional gradient-based optimizers converge to a local minimum, which is generally dependent upon the starting point.[23] In contrast, nonlocal optimization methods can reveal a variety of low-energy minima[24] and are far more explorative than traditional gradient-based geometry optimizers. Therefore, nonlocal optimizers are a good choice for de novo structure prediction of complex chemical systems. Metaheuristic algorithms nicely blend deterministic and stochastic aspects to generate diverse molecular structures. This is important if the chemical system is too complex to iterate through all possible conformations with a brute-force approach due to the high dimensionality of the search space. Nonlocal metaheuristic optimizers are generally not based on the gradient, and an entirely new solution is typically built at each metacycle of the algorithm, and this means that the energies do not need to be continuous. The traditional rigid partitioning of the QM and MM regions cannot be easily implemented in conjunction with metaheuristic methods. This is because neighboring fragments in any given cycle are not necessarily neighboring fragments from either the previous or subsequent optimization cycles. This is due to the highly explorative nature of the optimizer, which can place a fragment anywhere within a given search space, rather than the sequential application of short step lengths that are employed in conventional gradient-based optimizers (see Figure 1). Hence, there is a need for a method that can analyze the entire system at every metacycle, and this can be accomplished by employing an on-the-fly partitioning scheme (i.e. nonrigid partitioning of the entire system into either a QM or MM region). In general, partitioning the system into QM and MM regions is reliant upon a criterion, and to date, adaptive schemes can be categorized as either distance-based or ChemPhysChem 0000, 00, 1 – 9
&1&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
www.chemphyschem.org
zero, in the regions of covalent bonding and noncovalent interactions. These two regions can further be distinguished by using the density values; low density 1!0 combined with low s(1) represents the noncovalent interactions. As a consequence, this allows a spatial region W to be mapped from s(1) for constructing a visualizable RDG isosurface, V(r), in real space, which indicates the presence of noncovalent interactions between molecules (or molecular fragments). The spatial Figure 1. Schematic of a nonlocal optimizer (green) and a local optimizer (red) for the water dimer that share the extent of this RDG isosurface is same starting point. considered as a numerical indicator of the strength of the NCI between two interacting systems. This merit can be obtained number-based adaptive schemes (see Figure 2). The former from the integrated volume of RDG, defined as [Eq. (2)]: scheme defines the distance cutoff between the position of each molecule and the center of the QM core region, and the Z smax Z 1max latter one is defined by predefining the total number of molef ¼ dW ð2Þ 0 0 cules that are allowed in the QM region.[25] Setting an appropriate QM/MM boundary is crucial for successfully modeling moin which smax and 1max are set to threshold values, and f is the lecular systems with QM/MM methods. If the boundary to be integrated volume of the spatial area of the RDG isosurface set is too close to the active site, the chemical realism of the V(r). Hence, one can determine the presence or absence of model will be reduced. In contrast, if one positions the bounnoncovalent interactions by analyzing the intermolecular dendary too far away from the active site, the computational exsity using the NCI index. pense of the QM calculation may become problematic, possibly including difficulties in self-consistent field and geometry In summary, two QM/MM partitioning criteria were previousconvergence. Each molecular system is different in this respect, ly proposed (distance and number based), and herein we exand one must perform thorough tests to validate the appropriplored an alternative possibility of employing a nonempirical ateness of the particular QM/MM partitioning criteria used. partitioning criterion based on the NCI concept (see Figure 2). In 2010, Johnson et al.[26] introduced the concept of utilizing reduced charge density gradients (RDGs) to visually analyze Computational Section noncovalent interactions (NCI), and the NCIPLOT program was The DBA-QM/MM method was implemented into our own batchsubsequently released by Contreras-Garca et al.[27] The NCI based framework: JACOB.[32] concept was further developed to incorporate more quantifiable analysis of noncovalent interactions in fluctuating environPartitioning ments.[28] Experimental charge-density studies with the use of [29] [30, 31] models have independently the MOON and Multipole In principle, a density-based partitioning scheme requires a model confirmed the presence of RDGs in the crystalline environof the density for both the QM and the MM regions. The simplest ment. approach, and therefore the one adopted herein, is to assume atom-centered spherical densities for both regions. This spherical In short, according to the NCI concept, if two noncovalently density is also referred to as a promolecular density,[33] and it is interacting molecules are brought together, a finite RDG based a superposition of typically gas-phase neutral atomic densities. Of on Equation (1):
sð1Þ ¼
jr1j 4 2ð3p Þ 13 1
2
1 3
ð1Þ
where 1 is the total electron density and r1 is the electron density gradient, appears in the intermolecular region in which the total electronic density 1!0. To identify the intermolecular interactions, one can examine the plots of s(1) versus 1. The s(1) plots were shown to have very small values, approaching 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
course, the calculated/tabulated promolecular density for the MM region could be combined with an aspherical quantum mechanically derived charge density that naturally arises from the QM region. However, to expedite the numerical calculation, promolecular densities were used for both the QM and MM regions in the current implementation of our density-based partitioning scheme. The promolecular densities were shown to be capable of reproducing the topological features of a fully quantum mechanical calculation (i.e. aspherical density).[34, 35]
The NCIPLOT v1.0 program by Contreras-Garca et al.[27] was employed to analyze the intermolecular interactions by using the moChemPhysChem 0000, 00, 1 – 9
&2&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
www.chemphyschem.org Therefore, one can specifically investigate only the intermolecular interactions within supramolecular assemblies, and this approach was adopted for all systems in our present study. In our DBA-QM/MM method, any constituent of the system that must be treated quantum mechanically, for example, species involved in bond breaking/forming processes were placed into the static core region (QMcore). The region of space between the QMcore and each other molecule in the environment was analyzed for the presence of noncovalent interactions based on the NCI index. If, and only if, a noncovalent interaction was detected between the QMcore region and any given neighboring molecule, then that neighboring molecule was then included into the QM region. If no noncovalent interactions were observed to a given neighboring molecule, then this neighboring molecule was placed into the MM region, whereby the interaction between the QM and MM was treated by the classical electrostatic and van der Waals coupling terms. Repartitioning of the entire system throughout the optimization procedure was performed on the fly (see Figure 3).
Figure 3. The partitioning procedure constructs pairwise mutually independent (trivially parallelizable) calculations for NCIPLOT.[27]
Hybrid QM/MM The two-layer ONIOM method[36] within Gaussian 09[37] was used for the QM/MM single-point energy evaluations at the B3LYP[38, 39]/ 6-31G* or 6-311 + + G(2d,2p):[40–42] Amber[43] level of theory. These levels of theory were chosen, as they were utilized in the reference calculations, and they also represent common choices in the computational chemistry community. Furthermore, the same functional was used throughout this work because we were solely interested in investigating the effects of the partitioning scheme, which are not functional dependent. No healing region (alternatively called transition zone) was employed in this current study, and this had the distinct advantage of being computationally faster than methods that employ such a transition. The application of a nonlocal optimization technique in this current study meant that the healing region was unnecessary, see below.
Nonlocal Optimization Figure 2. Schematic partitioning of a small water cluster: QM (vdW radii)/ transition (ball and stick)/MM (thin wire). The green isosurfaces indicate the presence of a noncovalent interaction according to the NCI index.
lecular densities that were constructed for each monomer (i.e. 1molecular). The INTERMOLECULAR keyword was used so that only contributions from the intermolecular interactions would remain.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Our own Hybrid Metaheuristic Optimizer (HMO) method[44] was used to perform nonlocal optimization[24] of the molecular systems. The HMO method hybridizes a particle swarm optimizer (PSO)[45] with an ant colony optimizer (ACO),[46, 47] whereby the PSO optimizes the parameter space of the underlying ACO (see Table 1), whereas the ACO itself optimizes the conformational space of the ChemPhysChem 0000, 00, 1 – 9
&3&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
www.chemphyschem.org
Table 1. Parameters (optimized by using HMO) used for the ACO that were fitted to a water trimer by employing the PM6-DH + [48] Hamiltonian in MOPAC2012.[49] Parameter[a]
Optimized value
number of ants n pheromone resistance probability to choose global best global best percentage explorative/exploitive search controller trail width
20 0.25 0.30 0.30 1.00 0.01
[a] The presented parameters were used throughout this entire study.
molecular system. A supramolecular search space was used for each of the solvent water molecules (six degrees of freedom for each solvent molecule, three translational and three rotational) and the dihedral search space was used for the l-alanyl–l-alanine dipeptide (two degrees of freedom that is, f and y). The overall procedure for the DBA-QM/MM is depicted in the sequence diagram shown in Figure 4. The three-step process was run until the geometry converged or the maximum number of predetermined metacycles was reached.
Figure 5. A scan of the intercentroid distance of the water dimer showing the integrated volume from the NCI analysis of the promolecular density. Water dimers are overlaid schematically.
The orientational dependence of the intermolecular noncovalent interaction for the water dimer is plotted in Figure 6. Importantly, this particular dependence cannot be recovered 2. Results and Discussion by the distance-based adaptive scheme, and this nicely highWater in Water lights an additional advantage of our density-based approach for investigating nonspherically symmetric systems, that is, A correct description of water and solute–water clusters is critibeyond rare gas matrices. cal for the accurate prediction of structural and thermodynamFurthermore, if one investigates the atom dependency of ic properties. Clusters are formed due to attractive noncovalent the RDG isosurface, one sees that the trends in volumes of the interactions, and they are intermediate between the gas phase RDG reflects the same trends that are seen in the sum of the and the condensed phase; consequently, the interest in investivan der Waals (vdW) radii: O···O (RDG = 15 937 a.u., vdW = gating clusters has increased tremendously over the years. Al3.04 ), O···H (RDG = 14 702 a.u., vdW = 2.72 ), and H···H = though water clusters have been observed experimentally, in(RDG = 13 605 a.u., vdW = 2.40 ); all values were computed at vestigating supramolecular structures in water is difficult due a distance of 3.00 . to continuous H-bond breaking and forming.[50] To systematically examine the performance of the DBA-QM/ The aim of this section is to show a proof of principle for MM partitioning scheme, we optimized the water dimer, water the prototypical solvent, that is, water. An NCI analysis was trimer, and water tetramer. The QM/MM optimizations were performed for the water dimer, whereby the NCI isosurface disperformed over 20 independent runs to assess the reproduciappears if the intermolecular (center of masses) distance bebility. The search space was constructed in such a way that the comes too large (> 4 , see Figure 5). maximum distance between the centers of masses of molecules could be 14 and the search space granularity was set to 0.001 . All possible configurations from the density-based partitioning for the investigated water in water systems are shown in Figure 7. One can see that there are n possible configurations, for which n is the number of non-QMcore molecules. The distances between the centers of masses of the water molecules tended to decrease Figure 4. A sequence diagram showing the key features of the DBA-QM/MM method by using wrappers for NCIduring the nonlocal optimization PLOT[27] and G09.[37] 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
ChemPhysChem 0000, 00, 1 – 9
&4&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
www.chemphyschem.org l-Alanyl–l-Alanine (LALA) Peptide
Figure 6. A scan of the three rotational degree of freedom (a, b, g) for the water dimer by using an NCI analysis of the promolecular density at a constant intercentroid distance of 3 .
The effect of hydration on small peptides and amino acids is not completely understood because these effects are hard to verify with experiments.[52] Therefore, modeling can be used to increase our collective understanding. The zwitterionic form of the l-alanyl–l-alanine (LALAzw) peptide was studied as a model system, because the conformation of the molecule can be characterized by the f and y dihedral angles (see Figure 9).[53] Furthermore, water was shown to play an important role in the conformational changes of LALAzw due to the charged ends of the peptide. Finally, the predominant structure of this particular dipeptide in the pres-
Figure 7. Configurations possible for the dimer, trimer, and tetramer shown top to bottom, respectively.
procedure. The relative energies (DE, kcal mol1) of the pureQM- and DBA-QM/MM HMO optimized structures for each system with respect to the reference values[51] are shown in Figure 8. As the search space increased, so did the diversity of nonlocal minima for both the pure-QM reference and the DBAQM/MM method. Upon comparing the optimized structures, the water molecules were optimized into a nearly isoenergetic conformation from the previous reported structure. The results are very encouraging that our DBA-QM/MM method is capable of exploring multiple minima of comparable energy to the full QM reference optimizations on the explicitly solvated potential energy surface. 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Figure 8. Relative energies compared to those obtained by Bryantsev et al.[51] (DE, kcal mol1) computed for the pure-QM and DBA-QM/MM HMO optimized water dimer, trimer, and tetramer structures.
Figure 9. The zwitterionic form of the l-alanyl–l-alanine dipeptide.
ChemPhysChem 0000, 00, 1 – 9
&5&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
Figure 10. Schematic of LALAzw with 14 water molecules, showing a possible configuration of the entire system.
ence of water is known. It has been shown that LALAzw (in explicit water) predominantly adopts a conformation with the two dihedral angles f 1258, y 1438 during classical molecular dynamics simulations.[53] Furthermore, the dihedral angles f= 1538 and y = 1478 are reported for the LALAzw structure in implicit solvation models.[54, 55] Knapp-Mohammady et al. showed that at least 4 water molecules are required to solvate the end groups of the zwitterionic form of LALA, and dihedral angles for LALAzw with 4, 6, 7, 10, 11, 12, 13, and 14 water molecules have been reported.[56] Therefore, these same systems (i.e. LALAzw with 4–14 explicit water molecules) were adopted to validate the DBA-QM/MM method (see Figure 10). Hybrid QM/MM calculations of the solvated LALAzw were performed at the B3LYP/6-31G* level of theory combined with the Amber force field. The calculations were performed by using the same basis set and functional as those used for the reference structures.[55] We compared our density-based partitioning approach with both the distance- and number-based partitioning approaches for each of the solvated LALAzw systems. It was recently shown that 8 1 water molecules directly interact with LALAzw.[57] Considering that the experimental hydration number is known, the number of solvent molecules around the solute was set to nine in the number-based-QM/ MM method. In the distance-based-QM/MM method, the boundary cutoff was set to 5 . Importantly, no transition zones were used for any of the adaptive partitioning schemes in this present study to directly compare each partitioning method. The median and standard deviations for f and y obtained for each method (i.e. all using the same nonlocal optimizer) were used to plot the normal distributions (see Figure 11). All of the dihedral angles obtained from the different adaptive schemes with different numbers of water molecules (i.e. varying from 4 to 14) are given in the Supporting Information. The DBA-QM/MM method performed in a manner similar to that of the other methods for the f dihedral angle. In the case of the y dihedral angle, the DBA-QM/MM method underestimated the angle to the same extent as the number-based method. The distance-based method performed slightly better 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.chemphyschem.org
Figure 11. Polar graph of the normal distributions of the dihedral angles (f and Y) for LALAzw by using different partitioning schemes over LALAzw plus 4–14 solvent molecules. Reference values are B3LYP/6-311 + + G**/ CPCM (star)[54] and BPW91/6-311 + + G**/PCM (square).[55]
than the DBA or number-based scheme for the y dihedral angle. In general, all three partitioning methods showed similar standard deviations. The pure QM method itself showed deviation across the dihedral angles over the different numbers of explicit water molecules included in the system. Overall, the similar results indicate that all of the methods are capable of returning more or less similar structures for an explicitly solvated LALAzw system. To ensure that our validation of the DBA-QM/MM method was robust with respect to traditional gradient-based optimizers, we also compared the nonlocal optimization results by using different adaptive methods to the reference dihedral angles that were taken from the literature[54, 55] (see Figure 11). Similar values for these dihedral angles were observed across the implicitly solvated gradient-based optimized structures; therefore, the validation of the DBA-QM/MM method is not biased by the choice of reference. In summary, the DBA-QM/MM method performs quite well with respect to the two other well-established methods of adaptive QM/MM partitioning for the explicitly solvated LALAzw system. This is encouraging, as our aim of this current study was not to present a method that necessarily outperformed the two pre-existing methods, but that, more importantly, was devoid of empirical parameters. This means that DBA-QM/MM can in principle be a black-box method that avoids laborious parameterization. The extension of the DBAQM/MM method to other solvated systems (i.e. different solutes and solvents) or to other environments such as surface phenomena is the focus of future work and will be reported in due course.
3. Conclusions A density-based adaptive hybrid QM/MM method able to optimize molecular systems was presented. In general, we showed that analysis of the NCI index was a valid approach to select ChemPhysChem 0000, 00, 1 – 9
&6&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES the boundary cut off for QM/MM calculations by using a number of test systems. The DBA partitioning scheme was shown to adequately describe the highly changing partitioning of the system in combination with the ONIOM subtractive QM/ MM Hamiltonian. Although a subtractive QM/MM Hamiltonian was employed within this work, the DBA-QM/MM partitioning scheme was equally applicable to an additive QM/MM scheme. The DBA QM/MM method was combined with a nonlocal optimization procedure, because continuous energies and gradients were not needed; hence, no healing regions were needed to remove discontinuities. In principle, the DBA-QM/MM method could alternatively be applied with molecular dynamics calculations, if smoothing functions were employed; however, the computational requirements would be significant. The major advantage of our density-based adaptive partitioning is that the method is general, and therefore, unlike the distance- and number-based methods, it does not need to be fitted upon investigating new chemical systems. This is particularly important if the hydration number of a system is not known a priori. Therefore, our nonparameterized density-based adaptive QM/MM method is applicable to a wide range of complex chemical systems.
Acknowledgements Generous financial support from the Deutsche Forschungsgemeinschaft is gratefully acknowledged. Keywords: adaptive · molecular modeling · noncovalent interactions · quantum chemistry · topology [1] N. Lpez, G. Pacchioni, F. Maseras, F. Illas, Chem. Phys. Lett. 1998, 294, 611. [2] E. Eichinger, P. Tavan, J. Hutter, M. Parrinello, J. Chem. Phys. 1999, 110, 10452 – 10467. [3] J. Gao, J. Chem. Phys. 1998, 109, 2346. [4] A. Warshel, M. Levitt, J. Mol. Biol. 1976, 103, 227 – 249. [5] H. M. Senn, W. Thiel, Angew. Chem. 2009, 121, 1220 – 1254; Angew. Chem. Int. Ed. 2009, 48, 1198 – 1229. [6] S. Kumbhar, S. Johannsen, R. K. O. Sigel, M. P. Waller, J. Mller, J. Inorg. Biochem. 2013, 127, 203 – 210. [7] T. Dresselhaus, N. D. Weikart, H. D. Mootz, M. P. Waller, RSC Adv. 2013, 3, 16122 – 16129. [8] T. Kerdcharoen, K. R. Liedl, B. M. Rode, Chem. Phys. 1996, 211, 313 – 323. [9] P. Chieux, H. Bertagnolli, J. Phys. Chem. 1984, 88, 3726 – 3730. [10] T. Kerdcharoen, K. Morokuma, Chem. Phys. Lett. 2002, 355, 257 – 262. [11] T. Kerdcharoen, K. Morokuma, J. Chem. Phys. 2003, 118, 8856 – 8862. [12] G. Csnyi, T. Albaret, M. C. Payne, A. De Vita, Phys. Rev. Lett. 2004, 93, 175503. [13] A. Heyden, H. Lin, D. G. Truhlar, J. Phys. Chem. B 2007, 111, 2231 – 2241. [14] B. Ensing, S. O. Nielsen, P. B. Moore, M. L. Klein, M. Parrinello, J. Chem. Theory Comput. 2007, 3, 1100 – 1105. [15] R. E. Bulo, B. Ensing, J. Sikkema, L. Visscher, J. Chem. Theory Comput. 2009, 5, 2212 – 2221. [16] H. Wang, C. Schtte, L. Delle Site, J. Chem. Theory Comput. 2012, 8, 2878 – 2887. [17] K. Park, A. W. Gçtz, R. C. Walker, F. Paesani, J. Chem. Theory Comput. 2012, 8, 2868 – 2877. [18] C. Rowley, B. Roux, J. Chem. Theory Comput. 2012, 8, 3526 – 3535. [19] B. Lev, B. Roux, S. Y. Noskov, J. Chem. Theory Comput. 2013, 9, 4165 – 4175.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.chemphyschem.org [20] N. Bernstein, C. Vrnai, I. Solt, S. A. Winfield, M. C. Payne, I. Simon, M. Fuxreiter, G. Csnyi, Phys. Chem. Chem. Phys. 2012, 14, 646 – 656. [21] S. Pezeshki, H. Lin, J. Chem. Theory Comput. 2011, 7, 3625 – 3634. [22] R. E. Bulo, C. Michel, P. Fleurat-Lessard, P. Sautet, J. Chem. Theory Comput. 2013, 9, 5567 – 5577. [23] H. B. Schlegel, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 790 – 809. [24] B. Hartke, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 879 – 887. [25] N. Takenaka, Y. Kitamura, Y. Koyano, M. Nagaoka, Chem. Phys. Lett. 2012, 524, 56 – 61. [26] E. R. Johnson, S. Keinan, P. Mori-Snchez, J. Contreras-Garca, A. J. Cohen, W. Yang, J. Am. Chem. Soc. 2010, 132, 6498 – 6506. [27] J. Contreras-Garca, W. Yang, E. R. Johnson, J. Phys. Chem. A 2011, 115, 12983 – 12990. [28] P. Wu, R. Chaudret, X. Hu, W. Yang, J. Chem. Theory Comput. 2013, 9, 2226 – 2234. [29] J. Yang, M. P. Waller, J. Comput. Chem. 2013, 34, 466 – 470. [30] G. Saleh, C. Gatti, L. Lo Presti, J. Contreras-Garca, Chem. Eur. J. 2012, 18, 15523 – 15536. [31] G. Saleh, C. Gatti, L. Lo Presti, Comput. Theor. Chem. 2012, 998, 148 – 163. [32] M. P. Waller, T. Dresselhaus, J. Yang, J. Comput. Chem. 2013, 34, 1420 – 1428. [33] P. Coppens, IUCr Texts on Crystallography, Vol. 4, Oxford University Press, Oxford, 1997. [34] M. A. Spackman, Chem. Phys. Lett. 1999, 301, 425 – 429. [35] C. Gatti, E. May, R. Destro, F. Cargoni, J. Phys. Chem. A 2002, 106, 2707 – 2720. [36] M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber, K. Morokuma, J. Phys. Chem. 1996, 100, 19357 – 19363. [37] Gaussian 09, Revision D.01, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, . Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2009.. [38] A. D. Becke, J. Chem. Phys. 1993, 98, 5648 – 5652. [39] C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785 – 789. [40] R. Krishnan, J. S. Binkley, R. Seeger, J. A. Pople, J. Chem. Phys. 1980, 72, 650. [41] M. J. Frisch, J. A. Pople, J. S. Binkley, J. Chem. Phys. 1984, 80, 3265. [42] T. Clark, J. Chandrasekhar, G. W. Spitznagel, P. v. R. Schleyer, J. Comput. Chem. 1983, 4, 294. [43] W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, Jr., D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, P. A. Kollman, J. Am. Chem. Soc. 1995, 117, 5179 – 5197. [44] T. Dresselhaus, J. Yang, S. Kumbhar, M. P. Waller, J. Chem. Theory Comput. 2013, 9, 2137 – 2149. [45] T. Huber, W. F. van Gunsteren, J. Phys. Chem. A 1998, 102, 5937 – 5943. [46] A. Shmygelska, H. H. Hoos, BMC Bioinf. 2005, 6, 30. [47] S. Fidanova, I. Lirkov, Proceedings of the International Multiconference on Computer Science and Information Technology, IMCSIT 2008, 887 – 891. [48] M. Korth, J. Chem. Theory Comput. 2010, 6, 3808 – 3816. [49] J. J. P. Stewart, J. Mol. Model. 2013, 19, 1 – 32. [50] J. D. Smith, C. D. Cappa, K. R. Wilson, R. C. Cohen, P. L. Geissler, R. J. Saykally, Proc. Natl. Acad. Sci. USA 2005, 102, 14171 – 14174. [51] V. S. Bryantsev, M. S. Diallo, A. C. T. van Duin, W. A. Goddard III, J. Chem. Theory Comput. 2009, 5, 1016 – 1026. [52] A. Ben-Naim, Statistical Thermodynamics for Chemists and Biochemists, Plenum, New York, 1992. [53] D. Nerukh, J. Mol. Liq. 2012, 176, 65 – 70. [54] J. Sˇebek, B. Gyurcsik, J. Sˇebestk, Z. Kejk, L. Bednrov, P. Bourˇ, J. Phys. Chem. A 2007, 111, 2750 – 2760.
ChemPhysChem 0000, 00, 1 – 9
&7&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES [55] V. Sychrovsky´, M. Budeˇsˇnsky´, L. Benda, V. Sˇpirko, Z. Vokcˇov, J. Sˇebestk, P. Bourˇ, J. Phys. Chem. B 2008, 112, 1796 – 1805. [56] M. Knapp-Mohammady, K. J. Jalkanen, F. Nardi, R. C. Wade, S. Suhai, Chem. Phys. 1999, 240, 63 – 77.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.chemphyschem.org [57] P. G. Takis, K. D. Papavasileiou, L. D. Peristeras, V. S. Melissas, A. N. Troganis, Phys. Chem. Chem. Phys. 2013, 15, 7354 – 7362. Received: March 7, 2014 Published online on && &&, 2014
ChemPhysChem 0000, 00, 1 – 9
&8&
These are not the final page numbers! ÞÞ
ARTICLES M. P. Waller,* S. Kumbhar, J. Yang && – && A Density-Based Adaptive Quantum Mechanical/Molecular Mechanical Method Mechanical breakdown: A densitybased partitioning scheme is introduced for adaptive quantum mechanics/molecular mechanics (QM/MM). Unlike distance- and number-based methods, this method does not contain any fitting pa-
rameters, which is particularly important if the hydration number of a system is not known a priori. Thus, this densitybased QM/MM method is applicable to a wide range of modeling areas.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
ChemPhysChem 0000, 00, 1 – 9
&9&
These are not the final page numbers! ÞÞ
A Density-Based Adaptive Quantum Mechanical/Molecular Mechanical Method Mark P. Waller,* Sadhana Kumbhar, and Jack Yang[a] We present a density-based adaptive quantum mechanical/molecular mechanical (DBA-QM/MM) method, whereby molecules can switch layers from the QM to the MM region and vice versa. The adaptive partitioning of the molecular system ensures that the layer assignment can change during the optimization procedure, that is, on the fly. The switch from a QM molecule to a MM molecule is determined if there is an absence of noncovalent interactions to any atom of the QM core
region. The presence/absence of noncovalent interactions is determined by analysis of the reduced density gradient. Therefore, the location of the QM/MM boundary is based on physical arguments, and this neatly removes some empiricism inherent in previous adaptive QM/MM partitioning schemes. The DBA-QM/MM method is validated by using a water-in-water setup and an explicitly solvated l-alanyl–l-alanine dipeptide.
1. Introduction Modeling complex chemical systems in the gas or liquid phase, in the solid state, or on surfaces typically requires advanced mutliscale techniques.[1–3] Hybrid quantum mechanical/ molecular mechanical (QM/MM) methods partition a system into an active region that is treated with a high-level electronic structure method, and the surrounding environment can be inexpensively treated by a lower-level force field.[4, 5] One can therefore achieve good accuracy, while retaining efficiency, for studying complex chemical systems. For this reason, QM/MM methods have emerged as the method of choice for studying reaction mechanisms in large biomolecular systems. Furthermore, if certain parts of a system are either missing from the experimental data or if significant parts of the system have been modified, that is, by using chemical biology tools, then the whole system can be refined/remodeled with QM/MM methods.[6, 7] Multiscale methods for complex chemical systems are traditionally based on the rigid partitioning of the active region and the environment. To overcome this shortcoming and to investigate nonlocalized active regions more accurately, on-the-fly switching (i.e. adaptive QM/MM) is of great interest to the modeling community.[8–22] In adaptive schemes, the QM and MM assignments of molecules are not rigid, and molecules can therefore freely change their layers by repartitioning the system. Accurate modeling of complex chemical systems is largely reliant upon a solid experimental foundation, that is, X-ray [a] Dr. M. P. Waller, S. Kumbhar, Dr. J. Yang Theoretische Organische Chemie, Organisch-Chemisches Institut and Center for Multiscale Theory and Computation Westfalische Wilhelms-Universitt Mnster Corrensstrasse 40, 48149 Mnster (Germany) Fax: (+ 49) 0251 83 33202 E-mail: [email protected] Supporting Information for this article is available on the WWW under http://dx.doi.org/10.1002/cphc.201402105.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
crystallography or elegant NMR spectroscopy experiments to provide a starting structure. Traditional gradient-based optimizers converge to a local minimum, which is generally dependent upon the starting point.[23] In contrast, nonlocal optimization methods can reveal a variety of low-energy minima[24] and are far more explorative than traditional gradient-based geometry optimizers. Therefore, nonlocal optimizers are a good choice for de novo structure prediction of complex chemical systems. Metaheuristic algorithms nicely blend deterministic and stochastic aspects to generate diverse molecular structures. This is important if the chemical system is too complex to iterate through all possible conformations with a brute-force approach due to the high dimensionality of the search space. Nonlocal metaheuristic optimizers are generally not based on the gradient, and an entirely new solution is typically built at each metacycle of the algorithm, and this means that the energies do not need to be continuous. The traditional rigid partitioning of the QM and MM regions cannot be easily implemented in conjunction with metaheuristic methods. This is because neighboring fragments in any given cycle are not necessarily neighboring fragments from either the previous or subsequent optimization cycles. This is due to the highly explorative nature of the optimizer, which can place a fragment anywhere within a given search space, rather than the sequential application of short step lengths that are employed in conventional gradient-based optimizers (see Figure 1). Hence, there is a need for a method that can analyze the entire system at every metacycle, and this can be accomplished by employing an on-the-fly partitioning scheme (i.e. nonrigid partitioning of the entire system into either a QM or MM region). In general, partitioning the system into QM and MM regions is reliant upon a criterion, and to date, adaptive schemes can be categorized as either distance-based or ChemPhysChem 0000, 00, 1 – 9
&1&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
www.chemphyschem.org
zero, in the regions of covalent bonding and noncovalent interactions. These two regions can further be distinguished by using the density values; low density 1!0 combined with low s(1) represents the noncovalent interactions. As a consequence, this allows a spatial region W to be mapped from s(1) for constructing a visualizable RDG isosurface, V(r), in real space, which indicates the presence of noncovalent interactions between molecules (or molecular fragments). The spatial Figure 1. Schematic of a nonlocal optimizer (green) and a local optimizer (red) for the water dimer that share the extent of this RDG isosurface is same starting point. considered as a numerical indicator of the strength of the NCI between two interacting systems. This merit can be obtained number-based adaptive schemes (see Figure 2). The former from the integrated volume of RDG, defined as [Eq. (2)]: scheme defines the distance cutoff between the position of each molecule and the center of the QM core region, and the Z smax Z 1max latter one is defined by predefining the total number of molef ¼ dW ð2Þ 0 0 cules that are allowed in the QM region.[25] Setting an appropriate QM/MM boundary is crucial for successfully modeling moin which smax and 1max are set to threshold values, and f is the lecular systems with QM/MM methods. If the boundary to be integrated volume of the spatial area of the RDG isosurface set is too close to the active site, the chemical realism of the V(r). Hence, one can determine the presence or absence of model will be reduced. In contrast, if one positions the bounnoncovalent interactions by analyzing the intermolecular dendary too far away from the active site, the computational exsity using the NCI index. pense of the QM calculation may become problematic, possibly including difficulties in self-consistent field and geometry In summary, two QM/MM partitioning criteria were previousconvergence. Each molecular system is different in this respect, ly proposed (distance and number based), and herein we exand one must perform thorough tests to validate the appropriplored an alternative possibility of employing a nonempirical ateness of the particular QM/MM partitioning criteria used. partitioning criterion based on the NCI concept (see Figure 2). In 2010, Johnson et al.[26] introduced the concept of utilizing reduced charge density gradients (RDGs) to visually analyze Computational Section noncovalent interactions (NCI), and the NCIPLOT program was The DBA-QM/MM method was implemented into our own batchsubsequently released by Contreras-Garca et al.[27] The NCI based framework: JACOB.[32] concept was further developed to incorporate more quantifiable analysis of noncovalent interactions in fluctuating environPartitioning ments.[28] Experimental charge-density studies with the use of [29] [30, 31] models have independently the MOON and Multipole In principle, a density-based partitioning scheme requires a model confirmed the presence of RDGs in the crystalline environof the density for both the QM and the MM regions. The simplest ment. approach, and therefore the one adopted herein, is to assume atom-centered spherical densities for both regions. This spherical In short, according to the NCI concept, if two noncovalently density is also referred to as a promolecular density,[33] and it is interacting molecules are brought together, a finite RDG based a superposition of typically gas-phase neutral atomic densities. Of on Equation (1):
sð1Þ ¼
jr1j 4 2ð3p Þ 13 1
2
1 3
ð1Þ
where 1 is the total electron density and r1 is the electron density gradient, appears in the intermolecular region in which the total electronic density 1!0. To identify the intermolecular interactions, one can examine the plots of s(1) versus 1. The s(1) plots were shown to have very small values, approaching 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
course, the calculated/tabulated promolecular density for the MM region could be combined with an aspherical quantum mechanically derived charge density that naturally arises from the QM region. However, to expedite the numerical calculation, promolecular densities were used for both the QM and MM regions in the current implementation of our density-based partitioning scheme. The promolecular densities were shown to be capable of reproducing the topological features of a fully quantum mechanical calculation (i.e. aspherical density).[34, 35]
The NCIPLOT v1.0 program by Contreras-Garca et al.[27] was employed to analyze the intermolecular interactions by using the moChemPhysChem 0000, 00, 1 – 9
&2&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
www.chemphyschem.org Therefore, one can specifically investigate only the intermolecular interactions within supramolecular assemblies, and this approach was adopted for all systems in our present study. In our DBA-QM/MM method, any constituent of the system that must be treated quantum mechanically, for example, species involved in bond breaking/forming processes were placed into the static core region (QMcore). The region of space between the QMcore and each other molecule in the environment was analyzed for the presence of noncovalent interactions based on the NCI index. If, and only if, a noncovalent interaction was detected between the QMcore region and any given neighboring molecule, then that neighboring molecule was then included into the QM region. If no noncovalent interactions were observed to a given neighboring molecule, then this neighboring molecule was placed into the MM region, whereby the interaction between the QM and MM was treated by the classical electrostatic and van der Waals coupling terms. Repartitioning of the entire system throughout the optimization procedure was performed on the fly (see Figure 3).
Figure 3. The partitioning procedure constructs pairwise mutually independent (trivially parallelizable) calculations for NCIPLOT.[27]
Hybrid QM/MM The two-layer ONIOM method[36] within Gaussian 09[37] was used for the QM/MM single-point energy evaluations at the B3LYP[38, 39]/ 6-31G* or 6-311 + + G(2d,2p):[40–42] Amber[43] level of theory. These levels of theory were chosen, as they were utilized in the reference calculations, and they also represent common choices in the computational chemistry community. Furthermore, the same functional was used throughout this work because we were solely interested in investigating the effects of the partitioning scheme, which are not functional dependent. No healing region (alternatively called transition zone) was employed in this current study, and this had the distinct advantage of being computationally faster than methods that employ such a transition. The application of a nonlocal optimization technique in this current study meant that the healing region was unnecessary, see below.
Nonlocal Optimization Figure 2. Schematic partitioning of a small water cluster: QM (vdW radii)/ transition (ball and stick)/MM (thin wire). The green isosurfaces indicate the presence of a noncovalent interaction according to the NCI index.
lecular densities that were constructed for each monomer (i.e. 1molecular). The INTERMOLECULAR keyword was used so that only contributions from the intermolecular interactions would remain.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Our own Hybrid Metaheuristic Optimizer (HMO) method[44] was used to perform nonlocal optimization[24] of the molecular systems. The HMO method hybridizes a particle swarm optimizer (PSO)[45] with an ant colony optimizer (ACO),[46, 47] whereby the PSO optimizes the parameter space of the underlying ACO (see Table 1), whereas the ACO itself optimizes the conformational space of the ChemPhysChem 0000, 00, 1 – 9
&3&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
www.chemphyschem.org
Table 1. Parameters (optimized by using HMO) used for the ACO that were fitted to a water trimer by employing the PM6-DH + [48] Hamiltonian in MOPAC2012.[49] Parameter[a]
Optimized value
number of ants n pheromone resistance probability to choose global best global best percentage explorative/exploitive search controller trail width
20 0.25 0.30 0.30 1.00 0.01
[a] The presented parameters were used throughout this entire study.
molecular system. A supramolecular search space was used for each of the solvent water molecules (six degrees of freedom for each solvent molecule, three translational and three rotational) and the dihedral search space was used for the l-alanyl–l-alanine dipeptide (two degrees of freedom that is, f and y). The overall procedure for the DBA-QM/MM is depicted in the sequence diagram shown in Figure 4. The three-step process was run until the geometry converged or the maximum number of predetermined metacycles was reached.
Figure 5. A scan of the intercentroid distance of the water dimer showing the integrated volume from the NCI analysis of the promolecular density. Water dimers are overlaid schematically.
The orientational dependence of the intermolecular noncovalent interaction for the water dimer is plotted in Figure 6. Importantly, this particular dependence cannot be recovered 2. Results and Discussion by the distance-based adaptive scheme, and this nicely highWater in Water lights an additional advantage of our density-based approach for investigating nonspherically symmetric systems, that is, A correct description of water and solute–water clusters is critibeyond rare gas matrices. cal for the accurate prediction of structural and thermodynamFurthermore, if one investigates the atom dependency of ic properties. Clusters are formed due to attractive noncovalent the RDG isosurface, one sees that the trends in volumes of the interactions, and they are intermediate between the gas phase RDG reflects the same trends that are seen in the sum of the and the condensed phase; consequently, the interest in investivan der Waals (vdW) radii: O···O (RDG = 15 937 a.u., vdW = gating clusters has increased tremendously over the years. Al3.04 ), O···H (RDG = 14 702 a.u., vdW = 2.72 ), and H···H = though water clusters have been observed experimentally, in(RDG = 13 605 a.u., vdW = 2.40 ); all values were computed at vestigating supramolecular structures in water is difficult due a distance of 3.00 . to continuous H-bond breaking and forming.[50] To systematically examine the performance of the DBA-QM/ The aim of this section is to show a proof of principle for MM partitioning scheme, we optimized the water dimer, water the prototypical solvent, that is, water. An NCI analysis was trimer, and water tetramer. The QM/MM optimizations were performed for the water dimer, whereby the NCI isosurface disperformed over 20 independent runs to assess the reproduciappears if the intermolecular (center of masses) distance bebility. The search space was constructed in such a way that the comes too large (> 4 , see Figure 5). maximum distance between the centers of masses of molecules could be 14 and the search space granularity was set to 0.001 . All possible configurations from the density-based partitioning for the investigated water in water systems are shown in Figure 7. One can see that there are n possible configurations, for which n is the number of non-QMcore molecules. The distances between the centers of masses of the water molecules tended to decrease Figure 4. A sequence diagram showing the key features of the DBA-QM/MM method by using wrappers for NCIduring the nonlocal optimization PLOT[27] and G09.[37] 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
ChemPhysChem 0000, 00, 1 – 9
&4&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
www.chemphyschem.org l-Alanyl–l-Alanine (LALA) Peptide
Figure 6. A scan of the three rotational degree of freedom (a, b, g) for the water dimer by using an NCI analysis of the promolecular density at a constant intercentroid distance of 3 .
The effect of hydration on small peptides and amino acids is not completely understood because these effects are hard to verify with experiments.[52] Therefore, modeling can be used to increase our collective understanding. The zwitterionic form of the l-alanyl–l-alanine (LALAzw) peptide was studied as a model system, because the conformation of the molecule can be characterized by the f and y dihedral angles (see Figure 9).[53] Furthermore, water was shown to play an important role in the conformational changes of LALAzw due to the charged ends of the peptide. Finally, the predominant structure of this particular dipeptide in the pres-
Figure 7. Configurations possible for the dimer, trimer, and tetramer shown top to bottom, respectively.
procedure. The relative energies (DE, kcal mol1) of the pureQM- and DBA-QM/MM HMO optimized structures for each system with respect to the reference values[51] are shown in Figure 8. As the search space increased, so did the diversity of nonlocal minima for both the pure-QM reference and the DBAQM/MM method. Upon comparing the optimized structures, the water molecules were optimized into a nearly isoenergetic conformation from the previous reported structure. The results are very encouraging that our DBA-QM/MM method is capable of exploring multiple minima of comparable energy to the full QM reference optimizations on the explicitly solvated potential energy surface. 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Figure 8. Relative energies compared to those obtained by Bryantsev et al.[51] (DE, kcal mol1) computed for the pure-QM and DBA-QM/MM HMO optimized water dimer, trimer, and tetramer structures.
Figure 9. The zwitterionic form of the l-alanyl–l-alanine dipeptide.
ChemPhysChem 0000, 00, 1 – 9
&5&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES
Figure 10. Schematic of LALAzw with 14 water molecules, showing a possible configuration of the entire system.
ence of water is known. It has been shown that LALAzw (in explicit water) predominantly adopts a conformation with the two dihedral angles f 1258, y 1438 during classical molecular dynamics simulations.[53] Furthermore, the dihedral angles f= 1538 and y = 1478 are reported for the LALAzw structure in implicit solvation models.[54, 55] Knapp-Mohammady et al. showed that at least 4 water molecules are required to solvate the end groups of the zwitterionic form of LALA, and dihedral angles for LALAzw with 4, 6, 7, 10, 11, 12, 13, and 14 water molecules have been reported.[56] Therefore, these same systems (i.e. LALAzw with 4–14 explicit water molecules) were adopted to validate the DBA-QM/MM method (see Figure 10). Hybrid QM/MM calculations of the solvated LALAzw were performed at the B3LYP/6-31G* level of theory combined with the Amber force field. The calculations were performed by using the same basis set and functional as those used for the reference structures.[55] We compared our density-based partitioning approach with both the distance- and number-based partitioning approaches for each of the solvated LALAzw systems. It was recently shown that 8 1 water molecules directly interact with LALAzw.[57] Considering that the experimental hydration number is known, the number of solvent molecules around the solute was set to nine in the number-based-QM/ MM method. In the distance-based-QM/MM method, the boundary cutoff was set to 5 . Importantly, no transition zones were used for any of the adaptive partitioning schemes in this present study to directly compare each partitioning method. The median and standard deviations for f and y obtained for each method (i.e. all using the same nonlocal optimizer) were used to plot the normal distributions (see Figure 11). All of the dihedral angles obtained from the different adaptive schemes with different numbers of water molecules (i.e. varying from 4 to 14) are given in the Supporting Information. The DBA-QM/MM method performed in a manner similar to that of the other methods for the f dihedral angle. In the case of the y dihedral angle, the DBA-QM/MM method underestimated the angle to the same extent as the number-based method. The distance-based method performed slightly better 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.chemphyschem.org
Figure 11. Polar graph of the normal distributions of the dihedral angles (f and Y) for LALAzw by using different partitioning schemes over LALAzw plus 4–14 solvent molecules. Reference values are B3LYP/6-311 + + G**/ CPCM (star)[54] and BPW91/6-311 + + G**/PCM (square).[55]
than the DBA or number-based scheme for the y dihedral angle. In general, all three partitioning methods showed similar standard deviations. The pure QM method itself showed deviation across the dihedral angles over the different numbers of explicit water molecules included in the system. Overall, the similar results indicate that all of the methods are capable of returning more or less similar structures for an explicitly solvated LALAzw system. To ensure that our validation of the DBA-QM/MM method was robust with respect to traditional gradient-based optimizers, we also compared the nonlocal optimization results by using different adaptive methods to the reference dihedral angles that were taken from the literature[54, 55] (see Figure 11). Similar values for these dihedral angles were observed across the implicitly solvated gradient-based optimized structures; therefore, the validation of the DBA-QM/MM method is not biased by the choice of reference. In summary, the DBA-QM/MM method performs quite well with respect to the two other well-established methods of adaptive QM/MM partitioning for the explicitly solvated LALAzw system. This is encouraging, as our aim of this current study was not to present a method that necessarily outperformed the two pre-existing methods, but that, more importantly, was devoid of empirical parameters. This means that DBA-QM/MM can in principle be a black-box method that avoids laborious parameterization. The extension of the DBAQM/MM method to other solvated systems (i.e. different solutes and solvents) or to other environments such as surface phenomena is the focus of future work and will be reported in due course.
3. Conclusions A density-based adaptive hybrid QM/MM method able to optimize molecular systems was presented. In general, we showed that analysis of the NCI index was a valid approach to select ChemPhysChem 0000, 00, 1 – 9
&6&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES the boundary cut off for QM/MM calculations by using a number of test systems. The DBA partitioning scheme was shown to adequately describe the highly changing partitioning of the system in combination with the ONIOM subtractive QM/ MM Hamiltonian. Although a subtractive QM/MM Hamiltonian was employed within this work, the DBA-QM/MM partitioning scheme was equally applicable to an additive QM/MM scheme. The DBA QM/MM method was combined with a nonlocal optimization procedure, because continuous energies and gradients were not needed; hence, no healing regions were needed to remove discontinuities. In principle, the DBA-QM/MM method could alternatively be applied with molecular dynamics calculations, if smoothing functions were employed; however, the computational requirements would be significant. The major advantage of our density-based adaptive partitioning is that the method is general, and therefore, unlike the distance- and number-based methods, it does not need to be fitted upon investigating new chemical systems. This is particularly important if the hydration number of a system is not known a priori. Therefore, our nonparameterized density-based adaptive QM/MM method is applicable to a wide range of complex chemical systems.
Acknowledgements Generous financial support from the Deutsche Forschungsgemeinschaft is gratefully acknowledged. Keywords: adaptive · molecular modeling · noncovalent interactions · quantum chemistry · topology [1] N. Lpez, G. Pacchioni, F. Maseras, F. Illas, Chem. Phys. Lett. 1998, 294, 611. [2] E. Eichinger, P. Tavan, J. Hutter, M. Parrinello, J. Chem. Phys. 1999, 110, 10452 – 10467. [3] J. Gao, J. Chem. Phys. 1998, 109, 2346. [4] A. Warshel, M. Levitt, J. Mol. Biol. 1976, 103, 227 – 249. [5] H. M. Senn, W. Thiel, Angew. Chem. 2009, 121, 1220 – 1254; Angew. Chem. Int. Ed. 2009, 48, 1198 – 1229. [6] S. Kumbhar, S. Johannsen, R. K. O. Sigel, M. P. Waller, J. Mller, J. Inorg. Biochem. 2013, 127, 203 – 210. [7] T. Dresselhaus, N. D. Weikart, H. D. Mootz, M. P. Waller, RSC Adv. 2013, 3, 16122 – 16129. [8] T. Kerdcharoen, K. R. Liedl, B. M. Rode, Chem. Phys. 1996, 211, 313 – 323. [9] P. Chieux, H. Bertagnolli, J. Phys. Chem. 1984, 88, 3726 – 3730. [10] T. Kerdcharoen, K. Morokuma, Chem. Phys. Lett. 2002, 355, 257 – 262. [11] T. Kerdcharoen, K. Morokuma, J. Chem. Phys. 2003, 118, 8856 – 8862. [12] G. Csnyi, T. Albaret, M. C. Payne, A. De Vita, Phys. Rev. Lett. 2004, 93, 175503. [13] A. Heyden, H. Lin, D. G. Truhlar, J. Phys. Chem. B 2007, 111, 2231 – 2241. [14] B. Ensing, S. O. Nielsen, P. B. Moore, M. L. Klein, M. Parrinello, J. Chem. Theory Comput. 2007, 3, 1100 – 1105. [15] R. E. Bulo, B. Ensing, J. Sikkema, L. Visscher, J. Chem. Theory Comput. 2009, 5, 2212 – 2221. [16] H. Wang, C. Schtte, L. Delle Site, J. Chem. Theory Comput. 2012, 8, 2878 – 2887. [17] K. Park, A. W. Gçtz, R. C. Walker, F. Paesani, J. Chem. Theory Comput. 2012, 8, 2868 – 2877. [18] C. Rowley, B. Roux, J. Chem. Theory Comput. 2012, 8, 3526 – 3535. [19] B. Lev, B. Roux, S. Y. Noskov, J. Chem. Theory Comput. 2013, 9, 4165 – 4175.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.chemphyschem.org [20] N. Bernstein, C. Vrnai, I. Solt, S. A. Winfield, M. C. Payne, I. Simon, M. Fuxreiter, G. Csnyi, Phys. Chem. Chem. Phys. 2012, 14, 646 – 656. [21] S. Pezeshki, H. Lin, J. Chem. Theory Comput. 2011, 7, 3625 – 3634. [22] R. E. Bulo, C. Michel, P. Fleurat-Lessard, P. Sautet, J. Chem. Theory Comput. 2013, 9, 5567 – 5577. [23] H. B. Schlegel, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 790 – 809. [24] B. Hartke, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 879 – 887. [25] N. Takenaka, Y. Kitamura, Y. Koyano, M. Nagaoka, Chem. Phys. Lett. 2012, 524, 56 – 61. [26] E. R. Johnson, S. Keinan, P. Mori-Snchez, J. Contreras-Garca, A. J. Cohen, W. Yang, J. Am. Chem. Soc. 2010, 132, 6498 – 6506. [27] J. Contreras-Garca, W. Yang, E. R. Johnson, J. Phys. Chem. A 2011, 115, 12983 – 12990. [28] P. Wu, R. Chaudret, X. Hu, W. Yang, J. Chem. Theory Comput. 2013, 9, 2226 – 2234. [29] J. Yang, M. P. Waller, J. Comput. Chem. 2013, 34, 466 – 470. [30] G. Saleh, C. Gatti, L. Lo Presti, J. Contreras-Garca, Chem. Eur. J. 2012, 18, 15523 – 15536. [31] G. Saleh, C. Gatti, L. Lo Presti, Comput. Theor. Chem. 2012, 998, 148 – 163. [32] M. P. Waller, T. Dresselhaus, J. Yang, J. Comput. Chem. 2013, 34, 1420 – 1428. [33] P. Coppens, IUCr Texts on Crystallography, Vol. 4, Oxford University Press, Oxford, 1997. [34] M. A. Spackman, Chem. Phys. Lett. 1999, 301, 425 – 429. [35] C. Gatti, E. May, R. Destro, F. Cargoni, J. Phys. Chem. A 2002, 106, 2707 – 2720. [36] M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber, K. Morokuma, J. Phys. Chem. 1996, 100, 19357 – 19363. [37] Gaussian 09, Revision D.01, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, . Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2009.. [38] A. D. Becke, J. Chem. Phys. 1993, 98, 5648 – 5652. [39] C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785 – 789. [40] R. Krishnan, J. S. Binkley, R. Seeger, J. A. Pople, J. Chem. Phys. 1980, 72, 650. [41] M. J. Frisch, J. A. Pople, J. S. Binkley, J. Chem. Phys. 1984, 80, 3265. [42] T. Clark, J. Chandrasekhar, G. W. Spitznagel, P. v. R. Schleyer, J. Comput. Chem. 1983, 4, 294. [43] W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, Jr., D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, P. A. Kollman, J. Am. Chem. Soc. 1995, 117, 5179 – 5197. [44] T. Dresselhaus, J. Yang, S. Kumbhar, M. P. Waller, J. Chem. Theory Comput. 2013, 9, 2137 – 2149. [45] T. Huber, W. F. van Gunsteren, J. Phys. Chem. A 1998, 102, 5937 – 5943. [46] A. Shmygelska, H. H. Hoos, BMC Bioinf. 2005, 6, 30. [47] S. Fidanova, I. Lirkov, Proceedings of the International Multiconference on Computer Science and Information Technology, IMCSIT 2008, 887 – 891. [48] M. Korth, J. Chem. Theory Comput. 2010, 6, 3808 – 3816. [49] J. J. P. Stewart, J. Mol. Model. 2013, 19, 1 – 32. [50] J. D. Smith, C. D. Cappa, K. R. Wilson, R. C. Cohen, P. L. Geissler, R. J. Saykally, Proc. Natl. Acad. Sci. USA 2005, 102, 14171 – 14174. [51] V. S. Bryantsev, M. S. Diallo, A. C. T. van Duin, W. A. Goddard III, J. Chem. Theory Comput. 2009, 5, 1016 – 1026. [52] A. Ben-Naim, Statistical Thermodynamics for Chemists and Biochemists, Plenum, New York, 1992. [53] D. Nerukh, J. Mol. Liq. 2012, 176, 65 – 70. [54] J. Sˇebek, B. Gyurcsik, J. Sˇebestk, Z. Kejk, L. Bednrov, P. Bourˇ, J. Phys. Chem. A 2007, 111, 2750 – 2760.
ChemPhysChem 0000, 00, 1 – 9
&7&
These are not the final page numbers! ÞÞ
CHEMPHYSCHEM ARTICLES [55] V. Sychrovsky´, M. Budeˇsˇnsky´, L. Benda, V. Sˇpirko, Z. Vokcˇov, J. Sˇebestk, P. Bourˇ, J. Phys. Chem. B 2008, 112, 1796 – 1805. [56] M. Knapp-Mohammady, K. J. Jalkanen, F. Nardi, R. C. Wade, S. Suhai, Chem. Phys. 1999, 240, 63 – 77.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.chemphyschem.org [57] P. G. Takis, K. D. Papavasileiou, L. D. Peristeras, V. S. Melissas, A. N. Troganis, Phys. Chem. Chem. Phys. 2013, 15, 7354 – 7362. Received: March 7, 2014 Published online on && &&, 2014
ChemPhysChem 0000, 00, 1 – 9
&8&
These are not the final page numbers! ÞÞ
ARTICLES M. P. Waller,* S. Kumbhar, J. Yang && – && A Density-Based Adaptive Quantum Mechanical/Molecular Mechanical Method Mechanical breakdown: A densitybased partitioning scheme is introduced for adaptive quantum mechanics/molecular mechanics (QM/MM). Unlike distance- and number-based methods, this method does not contain any fitting pa-
rameters, which is particularly important if the hydration number of a system is not known a priori. Thus, this densitybased QM/MM method is applicable to a wide range of modeling areas.
2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
ChemPhysChem 0000, 00, 1 – 9
&9&
These are not the final page numbers! ÞÞ
