An extensible interface for QM/MM molecular dynamics simulations with AMBER

We present an extensible interface between the AMBER molecular dynamics (MD) software package and electronic structure software packages for quantum mechanical (QM) and mixed QM and classical molecular mechanical (MM) MD simulations within both mechanical and electronic embedding schemes. With this interface, ab initio wave function theory and density functional theory methods, as available in the supported electronic structure software packages, become available for QM/MM MD simulations with AMBER. The interface has been written in a modular fashion that allows straight forward extensions to support additional QM software packages and can easily be ported to other MD software. Data exchange between the MD and QM software is implemented by means of files and system calls or the message passing interface standard. Based on extensive tests, default settings for the supported QM packages are provided such that energy is conserved for typical QM/MM MD simulations in the microcanonical ensemble. Results for the free energy of binding of calcium ions to aspartate in aqueous solution comparing semiempirical and density functional Hamiltonians are shown to demonstrate features of this interface. © 2013 Wiley Periodicals, Inc.     

QuanPol: A full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree–Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self‐consistent field method, Møller–Plesset perturbation theory method, and time‐dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self‐consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations‐All Atom (OPLS‐AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented. © 2013 Wiley Periodicals, Inc.     

Reaching multi‐nanosecond timescales in combined QM/MM molecular dynamics simulations through parallel horsetail sampling

We report an enhanced sampling technique that allows to reach the multi‐nanosecond timescale in quantum mechanics/molecular mechanics molecular dynamics simulations. The proposed technique, called horsetail sampling, is a specific type of multiple molecular dynamics approach exhibiting high parallel efficiency. It couples a main simulation with a large number of shorter trajectories launched on independent processors at periodic time intervals. The technique is applied to study hydrogen peroxide at the water liquid–vapor interface, a system of considerable atmospheric relevance. A total simulation time of a little more than 6 ns has been attained for a total CPU time of 5.1 years representing only about 20 days of wall‐clock time. The discussion of the results highlights the strong influence of the solvation effects at the interface on the structure and the electronic properties of the solute. © 2017 Wiley Periodicals, Inc.     

Methodological aspects of QM/MM calculations: A case study on matrix metalloproteinase‐2

We address methodological issues in quantum mechanics/molecular mechanics (QM/MM) calculations on a zinc‐dependent enzyme. We focus on the first stage of peptide bond cleavage by matrix metalloproteinase‐2 (MMP‐2), that is, the nucleophilic attack of the zinc‐coordinating water molecule on the carbonyl carbon atom of the scissile fragment of the substrate. This step is accompanied by significant charge redistribution around the zinc cation, bond cleavage, and bond formation. We vary the size and initial geometry of the model system as well as the computational protocol to demonstrate the influence of these choices on the results obtained. We present QM/MM potential energy profiles for a set of snapshots randomly selected from QM/MM‐based molecular dynamics simulations and analyze the differences in the computed profiles in structural terms. Since the substrate in MMP‐2 is located on the protein surface, we investigate the influence of the thickness of the water layer around the enzyme on the QM/MM energy profile. Thin water layers (0–2 Å) give unrealistic results because of structural reorganizations in the active‐site region at the protein surface. A 12 Å water layer appears to be sufficient to capture the effect of the solvent; the corresponding QM/MM energy profile is very close to that obtained from QM/MM/SMBP calculations using the solvent macromolecular boundary potential (SMBP). We apply the optimized computational protocol to explain the origin of the different catalytic activity of the Glu116Asp mutant: the energy barrier for the first step is higher, which is rationalized on structural grounds. © 2016 Wiley Periodicals, Inc.     

Application of QM simulations and multivariate analysis in the study of alkene reactivity in the zeolite H‐ZSM5

Reported herein are the results of an investigation into the effect of the extended framework of the zeolite ZSM‐5 on the reaction energetics and structures of (a) the physisorbed complex formed between the zeolite and six alkenes, (b) the corresponding chemisorbed alkoxide intermediate and (c) the transition states (TS) connecting the two. For this, quantum mechanical (QM) simulations of ZSM‐5 in the presence and absence of the zeolite framework have been employed. A 46T density functional theory (DFT) cluster model and a 3T:46T DFT:UFF ONIOM model are used to represent the former scenario and a simple 3T DFT cluster model for the latter. The structural implications of neglecting the zeolite framework have been rigorously compared using the multivariate statistical method principal components analysis (PCA). This method allows one to assess the correlated nature of the changes in structure along the reaction coordinate, for multiple different alkenes, in a facile, reliable way. Copyright © 2008 John Wiley & Sons, Ltd.     

An embedded QM/MM study for different SiO2/Al2O3 ratios of the HZSM‐5 zeolite and for their interaction with n‐heptane

We carried out a theoretical study of the HZSM‐5 zeolite, for different SiO₂/Al₂O₃ ratios, that interacts with the n‐heptane molecule. The study was performed using a QM/MM (quantum mechanics/molecular mechanics) methodology. For the QM part, we have chosen a hybrid Hartree‐Fock density functional theory (DFT). The hybrid ACM/DZP approach, as implemented in Turbomole, was used for the treatment of the QM cluster containing 84 atoms that represents a ring structure model of the zeolite‐n‐heptane interacting system. The MM part was represented by means of an electrostatic forcefield (ESFF), which assesses the electronic embedding. The chosen QM/MM silicalite base model contains 3862 atoms. The studied SiO₂/Al₂O₃ ratios were 2300, 573.5, 287.7, and 189.83, containing 1, 4, 8, and 12 Al atoms, respectively. For the first ratio, the site for the substitution of Al for Si was that of minimum QM total energy value, because this replacement was done in the QM region. For the other SiO₂/Al₂O₃ ratios, the Al atoms were randomly spread through the MM region in accordance with the Lowenstein substitution rule. These results show the importance of the environment on the electronic properties in the QM region, where the active site lies, and their effects on the earlier steps on the activation experienced by the n‐heptane moiety. A minimal content of 12 Al atoms produces significant effects of the environment on the electronic structure of the QM region. Moreover, the carbocationic character of n‐heptane increases with the aluminum content. © 2002 Wiley Periodicals, Inc. Int J Quantum Chem, 2002     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

应用QM与ABEEM/MM结合的方法研究NAMI-A的结构性质

应用量子力学(QM)与ABEEM浮动电荷力场(ABEEM/MM)相结合的方法研究了抗癌药物NAMI-A在水溶液中的结构性质. 所有的结构优化都是在DFT的B3LYP方法下采用6-31G(d,p)和LanL2DZ基组完成的, 没有加入任何限制性条件. 结果表明, 优化得到的NAMI-A构型受不同环境及方法的影响均有变化. 与气相中得到的构型相比, QM/MM迭代优化得到构型要比PCM的构型变化更明显. QM/MM (ABEEM/MM)迭代优化得到的NAMI-A构型比QM/MM (OPLS-AA)的变化要小. 总之, 溶剂通过极化效应对NAMI-A结构、电荷分布及径向分布函数等性质均有影响, 客观地处理极化效应才能正确地反映QM区的性质.     

On the effect of a variation of the force field, spatial boundary condition and size of the QM region in QM/MM MD simulations

During the past years, the use of combined quantum-classical, QM/MM, methods for the study of complex biomolecular processes, such as enzymatic reactions and photocycles, has increased considerably. The quality of the results obtained from QM/MM calculations is largely dependent on five aspects to be considered when setting up a molecular model: the QM Hamiltonian, the MM Hamiltonian or force field, the boundary and coupling between the QM and MM regions, the size of the QM region and the boundary condition for the MM region. In this study, we systematically investigate the influence of a variation of the molecular mechanics force field and the size of the QM region in QM/MM MD simulations on properties of the photoactive part of the blue light photoreceptor protein AppA. For comparison, we additionally performed classical MD simulations and studied the effect of a variation of the type of spatial boundary condition. The classical boundary conditions and the force field used in a QM/MM MD simulation are shown to have non-neglegible effects upon the structural and energetic properties of the protein which makes it advisable to minimize computational artifacts in QM/MM MD simulations by application of periodic boundary conditions and a thermodynamically calibrated force field. A comparison of the structural and energetic properties of MD simulations starting from two alternative, different X-ray structures for the blue light utilizing flavin protein in its dark state indicates a slight preference of the two force fields used for the so-called Anderson structure over the Jung structure.     

Molecular insight into the interaction mechanisms of amino‐2H‐imidazole derivatives with BACE1 protease: A QM/MM and QTAIM study

In this study, we described quantitatively the interactions between two new amino‐2H‐imidazole inhibitors ((R)‐1t and (S)‐1m) and BACE1 using a hybrid quantum mechanics‐molecular mechanical (QM/MM) method together with a quantum theory of atoms In Molecules (QTAIM) analysis. Our computational calculations revealed that the binding affinity of these compounds is mostly related to the amino‐2H‐imidazole core, which interact tightly with the aspartate dyad of the active site. The interactions were stronger when the inhibitors presented a bulky substituent with a hydrogen bond acceptor motif pointing toward Trp76, such as the 3,5‐dimethyl‐4‐methoxyphenyl group of compound (S)‐1m. Furthermore, the QTAIM analysis revealed that many hydrophobic interactions complement cooperatively the hydrogen bond which is not present when compound (R)‐1t is bound to the enzyme. The combined QM/MM‐QTAIM analysis allows identifying the interactions that account for the activity difference between compounds, even at a nanomolar range.     

The implementation of a fast and accurate QM/MM potential method in Amber

Version 9 of the Amber simulation programs includes a new semi-empirical hybrid QM/MM functionality. This includes support for implicit solvent (generalized Born) and for periodic explicit solvent simulations using a newly developed QM/MM implementation of the particle mesh Ewald (PME) method. The code provides sufficiently accurate gradients to run constant energy QM/MM MD simulations for many nanoseconds. The link atom approach used for treating the QM/MM boundary shows improved performance, and the user interface has been rewritten to bring the format into line with classical MD simulations. Support is provided for the PM3, PDDG/PM3, PM3CARB1, AM1, MNDO, and PDDG/MNDO semi-empirical Hamiltonians as well as the self-consistent charge density functional tight binding (SCC-DFTB) method. Performance has been improved to the point where using QM/MM, for a QM system of 71 atoms within an explicitly solvated protein using periodic boundaries and PME requires less than twice the cpu time of the corresponding classical simulation.     

How Chorismatases Regulate Distinct Reaction Channels in a Single Conserved Active Pocket: Mechanistic Analysis with QM/MM (ONIOM) Investigations

The FkbO and Hyg5 subfamilies of chorismatases share the same active-site architectures, but perform distinct reaction mechanisms, that is, FkbO employs a hydrolysis reaction whereas Hyg5 proceeds through an intramolecular mechanism. Despite extensive research efforts, the detailed mechanism of the product selectivity in chorismatases need to be further unmasked. In this study, the effects of the A/G residue group (A244^FkbO/G240^Hyg5) and the V/Q residue group (V209^FkbO/Q201^Hyg5) on the catalytic mechanisms are investigated by employing molecular dynamics simulations and hybrid quantum mechanical/molecular mechanical (QM/MM) calculations of the two wild-type models (FkbO/CHO and Hyg5/CHO; CHO=chorismate) and four mutants models (A244G-FkbO/CHO and G240A-Hyg5/CHO; V209Q-FkbO/CHO and Q201V-Hyg5/CHO). Our results showed that the A/G residue group mentioned by previous works would cause changes in the binding states of the substrate and the orientation of the catalytic glutamate, but only these changes affect the product selectivity in chorismatases limitedly. Interestingly, the distal V/Q residue group, which determines the internal water self-regulating ability at the active site, has significant impact on the selectivity of the catalytic mechanisms. The V/Q residue group is suggested to be an important factor to control the catalytic activities in chorismatases. The results are consistent with biochemical and structural experiments, providing novel insight into the mechanism of product selectivity in chorismatases.     

Contributions of pauli repulsions to the energetics and physical properties computed in QM/MM methods

Conventional combined quantum mechanical/molecular mechanical (QM/MM) methods lack explicit treatment of Pauli repulsions between the quantum‐mechanical and molecular‐mechanical subsystems. Instead, classical Lennard‐Jones (LJ) potentials between QM and MM nuclei are used to model electronic Pauli repulsion and long‐range London dispersion, despite the fact that the latter two are inherently of quantum nature. Use of the simple LJ potential in QM/MM methods can reproduce minimal geometries and energies of many molecular clusters reasonably well, as compared to full QM calculations. However, we show here that the LJ potential cannot correctly describe subtle details of the electron density of the QM subsystem because of the neglect of Pauli repulsions between the QM and MM subsystems. The inaccurate electron density subsequently affects the calculation of electronic and magnetic properties of the QM subsystem. To explicitly consider Pauli interactions with QM/MM methods, we propose a method to use empirical effective potentials on the MM atoms. The test case of the binding energy and magnetic properties of a water dimer shows promising results for the general application of effective potentials to mimic Pauli repulsions in QM/MM calculations. © 2013 Wiley Periodicals, Inc.     

Fast generation of broken‐symmetry states in a large system including multiple iron–sulfur assemblies: Investigation of QM/MM energies,clusters charges,and spin populations

A density functional theory study is presented regarding the energetics and the Mulliken population analyses of a quantum mechanical/molecular mechanical (QM/MM) system including multiple iron–sulfur clusters in the QM region. The [FeFe]‐hydrogenase from Desulfovibrio desulfuricans was studied, and both the active site (an Fe₆S₆ assembly generally referred to as the H‐cluster) and an ancillary Fe₄S₄ site were treated at the BP86‐RI/TZVP level. The antiferromagnetic coupling that characterizes both sites was modeled using the broken‐symmetry (BS) approach. For such a QM system, 36 different BS couplings can be defined, depending on the localization of spin excess on the various spin centers. All the BS states were obtained by means of an effective and simple method for spin localization, that is here described and compared with more sophisticated approaches already available in literature. The variation of the QM/MM energy among the various geometry‐optimized protein models was found to be less than 25 kJ mol^–1. This energy variation almost doubles if no geometry optimization is performed. A detailed analysis of the additive nature of these variations in QM/MM energy is reported. The Mulliken charges show very small variations among the 36 BS states, whereas the Mulliken spin populations were found to be somewhat more variable. The relevance of such variations is discussed in light of the available Mössbauer and Electron Paramagnetic Resonance (EPR) spectroscopic data for the enzyme. Finally, the influence of the basis set on the spin populations, charges, and structural parameters of the models was investigated, by means of QM/MM computations on the same system at the BP86‐RI/SVP level. © 2010 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

QM/MM study of catalytic methyl transfer by the N5-glutamine SAM-dependent methyltransferase and its inhibition by the nitrogen analogue of coenzyme

The combined density functional quantum mechanical/molecular mechanical (QM/MM) approach has been used to investigate methyl-transfer reactions catalyzed by the N(5)-glutamine S-adenosyl-L-methionine (SAM)-dependent methyltransferase (HemK) and the coenzyme-modified HemK with the replacement of SAM by a nitrogen analogue. Calculations reveal that the catalytic methyl transfer by HemK is an energy-favored process with an activation barrier of 15.7 kcal/mol and an exothermicity of 12.0 kcal/mol, while the coenzyme-modified HemK is unable to catalyze the methyl transfer because of a substantial barrier of 20.6 kcal/mol and instability of the product intermediate. The results lend support to the experimental proposal that the nitrogen analogue of the SAM coenzyme should be a practicable inhibitor for the catalytic methyl transfer by HemK. Comparative QM/MM calculations show that the protein environment, especially the residues Asn197 and Pro198 in the active site, plays a pivotal role in stabilizing the transition state and regulating the positioning of reactive groups.