{'fraction_non_alphanumeric': 0.03158707752949636, 'fraction_numerical': 0.013476442858284284, 'mean_word_length': 4.70985401459854, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spherical Channels for Modeling Atomic Interactions

Modeling the energy and forces of atomic systems is a fundamental problem in computational chemistry with the potential to help address many of the world's most pressing problems, including those related to energy scarcity and climate change. These calculations are traditionally performed using Density Functional Theory, which is computationally very expensive. Machine learning has the potential to dramatically improve the efficiency of these calculations from days or hours to seconds. We propose the Spherical Channel Network (SCN) to model atomic energies and forces. The SCN is a graph neural network where nodes represent atoms and edges their neighboring atoms. The atom embeddings are a set of spherical functions, called spherical channels, represented using spherical harmonics. We demonstrate, that by rotating the embeddings based on the 3D edge orientation, more information may be utilized while maintaining the rotational equivariance of the messages. While equivariance is a desirable property, we find that by relaxing this constraint in both message passing and aggregation, improved accuracy may be achieved. We demonstrate state-of-the-art results on the large-scale Open Catalyst 2020 dataset in both energy and force prediction for numerous tasks and metrics.


Introduction
Modeling the properties of atomic systems is a foundational challenge in computational chemistry and critical to advancing technologies across numerous application domains. Notable applications include drug discovery [31,38] and the design of new catalysts for renewable energy storage to help in addressing climate change [47,32]. For catalyst discovery, new materials are currently evaluated using Density Functional Theory (DFT) that can estimate atomic energies and forces but is computationally very expensive; taking hours or days to evaluate a single material. Machine Learning (ML) has the potential to approximate DFT and dramatically speed up these calculations, allowing for high throughput screening of new materials to help address some of the world's most pressing challenges.
Our goal is to approximate DFT calculations using ML. An ML model takes as input a set of atom positions and their atomic numbers. As outputs, a model calculates the structure's energy (or other properties) and the per-atom forces, i.e., the forces exerted on each atom by the other atoms. A common approach to this problem is to use Graph Neural Networks (GNNs) [16,46] where each node represents an atom and the set of nearby atoms as edges [36,14,21,34,37,44,30,23].
A key challenge in network design is balancing the use of model constraints. When modeling atomic systems, a common constraint is SO(3) rotation equivariance [42,3,1,41,33,35], i.e., if the atomic system is rotated, the energies should remain constant and the atomic forces should similarly rotate. While this provides strong priors on the model to help in generalization, especially for smaller datasets [7,31], it can result in limiting the expressiveness of the network due to restrictions on non-linear transformations for equivariant models [3,41,33,35], or limiting interactions to pairs [34,37,44], triplets [23] or quadruplets of atoms [27,24,12] for invariant models. Alternatively, a non-equivariant model [20] can provide more freedom to the model, but lead to the model needing to learn approximate equivariance through methods such as data augmentation [39]. To draw an analogy with image detection, the use of a CNN [26] provides translation equivariance and removes the need for the network to learn how to detect the same object at different locations. However, most CNNs are not equivariant to scale or rotation [43,42], but are still effective in learning approximate equivariance through rotation and scale diversity in the training data. What is the analogous balance of constraints for modeling atomic systems?
In this paper, we propose a GNN [16,46] that balances the use of model constraints to aid in generalization while providing the network with the flexibility to learn accurate representations. We introduce the Spherical Channel Network (SCN) that explicitly models relative orientations of all neighboring atoms; information that is critical to accurately predicting atomic properties. Each node's embedding is a set of functions defined on the surface of a sphere (S 2 → R). The functions are represented using spherical harmonics, similar to approaches that build strictly equivariant models [41,1,3,10]. During message passing, angular information between atoms is conveyed by rotating or steering [5] the embeddings based on each edge's orientation [40,27]. We identify an expanded set of spherical harmonic coefficients that are invariant to rotation, which can provide rich information while maintaining a message's rotation equivariance. In addition, we demonstrate that if the equivariance constraint is relaxed, improved performance can be achieved by using additional coefficients. We further improve the expressivity of the network by performing a pointwise non-linear function, which is only approximately equivariant, on the embeddings during message aggregation.
We demonstrate our Spherical Channel Network on the large-scale Open Catalyst 2020 (OC20) dataset [6], which contains atomic structures useful for numerous applications important to addressing climate change. State-of-the-art results are achieved for atomic force and initial structure to relaxed energy prediction with improvements of 8% − 11%. In addition, we demonstrate our model is more sample efficient as compared to other state-of-the-art models.

Approach
Given an atomic structure with n atoms, our goal is to predict the structure's energy E and the per-atom forces f i for each atom i ∈ n. These values are estimated using a Graph Neural Network (GNN) [16,46] where each node represents an atom and edges represent nearby atoms. As input, the network is given the distance d ij between atoms i and j, and each atom's atomic number a i . The neighbors N i for an atom i are determined using a fixed distance threshold, or by picking a fixed number of closest atoms.

Node Embeddings
The angular or relative orientation information between atoms is critical to accurately modeling atomic properties such as energies and forces [23,24]. Inspired by this, our node embedding models the angular information from all neighboring atoms using spherical functions. Each node i's embedding s i is a set of C functions or channels represented on a sphere (S 2 → R), whose argument is a 3D unit vector indicating the orientation. That is, s ic (r) is the value of channel c for node i for some orientationr ∈ R 3 . Since the spherical channels contain orientation information over the entire sphere, the network may reason about directional information for all neighboring atoms and not just atom pairs, triplets, etc. Each spherical channel c may be represented using several different approaches, such as discrete 2D grids sampled over the sphere, or spherical harmonics. As we describe later, we use spherical harmonics due to their property of being SO(3)-equivariant (3D rotation equivariant).
Using spherical harmonics, a function s c (r) defined on the sphere is represented by a set of weighted spherical harmonics Y m l with s c (r) = l,m x lmc Y m l (r) where l and m are the degree and order of the basis functions with m ∈ [−l, l] and l ∈ L. x lmc are the spherical harmonic coefficients, which we refer to as the coefficients of the spherical channel. For every degree l there exists 2l + 1 spherical harmonics (see Figure 1(a) for an illustration of functions for l ≤ 4), which results in a total of (L + 1) 2 basis functions and coefficients up to degree L. Therefore for a maximum degree of L with C channels, x has size (L + 1) 2 × C.
The spherical channels are updated by the GNN through message passing for K layers to obtain the final node embeddings S (K) . S (0) is initialized from an embedding based on the atom's atomic number a i for l = 0 coefficients and the l = 0 coefficients are set to zero. The nodes' embeddings are updated by first calculating a set of messages m ij for each edge, which are then aggregated at each node. Finally, the energy and forces are estimated from S (K) . We describe each of these steps in turn.

Message Passing
Given a target node t and its neighbors s ∈ N t we want to update the embeddings x (k) t at iteration k ∈ K. The embeddings x (k) t are a set of spherical harmonic coefficients indexed by their degree l, order m and channel c. An important and useful property of spherical harmonics is the function represented by the coefficients is steerable [11,5], i.e., it can be rotated in 3D space using a linear transformation of the coefficients. Specifically, for a 3D rotation matrix R there exists a matrix called a Wigner D-matrix D l of size (2l + 1 × 2l + 1) that rotates the coefficients of degree l by the rotation R. If x ilc are all coefficients of degree l for node i and channel c, then for any orientationr and rotation R there exists a D l such that: (1) When calculating the message m st from atom s to atom t, we want to use the information contained in both x (k) given the context of the edge's orientationd st ,d st = d st /|d st |. We do this by rotating the embeddings by R st for which R stdst = [0, 0, 1] , i.e., the direction of the edge's unit vectord st is aligned with the z-axis. Thus when the messages are calculated, the problem is simplified since the network can assume the atoms are aligned along the z-axis. After calculating the messages using a neural network F e , they are rotated back to the global coordinate frame: where the matrices D st and D −1 st perform the rotation on the embeddings' coefficients corresponding to the 3D rotation matrices R st and R −1 st respectively. To simplify notation, we assume the coefficients in x s across degree l and order m are flattened, and is a block diagonal matrix constructed from the set of Wigner D-matrices across degrees, l ∈ L.
In addition to the rotated embeddings, F e is also provided input edge information that is invariant to rotations; |d st | is the magnitude of the distance between atoms s and t, and a s and a t their respective atomic numbers. Note that the rotation matrix R st is not unique, since the roll rotation around the vector d st is not specified and is randomly sampled during training. The implications of this are discussed below.
The message function F e is computed using a neural network as illustrated in Figure 2(a). The atomic numbers a s and a t are used to look up two independent embeddings of size E = 128, and a set of 1D basis functions are used to represent |d st | using equally spaced Gaussians every 0.02 Å from 0 to the 8 Å with σ = 0.04 followed by a linear layer with E outputs. Their values are added together and passed through a neural network to produce a vector of size H. The rotated embeddings D st x (k) s and D st x (k) t are concatenated and passed through a single layer neural network to produce another vector of size H. These two are multiplied together to combine the edge information with the information contained in the rotated embeddings. Two more fully connected layers are performed with each followed by a SiLU non-linear activation function [19]. Finally, a single linear layer is used to expand the output to the original size of the embeddings.

Message Equivariance
For our message function (Equation (2)) to be equivariant, the output of F e needs to be invariant to global rotations. In general, this is not the case since the matrix D st is not unique, i.e., the roll rotation around the vector d st is not specified, and in practice it is randomly chosen. The roll rotation corresponds to a rotation about the z-axis after rotating x (k) s and x (k) t by D st . Due to this, all of the rotated coefficients will vary based on the random rotation chosen, except the m = 0 coefficients that are symmetric about the z-axis; see Appendix C and the m = 0 bases highlighted in green in Figure 1(a). If rotation equivariance is desired for messages, the input and output coefficients to F e can be restricted to only those for which m = 0. Since the other inputs |d st |, a s and a t to F e are also invariant to rotations, F e is invariant to rotations if only m = 0 coefficients are used. The resulting messages m st are equivariant to rotations once they are rotated back to the global coordinate frame using D −1 st . While equivariance is a desirable property, as we demonstrate later, using m ∈ [−1, 1] coefficients can help improve performance even though equivariance is violated.
If we choose to use coefficients beyond just m = 0 for F e , two approaches may be taken. The first is to simply add the coefficients to the inputs and outputs of F e and assume the neural network will learn a roughly equivariant mapping through the random sampling of the roll rotation. As we demonstrate in Appendix A, the learned functions are indeed roughly equivariant. The second strategy takes a more direct approach to encouraging the network to learn equivariant mappings by taking advantage of the proprieties of the coefficients as they're rotated by φ about the z-axis. If D φ rotates a set x of coefficients about the z-axis by φ, the m = 0 coefficients are constant as a function of φ, while the m ∈ {−1, 1} coefficients are sine and cosine functions of φ, see Appendix C: for some set of vectors α, β and γ of size L. x (0) are the m = 0 coefficients and similarly for x (−1) and x (1) for the m = −1 and m = 1 coefficients respectively. Similar proprieties hold for m > 1 and m < −1. We encourage the output of the message block towards equivariance by computing F e at multiple rotations φ about the z-axis and combining them based on m: For m ∈ [−1, 1], we compute four samples or "taps" at φ ∈ {0, 1 2 π, π, 3 2 π}, and combine them using: wherem k+1 st(0) and F φ e(0) are the output m = 0 coefficients, and similarly form k+1 st (1) and F If higher m are desired, similar calculations can be performed. However, a larger number of taps will be needed, e.g., m ∈ [−2, 2] requires 8 taps separated by 1 4 π radians. The number of coefficients used during message passing when calculating F e may be varied based on the edge properties to reduce the memory required by the network. For instance, atoms that are far away from each other may not need the same resolution of the spherical functions as those close together, so a lower maximum degree L may be used. By training several message passing networks with varying degrees for different ranges of edge distances, memory usage can be reduced while not impacting accuracy. Similarly, the range of m can typically be reduced to [−1, 1] or [−2, 2] since m for values greater than 2 (or less than -2) don't vary as much in the direction of the z-axis, reducing their observed utility to the network. See Figure 1(b) for examples of the spherical channels when the range of m is reduced. Note that while the number of coefficients used by F e may be reduced, the node embeddings maintain their original size based on the maximum L, so they may be able to aggregate information from all neighboring atoms at different distances and angles.

Message Aggregation
At this point in the network, the per edge messages are computed and rotated back to the global coordinate frame. Ideally, to allow for complex interactions between the messages, a non-linear function would be applied after summing all messages directed at a target node. If we did this by naively passing the summed coefficients m through a fully connected neural network, the network may have difficulties learning representations that are approximately rotation equivariant (a fundamental property of atom forces). Another option is to place constraints on the non-linearities performed to enforce equivariance using Clebsch-Gordan Nets [25] or steerable MLPs [5].
We propose applying an unconstrained non-linear pointwise function on the sphere, i.e., a function that is applied at every orientation without knowledge of the orientation from which the point was sampled [9]. For our function, we use a non-linear neural network F c (R C → R C ) that combines information across channels. In practice, F c is applied at a discrete number of orientations, which results in a transformation that we numerically demonstrate (see Appendix A) is approximately equivariant to rotations: are converted from a spherical harmonic representation to a spherical grid representation using G. The channels are concatenated and passed through a 3 layer 1 × 1 CNN. Each layer is followed by a SiLU activation. Finally, the channels are converted back to a spherical harmonic representation using G −1 and added to x Figure 4: Illustration of 4 spherical channels as they progress through the network's layers. The model has 12 layers and L = 6. Note how the spherical channels typically change slowly during each layer, but in some cases change more dramatically. Spherical channels are shown projected onto a 2D grid (polar and azimuthal).
where G is a function that converts a spherical function represented by spherical harmonic coefficients, to one represented by point samples on a sphere. For G we use a spherical coordinate system (polar and azimuthal) to generate a 2D discrete representation of the spherical functions, see Figure 3. A 1 × 1 convolutional neural network F c is applied to each discrete sample containing C channels, and the result is then converted back to the spherical harmonic representation using G −1 . This is analogous to having a 2D image represented in the frequency domain, converting it to the spatial domain using an inverse discrete Fourier transform, applying a transformation and converting back to the frequency domain. Since the same operation is applied to all orientations on the sphere, a transformation that is approximately equivariant (up to discrete sampling limitations) to rotations may be learned, see Appendix A for a more detailed discussion.
A diagram of the network is shown in Figure 3. To provide additional information to the network, x (k) t is also converted to a spherical grid representation, concatenated with G(m (k+1) t ) and provided as input to the neural network to compute the final updates to the coefficients x (k+1) t at iteration k + 1. The network can learn circular functions on the sphere if varying resolutions of x (k) are provided to the network, similar to Difference of Gaussian filters for 2D images [29]. The resolution of x (k) can be reduced by lowering the degree L before transforming by G. In our experiments, if multiple resolutions or bands are used, we use resolutions of degree L and L − 1. Three layers of 1 × 1 convolutions on 2C channels (or 4C if two bands are used) are performed with each followed by a SiLU non-linear activation function. To avoid aliasing, the spherical grid is sampled at a resolution of 2 * (L + 1). An illustration of the spherical channels as they progress through K = 12 layers is shown in Figure 4.

Energy and Forces Estimation
We compute the energy E by first estimating a per-atom energy using a pointwise function on the sphere and taking its integral over all possible orientations. The per-atom energy estimates are summed to obtain the overall system's energy: where F energy is a three layer fully connected neural network (R C → R 1 ) with SiLU activation functions. We calculate forces similarly by estimating a force magnitude |f | = F f orce (s (K) i (r)) in every directionr over the sphere. Integration is performed after multiplying the magnitude by the orientationr to obtain directional vectors: where F f orce is a three layer fully connected neural network (R C → R 1 ). In practice, a discrete approximation of the integral is performed for equations (11) and (12) using a set of 128 evenly distributed points on the sphere (see Appendix B.2). Note that the application of a pointwise function of an arbitrary neural network at a finite number of discrete orientations is not invariant (energies) or equivariant (forces) to rotations. However in practice, numerically this operation does closely approximate these properties (see Appendix A.2).

Experiments
We present results on the Open Catalyst 2020 (OC20) dataset [6]    important applications [47]. This dataset is a popular benchmark for the ML community. We begin by comparing results across all tasks on the test set. Since several papers only report results on the IS2RE task using approaches that directly predict relaxed energies, we also train a model specifically for this task and approach, and compare. Finally, we show numerous ablation studies comparing model variations on the smaller OC20 2M dataset.

Implementation Details
The implementation of spherical harmonics and their transformations uses the code provided by Euclidean neural networks (e3nn) [13].

OC20 Tasks
We compare against numerous models trained on the OC20 2M, All and MD datasets across the Structure to Energy and Forces (S2EF), Initial Structure to Relaxed Structure (IS2RS) and Initial Structure to Relaxed Energy (IS2RE) tasks [6]. OC20 2M contains 2 million, OC20 All has 133 million and OC20 MD has 38 million training examples. The 12 and 16 layer SCN models significantly outperform state-of-the-art GemNet-OC [12] on S2EF force prediction (≈ 14% improvement) and IS2RE energy predictions (≈ 10% improvement) when trained on OC20 2M, Table 1. Comparable results are found for energy MAE on S2EF with GemNet-OC. We also report throughput efficiency. Note that while the SCN models process fewer samples per second, they are much more sample efficient than GemNet as shown in Figure 2(b). Improved efficiency has also been shown for fully equivariant models [3].
We compare our largest SCN models trained on the All and the All+MD datasets (133M + 38M examples) across all OC20 tasks in Table 2  MAE for S2EF. On the challenging relaxation based IS2RE task that requires models to predict both accurate forces and energies, SCN outperforms GemNet-OC by over 7%. Note that we only train a single model to predict both energy and forces, where GemNet-OC-L-F+E + MD is the combination of two models, one trained specifically for forces and the other for energy.

OC20 IS2RE Direct
The task of Initial Structure to Relaxed Energy (IS2RE) may be accomplished using two approaches: 1) Use an S2EF model to relax the positions of the atoms (find local energy minima) and output the energy at the relaxed structure, or 2) directly predict the relaxed energy from the initial structure without performing the relaxation. Empirically, the first approach reported previously in the paper achieves higher accuracy (see Tables 1 and 2), while the second approach is more efficient both during training and inference. Although solving the same problem, the second approach of direct prediction is a fundamentally different problem than the first. The direct approach cannot take advantage of detailed position and angular information, since the initial structure only contains a rough placement of the atoms. Instead it must rely more on global knowledge of how the atoms are arranged. Thus, an approach that works well for relaxation based approaches to IS2RE might not work well for direct approaches.
As shown in Table 3, the SCN model achieves state-of-the-art results for energy MAE across all test splits for models that do not use ancillary loses [45,15]. Energy within Threshold (EwT) is comparable to SEGNN [5]. Improvements for SCN are most pronounced for the out-of-domain splits, demonstrating its increased ability to generalize to unseen data.

Ablation studies
We explore numerous model variations in Table 1 trained on the OC20 2M dataset and evaluated on the OC20 validation set. The first set of ablation experiments explores reducing the complexity of our model. If messages are not rotated before and after applying F m in Equation 2, i.e., A st is an identity matrix, the results are significantly worse for energies and even more so for forces since forces are highly dependent on angular information. When scaling the network, increasing L has a more significant impact on force estimation, while depth improves energy estimation. This may be due to higher L resulting in greater angular resolution, which is important to force estimation. Greater depth, which allows information to be shared more often between atoms, may lead to better energy predictions since energy is function of the entire atomic structure.
Sampling multiple z-axis rotations when computing F e with 4-taps (Equations 7, 8 and 9) produces significantly improved force predictions, while using two bands (L and L − 1) when aggregating messages improves both energy and force prediction. The use of 2 bands reduces model throughput a small amount, while sampling with 4-taps reduces throughput by over 2×. However, sample efficiency is further increased with 4-taps, resulting in similar training times for the same accuracy, Figure 2(b).

Related work
Molecular and atomic property prediction has made significant recent progress. Models relying on hand-crafted representations like Behler-Parrinello [4] and sGDML [8] have recently been surpassed by learned feature representations using GNNs [3,22,24,35]. Early GNN developments focused only on invariant representations. CGCNN [44] and SchNet [34] make energy and force estimations using only atomic species and pair-wise distance information. DimeNet [22,23], SphereNet [27], and GemNet [24,12] extend this to explicitly capture triplet and quadruplet angles, which are scalars. Utilizing invariant representations offers flexibility in model architecture design since the node or edge features are inherently invariant and model constraints are not needed to maintain equivariance. The enumeration of triplets and quadruplets of atoms to calculate functions of relative or dihedral angles requires careful model design to maintain efficiency. More recently, equivariant models, capturing both scalar and equivariant features, have outperformed traditional GNNs on small molecular datasets including MD17 [7] and QM9 [31]. While not strictly equivariant, our work falls into this later class of GNNs for atomic property predictions.
Models that strictly enforce equivariance to the SO(3) (3D rotations) or E(3) group (3D rigid transformations) share similarities to our use of spherical channels [35,1,3,9,5]. Both represent atoms by higher-order tensors and not just scalar information. For instance, Tensor Field Networks [41], Cormorant [1], NequIP [3], SEGNN [5] and BOTNet [2] use the same spherical harmonic representations as SCNs, while PaiNN [35] is an l = 1 variation. All models place constraints on the operations that can be performed to ensure equivarance, such as restricting non-linear operators to invariant scalar inputs, e.g., pair-wise atom distances. In networks like TFN, geometric information of higher degree l are mixed into the l = 0 features via tensor products on which a gated non-linear function is applied. Similar to our approach, the recently proposed SEGNN uses the spherical harmonic representation to enable steerable features based on edge orientations. SEGNN introduces steerable MLPs that enforce equivariance by restricting the learnable parameters to those that scale the Clebsch-Gordan tensor products followed by gated non-linearities [42]. In our approach, we increase the expressivity of the model by identifying a set of m = 0 coefficients across degrees l that are invariant to rotations given an edge's orientation to which unconstrained functions may be applied. We also demonstrate that improved performance can be achieved by relaxing the requirement for strict message equivariance by including a broader set of coefficients beyond m = 0. Furthermore, we project the spherical channels onto a grid and perform a pointwise 1 × 1 convolution followed by a non-linearity. Although not strictly equivariant, the pointwise operation allows for complex mixing of different degrees of spherical harmonics resulting in rich geometric descriptions.

Discussion
A limitation of SCN is that it directly predicts forces, so energy conservation is not enforced. This may limit its use in chemistry applications such as molecular dynamics. Non-energy conserving models are still practically useful for applications such as structure relaxations and transition state searches. The SCN model scales O(n 2 ) with respect to L, which may limit using L larger than 8. Using a larger set of m coefficients during message passing also faces challenges, since the network may find it more difficult to learn approximately equivariant mappings. In Table 1, we saw that increased depth improves results, especially to energy prediction. It remains an open question whether deeper networks will see further improvement.
While this work is motivated by problems we face in addressing climate change, advances in chemistry may have numerous use cases. Some of these use cases may have unintended consequences, such as the development of the Haber-Bosch process for ammonia fertilizer production that enabled the world to feed a growing population, but over fertilization has recently led to ocean "dead zones" and its production is very carbon intensive. More alarmingly, the knowledge gained from fertilizer production was used to create explosives during wartime [17]. We hope to steer research in this area towards beneficial applications by utilizing datasets such as OC20.
In conclusion, we propose a GNN that uses spherical harmonics to represent channel embeddings that contain explicit orientation information. While the network is encouraged to learn output mappings that are equivariant to rotations, best performance is achieved by relaxing the equivariant constraints and allowing for more expressive non-linear transformations. We demonstrate stateof-the-art results across numerous tasks on the large-scale OC20 dataset for modeling catalyst for applications addressing climate change. Appendix:

A Empirical measurement of equivariance to rotations
The SCN is not strictly equivariant to rotations, but depending on the design choices approximate equivarance may be achieved. We begin by empirically measuring the network's invariance and equivariance to rotation for energy and forces respectively. We accomplish this by computing the absolute difference between pairs of results computed from an input and a randomly rotated version of the input: where d is the 3D positional difference between pairs of atoms, a are their atomic numbers, and R is a random 3D rotation matrix. SCN energy and SCN f orce are the network's predictions for energy and forces respectively.
Mean Absolute Difference (MAD) results for various model choices are shown in Table 4 for models with 12 layers and L = 6. Differences are averaged over a model's outputs for a random 1,000 atomic structures. There are four sources that may lead the network to predict different values for rotated versions of the input: 1) the use of m = 0 coefficients during message passing, 2) non-linear message aggregation, Equation (3), 3) the energy and force output blocks, Equations (4,5), and 4) limits to numerical precision especially when using Automatic Mixed Precision (AMP).
In Table 4, we observe that the MAD increases when m ∈ [−1, 1] coefficients are used during message passing as compared to the m = 0 message passing that is equivariant. However, the MAD of the m ∈ [−1, 1] model that uses 4 taps during message passing is nearly identical to m = 0 when not using AMP. This shows that the 4-tap model is close to producing equivariant results. Note MAD is non-zero even for m = 0 due to the other sources of errors.

A.1 Message aggregation
Next, we discuss the degree to which our message aggregation Equation (3) is empirically equivariant. If Equation (3) was equivariant, the following would hold: where D is a block diagonal matrix containing Wigner-D matrices that perform a rotation of the coefficients. That is, the same result should be achieved whether the input or output is rotated. Since the transformation G uses a discrete sampling of the sphere, Equation (3) is clearly not strictly equivariant. The neural network F c may introduce non-linearities that result in frequencies higher than those able to be represented by the maximum degree L used by the spherical harmonics. As a result, the sampling used by G may be below the Nyquist rate needed to avoid aliasing. This may be minimized but not eliminated by utilizing smooth activation functions such as SiLU. However, in practice we find that the use of Equation (3) with sampling resolution of 2(L + 1) or higher does in practice lead to results that are approximately equivariant. In Table 5, we show several empirical results with different depths and activation functions. If only a linear layer is used, the function is equivariant and there is no difference in computed values. For non-linear activations, such as ReLU and SiLU, differences of 4% − 7% are seen between the rotated versions. Percentages are the Mean Absolute Difference (MAD) divided by the mean absolute value of the outputs. MAD is computed using: Minimal differences are seen between using SiLU or ReLU. However, higher percentage differences are observed when using two layers of non-linearities when compared to one. We use two-layers followed by a linear layer in the models for this paper.

A.2 Energy and Force Output Blocks
Our final set of experiments measure the empirical invariance and equivariance to rotations for the output energy and force blocks. The differences between the not rotated and rotated versions is computed using: i rr where R is a random 3D rotation matrix and D its corresponding block diagonal matrix containing the Wigner-D matrices. Results are shown in Table 6. For both energy and forces the results are nearly equivariant, and demonstrate the outputs blocks have negligible negative impact on the overall network's equivariance. This is likely due to the force and energy output blocks integrating over the entire sphere, which may result in any differences negating each other.  Table 6: Mean Absolute Differences (MAD) for 3-layer energy and force output blocks when the inputs are rotated, Equations (17,18). The MADs are small (< 1%) when compared to the energy and force MAEs, and significantly smaller than other blocks.

B Implementation details
In this section, we describe additional implementation details not in the main paper.

B.1 Varying spherical channel resolution
For atoms that are far away from each other, the messages may be computed using lower resolution spherical channels to reduce memory usage while maintaining similar accuracy. For all results in this paper, the 12 closest neighbors to each atom use the settings described in the paper. For neighbors ranked 13 to 40 by distance, L is reduced by 2 with a minimum value of 4, and H is reduced by a factor of 4. Messages are not calculated for neighboring atoms with ranked distances greater than 40 or a cutoff distance greater than 8 Å.

B.2 Point sampling on the sphere for output blocks
For both the energy and force output blocks the integrals in Equations 4 and 5 in the main paper are approximated using points samples on the sphere. Finding a set of evenly distributed points on a sphere for an arbitrary number of points remains an open problem. We use points that are approximately evenly distributed using spherical Fibonacci point sets [18]. Since the density of points does vary slightly across the sphere, we weight each point using a Gaussian weighted counting of nearby points on the unit sphere with σ = 0.5.

B.3 PaiNN baseline
The PaiNN baseline is our reimplementation of Schütt et al. [35] with the difference that forces are predicted directly from vectorial features via a gated equivariant block instead of gradients of the energy output. This breaks energy conservation but is essential for good performance on OC20.

C Note on spherical harmonics properties
In the paper it is stated that m = 0 spherical harmonics are invariant to rotations about the z-axis. This is easily seen by looking at the equations for the spherical harmonics. We show the Y m l equations for up to l = 2 below, parameterized by θ (polar) and φ (longitudinal rotation about the z-axis). Notice that all m = 0 spherical harmonics are only a function of θ, and thus invariant to changes in φ.   The validation error is calculated on a 30k subset of the validation ID dataset. Note significant overfitting is seen for energies but not forces. Errors are smoothed using a 0.9 exponential moving average. Plot is generated from training run in Table 2 Figure 5 shows the train and validation errors during training. Interestingly, we see very little overfitting on forces, but energy has significant overfitting. This may be due to the energies being a per-structure property, while the forces have more examples since they are a per-atom property. However, overfitting on the energies is still surprising given the large number of examples in the OC20 All + MD training dataset. This indicates that the SCN model can fit complex functions but could improve on its ability to better generalize through model improvements or data augmentation. Similar trends are found when trained on the OC20 2M dataset.