Before delving into development of complex deep learning-based regression models, it is worthwhile to investigate whether simple linear regression models, crafted with physically-motivated features, may provide sufficient accuracy and to provide a benchmark against more complex models. Figure indicates that single variable linear regression models based on simple descriptors can be derived with R 2 as high as ∼0.3. Here, max(R ionic-ov ) denotes the maximum in ionic radii overlap (using default pymatgen 44 values) between any atom in the crystal structure and the diffusing O atom along the vector of the linearly interpolated path, p = rs + η(r e -rs ), η ∈ [0, 1], between end and start site coordinates, re and rs , respectively, with ||p|| 2 denoting the path length. Some interesting qualitative insights are observed, such as an intuitive lower bound to E A for increasing max(R ionic-ov ) and ||p|| 2 . However, these simple models are generally not accurate enough (all have R 2 ≤ 0.35) for downstream diffusivity calculations in materials screening applications. d 2 GNN surrogate model. We now seek a better surrogate model for Equation (2), i.e., a single-shot predictor of all NEB energies for a vacancy migration event from one crystallographic site to another. The model should facilitate rapid inference of NEB energies for all symmetrically inequivalent paths in not just one structure but the potentially tens of thousands oftentimes examined in high-throughput screening studies. We therefore seek to construct a model, parameterized by weights θ, that relies only on the relaxed host crystal structure, C h , and the path tuple p as input,