Spectral tuning in visual pigments: an ONIOM(QM:MM) study on bovine rhodopsin and its mutants

We have investigated geometries and excitation energies of bovine rhodopsin and some of its mutants by hybrid quantum mechanical/molecular mechanical (QM/MM) calculations in ONIOM scheme, employing B3LYP and BLYP density functionals as well as DFTB method for the QM part and AMBER force field for the MM part. QM/MM geometries of the protonated Schiff-base 11- cis-retinal with B3LYP and DFTB are very similar to each other. TD-B3LYP/MM excitation energy calculations reproduce the experimental absorption maximum of 500 nm in the presence of native rhodopsin environment and predict spectral shifts due to mutations within 10 nm, whereas TD-BLYP/MM excitation energies have red-shift error of at least 50 nm. In the wild-type rhodopsin, Glu113 shifts the first excitation energy to blue and accounts for most of the shift found. Other amino acids individually contribute to the first excitation energy but their net effect is small. The electronic polarization effect is essential for reproducing experimental bond length alternation along the polyene chain in protonated Schiff-base retinal, which correlates with the computed first excitation energy. It also corrects the excitation energies and spectral shifts in mutants, more effectively for deprotonated Schiff-base retinal than for the protonated form. The protonation state and conformation of mutated residues affect electronic spectrum significantly. The present QM/MM calculations estimate not only the experimental excitation energies but also the source of spectral shifts in mutants.     

Optimization of the Lennard-Jones parameters for a combined ab initio quantum mechanical and molecular mechanical potential using the 3-21G basis set

A combined ab initio quantum mechanical and molecular mechanical (AI-QM/MM) potential for use in molecular modeling and simulation has been described. In this article, we summarize a procedure for deriving the empirical parameters embedded in a combined QM/MM model and suggest a set of Lennard-Jones parameters for the combined ab initio 3-21G and MM OPLS-TIP3P (AI-3/MM) potential. Interaction energies and geometrical parameters predicted with the AI-3/MM model for over 80 hydrogen-bonded complexes of organic compounds with water were found to be in good accord with ab initio 6-31G(d) results. We anticipate that the AI-3/MM potential should be reasonable for use in condensed phase simulations. © 1996 by John Wiley & Sons, Inc.     

Ab initio and QM/MM study of electron addition on the disulfide bond in thioredoxin

Thioredoxin controls the intracellular redox potential through a disulfide/dithiol couple. Under conditions of oxidative stress, this protein functions via one-electron exchange, in which formation of the disulfide radical anion occurs. Combined quantum mechanical (QM) and molecular mechanical (MM) calculations using two- and three-level ONIOM schemes were performed on the thioredoxin (Trx) protein of Chlamydomonas reinhardtii in its oxidized-disulfide and one-electron-reduced forms. In both cases, the active site disulfide moiety was described at the MP2(fc)/6-31+G(d) level, and larger regions of varying sizes around the active site were described at the B3LYP/6-31+G(d) level. The remainder of the 112 residues and 33 water molecules of the crystal structure (PDB entry 1EP7) were described by the AMBER force field. Adiabatic electron affinities were calculated for the disulfide bond in all systems. Separate QM or QM/QM calculations were performed on the QM regions to establish the role of the remainder of the protein on the active site properties. The radical anion species becomes more stable as the number of amide groups in the vicinity increases. One-electron reduction potentials were calculated for the small molecule models, and approximated for the protein for which the values are similar to the experimental one (approximately 0 V). This high reduction potential is due to interaction with the charged end of Lys40, as indicated by mutation in silico to norleucine. The inclusion of the protonated Asp30 side chain and a water molecule in the QM region leads to an increase in the electron affinity. Proton transfer from the Asp30 side chain to the Cys39 sulfur in the radical anion species is strongly disfavored. The radical anion is more stable than the protonated form, which is consistent with experimental results.     

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.     

Toward accurate barriers for enzymatic reactions: QM/MM case study on p-hydroxybenzoate hydroxylase

The hydroxylation reaction catalyzed by p-hydroxybenzoate hydroxylase has been investigated by quantum mechanical/molecular mechanical (QM/MM) calculations at different levels of QM theory. The solvated enzyme was modeled (approximately 23,000 atoms in total, 49 QM atoms). The geometries of reactant and transition state were optimized for ten representative pathways using semiempirical (AM1) and density functional (B3LYP) methods as QM components. Single-point calculations at B3LYP/MM optimized geometries were performed with local correlation methods [LMP2, LCCSD(T0)] and augmented triple-zeta basis sets. A careful validation of the latter approach with regard to all computational parameters indicates convergence of the QM contribution to the computed barriers to within approximately 1 kcal mol(-1). Comparison with the available experimental data supports this assessment.     

Molecular dynamics and free energy perturbation study of hydride-ion transfer step in dihydrofolate reductase using combined quantum and molecular mechanical model

We used molecular dynamics simulation and free energy perturbation (FEP) methods to investigate the hydride-ion transfer step in the mechanism for the nicotinamide adenine dinucleotide phosphate (NADPH)-dependent reduction of a novel substrate by the enzyme dihydrofolate reductase (DHFR). The system is represented by a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1 semiempirical molecular orbital method for the reacting substrate and NADPH cofactor fragments, the AMBER force field for DHFR, and the TIP3P model for solvent water. The FEP calculations were performed for a number of choices for the QM system. The substrate, 8-methylpterin, was treated quantum mechanically in all the calculations, while the larger cofactor molecule was partitioned into various QM and MM regions with the addition of “link” atoms (F, CH₃, and H). Calculations were also carried out with the entire NADPH molecule treated by QM. The free energies of reaction and the net charges on the NADPH fragments were used to determine the most appropriate QM/MM model. The hydride-ion transfer was also carried out over several FEP pathways, and the QM and QM/MM component free energies thus calculated were found to be state functions (i.e., independent of pathway). A ca. 10 kcal/mol increase in free energy for the hydride-ion transfer with an activation barrier of ca. 30 kcal/mol was calculated. The increase in free energy on the hydride-ion transfer arose largely from the QM/MM component. Analysis of the QM/MM energy components suggests that, although a number of charged residues may contribute to the free energy change through long-range electrostatic interactions, the only interaction that can account for the 10 kcal/mol increase in free energy is the hydrogen bond between the carboxylate side chain of Glu30 (avian DHFR) and the activated (protonated) substrate. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 977–988, 1998     

Mechanism of proteolysis in matrix metalloproteinase‐2 revealed by QM/MM modeling

The mechanism of enzymatic peptide hydrolysis in matrix metalloproteinase‐2 (MMP‐2) was studied at atomic resolution through quantum mechanics/molecular mechanics (QM/MM) simulations. An all‐atom three‐dimensional molecular model was constructed on the basis of a crystal structure from the Protein Data Bank (ID: 1QIB), and the oligopeptide Ace‐Gln‐Gly～Ile‐Ala‐Gly‐Nme was considered as the substrate. Two QM/MM software packages and several computational protocols were employed to calculate QM/MM energy profiles for a four‐step mechanism involving an initial nucleophilic attack followed by hydrogen bond rearrangement, proton transfer, and C? N bond cleavage. These QM/MM calculations consistently yield rather low overall barriers for the chemical steps, in the range of 5–10 kcal/mol, for diverse QM treatments (PBE0, B3LYP, and BB1K density functionals as well as local coupled cluster treatments) and two MM force fields (CHARMM and AMBER). It, thus, seems likely that product release is the rate‐limiting step in MMP‐2 catalysis. This is supported by an exploration of various release channels through QM/MM reaction path calculations and steered molecular dynamics simulations. © 2015 Wiley Periodicals, Inc.     

An examination of a density functional/molecular mechanical coupled potential

In this article we describe the coupling of a density functional (DF) Hamiltonian with the molecular mechanics (MM) potential function AMBER. We examine a series of test cases in which we compare the binding energies and geometries of the complexes predicted by this coupled potential with those predicted by other theoretical methods and experiment to establish the relative accuracy of the DF/MM coupled potential. We find that the DF/MM coupled potential performs well in most cases studied and, in general, outperforms the semiempirical/MM approach. The interaction energies and structures obtained using this method appear to be insensitive to the use of nonlocal (NL) corrections to the DF method. The is fortuitous because the NL treatment is significantly more computationally expensive than the local treatment. However, NL corrections may be required to predict accurately the shape of potential energy surfaces that involve bond breaking and formation. The DF/MM method has also been applied to the determination of the solvation free energy for a series of ions using free-energy perturbation methods. The results obtained are good and can be improved by a simple scaling of the Lennard-Jones parameters for the ion in question. © 1995 by John Wiley & Sons, Inc.     

YinYang atom: a simple combined ab initio quantum mechanical molecular mechanical model

A simple interface is proposed for combined quantum mechanical (QM) molecular mechanical (MM) calculations for the systems where the QM and MM regions are connected through covalent bonds. Within this model, the atom that connects the two regions, called YinYang atom here, serves as an ordinary MM atom to other MM atoms and as a hydrogen-like atom to other QM atoms. Only one new empirical parameter is introduced to adjust the length of the connecting bond and is calibrated with the molecule propanol. This model is tested with the computation of equilibrium geometries and protonation energies for dozens of molecules. Special attention is paid on the influence of MM point charges on optimized geometry and protonation energy, and it is found that it is important to maintain local charge-neutrality in the MM region in order for the accurate calculation of the protonation and deprotonation energies. Overall the simple YinYang atom model yields comparable results to some other QM/MM models.     

Importance of selecting proper basis set in quantum mechanical studies of potential energy surfaces of carbohydrates

An extensive quantum mechanical study of a water dimer suggests that the introduction of a diffuse function into the basis set, which significantly reduces the basis set superposition error (BSSE) in the hydrogen bonding energy calculation, is the key to better calculations of the potential energy surfaces of carbohydrates. This article examines the potential energy surfaces of selected d -aldo- and d -ketohexoses (a total of 82 conformers) by quantum mechanics (QM) and molecular mechanics (MM) methods. In contrast to the results with a smaller basis set (B3LYP/6-31G** 5d), we found at the higher level calculation (B3LYP/6-311++G(2d,2p)//B3LYP/6-31G** 5d) that, in most cases, the furanose forms are less stable than the pyranose forms. These discrepancies are mainly due to the fact that intramolecular hydrogen bonding energies are overestimated in the lower level calculations. The higher level QM calculations of the potential energy surfaces of d -aldo- and d -ketohexoses now are more comparable to the MM3 results. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1593–1603, 1999     

A minimal implementation of the AMBER-GAUSSIAN interface for ab initio QM/MM-MD simulation

For applying to a number of theoretical methodologies based on an ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics method connecting AMBER9 with GAUSSIAN03, we have developed an AMBER-GAUSSIAN interface (AG-IF), which can be one of the simplest architectures. In the AG-IF, only a few subroutines addition is necessary to retrieve the QM/MM energy and forces, obtained by GAUSSIAN, for solving a set of Newtonian equations of motion in AMBER. It is, therefore, easy to be modified for individual applications since AG-IF utilizes most of those functions originally equipped not only in AMBER but also in GAUSSIAN. In the present minimal implementation, only AMBER is modified, whereas GAUSSIAN is left unchanged. Moreover, a different method of calculating electrostatic forces of MM atoms interacting with QM region is proposed. Using the AG-IF, we also demonstrate three examples of application: (i) the QM versus MM comparison in the radial distribution function, (ii) the free energy gradient method, and (iii) the charge from interaction energy and forces.     

Simulations of liquid ammonia based on the combined quantum mechanical/molecular mechanical (QM/MM) approach

Two combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations, namely, HF/MM and B3LYP/MM, have been performed to investigate the local structure and dynamics of liquid ammonia. The most interesting region, a sphere containing a central reference molecule and all its nearest surrounding molecules (first coordination shell), was treated by the Hartree-Fock (HF) and hybrid density functional B3LYP methods, whereas the rest of the system was described by the classical pair potentials. On the basis of both HF and B3LYP methods, it is observed that the hydrogen bonding in this peculiar liquid is weak. The structure and dynamics of this liquid are suggested to be determined by the steric packing effects, rather than by the directional hydrogen bonding interactions. Compared to previous empirical as well as Car-Parrinello (CP) molecular dynamics studies, our QM/MM simulations provide detailed information that is in better agreement with experimental data.     

Structural forms of green fluorescent protein by quantum mechanics/molecular mechanics calculations

The molecular modeling of structural forms of the green fluorescent protein (GFP) with the Ser65Thr single-site mutation was performed by the quantum mechanics/molecular mechanics (QM/MM) method. Two model systems were constructed based on the crystallographic structure from the Protein Data Bank (PDB entry code 1EMA.) The model systems differ in the initial protonation state of the side chain of the amino acid residue Glu222 near the chromophore. The atomic coordinates of the protein macromolecule corresponding to the equilibrium geometric configurations were determined by total energy minimization using the QM/MM method within the density functional theory approximation PBE0/cc-pVDZ for the quantum subsystem that consists of the chromophore, a water molecule, and the side chains of Arg96, Glu222, and Ser205, and with the parameters of the AMBER force field for the molecular mechanics subsystem. In the analysis of the results, particular attention was given to the hydrogen bond redistribution in the chromophore-containing region of the protein caused by a change in the protonation state of the chromophore. The results obtained from the model containing the initially protonated side chain of Glu222 suggest a new interpretation of the photophysical processes in the green fluorescent protein.     

Molecular mechanics (MM3) calculations on lithium amide compounds

The MM3 force field has been extended to deal with the lithium amide molecules that are widely used as efficient catalysts for stereoselective asymmetric synthesis. The MM3 force field parameters have been determined on the basis of the ab initio MP2/6-31G* and/or DFT (B3LYP/6-31G*, B3-PW91/6-31G*) geometry optimization calculations. To evaluate the electronic interactions specific to the lithium amides derived from the diamine molecules properly, the Lewis bonding potential term for the interaction between the lithium atom and the nonbonded adjacent electronegative atom such as nitrogen was introduced into the MM3 force field. The bond dipoles were evaluated correctly from the electronic charges on the atoms calculated by fitting to the electrostatic potential at points selected. The MM3 results on the molecular structures, conformational energies, and vibrational spectra show good agreement with those from the quantum mechanical calculations.     

Combined quantum mechanical/molecular mechanical study on the pentacoordinated ferric and ferrous cytochrome P450cam complexes

The pentacoordinated ferric and ferrous cytochrome P450(cam) complexes have been investigated by combined quantum mechanical/molecular mechanical (QM/MM) calculations in the presence of a protein/solvent environment and by QM calculations on the isolated QM regions with use of density functional theory. The B3LYP functional has been found more reliable than the BLYP and BHLYP functionals for estimating the relative state energies. The B3LYP/CHARMM calculations with an all-electron basis set for iron give high-spin ground states for the title complexes, in agreement with experiment. The comparison of the B3LYP/CHARMM results of the entire protein system with the B3LYP calculations on the naked QM regions shows that the amount of stabilization by the protein environment is largest for the intermediate-spin states, followed by the high-spin states of the complexes. The calculation of M?ssbauer parameters in the presence of the enzyme environment confirms the double occupation of the d(xz) orbital in the quintet spin state of the ferrous complex, consistent with the computed QM/MM energies in the enzyme environment, while the d(x)2(-)(y)2 orbital is doubly occupied in the gas-phase quintet state.     

Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations

Combined ab initio quantum mechanical and molecular mechanical calculations have been widely used for modeling chemical reactions in complex systems such as enzymes, with most applications being based on the determination of a minimum energy path connecting the reactant through the transition state to the product in the enzyme environment. However, statistical mechanics sampling and reaction dynamics calculations with a combined ab initio quantum mechanical (QM) and molecular mechanical (MM) potential are still not feasible because of the computational costs associated mainly with the ab initio quantum mechanical calculations for the QM subsystem. To address this issue, a reaction path potential energy surface is developed here for statistical mechanics and dynamics simulation of chemical reactions in enzymes and other complex systems. The reaction path potential follows the ideas from the reaction path Hamiltonian of Miller, Handy and Adams for gas phase chemical reactions but is designed specifically for large systems that are described with combined ab initio quantum mechanical and molecular mechanical methods. The reaction path potential is an analytical energy expression of the combined quantum mechanical and molecular mechanical potential energy along the minimum energy path. An expansion around the minimum energy path is made in both the nuclear and the electronic degrees of freedom for the QM subsystem internal energy, while the energy of the subsystem described with MM remains unchanged from that in the combined quantum mechanical and molecular mechanical expression and the electrostatic interaction between the QM and MM subsystems is described as the interaction of the MM charges with the QM charges. The QM charges are polarizable in response to the changes in both the MM and the QM degrees of freedom through a new response kernel developed in the present work. The input data for constructing the reaction path potential are energies, vibrational frequencies, and electron density response properties of the QM subsystem along the minimum energy path, all of which can be obtained from the combined quantum mechanical and molecular mechanical calculations. Once constructed, it costs much less for its evaluation. Thus, the reaction path potential provides a potential energy surface for rigorous statistical mechanics and reaction dynamics calculations of complex systems. As an example, the method is applied to the statistical mechanical calculations for the potential of mean force of the chemical reaction in triosephosphate isomerase.     

A combined quantum chemical/molecular mechanical study of hydrogen-bonded systems

A computational approach, which involves the combination of the OPLS force field and molecular orbital MNDO , AM 1, and PM 3 methods, has been developed to describe the effects of a large, molecular mechanically simulated environment on the Hamiltonian of a quantum chemical system. To test the validity of the combined quantum mechanical/molecular mechanical (QM /MM ) potential, a systematic study of the structures and energies of neutral and charged hydrogen-bonded complexes has been carried out, including comparisons with pure semiempirical calculations and available experimental and ab initio data. It is shown that, in many cases, the hybrid QM /MM potential behaves better than do related MNDO /M , AM 1, and PM 3 methods. As a case in point, the draw-back of AM 1 favoring bifurcated H-bonded structures over single ones is not presented in the combined AM 1/OPLS scheme. Possible ways of improvement of the combined QM /MM potential are discussed. © 1992 John Wiley & Sons, Inc.     

An assessment of pure,hybrid, meta,and hybrid‐meta GGA density functional theory methods for open‐shell systems: The case of the nonheme iron enzyme 8R–LOX

The performance of a range density functional theory functionals combined in a quantum mechanical (QM)/molecular mechanical (MM) approach was investigated in their ability to reliably provide geometries, electronic distributions, and relative energies of a multicentered open‐shell mechanistic intermediate in the mechanism 8R–Lipoxygenase. With the use of large QM/MM active site chemical models, the smallest average differences in geometries between the catalytically relevant quartet and sextet complexes were obtained with the B3LYP^* functional. Moreover, in the case of the relative energies between ⁴II and ⁶II , the use of the B3LYP^* functional provided a difference of 0.0 kcal mol^–1. However, B3LYP^± and B3LYP also predicted differences in energies of less than 1 kcal mol^–1. In the case of describing the electronic distribution (i.e., spin density), the B3LYP^*, B3LYP, or M06‐L functionals appeared to be the most suitable. Overall, the results obtained suggest that for systems with multiple centers having unpaired electrons, the B3LYP^* appears most well rounded to provide reliable geometries, electronic structures, and relative energies. © 2012 Wiley Periodicals, Inc.     

Flexible effective fragment QM/MM method: validation through the challenging tests

A new version of the QM/MM method, which is based on the effective fragment potential (EFP) methodology [Gordon, M. et al., J Phys Chem A 2001, 105, 293] but allows flexible fragments, is verified through calculations of model molecular systems suggested by different authors as challenging tests for QM/MM approaches. For each example, the results of QM/MM calculations for a partitioned system are compared to the results of an all-electron ab initio quantum chemical study of the entire system. In each case we were able to achieve approximately similar or better accuracy of the QM/MM results compared to those described in original publications. In all calculations we kept the same set of parameters of our QM/MM scheme. A new test example is considered when calculating the potential of internal rotation in the histidine dipeptide around the C(alpha)bond;C(beta) side chain bond.