Incorporating Gibbs free energy into interatomic potential fitting

This paper presents a novel method for fitting interatomic potentials by incorporating high-temperature Gibbs free energy data via Hamiltonian thermodynamic integration, a strategy validated on Ni and Fe-O systems that effectively complements conventional fitting techniques for structural and elastic properties.

The Gist

Here is an explanation of the paper, translated into everyday language with some creative analogies.

The Big Picture: Tuning the "Virtual Universe"

Imagine you are a video game developer trying to create a realistic simulation of a world made of atoms. To make the game run fast enough for a computer to handle, you can't use the super-accurate, super-slow physics engine that nature actually uses (which is like trying to calculate the trajectory of every single raindrop in a storm). Instead, you use a "shortcut" engine—a set of rules called an Interatomic Potential.

Think of these rules like a recipe. If the recipe is perfect, your virtual atoms will behave exactly like real atoms: they will melt at the right temperature, stretch like rubber, or snap like glass. If the recipe is off by even a tiny bit, your virtual world falls apart, or the metal melts at the wrong temperature.

The Problem:
For a long time, scientists could only tune these recipes to get the "room temperature" behavior right (like how hard a metal is). But getting the recipe right for high temperatures (like the inside of a star or the Earth's core) was incredibly hard. It's like trying to tune a car engine to run perfectly at 200 mph, but you only have data on how it runs at 20 mph.

The Solution:
The authors of this paper developed a new "tuning knob" that lets them adjust the recipe specifically to match the Gibbs Free Energy.

• What is Gibbs Free Energy? Think of it as the "happiness score" of a system. Nature always wants to be as "happy" (low energy) as possible. If you want to know if a solid block of metal will melt into a liquid, you just compare the "happiness scores" of the solid and the liquid. Whichever is happier wins.
• The Goal: They wanted to tweak their recipe so that the "happiness score" of their virtual atoms matched the "happiness score" of real atoms (calculated by super-computers) at extreme heat and pressure.

---

How They Did It: The "Thermodynamic Integration" Hike

The authors used a clever mathematical trick called Hamiltonian Thermodynamic Integration (HTI). Let's break this down with an analogy.

Imagine you are standing at the bottom of a foggy mountain (your current, imperfect recipe). You want to reach the summit (the perfect recipe that matches reality). You can't see the top, but you have a compass that tells you exactly which direction is "uphill" (or in this case, which way to tweak the numbers to get closer to the target).

1. The Gradient (The Compass): Instead of guessing, they calculate the "slope" of the mountain. They ask: "If I change this specific number in my recipe by a tiny bit, how much does the 'happiness score' (Free Energy) change?"
2. The Iteration (The Hike): They take a step in the direction that improves the score. Then they stop, recalculate the slope, and take another step.
3. The Result: After just a few steps (iterations), they are right at the summit. Their virtual atoms now have the exact same "happiness score" as the real atoms at high temperatures.

The "RECAL" Shortcut (Don't Re-Run the Movie)

One of the smartest parts of their method is something they call the RECAL protocol.

Usually, to test a new recipe, you have to run a massive, time-consuming movie simulation of the atoms moving around. If you have to do this 100 times to tune the recipe, it takes forever.

The authors found a shortcut. They realized they didn't need to re-watch the whole movie. They could just take the "frames" (the snapshots of where atoms were) from a previous run and simply re-calculate the energy using the new numbers. It's like having a photo album of a party; you don't need to throw the party again to see how the lighting would look if you changed the bulbs. You just look at the photos and imagine the new lighting. This made the process incredibly fast.

The Proof: Testing on Three Levels

To prove their method works, they tested it on three different "levels" of difficulty:

1. The Toy Model (The Uhlenbeck-Ford Model):

   • Analogy: Testing a new steering wheel on a remote-controlled car in a parking lot.
   • Result: It worked perfectly. They could tune the "happiness score" to match a known mathematical target almost instantly.

2. The Nickel System (Pure Metal):

   • Analogy: Tuning the engine of a real car to run at 200 mph.
   • Context: They looked at Nickel under the crushing pressure of the Earth's core.
   • Result: They adjusted the recipe so that the virtual Nickel melted at the exact same temperature as real Nickel. Before this, the virtual Nickel was melting hundreds of degrees too early or too late.

3. The Iron-Oxygen Soup (Binary Liquid):

   • Analogy: Tuning a recipe for a complex soup where you can change the ratio of salt to pepper.
   • Context: They looked at a mix of Iron and Oxygen (like the stuff inside the Earth's core).
   • Result: They successfully tuned the recipe so that the "happiness" of the mixture matched reality, no matter how much oxygen was in the mix.

Why This Matters

This paper is a game-changer for materials science because:

• It's Fast: It doesn't require super-computers to run for years.
• It's Accurate: It allows scientists to create "virtual materials" that behave correctly even in extreme environments (like inside a star or a nuclear reactor).
• It's Flexible: It can be combined with other tests (like how stretchy the metal is) to make the recipe perfect in every way.

In a nutshell: The authors built a "GPS" for tuning atomic recipes. Instead of blindly guessing how to make virtual atoms behave like real ones at high heat, they gave scientists a direct path to the target, making it possible to simulate the most extreme conditions in the universe with high precision.

Technical Summary

Here is a detailed technical summary of the paper "Incorporating Gibbs free energy into interatomic potential fitting" by Wei and Sun.

1. Problem Statement

The development of accurate interatomic potentials (IPs) for molecular dynamics (MD) simulations is critical for predicting material properties, particularly under extreme conditions (high temperature and pressure).

• Limitations of Current Methods: While ab initio molecular dynamics (AIMD) offers high accuracy, it is computationally expensive and limited in scale. Classical MD relies on semi-empirical potentials (e.g., Embedded Atom Method, EAM), which are efficient but often struggle to reproduce high-temperature thermodynamic properties.
• The Gap: Existing fitting protocols typically target 0 K properties (lattice constants, elastic tensors) or structural data (pair distribution functions). Directly fitting high-temperature Gibbs free energy ($G$) is difficult because:
  1. $G$ is not a direct output of standard MD simulations.
  2. Historically, calculating high-temperature free energy from first principles was computationally prohibitive.
  3. Without accurate $G$, potentials fail to predict phase diagrams, melting points, and nucleation driving forces correctly.
• Objective: To develop a robust method to incorporate high-temperature Gibbs free energy data (derived from AIMD or deep learning) directly into the fitting process of semi-empirical interatomic potentials.

2. Methodology

The authors propose an iterative fitting framework based on Hamiltonian Thermodynamic Integration (HTI) and the Newton-Raphson optimization method.

A. Theoretical Framework

• Thermodynamic Integration: The difference in Gibbs free energy between an initial potential parameter set ($\mathbf{X}_0$) and a target set ($\mathbf{X}_{target}$) is expressed as an integral of the potential energy gradient along a path in parameter space:
  $$G(\mathbf{X}) - G(\mathbf{X}_0) = \int_{\mathbf{X}_0}^{\mathbf{X}} \left\langle \frac{\partial U}{\partial \mathbf{X}'} \right\rangle_{\mathbf{X}'} d\mathbf{X}'$$
• Iterative Update (Newton-Raphson): Instead of solving the integral directly, the authors use an iterative scheme to update parameters $\mathbf{X}_n$ to approach the target free energy $G_{target}$:
  $$\mathbf{X}_{n+1} = \mathbf{X}_n - \frac{G(\mathbf{X}_n) - G_{target}}{\nabla G(\mathbf{X}_n)}$$
  where $\nabla G(\mathbf{X}_n) = \langle \frac{\partial U}{\partial \mathbf{X}} \rangle_{\mathbf{X}_n}$.
• The RECAL Protocol: A key efficiency feature is the "RECAL" (RECALculate) protocol. Instead of running new MD simulations for every parameter update, the ensemble of atomic configurations generated by the current potential is reused. The energy gradients ($\partial U / \partial \mathbf{X}$) are computed by re-evaluating the energy of these existing configurations under the new parameter set.
• Multi-Property Fitting: The method extends to fit multiple targets simultaneously (e.g., free energy, elastic constants, liquid structure) by constructing a Jacobian matrix of derivatives for all properties and solving a linear system to update parameters.

B. Implementation Details

• Potential Form: The method was applied to Embedded Atom Method (EAM) potentials, specifically using polynomial basis functions for the pair potential and embedding energy. This allows for analytical derivatives of the energy with respect to parameters.
• Simulation Setup:
  • Code: LAMMPS.
  • Ensembles: NPT (isothermal-isobaric) for free energy calculations; NVT for specific property checks.
  • Algorithms: Nosé-Hoover thermostat/barostat; Solid-Liquid Coexistence (SLC) for melting points.
• Target Data: High-accuracy Gibbs free energy data from previous DFT studies (specifically for Ni and Fe-O systems under Earth's inner core conditions).

3. Key Contributions

1. Direct Free Energy Fitting: The paper introduces a practical, iterative algorithm to directly minimize the error between a semi-empirical potential's Gibbs free energy and a target ab initio value.
2. Computational Efficiency: By utilizing the RECAL protocol and analytical derivatives (for polynomial-based EAM), the method avoids the need for expensive re-sampling of phase space during every iteration.
3. Generalizability: The framework is shown to be compatible with standard fitting targets (elasticity, structure) and can be adapted to other potential forms (including machine learning potentials) via perturbation methods if analytical derivatives are unavailable.
4. Robustness: The method is demonstrated to be robust against variations in initial parameter guesses.

4. Results

The authors validated the method on three distinct systems:

• Case 1: Uhlenbeck-Ford Model (Toy Model)

  • Goal: Fit a single parameter ($\sigma$) to match an analytical free energy.
  • Outcome: The iterative scheme converged rapidly (within 3 iterations) to the target free energy, validating the mathematical derivation and the RECAL protocol.

• Case 2: Pure Nickel (Ni) at High Pressure (323 GPa)

  • Goal: Fit an EAM potential for Ni to reproduce ab initio free energy differences between liquid and solid phases (bcc, hcp, fcc) at 6000–7000 K.
  • Process: Started with a potential fitted only to elastic constants and liquid structure (which had a ~200 K error in melting point).
  • Outcome: After incorporating free energy targets (at $T_m$ and $T_m - 400$ K), the potential converged in 2 iterations.
    • The RMS error for free energy dropped to 2.5 meV/atom.
    • The predicted melting temperature for fcc-Ni was within 40 K of the DFT target.
    • Elastic constants and liquid structure ($g(r)$) remained well-fitted.

• Case 3: Binary Iron-Oxygen Liquid (Fe$_{1-x}$O$_x$) at 323 GPa, 5500 K

  • Goal: Fit the composition-dependent mixing free energy and liquid structure.
  • Process: Used an existing EAM potential as a starting point, which showed significant deviation in mixing free energy.
  • Outcome: The method successfully fitted the mixing free energy for various oxygen concentrations (4%–20%).
    • Converged in 3 iterations with a final mixing free energy RMSE of 0.4 meV/atom.
    • The liquid pair distribution functions ($g(r)$) matched ab initio data accurately.
    • Robustness Test: Starting from four different initial parameter sets with large initial errors (2–27 meV/atom), the method converged to the same solution within one iteration, demonstrating insensitivity to initial conditions.

5. Significance

• Bridging the Accuracy-Efficiency Gap: This work enables the development of semi-empirical potentials that retain the computational speed of classical MD while achieving the thermodynamic accuracy of ab initio methods for high-temperature applications.
• Planetary Science Applications: The specific application to Ni and Fe-O under Earth's inner core conditions (high P, high T) highlights the method's utility for geophysics, where experimental data is scarce and accurate phase diagrams are crucial.
• Future-Proofing: The authors note that while demonstrated on EAM, the logic applies to Machine Learning Potentials (MLPs). Since MLPs are currently trained primarily on forces and energies (0 K), this method provides a pathway to train MLPs on thermodynamic data, significantly improving their reliability for phase transition and melting simulations.
• Methodological Shift: It moves the field away from indirect fitting (e.g., adjusting melting points via Gibbs-Duhem integration) toward direct, rigorous fitting of the fundamental thermodynamic potential ($G$).