Thermodynamic Integration Method Research Articles

Binding of Phosphate Species to Ca2+ and Mg2+ in Aqueous Solution.

Phosphate derivatives and their interaction with metal cations are involved in many important biological phenomena, so an accurate characterization of the phosphate-metal interaction is necessary to properly understand the role of phosphate-metal contacts in mediating biological function. Herein, we improved the standard 12-6 Lennard-Jones (LJ) potential via the usage of the 12-6-4 LJ model, which incorporates ion-induced dipole interactions. Via parameter scanning, we fine-tuned the 12-6-4 LJ polarizability values to obtain accurate absolute binding free energies for the phosphate anions H2PO4-, HPO42-, PO43- coordinating with Ca2+ and Mg2+. First, we modified the phosphate 12-6-4 LJ parameters to reproduce the solvation free energies of the series of phosphate anions using the thermodynamic integration (TI) method. Then, using the potential mean force (PMF) method, the polarizability of the metal-phosphate interaction was obtained. We show that the free energy profiles of phosphate ions coordinated to Ca2+ and Mg2+ generally show similar trends at longer metal-phosphate distances, while the absolute binding energy values increased with deprotonation. The resulting parameters demonstrate the flexibility of the 12-6-4 LJ-type nonbonded model and its usefulness in accurately describing cation-anion interactions.

Zirconium phase diagram from ab initio molecular dynamics

The paper presents analysis of precision of a phonon-based free energy calculation technique (Temperature Dependent Effective Potential method) in comparison with the thermodynamic integration method within classical molecular dynamics. It is shown that the use of the second order force constants is insufficient for accurate calculations of free energies at elevated temperatures due to noticeable anharmonic contributions to lattice vibrations. On the other hand, reliable calculations of the crystal phase equilibrium lines requires free energy accuracy better than 1–3 meV per atom. Based upon the obtained results the Temperature Dependent Effective Potential technique is applied to the detailed study of the solid-state Zr phase diagram from ab initio molecular dynamics that may be achieved if force constants up to the fourth order are used. The calculated zirconium phase diagram is in good agreement with experiment. The α−β−ω triple point was determined at P≈6.74 GPa, T≈837 K.

Atomistic investigation of phase transformations in NiTiCu shape memory alloys

NiTiCu shape memory alloys have attracted much attention for their ability to reduce hysteresis, which is required fast responses in actuator and sensor applications, but the correlation between the properties and specific process of phase transformation is unclear. In this work, an atomic-level understanding of the detailed phase transformation behavior of NiTiCu shape-memory alloys is achieved based on a newly developed interatomic potential for the Ni–Ti–Cu ternary system. Free energy calculations based on the thermodynamic integration method faithfully reproduce reported experimental trends, such as the compositional dependency of different phase transformation paths. Further, large-scale molecular dynamics simulations performed on nanocrystalline alloys provide visual evidence of the detailed phase transformation process involving different martensitic transformation paths. Our results suggest that although both the composition and grain size can be adjusted to realize slim hysteresis, the effect of increasing the Cu amount is only significant if changes in the phase transformation modes are involved at relatively low Cu contents. When the material exhibits only a single phase transformation mode at relatively high Cu contents, the additional Cu is not highly effective in reducing the hysteresis.

Nonbonded Force Field Parameters from MBIS Partitioning of the Molecular Electron Density Improve Binding Affinity Predictions of the T4-Lysozyme Double Mutant.

The use of computer simulation for binding affinity prediction is growing in drug discovery. However, its wider use is constrained by the accuracy of the free energy calculations. The key sources of error are the force fields used to depict molecular interactions and insufficient sampling of the configurational space. To improve the quality of the force field, we developed a Python-based computational workflow. The workflow described here uses the minimal basis iterative stockholder (MBIS) method to determine atomic charges and Lennard-Jones parameters from the polarized molecular density. This is done by performing electronic structure calculations on various configurations of the ligand when it is both bound and unbound. In addition, we validated a simulation procedure that accounts for the protein and ligand degrees of freedom to precisely calculate binding free energies. This was achieved by comparing the self-adjusted mixture sampling and nonequilibrium thermodynamic integration methods using various protein and ligand conformations. The accuracy of predicting binding affinity is improved by using MBIS-derived force field parameters and a validated simulation procedure. This improvement surpasses the chemical precision for the eight aromatic ligands, reaching a root-mean-square error of 0.7 kcal/mol.

Acceleration of Solvation Free Energy Calculation via Thermodynamic Integration Coupled with Gaussian Process Regression and Improved Gelman-Rubin Convergence Diagnostics.

The determination of the solvation free energy of ions and molecules holds profound importance across a spectrum of applications spanning chemistry, biology, energy storage, and the environment. Molecular dynamics simulations are powerful tools for computing this critical parameter. Nevertheless, the accurate and efficient calculation of the solvation free energy becomes a formidable endeavor when dealing with complex systems characterized by potent Coulombic interactions and sluggish ion dynamics and, consequently, slow transition across various metastable states. In the present study, we expose limitations stemming from the conventional calculation of the statistical inefficiency g in the thermodynamic integration method, a factor that can hinder the determination of convergence of the solvation free energy and its associated uncertainty. Instead, we propose a robust scheme based on Gelman-Rubin convergence diagnostics. We leverage this improved estimation of uncertainties to introduce an innovative accelerated thermodynamic integration method based on the Gaussian Process regression. This methodology is applied to the calculation of the solvation free energy of trivalent rare-earth elements immersed in ionic liquids, a scenario in which the aforementioned challenges render standard approaches ineffective. The proposed method proves to be effective in computing solvation free energy in situations where traditional thermodynamic integration methods fall short.

Calculation of Melting Temperature Using Nonequilibrium Thermodynamic Integration Methods

AbstractMelting temperature is a fundamental material property and is defined as the temperature at which the solid and liquid phases have the same free energy. However, there is no systematic study employing atomic simulations to calculate melting temperature using this definition. Here, molecular dynamics simulations and nonequilibrium thermodynamic integration methods are combined to calculate the melting temperature of Al and Cu. Results show that to accurately obtain the melting temperature, the model size should be considered carefully because the free energies of both solid and liquid phases are inversely proportional to the model size, causing a model size dependence on the calculated melting temperature. In addition, the melting temperature for various (semi‐) empirical potential models for Al and Cu is calculated and verified against experimental values to provide guidelines for the choice of potential models for simulation‐based problems relevant to the solid–liquid phase transformation.

Harmonic models and molecular dynamics simulations of isomorph behavior of Lennard-Jones fluids: Excess entropy and high temperature limiting behavior

Henchman’s approximate harmonic model of liquids is extended to predict the thermodynamic behavior along lines of constant excess entropy (“isomorphs”) in the liquid and supercritical fluid regimes of the Lennard-Jones (LJ) potential phase diagram. Simple analytic expressions based on harmonic cell models of fluids are derived for the isomorph lines, one accurate version of which only requires as input parameters the average repulsive and attractive parts of the potential energy per particle at a single reference state point on the isomorph. The new harmonic cell routes for generating the isomorph lines are compared with those predicted by the literature molecular dynamics (MD) methods, the small step MD method giving typically the best agreement over a wide density and temperature range. Four routes to calculate the excess entropy in the MD simulations are compared, which includes employing Henchman’s formulation, Widom’s particle insertion method, thermodynamic integration, and parameterized LJ equations of state. The thermodynamic integration method proves to be the most computationally efficient. The excess entropy is resolved into contributions from the repulsive and attractive parts of the potential. The repulsive and attractive components of the potential energy, excess Helmholtz free energy, and excess entropy along a fluid isomorph are predicted to vary as ∼T−1/2 in the high temperature limit by an extension of classical inverse power potential perturbation theory statistical mechanics, trends that are confirmed by the MD simulations.

Highly Accurate Densities and Isobaric and Isochoric Heat Capacities of Compressed Liquid Water Derived from New Speed of Sound Measurements

Comprehensive and accurate measurements of the speed of sound in liquid water are reported. The measurements were carried out by a double-path-length pulse-echo technique and cover the temperature range from 273.65 K to 368.15 K with pressures up to 100 MPa. The relative expanded (k=2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${ k}=2$$\\end{document}) uncertainties are 2.1 mK in temperature, 45 parts-per-million (ppm) in pressure, and between 40 ppm and 70 ppm in speed of sound. Furthermore, values for the density and specific isobaric and isochoric heat capacities were derived from the speed of sound data in the measured temperature range up to 100 MPa by the method of thermodynamic integration. Very accurate values for the derived properties were obtained by using density data of Takenaka and Masui (Metrologia 27:165–171, 1990) and isobaric heat capacity data of Osborne et al. (J Res Natl Bur Stand 23:197, 1939) at ambient pressure as initial values in combination with an accurate correlation of our speed of sound. The relative expanded (k=2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$${ k}=2$$\\end{document}) uncertainties of the derived properties are 2 ppm in density, 0.11% in isobaric heat capacity, and 0.12% in isochoric heat capacity. The experimental speeds of sound and derived properties are compared with experimental data of other authors from the literature, the IAPWS-95 (International Association of the Properties of Water and Steam) formulation for the thermodynamic properties of water, and a recent equation of state for supercooled water.

Atomistic simulations and machine learning of solute grain boundary segregation in Mg alloys at finite temperatures

Understanding solute segregation thermodynamics is the first step in investigating grain boundary (GB) properties, such as strong yttrium (Y) effects on grain growth and texture evolution in micro-scale polycrystalline magnesium (Mg) alloys. To estimate the average GB segregation behavior in low-solute-concentration Mg alloys (e.g., 2 at.% Y), a state-of-the-art spectral approach is applied based on a per-site segregation energy spectrum for Y solute atoms at zero K obtained from molecular statistics (MS) simulations of ∼104 GB sites in Mg symmetric tilt GBs (STGBs). Although selected MS simulation results are consistent with verification by density functional theory (DFT) calculations, estimates of average segregation tendency based on the zero-K energy spectrum deviate from experimental observations. To resolve this problem, thermodynamic integration (TI) methods based on molecular dynamics (MD) simulations are used to determine the per-site segregation free energies of Y at representative GB sites, which show contributions beyond harmonic approximations can be important for certain GB sites at high temperatures. A surrogate model of per-site segregation free energy is constructed from a small set of TI data points using stacking cross-validation regressors and physics-informed descriptors. This model is applied to predict the Y segregation free energy spectra for thousands of GB sites in Mg STGBs with uncertainty quantification. The average segregation tendency predicted by the spectral approach based on the free energy spectra agrees well (within the uncertainty range) with experimental observations for micro-scale polycrystalline Mg alloys at typical thermomechanical processing temperatures (500 ∼ 800 K), where thermodynamic equilibrium states are likely to be achieved due to fast diffusion.

Development of a model for granule-bound starch synthase I activity using free-energy calculations

Starch is a branched polymer of glucose with two components, both of which have (1 → 4)-α linear links and (1 → 6)-α branch points: amylopectin, of high molecular weight with many short branches, and amylose, of lower molecular weight and only a few long-chain branches. Granule-bound starch synthase I (GBSSI) is one of the main enzymes controlling amylose synthesis and chain-length distribution. As production of different GBSSI mutants is time-consuming and laborious, molecular dynamics (MD) simulations are used here to predict the binding of different GBSSI mutants to a representative amylose fragment. The simulations were atomistic, with explicit solvent and docking, a method successfully used to understand the binding of wild-type GBSSI to amylose fragments. The binding of GBSSI to G5 (a pentasaccharide amylose fragment) is combined with free-energy calculations employing a thermodynamic integration method to predict the effects of mutations on enzyme activity. Ten GBSSI mutants with different enzyme activities were analyzed to find the structural and energy changes among different single amino-acid mutants and their possible relationship to starch characteristics. Comparing the structural changes and the relative binding free energy of G5 to the wild type GBSSI and GBSSI mutants, it was found that mutants with negative binding energy (lower than −2.0 kcal/mol) are more likely to have higher enzyme activity and amylose content compared to the wild type. This theoretical paper used simulations and robust free energy calculations to interpret in planta data with potential predictions as to what mutants might be generated to give desired properties. This study can be used to help develop grains with improved functional properties.

Constrained Hybrid Monte Carlo Sampling Made Simple for Chemical Reaction Simulations.

Most electrochemical reactions should be studied under a grand canonical ensemble condition with a constant potential and/or a constant pH value. Free energy profiles provide key insights into understanding the reaction mechanisms. However, many molecular dynamics (MD)-based theoretical studies for electrochemical reactions did not employ an exact grand canonical ensemble sampling scheme for the free energy calculations, partially due to the issues of discontinuous trajectories induced by the particle-number variations during MD simulations. An alternative statistical sampling approach, the Monte Carlo (MC) method, is naturally appropriate for the open-system simulations if we focus on the thermodynamic properties. An advanced MC scheme, the hybrid Monte Carlo (HMC) method, which can efficiently sample the configurations of a system with large degrees of freedom, however, has limitations in the constrained-sampling applications. In this work, we propose an adjusted constrained HMC method to compute free energy profiles using the thermodynamic integration (TI) method. The key idea of the method for handling the constraint in TI is to integrate the reaction coordinate and sample the rest degrees of freedom by two types of MC schemes, the HMC scheme and the Metropolis algorithm with unbiased trials (M(RT)2-UB). We test the proposed method on three different systems involving two kinds of reaction coordinates, which are the distance between two particles and the difference of particles' distances, and compare the results to those generated by the constrained M(RT)2-UB method serving as benchmarks. We show that our proposed method has the advantages of high sampling efficiency and convenience of implementation, and the accuracy is justified as well. In addition, we show in the third test system that the proposed constrained HMC method can be combined with the path integral method to consider the nuclear quantum effects, indicating a broader application scenario of the sampling method reported in this work.

Anharmonicity in grain boundary energy for Al: Thermodynamic integration with artificial-neural-network potential

Anharmonicity in grain boundary (GB) energies for Al is elucidated by employing a thermodynamic integration method with a high-accuracy artificial-neural-network potential. Symmetric tilt GBs with the [001] and [11¯0] rotational axes are examined. It is found that anharmonic components tend to be small at low temperatures and increasingly contribute to reduction of GB free energies at elevated temperatures. The magnitude of anharmonicity is roughly correlated with excess volume at GBs. This is presumably because GB atoms tend to have longer bond lengths than the bulk atom, and hence the inclination of the potential energy surface at GBs is more moderate than in the bulk. On the other hand, the Σ13(510) GB takes non-negligible values of the anharmonic component even at intermediate temperatures. This GB contains a distinct atom around which the potential-energy curve significantly deviates from that of the harmonic approximation. Similar trends are also observed for other [001] GBs.

Entropy-driven phases at high coverage adsorption of straight rigid rods on three-dimensional cubic lattices.

Combining Monte Carlo simulations and thermodynamic integration method, we study the configurational entropy per site of straight rigid rods of length k (k-mers) adsorbed on three-dimensional (3D) simple cubic lattices. The process is monitored by following the dependence of the lattice coverage θ on the chemical potential μ (adsorption isotherm). Then, we perform the integration of μ(θ) over θ to calculate the configurational entropy per site of the adsorbed phase s(k,θ) as a function of the coverage. Based on the behavior of the function s(k,θ), different phase diagrams are obtained according to the k values: k≤4, disordered phase; k=5,6, disordered and layered-disordered phases; and k≥7, disordered, nematic and layered-disordered phases. In the limit of θ→1 (full coverage), the configurational entropy per site is determined for values of k ranging between 2 and 8. For k≥6, MC data coincide (within the statistical uncertainty) with recent analytical predictions [D. Dhar and R. Rajesh, Phys. Rev. E 103, 042130 (2021)2470-004510.1103/PhysRevE.103.042130] for very large rods. This finding represents the first numerical validation of the expression obtained by Dhar and Rajesh for d-dimensional lattices with d>2. In addition, for k≥5, the values of s(k,θ→1) for simple cubic lattices are coincident with those values reported in [P. M. Pasinetti et al., Phys. Rev. E 104, 054136 (2021)2470-004510.1103/PhysRevE.104.054136] for two-dimensional (2D) square lattices. This is consistent with the picture that at high densities and k≥5, the layered-disordered phase is formed on the lattice. Under these conditions, the system breaks to 2D layers, and the adsorbed phase becomes essentially 2D. The 2D behavior of the fully covered lattice reinforces the conjecture that the large-k behavior of entropy per site is superuniversal, and holds on d-dimensional hypercubical lattices for all d≥2.

Different metal coordination in sub- and super-critical fluids: Do molybdenum(IV) chloride complexes contribute to mass transfer in magmatic systems?

Molybdenum (Mo) is extracted predominantly from magmatic-hydrothermal ore deposits, and knowledge of the thermodynamic properties of aqueous Mo complexes over a wide range in pressures (P) and temperatures (T) is required to predict the transport and deposition of Mo in ore-forming environments. These properties are usually derived from experimental studies conducted over a limited P-T range, with uncertainties increasing as temperature and pressure increase, because of inherent technical difficulties and decreasing number of studies. There is a discrepancy between numerous experimental studies at T < 400 °C and a handful of studies at magmatic-hydrothermal conditions. The low T (<400 °C) experiments suggested that Mo(VI) chloride complexes are unlikely to play a significant role in nature, whereas some studies in magmatic-hydrothermal conditions showed the predominance of Mo(VI) chloride complexes. In an attempt to resolve these inconsistencies, we investigated Mo(VI)-Cl complexations from 350 to 750 °C (500 to 2000 bar) using ab initio molecular dynamics (MD) simulations. The simulations reveal that in high-temperature magmatic-hydrothermal fluids (750 °C, 2000 bar), a tetrahedral Mo(VI) chloride complex, MoO2(OH)Cl(aq), becomes stable; this coordination differs from the octahedral Mo(VI)-chloride complexes identified in highly acidic sub-critical fluids. The ab initio thermodynamic integration method was employed to calculate the equilibrium constants of the reaction.MoO2(OH)Cl(aq) + H2O = MoO2(OH)2(aq) + HCl(aq)Thermodynamic modelling of the speciation using the calculated equilibrium constants shows that at lower temperatures (350 °C, 500 bar), Mo(VI)-chloro complexes are only important at highly acidic conditions (pH < 1), in accordance with previous studies. However, the MD results show that at higher temperatures, Mo-chloro-complexes may be important in transporting Mo in mildly acidic conditions up to the quartz-muscovite-K-feldspar buffered pH. Hence, Mo(VI) chloro-complexes could play an important role in the formation of magmatic-hydrothermal Mo deposits (e.g., Climax-type Mo deposits). We conclude that the mode of metal transport can change dramatically at high temperatures as a result of coordination changes caused by changes in the solvation properties of water; and that quantitative MD is an important tool to support the interpretation of scant experimental results available for magmatic-hydrothermal systems.

Efficient Training of the Machine-Learning Interatomic Potential Based on an Artificial Neural Network for Estimating the Helmholtz Free Energy of Alkali Metals

Although the Helmholtz free energy can be accurately calculated using the thermodynamic integration (TI) method based on first-principles molecular dynamics (FPMD) simulations, the computational cost of this method is very high. By training the data obtained by FPMD simulations based on an artificial neural network (ANN), an interatomic potential with a very low computational cost can be constructed. To learn the FPMD data efficiently, it is important to elucidate the relation between the number of training data and accuracy of the free energy obtained by the TI method. We statistically evaluated the minimum number of training data by constructing multiple ANN potentials to accurately estimate the Helmholtz free energy of rubidium. The root mean square errors (RMSEs) of energy, atomic force, and pressure between the ANN and FPMD converged as the number of training data increased. The radial distribution function of the FPMD can be reproduced, even if the RMSEs do not converge. However, the number of training data for which the RMSEs converge is required to calculate the free energy with sufficient accuracy. This study demonstrated that the quality of an ANN potential can be estimated by statistically investigating the dependence of the RMSEs on the number of training data.

Temperature Transferable and Thermodynamically Consistent Coarse-Grained Model for Binary Polymer Systems

For the simulations of polymer systems, coarse-grained (CG) models are often developed to tackle the length and time scale limitations that are not feasible by all-atom molecular dynamic simulations. However, due to the necessary simplification or reduction in atomistic degrees of freedom, CG models usually have poor transferability over various temperatures, compositions, or thermodynamic states. In this work, the structure-based iterative Boltzmann inversion (IBI) method is further developed to obtain temperature transferable and thermodynamically consistent CG potentials for binary polymer systems. In particular, the conventional IBI procedure is first utilized to optimize single-component CG potentials from homopolymer systems while cross-interaction potentials from a random block copolymer system. Afterward, a thermodynamic integration method is applied to refine the cross-interaction potential between two constituent components in binary systems until the Flory–Huggins parameter (χ) between two components is matched with experimental value. We apply the above optimization procedure for the polystyrene-block-poly(methyl methacrylate) (PS-b-PMMA) as an example, which is widely studied experimentally, and the χ value can be easily obtained in the literature. To validate the transferability of the obtained CG potential, systems of both PS and PMMA homopolymers, random block copolymers, diblock copolymers, and PS/PMMA blends with various compositions are simulated. Many properties, including the mass density, radius distribution function, and bonded distribution functions, between CG beads generated from CG simulations match very well with those from atomistic simulations in the simulated temperature range, i.e., from 453 to 500 K. More importantly, phase morphologies of diblock copolymers and phase diagram of PS/PMMA binary blends over wide temperature and compositional ranges generated from CG simulations using our model are in good agreement with experimental results and predictions from self-consistent field theory.

Noble gas (He, Ne, and Ar) solubilities in high-pressure silicate melts calculated based on deep-potential modeling

Noble gas solubility in silicate melts has provided important insights into how these gases are incorporated into the Earth's interior and how they can be used as tracers of geochemical processes. However, previous studies have shown conflicting results regarding the noble gas solubility in silicate melts at high pressures. Some suggest that a drastic drop in solubility occurs above 5 GPa, while others indicate that maximum solubility is reached at this pressure but the solubility keeps unchanged when the pressure further increases. To clarify this discrepancy and gain a better understanding of the underlying mechanisms, we conducted a solubility investigation using the thermodynamic integration (TI) method and molecular dynamics (MD) simulations based on a machine-learning based force-field called deep potential (DP), which was trained by density functional theory (DFT) calculations. The effects of temperature, pressure and melt composition were investigated up to ∼30 GPa, 2500 K for the melt of model basalt, olivine and enstatite. Our results showed that the solubility of all noble gases (He, Ne and Ar) increased with temperature, while they increased with pressure until maximum solubility was reached at around 5 GPa. Beyond this pressure, solubility gradually decreased as pressure continued to increase. The structural and dynamic analysis of silicate melts revealed an abnormal increase in the self-diffusivity of oxygen atoms in silicate melts from ambient pressure to ∼5–8 GPa, which can be attributed to the formation of Si/Al-O5-6 in the melts upon compression. Interestingly, the upper limit pressure (5–8 GPa) for this abnormal increase in oxygen diffusivity was also the pressure at which noble gas solubilities reached their maximum. Our findings suggest that the maximum noble gas solubility at around 5 GPa may be a universal phenomenon in silicate melts. However, the specific values and changing trends of solubility depend on the specific melt compositions and noble gases. The complexity of noble gas solubility in silicate melts at high temperatures and pressures suggests that there is significant potential for improving the high-pressure noble gas partition coefficients between minerals, silicate melts, and even the Earth's core. Such improvements would extend our understanding of the behavior of noble gases in the Earth's interior and their role in tracing geodynamic processes.

Computational study on the encapsulation of glucosamine anomers by cucurbit[6]uril and cucurbit[8]uril in aqueous solution

Recently, we have investigated the cucurbit[7]uril (CB7) recognition of the α- and β-anomers of neutral, protonated, and acetylated forms of glucosamine in water. In the present work, we employed molecular dynamics (MD) and thermodynamic integration methods (TI) to investigate the recognition of these molecules by cucurbit[6]uril (CB6) and cucurbit[8]uril (CB8). MD revealed the formation of stable 1:1 inclusion complexes by all studied molecules with cucurbit[n]urils (CBn), and 2:1 complex by the α-anomer of the acetylated form of glucosamine with the large homologue CB8. CB6 forms roughly twice as many hydrogen bonds with the guest molecules as CB8. MM-PBSA results indicated that the electrostatic contribution to the binding free energy of each guest:CB complex was larger for CB6 than for CB8, and that CB6 and CB8 have lower affinity toward the different forms of glucosamine compared to CB7. Furthermore, TI was used to estimate the relative affinities of CB6 and CB8 toward the α- and β-anomers for each form of the studied glucosamine and compare with CB7.

Alchemical free energy to probe RNA/ligand interactions.

Alchemical free energy calculations are increasingly being applied in academic research as well as in the pharmaceutical industry. These methods originate in theories of Kirkwood and Zwanzig, where Kirkwood introduced a coupling parameter approach, and Zwanzig showed that the free energy difference between two states can be computed via an appropriate exponential average of energy differences over an ensemble of configurations. In recent years, alchemical free energy calculations have benefited from improved force fields, new sampling algorithms, and the emergence of parallel computing, which have been shown to yield accurate results. These calculations are routinely applied to large-scale free energy calculations for ligand modification, and for testing mutations in amino acids and nucleotides. Here, we report alchemical calculations on RNA/peptide and RNA/small-molecule complexes, which are conducted using free energy perturbation (FEP) and thermodynamic integration (TI) methods. We discovered that the peptides employ base-specific hydrogen bonding, sugar-phosphate backbone interactions, and cation-pi interactions for the recognition of RNA molecules. In contrast, the small molecules recognize the RNA via base-specific and pi-pi stacking interactions.

Insight into the role of excess hydroxide ions in silicate condensation reactions.

The formation of silicate oligomers in the early stages is key to zeolite synthesis. The pH and the presence of hydroxide ions are important in regulating the reaction rate and the dominant species in solutions. This paper describes the formation of silicate species, from dimers to 4-membered rings, using ab initio molecular dynamics simulations in explicit water molecules with an excess hydroxide ion. The thermodynamic integration method was used to calculate the free energy profile of the condensation reactions. The hydroxide group's role is not only to control the pH of the environment, but also to actively participate in the condensation reaction. The results show that the most favorable reactions are linear-tetramer and 4-membered-ring formation, with overall barriers of 71 kJ mol-1 and 73 kJ mol-1, respectively. The formation of trimeric silicate, with the largest free-energy barrier of 102 kJ mol-1, is the rate-limiting step under these conditions. The excess hydroxide ion aids in the stabilization of the 4-membered-ring structure over the 3-membered-ring structure. Due to a relatively high free-energy barrier, the 4-membered ring is the most difficult of the small silicate structures to dissolve in the backward reaction. This study is consistent with the experimental observation that silicate growth in zeolite synthesis is slower in a very-high-pH environment.