Energy decomposition analysis in solution based on the fragment molecular orbital method

We develop the pair interaction energy decomposition analysis (PIEDA) in solution by combining the fragment molecular orbital (FMO) method with the polarizable continuum model (PCM). The solvent screening of the electrostatic interaction and the desolvation penalty in complex formation are described by this approach from ab initio calculations of fragments and their pairs. The applications to the complex of solvated sodium and chlorine ions, as well as to lysine and aspartic acid, show how the analysis helps reveal the physical picture. The PIEDA/PCM method is also applied to a small protein chignolin (PDB: 1UAO), and the solvent screening of the pair interactions is discussed.     

An energy decomposition analysis for intramolecular non-covalent interaction in solvated environment

In this work, the intra-EDA method, which is a recently developed energy decomposition analysis scheme for intramolecular non-covalent interaction is extended from gas phase to solvated environment. It is the first analysis scheme that performs analysis for intramolecular interaction in solution. By fragmentation scheme, a molecule is divided into intramolecular interacting fragments and environmental fragments via single bond homolysis breaking. The solvent effect is taken into account by implicit solvation model. Intramolecular interaction free energy is estimated as the separated treatment of inter-fragment interactions in dielectric environment. The analysis results highlight the importance of solvent effects to intramolecular non-covalent interaction.     

Energy decomposition analysis for intramolecular non-covalent interaction in solvated environment

In this work, the intra-EDA method, which is a recently developed energy decomposition analysis scheme for intramolecular non-covalent interaction is extended from gas phase to solvated environment. It is the first analysis scheme that performs analysis for intramolecular interaction in solution. By fragmentation scheme, a molecule is divided into intramolecular interacting fragments and environmental fragments via single bond homolysis breaking. The solvent effect is taken into account by implicit solvation model. Intramolecular interaction free energy is estimated as the separated treatment of inter-fragment interactions in dielectric environment. The analysis results highlight the importance of solvent effects to intramolecular non-covalent interaction.     

Atomic decomposition of the protein solvation free energy and its application to amyloid-beta protein in water

We report the development of an atomic decomposition method of the protein solvation free energy in water, which ascribes global change in the solvation free energy to local changes in protein conformation as well as in hydration structure. So far, empirical decomposition analyses based on simple continuum solvation models have prevailed in the study of protein-protein interactions, protein-ligand interactions, as well as in developing scoring functions for computer-aided drug design. However, the use of continuum solvation model suffers serious drawbacks since it yields the protein free energy landscape which is quite different from that of the explicit solvent model and since it does not properly account for the non-polar hydrophobic effects which play a crucial role in biological processes in water. Herein, we develop an exact and general decomposition method of the solvation free energy that overcomes these hindrances. We then apply this method to elucidate the molecular origin for the solvation free energy change upon the conformational transitions of 42-residue amyloid-beta protein (Aβ42) in water, whose aggregation has been implicated as a primary cause of Alzheimer's disease. We address why Aβ42 protein exhibits a great propensity to aggregate when transferred from organic phase to aqueous phase.     

Enthalpies of Adduct Formation between Boron Trifluoride and Selected Organic Bases in Solution: Toward an Accurate Theoretical Entry to Lewis Basicity

The Lewis basicity of selected organic bases, modeled by the enthalpies of adduct formation between gaseous BF₃ and bases in dichloromethane (DCM) solution, is critically examined. Although experimental enthalpies for a large number of molecules have been reported in the literature, it may be desirable to estimate missing or uncertain data for important Lewis bases. We decided to use high-level ab initio procedures, combined with a polarized continuum solvation model, in which the solvated species were the clusters formed by specific hydrogen bonding of DCM with the Lewis base and the Lewis base/BF₃ adduct. This mode of interaction with DCM corresponds to a specific solvation model (SSM). The results essentially showed that the enthalpy of BF₃ adduct formation in DCM solution was clearly influenced by specific interactions, with DCM acting as hydrogen-bonding donor (HBD) molecule in two ways: base/DCM and adduct/DCM, confirming that specific solvation is an important contribution to experimentally determined Lewis basicity scales. This analysis allowed us to conclude that there are reasons to suspect some gas-phase values to be in error by more than the stated experimental uncertainty. Some experimental values in DCM solution that were uncertain for identified reasons could be complemented by the computed values.     

Combining a polarizable force-field and a coarse-grained polarizable solvent model. II. Accounting for hydrophobic effects

A revised and improved version of our efficient polarizable force-field/coarse grained solvent combined approach (Masella, Borgis, and Cuniasse, J. Comput. Chem. 2008, 29, 1707) is described. The polarizable pseudo-particle solvent model represents the macroscopic solvent polarization by induced dipoles placed on mobile pseudo-particles. In this study, we propose a new formulation of the energy term handling the nonelectrostatic interactions among the pseudo-particles. This term is now able to reproduce the energetic and structural response of liquid water due to the presence of a hydrophobic spherical cavity. Accordingly, the parameters of the energy term handling the nonpolar solute/solvent interactions have been refined to reproduce the free-solvation energy of small solutes, based on a standard thermodynamic integration scheme. The reliability of this new approach has been checked for the properties of solvated methane and of the solvated methane dimer, as well as by performing 10 × 20 ns molecular dynamics (MD) trajectories for three solvated proteins. A long-time stability of the protein structures along the trajectories is observed. Moreover, our method still provides a measure of the protein solvation thermodynamic at the same accuracy as standard Poisson-Boltzman continuum methods. These results show the relevance of our approach and its applicability to massively coupled MD schemes to accurately and intensively explore solvated macromolecule potential energy surfaces.     

Performance of the DFT-D method, paired with the PCM implicit solvation model, for the computation of interaction energies of solvated complexes of biological interest

In this work we investigate the performance of the DFT method, augmented with an empirical dispersion function (DFT-D), paired with the PCM implicit solvation model, for the computation of noncovalent interaction energies of biologically-relevant, solvated model complexes. It is found that this method describes intermolecular interactions within water and ether (protein-like) environments with roughly the same accuracy as in the gas phase. Another important finding is that, when environmental effects are taken into account, the empirical dispersion term associated with the DFT-D method need be modified very little (or not at all), in order to obtain the optimum, most well balanced, performance.     

Solvent effects on glycine. I. A supermolecule modeling of tautomerization via intramolecular proton transfer

The relative stabilities of glycine tautomers involved in the intramolecular proton transfer are investigated computationally by considering glycine-water complexes containing up to five water molecules. The supermolecule results are compared with continuum calculations. Specific solute-solvent interactions and solvent induced changes in the solute wave function are considered using the natural bond orbitals (NBO) method. The stabilization of the zwitterion upon solvation is explained by the changes in the wave functions localized on the forming and breaking bonds as well as by the different interaction energies in the zwitterionic and neutral clusters. Only the neutral species exist in mono- and dihydrated clusters and in the gas phase. In the smaller clusters, zwitterions are mainly stabilized by conformational effects, whereas in larger clusters, in particular when glycine is solvated on both sides of its heavy atom backbone, polarization effects dominate the stability of a given tautomer. Generally, the strength of the solute-solvent interactions is governed by the intermolecular charge transfer interactions. As the solvation progresses, the hypothetical gaseous zwitterion is better solvated than the gaseous neutral, making zwitterion to neutral tautomerization progressively less exothermic for clusters containing up to three water molecules, and endothermic for larger clusters. The neutral isomer does not exist for some solvent arrangements with five water molecules. Only solvent arrangements in which water molecules do not interact with the reactive proton are considered. Hence, the experimentally observed double well potential energy surface may be due to such an interaction or to a different reaction mechanism.     

A new method to determine electrostatic potential around a macromolecule in solution from molecular wave functions

The three-dimensional reference interaction site model integral equation theory (3D-RISM) combined with the ab initio molecular orbital method (3D-RISM-SCF) is applied to a solvated macromolecular system. The solvation structure around a solute molecule is obtained from the 3D-RISM integral equation under the electrostatic potential of the solute molecule, calculated by the ab initio molecular orbital theory. The electrostatic potential should be calculated on each grid point in the three-dimensional real space. Therefore, the calculation of the electrostatic potential is the most time consuming part in this method. In this article, we propose a new procedure to save the computational cost for calculating the electrostatic potential and the solvated fock matrix. The strategy of this procedure is to evaluate the electrostatic potential and the solvated fock matrix in different ways, depending on the distance between solute and solvent. Inside the repulsive cores of solute atoms, it is possible to avoid the calculation of electrostatic potential and solvated Fock matrix by assuming the potential to be infinity. In the region sufficiently far from solute, they are evaluated classically by putting the effective point charge on each atom. In the intermediate region, the electrostatic potential is evaluated directly by integrating the molecular orbitals of the solute molecule. The electronic structure and the energy gradient of Methionine-Enkephalin and solvation structure are estimated by using this procedure in aqueous solution, and are compared with the results from other procedures. The results are compared also with those from the continuum model.     

A study on orientation and absorption spectrum of interfacial molecules by using continuum model

In this work, a numerical procedure based on the continuum model is developed and applied to the solvation energy for ground state and the spectral shift against the position and the orientation of the interfacial molecule. The interface is described as a sharp boundary separating two bulk media. The polarizable continuum model (PCM) allows us to account for both electrostatic and nonelectrostatic solute-solvent interactions when we calculate the solvation energy. In this work we extend PCM to the interfacial system and the information about the position and orientation of the interfacial molecule can be obtained. Based on the developed expression of the electrostatic free energy of a nonequilibrium state, the numerical procedure has been implemented and used to deal with a series of test molecules. The time-dependent density functional theory (TDDFT) associated with PCM is used for the electron structure and the spectroscopy calculations of the test molecules in homogeneous solvents. With the charge distribution of the ground and excited states, the position- and orientation-dependencies of the solvation energy and the spectrum have been investigated for the interfacial systems, taking the electrostatic interaction, the cavitation energy, and the dispersion-repulsion interaction into account. The cavitation energy is paid particular attention, since the interface portion cut off by the occupation of the interfacial molecule contributes an extra part to the stabilization for the interfacial system. The embedding depth, the favorable orientational angle, and the spectral shift for the interfacial molecule have been investigated in detail. From the solvation energy calculations, an explanation has been given on why the interfacial molecule, even if symmetrical in structure, tends to take a tilting manner, rather than perpendicular to the interface.     

How does halogen bonding behave in solution? A theoretical study using implicit solvation model

A systematic study of halogen bonding interactions in gas phase and in solution was carried out by means of quantum chemical DFT/B3LYP method. Three solvents with different polarities (chloroform, acetone, and water) were selected, and solvation effects were considered using the polarized continuum model (PCM). For charged halogen-bonded complexes, the strength of the interactions tends to significantly weaken in solution, with a concomitant elongation of intermolecular distances. For neutral systems, halogen bond distances are shown to shorten and the interaction energies change slightly. Computations also reveal that in the gas phase the binding affinities decrease in the order Cl(-) > Br(-) > I(-), while in solution the energy gaps of binding appear limited for the three halide anions. According to free energy results, many systems under investigation are stable in solution. Particularly, calculated free energies of formation of the complexes correlate well with halogen-bonding association constants determined experimentally. The differences of the effects of solvent upon halogen and hydrogen bonding were also elucidated. This study can establish fundamental characteristics of halogen bonding in media, which would be very helpful for applying this noncovalent interaction in medicinal chemistry and material design.     

Computational approach to nuclear magnetic resonance in 1-Alkyl-3-methylimidazolium ionic liquids

A quantum-chemical computational approach to accurately predict the nuclear magnetic resonance (NMR) properties of 1-alkyl-3-methylimidazolium ionic liquids has been performed by the gauge-including atomic orbitals method at the B3LYP/6-31++G** level using different simulated ionic liquid environments. The first molecular model chosen to describe the ionic liquid system includes the gas-phase optimized structures of ion pairs and separated ions of a series of imidazolium salts containing methyl, butyl, and octyl substituents and PF6-, BF4-, and Br- anions. In addition, a continuum polarizable model of solvation has been applied to predict the effects of the medium polarity on the molecular properties of 1,3-dimethylimidazolium hexafluorophosphate (MmimPF6). Furthermore, the specific acidic and basic solute-solvent interactions have been simulated by a discrete solvation model based on molecular clusters formed by MmimPF6 species and a discrete number of water molecules. The computational prediction of the NMR spectra allows a consistent interpretation of the dispersed experimental evidence in the literature. The following are main contributions of this work: (a) Theoretical results state the presence of a chemical equilibrium between ion-pair aggregates and solvent-separated counterions of 1-alkyl-3-methylimidazolium salts which is tuned by the solvent environment; thus, strong specific (acidic and basic) and nonspecific (polarity and polarizability) solvent interactions are predicted favoring the dissociated ionic species. (b) The calculated 1H and 13C NMR properties of these ionic liquids are revealed as highly dependent on the nature of solute-solvent interactions. Thus, the chemical shift of the hydrogen atom in position two of the imidazolium ring is deviated to high values by the specific interactions with water molecules, whereas nonspecific interaction with water (as a solvent) affects, in the opposite direction, this 1H NMR parameter. (c) Last, current calculations support the presence of hydrogen bonding between counterions, suggesting the importance of this interaction in the properties of the solvent in the 1-alkyl-3-methylimidazolium ionic liquids.     

Theoretical study of metal-ligand interaction in Sm(III), Eu(III), and Tb(III) complexes of coumarin-3-carboxylic acid in the gas phase and solution

The interaction of lanthanide(III) cations (Ln(III) = Sm(III), Eu(III), and Tb(III)) with the deprotonated form of the coumarin-3-carboxylic acid (cca-) has been investigated by density functional theory (DFT/B3LYP) and confirmed by reference MP2 and CCSD(T) computations. Solvent effects on the geometries and stabilities of the Ln(III) complexes were computed using a combination of water clusters and a continuum solvation model. The following two series of systems were considered: (i) Ln(cca)2+, Ln(cca)2+, Ln(cca)3 and (ii) Ln(cca)(H2O)2Cl2, Ln(cca)2(H2O)2Cl, Ln(cca)3. The strength and character of the Ln(III)-cca- bidentate bonding were characterized by calculated Ln-O bond lengths, binding energies, ligand deformation energies, energy partitioning analysis, sigma-donation contributions, and natural population analyses. The energy decomposition calculations predicted predominant electrostatic interaction terms to the Ln-cca bonding (ionic character) and showed variations of the orbital interaction term (covalent contributions) for the Ln-cca complexes studied. Electron distribution analysis suggested that the covalent contribution comes mainly from the interaction with the carboxylate moiety of cca-.     

Role of solvation in the energy stabilization inside the hydrophobic core of the protein rubredoxin

There are many forces that contribute to the stability of a protein; among these are dispersion interactions, hydrogen bonding, and solvation effects. In a recent work, Vondrasek et al. estimated the in vacuo stabilization energy of the hydrophobic core of the protein rubredoxin using high level ab initio methods (Vondrasek, J.; et al. J. Am. Chem. Soc. 2005, 127, 2615). In this work, we evaluate the effects of solvation on the stability of the hydrophobic core of this protein. Solvation calculations are made using the polarizable continuum method at the MP2/aug-cc-pVDZ level of theory. It is found that, in a protein-like environment (mimicked by a continuum solvent with a dielectric constant of approximately 4), the stability of rubredoxin's hydrophobic core is decreased by 40-50%. We also observed that the stabilization energy of the hydrophobic core is only slightly lower in a protein-like medium than in an aqueous one (DeltaGether-DeltaGwater approximately 1.0-3.5 kcal/mol).     

Application of the frozen atom approximation to the GB/SA continuum model for solvation free energy

The generalized Born/surface area (GB/SA) continuum model for solvation free energy is a fast and accurate alternative to using discrete water molecules in molecular simulations of solvated systems. However, computational studies of large solvated molecular systems such as enzyme-ligand complexes can still be computationally expensive even with continuum solvation methods simply because of the large number of atoms in the solute molecules. Because in such systems often only a relatively small portion of the system such as the ligand binding site is under study, it becomes less attractive to calculate energies and derivatives for all atoms in the system. To curtail computation while still maintaining high energetic accuracy, atoms distant from the site of interest are often frozen; that is, their coordinates are made invariant. Such frozen atoms do not require energetic and derivative updates during the course of a simulation. Herein we describe methodology and results for applying the frozen atom approach to both the generalized Born (GB) and the solvent accessible surface area (SASA) parts of the GB/SA continuum model for solvation free energy. For strictly pairwise energetic terms, such as the Coulombic and van-der-Waals energies, contributions from pairs of frozen atoms can be ignored. This leaves energetic differences unaffected for conformations that vary only in the positions of nonfrozen atoms. Due to the nonlocal nature of the GB analytical form, however, excluding such pairs from a GB calculation leads to unacceptable inaccuracies. To apply a frozen-atom scheme to GB calculations, a buffer region within the frozen-atom zone is generated based on a user-definable cutoff distance from the nonfrozen atoms. Certain pairwise interactions between frozen atoms in the buffer region are retained in the GB computation. This allows high accuracy in conformational GB comparisons to be maintained while achieving significant savings in computational time compared to the full (nonfrozen) calculation. A similar approach for using a buffer region of frozen atoms is taken for the SASA calculation. The SASA calculation is local in nature, and thus exact SASA energies are maintained. With a buffer region of 8 A for the frozen-atom cases, excellent agreement in differences in energies for three different conformations of cytochrome P450 with a bound camphor ligand are obtained with respect to the nonfrozen cases. For various minimization protocols, simulations run 2 to 10.5 times faster and memory usage is reduced by a factor of 1.5 to 5. Application of the frozen atom method for GB/SA calculations thus can render computationally tractable biologically and medically important simulations such as those used to study ligand-receptor binding conformations and energies in a solvated environment.     

有机π共轭配体溶剂化效应与分子间相互作用的理论研究

采用可极化的连续介质模型(PCM), 运用密度泛函理论(DFT), 在B3LYP/6-31+G^**水平下研究了溶剂极性对有机π共轭配体N,N'-Bis-(3-pyridyl)ethylene-bis-urea(BPEBU)中syn-anti构象的分子几何和电子结构的影响, 并借助分子动力学模拟的方法, 采用明确溶剂模型研究了溶质-溶剂分子间的相互作用. 密度泛函理论计算结果表明, 随着溶剂极性的增强, BPEBU中尿素基上的CO键和N-H键以及吡啶环上的C-N键被明显极化, 使羰基氧原子和吡啶氮原子的电负性明显增强, 尿素基的N-H键上氢原子的正电荷也显著增加. 分子动力学模拟统计的结果表明, 在极性较强的乙醇溶液中, 有明确的O…H-O, N…H-O和N-H…O等3种氢键作用存在, 而在丙酮溶液中, 只有N…H-O一种氢键作用存在, 而且与乙醇溶液中的N…H-O作用相比要弱些. 另外, 采用密度泛函理论方法结合连续/明确的混合溶剂模型, 优化得到了溶质-溶剂三聚体的超分子簇结构, 与分子动力学模拟的第一溶剂层中的超分子结构相比, 两者定性一致.     

Shells theory of solvation and the long-range Born correction

A new solvation model, named shells theory of solvation, is proposed. In this approach, the solvent is divided in two regions, the S₁ shell, close to the solute and describing specific solute–solvent interactions, and the S₂ shell, representing the remain solvent and accounting for the long-range interaction contribution. A simple theoretical equation can be derived which allows the computation of the solvation free energy using two-point thermodynamic integration and configurations generated from molecular dynamics simulation. The discrete/continuum version of this theory provides rigorous theoretical foundations for the popular long-range Born correction and presents a new reliable expression for including this contribution. Further, it converges to the full discrete representation of the solvent when the number of solvent molecules goes to infinity. The method can be easily applied when the solute–solvent interaction (S₁ shell) is treated by full quantum mechanics, while the S₂ shell is described by a dielectric continuum solvation method. A simple test of the theory was done for solvation of fluoride ion in benzene solution. The S₁ shell was composed of the fluoride ion plus 32 benzene molecules, and the interaction with the S₂ shell was calculated at Hartree–Fock level with the MINI basis set and using the polarizable continuum model.     

Solvation of excess electrons in LiF ionic pair matrix: evidence for a solvated dielectron from ab initio molecular dynamics simulations and calculations

Ab initio molecular dynamics simulations and first-principles calculations reveal the existence of a solvated dielectron species, (2e)s, in an LiF ionic matrix. The nature of the solvation mechanism and the stability of the species was explored. In addition to electrostatic interactions, a hole-orbital coupling among solvent molecules may significantly enhance the stability of the solvated electrons and govern the extent of electron solvation. This hole-orbital coupling is different from either an electrostatic coupling or conventional chemical bonding, and it may be described as a transition between them.