《1. Introduction》

1. Introduction

It is somewhat ironic that while we understand a great deal about our solar system and about the surrounding stars and galaxies up to a certain distance, exact knowledge on the constituents of our own planet, and on their interplay, is rather limited. The deepest borehole ever drilled only penetrates 12 km into the crust [1]; thus, in principle, we do not have direct evidence of the materials buried further down and of their behaviors. Most of our information on the composition of the Earth’s mantle and core has been obtained indirectly from sources such as seismology, extruded minerals due to volcanic activities, or minerals embedded in diamonds. In this regard, precise knowledge of the physics and chemistry of materials under extreme conditions is paramount in order to construct coherent, logical, and reliable interpretations based on these findings. In the early days, the elastic constants of various rock types were used to interpret the density profile as a function of depth. Since then, theory and modeling have been indispensable in geological research [2]. The amalgamation of basic physics and chemistry principles has led to many successes and discoveries [3]. One example is the theory of chemical bonds, which was developed in the late 19th century, and has since been refined by many later scientists; this theory has led to predictive principles and rules that can be used to interpret the crystal structure of many minerals under ambient pressure, requiring only information on several key structural parameters and physical properties [4,5]. Arguably the greatest improvement in this theoretical arena that has been made in the last 20 years is that many quantities can now be computed reliably from first principles [6,7]; the method is based on quantum theory and involves solving the Schrödinger equation with only knowledge of the constituent elements as input. Recently, even crystal structures can be predicted reliably. Computer modeling has become a new paradigm that opens up a new dimension in the study of the behavior of materials under extreme pressure and temperature that are not easily amendable to experimentation. This article reviews several computational approaches with an emphasis on applications in the investigation of chemical reactions and transport properties. It is not the intention of this short review to provide a comprehensive survey of the state of the art in this field, or to gather all relevant articles. Rather, the goal is to provide the general reader with a synopsis of the current status of this field, drawing examples from the author’s recent work in geoscience.

The property of a material is solely determined by the interatomic forces between its atoms. Once the interatomic forces are defined, in principle, all properties can be computed. For example, static properties such as the crystal structure correspond to the lowest energy arrangement of the atoms in the unit cell; the elasticity is the response of the system to an external stress. Dynamic and transport properties at finite temperature can also be computed using classical molecular dynamics (MD) in conjunction with statistical mechanical theory [8]. For example, thermal conductivity is the transport of energy through a material, and the shear viscosity of a liquid is the elastic (viscous) stress arising from the internal forces due to deformation related to the strain rate. The calculation of material properties from numerical simulations is not new. In fact, almost all theoretical formalisms have already been established [9–12] since the MD method was conceived and implemented in the 1960s. Thus, once the interatomic forces for a system are defined by fitting the interatomic potential energy functions either to experimental observables (empirical) or to high-level quantum chemistry calculations on atomic or molecular clusters (non-empirical), a variety of properties can be predicted. For example, several highly successful interatomic potentials for Si–O, O–O and Si–Si interactions have been constructed [13,14] and have been shown to reproduce the experimental structure, elastic constant, vibrational spectrum, and structural phase transition for a family of crystalline silica and amorphous silica [15] at low pressures. Potential functions for multicomponent systems are also available and have provided useful information on the structure and dynamics of many silicate and aluminum–silicate solids and melts [16]. Classical MD calculation with reliable interatomic potentials is still a powerful predictive tool using today’s computer technology; the simulation of time duration on the order of 10 to 100 nanoseconds can be performed on millions of atoms. Very long simulation is often required for the study of slowly relaxing liquids and glasses. A major drawback of the predefined interatomic potential is that it can only be applied to the environment within the regime of the parameterization. Therefore, most potential parameters are not transferable. This, of course, will affect the reliability of the prediction. Thus, the behaviors of minerals under Earth’s conditions over a broad temperature and pressure range are not guaranteed using atomic potentials. To overcome this difficulty, it is necessary to resort to the first-principles quantum mechanical method without adjustable parameters. This is the subject of this mini-review.

《2. Theoretical background》

2. Theoretical background

According to the postulate of quantum mechanics [7,17], the electronic state of a many-particle system (electrons, nuclei, etc.) is uniquely defined by a wave function  which depends explicitly onthe electron and nucleus  coordinates. The total energy of the system is the expectation value of the energy Hamiltonian,  the total energy Hamiltonian, is simply the sum of the kinetic energy and potential energy. If the electron and nuclear motions (Born–Oppenheimer’s approximation) can be separated, the Hamiltonian for the electrons can be written as follows:

where m is the mass of an electron and the potential energy U is the sum of the electron–electron, electron–nuclei and nuclei–nuclei Coulomb interactions. Replacing the momentum operator by its quantum representation, the energy Hamiltonian becomes:

This is the Schrödinger equation of the system. Except for the hydrogen atom, the exact form of the electronic wave function  is unknown, and can be very complicated for many-particle systems. One method is to assume that each electron in the system only experiences a mean-field of the other electrons, and that each electron has its own distinct wave function  (the independent particle approximation). Then, the manyelectron total wave function ，(antisymmetrized due to Pauli’s exclusion principle) can be approximated by a determinant of the product of the individual and independent orbital wave functions. Substituting the total wave function into Eq. (2) and employing the variation principle to minimize the total energy with respect to the electron orbitals, the ground state wave function and energy can be obtained. This procedure results in a set of coupled differential equations that can be solved in an iterative manner. This approach, known as the Hartree–Fock method, has been shown to be quite efficient for molecules; however, it is laborious for extended systems due to the need to evaluate a large number of multi-center two-electron Coulomb, ，and exchange, integrals.

A different formulation of the electronic structure theory that is solely based on the electron density without invoking the wave function explicitly was formulated by Hohenberg and Kohn in 1964 [18]. The density function theory (DFT) proved that electron density is the primary variable and that the total energy of an interacting, many-electron system is a functional of the electron density. Later, Sham and Kohn [19] proposed a practical scheme to compute the total energy of an electronic system. The DFT theory eliminates the burden of the calculation of the exchange integrals, and replaces the electron exchange and correlation interactions with a universal functional that is a function of density, although the exact form of this functional is unknown. The Kohn– Sham (KS) equation for the total energy of a system can be written as follows:

where n is the electron density at. The functional is the sum of the potential energies of the electron–electron and electron– nuclei interactions, and the electron kinetic energy. Once again, the exact form of   is not known. It is easy to separate the contribution of the electronic Coulomb interaction from the total electronic interactions;  can be further decomposed into the following:

where is still a functional containing the kinetic energy of a system of non-interacting electrons with density n [19], and the remaining electronic interaction terms are grouped into , which is known as the exchange–correlation (XC) term.

According to the variational theorem, the ground state total energy is obtained from minimizing the total energy against the electron density subject to the constraint that the number particle, N, in the system is conserved; that is, . Thus, taking the differential of Eq. (3) yields

Here, μ is the chemical potential of a uniform gas with electron density n. The XC potential is defined as .

Since the exact kinetic energy functional for an interacting electron system is unknown, it is necessary to approximate the kinetic energy from the many-body wave function as a product of the oneelectron wave functions, as in the Hartree–Fock method: Each electron in the interacting electronic system satisfies the one-electron KS equation:

where the effective potential is defined as follows:

To apply the KS DFT to practical problems, an approximation must be made for the XC functional . The first approximation is the local-density approximation (LDA) in which the XC potential [19] at each point in the non-uniform real system is set to the exchange–correlation potential of a uniform electron gas calculated by the quantum Monte Carlo at the same electron density [20]. The simple LDA approximation neglects the inhomogeneity of the electron densities and, therefore, tends to overestimate the strengths of chemical bonds near equilibrium and underestimate the interatomic distances. Consideration of the curvature of the electron density by including the electron density gradient (generalized gradient approximation, GGA) [21] improves the predicted structure or sometimes the total energy, but not both. This problem is associated with the fact that the GGA does not satisfy all the known constraints to a semi-local functional. To alleviate this problem, the kinetic energy density is introduced to account for the slow variation of the electron density for metallic bonds and the tighter density for covalent bonds. This type of functional is known as a meta-GGA. Very recently, a ‘‘strongly strained and approximately normed” (SCAN) meta-GGA function was constructed to satisfy all 17 known exact constraints required for the exchange–correlation energy of a semi-local functional [22]. Since the orbitals φ_i have already been computed, the inclusion of the non-interacting kinetic energy term into the semilocal exchange–correlation energy adds little to the computational cost. The SCAN functional has been demonstrated to give superior geometry and total energy in comparison with all other available functionals [23]. It is by far the most efficient and reliable functional for main-group compounds.

A key ingredient of a DFT calculation is the evaluation of the forces acting on the atoms. This information is essential for geometry optimization and MD simulation. The force acting on an atom, F_i, is the negative derivative of the energy with respect to the displacement λ_j, and can be evaluated conveniently via the Hellman–Feynman theorem [24]:

in which the cumbersome integrals over the derivatives of the wave functions are simplified to a more manageable expression on the derivative of the energy Hamiltonian.

Once the functional and the basis set(s) have been chosen [25], the structure and energy for isolated or extended systems can be computed. The application of DFT in electronic structure calculations has achieved great success for atoms, molecules, and condensed matter since the 1980s.Once the forces acting on the atoms and the model supercell are determined, the atoms’ trajectories can be sampled according to the Boltzmann statistic in the phase space as a function of time at any chosen temperature and pressure by solving the classical Newtonian equation of motion [26], , where m_i is the atomic mass and is the temporal velocity. Useful results obtained from an MD calculation are the temporal positions, velocity, and forces of each atom. From these quantities, both the equilibrium and the thermal transport properties can be calculated [8,26]. Equilibrium thermodynamics variables (e.g., A) of a given ensemble H at temperature T are computed from statistical theory applying the ergodic hypothesis—that is, that the thermodynamic average can be approximate by the time average of the calculated quantity is the Boltzmann constant) [27]:

Thermal transport properties can also be computed from equilibrium MD simulations based on the time correlation formalism pioneered by Green and Kubo [8–12]. In essence, the theory assumes that the dynamics of a transport process is stochastic (Markov’s process) and the perturbation from the thermal equilibrium state is small. The rate of dissipation of a random process (A) is characterized by a decaying time correlation function (if A is a complex dynamical quantity). For example, the diffusion rate (D) is the mean-squared (position) displacement of an atom (diffusion) as a function of time. It can be shown that as the time , the integral of the velocity autocorrelation gives the diffusion coefficient [8]:

By analogy, a time correlation function can be constructed for the viscosity. As mentioned above, viscosity is the response of liquid to elastic deformation. The temporal atom position,  of a particle in a deformed model system is related to the reference,  ，via a linear transformation of the deformation matrix h, ；h is deformed from the reference system by an infinitesimal strain ε ， and I is the identity matrix. At t ¼ 0, the strain is switched off, and the system will experience a δ-function spike in the shear rate, . The timedependence of the shear stress, Pxyð Þt , to the sudden change is the following [8]:

The stress tensor (σ) describes the variation of the total energy under an infinitesimal distortion of the basis vectors  under a strain  .Applying the Hellman– Feynman theorem (Eq. (9)), the stress components can be evaluated:  The steady-state stress (e.g., P_xy) resulting from the steady shear is just the integration over time, and the shear viscosity , and the shear viscosity  is given by the following [8]:

Another useful transport property for geoscience is the thermal conductivity . The thermal conductivity is the accumulated vectorial heat flux ¯the rate of energy momenta transfer per unit area and unit time,  The expression for the thermal conductivity in the form of the heat flux autocorrelation function becomes [8,28]:

Note that the information required for the calculation of the diffusion rate, viscosity, and thermal conductivity is the temporal positions, shear stresses on the model cell, and particle energies. All these quantities can be obtained from MD calculations. The time correlation function is a powerful technique in which many dynamic properties, such as the vibrational density of states (VDOS, infrared (IR) and Raman spectra), can be computed from the atomic trajectories of MD calculations [26]. In the following discussion, we illustrate applications of these methods to a few selected problems in geoscience under extreme conditions.

《3. Illustrative examples》

3. Illustrative examples

In this section, the application of the methodologies presented above to several geologically related problems is illustrated. First, ‘‘static” DFT calculations are used to investigate structural changes in silicon (Si) and germanium oxide glasses [29]. It will be shown that the conventional rigid-sphere concept of definite atomic sizes is no longer applicable at high pressure. In glass structures, at any given pressure, there is no unique Si–O or Ge–O coordination. Therefore, any attempt to relate the coordination number (CN) to the properties of the glass is invalid. This author’s research team has previously described how DFT calculations, in conjunction with modern-day structural search techniques, can be used to predict stable polymorphs of iron (Fe) halides at high pressure [30]. The results show that Fe-rich halides can exist in the Earth’s core; this suggests that heavy halogens may be sequestered into the core, providing a plausible explanation for the depletion of terrestrial volatile halogens. First-principles molecular dynamics (FPMD) calculations were employed to study possible reactions between carbon dioxide (CO₂) and crystalline silica (SiO₂) [31]. It was found that chemical reactions with SiO₂ will occur at relatively low pressure and temperature for CO₂ trapped in the pores of a zeolite. However, if liquid CO₂ is brought into contact with the interface of crystalline quartz and stishovite, no reaction is expected, up to the pressure and temperature conditions of the transition zone. In another study, we found that H2 can react readily with quartz at low pressure and mild temperature to form water [32]. The remaining discussion is devoted to calculations of transport properties. We examined the viscosity of pure calcium carbonate (CaCO₃) melts and the viscosity of CaCO₃ melts in the presence of a small amount of water [33], and found that the viscosity is lower in the hydrated melt than in the neat liquid. We have also discussed a new efficient approximation to compute the thermal conductivity of ordered and disordered solids at high pressure and temperature [34].

《3.1. Structure of SiO2 and GeO2 glasses》

3.1. Structure of SiO₂ and GeO₂ glasses

Since the properties of liquid minerals are difficult to measure at high pressure and temperature, quenched glasses are often used as the first approximation to interpret the structure of the liquid. As mentioned above (supra infra), the concepts of atomic sizes and valence, which were originally introduced by Pauling [4] and O’keeffee and Navrotsky [5], have been a very powerful tool to elucidate the crystal structures of a variety of minerals at ambient pressure. The extension of this powerful concept to disordered glasses is not straightforward. Nevertheless, it is tempting to classify the properties of high-pressure liquids and glasses according to certain observable parameters. An attempt was made recently to correlate the bulk properties of a variety of oxide glasses to the oxygen-packing fraction [35]. This idea does not seem unreasonable, as at high pressure, the bulkier oxygen (O) atoms (oxide anions) in oxide glasses are pressed close to each other, which changes the local environment and should thus affect certain bulk properties such as compressibility and viscosity. The difficulty lies with the assignment of a unique O ionic radius, which is needed to compute the oxygen-packing fraction in a highpressure glass. For this purpose, it was argued that the ionic radius is related to the cation-oxygen CN and the bond length, and that it can be derived from a set of empirical relationships [35]. Since the local geometry (i.e., the CN) is expected to change with pressure, different formulas must be used in different pressure regimes. To apply the method, it is necessary to choose, a priori, the local geometry (i.e., 4-, 5- or 6-coordination) in the oxide glass. This raises the question of the validity of assigning a single CN in oxide glasses at a given pressure. Moreover, there is no rationale to assume a constant atomic radius ratio for a given bond type over a broad pressure range, as the ‘‘size” of individual atom is influenced by the surroundings (i.e., packing) and should be different from that under ambient pressure. To address this problem, the structure and relevant geometric parameters of SiO₂ [36] and germanium dioxide (GeO₂) [29] glass at different pressures were analyzed using theoretical structures obtained from FPMD calculations [26]. The definition of CN in a glass is somewhat arbitrary, as it is dependent on the choice of the cutoff separation. In the study, the distance of the first minimum in the Si–O and Ge–O radial distribution function was chosen as the cutoff separation. The results for GeO₂ are summarized succinctly in Figs. 1 and 2. Similar results were found for the SiO₂ glass. Obviously, except at the lowest pressure, no single unique CN can be ascribed to the glass structure. The structure always consists of a mixture of CNs. Fig. 1(a) shows the proportion of different CNs as a function of pressure in GeO₂ glass. From 0 to 15 GPa, there is a steady decrease in four-fold coordination along with increasing 5- and 6-fold coordination. Between 25–70 GPa, the glass has mixed 5-, 6-, and 7-coordination. The amount of CN greater than 7-coordination increases with a further increase in pressure. These theoretical results show that there is no basis to assume a unique CN for a glass at any pressure. The next question is, how does atomic size change with pressure? Is the change in the atomic size of Ge(Si) and O atoms uniform regardless of the pressure? Will the ratio of the Ge(Si)/O atomic size still be a constant as the Ge(Si)–O distance changes with pressure? Even in crystalline solids, due to multiple different bond lengths, the assignment of the size of an atom in an extended solid is ambiguous. To remove the arbitrariness, Gibbs et al. [37] proposed computing the atomic size from a mathematical analysis of the electron density distribution based on the atom-in-molecule (AIM) quantum theory pioneered by Bader [38]. According to the AIM theory, atomic regions in a molecule or in a solid are defined uniquely from the electron topology. The Bader volume of an atom is the space enclosed by the zero electron flux surface encompassing the atom. This rigorous definition based on quantum mechanical theory removes the arbitrariness and ambiguity of assigning a spherical radius to an otherwise severely distorted atomic density of glassy solids under extreme pressure. For example, this method of partitioning the atom volumes has been shown to be successful in explaining changes in the structure to the electronic property and structural transformation in metallic glasses [39]. A comparison of the Bader volumes for germanium (Ge) and O in GeO₂ glass as a function of pressure is shown in Fig. 1(b). As expected, the ‘‘atom” volume decreases with increasing pressure. The average O Bader volume at 0 GPa is 12.5 ų and decreases to 7.5 ų at 120 GPa, a 40% decrease. In comparison, the change in the Ge Bader volume is much smaller: from 6.2 ų at 0 GPa to 5.2 ų at 160 GPa, a mere 16% decrease. Most significantly, the ratio of the Bader volume changes with pressure and is not a constant. The theoretical results refute the oxygenpacking fraction hypothesis.

《Fig. 1》

Fig. 1. (a) Pressure dependence of the percentage contribution from varying Ge–O coordination of a GeO₂ glass; (b) variation of the Ge (black dots) and O (red dots) Bader volume in GeO₂ glass as a function of pressure.

《Fig. 2》

Fig. 2. (a) Convex-Hull plot for the stability of binary Fe-iodides polymorphs with different stoichiometries at selected pressures; (b) crystal structure of the most stable high-pressure phases revealed from a structure search.

《3.2. The missing halogen paradox》

3.2. The missing halogen paradox

The depletion of volatile elements is related to the nebular volatility. Carbonaceous (CI) chondrites are a group of primitive meteorites possessing an elemental distribution that closely resembles the nebula composition. A comparison of the concentration of elements in Earth with the abundance in CI chondrites provides valuable information on the early stage of Earth’s accretion. Surprisingly, the abundances of heavy halogens (bromine (Br) and iodine (I)) on Earth were found to be an order of magnitude less than their abundances in CI chondrites. This deficit is relevant to Earth’s formation and compositional evolution, and may affect the salinity of the present-day oceans. A hypothesis for this ‘‘missing halogen paradox” is that the halogens may have dissolved in liquid Fe in the early Earth, and subsequently been sequestered into the Earth’s core after crystallization. To constrain the halogens’ concentrations in the core, it is important to know whether stable Fe halide compounds can be formed at the pressures of Earth’s mantle and core. For this purpose, the phase diagrams of Fe halides have been explored using the intelligence-based structural search technique. This approach employs first-principles DFT calculations in conjunction with a structural search technique [40,41], and is a very powerful strategy to locate the most stable crystal structure in an unbiased fashion, given only the knowledge of the constituent atoms at a given thermodynamic condition. This kind of calculation has been shown to provide new insight into crystal types and chemical bonding for materials under pressure conditions not amendable to experiment. The search procedure is as follows: An initial set of structures with different Fe/halogen stoichiometries is generated randomly. The trial structures are first optimized using DFT calculations, and then improved using the particleswarm optimization algorithm [42]. Those structures with better energies are selected as the trial structures for the next cycle of optimization. The procedure is repeated until no new energetically favorable structure is found. In the case of the Fe halides case, this method successfully reproduced the known ambient structure of FeCl₃ (molysite), FeBr₃, FeCl₂ (lawrencite), FeBr₂, and FeI₂. To investigate the structure and stability of the high-pressure polymorphs, calculations were performed on Fe halides with an Fe:halide stoichiometry from 4:1 to 1:4 at the pressures of the transition zone (25 GPa), mantle core boundary (70 GPa), and deep core (135 and 360 GPa). The relative stability of the polymorphs found in the structural search is depicted in a convex-hull plot. The convexhull lines connect the predicted structures of different stoichiometries that are stable against decomposition into their constituent elements or mixtures of compounds at a given pressure. To illustrate the results, the convex-hull plot for the Fe–I system is shown in Fig. 2. A series of stable Fe–I compounds of different stoichiometries was found by the calculations. At low pressure, iodine-rich binary compounds, such as FeI₄ and FeI₃, are found to be most stable. The 1:1 FeI becomes most stable at 125 GPa. At high pressure, Fe-rich compounds are preferred, and Fe₂I is most stable at 320 GPa. The same sequence in the change of the stoichiometry of stable Fe-halides at high pressure is followed in the chemical evolution of the binary Fe–Cl and Fe–Br compounds. At low pressures (≤ 135 GPa for the Fe–Cl and Fe–Br systems and ≤ 70 GPa for the Fe–I system), halogen-rich compounds are more stable. Stable Fe-halides are only formed under Fe-rich conditions at high pressures. Furthermore, the enthalpy differences between the Fe halides and the constituent elements, solid Fe and halogens, were found to be larger at higher pressure. These results suggest that binary Fe-halides, particularly the iodides, can be stabilized in the Earth’s inner core. Significantly, the theoretical results confirm that stable Fe halides can also exist in the form of halogenrich phases at pressures less than 100 GPa—a condition close to that of the early Earth. Therefore, a plausible scenario to describe the depletion of the volatile halides is that Fe halides were formed during the early stage of accretion, and then sequestered into the core as the Earth grew and the pressure increased. The halogenrich Fe-halides that formed at low pressure during the early stage might have undergone chemical evolution with the release of excess halogens to satisfy stoichiometric requirements for chemical stability at high pressure. Remains of Fe-rich halides may still exist in the present-day Earth’s inner core.

《3.3. Chemical behavior of CO2 and SiO2 under mantle conditions》

3.3. Chemical behavior of CO₂ and SiO₂ under mantle conditions

CO₂ and SiO₂, both group IV oxides, are abundant on Earth. Under normal conditions, SiO₂ will not react with CO₂. It is well known that the nature of chemical bonding is strongly influenced by pressure. For example, Si is a semiconductor at room pressure but transforms to a metallic superconductor at 11 GPa! A more dramatic example is sodium, a prototypical model of a free electron metal, which becomes a large band gap insulator under extreme pressure. More relevant to the subject here, when compressed to 40 GPa and heated to 2000 K, solid CO₂ polymerizes into a fully 4-coordinated extended solid, and then transforms into a metallic liquid at 60 GPa and 4000 K. Therefore, at high pressure and temperature, the chemical nature of CO₂ and SiO₂ can be modified to promote chemical reactions. Furthermore, under these conditions, CO₂ is a supercritical fluid and can percolate into porous materials. If reactions do occur, this knowledge may have a profound impact on our understanding of the deep carbon cycle. Santoro et al. [43] first observed the formation of silicon carbonates from a CO₂-loaded microporous SiO₂ (zeolite) sample that was compressed to 18–26 GPa and heated to 600 and 980 K. Spectroscopic evidence unequivocally showed that reactions had occurred. However, no detailed structural information or mechanistic analysis was presented [43]. In view of the potential importance and relevance of these reactions to Earth science, a theoretical study was performed on a model 1:1 CO₂-filled zeolite SSZ-56 system using isothermalisobaric ensemble (constant pressure–constant temperature, NPT) MD calculations. The motivation to use MD lies in its ability to mimic the experiment by applying external pressure at different temperatures to the model system, and to monitor the interactions between CO₂ and SiO₂ before eventual thermodynamic equilibrium is reached (i.e., all reactions have ceased). Analysis of the equilibrium atomic trajectories will provide information on the vibrational spectra of the species formed that can be compared directly with experimental IR and Raman measurements [43].

In the initial thermal equilibrium period at 0 GPa and 300 K, no reaction between the SiO₂ framework and the trapped CO₂ was observed. Upon compression to 26 GPa and above, the SiO₂ framework started to shrink and became distorted, as indicated by the shortening of the crystal cell axes, particularly in the crystallographic b direction. When the compressed structure was heated to 1000 K at 26 GPa, chemical reactions were initiated between CO₂ and the exposed SiO₂ on the channel surface (Figs. 3(a) and (b)). Consequently, unidentate, bidentate, and bridged carbonates (Fig. 3) can be identified from the products. The predicted species are consistent with the suggestions derived from analysis of the experimental vibrational spectra. After 8 ps, no further reaction was observed, although the calculation was continued to 25 ps. The exposed SiO₂ surface is now completely covered by the products from the reaction with CO₂, and the remaining free CO₂ in the cavity can no longer interact with it. The model system has reached thermal equilibrium and the additional equilibrium atomic trajectories can then be used to compute the dynamic properties of the system. Inspection of the final structure of the model reveals that the reaction products are mainly composed of polymeric (–O–(C=O)–O–) chains of 3-coordinated carbonate (CO₃) situated in the large channels of the zeolite, together with some unreacted CO₂ molecules. Surprisingly, the zeolite framework is only slightly distorted after the reactions have completed. This is because the chemical reactions occur only on the surface of the empty channels, and the overall framework is not severely affected, even after pressurization and thermal treatments [43].

《Fig. 3》

Fig. 3. Structure of the SiO₂–CO₂ model after NPT simulation at 1000 K and 26 GPa, showing the presence of (a) unidentate, (b) bidentate, and (c) bridged carbonate fragments.

A distinct advantage of atomistic dynamic simulation is that the reaction mechanism(s) can be unraveled from an examination of the atom trajectories. In this case (Fig. 4), the analysis reveals that the chemical reactions proceeded in three steps. Step i: At high temperature, the CO₂ molecules are thermally excited and execute large-amplitude O–C–O bending motions. Step ii: The reduction in the volume of the model by pressure brings the CO₂ and SiO₂ closer together. The bent CO₂ weakens the C@O bonds and initiates the reaction with the Si of a SiO4 unit forming the surface of the zeolite channels. In a frontier orbital description [44], the lone pair electrons on the O atom of CO₂ are transferred to the empty σ* orbital of SiO₂ (Fig. 5) and concomitantly weaken the Si–O bond. Subsequently, the carbon (C) bonds to the bridging Si–O–Si and forms CO₃. Step iii: Consecutive reactions with other free CO₂ molecules available in the cavity produce the polymeric structure.

《Fig. 4》

Fig. 4. Scheme of the three primary reaction steps.

《Fig. 5》

Fig. 5. Frontier orbital interaction diagram depicting the insertion of the lone pair electrons from CO₂ to the empty Si–O σ* orbital, initiating the chemical reaction between CO₂ and SiO2. Reproduced form Ref. [44] with permission of Academic Press, 1971.

The theoretical results presented above illustrate the versatility of the MD method. In addition to the study of equilibrium properties, MD can be used to study non-equilibrium chemical processes and, from the analysis of atom trajectories, to uncover the atomistic mechanism(s) of the reaction(s).

《3.4. Water formation in Earth’s outer mantle》

3.4. Water formation in Earth’s outer mantle

SiO₂, in the form of SiO₂ and silicates, is one of the most abundant materials on Earth and other terrestrial planets. Quartz, the most common crystalline form, is chemically inert under ambient conditions. However, Shinozaki et al. [45] recently reported the dissolution of SiO₂ in gaseous hydrogen (H₂) under relatively mild pressure and temperature conditions corresponding to those of the Earth’s upper mantle. Furthermore, in situ Raman and IR spectroscopy show evidence of the formation of silane (SiH4) and water (H₂O). The observation of water under this thermodynamic condition is exciting, as it may have an important implication regarding the origin of water on Earth. Moreover, the possibility of the presence of trace water may affect the mantle dynamics that trigger and nucleate deep earthquakes in the continental mantle lithosphere. To understand this novel phenomenon, MD was performed on a model built from a slab of α-quartz containing 1296 SiO₂ with an empty region filled with 3294 H₂ molecules conforming to a gas pressure of 1 GPa at 100 K. It is not feasible at the present time to perform first-principles DFT MD on such a large system. Instead, the reactive force field ReaxFF [46,47], a many-body potential capable of mimicking chemical changes, was used. The ReaxFF calculates bond orders and atom valences dynamically based on the atomic charges determined by the charge-equilibration (QEq) algorithm. The interatomic interactions are adjusted on the fly in response to chemical and structural changes during the process. The model system was slowly brought to the experimental conditions at 1700 K and 2 GPa. Initially, hydrogen molecules were found to diffuse into the subsurface region of the quartz slab (Fig. 6). As the channels in quartz are relatively spacious, the H₂ molecules could penetrate deep into the slab; eventually, the trapped hydrogen led to a semi-uniform concentration profile (Figs. 6(b)–(d)) throughout the slab. At this stage, the hydrogen maintained its molecular structure but the H–H bond was lengthened from a molecular bond length of 0.76 Å to 1.0–1.1 Å (Figs. 6(b)–(d)). Subsequently, the weakened hydrogen molecules dissociated (Fig. 6(e)) into atomic hydrogen and interacted chemically with the O atoms of the quartz. The reaction weakened the Si–O bond and led to the formation of water molecules inside the slab. This process continued until a pocket of water was confined inside the slab (Fig. 6(e)). The reacted quartz was transformed into amorphous SiO_xH_y. In summary, the SiO₂/H₂ dissolution reaction proceeded in three distinctive stages (Fig. 7). First, H₂ penetrated into the porous channels of quartz. Next, the H2 interacted with one of the corner-sharing O atoms and with a Si atom. The reaction led to the dissociation of H₂, forming an O–H and a Si–H bond. A second H2 then reacted with the Si of a neighboring SiO₂ molecule in a similar fashion, leading to the formation of H₂O and the breaking of the SiO₂ tetrahedral linkage. The MD calculations successfully reproduced all the features observed in the experiments. This study supports the suggestion that a reaction between SiO₂ (and perhaps silicates) and H2 in the mantle may be a possible indigenous source for the production of water on Earth.

《Fig. 6》

Fig. 6. Snapshots of the SiO₂ slab model at different stages of the performed MD simulation: (a) initial structure; (b) initial model equilibrated to 100 K and 1 GPa; (c) structure at 580 K and 1.3 GPa; (d) structure at 1220 K and 1.7 GPa; (e) final structure at 1700 K and 2 GPa. All structures are oriented to the (100) direction. Si atoms are shown in green, O in red, and H in light blue.

《Fig. 7》

Fig. 7. Water formation from the reaction of hydrogen with quartz. From left to right: A H₂ molecule first adsorbs to a Si–O pair in the crystal lattice, dissociates to atomic hydrogen, and then forms a water molecule with the remaining hydrogen atoms attached to the Si atom.

Analysis of the structure of the confined water revealed a surprise. The O–O radial distribution function, , (Fig. 8(a)) in the water region is remarkably similar to that in high-pressure bulk water with a density of 1.95 g·cm^-3 reported in a previous FPMD calculation [48]. In particular, the absence of a clear first coordination shell—that is, the first minimum in the O–H radial distribution at distance r = 1.25 Å is nonzero—is an indication for the dissociation of the water molecules [48]. Visual inspection of the water structure and analysis of the chemical species present in the water region confirmed this suggestion. The molecular species found are summarized graphically in Fig. 8(b). The presence of OH^– , H₂O, H₃O⁺ , and H₄O²⁺ in the water region can be clearly seen. These theoretical results are unexpected and may have some significant consequences. The presence of a large amount of ionic species, such as H₃O⁺ , H₄O₂⁺ , and OH, from the dissociation of water under the anomalously high confining pressure may alter not only its chemical properties, such as enhanced reactivity, but also its physical properties (e.g., giving it a higher thermal diffusivity and electrical conductivity than its pure counterpart under similar conditions); in turn, this would exert a profound effect on diverse geological processes, from mantle metasomatism to partial melting and magmatism.

《Fig. 8》

Fig. 8. (a) The O–O, O–H and H–H radial distribution functions of the water layer confined in the SiO2 slab; (b) chemical species identified (OH^-, H₂O, H₃O⁺ , and H₄O²⁺ compounds) in the confined water layer. ID: identified.

《3.5. Viscosity of CaCO3 melts》

3.5. Viscosity of CaCO₃ melts

CaCO₃ is a common chemical component of rock minerals. In nature, it exists in its pure form as calcite and aragonite. A large amount of carbonate is also found in limestone, chalk, and marble. Therefore, carbonate minerals are potential carbon reservoirs. Knowledge on the mobility of carbonate-rich melts in the Earth’s mantle is important in understanding the deep carbon cycle and related geochemical and geophysical phenomena. Viscosity of mineral melts is a key parameter in the quantification of volcanic processes and magma flow. However, it is challenging to measure viscosity at high pressure and high temperature. Therefore, few experimental viscosities are known under the conditions of the transition zone or deeper, which is an active region of deep earthquakes. Recently, ultrafast radiographic imaging using synchrotron X-rays was used to measure the viscosity of CaCO₃ melt to 6.2 GPa at 2000–2500 K. It was found that the viscosity of carbonate melts is comparable to that of liquid water under ambient conditions, and is about one order of magnitude less than that of calcium silicate melts [49]. To investigate the origin of this surprisingly low viscosity and to obtain information relevant to the Earth’s transition zone, FPMD calculations were performed on a CaCO₃ melt up to 35 GPa at 2000 and 3000 K. The viscosity was computed from the Green– Kubo time autocorrelation function P_xy—that is, the off-diagonal component of the instantaneous stress tensor (Eq. (13)) from the MD calculations in the constant volume–constant temperature (NVT) canonical ensemble. One of the inconveniences of the correlation function method is the persistent small oscillations occur in the correlation function over long time. Therefore, an accurate estimation of the viscosity requires equilibrium MD simulation greater than 50 ps. Even longer simulation time is required when the melt is close to the liquid–glass boundary. The shear stress time autocorrelation function of a CaCO₃ melt calculated at 32.6 GPa and 3000 K, as shown in Fig. 9(a), clearly shows the long oscillatory tail. For this system, the viscosity is estimated to be (14 ± 1) mPas. Similar viscosity calculations were performed at selected pressures at 2000 and 3000 K. The results are summarized in Fig. 9(b). At 4.1 GPa and 2000 K, the theoretical viscosity of 5.8 mPa·s is in reasonable agreement with the experimental value of about 6 mPas measured at 2063 K and 6.2 GPa. Above 8 GPa at 2000 K, the CaCO₃ melt becomes a viscous liquid with a very low self-diffusion coefficient. It is no longer practical to compute the viscosity, as a very long simulation time is needed. On the other hand, the CaCO₃ melt at 3000 K is fluid-like, and the viscosities can be computed up to 32.6 GPa. At 0.6 GPa and 3000 K, the calculated viscosity is 2 mPa·s. This value is an order of magnitude lower than that of most silicate melts under similar thermodynamic conditions and, incidentally, is close to the viscosity of liquid water at ambient conditions. This comparison, however, is fortuitous. Analysis of the structure of the melts shows no hint of polymerization among the carbonate anions. At low pressure, the CaCO₃ melt behaves like an ideal molecular liquid. In fact, the directly computed viscosities are in reasonable agreement with estimates from the Stoke–Einstein equation using the computed diffusion coefficients. Similar to the experimental trend [49], the calculated viscosity of the melt at 3000 K only increases slowly at less than 10 GPa. A potentially significant prediction is that the viscosity increases abruptly at pressures greater than 10 GPa, indicating a change in the short-range order in the liquid. This prediction is confirmed from visual examination of the temporal atom configurations, which show that there are short-lived interactions between the carbonate anions and the Ca²⁺ ions at pressures greater than 10 GPa. This interaction can be quantified from the calculations of the velocity cross-correlation (VCC) between calcium (Ca) and the C and O atoms of the carbonates ( and ). The calculated Ca–O and Ca–C VCCs at 3000 K at several pressures are compared in Figs. 10(a) and (b). At pressures up to 7.1 GPa, the rapid and smooth decay of the correlation function at short simulation time and the absence of weak oscillations at long time in the VCC profiles are indicative of a normal diffusive fluid behavior. However, at pressures greater than 11.2 GPa, oscillatory features between 100 and 350 fs start to appear in the VCC. These features are due to correlated Ca...O and Ca...C motions. These oscillatory features are expected to develop a distinct vibrational band in the power spectrum. The single-particle VDOS can be computed from the Fourier transform (power spectrum) of the velocity time autocorrelation function,  dt. Indeed, a new band at 400–500 cm^-1 due to the Ca...CO₃ vibrations emerges in the Ca projected VDOS (Fig. 10(c)) at pressures greater than 10 GPa. The vibrational energy matches nicely with the oscillations in the VCC, with a period of about 71 fs or 470 cm^-1 . It can be concluded without any ambiguity that at high pressure (> 4 GPa at 2000 K and > 10 GPa at 3000 K), the Ca²⁺ and CO^2-₃ are compressed against each other and their motions are correlated. As a result, the atom diffusion rate decreases and the viscosity of the melt increases suddenly. In addition, the increase in the viscosity affects the heat transport (thermal conductivity). These effects may have a significant effect on the mobility of the melt in the Earth’s upper mantle. We have also studied the effect of water impurity on the viscosity of CaCO₃ melt. NVT MD simulation was performed on a model system with 15% mole-percent of water in CaCO₃ melt at 30 GPa and 3000 K. The calculated viscosity was 10 mPa·s, which is 10% lower than that of pure CaCO₃ melt under the same thermodynamic conditions.

《Fig. 9》

Fig. 9. (a) Shear stress autocorrelation function (black solid line), cumulative integral (blue solid line), and viscosity (red dashed line) of CaCO₃ at 3000 K and 32.6 GPa; (b) pressure dependence of the viscosity of molten CaCO₃ at 2000 and 3000 K, and of 15% water-CaCO₃ solution at 30 GPa. exp.: experimental value.

《Fig. 10》

Fig. 10. Pressure-dependent VCC functions of (a) Ca–C and (b) Ca–O in CaCO₃ melts. (c) Ca projected velocity density of states at 3000 K.

《3.6. Thermal conductivity of periclase (MgO)》

3.6. Thermal conductivity of periclase (MgO)

Thermal conductivity is an important thermophysical property of a material. Knowledge on the thermal conductivity of minerals at high pressure and high temperature is relevant to our understanding of the rate of heat loss in the core and the evolution of the Earth. Direct measurement of the thermal conductivity under extreme conditions is a formidable task. In recent years, many efforts have been made to develop efficient computational schemes to calculate the thermal conductivity from first-principles electronic calculations [50,51]. For systems that can be adequately described by many-body pairwise additive potentials, the thermal conductivity is most conveniently calculated from the heat flux time correlation function using equilibrium MD employing the Green–Kubo formalism (Eq. (14)) [28]. For the pairwise additive potential, the kinetic and potential energy term of the heat flux can be partitioned to each individual atom contribution. This method, however, cannot be extended to first-principles electronic structure-based MD, as there is no valid way to define the potential energy of a single atom in a quantum mechanical system. Recently [34], it was shown that if the convection term of the thermal conductivity is ignored, such as in solids (crystalline, disordered, and amorphous), the thermal conductivity can be estimated from the sum of the energy momentum of individual atoms (see Eq. (14) [28]). The reasoning is as follows: The Fourier law of heat conduction states that the thermal conductivity of a material k is the ratio of the heat flux q transported to the thermal gradient ():

The heat current J is the sum of the rate of the energy of each particle ε_i transport through a surface of cross-sectional area A (normal to the heat current) per unit time:

The energy associated with each particle is the sum of the kinetic and potential (u) energy:

Substituting into Eq. (16), the heat flux can be written as follows:

The first term in Eq. (17) is the thermal conductivity due to convection, and the second term is the thermal conductivity due to conduction. In a solid, the convection term is zero and can be neglected. The energy momenta R for a solid can be defined as follows:

Note that the heat current J is the time derivative of R—that is, . The energy momentum for a particle k is:

Instead of using the Green–Kubo formula [8] to compute the thermal conductivity k by integrating the time autocorrelation of the heat current J, this quantity can be computed using an equivalent Einstein relationship for the diffusion of the energy momenta—that is, the mean-squared displacement of the energy momenta [52]:

In principle, the instantaneous energy momenta (Eq. (20)) can be computed from the knowledge of the temporal positions, velocities, and net forces on each particle. These quantities are available in an MD calculation.

The simplified approach (Eqs. (20) and (21)) for the prediction of the thermal conductivity has been validated by investigations on a diverse class of ordered and disordered solids [34]. The mineral periclase (MgO) is a common component in the Earth’s lower mantle. In view of its importance, its simple cubic structure, and the availability of accurate thermal conductivity data on the single crystal [53] at high temperature and on powder samples of lesser quality at high pressure [54,55], periclase is often chosen as a test case to validate new computational schemes. The thermal conductivity of MgO at 0 GPa and 100, 300, 600, 900, and 1200 K was computed using the above-described method with NVT MD. The numerical results obtained using a 4 × 4 × 4 supercell are compared with experimental data and previous theoretical anharmonic lattice dynamics calculations in Fig. 11(a). The thermal conductivity of MgO at 300 K computed using the new method, that is, (70.3 ± 8.9) W·K^-1·m^-1 , is consistent with the rather scattered measured data, which ranges from 36 to 70 W·K^-1·m^-1 . This value is also in agreement with other theoretical estimates of 65– 111 W·K^-1·m^-1 . It should be noted that the imperfections in the density functional used in the calculations and/or the natural distribution of magnesium (Mg) and O isotopes in the measured samples have been ignored. At a given pressure, the thermal conductivity should decrease with temperature due to predominant Umklapp phonon–phonon scattering processes [56]. This trend was reproduced by the calculations, with the thermal conductivity decreasing from 153 W·K^-1·m^-1 at 100 K to 13.7 W·K^-1·m^-1 at 1200 K. The absolute magnitudes of the calculated thermal conductivity are also in good agreement with experiments over a broad temperature range. An added advantage of the MD approach is that the statistical error associated with the calculated thermal conductivity can be estimated as the standard deviation from the average value computed using time segments of the MD trajectory at different time origins but with the same time length. As shown in Fig. 11(a), the statistical errors become smaller at higher temperature. This is an indication that the sampling of the configuration space becomes more efficient at high temperature. Furthermore, for a large enough model, the MD result is expected to be more accurate at high temperature, as the anharmonic large-amplitude atomic motions that enhance the Umklapp scattering processes [56] are explicitly included. This is contrary to the anharmonic lattice dynamics method, in which the effect of temperature only changes the phonon population via the Bose–Einstein statistics, without corrections to the anharmonic force constants and the harmonic phonon frequencies, both of which are computed at 0 K. Indeed, compared with experimental and MD results, the thermal conductivities obtained from anharmonic lattice dynamics calculations decrease much faster at high temperatures [52,53]. Near the melting point, higher order force constants will become important, and the use of the harmonic frequencies is no longer valid. In principle, all orders of the anharmonic effects are included in the MD method. To test the effect of long-wavelength phonons and long-range interactions—that is, the convergence of the predicted thermal conductivity versus system size—another calculation was made at 300 K using a larger 5×5 ×5 supercell model. The computed thermal conductivity of (73.3 ± 9.2) W·K^-1·m^-1 for the larger model falls within the error estimate of the smaller 4 × 4 × 4 model of (70.3 ± 8.9) W·K^-1·m^-1 (vide infra). Therefore, at least the system size dependence in the MgO test case is not too severe.

《Fig. 11》

Fig. 11. Comparison of experimental and theoretical thermal conductivities at (a) high temperature and ambient pressure and (b) high pressure and 300 K. (a) Reproduced from Ref. [53] with permission of National Academy of Sciences, 2010; (b) reproduced from Ref. [54] with permission of Springer Nature Publishing AG,
2013 and Ref. [55] with permission of American Geophysical Union, 2014.

Directin situ experimental determination of the absolute thermal conductivity of a material at high pressure and high temperature is not feasible. Under these extreme conditions, most information was obtained indirectly from the measurement of the thermal diffusivities from powder samples, and then converted into thermal conductivities by considering the heat capacity at constant pressure. The heat capacity is often estimated from the extrapolation of the equation of states in combination with other known elastic properties. Recently, attempts have been made to measure absolute thermal conductivity in a diamond anvil cell using transient thermoreflectance spectroscopy [54]. A potential shortcoming of this method is that (semi)-empirical corrections must be applied to correct for the contribution from the sample environment. Using this method, the thermal conductivity of a single-crystal MgO was reported up to 60 GPa at room temperature [54]. Using the same 4 × 4 × 4 supercell model, the thermal conductivity at 3300 K and at selected pressures was computed. The theoretical results are compared to the experimental values in Fig. 11(b). The theory correctly predicted the observed trend in which the thermal conductivity increases with increasing pressure. At the highest pressure of 137 GPa, the calculated thermal conductivity is (357 ± 41) W·K^-1·m^-1 , which is almost five times larger than the ambient pressure value of (70.3 ± 8.9) W·K^-1·m^-1 . In general, the predicted thermal conductivity is slightly higher than the observed value. For example, the calculated values at 47 and 65 GPa of (178 ± 25) W·K^-1·m^-1 and (211 ± 30) W·K^-1·m^-1 are larger than the measured values of 140 and 162 W·K^-1·m^-1 , respectively. Taking into consideration the fairly large experimental errors at high pressure (±20 W·K^-1·m^-1 at 47 GPa and ± 30 W·K^-1·m^-1 at 60 GPa), the theoretical thermal conductivities calculated using the new MD method are not unreasonable. The system size dependency is also small, as the thermal conductivity of (218 ± 31) W·K^-1·m^-1 that was calculated using the larger 5 × 5 × 5 supercell at 65 GPa is within the error bound of (211 ± 39) W·K^-1·m^-1 computed from the smaller 4 × 4 × 4 model. In contrast to the situation at high temperature, the statistical errors become larger at higher pressure. This is because when the atoms are compressed against each other, thermal sampling of the configuration space is restricted. A longer MD run is required to achieve higher numerical accuracy.

Although the pressure trend of the thermal conductivity is correctly reproduced by the calculations, the absolute values are consistently overestimated, and the deviation from the experimental results increases with increasing pressure [55]. However, the MD results are very similar to earlier lattice anharmonic dynamics calculations [53]. The cause of the discrepancy is unlikely to be a deficiency in the theory. Since the experimentally measured quantity was thermal diffusivity [55], information on the heat capacity at constant pressure, C_p, is needed to convert this quantity into thermal conductivity. Although C_p at high pressure is not a measureable quantity, it can be estimated from the known experimental equation of states, using information on elastic properties at low pressures [57]. It is possible that the extrapolation of this quantity to high pressure may not be valid. This speculation is confirmed by a comparison of the empirical heat capacity at C_p of MgO obtained from the equation of states [54] with the results of quasi-harmonic (QHA) lattice dynamics calculations at several pressures (Fig. 12). It is shown that the estimated C_p and the calculated C_p at ambient pressure are in agreement, but that the C_p estimated from the equation of states (‘‘empirical”) is lower than the values computed by the quasianharmonic theory at higher pressures. The difference becomes even larger at higher pressure. Therefore, the ‘‘experimental” thermal conductivity is consistently lower than the theoretical value. If the ‘‘experimental” thermal conductivities were scaled with the appropriate (QHA/empirical) C_p ratios (inset, Fig. 12), the absolute agreements with the theoretical results would be greatly improved (Fig. 11(b)).

《Fig. 12》

Fig. 12. Comparison of the heat capacity at C_p of MgO from theoretical quasiharmonic approximation (QHA) (red line) and estimated from the experimental equation of states (Exp. EOS); (black line) and a second-order polynomial fit (blue line) reported in the Supplementary Materials of Ref. [54].

The computational scheme neglects the contribution from convection and is not applicable to liquids. In a full quantum mechanical theory, it is not possible to partition the energy momenta into individual atomic components. The very reasonable agreement between the predicted and observed thermal conductivity for a large class of solids [34] may be rationalized by the following intuitive argument: Starting from the equilibrium potential energy surface, the temporal changes in the total energy of the system can be approximated as a sum of the perturbations on individual atoms. The energy change associated with each atom is simply the displacement times the net force acting on that atom, with the forces computed via the Hellmann Feynman theorem. The total of the energy momenta is the sum of the individual atom contributions. For a solid, Eq. (21) is a reasonable first approximation. Despite the inexactness, empirical evidence supports the accuracy of the MD method and indicates that this method can be applied to high temperatures, at which high-order anharmonic effects cannot be neglected.

《4. Outlook》

4. Outlook

Theoretical modeling has been and will continue to be an indispensable tool for high-pressure/high-temperature research. In this mini-review, only a small and restricted fraction of possible applications were illustrated. Using the limited examples presented here, it was demonstrated that reliable thermophysical transport properties such as viscosity and thermal conductivity can be computed using FPMD. The advantage of the MD approach is that the temperature effect and the associate anharmonic atomic vibrations are implicitly included in the algorithm. In principle, DFT lattice dynamics calculations correspond to the athermal (i.e., 0 K) condition. Correction of the thermal effect using perturbative treatments, such as using the quasi-harmonic theory, [58] is needed in order to make comparisons between theoretical results and experimental measurements. These corrections are valid only if the temperature does not far exceed the Debye temperature of the solid. Unfortunately, many interesting anomalies in minerals occur near the melting point, which is a temperature much higher than the Debye temperature. Significant existing challenges in experiments under extreme temperature and pressure conditions prohibit the direct measurement of viscosity and thermal conductivity. Therefore, FPMD may be the preferred computational method. For example, a potential application of this method is in the study of the shear weakening close to the melting of hexagonal close-packed Fe (hcp-Fe) under inner-core conditions [59]. Although not mentioned here, the elastic moduli can be computed on the fly from the fluctuations of the strains in the isobaric– isoenthalpy (p, H, N) ensemble [60–62]. In this brief review, it was also shown that FPMD can be used to explore novel chemical processes and to reveal atomistic mechanisms that may not be amendable by current experiment under extreme conditions.

The FPMD calculations illustrated here require substantial computational resources. To extract reliable statistics on the properties computed using the time correlation function method, lengthy simulations are needed. The time steps required for data analysis after thermodynamic equilibrium may be reduced using the corresponding but equivalent Einstein diffusion equations, as was done for the calculation of thermal conductivity. In the examples described here, FPMD simulation of several hundred atoms at one thermodynamic data point using current computer technology would consume 1–2 months of computing time on a massively parallel computer employing 100–200 cores. This author is optimistic that the MD approach will be widely used in the near future, in view of rapid advancements in computational capacity and capability, as well as active developments in more reliable and efficient algorithms for large-scale electronic structure calculations.

《Acknowledgements》

Acknowledgements

The author wishes to acknowledge his collaborators, Drs. N.J. English, Y. Pan, and T. Iitaka, and previous postdocs and students, Drs. X. Du, M. Wu, and X. Yong, for their contributions to the studies reported here. He also thanks WestGrid (Canada), RIKEN (Japan), and HPSTAR (China) for the allocation of computational resources.