A Many-Body, Fully Polarizable Approach to QM/MM Simulations

PyDFT-QMMM: A modular, extensible software framework for DFT-based QM/MM molecular dynamics

PyDFT-QMMM is a Python-based package for performing hybrid quantum mechanics/molecular mechanics (QM/MM) simulations at the density functional level of theory. The program is designed to treat short-range and long-range interactions through user-specified combinations of electrostatic and mechanical embedding procedures within periodic simulation domains, providing necessary interfaces to external quantum chemistry and molecular dynamics software. To enable direct embedding of long-range electrostatics in periodic systems, we have derived and implemented force terms for our previously described QM/MM/PME approach [Pederson and McDaniel, J. Chem. Phys. 156, 174105 (2022)]. Communication with external software packages Psi4 and OpenMM is facilitated through Python application programming interfaces (APIs). The core library contains basic utilities for running QM/MM molecular dynamics simulations, and plug-in entry-points are provided for users to implement custom energy/force calculation and integration routines, within an extensible architecture. The user interacts with PyDFT-QMMM primarily through its Python API, allowing for complex workflow development with Python scripting, for example, interfacing with PLUMED for free energy simulations. We provide benchmarks of forces and energy conservation for the QM/MM/PME and alternative QM/MM electrostatic embedding approaches. We further demonstrate a simple example use case for water solute in a water solvent system, for which radial distribution functions are computed from 100 ps QM/MM simulations; in this example, we highlight how the solvation structure is sensitive to different basis-set choices due to under- or over-polarization of the QM water molecule's electron density.

Seamless integration of GEM, a density based-force field, for QM/MM simulations via LICHEM, Psi4, and Tinker-HP

Hybrid quantum mechanics/molecular mechanics (QM/MM) simulations have become an essential tool in computational chemistry, particularly for analyzing complex biological and condensed phase systems. Building on this foundation, our work presents a novel implementation of the Gaussian Electrostatic Model (GEM), a polarizable density-based force field, within the QM/MM framework. This advancement provides seamless integration, enabling efficient and optimized QM/GEM calculations in a single step using the LICHEM Code. We have successfully applied our implementation to water dimers and hexamers, demonstrating the ability to handle water systems with varying numbers of water molecules. Moreover, we have extended the application to describe the double proton transfer of the aspartic acid dimer in a box of water, which highlights the method's proficiency in investigating heterogeneous systems. Our implementation offers the flexibility to perform on-the-fly density fitting or to utilize pre-fitted coefficients to estimate exchange and Coulomb contributions. This flexibility enhances efficiency and accuracy in modeling molecular interactions, especially in systems where polarization effects are significant.

Assessing Implicit and Explicit Polarizable Solvation Models for Nuclear-Electronic Orbital Systems: Quantum Proton Polarization and Solvation Energetics

Accurate simulations of many chemical processes require the inclusion of both nuclear quantum effects and a solvent environment. The nuclear-electronic orbital (NEO) approach, which treats electrons and select nuclei quantum mechanically on the same level, combined with a polarizable continuum model (PCM) for the solvent environment, addresses this challenge in a computationally practical manner. In this work, the NEO-PCM approach is extended beyond the IEF-PCM (integral equation formalism PCM) and C-PCM (conductor PCM) approaches to the SS(V)PE (surface and simulation of volume polarization for electrostatics) and ddCOSMO (domain decomposed conductor-like screening model) approaches. IEF-PCM, SS(V)PE, C-PCM, and ddCOSMO all exhibit similar solvation energies as well as comparable nuclear polarization within the NEO framework. The calculations show that the nuclear density does not leak out of the molecular cavity because it is much more localized than the electronic density. Finally, the polarization of quantized protons is analyzed in both continuum solvent and explicit solvent environments described by the polarizable MB-pol model, illustrating the impact of specific hydrogen-bonding interactions captured only by explicit solvation. These calculations highlight the relationship among solvation formalism, nuclear polarization, and energetics.

The application of QM/MM simulations in heterogeneous catalysis

The QM/MM simulation method is provenly efficient for the simulation of biological systems, where an interplay of extensive environment and delicate local interactions drives a process of interest through a funnel on a complex energy landscape. Recent advances in quantum chemistry and force-field methods present opportunities for the adoption of QM/MM to simulate heterogeneous catalytic processes, and their related systems, where similar intricacies exist on the energy landscape. Herein, the fundamental theoretical considerations for performing QM/MM simulations, and the practical considerations for setting up QM/MM simulations of catalytic systems, are introduced; then, areas of heterogeneous catalysis are explored where QM/MM methods have been most fruitfully applied. The discussion includes simulations performed for adsorption processes in solvent at metallic interfaces, reaction mechanisms within zeolitic systems, nanoparticles, and defect chemistry within ionic solids. We conclude with a perspective on the current state of the field and areas where future opportunities for development and application exist.

DFT-based QM/MM with Particle-Mesh Ewald for Direct, Long-Range Electrostatic Embedding

We present a DFT-based, QM/MM implementation with long-range electrostatic embedding achieved by direct real-space integration of the particle mesh Ewald (PME) computed electrostatic potential. The key transformation is the interpolation of the electrostatic potential from the PME grid to the DFT quadrature grid, from which integrals are easily evaluated utilizing standard DFT machinery. We provide benchmarks of the numerical accuracy with choice of grid size and real-space corrections, and demonstrate that good convergence is achieved while introducing nominal computational overhead. Furthermore, the approach requires only small modification to existing software packages, as is demonstrated with our implementation in the OpenMM and Psi4 software. After presenting convergence benchmarks, we evaluate the importance of long-range electrostatic embedding in three solute/solvent systems modeled with QM/MM. Water and BMIM/BF4 ionic liquid were considered as ``simple' and ``complex' solvents respectively, with water and p-phenylenediamine (PPD) solute molecules treated at QM level of theory. While electrostatic embedding with standard real-space truncation may introduce negligible error for simple systems such as water solute in water solvent, errors become more significant when QM/MM is applied to complex solvents such as ionic liquids. An extreme example is the electrostatic embedding energy for oxidized PPD in BMIM/BF4 for which real-space truncation produces severe error even at 2-3 nm cutoff distances. This latter example illustrates that utilization of QM/MM to compute redox potentials within concentrated electrolytes/ionic media requires carefully chosen long-range electrostatic embedding algorithms, with our presented algorithm providing a general and robust approach.

QM/MM Energy Decomposition Using the Interacting Quantum Atoms Approach

The interacting quantum atoms (IQA) method decomposes the quantum mechanical (QM) energy of a molecular system in terms of one- and two-center (atomic) contributions within the context of the quantum theory of atoms in molecules. Here, we demonstrate that IQA, enhanced with molecular mechanics (MM) and Poisson-Boltzmann surface-area (PBSA) solvation methods, is naturally extended to the realm of hybrid QM/MM methodologies, yielding intra- and inter-residue energy terms that characterize all kinds of covalent and noncovalent bonding interactions. To test the robustness of this approach, both metal-water interactions and QM/MM boundary artifacts are characterized in terms of the IQA descriptors derived from QM regions of varying size in Zn(II)- and Mg(II)-water clusters. In addition, we analyze a homologous series of inhibitors in complex with a matrix metalloproteinase (MMP-12) by carrying out QM/MM-PBSA calculations on their crystallographic structures followed by IQA energy decomposition. Overall, these applications not only show the advantages of the IQA QM/MM approach but also address some of the challenges lying ahead for expanding the QM/MM methodology.

Flexible boundary layer using exchange for embedding theories. I. Theory and implementation

Embedding theory is a powerful computational chemistry approach to exploring the electronic structure and dynamics of complex systems, with Quantum Mechanical/Molecular Mechanics (QM/MM) being the prime example. A challenge arises when trying to apply embedding methodology to systems with diffusible particles, e.g., solvents, if some of them must be included in the QM region, for example, in the description of solvent-supported electronic states or reactions involving proton transfer or charge-transfer-to-solvent: without a special treatment, inter-diffusion of QM and MM particles will eventually lead to a loss of QM/MM separation. We have developed a new method called Flexible Boundary Layer using Exchange (FlexiBLE) that solves the problem by adding a biasing potential to the system that closely maintains QM/MM separation. The method rigorously preserves ensemble averages by leveraging their invariance to an exchange of identical particles. With a careful choice of the biasing potential and the use of a tree algorithm to include only important QM and MM exchanges, we find that the method has an MM-forcefield-like computational cost and thus adds negligible overhead to a QM/MM simulation. Furthermore, we show that molecular dynamics with the FlexiBLE bias conserves total energy, and remarkably, sub-diffusional dynamical quantities in the inner QM region are unaffected by the applied bias. FlexiBLE thus widens the range of chemistry that can be studied with embedding theory.

General Many-Body Framework for Data-Driven Potentials with Arbitrary Quantum Mechanical Accuracy: Water as a Case Study

We present a general framework for the development of data-driven many-body (MB) potential energy functions (MB-QM PEFs) that represent the interactions between small molecules at an arbitrary quantum-mechanical (QM) level of theory. As a demonstration, a family of MB-QM PEFs for water is rigorously derived from density functionals belonging to different rungs across Jacob's ladder of approximations within density functional theory (MB-DFT) and from Møller-Plesset perturbation theory (MB-MP2). Through a systematic analysis of individual MB contributions to the interaction energies of water clusters, we demonstrate that all MB-QM PEFs preserve the same accuracy as the corresponding ab initio calculations, with the exception of those derived from density functionals within the generalized gradient approximation (GGA). The differences between the DFT and MB-DFT results are traced back to density-driven errors that prevent GGA functionals from accurately representing the underlying molecular interactions for different cluster sizes and hydrogen-bonding arrangements. We show that this shortcoming may be overcome, within the MB formalism, by using density-corrected functionals (DC-DFT) that provide a more consistent representation of each individual MB contribution. This is demonstrated through the development of a MB-DFT PEF derived from DC-PBE-D3 data, which more accurately reproduce the corresponding ab initio results.

Development and application of quantum mechanics/molecular mechanics methods with advanced polarizable potentials

Quantum mechanics/molecular mechanics (QM/MM) simulations are a popular approach to study various features of large systems. A common application of QM/MM calculations is in the investigation of reaction mechanisms in condensed-phase and biological systems. The combination of QM and MM methods to represent a system gives rise to several challenges that need to be addressed. The increase in computational speed has allowed the expanded use of more complicated and accurate methods for both QM and MM simulations. Here, we review some approaches that address several common challenges encountered in QM/MM simulations with advanced polarizable potentials, from methods to account for boundary across covalent bonds and long-range effects, to polarization and advanced embedding potentials.

Biomolecular QM/MM Simulations: What Are Some of the "Burning Issues"?

QM/MM simulations have become an indispensable tool in many chemical and biochemical investigations. Considering the tremendous degree of success, including recognition by a 2013 Nobel Prize in Chemistry, are there still "burning challenges" in QM/MM methods, especially for biomolecular systems? In this short Perspective, we discuss several issues that we believe greatly impact the robustness and quantitative applicability of QM/MM simulations to many, if not all, biomolecules. We highlight these issues with observations and relevant advances from recent studies in our group and others in the field. Despite such limited scope, we hope the discussions are of general interest and will stimulate additional developments that help push the field forward in meaningful directions.