Atomistic Calculations using GNN

Disclaimer

I'm still very much working through M3GNet, CHGNet, and MACE papers. As a result I may get some of the GNN concepts and descriptins incorrect. Apologies in advance. I will provide updates as I continue learning.

I'm in the process of better understanding the recent release of graph neural network (GNN) interatomic potentials. I'm primarily focused on M3GNet and CHGNet potentials which include 3-body effects by updating the GNN structure [1,2]. The CHGNet adds the feature that it captures dynamic charge transfer effects, although the predictive improvement seems marginal over M3GNet/MEGNet based on a recent matbench-discovery preprint [3]. I will say that in general, I'm fairly optimistic about this class of ML potential and its use as an initial effort to do MD simulations across the periodic table. For reference here are some figure-of-merits for the GNN from ref. [3]:

Table showing a comparison of GNN using Matbench Discovery [3].

Based on the results it's clear that M3GNet and CHGNet perform better overall. CHGNet is at the top and especially does better in the discovery acceleration factor (DAF) which is a measure of how well the model predicts a stable structure (as compared to a random/dummy model).

So what is allowing these potentials to work in the first place? The answer is intricate GNNs. So I'm going to do my best here to try and describe a GNN in my understanding. I'll probably get a lot of points wrong but here's my go at it.

We first need to define what a graph is. A graph is a connected network that contains 3 major components. A set of vertices or nodes, $v \in V$, a set of edges, $e \in E | e; \text{connects to}\; v$, and global state, $U$. Okay, this seems reasonable and at first glance, you can see how an atomic representation of a material can be mapped to a graph. The vertices/nodes are the atoms, the edges are the bonds, and the global state are system properties like density, pressure, and temperature. There are however two additional needs for an atomic system represented by a graph:

1. Atoms have defined features: mass, atomic number, charge, and valency.
2. Bonds can be characterized by binding energy and lengths but also bond-pairs matter because bonding angles and dihedrals are characteristic of the chemistry.

So what do you do, you assign features to the vertices and edges corresponding to these. If you break this down then you really have a large set of graphs where each corresponds to a graph with a scalar value assigned to the nodes and edges; these are like weights and bias maybe? Great, I think, but now you need to feed this into a neural network, and for potential, you need it to predict energies, forces, and stresses. Additionally, because physics is invariant to the orientation of the system, meaning if we rotate the entire system by $\pi/2$ we shouldn't see anything change, so our NN needs to handle such inputs (i.e. if we rotate a molecule the prediction is the same). I'll touch on this invariance condition at the end.

Alright, so we have some additional details about the graph or an atomic system and constraints for the NN. We then may want to ask about the feature dimensionality of the graph; or collection of graphs if we want to think like that. This can easily become very large and so one may want to find a latent space representation that is easier to train on to match the energies and forces.  To do this we can use a convolutional NN that finds a single graph representation. Then we can feed this graph into a multi-layer perceptron network to make predictions for the energies, forces, and stresses. 

So what does M3GNet and CHGNet do here? Based on how I understand the papers, they add something that is akin to an attention mechanism to CNN. What is attention? I'm still learning this as well but essentially it is a way to have the CNN consider what the environment looks like and use that information to update the subsequent graph. In other words, say an edge corresponding to a bond is passed through the NN but doesn't know anything about the other edges or nodes, then it will simply learn about what it means to be that type of edge in the grand scheme of the network. However, if the edge is informed, via an attention mechanism (?), about the environment then one updates what it means to be an edge/bond in that context. 

There is also the addition of some more physics-informed aspects of informing/updating the edge features (i.e. bonds) in the graph. This is done by the inclusion of many-body interactions which give sets of nodes and edges, then edge features can be updated. This is very similar to the Tersoff-style interatomic potential which uses a bond-order parameter to modify the pairwise interaction between atoms. In M3GNet this step seems to occur before the update (attention) in the CNN. Basically, a GNN is trained to take in nodes and edges where the edges have features representing the bond distance as an expansion in radial basis functions and these are what gets updated and pass on to the graph CNN. Again this is what I think is going on, but could be wrong.

If I get anything wrong here please leave a comment to correct me. Also on the topic of invariance, I believe that graph CNN can be structured so that they are resistant to translation and rotation of inputs [4].

References

[1] C. Chen, S.P. Ong, A universal graph deep learning interatomic potential for the periodic table, Nat Comput Sci. 2 (2022) 718–728. https://doi.org/10.1038/s43588-022-00349-3.

[2] B. Deng, P. Zhong, K. Jun, K. Han, C.J. Bartel, G. Ceder, CHGNet: Pretrained universal neural network potential for charge-informed atomistic modeling, (2023). https://doi.org/10.48550/arXiv.2302.14231.

[3] Riebesell, J., Goodall, R., Jain, A., Persson, K., & Lee, A. (Date TBD). Can machine learning identify stable crystals? [Preprint]. Matbench Discovery. Retrieved from https://matbench-discovery.materialsproject.org/preprint.

[4] N. Keriven, G. Peyré, Universal Invariant and Equivariant Graph Neural Networks, (2019). https://doi.org/10.48550/arXiv.1905.04943.

Reuse and Attribution

Why not just use a CPU

Modern science and engineering uses computing for a vast array of tasks. Some are basic such as text or image processing, while others are more specific to solving difficult math and physics equations using numerical methods¹. For much of computing the central processing unit (CPU) was the most important component for handling arithmetic, logic, controlling, and input/output (I/O) operations that are used to solve computational problems. The reason is tha the CPU is a very general purpose device and is very flexible in use. However, this comes with a downside: operations that might be simple to execute could actually be quite time consuming on a CPU.

So what is there to compare to, in other words, what other type of computing device could be used? A processing device that is tailored for specific type of operations. There are application specific integrated circuits(ASICs) which are gate arrays of transistors that can be programmed to behave in a specific manner that is enhanced for a specific type of application or task. Another is the Graphical Processing Unit (GPU) which is what I want to talk about today.

Ancient 🗞️

I'm writing this post because I wanted to capture the difference between GPUs and CPUs, but this post is probably 10 years too late, in other words, many have written post like this before.

The GPU vs. CPU

The fundamental difference between GPUs and CPUs is in their architecture and design philosophy [1-2]. The CPU is considered "general purpose" and is optimized for sequential processing. They typically have a few physical cores (typically 4-16 in consumer models) that are highly sophisticated, capable of handling a wide variety of operations efficiently.

The GPU, on the other hand, was originally designed for rendering computer graphics, which involves performing many similar calculations simultaneously [1-2]. This led to an architecture with hundreds or thousands of simpler cores, each capable of executing the same instruction on different data points simultaneously. This approach, known as Single Instruction, Multiple Data (SIMD), makes GPUs extraordinarily efficient at parallel processing [1].

Why can GPUs be better for Scientific Computing?

We can go through various points to understand why a GPU is better for many scientific² computing applications:

1. Parallel Processing: Many scientific algorithms, such as matrix operations, Monte Carlo simulations, and numerical integrations, can be parallelized. GPUs excel at these tasks, potentially offering orders of magnitude speedup compared to CPUs.

2. Memory Bandwidth: GPUs typically have much higher memory bandwidth than CPUs. This is crucial for data-intensive computations where the bottleneck is often moving data rather than performing calculations.

3. Floating-Point Performance: Modern GPUs are optimized for floating-point arithmetic, which is essential in many scientific computations. They can perform these operations much faster than CPUs, especially when dealing with single-precision numbers.

4. Energy Efficiency: For the same computational power, GPUs often consume less energy than CPUs, which can be a significant factor in large-scale scientific computing environments.

5. Cost-Effectiveness: In terms of raw computational power per dollar, GPUs often outperform CPUs, making them an attractive option for budget-constrained research projects.

When is a GPU not better?

Because GPUs are specialized computing devices, they have some downsides:

1. Sequential Tasks: Algorithms that can't be easily parallelized often perform better on CPUs due to their higher clock speeds and more sophisticated cores.

2. Complex Control Flow: Tasks involving many conditional statements or complex branching logic are generally better suited to CPUs.

3. Small Data: The overhead of transferring data to and from the GPU can outweigh the benefits for small data.

4. I/O-Bound Tasks: For operations that are limited by input/output speeds rather than computational power, the GPU's advantages may be negated.

Some challenges with GPU programming

For at least 10 years I've tried to get reasonable familiar with GPU programming, mostly with CUDA, which is Nvidia's GPU programming language. There are other GPU programming languages such as OpenCL and HIP that are similar to CUDA. The issue is they are less adopted and seem more difficult. But in general I have found GPU programming to be a tour-de-force, that is it is always more difficult then I anticipate. Here are some reasons:

1. Different Programming Model: GPU programming requires thinking in terms of parallel execution, which can be counterintuitive or unnatural.

2. Memory Management: GPU memory is limited and efficient GPU programming requires careful management memory types and data transfer between CPU and GPU.

3. Debugging Tools: I've never been good with debugging tools, but GPU code is much harder to debug.

4. Optimization: Achieving peak performance usually requires understanding specific GPU architectures.

Example PyCUDA code

To illustrate GPU programming in a context relevant to materials science, let's consider a simple example of calculating the Lennard-Jones potential, which is commonly used to model interatomic interactions in molecular dynamics simulations. Here's a CUDA implementation using Python and PyCUDA [3]:

import numpy as np
import pycuda.autoinit
import pycuda.driver as drv
from pycuda.compiler import SourceModule


# GPU implementation
cuda_module = SourceModule("""
global void lennard_jones_gpu(float potential, 
                                  float positions, 
                                  int num_atoms, 
                                  float epsilon, 
                                  float sigma)
{
    const int i = threadIdx.x + blockDim.x * blockIdx.x;
    if (i < num_atoms) {
        float pote = 0.0f;
        float3 pos_i = make_float3(positions[3i], 
                                   positions[3i+1], 
                                   positions[3i+2]);


        float sigma6 = powf(sigma, 6);
        float sigma12 = sigma6 * sigma6;
        for (int j = 0; j < num_atoms; j++) {
            if (i != j) {
                float3 pos_j = make_float3(positions[3j], 
                                           positions[3j+1], 
                                           positions[3j+2]);
                float3 r = make_float3(pos_i.x - pos_j.x, 
                                       pos_i.y - pos_j.y, 
                                       pos_i.z - pos_j.z);
                float r2 = r.xr.x + r.yr.y + r.zr.z;
                float r6 = r2r2r2;
                float r12 = r6r6;


                pote += 4.0f * epsilon * ((sigma12/r12) - (sigma6/r6));
            }
        }
        potential[i] = pote;
    }
}
""")


lennard_jones_gpu = cuda_module.get_function("lennard_jones_gpu")


# Parameters
num_atoms = 1000
epsilon = 1.0

sigma = 1.0 



# Arrays
positions = np.random.rand(num_atoms * 3).astype(np.float32)
potential = np.zeros(num_atoms).astype(np.float32)


# Call the kernel
lennard_jones_gpu(
    drv.Out(potential), drv.In(positions),
    np.int32(num_atoms), np.float32(epsilon), np.float32(sigma),
    block=(256, 1, 1), grid=((num_atoms + 255) // 256, 1)
)

The keyword args block and grid specify to the kernel how to perform the parallelization (i.e., how to decompose the arrays into parallel computes). This example demonstrates how we can use GPU computing to calculate the Lennard-Jones potential for a system of atoms, which is a fundamental operation in many molecular dynamics simulations. The GPU's parallel processing capabilities allow us to compute the potential for many atoms simultaneously, which can lead to significant speedups compared to CPU-based calculations, especially for larger systems.

Benchmarking CPU vs GPU: Lennard-Jones Potential Calculation

To truly appreciate the performance gains in using GPUs for computational materials science, I put together a jupyter notebook which compares the execution time of Lennard-Jones potential calculations on both CPU and GPU. I think it helps demonstrate how much speed-up can be had with GPU hardware for the right problem.

Number of Atoms	Single CPU Time (s)	GPU Time (s)	Speedup
100	0.1490	0.0008	~1830x
500	2.7130	0.0009	~3100x
1000	10.8806	0.0009	~11000x
2000	43.7833	0.0127	~3400x
5000	274.5926	0.0135	~20000x

One pecularity is that the 2000 atom system shows a drop in speed-up, not sure what causes this. It could be the block and grid settings, but then for 5000 atoms we see the speed-up factor jumps again.

Attention

The speed-up will be GPU hardware dependent and I'm not comparing it to multi-core CPU performance. The hardware used to generate the speed-up table was a RTX A4000 GPU so other cards may be better or worse, just depends. It's also important to note that this is a simplified example. Practical atomistic simulations often include additional routines and more complex force fields.

Summary

I think this benchmark demonstrates why GPUs can be so powerful for scientific computing tasks like MD simulations. As the number of atoms in the system increases, the GPU's ability to perform many calculations in parallel leads to increasingly dramatic performance improvements over an single CPU. Yes MPI and OpenMP and can be used to provide great speed-ups.

As I continue to explore the use of GPUs in materials science and AI, I'm pretty excited by the possibilities these performance improvements open up. With GPUs, we can simulate larger systems for longer timescales, potentially leading to new insights in materials research and drug discovery.

Footnotes

---

1. In numerical methods, the problem at hand is represented in a discretized fashion that is solved using basic arithmetic operations. However, it is possible to represent the problem using symbolic formalism and then the symbolic form is manipulated based on rules and conditions to arrive at a solution. ↩

2. I'm using scientific computing applications as a stand in for both physical science and ML/AI computing. ↩

References

[1] D.B. Kirk, W.W. Hwu, Programming Massively Parallel Processors, Third Edition: A Hands-On Approach, CreateSpace Independent Publishing Platform, 2017. URL.

[2] J. Nickolls, W.J. Dally, The GPU Computing Era, IEEE Micro 30 (2010) 56–69. DOI.

[3] A. Klöckner, N. Pinto, Y. Lee, B. Catanzaro, P. Ivanov, A. Fasih, PyCUDA and PyOpenCL: A scripting-based approach to GPU run-time code generation, Parallel Computing 38 (2012) 157–174. DOI.

---

Reuse and Attribution

I now know why Boron Carbide shatters

Hard high-temperature ceramics is an area I've been involved with throughout my career. One of these materials is Boron Carbide. The B₄C polymorph is one of the hardest known materials and is incredibly lightweight due to its constituents and density. These attributes make B₄C a leading candidate for personal body armor, vehicle shielding, and other high-impact applications. Its remarkable properties originate from its rhombohedral crystal structure that consists of a 12-atom icosahedra that are connected by stiff three-atom chains.

Figure 1. Unit cell of B₄C with 12-atom icosahedra & C-B-C fragment (Materialscientist, CC BY-SA)

Despite its exceptional hardness boron carbide exhibits a critical weakness under high-speed impact: it can fail catastrophically through amorphization, where the ordered crystal collapses into a glassy/disordered state. This phase change results in loss of outstanding hardness and strength. For a long time the atomic-scale mechanisms behind this phenomenon weren't really well characterized or understood. I remember shock physics molecular dynamics studies trying to look at this and it was never quite clear what the incipient (i.e. onset/preceding) steps were to amorphiziation. From a modeling perspective, a lot of it is related to the ability of interatomic potentials to capture the atomic-scale process. That seems to have changed. A recent study by Ghaffari et al. [1] used machine-learned interatomic potentials (MLIPs) and MD simulations to reveal how the crystal orientation, polymorphism, and velocity impact the amorphization failure process.

One of the nails in the coffin, so to speak, is shown in Figure 2, which reveals thin amorphous bands forming at ~45° relative to the impact direction, demonstrating the shear-driven nature of failure. Previous MD simulations could not reproduce this experimentally observed phenomenon.

Figure 2. Amorphous bands at ~45° relative to impact (Fig. 4 from Ghaffari et al. 2025)

The Paradox of Boron Carbide

The exceptional hardness of most carbides, like B₄C (Vickers ≥ 30 GPa)¹, originate from the covalent bonds, and in the case of B₄C its icosahedra and three-atom chains. However, under dynamic loading beyond its Hugoniot Elastic Limit (HEL), which is the transition condition where a material goes from elastic to plastic deformation under shock loading and for B₄C it is experimentally around 20 to 25 GPa. So under shock loading B₄C softens instead of hardening, exhibiting glass-like mechanical behavior [1]. Amorphization degrades the structure of the material by breaking up icosahedral cages, resulting in loss of long-range order, and therefore sharp reduction in shear strength and hardness. In experiments they have observed narrow amorphous bands forming along shear planes [2-4], but these of course lacked temporal and spatial resolution to identify the atomic processes responsible for the amorphization.

New kid on the block: MLIPs

I've posted a bit on this blog about machine-learned interatomic potentials (MLIPs) because they have shifted the type of classical atomistic simulations that can be now, mostly in terms of chemistries simulated, but with the work by Ghaffari et al. [1] its pushed this to new physics regimes and mechanisms. From what I recall, simulating shock physics amorphization in carbides has been a challenge with classical interatomic potentials. I'm sure some have tried AIMD but the spatial scale would always be the limiting factor, i.e., nearly impossible to capture the amorphous band formation. Force fields such as ReaxFF have been tried [5] but fail to reproduce the observed experimental features. Ghaffari et al. [1] address this by developing a machine-learned interatomic potential trained on DFT-calculated structures for four dominant polymorphs, amorphous phases, and high-pressure states. Using DeePMD-kit, the resulting model achieved high accuracy and enables MD simulations of systems scaled up to ~450K atoms.

tl;dr

Simulations using their trained boron carbide MLIP model successfully capture shear-induced amorphous band formation for the first time, directly matching experimental observations and providing atomistic process responsible for the amorphization.

So what did they find?

I mentioned earlier in Figure 2 that the amorphous bands formation being directly observed in the MD simulations. This was one of the key results that they were able to achieve, but some other important details were also revealed.

1. Crystal Orientation Strongly Influences HEL

The HEL depends strongly on the orientation of the three-atom chain relative to the impact direction. The lowest HEL values, between 18 and 22 GPa, occurred when the chains were either aligned parallel or perpendicular to the impact. In contrast, orientations such as 45° and 60° showed the highest HEL values, between 32 and 38 GPa [1]. This seems to overturn the earlier assumptions from ReaxFF-based simulations [5] that aligning the orientation with the least compliance with the impact direction would maximize strength. But what the results in the paper indicate is the failure is driven primarily by shear stress rather than compression. Amorphous bands consistently form at planes oriented approximately 45° from the impact direction, corresponding to the planes of maximum shear stress.

2. Polymorphism Governs Toughness

Boron carbide polymorphs differ subtly in atomic arrangement but behave very differently under shock loading, as highlighted in Figure 3. The paper directly compares B₁₂(CBC), which contains a C–B–C chain, and B₁₂(CCC), which has a C–C–C chain, revealing how polymorphism controls toughness. B₁₂(CBC) exhibits the highest Hugoniot Elastic Limit (HEL) and the greatest resistance to amorphization, surviving shock by forming amorphous bands that expand in volume. In contrast, B₁₂(CCC) has a much lower HEL and fails via irreversible structural collapse and significant volume shrinkage. This difference is directly linked to the atomic arrangement: in B₁₂(CBC), the middle boron atom temporarily migrates into a cage space during compression and returns upon unloading, resulting in reversible chain bending and recovery of the lattice structure. For B₁₂(CCC), the central carbon atom instead forms permanent bonds with nearby icosahedra due to the higher bond dissociation energy of C–B bonding (448 KJ/mol) relative to B–B bonding (297 KJ/mol), leading to irreversible bonding, collapse, and densification. Figure 4 provides side-by-side atomistic snapshots of these chain mechanics, illustrating the reversible bending and recovery in CBC versus the irreversible bonding and collapse in CCC. Figure 5 further strengthens this mechanistic explanation by showing the cage-space migration pathway that enables reversibility in CBC.

Figure 3. Reversible vs Irreversible Chain Bending (Fig. 8 from Ghaffari et al. 2025)

Figure 4. Reversible vs Irreversible Chain Bending (Fig. 10 from Ghaffari et al. 2025)

Figure 5. Cage-Space Migration (Fig. 11 from Ghaffari et al. 2025)

3. Two Distinct Regimes of Amorphization

The study also identifies two fundamentally different deformation regimes controlled by impact velocity [1]. In Figure 6 the response to 2.5, 3.0, and 4.0 km/s shocks is shown and it establishes distinct band-dominated versus collapse-dominated regimes. At 2.5 km/s, localized amorphous bands dominate, and the volume inside these bands expands by ~6% relative to the surrounding crystal lattice, which then suppresses crack initiation (i.e. extrinsic toughening). At 3.0 km/s, amorphous banding and partial structural collapse both coexist. But at 4.0 km/s and above, widespread collapse of icosahedral units occurs, producing significant volume shrinkage and creating voids that promote crack growth. The preservation or destruction of icosahedral motifs determines whether a region undergoes expansion, local compressive shielding, or irreversible densification.

Figure 6. Atomic Volume Variation

Implications for Materials Design

The results from Ghaffari et al. [1] provide an atomic-level framework for designing boron carbide that is more resistant to amorphization. Orienting grains to align their stiffest chains with the planes of maximum shear resistance, approximately 45°, increases HEL and delays failure. Favoring polymorphs containing C–B–C chains during synthesis improves toughness by enabling reversible chain bending and reducing susceptibility to catastrophic collapse. Furthermore, the demonstrated accuracy of machine-learned interatomic potentials shows that AI-driven simulations are now a critical tool for predicting material behavior under extreme conditions and guiding the design of next-generation superhard ceramics.

Footnotes

---

1. In contrast to Boron Carbide, Tungsten carbide is another hard and refractory material, but Tungsten carbide density is ~7x (15.7 g/cm³) that of Boron Carbide. I believe the hardness of Boron Carbide is also slightly higher than that of Tungsten carbide. ↩

References

[1] K. Ghaffari, S. Bavdekar, D.E. Spearot, G. Subhash, Influence of Crystal Orientation and Polymorphism on the Shock Response of Boron Carbide, SSRN Preprint, (2025). https://ssrn.com/abstract=5186595.

[2] G. Parsard, G. Subhash, P. Jannotti, Amorphization-Induced Volume Change and Residual Stresses in Boron Carbide, Journal of the American Ceramic Society. 101 (2018) 2606–2615. https://doi.org/10.1111/jace.15442.

[3] M. Chen, J.W. McCauley, K.J. Hemker, Shock-Induced Localized Amorphization in Boron Carbide, Science. 299 (2003) 1563–1566. https://doi.org/10.1126/science.1080067.

[4] K.M. Reddy, P. Liu, A. Hirata, T. Fujita, M.W. Chen, Atomic Structure of Amorphous Shear Bands in Boron Carbide, Nature Communications. 4 (2013) 2483. https://doi.org/10.1038/ncomms3483.

[5] A.A. Cheenady, A. Awasthi, M. DeVries, C. Haines, G. Subhash, Shock Response of Single-Crystal Boron Carbide Along Orientations with the Highest and Lowest Elastic Moduli, Physical Review B 104 (2021) 184110. https://doi.org/10.1103/PhysRevB.104.184110

---

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

Dual Numbers

A while back when I was doing some exploration of writing a simple NN code to improve my understanding of neural networks and deep learning in general, I came across dual numbers. They're a type of number that generalizes the concept of real and complex numbers. But what makes them so interesting is that they can encode both a function value and its derivative in a single number. This means that we can use them to simplify the calculation of derivatives and solve complex problems efficiently.

So how does one think of dual numbers? What's the difference between a dual number and a complex number? One way to think about dual numbers is that they consist of two parts: a scalar part and a skew part. The scalar part is just a regular real number, while the skew part is a multiple of a new number, often denoted as $\epsilon$, that satisfies the property $\epsilon^2=0$. This means that every dual number can be written as $a+b\epsilon$, where $a$ and $b$ are real numbers.

What's most interesting is that the skew part of a dual number is that it provides an approximation of the first derivative of a function evaluated at a particular point. By using the dual number representation of the function at that point, one can calculate both the function value and its derivative in one shot.  One reason dual numbers have applications in deep learning is that algebra on dual numbers provides the chain rule for calculus, therefore they can be used to compute derivatives of complicated functions involving multiple variables and interdependencies.

As an example, say I want to evaluate the function $f(x)=x^2+2x$ at $x=3$. The dual number representation of $f(3)$ is $f(3+\epsilon)=f(3)+f'(3)\epsilon$, where $f'(x)=\frac{df(x)}{dx}$. We can compute $f(3)$ directly as $f(3)=3^2+2\cdot3=9+6=15$. To compute $f'(3)$, we can take the derivative of $f$ with respect to $x$: $f'(x)=2x+2$. Evaluating this at $x=3$, we get $f'(3)=2\cdot3+2=8$. Therefore, the dual number representation of $f$ at $x=3$ is $15+8\epsilon$.

One of the benefits of dual numbers is the derivative of the composition of two functions, $f(g(x))$ requires only the derivatives of the individual functions. Specifically, if $f(x)$ and $g(x)$ are two functions, then the dual number representation of their composition $f(g(x))$ is $(f(g(x)), f'(g(x))g'(x))$. This is especially useful when dealing with complex functions involving multiple variables and complicated interdependencies.

Dual numbers are actually a useful mathematical concept because they have practical applications in a wide range of fields. It's pretty cool that one can encode function values and derivatives in a single number, which makes it possible to simplify the calculation of derivatives and solve complex problems efficiently. On my computational blog, I have an example using dual numbers to calculate the derivative of an interatomic potential, Dual Numbers Pluto blog.

Reuse and Attribution

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

DFT is gettings its share of AI

I made a post on LinkedIn the other day sharing a recent preprint[1] by Microsoft AI that introduces a significant advancement in exchange-correlation functionals (XC) for DFT. This, to me, seemed like an important piece of news given that much of the data used for training the recent explosion in foundational interatomic potentials for atomistics uses massive DFT datasets.

The biggest achievement seems to be the ability for Microsoft AI's DFT XC functional, called Skala, to achieve near chemical accuracy (i.e., 1 kcal/mol) for molecular systems while doing it with reasonable computational cost. What this means is that we may get a new rung added to the chemistry Jacob's ladder of accuracy for DFT. I wanted to illusrate this so I made a diagram of the chemistry Jacob's ladder with the new Skala XC added as shown in Figure 1. I got the idea to annotate it from Keiran Rowell's blog[2].

The idea of Jacobs Ladder in chemistry, I think, was first introduced by J. Perdew in 2013(?) [3], I believe. The primary idea is that like the in the biblical story, Jacob dreams of a ladder reaching from earth to heaven and God standing above it promising Jacob divine blessing and protection. In the spirit of this, we the scientist dream of reaching chemical accuracy using our ladder of DFT approximations. Why is this important? Well, chemical accuracy with DFT would enable accurate enough predictions in reasonable computational time to improve discovery, design, and understanding of chemicals and materials. In particular understanding reaction mechanisms and kinetics would be much more tractable if chemical accuracy is reached.

Figure 1. Computational chemistry Jacob's ladder for DFT with the Microsoft AI Skala XC added.

DFT XC Overview

First, let me briefly summarize some basics of density functional theory (DFT) to explain why the exchange-correlation (XC) functional is so important. The Nobel-worthy Hohenberg–Kohn theorems and the formulation of the Kohn–Sham equations established that the total ground-state (i.e., 0 K) energy of a many-body electronic system can be written as:

$ \begin{equation} E[\rho] = T_s[\rho] + V_{ext}[\rho] + J[\rho] + E_{xc}[\rho] \end{equation} \label{eq:dft_energy} $

The first three terms are kinetic energy, external potential energy, and Coulombic repulsion energy. The last term is the XC energy which is the only term we don't know exactly and is crucial because it encapsulates the many-body quantum mechanical effects of the electrons whereas the others are single-particle terms.

Throughout the successful history of DFT, different approximations have been used for the XC term. The first was local density approximation (LDA), with the idea being that the XC energy is simply a function of the local density of the electrons. For all its simplicity, LDA did work reasonably for some systems but shortcomings were quickly discovered.

Then the generalized gradient approximation (GGA) was introduced, where the XC energy is now not only a function of the local density but also the gradient of the density. This introduced the idea of non-local XC functionals (i.e., the energy now depends on how the density changes with position, rather than just at that position).

After GGA, what followed was a "flavor-of-the-month" approach toward improving XC functions. There were hybrid functionals that mixed DFT XC with the Hartree-Fock exchange. Then meta-GGA was introduced with the idea of adding higher-order derivatives to the density. Add to that expensive XC functionals that include perturbation quantum chemistry methods. All have their improvements and shortfalls. You can see all the major players in the rungs of the ladder in Figure 1 showing how different XC functionals improve the range of accuracy as you climb up the ladder.

These are just the formalisms; there are also ad-hoc tweaks and specializations that get added based on your domain use and expertise. The limitations of such approaches are:

1. Handcrafted Features: Most functionals rely on fixed analytic forms and known constraints.
2. Non-locality Deficiency: Electron correlation is inherently non-local; local/semi-local functionals cannot fully capture this.
3. Empirical Tuning vs. Physical Justification: Empirical functionals may generalize poorly or violate constraints.
4. Slow Progress at Higher Rungs: Despite decades of work, hybrid/double-hybrid methods still fall short of universal chemical accuracy (error < 1 kcal/mol).

Skala Neural XC

The efforts by Microsoft AI represent a significant leap in XC functional design. Unlike traditional functionals that have relied on hand-crafted mathematical forms, Skala incorporates deep learning (i.e., Neural Networks) to achieve near-chemical accuracy while maintaining computational efficiency comparable to meta-GGA functionals[1]. Skala is particularly impressive because it navigates the trade-off between accuracy and computational cost. The Microsoft AI team has designed it to bridge the gap between semi-local functionals (fast but less accurate) and hybrid/double-hybrid functionals (accurate but computationally expensive).

Architecture and Design Philosophy

As mentioned, Skala XC is a neural network architecture trained to learn the non-local electron density interactions without requiring the full computational burden of exact exchange calculations. Like meta-GGA, it starts with seven semi-local density-derived features, but the design then employs what the authors call a "coarse-fine grid structure" that captures long-range density correlations consistent with multipole-like behavior.

For the exchange, the authors use LDA but incorporate a neural enhancement factor $f_\theta$:

$ \begin{equation} E_{xc}^\theta[\rho] = -\frac{3}{4}\left(\frac{6}{\pi}\right)^{1/3} \int \left(\rho_\uparrow(r)^{4/3} + \rho_\downarrow(r)^{4/3}\right) f_\theta[x[\rho]](r) \, dr \label{eq:skala_xc} \end{equation} $

Not knowing too much about XC design, this does seem clever though because $f_\theta$ operates on a feature vector $x[\rho]$ that encodes both local and non-local density information while maintaining computational efficiency (would like to understand this but above my head for now). This is reminiscent of delta-learning in atomistic ML, though Skala performs a direct functional approximation rather than a residual correction atop an existing functional.

What data was used?

The Microsoft AI team used a quantum chemistry dataset that seems only they could have curated, something like 150,000 data points spanning thermochemistry, conformational energies, noncovalent interactions, and ionization potentials. This was all done with gold-standard wavefunction methods like CCSD(T) to generate the data, 🤯. Their MSR-ACC/TAE dataset alone includes about 80,000 total atomization energies with errors below 1 kcal/mol, which is just staggering when you consider the computational expense of generating such reference data. The datasets alone might be useful for others.

Peformance Metrics

I'm not too familiar with benchmarks in this space, but they mention the challenging W4-17 benchmark and Skala achieves a mean absolute error (MAE) of around 1.0 kcal/mol, outperforming more established functionals. On the GMTKN55 benchmark, it scores a WTMAD-2 of 3.89 kcal/mol, which puts it in competition with the best hybrid functionals while requiring significantly less computational resources.

Another thing that is super interesting is that for large atomic numbers that are out of distribution for the trained Scala XC, They are showing that it maintains very good accuracy. This is usually a failure point for data-driven ML models where out of distribution data is not well handled particularly well [4]. The authors also method that the self-consistent field convergence is stable, crucial for practical application in DFT --- since, unlike post-SCF models, Skala is trained and used self-consistently.

Redefining Jacob's Ladder?

So going back to the Jacob's ladder of DFT functionals in Figure 1. Whether Skala constitutes a sixth rung or a bypass of the traditional ladder remains obviously open to interpretation, but it clearly represents a shift. I'm not a veteran in this space so hard for me to be definitive, but it probably represents some departure to non-locality from semi-local features to effectively bypass the traditional way of doing things.

It will be interesting to see from the "Kings of DFT" to see how they view Skala XC. Based on the Microsoft AI preprint this functional(s) is showing systematic improvement through use of high-quality data and training procedures, and the authors claim it can encode known physics through appropriate constraints and design. On the other hand, they lack the theoretical transparency and interpretability that traditional functional forms offer will probably be a critiqued, although existing functionals are not perfect either.

I'm just interested to see how this gets used for the other down stream applications in solid-state physics and materials science. Will we get DFT training sets for Materials Project structures that range in the 1 million of structures and all have chemical accuracy? If so that could in turn make these foundation models even better for materials discovery and classical atomistic modeling.

---

References

[1] G. Luise, C.-W. Huang, T. Vogels, D.P. Kooi, S. Ehlert, S. Lanius, K.J.H. Giesbertz, A. Karton, D. Gunceler, M. Stanley, W.P. Bruinsma, L. Huang, X. Wei, J.G. Torres, A. Katbashev, B. Máté, S.-O. Kaba, R. Sordillo, Y. Chen, D.B. Williams-Young, C.M. Bishop, J. Hermann, R. van den Berg, P. Gori-Giorgi, Accurate and scalable exchange-correlation with deep learning, (2025). DOI.

[2] K. Rowell, An Ersatz Ansatz, Blog (2023). https://keiran-rowell.github.io/guide/2023-04-12-compchem-methods-basics (accessed June 21, 2025).

[3] J.P. Perdew, Climbing the ladder of density functional approximations, MRS Bull. 38 (2013) 743–750. DOI.

[4] K. Li, A.N. Rubungo, X. Lei, D. Persaud, K. Choudhary, B. DeCost, A.B. Dieng, J. Hattrick-Simpers, Probing out-of-distribution generalization in machine learning for materials, Commun Mater 6 (2025) 1–10. DOI.

---

Reuse and Attribution