Gradient-based grand canonical optimization enabled by graph neural networks with fractional atomic existence

TL;DR

This work addresses the challenge of exploring configurational space with variable atom counts by introducing fractional atomic existence as a continuous variable in graph neural network interatomic potentials. It derives a differentiable Gibbs free energy $ΔG(x,q)$ and develops a gradient-based grand canonical optimization workflow, enabled by a SchNet-based message-passing architecture augmented with existence variables. Key contributions include formalizing existence conditions (embedding/readout equivalence and continuity), implementing fractional existence in MP, and demonstrating targeted configuration generation on Cu(110) copper-oxide surfaces with validation against DFT. The approach is architecture-agnostic and compatible with pretrained models, offering a pathway to generative, chemistry-aware exploration under chemical potentials and potentially extending to diffusion-type generative models.

Abstract

Machine learning interatomic potentials have become an indispensable tool for materials science, enabling the study of larger systems and longer timescales. State-of-the-art models are generally graph neural networks that employ message passing to iteratively update atomic embeddings that are ultimately used for predicting properties. In this work we extend the message passing formalism with the inclusion of a continuous variable that accounts for fractional atomic existence. This allows us to calculate the gradient of the Gibbs free energy with respect to both the Cartesian coordinates of atoms and their existence. Using this we propose a gradient-based grand canonical optimization method and document its capabilities for a Cu(110) surface oxide.

Figures (8)

• Figure 1: Illustration of the different steps involved in message-passing. (a) Graph construction and neighbourhood definition with two node embeddings shown. (b) Message passing from neighbouring nodes with two messages shown. (c) Aggregation of messages to compute the full message, with adaptations to account for node existence shown in red. (d) Node embedding update based on the message and the previous embedding. (e) Readout layer using node predictions to compute a total property modified to take into account the existence of nodes.
• Figure 2: Configuration snapshots from grand canonical optimization trajectory. The area of the atoms is linearly proportional to their existence and lighter colors indicate atoms that are further from the surface slab for which the copper atoms are depicted in gray. Initially, the configurations involve many atoms of various existences inserted at random and, as the optimization progresses only those atoms that provide stability at the chosen chemical potentials survive with the rest going out of existence.
• Figure 3: The conditions on the inclusion of existence in a graph neural network, tested with the trained network described in Section \ref{['sec:training_dataset']}. In (a) two dimensions of the learned atomic embeddings are plotted for a CuO3 cluster which can be compared to (b) where the same embeddings are shown for the same configuration with an additional oxygen atom with null existence, this confirms that embedding equivalence is upheld. In (c) the existence of the additional oxygen atom is increased to 1/2 and in (d) it is increased 1, this shows how the embeddings smoothly sweep across the embedding space as a function of the existence. Finally in (e) both read-out equivalence and continuity are tested, showing the energy as a function of the existence of an atom in a small cluster along side the energy of the cluster without and with the extra atom.
• Figure 4: Left: The $c(6\times 2)$ copper slab along with the most favourable configuration found at the five different chemical potentials targeted with the grand canonical optimization scheme. Right: Gibbs free energy diagram calculated using the model and the configurations generated by the grand canonical optimization procedure. The small discrepancy of the Gibbs free energy of the Cu6O0 configuration, not being equal to zero, is caused by a slight relaxation of by the added layer compared to the bare slab.
• Figure 5: (a-c): Distributions of Gibbs free energies calculated at different chemical potentials (in eV), each color corresponds to configurations identified by our optimization procedure with a specific chemical potential for oxygen. In each case the configurations that were generated at the desired chemical potential have lower Gibbs energies, showings the methods capability at targeting specific conditions. (d): Distributions of the number of oxygen in configurations generated at different chemical potentials. (e): The number of oxygen and copper in configurations generated at $\Delta \mu_\mathrm{O} = -1.25$ eV.