Machine learning potentials for complex aqueous systems made simple

The implementation of our MLP framework relies on committee neural network potentials (C-NNPs) Schran2020/10.1063/5.0016004, in our case built from Behler-Parrinello NNPs Behler2007/10.1103/PhysRevLett.98.146401; Behler2017/10.1002/anie.201703114. While these models only consider the local environment of each atom explicitly within a finite range, long-range contributions can in principle be incorporated in MLPs as demonstrated in the literature Artrith2011/10.1103/PhysRevB.83.153101; Grisafi2019/10.1063/1.5128375; Ko2021/10.1038/s41467-020-20427-2; Yue2021/10.1063/5.0031215. The main idea behind the current approach is the combination of multiple NNPs in a “committee model”, where the committee members are separately trained from independent random initial conditions to a subset of the total training set. The resulting committee model has multiple benefits over its individual members. While the committee average, which is used as the prediction of the whole model, has been shown to provide better performance than the individual NNPs, the committee disagreement, which is defined as the standard deviation between the committee member predictions, grants access to an estimate of the error of the model. If scaled by a constant obtained by comparison to the true validation error Imbalzano2021/10.1063/5.0036522, the committee disagreement provides an objective measure of the accuracy of the underlying model. To construct a training set of such a model in an automated and data-driven way, new configurations with the highest disagreement can be added to the training set. This is an active learning strategy called query by committee and can be used to systematically improve a machine learning model, while making efficient use of the limited data available which provides an important advantage for the current application. Further details on the C-NNP methodology can be found in Ref. Schran2020/10.1063/5.0016004.

These concepts can be used for the rapid development of MLPs as described in the following and schematically depicted in the top panel of Fig. 1. Initially, a small-scale reference AIMD simulation is performed to sample the system of interest under the selected thermodynamic condition. Small-scale refers here mostly to the envisaged time and length scales of the final application, while they might still be considered large from the ab initio perspective. In practice, AIMD trajectories with a length of 30 ps have been sufficient for this purpose with system sizes that are large enough to cover all required local chemical environments. Given this trajectory, which can also come from existing previous work, we then construct the training set of the C-NNP model by selecting the most representative configurations for the system and condition of interest in an iterative active learning procedure. In the beginning, a small set of 20 structures is randomly selected from the reference trajectory in order to train an initial C-NNP model. Next, query by committee is used to actively select 20 configurations based on the highest committee disagreement in the atomic forces from the set of candidate structures provided by the reference trajectory. These points are added to the training set and an improved C-NNP model is trained to the extended training set. No additional ab initio reference calculations are required in this process, as the potential energy and atomic forces are already available from the reference AIMD simulation. Such iterations are repeated until convergence of the committee disagreement is observed, indicated by marginal improvements of the disagreement in subsequent iterations and no substantial difference in disagreement between the selected points and those already present in the training set. This implies that a sufficient variety of structures has been added to the training set to yield an accurate and robust C-NNP model.

The uncomplicated construction of MLPs for new systems requires a general set of atomic descriptors. In the case of the atom-centered symmetry functions Behler2011/10.1063/1.3553717, employed here, we make that possible by generating a systematic set of radial and angular functions. This set consists of ten equidistantly shifted radial functions with a fixed width and four angular functions applied to each pair and triple of atoms, respectively. In addition, we apply the same hyperparameters, such as number of committee members, hidden layers and nodes, as well as neural network optimization parameters, to every system to remove as much user input from our procedure as possible. Thanks to the active learning procedure that adapts the training set to the flexibility of a particular model, we have in practice observed no limitations of such an application of a general set of parameters and settings. The remaining task of the user is to provide the reference trajectory for the system of interest under the chosen simulation conditions. Given that input, a C-NNP model can be obtained in practice without further adjustments and in a short amount of time.

Once a new model is trained, it can be applied to the system of interest close to the same state point sampled by the original reference trajectory, reducing the computational cost compared to the ab initio reference by orders of magnitude. During such simulations, the committee disagreement of the C-NNP model is a valuable tool to gauge the validity of the model. It provides an intrinsic error estimate and thus can be monitored over the course of the simulation and compared to the training process. These concepts enable the straightforward extension of both time and length scales in molecular simulations. Per design, the construction in the above fashion will be state dependent. The heart of the matter is that this lack of generality is compensated by the speed and simplicity in its fitting.

Besides this simple and robust framework for the generation of MLPs for complex systems, special emphasis lies on the selection of a suitable electronic structure reference method. MLPs will always only be as good as their underlying reference method and a user has to make a careful choice for each system of interest. In this context, density functional theory (DFT) has become indispensable for the investigation of aqueous systems, while careful benchmark studies Grossman2003/10.1063/1.1630560; Distasio2014/10.1063/1.4893377; Forster-Tonigold2014/10.1063/1.4892400; Gillan2016/10.1063/1.4944633; Chen2017/10.1073/pnas.1712499114; Marsalek2017/10.1021/acs.jpclett.7b00391; Brandenburg2019/10.1063/1.5121370; Schienbein2020/10.1002/anie.202009640 can guide the selection of suitable functionals. In addition, there has been promising steps to better understand remaining limitations of existing functionals and provide potential solutions in recent studies Sharkas2020/10.1073/pnas.1921258117; Wagle2021/10.1063/5.0041620; Duignan2021/10.1021/acs.accounts.1c00107. Combined with the inexpensive yet reliable representation of interactions, as shown in this work, this opens up the possibility for the uncomplicated but accurate investigation of many technologically relevant systems.

Following our active learning workflow, we have obtained C-NNP models for six different aqueous phase systems, which are shown in the bottom panel of Fig. 1. These diverse systems involve various types of interactions, feature bulk, interface and confinement regions and go up to a chemical composition of four elements. In all cases, the training set of the model is exclusively based on a set of structures generated by ab initio molecular dynamics simulations, as described in detail in the Materials and Methods section. Given that input, all six models have been generated without further adjustments. For all these systems, convergence was achieved with a compact training set of roughly 300 structures, highlighting the advantages of the active learning procedure.

An automated training procedure calls for an efficient and robust validation protocol. Through extensive comparisons of our models to the underlying ab initio reference trajectories we have identified a general set of properties that serve to provide a stern test of our models. These can be evaluated for any system of interest and provide a broad overview of the performance of a MLP, while further tests are included in the Supporting Information. The selected properties exemplify the performance of the models for structural and dynamical properties as well as the precision of the force prediction. Specifically, the performance for structural properties is assessed by the match of the radial distribution functions (RDF) for all involved species comparing the ab initio reference to the model prediction based on molecular dynamics simulations. All RDFs of a given system provide a comprehensive summary of the structural arrangement of the system of interest and are thus ideal to evaluate the performance of the machine learning model for thermodynamic properties. Dynamical properties are validated by comparing the species-resolved vibrational density of states (VDOS) obtained with the model and the reference, which gives a comprehensive overview of inter- and intramolecular motions. Finally, the force prediction of the model is validated by the force root mean square error (RMSE) of a randomly selected subset of structures from the ab initio reference simulation. This quantity is chosen since the forces ultimately drive the molecular dynamics simulations when using the model. In order to make these properties comparable for all systems, they are reduced into a score by suitable difference measurements and subsequent averaging over the involved species as described in detail in the Supporting Information. The entire testing protocol functions in an automated manner and efficiently provides a condensed summary of the accuracy of each model.

The resulting summary of the quality assessment for all six studied systems is shown in Fig. 2. From this analysis it is clear that all models reproduce the three selected properties with rather high precision, where the RDF score ranges between 100 to 98%, the VDOS score between 98 to 96%, and the force score between 95 to 86%. To illustrate the meaning of these values, we depict the individual functions (RDF and VDOS) and the force correlation for the solvated fluoride ion C-NNP model along with the total scores in Fig. 2, while all other properties for the remaining models are compiled in the Supporting Information. This comparison shows that the selected properties are reproduced with good agreement to the ab initio reference method by our six C-NNP models. In addition, all C-NNP results included in this performance summary are based on substantially extended simulation times compared to the AIMD references as described in detail in the Supporting Information. This highlights the robust nature of our models, enabling reliable predictions over long time scales. We note that statistical fluctuations due to the more converged nature of our C-NNP simulations accounts only for very minor changes of the final property scores on the order of 0.5%.

Besides the three sets of properties that we have quantitatively validated here, we have performed additional performance tests, in particular for the more complex systems of water confined in nanotubes and between \chMoS2 sheets as well as water on \chTiO2. These tests include the detailed analysis of the global structure of the solid and liquid subsystems, as encoded by the density profiles, the hydrogen bonding of water in the various systems, and the orientation of water with respect to the involved interfaces. For all these tests, which are presented in detail in the Supporting Information, we observe good agreement between our C-NNP results and the AIMD reference simulations within the statistics of those shorter AIMD runs. We are therefore confident that our performance overview, presented in Fig. 2, underlines the high quality of our C-NNP models.

The quality assessment for the six different systems included in this work clearly highlights that our simple and straightforward process to develop machine learning potentials is able to provide robust and accurate models for the selected thermodynamic condition. Compared to the typical DFT setups employed here, the evaluation of the potential energy and atomic forces is usually four to five orders of magnitude faster with the C-NNP model. As a consequence, all chosen systems could now be studied in detail using exhaustive simulations that are accessible with the developed models. Given the focus on general properties in the testing protocol, we expect that it could prove useful for the development of potentials for various other solid-liquid systems of technological and/or scientific interest.

Let us finally showcase the potential of the presented methodology to extend the length and time scales of molecular simulations and thus further the understanding of a system of interest. For that purpose we investigate structural and dynamical properties of water in contact with rutile \chTiO2(110). This system is of scientific and technological importance due to the application of \chTiO2 for example in photocatalysis or self-cleaning coatings and sensors. In addition, it is an established prototypical oxide system in surface science Pang2013/10.1021/cr300409r and a rather controversial benchmark system both for theory and experiment Diebold2017/10.1063/1.4996116. For example, the extent of the mobility of water in the contact layers, relevant e.g. for a detailed understanding of catalytic processes, has been the focus of substantial research Predota2007/10.1021/jp065165u; Liu2010/10.1103/PhysRevB.82.161415; Spencer2009/10.1021/jp8109918; English2012/10.1080/00268976.2012.683888; Agosta2017/10.1063/1.4991381.

In order to shed light on these questions, we have used the developed C-NNP model to simulate rutile \chTiO2(110) in contact with water. The model was constructed from a 30 ps AIMD simulation at 300 K with the optB88-vdW functional Klimes2010/10.1088/0953-8984/22/2/022201, involving a four O-Ti-O trilayer slab in contact with 80 water molecules, forming a 1.5 nm water film on the surface. After the development and benchmarking of the C-NNP model as shown in the previous section, we made a 2$\times$2 model of the interface (resulting in a \chTiO2/H2O setup with 1728 atoms) and ran 5 ns of MD. Reaching such length and time scales with AIMD simulations would represent an enormous computational burden, while they can be routinely performed with the C-NNP models. Further details of these simulations can be found in the Materials and Methods section.

First, we analyze the water structuring by looking at the density profile of water on the \chTiO2 surface shown in panel b) of Fig. 3. The first important observation is the close match between the crude density profile obtained for our AIMD simulation and the statistically converged profile obtained with the significantly extended C-NNP simulation. Overall, we observe a highly structured arrangement of water in the first two layers, which correspond to the water adsorbed on the 5-coordinated titanium site for the first peak and the water around the 6-coordinated titanium site for the second peak, as shown in the snapshot in the top panel. This density profile is substantially more structured than for AIMD simulations of water on rutile Liu2010/10.1103/PhysRevB.82.161415 with the PBE functional, which highlights the complex dependence of interfacial properties on the chosen functional. Given the improved understanding of the dependence of the properties of water on the DFT functional, in particular regarding the inclusion of dispersion interactions Gillan2016/10.1063/1.4944633, we conclude that modern DFT approaches predict a highly structured arrangement of water on the rutile surface reaching up to about 1 nm into the liquid. This is also in good agreement with the density profiles obtained from MLP simulations of water on the anatase (101) surface Andrade2020/10.1039/c9sc05116c that used the SCAN functional as the reference method.

In a next step, we evaluate the water diffusion coefficient resolved by the distance from the \chTiO2 surface. Specifically, we make use of the mean square displacement, which we spatially decompose based on the position of each water molecule at zero delay Pluharova2019/10.1021/acs.jctc.8b00111 — an approach made feasible by the extensive statistics provided by our C-NNP model. We then obtain a local estimate of the diffusion coefficient by the well-known Einstein relation, which can be evaluated separately for the parallel (xy) and perpendicular (z) directions with respect to the interface. Panel c) in Fig. 3 depicts the resulting water diffusion constant D_xy and D_z as a function of the distance from the \chTiO2 surface. As anticipated from the very structured density profile, the mobility close to the surface is reduced substantially, where essentially no diffusion is observed in the strongly adsorbed contact layer. Beyond the first and second layer, the diffusivity in the xy-direction increases steadily up to the water-vacuum interface. At the same time, D_z features a plateau around 1 nm from the surface and a substantial increase in diffusion towards the vapor interface. Overall, this analysis highlights the strong influence of the \chTiO2 interface on the water diffusion stretching more than 1 nm into the liquid.

Next, we address the accuracy of our extended C-NNP simulations by analyzing the intrinsic error estimate of our model, given by the atomic force committee disagreement $\sigma_{\text{F}}$. Panel d) in Fig. 3 resolves this local error estimate of the atomic force components for all water atoms as a function of the distance from the surface. This error estimate is a product of the committee disagreement $\sigma_{\text{F}}$, as directly provided by our C-NNP simulations, and a scaling factor $\alpha$ determined to match the force RMSE of a validation set as proposed in Ref. Imbalzano2021/10.1063/5.0036522. In addition, we have evaluated this error estimate with our C-NNP model for all configurations of the original AIMD simulation to assess if the increased system size or C-NNP generated structures lead to higher errors. Overall, we observe error estimates between 40 to 80 meV/Å over the entire water region with only slightly higher values close to the \chTiO2 surface, indicating the increased complexity of the involved interactions in this inhomogeneous region. Furthermore, the error (averaged over all water atoms) does not deteriorate over the course of the 5 ns long trajectory, fluctuating around an average of 60 meV/Å, as shown in the inset of panel d). Such atomic force errors are similar or even smaller than those reported for other developed MLPs e.g. for pure water Morawietz2016/10.1073/pnas.1602375113; Cheng2019/10.1073/pnas.1815117116; Gartner2020/10.1073/pnas.2015440117; Schran2020/10.1063/5.0016004. At the same time, the error estimate obtained for the AIMD configurations features essentially the same distance resolved profile, which reveals that our C-NNP simulations are able to conserve their predictive power, while substantially extending both time and length scales of the simulations.

Finally, we analyze the free energy profile of water adsorbed in the first two contact layers on the \chTiO2 surface, as depicted in panels e) and f) in Fig. 3 for the AIMD and C-NNP simulation, respectively. From the direct comparison between the AIMD and C-NNP results it is clear that the limited statistics of the shorter AIMD simulation is insufficient to provide reliable insight into this property. Only with the extensive sampling enabled by the C-NNP model, the free energy profile can be fully converged. The C-NNP free energy profile clearly underlines the strong preference for water adsorption above the fivefold coordinated titanium sites in the first contact layer and the slightly weaker adsorption of water around the sixfold coordinated titanium and threefold coordinated oxygen sites in the second contact layer. In between these adsorption sites, substantial free energy barriers are observed, highlighting the immobile nature of the two contact layers as also revealed by the analysis of the water diffusion.

In summary, the extensive simulations with our C-NNP model, obtained in a straightforward and efficient workflow, provide detailed insight into the properties of water on the rutile surface. We observe a pronounced water layering effect with strong density fluctuations and clear evidence of a highly structured arrangement of water in the first adsorption layers. In addition, our analysis of the water dynamics reveals almost no water diffusion close to the interface and a strong influence on the diffusion stretching more than 1 nm into the liquid. The treatment of a complex interfacial system such as this one requires an accurate description of the binding at the various adsorption sites as well as long-time sampling of the dynamics. The C-NNP model developed here delivers on both fronts, which highlights the potential of our approach to deepen understanding of technologically relevant solid-liquid systems.

In this work, we have presented a machine learning framework that makes the generation of MLPs simple. We have also showcased its versatility for the description of a range of complex aqueous systems. Making use of committee neural network potentials, we have shown how MLPs can be obtained in a straightforward and robust process from a single reference simulation. By essentially removing the need to adjust any hyperparameters, a new system of interest can be tackled in a direct, data-driven way. We have demonstrated the potential of this approach employing it for six complex liquid and solid-liquid systems and have evaluated the quality of the resulting models in detail for various properties underlining the high accuracy of our models. This important final step is realized with an automated validation protocol that is fully integrated into our framework. These developments are directly accessible to the community as they build exclusively on open-source solutions and we make our underlying software package and all templates available.

In its spirit similar to on-the-fly learning techniques, we depart from the goal of a high degree of transferability or generality to concentrate exclusively on the thermodynamic condition relevant for the chosen scientific question. Under these constraints, we have shown how a robust and accurate machine learning potential can be obtained with limited user input in an uncomplicated process. We note that we have also explored the application of our models for elevated and lowered temperatures, which was possible without problems in a $\pm$30 K regime. During such simulations, the intrinsic error estimate of our approach is a very useful tool to gauge the validity of the results. Due to these promising signs, we plan to systematically validate the robustness of our models to venture beyond the chosen thermodynamic state point, for example to describe variations in pressure.

We note that we built our development on established components, such as an efficient ML structure-energy representation and active learning concepts. At the same time, we expect that the concepts laid out in this work are transferable to other MLPs and active learning approaches, making it possible to achieve similar results. The novelty of our work lies in a change in perspective, where relatively little effort is put into creating an intermolecular potential which, in spite of concentrating on a sub-section of phase space, is still robust and accurate enough to be used to describe the non-trivial behavior of complex molecular systems. The importance clearly is not in the individual pieces but rather in the end-to-end framework and the broad range of applications made possible. This work therefore enables simulations that were not possible a short time ago, pushing forward the straightforward and reliable application of MLPs.

Looking to the future, we have limited ourselves to aqueous systems with four species in this work. However, we are confident that systems with more element types can also be tackled with our approach. In addition, we have not explored reactive processes in this study. Since MLPs are able to describe bond breaking and bond making events by design, we again anticipate the straightforward application of our methodology to such situations. This will be especially important to investigate interesting surface reaction phenomena, such as water dissociation on reactive surfaces or the recently reported reversible hydrolysis of zeolites in contact with water Heard2019/10.1038/s41467-019-12752-y. Key to a successful application in such situations will be the sampling of the relevant reactive process by the initial AIMD simulation used as input to our active learning protocol. Finally, we are currently exploring a training protocol in which structures are generated by classical molecular dynamics as input for our active learning protocol to minimize the need for expensive ab initio calculations. In this approach the expensive quantum computing engine is only used to obtain the ab initio potential energies and atomic forces for those configurations identified to be most important for the generation of the model. Such an approach has potential for additional cost savings over the one presented here and does not require expertise in AIMD simulations, thus more readily opening up the approach to researchers from the classical force field community.

In our six showcase applications of complex aqueous systems we have developed models with different DFT functionals, all of which represent reasonable choices for the aqueous systems studied. However, this illustrates a broader issue which is that there is currently no ‘perfect’ DFT functional for water and complex aqueous interfaces. We believe that the approach developed here could become a valuable tool in this long-standing quest to find suitable DFT functionals. Since our procedures provide the ability to reveal the true converged thermal performance of any given functional for a realistic system, at a modest cost, the systematic exploration of the performance of DFT methods for complex disordered systems becomes feasible. This makes it possible to go beyond the usual energetic benchmarks of relatively small systems in the absence of temperature and thus facilitates direct comparison with experiment. Due to the moderate size of our training sets, our framework is also expected to be easily extendable to more expensive ab initio methods, e.g. at the hybrid DFT level or considering explicit electron correlation, thus making these methods available for the realistic simulation of complex interfacial systems.

Overall, the developments reported herein will enable the investigation of complex aqueous processes such as water structuring in contact with interfaces and wetting or ice formation on surfaces in a straightforward manner. Although here we applied it to aqueous systems, we believe that the methodology will also prove useful for other materials and liquids in contact with solids, as well as general solvation phenomena, enabling the fast screening of different materials at ab initio accuracy. It will also be particularly useful for situations where long sampling is required as for the exploration of free energy surfaces, or calculations of dynamical properties, such as the friction or viscosity of liquids in contact with interfaces. In summary, this work outlines a straightforward strategy for the uncomplicated yet accurate investigation of many technologically and scientifically relevant systems by molecular simulations.

The introduced machine learning framework has been implemented in the AML Python package, which interleaves the required simulation packages and data manipulation steps in a user friendly environment. The AML package is available free of charge at https://github.com/MarsalekGroup/aml and enables the straightforward generation of a C-NNP model given a reference trajectory as input. With this code all six C-NNP models were developed for the various aqueous phase systems studied here.

NNP optimizations were performed with the open-source n2p2 code Singraber2019/10.1021/acs.jctc.8b01092 using the optimization parameters and symmetry functions as provided in the template file in the associated data repository for this paper. All additional information on the C-NNP fitting procedure can be found in the SI, while all training input files, training sets and parameters of the final models are publicly available at https://github.com/water-ice-group/simple-MLP.

The reference AIMD simulations used as the starting point for our C-NNP models employed quite different DFT settings, while all having been performed with the CP2K simulation package Kuhne2020/10.1063/5.0007045. We provide full detail about these reference simulations in the SI, but the typical simulation setups consist of 64 to 110 water molecules reaching simulation times between 30 to 130 ps.

MD simulations using the C-NNP models were also performed with CP2K, which features an open-source implementation of the C-NNP methodology since release 8.1. All C-NNP simulations for our validation protocol were propagated for at least 0.5 ns to allow for the converged computation of the RDF and VDOS. Further details of the validation protocol and the associated simulations can be found in the SI, while the validation can be performed with the AML Python package. The C-NNP simulations of the \chTiO2 water system were propagated for 5 ns for a 2$\times$2$\times$1 supercell of the original AIMD setup, resulting in total in 1728 atoms making up 320 water molecules on four O-Ti-O trilayers in a 23.672, 25.988, 42.0 Å periodic box. This simulation employed a molecular dynamics time step of 1 fs, while using deuterium masses for the hydrogen atoms. The temperature of 300 K was maintained with a canonical sampling through velocity rescaling thermostat Bussi2007/10.1063/1.2408420.

In order to validate the six models developed by our rapid machine learning framework, we have established a validation protocol that makes it possible to compare the performance of the various models for different systems in a direct and condensed manner. For that purpose we have selected three main categories for structural and dynamical properties as well as the precision of the force prediction. The main idea behind this validation protocol is that it probes the performance of a given model for the thermodynamic condition it is developed for. Thus, all three categories are compared directly against the AIMD simulation that was used as the starting point of the development of the models. This enables the straight-forward evaluation of the models without the need for additional benchmark simulations with the expensive reference method.

Structural properties for complex liquid-solid systems are directly probed by the radial distribution functions (RDFs) of the various species, which provide detailed insight into the two-component structural arrangement. For the scoring of the models, we compute all RDFs for a given system, both for the AIMD reference simulation and for independent C-NNP simulations using the developed model. For an $N$ component system this results in $\binom{N+1}{2}$ RDFs which can be directly compared between the AIMD and C-NNP results by using a suitable norm $d^{\text{RDF}}$

which provides a measure of the similarity of two RDFs ranging from 0 (for most different) to 1 (for identical). Averaging over all $\binom{N+1}{2}$ norms $d^{\text{RDF}}$ finally yields a single number that can be converted into percent to provide a condensed score for the performance of the C-NNP model for structural properties.

Dynamical properties are directly encoded by the vibrational density of states (VDOS), which is obtained by Fourier transform of the velocity autocorrelation function. The VDOS can be computed separately for all components in a system of interest and thus provides detailed insight into the dynamical properties of the system, probing vastly different processes over a broad frequency range. For an $N$ component system $N$ species-resolved VDOS are computed, both for the AIMD and C-NNP simulations. Using a similar norm $d^{\text{VDOS}}$ as for the RDFs

the AIMD and C-NNP results can be condensed into a measure of the similarity between the different functions, which after averaging over the different species and conversion into percent provides the score of a C-NNP model for dynamical properties.

Finally, the precision of the C-NNP model for the prediction of the forces is evaluated to provide another property score. The forces are what drives the dynamics of the system of interest and are thus of fundamental importance for an accurate description. We generated a test set for the evaluation of the force performance, by selecting a large subset of 1000 structures and associated forces from the original AIMD simulation. The root mean square error (RMSE), calculated separately for the $N$ species in each system

is used as a suitable measure for the force prediction. Since the magnitude of the forces can fluctuate strongly for different systems, but also within a given system (solid compared to liquid atoms), the RMSE is put into relation of the average force fluctuation of a given species

The resulting $N$ scaled force errors $F^{\text{Force}}=\frac{F^{\text{RMSE}}}{F^{\text{RMS}}}$ are averaged and converted into percent to provide the force score of the C-NNP model.