Next Article in Journal
Subcritical Water-Carbon Dioxide Pretreatment of Oil Palm Mesocarp Fiber for Xylooligosaccharide and Glucose Production
Next Article in Special Issue
Orientation of Laurdan in Phospholipid Bilayers Influences Its Fluorescence: Quantum Mechanics and Classical Molecular Dynamics Study
Previous Article in Journal
The Effect of Light Exposure at Night (LAN) on Carcinogenesis via Decreased Nocturnal Melatonin Synthesis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Parametrization of Combined Quantum Mechanical and Molecular Mechanical Methods: Bond-Tuned Link Atoms

Department of Chemistry, Chemical Theory Center, and Supercomputing Institute, University of Minnesota, Minneapolis, MN 55455-0431, USA
*
Authors to whom correspondence should be addressed.
Molecules 2018, 23(6), 1309; https://doi.org/10.3390/molecules23061309
Submission received: 25 April 2018 / Revised: 22 May 2018 / Accepted: 27 May 2018 / Published: 30 May 2018

Abstract

:
Combined quantum mechanical and molecular mechanical (QM/MM) methods are the most powerful available methods for high-level treatments of subsystems of very large systems. The treatment of the QM−MM boundary strongly affects the accuracy of QM/MM calculations. For QM/MM calculations having covalent bonds cut by the QM−MM boundary, it has been proposed previously to use a scheme with system-specific tuned fluorine link atoms. Here, we propose a broadly parametrized scheme where the parameters of the tuned F link atoms depend only on the type of bond being cut. In the proposed new scheme, the F link atom is tuned for systems with a certain type of cut bond at the QM−MM boundary instead of for a specific target system, and the resulting link atoms are call bond-tuned link atoms. In principle, the bond-tuned link atoms can be as convenient as the popular H link atoms, and they are especially well adapted for high-throughput and accurate QM/MM calculations. Here, we present the parameters for several kinds of cut bonds along with a set of validation calculations that confirm that the proposed bond-tuned link-atom scheme can be as accurate as the system-specific tuned F link-atom scheme.

Graphical Abstract

1. Introduction

The combined quantum mechanical and molecular mechanical (QM/MM) method is now well established for molecular simulations of large and complex chemical systems that are too big to be treated accurately by full QM methods [1,2,3,4,5,6]. The method combines the accuracy of QM methods and the efficiency of MM methods by treating a small-scale primary system at the QM level and a large-scale secondary system at MM level. Examples include catalytic systems where one needs to include effects beyond the active site and QM/MM methods have been widely applied to both enzyme kinetics [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22] and metal-organic frameworks (MOFs) [23,24,25,26,27,28,29,30,31].
The treatment of the QM−MM boundary is a challenging issue within the framework of QM/MM, especially when the boundary between the QM fragment and the MM surroundings passes through a covalent bond, which is unavoidable in treating enzymes and MOFs. A number of approaches have been proposed to cap the dangling bonds of the QM subsystem, such as a link atom or pseudoatom [32,33,34,35,36,37,38,39,40,41,42,43,44,45], or localized or generalized hybrid orbitals [46,47,48,49]. (Note that the H* method of Chaquin [39,41,42] was not used in the context of QM/MM calculations.) The link atom scheme is the subject of the present article.
Within the link-atom scheme, an H atom is the most widely used choice, with the QM−MM boundary usually set at C−C bonds because electrostatic balance is less of an issue at nonpolar bonds; nevertheless, a previous work [45] found that a F atom tuned for a particular structure of a particular system can perform better than H as a link atom in most cases, especially when the MM boundary atom is electronegative. The reason for choosing F as a tuned atom is that it is the most electronegative element and the preferred tuning procedure is to add a repulsive pseudopotential. Countervailing this success, the original tuned F link-atom scheme has some limitations. First, it is generally limited to QM/MM systems having just one kind of cut bond because there is no general scheme to tune F link atoms used to cap inequivalent cut bonds in the same fragment. Second, it is sometimes hard to choose a representative structure for the reaction of interest because multiple structures play an important role in the reaction. Third, the need for tuning is cumbersome for high-throughput screening. Here, we propose a new way to get the tuning parameters that overcomes all three of these drawbacks of the original scheme. The link atoms obtained by the new scheme are called bond-tuned link atoms. The new scheme is validated against a data set of 20 QM/MM systems involving 10 different cut bonds, as illustrated in the test set shown in Figure 1.

2. Method

A QM/MM calculation proceeds by dividing the entire system (ES) into two parts: a subsystem that will be treated by QM and a subsystem that will be treated by MM. In the present work, we assume that these subsystems are connected by a covalent bond. When the QM subsystem is pulled out of the ES, the bond is cut and the atom on the QM side of the bond (this atom is called Q1, see CO_1 in Figure 1 as an example) is unsaturated. A link atom is added to make a bond to Q1, and the QM subsystem combined with the link atom is called the capped QM system; the QM subsystem without the link atom is called the uncapped QM subsystem. The link atom is often taken to be a normal (nontuned) hydrogen atom, but in the present work, we take it to be a tuned fluorine atom. The definition of a tuned link atom is that it has a parameter varied so that the charge on the uncapped QM subsystem in a calculation on the capped QM system is the same as the charge on the uncapped QM subsystem in a calculation on the ES or a good model of the ES. The ES or model thereof that is used for tuning will be called the tuning system.
As previously described [45], to construct a tuned F atom (denoted as F*), the 1 s2 core of a real F atom is replaced by the sum of the CRENBL effective core potential (ECP) for F [52] plus the tuning pseudopotential U ( r ) . The tuning potential is [45]
U ( r ) = C   exp [ ( r / a 0 ) 2 ]
where a 0 is the Bohr radius and C (which will be given in hartrees) is the tuning parameter fitted to satisfy the definition stated above.
In applications (discussed below), the calculations on the capped QM system are carried out in the presence of background charges representing the electrostatic field due to the MM subsystem; this is called electronic embedding [53]. In previous work, we found that background charges have only a small effect on the final tuning parameter [50], and so the tuning was carried out without background charges. We followed the same protocol here.
In the tuned F link atom scheme, the tuning system is the specific system to be used in a particular QM/MM application [45]. An F atom tuned in this way is called a system-specific tuned F atom.
In the present work, we derive and validate the bond-tuned link atoms. Each of these tuned link atoms is parametrized for a particular kind of cut single bond at the QM−MM boundary, for example, where the Q1 atom is C and the other atom on the MM side of the cut bond (this atom is called M1, see CO_1 in Figure 1 as an example) is O. To derive the bond-tuned parameter for a certain kind of cut bond, we use (as the tuning system) a simple but representative molecule including that kind of bond. The present article derives parameters for 10 kinds of cut bonds, those in Figure 1, and the tuning system for each of these bond types is given in Table 1.
To evaluate the performance of the bond-tuned link atom scheme, we considered both H link atom and tuned F link atom schemes. The H and F* link atoms were placed along the axes of Q1−M1 bonds. The standard bond lengths listed in Table 2 were employed for Q1−H and Q1−F* link-atom bonds.
As described above, the tuning process requires calculating the charges on QM subsystems. We use the CM5 charge model [54], which is especially well suited for this purpose because of its good stability with respect to changing the basis set [51], that is, the tuning parameter is not strongly dependent on the basis set if using CM5 charges. In addition, CM5 charges can nicely reproduce the dipole moments of full quantum mechanical calculations [51]. A previous study [51] showed that “as compared to Mulliken charges, the CM5 charges describe the charge distributions in test molecules better, and they reproduce the dipole moments of full quantum mechanical calculations better.”

3. Details of the Validation Calculations

Although we turned background charges off during the tuning process, which is called mechanical embedding [53], background charges are actually considered in all of our QM/MM applications, i.e., the more realistic electronic embedding is used in the QM/MM calculations.
Electronic embedding calculations require specifying the method for treatment of boundary charges, i.e., the charges at the MM boundary. For the treatment of boundary charges in QM/MM calculations, we used the two previously recommended charge modification schemes [50], i.e., the balanced redistributed charge-2 (BRC2) scheme [45] and the balanced smeared redistributed charge (BSRC) scheme [50]. Here, we summarize these two charge modification schemes.
As the capped QM system has the same total charge as the original entire system, the charge on the entire QM/MM system (a capped QM system plus a MM subsystem) is not the same as the total charge on the original entire system. Therefore, the first step in both the BRC2 scheme and the BSRC scheme is to adjust the charge on the MM subsystem in order to conserve the total charge of the entire QM/MM system [45,50]. More specifically, the charge on the M1 atom is adjusted to make the total charge of the MM subsystem be zero [45,50].
Then, in the BRC2 scheme, the adjusted M1 charge is redistributed to all M2 atoms (i.e., the MM atoms directly bonded to M1 atoms, see CO_1 in Figure 1 as an example) evenly; in this scheme, all MM charges are point charges.
In the BSRC scheme, the adjusted M1 charge is redistributed to the midpoints of all M1−M2 bonds evenly. In this scheme, all MM charges are point charges except the redistributed charges, which are smeared as [50]
q R C * = q R C q R C   ( 1 + r / r 0 )   exp ( 2 r / r 0 )
where q R C is the redistributed charge, r is the distance of the charge density from the redistributed charge center, and r 0 is the smearing width of the charge density. Here, we use the previously recommended smearing width of 1 Å [50].
The additive QM/MM scheme with electronic embedding was adopted for QM/MM calculations; in this scheme the QM/MM energy of the entire system is given by [50,55]
E ( QM / MM ; E S ) = E ( Q M ; C Q M * * ) + [ E ( v a l ; E S ) E ( v a l ; C Q M ) ] + [ E ( v d W ; E S ) E ( v d W ; C Q M ) ] + E ( C o u l ; M M * * )
where E ( Q M ; C Q M * * ) is the quantum mechanical energy of the capped QM system embedded in the modified electrostatic field of the MM subsystem with adjusted M1 charge, the first bracketed energy difference is the difference in MM valence (val) energy terms between the entire system and the capped QM system, the second bracketed energy difference is the difference in MM van der Waals (vdW) energy terms between the entire system and the capped QM system, and E ( C o u l ; M M * * ) is the Coulomb (Coul) interaction energy of the MM subsystem with adjusted M1 charge. Note that the energy differences in brackets are independent of all decisions about electrostatics and charges.
To validate the robustness of the bond-tuned link atom scheme, we used the test suite that was used in previous works [50,51] plus five new and challenging molecules (i.e., CO_5, CO_6, CO_7, CO_8, and CC_4). The CO_5 and CO_6 test molecules were respectively constructed by replacing the two CH2 groups in the QM subsystem of CO_1 with two strong electron-withdrawing CF2 groups and by replacing the CH3 group in the MM subsystem of CO_1 with a strong electron-withdrawing CF3 group. The CO_7 and CO_8 test molecules have the same QM−MM boundary and MM subsystem as CO_6, but one (for CO_7) or two (for CO_8) more CH2 groups are included in the QM subsystem. The CC_4 and CO_7 test molecules have the same structure and they only differ in the position of QM−MM boundary (see Figure 1). All 20 test molecules are shown in Figure 1, and the deprotonation processes of these molecules are the test reactions used to evaluate the new scheme. We considered the deprotonation processes because they provide the most sensitive test of electrostatics because they involve a whole unit change in charge. If the method works well for the deprotonation processes, it should work well for other kinds of chemical processes involving less severe changes in charge distributions.
The expression for the deprotonation energy is
E ( deprotonation ) = E ( deprotonated   species ) E ( protonated   species )
The geometries of the protonated species were optimized at the full QM level (the Cartesian coordinates for these species can be obtained from the Supplementary Materials), while the deprotonated species were constructed by just removing a proton (the one marked with an asterisk in Figure 1) from the optimized protonated species without further optimization. Thus, we are considering what may be called vertical or sudden deprotonation, not adiabatic deprotonation. The reason for this is that the purpose of the present paper is to test the treatment of electrostatics for the demanding process of protonation/deprotonation in a way where the electrostatic effect is not obscured or partially cancelled or enhanced by geometry changes.
Determination of the tuning parameters for the bond-tuned scheme was carried out using the neutral molecules of Table 1. For testing the system-specific tuned F link atom scheme [45], tuning was performed using the protonated species of the molecules in the test suite; in addition, to be consistent with the tuning process of the bond-tuned link atom scheme, CM5 charges [54] were used and background charges were not considered while tuning. As stated above, a previous paper [50] reported that the effect of background charges on the tuning parameter is small.
All QM/MM calculations were carried out in our own QMMM 2017 program [56], which uses Gaussian 16 [57] as the QM engine and modified TINKER 6.3 [58] as the MM engine. All QM calculations were performed using Gaussian 16 [57]. The M06-2X density functional [59,60] and the 6-311G** [61] basis set were used for the calculations on the capped QM system and the entire system. The MMFF94 force field [62] was used for the valence and van der Waals terms of Equation (3). For the original MM charges (i.e., MM charges before applying the charge modification schemes), M06-2X/6-311G**/CM5 charges of the protonated species were used.
A previous paper in our group [51] has already tested the sensitivity of the tuning parameters and deprotonation energies to the basis set. That paper showed that the tuning parameters have at most minor changes with the variation of basis set; for different basis sets, the deprotonation energies are very similar. Thus, we only consider one basis set in the present paper.

4. Results and Discussion

4.1. Tuning Parameters

The tuning parameters C (in hartrees), obtained as described above, are listed in Table 3. In general, the system-specific tuning parameters are similar for molecules containing the same type of cut bond and show larger differences among molecules containing different types of cut bonds. The bond-tuned parameters are in qualitative agreement with the system-specific tuning parameters, although for CO_1, CO_4, CO_6, CO_7, CO_8, CN_1, CC_4, and OC_1 we observed relatively large deviations between the two sets of tuning parameters, with the system-specific tuning parameters typically smaller than those of the corresponding bond-tuned parameters. Equation (1) shows that for a tuning potential, a smaller tuning parameter C leads to a more electronegative tuned F atom; a positive (negative) tuning parameter means that a repulsive (attractive) tuning potential is added to a F atom. It is reasonable to infer that the molecules (CO_1, CO_4, CO_6, CO_7, CO_8, CN_1, CC_4, and OC_1) that have significantly smaller system-specific tuning parameters than bond-tuned parameters show this trend because they have strong electron-withdrawing groups (carbonyl and CF3) in the MM subsystems. Consistent with this hypothesis, we see that CO_5 has a stronger electron-withdrawing character than CO_1 in the QM region and has 0.2 hartree increase in the system-specific tuning parameter, while CO_6, which has stronger electron-withdrawing character in the MM region, has C decreased by 0.2 hartree. It is also consistent that CO_7 and CO_8, which have the same QM−MM boundary and MM subsystem as CO_6, have system-specific tuning parameters close to the corresponding parameter for CO_6. A final consistency check reported in Table 3 is that for the molecules that have C as the Q1 atom, but with different M1 atoms (O, N, C, S, Si), the bond-tuned parameters indicate that the electronegativities of the tuned F atoms have the order of C−Si < C−C < C−S < C−N < C−O. This order matches well with the electronegativity order of the MM boundary atoms, which lends support to the physicality of the derived bond-tuned parameters.

4.2. CM5 Charge Analyses

We computed the sum of the CM5 charges of the QM subsystems (i) of the entire molecules in the test suite; (ii) of the H-capped QM systems, and (iii) of the bond-tuned capped QM systems; the results are compared in Table 4. The table and Figure 2 show that the H-capped QM subsystem models fail to reproduce the total charges of the QM subsystems of the molecules, and the error can be as large as ~0.3 e. In contrast, we found that for most cases, the CM5 charges calculated from the QM subsystem capped with bond-tuned link atoms and those obtained from calculations on the entire system model are very similar, even for the molecules that have relatively large deviations of the tuning parameters (i.e., for CO_1, CO_4, CO_6, CO_7, CO_8, CN_1, CC_4, and OC_1, as discussed in Section 4.1). These results indicate that the bond-tuned link atoms perform significantly better than the H link atoms in terms of charges, and in fact, they accomplish the goal of making the charge on the QM subsystem realistic.

4.3. Deprotonation Energies

Table 5 shows the full QM deprotonation energies and the QM/MM signed errors and mean unsigned errors (MUEs) of the deprotonation energies (QM/MM deprotonation energies are given in the Supplementary Materials). For the QM/MM calculations, we found that for each link atom scheme, the two recommended charge modification schemes (i.e., BRC2 and BSRC) give very similar results. The key result in Table 5 is that the bond-tuned link atoms perform as well as the system-specific tuned F link atoms (MUE: 2.5 vs. 2.3 kcal/mol for the BRC2 scheme; 2.5 vs. 2.2 kcal/mol for the BSRC scheme) and much better than the H link atoms (MUE: 7.2 kcal/mol for the BRC2 scheme; 7.5 kcal/mol for the BSRC scheme). Moreover, for many molecules (e.g., OC_1), the bond-tuned link atoms even outperform the system-specific tuned F atoms in terms of QM/MM signed errors. These results show that using the bond-tuned link atom to cap the QM boundary atom can catch the main natures of the cut bond, and they indicate that the tuning parameter is transferable among molecules with the same type of cut bond. The bond-tuned link atom scheme is as straightforward and convenient as the popular H link atom scheme but as accurate as using the system-specific tuned F link atom scheme.
Table 5 also shows that the eight seemingly most challenging molecules (i.e., CO_1, CO_4, CO_6, CO_7, CO_8, CN_1, CC_4, and OC_1, as discussed above) all have similar errors between the system-specific and bond-tuned link atoms, with most (except CO_6) differences being 1.4 kcal/mol or less, though the two schemes have relatively large deviations on the tuning parameter for these molecules (see Table 3), which is very encouraging.
Compared to the system-specific tuned F link atom scheme, the existence of strong electron-withdrawing groups may introduce additional error in the bond-tuned link atom scheme. However, it is encouraging to see that the additional error is not very large (~1.2 kcal/mol in the case of CO_1, ~0.8 kcal/mol in the case of CO_4, and ~2.1 kcal/mol in the case of CO_6). In addition, the error of ~6 kcal/mol when using the bond-tuned link atom scheme in the case of CO_6 is still very small compared to the large deprotonation energy of 385.6 kcal/mol. Moreover, for CO_6, the bond-tuned link atom still performs significantly better than the H link atom. We also found that when moving the QM−MM boundary away from the reaction site (CO_6 → CO_7 → CO_8), the errors decrease significantly for all the link atom schemes. These results indicate that for QM/MM calculations employing the bond-tuned link atom scheme, increasing the size of the QM region is an effective way to reduce the error when strong electron-withdrawing (or donating) groups exist near the QM−MM boundary. We noticed that choosing the position of QM−MM boundary is also important. CC_4 and CO_7 have the same structure; they only differ in the position of QM−MM boundary, with a C−C (C−O) bond being cut for CC_4 (CO_7). Although CC_4 has a smaller QM region than CO_7, it has smaller errors than CO_7 for all the link atom schemes, which is mainly due to two reasons: (1) the strong electron-withdrawing groups are further away from the QM−MM boundary in CC_4 compared to CO_7; (2) usually setting the QM−MM boundary at a C−C bond instead of a C−O bond gives smaller errors.

5. Concluding Remarks

The introduction mentioned three limitations of the system-specific tuned F link-atom scheme. The bond-tuned scheme presented here overcomes all three of these limitations. It can be preparametrized with general parameters, and it can be applied to QM/MM systems including multiple cut bonds at the QM−MM boundary, for example, a MOF with the QM−MM boundary cutting through its inorganic node. A tuned F atom is just as convenient to use as the popular H link atom but is more accurate, and it has good stability [51] with respect to variation of basis sets because it is tuned with CM5 charges [54].
The overall performance (CM5 charges and deprotonation energies) of the proposed bond-tuned link-atom scheme is quite similar to the previous (system-specific) tuned F link-atom scheme [45]. The existence of strong electron-withdrawing groups near the QM−MM boundary may increase the error more when using the bond-tuned link atom scheme than when using the system-specific tuned F link-atom scheme, but the average error in the new method is still small.
Since the transferability of the tuning parameter among molecules with the same type of cut bond is validated in this work and the good stability of the tuning parameter with respect to basis set variations is validated in a previous work [51], future studies can focus on extending the database of bond-tuned parameters, which can serve for the high-throughput and accurate QM/MM calculations.

Supplementary Materials

The following are available online at https://www.mdpi.com/1420-3049/23/6/1309/s1. QM/MM deprotonation energies and Cartesian coordinates for molecules (optimized using M06-2X/6-311G**) in the test suite.

Author Contributions

Conceptualization, X.-P.W., L.G. and D.G.T.; Data curation, X.-P.W.; Funding acquisition, L.G. and D.G.T.; Investigation, X.-P.W.; Methodology, X.-P.W., L.G. and D.G.T.; Software, X.-P.W.; Supervision, L.G. and D.G.T.; Visualization, X.-P.W.; Writing—original draft, X.-P.W.; Writing—review & editing, L.G. and D.G.T.

Acknowledgments

This work was supported as part of the Inorganometallic Catalyst Design Center, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award DE-SC0012702. The authors acknowledge the Minnesota Supercomputing Institute (MSI) and National Energy Research Scientific Computing Center (NERSC) under Grant No. 90809 for providing computing resources.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Gao, J.; Thompson, M.A. (Eds.) Combined Quantum Mechanical and Molecular Mechanical Methods; ACS Symposium Series 712; American Chemical Society: Washington, DC, USA, 1998. [Google Scholar]
  2. Sherwood, P. Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; John von Neumann Institute for Computing: Jülich, Germany, 2000; p. 285. [Google Scholar]
  3. Lin, H.; Truhlar, D.G. QM/MM: What have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 2007, 117, 185–199. [Google Scholar] [CrossRef]
  4. Senn, H.M.; Thiel, W. QM/MM methods for biomolecular systems. Angew. Chem. Int. Ed. 2009, 48, 1198–1229. [Google Scholar] [CrossRef] [PubMed]
  5. Bernstein, N.; Kermode, J.R.; Csányi, G. Hybrid atomistic simulation methods for materials systems. Rep. Prog. Phys. 2009, 72, 026501. [Google Scholar] [CrossRef]
  6. Brunk, E.; Rothlisberger, U. Mixed quantum mechanical/molecular mechanical molecular dynamics simulations of biological systems in ground and electronically excited states. Chem. Rev. 2015, 115, 6217–6263. [Google Scholar] [CrossRef] [PubMed]
  7. Gao, J.; Truhlar, D.G. Quantum mechanical methods for enzyme kinetics. Annu. Rev. Phys. Chem. 2002, 53, 467–505. [Google Scholar] [CrossRef] [PubMed]
  8. Nam, K.; Gao, J.; York, D.M. Quantum mechanical/molecular mechanical simulation study of the mechanism of hairpin ribozyme catalysis. J. Am. Chem. Soc. 2008, 130, 4680–4691. [Google Scholar] [CrossRef] [PubMed]
  9. Lundberg, M.; Kawatsu, T.; Vreven, T.; Frisch, M.J.; Morokuma, K. Transition states in a protein environment—ONIOM QM:MM modeling of isopenicillin N synthesis. J. Chem. Theory Comput. 2009, 5, 222–234. [Google Scholar] [CrossRef] [PubMed]
  10. Rosta, E.; Woodcock, H.L.; Brooks, B.R.; Hummer, G. Artificial reaction coordinate ‘‘tunneling’’ in free-energy calculations: The catalytic reaction of RNase H. J. Comput. Chem. 2009, 30, 1634–1641. [Google Scholar] [CrossRef] [PubMed]
  11. Hu, H.; Yang, W. Development and application of ab initio QM/MM methods for mechanistic simulation of reactions in solution and in enzymes. J. Mol. Struct. Theochem. 2009, 898, 17–30. [Google Scholar] [CrossRef] [PubMed]
  12. Shaik, S.; Cohen, S.; Wang, Y.; Chen, H.; Kumar, D.; Thiel, W. P450 enzymes: Their structure, reactivity, and selectivity—Modeled by QM/MM calculations. Chem. Rev. 2010, 110, 949–1017. [Google Scholar] [CrossRef] [PubMed]
  13. Lonsdale, R.; Harvey, J.N.; Mulholland, A.J. Compound I reactivity defines alkene oxidation selectivity in cytochrome P450cam. J. Phys. Chem. B 2010, 114, 1156–1162. [Google Scholar] [CrossRef] [PubMed]
  14. Kanaan, N.; Ferrer, S.; Martí, S.; Garcia-Viloca, M.; Kohen, A.; Moliner, V. Temperature dependence of the kinetic isotope effects in thymidylate synthase: A theoretical study. J. Am. Chem. Soc. 2011, 133, 6692–6702. [Google Scholar] [CrossRef] [PubMed]
  15. Hou, G.; Cui, Q. QM/MM analysis suggests that alkaline phosphatase (AP) and nucleotide pyrophosphatase/phosphodiesterase slightly tighten the transition state for phosphate diester hydrolysis relative to solution: Implication for catalytic promiscuity in the AP superfamily. J. Am. Chem. Soc. 2012, 134, 229–246. [Google Scholar] [CrossRef] [PubMed]
  16. Layfield, J.P.; Hammes-Schiffer, S. Calculation of vibrational shifts of nitrile probes in the active site of ketosteroid isomerase upon ligand binding. J. Am. Chem. Soc. 2013, 135, 717–725. [Google Scholar] [CrossRef] [PubMed]
  17. Gómez, H.; Rojas, R.; Patel, D.; Tabak, L.A.; Lluch, J.M.; Masgrau, L. A computational and experimental study of O-glycosylation. catalysis by human UDP-GalNAc polypeptide:GalNAc transferase-T2. Org. Biomol. Chem. 2014, 12, 2645–2655. [Google Scholar] [CrossRef] [PubMed]
  18. Zhou, J.; Li, M.; Chen, N.; Wang, S.; Luo, H.-B.; Zhang, Y.; Wu, R. Computational design of a time-dependent histone deacetylase 2 selective inhibitor. ACS Chem. Biol. 2015, 10, 687–692. [Google Scholar] [CrossRef] [PubMed]
  19. Luk, L.Y.P.; Ruiz-Pernía, J.J.; Adesina, A.S.; Loveridge, E.J.; Tuñõn, I.; Moliner, V.; Allemann, R.K. Chemical ligation and isotope labeling to locate dynamic effects during catalysis by dihydrofolate reductase. Angew. Chem. Int. Ed. 2015, 54, 9016–9020. [Google Scholar] [CrossRef] [PubMed]
  20. Quesne, M.G.; Borowski, T.; De Visser, S.P. Quantum mechanics/molecular mechanics modeling of enzymatic processes: Caveats and breakthroughs. Chem. Eur. J. 2016, 22, 2562–2581. [Google Scholar] [CrossRef] [PubMed]
  21. Sousa, S.F.; Ribeiro, A.J.M.; Neves, R.P.P.; Brás, N.F.; Cerqueira, N.M.; Fernandes, P.A.; Ramos, M.J. Application of quantum mechanics/molecular mechanics methods in the study of enzymatic reaction mechanisms. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2017, 7, e1281. [Google Scholar] [CrossRef]
  22. Wang, B.; Johnston, E.M.; Li, P.; Shaik, S.; Davies, G.J.; Walton, P.H.; Rovira, C. QM/MM studies into the H2O2-dependent activity of lytic polysaccharide monooxygenases: Evidence for the formation of a caged hydroxyl radical intermediate. ACS Catal. 2018, 8, 1346–1351. [Google Scholar] [CrossRef]
  23. Choomwattana, S.; Maihom, T.; Khongpracha, P.; Probst, M.; Limtrakul, J. Structures and mechanisms of the carbonyl-ene reaction between MOF-11 encapsulated formaldehyde and propylene: An ONIOM study. J. Phys. Chem. C 2008, 112, 10855–10861. [Google Scholar] [CrossRef]
  24. Oxford, G.A.E.; Snurr, R.Q.; Broadbelt, L.J. Hybrid quantum mechanics/molecular mechanics investigation of (salen)Mn for use in metal−organic frameworks. Ind. Eng. Chem. Res. 2010, 49, 10965–10973. [Google Scholar] [CrossRef]
  25. Zheng, M.; Liu, Y.; Wang, C.; Liu, S.; Lin, W. Cavity-induced enantioselectivity reversal in a chiral metal-organic framework Brønsted acid catalyst. Chem. Sci. 2012, 3, 2623–2627. [Google Scholar] [CrossRef]
  26. Yu, D.; Yazaydin, A.O.; Lane, J.R.; Dietzel, P.D.C.; Snurr, R.Q. A combined experimental and quantum chemical study of CO2 adsorption in the metal-organic framework CPO-27 with different metals. Chem. Sci. 2013, 4, 3544–3556. [Google Scholar] [CrossRef]
  27. Hirao, H.; Ng, W.K.H.; Moeljadi, A.M.P.; Bureekaew, S. Multiscale model for a metal–organic framework: High-spin rebound mechanism in the reaction of the oxoiron(IV) species of Fe-MOF-74. ACS Catal. 2015, 5, 3287–3291. [Google Scholar] [CrossRef]
  28. Moeljadi, A.M.P.; Schmid, R.; Hirao, H. Dioxygen binding to Fe-MOF-74: Microscopic insights from periodic QM/MM calculations. Can. J. Chem. 2016, 94, 1144–1150. [Google Scholar] [CrossRef]
  29. Doitomi, K.; Xu, K.; Hirao, H. The mechanism of an asymmetric ring-opening reaction of epoxide with amine catalyzed by a metal-organic framework: Insights from combined quantum mechanics and molecular mechanics calculations. Dalton Trans. 2017, 46, 3470–3481. [Google Scholar] [CrossRef] [PubMed]
  30. Doitomi, K.; Hirao, H. Hybrid computational approaches for deriving quantum mechanical insights into metal–organic frameworks. Tetrahedron Lett. 2017, 58, 2309–2317. [Google Scholar] [CrossRef]
  31. Wu, X.-P.; Gagliardi, L.; Truhlar, D.G. Combined quantum mechanical and molecular mechanical method for metal–organic frameworks: Proton topologies of NU-1000. Phys. Chem. Chem. Phys. 2018, 20, 1778–1786. [Google Scholar] [CrossRef] [PubMed]
  32. Singh, U.C.; Kollman, P.A. A combined ab initio quantum mechanical and molecular mechanical method for carrying out simulations on complex molecular systems: Applications to the CH3Cl + Cl exchange reaction and gas phase protonation of polyethers. J. Comput. Chem. 1986, 7, 718–730. [Google Scholar] [CrossRef]
  33. Koga, N.; Morokuma, K. A simple scheme of estimating substitution or substituent effects in the ab initio MO method based on the shift operator. Chem. Phys. Lett. 1990, 172, 243–248. [Google Scholar] [CrossRef]
  34. Bakowies, D.; Thiel, W. Hybrid models for combined quantum mechanical and molecular mechanical approaches. J. Phys. Chem. 1996, 100, 10580–10594. [Google Scholar] [CrossRef]
  35. Zhang, Y.; Lee, T.-S.; Yang, W. A pseudobond approach to combining quantum mechanical and molecular mechanical methods. J. Chem. Phys. 1999, 110, 46–54. [Google Scholar] [CrossRef]
  36. Antes, I.; Thiel, W. Adjusted connection atoms for combined quantum mechanical and molecular mechanical methods. J. Phys. Chem. A 1999, 103, 9290–9295. [Google Scholar] [CrossRef]
  37. Alary, F.; Poteau, R.; Heully, J.-L.; Barthelat, J.-C.; Daudey, J.-P. A new method for modelling spectator chemical groups in ab initio calculations: Effective group potentials. Theor. Chem. Acc. 2000, 104, 174–178. [Google Scholar] [CrossRef]
  38. DiLabio, G.A.; Hurley, M.M.; Christiansen, P.A. Simple one-electron quantum capping potentials for use in hybrid QM/MM studies of biological molecules. J. Chem. Phys. 2002, 116, 9578–9584. [Google Scholar] [CrossRef]
  39. Dumont, E.; Chaquin, P. The H* method: Hydrogen atoms with a fictitious nuclear charge: A versatile theoretical tool for study of atom and group properties as substituents: Electronegativity and partition of σ and π contributions. J. Mol. Struct. THEOCHEM 2004, 680, 99–106. [Google Scholar] [CrossRef]
  40. Von Lilienfeld, O.A.; Tavernelli, I.; Rothlisberger, U.; Sebastiani, D. Variational optimization of effective atom centered potentials for molecular properties. J. Chem. Phys. 2005, 122, 014113. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. Dumont, E.; Chaquin, P. Diels-Alder reaction: A theoretical comprehensive study of substituent effects using the ‘H* method’. J. Mol. Struct. THEOCHEM 2006, 758, 161–167. [Google Scholar] [CrossRef]
  42. Dumont, E.; Chaquin, P. Investigation of pure inductive effects on benzene ring by 13C NMR chemical shifts: A theoretical study using fictitious nuclear charges of hydrogen atoms (‘H* method’). Chem. Phys. Lett. 2007, 435, 354–357. [Google Scholar] [CrossRef]
  43. Ohnishi, Y.Y.; Nakao, Y.; Sato, H.; Sakaki, S. Frontier orbital consistent quantum capping potential (FOC-QCP) for bulky ligand of transition metal complexes. J. Phys. Chem. A 2008, 112, 1946–1955. [Google Scholar] [CrossRef] [PubMed]
  44. Lewin, J.L.; Cramer, C.J. Modified carbon pseudopotential for use in ONIOM calculations of alkyl-substituted metallocenes. J. Phys. Chem. A 2008, 112, 12754–12760. [Google Scholar] [CrossRef] [PubMed]
  45. Wang, B.; Truhlar, D.G. Combined quantum mechanical and molecular mechanical methods for calculating potential energy surfaces: Tuned and balanced redistributed-charge algorithm. J. Chem. Theory Comput. 2010, 6, 359–369. [Google Scholar] [CrossRef] [PubMed]
  46. Théry, V.; Rinaldi, D.; Rivail, J.L.; Maigret, B.; Ferenczy, G.G. Quantum mechanical computations on very large molecular systems: The local self-consistent field method. J. Comput. Chem. 1994, 15, 269–282. [Google Scholar] [CrossRef]
  47. Gao, J.; Amara, P.; Alhambra, C.; Field, M.J. A generalized hybrid orbital (GHO) method for the treatment of boundary atoms in combined QM/MM calculations. J. Phys. Chem. A 1998, 102, 4714–4721. [Google Scholar] [CrossRef]
  48. Pu, J.; Gao, J.; Truhlar, D.G. Generalized hybrid orbital (GHO) method for combining ab initio Hartree−Fock wave functions with molecular mechanics. J. Phys. Chem. A 2004, 108, 632–650. [Google Scholar] [CrossRef]
  49. Monari, A.; Rivail, J.-L.; Assfeld, X. Theoretical modeling of large molecular systems. Advances in the local self consistent field method for mixed quantum mechanics/molecular mechanics calculations. Acc. Chem. Res. 2013, 46, 596–603. [Google Scholar] [CrossRef] [PubMed]
  50. Wang, B.; Truhlar, D.G. Geometry optimization using tuned and balanced redistributed charge schemes for combined quantum mechanical and molecular mechanical calculations. Phys. Chem. Chem. Phys. 2011, 13, 10556–10564. [Google Scholar] [CrossRef] [PubMed]
  51. Wang, B.; Truhlar, D.G. Tuned and balanced redistributed charge scheme for combined quantum mechanical and molecular mechanical (QM/MM) methods and fragment methods: Tuning based on the CM5 charge model. J. Chem. Theory Comput. 2013, 9, 1036–1042. [Google Scholar] [CrossRef] [PubMed]
  52. Pacios, L.F.; Christiansen, P.A. Ab initio relativistic effective potentials with spin-orbit operators. I. Li through Ar. J. Chem. Phys. 1985, 82, 2664–2671. [Google Scholar] [CrossRef]
  53. Antes, I.; Thiel, W. On the treatment of link atoms in hybrid methods. ACS Symp. Ser. 1998, 712, 50–65. [Google Scholar]
  54. Marenich, A.V.; Jerome, S.V.; Cramer, C.J.; Truhlar, D.G. Charge model 5: An extension of Hirshfeld population analysis for the accurate description of molecular interactions in gaseous and condensed phases. J. Chem. Theory Comput. 2012, 8, 527–541. [Google Scholar] [CrossRef] [PubMed]
  55. Lin, H.; Truhlar, D.G. Redistributed charge and dipole schemes for combined quantum mechanical and molecular mechanical calculations. J. Phys. Chem. A 2005, 109, 3991–4004. [Google Scholar] [CrossRef] [PubMed]
  56. Lin, H.; Zhang, Y.; Pezeshki, S.; Duster, A.W.; Wang, B.; Wu, X.-P.; Gagliardi, L.; Truhlar, D.G. QMMM 2017; University of Minnesota: Minneapolis, MN, USA, 2017. [Google Scholar]
  57. 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.; et al. Gaussian 16; Gaussian, Inc.: Wallingford, CT, USA, 2016. [Google Scholar]
  58. Ponder, J.W. TINKER, Version 6.3; Washington University: St. Louis, MO, USA, 2014. [Google Scholar]
  59. Zhao, Y.; Truhlar, D.G. Density functionals with broad applicability in chemistry. Acc. Chem. Res. 2008, 41, 157–167. [Google Scholar] [CrossRef] [PubMed]
  60. Zhao, Y.; Truhlar, D.G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: Two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215–241. [Google Scholar]
  61. Hehre, W.J.; Random, L.; Schleyer, P.v.R.; Pople, J.A. Ab Initio Molecular Orbital Theory; Wiley: New York, NY, USA, 1986. [Google Scholar]
  62. Halgren, T.A. Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. J. Comput. Chem. 1996, 17, 490–519. [Google Scholar] [CrossRef]
Sample Availability: Samples of the compounds are not available from the authors.
Figure 1. Test Suite. The asterisk (*) denotes a deprotonation site. The QM region is on the left of the cut bond, and the MM region is on its right. The Q1, M1, and M2 atoms of CO_1 are labeled. This test suite includes the whole test suite (15 molecules) that used in previous works [50,51] plus five new molecules, i.e., CO_5, CO_6, CO_7, CO_8, and CC_4.
Figure 1. Test Suite. The asterisk (*) denotes a deprotonation site. The QM region is on the left of the cut bond, and the MM region is on its right. The Q1, M1, and M2 atoms of CO_1 are labeled. This test suite includes the whole test suite (15 molecules) that used in previous works [50,51] plus five new molecules, i.e., CO_5, CO_6, CO_7, CO_8, and CC_4.
Molecules 23 01309 g001
Figure 2. Histogram for CM5 charge deviations of the capped-QM-system results from the entire-system results shown in Table 4.
Figure 2. Histogram for CM5 charge deviations of the capped-QM-system results from the entire-system results shown in Table 4.
Molecules 23 01309 g002
Table 1. Tuning molecules for various cut bonds (the atoms or fragments appearing before and after the dashes are in the QM and MM subsystems, respectively).
Table 1. Tuning molecules for various cut bonds (the atoms or fragments appearing before and after the dashes are in the QM and MM subsystems, respectively).
Cut bondC−OC−NC−CN−CO−C
MoleculeH3C−OHH3C−NH2H3C−CH3H2N−CH3HO−CH3
Cut bondC−SS−SS−CC−SiO−N
MoleculeH3C−SHHS−SHHS−CH3H3C−SiH3HO−NH2
Table 2. Standard bond lengths (Å) [45].
Table 2. Standard bond lengths (Å) [45].
BondC−HN−HO−HS−H
Distance1.091.010.951.34
BondC−F*N−F*O−F*S−F*
Distance1.331.411.411.65
Table 3. System-specific and bond-tuned parameters C (in hartrees) of the tuned F link atoms.
Table 3. System-specific and bond-tuned parameters C (in hartrees) of the tuned F link atoms.
System-SpecificBond-Tuned
CO_1−0.01850.2256
CO_20.20750.2256
CO_30.24340.2256
CO_4−0.03100.2256
CO_50.19540.2256
CO_6−0.21240.2256
CO_7−0.22830.2256
CO_8−0.23920.2256
CN_1−0.01220.3181
CC_10.68060.8463
CC_20.77150.8463
CC_30.73150.8463
CC_40.54990.8463
NC_11.07981.0966
OC_10.85991.0625
CS_10.60790.5888
SS_10.82560.8108
SC_10.97501.0658
CSi_10.75200.9001
ON_10.57660.6321
Table 4. Sum of partial atomic charges of all QM atoms using the CM5 charge model and deviations of the capped-QM-system results from the entire-system results.
Table 4. Sum of partial atomic charges of all QM atoms using the CM5 charge model and deviations of the capped-QM-system results from the entire-system results.
Entire SystemCapped QM SubsystemDeviations
H LinkBond-Tuned LinkH LinkBond-Tuned Link
CO_10.185−0.0830.130−0.27−0.06
CO_20.131−0.0850.127−0.220.00
CO_30.107−0.1030.111−0.210.00
CO_40.189−0.0840.131−0.27−0.06
CO_50.078−0.1220.071−0.20−0.01
CO_60.229−0.0820.131−0.31−0.10
CO_70.233−0.0830.132−0.32−0.10
CO_80.238−0.0800.135−0.32−0.10
CN_10.157−0.1120.081−0.27−0.08
CC_1−0.003−0.112−0.043−0.11−0.04
CC_20.003−0.085−0.014−0.09−0.02
CC_30.000−0.102−0.027−0.10−0.03
CC_40.054−0.085−0.015−0.14−0.07
NC_1−0.154−0.304−0.158−0.150.00
OC_1−0.109−0.342−0.157−0.23−0.05
CS_10.041−0.0840.046−0.130.01
SS_1−0.009−0.120−0.006−0.110.00
SC_1−0.049−0.120−0.073−0.07−0.02
CSi_10.007−0.086−0.028−0.09−0.04
ON_1−0.041−0.340−0.055−0.30−0.01
average −0.20−0.04
Table 5. QM deprotonation energies (DE), QM/MM signed errors and mean unsigned errors (MUEs) of deprotonation energies (all energies and errors in kcal/mol) for the test suite using the BRC2 and BSRC schemes with H link atoms and with system-specific and bond-tuned link atoms.
Table 5. QM deprotonation energies (DE), QM/MM signed errors and mean unsigned errors (MUEs) of deprotonation energies (all energies and errors in kcal/mol) for the test suite using the BRC2 and BSRC schemes with H link atoms and with system-specific and bond-tuned link atoms.
Molecule 1DEH LinkSystem-Specific Tuned F LinkBond-Tuned Link
BRC2BSRCBRC2BSRCBRC2BSRC
CO_1HOCH2CH2−OC(O)CH3393.711.411.62.93.14.14.3
CO_2HOCH2CH2−OCH2NH2399.68.49.40.61.70.71.8
CO_3HOOCCH2−OCHOHCH3365.85.06.6−1.00.6−1.00.5
CO_4HOCH2CH2CH2−OC(O)CH3396.87.37.41.61.82.42.6
CO_5HOCF2CF2−OC(O)CH3356.76.37.2−0.10.90.11.1
CO_6HOCH2CH2−OC(O)CF3385.612.512.44.03.86.15.9
CO_7HOCH2CH2CH2−OC(O)CF3391.08.28.12.42.23.73.6
CO_8HOCH2CH2CH2CH2−OC(O)CF3393.46.16.01.41.22.42.2
CN_1HOOCCH2−NHCOCH3354.79.09.4−1.1−0.50.61.1
CC_1HOOCCH2−CH2F360.84.64.7−1.9−1.8−0.9−0.9
CC_2HOCH2CH2−CHNH2CONH2404.26.26.20.30.30.80.7
CC_3HOCH2−CH2OH400.814.714.81.82.02.93.0
CC_4HOCH2CH2−CH2OC(O)CF3391.05.75.6−0.9−1.00.60.6
NC_1HOOCCH2NH−CH2CH2OH375.33.34.2−3.4−3.0−3.3−2.9
OC_1HOCH2CH2O−CH2CONH2397.83.64.0−7.1−7.0−5.7−5.6
CS_1HOCH2CH2−SCH3393.712.813.46.47.06.36.9
SS_1HOCH2CH2S−SCH3388.63.94.21.01.41.01.3
SC_1HOCH2CH2S−CH2CH2OH392.7−0.8−0.4−2.4−2.2−1.9−1.7
CSi_1HOCH2CH2−SiH2F394.57.76.41.60.22.41.1
ON_1HOCH2CH2O−N(CH3)2398.66.38.3−4.0−2.3−3.6−1.9
MUE 7.27.52.32.22.52.5
1 The QM subsystem is in bold font in the table.

Share and Cite

MDPI and ACS Style

Wu, X.-P.; Gagliardi, L.; Truhlar, D.G. Parametrization of Combined Quantum Mechanical and Molecular Mechanical Methods: Bond-Tuned Link Atoms. Molecules 2018, 23, 1309. https://doi.org/10.3390/molecules23061309

AMA Style

Wu X-P, Gagliardi L, Truhlar DG. Parametrization of Combined Quantum Mechanical and Molecular Mechanical Methods: Bond-Tuned Link Atoms. Molecules. 2018; 23(6):1309. https://doi.org/10.3390/molecules23061309

Chicago/Turabian Style

Wu, Xin-Ping, Laura Gagliardi, and Donald G. Truhlar. 2018. "Parametrization of Combined Quantum Mechanical and Molecular Mechanical Methods: Bond-Tuned Link Atoms" Molecules 23, no. 6: 1309. https://doi.org/10.3390/molecules23061309

APA Style

Wu, X. -P., Gagliardi, L., & Truhlar, D. G. (2018). Parametrization of Combined Quantum Mechanical and Molecular Mechanical Methods: Bond-Tuned Link Atoms. Molecules, 23(6), 1309. https://doi.org/10.3390/molecules23061309

Article Metrics

Back to TopTop