Interpolation and differentiation of alchemical degrees of freedom in machine learning interatomic potentials

TL;DR

This work introduces continuous alchemical degrees of freedom within graph-based interatomic potentials by augmenting the input graph with alchemical atoms weighted by a compositional vector $\bm{\lambda}$. The approach yields end-to-end differentiability with respect to composition, enabling gradient-based optimization and efficient nonequilibrium free-energy calculations via thermodynamic integration along alchemical pathways. The authors demonstrate the method on solid-solution representations, solid-solution disorder energetics, and alchemical free-energy calculations for vacancies and phase transformations, achieving accurate interpolation beyond Vegard’s law and faster convergence than traditional methods like Frenkel–Ladd. The framework offers a scalable avenue to model compositional disorder, optimize material compositions for targeted properties, and quantify composition-dependent thermodynamics using pre-trained universal MLIPs.

Abstract

Machine learning interatomic potentials (MLIPs) have become a workhorse of modern atomistic simulations, and recently published universal MLIPs, pre-trained on large datasets, have demonstrated remarkable accuracy and generalizability. However, the computational cost of MLIPs limits their applicability to chemically disordered systems requiring large simulation cells or to sample-intensive statistical methods. Here, we report the use of continuous and differentiable alchemical degrees of freedom in atomistic materials simulations, exploiting the fact that graph neural network MLIPs represent discrete elements as real-valued tensors. The proposed method introduces alchemical atoms with corresponding weights into the input graph, alongside modifications to the message-passing and readout mechanisms of MLIPs, and allows smooth interpolation between the compositional states of materials. The end-to-end differentiability of MLIPs enables efficient calculation of the gradient of energy with respect to the compositional weights. With this modification, we propose methodologies for optimizing the composition of solid solutions towards target macroscopic properties, characterizing order and disorder in multicomponent oxides, and conducting alchemical free energy simulations to quantify the free energy of vacancy formation and composition changes. The approach offers an avenue for extending the capabilities of universal MLIPs in the modeling of compositional disorder and characterizing the phase stability of complex materials systems.

Figures (6)

• Figure 1: Alchemical modification scheme for machine learning interatomic potentials. (a) Alchemical graph augmentation: The relevant original atoms are split into alchemical atoms with different elemental identities, which are associated with alchemical weights $\lambda_i$. (b) Alchemical message passing: At the message aggregation step (equation \ref{['eq:alchemical_message']}), each message contribution from neighboring atoms is weighted according to the asymmetric weighting scheme in equation \ref{['eq:edge_weight']}. Only the weights from alchemical to non-alchemical atoms are weighted according to the alchemical weights of the source atoms to ensure consistency with the message-passing scheme in the original graph. (c) Alchemical energy readout: The energy contributions from the alchemical atoms are weighted according to their respective alchemical weights.
• Figure 2: Lattice parameters for solid solutions. (a) The starting structures, CeO2 and BiSBr, for solid solutions. (b) Lattice parameter $a$ for Ce_1-xM_xO_2 (M = Zr, Sn) as a function of the compositional weight $x$. (c) Lattice parameters $a$, $b$, and $c$ for BiSX_1-xY_x (X, Y = Cl, Br, I) as a function of $x$. The upper panels are the result of the alchemically modified MACE-MP-0 medium model batatia2024foundation, and the lower panels are the experimental results from baidya2016understanding and schultz2014strong for (b) and (c), respectively. Arrows in the rightmost panels indicate the composition with the minimum value of $c$.
• Figure 3: Compositional optimization. (a) Lattice parameter optimization for solid solutions of LiCl, NaCl, and KCl. The left panel shows optimized lattice parameters as a color gradient, obtained by relaxing the cell geometry for each compositional weight. The right panel displays hydrostatic stress, with color intensity representing stress magnitude and arrows indicating gradient direction, calculated by fixing the cell dimensions to those of NaCl. Since the energy output is end-to-end differentiable with respect to the alchemical weights, the composition can be optimized to match target cell dimensions (minimizing stress) by following these gradients. The compositions with cell parameters matching NaCl (left) and those obtained by minimizing stress (right), indicated by the dotted lines, align in both figures. (b) The optimization for the lattice-matching condition for solid solutions Al_1-xSc_xN and Al_1-xY_xN with GaN. The most stable polymorph structures are shown on the left. The plot on the right shows the cell dimension $a$ obtained by optimizing for each compositional weight (Scan), calculated from the corresponding supercell (Supercell), and the compositional weights optimized by gradient descent to match the $a$ value for GaN (Optimized). All results are obtained using the alchemically modified MACE-MP-0 medium model batatia2024foundation.
• Figure 4: Disordered energetics in multicomponent perovskite oxides. (a) Crystal structure schematics for fully cation-disordered A2B$'$B$"$O6 perovskite oxide solid solutions, illustrating different alchemical supercell sizes and number of atoms, and representative 320-atom 4$\times$4$\times$4 SQS structures. (b) MACE-relaxed energy differences between cation-disordered $2 \times 2 \times 2$ alchemical cells and smaller or larger supercells, evaluated across 100 A2B$'$B$"$O6 compositions from peng2024learning. (c) Difference between the unrelaxed and MACE-relaxed structures for various alchemical cell sizes, including $4 \times 4 \times 4$ SQS structures, quantified by cosine distance between local structure fingerprints law2023upperpeng2024data. (d) Comparison of MACE-relaxed energies for $2 \times 2 \times 2$ alchemical cells versus $4 \times 4 \times 4$ SQS structures. (e) MACE-relaxed energy differences among $2 \times 2 \times 2$ alchemical cells, $4 \times 4 \times 4$ SQS structures, and all cation-ordered configurations with four B$'$ and four B$"$ cations on eight B sites in the $2 \times 2 \times 2$ supercell. Compositions are sorted by energy differences between $4 \times 4 \times 4$ SQS and lowest-energy ordered arrangements. Experimentally characterized ordered and disordered compositions peng2024data are marked in the upper and lower regions, respectively. (f) Receiver operating characteristic (ROC) curves for experimental order/disorder classification peng2024data based on relative energy values of $4 \times 4 \times 4$ SQS or $2 \times 2 \times 2$ alchemical cells in (e) with respect to the lowest-energy cation-ordered arrangements, with area under the curve (AUC) values shown. All results are derived using the alchemically modified MACE-MP-0 medium model batatia2024foundation.
• Figure 5: Free energy of vacancy formation in BCC iron. (a) Transformations used to determine the Gibbs free energy of the perfect crystal and the crystal with a defect. The alchemical pathway used here transforms the perfect crystal into the crystal with a defect and a single atom attached to a spring to avoid diffusion. (b) The intermediate state parameterized by $\lambda$ for the alchemical pathway in (a). The atom to be removed is assigned an alchemical weight of $1 - \lambda$, and the energy of the harmonic oscillator is scaled by $\lambda$. (c) The free energy of vacancy (equation \ref{['eq:vacancy']}) computed by the Frenkel--Ladd path and alchemical path. (d) Statistical efficiency for the Frenkel--Ladd paths and alchemical path at 100 K against the switching time. Upper panel shows the deviation of Gibbs free energy from the reference value at the longest switching time (60 ps), and the lower panel shows average dissipated energies (equation \ref{['eq:energy_diss']}).