Global structure searches under varying temperatures and pressures using polynomial machine learning potentials: A case study on silicon

This study proposes a robust methodology utilizing polynomial machine learning potentials to systematically enumerate crystal structures and evaluate phase stability for elemental silicon under high-pressure (up to 100 GPa) and finite-temperature (up to 1000 K) conditions.

The Gist

Imagine you are trying to find the perfect arrangement of furniture in a room, but the room is constantly changing size (pressure) and the temperature is fluctuating wildly. You want to find the most comfortable, stable layout for every possible scenario. Now, imagine doing this not for a living room, but for the atoms inside a piece of silicon, and you have to test millions of layouts to find the absolute best one.

That is essentially what this paper does, but with a high-tech twist. Here is the story of how the researchers solved this puzzle.

The Problem: The "Needle in a Haystack"

Silicon is the stuff computer chips are made of. We know it well at room temperature and normal pressure (it's the diamond-like structure we use in electronics). But if you squeeze silicon really hard (high pressure) or heat it up, it changes shape into different "polymorphs" (like a chameleon changing colors).

Scientists want to predict exactly which shape silicon will take under any combination of pressure and heat. To do this, they usually use Density Functional Theory (DFT). Think of DFT as a super-accurate, super-expensive GPS. It tells you exactly where you are, but it takes a long time to calculate a single step. To map out the entire landscape of silicon's shapes, you would need to take billions of steps. Doing this with the "super-accurate GPS" would take a supercomputer years to finish.

The Solution: The "Smart Shortcut"

The researchers developed a Machine Learning Potential (MLP). Think of this as a student who has studied the GPS maps for a while.

• The Teacher (DFT): Knows everything perfectly but is slow.
• The Student (MLP): Has learned the patterns. It's not quite as perfect as the teacher, but it is millions of times faster.

However, there was a catch. Previous "students" (MLPs) were good at learning the room at normal pressure, but when the room got squeezed (high pressure), they got confused and gave bad directions. They couldn't handle the stress.

The Innovation: Training the Ultimate Student

The authors created a new, super-smart student using a Polynomial Machine Learning Potential. Here is how they trained it:

1. The Curriculum (The Dataset): Instead of just showing the student the room at normal pressure, they showed it the room being squeezed, distorted, and heated up. They generated thousands of "what-if" scenarios where the silicon atoms were pushed to their limits.
2. The Hybrid Approach: They didn't just use one model. They built a Hybrid MLP, which is like having two experts working together. One expert looks at the immediate neighbors of an atom (short-range), and the other looks further out (long-range). They combine their opinions to get a perfect prediction.
3. The Iterative Loop (The "Study Hall"):
   • The student tries to find the best furniture arrangement (structure search).
   • Sometimes the student makes a mistake and picks a weird arrangement.
   • The researchers catch the mistake, ask the slow "Teacher" (DFT) for the real answer, and feed that new information back to the student.
   • The student learns from the mistake and gets better. They repeat this cycle until the student is nearly as good as the teacher but still super fast.

The Grand Tour: Mapping the Silicon Universe

Once they had this super-trained student, they sent it on a massive journey:

• The Search: They let the student explore millions of random atomic arrangements under pressures up to 100 GPa (that's about 1 million times the atmospheric pressure at sea level!) and temperatures up to 1000 K.
• The Filter: The student found the "local minimums" (the stable spots).
• The Vibration Check (SSCHA): Atoms aren't static; they vibrate. At high temperatures, these vibrations can make a structure unstable. The researchers used a method called SSCHA (Stochastic Self-Consistent Harmonic Approximation) to simulate these vibrations. Think of this as checking if a chair is stable not just when you sit still, but when you wiggle around in it.

The Results: A New Map

The result is a Pressure-Temperature Phase Diagram. It's a map that tells you:

• "If you are at 50 GPa and 500 K, silicon will be in this specific shape."
• "If you heat it up to 800 K, it will shift to that other shape."

They confirmed almost all the shapes scientists had already discovered in experiments. But they also found something new: a specific shape (called the α-La type) that becomes the most stable at very high pressures and temperatures, filling in a gap in our knowledge.

Why This Matters

This paper isn't just about silicon. It's about how we do science.

• Before: We had to choose between being accurate (but slow) or fast (but inaccurate).
• Now: This study shows we can have both. By training a smart AI model on diverse, high-pressure data, we can explore the universe of materials much faster than ever before.

In a nutshell: The researchers built a "super-fast, super-smart AI assistant" that learned how silicon behaves under extreme stress. They used this assistant to map out the entire landscape of silicon's shapes, creating a reliable guide for future materials science that was previously impossible to draw.

Technical Summary

1. Problem Statement

Predicting crystal structures and phase stability under extreme conditions (high pressure and finite temperature) is a significant challenge in materials science.

• Limitations of DFT: While Density Functional Theory (DFT) is the gold standard for accuracy, global structure searches (GSS) and finite-temperature stability analyses (requiring free energy calculations) are computationally prohibitive when applied to vast configurational spaces and numerous thermodynamic states.
• Limitations of Existing MLPs: Existing Machine Learning Potentials (MLPs) often lack the necessary accuracy across a broad range of pressures (0–100 GPa) and structural complexities (low-symmetry local minima) to reliably perform GSS. Furthermore, applying stochastic self-consistent harmonic approximation (SSCHA) calculations to a wide set of local minimum structures is computationally demanding, especially for low-symmetry systems where force constant (FC) optimization is difficult.
• The Gap: There is a need for a robust, integrated framework that combines high-accuracy MLPs with iterative search strategies and efficient anharmonic free energy calculations to generate reliable pressure-temperature ($P-T$) phase diagrams.

2. Methodology

The authors propose a comprehensive, iterative framework consisting of four main stages:

A. Dataset Construction and MLP Development

• Data Generation: The study starts with a base dataset of ~12,000 structures (Dataset 1) optimized at zero pressure. To address high-pressure accuracy, a new dataset (Dataset 2) is created by optimizing 86 prototype structures at 25, 50, 75, and 100 GPa, followed by random lattice expansions, distortions, and atomic displacements.
• Polynomial MLP Formulation: The study utilizes Polynomial Machine Learning Potentials (MLPs) based on polynomial rotational invariants derived from radial and spherical harmonic functions.
  • Hybrid Models: A "Hybrid Polynomial MLP" is introduced, combining two distinct MLP models with different input parameters (e.g., cutoff radii of 5.0 Å and 6.0 Å) to capture both short-range and long-range interactions effectively.
  • Training: The models are trained using weighted linear ridge regression on energy, force, and stress tensor components from DFT calculations (VASP code, PBE functional).

B. Iterative Global Structure Search (GSS)

• Random Structure Search (RSS): The authors employ RSS to generate random initial structures, which are then locally optimized using the MLP.
• Iterative Update: To ensure reliability, an iterative loop is implemented:
  1. Perform GSS with the current MLP.
  2. Identify local minimum structures.
  3. Perform single-point DFT calculations on these minima.
  4. Add these DFT results to the training dataset and retrain the MLP.
  5. Repeat until convergence.
• Scale: The search covers pressures from 0 to 100 GPa and identifies over 28,000 local minimum structures.

C. Finite-Temperature Phase Stability (SSCHA)

• Stochastic Self-Consistent Harmonic Approximation (SSCHA): To evaluate phase stability at finite temperatures (up to 1000 K), the authors perform SSCHA calculations. This method captures anharmonic vibrational effects, which are crucial for predicting dynamical stability in structures that appear unstable under the harmonic approximation.
• Two-Stage Workflow:
  1. General-Purpose MLP: Used to screen ~80 low-enthalpy structures and perform initial SSCHA calculations to identify candidates for high-accuracy analysis.
  2. Specialized MLPs: For the top 27 candidates (those with Gibbs free energies within 5 meV/atom of the global minimum), specialized MLPs are trained on-the-fly using DFT datasets specific to each structure. This ensures high accuracy for the specific local environments and large atomic displacements required by SSCHA.

D. Phase Diagram Construction

• Helmholtz free energies calculated via SSCHA are fitted to the Rose–Vinet Equation of State (EOS) to determine Gibbs free energies across the $P-T$ space.
• The resulting phase diagram is compared against experimental data and previous theoretical studies.

3. Key Contributions

• Robust MLP Framework: Development of a hybrid polynomial MLP capable of high-accuracy predictions (RMSE ~2.8 meV/atom for energy) across a wide pressure range (0–100 GPa) and diverse structural motifs, including low-symmetry local minima.
• Iterative Search Strategy: Demonstration of an iterative RSS-MLP-DFT loop that systematically improves MLP accuracy and ensures the discovery of global minimum structures that might be missed by static potentials.
• Scalable SSCHA Implementation: A novel procedure for performing SSCHA calculations on a large set of complex, low-symmetry structures by combining general-purpose screening with specialized, on-the-fly MLP training. This overcomes the computational bottleneck of traditional DFT-based SSCHA.
• Comprehensive Silicon Phase Diagram: The generation of a highly accurate $P-T$ phase diagram for elemental silicon (0–100 GPa, 0–1000 K) that resolves discrepancies in previous theoretical works.

4. Key Results

• Structural Discovery: The search successfully identified all experimentally observed silicon phases (Si-I, Si-II, Si-V, Si-VI, Si-VII, Si-X) and the metastable Si-XI phase. It also identified the $\alpha$-La type structure as a global minimum in the 83.8–87.7 GPa range at 0 K.
• Accuracy Validation:
  • The polynomial MLP outperformed traditional potentials (MEAM, Tersoff, SNAP) and even the Gaussian Approximation Potential (GAP) in predicting energies of local minima, with a Mean Absolute Error (MAE) of only 19 meV/atom compared to >260 meV/atom for GAP.
  • The MLP accurately reproduced DFT energy-volume curves and phonon density of states for high-pressure phases.
• Phase Stability Insights:
  • Anharmonic Effects: SSCHA calculations revealed that anharmonic effects significantly expand the stability fields of certain phases (e.g., Si-XI, Si-VII, $\alpha$-La) at higher temperatures compared to Quasi-Harmonic Approximation (QHA) results.
  • Comparison with Literature: The study's phase diagram contradicts previous theoretical works (Li et al., Paul et al.) that predicted the $\alpha$-La phase to be stable at lower pressures (<60 GPa). The current study places $\alpha$-La stability at 85–100 GPa (at 1000 K), which aligns better with experimental dynamic-compression data.
  • Experimental Agreement: The predicted stability ranges for the Simple Hexagonal (SH) and Hexagonal Close-Packed (HCP) phases at 1000 K show strong agreement with experimental shock-compression data.

5. Significance

• Methodological Advancement: This work establishes a generalizable, robust framework for predicting crystal structures and phase stability under extreme conditions. The combination of polynomial MLPs, iterative active learning, and specialized SSCHA workflows can be applied to other complex elemental and compound systems.
• Resolution of Discrepancies: By providing a more accurate theoretical phase diagram for silicon, the study clarifies conflicting results in the literature and offers a benchmark for future high-pressure physics research.
• Computational Efficiency: The approach demonstrates that high-accuracy anharmonic free energy calculations for hundreds of structures are feasible using MLPs, reducing the computational cost by orders of magnitude compared to pure DFT approaches while maintaining DFT-level accuracy.
• Material Design: The methodology provides a powerful tool for discovering new high-pressure polymorphs and understanding phase transitions in materials relevant to geophysics, planetary science, and advanced manufacturing.