An application of coupled reference interaction site model/molecular dynamics to the conformational analysis of the alanine dipeptide

We present an application of our recently proposed coupled reference interaction site model (RISM) molecular dynamics (MD) solvation free energy methodology [Freedman and Truong, Chem. Phys. Lett. 381, 362 (2003); J. Chem. Phys. 121, 2187 (2004)] to study the conformational stability of alanine dipeptide in aqueous solution. In this methodology, radial distribution functions obtained from a single MD simulation are substituted into a RISM expression for solvation free energy. Consequently, iterative solution of the RISM equation is not needed. The relative solvation free energies of seven different conformations of the alanine dipeptide in aqueous solution are calculated. Results from the coupled RISM/MD methodology are in good agreement with those from earlier simulations using the accurate free energy perturbation approach, showing that the alphaR conformation is most stabilized by solution. This study establishes a framework for applying this coupled RISM/MD method to larger biological systems.     

Improving the performance of the coupled reference interaction site model-hyper-netted chain (RISM-HNC)/simulation method for free energy of solvation

The coupled reference interaction site model-hyper-netted chain (RISM-HNC)/ simulation methodology determines solvation free energies as a function of the set of all radial distribution functions of solvent atoms about atomic solute sites. These functions are determined from molecular dynamics (MD) or Monte Carlo (MC) simulations rather than from solving the RISM and HNC equations iteratively. Previous applications of the method showed that it can predict relative free energies of solvation for small solutes accurately. However, the errors scale with the system size. In this study, we propose the use of the hard-sphere free energy as the reference and a linear response approximation to improve the performance, i.e., accuracy and robustness, of the method, particularly removing the size dependency of the error. The details of the new formalism are presented. To validate the proposed formalism, solvation free energies of N-methylacetamide and methylamine are computed using the new RISM-HNC-based expressions in addition to a linear response expression, which are compared to previous thermodynamic integration and thermodynamic perturbation results performed with the same force field. Additionally, free energies of solvation for cyclohexane, pyridine, benzene and derivatives, and other small organic molecules are calculated and compared to experimental values.     

Photo absorption of ‐coumaric acid in aqueous solution: RISM‐SCF‐SEDD theory approach

Photo absorption properties of p‐coumaric acid, the chromophore of photoactive yellow protein, in aqueous solution were investigated by means of reference interaction site model self‐consistent field with spatial electron density distribution (RISM‐SCF‐SEDD) method. RISM‐SCF‐SEDD is a combination methodology of electronic structure theory and statistical mechanics for molecular liquids. Here, time‐dependent density functional theory was coupled with RISM equation to study the electronic structure of p‐coumaric acid in aqueous system. Excitation energies of the chromophore in its neutral, two monoanionic and dianionic forms were computed to elucidate the effect of the deprotonation and solvation on the spectroscopic properties. We found that solvation strongly affects the excitation character of the chromophore, especially for phenolate anion and dianion. The free energy difference among the four protonation states is also discussed. © 2017 Wiley Periodicals, Inc.     

Efficient treatment of solvation shells in 3D molecular theory of solvation

We developed a technique to decrease memory requirements when solving the integral equations of three‐dimensional (3D) molecular theory of solvation, a.k.a. 3D reference interaction site model (3D‐RISM), using the modified direct inversion in the iterative subspace (MDIIS) numerical method of generalized minimal residual type. The latter provides robust convergence, in particular, for charged systems and electrolyte solutions with strong associative effects for which damped iterations do not converge. The MDIIS solver (typically, with 2 × 10 iterative vectors of argument and residual for fast convergence) treats the solute excluded volume (core), while handling the solvation shells in the 3D box with two vectors coupled with MDIIS iteratively and incorporating the electrostatic asymptotics outside the box analytically. For solvated systems from small to large macromolecules and solid–liquid interfaces, this results in 6‐ to 16‐fold memory reduction and corresponding CPU load decrease in MDIIS. We illustrated the new technique on solvated systems of chemical and biomolecular relevance with different dimensionality, both in ambient water and aqueous electrolyte solution, by solving the 3D‐RISM equations with the Kovalenko–Hirata (KH) closure, and the hypernetted chain (HNC) closure where convergent. This core–shell‐asymptotics technique coupling MDIIS for the excluded volume core with iteration of the solvation shells converges as efficiently as MDIIS for the whole 3D box and yields the solvation structure and thermodynamics without loss of accuracy. Although being of benefit for solutes of any size, this memory reduction becomes critical in 3D‐RISM calculations for large solvated systems, such as macromolecules in solution with ions, ligands, and other cofactors. © 2012 Wiley Periodicals, Inc.     

A coupled reference interaction site model/molecular dynamics study of the potential of mean force curve of the SN2 Cl- + CH3Cl reaction in water

An application of the coupled reference interaction site model (RISM)/simulation methodology to the calculation of the potential of mean force (PMF) curve in aqueous solution for the identity nucleophilic substitution reaction Cl(-) + CH(3)Cl is performed. The free energy of activation is calculated to be 27.1 kcal/mol which compares very well with the experimentally determined barrier height of 26.6 kcal/mol. Furthermore, the calculated PMF is almost superimposed with that previously calculated using the computationally rigorous Monte Carlo with importance sampling method (Chandrasekhar, J.; Smith, S. F.; Jorgensen, W. L. J. Am. Chem. Soc. 1985, 107, 154). Using the calculated PMF, a crude estimate of the solvated kinetic transmission coefficient also compares well with that of previous more accurate simulations. These results indicate that the coupled RISM/simulation method provides a cost-effective methodology for studying reactions in solution.     

Theoretical study on aquation reaction of cis-platin complex: RISM-SCF-SEDD, a hybrid approach of accurate quantum chemical method and statistical mechanics

The ligand exchange process of cis-platin in aqueous solution was studied using RISM-SCF-SEDD (reference interaction site model-self-consistent field with spatial electron density distribution) method, a hybrid approach of quantum chemistry and statistical mechanics. The analytical nature of RISM theory enables us to compute accurate reaction free energy in aqueous solution based on CCSD(T), together with the microscopic solvation structure around the complex. We found that the solvation effect is indispensable to promote the dissociation of the chloride anion from the complex.     

Hydration of ionic species studied by the reference interaction site model with a repulsive bridge correction

We have tested the reference interaction site model (RISM) for the case of the hypernetted chain (HNC) and the partially linearized hypernetted chain (PLHNC) closures improved by a repulsive bridge correction (RBC) for ionic hydrated species. We have analyzed the efficiency of the RISM/HNC+RBC and RISM/PLHNC+RBC techniques for decomposition of the electrostatic and the nonpolar hydration energies on the energetic and the enthalpic parts for polyatomic ions when the repulsive bridge correction is treated as a thermodynamic perturbation, and investigate the repulsive bridge effect on the electrostatic potential induced by solvent on solute atoms. For a number of univalent and bivalent atomic ions, molecular cations, and anions, the method provides hydration energies deviating only by several percents from the experimental data. In most cases, the enthalpic contributions to the free energies are also close to the experimental results. The above models are able to satisfactory predict the hydration energies as well as the electrostatic potential around the ionic species. For univalent atomic ions, they also provide qualitative estimates of the Samoilov activation energies.     

Coupled semiempirical quantum mechanics and molecular mechanics (QM/MM) calculations on the aqueous solvation free energies of ionized molecules

The aqueous solvation free energies of ionized molecules were computed using a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1, MNDO, and PM3 semiempirical molecular orbital methods for the solute molecule and the TIP3P molecular mechanics model for liquid water. The present work is an extension of our model for neutral solutes where we assumed that the total free energy is the sum of components derived from the electrostatic/polarization terms in the Hamiltonian plus an empirical “nonpolar” term. The electrostatic/polarization contributions to the solvation free energies were computed using molecular dynamics (MD) simulation and thermodynamic integration techniques, while the nonpolar contributions were taken from the literature. The contribution to the electrostatic/polarization component of the free energy due to nonbonded interactions outside the cutoff radii used in the MD simulations was approximated by a Born solvation term. The experimental free energies were reproduced satisfactorily using variational parameters from the vdW terms as in the original model, in addition to a parameter from the one-electron integral terms. The new one-electron parameter was required to account for the short-range effects of overlapping atomic charge densities. The radial distribution functions obtained from the MD simulations showed the expected H-bonded structures between the ionized solute molecule and solvent molecules. We also obtained satisfactory results by neglecting both the empirical nonpolar term and the electronic polarization of the solute, i.e., by implementing a nonpolarization model. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1028–1038, 1999     

Combination of molecular dynamics method and 3D-RISM theory for conformational sampling of large flexible molecules in solution

We have developed an algorithm for sampling the conformational space of large flexible molecules in solution, which combines the molecular dynamics (MD) method and the three-dimensional reference interaction site model (3D-RISM) theory. The solvent-induced force acting on solute atoms was evaluated as the gradient of the solvation free energy with respect to the solute-atom coordinates. To enhance the computation speed, we have applied a multiple timestep algorithm based on the RESPA (Reversible System Propagator Algorithm) to the combined MD/3D-RISM method. By virtue of the algorithm, one can choose a longer timestep for renewing the solvent-induced force compared with that of the conformational update. To illustrate the present MD/3D-RISM simulation, we applied the method to a model of acetylacetone in aqueous solution. The multiple timestep algorithm succeeded in enhancing the computation speed by 3.4 times for this model case. Acetylacetone possesses an intramolecular hydrogen-bonding capability between the hydroxyl group and the carbonyl oxygen atom, and the molecule is significantly stabilized due to this hydrogen bond, especially in gas phase. The intramolecular hydrogen bond was kept intact during almost entire course of the MD simulation in gas phase, while in the aqueous solutions the bond is disrupted in a significant number of conformations. This result qualitatively agrees with the behavior on a free energy barrier lying upon the process for rotating a torsional degree of freedom of the hydroxyl group, where it is significantly reduced in aqueous solution by a cancellation between the electrostatic interaction and the solvation free energy.     

Solvent effects in four-component relativistic electronic structure theory based on the reference interaction-site model

A combined method of the Dirac–Hartree–Fock (DHF) method and the reference interaction-site model (RISM) theory is reported; this is the initial implementation of the coupling of the four-component relativistic electronic structure theory and an integral equation theory of molecular liquids. In the method, the DHF and RISM equations are solved self-consistently, and therefore the electronic structure of the solute, including relativistic effects, and the solvation structure are determined simultaneously. The formulation is constructed based on the variational principle with respect to the Helmholtz energy, and analytic free energy gradients are also derived using the variational property. The method is applied to the iodine ion (I⁻), methyl iodide (CH₃I), and hydrogen chalcogenide (H₂X, where X = O–Po) in aqueous solutions, and the electronic structures of the solutes, as well as the solvation free energies and their component analysis, solvent distributions, and solute–solvent interactions, are discussed.     

New approach to free energy of solvation applying continuum models to molecular dynamics simulation

A new approach to the calculation of the free energy of solvation from trajectories obtained by molecular dynamics simulation is presented. The free energy of solvation is computed as the sum of three contributions originated at the cavitation of the solute by the solvent, the solute-solvent nonpolar (repulsion and dispersion) interactions, and the electrostatic solvation of the solute. The electrostatic term is calculated based on ideas developed for the broadly used continuum models, the cavitational contribution from the excluded volume by the Claverie-Pierotti model, and the Van der Waals term directly from the molecular dynamics simulation. The proposed model is tested for diluted aqueous solutions of simple molecules containing a variety of chemically important functions: methanol, methylamine, water, methanethiol, and dichloromethane. These solutions were treated by molecular dynamics simulations using SPC/E water and the OPLS force field for the organic molecules. Obtained free energies of solvation are in very good agreement with experimental data.     

Comparison of implicit solvent models for the simulation of protein-surface interactions

Empirical force field-based molecular simulations can provide valuable atomistic-level insights into protein-surface interactions in aqueous solution. While the implicit treatment of solvation effects is desired as a means of improving simulation efficiency, existing implicit solvent models were primarily developed for the simulation of peptide or protein behavior in solution alone, and thus may not be appropriate for protein interactions with synthetic material surfaces. The objective of this research was to calculate the change in free energy as a function of surface-separation distance for peptide-surface interactions using different empirical force field-based implicit solvation models (ACE, ASP, EEF1, and RDIE with the CHARMM 19 force field), and to compare these results with the same calculations conducted using density functional theory (DFT) combined with the self-consistent reaction field (SCRF) implicit solvation model. These comparisons show that distinctly different types of behavior are predicted with each implicit solvation method, with ACE providing the best overall agreement with DFT/SCRF calculations. These results also identify areas where ACE is in need of improvement for this application and provide a basis for subsequent parameter refinement.     

A solvation‐free‐energy functional: A reference‐modified density functional formulation

The three‐dimensional reference interaction site model (3D‐RISM) theory, which is one of the most applicable integral equation theories for molecular liquids, overestimates the absolute values of solvation‐free‐energy (SFE) for large solute molecules in water. To improve the free‐energy density functional for the SFE of solute molecules, we propose a reference‐modified density functional theory (RMDFT) that is a general theoretical approach to construct the free‐energy density functional systematically. In the RMDFT formulation, hard‐sphere (HS) fluids are introduced as the reference system instead of an ideal polyatomic molecular gas, which has been regarded as the appropriate reference system of the interaction‐site‐model density functional theory for polyatomic molecular fluids. We show that using RMDFT with a reference HS system can significantly improve the absolute values of the SFE for a set of neutral amino acid side‐chain analogues as well as for 504 small organic molecules. © 2015 Wiley Periodicals, Inc.     

The oxidation of tyrosine and tryptophan studied by a molecular dynamics normal hydrogen electrode

The thermochemical constants for the oxidation of tyrosine and tryptophan through proton coupled electron transfer in aqueous solution have been computed applying a recently developed density functional theory (DFT) based molecular dynamics method for reversible elimination of protons and electrons. This method enables us to estimate the solvation free energy of a proton (H(+)) in a periodic model system from the free energy for the deprotonation of an aqueous hydronium ion (H(3)O(+)). Using the computed solvation free energy of H(+) as reference, the deprotonation and oxidation free energies of an aqueous species can be converted to pK(a) and normal hydrogen electrode (NHE) potentials. This conversion requires certain thermochemical corrections which were first presented in a similar study of the oxidation of hydrobenzoquinone [J. Cheng, M. Sulpizi, and M. Sprik, J. Chem. Phys. 131, 154504 (2009)]. Taking a different view of the thermodynamic status of the hydronium ion, these thermochemical corrections are revised in the present work. The key difference with the previous scheme is that the hydronium is now treated as an intermediate in the transfer of the proton from solution to the gas-phase. The accuracy of the method is assessed by a detailed comparison of the computed pK(a), NHE potentials and dehydrogenation free energies to experiment. As a further application of the technique, we have analyzed the role of the solvent in the oxidation of tyrosine by the tryptophan radical. The free energy change computed for this hydrogen atom transfer reaction is very similar to the gas-phase value, in agreement with experiment. The molecular dynamics results however, show that the minimal solvent effect on the reaction free energy is accompanied by a significant reorganization of the solvent.     

烷基芳基磺酸盐的分子动力学模拟与自由能微扰计算

为研究不同结构的表面活性剂分子在溶液中胶束化能力的差异, 采用分子动力学方法模拟三种烷基芳基磺酸盐在真空和水溶液环境下的结构与相互作用. 利用自由能微扰(FEP)方法计算了水合自由能, 发现与用传统热力学表面张力法测定自制的烷基芳基磺酸盐结果一致. 研究表明: 烷基芳基磺酸盐在水溶液中的胶束化过程是自发进行的, 随着分子结构中芳环向长烷基链中间位置移动, 胶束化能力和胶束稳定性均下降; 疏水基周围水分子的“冰山结构”会影响胶束的稳定性, 而水分子中氢键的生存周期是反映冰山结构变化的重要指标; 同时, 亲水基与水分子间形成氢键的数目会增强或减弱分子脱离胶束体的趋势, 从而影响胶束结构的稳定性.     

A novel quantum mechanical/molecular mechanical approach to the free energy calculation for isomerization of glycine in aqueous solution

The free energy change associated with the isomerization reaction of glycine in water solution has been studied by a hybrid quantum mechanical/molecular mechanical (QM/MM) approach combined with the theory of energy representation (QM/MM-ER) recently developed. The solvation free energies for both neutral and zwitterionic form of glycine have been determined by means of the QM/MM-ER simulation. The contributions of the electronic polarization and the fluctuation of the QM solute to the solvation free energy have been investigated. It has been found that the contribution of the density fluctuation of the zwitterionic solute is estimated as -4.2 kcal/mol in the total solvation free energy of -46.1 kcal/mol, while that of the neutral form is computed as -3.0 kcal/mol in the solvation free energy of -15.6 kcal/mol. The resultant free energy change associated with the isomerization of glycine in water has been obtained as -7.8 kcal/mol, in excellent agreement with the experimental data of -7.3 or -7.7 kcal/mol, implying the accuracy of the QM/MM-ER approach. The results have also been compared with those computed by other methodologies such as the polarizable continuum model and the classical molecular simulation. The efficiency and advantage of the QM/MM-ER method has been discussed.     

Diabatic free energy curves and coordination fluctuations for the aqueous Ag+Ag2+ redox couple: a biased Born-Oppenheimer molecular dynamics investigation

Biased Born-Oppenheimer molecular dynamics simulations are performed to compute redox potential and free energy curves for the redox half reaction Ag(+)-->Ag(2+)+e(-) in aqueous solution. The potential energy surfaces of reactant and product state are linearly coupled and the system transferred from the reduced state to the oxidized state by variation of the coupling parameter from 0 to 1. The redox potential is obtained by thermodynamic integration of the average ionization energy of Ag(+). Diabatic free energy curves of reduced (R) and oxidized (O) states are obtained to good statistical accuracy by reweighting and combining the set of biased distributions of the ionization energy. The diabatic free energy curves of Ag(+) and Ag(2+) are parabolic over a wide range of the reaction coordinate in agreement with the linear response assumption that underlies Marcus theory. However, we observe deviations from parabolic behavior in the equilibrium region of Ag(+) and find different values for the reorganization free energy of R (1.4 eV) and O (0.9 eV). The computed reorganization free energy of Ag(2+) is in good agreement with the experimental estimate of 0.9-1.2 eV obtained from photoelectron spectroscopy. As suggested by our calculations, the moderate deviation from linear response behavior found for Ag(+) is likely related to the highly fluxional solvation shell of this ion, which exhibits water exchange reactions on the picosecond time scale of the present molecular dynamics simulation.     

Ab initio theoretical study of temperature and density dependence of molecular and thermodynamic properties of water in the entire fluid region: autoionization processes

The temperature and density dependence of the molecular and thermodynamic properties of water is investigated theoretically by means of the ab initio electronic structure theory combined with the reference interaction site model method, so-called RISM-SCF. We consider the autoionization process (H2O + H2O right harpoon over left harpoon H3O+ + OH-) by regarding H2O, H3O+, and OH- as "solute" molecules in an aqueous solution and evaluate molecular geometry, electronic structure, solvation structure, and the ionic product of water (pKw) of these species as functions of thermodynamic conditions. In our previous paper, we calculated these properties by using essentially the same method in a wide range of density values (0.6-1.4 g/cm3). However, the calculation was limited at rather higher density (>0.6 g/cm3) due to the difficulty of convergence, which is inherent to the hypernetted-chain (HNC) closure. The problem is overcome in this study by employing the Kovalenko-Hirata (KH) closure which hybridizes the HNC and the mean-spherical approximation (MSA). Here, we present the results for the thermodynamic range of densities from 0.025 to 1.0 g/cm3 and for temperatures from 300 to 800 K including the supercritical point.     

Free-energy analysis of the molecular binding into lipid membrane with the method of energy representation

A statistical-mechanical treatment of the molecular binding into lipid membrane is presented in combination with molecular simulation. The membrane solution is viewed as an inhomogeneous, mixed solvent system, and the free energy of solvation of a solute in membrane is computed with a realistic set of potential functions by the method of energy representation. Carbon monoxide, carbon dioxide, benzene, and ethylbenzene are adopted as model solutes to analyze the binding into 1,2-dimyristoyl-sn-glycero-3-phosphatidylcholine (DMPC) membrane. It is shown that the membrane inside is more favorable than bulk water and that the solute distribution is diffuse throughout the membrane inside. The membrane-water partition coefficient is then constructed with the help of the Kirkwood-Buff theory from the solvation free energy obtained separately in the hydrophobic, glycerol, headgroup, and aqueous regions. To discuss the role of repulsive and attractive interactions, the solvation free energy is partitioned into the DMPC and water contributions and the effect of water to stabilize the benzene and ethylbenzene solutes within the membrane is pointed out.