Thermodynamic Integration Method Research Articles

Interfacial segregation of interacting vacancies and their role on the wetting critical properties of the Blume-Emery-Griffiths model.

We study the wetting critical behavior of the three-state (s=±1,0) Blume-Emery-Griffiths model using numerical simulations. This model provides a suitable scenario for the study of the role of vacancies on the wetting behavior of a thin magnetic film. To this aim we study a system confined between parallel walls with competitive short-range surface magnetic fields (h_{L}=-|h_{1}|). We locate relevant critical curves for different values of the biquadratic interaction and use a thermodynamic integration method to calculate the surface tension as well as the interfacial excess energy and determine the wetting transition. Furthermore, we also calculate the local position of the interface along the film and its fluctuations (capillary waves), which are a measure of the interface width. To characterize the role played by vacancies on the interfacial behavior we evaluate the excess density of vacancies, i.e., the density difference between a system with and without interface. We also show that the temperature dependence of both the local position of the interface and its width can be rationalized in term of a finite-size scaling description, and we propose and successfully test the same scaling behavior for the average position of the center of mass of the vacancies and its fluctuations. This shows that the excess of vacancies can be associated to the presence of the interface that causes the observed segregation. This segregation phenomena is also evidenced by explicitly evaluating the interfacial free energy.

Free-Energy Landscape of a pH-Modulated G·C Base Pair Transition from Watson-Crick to Hoogsteen State in Duplex DNA.

In double-helical DNAs, the most stable Watson-Crick (WC) base pair (bp) can be in thermal equilibrium with much less abundant Hoogsteen (HG) bp by the spontaneous rotation of the glycosidic angle in purine bases. Previous experimental studies showed that in the case of a G·C bp, the population of the transient HG is enhanced as a protonated form (HG+) through the protonation of the cytosine base under weakly acidic conditions. Hence, pH is a key factor that can modulate this transition event from the WC to HG+ bp. In this study, to computationally probe the overall free-energy landscapes of this pH-modulated G·C HG breathing, a comprehensive classical molecular dynamics (MD) simulation protocol is proposed using an enhanced sampling MD in conjunction with the standard thermodynamic integration method. From this MD protocol proposed, the free-energy surfaces of the G·C bp transition from the WC to HG bp were constructed successfully at any pH range, producing pH-dependent free-energy quantities in close agreement with previously reported experimental results. The simulation protocol is expected to provide valuable atomistic insight into the DNA bp transition events coupled with protonation or tautomeric shift in a target bp.

Theoretical insights into C-H bond activation of methane by transition metal clusters: the role of anharmonic effects.

In heterogeneous catalysis, the determination of active phases has been a long-standing challenge, as materials' properties change under operational conditions (i.e. temperature (T) and pressure (p) in an atmosphere of reactive molecules). As a first step towards materials design for methane activation, we study the T and p dependence of the composition, structure, and stability of metal oxide clusters in a reactive atmosphere at thermodynamic equilibrium using a prototypical model catalyst having wide practical applications: free transition metal (Ni) clusters in a combined oxygen and methane atmosphere. A robust methodological approach is employed, where the starting point is systematic scanning of the potential energy surface (PES) to obtain the global minimum structures using a massively parallel cascade genetic algorithm (cGA) at the hybrid density functional level. The low energy clusters are further analyzed to estimate their thermodynamic stability at realistic T, pO2 and pCH4 using ab initio atomistic thermodynamics (aiAT). To incorporate the anharmonicity in the vibrational free energy contribution to the configurational entropy, we evaluate the excess free energy of the clusters numerically by a thermodynamic integration method with ab initio molecular dynamics (aiMD) simulation inputs. By analyzing a large dataset, we show that the conventional harmonic approximation miserably fails for this class of materials, and capturing the anharmonic effects on the vibration free energy contribution is indispensable. The latter has a significant impact on detecting the activation of the C–H bond, while the harmonic infrared spectrum fails to capture this, due to the wrong prediction of the vibrational modes.

Strong Sequestration of Hydrogen Into the Earth's Core During Planetary Differentiation

AbstractWe explore the partitioning behavior of hydrogen between coexisting metal and silicate melts at conditions of the magma ocean and the current core–mantle boundary with the help of density functional theory molecular dynamics. We perform simulations with the two‐phase and thermodynamic integration methods. We find that hydrogen is weakly siderophile at low pressure (20 GPa and 2,500 K), and becomes much more strongly so with pressure, suggesting that hydrogen is transported to the core in a significant amount during core segregation and is stable there. Based on our results, the core likely contains ~1 wt% H, assuming single‐stage formation and equilibration at 40 GPa. Our two‐phase simulations further suggest that silicon is entrained in the core‐forming metal, while magnesium remains in the silicate phase. This preferred incorporation of silicon hints at an explanation for the elevated Mg/Si ratio of the bulk silicate Earth relative to chondritic compositions.

Repulsive Soft-Core Potentials for Efficient Alchemical Free Energy Calculations.

In alchemical free energy (FE) simulations, annihilation and creation of atoms are generally achieved with the soft-core potential that shifts the interparticle separations. While this soft-core potential eliminates the numerical instability occurring near the two end states of the transformation, it makes the hybrid Hamiltonian vary nonlinearly with respect to the parameter λ, which interpolates between the Hamiltonians representing the two end states. This complicates FE estimation by Bennett acceptance ratio (BAR), free energy perturbation (FEP), and thermodynamic integration (TI) methods, thus reducing their calculation efficiency. In this work, we develop a new type of repulsive soft-core potential, called Gaussian soft-core (GSC) potential, with two parameters controlling its maximum and width. The main advantage of this potential is the linearity of the hybrid Hamiltonian with respect to λ, thus permitting the direct application of BAR, FEP, TI, and other variant FE methods. The accuracy and efficiency of the GSC potential are demonstrated by comparing the free energies of annihilation determined for 13 small molecules and an alchemical mutation of a protein side chain. In addition, in combination with a TI integrand (∂H/∂λ) estimation strategy, we show that GSC can considerably reduce the number of λ simulations compared to the commonly used separation-shifted soft-core potential.

Urea-water solvation of protein side chain models

Aqueous urea stabilizes the unfolded states of protein due to their ability to solvate both hydrophilic and hydrophobic residues favorably. The nature of interactions that stabilize different types of amino acid side chains in their solvent exposed state is still not understood. To gain insights into the molecular level details of urea interactions with proteins in their unfolded states, we have performed atomistic molecular dynamics simulations and free energy calculations using the thermodynamic integration method on model systems representing side chains of all amino acids in different solvent environments (water and varying concentrations of aqueous urea). A systematic analysis of structural, energetic and dynamic parameters has been done to understand the detailed atomistic mechanism. The main aim of the current study is to unravel the nature of urea-amino acid interactions by emphasizing on the chemical nature of amino acid side chain models. The preferential interactions of urea over water with each side chain and backbone model systems in various concentrations of aqueous urea were quantified using the two-domain model, and it is validated by mean lifetime calculations. Interestingly, almost all amino acids showed a preference for urea over water. The order of preferences depending on the chemical nature of the amino acids is obtained with the aromatic groups exhibiting the highest preferences followed by hydrophobic groups, followed by amides and basic groups, and the least by nucleophilic groups. The extensive energetic analysis revealed, these preferential interactions are enthalpically and entropically driven and are dominated by dispersion effects. Spatial density distributions and radial distribution analyses provide insights to understand the different modes and urea orientation towards preferred sites of interactions by which urea-protein interactions stabilize proteins in their unfolded states by forming favorable interactions with exposed amino acids side chains.

Magnetocaloric Properties of an Ising Antiferromagnet on a Kagome Lattice

Owing to a high degree of geometrical frustration an Ising antiferromagnet on a kagome lattice (IAKL) is known to exhibit no long-range ordering at any temperature, including the ground state. Nevertheless, at low temperatures it shows a strongly correlated, highly fluctuating regime known as a cooperative paramagnet or classical spin liquid. In the ground state it is characterized by a macroscopic degeneracy which translates to a relatively large value of the residual entropy. It has been shown that the presence of a macroscopic degeneracy associated with geometrical frustration below the saturation field can facilitate an enhanced magnetocaloric effect (MCE), which can exceed that of an ideal paramagnet with equivalent spin by more than an order of magnitude. In the present study we investigate magnetic and magnetocaloric properties of IAKL by Monte Carlo simulation. In particular, we calculate the entropy of the system using the thermodynamic integration method and evaluate quantities which characterize MCE, such as the isothermal entropy and adiabatic temperature changes in a varying magnetic field. It is found that IAKL shows the most interesting magnetocaloric properties at low temperatures and moderate magnetic fields, suggesting its potential to be used in technological applications for low-temperature magnetic refrigeration.

Yttrium complexation and hydration in chloride-rich hydrothermal fluids: A combined ab initio molecular dynamics and in situ X-ray absorption spectroscopy study

Accurate knowledge of rare earth elements (REE) speciation in high pressure – high temperature fluids is required to model REE transport and precipitation in subduction zones and magmatic-hydrothermal environments, and the formation of rare metal deposits. Recent experiments (lanthanum, ytterbium, erbium) have demonstrated that REE chloride complexes are the main REE form in many hydrothermal fluids (Migdisov et al., 2016). However, the speciation of yttrium (Y(III)), a cation with an ionic radius similar to that of Ho(III), remains poorly constrained in chloride-rich hydrothermal solutions.We used ab initio molecular dynamics (MD) simulations to calculate the nature of Y(III)-Cl complexes and the thermodynamic properties of these species at temperatures up to 500 °C and pressures of 800 bar and 1000 bar. The MD results were complemented by in situ X-ray absorption spectroscopy (XAS) measurements. Our results indicate that at temperatures below 200 °C, chloro-complexes do not form readily, even in highly concentrated brines. At ambient condition, the Y(III) aqua ion binds to eight water molecules in a square antiprism geometry, which is consistent with previous ab initio studies (Ikeda et al., 2005a). The thermodynamic integration method was employed to calculate the formation constants (logKΘ) of Y(III)-Cl− complexes in two simulation boxes containing different Y:Cl ratios; we obtained very consistent results of the standard logKΘ of the individual complexes from the two independent calculations, which confirms that the thermodynamic integration method is reliable and not significantly affected by technical limitations in box size, box composition, or simulation time.Based on the derived formation constants, we fit modified Ryzhenko–Bryzgalin (MRB) equation of state parameters, which enable extrapolation of the formation constants at elevated temperature and pressure. The results are consistent with the XAS data, and show that the stability of Y(III)-Cl complexes increases with increasing temperature, Y(III) forming high order Cl− complexes (up to YCl4−) in high salinity solutions at high temperature and pH = 3. We also compare the extrapolated logKΘ with the available data for other REE at 150 °C, 200 °C and 250 °C. At 200 °C, yttrium behaves more like a heavy REE, but from 200 °C to 250 °C, the formation constants of Y(III)-Cl complexes increase dramatically and behave more like the light REE. The difference of Cl− dominant species between Ho(III) (HoCl2+) and Y(III) (YCl2+) may account for the formation of anomalous Y/Ho ratios in some hydrothermal environments.

Free energy profile of permeation of Entecavir through Hepatitis B virus capsid studied by molecular dynamics calculation

Abstract Entecavir, triphosphorylated in liver cells, is an antiviral reagent against Hepatitis B virus (HBV). The reagent inhibits reverse transcription of RNA inside the virus capsid. In the present study, free energy profile of an Entecavir triphosphate (ETVTP) molecule has been calculated when it passes through pores of the capsid along two- and three-fold rotational symmetry axes in order to investigate permeation pathway of the reagent to the inside of the capsid. The calculations have been done based on thermodynamic integration (TI) method combined with all-atomistic molecular dynamic (MD) calculations. A free energy minimum of −19 kJ/mol was found at the entrance of the pore from the outside along the three-fold symmetry axis. This stabilization is from the interaction of negatively charged ETVTP with positively charged capsid methionine residues. This excess free energy concentrates of the reagent at the entrance of the pore by a factor of about 2000. A free energy barrier of approximately 13 kJ/mol was also found near the exit of the pore to the inside of the capsid due to narrow space of the pore surrounded by hydrophobic wall made by proline residues and negatively charged wall by aspartic acid residues. There, ETVTP is partially dehydrated in order to pass through the narrow space, which causes the great free energy loss. Further, the negatively charged residues produce repulsive forces on the ETVTP molecule. In contrast, in the case of the pore along the two-fold symmetry axis, the calculated free energy profile showed shallower free energy minimum, −4 kJ/mol at the entrance in spite of the similarly high barrier, 7 kJ/mol, near the exit of the pore.

A computational study of excess properties for mW potential model of water in supercooled region

In this work, we studied the liquid–crystal transition of water using monatomic potential model (mW) (Molinero and Moore, 2009). To simulate water, the NPT-MC simulation were performed and obtained the average properties as potential energy (〈ϕ′〉) and density (〈ρ〉) at various temperatures which are in the very good agreement (i.e. 〈ϕ′〉=−1.766 and 〈ρ〉=0.448 at 202 K) with the literature values (Hujo et al., 2011). The free energy difference between the two liquid phases of water was computed by thermodynamic integration along with a reversible path at zero pressure and various temperatures. The recently developed thermodynamic integration (TDI) method was used to computed the excess Gibbs free energy of high density liquid (HDL) phase with respect to the crystalline phase at different temperatures in the supercooled region of mW-water model (Molinero and Moore, 2009). Based on the slope of excess Gibbs free energy with respect to temperature the excess entropy of the high density liquid (HDL) phase was obtained, which shows an anomalous behavior at or near the liquid–liquid transition temperature of Tll=202 K (Molinero and Moore, 2009). The excess entropy of the liquid phase in the supercooled mW water shows a sharp decline at or near the liquid–liquid phase transition temperature (Tll). The excess entropy of the HDL phase decreases as the temperature decreases to 197 K and increases at 192 K. Our results are in good agreement with the previous observation of a anomalous dependency of the density on temperature in MD simulations starting in the HDL phase at a temperature just above Tll (Moore and Molinero, 2009).

Evaluation of the 1-octanol/water partition coefficient of quinolones and quinolones compounds via free energy estimated in quantum chemical calculations

A rapid screening assay for the measurement of octanol-water partition coefficients (logP) of pharmaceuticals such as quinolines, quinolones and fluoroquinolones compounds were developed by using quantum chemical calculations. In this paper, we report the logP of the pharmaceutical derivatives compounds that determined by the calculation of the difference in the solvation free energies in water-octanol using the thermodynamic integration method. The logP values calculated in this method are closer to the experimental values compared to other ab initio methods. Solvation free energy in water and octanol, free energy of cavity formation in water and Henry’s constants, and some other parameters are determined at the density functional theory (DFT) and Hartree-Fock (HF) level with 6-31++G** basis set. Moreover, surface area, mass, refractivity, volume, polarizability have been calculated for some of the compounds with the same levels.

Monte Carlo simulation and thermodynamic integration applied to protein charge transfer.

We introduce a combination of Monte Carlo simulation and thermodynamic integration methods to address a model problem in free energy computations, electron transfer in proteins. The feasibility of this approach is tested using the ferredoxin protein from Clostridium acidurici. The results are compared to numerical solutions of the Poisson-Boltzmann equation and data from recent molecular dynamics simulations on charge transfer in a protein complex, the NrfHA nitrite reductase of Desulfovibrio vulgaris. Despite the conceptual and computational simplicity of the Monte Carlo approach, the data agree well with those obtained by other methods. A link to experiments is established via the cytochrome subunit of the bacterial photosynthetic reaction center of Rhodopseudomonas viridis.

Global Thermodynamic Properties of Complex Spin Systems Calculated from Density of States and Indirectly by Thermodynamic Integration Method

Evaluation of global thermodynamic properties such as the entropy or the free energy of complex systems featuring a high degree of frustration or disorder is often desirable. Nevertheless, they cannot be measured directly in standard Monte Carlo simulation. Therefore, they are either evaluated indirectly from the directly measured quantities, for example by the thermodynamic integration method (TIM), or by applying more sophisticated simulation methods, such as the Wang-Landau (WL) algorithm, which can directly sample density of states. In the present investigation we compare the performance of the WL and TIM methods for the calculation of the entropy of an Ising antiferromagnetic system on a Kagome lattice – a typical example of a complex spin system with high geometrical frustration resulting in a non-zero residual entropy the value of which is exactly known. It is found that the easier to implement TIM can yield results of comparable accuracy with that of the more involved WL method.

Atomic insights into the effects of pathological mutants through the disruption of hydrophobic core in the prion protein

Destabilization of prion protein induces a conformational change from normal prion protein (PrPC) to abnormal prion protein (PrPSC). Hydrophobic interaction is the main driving force for protein folding, and critically affects the stability and solvability. To examine the importance of the hydrophobic core in the PrP, we chose six amino acids (V176, V180, T183, V210, I215, and Y218) that make up the hydrophobic core at the middle of the H2-H3 bundle. A few pathological mutants of these amino acids have been reported, such as V176G, V180I, T183A, V210I, I215V, and Y218N. We focused on how these pathologic mutations affect the hydrophobic core and thermostability of PrP. For this, we ran a temperature-based replica-exchange molecular dynamics (T-REMD) simulation, with a cumulative simulation time of 28 μs, for extensive ensemble sampling. From the T-REMD ensemble, we calculated the protein folding free energy difference between wild-type and mutant PrP using the thermodynamic integration (TI) method. Our results showed that pathological mutants V176G, T183A, I215V, and Y218N decrease the PrP stability. At the atomic level, we examined the change in pair-wise hydrophobic interactions from valine-valine to valine-isoleucine (and vice versa), which is induced by mutation V180I, V210I (I215V) at the 180th–210th (176th–215th) pair. Finally, we investigated the importance of the π-stacking between Y218 and F175.

Investigation on the efficiency and accuracy of methods for calculating melting temperature by molecular dynamics simulation

Finding the efficient and accurate schemes to calculate the melting temperature is important for molecular dynamics simulation. We propose a modified void method based on heterogeneous nucleation to validly calculate the melting temperature from solid-liquid coexistence state. It decreases the degree of superheating in the conventional voids method. The efficiency and accuracy of five methods (i.e., hysteresis method, two phase coexistence method, the interface pinning method, the Frenkel-Ladd path method and the modified void method) have been discussed in detail. The calculated results of hysteresis method strongly depend on heating/cooling rate. For the nonequilibrium thermodynamic integration method, the calculated accuracy is consistent with the interface pinning method for monoatomic Cu. For the modified void method the calculated results agree with the nonequilibrium thermodynamic integration method.

Steering a solute between coexisting solvation states: Revisiting nonequilibrium work relations and the calculation of free energy differences.

By analogy with single-molecule pulling experiments, we present a computational framework to obtain free energy differences between complex solvation states. To illustrate our approach, we focus on the calculation of solvation free energies (SFEs). However, the method can be readily extended to cases involving more complex solutes and solvation conditions as well as to the calculation of binding free energies. The main idea is to drag the solute across the simulation box where atomistic and ideal gas representations of the solvent coexist at constant temperature and chemical potential. At finite pulling speeds, the resulting work allows one to extract SFEs via nonequilibrium relations, whereas at infinitely slow pulling speeds, this process becomes equivalent to the thermodynamic integration method. Results for small molecules well agree with literature data and pave the way to systematic studies of arbitrarily large and complex molecules.

Development of a Robust Indirect Approach for MM → QM Free Energy Calculations That Combines Force-Matched Reference Potential and Bennett's Acceptance Ratio Methods.

We use the PBE0/6-31G* density functional method to perform ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations under periodic boundary conditions with rigorous electrostatics using the ambient potential composite Ewald method in order to test the convergence of MM → QM/MM free energy corrections for the prediction of 17 small-molecule solvation free energies and eight ligand binding free energies to T4 lysozyme. The "indirect" thermodynamic cycle for calculating free energies is used to explore whether a series of reference potentials improve the statistical quality of the predictions. Specifically, we construct a series of reference potentials that optimize a molecular mechanical (MM) force field's parameters to reproduce the ab initio QM/MM forces from a QM/MM simulation. The optimizations form a systematic progression of successively expanded parameters that include bond, angle, dihedral, and charge parameters. For each reference potential, we calculate benchmark quality reference values for the MM → QM/MM correction by performing the mixed MM and QM/MM Hamiltonians at 11 intermediate states, each for 200 ps. We then compare forward and reverse application of Zwanzig's relation, thermodynamic integration (TI), and Bennett's acceptance ratio (BAR) methods as a function of reference potential, simulation time, and the number of simulated intermediate states. We find that Zwanzig's equation is inadequate unless a large number of intermediate states are explicitly simulated. The TI and BAR mean signed errors are very small even when only the end-state simulations are considered, and the standard deviations of the TI and BAR errors are decreased by choosing a reference potential that optimizes the bond and angle parameters. We find a robust approach for the data sets of fairly rigid molecules considered here is to use bond + angle reference potential together with the end-state-only BAR analysis. This requires QM/MM simulations to be performed in order to generate reference data to parametrize the bond + angle reference potential, and then this same simulation serves a dual purpose as the full QM/MM end state. The convergence of the results with respect to time suggests that computational resources may be used more efficiently by running multiple simulations for no more than 50 ps, rather than running one long simulation.

Effect of the Solvent on the Oxygen Evolution Reaction at the TiO2–Water Interface

We investigate the solvation effect of water on the overpotentials of the oxygen evolution reaction on rutile TiO2 by applying the thermodynamic integration method on atomistic model interfaces with and without the water molecules. We compare the results at the vacuum interface with the commonly used computational hydrogen electrode method, finding overall good agreement. The effect of the solvent is found to be twofold. First, the explicit treatment of the solvent can lead to equilibrium configurations differing from the relaxed structures without solvent. Second, the overpotentials can be affected by up to 0.5 eV. The energetics are subject to electrostatic effects at the interface rather than to modifications in the hydrogen bond network. These results provide a promising perspective on the use of implicit models for treating the solvent.

Energetic, Structural and Dynamic Properties of Nucleobase-Urea Interactions that Aid in Urea Assisted RNA Unfolding

Understanding the structure-function relationships of RNA has become increasingly important given the realization of its functional role in various cellular processes. Chemical denaturation of RNA by urea has been shown to be beneficial in investigating RNA stability and folding. Elucidation of the mechanism of unfolding of RNA by urea is important for understanding the folding pathways. In addition to studying denaturation of RNA in aqueous urea, it is important to understand the nature and strength of interactions of the building blocks of RNA. In this study, a systematic examination of the structural features and energetic factors involving interactions between nucleobases and urea is presented. Results from molecular dynamics (MD) simulations on each of the five DNA/RNA bases in water and eight different concentrations of aqueous urea, and free energy calculations using the thermodynamic integration method are presented. The interaction energies between all the nucleobases with the solvent environment and the transfer free energies become more favorable with respect to increase in the concentration of urea. Preferential interactions of urea versus water molecules with all model systems determined using Kirkwood-Buff integrals and two-domain models indicate preference of urea by nucleobases in comparison to water. The modes of interaction between urea and the nucleobases were analyzed in detail. In addition to the previously identified hydrogen bonding and stacking interactions between urea and nucleobases that stabilize the unfolded states of RNA in aqueous solution, NH-π interactions are proposed to be important. Dynamic properties of each of these three modes of interactions have been presented. The study provides fundamental insights into the nature of interaction of urea molecules with nucleobases and how it disrupts nucleic acids.

Investigation of Chitosan, Quercetin and Water Interactions using Free Energy Calculations and Molecular Docking

Predicting the preferred orientation of a ligand to a receptor and the free energy (ΔG AB ) of these molecules in a solution is very useful for pharmaceutical applications. In this work, ΔG AB between chitosan and quercetin molecules with water were obtained by the thermodynamic integration (TI) method and the chitosan+quercetin complex was investigated by Molecular Docking (MDK). The topology of the molecules was compatible with the CHARMM force-field. Both molecules were solvated with the GROMACS SPC water model. The ΔG AB values of the molecule+water systems with the thermodynamic integration (TI) method was done in a 10 ns simulation. The complex was achieved in a 2.5 nm x 2.5 nm x 2.5 nm box. The ΔG AB values for the molecules+water systems were 14.36 +/- 0.46 and 0.19 +/- 0.53 kJ/mol, respectively. For the complex, nine conformations were obtained for the quercetin around the chitosan. The best conformation showed an affinity energy of -3.8 kcal/mol. The ΔG AB results for the systems confirm the low or no solubility of the molecules in the water under normal conditions, as already demonstrated experimentally. The docking showed that the interaction between chitosan and quercetin molecules involves hydrogen bonds in a specific position. DOI: http://dx.doi.org/10.17807/orbital.v11i2.1344