# Force Field Analysis Software and Tools (FFAST): Assessing Machine Learning
Force Fields Under the Microscope

Gregory Fonseca  Department of Physics and Materials Science, University of
Luxembourg, L-1511 Luxembourg City, Luxembourg  Igor Poltavsky  Department of
Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg
City, Luxembourg  Alexandre Tkatchenko<EMAIL_ADDRESS>Department
of Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg
City, Luxembourg

(11 August 2023)

###### Abstract

As the sophistication of Machine Learning Force Fields (MLFF) increases to
match the complexity of extended molecules and materials, so does the need for
tools to properly analyze and assess the practical performance of MLFFs. To go
beyond average error metrics and into a complete picture of a model’s
applicability and limitations, we develop FFAST (Force Field Analysis Software
and Tools): a cross-platform software package designed to gain detailed
insights into a model’s performance and limitations, complete with an easy-to-
use graphical user interface. The software allows the user to gauge the
performance of many popular state-of-the-art MLFF models on various popular
dataset types, providing general prediction error overviews, outlier detection
mechanisms, atom-projected errors, and more. It has a 3D visualizer to find
and picture problematic configurations, atoms, or clusters in a large dataset.
In this paper, the example of the MACE and Nequip models are used on two
datasets of interest – stachyose and docosahexaenoic acid (DHA) – to
illustrate the use cases of the software. With it, it was found that carbons
and oxygens involved in or near glycosidic bonds inside the stachyose molecule
present increased prediction errors. In addition, prediction errors on DHA
rise as the molecule folds, especially for the carboxylic group at the edge of
the molecule. We emphasize the need for a systematic assessment of MLFF models
for ensuring their successful application to study the dynamics of molecules
and materials.

##  I Introduction

Many of the cornerstones of modern advancements in industry and technology
rely on the introduction of novel or optimized materials, newly designed
drugs, and the understanding of physico-chemical phenomena on a molecular
level. Here, computer simulations and predictions of properties play a pivotal
role in enabling, enhancing, or accelerating research and development. Lee
_et al._ (2019); Ramos  _et al._ (2019); Mansbach  _et al._ (2019); Jing  _et
al._ (2019); Ekins (2016); Smalley (2017); Elton  _et al._ (2018, 2019); Faber
_et al._ (2016); Naserifar  _et al._ (2021); Cova and Pais (2019)

Amidst the computational tools, Machine Learning Force Fields (MLFF) are
steadily rising in popularity, with a wide variety of models achieving
remarkable predictive accuracy Batzner  _et al._ (2022); Geiger and Smidt
(2022); Batatia  _et al._ (2022a, b); Bartók  _et al._ (2010); Bartók, Kondor,
and Csányi (2013); Behler (2016); Behler and Parrinello (2007); Schütt  _et
al._ (2017a, b, 2018, 2019); Rupp  _et al._ (2012); Bartók  _et al._ (2017);
Unke and Meuwly (2019); Unke  _et al._ (2021a, b); Chmiela  _et al._ (2018,
2019); Sauceda  _et al._ (2022); Qiao  _et al._ (2020); Smith, Isayev, and
Roitberg (2017); Devereux  _et al._ (2020); Deringer  _et al._ (2018); Ryczko
_et al._ (2018); Frank, Unke, and Müller (2023); Christensen  _et al._ (2020);
Wang  _et al._ (2019); Ko  _et al._ (2021a); Kocer, Ko, and Behler (2022);
Gasteiger, Groß, and Günnemann (2022); Gasteiger, Becker, and Günnemann
(2022); Deringer  _et al._ (2021); Bartók  _et al._ (2018); Thomas  _et al._
(2018); Khorshidi and Peterson (2016); Li, Kermode, and De Vita (2015); Keith
_et al._ (2021); Ko  _et al._ (2021b); Yao  _et al._ (2018); Grisafi and
Ceriotti (2019) for various molecules and materials of ever-increasing sizes
and complexity Hoja  _et al._ (2021); Blum and Reymond (2009); Ramakrishnan
_et al._ (2014); Ruddigkeit  _et al._ (2012); Chmiela  _et al._ (2017, 2023);
Smith  _et al._ (2020). The ultimate goal of MLFF is to bring the quality of
expensive ab initio methods to systems of larger scales, which is currently
mostly feasible for efficient but significantly less accurate empirical
mechanistic force fields. MLFFs aim to advance the FF development by learning
highly accurate data generated by state-of-the-art quantum chemistry methods
to then reproduce their results in a fraction of the computational time.
However, in order to use these models in practice, it is crucial to understand
which tasks a given MLFF is appropriate for. ML models are very sensitive to
their training data, the training procedure, hyperparameters and many other
details: two models with similar overall accuracy but trained on different
data, hyperparameters, or architectures can still present fundamental
differences in actual applications (for example molecular dynamics) that are
not captured by simple error metrics Fonseca  _et al._ (2021); Fu  _et al._
(2022).

For practical applications, it is important to know that our model is accurate
and when it is likely to fail. This information is crucial, given that a
single configuration that the model predicts inaccurately can move the system
under study into an improbable or unphysical state, which could affect the
remaining simulation. The probability of such a misprediction is steadily
growing with increasing size, chemical, and structural complexity of the
systems under study. While no silver bullets exist to make ML models perfectly
stable, analyzing models beyond usual error metrics can avoid many pitfalls
and limitations. This is why we developed FFAST (Force Field Analysis Software
and Tools): a cross-platform software package designed to gain detailed
insights into a model’s performance and limitations, complete with an easy-to-
use graphical user interface. Note that while this paper’s main focus is
MLFFs, all the described tools are equally applicable to empirical force
fields or – in fact – force fields of any kind.

The article is organized as follows: in the Software section, the software and
its components are listed and briefly explained. In the Workflow section, part
of those components are explored in a realistic use-case from start to finish.
The Applications section applies the above workflow to actual datasets and
models and analyses the results. Conclusions and outlooks can be found in the
final section.

##  II Software

Out of the box, FFAST can load a variety of MLFF models (currently supported
sGDML Chmiela  _et al._ (2018), SchNet Schütt  _et al._ (2017a), Nequip
Batzner  _et al._ (2022); Geiger and Smidt (2022), MACE Batatia  _et al._
(2022a, b), SpookyNet Unke  _et al._ (2021a), and pre-predicted
forces/energies) as well as dataset formats (currently supported .xyz, .npz,
and .db). The ML models are used to generate predictions of the energy and
forces on the loaded datasets. To prevent excessive loading times when
handling large datasets or expensive models, this step can be accomplished in
headless mode (i.e., without a graphical interface) to pre-compute the
predictions externally (e.g., on a high-performance computer).

Once loaded, various analysis tools are available in a user-friendly
interface.

  * •

Error distributions on both energy and force predictions are available, able
to visualize multiple dataset/model combinations for easy comparisons. By
default, mean average errors (MAE) are used, and the plots visualize a
gaussian kernel-density estimate.

  * •

Error timelines show the MAE across the given time-ordered dataset (e.g., a
dynamic). An adjustable smoothing factor is provided to disregard fluctuations
by averaging each point over the given number of neighbors.

  * •

Plots of cluster errors provide a way to find problematic regions of
configurational space. By default, agglomerative clustering on pairwise
interatomic distances is followed up by KMeans on energies for a total of 40
clusters.

  * •

Error scatter plots are available to find outliers quickly, both for energies
and forces.

  * •

Distribution of atomic errors can be used to determine the difference in
prediction between chosen elements. The procedure is otherwise identical to
error distributions.

The user can zoom in on most plots and choose to create subdatasets. These are
subsets of loaded datasets to be further analyzed by the software. Subsets and
datasets can be saved in any format that is provided by the software in case
external editing or analysis is needed.

The program comes equipped with a 3D visualization tool. The molecular
structures inside of a dataset or subset can be viewed in one or more
interactive windows. The visualizer comes with convenience features like
aligning geometries along a chosen plane for easier visualization, choosing
the bonds to be visualized, or animating the molecule throughout the selected
dataset. One can extract information on the geometry such as atomic distances,
angles or dihedrals. The user can enable the plotting of force prediction
errors on each atom of the molecule, either at the given timeframe or averaged
throughout the dataset/subset. Finally, atoms of interest can be selected and
analyzed independently in all analysis tools above.

One can also use the above-described tools to analyze the underlying reference
dataset by creating a dummy model. In this case, FFAST will operate on the
reference energy and forces instead of the difference between the reference
data and ML models’ prediction.

The software is designed to be modular, allowing users comfortable with Python
to add features they need for their workflow. The workflow is streamlined such
that swapping between software codes and separate scripts for different models
is eliminated as much as possible. The only dependencies are those of the
external ML models themselves (subject to their own installation process) and
readily available Python packages. As such, FFAST provides a platform that
beginners and experts alike can use to quickly and precisely assess their
force field performance in an unbiased way.

##  III Workflow

To illustrate the usefulness of the program, an example workflow is presented
in this chapter. It is important to note that this use-case serves to guide
and highlight important points of the software, not to be an all-inclusive
exposition of all its features.

As a first step, one or multiple ML models and datasets of interest can be
loaded. If the datasets are large or the MLFFs are expensive, the program can
be used in headless mode on a supercomputer to pre-compute the forces and
energies. This allows the user to avoid unnecessary recalculations and
prolonged loading times while using the software. For all applications in this
paper, data for molecular systems are taken from the MD22 dataset Chmiela  _et
al._ (2023). This chapter specifically focuses on its stachyose trajectory as
an example. The particular selection of models used as illustrations in this
paper are Nequip and MACE. The parameters used for all models reflect those
recommended in the original paper or official software description of the
respective packages. In all cases, all models are trained on the same 1000
training points for each dataset, respectively.

Prediction error overview: The first screen to greet the user contains basic
information providing an overview of a model’s overall performance. This
overview includes basic error metrics such as mean average errors (MAEs), root
mean squared errors (RMSEs), error distributions, and error timelines
throughout a dataset. An example is given in Fig. 1. Error distributions are
one of the fastest ways to determine a model’s performance. From a glance,
they provide a rough estimate of the accuracy across a given dataset and what
to expect for the low and high end of the error spectrum. Furthermore, one can
look at the distribution patterns (such as deviation from a normal curve) to
understand if a systematic error occurs on specific data subsets.

Figure 1: Basic error screen including timelines and distributions for both
energy and forces mean average errors across an entire stachyose dataset. In
blue, a Nequip model trained on 1000 points of the dataset is used.

Error timelines show the prediction errors throughout the dataset’s lifecycle.
This is mostly useful if the data is generated from molecular dynamic
simulations, as it allows the user to observe general trends or rough moments
in time where an MLFF model’s performance might have been lacking. The average
window size can be adjusted manually to find a good balance between coarseness
and precision.

Outlier detection: Many tools are provided to find outliers in the dataset or
in the prediction error profile. While average errors across entire datasets
are helpful, they do not provide any indication as to why and where the FF
fails. As such, it is essential to determine which configurations cause issues
and whether those are outliers in the dataset or important representative
geometries.

Correlation scatters are scatter plots of the predicted values against the
true values. While simple, they allow users to visually determine which
points, likely far from the expected linear correlation, are worth a second
look. Besides, the general shape of a correlation plot gives a good indication
of the stability of the model, as in practice, a single bad prediction during
a time step can result in the entire dynamics steering into the non-physical
territory.

A more in-depth way to visualize a model’s stability across areas of
configurational space available in a given dataset is through the usage of
clustering algorithms Fonseca  _et al._ (2021), see example on Fig. 2. In
essence, the dataset is split into a select number of distinct clusters, each
containing a varying amount of configurations. The groups are chosen such as
to maximize the similarity between configurations of the same group, thereby
splitting the entire datasets into a digestible number of qualitatively
different configurations. Error analysis on those different clusters provides
an overview akin to an error distribution; however, its discrete nature
focuses on types of configuration rather than single points. Thus the user can
apply the chemical intuition and knowledge of the molecule at hand to
potentially find a) why a model fails on certain clusters but not others, b)
what the application range of the model is, or c) outliers in the dataset.

Figure 2: Mean average error of force predictions of a Nequip model trained on
1000 points for 50 different clusters (ordered by ascending error) of a
stachyose dataset.

Atomic errors: While the previous plots focus on finding single or groups of
configurations with particular prediction error patterns, there is also merit
to focusing on models’ performance for different atom types. With the
increasing chemical composition and structural complexity of the systems of
interest for MLFFs, atoms of the same type interact differently with their
environment based on their environment composition. Thus, one can expect that
ML model predictions vary across atom types. FFAST provides the advantage of
visualizing a given model’s force prediction error distribution across
different chemical elements constituting the system. Alternatively, one can
view the error distribution for a selected chemical element across multiple ML
models or datasets. An example can be found in Fig. 3

Figure 3: Mean average error density of the force predictions of a Nequip
model trained on 1000 points of a stachyose dataset. The densities are done
separately for each atom type. Colors correspond to elements (Hydrogen: white,
Oxygen: red, Carbon: grey)

3D visualization : All the previous features offered a way to get a general
understanding of a FF performance and find outliers for the models. However,
to extract chemical intuition from these results, it is paramount to visualize
them in 3D space. FFAST allows most of its plots mentioned above to be used
for visualization purposes. Zoom into any subarea of interest and visualize
only configurations within that region. Beyond giving form to the points,
areas, and clusters of interest, the 3D visualization tool allows for
configuration-by-configuration and atom-by-atom basis analysis. One can view
the force prediction error for each atom on an individual geometry, or errors
averaged throughout all selected configurations. That way, if a specific
mechanism is at the root of poor predictions, combining outlier detection (to
find the geometries involved) and atomic force prediction error (to see the
atoms involved) can provide a path to revealing its cause.

Figure 4: Cutout of the 3D visualiser on the stachyose molecule. (Left) atoms
colored according to their element. (Right) atoms colored according to their
respective mean average force error as predicted by a Nequip model trained on
1000 points.

##  IV Applications

After demonstrating the main features of the FFAST, we will apply it to
analyze in detail the capabilities of two state-of-the-art ML architectures,
namely Nequip and MACE, to reconstruct PESs and FFs of flexible organic
molecules. While both ML models can easily reach below chemical accuracy MAEs
on precomputed test datasets, obtaining stable long-time dynamics is a
nontrivial task. Therefore, understanding the origin of prediction errors on
an atomistic level is crucial for accurate applications and model prediction
improvements.

###  IV.1 Stachyose

This subsection focuses on the stachyose molecule (C24H42O21), consisting of
87 carbon, oxygen, and hydrogen atoms. The primary ML architecture is the
Nequip model trained on the abovementioned dataset and analyzed using the
FFAST software. The basic errors screen (see Fig. 1) reveals some noteworthy
trends. Most of the configurations in the dataset fall within chemical
accuracy for the chosen model, as indicated by the energy MAE distribution’s
main contributor lying below 1 $kcal/mol$. This also applies to the force MAE
distribution with the overall MAE of 0.58 $kcal/(mol\,\mbox{\AA})$.
Interestingly, the force error distribution demonstrates a well-defined double
peak shape, suggesting non-equal force reconstruction for different system
components. Also, the forces and energy MAE timelines shows oscillating
behavior clearly distinguishing sets of configurations, reflecting that the
molecular dynamics that generated the dataset explored qualitatively different
molecular structures throughout the simulation instead of oscillating in or
around an equilibrium state.

The atomic error plot in Fig. 3 expands on the overall distributions and
distinctly shows a fundamental difference between hydrogen and other elements.
The overall MAE and RMSE of each atom type can be found in Table 1.
Importantly, the well-predicted hydrogen atoms, making up almost half of the
molecule, highly influence the overall prediction errors. This also largely
explains the double peak structure found in the overall force error
distribution discussed above. Nevertheless, the secondary peaks in the
distributions are also visible on a per-atom basis, notably for carbon but
also in oxygen.

| H | C | $\text{C}_{b}$ | $\text{C}_{r}$ | $\text{C}_{s}$ | O | $\text{O}_{b}$ | $\text{O}_{r}$ | $\text{O}_{s}$ | All
---|---|---|---|---|---|---|---|---|---|---
MAE | 0.40 | 0.83 | 0.92 | 0.78 | 0.77 | 0.66 | 0.89 | 0.82 | 0.57 | 0.58
RMSE | 0.55 | 1.10 | 1.23 | 1.04 | 1.02 | 0.90 | 1.20 | 1.10 | 0.76 | 0.82
Table 1: Overall force MAE and RMSE in $kcal/(mol\,\mbox{\AA})$ for each atom
type in the Stachyose dataset, as predicted by a Nequip model of 1000 training
points. The subscripts $b$, $r$, and $s$ indicate filtering for only atoms
touching a glycosidic bond, inside the rest of the ring, or in a side chain,
respectively.

The reason for the double-peak structure of carbon and oxygen force error
distributions is revealed when looking at the average force prediction error
for every individual atom in the 3D visualizer; see Fig. 4 (right). The
molecule consists of a main chain of one furanose (5-membered carbohydrate
ring with oxygen) and three pyranoses (6-membered carbohydrate ring with
oxygen). Upon inspection, one can notice that the carbons involved in or
directly touching the glycosidic bonds between the rings have notably worse
errors than those elsewhere on the rings (by about $\sim 18\%$). A similar
situation is found for the oxygens inside said bonds, which have worse
predictions than those in the rings or, more significantly, than those in
side-chain hydroxyl groups (by about $\sim 56\%$, see Table 1). Remembering
that a model’s actual simulation performance is largely defined by its most
problematic atoms rather than an overall mean error is important. Hence, it
would be wise to consider this particular model’s practical MAE and RMSE to be
around $\sim$ 0.9 $kcal/(mol\,\AA)$ and $\sim$ 1.2 $kcal/(mol\,\mbox{\AA})$
respectively. Also, based on this analysis, one could consider introducing
atom types based on their chemical environment rather than purely by its
nuclear charge, akin to methods employed for empirical force fields with
chemically diverse systems Ewig  _et al._ (2001).

Figure 5: Mean average error density of the force predictions of a MACE model
trained on 1000 points of a stachyose dataset. The densities are done
separately for each atom type. Colors correspond to elements (Hydrogen: white,
Oxygen: red, Carbon: grey).

All the prior discussions can be repeated with a separate model, e.g., a MACE
model created with the same training set. In Fig. 5, one can see that a very
similar trend can be observed, with even more distinguishable peaks in the
error distributions. Also, all the observations made in the 3D visualizer for
the MACE model are qualitatively identical to those made for the Nequip model
showing significant similarity in both models’ performance for the considered
system.

###  IV.2 Docosahexaenoic acid

A second example illustrating the use-case of the FFAST is that of the
docosahexaenoic acid (DHA) molecule, a fatty acid with the chemical formula
C22H32O2 for a total of 56 atoms. The molecule comprises a 22-carbon chain
with six double bonds and a carboxylic head. DHA is rather flexible due to its
long hydrocarbon tail, meaning the molecule can visit various extended and
folded states at ambient conditions. Therefore, understanding how well the
folding and unfolding processes are represented in the reference dataset is
crucial for constructing a reliable MLFF. FFAST provides a simple tool to
visualize the folding/unfolding processes by tracking the molecule’s radius of
gyration (or gyradius), as shown in Fig. 6. This does not require any ML model
to be loaded. Fig. 6 shows that, in total, the dataset consists of six compact
and six extended states. Moreover, the radius of gyration correlates with the
molecule’s potential energy after an equilibration period of approximately 10k
simulation steps. One region where the correlation is broken is the fourth
folded state (around 40k step). While the molecule is compact, its potential
energy is relatively high for those configurations, suggesting the occurrence
of some non-trivial interaction patterns demanding further investigation.
Interestingly, both ML models trained in this subsection would demonstrate the
largest energy MAEs and noticeable force MAEs strictly for this set of
structures (see Fig. 8), fortifying the previous conclusion.

Figure 6: Gyradius (white) and energy (orange) throughout the trajectory of a
DHA molecule. Values are averaged over a window of 2500 points to smooth out
noise.

Another essential analysis that can quickly be performed using FFAST and is a
prerequisite for training reliable MLFF for flexible molecules: ensuring the
quality of the training set. Here, we compare the distributions of gyration
radius, forces, and energies within the given (complete) dataset and the
generated training (sub)set of $1000$ points, see Figure 7. The comparison
verifies that the training set is indeed representative of the dataset and the
obtained ML model should not unexpectedly enter interpolation regimes while
within the reference’s configurational space.

Figure 7: Force distribution (left), energy distribution (middle) and gyradius
distribution (right) of the entire DHA dataset (white) against a training set
of $1000$ points (cyan). The training set was generated using the sGDML
training point selection method.

With such an insight into the content of the given reference and training
datasets, we can build MLFFs aiming at the description of folded and unfolded
DHA geometries and the transitions between them. Below, we will compare the
performance of two competing state-of-the-art ML architectures under the same
conditions. In Fig. 8, one can see Nequip (orange) and MACE (red), both
generated using the same training and validation set of 1000 points each. One
can observe that in this case, the MACE model has better general accuracy than
the Nequip equivalent, with each presenting an overall MAE on forces of 0.20
$kcal/(mol\,\AA)$ and 0.33 $kcal/(mol\,\mbox{\AA})$, respectively.
Furthermore, the error distribution peaks for the MACE model are narrower for
both energy and forces. Interestingly, the difference in accuracy does not
come from a better prediction of forces acting on a specific type of atoms
within the MACE model. A similar 50 % decrease in the MAEs is observed for
hydrogens, oxygens, and carbons, as can be seen by comparing atomic errors
(similar to the analyses shown in Fig. 3). This hints at the overall better
performance of the MACE architecture for reconstructing the PES of the DHA
molecule.

Figure 8: Basic error screen including timelines and distributions for both
energy and forces mean average errors across an entire stachyose dataset. In
orange and red respectively, Nequip and MACE models trained on 1000 points of
the dataset are used.

Notably, distinct loosely periodic peaks can be observed both in the energy
and force MAE timeline. These peaks match between Nequip and MACE, meaning
that they are unlikely to be due to an artifact of the models’ predictions but
rather a fundamental change in geometry throughout those time steps. Zooming
in on the valleys/peaks and displaying the configurations in the 3D visualizer
reveals that the low error regions correspond to extended geometries while the
high error peaks contain folded geometries. An example for each can be seen in
Fig. 9

Figure 9: (Left) Example of an extended DHA molecule, as found in the valleys
of the MAE timeline and low-error clusters. (Right) Example of a folded DHA
molecule, as found in the peaks of the MAE timeline and high-error clusters.
The geometries were chosen from the lowest and highest force prediction error
cluster, respectively.

The mean force prediction cluster errors (see Fig. 10) show a relatively flat
profile, i.e., the highest average force prediction error on a group of
similar geometries is less than two times higher than that of the lowest.
Nevertheless, similarly to the error timeline plots, visualizing the high
error clusters consistently displays folded configurations and vice-versa. One
can find an example of the lowest and highest error clusters in Fig. 9.

Figure 10: Mean average error of force predictions for 50 different clusters
of a stachyose dataset. The models used are a Nequip model (red) and a MACE
model (orange) trained on 1000 points. The clusters are ordered in ascending
prediction error of the MACE model.

Finally, additional useful information can be obtained by visualizing the
average force prediction error on every atom in the 3D visualizer, as seen in
Fig. 11. Once again, hydrogens are very well predicted across the board,
suggesting that a significant contribution to the low overall prediction error
comes from the 32 hydrogen atoms in the molecule. In contrast, one can see
that the carbon atoms inside a carboxylic group and those closest to it (on
the left side of the figure) are significantly worse predicted than the rest
of the chain. While not shown here, the Nequip model presents a similar
overall trend. One can easily explain the observed nonuniform force
reconstruction for carbon atoms (both for Nequip and MACE) by combining the
information obtained from the FFAST software and our physical/ML intuition.
The carbon atoms near the head of the DHA molecule containing oxygen atoms are
exposed to a more complex chemical composition of their neighborhood than
those in other parts of the molecule. Thus, the descriptor space for such
carbon atoms contains more possible states, especially in folded-like
configurations. This leads to noticeably more challenging training tasks for
the ML models resulting in larger prediction errors. This again hints at the
need to reconsider the definition of atomic types based on the atom’s
environment composition rather than its nuclear weight when we move from
simple periodic systems or small molecules to more challenging systems.

Figure 11: Extended configuration of a DHA molecule. Atoms are colored
according to their respective mean average force error as predicted by a MACE
model.

All in all, our analysis shows that beyond the insignificant for most
practical applications error difference in less than 0.25 $kcal/mol$ for the
energy and 0.20 $kcal/(mol\,\mbox{\AA})$ for the forces, both the Nequip and
MACE models demonstrate equivalent results for the given test problem.
Therefore, one can in principle apply any of these two models to study
molecules similar to DHA or stachyose examples.

##  V Outlook

The applications in this manuscript (as well as many other examples omitted)
demonstrate that increased chemical and structural complexity of mid and
large-size molecules or periodic systems require a detailed analysis of MLFFs’
stability, reliability, and performance. For such cases, standard metrics such
as MAE and RMSE are only one piece of the puzzle. High overall accuracies for
reconstructing a PES or FF can largely stem from the good performance of an ML
model on specific (albeit large) parts of the system; see the case of
hydrogens in the examples above. At the same time, forces for other vital
system components might be predicted at a significantly lower accuracy. For
both the examples of DHA and stachyose, the largest errors for carbon atoms
lie in a region where the overall force error distribution shows negligible
probability. This effect is clearly observed despite the molecules containing
less than a hundred atoms. For more complex systems, such as molecules on
surfaces, periodic structures containing rigid crystal frames and molecular
interstitials, etc., the difference between mean errors and region-specific
errors could become even more prominent. When looking to enhance models and go
beyond current capabilities, all these effects can heavily influence the next
directions to undertake. To create a reliable way to systematically improve ML
models, one first needs to agree on exactly what the current limitations are
and which areas need targeted improvements. Furthermore, reaching certain
limitations does not necessarily mean the employed ML architectures cannot
handle the given tasks. With a detailed insight into the models’ performance,
one can often significantly increase its applicability by adjusting the
training process or the underlying data. Therefore, FFAST or similar software
packages should become a staple tool for ML model developers and users.

Another challenge where FFAST can be indispensable is assessing the stability
of MLFFs MD trajectories. Providing stable long-time dynamics for complex
systems is not at all guaranteed, even for state-of-the-art MLFFs. A long
enough simulation eventually reaches a region where extrapolation is
necessary, which is likely to generate unphysical forces and heavily influence
the rest of the dynamics. These effects can sometimes be easily detected (e.g.
breaking of bonds, dissociated hydrogen atoms, etc), but less obvious
consequences can also follow. Namely, unphysical behavior can occur without
showing apparent damage to the molecular structure within the MD run.
Physical/chemical property computed thereafter would be affected without clear
evidence to the user that extrapolation regions were reached. While some
techniques, such as predicting the model uncertainty, can partially alleviate
this problem, the fact that the given ML model is confident in its prediction
does not necessarily mean that the prediction is correct. The outlier
detection tools and interactive visualization implemented in FFAST provide an
efficient alternative instrument to detect such situations without significant
human effort. Considering that MLFFs by construction focus first on the most
typical interaction patterns, such analyses also highlight non-trivial physics
and chemistry that might be present in the system. Therefore, combining FFAST
with modern MLFF architectures should also enable a better understanding of
the interplay between different types of interactions for large and complex
systems.

##  VI Conclusion

With increased complexities of systems that can be simulated using novel MLFFs
comes a need to provide insightful analysis for these models. Even the most
advanced MLFFs demonstrate highly heterogeneous predictions across
configurational space (CS), invisible to overall error metrics. Thus, the
detailed breakdowns of the ML models’ performance on different domains of CS
and various parts of the system under study should become standard practice.
FFAST is an actively developing tool allowing experts and non-experts to get
an in-depth insight into the performance of MLFFs or any other force field of
choice on any dataset. In this paper, a list of provided features was given
and explained, including the example of a best practice scenario. The features
were showcased using two datasets of medium-sized flexible organic molecules:
docosahexaenoic acid and stachyose.

Two state-of-the-art ML models (Nequip and MACE) were trained for each
dataset. With those, a general FFAST workflow was showcased and the following
conclusions followed. Hydrogens are generally significantly better predicted
for both molecules than other atoms, considerably lowering the overall ML
models’ MAEs and RMSEs in hydrogen-rich systems. In contrast, forces acting on
carbons and oxygens are unevenly reconstructed. Atoms with chemically or
structurally more diverse environments present significantly higher prediction
errors compared to the other atoms of the same type. For instance, the
prediction errors on DHA rise as the molecule folds, and the primary source of
the errors comes from the carboxylic group.

Beyond a comprehensive analysis of the stability and reliability of ML models,
FFAST provides the means to analyze reference data. For instance, one can
easily find the number of folding/unfolding processes in a given MD trajectory
and the correlation of that with the molecule’s potential energy. Another
possible application is comparing energy/force distributions in different
datasets (training, validation, and test), which should precede any ML model
training. As FFAST offers detailed reference datasets and MLFF error analysis,
field-specialized knowledge can be applied to enhance MLFFs after identifying
the remaining challenges. The software is provided as an open-source project
on github.com/fonsecag/FFAST and is designed to be modulable and expandable to
adapt to the user’s needs.

###### Acknowledgements.

We acknowledge financial support from the Luxembourg National Research (FNR)
under the AFR project 14593813, FNR DTU-PRIDE MASSENA, and the European
Research Council (ERC-AdG grant FITMOL).

##  VII Bibliography

## References

  * Lee  _et al._ (2019) A. C. Lee, J. L. Harris, K. K. Khanna, and J.-H. Hong, “A Comprehensive Review on Current Advances in Peptide Drug Development and Design,” Int. J. Mol. Sci. 20, 2383 (2019).
  * Ramos  _et al._ (2019) J. Ramos, J. Muthukumaran, F. Freire, J. Paquete-Ferreira, A. R. Otrelo-Cardoso, D. Svergun, A. Panjkovich, and T. Santos-Silva, “Shedding Light on the Interaction of Human Anti-Apoptotic Bcl-2 Protein with Ligands through Biophysical and in Silico Studies,” Int. J. Mol. Sci. 20, 860 (2019).
  * Mansbach  _et al._ (2019) R. A. Mansbach, T. Travers, B. H. McMahon, J. M. Fair, and S. Gnanakaran, “Snails In Silico: A Review of Computational Studies on the Conopeptides,” Mar. Drugs 17, 145 (2019).
  * Jing  _et al._ (2019) Z. Jing, C. Liu, S. Y. Cheng, R. Qi, B. D. Walker, J.-P. Piquemal, and P. Ren, “Polarizable Force Fields for Biomolecular Simulations: Recent Advances and Applications,” Annu. Rev. Biophys. 48, 371–394 (2019).
  * Ekins (2016) S. Ekins, “The Next Era: Deep Learning in Pharmaceutical Research,” Pharm. Res. 33, 2594–2603 (2016).
  * Smalley (2017) E. Smalley, “AI-powered drug discovery captures pharma interest,” Nat. Biotechnol. 35, 604–605 (2017).
  * Elton  _et al._ (2018) D. C. Elton, Z. Boukouvalas, M. S. Butrico, M. D. Fuge, and P. W. Chung, “Applying machine learning techniques to predict the properties of energetic materials,” Sci. Rep. 8, 9059 (2018).
  * Elton  _et al._ (2019) D. C. Elton, Z. Boukouvalas, M. D. Fuge, and P. W. Chung, “Deep learning for molecular design—a review of the state of the art,” Mol. Syst. Des. Eng. 4, 828–849 (2019).
  * Faber  _et al._ (2016) F. A. Faber, A. Lindmaa, O. A. von Lilienfeld, and R. Armiento, “Machine Learning Energies of 2 Million Elpasolite $(AB{}C{}_{2}{}D{}_{6})$ Crystals,” Phys. Rev. Lett. 117, 135502 (2016).
  * Naserifar  _et al._ (2021) S. Naserifar, Y. Chen, S. Kwon, H. Xiao, and W. A. Goddard, III, “Artificial Intelligence and QM/MM with a Polarizable Reactive Force Field for Next-Generation Electrocatalysts,” Matter 4, 195–216 (2021).
  * Cova and Pais (2019) T. F. G. G. Cova and A. A. C. C. Pais, “Deep Learning for Deep Chemistry: Optimizing the Prediction of Chemical Patterns,” Front. Chem. 7 (2019).
  * Batzner  _et al._ (2022) S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, “E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials,” Nat. Commun. 13, 2453 (2022).
  * Geiger and Smidt (2022) M. Geiger and T. Smidt, “E3nn: Euclidean neural networks,” (2022), arxiv:2207.09453 [cs.LG] .
  * Batatia  _et al._ (2022a) I. Batatia, D. P. Kovacs, G. N. C. Simm, C. Ortner, and G. Csanyi, “MACE: Higher order equivariant message passing neural networks for fast and accurate force fields,” in  _NeurIPS_ , edited by A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho (2022).
  * Batatia  _et al._ (2022b) I. Batatia, S. Batzner, D. P. Kovács, A. Musaelian, G. N. C. Simm, R. Drautz, C. Ortner, B. Kozinsky, and G. Csányi, “The design space of E(3)-equivariant atom-centered interatomic potentials,” (2022b), arxiv:2205.06643 [stat.ML] .
  * Bartók  _et al._ (2010) A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, “Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons,” Phys. Rev. Lett. 104, 136403 (2010).
  * Bartók, Kondor, and Csányi (2013) A. P. Bartók, R. Kondor, and G. Csányi, “On representing chemical environments,” Phys. Rev. B 87, 184115 (2013).
  * Behler (2016) J. Behler, “Perspective: Machine learning potentials for atomistic simulations,” J. Chem. Phys. 145, 170901 (2016).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, “Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces,” Phys. Rev. Lett. 98, 146401 (2007).
  * Schütt  _et al._ (2017a) K. T. Schütt, P.-J. Kindermans, H. E. Sauceda, S. Chmiela, A. Tkatchenko, and K.-R. Müller, “SchNet: A continuous-filter convolutional neural network for modeling quantum interactions,” (2017a), arxiv:1706.08566 [stat.ML] .
  * Schütt  _et al._ (2017b) K. T. Schütt, F. Arbabzadah, S. Chmiela, K. R. Müller, and A. Tkatchenko, “Quantum-chemical insights from deep tensor neural networks,” Nat. Commun. 8, 13890 (2017b).
  * Schütt  _et al._ (2018) K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, and K.-R. Müller, “SchNet – A deep learning architecture for molecules and materials,” J. Chem. Phys. 148, 241722 (2018).
  * Schütt  _et al._ (2019) K. T. Schütt, P. Kessel, M. Gastegger, K. A. Nicoli, A. Tkatchenko, and K.-R. Müller, “SchNetPack: A Deep Learning Toolbox For Atomistic Systems,” J. Chem. Theory Comput. 15, 448–455 (2019).
  * Rupp  _et al._ (2012) M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, “Fast and accurate modeling of molecular atomization energies with machine learning,” Phys. Rev. Lett. 108, 058301 (2012).
  * Bartók  _et al._ (2017) A. P. Bartók, S. De, C. Poelking, N. Bernstein, J. R. Kermode, G. Csányi, and M. Ceriotti, “Machine learning unifies the modeling of materials and molecules,” Sci. Adv. 3, 1701816 (2017).
  * Unke and Meuwly (2019) O. T. Unke and M. Meuwly, “PhysNet: A Neural Network for Predicting Energies, Forces, Dipole Moments, and Partial Charges,” J. Chem. Theory Comput. 15, 3678–3693 (2019).
  * Unke  _et al._ (2021a) O. T. Unke, S. Chmiela, M. Gastegger, K. T. Schütt, H. E. Sauceda, and K.-R. Müller, “SpookyNet: Learning force fields with electronic degrees of freedom and nonlocal effects,” Nat. Commun. 12, 7273 (2021a).
  * Unke  _et al._ (2021b) O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko, and K.-R. Müller, “Machine Learning Force Fields,” Chem. Rev. 121, 10142–10186 (2021b).
  * Chmiela  _et al._ (2018) S. Chmiela, H. E. Sauceda, K.-R. Müller, and A. Tkatchenko, “Towards exact molecular dynamics simulations with machine-learned force fields,” Nat. Commun. 9, 3887 (2018).
  * Chmiela  _et al._ (2019) S. Chmiela, H. E. Sauceda, I. Poltavsky, K.-R. Müller, and A. Tkatchenko, “sGDML: Constructing accurate and data efficient molecular force fields using machine learning,” Comput. Phys. Commun. 240, 38–45 (2019).
  * Sauceda  _et al._ (2022) H. E. Sauceda, L. E. Gálvez-González, S. Chmiela, L. O. Paz-Borbón, K.-R. Müller, and A. Tkatchenko, “BIGDML—Towards accurate quantum machine learning force fields for materials,” Nat. Commun. 13, 3733 (2022).
  * Qiao  _et al._ (2020) Z. Qiao, M. Welborn, A. Anandkumar, F. R. Manby, and T. F. Miller, III, “OrbNet: Deep learning for quantum chemistry using symmetry-adapted atomic-orbital features,” J. Chem. Phys. 153, 124111 (2020).
  * Smith, Isayev, and Roitberg (2017) J. S. Smith, O. Isayev, and A. E. Roitberg, “ANI-1: An extensible neural network potential with DFT accuracy at force field computational cost,” Chem. Sci. 8, 3192–3203 (2017).
  * Devereux  _et al._ (2020) C. Devereux, J. S. Smith, K. K. Huddleston, K. Barros, R. Zubatyuk, O. Isayev, and A. E. Roitberg, “Extending the Applicability of the ANI Deep Learning Molecular Potential to Sulfur and Halogens,” J. Chem. Theory Comput. 16, 4192–4202 (2020).
  * Deringer  _et al._ (2018) V. L. Deringer, N. Bernstein, A. P. Bartók, M. J. Cliffe, R. N. Kerber, L. E. Marbella, C. P. Grey, S. R. Elliott, and G. Csányi, “Realistic Atomistic Structure of Amorphous Silicon from Machine-Learning-Driven Molecular Dynamics,” J. Phys. Chem. Lett. 9, 2879–2885 (2018).
  * Ryczko  _et al._ (2018) K. Ryczko, K. Mills, I. Luchak, C. Homenick, and I. Tamblyn, “Convolutional neural networks for atomistic systems,” Comput. Mater. Sci. 149, 134–142 (2018).
  * Frank, Unke, and Müller (2023) J. T. Frank, O. T. Unke, and K.-R. Müller, “So3krates: Equivariant attention for interactions on arbitrary length-scales in molecular systems,” (2023), arxiv:2205.14276 [cs.LG] .
  * Christensen  _et al._ (2020) A. S. Christensen, L. A. Bratholm, F. A. Faber, and O. Anatole von Lilienfeld, “FCHL revisited: Faster and more accurate quantum machine learning,” J. Chem. Phys. 152, 044107 (2020).
  * Wang  _et al._ (2019) J. Wang, S. Olsson, C. Wehmeyer, A. Pérez, N. E. Charron, G. de Fabritiis, F. Noé, and C. Clementi, “Machine Learning of Coarse-Grained Molecular Dynamics Force Fields,” ACS Cent. Sci. 5, 755–767 (2019).
  * Ko  _et al._ (2021a) T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, “General-Purpose Machine Learning Potentials Capturing Nonlocal Charge Transfer,” Acc. Chem. Res. 54, 808–817 (2021a).
  * Kocer, Ko, and Behler (2022) E. Kocer, T. W. Ko, and J. Behler, “Neural Network Potentials: A Concise Overview of Methods,” Annu. Rev. Phys. Chem. 73, 163–186 (2022).
  * Gasteiger, Groß, and Günnemann (2022) J. Gasteiger, J. Groß, and S. Günnemann, “Directional message passing for molecular graphs,” (2022), arxiv:2003.03123 [cs.LG] .
  * Gasteiger, Becker, and Günnemann (2022) J. Gasteiger, F. Becker, and S. Günnemann, “GemNet: Universal directional graph neural networks for molecules,” (2022), arxiv:2106.08903 [physics.comp-ph] .
  * Deringer  _et al._ (2021) V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti, and G. Csányi, “Gaussian Process Regression for Materials and Molecules,” Chem. Rev. 121, 10073–10141 (2021).
  * Bartók  _et al._ (2018) A. P. Bartók, J. Kermode, N. Bernstein, and G. Csányi, “Machine Learning a General-Purpose Interatomic Potential for Silicon,” Phys. Rev. X 8, 041048 (2018).
  * Thomas  _et al._ (2018) N. Thomas, T. Smidt, S. Kearnes, L. Yang, L. Li, K. Kohlhoff, and P. Riley, “Tensor field networks: Rotation- and translation-equivariant neural networks for 3D point clouds,” (2018), arxiv:1802.08219 [cs.LG] .
  * Khorshidi and Peterson (2016) A. Khorshidi and A. A. Peterson, “Amp: A modular approach to machine learning in atomistic simulations,” Comput. Phys. Commun. 207, 310–324 (2016).
  * Li, Kermode, and De Vita (2015) Z. Li, J. R. Kermode, and A. De Vita, “Molecular Dynamics with On-the-Fly Machine Learning of Quantum-Mechanical Forces,” Phys. Rev. Lett. 114, 096405 (2015).
  * Keith  _et al._ (2021) J. A. Keith, V. Vassilev-Galindo, B. Cheng, S. Chmiela, M. Gastegger, K.-R. Müller, and A. Tkatchenko, “Combining Machine Learning and Computational Chemistry for Predictive Insights Into Chemical Systems,” Chem. Rev. 121, 9816–9872 (2021).
  * Ko  _et al._ (2021b) T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, “A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer,” Nat. Commun. 12, 398 (2021b).
  * Yao  _et al._ (2018) K. Yao, J. E. Herr, D. W. Toth, R. Mckintyre, and J. Parkhill, “The TensorMol-0.1 model chemistry: A neural network augmented with long-range physics,” Chem. Sci. 9, 2261–2269 (2018).
  * Grisafi and Ceriotti (2019) A. Grisafi and M. Ceriotti, “Incorporating long-range physics in atomic-scale machine learning,” J. Chem. Phys. 151, 204105 (2019).
  * Hoja  _et al._ (2021) J. Hoja, L. Medrano Sandonas, B. G. Ernst, A. Vazquez-Mayagoitia, R. A. DiStasio Jr., and A. Tkatchenko, “QM7-X, a comprehensive dataset of quantum-mechanical properties spanning the chemical space of small organic molecules,” Sci. Data 8, 43 (2021).
  * Blum and Reymond (2009) L. C. Blum and J.-L. Reymond, “970 Million Druglike Small Molecules for Virtual Screening in the Chemical Universe Database GDB-13,” J. Am. Chem. Soc. 131, 8732–8733 (2009).
  * Ramakrishnan  _et al._ (2014) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, “Quantum chemistry structures and properties of 134 kilo molecules,” Sci. Data 1, 140022 (2014).
  * Ruddigkeit  _et al._ (2012) L. Ruddigkeit, R. van Deursen, L. C. Blum, and J.-L. Reymond, “Enumeration of 166 Billion Organic Small Molecules in the Chemical Universe Database GDB-17,” J. Chem. Inf. Model. 52, 2864–2875 (2012).
  * Chmiela  _et al._ (2017) S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, “Machine learning of accurate energy-conserving molecular force fields,” Sci. Adv. 3, 1603015 (2017).
  * Chmiela  _et al._ (2023) S. Chmiela, V. Vassilev-Galindo, O. T. Unke, A. Kabylda, H. E. Sauceda, A. Tkatchenko, and K.-R. Müller, “Accurate global machine learning force fields for molecules with hundreds of atoms,” Sci. Adv. 9, adf0873 (2023).
  * Smith  _et al._ (2020) J. S. Smith, R. Zubatyuk, B. Nebgen, N. Lubbers, K. Barros, A. E. Roitberg, O. Isayev, and S. Tretiak, “The ANI-1ccx and ANI-1x data sets, coupled-cluster and density functional theory properties for molecules,” Sci. Data 7, 134 (2020).
  * Fonseca  _et al._ (2021) G. Fonseca, I. Poltavsky, V. Vassilev-Galindo, and A. Tkatchenko, “Improving molecular force fields across configurational space by combining supervised and unsupervised machine learning,” J. Chem. Phys. 154, 124102 (2021).
  * Fu  _et al._ (2022) X. Fu, Z. Wu, W. Wang, T. Xie, S. Keten, R. Gomez-Bombarelli, and T. Jaakkola, “Forces are not enough: Benchmark and critical evaluation for machine learning force fields with molecular simulations,” (2022), arxiv:2210.07237 [physics.comp-ph] .
  * Ewig  _et al._ (2001) C. S. Ewig, R. Berry, U. Dinur, J.-R. Hill, M.-J. Hwang, H. Li, C. Liang, J. Maple, Z. Peng, T. P. Stockfisch, T. S. Thacher, L. Yan, X. Ni, and A. T. Hagler, “Derivation of class II force fields. VIII. Derivation of a general quantum mechanical force field for organic compounds,” J. Comput. Chem. 22, 1782–1800 (2001).