Variational and perturbative formulations of quantum mechanical/molecular mechanical free energy with mean-field embedding and its analytical gradients

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.

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.

Ab Initio Evaluation of the Redox Potential of Cytochrome c

Various biochemical activities of metabolism and biosynthesis are fulfilled by redox processes with explicit electron exchange, which furnish redox enzymes with high chemical reactivity. However, theoretical investigation of a redox process, which simultaneously involves a complex electronic change at a redox metal center and conformational reorganization of the surrounding protein environment coupled to the electronic change, requires computationally conflicting approaches, highly accurate quantum chemical calculations, and long-time molecular dynamics (MD) simulations, limiting the physicochemical understanding of biological redox processes. Here, we theoretically examined a redox process of cytochrome c by means of a hybrid molecular simulation technique, which enables one to consistently treat the redox center at the ab initio quantum chemistry level of theory and the protein reorganization with long-time MD simulations on the microsecond timescale. The calculations successfully evaluated a large absolute redox potential, 4.34 eV, with errors of only 0.03 to 0.34 eV to the experimental ones without any problem-specific empirical parameters. Through the long-time MD sampling, large and nonlinear reorganization of the protein environment was unveiled and the molecular determinants for the redox potential were identified. The present ab initio approach significantly expands the applicability of theoretical investigation to biological redox systems with more electronically complicated redox centers such as polynuclear transition metal complexes.

Quantum Chemical and QM/MM Models in Biochemistry

Quantum chemical (QC) calculations provide a basis for deriving a microscopic understanding of enzymes and photobiological systems. Here we describe how QC models can be used to explore the electronic structure, dynamics, and energetics of biomolecules. We introduce the hybrid quantum mechanics/classical mechanics (QM/MM) approach, where a quantum mechanically described system of interest is embedded in a classically described force field representation of the biochemical surroundings. We also discuss the QM cluster model approach, as well as embedding theories, that provide complementary methodologies to model quantum mechanical effects in biomolecules. The chapter also provides some practical guides for building quantum biochemical models using the quinone reduction catalysis in respiratory complex I and a model reaction in solution as examples.

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.

Free-energy analysis of physisorption on solid-liquid interface with the solution theory in the energy representation

Physisorption of urea on its crystal in contact with water was subject to energetics analysis with all-atom molecular dynamics simulation. The transfer free energy of urea to an adsorption site was treated in the framework of the energy-representation theory of solutions, which allows a fast computation of the free energy in an inhomogeneous environment with solid-liquid interface. The preference of adsorption was then compared between the (001) and (110) faces, and it was found that the physisorption is more favorable on (001) than on (110) in correspondence to the hydrogen bonding between the adsorbed urea and the crystal urea. Among the terrace configurations of adsorption, the attractive interaction governs the preferable site with a minor role of the repulsive interaction. The effect of an edge was also treated by examining the terrace and step and was shown to be strongly operative on the (110) face when the CO group of the adsorbed urea points toward the edge. The present work demonstrates that the solution theory can be a framework for analyzing the energetics of physisorption and addressing the roles of the crystal and liquid at the interface through the systematic decomposition of free energy.

Structure and permeability of ionomers studied by atomistic molecular simulation combined with the theory of solutions in the energy representation

Ionomers play a key role in forming the catalyst layer of polymer electrolyte fuel cells. In the present work, we performed atomistic molecular dynamics simulations and free-energy calculations with the energy-representation method for sulfonated polyethersulfone (SPES) and its derivatives toward the rational design of ionomers for carbon alloy catalysts. It was observed that H₂O aggregates strongly in the branched SPES systems with fluorocarbons and is located homogeneously in the systems without fluorocarbons. The O₂ permeability was then examined within the framework of the solubility-diffusion mechanism. The permeability was seen to be large for the branched SPES with fluorocarbons, indicating that the performance of ionomers as a permeation medium for O₂ may be tuned by the flexibility and branching of the polymer chain.

Probabilistic analysis for identifying the driving force of protein folding

Toward identifying the driving force of protein folding, energetics was analyzed in water for Trp-cage (20 residues), protein G (56 residues), and ubiquitin (76 residues) at their native (folded) and heat-denatured (unfolded) states. All-atom molecular dynamics simulation was conducted, and the hydration effect was quantified by the solvation free energy. The free-energy calculation was done by employing the solution theory in the energy representation, and it was seen that the sum of the protein intramolecular (structural) energy and the solvation free energy is more favorable for a folded structure than for an unfolded one generated by heat. Probabilistic arguments were then developed to determine which of the electrostatic, van der Waals, and excluded-volume components of the interactions in the protein-water system governs the relative stabilities between the folded and unfolded structures. It was found that the electrostatic interaction does not correspond to the preference order of the two structures. The van der Waals and excluded-volume components were shown, on the other hand, to provide the right order of preference at probabilities of almost unity, and it is argued that a useful modeling of protein folding is possible on the basis of the excluded-volume effect.

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.

Computing conformational free energy differences in explicit solvent: An efficient thermodynamic cycle using an auxiliary potential and a free energy functional constructed from the end points

Many biomolecules undergo conformational changes associated with allostery or ligand binding. Observing these changes in computer simulations is difficult if their timescales are long. These calculations can be accelerated by observing the transition on an auxiliary free energy surface with a simpler Hamiltonian and connecting this free energy surface to the target free energy surface with free energy calculations. Here, we show that the free energy legs of the cycle can be replaced with energy representation (ER) density functional approximations. We compute: (1) The conformational free energy changes for alanine dipeptide transitioning from the right-handed free energy basin to the left-handed basin and (2) the free energy difference between the open and closed conformations of β-cyclodextrin, a "host" molecule that serves as a model for molecular recognition in host-guest binding. β-cyclodextrin contains 147 atoms compared to 22 atoms for alanine dipeptide, making β-cyclodextrin a large molecule for which to compute solvation free energies by free energy perturbation or integration methods and the largest system for which the ER method has been compared to exact free energy methods. The ER method replaced the 28 simulations to compute each coupling free energy with two endpoint simulations, reducing the computational time for the alanine dipeptide calculation by about 70% and for the β-cyclodextrin by > 95%. The method works even when the distribution of conformations on the auxiliary free energy surface differs substantially from that on the target free energy surface, although some degree of overlap between the two surfaces is required. © 2016 Wiley Periodicals, Inc.

Free-energy analysis of protein solvation with all-atom molecular dynamics simulation combined with a theory of solutions

The structure of a protein is strongly influenced by solvation. In the present review, the solvation effect is treated within the framework of statistical thermodynamics. The key quantity is the solvation free energy of protein and a fast method for its computation with explicit solvent is introduced. The applications of the method to the structural fluctuation of protein in water, structure determination of protein-protein complex, and urea effect on protein structure are then described, and the roles of solvent water and cosolvent are discussed from the standpoint of energetics.

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.

Correlation analysis for heat denaturation of Trp-cage miniprotein with explicit solvent

Energetics was analyzed for Trp-cage miniprotein in water to elucidate the solvation effect in heat denaturation. The solvation free energy was computed for a set of protein structures at room and high temperatures with all-atom molecular dynamics simulation combined with the solution theory in the energy representation, and its correlations were investigated against the intramolecular (structural) energy of the protein and the average interaction energy of the protein with the solvent water. It was observed both at room and high temperatures that the solvation free energy is anticorrelated to the structural energy and varies in parallel to the electrostatic component of the protein-water interaction energy without correlations to the van der Waals and excluded-volume components. When the set of folded structures sampled at room temperature was compared with the set of unfolded ones at high temperature, it was found that the preference order of the two sets is in correspondence to the van der Waals and excluded-volume components in the sum of the protein intramolecular and protein-water intermolecular interactions and is not distinguished by the electrostatic component.

Multiple environment single system quantum mechanical/molecular mechanical (MESS-QM/MM) calculations. 1. Estimation of polarization energies

In combined quantum mechanical/molecular mechanical (QM/MM) free energy calculations, it is often advantageous to have a frozen geometry for the quantum mechanical (QM) region. For such multiple-environment single-system (MESS) cases, two schemes are proposed here for estimating the polarization energy: the first scheme, termed MESS-E, involves a Roothaan step extrapolation of the self-consistent field (SCF) energy; whereas the other scheme, termed MESS-H, employs a Newton–Raphson correction using an approximate inverse electronic Hessian of the QM region (which is constructed only once). Both schemes are extremely efficient, because the expensive Fock updates and SCF iterations in standard QM/MM calculations are completely avoided at each configuration. They produce reasonably accurate QM/MM polarization energies: MESS-E can predict the polarization energy within 0.25 kcal/mol in terms of the mean signed error for two of our test cases, solvated methanol and solvated β-alanine, using the M06-2X or ωB97X-D functionals; MESS-H can reproduce the polarization energy within 0.2 kcal/mol for these two cases and for the oxyluciferin–luciferase complex, if the approximate inverse electronic Hessians are constructed with sufficient accuracy.

Exploring solvent effects upon the Menshutkin reaction using a polarizable force field

The energetics of the Menshutkin reaction between triethylamine and ethyl iodide have been computed using B3LYP and MP2 with the LANL2DZ, LANL2DZd, SVP, MIDI!, 6-311G(d,p), and aug-cc-PVTZ basis sets. Small- and large-core energy-consistent relativistic pseudopotentials were employed. Solvent effect corrections were computed from QM/MM Monte Carlo simulations utilizing free-energy perturbation theory, PDDG/PM3, and both a nonpolarizable OPLS and polarizable OPLS-AAP force field. The B3LYP/MIDI! theory level provided the best DeltaG(++) values with a mean absolute error (MAE) of 4.9 kcal/mol from experiment in cyclohexane, CCl(4), THF, DMSO, acetonitrile, water, and methanol. However, the relative rates in cyclohexane, and to a certain extent CCl(4), were determined to be greatly underestimated when using the nonpolarizable OPLS force field. An overall reduction in the MAE to 3.1 kcal/mol using B3LYP/MIDI!/OPLS-AAP demonstrated the need for a fully polarizable force field when computing solvent effects for highly dipolar transition structures in low-dielectric media. The MAEs obtained with PDDG/PM3/OPLS and OPLS-AAP of 5.3 and 3.8 kcal/mol, respectively, provided comparable results to B3LYP at a fraction of the computational resources. The large rate accelerations observed in the reaction were correlated to an increased stabilization of the emerging charge separation at the transition state via favorable solute-solvent interactions.

A wave-function based approach for polarizable charge model: Systematic comparison of polarization effects on protic, aprotic, and ionic liquids

We first describe a wave-function based formalism of polarizable charge model by starting from the Hartree product ansatz for the total wave function and making the second-order expansion of individual molecular energies with the use of partial charge operators. The resulting model is shown to be formally equivalent to the charge response kernel model that starts from the linear-response approximation to partial charges, and also closely related to a family of fluctuating charge models that are based on the electronegativity equalization principle. We then apply the above model to a systematic comparison of polarization effects on qualitatively different liquids, namely, protic solvents (water and methanol), an aprotic polar solvent (acetonitrile), and imidazolium-based ionic liquids. Electronic polarization is known to decelerate molecular motions in conventional solvents while it accelerates them in ionic liquids. To obtain more insights into these phenomena, we consider an effective decomposition of total polarization energy into molecular contributions, and show that their statistical distribution is well-correlated with the acceleration/deceleration of molecular motions. In addition, we perform effective nonpolarizable simulations based on mean polarized charges, and compare them with fully polarizable simulations. The result shows that the former can reproduce structural properties of conventional solvents rather accurately, while they fail qualitatively to reproduce acceleration of molecular motions in ionic liquids.