The Everything Bagel Phase

I was thinking about how universal these foundation MLIPs are, or more specifically, what universality actually means. Some use "foundation" or "pre-trained," which is probably more accurate, i.e. a training recipe of broad data and then deploy, whereas universal would be a potential energy surface (PES) that still behaves for any chemistry regardless if it was captured in the training data. The word fight is mostly vocabulary. Still, the blunt question is, what does "universal" get you once the chemistry no longer resembles anything in the training set, where out-of-distribution means that local neighborhood was not an environment the fit saw.

So here is the stunt I came up with: take one cell and every element from hydrogen to francium, $1 \leq Z \leq 87$, put them in a periodic box at equal atomic fractions, attach a foundation MLIP with coverage, and ask how sampling different configurations via swapping sites and relaxing the cell samples the PES.

This is not some insane equiatomic high-entropy alloy ... although it is an insane chemistry ... but it mixes metals, nonmetals, halogens, noble gases, lanthanides, and radioactive species that barely exist long enough to sit in a box. The chance anyone synthesizes anything close is next to nil. For fun let me call this fictitious initial state the Everything Bagel Phase. Phase is loose here, it names the starting configuration, not a claim that 87 elements form a stable homogeneous thermodynamic phase.

The digital experiment

Equal mole fractions, $x_i=1/87$, so $N=87m$ with $m$ copies of each element. The probe³ below uses $m=1$ ($N=87$), which may introduce finite-size noise effects. A serious run wants tens or hundreds of atoms per species but not trying to spend more than few minutes on this with my compute.

I kept the composition closed ... and of course we can't capture nuclear decay since we have classical nuclei on a learned PES, no explicit electrons or plasma chemistry. Also we are working with periodic boundaries. The search below is a labeling on this PES, so it is not a claim that the cell found the thermodynamic ground state of an 87-component mixture.

The oracles here are just whichever few MLIP ASE calculators will take the whole cell. The calculator has to take every $Z$ from 1 through 87, and the thing to note on the Matbench Discovery diatomics page is that Po, At, Rn, Fr, and Ra are already off the scored curves so the metric stays compatible with Materials Project compositions, the same cut MACE-MPA-0 [1] and M3GNet [2] follow when their element tables skip $Z=84$ through $88$, which puts those models out on this cell. I stayed on the three that accept the full cell, mattersim_v1_5m, orb_v3 (orb-v3-conservative-inf-mpa), and orb_v3_direct_inf (orb-v3-direct-inf-mpa) [3], with orb-models 0.5.4. Both ORB models one-hot atomic number out to 117, and their reference-energy tables are filled through Ra, with Po, At, Rn, and Fr in the table and a filler near $-4.9$ eV from Ac onward. Open MatterSim v1.0.0-5M [4] is one of the three, its embedding runs to $Z=94$ and its paper describes support for the first 89 elements, and the H through Fr cell has a finite energy on all three. The swap and cell-relax loops are in the footnotes.³

Random packing at some guessed "average density" mostly tests overlapping cores. There is no meaningful average density for H, Xe, W, Hg, F, and Fr together I can think of. A less pathological path is a radius-aware pack at low density, then a soft-repulsive cleanup, then the MLIP. Figure 1 is one instance of H through Fr ($N=87$) in a cubic cell of side $13.18$ Å ($V=2290.6$ Å$^3$), density from a covalent-radius volume proxy at packing fraction $0.64$ ($\rho \approx 6.64$ g/cm$^3$), and the soft-overlap push got the shortest pair to about $1.44$ Å.

Figure 1. One atom of each element H through Fr in a soft-overlap-relaxed periodic cell. Left, colored by atomic number $Z$. Right, the same cell colored by chemical class.

Local minima searching

Starting with a trial, which is one proposal of two random atom pairs swapping labels, then FIRE relaxes the cell and the positions together on a FrechetCellFilter at $P=0$, and the state updates only when $E/N$ drops. The shortest pair permitted is at least $1.0$ Å, and the volume is constrained to stay between $0.45$ and $2.2$ times the initial volume.

There are sixty trials on one seed for each of mattersim_v1_5m, orb_v3 (conservative, unlimited neighbors), and orb_v3_direct_inf (direct forces, same cutoff), so one hundred twenty pair swaps were proposed per model, of which MatterSim kept fourteen, ORB v3 conservative kept thirteen, and ORB v3 direct kept fifteen. Rather than compare absolute $E/N$ across models, in Figure 2 the drop from each model's own start against accepted trials is shown. MatterSim finishes $0.20$ eV/atom down, its largest single drop being Mg with S and Re with Ni ($-0.042$), ORB v3 conservative finishes $0.19$ eV/atom down, its largest single drop being Pb with Y and Dy with W ($-0.043$), and ORB v3 direct finishes $0.16$ eV/atom down, its largest single drop being Mo with As and Sn with Ni ($-0.031$), with a short further relax of those three landing at $-4.949$, $-5.191$, and $-5.137$ eV/atom.

The first trial is the same on all three, Ti with At and Eu with Ni, kept by MatterSim and by the direct model, and the conservative model leaves the state unchanged. Density barely changes, and those values and shortest/largest first-neighbor distances are shown in the table below ... just as point of reference, stainless steel is about $8$ g/cm$^3$.

	Packed	MatterSim	ORB cons.	ORB direct
Density (g/cm$^3$)	6.64	6.66	6.68	6.60
Shortest pair (Å)	1.77	1.75	1.79	1.72
Largest first-neighbor pair (Å)	3.33	3.40	3.43	3.47

Because every element appears once, a partial radial distribution function (RDF) by element is a single distance, so the structure info is shown in Figure 3, each entry is how many atoms of one class have the other class as their closest neighbor.² MatterSim piles metals next to metals ($11$ to $16$), ORB conservative sits in between, and ORB direct does not follow that metal-metal rise, while lanthanide-halogen contacts fall from $5$ on the packed cell toward $1$ under ORB direct. Alkali-halogen stays rare, and metal-halogen is zero for every model. With one atom per element these counts are small, so this is a packing-and-disagreement check, and not so much a motif one.

The final structure in Figure 4 is the MatterSim cell after those sixty trials, $-4.949$ eV/atom at $\rho \approx 6.66$ g/cm$^3$. Sixty trials are a local sample of the labeling, not a ground state.¹

Figure 2. Sixty trials on the packed $N=87$ cell, plotted as accepted trials only. Each trial swaps two random pairs and FIRE-relaxes cell and positions. The curves are the energy drop from each model's own start for MatterSim, ORB v3 conservative, and ORB v3 direct. Leader labels mark the first kept trial and each model's largest single drop.

Figure 3. Nearest-neighbor class pairing after sixty trials. Each panel is a count matrix, row is the center atom's class, column is its closest neighbor's class, for the packed start and each model's incumbent. Class abbreviations marked on the axes are defined in the footnotes.

Figure 4. MatterSim cell after sixty trials, $E/N\approx -4.949$ eV/atom, $\rho\approx 6.66$ g/cm$^3$. Left, colored by atomic number $Z$. Right, the same cell colored by chemical class.

What the swaps say

Of course trying to synthesize the Everything Bagel is a stunt rather than providing utility. My main interest was just to see what foundation MLIPs would do if handed some ridiculous chemical system cell. I also was curious about the agreement between them, although in this case its just really two types and one of a different flavor. The MLIPs did return finite energy, and under the volume constraint the drop from each model's own start seemed reasonable and controlled. While the kept atom pairs are not the same, and after the first trial those three searches are already on different structures, so there is a split in the search/acceptance.

I'm not entirely sure this answered my initial question about how well these would do with something so out of distribution, although truth is my question may have actually not really been well posed because there is nothing to validate against. But it is an interesting exercise about MLIPs and their behavior under extreme conditions.

Footnotes

---

1. You can view and download the same final cell as Figure 4, rotatable in MatterViz. ↩

2. In Figure 3, Ae* is alkaline earth, Ln† is lanthanide, and MG‡ is other main-group (B, Al, Si, P, S, Ga, Ge, As, Se, In, Sn, Sb, Te, Tl, Pb, Bi, Po). ↩

3. Core ASE loops (calculator attach omitted) code block below. Each swap trial relaxes the cell and the positions together. ↩↩

   ASE Loops

   
             from ase.optimize import FIRE
             from ase.filters import FrechetCellFilter
   
             def relax(atoms, fmax=1.0, steps=40):
                 """FIRE on cell and positions at zero pressure."""
                 filt = FrechetCellFilter(atoms, scalar_pressure=0.0)
                 FIRE(filt, 
                      logfile=None,
                      maxstep=0.04).run(fmax=fmax, steps=steps)
                 return atoms
   
             def swap(numbers, rng, n=2):
                 """Swap n random pairs of species labels."""
                 out = numbers.copy()
                 for _ in range(n):
                     i, j = rng.choice(len(out), 2, replace=False)
                     out[i], out[j] = out[j], out[i]
                 return out
           

References

[1] I. Batatia, P. Benner, Y. Chiang, others, A Foundation Model for Atomistic Materials Chemistry, The Journal of Chemical Physics 163 (2025) 184110. https://doi.org/10.1063/5.0297006.

[2] C. Chen, S.P. Ong, A Universal Graph Deep Learning Interatomic Potential for the Periodic Table, Nature Computational Science 2 (2022) 718–728. https://doi.org/10.1038/s43588-022-00349-3.

[3] B. Rhodes, S. Vandenhaute, V. {\v{S}}imkus, others, {Orb-v3}: Atomistic Simulation at Scale, arXiv Preprint arXiv:2504.06231 (2025). https://doi.org/10.48550/arXiv.2504.06231.

[4] H. Yang, C. Hu, others, {MatterSim}: A Deep Learning Atomistic Model Across Elements, Temperatures and Pressures, arXiv Preprint arXiv:2405.04967 (2024). https://doi.org/10.48550/arXiv.2405.04967.

Reuse and Attribution

Making KMC Cool

If you've ever tried to simulate a deposition process for film growth, you'll quickly find that a lot of different events take place. Moreover, these events occur at different timescales. In a deposition process, for example, the particles/species are arriving at the surface at a much faster timescale than, say, surface diffusion. To put it in context, say we have a 1 mm$^2$ area, and our particle flux is $10^{21} \frac{\#}{\text{m}^2 \cdot \text{s}}$, while the diffusion coefficient is $10^{-4} \frac{\text{cm}^2}{\text{s}}$. If we assume some desired number of particles to deposit to create a film on this area, then the timescales to see events are something like microseconds for deposition and thousands of seconds for diffusion.

These timescales are vastly different and hence will need to be handled differently. In the deposition process, we could treat it directly using atomistic simulations, but it turns out that for most practical fluxes, the timescales are just out of reach. To solve this problem, we can turn to a technique called kinetic Monte Carlo (KMC). For diffusion, it's clear that we need to use KMC or solve mass transport equations.

Kinetic Monte Carlo (KMC)

Most people are familiar with Monte Carlo methods. These are stochastic approaches to solving problems where, instead of computing the solutions directly, we sample from probable events. The probability is determined by the probability mass or density function. Monte Carlo methods can be used to solve all kinds of problems, for example performing integration. Say we have a function we want to integrate but don't know a deterministic algorithm or analytical solution. We can use Monte Carlo to randomly sample the area under the curve and estimate a definite integral. The ultimate accuracy depends on how much sampling you decide to perform.

Monte Carlo methods let us perform random sampling to estimate the behavior of a system or process; if we think back to our integral example, we can think of the system or process being described by a high-dimensional probability density or mass function that we are trying to integrate to get the expected behavior (i.e., value).

So what does it mean when we add the prefix "kinetic" to Monte Carlo? It simply designates that the process has a time component to the evolution. Think of a chemical reaction or diffusive dynamics. In KMC, the events that are described are derived from rates, $\nu$, which are usually in inverse time units. The probability of an event is determined by the rate, and this in turn informs the elapsed time. To expand, let's look at the math.

Say we have some event, $\Omega_i$, that corresponds to a diffusive hop. There is some energy required to overcome for the diffusive hop to occur, given as $E_{a}^i$. We can describe the rate as $\nu_i = \nu_0^i \exp\left(-\frac{E_{a}^i}{k_B T}\right)$, where $k_B$ is the Boltzmann constant, $T$ is the temperature, and $\nu_0^i$ is the prefactor (attempt frequency) for event $i$. So now assume we have $N$ events that describe some diffusive hop to a different site/location. We can describe the probability of any one event occurring by:

\begin{equation} P(\Omega_i; E_a^i, \nu_0^i) = \frac{\nu_0^i \exp\left(-\dfrac{E_{a}^i}{k_B T}\right)}{\sum\limits_{j=1}^N \nu_0^j \exp\left(-\dfrac{E_{a}^j}{k_B T}\right)} \label{eq:probability} \end{equation}

I intentionally include that the probability is parameterized based on selecting the activation energy and the prefactor. However, we could easily consider this as a conditional probability if we have some knowledge of the range of values for these terms.

Thus, you have a way to sample from the event catalog now. So how do you define time based on the sampling? If an event $\Omega_i$ is selected, then we can calculate the time increment associated with this event. In kinetic Monte Carlo, we generate a second random number to calculate the time increment. Specifically, the time increment $\Delta t$ is given by:

\begin{equation} \Delta t = -\dfrac{1}{R_N} \ln r \end{equation}

where $r$ is a uniform random number between 0 and 1, and $R_N$ is the cumulative rate for all events/processes, i.e., $R_N = \sum\limits_{j=1}^N \nu_j$. This time increment reflects the stochastic time elapsed before the occurrence of event $\Omega_i$.

It's a pretty simple formalism but very powerful to simulate kinetic processes at the right timescales. The challenge is always two-fold:

1. Treating a global event table with all process timescales and rates can lead to the need for very large sampling trials. In other words, if we treat the deposition and diffusive processes on the same event table, diffusive events will have a lower chance of being sampled. The easiest way to address this is simply by decoupling them and having two KMC times simultaneously evolve.

2. A priori knowledge of what should happen, i.e., this atom deposits at this site and reaction $\alpha$ occurs if it has species A but does reaction $\beta$ if it has species B. For well-understood systems that have simple species, this is manageable, but it should be clear that the combinatorics of events for multi-species deposition and diffusion becomes prohibitive, and more importantly, we have many model parameters to decide on.

Lattice KMC

One way to simplify the KMC approach and to simulate crystalline materials is to use a lattice formalism. This means we can describe lattice points where events can occur and thus it somewhat restricts what can happen. In particular, diffusion hops are constrained to occur between neighboring lattice sites. This works well for crystalline materials with a single basis, but one may need to fall back on mesh methods to treat more complicated crystals.

🧠 Thinking out Loud

I haven't really thought too much about meshing approaches, but the idea would be to define nested lattices where each basis point is defined. One could then define the events that are allowed within a lattice and between lattices. Take a Cesium Chloride crystal; we would have two interpenetrating simple cubic lattices.

Some of the challenges with lattice KMC are that you can't treat amorphous materials and grain boundaries. And it gets complicated for multiple species in my view.

You're also still plagued with the fact that you need to enumerate all the events and assign the rates. Maybe not too difficult for a monoatomic species, where you can find a good amount of information for the surface deposition/diffusion and bulk diffusion.

On-the-Fly Off-Lattice KMC

First, imagine we remove the lattice constraint and let points (representing atoms or molecules) occupy any point in space. Then, say we have a way to determine for a species the value of $E_a^i$ for some motion it could have in space. Think of a surface atom hopping on top of another surface atom. If you could find and enumerate all the possible motions such an atom could have, you would come up with a list of activation energies and therefore rates. Thus, you would have cataloged an event table you could then sample from to perform KMC. This is on-the-fly KMC.

Basically, you need to have a model/description for the energy barrier associated with different dynamical processes. We do this through the potential energy surface (PES) model, which is the functional surface describing, for a given configuration of atoms in space, what is the potential energy and the associated gradients (i.e., forces).

Once you have a way to describe the PES, then we need a methodology to sample different pathways (e.g., reaction coordinates) to determine what pathways through space an atom(s) could take and the associated energetic barrier. Then, through transition state theory, we can determine the rate of each pathway; that is, we get an Arrhenius expression for the rate of a reaction coordinate:

\begin{equation} \nu_i = \nu_0 \exp\left(-\dfrac{E_{a}^i}{k_B T}\right) \end{equation}

Challenges with Transition State Theory

The two major hurdles with transition state theory are:

1. You need to have an accurate PES description for the chemical system you're trying to describe. This usually requires force-field or interatomic potential development (see my posts discussing such topics).

2. The pathway along each reaction coordinate either needs an intelligent way to search around an initial point or requires an initial and final point to be connected. The former are called single-ended methods, and techniques like the dimer method or growing string method are used. These essentially find the saddle points of the PES, and thus the final state is just the "downhill" path from the saddle point. The other method is dual-ended, where the initial and final coordinates are known, but possible paths are not. We then can find the possible barriers and usually only tabulate the one that is the lowest.

Most of the effort is always centered around proper PES description and setting up the TST calculations, which can be truly challenging. This is why we probably don't see many on-the-fly KMC codes in practice. The only one that I came across was EON from the Henkelman group at UT Austin [1]. It uses a client-server approach where the server runs the TST calculations and then hands the event table back to the KMC client.

Putting It All Together

In practice, I think the outline for an on-the-fly off-lattice KMC code that simulates deposition and diffusion would be something like:

flowchart TD Start{Start Simulation} --> Init[Initialize Deposition Events] Init --> OuterLoop{Max steps?} OuterLoop -->|No| Deposition([Sample deposition event]) Deposition --> InnerKMC[Perform inner KMC loop] InnerKMC --> Update[Update KMC times] Update --> Stats[(Store Statistics)] Stats --> OuterLoop OuterLoop -->|Yes| End(((Terminate))) subgraph InnerKMC [Inner KMC Loop] TST(Transition State Calculation) TST --> SampleTS(Sample from transition state table) end

As this blog was more a discussion on KMC, just so I can jot down my thoughts, I'll leave it at this for now. If you want to read more on lattice KMC, you can read the chapter in the book by LeSar on it; it's a pretty good starting point [2].

References

[1] S.T. Chill, M. Welborn, R. Terrell, L. Zhang, J.-C. Berthet, A. Pedersen, H. Jónsson, G. Henkelman, EON: software for long time simulations of atomic scale systems, Modelling Simul. Mater. Sci. Eng. 22 (2014) 055002. DOI.

[2] R. LeSar, Introduction to Computational Materials Science: Fundamentals to Applications, Cambridge University Press, Cambridge; New York, 2013. DOI

---

Reuse and Attribution