Efficient Calculation of QM/MM Frequencies with the Mobile Block Hessian

Efficient Method for Numerical Calculations of Molecular Vibrational Frequencies by Exploiting Sparseness of Hessian Matrix

Molecular vibrational frequency analysis plays an important role in theoretical and computational chemistry. However, in many cases, the analytical frequencies are unavailable, whereas frequency calculations using conventional numerical methods are very expensive. In this work, we propose an efficient method to numerically calculate the frequencies. Our main strategies are to exploit the sparseness of the Hessian matrix and to construct the N-fold two-variable potential energy surfaces to fit the parabola parameters, which are later used for the construction of Hessian matrices. A set of benchmark calculations is performed for typical molecules of different sizes and complexities using the proposed method. The obtained frequencies are compared to those calculated with the analytical methods and conventional numerical methods. It is shown that the results yielded with the new method are in very good agreement with corresponding accurate values (with a maximum error of ∼20 cm-1), while the required computation resource is largely reduced compared to that required by conventional numerical methods. For medium-sized molecules, the calculational scaling is lowered to O(N1.6) (this work) from that of O(N2) (conventional numerical methods). For even larger molecules, more computational savings can be achieved, and the scaling is estimated to be quasilinear with respect to the molecular size.

AutoParams: An Automated Web-Based Tool To Generate Force Field Parameters for Molecular Dynamics Simulations

Many research questions benefit from molecular dynamics simulations to observe the motions and conformations of molecules over time, which rely on force fields that describe sets of common molecules by category. With the increase of importance for large data sets used in machine learning and growing computational efficiency, the ability to rapidly create large numbers of force field inputs is of high importance. Unusual molecules, such as nucleotide analogues, functionalized carbohydrates, and modified amino acids, are difficult to describe consistently using standard force fields, requiring the development of custom parameters for each unique molecule. While these parameters may be created by individual users, the process can become time-consuming or may introduce errors that may not be immediately apparent. We present an open-source automated parameter generation service, AutoParams, which requires minimal input from the user and creates useful Amber force field parameter sets for most molecules, particularly those that combine molecular types (e.g., a carbohydrate functionalized with a benzene). We include hierarchical atom-typing logic that makes it straightforward to expand with additional force fields and settings, and options for creating monomers in polymers, such as functionalized amino acids. It can be straightforwardly linked to any charge generation program and currently has interfaces to Psi4, PsiRESP, and TeraChem. It is open source and is available via GitHub. It includes error checking and testing protocols to ensure the parameters will be sufficient for subsequent molecular dynamics simulations and streamlines the creation of force field databases.

The effect of hydrogen-bonding on flavin's infrared absorption spectrum

Cluster and continuum solvation computational models are employed to model the effect of hydrogen bonding interactions on the vibrational modes of lumiflavin. Calculated spectra were compared to experimental Fourier-transform infrared (FTIR) spectra in the diagnostic 1450-1800 cm^-1 range, where intense ν_C=C, ν_C=N, [Formula: see text] , and [Formula: see text] stretching modes of flavin's isoalloxazine ring are found. Local mode analysis is used to describe the strength of hydrogen-bonding in cluster models. The computations indicate that ν_C=C and ν_C=N mode frequencies are relatively insensitive to intermolecular interactions while the [Formula: see text] and [Formula: see text] modes are sensitive to direct (and also indirect for [Formula: see text] ) hydrogen-bonding interactions. Although flavin is neutral, basis sets without the diffuse functions provide incorrect relative frequencies and intensities. The 6-31+G* basis set is found to be adequate for this system, and there is limited benefit to considering larger basis sets. Calculated vibrational mode frequencies agree with experimentally determined frequencies in solution when cluster models with multiple water molecules are used. Accurate simulation of relative FTIR band intensities, on the other hand, requires a continuum (or possibly quantum mechanical/molecular mechanical) model that accounts for long-range electrostatic effects. Finally, an experimental peak at ca. 1624 cm^-1 that is typically assigned to the [Formula: see text] vibrational stretching mode has a complicated shape that suggests multiple underlying contributions. Our calculations show that this band has contributions from both the C₆-C₇ and C₂ = O stretching vibrations.

Anharmonic Vibrational Calculations Based on Group-Localized Coordinates: Applications to Internal Water Molecules in Bacteriorhodopsin

An efficient anharmonic vibrational method is developed exploiting the locality of molecular vibration. Vibrational coordinates localized to a group of atoms are employed to divide the potential energy surface (PES) of a system into intra- and inter-group contributions. Then, the vibrational Schrödinger equation is solved based on a PES, in which the inter-group coupling is truncated at the harmonic level while accounting for the intra-group anharmonicity. The method is applied to a pentagonal hydrogen bond network (HBN) composed of internal water molecules and charged residues in a membrane protein, bacteriorhodopsin. The PES is calculated by the quantum mechanics/molecular mechanics (QM/MM) calculation at the level of B3LYP-D3/aug-cc-pVDZ. The infrared (IR) spectrum is computed using a set of coordinates localized to each water molecule and amino acid residue by second-order vibrational quasi-degenerate perturbation theory (VQDPT2). Benchmark calculations show that the proposed method yields the N-D/O-D stretching frequencies with an error of 7 cm-1 at the cost reduced by more than five times. In contrast, the harmonic approximation results in a severe error of 150 cm-1. Furthermore, the size of QM regions is carefully assessed to find that the QM regions should include not only the pentagonal HBN itself but also its HB partners. VQDPT2 calculations starting from transient structures obtained by molecular dynamics simulations have shown that the structural sampling has a significant impact on the calculated IR spectrum. The incorporation of anharmonicity, sufficiently large QM regions, and structural samplings are of essential importance to reproduce the experimental IR spectrum. The computational spectrum paves the way for decoding the IR signal of strong HBNs and helps elucidate their functional roles in biomolecules.

Harmonic Infrared and Raman Spectra in Molecular Environments Using the Polarizable Embedding Model

We present a fully analytic approach to calculate infrared (IR) and Raman spectra of molecules embedded in complex molecular environments modeled using the fragment-based polarizable embedding (PE) model. We provide the theory for the calculation of analytic second-order geometric derivatives of molecular energies and first-order geometric derivatives of electric dipole moments and dipole-dipole polarizabilities within the PE model. The derivatives are implemented using a general open-ended response theory framework, thus allowing for an extension to higher-order derivatives. The embedding-potential parameters used to describe the environment in the PE model are derived through first-principles calculations, thus allowing a wide variety of systems to be modeled, including solvents, proteins, and other large and complex molecular environments. Here, we present proof-of-principle calculations of IR and Raman spectra of acetone in different solvents. This work is an important step toward calculating accurate vibrational spectra of molecules embedded in realistic environments.

UV-visible absorption spectrum of FAD and its reduced forms embedded in a cryptochrome protein

Cryptochromes are a class of flavoproteins proposed as candidates to explain magnetoreception of animals, plants and bacteria. The main hypothesis is that a biradical is formed upon blue-light absorption by flavin adenine dinucleotide (FAD). In a protein milieu, the oxidized form of FAD can be reduced, leading to four redox derivative forms: anionic and neutral semi-reduced radicals, and anionic and neutral fully reduced forms. All these forms have a characteristic electronic absorption spectrum, with a strong vibrational resolution. Here, we carried out a normal mode analysis at the electrostatic embedding QM/MM level of theory to compute the vibrationally resolved absorption spectra of the five redox forms of FAD embedded in a plant cryptochrome. We show that explicitly accounting for vibrational broadening contributions to electronic transitions is essential to reproduce the experimental spectra. In the case of the neutral radical form of FAD, the absorption spectrum is reproduced only if the presence of a tryptophan radical is considered.

Efficient Analytic Second Derivative of Electrostatic Embedding QM/MM Energy: Normal Mode Analysis of Plant Cryptochrome

Analytic second derivatives of electrostatic embedding (EE) quantum mechanics/molecular mechanics (QM/MM) energy are important for performing vibrational analysis and simulating vibrational spectra of quantum systems interacting with an environment represented as a classical electrostatic potential. The main bottleneck of EE-QM/MM second derivatives is the solution of coupled perturbed equations for each MM atom perturbation. Here, we exploit the Q-vector method [J. Chem. Phys., 2019, 151, 041102] to workaround this bottleneck. We derive the full analytic second derivative of the EE-QM/MM energy, which allows us to compute QM, MM, and QM-MM Hessian blocks in an efficient and easy to implement manner. To show the capabilities of our method, we compute the normal modes for the full Arabidopsis thaliana plant cryptochrome. We show that the flavin adenine dinucleotide vibrations (QM subsystem) strongly mix with protein modes. We compute approximate vibronic couplings for the lowest bright transition, from which we extract spectral densities and the homogeneous broadening of FAD absorption spectrum in protein using vibrationally resolved electronic spectrum simulations.

Projector-Based Embedding Eliminates Density Functional Dependence for QM/MM Calculations of Reactions in Enzymes and Solution

Combined quantum mechanics/molecular mechanics (QM/MM) methods are increasingly widely utilized in studies of reactions in enzymes and other large systems. Here, we apply a range of QM/MM methods to investigate the Claisen rearrangement of chorismate to prephenate, in solution, and in the enzyme chorismate mutase. Using projector-based embedding in a QM/MM framework, we apply treatments up to the CCSD(T) level. We test a range of density functional QM/MM methods and QM region sizes. The results show that the calculated reaction energetics are significantly more sensitive to the choice of density functional than they are to the size of the QM region in these systems. Projector-based embedding of a wave function method in DFT reduced the 13 kcal/mol spread in barrier heights calculated at the DFT/MM level to a spread of just 0.3 kcal/mol, essentially eliminating dependence on the functional. Projector-based embedding of correlated ab initio methods provides a practical method for achieving high accuracy for energy profiles derived from DFT and DFT/MM calculations for reactions in condensed phases.

Frequency Range Selection Method for Vibrational Spectra

Theoretical calculations of vibrational properties are widely used to explain and predict experimental spectra. However, with standard quantum chemical methods all molecular motions are considered, which is rather time-consuming for large molecules. Because typically only a specific spectral region is of experimental interest, we propose here an efficient method that allows calculation of only a selected frequency interval. After a computationally cheap low-level estimate of the molecular motions, the computational time is proportional to the number of normal modes needed to describe this frequency range. Results for a medium-sized molecule show a reduction in computational time of up to 1 order of magnitude with negligible loss in accuracy. We also show that still larger computational savings are possible by using an additional intensity-selection procedure.

Ab Initio QM/MM Modeling of the Rate-Limiting Proton Transfer Step in the Deamination of Tryptamine by Aromatic Amine Dehydrogenase

Aromatic amine dehydrogenase (AADH) and related enzymes are at the heart of debates on the roles of quantum tunneling and protein dynamics in catalysis. The reaction of tryptamine in AADH involves significant quantum tunneling in the rate-limiting proton transfer step, shown by large H/D primary kinetic isotope effects (KIEs), with unusual temperature dependence. We apply correlated ab initio combined quantum mechanics/molecular mechanics (QM/MM) methods, at levels up to local coupled cluster theory (LCCSD(T)/(aug)-cc-pVTZ), to calculate accurate potential energy surfaces for this reaction, which are necessary for quantitative analysis of tunneling contributions and reaction dynamics. Different levels of QM/MM treatment are tested. Multiple pathways are calculated with fully flexible transition state optimization by the climbing-image nudged elastic band method at the density functional QM/MM level. The average LCCSD(T) potential energy barriers to proton transfer are 16.7 and 14.0 kcal/mol for proton transfer to the two carboxylate atoms of the catalytic base, Asp128β. The results show that two similar, but distinct pathways are energetically accessible. These two pathways have different barriers, exothermicity and curvature, and should be considered in analyses of the temperature dependence of reaction and KIEs in AADH and other enzymes. These results provide a benchmark for this prototypical enzyme reaction and will be useful for developing empirical models, and analyzing experimental data, to distinguish between different conceptual models of enzyme catalysis.

Multiscale Quantum Mechanics/Molecular Mechanics Simulations with Neural Networks

Molecular dynamics simulation with multiscale quantum mechanics/molecular mechanics (QM/MM) methods is a very powerful tool for understanding the mechanism of chemical and biological processes in solution or enzymes. However, its computational cost can be too high for many biochemical systems because of the large number of ab initio QM calculations. Semiempirical QM/MM simulations have much higher efficiency. Its accuracy can be improved with a correction to reach the ab initio QM/MM level. The computational cost on the ab initio calculation for the correction determines the efficiency. In this paper we developed a neural network method for QM/MM calculation as an extension of the neural-network representation reported by Behler and Parrinello. With this approach, the potential energy of any configuration along the reaction path for a given QM/MM system can be predicted at the ab initio QM/MM level based on the semiempirical QM/MM simulations. We further applied this method to three reactions in water to calculate the free energy changes. The free-energy profile obtained from the semiempirical QM/MM simulation is corrected to the ab initio QM/MM level with the potential energies predicted with the constructed neural network. The results are in excellent accordance with the reference data that are obtained from the ab initio QM/MM molecular dynamics simulation or corrected with direct ab initio QM/MM potential energies. Compared with the correction using direct ab initio QM/MM potential energies, our method shows a speed-up of 1 or 2 orders of magnitude. It demonstrates that the neural network method combined with the semiempirical QM/MM calculation can be an efficient and reliable strategy for chemical reaction simulations.

A radical rebound mechanism for the methane oxidation reaction promoted by the dicopper center of a pMMO enzyme: a computational perspective

Hydrogen Atom Transfer (HAT) promoted by a triplet state of the bis-oxoCu₂(iii) core generates a new radical rebound mechanism for the hydroxylation of methane catalyzed by the binuclear copper site of a pMMO enzyme.

Multiple active zones in hybrid QM/MM molecular dynamics simulations for large biomolecular systems

A new QM/MM molecular dynamics approach that can deal with the dynamics of large real systems involving several simultaneous active zones is presented.

Efficiently computing pathway free energies: New approaches based on chain-of-replica and Non-Boltzmann Bennett reweighting schemes

BACKGROUND

Accurately modeling condensed phase processes is one of computation's most difficult challenges. Include the possibility that conformational dynamics may be coupled to chemical reactions, where multiscale (i.e., QM/MM) methods are needed, and this task becomes even more daunting.

METHODS

Free energy simulations (i.e., molecular dynamics), multiscale modeling, and reweighting schemes.

RESULTS

Herein, we present two new approaches for mitigating the aforementioned challenges. The first is a new chain-of-replica method (off-path simulations, OPS) for computing potentials of mean force (PMFs) along an easily defined reaction coordinate. This development is coupled with a new distributed, highly-parallel replica framework (REPDstr) within the CHARMM package. Validation of these new schemes is carried out on two processes that undergo conformational changes. First is the simple torsional rotation of butane, while a much more challenging glycosidic rotation (in vacuo and solvated) is the second. Additionally, a new approach that greatly improves (i.e., possibly an order of magnitude) the efficiency of computing QM/MM PMFs is introduced and compared to standard schemes. Our efforts are grounded in the recently developed method for efficiently computing QM-based free energies (i.e., QM-Non-Boltzmann Bennett, QM-NBB). Again, we validate this new technique by computing the QM/MM PMF of butane's torsional rotation.

CONCLUSIONS

The OPS-REPDstr method is a promising new approach that overcomes many limitations of standard pathway simulations in CHARMM. The combination of QM-NBB with pathway techniques is very promising as it offers significant advantages over current procedures.

GENERAL SIGNIFICANCE

Efficiently computing potentials of mean force is a major, unresolved, area of interest. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Correcting for the free energy costs of bond or angle constraints in molecular dynamics simulations

BACKGROUND

Free energy simulations are an important tool in the arsenal of computational biophysics, allowing the calculation of thermodynamic properties of binding or enzymatic reactions. This paper introduces methods to increase the accuracy and precision of free energy calculations by calculating the free energy costs of constraints during post-processing. The primary purpose of employing constraints for these free energy methods is to increase the phase space overlap between ensembles, which is required for accuracy and convergence.

METHODS

The free energy costs of applying or removing constraints are calculated as additional explicit steps in the free energy cycle. The new techniques focus on hard degrees of freedom and use both gradients and Hessian estimation. Enthalpy, vibrational entropy, and Jacobian free energy terms are considered.

RESULTS

We demonstrate the utility of this method with simple classical systems involving harmonic and anharmonic oscillators, four-atomic benchmark systems, an alchemical mutation of ethane to methanol, and free energy simulations between alanine and serine. The errors for the analytical test cases are all below 0.0007kcal/mol, and the accuracy of the free energy results of ethane to methanol is improved from 0.15 to 0.04kcal/mol. For the alanine to serine case, the phase space overlaps of the unconstrained simulations range between 0.15 and 0.9%. The introduction of constraints increases the overlap up to 2.05%. On average, the overlap increases by 94% relative to the unconstrained value and precision is doubled.

CONCLUSIONS

The approach reduces errors arising from constraints by about an order of magnitude. Free energy simulations benefit from the use of constraints through enhanced convergence and higher precision.

GENERAL SIGNIFICANCE

The primary utility of this approach is to calculate free energies for systems with disparate energy surfaces and bonded terms, especially in multi-scale molecular mechanics/quantum mechanics simulations. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Molecular properties of excited electronic state: formalism, implementation, and applications of analytical second energy derivatives within the framework of the time-dependent density functional theory/molecular mechanics

This work extends our previous works [J. Liu and W. Z. Liang, J. Chem. Phys. 135, 014113 (2011); J. Liu and W. Z. Liang, J. Chem. Phys. 135, 184111 (2011)] on analytical excited-state energy Hessian within the framework of time-dependent density functional theory (TDDFT) to couple with molecular mechanics (MM). The formalism, implementation, and applications of analytical first and second energy derivatives of TDDFT/MM excited state with respect to the nuclear and electric perturbations are presented. Their performances are demonstrated by the calculations of adiabatic excitation energies, and excited-state geometries, harmonic vibrational frequencies, and infrared intensities for a number of benchmark systems. The consistent results with the full quantum mechanical method and other hybrid theoretical methods indicate the reliability of the current numerical implementation of developed algorithms. The computational accuracy and efficiency of the current analytical approach are also checked and the computational efficient strategies are suggested to speed up the calculations of complex systems with many MM degrees of freedom. Finally, we apply the current analytical approach in TDDFT/MM to a realistic system, a red fluorescent protein chromophore together with part of its nearby protein matrix. The calculated results indicate that the rearrangement of the hydrogen bond interactions between the chromophore and the protein matrix is responsible for the large Stokes shift.

Kinetic Isotope Effects from QM/MM Subset Hessians: "Cut-Off" Analysis for SN2 Methyl Transfer in Solution

Isotopic partition-function ratios and kinetic isotope effects for reaction of S-adenosylmethionine with catecholate in water are evaluated using a subset of 324 atoms within its surrounding aqueous environment at the AM1/TIP3P level. Two alternative methods for treating motion in the six librational degrees of freedom of the subset atoms relative to their environment are compared. A series of successively smaller subset Hessians are generated by cumulative deletion of rows and columns from the initial 972 × 972 Hessian. We find that it is better to treat these librations as vibrations than as translations and rotations and that there is no need to invoke the Teller-Redlich product rule. The validity of "cut-off" procedures for computation of isotope effects with truncated atomic subsets is assessed: to ensure errors in ln(KIE) < 1% (or 2% for the quantum-corrected KIE) for all isotopic substitutions considered, it is necessary to use a less-restrictive procedure than is suggested by the familiar two-bond cutoff rule.

MSCALE: A General Utility for Multiscale Modeling

The combination of theoretical models of macromolecules that exist at different spatial and temporal scales has become increasingly important for addressing complex biochemical problems. This work describes the extension of concurrent multiscale approaches, introduces a general framework for carrying out calculations, and describes its implementation into the CHARMM macromolecular modeling package. This functionality, termed MSCALE, generalizes both the additive and subtractive multiscale scheme (e.g. QM/MM ONIOM-type), and extends its support to classical force fields, coarse grained modeling (e.g. ENM, GNM, etc.), and a mixture of them all. The MSCALE scheme is completely parallelized with each subsystem running as an independent, but connected calculation. One of the most attractive features of MSCALE is the relative ease of implementation using the standard MPI communication protocol. This allows external access to the framework and facilitates the combination of functionality previously isolated in separate programs. This new facility is fully integrated with free energy perturbation methods, Hessian based methods, and the use of periodicity and symmetry, which allows the calculation of accurate pressures. We demonstrate the utility of this new technique with four examples; (1) subtractive QM/MM and QM/QM calculations; (2) multi-force field alchemical free energy perturbation; (3) integration with the SANDER module of AMBER and the TINKER package to gain access to potentials not available in CHARMM; and (4) mixed resolution (i.e. coarse grain / all-atom) normal mode analysis. The potential of this new tool is clearly established and in conclusion an interesting mathematical problem is highlighted and future improvements are proposed.