[3]Liyang Zhou [1,2]Jian Jiang

1]Beijing National Laboratory for Molecular Sciences, State Key Laboratory of
Polymer Physics and Chemistry, Institute of Chemistry, Chinese Academy of
Sciences, Beijing, 100190, P. R. China

2]University of Chinese Academy of Sciences, Beijing, 100049, P. R. China

3]Juhua Group Co., Ltd, Quzhou, 324004, P. R. China

# Efficient Machine Learning Force Field for Large-Scale Molecular Simulations
of Organic Systems

Junbao Hu<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [

###### Abstract

To address the computational challenges of ab initio molecular dynamics and
the accuracy limitations of empirical force fields, the introduction of
machine learning force fields has proven effective in various systems
including metals and inorganic materials. However, in large-scale organic
systems, the application of machine learning force fields is often hindered by
impediments such as the complexity of long-range intermolecular interactions
and molecular conformations, as well as the instability in long-time molecular
simulations.Therefore, we propose a universal multiscale higher-order
equivariant model combined with active learning techniques, efficiently
capturing the complex long-range intermolecular interactions and molecular
conformations. Compared to existing equivariant models, our model achieves the
highest predictive accuracy, and magnitude-level improvements in computational
speed and memory efficiency. In addition, a bond length stretching method is
designed to improve the stability of long-time molecular simulations.
Utilizing only 901 samples from a dataset with 120 atoms, our model
successfully extends high precision to systems with hundreds of thousands of
atoms. These achievements guarantee high predictive accuracy, fast simulation
speed, minimal memory consumption, and robust simulation stability, satisfying
the requirements for high-precision and long-time molecular simulations in
large-scale organic systems.

## Main

Molecular dynamics (MD) simulation has gained significant attention in recent
years across various disciplines, spanning physics, chemistry, biology, and
materials science. This cutting-edge technology, by simulating interactions
between molecules or atoms, offers researchers a means to investigate the
microstructure of substances and understand their macroscopic properties. This
technique has proved invaluable for experimental design, development of new
materials, and advances in biomedical research [1, 2, 3, 4].

Traditional molecular simulations grapple with the dilemma of balancing the
high computational cost of ab initio molecular dynamics (AIMD) against the low
precision of empirical force fields. A resolution to this challenge is found
in the application of machine learning, leveraging its powerful fitting
capabilities. The fundamental idea of machine learning force field (MLFF) is
to establish a mapping from molecular coordinates to the labels of high-
precision quantum chemistry data, including potential energy and forces. This
approach eliminates the need to solve the intricate Schrödinger equation,
resulting in a significant acceleration and achieving a balance between
prediction precision and simulation speed [5].

Since the advent of BPNN in 2007[6], numerous MLFF models have been proposed
to improve prediction accuracy, and their performance has been systematically
investigated in public datasets (MD17[7], MD22[8], OC22[9]). From the
perspective of tensor order (denoted by $l$), existing 3D molecular
representation learning models can be categorized into two main classes: one
is the invariant graph neural networks with only scalar features (i.e., $l=0$)
, including SchNet[10], DeePMD[11], DTNN[12], PhysNet[13], ComENet[14],
SphereNet[15]; the other is the equivariant graph neural networks with vector
features (i.e., $l=1$) including EGNN[16], PaiNN[17], GVP-GNN[18], EQGAT[19],
as well as higher-order features (i.e., $l>1$) such as TFN[20], Cormorant[21],
SEGNN[22], Equiformer[23], NequlP[24], Allegro[25], BotNet[26], MACE[27].
Specifically, invariant methods directly use invariant geometric features such
as distance and angles as input, ensuring invariance to the rotation and
translation transformations on input molecules. In contrast, the equivariant
model can maintain certain symmetries or properties under specific
transformations such as rotation and translation.[28]. Models with
equivariance properties for $l\geq 1$ generally outperform invariant ones on
various public datasets and tests[16, 24, 29, 30]. Furthermore, in terms of
the expressive power of geometric graph neural networks, higher-order
equivariant models demonstrate more expressive capabilities compared to the
first-order equivariant models[31]. Increasing the equivariant order often
leads to improved model accuracy[32, 27, 33]. Therefore, higher-order
equivariant force field models such as NequIP, Allegro, and MACE achieve
excellent performance on multiple MD simulation metrics [30] with nearly a
1000-fold improvement in data efficiency [24, 25].

The high accuracy MLFFs (the errors in atom energy and forces are around 1
meV/atom and 50 meV/Å, respectively [34, 35]) have been successfully applied
in diverse systems such as metallic and non-metallic inorganic materials[36].
However, there is limited research to date on the application of MLFF in
organic systems. In fact, due to the weak generalization ability of machine
learning models[37], the application of MLFF in large-scale organic systems is
usually challenged by several impediments such as the complexity of long-range
intermolecular interactions and molecular conformations, as well as the
instability in long-time molecular simulations. Therefore, molecular
simulations based on MLFF often fail (e.g., loss of accuracy, bond breaking,
and atomic overlap) when encountering scenarios not present in the training
set [30, 38] .

Expanding the receptive field of neural network models by appropriately
increasing the number of interaction layers has been proven to enhance the
predictive accuracy of MLFF in multi-molecular interaction systems[39].
However, increasing the number of interaction layers results in higher
computational costs and can lead to over-smoothing issues that diminish the
expressive capability of the machine learning model[40]. In fact, expanding
the model’s receptive field by increasing the number of interaction layers
does not guarantee accurate characterization of long-range intermolecular
interactions[29].

The key to accurately capturing long-range intermolecular interactions is
increasing the cutoff radius of the neural network model’s hyperparameter.
However, the number of neighboring atoms increases cubically with increasing
cutoff radius, resulting in a significant increase in simulation time and
memory consumption. This limitation largely hampers the application of MLFF in
long-time MD simulations for large-scale organic systems. Consequently, many
models adopt modular methods to consider long-range intermolecular
interactions. For instance, a fragment-based method[39] was used to
demonstrate the long-range intermolecular interactions. However, this
fragment-based approch lacks generality and may require specific fragmentation
schemes for certain molecules. The DeePMD-LR method[41] directly incorporates
existing empirical equations to account for the long-range intermolecular
interactions, but the limited applicability and the requirement for expert
knowledge in empirical force fields restrict the use of this method. The
fourth generation HDNNP network[42] employs a separate neural network module
to predict long-range interactions, which results in an overly complex model
with significant computational costs, limiting its application in long-time MD
simulations. The MFN model[43] considers the long-range intermolecular
interactions by adopting a new architecture that incorporates analytic matrix-
valued functions. However, the computational complexity and the memory
consumption scale as $O(N^{2})$ and $O(N^{4/3})$, respectively, where N is the
number of atoms. Therefore, there is a pressing need for an efficient end-to-
end method to capture the long-range intermolecular interactions in the
organic systems.

In addition, the extensive conformational variability of organic compounds
requires that the training set ensures not only quantity but also
diversity[44]. Unfortunately, AIMD proves to be prohibitively expensive to
generate a large number of diverse conformations. Consequently, an efficient
method for collecting training sets becomes imperative, emphasizing the need
for non-redundant datasets to mitigate the labeling costs associated with
quantum-mechanical calculations such as density functional theory (DFT).
Moreover, the impracticality of preparing the training set from high-precision
quantum DFT calculations for large organic systems further complicates the
labeling process. As a result, the labeled datasets are typically derived from
smaller systems[45]. This limitation underscores the necessity of
investigating whether MLFF methods can maintain high-precision generalization
when applied to larger systems.

Finally, to address the issue of instability in long-time molecular
simulations, various methods have been proposed. For example, directly
incorporating configurations not present in the training set can enhance the
stability of the simulation [38], but it is impractical to exhaustively
consider all possible configurations. Reducing the time step in MD is another
strategy to decrease the incidence of simulation crashes [24, 25]. However,
this approach comes at the cost of a proportional increase in simulation time.
Although larger models with millions or more parameters may alleviate the
issues mentioned above, they do not guarantee simulation stability [30].
Therefore, a more efficient and straightforward MLFF model is required to
ensure stability during long-time MD simulations.

To address the aforementioned challenges, we propose an end-to-end, highly
efficient, and broadly applicable higher-order equivariant model to
effectively handle the long-range intermolecular interactions. This method
exhibits excellent competitiveness in various aspects, including prediction
accuracy, simulation speed, memory consumption, and MD results. Moreover, the
incorporation of a bond length stretching technique significantly boosts the
simulation stability of MLFF models in organic systems. In addition, an active
learning approach based on committee queries is adopted to efficiently collect
datasets. Notably, this study successfully extends high precision from a
limited dataset comprising 901 samples with 120 atoms to larger systems
encompassing ten thousand atoms. This achievement plays a pivotal role in
facilitating long-time MD simulations in large-scale organic systems.

## Results

### Efficient and universal multiscale higher-order equivariant modeling
architecture

Figure 1: |The framework of the universal multiscale higher-order equivariant
model. The model comprises an embedding layer, an interaction layer, and an
output layer. The interaction layer may consist of multiple layers that are
responsible for considering interactions between the central atom and its
neighboring atoms. The yellow module is designed to efficiently capture the
long-range intermolecular interactions. This module achieves efficient
processing of long-range message by reducing the channel numbers of nodes and
lowering expansion order of direction information. After the aggregation of
the long-range message, the feature dimensions of the long-range nodes and the
local nodes are aligned, and then they are added together. Here, $\sum$
denotes summation, $\parallel$ represents feature concatenation, TP denotes
tensor product, MLP stands for multilayer perceptron, and Attention indicates
the attention mechanism.

Using the Quantum Mechanics (QM) method to compute large-scale systems remains
a challenging task[46, 47]. However, certain studies focus only on specific
regions within the system, such as the active sites of enzymes. Therefore, the
Quantum Mechanics/Molecular Mechanics (QM/MM) hybrid method, also known as a
multiscale method, is widely employed to address this challenge. The key
concept of this method is that for the target regions, a high computational
cost QM method is used, and for nontarget regions, a low computational cost
Molecular Mechanics (MM) method is employed; thereby achieving efficient
computations for large-scale systems[47].

Most of the MLFF models are based on the locality hypothesis of atoms[48, 49].
This assumption claims that atoms only interact with their neighbors and is
physically justified by the nearsightedness principle of electronic matter.
[50] Furthermore, this assumption infers that interactions between atoms that
are far apart are negligible. In order to take into account long-range
intermolecular interactions, the MLFF models must increase receptive field. In
this work, inspired by the QM/MM method, we propose a universal multiscale
higher-order equivariant model, in which a high-cost module is used to
demonstrate the short-range interactions related to the central atom, while a
low-cost module is employed to capture the corresponding long-range
interactions. This method adopts a multiscale strategy that empowers MLFF
models to attain optimal precision and efficiency in handling intermolecular
long-range interactions. Specifically, the potential energy of each atom is
divided into short-range and long range parts, and the total potential energy
of the system is obtained by summing the contributions from all atoms. The
atomic forces are the negative gradients of the total potential energy with
respect to the coordinates, ensuring the energy conservation of the model[51]:

$E_{pot}=\sum_{i\in
N_{\text{atoms}}}(E_{i}^{\text{short}}+E_{i}^{\text{long}})$
$\overrightarrow{F_{i}}=-\triangledown_{i}E_{pot}$

As mentioned above, in terms of fitting atomic potential energy, higher-order
equivariant models significantly outperform the first-order equivariant models
and invariant models. Therefore, in this work, we only focus on higher-order
equivariant models. Various higher-order equivariant models have been proposed
through the construction of different nonlinear functions, convolution
filters, message functions, aggregation functions, and update functions.
Examples of these models include TFN [20], Equiformer [23], Cormorant [21],
MACE [27], NequIP [52], SEGNNs [53] and so on. Our universal and efficient
multiscale higher-order equivariant model is proposed based on these
groundworks (Fig. 1). This multiscale model achieves fast processing of long-
range messages by reducing the channel dimensions of nodes through equivariant
linear transformations and lowering the expansion order of the embedding
direction information. After the long-range messages are aggregated, another
linear transformation is used to align the long-range and short-range
features. Finally, the summation of these two features gives the entire
multiscale equivariant model. We have to note that this multiscale model is
general. Any higher-order equivariant model can be combined with this model.
In this work, our multiscale higher-order equivariant model is constructed
based on the MACE model. Hereafter, it will be referred to as “Our” or “Our
(MS-MACE)”.

### Intermolecular interaction

Figure 2: |The potential energy predictions from three different machine
learning models. a The schematic diagram of two formaldehyde molecules. The
shortest and farthest distances between the right molecule and the carbon atom
of the left molecule are 8 Å and 9 Å, respectively. b the potential energy
predictions on the test dataset. The blue dash line denotes the model with
three interaction layers (n=3) and a cutoff radius of 8 Å, the green dash line
represents the model with one interaction layer (n=1) and a cutoff radius of 9
Å, the red dash line represents our multiscale model, and the black line
represents the GFN2-xTB results.

Increasing the number of interaction layers generally enhances the receptive
field of the model [54]. However, for long-range intermolecular interactions,
there may be scenarios that lack intermediary atoms, which causes an
interruption in propagation, consequently leading to the ineffectiveness of
this model. On the contrary, increasing the cutoff radius of the neural
network model’s hyperparameter proves to be the key to accurately capturing
long-range intermolecular interactions. For example, we consider the total
potential energy of two formaldehyde molecules that are arranged in a face-to-
face parallel alignment in space. Taking the carbon atom of the left molecule
as the center atom, the shortest and longest distances to the right molecule
are 8Å and 9 Å, respectively, as shown in Fig. 2a. To prepare the dataset, we
rotate one of the molecules around the central axis of two carbon atoms by 1
degree each time, and collect 360 samples. In this dataset, we select one
sample every 3 degrees to include in the training dataset, while the remaining
configurations serve as the test dataset. The energy and forces for each
sample are calculated using the semi-empirical GFN2-xTB level of theory[55].
Three different machine learning models based on the MACE model are
considered: (i) three interaction layers (n = 3) with a cutoff radius of 8 Å;
(ii) one interaction layer (n = 1) with a cutoff radius of 9 Å; (iii) our
multiscale model with one interaction layer (n = 1) and cutoff radii of 3 Å
for short-range interactions and 9 Å for long-range interactions. To guarantee
sufficient fitting capacity for all models, in the first two models and the
short-range part of the third model, the node features are set to 256x0e, the
direction information is embedded using a third-order expansion, the cutoff
polynomial envelope is set to 60, and the batch size is set to 2. In the long-
range part of the third model (our MS-MACE model), the node features are set
to 32x0e with only first order expansion for direction information embeddings.
All other hyperparameters remain the same.

Due to the lack of intermediary atoms between these two formaldehyde
molecules, the transmission of the message is interrupted. Consequently,
merely increasing the number of interaction layers without enlarging the
cutoff radius is insufficient to capture long-range intermolecular
interactions. Therefore, according to Fig. 2(b), the first model (n = 3,
$r_{c}$ = 8 Å) can only predict the average energy (the blue dash line). On
the other hand, although the number of interaction layer is reduced to one (n
= 1), the variations of the potential energy with respect to angles can be
effectively captured by increasing the cutoff radius to 9 Å (the green dash
line). However, a longer cutoff radius incurs significant computational cost
and memory consumption, especially for large-scale systems. In our multiscale
model, inspired by the QM/MM method, we consider the short-range interactions
using a very short cutoff radius ($r_{c}$ = 3 Å) but a larger channel
dimensions (256x0e) and a higher-order expansion for the direction
information, while the long-range intermolecular interactions are captured by
a large cutoff radius ($r_{c}$ = 9 Å) but a smaller channel dimensions (32x0e)
and only first-order expansion for the embedded direction information.
According to Fig. 2(b), our model gives the most accurate predictions (red
dash line). We have to note that our model not only exhibits the highest
accuracy, but also achieves magnitude-level improvements in computational
speed and memory efficiency. That is to say, to consider long-range
intermolecular interactions, the MLFF models should increase the cutoff radius
rather than the number of interaction layers, the channel dimensions of node
features, and the expansion order for the embedded direction information. The
accuracy and efficiency of our multiscale model will be discussed in detail
below and in Fig. S1 of the Supporting Materials.

### Stability of long-time simulation

Figure 3: |Long-time molecular simulation stability. a, The violin plots of
the C-F bond length distributions for perfluorotri-n-butylamine in the
training datasets are collected using GFN2-xTB level of theory at different
temperatures. In addition, the bond length distribution of the training
dataset for 300 K, improved by the bond length stretching method, is also
presented (the red shape). b and c show the average absolute errors of energy
and forces for different machine learning models on the test dataset of
various C-F bond lengths. The plus symbol denotes the MS-NequIP model, while
the dot symbol denotes the MS-MACE model. The blue color represents the
original training dataset derived from the 300 K trajectories (hereafter, this
training dataset is referred to as “300K”), while the red color represents the
training dataset improved by the bond length stretching method (hereafter,
this training dataset is referred to as “300K-Stretch”). The black line
indicates the GFN2-xTB result and the blue dash line is the zero line. d, The
violin plots for the stability of MD simulations. In e and f, the curves
describe the C-F bond lengths change in MD iterations at 600 K with initial
bond lengths of 1.9 Å and 2.1 Å, respectively. The blue and green lines
represent the MS-NequIP model trained from the “300K” dataset and the
“300K-Stretch” dataset, respectively. Similarly, the yellow and red lines
represent the MS-MACE model trained from the “300K” and “300K-Stretch”
datasets, respectively.

The existing MLFF models are easy to form excessively long or short bonds
during long-time MD simulations in organic systems, leading to the collapse of
simulations[38, 56, 30]. This is due to the fact that most of the samples in
the training dataset are obtained from sampling near the equilibrium
state[57]], and the proportion of molecular conformations with abnormal bond
lengths in the dataset is very low. Consequently, the weak generalization
ability of machine learning models leads to the collapse of MD simulations.
However, it is challenging to directly collect conformations with extremely
long or short bond lengths from MD trajectories, because the MD process would
collapse instantaneously in this scenario. Therefore, we propose a data
augmentation method to enhance the stability of the model during long-time MD
simulations. We first collect conformations from normal MD trajectories and
then use a method of stretching bond lengths to expand the training dataset.

To demonstrate this bond length stretching method, we take an organic molecule
as example, i.e., perfluorotri-n-butylamine (molecular formula
$\text{C}_{12}\text{F}_{27}\text{N}$, CAS No. 86508-42-1), which is a
nonconducting, well-thermally and chemically stabilized dielectric liquid and
is widely used as a fully immersible electronics coolant. Molecular
simulations of 300 ps are performed at 300 K, 600 K, and 800 K using AIMD
simulations based on semi-empirical GFN2-xTB level of theory under the NVT
ensemble. Fig.3a show the distribution of C-F bond length at different
temperatures. The results indicate that higher temperatures lead to broader
distributions. However, even at 800 K, the distribution of bond lengths
remains too narrow to guarantee the stability of long-time MD simulations.
Inspired by this fact, the bond length stretching method is developed as
follows. First, the training dataset is constructed by randomly selecting 600
samples of conformations from the trajectories at 300 K. Then, in every sample
within this training dataset, we randomly select a chemical bond and multiply
the length of this bond by a factor of $\lambda$, where $\lambda$ is a random
variable in the range of [0.85, 2.0]. This range is chosen to avoid
overstretching and overcompression of the bonds. For comparison, the bond
length distribution of the dataset with stretched bonds is also added to Fig.
3a. The results indicate that the stretching method can significantly improve
the bond length distributions.

To verify the effectiveness of this bond stretching method, we first combine
the two representative equivariant models, i.e., the MACE and NequIP models,
with our multiscale model, resulting in the MS-MACE and MS-NequIP models,
respectively. Then these two multiscale models are trained on datasets that
are both processed and unprocessed by the bond length stretching method.
According to Figs. 3b and 3c, our bond length stretching method significantly
improves the prediction capabilities of potential energy and forces. As shown
in Fig. 3c, without employing the bond length stretching method, these two
MLFF models would even yield unphysical predictions of forces for extremely
short or long bond lengths. Furthermore, 150 conformations are randomly
sampled from the MD trajectories at 800 K to serve as the initial
configurations. Then 50000 steps of Langevin dynamics simulations based on the
above MLFF models are performed at 800 K with a time step of 3 fs. The MD
simulation stability of MLFF models can be analyzed by counting the number of
frames in which collapse occurs (the criterion for collapse is that there
exists a bond length deviating by more than 5 Å from its equilibrium
state)[38, 30]. According to Fig. 3d, without employing the bond length
stretching method, the stability of the MS-MACE model is superior to that of
the MS-NequIP model, with a success rate of 88.67% for the former and only
2.0% for the latter, further emphasizing the excellent simulation stability of
the MACE framework [29]. With employing the bond length stretching method, the
simulation stability of both models is greatly enhanced, with a success rate
of 100%. This suggests that the bond length stretching method exhibits high
generality and effectiveness.

Figure 4: |Active learning workflow and results. a The schematic diagram of
active learning techniques within the framework of the multiscale higher-order
equivariant model combined with the bond length stretching method. b, c, d and
e respectively show the distribution of the maximum force variance
($\sigma_{f}^{\text{max}}$) in the first, second, third, and tenth rounds of
the data collecting process (or exploring process). Different colors represent
exploration at different temperatures, where the blue, green, yellow, and red
lines correspond to 300 K, 500 K, 700 K, and 900 K, respectively.

In addition, according to the inset of Fig. 3c, without employing the bond
length stretching method, MS-NequIP and MS-MACE models give unphysical
preditions (i.e., repulsive forces) when the bond lengths are larger than 1.9
Å and 2.0 Å, respectively. This seems to suggest that the threshold bond
length for unphysical predictions determines the stability of the long-time MD
simulations. Therefore, we infer that the simulation stability of the MLFF
model is highly correlated with the predictive capability of extreme bond
lengths. To further examine the simulation stability of each model, we
investigate the change of C-F bond length over the simulation time in Langevin
dynamics simulations with a time step of 1 fs at 600 K, where the initial bond
lengths are 1.9 Å (Fig. 3e) and 2.1 Å (Fig. 3f), respectively. According to
Fig. 3c, in the case where the bond length is 1.9 Å (larger than the
equilibrium bond length), the MS-MACE predicts a negative force (i.e., elastic
rebound force). Therefore, the elongated chemical bond can be restored to its
equilibrium state (see the yellow line in Fig. 3e). However, in this case, the
MS-NequIP model gives an unphysical prediction (a repulsive force) (Fig. 3c),
resulting an abnormal elongation of the C-F bond during the MD simulations
(see the blue line in Fig. 3e), and ultimately leading to the collapse of
simulations. At a bond length of 2.1 Å, without employing the bond length
stretching method, both models give unphysical predictions (i.e., repulsive
forces), resulting in incorrect elongations of the C-F bonds (see the yellow
and blue lines in Fig. 3f). On the other hand, by the application of our bond
length stretching method, both models can make correct predictions for the
cases where the initial bond lengths are 1.9 Å and 2.1 Å, respectively (see
the green and red curves in Figs. 3e and 3f). In summary, the bond length
stretching method proves to be an effective and universally applicable method
to enhance the stability of long-time MD simulations.

### Efficient collection of training datasets

Due to the limited generalization capability of neural networks, the
performance quality of MLFF model is heavily dependent on the quality of the
training dataset. Moreover, labeling the dataset requires the use of quantum-
mechanical calculations, which involve a high computational cost. Therefore, a
key task prior to training MLFF model is to construct a dataset that is as
comprehensive and non-redundant as possible [34, 58]. Active learning
techniques have emerged as a prevalent strategy for collecting datasets, with
committee querying methods being particularly prominent in the field of MLFF
[59, 45, 49].

Table 1: The percentage of accurate ($\sigma_{f}^{\text{max}}<\sigma_{\text{low}}$), effective candidate ($\sigma_{\text{low}}\leq\sigma_{f}^{\text{max}}\leq\sigma_{\text{high}}$), and failed ($\sigma_{f}^{\text{max}}>\sigma_{\text{high}}$) samples in the $i^{th}$ iteraction. i | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
Accurate | 73.20 | 97.42 | 98.74 | 98.68 | 98.15 | 99.49 | 99.46 | 99.51 | 99.67 | 99.68
Candidate | 26.72 | 2.58 | 1.25 | 1.32 | 1.85 | 0.51 | 0.54 | 0.49 | 0.33 | 0.32
Failed | 0.08 | 0.00 | 0.01 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
|  |  |  |  |  |  |  |  |  |

The fundamental idea of an active learning technique is to use the
disagreement of the committee to quantify the generalization error [60].
Specifically, if a training dataset is of sufficiently high quality, a MLFF
model, operating under various random seed conditions, will yield accurate
predictions with low variance for normal samples that are not present in the
training set. Otherwise, predictions with high variance indicate that the
training dataset is incomplete and requires the addition of more effective
samples for further improvement. Typically, the maximum variance in force
predictions made by the MLFF model under different random seeds is used as a
criterion to determine whether a sample should be added to the training set.
This maximum variance is calculated based on the following formula:

$\sigma_{f}^{\text{max}}=\mathop{\text{max}}\limits_{i}\left\\{\sqrt{\langle\|F_{i}(R_{t})-\langle
F_{i}(R_{t})\rangle\|^{2}\rangle}\right\\}$ (1)

where $i$ is the index of the atom in the candidate sample. If the maximum
variance satisfies the inequality
$\sigma_{\text{low}}\leq\sigma_{f}^{\text{max}}\leq\sigma_{\text{high}}$, then
the sample will be added to the training dataset, where the lower and upper
limits are set to $\sigma_{\text{low}}$ = 100 meV/Å and $\sigma_{\text{high}}$
= 300 meV/Å, respectivley.[45] If the maximum variance exceeds the upper
limit, it suggests that the corresponding sample may be abnormal. Using it as
labeled data could potentially be detrimental to the training effectiveness of
the MLFF model. An extremely low maximum variance ($<\sigma_{\text{low}}$)
indicates that the existing training dataset is sufficiently comprehensive for
the MLFF model to make correct predictions about the corresponding sample.
That is to say that this sample does not need to be added to the training
dataset.

In this section, similar to the previous section, we will take perfluorotri-n-
butylamine molecules as an example to demonstrate the powerful capabilities of
active learning technique for data collection within the framework of our
multiscale higher-order equivariant model. To construct a dataset that is as
comprehensive and non-redundant as possible, the workflow of the active
learning method involves a series of continuous iterations. First, an
annealing simulation is conducted using the general AMBER force field
(GAFF)[61] on a system consisting of 200 perfluorotri-n-butylamine molecules
in the NPT ensemble. During the annealing simulation, the system cools down
from 800 K to 280 K in 60 ns, and 300 frames of trajectory are randomly
extracted with a minimum sampling interval of 10 ps. In each frame, the three
molecules closest to the center of the simulation box are selected as a sample
(120 atoms) and added to the initial dataset. To enhance the simulation
stability during data collection, the initial dataset is improved using the
bond length stretching method proposed in the previous section. Then the
improved dataset is labeled using the semi-empirical GFN2-xTB level of theory.
Second, four MLFF models based on our multiscale higher-order equivariant
model (i.e., MS-MACE) with different random seeds and a batch size of 64 are
trained on the initial dataset (300 samples). Then high precision MD
simulations are conducted in the NVT ensemble at 300 K, 500 K, 700 K, and 900
K, respectively, using the four trained MLFF models with a step length of 1 fs
and a total simulation time of 45 ps. In the meantime, a certain number (10,
20, 40, and 50 for 300 K, 500 K, 700 K, and 900 K, respectively.) of new
samples are randomly selected as the candidate samples from the MD
trajectories. Finally, the candidate samples are assessed using Eq. (1), and
the effective candidates are labeled and added to the training dataset. The
above process is repeated until the proportion of effective candidates
generated from the pool of candidate samples converges to a value less than
0.5%.

Figure 5: |a and b are the changes of the average absolute errors of the
energy and forces per atom with the increase of the atom numbers,
respectively. c and d are the changes of time cost and memory consumption of
simulations with the increase of atom numbers. The straight lines are obtained
by linear fittings to the data points using the least squares method. The
blue, green, yellow, and black lines are from the Allegro, NequIP, MACE, and
our MS-MACE models, respectively.

Fig. 4b-e depict the distribution of $\sigma_{f}^{\text{max}}$ at different
temperatures in the 1st, 2nd, 3rd, and 10th iterations. Table 1 presents the
percentages of accurate ($\sigma_{f}^{\text{max}}<\sigma_{\text{low}}$),
effective candidate
($\sigma_{\text{low}}\leq\sigma_{f}^{\text{max}}\leq\sigma_{\text{high}}$),
and failed ($\sigma_{f}^{\text{max}}>\sigma_{\text{high}}$) samples in each
iteration. In the 1st iteration, 26.80% of samples exceed the lower limit of
the force variance (i.e., $\sigma_{f}^{\text{max}}\geq\sigma_{\text{low}}$),
and the density distribution curves of $\sigma_{f}^{\text{max}}$ at different
temperatures are relatively broad (Fig. 4b). In the first iteration, 26.72% of
the candidates become effective candidates and are added to the initial
training dataset. After just one iteration, the percentage of accurate samples
increased to 97.42%, and the density distribution of $\sigma_{f}^{\text{max}}$
is primarily concentrated in the range between 0 and $\sigma_{\text{low}}$
(Fig. 4c). The proportion of effective candidate samples drops to 2.58%. After
10 steps of iteration, the percentage of accurate samples approaches 99.7%,
and only 0.32% of samples need to be added to the dataset. Referring to
pertinent literature [59, 45], we conclude that the iterations of data
collection using an active learning technique have converged. In the end, we
collect a training dataset consisting of 901 samples.

### Prediction accuracy and molecular simulation efficiency

Figure 6: |Comparisons with AIMD simulations. a is the RDF of F-F atoms. b is
the angle distributions of F-C-F bonds, and c is the dihedral distributions of
F-C-C-F. The blue lines are the results from AIMD simulations and the red
dotted lines are the results from our multiscale higher-order equivariant
model.

In this section, we compare our multiscale model (referred as “MS-MACE” model)
with three other representative equivariant models, including MACE, NequIP,
and Allegro, in terms of prediction accuracy, simulation speed, and GPU memory
consumption. It is important to note that achieving complete consistency in
hyperparameters for all models is impractical due to differences in model
architecture. For example, MACE and NequIP models utilize atoms for message
passing, while the Allegro model employs edges for message passing[28].
However, to ensure the validity of the comparisons, the shared parameters of
the models are set identically. The detailed information on the parameter
settings can be found in the Supporting Materials and the section of Method.
In addition, it merits emphasis that comparisons with invariant models have
not been conducted, because equivariant models typically exhibit superior
performance over invariant models in the construction of MLFF.

In the previous section, we have obtained a dataset comprising a total of 901
samples, each containing three perfluorotri-n-butylamine molecules. This
dataset is randomly divided into training and validation sets, comprising 801
and 100 samples, respectively. Subsequently, we will construct seven test
sets, each comprising 100 samples. Within a single test set, each sample
contains an equal number of molecules with different conformations. Across the
first to the seventh test sets, each sample will contain 20, 40, 60, 80, 100,
120, and 140 molecules, respectively. In order to obtain these seven sets, we
conduct all-atom MD simulations on a system with 1000 molecules in the NPT
ensemble at 600 K using GAFF force field. All of the samples in the seven test
sets are randomly extracted from the MD trajectories, and then labeled using
the semi-empirical GFN2-xTB level of theory. In addition, aiming to reduce the
impact of randomness, the average prediction errors for each model are
obtained from the model with five different random number seeds.

Fig. 5 shows the results of the comparisons among different equivariant
models. All models exhibit high precisions in predicting energy and forces on
the test sets. Moreover, the time cost and GPU memory consumption of
simulations change linearly with the increase of atom number in the system.
This could be attributed to the fact that all models are based on a local
assumption [62, 48, 49].

Compared to other models, our MS-MACE model achieves the best performance
across various metrics, including accuracy, simulation speed, and GPU memory
consumption. Specifically, according to Figs. 5a and 5b, our MS-MACE model
gives the most accurate predictions for both energy and force. In particular,
the average error in force approaches an asymptotic value of 10.5 meV/Å,
indicating that even in a large-scale system comprising several hundred
thousand atoms, the force error remains around 10.5 meV/Å. That is to say,
although each sample in the training dataset only contains 120 atoms, our
multiscale MLFF model can successfully extend high precision to large-scale
systems with hundred thousand atoms. In addition, according to Figs. 5c and
5d, our MS-MACE model achieves magnitude-level improvements in computational
speed and memory efficiency. Furthermore, due to the linear characteristics of
the time cost and memory consumption changes, it is easy to speculate that
when the number of atoms in the system exceeds 100000, the computational speed
and memory efficiency will see orders of magnitude improvement. In summary,
our MS-MACE model demonstrates exceptional competitiveness in energy and force
predicting, simulation speed, and GPU memory consumption.

### Compared with AIMD simulations

In this section, by comparing with the results from AIMD simulations, we
validate the capability of our multiscale higher-order equivariant model. We
conduct AIMD simulations in the NVT ensemble at 600 K for 500 ps. The system
comprises five perfluorotri-n-butylamine molecules, where periodic boundary
conditions are used in all three directions. Under the same conditions, MD
simulations based on MS-MACE model are performed. To reduce the influence of
randomness, all simulations are repeated three times independently. Fig. 6
shows the comparisons of radial distribution functions (RDFs) for F-F atoms,
the angle distributions for F-C-F bonds, and the dihedral distributions for
F-C-C-F bonds. The results indicate that the MS-MACE model can perfectly
reproduce the results from AIMD simulations, highlighting the accuracy of our
multiscale higher-order equivariant model.

## Discussion

Although the polarity of the perfluorotri-n-butylamine molecule is extremely
weak, a cutoff radius of 8 Å is still required to accurately capture the
intermolecular long-range interactions. However, a large cutoff radius is a
disaster for the MLFF model. According to Figs. 5c and 5d, the cost of
simulation time and memory consumption is enormous for a general MLFF model,
which is unaffordable for simulating large-scale systems. Moreover, the cutoff
radius in systems with dipole-dipole or Coulomb interactions will be even
larger. Therefore, a multiscale higher-order equivariant model with
siginificant advantages in prediction accuracy, simulation speed, and memory
consumption is necessary for large-scale organic systems.

In addition, long-time simulation stability is one of the most important and
challenging issues in the development of the MLFF model. To address this
issue, we proposed a bond length stretching method, in which we randomly
select chemical bonds from the training dataset and randomly stretch or
compress them within a reasonable range. This method achieves long-time
simulation stability comparable to AIMD simulations with almost no additional
computational cost. Moreover, this bond length stretching method can
effectively improve the long-time simulation stability of all MLFF models.

Finally, a training set that is as comprehensive and non-redundant as possible
is also extremely important for the MLFF model, as it directly determines the
prediction accuracy and training efficiency of a MLFF model. The active
learning technique provides a method for constructing a comprehensive training
set; however, active learning is merely a concept. Its effectiveness depends
on the corresponding MLFF model. In this work, the active learning concept is
combined with the bond length stretching method within the framework of our
multiscale higher-order equivariant model. The results show that this method
can efficiently and stably construct a comprehensive and low-redundancy
training set. For example, we built a training set containing only 901
samples, each with merely 120 atoms. However, based on this dataset, our
multiscale model can achieve high-precision, high-efficiency, and low-memory
consumption in long-time stable simulations of systems with hundreds of
thousands of atoms.

In summary, in this work we proposed a universal multiscale higher-order
equivariant framework for constructing a MLFF model. Moreover, within this
framework, a bond length stretching method and an active learning workflow
have been designed to realize the long-time stability of MD simulations and
efficient collection of a comprehensive, low-redundancy training set. By
employing our multiscale model, one can achieve high-precision, high-speed,
and low GPU memory consumption in long-time stable simulations of large-scale
organic systems. In our subsequent work, we will explore the application of
this multiscale model in systems with polar or charged molecules.

## Method

### Equivariant model

Incorporating the data symmetry in machine learning models can improve the
efficiency of data collection and the generalization capability of the models
[63]. For atomic systems, if the coordinates of the system rotate, quantities
like forces and dipole moments should rotate accordingly. More strictly
speaking, if the function $f:X\rightarrow Y$ is equivariant under the action
of a group of transformation $G$, the the following equation holds:

$f(D_{X}[g]x)=D_{Y}[g]f(x).$ (2)

In Eq. (2), $X$ and $Y$ are two vector spaces, $x$, $y$, and $g$ are elements
in $X$, $Y$, and $G$, respectively, and $D_{X}[g]$ and $D_{Y}[g]$ are the
transformation matrices parametrized by g in $X$ and $Y$. A natural and
effective method to ensure the equivariance of the transformation is to impose
constraints on the transformation matrices to consider the symmetric
properties of the data.[26, 64]

Han[65] categorizes the equivariant models into three types: vector-based, Lie
group-based, and irreducible representation-based. Among these three types of
models, the irreducible representation-based approach , which takes advantage
of the transformation properties of spherical harmonics $Y_{m}^{l}$, exhibits
higher-order equivariant capabilities and excels in various force field tasks
[30]. Currently, this method can be unified within the e3nn framework [32].

### Multi-scale higher-order equivariant model

Taking the MACE model as an example, we will illustrate the construction of
the multi-scale higher-order equivariant model as shown in Fig. 1. First, the
atomic numbers are initially mapped to a one-dimensional vector
$\delta_{z_{i}}$ through one-hot encoding. Subsequently, a linear
transformation initializes the node as $h_{i,c00}^{(0)}$, where the index “00”
indicates that the current node contains only scalar information. In the
subsequent network framework, node features are denoted by
$\bar{h}^{(t,\text{short/long})}_{i,cl_{2}m_{2}}$, where $i$ represents the
atomic number, $l_{2}m_{2}$ signifies the specific spherical harmonic
features, $c$ denotes the channel count, $t$ indicates the interaction layer,
and “short/long” is a shorthand to avoid repetition. This is because the
operations for short-range and long-range interactions in the interaction
layer are similar and implemented through Eqs. (4) to (8). The key difference
lies in the fact that the number of neighboring atoms to handle long-range
interactions is larger. Therefore, the number of channels $c$ and the order
$l$ for long-range interactions should be smaller than those for short-range
interactions to reduce computational complexity. Finally, short-range and
long-range information is summed up using Eq. (9).

Eq. (4) represents a linear transformation that satisfies equivariance
requirements, which only occurs between the same orders. Therefore, it is
necessary to specify the current $L$-equivariant information. Eq. (Multi-scale
higher-order equivariant model) denotes the radial embedding function, using
Bessel basis functions and polynomial smooth truncation functions.[66] Here,
$n$ represents the embedding dimension, and the truncation function for long-
range interactions should be greater than that for short-range interactions.
Eq. (6) expands the radial information to a specified dimension using a
learnable multilayer perceptron (MLP), which is related to the tensor product
in Eq. (Multi-scale higher-order equivariant model).

$\displaystyle h_{i,c00}^{(0)}$ $\displaystyle={\sum_{z}}W_{cz}\delta_{z_{i}}$
(3) $\displaystyle\bar{h}^{(t,\text{short/long})}_{i,cl_{2}m_{2}}$
$\displaystyle=\sum_{\tilde{c}}W_{c\tilde{c}l_{2}}^{(t)}h^{(t)}_{i,\tilde{c}l_{2}m_{2}}$
(4) $\displaystyle j^{\text{short/long}}_{n}(r_{ij})$
$\displaystyle=\sqrt{\frac{2}{r^{\text{short/long}}_{\text{cut}}}}\frac{\sin{\left(n\pi\frac{r_{ij}}{r^{\text{short/long}}_{\text{cut}}}\right)}}{r_{ij}}$
$\displaystyle\qquad\qquad\times f^{\text{short/long}}_{\text{cut}}(r_{ij})$
(5) $\displaystyle
R_{c\eta_{1}l_{1}l_{2}l_{3}}^{(t,\text{short/long})}(r_{ij})$
$\displaystyle={\rm
MLP}\left(\left\\{{j^{\text{short/long}}_{n}}(r_{ij})\right\\}\right)$ (6)
$\displaystyle\phi_{ij,c\eta_{1}l_{3}m_{3}}^{(t,\text{short/long})}$
$\displaystyle=\sum_{l_{1}l_{2}m_{1}m_{2}}C_{\eta_{1},l_{1}m_{1}l_{2}m_{2}}^{l_{3}m_{3}}$
$\displaystyle\times
R_{c\eta_{1}l_{1}l_{2}l_{3}}^{(t,\text{short/long})}(r_{ij})Y^{m_{1}}_{l_{1}}(\boldsymbol{\hat{r}}_{ij})\bar{h}^{(t,\text{short/long})}_{j,cl_{2}m_{2}}$
(7) $\displaystyle A_{i,cl_{3}m_{3}}^{(t,\text{short/long})}$
$\displaystyle=\sum_{\tilde{k},\eta_{1}}W_{c\tilde{c}\eta_{1}l_{3}}^{(t)}\sum_{j\in\mathcal{N}(i)}\phi_{ij,\tilde{k}\eta_{1}l_{3}m_{3}}^{(t,\text{short/long})}$
(8) $\displaystyle A_{i,cl_{3}m_{3}}^{(t)}$
$\displaystyle=(A_{i,cl_{3}m_{3}}^{(t,\text{short})}+A_{i,cl_{3}m_{3}}^{(t,\text{long})})/2$
(9) $\displaystyle{B}^{(t),\nu}_{i,\eta_{\nu}cLM}$
$\displaystyle=\sum_{{l}{m}}\mathcal{C}^{LM}_{\eta_{\nu}lm}\prod_{\xi=1}^{\nu}A_{i,cl_{\xi}m_{\xi}}^{(t)}$
(10) $\displaystyle m_{i,cLM}^{(t)}$
$\displaystyle=\sum_{\nu}\sum_{\eta_{\nu}}W_{z_{i}\eta_{\nu}cL}^{(t),\nu}{B}^{(t),\nu}_{i,\eta_{\nu}cLM}$
(11) $\displaystyle h^{(t+1)}_{i,cLM}$
$\displaystyle=\sum_{\tilde{c}}W_{cL,\tilde{c}}^{(t)}m_{i,\tilde{c}LM}^{(t)}+\sum_{\tilde{c}}W_{z_{i}cL,\tilde{c}}^{(t)}h^{(t)}_{i,\tilde{c}LM}$
(12) $\displaystyle E_{i}$
$\displaystyle=\sum_{t}\sum_{c}W^{(t)}_{c}h^{(t)}_{i,cLM}$ (13)
$\displaystyle{F}$ $\displaystyle=-\nabla\sum_{i}E_{i}$ (14)

Figure 7: |Equation calculation flow diagram.A single interaction layer
calculation is completed from $h^{(t)}_{i,cLM}$ to $h^{(t+1)}_{i,cLM}$. The
final energy is obtained by summing the linearly read-out results from each
layer’s $h^{(t+1)}_{i,cLM}$.The forces come from the negative gradient of the
final energy with respect to the coordinates.

Eq. (Multi-scale higher-order equivariant model) is the key part of the
equivariant network, merging neighboring atomic features
${h}^{(t,\text{short/long})}_{j,cl_{2}m_{2}}$ with the radial information
$R_{c\eta_{1}l_{1}l_{2}l_{3}}^{(t,\text{short/long})}(r_{ij})$ and directional
information $Y^{m_{1}}_{l_{1}}$ through convolutional filtering operations.
Eq. (8) involves pooling operations and linear transformations. It is worth
noting that the handling of long-range information modules may require zero-
padding to align with short-range node information, facilitating the final
summation in Eq. (9).

Eqs. (10) and (11) outline the process of efficiently calculating many-body
interactions in the MACE framework. Details can be found in the literature.
[27, 29] Eq. (12) involves residual connections [67], designed to update node
features. Fig. (7) illustrates the calculation process.

### Hyperparameter settings

All models were trained on a NVIDIA RTX 4090 GPU in single-GPU training using
float32 precision. Unless explicitly stated, the default hyperparameter
settings for all models in this paper are as follows: the embedding dimension
of the radial basis function is 8; smooth truncation is set to be a polynomial
envelope function with p=6; the radial MLP is [64, 64, 64]; the dimension of
the readout layer is 16; node features are represented as 64x0e; the number of
layers in the interaction layer is 2; directional information is expanded to
the 3rd order; and the cutoff radius is set to 8 Å.

For multiscale higher-order equivariant models, the hyperparameter settings
for the short-range node features, the expansion order of short-range
directional information, and the number of interaction layers remain the same
as described above. However, the long-range node features are represented as
8x0e and the long-range directional information is expanded to the first
order. The cutoff radius for short-range and long-range are set to 3 Å and 8
Å, respectively. For the Allegro model, the env_embed_multiplicity is set to
8, and latent_mlp_latent_dimensions is set to 256.

The training hyperparameters include an initial learning rate of 0.01 and a
ReduceLROnPlateau scheduler, which reduces the learning rate when the
validation loss does not improve over a certain number of epochs. To update
the evaluation and final model weights of the validation dataset, an
exponential moving average with a weight of 0.99 is applied. The optimizer is
Adam, and the total number of training epochs is set to 4000. Energy is
normalized by the moving average of the potential energy. The loss function is
as follows:

$\displaystyle\mathcal{L}=$
$\displaystyle\frac{\lambda_{E}}{B}\sum_{b=1}^{B}\left(\frac{E_{b}-\hat{E}_{b}}{N_{b}}\right)^{2}$
$\displaystyle+\frac{\lambda_{F}}{3B}\sum_{b=1}^{B}\frac{1}{N_{b}}\sum_{i_{b},\alpha=1}^{N_{b},3}\left(-\frac{\partial
E_{b}}{\partial r_{i_{b},\alpha}}-\hat{F}_{i_{b},\alpha}\right)^{2}$ (15)

The initial weight values of the force and energy follow common settings [27,
29], where the force weight is set to 1000, and the energy weight is the
number of atoms of the system. For MACE and MS-MACE models, Stochastic Weight
Averaging (SWA) [68, 69] is enabled at 75% of the total iteration count. After
initiating SWA, as in [29], the energy weights and force weights of the loss
function are reset to 1000 and 10, respectively.

### DFT settings

Since this paper requires calculations with thousands of atoms and the need to
perform AIMD, the energy and forces of all data sets are calculated using the
lower computational cost but good accuracy method, i.e., the semi-empirical
GFN2-xTB level of theory[55]. This choice is made to balance computational
cost and accuracy. The convergence level follows the default settings of the
xtb software [70].

It is essential to note that in the process of data collection, a pretrained
MLFF based on the multiscale higher-order equivariant model is used to perform
high-precision MD, which is significantly faster than the AIMD simulation.
Therefore, the data sampling stage is not computationally intensive. In
practice, the primary cost of generating data labeling still lies in labeling
using DFT. More accurate quantum-mechanical calculations can be employed based
on specific requirements.

### Molecular dynamics settings

Periodic AIMD simulations based on the semi-empirical GFN2-xTB level of theory
are implemented using the DFTB+ software [71]. The MD simulations based on the
MLFF models are executed using version 1.0 of the OpenMM-Torch plugin and
version 8.0.0 of the OpenMM software.

### Test of simulation speed and memory consumption

Before benchmarking, it is essential to preheat the GPU using partial data.
The timer employs torch.cuda.Event(), and the dataset with different numbers
of atoms should contain at least 20 samples. Each sample is run 50 times, and
the average value is taken as the final result. Similarly,
torch.cuda.max_memory_allocated() from the official torch is used to record
the peak GPU memory consumption. Before each recording round,
torch.cuda.reset_peak_memory_stats() and torch.cuda.empty_cache() are employed
to reset information and release excess cache. The GPU of the test platform is
a 40G A100.

## Code availability

The code and sample scripts will be released after review.

## Data availability

The dataset will be released after review.

## Competing interests

The authors declare no competing interests.

## Acknowledgements

This work was supported by the National Key R&D Program of China(No.
2021YFB3803200) and the National Natural Science Foundation of China under
Grant No. 22273112.

## References

  *   * Zhang et al. [2018] Zhang, J., Wang, X., Zhu, Y., Shi, T., Tang, Z., Li, M., Liao, G.: Molecular dynamics simulation of the melting behavior of copper nanorod. Computational Materials Science 143, 248–254 (2018)
  * Zhong et al. [2023] Zhong, Z., Du, G., Wang, Y., Jiang, J.: Phase behaviors of dialkyldimethylammonium bromide bilayers. Langmuir 39(31), 11081–11089 (2023)
  * Ma et al. [2023] Ma, L., Zhong, Z., Hu, J., Qing, L., Jiang, J.: Long-lived weak ion pairs in ionic liquids: An insight from all-atom molecular dynamics simulations. The Journal of Physical Chemistry B (2023)
  * Perilla and Schulten [2017] Perilla, J.R., Schulten, K.: Physical properties of the hiv-1 capsid from all-atom molecular dynamics simulations. Nature communications 8(1), 15959 (2017)
  * Mouvet et al. [2022] Mouvet, F., Villard, J., Bolnykh, V., Rothlisberger, U.: Recent advances in first-principles based molecular dynamics. Accounts of Chemical Research 55(3), 221–230 (2022)
  * Behler and Parrinello [2007] Behler, J., Parrinello, M.: Generalized neural-network representation of high-dimensional potential-energy surfaces. Physical review letters 98(14), 146401 (2007)
  * Chmiela et al. [2017] Chmiela, S., Tkatchenko, A., Sauceda, H.E., Poltavsky, I., Schütt, K.T., Müller, K.-R.: Machine learning of accurate energy-conserving molecular force fields. Science Advances 3(5), 1603015 (2017) https://doi.org/10.1126/sciadv.1603015
  * Chmiela et al. [2023] Chmiela, S., Vassilev-Galindo, V., Unke, O.T., Kabylda, A., Sauceda, H.E., Tkatchenko, A., Müller, K.-R.: Accurate global machine learning force fields for molecules with hundreds of atoms. Science Advances 9(2), 0873 (2023)
  * Tran et al. [2023] Tran, R., Lan, J., Shuaibi, M., Wood, B.M., Goyal, S., Das, A., Heras-Domingo, J., Kolluru, A., Rizvi, A., Shoghi, N., et al.: The open catalyst 2022 (oc22) dataset and challenges for oxide electrocatalysts. ACS Catalysis 13(5), 3066–3084 (2023)
  * Schütt et al. [2017] Schütt, K., Kindermans, P.-J., Sauceda Felix, H.E., Chmiela, S., Tkatchenko, A., Müller, K.-R.: Schnet: A continuous-filter convolutional neural network for modeling quantum interactions. Advances in neural information processing systems 30 (2017)
  * Wang et al. [2018] Wang, H., Zhang, L., Han, J., Weinan, E.: Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics. Computer Physics Communications 228, 178–184 (2018)
  * Schütt et al. [2017] Schütt, K.T., Arbabzadah, F., Chmiela, S., Müller, K.R., Tkatchenko, A.: Quantum-chemical insights from deep tensor neural networks. Nature communications 8(1), 13890 (2017)
  * Unke and Meuwly [2019] Unke, O.T., Meuwly, M.: Physnet: A neural network for predicting energies, forces, dipole moments, and partial charges. Journal of chemical theory and computation 15(6), 3678–3693 (2019)
  * Wang et al. [2022] Wang, L., Liu, Y., Lin, Y., Liu, H., Ji, S.: Comenet: Towards complete and efficient message passing for 3d molecular graphs. Advances in Neural Information Processing Systems 35, 650–664 (2022)
  * Liu et al. [2022] Liu, Y., Wang, L., Liu, M., Lin, Y., Zhang, X., Oztekin, B., Ji, S.: Spherical message passing for 3d molecular graphs. In: International Conference on Learning Representations (2022). https://openreview.net/forum?id=givsRXsOt9r
  * Satorras et al. [2021] Satorras, V.G., Hoogeboom, E., Welling, M.: E (n) equivariant graph neural networks. In: International Conference on Machine Learning, pp. 9323–9332 (2021). PMLR
  * Schütt et al. [2021] Schütt, K., Unke, O., Gastegger, M.: Equivariant message passing for the prediction of tensorial properties and molecular spectra. In: International Conference on Machine Learning, pp. 9377–9388 (2021). PMLR
  * Jing et al. [2021] Jing, B., Eismann, S., Suriana, P., Townshend, R.J.L., Dror, R.: Learning from protein structure with geometric vector perceptrons. In: International Conference on Learning Representations (2021). https://openreview.net/forum?id=1YLJDvSx6J4
  * Le et al. [2022] Le, T., Noé, F., Clevert, D.-A.: Equivariant Graph Attention Networks for Molecular Property Prediction (2022)
  * Thomas et al. [2018] Thomas, N., Smidt, T., Kearnes, S., Yang, L., Li, L., Kohlhoff, K., Riley, P.: Tensor field networks: Rotation- and translation-equivariant neural networks for 3D point clouds (2018)
  * Anderson et al. [2019] Anderson, B., Hy, T.S., Kondor, R.: Cormorant: Covariant molecular neural networks. Advances in neural information processing systems 32 (2019)
  * Brandstetter et al. [2022] Brandstetter, J., Hesselink, R., Pol, E., Bekkers, E.J., Welling, M.: Geometric and physical quantities improve e(3) equivariant message passing. In: International Conference on Learning Representations (2022). https://openreview.net/forum?id=_xwr8gOBeV1
  * Liao and Smidt [2023] Liao, Y.-L., Smidt, T.: Equiformer: Equivariant graph attention transformer for 3d atomistic graphs. In: The Eleventh International Conference on Learning Representations (2023). https://openreview.net/forum?id=KwmPfARgOTD
  * Batzner et al. [2022] Batzner, S., Musaelian, A., Sun, L., Geiger, M., Mailoa, J.P., Kornbluth, M., Molinari, N., Smidt, T.E., Kozinsky, B.: E (3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature communications 13(1), 2453 (2022)
  * Musaelian et al. [2023] Musaelian, A., Batzner, S., Johansson, A., Sun, L., Owen, C.J., Kornbluth, M., Kozinsky, B.: Learning local equivariant representations for large-scale atomistic dynamics. Nature Communications 14(1), 579 (2023)
  * Batatia et al. [2022a] Batatia, I., Batzner, S., Kovács, D.P., Musaelian, A., Simm, G.N.C., Drautz, R., Ortner, C., Kozinsky, B., Csányi, G.: The Design Space of E(3)-Equivariant Atom-Centered Interatomic Potentials (2022)
  * Batatia et al. [2022b] Batatia, I., Kovacs, D.P., Simm, G., Ortner, C., Csanyi, G.: Mace: Higher order equivariant message passing neural networks for fast and accurate force fields. In: Koyejo, S., Mohamed, S., Agarwal, A., Belgrave, D., Cho, K., Oh, A. (eds.) Advances in Neural Information Processing Systems, vol. 35, pp. 11423–11436 (2022)
  * Zhang et al. [2023] Zhang, X., Wang, L., Helwig, J., Luo, Y., Fu, C., Xie, Y., Liu, M., Lin, Y., Xu, Z., Yan, K., et al.: Artificial intelligence for science in quantum, atomistic, and continuum systems. arXiv preprint arXiv:2307.08423 (2023)
  * Kovács et al. [2023] Kovács, D.P., Batatia, I., Arany, E.S., Csányi, G.: Evaluation of the MACE force field architecture: From medicinal chemistry to materials science. The Journal of Chemical Physics 159(4), 044118 (2023) https://doi.org/10.1063/5.0155322
  * Fu et al. [2022] Fu, X., Wu, Z., Wang, W., Xie, T., Keten, S., Gomez-Bombarelli, R., Jaakkola, T.: Forces are not enough: Benchmark and critical evaluation for machine learning force fields with molecular simulations. arXiv preprint arXiv:2210.07237 (2022)
  * Joshi et al. [2023] Joshi, C.K., Bodnar, C., Mathis, S.V., Cohen, T., Lio, P.: On the expressive power of geometric graph neural networks. arXiv preprint arXiv:2301.09308 (2023)
  * Geiger and Smidt [2022] Geiger, M., Smidt, T.: e3nn: Euclidean neural networks. arXiv preprint arXiv:2207.09453 (2022)
  * Rackers et al. [2023] Rackers, J.A., Tecot, L., Geiger, M., Smidt, T.E.: A recipe for cracking the quantum scaling limit with machine learned electron densities. Machine Learning: Science and Technology 4(1), 015027 (2023)
  * Tokita and Behler [2023] Tokita, A.M., Behler, J.: Tutorial: How to train a neural network potential. arXiv preprint arXiv:2308.08859 (2023)
  * Klicpera et al. [2021] Klicpera, J., Becker, F., Günnemann, S.: Gemnet: Universal directional graph neural networks for molecules. In: Beygelzimer, A., Dauphin, Y., Liang, P., Vaughan, J.W. (eds.) Advances in Neural Information Processing Systems (2021). https://openreview.net/forum?id=HS_sOaxS9K-
  * Wen et al. [2022] Wen, T., Zhang, L., Wang, H., Weinan, E., Srolovitz, D.J.: Deep potentials for materials science. Materials Futures 1(2), 022601 (2022)
  * Gao et al. [2020] Gao, X., Ramezanghorbani, F., Isayev, O., Smith, J.S., Roitberg, A.E.: Torchani: A free and open source pytorch-based deep learning implementation of the ani neural network potentials. Journal of chemical information and modeling 60(7), 3408–3415 (2020)
  * Wang et al. [2023] Wang, Z., Wu, H., Sun, L., He, X., Liu, Z., Shao, B., Wang, T., Liu, T.-Y.: Improving machine learning force fields for molecular dynamics simulations with fine-grained force metrics. The Journal of chemical physics 159(3), 035101 (2023)
  * Li et al. [2023] Li, Y., Wang, Y., Huang, L., Yang, H., Wei, X., Zhang, J., Wang, T., Wang, Z., Shao, B., Liu, T.-Y.: Long-short-range message-passing: A physics-informed framework to capture non-local interaction for scalable molecular dynamics simulation. arXiv preprint arXiv:2304.13542 (2023)
  * Di Giovanni et al. [2023] Di Giovanni, F., Giusti, L., Barbero, F., Luise, G., Lio, P., Bronstein, M.M.: On over-squashing in message passing neural networks: The impact of width, depth, and topology. In: International Conference on Machine Learning, pp. 7865–7885 (2023). PMLR
  * Zhang et al. [2022] Zhang, L., Wang, H., Muniz, M.C., Panagiotopoulos, A.Z., Car, R., et al.: A deep potential model with long-range electrostatic interactions. The Journal of Chemical Physics 156(12) (2022)
  * Ko et al. [2021] Ko, T.W., Finkler, J.A., Goedecker, S., Behler, J.: A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer. Nature communications 12(1), 398 (2021)
  * Batatia et al. [2023] Batatia, I., Schaaf, L.L., Chen, H., Csányi, G., Ortner, C., Faber, F.A.: Equivariant matrix function neural networks. arXiv preprint arXiv:2310.10434 (2023)
  * Zuo et al. [2020] Zuo, Y., Chen, C., Li, X., Deng, Z., Chen, Y., Behler, J., Csányi, G., Shapeev, A.V., Thompson, A.P., Wood, M.A., et al.: Performance and cost assessment of machine learning interatomic potentials. The Journal of Physical Chemistry A 124(4), 731–745 (2020)
  * Huang et al. [2021] Huang, J., Zhang, L., Wang, H., Zhao, J., Cheng, J., et al.: Deep potential generation scheme and simulation protocol for the li10gep2s12-type superionic conductors. The Journal of Chemical Physics 154(9) (2021)
  * Chung et al. [2015] Chung, L.W., Sameera, W., Ramozzi, R., Page, A.J., Hatanaka, M., Petrova, G.P., Harris, T.V., Li, X., Ke, Z., Liu, F., et al.: The oniom method and its applications. Chemical reviews 115(12), 5678–5796 (2015)
  * Collins and Bettens [2015] Collins, M.A., Bettens, R.P.: Energy-based molecular fragmentation methods. Chemical reviews 115(12), 5607–5642 (2015)
  * Grisafi and Ceriotti [2019] Grisafi, A., Ceriotti, M.: Incorporating long-range physics in atomic-scale machine learning. The Journal of chemical physics 151(20) (2019)
  * Fedik et al. [2022] Fedik, N., Zubatyuk, R., Kulichenko, M., Lubbers, N., Smith, J.S., Nebgen, B., Messerly, R., Li, Y.W., Boldyrev, A.I., Barros, K., et al.: Extending machine learning beyond interatomic potentials for predicting molecular properties. Nature Reviews Chemistry 6(9), 653–672 (2022)
  * Prodan and Kohn [2005] Prodan, E., Kohn, W.: Nearsightedness of electronic matter. Proceedings of the National Academy of Sciences 102(33), 11635–11638 (2005) https://doi.org/10.1073/pnas.0505436102
  * Chmiela et al. [2019] Chmiela, S., Sauceda, H.E., Poltavsky, I., Müller, K.-R., Tkatchenko, A.: sgdml: Constructing accurate and data efficient molecular force fields using machine learning. Computer Physics Communications 240, 38–45 (2019)
  * Batzner et al. [2022] Batzner, S., Musaelian, A., Sun, L., Geiger, M., Mailoa, J.P., Kornbluth, M., Molinari, N., Smidt, T.E., Kozinsky, B.: E (3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature communications 13(1), 2453 (2022)
  * Brandstetter et al. [2021] Brandstetter, J., Hesselink, R., Pol, E., Bekkers, E.J., Welling, M.: Geometric and physical quantities improve e (3) equivariant message passing. arXiv preprint arXiv:2110.02905 (2021)
  * Lubbers et al. [2018] Lubbers, N., Smith, J.S., Barros, K.: Hierarchical modeling of molecular energies using a deep neural network. The Journal of chemical physics 148(24) (2018)
  * Bannwarth et al. [2019] Bannwarth, C., Ehlert, S., Grimme, S.: Gfn2-xtb—an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions. Journal of chemical theory and computation 15(3), 1652–1671 (2019)
  * Stocker et al. [2022] Stocker, S., Gasteiger, J., Becker, F., Günnemann, S., Margraf, J.T.: How robust are modern graph neural network potentials in long and hot molecular dynamics simulations? Machine Learning: Science and Technology 3(4), 045010 (2022)
  * Liu et al. [2023] Liu, Y., He, X., Mo, Y.: Discrepancies and the error evaluation metrics for machine learning interatomic potentials. npj Computational Materials 9(1), 174 (2023)
  * Schütt et al. [2020] Schütt, K.T., Chmiela, S., Von Lilienfeld, O.A., Tkatchenko, A., Tsuda, K., Müller, K.-R.: Machine learning meets quantum physics. Lecture Notes in Physics (2020)
  * Zhang et al. [2020] Zhang, Y., Wang, H., Chen, W., Zeng, J., Zhang, L., Wang, H., Weinan, E.: Dp-gen: A concurrent learning platform for the generation of reliable deep learning based potential energy models. Computer Physics Communications 253, 107206 (2020)
  * Schran et al. [2020] Schran, C., Brezina, K., Marsalek, O.: Committee neural network potentials control generalization errors and enable active learning. The Journal of Chemical Physics 153(10) (2020)
  * Wang et al. [2004] Wang, J., Wolf, R.M., Caldwell, J.W., Kollman, P.A., Case, D.A.: Development and testing of a general amber force field. Journal of computational chemistry 25(9), 1157–1174 (2004)
  * Huguenin-Dumittan et al. [2023] Huguenin-Dumittan, K.K., Loche, P., Haoran, N., Ceriotti, M.: Physics-inspired equivariant descriptors of nonbonded interactions. The Journal of Physical Chemistry Letters 14, 9612–9618 (2023)
  * Cohen and Welling [2016] Cohen, T., Welling, M.: Group equivariant convolutional networks. In: International Conference on Machine Learning, pp. 2990–2999 (2016). PMLR
  * Darby et al. [2023] Darby, J.P., Kovács, D.P., Batatia, I., Caro, M.A., Hart, G.L., Ortner, C., Csányi, G.: Tensor-reduced atomic density representations. Physical Review Letters 131(2), 028001 (2023)
  * Han et al. [2022] Han, J., Rong, Y., Xu, T., Huang, W.: Geometrically equivariant graph neural networks: A survey. arXiv preprint arXiv:2202.07230 (2022)
  * Gasteiger et al. [2020] Gasteiger, J., Groß, J., Günnemann, S.: Directional message passing for molecular graphs. In: International Conference on Learning Representations (ICLR) (2020). https://arxiv.org/abs/2003.03123
  * He et al. [2016] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 770–778 (2016)
  * Izmailov et al. [2018] Izmailov, P., Podoprikhin, D., Garipov, T., Vetrov, D., Wilson, A.G.: Averaging weights leads to wider optima and better generalization. arXiv preprint arXiv:1803.05407 (2018)
  * Athiwaratkun et al. [2018] Athiwaratkun, B., Finzi, M., Izmailov, P., Wilson, A.G.: There are many consistent explanations of unlabeled data: Why you should average. arXiv preprint arXiv:1806.05594 (2018)
  * Bannwarth et al. [2021] Bannwarth, C., Caldeweyher, E., Ehlert, S., Hansen, A., Pracht, P., Seibert, J., Spicher, S., Grimme, S.: Extended tight-binding quantum chemistry methods. Wiley Interdisciplinary Reviews: Computational Molecular Science 11(2), 1493 (2021)
  * Hourahine et al. [2020] Hourahine, B., Aradi, B., Blum, V., Bonafé, F., Buccheri, A., Camacho, C., Cevallos, C., Deshaye, M., Dumitrică, T., Dominguez, A., et al.: Dftb+, a software package for efficient approximate density functional theory based atomistic simulations. The Journal of chemical physics 152(12) (2020)