Acceleration of Atomistic Simulations of Lower Mantle Silicate Melts by Combination of Monte Carlo and Molecular Dynamics

Silicate phases are the main constituents of rocky planets and hence their properties determine to a major extent the differentiation and geodynamics of such planets. However, the study of silicate melts in the lower mantle is limited by the difficulty and cost of laboratory experiments at simultaneous high pressure and temperature. Molecular dynamics (MD) simulations offer a powerful alternative tool to study melt properties over a wide range of composition and imposed conditions. Increasing computer power has led to several MD-based studies of silicate melts under pressure in the past decade.

To obtain the properties of a system by MD, it is critical to efficiently achieve (and demonstrate the achievement of) equilibrium states. However, multicomponent silicate melts are very viscous, have complex intermediate-range order phenomena, and can easily demonstrate glassy or metastable behavior on simulation timescales at lower mantle pressure and temperature. As a result, the calculations are expensive, the number of atoms and length of runs are limited, and the determination of equilibrium is sometimes too subtle to capture by traditional criteria.

To ensure data quality and reduce simulation cost, we need accelerated simulation techniques and more robust equilibration criteria. In this study, we combine traditional MD simulations with Monte Carlo (MC) swaps of cation positions to explore the configuration space more quickly than does ordinary diffusion in a viscous melt system. The goal of this preliminary study using empirical force-field MD is to identify the schedule of MD and MC steps that minimizes simulation wall time. We expect this technique to save major computational cost when ported to ab initio MD (AIMD) simulations of complex natural systems, reducing the number of calculations of electron density distribution. We also propose sensitive new criteria for achievement of equilibrium melt states based on the local structural environments around the atoms and on the success rates of the proposed MC cation swaps. Compared to traditional equilibration criteria (fluctuation of thermodynamic quantities, mean-squared displacements, etc.), we conclude that the new criteria can reveal more subtle equilibration properties such as the ordering of 2^nd nearest neighbors and radial charge distributions.