Efficient Approach to Reactive Molecular Dynamics with Accurate Forces

Connecting the quantum and classical mechanics simulation world: Applications of reactive step molecular dynamics simulations

This article presents the application of the reactive step molecular dynamics simulation method [M. Biedermann, D. Diddens, and A. Heuer, J. Chem. Theory Comput. 17, 1074 (2021)] toward two different atomistic, chemically reactive systems. During reactive steps, transitions from reactant to product molecules are modeled according to physically correct transition probabilities based on quantum chemical information about the reactions such as molecular reaction rates via instant exchange of the employed force field and a subsequent, short relaxation of the structure. In the first application, we study the follow-up reactions of singly reduced ethylene carbonate (EC) radicals in EC solution, first, via extensive ab initio molecular dynamics simulations and, second, with the reactive step algorithm. A direct comparison of both simulation methods shows excellent agreement. Then, we employ the reactive step algorithm to simulate the enolate formation of 2-methylcyclopropanone with the base lithium diisopropylamine. Thereby, we can demonstrate that the reactive step algorithm is also capable of capturing effects from kinetic vs thermodynamic control of chemical reactions during simulation.

Theoretical approaches for dynamical ordering of biomolecular systems

BACKGROUND

Living systems are characterized by the dynamic assembly and disassembly of biomolecules. The dynamical ordering mechanism of these biomolecules has been investigated both experimentally and theoretically. The main theoretical approaches include quantum mechanical (QM) calculation, all-atom (AA) modeling, and coarse-grained (CG) modeling. The selected approach depends on the size of the target system (which differs among electrons, atoms, molecules, and molecular assemblies). These hierarchal approaches can be combined with molecular dynamics (MD) simulation and/or integral equation theories for liquids, which cover all size hierarchies.

SCOPE OF REVIEW

We review the framework of quantum mechanical/molecular mechanical (QM/MM) calculations, AA MD simulations, CG modeling, and integral equation theories. Applications of these methods to the dynamical ordering of biomolecular systems are also exemplified.

MAJOR CONCLUSIONS

The QM/MM calculation enables the study of chemical reactions. The AA MD simulation, which omits the QM calculation, can follow longer time-scale phenomena. By reducing the number of degrees of freedom and the computational cost, CG modeling can follow much longer time-scale phenomena than AA modeling. Integral equation theories for liquids elucidate the liquid structure, for example, whether the liquid follows a radial distribution function.

GENERAL SIGNIFICANCE

These theoretical approaches can analyze the dynamic behaviors of biomolecular systems. They also provide useful tools for exploring the dynamic ordering systems of biomolecules, such as self-assembly. This article is part of a Special Issue entitled "Biophysical Exploration of Dynamical Ordering of Biomolecular Systems" edited by Dr. Koichi Kato.

Quantum mechanical fragment methods based on partitioning atoms or partitioning coordinates

Conspectus The development of more efficient and more accurate ways to represent reactive potential energy surfaces is a requirement for extending the simulation of large systems to more complex systems, longer-time dynamical processes, and more complete statistical mechanical sampling. One way to treat large systems is by direct dynamics fragment methods. Another way is by fitting system-specific analytic potential energy functions with methods adapted to large systems. Here we consider both approaches. First we consider three fragment methods that allow a given monomer to appear in more than one fragment. The first two approaches are the electrostatically embedded many-body (EE-MB) expansion and the electrostatically embedded many-body expansion of the correlation energy (EE-MB-CE), which we have shown to yield quite accurate results even when one restricts the calculations to include only electrostatically embedded dimers. The third fragment method is the electrostatically embedded molecular tailoring approach (EE-MTA), which is more flexible than EE-MB and EE-MB-CE. We show that electrostatic embedding greatly improves the accuracy of these approaches compared with the original unembedded approaches. Quantum mechanical fragment methods share with combined quantum mechanical/molecular mechanical (QM/MM) methods the need to treat a quantum mechanical fragment in the presence of the rest of the system, which is especially challenging for those parts of the rest of the system that are close to the boundary of the quantum mechanical fragment. This is a delicate matter even for fragments that are not covalently bonded to the rest of the system, but it becomes even more difficult when the boundary of the quantum mechanical fragment cuts a bond. We have developed a suite of methods for more realistically treating interactions across such boundaries. These methods include redistributing and balancing the external partial atomic charges and the use of tuned fluorine atoms for capping dangling bonds, and we have shown that they can greatly improve the accuracy. Finally we present a new approach that goes beyond QM/MM by combining the convenience of molecular mechanics with the accuracy of fitting a potential function to electronic structure calculations on a specific system. To make the latter practical for systems with a large number of degrees of freedom, we developed a method to interpolate between local internal-coordinate fits to the potential energy. A key issue for the application to large systems is that rather than assigning the atoms or monomers to fragments, we assign the internal coordinates to reaction, secondary, and tertiary sets. Thus, we make a partition in coordinate space rather than atom space. Fits to the local dependence of the potential energy on tertiary coordinates are arrayed along a preselected reaction coordinate at a sequence of geometries called anchor points; the potential energy function is called an anchor points reactive potential. Electrostatically embedded fragment methods and the anchor points reactive potential, because they are based on treating an entire system by quantum mechanical electronic structure methods but are affordable for large and complex systems, have the potential to open new areas for accurate simulations where combined QM/MM methods are inadequate.

Concerted hydrogen atom and electron transfer mechanism for catalysis by lysine-specific demethylase

We calculate the free energy profile for the postulated hydride transfer reaction mechanism for the catalysis of lysine demethylation by lysine-specific demethylase LSD1. The potential energy surface is obtained by using combined electrostatically embedded multiconfiguration molecular mechanics (EE-MCMM) and single-configuration molecular mechanics (MM). We employ a constant valence bond coupling term to obtain analytical energies and gradients of the EE-MCMM subsystem, which contains 45 quantum mechanics (QM) atoms and which is parametrized with density functional calculations employing specific reaction parameters obtained by matching high-level wave function calculations. In the MM region, we employ the Amber ff03 and TIP3P force fields. The free energy of activation at 300 K is calculated by molecular dynamics (MD) umbrella sampling on a system with 102,090 atoms as the maximum of the free energy profile along the reaction coordinate as obtained by the weighted histogram analysis method with 17 umbrella sampling windows. This yields a free energy of activation of only 10 kcal/mol, showing that the previously postulated direct hydride transfer reaction mechanism is plausible, although we find that it is better interpreted as a concerted transfer of a hydrogen atom and an electron.

Multiscale reactive molecular dynamics

Many processes important to chemistry, materials science, and biology cannot be described without considering electronic and nuclear-level dynamics and their coupling to slower, cooperative motions of the system. These inherently multiscale problems require computationally efficient and accurate methods to converge statistical properties. In this paper, a method is presented that uses data directly from condensed phase ab initio simulations to develop reactive molecular dynamics models that do not require predefined empirical functions. Instead, the interactions used in the reactive model are expressed as linear combinations of interpolating functions that are optimized by using a linear least-squares algorithm. One notable benefit of the procedure outlined here is the capability to minimize the number of parameters requiring nonlinear optimization. The method presented can be generally applied to multiscale problems and is demonstrated by generating reactive models for the hydrated excess proton and hydroxide ion based directly on condensed phase ab initio molecular dynamics simulations. The resulting models faithfully reproduce the water-ion structural properties and diffusion constants from the ab initio simulations. Additionally, the free energy profiles for proton transfer, which is sensitive to the structural diffusion of both ions in water, are reproduced. The high fidelity of these models to ab initio simulations will permit accurate modeling of general chemical reactions in condensed phase systems with computational efficiency orders of magnitudes greater than currently possible with ab initio simulation methods, thus facilitating a proper statistical sampling of the coupling to slow, large-scale motions of the system.

Interpolated mechanics-molecular mechanics study of internal rotation dynamics of the chromophore unit in blue fluorescent protein and its variants

The excited state dynamics of the blue fluorescent protein (BFP) and its variants, azurite, EBFP1.2, and EBFP2.0, are studied using molecular dynamics (MD) simulations on potential energy surfaces (PESs) generated with the interpolated mechanics-molecular mechanics (IM/MM) scheme. This IM/MM strategy adopts the interpolated PES for an important area of the complex and the conventional force field for the remaining part. We focus on the internal rotation dynamics of the chromophore unit, which is directly related to its fluorescence property, and analyze the time evolutions of the nonrotated chromophore fractions based on trajectories over 10 μs of aggregate simulation time. The characteristics obtained from the calculated time progresses of the nonrotated chromophore fractions in BFP and other variants agree well with experimentally observed properties. The results show that the MD simulation with an IM/MM potential is an attractive approach for studying excited state dynamics of fluorescent proteins in consideration of its efficiency and reliability. We also attempt to investigate the detailed roles that the mutated residues play in delaying the excited state chromophore twisting and thus improving the fluorescence property, and discuss the contributions by the Coulombic and the steric interactions between the chromophore and the mutated residues.

The curious case of the hydrated proton

Understanding the hydrated proton is a critically important problem that continues to engage the research efforts of chemists, physicists, and biologists because of its involvement in a wide array of phenomena. Only recently have several unique properties of the hydrated proton been unraveled through computer simulations. One such process is the detailed molecular mechanism by which protons hop between neighboring water molecules, thus giving rise to the anomalously high diffusion of protons relative to other simple cations. Termed Grotthuss shuttling, this process occurs over multiple time and length scales, presenting unique challenges for computer modeling and simulation. Because the hydrated proton is in reality a dynamical electronic charge defect that spans multiple water molecules, the simulation methodology must be able to dynamically readjust the chemical bonding topology. This reactive nature of the chemical process is automatically captured with ab initio molecular dynamics (AIMD) simulation methods, where the electronic degrees of freedom are treated explicitly. Unfortunately, these calculations can be prohibitively expensive for more complex proton solvation and transport phenomena in the condensed phase. These AIMD simulations remain extremely valuable, however, in validating empirical models, verifying results, and providing insight into molecular mechanisms. In this Account, we discuss recent progress in understanding the solvation and transport properties of the hydrated excess proton. The advances are based on results obtained from reactive molecular dynamics simulations using the multistate empirical valence bond (MS-EVB) methodology. This approach relies on a dynamic linear combination of chemical bond topologies to model charge delocalization and dynamic bonding environments. When parametrized via a variational force-matching algorithm from AIMD trajectories, the MS-EVB method can be viewed as a multiscale bridging of ab initio simulation results to a simpler and more efficient representation. The process allows sampling of longer time and length scales, which would normally be too computationally expensive with AIMD alone. With the MS-EVB methodology, the statistically important components of the excess proton solvation and hopping mechanisms in liquid water have been identified. The most likely solvation structure for the hydrated proton is a distorted Eigen-type complex (H(9)O(4)(+)). In this state, the excess proton charge defect rapidly resonates between three possible distorted Eigen-type structures until a successful proton hop occurs. This process, termed the "special-pair dance", serves as a kind of preparatory phase for the proton hopping while the neighboring water hydrogen-bonding network fluctuates and ultimately rearranges to facilitate a proton hop. The modifications of the solvation structure and transport properties of the excess proton in concentrated acid solutions were further investigated. The Eigen-type solvation structure also possesses both "hydrophilic" and "hydrophobic" sides, which accounts for the affinity of the hydrated proton for the air-water interface. This unusual "amphiphilic" character of the hydrated proton further leads to the metastable formation of contact ion pairs between two hydrated protons. It also engenders a surprisingly constant degree of solubility of hydrophobic species as a function of acid concentration, which contrasts with a markedly variable solubility as a function of salt (such as NaCl or KCl) concentration.

Gradient-based multiconfiguration Shepard interpolation for generating potential energy surfaces for polyatomic reactions

This paper describes and illustrates a way to construct multidimensional representations of reactive potential energy surfaces (PESs) by a multiconfiguration Shepard interpolation (MCSI) method based only on gradient information, that is, without using any Hessian information from electronic structure calculations. MCSI, which is called multiconfiguration molecular mechanics (MCMM) in previous articles, is a semiautomated method designed for constructing full-dimensional PESs for subsequent dynamics calculations (classical trajectories, full quantum dynamics, or variational transition state theory with multidimensional tunneling). The MCSI method is based on Shepard interpolation of Taylor series expansions of the coupling term of a 2 x 2 electronically diabatic Hamiltonian matrix with the diagonal elements representing nonreactive analytical PESs for reactants and products. In contrast to the previously developed method, these expansions are truncated in the present version at the first order, and, therefore, no input of electronic structure Hessians is required. The accuracy of the interpolated energies is evaluated for two test reactions, namely, the reaction OH+H(2)-->H(2)O+H and the hydrogen atom abstraction from a model of alpha-tocopherol by methyl radical. The latter reaction involves 38 atoms and a 108-dimensional PES. The mean unsigned errors averaged over a wide range of representative nuclear configurations (corresponding to an energy range of 19.5 kcal/mol in the former case and 32 kcal/mol in the latter) are found to be within 1 kcal/mol for both reactions, based on 13 gradients in one case and 11 in the other. The gradient-based MCMM method can be applied for efficient representations of multidimensional PESs in cases where analytical electronic structure Hessians are too expensive or unavailable, and it provides new opportunities to employ high-level electronic structure calculations for dynamics at an affordable cost.