Accurate NMR Shieldings with σ-Functionals.

In recent years, density-functional methods relying on a new type of fifth-rung correlation functionals called σ-functionals have been introduced. σ-Functionals are technically closely related to the random phase approximation and require the same computational effort but yield distinctively higher accuracies for reaction and transition state energies of main group chemistry and even outperform double-hybrid functionals for these energies. In this work, we systematically investigate how accurate σ-functionals can describe nuclear magnetic resonance (NMR) shieldings. It turns out that σ-functionals yield very accurate NMR shieldings, even though in their optimization, exclusively, energies are employed as reference data and response properties such as NMR shieldings are not involved at all. This shows that σ-functionals combine universal applicability with accuracy. Indeed, the NMR shieldings from a σ-functional using input orbitals and eigenvalues from Kohn-Sham calculations with the exchange-correlation functional of Perdew, Burke and Ernzerhof (PBE) turned out to be the most accurate ones among the NMR shieldings calculated with various density-functional methods including methods using double-hybrid functionals. That σ-functionals can be used for calculating both reliable energies and response properties like NMR shieldings characterizes them as all-purpose functionals, which is appealing from an application point of view.

QM/MM Free Energy Calculations of Long-Range Biological Protonation Dynamics by Adaptive and Focused Sampling.

Water-mediated proton transfer reactions are central for catalytic processes in a wide range of biochemical systems, ranging from biological energy conversion to chemical transformations in the metabolism. Yet, the accurate computational treatment of such complex biochemical reactions is highly challenging and requires the application of multiscale methods, in particular hybrid quantum/classical (QM/MM) approaches combined with free energy simulations. Here, we combine the unique exploration power of new advanced sampling methods with density functional theory (DFT)-based QM/MM free energy methods for multiscale simulations of long-range protonation dynamics in biological systems. In this regard, we show that combining multiple walkers/well-tempered metadynamics with an extended system adaptive biasing force method (MWE) provides a powerful approach for exploration of water-mediated proton transfer reactions in complex biochemical systems. We compare and combine the MWE method also with QM/MM umbrella sampling and explore the sampling of the free energy landscape with both geometric (linear combination of proton transfer distances) and physical (center of excess charge) reaction coordinates and show how these affect the convergence of the potential of mean force (PMF) and the activation free energy. We find that the QM/MM-MWE method can efficiently explore both direct and water-mediated proton transfer pathways together with forward and reverse hole transfer mechanisms in the highly complex proton channel of respiratory Complex I, while the QM/MM-US approach shows a systematic convergence of selected long-range proton transfer pathways. In this regard, we show that the PMF along multiple proton transfer pathways is recovered by combining the strengths of both approaches in a QM/MM-MWE/focused US (FUS) scheme and reveals new mechanistic insight into the proton transfer principles of Complex I. Our findings provide a promising basis for the quantitative multiscale simulations of long-range proton transfer reactions in biological systems.

Efficient Exploitation of Numerical Quadrature with Distance-Dependent Integral Screening in Explicitly Correlated F12 Theory: Linear Scaling Evaluation of the Most Expensive RI-MP2-F12 Term.

We present a linear scaling atomic orbital based algorithm for the computation of the most expensive exchange-type RI-MP2-F12 term by employing numerical quadrature in combination with CABS-RI to avoid six-center-three-electron integrals. Furthermore, a robust distance-dependent integral screening scheme, based on integral partition bounds [Thompson, T. H.; Ochsenfeld, C. J. Chem. Phys. 2019, 150, 044101], is used to drastically reduce the number of the required three-center-one-electron integrals substantially. The accuracy of our numerical quadrature/CABS-RI approach and the corresponding integral screening is thoroughly assessed for interaction and isomerization energies across a variety of numerical integration grids. Our method outperforms the standard density fitting/CABS-RI approach with errors below 1 µEh even for small grid sizes and moderate screening thresholds. The choice of the grid size and screening threshold allows us to tailor our ansatz to a desired accuracy and computational efficiency. We showcase the approach's effectiveness for the chemically relevant system valinomycin, employing a triple-ζ F12 basis set combination (C54H90N6O18, 5757 AO basis functions, 10,266 CABS basis functions, 735,783 grid points). In this context, our ansatz achieves higher accuracy combined with a 135× speedup compared to the classical density fitting based variant, requiring notably less computation time than the corresponding RI-MP2 calculation. Additionally, we demonstrate near-linear scaling through calculations on linear alkanes. We achieved an 817-fold acceleration for C80H162 and an extrapolated 28,765-fold acceleration for C200H402, resulting in a substantially reduced computational time for the latter─from 229 days to just 11.5 min. Our ansatz may also be adapted to the remaining MP2-F12 terms, which will be the subject of future work.

Exploring Chemical Space Using Ab Initio Hyperreactor Dynamics.

In recent years, first-principles exploration of chemical reaction space has provided valuable insights into intricate reaction networks. Here, we introduce ab initio hyperreactor dynamics, which enables rapid screening of the accessible chemical space from a given set of initial molecular species, predicting new synthetic routes that can potentially guide subsequent experimental studies. For this purpose, different hyperdynamics derived bias potentials are applied along with pressure-inducing spherical confinement of the molecular system in ab initio molecular dynamics simulations to efficiently enhance reactivity under mild conditions. To showcase the advantages and flexibility of the hyperreactor approach, we present a systematic study of the method's parameters on a HCN toy model and apply it to a recently introduced experimental model for the prebiotic formation of glycinal and acetamide in interstellar ices, which yields results in line with experimental findings. In addition, we show how the developed framework enables the study of complicated transitions like the first step of a nonenzymatic DNA nucleoside synthesis in an aqueous environment, where the molecular fragmentation problem of earlier nanoreactor approaches is avoided.

Systematic QM/MM Study for Predicting 31P NMR Chemical Shifts of Adenosine Nucleotides in Solution and Stages of ATP Hydrolysis in a Protein Environment.

NMR (nuclear magnetic resonance) spectroscopy allows for important atomistic insights into the structure and dynamics of biological macromolecules; however, reliable assignments of experimental spectra are often difficult. Herein, quantum mechanical/molecular mechanical (QM/MM) calculations can provide crucial support. A major problem for the simulations is that experimental NMR signals are time-averaged over much longer time scales, and since computed chemical shifts are highly sensitive to local changes in the electronic and structural environment, sufficiently large averages over representative structural ensembles are essential. This entails high computational demands for reliable simulations. For NMR measurements in biological systems, a nucleus of major interest is 31P since it is both highly present (e.g., in nucleic acids) and easily observable. The focus of our present study is to develop a robust and computationally cost-efficient framework for simulating 31P NMR chemical shifts of nucleotides. We apply this scheme to study the different stages of the ATP hydrolysis reaction catalyzed by p97. Our methodology is based on MM molecular dynamics (MM-MD) sampling, followed by QM/MM structure optimizations and NMR calculations. Overall, our study is one of the most comprehensive QM-based 31P studies in a protein environment and the first to provide computed NMR chemical shifts for multiple nucleotide states in a protein environment. This study sheds light on a process that is challenging to probe experimentally and aims to bridge the gap between measured and calculated NMR spectroscopic properties.

Improved Sampling of Adaptive Path Collective Variables by Stabilized Extended-System Dynamics.

Because of the complicated multistep nature of many biocatalytic reactions, an a priori definition of reaction coordinates is difficult. Therefore, we apply enhanced sampling algorithms along with adaptive path collective variables (PCVs), which converge to the minimum free energy path (MFEP) during the simulation. We show how PCVs can be combined with the highly efficient well-tempered metadynamics extended-system adaptive biasing force (WTM-eABF) hybrid sampling algorithm, offering dramatically increased sampling efficiency due to its fast adaptation to path updates. For this purpose, we address discontinuities of PCVs that can arise due to path shortcutting or path updates with a novel stabilization algorithm for extended-system methods. In addition, we show how the convergence of simulations can be further accelerated by utilizing the multistate Bennett's acceptance ratio (MBAR) estimator. These methods are applied to the first step of the enzymatic reaction mechanism of pseudouridine synthases, where the ability of path WTM-eABF to efficiently explore intricate molecular transitions is demonstrated.

Highly accurate σ- and τ-functionals for beyond-RPA methods with approximate exchange kernels.

σ-Functionals are promising new developments for the Kohn-Sham correlation energy based upon the direct Random Phase Approximation (dRPA) within the adiabatic connection formalism, providing impressive improvements over dRPA for a broad range of benchmarks. However, σ-functionals exhibit a high amount of self-interaction inherited from the approximations made within dRPA. Inclusion of an exchange kernel in deriving the coupling-strength-dependent density-density response function leads to so-called τ-functionals, which - apart from a fourth-order Taylor series expansion - have only been realized in an approximate fashion so far to the best of our knowledge, most notably in the form of scaled σ-functionals. In this work, we derive, optimize, and benchmark three types of σ- and τ-functionals including approximate exchange effects in the form of an antisymmetrized Hartree kernel. These functionals, based on a second-order screened exchange type contribution in the adiabatic connection formalism, the electron-hole time-dependent Hartree-Fock kernel (eh-TDHF) otherwise known as RPA with exchange (RPAx), and an approximation thereof known as approximate exchange kernel (AXK), are optimized on the ASCDB database using two new parametrizations named A1 and A2. In addition, we report a first full evaluation of σ- and τ-functionals on the GMTKN55 database, revealing our exchange-including functionals to considerably outperform existing σ-functionals while being highly competitive with some of the best double-hybrid functionals of the original GMTKN55 publication. In particular, the σ-functionals based on AXK and τ-functionals based on RPAx with PBE0 reference stand out as highly accurate approaches for a wide variety of chemically relevant problems.

Analytical Second-Order Properties for the Random Phase Approximation: Nuclear Magnetic Resonance Shieldings.

A method for the analytical computation of nuclear magnetic resonance (NMR) shieldings within the direct random phase approximation (RPA) is presented. As a starting point, we use the RPA ground-state energy expression within the resolution-of-the-identity approximation in the atomic-orbital formalism. As has been shown in a recent benchmark study using numerical second derivatives [Glasbrenner, M. J. Chem. Theory Comput. 2022, 18, 192], RPA based on a Hartree-Fock reference shows accuracies comparable to coupled cluster singles and doubles (CCSD) for NMR chemical shieldings. Together with the much lower computational cost of RPA, it has emerged as an accurate method for the computation of NMR shieldings. Therefore, we aim to extend the applicability of RPA NMR to larger systems by introducing analytical second-order derivatives, making it a viable method for the accurate and efficient computation of NMR chemical shieldings.

Postsynthetic Transformation of Imine- into Nitrone-Linked Covalent Organic Frameworks for Atmospheric Water Harvesting at Decreased Humidity.

Herein, we report a facile postsynthetic linkage conversion method giving synthetic access to nitrone-linked covalent organic frameworks (COFs) from imine- and amine-linked COFs. The new two-dimensional (2D) nitrone-linked covalent organic frameworks, NO-PI-3-COF and NO-TTI-COF, are obtained with high crystallinity and large surface areas. Nitrone-modified pore channels induce condensation of water vapor at 20% lower humidity compared to their amine- or imine-linked precursor COFs. Thus, the topochemical transformation to nitrone linkages constitutes an attractive approach to postsynthetically fine-tune water adsorption properties in framework materials.

Efficient Method for the Computation of Frozen-Core Nuclear Gradients within the Random Phase Approximation.

A method for the evaluation of analytical frozen-core gradients within the random phase approximation is presented. We outline an efficient way to evaluate the response of the density of active electrons arising only when introducing the frozen-core approximation and constituting the main difficulty, together with the response of the standard Kohn-Sham density. The general framework allows to extend the outlined procedure to related electron correlation methods in the atomic orbital basis that require the evaluation of density responses, such as second-order Møller-Plesset perturbation theory or coupled cluster variants. By using Cholesky decomposed densities─which reintroduce the occupied index in the time-determining steps─we are able to achieve speedups of 20-30% (depending on the size of the basis set) by using the frozen-core approximation, which is of similar magnitude as for molecular orbital formulations. We further show that the errors introduced by the frozen-core approximation are practically insignificant for molecular geometries.

Fully Automated Generation of Prebiotically Relevant Reaction Networks from Optimized Nanoreactor Simulations.

The nanoreactor approach first introduced by the group of Martínez [Wang et al. Nat. Chem. 2014, 6, 1044-1048] has recently attracted much attention because of its ability to accelerate the discovery of reaction pathways. Here, we provide a comprehensive study of various simulation parameters and present an alternative implementation for the reactivity-enhancing spherical constraint function, as well as for the detection of reaction events. In this context, a fully automated postsimulation evaluation procedure based on RDKit and NetworkX analysis is introduced. The chemical and physical robustness of the procedure is examined by investigating the reactivity of selected homogeneous systems. The optimized procedure is applied at the GFN2-xTB level of theory to a system composed of HCN molecules and argon atoms, acting as a buffer, yielding prebiotically plausible primary and secondary precursors for the synthesis of RNA. Furthermore, the formose reaction network is explored leading to numerous sugar precursors. The discovered compounds reflect experimental findings; however, new synthetic routes and a large collection of exotic, highly reactive molecules are observed, highlighting the predictive power of the nanoreactor approach for unraveling the reactive manifold.

From free-energy profiles to activation free energies.

Given a chemical reaction going from reactant (R) to the product (P) on a potential energy surface (PES) and a collective variable (CV) discriminating between R and P, we define the free-energy profile (FEP) as the logarithm of the marginal Boltzmann distribution of the CV. This FEP is not a true free energy. Nevertheless, it is common to treat the FEP as the "free-energy" analog of the minimum potential energy path and to take the activation free energy, ΔFRP ‡, as the difference between the maximum at the transition state and the minimum at R. We show that this approximation can result in large errors. The FEP depends on the CV and is, therefore, not unique. For the same reaction, different discriminating CVs can yield different ΔFRP ‡. We derive an exact expression for the activation free energy that avoids this ambiguity. We find ΔFRP ‡ to be a combination of the probability of the system being in the reactant state, the probability density on the dividing surface, and the thermal de Broglie wavelength associated with the transition. We apply our formalism to simple analytic models and realistic chemical systems and show that the FEP-based approximation applies only at low temperatures for CVs with a small effective mass. Most chemical reactions occur on complex, high-dimensional PES that cannot be treated analytically and pose the added challenge of choosing a good CV. We study the influence of that choice and find that, while the reaction free energy is largely unaffected, ΔFRP ‡ is quite sensitive.

Accelerating Hybrid Density Functional Theory Molecular Dynamics Simulations by Seminumerical Integration, Resolution-of-the-Identity Approximation, and Graphics Processing Units.

The computationally very demanding evaluation of the 4-center-2-electron (4c2e) integrals and their respective integral derivatives typically represents the major bottleneck within hybrid Kohn-Sham density functional theory molecular dynamics simulations. Building upon our previous works on seminumerical exact-exchange (sn-LinK) [Laqua, H., Thompsons, T. H., Kussmann, J., Ochsenfeld, C., J. Chem. Theory Comput. 2020, 16, 1465] and resolution-of-the-identity Coulomb (RI-J) [Kussmann, J., Laqua, H., Ochsenfeld, C., J. Chem. Theory Comput. 2021, 17, 1512], the expensive 4c2e integral evaluation can be avoided entirely, resulting in a highly efficient electronic structure theory method, allowing for fast ab initio molecular dynamics (AIMD) simulations even with large basis sets. Moreover, we propose to combine the final self-consistent field (SCF) step with the subsequent nuclear forces evaluation, providing the forces at virtually no additional cost after a converged SCF calculation, reducing the total runtime of an AIMD simulation by about another 25%. In addition, multiple independent MD trajectories can be computed concurrently on a single node, leading to a greatly increased utilization of the available hardware─especially when combined with graphics processing unit acceleration─improving the overall throughput by up to another 5 times in this way. With all of those optimizations combined, our proposed method provides nearly 3 orders of magnitude faster execution times than traditional 4c2e integral-based methods. To demonstrate the practical utility of the approach, quantum-mechanical/molecular-mechanical dynamics simulations on double-stranded DNA were performed, investigating the relative hydrogen bond strength between adenine-thymine and guanine-cytosine base pairs. In addition, this illustrative application also contains a general accuracy assessment of the introduced approximations (integration grids, resolution-of-the-identity) within AIMD simulations, serving as a protocol on how to apply these new methods to practical problems.

Tensor-Hypercontracted MP2 First Derivatives: Runtime and Memory Efficient Computation of Hyperfine Coupling Constants.

We employ our recently introduced tensor-hypercontracted (THC) second-order Møller-Plesset perturbation theory (MP2) method [Bangerter, F. H., Glasbrenner, M., Ochsenfeld, C. J. Chem. Theory Comput. 2021, 17, 211-221] for the computation of hyperfine coupling constants (HFCCs). The implementation leverages the tensor structure of the THC factorized electron repulsion integrals for an efficient formation of the integral-based intermediates. The computational complexity of the most expensive and formally quintic scaling exchange-like contribution is reduced to effectively subquadratic, by making use of the intrinsic, exponentially decaying coupling between tensor indices through screening based on natural blocking. Overall, this yields an effective subquadratic scaling with a low prefactor for the presented THC-based AO-MP2 method for the computation of isotropic HFCCs on DNA fragments with up to 500 atoms and 5000 basis functions. Furthermore, the implementation achieves considerable speedups with up to a factor of roughly 600-1000 compared to previous implementations [Vogler, S., Ludwig, M., Maurer, M., Ochsenfeld, C. J. Chem. Phys. 2017, 147, 024101] for medium-sized organic radicals, while also significantly reducing storage requirements.

Development of hetero-triaryls as a new chemotype for subtype-selective and potent Sirt5 inhibition.

In contrast to other sirtuins (NAD+-dependent class III lysine deacylases), inhibition of Sirt5 is poorly investigated, yet. Our present work is based on the recently identified Sirt5 inhibitor balsalazide, an approved drug with negligible bioavailability after oral administration. After gaining first insights into its structure-activity relationship in previous work, we were able to now develop heteroaryl-triaryls as a novel chemotype of drug-like, potent and subtype-selective Sirt5 inhibitors. The unfavourable azo group of the lead structure was modified in a systematic and comprehensive manner, leading us to a few open-chained and, most importantly, five-membered heteroaromatic substitutes (isoxazole CG_209, triazole CG_220, pyrazole CG_232) with very encouraging in vitro activities (IC50 on Sirt5 in the low micromolar range, <10 µM). These advanced inhibitors were free of cytotoxicity and showed favourable pharmacokinetic properties, as confirmed by permeability into mitochondria using live cell imaging experiments. Furthermore, results from calculations of the relative free binding affinities of the analogues compared to balsalazide as reference compound agreed well with the trends for inhibitory activities obtained in the in vitro experiments. Therefore, this method can be used to predict the affinity of closely related future potential Sirt5 inhibitors. Encouraged by our findings, we employed chemoproteomic selectivity profiling to confirm Sirt5 as main target of balsalazide and one of its improved analogues. An immobilised balsalazide-analogue specifically pulled down Sirt5 from whole cell lysates and competition experiments identified glutaryl-CoA dehydrogenase (GCDH) and nucleotide diphosphate kinase (NME4) as potential off-targets, once again confirming the selectivity of the novel balsalazide-derived Sirt5 inhibitors. In summary, a combination of targeted chemical synthesis, biological work, and computational studies led to a new generation of tailored Sirt5 inhibitors, which represent valuable chemical tools for the investigation of the physiological role of Sirt5, but could also serve as advanced lead structures for drug candidates for systemic use.

Statistically optimal analysis of the extended-system adaptive biasing force (eABF) method.

The extended-system adaptive biasing force (eABF) method and its newer variants offer rapid exploration of the configuration space of chemical systems. Instead of directly applying the ABF bias to collective variables, they are harmonically coupled to fictitious particles, which separates the problem of enhanced sampling from that of free energy estimation. The prevalent analysis method to obtain the potential of mean force (PMF) from eABF is thermodynamic integration. However, besides the PMF, most information is lost as the unbiased probability of visited configurations is never recovered. In this contribution, we show how statistical weights of individual frames can be computed using the Multistate Bennett's Acceptance Ratio (MBAR), putting the post-processing of eABF on one level with other frequently used sampling methods. In addition, we apply this formalism to the prediction of nuclear magnetic resonance shieldings, which are very sensitive to molecular geometries and often require extensive sampling. The results show that the combination of enhanced sampling by means of extended-system dynamics with the MBAR estimator is a highly useful tool for the calculation of ensemble properties. Furthermore, the extension of the presented scheme to the recently published Gaussian-accelerated molecular dynamics eABF hybrid is straightforward and approximation free.

Highly Efficient and Accurate Computation of Multiple Orbital Spaces Spanning Fock Matrix Elements on Central and Graphics Processing Units for Application in F12 Theory.

We employ our recently published highly efficient seminumerical exchange (sn-LinK) [Laqua, H.; Thompson, T. H.; Kussmann, J.; Ochsenfeld, C. J. Chem. Theory Comput. 2020, 16, 1456-1468] and integral-direct resolution of the identity Coulomb (RI-J) [Kussmann, J.; Laqua, H.; Ochsenfeld, C. J. Chem. Theory Comput. 2021, 17, 1512-1521] methods to significantly accelerate the computation of the demanding multiple orbital spaces spanning Fock matrix elements present in R12/F12 theory on central and graphics processing units. The errors introduced by RI-J and sn-LinK into the RI-MP2-F12 energy are thoroughly assessed for a variety of basis sets and integration grids. We find that these numerical errors are always below "chemical accuracy" (∼1 mH) even for the coarsest settings and can easily be reduced below 1 µH by employing only moderately large integration grids and RI-J basis sets. Since the number of basis functions of the multiple orbital spaces is notably larger compared with conventional Hartree-Fock theory, the efficiency gains from the superior basis scaling of RI-J and sn-LinK (O(Nbas2) instead of O(Nbas4) for both) are even more significant, with maximum speedup factors of 37 000 for RI-J and 4500 for sn-LinK. In total, the multiple orbital spaces spanning Fock matrix evaluation of the largest tested structure using a triple-ζ F12 basis set (5058 AO basis functions, 9267 CABS basis functions) is accelerated over 1575× using CPUs and over 4155× employing GPUs.

Exponential averaging versus umbrella sampling for computing the QM/MM free energy barrier of the initial step of the desuccinylation reaction catalyzed by sirtuin 5.

The computational characterization of enzymatic reactions poses a great challenge which arises from the high dimensional and often rough potential energy surfaces commonly explored by static QM/MM methods such as adiabatic mapping (AM). The present study highlights the difficulties in estimating free energy barriers via exponential averaging over AM pathways. Based on our previous study [von der Esch et al., J. Chem. Theory Comput., 2019, 15, 6660-6667], where we analyzed the first reaction step of the desuccinylation reaction catalyzed by human sirtuin 5 (SIRT5) by means of QM/MM adiabatic mapping and machine learning, we use, here, umbrella sampling to compute the free energy profile of the initial reaction step. The computational investigations show that the initial step of the desuccinylation reaction proceeds via an SN2-type reaction mechanism in SIRT5, suggesting that the first step of the deacylation reactions catalyzed by sirtuins is highly conserved. In addition, the direct comparison of the extrapolated free energy barrier from minimal energy paths and the computed free energy path from umbrella sampling further underlines the importance of extensive sampling.

How to obtain reaction free energies from free-energy profiles.

For chemical reactions that occur via the rearrangement of atoms from a configuration about one minimum (reactant, R) of the potential energy surface (PES) to a configuration about another minimum (product, P), an exact relation between the Helmholtz reaction free energy (ΔFRP) and the free-energy profile (FEP) can be derived. Since the FEP assumes a form similar to that of the PES along the minimum energy path between R and P, there is an unfortunate tendency to regard the FEP as the "free-energy" analog of the minimum energy path and consequently to equate ΔFRP to the difference between the values of the FEP at the minima corresponding to R and P. Analytic treatments of one- and two-dimensional models are presented that show how this mistaken idea leads to errors. In effect, treating the FEP by analogy with the minimum energy path neglects the role of entropy. The FEP is a function of a collective variable (CV), which must be chosen to describe the course of the rearrangement consistently with the exact relation between ΔFRP and the FEP. For large systems of common interest, the PES is often so complex that a straightforward way of choosing a CV is lacking. Consequently, one is forced to make an educated guess. A criterion for judging the quality of the guess is proposed and applied to a two-dimensional model.

Factors That Determine the Variation of Equilibrium and Kinetic Properties of QM/MM Enzyme Simulations: QM Region, Conformation, and Boundary Condition.

To analyze the impact of various technical details on the results of quantum mechanical (QM)/molecular mechanical (MM) enzyme simulations, including the QM region size, catechol-O-methyltransferase (COMT) is studied as a model system using an approximate QM/MM method (DFTB3/CHARMM). The results show that key equilibrium and kinetic properties for methyl transfer in COMT exhibit limited variations with respect to the size of the QM region, which ranges from ∼100 to ∼500 atoms in this study. With extensive sampling, local and global structural characteristics of the enzyme are largely conserved across the studied QM regions, while the nature of the transition state (e.g., secondary kinetic isotope effect) and reaction exergonicity are largely maintained. Deviations in the free energy profile with different QM region sizes are similar in magnitude to those observed with changes in other simulation protocols, such as different initial enzyme conformations and boundary conditions. Electronic structural properties, such as the covariance matrix of residual charge fluctuations, appear to exhibit rather long-range correlations, especially when the peptide backbone is included in the QM region; this observation holds when a range-separated DFT approach is used as the QM region, suggesting that delocalization error is unlikely the origin. Overall, the analyses suggest that multiple simulation details determine the results of QM/MM enzyme simulations with comparable contributions.