Combined QM/MM and path integral simulations of kinetic isotope effects in the proton transfer reaction between nitroethane and acetate ion in water

An integrated Feynman path integral-free energy perturbation and umbrella sampling (PI-FEP/UM) method has been used to investigate the kinetic isotope effects (KIEs) in the proton transfer reaction between nitroethane and acetate ion in water. In the present study, both nuclear and electronic quantum effects are explicitly treated for the reacting system. The nuclear quantum effects are represented by bisection sampling centroid path integral simulations, while the potential energy surface is described by a combined quantum mechanical and molecular mechanical (QM/MM) potential. The accuracy essential for computing KIEs is achieved by a FEP technique that transforms the mass of a light isotope into a heavy one, which is equivalent to the perturbation of the coordinates for the path integral quasiparticle in the bisection sampling scheme. The PI-FEP/UM method is applied to the proton abstraction of nitroethane by acetate ion in water through molecular dynamics simulations. The rule of the geometric mean and the Swain-Schaad exponents for various isotopic substitutions at the primary and secondary sites have been examined. The computed total deuterium KIEs are in accord with experiments. It is found that the mixed isotopic Swain-Schaad exponents are very close to the semiclassical limits, suggesting that tunneling effects do not significantly affect this property for the reaction between nitroethane and acetate ion in aqueous solution.     

QM/MM calculation of solvent effects on absorption spectra of guanine

Electronic spectra of guanine in the gas phase and in water were studied by quantum mechanical/molecular mechanical (QM/MM) methods. Geometries for the excited‐state calculations were extracted from ground‐state molecular dynamics (MD) simulations using the self‐consistent‐charge density functional tight binding (SCC‐DFTB) method for the QM region and the TIP3P force field for the water environment. Theoretical absorption spectra were generated from excitation energies and oscillator strengths calculated for 50 to 500 MD snapshots of guanine in the gas phase (QM) and in solution (QM/MM). The excited‐state calculations used time‐dependent density functional theory (TDDFT) and the DFT‐based multireference configuration interaction (DFT/MRCI) method of Grimme and Waletzke, in combination with two basis sets. Our investigation covered keto‐N7H and keto‐N9H guanine, with particular focus on solvent effects in the low‐energy spectrum of the keto‐N9H tautomer. When compared with the vertical excitation energies of gas‐phase guanine at the optimized DFT (B3LYP/TZVP) geometry, the maxima in the computed solution spectra are shifted by several tenths of an eV. Three effects contribute: the use of SCC‐DFTB‐based rather than B3LYP‐based geometries in the MD snapshots (red shift of ca. 0.1 eV), explicit inclusion of nuclear motion through the MD snapshots (red shift of ca. 0.1 eV), and intrinsic solvent effects (differences in the absorption maxima in the computed gas‐phase and solution spectra, typically ca. 0.1–0.3 eV). A detailed analysis of the results indicates that the intrinsic solvent effects arise both from solvent‐induced structural changes and from electrostatic solute–solvent interactions, the latter being dominant. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

LICHEM: A QM/MM program for simulations with multipolar and polarizable force fields

We introduce an initial implementation of the LICHEM software package. LICHEM can interface with Gaussian, PSI4, NWChem, TINKER, and TINKER–HP to enable QM/MM calculations using multipolar/polarizable force fields. LICHEM extracts forces and energies from unmodified QM and MM software packages to perform geometry optimizations, single‐point energy calculations, or Monte Carlo simulations. When the QM and MM regions are connected by covalent bonds, the pseudo‐bond approach is employed to smoothly transition between the QM region and the polarizable force field. A series of water clusters and small peptides have been employed to test our initial implementation. The results obtained from these test systems show the capabilities of the new software and highlight the importance of including explicit polarization. © 2016 Wiley Periodicals, Inc.     

Infrared spectroscopy of N-methylacetamide in water from high-level QM/MM calculations简

The infrared(IR) spectra of the N-methylacetamide molecule in water are calculated by using the MD simulation with high-level QM/MM corrections. The B3LYP and MP2 levels with 6-311++G** basis set are used for the QM region, respectively. Our results show all IR spectra at the B3LYP level are well consistent with the corresponding MP2 results. A dynamical charge fluctuation is observed for each atom along the simulation trajectories due to the electrostatic polarization(EP) effects from surrounding solvent environment. We find that the QM/MM corrected IR spectra satisfactorily reproduce the experimental vibrational features of amide I–III modes.     

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.     

Optimizing efficiency of perturbative Monte Carlo method

We introduce error weighting functions into the perturbative Monte Carlo method for use with a hybrid ab initio quantum mechanics/molecular mechanics (QM/MM) potential. The perturbative Monte Carlo approach introduced earlier provides a means to reduce the number of full SCF calculations in simulations using a QM/MM potential by evoking perturbation theory to calculate energy changes due to displacements of an MM molecule. The use of weighting functions, introduced here, allows an optimal number of MM molecule displacements to occur between the performance of the full self-consistent field calculations. This will allow the ab initio QM/MM approach to be applied to systems that require more accurate treatment of the QM and/or MM regions. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1632–1638, 1998     

应用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区的性质.     

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.     

Combined QM(DFT)/MM molecular dynamics simulations of the deamination of cytosine by yeast cytosine deaminase (yCD)

Extensive combined quantum mechanical (B3LYP/6‐31G*) and molecular mechanical (QM/MM) molecular dynamics simulations have been performed to elucidate the hydrolytic deamination mechanism of cytosine to uracil catalyzed by the yeast cytosine deaminase (yCD). Though cytosine has no direct binding to the zinc center, it reacts with the water molecule coordinated to zinc, and the adjacent conserved Glu64 serves as a general acid/base to shuttle protons from water to cytosine. The overall reaction consists of several proton‐transfer processes and nucleophilic attacks. A tetrahedral intermediate adduct of cytosine and water binding to zinc is identified and similar to the crystal structure of yCD with the inhibitor 2‐pyrimidinone. The rate‐determining step with the barrier of 18.0 kcal/mol in the whole catalytic cycle occurs in the process of uracil departure where the proton transfer from water to Glu64 and nucleophilic attack of the resulting hydroxide anion to C2 of the uracil ring occurs synchronously. © 2016 Wiley Periodicals, Inc.     

Geometry optimization with QM/MM,ONIOM, and other combined methods. I. Microiterations and constraints

Hybrid energy methods such as QM/MM and ONIOM, that combine different levels of theory into one calculation, have been very successful in describing large systems. Geometry optimization methods can take advantage of the partitioning of these calculations into a region treated at a quantum mechanical (QM) level of theory and the larger, remaining region treated by an inexpensive method such as molecular mechanics (MM). A series of microiterations can be employed to fully optimize the MM region for each optimization step in the QM region. Cartesian coordinates are used for the MM region and are chosen so that the internal coordinates of the QM region remain constant during the microiterations. The coordinates of the MM region are augmented to permit rigid body translation and rotation of the QM region. This is essential if any atoms in the MM region are constrained, but it also improves the efficiency of unconstrained optimizations. Because of the microiterations, special care is needed for the optimization step in the QM region so that the system remains in the same local valley during the course of the optimization. The optimization methodology with microiterations, constraints, and step-size control are illustrated by calculations on bacteriorhodopsin and other systems.     

QM/MM study of the taxadiene synthase mechanism

Combined quantum mechanics/molecular mechanics (QM/MM) calculations were used to investigate the reaction mechanism of taxadiene synthase (TXS). TXS catalyzes the cyclization of geranylgeranyl diphosphate (GGPP) to taxadiene (T) and four minor cyclic products. All these products originate from the deprotonation of carbocation intermediates. The reaction profiles for the conversion of GGPP to T as well as to minor products were calculated for different configurations of relevant TXS carbocation complexes. The QM region was treated at the M06-2X/TZVP level, while the CHARMM27 force field was used to describe the MM region. The QM/MM calculations suggest a reaction pathway for the conversion of GGPP to T, which slightly differs from previous proposals regarding the number of reaction steps and the conformation of the carbocations. The QM/MM results also indicate that the formation of minor products via water-assisted deprotonation of the carbocations is highly exothermic, by about −7 to −23 kcal/mol. Curiously, however, the computed barriers and reaction energies indicate that the formation of some of the minor products is more facile than the formation of T. Thus, the present QM/MM calculations provide detailed insights into possible reaction pathways and into the origin of the promiscuity of TXS, but they do not reproduce the product distribution observed experimentally. © 2019 Wiley Periodicals, Inc.     

Ab initio QM/MM study of excited state electron transfer between pyrene and 4,4′‐bis(dimethylamino)‐diphenylmethane with different solvent systems: Role of hydrogen bonding within solvent molecules

The exciplex is a charge transfer species formed in the process of electron transfer between an electron donor and an electron acceptor and hence is very sensitive to solvent polarity. In order to understand the role of solvent in exciplex formation between pyrene (PY) and 4,4′‐bis(dimethylamino)diphenylmethane (DMDPM), we used two types of solvent approximations: an implicit solvent model and an explicit solvent model. The difference in energies between the excited and the meta‐stable Frank–Condon state (ΔE) of the structures were assumed to correspond to the emission maximum of the exciplex in different solvents. The ΔE values show the trend of stabilization of the exciplex with an increase in solvent polarity. This trend in stabilization is substantially more prominent in the explicit solvent model than that with the implicit solvent model. The ΔE value obtained in methanol reflects equal stabilization compared to that in a more polar solvent, N,N‐dimethylformamide. This extra stabilization of the exciplex may be explained on the basis of the H‐bonding capability of the protic solvent, methanol. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2005     

Extension of QM/MM docking and its applications to metalloproteins

To overcome the limitation of conventional docking methods which assume fixed charge model from force field parameters, combined quantum mechanics/molecular mechanics (QM/MM) method has been applied to docking as a variable charge model and shown to exhibit improvement on the docking accuracy over fixed charge based methods. However, it has also been shown that there are a number of examples for which adoption of variable‐charge model fails to reproduce the native binding modes. In particular, for metalloproteins, previously implemented method of QM/MM docking failed most often. This class of proteins has highly polarized binding sites at which high‐coordinate‐numbered metal ions reside. We extend the QM/MM docking method so that protein atoms surrounding the binding site along with metal ions are included as quantum region, as opposed to only ligand atoms. This extension facilitates the required scaling of partial charges on metal ions leading to prediction of correct binding modes in metalloproteins. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Solvent effects on the electronic absorption spectrum of formamide studied by a sequential Monte Carlo/quantum mechanical approach

Sequential Monte Carlo/quantum mechanical calculations are performed to study the solvent effects on the electronic absorption spectrum of formamide (FMA) in aqueous solution, varying from hydrogen bonds to the outer solvation shells. Full quantum-mechanical intermediate neglect of differential overlap/singly excited configuration interaction calculations are performed in the supermolecular structures generated by the Monte Carlo simulation. The largest calculation involves the ensemble average of 75 statistically uncorrelated quantum mechanical results obtained with the FMA solute surrounded by 150 water solvent molecules. We find that the n → π* transition suffers a blueshift of 1,600 cm⁻¹ upon solvation and the π → π* transition undergoes a redshift of 800 cm⁻¹. On average, 1.5 hydrogen bonds are formed between FMA and water and these contribute with about 20% and about 30% of the total solvation shifts of the n → π* and π → π* transitions, respectively. The autocorrelation function of the energy is used to sample configurations from the Monte Carlo simulation, and the solvation shifts are shown to be converged values. Received: 14 March 2002 / Accepted: 3 April 2002 / Published online: 24 June 2002     

量子力学和分子力学组合方法及其应用

QM/MM组合方法在研究凝聚态中的化学反应及生物大分子的结构和活性之间的关系等方面已取得重要进展。这一方法的要点在于将大体系配分成几部分,根据需要对不同部分进行不同级别的处理,因此既利用了量子力学的精确性,又利用了分子力学的高效性。对QM/MM组合理论及其一些最新进展作一简单介绍,并以最近进行了几个工作为例说明QM、MM组合方法的应用。     

Hybrid QM/MM studies on energetics of malonaldehyde in condensed phase

Energetics of eight enol isomers of a malonaldehyde (MA) at the ground state in CCl₄ solvent environment have been investigated using a hybrid quantum mechanical (QM)/molecular mechanics (MM) method. It is found that relative energies of the isomers slightly change due to interactions between MA and surrounding atoms. In an isolated environment all eight isomers have stable planar structure. On the other hand, most of the isomers have nonplanar structure in CCl₄ whose interactions with the solute molecules are, however, weak. Mainly, structural changes are found in the H? O? C?C dihedral angle ?_{H? O? C?C}, i.e., a hydrogen atom that is connected directly with an oxygen atom, is located at a nonplanar position, and other atoms remain almost planar. Vertical excitation energies of low‐lying excited states at the resultant optimized structure of each isomer are evaluated. The vertical excitation energies in CCl₄ are almost the same as those in the isolated environment, but some changes were found in the triplet excitation states. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2004     

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.     

Acceleration of a QM/MM‐QMC simulation using GPU

We accelerated an ab initio molecular QMC calculation by using GPGPU. Only the bottle‐neck part of the calculation is replaced by CUDA subroutine and performed on GPU. The performance on a (single core CPU + GPU) is compared with that on a (single core CPU with double precision), getting 23.6 (11.0) times faster calculations in single (double) precision treatments on GPU. The energy deviation caused by the single precision treatment was found to be within the accuracy required in the calculation, ～10^?5 hartree. The accelerated computational nodes mounting GPU are combined to form a hybrid MPI cluster on which we confirmed the performance linearly scales to the number of nodes. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011