Machine Learning for Multi-fidelity Scale Bridging and Dynamical Simulations of Materials

Ever since the first successful demonstration of molecular dynamics (MD) to determine the diffusion coefficient of liquid argon more than 50 years ago [1], MD has become an indispensable tool to understand and predict materials properties for a broad range of applications, including drug discovery[2, 3, 4, 5], materials design[6, 7, 8] and defect chemistry[9, 10, 11]. Its generality, i.e., being applicable to a diverse range of materials (metals, ceramics, amorphous glasses, polymers, or biomolecules), along with its versatility, i.e, the ability to capture thermodynamic, mechanical, electrical and chemical behavior of materials, make it a pervasive and vital theoretical tool, as highlighted in Figure 1a. Increasingly, empirical findings are supported with analogous MD simulations to explain materials behavior at the electronic, atomistic or molecular level [12, 13]. Beyond its scientific merit, the popularity of MD simulations is fueled by the advancements in the computational power, and the availability of efficient software packages (VASP[14], LAMMPS [15], ASE[16], NAMD [17], CHARMM [18], GROMACS [19], etc.) that allow even a novice user (with limited background in MD simulations) to quickly simulate their experiments virtually. However, depending on the the phenomenon or material (metals, polymers, biomolecules or proteins) being modeled, smaller MD sub-communities have emerged to address specific modeling needs. Machine learning (ML) based force fields (FF), which is the focus of this perspective, provide an opportunity to bridge distinction between these different sub-communities, and provide a common and flexible methodology applicable across diverse materials domain.

At its core, classical MD simulations uses a physical model typically termed ‘force-fields’ to iteratively compute inter-atomic forces, which when used to integrate Newton’s equations of motion using numerical methods such as velocity Verlet for all atoms of a system, can reveal its dynamics. Depending on the flavor of the model used to obtain these inter-atomic forces, i.e., ab-initio, classical or coarse-grained, the nature of the simulations can vary drastically. First-principles methods, such as density functional theory (DFT) [20, 21], provide the most general and reliable force estimate but are computationally too exhaustive to be able to simulate phenomenon that involves larger system sizes (tens of nanometers to microns and beyond for understanding mechanical failure, microstructure evolution, to name a few) and/or occurs over longer time scales (nanoseconds and beyond for phenomena such as nucleation and growth, fatigue, to name a few). The empirical FFs or classical inter-atomic potentials community have successfully surmounted this challenge by developing exclusive analytical functional forms, involving tens to hundreds of tunable parameters, to approximate the inter-atomic force (or energy) expression. Once these potentials are fit to a known set of thermodynamic, mechanical or physical properties, they can be utilized to simulate materials dynamics under different conditions of interest. To efficiently treat even larger system sizes, such as biomolecules or proteins, coarse-grained representations are necessary, wherein rather than an all atomistic treatment, molecular motifs form the basic building block. In classical atomistic molecular dynamics, it should be noted that the terms force-fields and inter-atomic potentials are often interchangeably used to describe the functional form that captures the interactions between the atoms in the system. For coarse-grained simulations, the term force-field is used since no explicit atoms are considered. Figure 1b illustrates the different time/length scale regimes treated using aforementioned theoretical models.

Owing to the simplicity of the force expression, empirical FFs or coarse-grained models are several orders of magnitude faster than ab-initio molecular dynamics (AIMD)[22]. However, both of these former approaches suffer from limited accuracy and transferability when applied for cases far from those considered during the parameter fitting process. There thus exists a significant gulf between the various flavors of molecular dynamics ranging from the high-fidelity yet high computational cost ab-initio MD to low-fidelity yet efficient coarse-grained MD. Perhaps, an even bigger deterrent for FFs with fixed functional forms is the years of human effort required to find a respectable set of parameters, referred to as FF parameterization, that reliably describes the physics of a specific material system.

The advent of big data analytics and machine learning has brought to the forefront a new generation of powerful atomistic and molecular models [23, 24, 25, 26, 27, 28, 29]. This has recently propelled the development of AI-assisted or machine learning based FF that have displayed the promise of bridging the accuracy(fidelity)/flexibility gap between the AIMD and classical/coarse-grained models, while being sufficiently rapid to reach larger length and time scales. Machine learning and AI are playing a critical role towards addressing this challenge in multi-fidelity scale bridging. Furthermore, these FFs can be built relatively quickly and in an autonomous manner, without the need of human intervention.

In this perspective, we discuss two general directions in which the ML methods are emerging as indispensable tools for the MD community: 1) efficient and autonomous parameterization of FFs with known or pre-defined functional forms, and 2) directly learning the FF functional from the available high-fidelity AIMD or experimental datasets (i.e., establishing a configuration-to-energy/force mapping). We discuss several examples from each areas, along with their advantages and limitations, and finally conclude with some research avenues that remain to be addressed in the future.

Historically, FFs or potentials have been constructed by fitting an approximate atomic energy (or force) expression to reproduce tabulated empirical (lattice parameters, elastic constants, melting point, etc.) or quantum mechanical (phase energetics, cluster and interface energies, atomic forces, etc.) data. The complexity or functional form of the atomic energy expression is chosen based on the types of atomic interactions (ionic, covalent, dispersive, hydrogen bonding, and electrostatic) believed to be dominant in the considered material system, with the number of fitting parameters varying from two—Lennard-Jones potential [30] for simpler noble elements—to nearly 100—ReaxFF[31] and COMB [32] potentials for complex multi-elemental or reactive systems. The ingenuity of a FF developer lies in not only choosing an appropriate functional form for the material system of interest, but perhaps more importantly, in tuning the associated parameters to reproduce the known material behavior. Unfortunately, little guidance is available to solve this high-dimensional FF parameterization problem, with developers mostly relying on “chemical intuition” (in simple words - the expertise of the computational chemist) or traditional heuristics. Moreover, for complex systems (glasses, reactions, etc.) the number of tunable parameters increases dramatically and such heuristics become difficult to implement. Thus, there is no surprise that FF parameterization is considered to be a very painful exercise, often entailing several years of research efforts, and is pursued by only a limited number of scientists.

ML and AI based tools, however, offer an efficient and autonomous alternative pathway for FF parameterization, especially for complex systems involving large number of tuning parameters. Two methods that have particularly shown promise in this area are the Bayesian optimization (BO) [33] and the genetic algorithm (GA) [34]. BO can be used for FF parameterization by coupling it with a supervised ML model (e.g. Gaussian process regression, GPR [35]) that predicts the performance of a given set of FF parameters with some uncertainty. The main idea behind the use of BO is to find an optimal set of parameters that offers the best “exploration versus exploitation trade-off,” that is, the best balance between exploring new parameter space where the ML model is uncertain about the FF performance, and on the other hand focusing the search in the region where the ML model predicts FF parameters to have low errors. By iteratively following this search strategy, FF parameters can be found that minimizes the parameterization cost function, as illustrated in Figure 2a. Bauchy and co-workers used this approach to develop FFs for glassy silica [36], a material system known to be difficult to parameterize using traditional methods. They combined BO with GPR to tune interatomic potential of the Buckingham form that involves ten independent parameters and pose a difficult high-dimensional optimization problem. FF performance was evaluated using a cost function measuring the difference in the partial pair density functions (PDF) of Si and O obtained from the FF and the baseline AIMD simulations. To initiate the BO process, they first trained a GPR model that predicts the cost function, given the FF parameters. Then using the expected improvement [37] scheme they iteratively selected the new set of FF parameters that could potentially lead to reduction in the cost function. After only 600 iterations they were able to find a good FF, which not only reproduced the correct Si and O PDFs, but was significantly different from those of the baseline BKS potential—which illustrates the roughness of the cost function, and the challenge with traditional FF development approach. Notably, this scheme is general and can be extended to develop FFs for other material systems.

Another approach that has shown promise to parameterize FFs are evolutionary optimization methods such as GA, which solves the optimization problems using the principle of natural selection. In analogy with how nature uses the basic steps of crossover, mutation and selection for evolution of a specie, GA can be used to evolve a generation of FF parameters to reduce the errors in the parameterization cost function. In this approach, starting from a randomly selected generation of (nearly 100) FF parameter candidates, 1) crossover—mixing parameter values among different sets—and mutation—randomly perturbing parameter values—operations are performed to obtain new candidate sets; 2) MD computations are then performed to evaluate performance of these candidates using the parameterization cost function; 3) top candidates with low cost evaluations are then retained as parents for the next generation of parameter sets, with the above steps iterated until a candidate parameter set with desired performance is obtained. Although the GA approach offers no theoretical convergence guarantees [38, 39], empirically it has been shown to perform extremely well, especially for problems involving rugged cost function with several deceptive low-value local minimas. For example, GA has been used to fit ReaxFF for SiOH system [40].

A major bottleneck in the use of GA (or another optimization scheme) for FF parameterization is the evaluation of the cost function (step 2) for each candidate parameter set, which involves performing several expensive MD computations (to estimate correct phase energetics, high-temperature high-pressure stability, etc.). To overcome this problem, Chan and co-workers used a hierarchical objective function scheme, termed HOGA, presented in Figure 2b [41, 42, 43]. In this approach, the desired set of property objectives are segregated into different hierarchical classes based on their computational cost. During the cost function evaluation (i.e. step 2), the processes is truncated if a candidate parameter set leads to large errors in a hierarchical classes, with the associated class penalty being assigned to the candidate set. Thus, the HOGA scheme accelerates the traditional evolutionary search by (i) efficiently sampling the parameter landscape in a given generation, and (ii) overcomes the limitation of a single objective that relies on assigning arbitrary weights. Further, it allows multi-fidelity scale bridging by incorporating training datasets at various length scales, i.e., phase energetics and structural properties from first principles or atomistic models, experimental measurements, and even on-the-fly MD based property evaluations. The sequencing and selection of hierarchical classes is at the discretion of the user, which can be modified based on the target property of interest. The success of the HOGA scheme was demonstrated by developing FFs for H₂O [41] and WSe₂ [42] systems using bond order potential model with 13 independent parameters. The coarse-grained (CG) models for water accurately captures its density anomaly, which has been a major limitation for other CG models. Similarly, WSe₂ FFs were used to study thermal transport behavior of diverse nanostructures at different temperatures using non-equilibrium MD simulations. In contrast to the traditional local optimization approaches, these evolutionary optimization methods allow for global optimization of the objective function. The success of these evolutionary approaches stems from their ability to more efficiently navigate the high-dimensional parameter space than conventional local optimizers or gradient based approaches [39].

Although BO and GA provide efficient ways to parameterize FFs, a key aspect that remains missing from the whole process of FF development is the generation of training data. In this regards, Chan and co-workers developed an autonomous framework termed BLAST—Bridging Length/Time scales via Atomistic Simulation Toolkit—that allows users to create their own potentials/FFs by following through each stage of FF development, i.e., generating appropriate training data sets, optimizing FF parameters, and finally cross-validating their FF predictions [23]. It employs a data-driven and ML based FF development approach that significantly deviates from the traditional approaches in several key aspects, including the nature of the training dataset, the use of ML algorithms for advanced sampling and FF parameter optimization (e.g., genetic algorithms, multi-objective optimization, neural networks, etc.), the available form of FFs, and cross validation and iterative improvement of ML based potentials. This approach has been successfully demonstrated for diverse material systems, including metal nanoclusters, oxides, nitrides, and 2D materials (stanene and WSe2)[44, 45, 46, 23, 42, 47].

Finally, we note the following in the area of FF parameter optimization using machine learning. Although optimization procedures such as BO and GA remain attractive to optimize the coefficients or parameters of any given FF, we envision the use of other emerging AI techniques, for instance, Monte Carlo Tree Search (MCTS) [48] with reinforcement learning to derive globally optimal solutions for these parameter search problems. We expect burgeon of automated workflows (e.g. BLAST - Bridging Length/timescales via Atomistic Simulation Toolkit [23]) that allow non-expert users to carry out their own potential model development depending on the target system and phenomenon. The term ’workflow’ describes the use of scripting languages to interface with various domain-specific simulation codes and/or execution of simulation campaigns across multiple HPC and cloud resources. These can be achieved via python or dedicated parallel scripting languages, such as SWIFT [49]. Likewise, other dataflow systems like Decaf developed by Tom Peterka allow for the parallel communication of coupled tasks in an HPC workflow [50]. This new development democratizes the FF development process and reduces the dependency of the MD community on a handful of FF developers for optimization of FF parameters. Furthermore, this approach of FF development is believed to remain popular as it returns a FF functional form which is physically meaningful, interpretable and much faster than a complete ML based FF, which is discussed next.

Despite the popularity of pre-defined functional forms, it has become imperative to address the limitations arising from their lack of flexibility. The use of pre-defined functions to represent the inter-atomic interactions inherently imposes a limit on its predictive power. This has led to the search for alternate more flexible, accurate and efficient models. Rather than using ML to find parameters of a FF with fixed functional form (bond order potential (BOP) [51], reactive force-fields ReaxFF [31], etc.), therefore, a different class of ML FFs have emerged in the last decade that aim at directly learning the functional form of the FF itself, i.e. establish a direct mapping from the atomic neighborhood, or fingerprint, to atomic energies—and their sum to total potential energy—using the reference quantum mechanical data. As no prior restriction is imposed on the functional form of these FFs, this methodology is quite general and can be used to learn atomic energy functionals of diverse materials (metals, ceramics, alloys, polymers) involving different atomic interactions with minimal human interference. This also precludes the limitations of the traditional FFs with fixed functional form designed for capturing specific inter-atomic interactions. This approach has been successfully used to develop ML FFs for numerous systems, including elemental bulk (e.g. Al [52, 53, 54], C [55, 56, 57], Li [58], Si[59, 60, 53, 61], Fe [62], Zr [63]), alloys[64], metallic clusters [65, 66], semiconductors [67], oxides [68, 69], water [70, 71], organic molecules [72, 73], among others, and study complex applications, such as surface diffusion [74], liquids, phase equilibria, etc. Different software packages (RuNNer [59], GAP suite [75], AMP [76], AGNI [52], DeepMD [77], AENet [69], MTP [78], SchNetPack [79], N2P2 [80], SNAP [81], etc.) that allow users to train and validate their own ML FF for a particular system of interest have also been released. Despite its popularity and success over the past several years, it is still an active area of research with advancements in atomic fingerprinting scheme or the supervised ML task (atomic positions to energy functional mapping) being reported frequently. Thus, below we provide only an essence of this approach, and refer the reader to detailed reviews on this topic for more information[82, 83, 84, 85, 86, 87, 88, 89, 25, 23].

Figure 3 shows a general schema typically used in the construction of these ML based FFs that consists of three stages: database generation, fingerprinting and supervised learning. In the first stage a reference high-fidelity first-principles (e.g. AIMD) database spanning diverse material configurational space and associated energetics, stresses and atomic forces is constructed. In the next step, a fingerprinting scheme is utilized to numerically represent the atomic configurations present within the reference dataset, which are then mapped to the associated atomic energies—or more formally, their sum to the total potential energy—in the last step using a supervised ML algorithm. Within each of these three stages, numerous methods have been proposed to improve their accuracy or speed, thereby creating different flavors of ML FFs in the community.

Perhaps the most critical component in the construction of these ML based FFs has been the development of a suitable atomic fingerprint in stage 2, as these fingerprints are required to be strictly invariant with respect to arbitrary translations, rotations, and exchange of like atoms, in addition to being continuous and differentiable with respect to small variations in atomic positions. Several candidates, including those based on atom-centered symmetry functions (ACSF) [59], smooth overlap of atomic positions (SOAP)[75, 90, 91], moment tensor representation [78], spectrum of London and Axilrod–Teller–Muto (SLATM) [92], the Faber–Christensen–Huang–Lilienfeld (FCHL) representations [93, 94], bispectra of neighborhood atomic densities [81], Coulomb matrices (and its variants) [95, 96], and others, have been proposed. A common theme across all of these fingerprinting approaches is to numerically capture atomic neighborhood around a reference atom using some form of distribution functions as qualitatively illustrated in Figure 3a. The accuracy vs speed trade-off among these different descriptors is being presently studied [97], but the ACSF and SOAP representations have been most successful in terms of number of applications. We, however, note that other descriptors have also been shown to produce reliable FFs, hinting that many of these representations maybe accurate enough for ML FF generation.

In terms of the first data generation step, researchers have proposed different MD or sampling techniques to generate high-quality diverse database with limited computational budget [98, 99, 100, 53]. Learning algorithms, ranging from neural networks [59], support vector machines [101], kernel ridge regression [53, 95], GPR [75], etc. have been utilized to carry out the supervised learning task in stage 3. Although it is not clear which supervised learning algorithm performs best, deep neural networks are expected to be a top choice because of the data intensive nature of this problem.

A relatively new development within the area of ML FFs is to learn and predict the atomic forces directly [102, 52, 103], as this is the fundamental quantity required to perform MD simulations. These approaches are inspired by the idea that it should be possible to predict atomic forces given just the atomic configuration, without going through the agency of the total potential energy. A major advantage of this approach is that atomic force can be uniquely assigned to an individual atom, while the potential energy is a global property of the entire system; partitioning the potential energy to atomic contributions does not have a formal basis, but is often assumed in many of the ML FFs mentioned earlier. Mapping atomic fingerprints to purely atomic properties can thus lead to powerful and accurate prescriptions. One such methodology called AGNI has been shown work well for a variety of elements (Al, C, Pt, etc.) for a variety of mechanical and thermal properties [52, 104, 105, 106].

Another promising development in this area has been the use of graph networks [107, 108, 109, 110], which by construction overcome the limitations of manually constructed atomic fingerprints (or feature vectors). A graph representation of molecules or crystals consists of feature vectors capturing information about all the constituting atoms (nodes/vertex), and bonds or atomic distances (edge). This graph representation can be fed into a ML model (neural networks) and can be mapped to material property of interest. So far this methodology has not been used for FF development, but the unprecedented level accuracy achieved by such models on the Materials Project [111] and the QM9 databases [112] to predict diverse properties, such as energy, polarization, vibrational frequency, band gap, bulk modulus, etc., offers great promise.

Previously, we discussed how GA can be used to effectively find parameters of a FF with known functional form. However, GA can also be utilized to identify new functional forms themselves using symbolic regression [113]. The key to this approach is to search within a physically meaningful hypothesis space. Analysis of interatomic potentials derived over the past several decades from physical principles, reveals that only a few fundamental operations are necessary to recover functional form of these classical potentials. These operations include addition, subtraction, multiplication, division, and power operators; constant values and distances between atoms; and an operator that performs a sum over functions of distances between a given atom and all neighbors within a given cutoff radius. Thus, if the GA is used to search within a combinatorial space of these physically meaningful operations, new FF functional forms can be learned that better fit the available quantum mechanical or empirical data. Symbolic regression is another crucial piece of this methodology which transforms a given functional form into a tree structure, as shown in Figure 3c, to make it amenable to basic GA crossover and mutation operations. Mueller and co-workers [114] successfully used this scheme to not only recover the exact form of famous Lennard-Jones and Sutton-Chen (SC) EAM potential, but also build ML potentials for Cu using DFT training data, with an exemplary potential being $E_{i}=\sum(r^{10.21-5.47r}-0.21^{r})f(r)+0.97\big{(}\sum 0.33^{r}f(r)\big{)}^{-1}$.

A key limitation of any FF, whether traditional or ML based, is that they may not work for applications (simulations) involving configurational domains that are different from their training domain [115]. As such the transferability of the potential function or force-field is an important measure of its predictive power. For example, FFs are known to be unsuitable for simulations far from equilibrium, as they are mostly trained to reproduce equilibrium properties and configurations. Likewise, neural network based potentials are interpolative and as such tend to fail miserably when encountering scenarios outside of their training set [116, 117]. Despite these challenges, the ML approach does provide an opportunity to systematically improve versatility and transferability of FFs over the traditional pre-defined functional forms.

The interpolative nature of the ML FF, especially those based on NN, demands a large number of training data to ensure that all possible physical scenarios are adequately captured. Their training, therefore, often requires large quantities ($\sim$10⁴ or greater) of data to ensure that the model adequately samples the energy landscape both near and far-from-equilibrium. A highly desirable goal is to minimize the number of training examples, especially if the underlying reference model is first-principles based and hence computationally expensive. This has led to the emergence of ML techniques such as ‘active learning’ and ‘transfer learning’ that focus on efficient sampling of the training data. The overarching goal is to reduce the number of training data to the absolute minimum without compromising on the adequate representation of various areas in a potential energy surface.

There have been many recent studies aimed at addressing the sampling challenge in generation of training data. Smith et al. [118] employed an active learning (AL) strategy based on the Query by Committee (QBC) scheme. QBC uses the disagreement between an ensemble of ML potentials to infer the reliability of the ensemble’s prediction. QBC allowed for automatic sampling of the regions of chemical space where the potential energy was not accurately described by the ML potential. Their AL approach was validated on a test set consisting of a diverse set of organic molecules and their results showed that one requires only 10% to 25% of the data to accurately represent the chemical space of these molecules. Similarly, Zhang et al. [119] introduced an AL scheme (deep potential generator (DP-GEN)) that constructs ML models for simulating materials at the molecular scale. Their procedure involve exploration, generation of accurate reference data, and training. They used Al and Al-Mg as representative cases and showed that ML models can be trained with minimum number of reference data. Another interesting body of work in this direction is on-the-fly learning to accelerate AIMD. These methods rely on building FFs that not only provide point estimates of energy and atomic forces, but additionally outputs associated prediction uncertainty. Ceriotti and co-workers have developed an inexpensive yet reliable way of estimating this uncertainty associated with ML model predictions of atomic and molecular properties. Their scheme is based on resampling wherein they generate multiple models based on subsampling of the same training data. They benchmark the accuracy of the uncertainty prediction by maximum likelihood estimation which in turn can correct for correlations between resampled models and improve performance of the uncertainty estimation via a cross-validation procedure [120]. By tracking model uncertainty during the MD simulation, a call to the high-fidelity (but expensive) DFT calculation can made when the system drifts to a configuration where the model is uncertain in the energy/force prediction beyond a certain threshold [102, 121, 122]. Vandermause et al. [123] sampled structures on-the-fly from AIMD and used an adaptive Bayesian inference method to automate the training of low-dimensional multiple element interatomic force fields. Their AL framework uses internal uncertainty of a GPR model to decide acceptance of model prediction or the need to augment training data. In all of the above studies, the overarching aim is to minimize the ab-initio training data required to describe the potential energy surface.

In another exemplary work by Huan and co-workers [104], “failed” data augmentation was used to systematically improve the domain of applications (and configurations) that can be handled using the AGNI ML FF. They systematically tested their ML FFs on new applications (stress/strain analysis, stacking fault energetics, and melting simulations), identified configurations where the present version of FF failed, appended those configurations to the original training data, and retrained the ML model on the updated training data, thereby increasing the domain of applicability of their FFs. Another crucial information that was transferred from the previous FF version during the retraining process was all the associated ML parameters, e.g., the fingerprint parameters, the learning algorithm, and the previously selected training data. This helped transfer the current performance/knowledge of the FF to the next generation which was then trained using the augmented dataset. More recently, Loeffler et al. [124] have developed an active learning approach that iteratively trains an ANN model to faithfully reproduce the coarse-grained potential energy surface of water clusters. They showed that one can train a NN potential with only 426 total structures in its training data. Their AL workflow initiates with a sparse training data set and is continually updated via a Nested Ensemble Monte Carlo scheme that sparsely queries the energy landscape and tests the network performance. The network is retrained with an updated training set that includes failed configurations/energies from previous iteration until convergence is attained. Such a trained network adequately reproduces the energies (within mean absolute error (MAE) of 2 meV/molecule) and forces (MAE 40 meV/Å) compared to the reference model. Such studies have laid the groundwork on for workflows that allow on-demand generation of training data.

With the recent advances in computing power, AI, ML and big data analytics are emerging as powerful tools for materials modeling. The flexibility, accuracy and computational efficiency of ML models can accelerate materials discovery and design and can also be leveraged for advanced nanoscale materials characterization. The ML models are primarily geared towards addressing one of the most critical aspects of MD i.e. its predictive power which hinges strongly on the nature of the interatomic potential used to describe the atomistic interactions in the system. The various examples presented in this perspective article highlight the advantages of using ML and data driven approaches for training new models. There is, however, still a lot of room for improvement and we believe that both supervised and unsupervised ML methods can play a key role towards the development of next-generation atomistic, molecular or coarse-grained models. Some of these areas are briefly discussed below.

MD with pre-defined functional forms continue to remain attractive owing to their computational efficiency. One advantage of pre-defined force-fields is that they are usually based on physical models that attempt to reliably capture atomistic interactions. This physical basis imparts predictability and transferability for such models. However, most functional forms are restricted to first order treatments of the physical models and ignore higher level terms owing to the increase in computational cost. This is especially true for polarizable or reactive systems, where the treatment of electrostatics still relies on point charges on atoms. To adequately describe a molecular dipole moment and its electrostatic field, one needs to include dipole moments of the atomic charge distributions. While this was established decades ago by Stone AJ [125], such terms are typically not included owing to their computational complexity, and therefore the accompanying polarization effects directly described. Likewise, the charge transfer dynamics in the case of reactive systems is still described by methods such as electronegativity equalization method (EEM) [126] and charge equilibration (Qeq) [127]. The charge calculations represent the most compute intensive part of reactive potential models. These charge transfer and polarizability effects can have strong significance for organic systems, and for glasses and ceramics, especially at the surfaces and interfaces. The advances in ML and the automated training processes should enable us to revisit and efficiently describe effects arising from higher order terms to adequately capture electrostatic contributions in polarizable and reactive systems.

While pre-defined functional forms are efficient, the choice of a pre-defined function limits their ability to capture the underlying chemistry and physics being modeled. Despite a lot of advances in sampling and generation of training data as well as improvements in global optimization algorithms, there is always going to be a ceiling limit for pre-defined functional forms. Flexibility in the functional form is one of the critical challenges. The re-emergence of neural network models in the post-AI winter era is aimed at providing more flexibility to users interested in modeling complex reactive systems that involve multiple bonding characteristics (e.g. metallic, covalent and ionic). However, the advantage offered by such flexible ML potentials can also become their limitation; when such ML potentials are used outside their domain of applicability (far from training data) they result in arbitrary unphysical behavior. Thus, more research on constraining these models to basic physical relations, such as the work on physics informed neural networks [116], will be highly useful.

All ML frameworks require carefully curation and sampling of training data. In the context of materials models, one should ensure that the training data used is sufficiently diverse (for instance, training set that consists of configurations spanning broad range of energies from near equilibrium to highly non-equilibrium). For the model to be robust, one has to ensure ample representation of the different parts of the potential energy surface i.e. both near equilibrium and far-from-equilibrium configurations in their training data. The availability of first-principles training datasets via Materials Project and other thermodynamic databases puts us in an enviable position to incorporate robustness in the ML developed models. However, it still is not properly understood what is the minimum diversity and size of training dataset required to build robust FFs. On the other hand, advances in workflows have also allowed for direct interfacing with electronic structure codes (e.g. VASP) allowing users to generate training data on-the-fly. This is an important emerging area and recent studies, especially on active and transfer learning strategies, highlight the importance of sampling and incorporation of diversity in training data.

In this respect, one of the significant departures in MD potential development compared to conventional fitting is the development of automated training workflows that explicitly include temperature dependent properties in its objective function. Traditionally, one would define objective functions that would primarily include the static properties such as configurational energies, cohesive energies, lattice constants to name a few. Such models, however, lack predictive power to capture dynamical and transport properties. Our recent work has introduced workflows that allow temperature dependent properties derived from on-the-fly molecular dynamics to be included as part of the training procedure. The ability to directly train MD potentials that capture transport properties and other temperature/pressure dependent properties is of tremendous significance to capture dynamical phenomena near phase boundaries. In this regard, Chan et al. trained a bond order potential for WSe₂ which was used to study thermal transport behavior of diverse nanostructures at different temperatures using non-equilibrium MD simulation[42]. The potential was fit to not only first principles data, but the temperature dependent stability of WSe₂ structures was also among the fitting criteria. The idea is to run on-the-fly MD simulations during the potential fitting process to assess the stability of various structures at different target temperatures[41, 42].

It should be noted that MD simulations with classical Newtonian mechanics ignore quantum zero-point energy contributions to structural fluctuations. The quantization of the energy of the vibration modes cannot be accounted for by using standard MD because it is based on classical statistics. In this regard, Dammak et al. [128] introduced a quantum thermal bath (QTB) that accounts for quantum statistics while using standard MD. The basic idea of the QTB is to use a Langevin-type approach. Barrat et al. implemented QTB in a flexible manner and at a low computational cost by synthesizing the corresponding noise ‘on the fly’ [129]. These corrections are especially important when evaluating thermal and other transport properties. We note that the effect for metals (especially transition elements) is generally small, but it becomes more important for glasses and ceramics, and very significant for organics. Thus, depending on the material type, such corrections should not be ignored and become important when one includes temperature dependent quantities in future training procedures.

Multi-fidelity scale bridging from the electronic structure to coarse-grained scales emphasizes a greater need to estimate the extent and nature of error propagation from the high-fidelity first-principles scale to coarse-grained mesoscopic scale. With the increasing use of higher-level theory (CCSD or QMC) to generate the training data and reduce the error necessitates quantification of uncertainties in the predictions at various scales. In this regard, Functional Uncertainty Quantification method (FunUQ) [130] to assess model form errors, such as interatomic potentials for MD simulations represents an important future direction. Cross-validation, sensitivity analysis and uncertainty quantification will be critical to test the robustness of their developed models. Preliminary work on the use of uncertainty quantification has also lead to autonomous MD workflows such as FLARE [123] that provide a merger between high-fidelity quantum mechanical and low-fidelity machine learning potentials.

We note that there are limitations to DFT—for instance, vibrational frequencies/phonon spectra calculated by DFT can have errors of $\sim$10% and this could represent a source of error in the training set. One therefore requires some sort of calibration to bring the calculated data in closer agreement with experiments (reality). This could be achieved by performing more higher level quantum calculations (CCSD or QMC) to capture the errors in DFT and establish some scaling relationship or via transfer learning where the weights of a neural network are first trained using DFT and then correctly by performing a limited number of higher level DFT calculations.

ML global optimization algorithms and multi-stage optimization strategies clearly represent a major advance over local optimization procedures that have been traditionally employed for optimizing parameters to describe the potential energy surface. Evolutionary optimization procedures such as genetic algorithms remain attractive but one can easily envision the use of other emerging AI techniques such as Monte Carlo Tree Search (MCTS) [48] to derive globally optimal solutions. The use of decision trees and dynamic learning to navigate more effectively through the potential energy landscape is expected to address some of the bottlenecks with the sluggishness of global optimization schemes.

A key aspect of optimization is to appropriately define objective functions that go beyond the simple sum of square difference. The selection of the appropriate weights poses major problems for optimization that focus on single objective functions. The weight selection in a single objective function is typically dependent on the expertise of the force-field developer. The goodness of fit (or the predictive capability of the force-field) is measured by a weighted sum of errors in predictions of a set of target observables. The general rule of thumb is that the weights should be selected such that the contribution to the overall error for each of the individual properties in the objective function is approximately the same. Nonetheless, the weight selection is quite challenging and typically, the force-fields developed tend to do well in a few properties and rather poorly in others—even though the overall objective function has reached a minimum. To circumvent the issues with weight selection, there have been developments in multi-objective strategies where each of the objectives is defined separately. One such strategy is based on Pareto optimization, where one retains a set of optimal solutions and the user can choose a solution depending on the properties of interest. A Pareto-optimal solution is one in which a small variation of parameters that leads to an improvement in any of the objectives will lead to at least one of the other objectives being less optimal; there can be curves or surfaces in the parameter space corresponding to such solutions and these are termed Pareto fronts. Such multi-objective evolutionary approaches are especially helpful if one is dealing with conflicting properties (e.g. elastic constants vs. defect formation energies). Other multi-objective approaches that have recently emerged are those that treat a hierarchy of objectives during an evolutionary search[41]. In this case, the quality of a proposed parameter set is evaluated based on a hierarchical objective function. In this evolutionary scheme, we truncate the evaluation of any parameter set which leads to large errors in hierarchical property classes and assign it a penalty depending on which class it fails at. This approach aids in an accelerated evolutionary search by (i) efficiently sampling the parameter landscape within a given generation, and (ii) overcoming the limitation of assigning arbitrary weights within a single objective thereby ensuring that all the properties (static or dynamic) are equally well described. Such an approach has been used recently to train a coarse-grained model for water. This model is able to successfully capture the various thermodynamic anomalies of water and outperforms most other existing atomistic and coarse-grained models in their predictive power[41]. Multi-objective optimization such as Pareto optimization and the recent success of HOGA (hierarchical objective genetic algorithm) [41] to find optimal solutions that represent a compromise between the various desired objectives is a step in the right direction. One needs to exercise caution in defining the objective appropriately otherwise even the best global optimization strategy may not yield the best results.

Finally, while a lot of work has focused on length-scale bridging, timescale challenges remain unaddressed. All-atom simulations typically employ 0.1-1 femtosecond timesteps and as such are suited for modeling phenomena in the nanosecond timescales. Even with exascale computing, one would primarily gain in terms of spatial scales that can be accessed. The clock-speeds and bandwidths are not going to change significantly and timescale challenges will remain. Coarse-graining allows for larger time-step during MD integration and thereby explore longer time-scales. In particular, coarse-grained models typically allow 10-50 femtosecond timesteps and represent one route to address this timescale challenge albeit typically at the cost of accuracy. The use of coarse-grained models along with advanced sampling techniques is becoming very popular among the community for exploring longer time-scales, exploring phase-transitions such as Lower/Upper Critical Solution Temperature LCST/UCST in polymeric solutions and brushes. In particular, the MARTINI coarse grain model [131] has gained in popularity and is increasingly being used to understand drug uptakes and transport across lipid bilayer membrane and membrane based dynamics. The future looks promising with the advent of quantum computing and the use of physics-based polarizable models coupled with path integral calculations for accurately determining the free energy of binding of drugs and also exploring of time-scale where rare events such as allosteric modulation, molecular recognition, protein folding and self-assembly. To explore such time-scales, chemists are resorting to polarizable united-atom and coarse-grained models that can more adequately treat electrostatics.

One of the major challenges with coarse-grained models is that the gain in efficiency comes at the cost of accuracy. This is particularly true for water, which is one of the most important solvent used in soft-matter simulations of polymers, proteins and other biomolecules. Traditionally, models to describe water range from fixed charge dispersion/repulsion 12-6 model, Buckingham potential, Stockmayer potential, ST2, force fields due to fluctuating charge, and the inclusion of van-der Waal parameters for explicit hydrogens have been explored in the past 80 years [132]. Further refinements include the introduction of sites on the HOH angle bisector, multiple sites $>$3 for water based on the TIPnP framework. The target data for water includes the phase diagram, ice-formation abilities, structural properties (radial distribution functions), transport properties including diffusion, thermal conductivity and viscosity not to forget the anomalous density variation with temperature of dense water. ML strategies can allow for development of coarse-grained surrogates that address the timescale challenges without sacrificing accuracy. ML CG water models developed by Chan et al. [41] that allow for on-the-fly sampling of temperature dependent properties have succeeded in capturing the various thermodynamic properties and anomalies of water at fraction of the computational cost of its atomistic counterpart. Such coarser surrogate model development augurs well for capturing hitherto unexplored dynamical phenomena in the mesoscopic regime.