Thermodynamic Integration Technique Research Articles

How Solvation Alters the Thermodynamics of Asymmetric Bond-Breaking: Quantum Simulation of NaK+ in Liquid Tetrahydrofuran.

Gas-phase potential energy surfaces (PESs) are often used to provide an intuitive understanding of molecular chemical reactivity. Most chemical reactions, however, take place in solution, and it is unclear whether gas-phase PESs accurately represent chemical processes in solvent environments. In this work we use quantum simulations to investigate the dissociation energetics of NaK+ in liquid tetrahydrofuran (THF) to understand the degree to which solvent interactions alter the gas-phase picture. Using umbrella sampling and thermodynamic integration techniques, we construct condensed-phase free energy surfaces of NaK+ on THF in both the ground and electronic excited states. We find that solvation by THF completely alters the nature of the NaK+ bond by reordering the thermodynamic dissociation products. Reaching the thermodynamic dissociation limit in THF also requires a long-range charge transfer process that has no counterpart in the gas phase. Gas-phase PESs, even with perturbations, cannot adequately describe the reactivity of simple asymmetric molecules in solution.

Deep‐Learning‐Based Prediction of the Tetragonal → Cubic Transition in Davemaoite

AbstractDavemaoite, that is, CaSiO3 perovskite (CaPv), is the third most abundant phase in the lower mantle and exhibits a tetragonal‐cubic phase transition at high pressures and temperatures. The phase boundary in CaPv has recently been proposed to be close to the cold slab adiabat and cause mid‐mantle seismic wave speed anomalies (Thomson et al., 2019, https://doi.org/10.1038/s41586‐019‐1483‐x). This study utilized accurate deep‐learning‐based simulations and thermodynamic integration techniques to compute free energies at temperatures ranging from 300 to 3,000 K and pressures up to 130 GPa. Our results indicate that CaPv exhibits a single cubic phase throughout lower‐mantle conditions. This suggests that the phase diagram proposed by Thomson et al. requires revision, and mid‐mantle seismic anomalies are likely attributable to other mechanisms.

A fast atomistic approach to finite-temperature surface elasticity of crystalline solids

Surface energies and surface elasticity largely affect the mechanical response of nanostructures as well as the physical phenomena associated with surfaces such as evaporation and adsorption. Studying surface energies at finite temperatures is therefore of immense interest for nanoscale applications. However, calculating surface energies and derived quantities from atomistic ensembles is usually limited to zero temperature or involves cumbersome thermodynamic integration techniques at finite temperature. Here, we illustrate a computational technique to identify the energy and elastic properties of surfaces of solids at non-zero temperature based on a Gaussian phase packets (GPP) approach (which in the isothermal limit coincides with a maximum-entropy formulation). Using this technique, we investigate the effect of temperature on the surface properties of different crystal faces for six pure metals – copper, nickel, aluminium, iron, tungsten and vanadium – thus covering both FCC and BCC lattice structures. While the obtained surface energies and stresses usually show a decreasing trend with increasing temperature, the elastic constants do not show such a consistent trend across the different materials and are quite sensitive to temperature changes. Validation is performed by comparing the obtained surface energy densities of selected BCC and FCC materials to those calculated via molecular dynamics.

The phase diagrams of beryllium and magnesium oxide at megabar pressures

We perform ab initio simulations of beryllium (Be) and magnesium oxide (MgO) at megabar pressures and compare their structural and thermodynamic properties. We make a detailed comparison of our two recently derived phase diagrams of Be (Wu et al 2021 Phys. Rev. B 104 014103) and MgO (Soubiran and Militzer 2020 Phys. Rev. Lett. 125 175701) using the thermodynamic integration technique, as they exhibit striking similarities regarding their shape. We explore whether the Lindemann criterion can explain the melting temperatures of these materials through the calculation of the Debye temperature at high pressure. From our free energy calculations, we find that the melting line of both materials is well represented by the Simon–Glazel fit T m(P) = T 0(1 + P/a)1/c , where T 0 = 1564 K, a = 15.8037 GPa and c = 2.4154 for Be, while T 0 = 3010 K, a = 10.5797 GPa and c = 2.8683 for the MgO in the B1. For the B2 phase, we use the values a = 26.1163 GPa and c = 2.2426. Both materials exhibit negative Clapeyron slopes on the boundaries between the two solid phases that are strongly affected by anharmonic effects, which also influence the location of the solid–solid–liquid triple point. We find that the quasi-harmonic approximation underestimates the stability range of the low-pressure phases, namely hcp for Be and B1 for MgO. We also compute the phonon dispersion relations at low and high pressure for each of the phases of these materials, and also explore how the phonon density of states is modified by temperature. Finally, we derive secondary shock Hugoniot curves in addition to the principal Hugoniot curve for both materials, and study their offsets in pressure between solid and liquid branches.

Molecular simulation of pentaethylene glycol monooctyl ether micelles in water. Prediction of the micellization Gibbs energy

The Gibbs energy (Δg0) and the enthalpy (Δh0) of the micellization of pentaethylene glycol monooctyl ether (C8E5) were predicted by molecular dynamics (MD) simulations of preformed micelles of various sizes at 298.15 K. Macroscopically, the equilibrium size distributions of the micellar aggregates in amphiphile solutions were modeled using a multiple chemical equilibrium approach, under the assumption that the micelles are spheres or sphero-cylinders. In the model, Δg0 is a function of the micelle size, and the parameters for its determination are related to the contributions of the amphiphiles on the micelle’s spherical caps and cylindrical body. To estimate these contributions, the thermodynamic integration technique was employed over the MD simulations. For Δg0, the predicted parameters are very similar to those determined by fitting the macroscopic model to measured tensiometric and calorimetric data. However, Δh0 was underpredicted compared to experimental data. The parameters calculated via MD allowed the prediction of the critical micelle concentration and of the mean aggregation number, in very good agreement with those previously published from experimental calorimetric and tensiometric data.

Ab initio surface free energies of tungsten with full account of thermal excitations

The surface free energies of seven different facets of tungsten (W) are obtained up to the melting point with full account of all the relevant thermal excitations; in particular, thermal atomic vibrations, electronic excitations, and their mutual coupling. The latter is done using ab initio molecular dynamics simulations coupled with the thermodynamic integration technique. In this way, the calculations contain almost no error but the one related to the used exchange-correlation functional, which makes the results truly first principles. The obtained results are compared with previous quasiharmonic calculations for the surface free energies of W and experimental data. The anharmonic contribution is, as expected, important for open surfaces at high temperatures, which leads to a temperature dependence of the surface energy anisotropy. The calculated Wulff shapes and surface energies are in excellent agreement with experimental data close to the melting point, where the crystalline structure of the surface layers is destroyed by a dramatic mobility of the atoms there.

The Earth’s core as a reservoir of water

Current estimates of the budget and distribution of water in the Earth have large uncertainties, most of which are due to the lack of information about the deep Earth. Recent studies suggest that the Earth could have gained a considerable amount of water during the early stages of its evolution from the hydrogen-rich solar nebula, and that a large amount of the water in the Earth may have partitioned into the core. Here we calculate the partitioning of water between iron and silicate melts at 20–135 GPa and 2,800–5,000 K, using ab initio molecular dynamics and thermodynamic integration techniques. Our results indicate a siderophile nature of water at core–mantle differentiation and core–mantle boundary conditions, which weakens with increasing temperature; nevertheless, we found that water always partitions strongly into the iron liquid under core-formation conditions for both reducing and oxidizing scenarios. The siderophile nature of water was also verified by an empirical-counting method that calculates the distribution of hydrogen in an equilibrated iron and silicate melt. We therefore conclude that the Earth’s core may act as a large reservoir that contains most of the Earth’s water. In addition to constraining the accretion models of volatile delivery, the findings may partially account for the low density of the Earth’s core implied by measured seismic velocities. The Earth’s core may host most of the planet’s water inventory, according to calculations of the partitioning behaviour of water at conditions of core formation.

QM/MM-MD dissociation of Ag+ and H+ mediated cytosine pairs: Monomers and dimers

We report a biased QM/MM molecular dynamics study of the dissociation process of cytosine-cytosine complexes mediated by Ag+ or H+ ions in monomeric and dimeric forms. We performed calculations under explicit solvent conditions and obtained the free energy profiles by thermodynamic integration technique to give deep insights on the dissociation process. For all geometries corresponding to key points on energy profiles the noncovalent interaction descriptors were calculated and the details of dissociation mechanism were revealed. Our findings by means of dissociation energy barrier analysis for Ag+-mediated cytosines suggested more favorable cis-configuration of C(Ag+)C with the energy of 17.8 kcal/mol over trans-configuration with the energy of 12.2 kcal/mol in monomeric form, in contrast to H+ mediated cytosines. Also, it is shown that the QM/MM-MD global minimum of dissociation free energy profiles for both monomer isomers (cis- and trans-) doesn’t correspond to the geometries found in crystal data and previous theoretical studies of equilibrium geometries. But the addition of an adjacent silver mediated cytosine pair, in which additional stabilizing factors arise such as interplane H-bonding, π-π stacking or argentophilic interaction, leads to a change in the dissociation profiles and improving the agreement with experiment. The trans-configuration of (C(Ag+)C)2 with the energy barrier of 16.6 kcal/mol is more favorable in the dissociation process then cis-configuration with the energy of 15.6 kcal/mol.

Carbon Partitioning Between the Earth's Inner and Outer Core

AbstractKnowledge of the abundance and distribution of light elements in the core is fundamental to the understanding of the Earth and other planetary systems. Recent studies (Li et al., 2018; Mashino et al., 2019) suggest the particular importance of carbon for explaining core properties, yet knowledge of carbon partitioning between the outer and inner core is unknown. By using the quasiharmonic approximation, ab initio molecular dynamics, and thermodynamic integration techniques, we have computed the chemical potential of carbon in liquid Fe and solid hcp‐Fe at core conditions. We find that substitutional carbon is more stable than interstitial carbon and other carbon defect cluster structures in solid Fe. Lattice strain and overcoordination effects lead to a high chemical potential of C in solid Fe compared to the liquid, and consequently carbon partitions almost completely into the liquid. We find that carbon can account for most of the density jump at the inner‐core boundary. This provides an alternative mechanism to the necessity of an oxygen‐rich outer core and may have significant implications for the composition and structure of the deep Earth.

TI-Stan: Model Comparison Using Thermodynamic Integration and HMC

We present a novel implementation of the adaptively annealed thermodynamic integration technique using Hamiltonian Monte Carlo (HMC). Thermodynamic integration with importance sampling and adaptive annealing is an especially useful method for estimating model evidence for problems that use physics-based mathematical models. Because it is based on importance sampling, this method requires an efficient way to refresh the ensemble of samples. Existing successful implementations use binary slice sampling on the Hilbert curve to accomplish this task. This implementation works well if the model has few parameters or if it can be broken into separate parts with identical parameter priors that can be refreshed separately. However, for models that are not separable and have many parameters, a different method for refreshing the samples is needed. HMC, in the form of the MC-Stan package, is effective for jointly refreshing the ensemble under a high-dimensional model. MC-Stan uses automatic differentiation to compute the gradients of the likelihood that HMC requires in about the same amount of time as it computes the likelihood function itself, easing the programming burden compared to implementations of HMC that require explicitly specified gradient functions. We present a description of the overall TI-Stan procedure and results for representative example problems.

Phase diagram of uranium fromab initiocalculations and machine learning

Experimental studies of materials at extreme conditions are challenging, and as a consequence, P-T phase diagrams are still unknown for many elements and materials. In this work, we present the P-T phase diagram of uranium calculated up to extreme conditions. First, we searched for possible crystal structures using the evolutionary algorithm USPEX. Their free energies were then calculated using thermodynamic integration (TI) and temperature-dependent effective potential techniques. TI was performed using molecular dynamics, employing a machine learning (ML) force field trained on energies and forces from density-functional calculations at the generalized gradient approximation level. The prediction error of the ML force field for the energy was less than 10 meV/atom. Using thermodynamic perturbation theory (including first and second order corrections), from the free energies of the ML force field, we obtained free energies and phase diagram at the level of quality of the underlying density-functional calculations at pressures up to 800 GPa and temperatures up to 16 000 K.

Closure for the Ornstein-Zernike equation with pressure and free energy consistency.

The Ornstein-Zernike (OZ) integral equation theory is a powerful approach to simple liquids due to its low computational cost and the fact that, when combined with an appropriate closure equation, the theory is thermodynamically complete. However, approximate closures proposed to date exhibit pressure or free energy inconsistencies that produce inaccurate or ambiguous results, limiting the usefulness of the Ornstein-Zernike approach. To address this problem, we combine methods to enforce both pressure and free energy consistency to create a new closure approximation and test it for a single-component Lennard-Jones fluid. The closure is a simple power series in the direct and total correlation functions for which we have derived analytical formulas for the excess Helmholtz free energy and chemical potential. These expressions contain a partial molar volumelike term, similar to excess chemical potential correction terms recently developed. Using our bridge approximation, we have calculated the pressure, Helmholtz free energy, and chemical potential for the Lennard-Jones fluid using the Kirkwood charging, thermodynamic integration techniques, and analytic expressions. These results are compared with those from the hypernetted chain equation and the Verlet-modified closure against Monte Carlo and equations-of-state data for reduced densities of ρ^{*}<1 and temperatures of T^{*}=1.5, 2.74, and 5. Our closure shows consistency among all thermodynamic paths, except for one expression of the Gibbs-Duhem relation, whereas the hypernetted chain equation and the Verlet-modified closure exhibit consistency between only a few relations. Accuracy of the closure is comparable to the Verlet-modified closure and a significant improvement to results obtained from the hypernetted chain equation.

Ab Initio Prediction of Potassium Partitioning Into Earth's Core

AbstractPotassium partitioning between molten silicates and liquid iron alloys is the fundamental process determining its incorporation into the Earth's core. In this study, it is investigated using the method of the ab initio molecular dynamics simulation combined with the thermodynamic integration technique. Results suggest that the potassium incorporation into iron alloys positively depends on temperature, while the effect of pressure is insignificant. Moreover, the existence of oxygen in liquid iron alloys significantly enhances the potassium solubility therein, whereas sulfur and silicon only have negligible effects. Electronic structure analyses reveal that potassium remains alkali‐metallic in liquid iron alloy systems under all conditions in this study, which is distinct from the characteristics reported for solid potassium. Atomic structure analyses indicate that the oxygen coordination number around potassium atom increases with oxygen concentration in liquid iron alloys, supporting the oxygen concentration dependence of the potassium partitioning between molten silicates and liquid iron alloys. Using the obtained partitioning coefficients combined with geochemical property, the maximum potassium concentration in the Earth's core is estimated to be ~30 ppm, which would be unlikely to affect the thermal evolution of the Earth's core significantly.

Role of a reaction coordinate in free energy calculations

The free energy calculation method emerges as a viable technique for ‘in-silico’ calorimetry. Efficient sampling techniques and the good choice of a reaction path connecting the reactant and the product state enable accurate computations of the free energy differences. We argue that in many cases the thermodynamic integration technique has the lowest variance when the transformation between the reactant and the product state proceeds along the natural path of the studied chemical reaction. We provide examples of free energy calculations for the fragmentation of the charged clusters and the swapping reaction of oligomer formation in proteins that follow a tentative reaction mechanism.

Free energy and entropy of a dipolar liquid by computer simulations.

Thermodynamic properties for a system composed of dipolar molecules are computed. Free energy is evaluated by means of the thermodynamic integration technique, and it is also estimated by using a perturbation theory approach, in which every molecule is modeled as a hard sphere within a square well, with an electric dipole at its center. The hard sphere diameter, the range and depth of the well, and the dipole moment have been calculated from properties easily obtained in molecular dynamics simulations. Connection between entropy and dynamical properties is explored in the liquid and supercooled states by using instantaneous normal mode calculations. A model is proposed in order to analyze translation and rotation contributions to entropy separately. Both contributions decrease upon cooling, and a logarithmic correlation between excess entropy associated with translation and the corresponding proportion of imaginary frequency modes is encountered. Rosenfeld scaling law between reduced diffusion and excess entropy is tested, and the origin of its failure at low temperatures is investigated.

Activity coefficients from molecular simulations using the OPAS method.

A method for determining activity coefficients by molecular dynamics simulations is presented. It is an extension of the OPAS (osmotic pressure for the activity of the solvent) method in previous work for studying the solvent activity in electrolyte solutions. That method is extended here to study activities of all components in mixtures of molecular species. As an example, activity coefficients in liquid mixtures of water and methanol are calculated for 298.15 K and 323.15 K at 1 bar using molecular models from the literature. These dense and strongly interacting mixtures pose a significant challenge to existing methods for determining activity coefficients by molecular simulation. It is shown that the new method yields accurate results for the activity coefficients which are in agreement with results obtained with a thermodynamic integration technique. As the partial molar volumes are needed in the proposed method, the molar excess volume of the system water + methanol is also investigated.

Calculating free energies of organic molecules on insulating substrates.

The challenges and limitations in calculating free energies and entropies of adsorption and interaction of organic molecules on an insulating substrate are discussed. The adhesion of 1,3,5-tri(4'-cyano-[1,1'-biphenyl]-4-yl)benzene (TCB) and 1,4-bis(4-cyanophenyl)-2,5-bis(decyloxy)benzene (CDB) molecules to step edges on the KCl(001) surface and the formation of molecular dimers were studied using classical molecular dynamics. Both molecules contain the same anchoring groups and benzene ring structures, yet differ in their flexibility. Therefore, the entropic contributions to their free energy differ, which affects surface processes. Using potential of mean force and thermodynamic integration techniques, free energy profiles and entropy changes were calculated for step adhesion and dimer formation of these molecules. However, converging these calculations is nontrivial and comes at large computational cost. We illustrate the difficulties as well as the possibilities of applying these methods towards understanding dynamic processes of organic molecules on insulating substrates.

Water around fullerene shape amphiphiles: A molecular dynamics simulation study of hydrophobic hydration.

Fullerene C60 sub-colloidal particle with diameter ∼1 nm represents a boundary case between small and large hydrophobic solutes on the length scale of hydrophobic hydration. In the present paper, a molecular dynamics simulation is performed to investigate this complex phenomenon for bare C60 fullerene and its amphiphilic/charged derivatives, so called shape amphiphiles. Since most of the unique properties of water originate from the pattern of hydrogen bond network and its dynamics, spatial, and orientational aspects of water in solvation shells around the solute surface having hydrophilic and hydrophobic regions are analyzed. Dynamical properties such as translational-rotational mobility, reorientational correlation and occupation time correlation functions of water molecules, and diffusion coefficients are also calculated. Slower dynamics of solvent molecules—water retardation—in the vicinity of the solutes is observed. Both the topological properties of hydrogen bond pattern and the "dangling" -OH groups that represent surface defects in water network are monitored. The fraction of such defect structures is increased near the hydrophobic cap of fullerenes. Some "dry" regions of C60 are observed which can be considered as signatures of surface dewetting. In an effort to provide molecular level insight into the thermodynamics of hydration, the free energy of solvation is determined for a family of fullerene particles using thermodynamic integration technique.

Atomistic simulations of liquid–liquid coexistence in confinement: comparison of thermodynamics and kinetics with bulk

We review a few simulation methods and results related to the structure and non-equilibrium dynamics in the coexistence region of immiscible symmetric binary fluids, in bulk as well as under confinement, with special emphasis on the latter. Monte Carlo methods to estimate interfacial tensions for flat and curved interfaces have been discussed. The latter, combined with a thermodynamic integration technique, provides contact angles for coexisting fluids attached to the wall. For such three-phase coexistence, results for the line tension are also presented. For the kinetics of phase separation, various mechanisms and corresponding theoretical expectations have been discussed. A comparative picture between the domain growth in bulk and confinement (including thin-film and semi-infinite geometry) has been presented from molecular dynamics simulations. Applications of finite-size scaling technique have been discussed in both equilibrium and non-equilibrium contexts.

High-temperature miscibility of iron and rock during terrestrial planet formation

The accretion of a terrestrial body and differentiation of its silicate/oxide mantle from iron core provide abundant energy for heating its interior to temperatures much higher than the present day Earth. The consequences of differentiation on the structure and composition of planets are typically addressed considering only the interaction of molten iron with an immiscible ‘rocky’ phase. We demonstrate that mixing in a representative system of liquid or solid MgO and liquid iron to a single homogeneous liquid occurs at sufficiently low temperature to be present in the aftermath of a giant impact. Applying the thermodynamic integration technique to density functional theory molecular dynamics simulations, we determine the solvus closure temperature for the Fe–MgO system for pressures up to 400 GPa. Solvus closure occurs at ∼4000 K at low pressure, and has a weak positive pressure dependence, such that its gradient with respect to depth is less steep than an adiabatic temperature profile. This predicts a new mode of core–mantle differentiation following the most energetic giant impacts, with exsolution of iron from the mixture beginning in the outer layers of the planet. We demonstrate that high-temperature equilibration results in delivery of nominally insoluble Mg-rich material to the early core. Since MgO is the least soluble major mantle component in iron at low temperatures, these results may represent an upper bound on temperature for mixing in terrestrial planets.