Assessing Ion-Water Interactions in the AMOEBA Force Field Using Energy Decomposition Analysis of Electronic Structure Calculations

Evaluating Polarizable Biomembrane Simulations against Experiments

Owing to the increase of available computational capabilities and the potential for providing a more accurate description, polarizable molecular dynamics force fields are gaining popularity in modeling biomolecular systems. It is, however, crucial to evaluate how much precision is truly gained with increasing cost and complexity of the simulation. Here, we leverage the NMRlipids open collaboration and Databank to assess the performance of available polarizable lipid models─the CHARMM-Drude and the AMOEBA-based parameters─against high-fidelity experimental data and compare them to the top-performing nonpolarizable models. While some improvement in the description of ion binding to membranes is observed in the most recent CHARMM-Drude parameters, and the conformational dynamics of AMOEBA-based parameters are excellent, the best nonpolarizable models tend to outperform their polarizable counterparts for each property we explored. The identified shortcomings range from inaccuracies in describing the conformational space of lipids to excessively slow conformational dynamics. Our results provide valuable insights for the further refinement of polarizable lipid force fields and for selecting the best simulation parameters for specific applications.

AMOEBA Force Field Predicts Accurate Hydrogen Bond Counts of Nitriles in SNase by Revealing Water-Protein Interaction in Vibrational Absorption Frequencies

Precisely quantifying the magnitude and direction of electric fields in proteins has long been an outstanding challenge in understanding biological functions. Nitrile vibrational Stark effect probes have been shown to be minimally disruptive to the protein structure and can be better direct reporters of local electrostatic field in the native state of a protein than other measures such as pKa shifts of titratable residues. However, interpretations of the connection between measured vibrational energy and electric field rely on the accurate molecular understanding of interactions of the nitrile group and its environment, particularly from hydrogen bonding. In this work, we compared the extent of hydrogen bonding calculated in two common force fields, the fixed charge force field Amber03 and polarizable force field AMOEBA, at 10 locations of cyanocysteine (CNC) in staphylococcal nuclease (SNase) against the experimental nitrile absorption frequency in terms of full width at half-maximum (FWHM) and frequency temperature line slope (FTLS). We observed that the number of hydrogen bonds correlated well in AMOEBA trajectories with respect to both the FWHM (r = 0.88) and the FTLS (r = -0.85), whereas the correlation of Amber03 trajectories was less reliable because the Amber03 force field predicted more hydrogen bonds in some mutants. Moreover, we demonstrated that contributions from the interactions between CNC and nearby water molecules were significant in AMOEBA trajectories but were not predicted by Amber03. We conclude that although the nitrile absorption peak shape could be qualitatively predicted by the fixed charge Amber03 force field, the detailed electrostatic environment measured by the nitrile probe in terms of the extent of hydrogen bonding could only be accurately observed in the AMOEBA trajectories, where the permanent dipole, quadrupole, and dipole-induced-dipole polarizable interactions were all taken into account. The significance of this finding to the goal of accurately predicting electric fields in complex biomolecular environments is discussed.

Computational Elucidation on the Conformational Control of Selectivity in Intramolecular Ring-Closing Metathesis vs Intermolecular Homometathesis

The ring-closing metathesis reaction of diene plays an important role in the construction of cyclic compounds. In this research, density functional theory (DFT) calculations were conducted to elucidate the mechanisms and origins of the selectivity of ring-closing metathesis and homometathesis. The computational results suggest that the selectivity is determined by the substrate conformation. For the ester-tethered substrate, the homometathesis is more favorable, due to the planar structure of ester facilitating the conjugative effect of the formed E-homometathesis product. For the amide-tethered substrate, the ring-closing metathesis product is the only observed product because the steric hindrance of N-substituents disfavors homometathesis.

Many-Body Effects in Aqueous Systems: Synergies Between Interaction Analysis Techniques and Force Field Development

Interaction analysis techniques, including the many-body expansion (MBE), symmetry-adapted perturbation theory, and energy decomposition analysis, allow for an intuitive understanding of complex molecular interactions. We review these methods by first providing a historical context for the study of many-body interactions and discussing how nonadditivities emerge from Hamiltonians containing strictly pairwise-additive interactions. We then elaborate on the synergy between these interaction analysis techniques and the development of advanced force fields aimed at accurately reproducing the Born-Oppenheimer potential energy surface. In particular, we focus on ab initio-based force fields that aim to explicitly reproduce many-body terms and are fitted to high-level electronic structure results. These force fields generally incorporate many-body effects through (a) parameterization of distributed multipoles, (b) explicit fitting of the MBE, (c) inclusion of many-atom features in a neural network, and (d) coarse-graining of many-body terms into an effective two-body term. We also discuss the emerging use of the MBE to improve the accuracy and speed of ab initio molecular dynamics.

Force Decomposition Analysis: A Method to Decompose Intermolecular Forces into Physically Relevant Component Contributions

Computational quantum chemistry can be more than just numerical experiments when methods are specifically adapted to investigate chemical concepts. One important example is the development of energy decomposition analysis (EDA) to reveal the physical driving forces behind intermolecular interactions. In EDA, typically the interaction energy from a good-quality density functional theory (DFT) calculation is decomposed into multiple additive components that unveil permanent and induced electrostatics, Pauli repulsion, dispersion, and charge-transfer contributions to noncovalent interactions. Herein, we formulate, implement, and investigate decomposing the forces associated with intermolecular interactions into the same components. The resulting force decomposition analysis (FDA) is potentially useful as a complement to the EDA to understand chemistry, while also providing far more information than an EDA for data analysis purposes such as training physics-based force fields. We apply the FDA based on absolutely localized molecular orbitals (ALMOs) to analyze interactions of water with sodium and chloride ions as well as in the water dimer. We also analyze the forces responsible for geometric changes in carbon dioxide upon adsorption onto (and activation by) gold and silver anions. We also investigate how the force components of an EDA-based force field for water clusters, namely MB-UCB, compare to those from force decomposition analysis.

Understanding Charge Dynamics in Dense Electronic Manifolds in Complex Environments

Photoinduced charge transfer (CT) excited states and their relaxation mechanisms can be highly interdependent on the environment effects and the consequent changes in the electronic density. Providing a molecular interpretation of the ultrafast (subpicosecond) interplay between initial photoexcited states in such dense electronic manifolds in condensed phase is crucial for improving and understanding such phenomena. Real-time time-dependent density functional theory is here the method of choice to observe the charge density, explicitly propagated in an ultrafast time domain, along with all time-dependent properties that can be easily extracted from it. A designed protocol of analysis for real-time electronic dynamics to be applied to time evolving electronic density related properties to characterize both in time and in space CT dynamics in complex systems is here introduced and validated, proposing easy to be read cross-correlation maps. As case studies to test such tools, we present the photoinduced charge-transfer electronic dynamics of 5-benzyluracil, a mimic of nucleic acid/protein interactions, and the metal-to-ligand charge-transfer electronic dynamics in water solution of [Ru(dcbpy)2(NCS)2]4-, dcbpy = (4,4'-dicarboxy-2,2'-bipyridine), or "N34-", a dye sensitizer for solar cells. Electrostatic and explicit ab initio treatment of solvent molecules have been compared in the latter case, revealing the importance of the accurate modeling of mutual solute-solvent polarization on CT kinetics. We observed that explicit quantum mechanical treatment of solvent slowed down the charge carriers mobilities with respect to the gas-phase. When all water molecules were modeled instead as simpler embedded point charges, the electronic dynamics appeared enhanced, with a reduced hole-electron distance and higher mean velocities due to the close fixed charges and an artificially increased polarization effect. Such analysis tools and the presented case studies can help to unveil the influence of the electronic manifold, as well as of the finite temperature-induced structural distortions and the environment on the ultrafast charge motions.

Inclusion of High-Field Target Data in AMOEBA's Calibration Improves Predictions of Protein-Ion Interactions

The reliability of molecular mechanics simulations to predict effects of ion binding to proteins depends on their ability to simultaneously describe ion-protein, ion-water, and protein-water interactions. Force fields (FFs) to describe protein-water and ion-water interactions have been constructed carefully and have also been refined routinely to improve accuracy. Descriptions for ion-protein interactions have also been refined, although in an a posteriori manner through the use of "nonbonded-fix (NB-fix)" approaches in which parameters from default Lennard-Jones mixing rules are replaced with those optimized against some reference data. However, even after NB-fix corrections, there remains a significant need for improvement. This is also true for polarizable FFs that include self-consistent inducible moments. Our recent studies on the polarizable AMOEBA FF suggested that the problem associated with modeling ion-protein interactions could be alleviated by recalibrating polarization models of cation-coordinating functional groups so that they respond better to the high electric fields present near ions. Here, we present such a recalibration of carbonyls, carboxylates, and hydroxyls in the AMOEBA protein FF and report that it does improve predictions substantially─mean absolute errors in Na⁺-protein and K⁺-protein interaction energies decrease from 8.7 to 5.3 and 9.6 to 6.3 kcal/mol, respectively. Errors are computed with respect to estimates from van der Waals-inclusive density functional theory benchmarked against high-level quantum mechanical calculations and experiments. While recalibration does improve ion-protein interaction energies, they still remain underestimated, suggesting that further improvements can be made in a systematic manner through modifications in classical formalism. Nevertheless, we show that by applying our many-body NB-fix correction to Lennard-Jones components, these errors are further reduced to 2.7 and 2.6 kcal/mol, respectively, for Na⁺ and K⁺ ions. Finally, we show that the recalibrated AMOEBA protein FF retains its intrinsic reliability in predicting protein structure and dynamics in the condensed phase.

A Combined DFT, Energy Decomposition, and Data Analysis Approach to Investigate the Relationship Between Noncovalent Interactions and Selectivity in a Flexible DABCOnium/Chiral Anion Catalyst System

Developing strategies to study reactivity and selectivity in flexible catalyst systems has become an important topic of research. Herein, we report a combined experimental and computational study aimed at understanding the mechanistic role of an achiral DABCOnium cofactor in a regio- and enantiodivergent bromocyclization reaction. It was found that electron-deficient aryl substituents enable rigidified transition states via an anion-π interaction with the catalyst, which drives the selectivity of the reaction. In contrast, electron-rich aryl groups on the DABCOnium result in significantly more flexible transition states, where interactions between the catalyst and substrate are more important. An analysis of not only the lowest-energy transition state structures but also an ensemble of low-energy transition state conformers via energy decomposition analysis and machine learning was crucial to revealing the dominant noncovalent interactions responsible for observed changes in selectivity in this flexible system.

Classical Exchange Polarization: An Anisotropic Variable Polarizability Model

Polarizability, or the tendency of the electron distribution to distort under an electric field, often depends on the local chemical environment. For example, the polarizability of a chloride ion is larger in gas phase compared to a chloride ion solvated in water. This effect is due to the restriction the Pauli exclusion principle places on the allowed electron states. Because no two electrons can occupy the same state, when a highly polarizable atom comes in close contact with other atoms or molecules, the space of allowed states can dramatically decrease. This constraint suggests that an accurate molecular mechanics polarizability model should depend on the radial distance between neighboring atoms. This paper introduces a variable polarizability model within the framework of the HIPPO (Hydrogen-like Intermolecular Polarizable Potential) force field, by damping the polarizability as a function of the orbital overlap of two atoms. This effectively captures the quantum mechanical exchange polarization effects, without explicit utilization of antisymmetrized wave functions. We show that the variable polarizability model remarkably improves the two-body polarization energies and three-body energies of ion-ion and ion-water systems. Under this model, no manual tuning of atomic polarizabilities for monatomic ions is required; the gas-phase polarizability can be used because an appropriate damping function is able to correct the polarizability at short range.

The Stannum-Ene Reactions of Benzyne and Cyclohexyne with Superb Chemoselectivity for Cyclohexyne

The stannum-ene reactions of both benzyne and cyclohexyne were realized, which is particularly suitable for cyclohexyne with a broad substrate scope and excellent chemoselectivity. Our DFT calculations via distortion/interaction analysis revealed that both stannum- and hydrogen-ene reactions with cyclohexyne have later transition states due to their higher distortion energies in the transition states than those in benzyne reactions, which lead to enhanced Pauli repulsion as the decisive factor in the interaction energy accompanied with enlarged energy gap between two types of ene reactions. Therefore, excellent chemoselectivity was disclosed in the cyclohexyne-ene reaction.

Systematic Evaluation of Ion Diffusion and Water Exchange

As a fundamental property of all fluids, diffusion plays myriad roles in both science and our daily lives. Diffusive properties of many liquids including water have been extensively studied both experimentally and theoretically, while for transition metal ions, there exist significant experimental data that have not been extensively studied theoretically. Hence, high-confidence predictions for challenging systems like radioactive ions that are biohazardous cannot be reliably generated. In this work, a workflow named ISAIAH (Ion Simulation using AMBER for dIffusion Action when Hydrated) was designed to accurately simulate the diffusion coefficients of 15 monoatomic ions with charges varying from -1 to +3 in four water models. As the results indicate, good agreement with experimental values was achieved, leading us to select ²³⁹Pu⁴⁺ (for which no experimental data are available) as a candidate ion to make a theoretical prediction of its diffusion coefficient in water. Among all the force field parameter sets, the ones parametrized using an augmented 12-6-4 Lennard-Jones (LJ) potential showed lower average unsigned errors (AUE) for ions of various radii and electron configurations relative to some 12-6 LJ parameters. This observation agrees well with the fact that diffusion is affected by both the hydration free energy (HFE) and the ion-oxygen distance (IOD) between solute and solvent molecules, both of which are handled well by the 12-6-4 model.

How the electron-deficient Cp ligand facilitates Rh-catalyzed annulations with alkynes

The dominant factors for the CpX ligand effects (Cp* versus CpE) on the reactivity for alkyne insertion into cationic and neutral rhodacycles are identified based on energy decomposition analysis.

Origin of Ligand Effects on Stereoinversion in Pd-Catalyzed Synthesis of Tetrasubstituted Olefins

The mechanism and origin of ligand effects on stereoinversion of Pd-catalyzed synthesis of tetrasubstituted olefins were investigated using DFT calculations and the approach of energy decomposition analysis (EDA). The results reveal that the stereoselectivity-determining steps are different when employing different phosphine ligands. This is mainly due to the steric properties of ligands. With the bulkier Xantphos ligand, the syn/anti-to-Pd 1,2-migrations determine the stereoselectivity. While using the less hindered P(o-tol)₃ ligand, the 1,3-migration is the stereoselectivity-determining step. The EDA results demonstrate that Pauli repulsion and polarization are the dominant factors for controlling the stereochemistry in 1,2- and 1,3-migrations, respectively. The origins of differences of Pauli repulsion and polarization between the two stereoselective transition states are further identified.

Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package

This article summarizes technical advances contained in the fifth major release of the Q-Chem quantum chemistry program package, covering developments since 2015. A comprehensive library of exchange-correlation functionals, along with a suite of correlated many-body methods, continues to be a hallmark of the Q-Chem software. The many-body methods include novel variants of both coupled-cluster and configuration-interaction approaches along with methods based on the algebraic diagrammatic construction and variational reduced density-matrix methods. Methods highlighted in Q-Chem 5 include a suite of tools for modeling core-level spectroscopy, methods for describing metastable resonances, methods for computing vibronic spectra, the nuclear-electronic orbital method, and several different energy decomposition analysis techniques. High-performance capabilities including multithreaded parallelism and support for calculations on graphics processing units are described. Q-Chem boasts a community of well over 100 active academic developers, and the continuing evolution of the software is supported by an "open teamware" model and an increasingly modular design.

Rational Design of Nonbonded Point Charge Models for Divalent Metal Cations with Lennard-Jones 12-6 Potential

Exploring a metal-involved biochemical process at a molecular level often requires a reliable description of metal properties in aqueous solution by classical nonbonded models. An additional C₄ term for considering ion-induced dipole interactions was previously proposed to supplement the widely used Lennard-Jones 12-6 potential (known as the 12-6-4 LJ-type model) with good accuracy. Here, we demonstrate an alternative to modeling divalent metal cations (M²⁺) with the traditional 12-6 LJ potential by developing nonbonded point charge models for use with 11 water models: TIP3P, SPC/E, SPC/E_b, TIP4P-Ew, TIP4P-D, and TIP4P/2005 and the more recent OPC3, TIP3P-FB, OPC, TIP4P-FB, and a99SB-disp. Our designed models simultaneously reproduce the experimental hydration free energy, ion-oxygen distance, and coordination number in the first hydration shell accurately for most of the metal cations, an accuracy equivalent to that of the complex 12-6-4 LJ-type and double exponential potential models. A systematic comparison with the existing M²⁺ models is presented as well in terms of effective ion radii, diffusion constants, water exchange rates, and ion-water interactions. Molecular dynamics simulations of metal substitution in Escherichia coli glyoxalase I variants show the great potential of our new models for metalloproteins.

A Quantum Chemical Topology Picture of Intermolecular Electrostatic Interactions and Charge Penetration Energy

Based on the Interacting Quantum Atoms approach, we present herein a conceptual and theoretical framework of short-range electrostatic interactions, whose accurate description is still a challenging problem in molecular modeling. For all the noncovalent complexes in the S66 database, the fragment-based and atomic decomposition of the electrostatic binding energies is performed using both the charge density of the dimers and the unrelaxed densities of the monomers. This energy decomposition together with dispersion corrections gives rise to a pairwise approximation to the total binding energy. It also provides energetic descriptors at varying distance that directly address the atomic and molecular electrostatic interactions as described by point-charge or multipole-based potentials. Additionally, we propose a consistent definition of the charge penetration energy within quantum chemical topology, which is mainly characterized in terms of the intramolecular electrostatic energy. Finally, we discuss some practical implications of our results for the design and validation of electrostatic potentials.

A Computational Mechanistic Study of Cp*Co(III)-Catalyzed Three-Component C-H Bond Addition to Terpenes and Formaldehydes: Insights into the Origins of Regioselectivity

Transition metal-catalyzed three-component reactions of arenes, dienes, and carbonyls enable the convergent synthesis of homoallylic alcohols. Controlling regioselectivity is a central challenge for the difunctionalization of substituted 1,3-dienes in which multiple unbiased C═C bonds exist. Here, the mechanisms of Cp*Co(III)-catalyzed three-component C-H bond addition to terpenes and formaldehydes were investigated by density functional theory calculations. The reaction proceeds via sequential C(sp²)-H activation, migratory insertion, β-hydride elimination, hydride reinsertion, and C-C bond formation to yield the final product. The migratory insertion is the rate- and regioselectivity-determining step of the overall reaction. We employed an energy decomposition approach to quantitatively dissect the contributions of different types of interactions to regioselectivity. For the 2-alkyl substituted 1,3-dienes, the orbital interactions in the 3,4-insertion are intrinsically more favorable as compared to that in the 4,3-insertion, while the stronger steric effects between metallacycle and 1,3-diene override the intrinsic electronic preference. However, the steric effects failed to rationalize the unfavorable 1,2-insertion that is analogous to 4,3-insertion and even bears smaller steric effects. The donor-acceptor interaction analysis indicates that orbital interactions between σ_Co-C and π_C═C decreased significantly in the 1,2-insertion transition state, which leads to higher activation energy barriers. These insights into the dominant effects controlling regioselectivity will enable rational design of new catalysts for selective functionalization of dienes.

Analysis of DNA interactions and GC content with energy decomposition in large-scale quantum mechanical calculations

GC content is a contributing factor to the stability of nucleic acids due to hydrogen bonding. More hydrogen bonding generally results in greater stability. Empirical evidence, however, has suggested that the GC content of a nucleic acid is a poor predictor of its stability, implying that there are sequence-dependent interactions besides what its GC content indicates. To examine how much such sequence-dependent interactions affect the interaction energies of double-stranded DNA (dsDNA) molecules, dsDNA molecules of different sequences are generated and examined in silico for variabilities in the interaction energies within each group of dsDNA molecules of the same GC content. Since the amount of hydrogen bonding depends on the GC content, holding the GC content fixed when examining the differences in interaction energies allows sequence-dependent interactions to be isolated. The nature of sequence-dependent interactions is then dissected using energy decomposition analysis (EDA). By using EDA, the components of the interactions that depend on the neighboring base pairs help explain some of the variability in the interaction energies of the dsDNA molecules despite having the same GC content. This work provides a new paradigm and tool for the study and analysis of the distributions of interaction components in dsDNA with the same GC content using EDA within large-scale quantum chemistry calculations.

From Intermolecular Interaction Energies and Observable Shifts to Component Contributions and Back Again: A Tale of Variational Energy Decomposition Analysis

Quantum chemistry in the form of density functional theory (DFT) calculations is a powerful numerical experiment for predicting intermolecular interaction energies. However, no chemical insight is gained in this way beyond predictions of observables. Energy decomposition analysis (EDA) can quantitatively bridge this gap by providing values for the chemical drivers of the interactions, such as permanent electrostatics, Pauli repulsion, dispersion, and charge transfer. These energetic contributions are identified by performing DFT calculations with constraints that disable components of the interaction. This review describes the second-generation version of the absolutely localized molecular orbital EDA (ALMO-EDA-II). The effects of different physical contributions on changes in observables such as structure and vibrational frequencies upon complex formation are characterized via the adiabatic EDA. Example applications include red- versus blue-shifting hydrogen bonds; the bonding and frequency shifts of CO, N₂, and BF bound to a [Ru(II)(NH₃)₅]^{2 +} moiety; and the nature of the strongly bound complexes between pyridine and the benzene and naphthalene radical cations. Additionally, the use of ALMO-EDA-II to benchmark and guide the development of advanced force fields for molecular simulation is illustrated with the recent, very promising, MB-UCB potential.

Electrostatic repulsion-controlled regioselectivity in nitrene-mediated allylic C–H amidations

The difference of repulsive electrostatic interactions between nitrene N and allyl carbon atoms is the dominant factor affecting the regioselectivity in metal nitrenoid-catalyzed allylic C–H amidations.

Origins of Lewis acid acceleration in nickel-catalysed C–H, C–C and C–O bond cleavage

The effects of charge transfer, Pauli repulsion and electrostatics/polarization are identified as dominant factors for Lewis acid accelerations in Ni-catalyzed C–X (X = H, C and O) bond cleavages.

Perylene bisimide cyclophanes as receptors for planar transition structures – catalysis of stereoinversions by shape-complementarity and noncovalent π–π interactions

Perylene bisimide cyclophanes, ideal receptors for planar aromatic compounds, act as π–π catalysts by stabilizing shape-complementary stereoinversion transition structures.

Consistent inclusion of continuum solvation in energy decomposition analysis: theory and application to molecular CO2 reduction catalysts

To facilitate computational investigation of intermolecular interactions in the solution phase, we report the development of ALMO-EDA(solv), a scheme that allows the application of continuum solvent models within the framework of energy decomposition analysis (EDA) based on absolutely localized molecular orbitals (ALMOs). In this scheme, all the quantum mechanical states involved in the variational EDA procedure are computed with the presence of solvent environment so that solvation effects are incorporated in the evaluation of all its energy components. After validation on several model complexes, we employ ALMO-EDA(solv) to investigate substituent effects on two classes of complexes that are related to molecular CO2 reduction catalysis. For [FeTPP(CO2-κC)]2- (TPP = tetraphenylporphyrin), we reveal that two ortho substituents which yield most favorable CO2 binding, -N(CH3)3 + (TMA) and -OH, stabilize the complex via through-structure and through-space mechanisms, respectively. The coulombic interaction between the positively charged TMA group and activated CO2 is found to be largely attenuated by the polar solvent. Furthermore, we also provide computational support for the design strategy of utilizing bulky, flexible ligands to stabilize activated CO2 via long-range Coulomb interactions, which creates biomimetic solvent-inaccessible "pockets" in that electrostatics is unscreened. For the reactant and product complexes associated with the electron transfer from the p-terphenyl radical anion to CO2, we demonstrate that the double terminal substitution of p-terphenyl by electron-withdrawing groups considerably strengthens the binding in the product state while moderately weakens that in the reactant state, which are both dominated by the substituent tuning of the electrostatics component. These applications illustrate that this new extension of ALMO-EDA provides a valuable means to unravel the nature of intermolecular interactions and quantify their impacts on chemical reactivity in solution.

Systematic Parametrization of Divalent Metal Ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB Water Models

Divalent metal ions play important roles in biological and materials systems. Molecular dynamics simulation is an efficient tool to investigate these systems at the microscopic level. Recently, four new water models (OPC3, OPC, TIP3P-FB, and TIP4P-FB) have been developed and better represent the physical properties of water than previous models. Metal ion parameters are dependent on the water model employed, making it necessary to develop metal ion parameters for select new water models. In the present work, we performed parameter scanning for the 12-6 Lennard-Jones nonbonded model of divalent metal ions in conjunction with the four new water models as well as four previous water models (TIP3P, SPC/E, TIP4P, and TIP4P-Ew). We found that these new three-point and four-point water models provide comparable or significantly improved performance for the simulation of divalent metal ions when compared to previous water models in the same category. Among all eight water models, the OPC3 water model yields the best performance for the simulation of divalent metal ions in the aqueous phase when using the 12-6 model. On the basis of the scanning results, we independently parametrized the 12-6 model for 24 divalent metal ions with each of the four new water models. As noted previously, the 12-6 model still fails to simultaneously reproduce the experimental hydration free energy (HFE) and ion-oxygen distance (IOD) values even with these new water models. To solve this problem, we parametrized the 12-6-4 model for the 16 divalent metal ions for which we have both experimental HFE and IOD values for each of the four new water models. The final parameters are able to reproduce both the experimental HFE and IOD values accurately. To validate the transferability of our parameters, we carried out benchmark calculations to predict the energies and geometries of ion-water clusters as well as the ion diffusivity coefficient of Mg²⁺. By comparison to quantum chemical calculations and experimental data, these results show that our parameters are well designed and have excellent transferability. The metal ion parameters for the 12-6 and 12-6-4 models reported herein can be employed in simulations of various biological and materials systems when using the OPC3, OPC, TIP3P-FB, or TIP4P-FB water model.

Probing radical-molecule interactions with a second generation energy decomposition analysis of DFT calculations using absolutely localized molecular orbitals

Intermolecular interactions between radicals and closed-shell molecules are ubiquitous in chemical processes, ranging from the benchtop to the atmosphere and extraterrestrial space. While energy decomposition analysis (EDA) schemes for closed-shell molecules can be generalized for studying radical-molecule interactions, they face challenges arising from the unique characteristics of the electronic structure of open-shell species. In this work, we introduce additional steps that are necessary for the proper treatment of radical-molecule interactions to our previously developed unrestricted Absolutely Localized Molecular Orbital (uALMO)-EDA based on density functional theory calculations. A "polarize-then-depolarize" (PtD) scheme is used to remove arbitrariness in the definition of the frozen wavefunction, rendering the ALMO-EDA results independent of the orientation of the unpaired electron obtained from isolated fragment calculations. The contribution of radical rehybridization to polarization energies is evaluated. It is also valuable to monitor the wavefunction stability of each intermediate state, as well as their associated spin density profiles, to ensure the EDA results correspond to a desired electronic state. These radical extensions are incorporated into the "vertical" and "adiabatic" variants of uALMO-EDA for studies of energy changes and property shifts upon complexation. The EDA is validated on two model complexes, H₂O˙F and FH˙OH. It is then applied to several chemically interesting radical-molecule complexes, including the sandwiched and T-shaped benzene dimer radical cation, complexes of pyridine with benzene and naphthalene radical cations, binary and ternary complexes of the hydroxyl radical with water (˙OH(H₂O) and ˙OH(H₂O)₂), and the pre-reactive complexes and transition states in the ˙OH + HCHO and ˙OH + CH₃CHO reactions. These examples suggest that this second generation uALMO-EDA is a useful tool for furthering one's understanding of both energetic and property changes associated with radical-molecule interactions.

Hierarchical phenomena in multicomponent liquids: simulation methods, analysis, chemistry

Complex, multicomponent, solutions have often been studied solely through the lens of specific applications of interest. Yet advances to both simulation methodologies (enhanced sampling, etc.) and analysis techniques (network analysis algorithms and others), are creating a trove of data that reveal transcending characteristics across vast compositional phase space. This perspective discusses technical considerations of the reliable and accurate simulations of complex solutions, followed by the advances to analysis algorithms that elucidate coupling of different length and timescale behavior (hierarchical phenomena). The different manifestations of hierarchical phenomena are presented across an array of solution environments, emphasizing fundamental and ongoing science questions. With a more advanced molecular understanding in hand, a quintessential application (solvent extraction) is discussed, where significant opportunities exist to re-imagine the technical scope of an established technology.

Developing an effective polarizable bond method for small molecules with application to optimized molecular docking

Electrostatic interaction plays an essential role in protein-ligand binding. Due to the polarization effect, electrostatic interactions are largely impacted by their local environments. However, traditional force fields use fixed point charge-charge interactions to describe electrostatic interactions but is unable to include the polarization effect. The lack of the polarization effect in the force field representation can result in substantial error in biomolecular studies, such as molecular dynamics and molecular docking. Docking programs usually employ traditional force fields to estimate the binding energy between a ligand and a protein for pose selection or scoring. The intermolecular interaction energy mainly consists of van der Waals and electrostatic interaction in the force field representation. In the current study, we developed an Effective Polarizable Bond (EPB) method for small organic molecules and applied this EPB method to optimize protein-ligand docking in computational tests for a variety of protein-ligand systems. We tested the method on a set of 38 cocrystallized structures taken from the Protein Data Bank (PDB) and found that the maximum error was reduced from 7.98 Å to 2.03 Å when using EPB Dock, providing strong evidence that the use of EPB charges is important. We found that our optimized docking approach with EPB charges could improve the docking performance, sometimes dramatically, and the maximum error was reduced from 12.88 Å to 1.57 Å in Optimized Docking (in the case of 1fqx). The average RMSD decreased from 2.83 Å to 1.85 Å. Further investigations showed that the use of the EBP method could enhance intermolecular hydrogen bonding, which is a major contributing factor to improved docking performance. Developed tools for the calculation of the polarized ligand charge from a protein-ligand complex structure with the EPB method are freely available on GitHub (https://github.com/Xundrug/EPB).

Polarisable force fields: what do they add in biomolecular simulations?

The quality of biomolecular simulations critically depends on the accuracy of the force field used to calculate the potential energy of the molecular configurations. Currently, most simulations employ non-polarisable force fields, which describe electrostatic interactions as the sum of Coulombic interactions between fixed atomic charges. Polarisation of these charge distributions is incorporated only in a mean-field manner. In the past decade, extensive efforts have been devoted to developing simple, efficient, and yet generally applicable polarisable force fields for biomolecular simulations. In this review, we summarise the latest developments in accounting for key biomolecular interactions with polarisable force fields and applications to address challenging biological questions. In the end, we provide an outlook for future development in polarisable force fields.

Electronic transitions for a fully polarizable QM/MM approach based on fluctuating charges and fluctuating dipoles: Linear and corrected linear response regimes

The fully polarizable Quantum Mechanics/Molecular Mechanics (QM/MM) approach based on fluctuating charges and fluctuating dipoles, named QM/FQFμ [T. Giovannini et al., J. Chem. Theory Comput. 15, 2233 (2019)], is extended to the calculation of vertical excitation energies of solvated molecular systems. Excitation energies are defined within two different solvation regimes, i.e., linear response (LR), where the response of the MM portion is adjusted to the QM transition density, and corrected-Linear Response (cLR) in which the MM response is adjusted to the relaxed QM density, thus being able to account for charge equilibration in the excited state. The model, which is specified in terms of three physical parameters (electronegativity, chemical hardness, and polarizability) is applied to vacuo-to-water solvatochromic shifts of aqueous solutions of para-nitroaniline, pyridine, and pyrimidine. The results show a good agreement with their experimental counterparts, thus highlighting the potentialities of this approach.

Energy Decomposition Analyses Reveal the Origins of Catalyst and Nucleophile Effects on Regioselectivity in Nucleopalladation of Alkenes

Nucleopalladation is one of the most common mechanisms for Pd-catalyzed hydro- and oxidative functionalization of alkenes. Due to the electronic bias of the π-alkene-palladium complexes, nucleopalladations with terminal aliphatic alkenes typically deliver the nucleophile to the more substituted sp² carbon to form the Markovnikov-selective products. The selective formation of the anti-Markovnikov nucleopalladation products requires the inherent electronic effects to be overridden, which is still a significant challenge for reactions with simple aliphatic alkenes. Because the interactions between the nucleophile and the alkene substrate are influenced by a complex combination of multiple types of steric and electronic effects, a thorough understanding of the interplay of these underlying interactions is needed to rationalize and predict the regioselectivity. Here, we employ an energy decomposition approach to quantitatively separate the different types of nucleophile-substrate interactions, including steric, electrostatic, orbital interactions, and dispersion effects, and to predict the impacts of each factor on regioselectivity. We demonstrate the use of this approach on the origins of catalyst-controlled anti-Markovnikov-selectivity in Hull's Pd-catalyzed oxidative amination reactions. In addition, we evaluated the regioselectivity in a series of nucleopalladation reactions with different neutral and anionic Pd catalysts and N- and O-nucleophiles with different steric and electronic properties. On the basis of these computational analyses, a generalized scheme is established to identify the dominant nucleophile-substrate interaction affecting the regioselectivity of nucleopalladations with different Pd catalysts and nucleophiles.

Probing Blue-Shifting Hydrogen Bonds with Adiabatic Energy Decomposition Analysis

The physical origin of blue-shifting hydrogen bonds remains a subject of debate, although many plausible explanations have been proposed. Using a molecular property decomposition analysis based on absolutely localized molecular orbitals, we investigated several representative F₃CH···Y (Y = H₂O, NH₃, Cl^-) complexes. We reveal that features of a blue-shifting H-bond already appear on the frozen surface where both polarization and charge transfer (CT) are "turned off", and that the final frequency shift observed depends on the strength of CT. Further decomposition of forces at the frozen level shows that Pauli repulsion is the only component that shortens the C-H bond in the short-range, while both permanent electrostatics and dispersion lengthen the bond. The effects of these forces from the medium to long-range are also discussed. Our analysis provides a complete picture for blue-shifting H-bonds and suggests two necessary conditions for their features to be observed at equilibrium structures: (i) stronger Pauli repulsion than the combination of electrostatic and dispersion forces; (ii) relatively weak CT that is insufficient to compensate for the blue-shifting effect of the frozen interaction.

Polarizable Force Fields for Biomolecular Simulations: Recent Advances and Applications

Realistic modeling of biomolecular systems requires an accurate treatment of electrostatics, including electronic polarization. Due to recent advances in physical models, simulation algorithms, and computing hardware, biomolecular simulations with advanced force fields at biologically relevant timescales are becoming increasingly promising. These advancements have not only led to new biophysical insights but also afforded opportunities to advance our understanding of fundamental intermolecular forces. This article describes the recent advances and applications, as well as future directions, of polarizable force fields in biomolecular simulations.

Mutually polarizable QM/MM model with in situ optimized localized basis functions

We extend our recently developed quantum-mechanical/molecular mechanics (QM/MM) approach [Dziedzic et al., J. Chem. Phys. 145, 124106 (2016)] to enable in situ optimization of the localized orbitals. The quantum subsystem is described with onetep linear-scaling density functional theory and the classical subsystem - with the AMOEBA polarizable force field. The two subsystems interact via multipolar electrostatics and are fully mutually polarizable. A total energy minimization scheme is employed for the Hamiltonian of the coupled QM/MM system. We demonstrate that, compared to simpler models using fixed basis sets, the additional flexibility offered by in situ optimized basis functions improves the accuracy of the QM/MM interface, but also poses new challenges, making the QM subsystem more prone to overpolarization and unphysical charge transfer due to increased charge penetration. We show how these issues can be efficiently solved by replacing the classical repulsive van der Waals term for QM/MM interactions with an interaction of the electronic density with a fixed, repulsive MM potential that mimics Pauli repulsion, together with a modest increase in the damping of QM/MM polarization. We validate our method, with particular attention paid to the hydrogen bond, in tests on water-ion pairs, the water dimer, first solvation shells of neutral and charged species, and solute-solvent interaction energies. As a proof of principle, we determine suitable repulsive potential parameters for water, K⁺, and Cl^-. The mechanisms we employed to counteract the unphysical overpolarization of the QM subsystem are demonstrated to be adequate, and our approach is robust. We find that the inclusion of explicit polarization in the MM part of QM/MM improves agreement with fully QM calculations. Our model permits the use of minimal size QM regions and, remarkably, yields good energetics across the well-balanced QM/MM interface.