Accurate and Efficient Treatment of Continuous Solute Charge Density in the Mean-Field QM/MM Free Energy Calculation

Accelerated constant-voltage quantum mechanical/molecular mechanical method for molecular systems at electrochemical interfaces

The structure and electronic properties of a molecule at an electrochemical interface are changed by interactions with the electrode surface and the electrolyte solution, which can be significantly modulated by an applied voltage. We present an efficient self-consistent quantum mechanics/molecular mechanics (QM/MM) approach to study a physisorbed molecule at a metal electrode-electrolyte interface under the constant-voltage condition. The approach employs a classical polarizable double electrode model, which enables us to study the QM/MM system in the constant-voltage ensemble. A mean-field embedding approximation is further introduced in order to overcome the difficulties associated with statistical sampling of the electrolyte configurations. The results of applying the method to a test system indicate that the adsorbed molecule is no less or slightly more polarized at the interface than in the bulk electrolyte solution. The geometry of the horizontally adsorbed molecule is modulated by their electrostatic interactions with the polarizable electrode surfaces and also the interactions with cations attracted toward the interface when the adsorbate is reduced. We also demonstrate that the approach can be used to quantitatively evaluate the reorganization energy of a one electron reduction reaction of a molecule in an electrochemical cell.

Thermodynamics of π–π Interactions of Benzene and Phenol in Water

The π–π interaction is a major driving force that stabilizes protein assemblies during protein folding. Recent studies have additionally demonstrated its involvement in the liquid–liquid phase separation (LLPS) of intrinsically disordered proteins (IDPs). As the participating residues in IDPs are exposed to water, π–π interactions for LLPS must be modeled in water, as opposed to the interactions that are often established at the hydrophobic domains of folded proteins. Thus, we investigated the association of free energies of benzene and phenol dimers in water by integrating van der Waals (vdW)-corrected density functional theory (DFT) and DFT in classical explicit solvents (DFT-CES). By comparing the vdW-corrected DFT and DFT-CES results with high-level wavefunction calculations and experimental solvation free energies, respectively, we established the quantitative credibility of these approaches, enabling a reliable prediction of the benzene and phenol dimer association free energies in water. We discovered that solvation influences dimer association free energies, but not significantly when no direct hydrogen-bond-type interaction exists between two monomeric units, which can be explained by the enthalpy–entropy compensation. Our comprehensive computational study of the solvation effect on π–π interactions in water could help us understand the molecular-level driving mechanism underlying the IDP phase behaviors.

Doubly Polarized QM/MM with Machine Learning Chaperone Polarizability

A major shortcoming of semiempirical (SE) molecular orbital methods is their severe underestimation of molecular polarizability compared with experimental and ab initio (AI) benchmark data. In a combined quantum mechanical and molecular mechanical (QM/MM) treatment of solution-phase reactions, solute described by SE methods therefore tends to generate inadequate electronic polarization response to solvent electric fields, which often leads to large errors in free energy profiles. To address this problem, here we present a hybrid framework that improves the response property of SE/MM methods through high-level molecular-polarizability fitting. Specifically, we place on QM atoms a set of corrective polarizabilities (referred to as chaperone polarizabilities), whose magnitudes are determined from machine learning (ML) to reproduce the condensed-phase AI molecular polarizability along the minimum free energy path. These chaperone polarizabilities are then used in a machinery similar to a polarizable force field calculation to compensate for the missing polarization energy in the conventional SE/MM simulations. Because QM atoms in this treatment host SE wave functions as well as classical polarizabilities, both polarized by MM electric fields, we name this method doubly polarized QM/MM (dp-QM/MM). We demonstrate the new method on the free energy simulations of the Menshutkin reaction in water. Using AM1/MM as a base method, we show that ML chaperones greatly reduce the error in the solute molecular polarizability from 6.78 to 0.03 ų with respect to the density functional theory benchmark. The chaperone correction leads to ∼10 kcal/mol of additional polarization energy in the product region, bringing the simulated free energy profiles to closer agreement with the experimental results. Furthermore, the solute-solvent radial distribution functions show that the chaperone polarizabilities modify the free energy profiles through enhanced solvation corrections when the system evolves from the charge-neutral reactant state to the charge-separated transition and product states. These results suggest that the dp-QM/MM method, enabled by ML chaperone polarizabilities, provides a very physical remedy for the underpolarization problem in SE/MM-based free energy simulations.

Recent Progress toward Ab Initio Modeling of Electrocatalysis

Electrode potential is the key factor for controlling electrocatalytic reactions at electrochemical interfaces, and moreover, it is also known that the pH and solutes (e.g., cations) of the solution have prominent effects on electrocatalysis. Understanding these effects requires microscopic information on the electrochemical interfaces, in which theoretical simulations can play an important role. This Perspective summarizes the recent progress in method development for modeling electrochemical interfaces, including different methods for describing the electrolytes at the interfaces and different schemes for charging up the electrode surfaces. In the final section, we provide an outlook for future development in modeling methods and their applications to electrocatalysis.

Reaction Path-Force Matching in Collective Variables: Determining Ab Initio QM/MM Free Energy Profiles by Fitting Mean Force

First-principles determination of free energy profiles for condensed-phase chemical reactions is hampered by the daunting costs associated with configurational sampling on ab initio quantum mechanical/molecular mechanical (AI/MM) potential energy surfaces. Here, we report a new method that enables efficient AI/MM free energy simulations through mean force fitting. In this method, a free energy path in collective variables (CVs) is first determined on an efficient reactive aiding potential. Based on the configurations sampled along the free energy path, correcting forces to reproduce the AI/MM forces on the CVs are determined through force matching. The AI/MM free energy profile is then predicted from simulations on the aiding potential in conjunction with the correcting forces. Such cycles of correction-prediction are repeated until convergence is established. As the instantaneous forces on the CVs sampled in equilibrium ensembles along the free energy path are fitted, this procedure faithfully restores the target free energy profile by reproducing the free energy mean forces. Due to its close connection with the reaction path-force matching (RP-FM) framework recently introduced by us, we designate the new method as RP-FM in collective variables (RP-FM-CV). We demonstrate the effectiveness of this method on a type-II solution-phase SN2 reaction, NH3 + CH3Cl (the Menshutkin reaction), simulated with an explicit water solvent. To obtain the AI/MM free energy profiles, we employed the semiempirical AM1/MM Hamiltonian as the base level for determining the string minimum free energy pathway, along which the free energy mean forces are fitted to various target AI/MM levels using the Hartree-Fock (HF) theory, density functional theory (DFT), and the second-order Møller-Plesset perturbation (MP2) theory as the AI method. The forces on the bond-breaking and bond-forming CVs at both the base and target levels are obtained by force transformation from Cartesian to redundant internal coordinates under the Wilson B-matrix formalism, where the linearized FM is facilitated by the use of spline functions. For the Menshutkin reaction tested, our FM treatment greatly reduces the deviations on the CV forces, originally in the range of 12-33 to ∼2 kcal/mol/Å. Comparisons with the experimental and benchmark AI/MM results, tests of the new method under a variety of simulation protocols, and analyses of the solute-solvent radial distribution functions suggest that RP-FM-CV can be used as an efficient, accurate, and robust method for simulating solution-phase chemical reactions.

Controlling potential difference between electrodes based on self-consistent-charge density functional tight binding

A proper understanding and description of the electronic response of the electrode surfaces in electrochemical systems are quite important because the interactions between the electrode surface and electrolyte give rise to unique and useful interfacial properties. Atomistic modeling of the electrodes requires not only an accurate description of the electronic response under a constant-potential condition but also computational efficiency in order to deal with systems large enough to investigate the interfacial electrolyte structures. We thus develop a self-consistent-charge density functional tight binding based method to model a pair of electrodes in electrochemical cells under the constant-potential condition. The method is more efficient than the (ab initio) density functional theory calculations so that it can treat systems as large as those studied in classical atomistic simulations. It can also describe the electronic response of electrodes quantum mechanically and more accurately than the classical counterparts. The constant-potential condition is introduced through a Legendre transformation of the electronic energy with respect to the difference in the number of electrons in the two electrodes and their electrochemical potential difference, through which the Kohn-Sham equations for each electrode are variationally derived. The method is applied to platinum electrodes faced parallel to each other under an applied voltage. The electronic response to the voltage and a charged particle is compared with the result of a classical constant-potential method based on the chemical potential equalization principle.

Enthalpy-Entropy Interplay in π-Stacking Interaction of Benzene Dimer in Water

Aromatic groups can engage in an interesting class of noncovalent interactions termed π-π interactions, which play a pivotal role in stabilizing a variety of molecular architectures, including nucleic acids, proteins, and supramolecular assemblies. When the aromatic compounds interact with each other in an aqueous environment, their association is facilitated by the hydrophobic effect-the trend of nonpolar solutes to aggregate in a polar solution. To develop an in-depth understanding of hydrophobic association, we investigate in the present work π-π interactions in water, employing as a paradigm the benzene dimer. Using DFT-CES, a mean-field QM/MM method recently developed by our group, we describe the benzene solute at a quantum-mechanical level. Full consideration of detailed solute-electron density enables an optimal description of the solute-solvent interactions, leading to an accurate prediction of hydration free energies. In π-stacking of benzene, we find an entropic stabilization associated with the shrinkage of the solvent-excluded volume, which agrees with the theory of hydrophobic effect at subnanoscales. However, at the equilibrium binding distance of the benzene dimer, we find that the entropic stabilization nearly cancels out due to the enthalpic cost required for dewetting the internal space. Such an enthalpy-entropy compensation leads the association free energy to be predominantly dictated by the solute-solute interaction enthalpy. The present work offers new insight into the mechanistic role of water and the primary thermodynamic driving force of hydrophobic association.

Hydration Thermodynamics of Non-Polar Aromatic Hydrocarbons: Comparison of Implicit and Explicit Solvation Models

The precise description of solute-water interactions is essential to understand the chemo-physical nature in hydration processes. Such a hydration thermodynamics for various solutes has been explored by means of explicit or implicit solvation methods. Using the Poisson-Boltzmann solvation model, the implicit models are well designed to reasonably predict the hydration free energies of polar solutes. The implicit model, however, is known to have shortcomings in estimating those for non-polar aromatic compounds. To investigate a cause of error, we employed a novel systematic framework of quantum-mechanical/molecular-mechanical (QM/MM) coupling protocol in explicit solvation manner, termed DFT-CES, based on the grid-based mean-field treatment. With the aid of DFT-CES, we delved into multiple energy parts, thereby comparing DFT-CES and PB models component-by-component. By applying the modified PB model to estimate the hydration free energies of non-polar solutes, we find a possibility to improve the predictability of PB models. We expect that this study could shed light on providing an accurate route to study the hydration thermodynamics for various solute compounds.

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.

Averaged Condensed Phase Model for Simulating Molecules in Complex Environments

The need for configurational sampling dramatically increases the cost of combined quantum mechanics/molecular mechanics (QM/MM) simulations of chemical processes in solution. We developed an averaged condensed phase environment (ACPE) model that constructs an effective polarizable environment directly from explicitly sampled molecular dynamics configurations via the K-means++ algorithm and a mathematically rigorous translation of the molecular mechanics parameters. The model captures detailed heterogeneous features in the environment that may be difficult to describe using a conventional polarizable continuum model. Instead of performing repeated QM/MM calculations for each new configuration of the environment, the ACPE approach allows one to perform a single QM calculation on an averaged configuration. Here, we demonstrate the model by computing electronic excitation energies for several small molecules in solution. The ACPE model predicts the excitation energies in excellent agreement with conventional configurational averaging yet with orders of magnitude of reduction in the computational cost.

Classical and semiclassical dynamics in statistical environments with a mixed dynamical and statistical representation

We present a basic theory to study real-time dynamics embedded in a large environment that is treated using a statistical method. In light of great progress in the molecular-level studies on time-resolved spectroscopies, chemical reaction dynamics, and so on, not only in the gas phase but also in condensed phases like liquid solvents and even in crowded environments in living cells, we need to bridge over a gap between statistical mechanics and microscopic real-time dynamics. For instance, an analogy to gas-phase dynamics in which molecules are driven by the gradient of the potential energy hyper-surfaces (PESs) suggests that particles in condensed phases should run on the free energy surface instead. The question is whether this anticipation is correct. To answer it, we here propose a mixed dynamics and statistical representation to treat chemical dynamics embedded in a statistical ensemble. We first define the entropy functional, which is a function of the phase-space position of the dynamical subsystem, being dressed with statistical weights from the statistical counterpart. We then consider the functionals of temperature, free energy, and chemical potential as their extensions in statistical mechanics, through which one can clarify the relationship between real-time microscopic dynamics and statistical quantities. As an illustrative example we show that molecules in the dynamical subsystem should run on the free-energy functional surface, if and only if the spatial gradients of the temperature functional are all zero. Otherwise, additional forces emerge from the gradient of the temperature functional. Numerical demonstrations are presented at the very basic level of this theory of molecular dissociation in atomic cluster solvents.

Solvatochromic Shift of Brooker's Merocyanine: Hartree-Fock Exchange in Time Dependent Density Functional Calculation and Hydrogen Bonding Effect

The Brooker's merocyanine exhibits a large hypsochromic shift from an apolar aprotic solvent to a polar protic solvent. Quantum chemical calculations have been performed to study the solvatochromism, but there remained a discrepancy between the calculated and experimental solvatochromic shifts. In this paper we evaluate quantum mechanically the excitation energies of the Brooker's merocyanine in water, methanol, acetonitrile, and dichloromethane to investigate what are important factors to accurately model the solvatochromism of the dye by using TDDFT in combination with implicit and explicit solvation models including the PCM, PCMSMD, RISM-SCF-SEDD, and mean-field QM/MM. The results severely depend on the density functional, especially on the amount of Hartree-Fock exchange included in the functional. Furthermore, an explicit description of the solute-solvent hydrogen bonds makes a non-negligible contribution to the shift. The experimental large solvatochromic shift can be accurately reproduced by the TDDFT/RISM-SCF-SEDD and mean-field QM/MM calculations with the LC-BOP functional, although the excitation energies in solutions are considerably overestimated. We also estimated the excitation energies and the solvatochromic shift at the SAC-CI/RISM-SCF-SEDD and mean-field QM/MM level, which are in very good agreement with the experimental values. These results indicate that if an explicit solvent model is used, TDDFT calculations using such a long-range corrected functional can accurately model the solvatochromism. However, an ab initio quantum chemical method including sufficient electron correlation effects is required to reproduce not only the solvatochromism but also the excitation energies in solutions.

Treating electrostatics with Wolf summation in combined quantum mechanical and molecular mechanical simulations

The Wolf summation approach [D. Wolf et al., J. Chem. Phys. 110, 8254 (1999)], in the damped shifted force (DSF) formalism [C. J. Fennell and J. D. Gezelter, J. Chem. Phys. 124, 234104 (2006)], is extended for treating electrostatics in combined quantum mechanical and molecular mechanical (QM/MM) molecular dynamics simulations. In this development, we split the QM/MM electrostatic potential energy function into the conventional Coulomb r(-1) term and a term that contains the DSF contribution. The former is handled by the standard machinery of cutoff-based QM/MM simulations whereas the latter is incorporated into the QM/MM interaction Hamiltonian as a Fock matrix correction. We tested the resulting QM/MM-DSF method for two solution-phase reactions, i.e., the association of ammonium and chloride ions and a symmetric SN2 reaction in which a methyl group is exchanged between two chloride ions. The performance of the QM/MM-DSF method was assessed by comparing the potential of mean force (PMF) profiles with those from the QM/MM-Ewald and QM/MM-isotropic periodic sum (IPS) methods, both of which include long-range electrostatics explicitly. For ion association, the QM/MM-DSF method successfully eliminates the artificial free energy drift observed in the QM/MM-Cutoff simulations, in a remarkable agreement with the two long-range-containing methods. For the SN2 reaction, the free energy of activation obtained by the QM/MM-DSF method agrees well with both the QM/MM-Ewald and QM/MM-IPS results. The latter, however, requires a greater cutoff distance than QM/MM-DSF for a proper convergence of the PMF. Avoiding time-consuming lattice summation, the QM/MM-DSF method yields a 55% reduction in computational cost compared with the QM/MM-Ewald method. These results suggest that, in addition to QM/MM-IPS, the QM/MM-DSF method may serve as another efficient and accurate alternative to QM/MM-Ewald for treating electrostatics in condensed-phase simulations of chemical reactions.

A hybrid framework of first principles molecular orbital calculations and a three-dimensional integral equation theory for molecular liquids: multi-center molecular Ornstein-Zernike self-consistent field approach

In this study, we reported the development of a new quantum mechanics/molecular mechanics (QM/MM)-type framework to describe chemical processes in solution by combining standard molecular-orbital calculations with a three-dimensional formalism of integral equation theory for molecular liquids (multi-center molecular Ornstein-Zernike (MC-MOZ) method). The theoretical procedure is very similar to the 3D-reference interaction site model self-consistent field (RISM-SCF) approach. Since the MC-MOZ method is highly parallelized for computation, the present approach has the potential to be one of the most efficient procedures to treat chemical processes in solution. Benchmark tests to check the validity of this approach were performed for two solute (solute water and formaldehyde) systems and a simple SN2 reaction (Cl(-) + CH3Cl → ClCH3 + Cl(-)) in aqueous solution. The results for solute molecular properties and solvation structures obtained by the present approach were in reasonable agreement with those obtained by other hybrid frameworks and experiments. In particular, the results of the proposed approach are in excellent agreements with those of 3D-RISM-SCF.

Recent advances in QM/MM free energy calculations using reference potentials

Background

Recent years have seen enormous progress in the development of methods for modeling (bio)molecular systems. This has allowed for the simulation of ever larger and more complex systems. However, as such complexity increases, the requirements needed for these models to be accurate and physically meaningful become more and more difficult to fulfill. The use of simplified models to describe complex biological systems has long been shown to be an effective way to overcome some of the limitations associated with this computational cost in a rational way.

Scope of review

Hybrid QM/MM approaches have rapidly become one of the most popular computational tools for studying chemical reactivity in biomolecular systems. However, the high cost involved in performing high-level QM calculations has limited the applicability of these approaches when calculating free energies of chemical processes. In this review, we present some of the advances in using reference potentials and mean field approximations to accelerate high-level QM/MM calculations. We present illustrative applications of these approaches and discuss challenges and future perspectives for the field.

Major conclusions

The use of physically-based simplifications has shown to effectively reduce the cost of high-level QM/MM calculations. In particular, lower-level reference potentials enable one to reduce the cost of expensive free energy calculations, thus expanding the scope of problems that can be addressed.

General significance

As was already demonstrated 40 years ago, the usage of simplified models still allows one to obtain cutting edge results with substantially reduced computational cost. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

•

We present some of the advances to accelerate high-level QM/MM calculations.

•

Quantitative limitations of low-level methods can be overcome by these approaches.

•

Reference potentials make free energy simulations feasible for large systems.

•

Automated fitting reduces the need of expensive sampling of high-level approaches.

•

Application of reference potentials can be extended to a wide range of processes.