Machine Learning Classical Interatomic Potentials for Molecular

Subscriber access provided by University of Kansas Libraries

Feature Article

Machine Learning Classical Interatomic Potentials for Molecular Dynamics from First-Principles Training Data Henry Chan, Badri Narayanan, Mathew J. Cherukara, Fatih G. Sen, Kiran Sasikumar, Stephen K. Gray, Maria K. Y. Chan, and Subramanian K.R.S. Sankaranarayanan J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b09917 • Publication Date (Web): 03 Jan 2019 Downloaded from http://pubs.acs.org on January 3, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Machine Learning Classical Interatomic Potentials for Molecular Dynamics from First-Principles Training Data Henry Chan1§, Badri Narayanan1, 4§, Mathew J. Cherukara2, Fatih G. Sen,1 Kiran Sasikumar, 1 Stephen K. Gray1, Maria K. Y. Chan1, Subramanian K.R.S. Sankaranarayanan1, 3* 1 Center

for Nanoscale Materials, Argonne National Laboratory, Argonne, IL 60439

2 X-ray 3 4

Science Division, Argonne National Laboratory, Argonne, IL 60439

Institute of Molecular Engineering, University of Chicago, Chicago, IL 60637

Mechanical Engineering Department, University of Louisville, Louisville, KY 40292

§ Equal contributions * Corresponding author – [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract: The ever-increasing power of modern supercomputers, along with the availability of highly scalable atomistic simulation codes, has begun to revolutionize predictive modeling of materials. In particular, molecular dynamics (MD) has led to breakthrough advances in diverse fields, including tribology, catalysis, sensing, and nanoparticle self-assembly. Furthermore, recent integration of MD simulations with X-ray characterization has demonstrated promise in real-time 3-D atomistic characterization of materials. The popularity of MD is driven by its applicability at disparate length/time-scales, ranging from ab initio MD (hundreds of atoms and tens of picoseconds) to all-atom classical MD (millions of atoms over microseconds), and coarse-grained (CG) models (microns and tens of micro-seconds). Nevertheless, a substantial gap persists between AIMD, which is highly accurate but restricted to extremely small sizes, and those based on classical force fields (atomistic and CG) with limited accuracy but access to larger length/time scales. The accuracy and predictive power of classical MD simulations is dictated by the empirical force fields, and their capability to capture the relevant physics. Here, we discuss some of our recent work on the use of machine learning (ML) to combine the accuracy and flexibility of electronic structure calculations with the speed of classical potentials. Our ML framework attempts to bridge the significant gulf that exists between the handful of research groups that develop new interatomic potential models (often requiring several years of effort) and the increasingly large user community from academia and industry that applies these models. Our data-driven approach represents significant departure from the status quo, and involves several steps including generation and manipulation of extensive training data sets through electronic structure calculations, defining novel potential functional forms, employing state-of-the-art ML algorithms to formulate highly optimized training procedures, and subsequently developing user-friendly workflow tool integrating these algorithms on high-performance computers (HPCs). Our ML approach showed marked success in developing force fields for a wide range of materials from metals, oxides, nitrides, heterointerfaces to two-dimensional (2-D) materials.

2


Page 2 of 47

Page 3 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Over the last decade, increases in high performance computing (HPC) power, along with algorithmic improvements, have significantly advanced the state-of-the-art in molecular understanding of various material interfaces1-3. These advances, in turn, have accelerated the development of new materials, or novel design concepts for various energy-related technologies, including sensing4, nano-electronics5-7, catalysis3, and energy storage8-10. Despite such marked developments in materials modeling, efforts to accurately capture atomic-scale dynamics, especially in low dimensional systems such as clusters and interfaces (commonly encountered in technological materials) are still in their infancy. The prime challenge in clusters and interfaces are their association with myriad inter-connected dynamical processes, including (a) chemical reactions, (b) transport of atoms and ions, (c) defect chemistry, (d) solvation dynamics. The complex interplay of these processes at nano- to meso- scales dictates micro-structural evolution, as well as stability and macroscopic performance of the devices (e.g., battery capacity and lifetime). Evidently, a clear understanding of dynamical processes at clusters and interfaces is critical for the design and development of novel functional materials across a broad range of energy applications. As such, achieving this knowledge, and subsequently, the ability to engineer dynamics at clusters and interfaces to obtain desired microstructure have emerged as key challenges in materials and interfacial science11. To address these issues in a meaningful way, an accurate description of atomic interactions in nanoscale clusters and across interfaces (e.g. solid-liquid) is crucial. Interatomic potential models that capture these interactions have been notoriously hard to develop due to widely different (and often rapidly changing) chemical environments (including ionic, covalent, dispersive, hydrogen bonding, and electrostatic forces) exhibited by such low dimensional systems.

3



a)

b)

Figure 1. Comparison of the computational cost associated with various flavors of MD simulation. (a) A typical first principles MD simulation – ion migration in perovskite nickelate5. With ~200 atoms, a 10 ps simulation consumes 100,000 core-hours of computer time. (b) A typical classical MD simulation – growth of silicene on Ir(111)6. With 100,000 atoms, a 1 ns simulation consumes ~10,000 core-hours. Density functional theory (DFT) can accurately capture diverse bonding environments within a single formulation, and is well-equipped to model clusters and interfaces (both reactive and non-reactive). Importantly, the power of current supercomputers, and the scalability of available DFT codes have made ab initio molecular dynamics (AIMD) simulations quite feasible. However, there is still a substantial gulf between the length/time-scales accessible to AIMD and classical molecular dynamics (MD) as illustrated by Figure 1. State-of-the-art AIMD simulations can simulate the trajectory of a few hundred atoms over tens of picoseconds; on the other hand, microsecond long trajectory for millions of atoms has been achieved using MD simulations powered by empirical force fields (FFs). Note that the atomic-scale processes underlying the functionality of low dimensional systems occur over angstroms to microns, and picosecond to hundreds of nanoseconds12. Although additional computing time can push the AIMD accessible timescales to ~100 ps, the dense linear algebra involved in DFT calculations renders scaling AIMD to much larger length scales (close to those accessible to classical MD) impossible. Evidently, the broad ranges of

4


Page 4 of 47

Page 5 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


time and spatial scales relevant to processes of interest in multi-component clusters and interfaces cannot be met without a new generation of highly accurate and robust, yet computationally efficient, interatomic potentials (force fields). Historically, only a limited number of research groups have developed interatomic potential models (often entailing several years of effort for each potential model), and subsequently disseminate these models to a broad user community from academia and industry. Typically, there is a disconnect between the force field developers and end users, which often results in inadvertent use of force fields to model scenarios where they performs poorly or are completely out of scope. Furthermore, even advanced users do not have the flexibility to adapt these pre-defined potential models to problems of their interest. The advent of bigdata analytics and easy access to HPC resources has brought powerful machine learning (ML) techniques to the forefront. These ML methods hold tremendous promise in resolving the disconnect between force field developers and the end-users, as well as empowering the users to develop new or tailor existing FFs to meet their needs. This feature article discusses some of the recent work carried out by our group at Argonne to overcome this barrier13-17. We are developing an automated framework that allows users to create their own models by generating training data sets, optimizing potential functions, and cross-validating their model predictions. Our framework employs a data-driven and ML based force field development approach that significantly deviate from the status quo in several key aspects: the nature of training dataset, the use of ML algorithms for advanced sampling and optimization schemes for fitting (e.g., genetic algorithms (GA), multi-objective optimization, neural networks to name a few), the form of the force fields, our treatment of charge transfer for reactivity, and cross validation and iterative improvement of ML based potential models.

2. Current state-of-the-art in empirical force-fields for clusters and interfaces A vast body of literature exists on development of FFs for a wide range of material classes with varying degrees of complexity. In the framework of most FF models, the potential energy of a system is treated as a sum of various short-range covalent (bonded) and long-range (non-bonded) interactions

5



between atoms. Typically, each of these partial contributions is written as a function of atomic positions, relative orientations, and charges. The contribution to the potential energy from bond interactions arise from bond stretching, angle bending, dihedral torsion, and over/under-coordination. These interactions are also classified as two-body (pair-wise), three-body (angle), four-body (conjugation) interaction terms. Note that these classifications are fairly arbitrary, and are a manifestation of varied philosophies among forcefield developers in different research domains (e.g., those working on soft-matter as opposed to inorganic systems). Non-bonded interactions, on the other hand, are primarily dispersive (e.g., van der Waals) or electrostatic in nature. Broadly, the mathematical formulation commonly employed for FFs for inorganic materials can be divided into two classes: (1) those that are composed of distinguishable two-body and three-body (and higher) energy terms, such as Stillinger-Weber (SW)18, and Vashishta potentials19, and (2) those that treat many-body interactions collectively, either as a function of atomic density (e.g., Finnis-Sinclair20, embedded-atom method (EAM)21, modified EAM (MEAM)22) or using a bond-order concept23-26 (e.g., Abell, Chelikowsky, Pettifor’s, Tersoff, and reactive empirical bond order potential (REBO)27). Despite the successes of these broad classes of FFs for wide range of material systems (including metals, covalent, or ionic solids), they are plagued with inherent limitations. The primary issue is the inability to allow for variable or dynamic charge transfer, and the overemphasis on equilibrium configurations, which can sometimes lead to large errors. For example, the MEAM potential fitted to equilibrium structures28 performs poorly in predicting the relative energies of far-from equilibrium Al-Cu alloy configurations sampled using rigorous evolutionary searches (Figure 2a).

6


Page 6 of 47

Page 7 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Limitations of existing force-fields for molecular simulations. (a) Formation enthalpy predicted using MEAM and DFT for the Al-Cu alloy system. Note that the MEAM errors are very large for the sampled non-equilibrium configurations especially for alloy compositions that were not included in the original fitting of Al-Cu. (b) Dynamical evolution of thermosensitive polymer undergoing coil-to-globule transition obtained from four different FFs (AMBER29, OPLS-AA30, CHARMM31, and PCFF32). Dynamics of conformational change of a thermosensitive polymer predicted by four popular FFs is seen to be widely different. Radius of gyration, RG, is a measure of chain conformation; RG ~ 18 Å signifies coil state whereas RG ~ 10 Å is the collapsed state. To describe re-distribution of atomic charges during the dynamic evolution of reactive interfaces, new formalisms have recently emerged, which integrate bond-order concepts with electronegativity equalization schemes. Two of such reactive force fields have gained a lot of popularity over the last decade, namely, ReaxFF (reactive FF) method developed by van Duin and co-workers, 33-34 and Charge-optimized many-body (COMB) potentials developed by Sinnott, Phillpot, and coworkers.35-36 In the ReaxFF formalism33-34, short-range bonded interactions are computed using instantaneous bond-orders around each unique pair of atoms; the bond-order for each atomic pair is strongly influenced by the local atomic, and accounts for multi-body effects. Instantaneous atomic charges are derived using an electronegativity equalization (EEM) or charge-equilibration (QEq) scheme37. Such a bond-order formalism coupled with dynamic charge re-distribution allows for an accurate description of bond formation and dissociation. The long-range non-bonded interactions are accounted for using van der Waals, and electrostatic terms. On the

7



other hand, COMB potentials are built upon the well-established Tersoff-formalism, while using chargeequilibration schemes to describe charge transfer. Regardless of the mathematical treatment employed by the different flavors of reactive force fields, they often possess a large number of independent parameters, most of which can neither be directly computed (via quantum chemistry) nor empirically measured. These parameters are commonly obtained by least-squares fitting against a set of structural, mechanical, and thermodynamic properties computed largely near equilibrium (e.g., cohesive energy, elastic constants, and structural properties of crystalline phases, or gas-phase configurations).38 For instance, the ReaxFF parameter set for CaCO3-water system contains 468 independent parameters that were obtained by fitting to the properties of four bulk condensed phases and one gas-phase anion. The over-reliance of the fitting protocol on equilibrium configurations often results in poor performance of even these sophisticated reactive potentials; e.g., the ReaxFF for CaCO3-water (with 468 independent parameters) fails at reproducing the structure of calcite/water interface as measured by xray reflectivity.39 In addition to static structural properties at reactive interfaces, existing FFs (depending on their training procedures) also exhibit a wide variability in describing dynamical phenomena. Figure 2b shows differences in dynamical evolution of four popular FFs for a model thermosensitive polymer undergoing thermally induced coil-to-globule transition40-42. While the thermodynamic end-points, i.e., coil and collapsed phases are reasonably well described by the four FFs, the dynamical evolution between the two phases is clearly different. In essence, the shortcomings of existing FFs can be broadly attributed to the following factors: 1. In a typical force field development, the training set is not balanced over the potential energy space. Commonly, equilibrium configurations are over-emphasized, while numerous far-from-equilibrium configurations and transition states are ignored. This often results in poor description of atomic-scale dynamics, which tends to be associated with non-equilibrium configurations. 2. Local optimization, and least squares are typical methods of choice for fitting the independent parameters. The success of these methods is largely controlled by initial guesses (or human intuition); in general, these methods are prone to convergence issues, and over-fitting especially in the large

8


Page 8 of 47

Page 9 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


parameter space required for reactive force field formulations, 3. Often there is substantial lack of quantitative cross-validation against known properties not used in the fitting procedure. 4. Many FFs employ fixed charges making it difficult to model nonstoichiometric defect formation, defect distributions charge and defect migration across interfaces and capture finite-size effects. 5. Most fitting procedures use single objective function, (i.e., a weighted sum of errors in prediction of target properties) to describe the quality of a given parameter set. The weights are arbitrarily chosen, primarily driven by intuition (or experience). This often results in inefficiencies given the large parameter space, and thwarts identification of good parameter sets within reasonable time.

3. Machine learning inter-atomic potentials to bridge the electronic and atomistic scales In this section, we discuss our recent efforts in the development of new, first-principles based, and robust inter-atomic potentials for accurate simulations of dynamical processes at reactive interfaces and low dimensional systems, such as clusters and molecules. As shown in Figure 3, our procedure involves (1) defining or selecting a functional form, (2) constructing an extensive training data set from electronic structure calculations, (3) formulating a fitting procedure, and (4) implementing these algorithms on HPCs. This new class of potentials enables accurate prediction of inter-atomic forces, and thereby, allows highfidelity dynamical and statistical simulations of low dimensional systems, properties and functionalities of new materials, as well as pathways and mechanisms of their synthesis and assembly. Using the developed methodology on computational facilities at Argonne, we have explored a range of target systems (Section 4), and gained crucial insights into the atomic scale processes in clusters and interfaces relevant to energy applications. Next, we describe each of the steps involved in developing an accurate, robust, and efficient FF.

9



Figure 3. An overview of our machine learning framework - BLAST (Bridging Length/Timescales via Atomistic Simulation Toolkit). Schematic illustrates the BLAST workflow used to train interatomic potentials from first-principles calculations.

3.1. Model Selection: The capability, transferability, and limitations of FFs are, to a great extent, dictated by their mathematical functional form describing the inter-atomic interactions. In popular MD codes (e.g. LAMMPS) for materials simulations, a wide range of pre-defined functional forms are available43; however, it is crucial to select one carefully such that it is suited for investigating the problem of interest, while being computationally as efficient as possible. For instance, simple pairwise potentials, e.g., Morse or Lennard-Jones forms, cannot describe directional bonding in covalent materials owing to the absence of energy contributions from 3-body and higher order interactions in their functional form 44. Apart from the material being studied, the choice of functional form is also strongly dependent on the phenomenon being explored. One of the example systems described below, namely Au clusters, provides a clear illustration of this dependence. As described in the section 4.1, Au clusters in the sub-nanometer size regime exhibit diverse structural motifs (e.g., planes, tubes, hollow cages, space-filled) encompassing a wide range of atomic coordination environments, with implications in catalysis applications. Although, conventionally, force-fields employing functions of local atomic density (based on pairwise distances), such as EAM, have been successfully employed to study the structural, elastic, and thermodynamic properties of bulk gold 45,

10


Page 10 of 47

Page 11 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the functional forms of these potentials are ill-equipped to describe the diverse geometries exhibited by Au nanoclusters. Indeed, even an exhaustive search for EAM parameters (fitting against a large training set containing Au clusters) could not yield a set that could reliably describe the structure and energetics of Au clusters in various geometric confirmations. This search failed because spherically symmetric potentials like EAM cannot predict open structures, (e.g., planar or cage-like) to be an energy minimum, and consequently, cannot be used to investigate size-dependent changes in structural motifs in Au clusters. In fact, such a study would require a FF that can account for bond-directionality (e.g., the Tersoff formalism23), and is contrary to conventional wisdom to employ EAM-type models for metals. Evidently, a careful assessment of the existing functional forms affirming its suitability for both the problem and material of interest represents the first key step in atomistic modeling. 3.2. Training Data Set: An ideal force-field should be able to describe the material structure and properties across all locations in the free energy landscape. To achieve this, a force-field should be trained with data that contains an inclusive set of structures, energies, elastic constants, surface energies, defect properties, and other properties relevant to the system of interest. We create a large training dataset encompassing wide range of bulk and surface structures, as well as cluster configurations, such that it provides an ample sampling of all possible atomic environments/coordination likely to be encountered in dynamical simulations. For bulk phases, we typically include lattice parameters, internal coordinates, elastic constants, and cohesive energies. Addition of surface configurations in the training set enables us to capture energetic differences between various coordination environments, model broken bonds, as well as describe the electrostatic interactions; this, in turn, greatly enhances the FF’s capability to describe defect properties. Interestingly, a random sampling approach to generate training configurations is problematic, since it tends to favor high energy structures (Figure 4). To address this issue, we have implemented an evolutionary structural search method based on GA46 and coupled it with DFT calculations17. This technique yields a uniform energetic distribution of training configurations, and provides a good sampling of both near- and far-from equilibrium configurations (Figure 4). The idea of GA is inspired by the idea of natural selection or survival of the fittest. It starts with a set of randomly generated configurations, which is defined

11



as the population. Next, the individuals (configurations) in the population are subjected to genetic operations, such as (a) crossover, and (b) mutation to obtain child structures, which constitutes evolving a generation. Crossover consists of selecting different motifs from each participating parent, and juxtaposing them together. Mutation, on the other hand, involves introducing a random structural distortion to the chosen parent. All the individuals (both parents and children) are then ranked according to their DFT energies; the ones with lower energies to defined to possess higher fitness. At the end of each generation, the structures with lowest fitness are eliminated such that the population size remains constant. Our search method is distinct from existing GA codes primarily in that it uses a spatial decomposition scheme17 for crossover with flexible number of parents, rather than the standard two-parent schemes. The number of parents can be increased until each parent contains the smallest meaningful atomic neighborhood. This approach speeds up the convergence of the GA run, and helps maintain diversity of the pool. In addition to crossover and mutation, other genetic operations, such as parent cloning and parent exchange are also applied. The rates to apply these operations are adjusted on the fly depending on the energy differences between generations. Consequently, a stagnant gene pool and early convergence during search is avoided. We applied this approach to gold nanoclusters (see Section 4.1) and successfully predicted the ground state structures for different cluster sizes, and built a good training set with unprecedented sampling of various atomic environments.

b)

a)

Figure 4. Comparison between the energy sampling of randomly generated (a) and evolutionary algorithm based (b) Au13 nanoclusters. As seen from the range of energetic distributions of the sampled

12


Page 12 of 47

Page 13 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


configurations, a GA-based sampling provides a more diverse set of training data compared to random sampling. Reproduced with permission from Ref (14). 3.3. Two stage optimization (Global optimization followed by local optimization) Optimizing forcefield parameters for complex functional forms requires a systematic training procedure, which allows efficient sampling of the entire parameter space. Traditionally, the procedure used to determine the optimal parameter set is almost exclusively the least squares (LS) method. There have been alternative fitting schemes such as the minimal maximum error (MME) procedure and the interpolative moving least squares (IMLS)47 technique, both of which are superior to the LS methodology in a number of respects. In these existing fitting procedures, the automation is realized through a built-in monitoring of a “local” (as distinct from the LS “global”) error indicator that triggers a response mechanism in the fitting procedure, which makes it difficult to fit to large and diverse training data sets. To explore the entire parameter space, an unknown number of starting guesses needs to be employed. Consequently, it is not surprising that most popular FFs that were derived using local optimization procedures often have poor predictive power. ML techniques48-52 are beginning to emerge as an alternative approach to construct very accurate and efficient potentials by rigorously training against an extensive training dataset derived from high-level electronic structure calculations. Global optimization methods provide efficient ways of exploring the parameter space, and are wellsuited to for force-field parameterization. Therefore, we employ evolutionary algorithms namely GA, particle swarm or differential evolution to perform an unbiased global search of the parameter space by maintaining a diverse population. We find that these techniques are quite successful at identifying promising regions of interest in the parameter landscape, much faster than random sampling or traditional “trial-and-error” methods 14, 17, 53. To further improve the efficiency of the search, we combine global search methods with local optimization; this ensures that all the sampled parameter sets are local minimum in the parameter landscape. For local optimization, we employ the commonly used simplex54 and LevenbergMarquardt55 methods that are known to perform reasonably well for many non-linear least-squares

13



problems. In addition to these common methods, we have also tested more sophisticated derivative-free methods for non-linear least squares problems as implemented in POUNDERS (Practical Optimization Using No Derivatives for sums of Squares) Interpolation in Trust-Regions)

57

56

and ORBIT (Optimization by Radial Basis Function

for the first time in FF fitting. These methods are overall less efficient

than global optimizers such as GA, but given good initial guesses, they can locate solutions that are better. For instance, we have compared various optimization algorithms including POUNDERS on a model IrO2 system58. For this model system, we do find that POUNDERS seems to perform better than GA i.e. requires fewer evaluations to converge to the optimal parameter set. Details are in Ref (58). More detailed studies on diverse classes of materials (and functional forms) are needed to evaluate the potential of methods such as POUNDERS for force-field optimization. Regardless of the global/local optimization methods employed, the goodness of fit (or the predictive capability of the FF) is measured by a weighted sum of errors in predictions of a set of target observables. As aforementioned, such weights are usually chosen at the whim of the developer, and in many instances, are a hindrance to achieve the best possible set. One way of circumventing the spurious effects of arbitrary weighting is to use multi-objective optimization. Using this route leads to a set of Pareto-optimal solutions59 (parameter-sets), and enables selection of FF parameters according to target applications (Figure 5). (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 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 as in the case of ZrN 15). Overall, our ML framework allows users to explore and use combinations of different supervised learning algorithms and optimization schemes including local, global, single- and multi-objective evolutionary algorithms (GA, particle swarm, and differential evolution) and linear/sparse regression methods (lasso) to train their FFs.

14


Page 14 of 47

ns ti o ra ne ge

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Average Error in Elastic Constants (1))

Page 15 of 47

Pa re

to fro nt Average Error in Energies (2))

Figure 5. Multi-objective evolutionary optimization of a potential model. Schematic shows the evolution of a typical Pareto front with generations during a typical multi-objective evolutionary optimization run for two objective functions, namely, the sum of squared errors in the prediction of structural and elastic properties Δ(1), and that for configurational energies Δ(2). 3.4. Automated Workflow: The general workflow of our ML framework (BLAST) for training a FF is illustrated in Figure 6. The idea is that a user will specify the form of the empirical potential, choose any of the ML or optimization algorithms, provide some general input parameters such as number of configurations, crystal type and symmetry etc., and the desired or target output. A backend (e.g., Swift60) will then launch the various user specified ML and optimization algorithms in parallel and directly port the input/output to/from either VASP for DFT or LAMMPS for MD. A user can employ a combination of global optimization schemes (GA, particle swarm, differential evolution, or neural network) as well as local optimization techniques (simplex, Levenberg-Marquardt) to generate the optimal set. Briefly, a user first generates a training set with necessary properties, e.g., structures, energies, elastic constants, defect energies

15



etc.; these properties are primarily evaluated by DFT calculations, but there will be provisions to include experimental data as well. After constructing the training set, a user will perform a global optimization (e.g., GA) of the FF parameters. In the case of GA optimization, first a population Np (~ tens to several hundreds) of parameter sets are generated randomly. For each parameter set (members) in the population, the target properties are evaluated using popular MD packages (e.g. LAMMPS). Next, our ML workflow will compute the weighted sum of errors in prediction for every member, which is defined as its objective function (OF). The members are then sorted in the increasing order of their OF – a member with lower value of OF is better (fitter) than those with high OF values. A certain number of m best members (low OF) among the total Np are chosen, which are then subjected to genetic operations (e.g., crossover, mutation etc.). These operations result in “offspring” members containing the traits (binary segments) from their parents. Out of the Np members, and the offsprings, only the best Np sets are retained for the next generation. Such an optimization routine ensures that only good parameter sets survive after each generation; upon repeating this workflow for sufficient generations, we sample diverse regions in the parameter space before the run converges, i.e., the values of OF for at least Np/4 members in the population are below the prescribed tolerance. Typically, tens of such GA runs starting with different random population are necessary to scour the parameter space; we choose ~10 best parameter sets from each of these GA runs, which give similar low OF values. Then, using these sets as initial guesses, we perform local optimization using simplex or Levenberg-Marquardt techniques. The optimized parameter sets further cross-validated to obtain the best one.

16


Page 16 of 47

Page 17 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Automated machine learning workflow to train first principles based forcefields for molecular dynamics simulations.

4. Case studies (model selection, nature of the training data, training procedure) 4.1. Metal clusters (planar to globular transition in gold) Nanoscale effects in gold result in exceptional chemical, optical, and electronic properties. This makes Au nanoclusters promising for applications in various technological sectors, including optoelectronics, biorecognition, and catalysis61-62. Their enticing potential has fueled numerous theoretical and experimental investigations 63-66; however, a fundamental understanding of the size dependent structural transformation in gold clusters remains a topic of intense interest. For instance, there is an open debate on the critical value of cluster size (i.e., number of atoms) beyond which globular isomers become energetically preferable over planar ones, with reported values ranging from 7 to 15 63-64, 66-67. Furthermore, since many isomers at a given

17



Page 18 of 47

cluster size are energetically close to each other (~20 meV/atom),63 it is possible that they may undergo structural transitions under the influence of external stimuli, e.g., temperature fluctuations; the knowledge of such transitions and the associated mechanisms is still in its infancy. For Au, there are four prominent FFs currently available in the literature, namely, Sutton-Chen (SC)68, Gupta69-70, Embedded Atom Method (EAM)

71,

and Reactive force field (ReaxFF)

72.

These

traditional FFs have, however, been parameterized by fitting against structural and thermodynamic properties of bulk face-centered cubic (FCC) configurations of gold. To test these FFs, we generated a test set of ~1000 Au13 (size: 13 atoms) nano-clusters and compared the relative energetics predicted by these potentials with those obtained from DFT calculations. All the traditional FFs describe Au nanoclusters poorly, as indicated by the large average error compared to DFT energies as well as well as significant scatter (Figure 7). Evidently, currently available classical FFs fail to describe the structural features and the diverse geometries of gold at the nanoscale; these deficiencies could possibly be related to (a) the use of spherically symmetric functional forms that cannot favor planar clusters, and (b) the overemphasis of nearequilibrium structures in the fitting procedure. In Ref (14), we used our ML framework to develop an interatomic potential capable of describing the diverse size dependent motifs of gold (planar, globular, pyramidal, icosahedral, cage etc.) as well as bulk. 4.1.1. Model selection: We have developed a potential formalism termed hybrid bond order potential (HyBOP), wherein the total potential energy is composed of partial energy contributions arising from both short-range and long-range atomic interactions. The potential also employs a scaling function that captures the size and distance dependence of the long-range contribution. In other words, for small sized clusters, the interactions are dominated by short-range interactions and the long-range interactions become increasingly significant as we approach the bulk. The HyBOP contains 17 independent parameters that need to be optimized. 4.1.2. Training data: Our training data include configurations that predict both near equilibrium and far from equilibrium structures, i.e., a continuous range of energies (vs. structures) from high to low. We use the structural GA as described in Section 3.2 and coupled it with DFT calculations to generate a diverse

18


Page 19 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


training set. A comparison between the DFT energy sampling of randomly generated, and GA-based Au13 nanoclusters is shown in Figure 7b. It is seen that randomly created nanoclusters tend to emphasize a narrower phase space far from equilibrium. With GA, a better sampling of the potential energy surface is obtained as shown by the broader range of energies spanned by the configurations. Overall, our training data set included ~1200 configurations of Au13 planar and globular structures as well as their energies. Bulk properties such as lattice constant, elastic constants, internal coordinates and cohesive energy were also included in the training data set. The structures used in the training set are sufficiently diverse and amply represent different parts of the potential energy surface, as detailed in the section 3.2. 4.1.3. Parameter optimization: We use our ML framework to optimize the HyBOP parameters. Global optimization is performed using the single/multi objective genetic algorithm toolbox73. The optimization procedure is initiated by randomly generating a population of Np = 60 parameter sets. For each HyBOP parameter set (i.e., 17 in total), we used the MD simulation package LAMMPS to compute the energies for all structures in the training set. The objective function is defined as the weighted average of the difference between HyBOP-predicted and DFT energies for the structures in the training set. We employ a single objective and the weights for the various observables (elastic, internal coordinates, lattice constants etc.) are chosen such that the total weighted sum of squared error for each property type is approximately equal. We choose ~ 10 different parameters from the GA run and run a local optimization using the simplex method until the error in successive steps is less than the tolerance (set to be 10–14).

19



Figure 7. A hybrid bond order potential (HyBOP) model developed in Ref (14) for capturing the diverse geometries of gold clusters. (a) Cohesive energy of 20 near global minimum Au13 configurations evaluated using DFT, Tersoff-inspired BOP, and our hybrid bond-order potential (HyBOP). (b) HyBOP predicts planar P1 structure for Au13 to be the global energy minimum, which agrees with DFT prediction. Other FFs incorrectly predict globular Ih to be the most stable structure for Au13. (c) Global energy minimum configurations of Au nanoclusters at different sizes as predicted by HyBOP. (d, e) Comparison of the DFT (PBE-GGA) predicted cohesive energies for 345 planar and 901 globular Au13 nanoclusters against HyBOP predictions and those by various popular FFs (Sutton–Chen, reactive force field, embedded atom method, and Tersoff-type BOP). Each panel provides the mean absolute errors (δ) in the predicted energies as compared to the DFT values. Reproduced with permission from Ref (14).

4.1.4. Model validation: Our HyBOP parameters are extremely robust and predict planar global-energy minimum configurations for smaller sizes of 4–12 atoms as well in excellent agreement with DFT overcoming the limitations of existing FFs (Figure 7a). Briefly, the ML trained HyBOP parameters accurately predict (i) global energy-minimum configurations for various cluster sizes as in Figure. 7b, (ii) evolution of structural motifs (planes, hollow cages, space-filled) with cluster size as shown in Figure 7c, and (iii) critical size for transition from planar to globular structures in Figure 7b. In addition, these BOP

20


Page 20 of 47

Page 21 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


parameters can capture the transition/intermediate states, which arise during various dynamical phenomena in Au nanoclusters. A few examples of such phenomena include formation of Au clusters either atom-byatom or via combination of smaller clusters, phase transition between different isomers at a given cluster size, and melting. 4.2. 2-D material such as Stanene (Thermal transport) The recent successful synthesis of 2-D tin or stanene on a Bi2Te3(111) substrate via molecular beam epitaxy (MBE) has opened up new avenues for exploiting the exotic electronic, thermal and mechanical properties of this fascinating 2-D material74. Despite these encouraging prospects, an atomistic-scale understanding of the finite temperature effects on structure and transport properties has not been achieved. A major roadblock has been the inability to perform finite temperature as well as long time and length scale simulations due to a lack of accurate FFs for emerging and interesting 2-D systems such as stanene. We employed our ML framework to develop the first interatomic potential for stanene (Figure 8a).16 4.2.1. Model selection: We adopt the bond order formalism outlined by Tersoff27 as the functional form of this force field. The Tersoff formalism captures many-body interactions, bond breaking and formation, and is computationally inexpensive compared to other bond order potentials28. Tersoff potentials have also been successful in describing atomistic interactions and thermal transport in other analogous 2D materials such as graphene29 and hexagonal boron nitride30 and hence was chosen to model stanene. 4.2.2. Training data: The potential parameters of this functional form were trained against an ab-initio derived dataset that included the equation of state, cohesive energy, lattice constants, buckling, and in-plane stiffness calculated using DFT. We also fit to the DFT phonon dispersion (Figure 8 b,c) to ensure that the potential accurately reproduces the phonon transport in the material as well as mechanical properties. 4.2.3. Parameter optimization: We use a two-stage optimization process to obtain the best-fit set of parameters to the training set. In the first stage, GA is used to perform a global parameter search. Ten isolated populations of 200 individuals (a total of 2000 individuals) are given random initial values for the 12 adjustable Tersoff parameters. We evaluate the fitness of each individual by comparing the predicted

21



value from MD vs. DFT calculations. Individuals are ranked by their fitness function before performing genetic operations i.e. crossover and mutation. The independent GA runs continue until the fitness of their best individuals have converged. We choose the parameter set with the lowest objective function among all GA runs and use these as input to a local optimization procedure using Nelder-Mead simplex. The best set obtained from this procedure captures the density of states and the phonon dispersion, which are in excellent agreement with ab initio calculations (Figure 8 b,c).

a)

c)

b)

ML-BOP

Figure 8. A machine learned bond order potential (BOP) derived from an ab initio derived training data set. (a) An MD simulation snapshot showing rippling in stanene. Comparison of the (b) phonon dispersions and (c) phonon density of states (DOS) of stanene computed from DFT calculations and the ML trained BOP. Reproduced with permission from Ref (16).

4.2.4. Model validation: As a representative case study, we use this ML trained model to perform MD simulations to study stanene’s structure and temperature-dependent thermal conductivity (Figure 8a). In comparison to 2-D materials such as graphene, stanene is highly rippled owing to its low in-plane stiffness (stanene: ~ 25 N/m; graphene: ~ 480 N/m). Hence, it shows stronger temperature dependence compared to that in graphene. Simulations performed using this ML model show that stanene nanostructures have significantly lower thermal conductivity compared to graphene based structures owing to their softness (i.e., low phonon group velocities) and high degree of anharmonicity. This ML model facilitates the exploration of stanene based low dimensional heterostructures for thermoelectric and thermal management applications.

22


Page 22 of 47

Page 23 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.3. Structure and Dynamics of Heterointerfaces (mixed metallic-covalent systems) Nanoscale hybrid architectures composed of metals and organic materials provide a unique opportunity to design next-generation of multifunctional materials (often with conflicting/exotic set of physical properties) for applications in flexible electronics, catalysis, and energy storage. Despite their great promises, a fundamental understanding of the inter-relationships between structure, morphology, atomic scale dynamics, chemistry, and physical properties of such mixed metallic-covalent systems has not been achieved. Large-scale MD simulations represent a viable route to attain such knowledge; however, a lack of robust, efficient, accurate, and predictive FFs remains a major impediment in the success of these simulations. The chemical complexity of the metal–organic interfaces (metallic in metal-rich regions and covalent in organic-rich regions) requires development of FFs that accurately describe the atomic interactions in disparate environments. It is often difficult to accurately yet efficiently treat such diverse chemical bonding characteristics within a single potential formalism. In the literature, it is quite common to use different models to resolve such problem (e.g., EAM for metallic, Tersoff for covalent and Buckingham for oxide interactions). Such an approach often leads to inaccuracies, in particular, at the interfacial region that exhibits mixed metal-covalent characteristics. In this context, ML algorithms allowed us to introduce a first-principles based bond-order potential (BOP) that utilizes a single framework that accurately describing metallic, covalent, and mixed ionic-covalent interactions. We choose a representative metal-organic system, namely cobalt-carbon to exemplify the suitability of our ML framework.7 4.3.1. Model selection: To model the interactions in a Co-C system, we choose the Tersoff-Brenner formalism44. This formalism has been successful in describing a diverse range of materials, including metals, ionic and covalent solids, as well as those with mixed bonding characteristics. The formalism is well-known for its accuracy in describing the structure, thermodynamics, mechanical, and surface properties of various forms of carbon.

23



4.3.2. Training data: We use DFT to derive our training set which contains (a) lattice parameters, (b) equations of state (c) elastic constants, and (d) cohesive energies of different crystalline phases of cobalt as well as Co2C. For cobalt, we consider four stable and metastable polymorphs, namely hexagonal close packed (hcp), face-centered cubic (fcc), body-centered cubic (bcc), and diamond cubic. For cobalt carbides, we employ a metastable Co2C phase with an orthorhombic symmetry (space group: Pmnm) and we retain the original Tersoff parameters optimized for carbon. 4.3.3. Parameter optimization: Our ML framework employs a combination of global optimization via GA and local minimization (simplex) to optimize the independent BOP parameters, for Co–Co, and Co–C interactions (i.e., a total of 22 adjustable parameters). Our ML trained model reproduces cohesive energies, lattice parameters and elastic constants for five different polymorphs of Co as well as structural and thermodynamic properties of the carbide phases, in excellent accordance with DFT calculations and experiments.

Figure 9. Multi-million atom simulations of a Co-C hybrid system modeled with ML trained potentials. MD simulations of self-assembly of Co515 clusters and C60 at 800 K. (a–c) Final configurations of the porous heterostructures at various C/Co ratios obtained from our MD simulations. The Co and C atoms are depicted as blue and grey spheres, respectively. (d–f) Chain/branched structures of Co515 clusters identified using a clustering algorithm. Individual clusters are labeled by different colors. Reproduced from Ref (7) with permission from Royal Society of Chemistry.

24


Page 24 of 47

Page 25 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.3.4. Model validation: We use our ML model for Co-C to perform massively parallel multi-million atom simulations on leadership class computing resources to investigate the atomic scale dynamics associated with self-assembly of porous Co/C hetero-structures with applications in flexible electronics (Figure 9). Our work represents an important step towards establishing a systematic ML protocol for training ab-initio based FFs for nanoscale hetero-interfaces. This example highlights how ML strategies can be used to (i) develop a single first principles-based, computationally robust but inexpensive framework to describe complex interactions in nanoscale hybrids and across diverse interfaces, and (ii) provide new insights into a wide range of bulk/surface phenomena as well as transport processes including segregation during their self-assembly. This approach is also generalizable to a broad class of organic/inorganic hybrid nanoarchitectures. 4.4. Nitrides (Multi-objective optimization) Zirconium nitride (ZrN) is a representative of an important class of nitride materials that exhibits exceptional mechanical, chemical, and electrical properties. This makes the material attractive for a wide range of technological applications, including wear-resistant coating, protection from corrosion, cutting/shaping tools, and nuclear breeder reactors. Despite its broad use, an atomic-scale understanding of the superior performance of ZrN, and its response to external stimuli, e.g., temperature, applied strain, etc., has not been achieved. For heterosystems such as nitrides, it is often difficult to develop interatomic potential models that accurately describe the interactions between Zr and N atoms. Moreover, it is often difficult to treat surface and bulk properties within a single objective function. In Ref (15), we demonstrated the effectiveness of multi-objective global optimization strategies to alleviate the problems associated with assigning weights to the surface vs. bulk properties during force field fitting. 4.4.1. Model selection: For ZrN, we choose the modified embedded-atom method (MEAM) interatomic potential since it can simultaneously describe a wide range of crystal structures (body-centered cubic (bcc), face-centered cubic (fcc), hexagonal close-packed (hcp), diamond and even gaseous elements).

25



4.4.2. Training data: Our training set is derived from DFT and contained (a) equation of state, (b) elastic constants, and (c) energies of low index surfaces [(001), (110), and (111)] of rocksalt ZrN, as well as lattice parameters and heat of formation of another commonly observed compound, Zr3N4. 4.4.3. Parameter optimization: In all global and local optimization schemes, the goodness of fit, and hence, the predictive capability of the force field, depends on the weighting scheme selected between different observables. The selection of weights is difficult because the errors often have different orders of magnitude and originate from different number of measurements. In order to obtain a more reliable and accurate fit, a more systematic way of weighting is required. Our framework circumvents the problem of weight selection by using a multi-objective GA (MOGA) optimization algorithm73. This gives us a set of solutions at the Pareto front and enables selection of force field parameters according to user desired properties (Figure 9a).

Figure 10. Multi-objective evolutionary optimization of an MEAM potential for ZrN system. (a) Evolution of the Pareto front with generation number g during a typical multi-objective evolutionary optimization run for two objective functions, namely, the sum of squared errors in the prediction of structural and elastic properties Δ(1), and that for surface energies Δ(2). The final optimized parameter set obtained in this study is represented by the orange star. (b) Atomistic representation of brittle fracture of nanopillar ZrN. The corresponding strains (with reference to the unstrained length of ZrN) are shown for each atomistic snapshot. The MD simulation for the above representative case was performed at a strain rate of 109 s–1. The atoms are colored by the per-atom potential energy. Unlike the bulk ZrN, significant amorphization is seen near the crack plane. After brittle fracture (at ∼10%), the crack propagation proceeds gradually with the formation of a narrow-connected neck and eventual separation of the cracked ZrN pillars

26


Page 26 of 47

Page 27 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


at ∼50% strain. Reproduced with permission from Ref (15). The evolution of the Pareto front is shown in Figure 10a as a function of generations during our multi-objective GA (MOGA) run. The Pareto front at each generation provides a collection of nondominated MEAM parameter sets in the objectives Δ(1) and Δ(2). Initially (g = 0), the population of MEAM parameter sets is randomly chosen within a physically allowable domain in the parameter space. This results in parameter sets that exhibit extremely high objective function values (poor suitability). With increasing MOGA generations, fitter MEAM parameter sets (with lower objective values) are sampled, which causes the Pareto front to advance gradually toward lower values of Δ(1) and Δ(2). Eventually, MOGA run converges and we notice that the Pareto front stops advancing at ∼400 generations; these represent the population of Np best candidate parameter sets. From the final generation, we pick an MEAM potential accurately predicts structure, thermodynamics, energetic ordering of bulk polymorphs, as well as elastic and surface properties of ZrN, in excellent agreement with DFT calculations and experiments. 4.4.4. Model validation: As a representative test case, we utilize our MEAM potential to study the atomic scale mechanisms underlying fracture of bulk and nanopillar ZrN under applied uniaxial strains, as well as the impact of strain rate on their mechanical behavior (Figure 10b). Our ML trained MEAM potential model represents the first FF for ZrN which captures both bulk and surface properties, making it suitable for studying structural and dynamical response to external stimuli such as, for example, temperature and applied strain. This example highlights the advantage of multi-objective evolutionary optimization strategies in developing first principles-based, computationally robust but inexpensive, models to investigate a wide range of bulk/surface phenomena in materials models.

4.5. Variable charge model for iridium oxide All the examples highlighted in sections 4.1−4.4 represent non-reactive systems. In many systems that involve charged entities like ions, there is a clear need to explicitly account for polarizability and charge transfer effects. Variable charge methods compute the electrostatic interactions dynamically based on the

27



connectivity and geometry of the atomic configurations. This allows to capture the chemical environment dependent properties to be estimated more accurately such as surfaces and point defects. Iridium oxide (IrO2) represents a prototypical system that demonstrates exciting electrical, electrochemical and spintronic properties. IrO2 is one of the most efficient pH sensors and IrO2 nanoparticles are currently amongst the best electrocatalyst for water splitting reaction in producing hydrogen from sunlight owing to their high activity, and exceptional stability 75-77. The catalytic properties of these clusters are sensitive to their size and structure

77.

Most of the catalytic and surface chemical

properties depend on the size, morphology, and structure of IrO2 surfaces, and nanoclusters. Therefore, a deeper understanding of the structure and dynamics of IrO2 surfaces and nanoclusters at the atomic scale is essential to develop improved solar fuel catalysts and electrodes made from IrO2. Interestingly, however, despite its tremendous technological importance, that there was a lack of an empirical potential model for simulating systems such as Ir-O. To address this need, we employed our ML framework to develop a variable charge potential model in Ref (13) to probe structure and surface stability of IrO2 (Figure 11). The key steps involved include: 4.5.1. Model selection: We choose a variable charge potential (MS-Q) function that couples a Morse functional form with a Coulomb term. The charge equilibration method (QEq) is used to compute the charges on the atoms and compute the electrostatic contribution to the total energy 37. 4.5.2. Training data: We generate an extensive training data set using DFT that includes structural, mechanical, and thermodynamic properties for various experimental and hypothetical crystal structures of IrO2 bulk structures. We considered ground state rutile (P42/mnm) structure78 and high pressure (>15 GPa) pyrite (Pa-3) phase79 of IrO2. All DFT calculations included Hubbard U=1.0 eV correction and spin orbit coupling (SOC) in PBE functional 80. The cohesive energies (Ecoh) were adjusted by 1.198 eV per O2, to correct the error in calculation of O2 reaction energies with GGA-PBE functional 81. In addition to these experimental structures, we also computed energies and structural properties for the hypothetical anatase (I41/amd), columbite (Pbcn) and brookite (Pbca) polymorphs of IrO2. Surface energies of (110) and (100) surfaces of rutile IrO2 were also included in the training set. The use of different polymorphs and two

28


Page 28 of 47

Page 29 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


surfaces in the training set provided a variety of local environments in order to capture variable charge effects in the fitted force field. 4.5.3. Parameter optimization: The training of Morse potential parameters De, α, and r0 for Ir-Ir, Ir-O and O-O pairs and QEq parameters χ, and J for Ir and O elements requires an efficient optimization scheme. We employed a combination of global optimization performed using GA followed by the simplex method. An unbiased, and global search performed by maintaining a diverse population allows us to potentially identify good regions of interest. This represents the starting point for the simplex method, which is performed to identify the best set of parameters. The computed lattice parameters for the various Ir-O polymorphs are within 2% error compared to respective DFT values. The discrepancy in predicted elastic constants with MS-Q are less than 20% compared to DFT values and the bulk modulus is predicted to within 5%. The ML trained MSQ potential also adequately captures the electrostatics despite the fact that atomic charges were not included in the training set. We obtained an excellent agreement between MS-Q and DFT Bader charges. For the rutile phase, the respective Ir and O charges of +1.685, and -0.843 matched very well with DFT Bader charges of +1.658 and -0.829, respectively. For other polymorphs, the MS-Q and Bader charges also matched well with little discrepancy in the amount of charges relative to rutile phase.

Figure 11. A first-principles-based, variable-charge force field that accurately predict bulk and nanoscale structural and thermodynamic properties of IrO2. Reproduced from Ref (13) with permission with permission from Royal Society of Chemistry.

29



For catalysis applications at nanoscale, predicting the surface properties accurately is important. To test the derived potential in predicting surface properties of IrO2; we compute the surface energies (γ) of low index surfaces of rutile IrO2 namely, (110), (101), (100) and (001) with the MS-Q potential. We compared these with our DFT results and earlier revPBE

82

calculations. We find a strong agreement

between the MS-Q γ values and the γ values obtained with DFT. In both methods IrO2 (110) surface was the most stable, and the relative stability of surfaces followed:  (110)   (101)   (100)   (001) , which matches well with the results for other rutile structures 83. The similarity of γ values between MS-Q and DFT is evident in the constructed the Wulff plot 84 shown in Figure 12a,b. This indicates that both MS-Q and DFT predict a similar surface structure and crystallite shape for the IrO2. Note that the MS-Q predicts slightly larger area coverage for (100) compared to DFT.

c)

Figure 12. Predictive power of the variable charge Morse potential model trained using the ML framework shown in Figure 3. Wulff shape of rutile IrO2 created using a) DFT and b) MS-Q computed low index surface energies. c) Adsorption sites considered on IrO2 nanocluster. Large atoms are Ir and small atoms are O. The adsorbed oxygen atoms at the corresponding site were numbered and colored differently. Reproduced from Ref (13) with permission from Royal Society of Chemistry.

30


Page 30 of 47

Page 31 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.5.4. Model validation: To understand the water splitting and oxygen evolution reaction (OER) activity on IrO2 surfaces, we constructed a 2.5 nm sized nanocluster of IrO2 based on the Wulff construction shown

Ead  EslabO  Eslab  EO in that composed of (110), (101) and (100) facets and studied the O adsorption energetics. The adsorption energy (Ead) of a single O atom on rutile IrO2 was calculated using where Eslab+O is the total energy of the (110) surface slab with the adsorbed O atom and Eslab is the total energy of the pristine (110) surface slab. Here, we assumed that the dissociation of O2 molecule reaction which has a large energy barrier to atomic O has already happened so that EO is the total energy of the isolated O atom, while in MS-Q EO is zero since the zero of the energy is set at the isolated atoms limit. The nanocrystal constructed has 785 atoms, i.e. 524 O and 261 Ir atoms. The nanostructure is near stoichiometric with an O:Ir ratio of 2.008. The O adsorption characteristics can be different at the edges and corners of the crystal so we considered different surface sites on the IrO2 nanocluster such as surface Ir top sites on (110) and (101) surfaces, edges and corners formed at the intersection of (110), (100), and (101) facets which are depicted in Figure 12c. Table 1 indicates that corner sites at the IrO2 nanocrystal surface have up to 1.4 eV lower binding energy for O compared to bare (110) and (101) surfaces. In a previous DFT study85, it was demonstrated that the activity of IrO2 for water splitting, #

Adsorption site

Ead (eV/O atom)

dIr-O (Å)

1

(110) Ir top

-3.22

1.98

2

(101) Ir top

-3.06

1.95

3

(101)||(101) edge bridge

-4.08

2.01

on IrO2 surfaces. Here, we find that O

4

(101)|(101) corner

-3.37

1.92

5

(110)||(100) edge

-3.14

1.93

binding can be significantly stronger at the

6

(110)||(101) edge 1

-3.05

1.97

7

(110)||(101) edge 2

-3.44

1.92

8

(110)||(101)||(101) corner 1

-4.54

1.87

9

(110)||(101)||(101) corner 2

-2.71

1.95

and oxygen evolution reactions is directly related to the binding energy of atomic O

corners of a nanocrystal, and one can engineer the catalytic activity of IrO2 nanocrystals by selectively creating facets

Table 1. Oxygen adsorption energies on IrO2 nanocrystal. Reproduced from Ref (13).

to maximize the catalytic performance. In

31



addition, the amount of charge transfer at adsorption sites, in conjunction with coordination number, was found to predict the MS-Q adsorption energy. The use of variable charge is therefore critical for modeling catalytic processes. MS-Q was shown to be suitable for studying the O binding on IrO2 surfaces and it can be used to further investigate shape and size effects on binding properties of IrO2 nanoclusters that will affect its catalytic activity.

4.6. Reactive simulations of Ta/TaOx hetero-structure There are a wide variety of materials chemistry problems where bond formation and bond breakage phenomena are intricately associated with the dynamical evolution of interfaces. For example, oxidation and oxide growth of metals involves formation of O-metal bonds, which induces significant positive charges on the metal atoms and significant negative charges on the oxygen atoms. The atomic charges during the course of the oxidation process vary dynamically and are therefore environment dependent. In a metal-oxide heterostructure, the charges change continuously from a zero value in a fully metallic region to their valency-determined maximum value in the stoichiometric oxide. The charge transfer dynamics and oxygen stoichiometry have a strong bearing on their functional properties. For example, thin film heterostructures of Ta and Ta2O5 are potential platforms for neuromorphic computing86, where an applied electric field, has been exploited to switch between the two thermodynamically stable states i.e. metallic Ta and insulating Ta2O5. Fixed charge and even variable charge models such as MS-Q13 are not capable of capturing the multiple oxidation states encountered in oxide/metal hetero-interfaces. In Ref (86), use our ML framework to introduce a charge transfer ionic potential (CTIP) model for Ta/TaOx system that can describe the thermo-physical, structural and surface properties of Ta, and the various Ta2O5 polymorphs, as well as their interfaces. 4.6.1. Model selection: We choose the charge transfer ionic interatomic potential developed by Zhou and Wadley87 where the interaction potential is divided into two parts: an electrostatic contribution (Ees) for ionic interactions dominant in the oxide, and a non-electrostatic energy (Em) based on the embedded atom

32


Page 32 of 47

Page 33 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


method to model metallic interactions. The non-electrostatic interactions between the metal atoms are described using the embedded atom method (EAM) whereas the electrostatic interactions are calculated based on charges computed using electronegativity equalization method (EEM). 4.6.2. Training data: We employ three experimentally reported crystalline polymorphs of Ta2O5

88-89,

namely, monoclinic (C2/c), δ-hexagonal (P6/mmm) and β-orthorhombic (Pmmm). For each of these polymorphs, the training set included lattice constants, cohesive energy, internal coordinates as well as equation of state (energy vs. volume) from DFT+U calculations. Our training set additionally consists of 13 independent elastic constants of the monoclinic phase. Cross-validation is carried out by calculating surface energies of bcc Ta and monoclinic Ta2O5, the bulk modulus of all three polymorphs as well as the structure and dynamics of the polymorphs at 300 K and comparing with those from experiments and DFT calculations. 4.6.3. Parameter optimization: There are a total of 26 parameters, which require to be trained for the Ta/TaOx CTIP model. Note that the objective functional landscape for a 26-dimensional parameter set can be quite complex and consists of multiple minima. The optimized parameters obtained using local optimization alone will depend heavily on the initial guess as the objective settles into the nearest local minimum. GA optimization is used to perform a global search; we begin with a random population of Np = 50 parameter sets each with 26 parameters. We perform the appropriate genetic operations such as crossover and mutation and assess the fitness of a particular parameter set (individual) by the accuracy of its predicted properties against the training set. We choose Np individuals with the lowest objective function values for the next iteration of genetic operations and the process is continued till the objective converges to within tolerance limits. A multi-start simplex method is performed using the promising parameter sets (15 sets with the lowest error in prediction) to obtain the best possible parameter set. We compare the predicted properties against both the training and test sets to assess the accuracy of our ML trained CTIP model. The ML model captures the structure, energy and elastic constant components with reasonable accuracy for the monoclinic, hexagonal and orthorhombic polymorphs. The

33



energy ordering between the three polymorphs is also accurately computed and the equation of state (energy vs. volume shown in Figure 13a-c) for all the polymorphs agree well with DFT calculations. The surface energies of bcc Ta and monoclinic Ta2O5 for the (001), (010) and (100) surfaces are also in excellent agreement with both experiments and first-principles calculations, despite not being included in the training set. (d)

Simulation Time

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 13. Performance of our ML trained CTIP-EAM potential model with DFT calculations. The equations of state for (a) monoclinic, (b) hexagonal, and (c) orthorhombic Ta2O5 calculated using the EAM+Qeq parameters developed in Ref (86) (solid blue squares) and DFT calculations (solid red circles). Solid lines correspond to the Murnaghan fit. The energies are relative to the cohesive energy of the crystal at equilibrium (E0) as evaluated by the corresponding level of theory. The crystal volumes are normalized by the equilibrium value in the framework of the corresponding level of theory. (d) Snapshots showing the oxidation of Ta surface studied using this ML trained CTIP-EAM model. Reproduced with permission from Ref (86). 4.6.4. Model validation: We use the ML trained CTIP to investigate the differences in oxidation kinetics during the initial stages of tantalum oxide growth as a function of both temperature and the atomic/molecular nature of the oxidizing species (Figure 13b). We study the formation mechanism and microstructure of these oxide films, the early stage oxidation growth kinetics as well as atomistic details of the composition, microstructure and the limiting thickness of the Ta oxide films. Our simulation results are in good agreement with experimental results on low-temperature Ta oxidation. These ML trained CTIP

34


Page 34 of 47

Page 35 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


models can be used to probe the atomistic mechanisms responsible for the sub-nanosecond timescale switching behavior of TaOx memristors, and other transport phenomena at Ta/TaOx interfaces under external stimuli. It should be noted that while all the examples shown above were focused on inorganic systems, the workflow and the methodology can be extended to soft-matter systems as well. In our recent work90, we used a similar workflow to train the parameters of a polarizable force field ab initio data from quantum mechanics (QM) calculations of molecular clusters at the MP2/6-31G(d,p), DFMP2(fc)/jul-cc-pVDZ, and DFMP2(fc)/jul-cc-pVTZ levels. The target properties included experimental condensed phase properties (i.e., density and heat of vaporization) and the optimized model was found to better than the original AMOEBA force field (amoeba09.prm), which was optimized empirically to match liquid properties (see Ref (90) for details). We thus envision wide usage of automated ML workflows for both inorganic and softmatter systems. 5. Future Directions and Perspective With continued increase in computing power, MD is emerging as a powerful tool for materials modeling as well as for advanced nanoscale materials characterization. The predictive power of MD hinges strongly on the interatomic potential used to describe the atomistic interactions in the system. While the ML framework and the examples presented above highlight the advantages of using data driven approaches for training new models, there is still a lot of room for improvement. Although the list below is not exhaustive, we believe that these are certainly areas that would require further focused studies in the near future from researchers worldwide. 

Most of the classical MD simulations employ pre-defined functional forms that can often limit the chemistry and physics that can be captured. While it appears that there can be significant improvements made by using data driven approaches that employ extensive training data sets and advanced optimization, there will always be a ceiling limit imposed by the use of pre-defined functional forms. Existing force-fields with pre-defined functional form are thus not quite flexible and cannot be

35



transferred easily from one material class to another (e.g. metals to oxides). There is clearly a need to incorporate flexibility in the functional form. In this regard, neural network models91 and potential models based on basis function expansions provide more flexibility. 

Regardless of the application domain or area, all ML frameworks require carefully prepared and sampled training data. In using ML for developing material models, the training data should be sufficiently diverse i.e. consist of structures, which span a wide range of energies. As highlighted in the example of Au clusters, there is a clear need to go beyond sampling near equilibrium configurations and generating a training data that includes far from equilibrium configurations as well. To develop a predictive model, the developer has to ensure ample representation from different parts of the potential energy surface in the training data for the model to be robust.



Classical MD is limited in the timescales it can access. For example, all-atom simulations that employ 0.1-1 femtosecond timesteps can only access nanosecond timescales. It is worth noting that even with the move towards exascale computing, the gains are going to be primarily in terms of spatial scales that can be accessed but given that the clockspeeds and bandwidths are not going to change significantly, timescale challenges will remain92. Accelerated MD techniques93 represent one way to address this time-scale challenge, but another option in some cases is to employ coarse-grained (CG) models94 that typically allow 10-50 femtosecond timesteps. The challenge with such models is that the gain in efficiency comes at the cost of accuracy. Many of the techniques highlighted in this work can be extended to train CG models which could potentially allow us to overcome some timescale challenges without sacrificing accuracy. In particular, there is a significant gap that exists between the atomistic models and accurate CG models that can access the mesoscopic regime.



Classical MD does well mostly in the case of static properties but often lacks predictive power when it comes to dynamical and transport properties. One way to address this challenge is to include forces in the training dataset. Going forward, we envision that temperature dependent properties derived from on-the-fly MD can also be included as part of the training procedure. This would allow us to directly

36


Page 36 of 47

Page 37 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


train MD potentials that can also capture transport properties and other dynamical or temperature dependent properties of interest. 

Cross validation and iterative improvement techniques are seldom used in the fitting of potentials. Even with a first principles based dataset (e.g. derived from DFT), there will always be errors that will be propagated to the atomistic potential model. While one can envision using higher-level theory95-98 (CCSD or QMC) to generate the training data and reduce the error, there is a clear need for quantification of uncertainties in the predictions at various scales99. Recent work by Strachan and coworkers on Functional Uncertainty Quantification method (FunUQ) to assess model form errors, such as interatomic potentials for MD simulations represents an important future direction100-101. Crossvalidation, sensitivity analysis and uncertainty quantification will be critical to improve the quality and predictive power of the interatomic potentials for MD.



Our examples clearly highlight the advantage of ML global optimization algorithms and multi-stage optimization strategies over local optimization procedures that were traditionally employed. Evolutionary optimization procedures such as GA remain attractive but one can easily envision the use of other techniques such as differential evolution to derive globally optimal solutions. Use of more recent derivative free optimization techniques such as POUNDerS56 either standalone or in conjunction with other global optimization methods may have advantages in training interatomic potentials for MD.



Finally, one needs to have objective definitions that go beyond the simple sum of squared difference. Single objective functions are often plagued by the selection of the appropriate weights. In this respect, multi-objective optimization73 such as Pareto optimization can provide solutions that represent a compromise between the various desired objectives. One needs to exercise caution in defining the objective appropriately otherwise even global minimum will not yield the best set. In summary, machine learning has the potential to seamlessly bridge the various disparate scales in materials modeling, thereby dramatically accelerating the materials design and discovery process. This can lead to transformative advances by allowing us to access time and length scales in materials modeling that were hitherto considered to be inaccessible to molecular dynamics.

37



Acknowledgements: The authors thank Alper Kinaci, and Troy Loeffler for their significant contributions to the development of the ML workflow for training first-principles based force-fields. Use of the Center for Nanoscale Materials, an Office of Science user facility, was supported by the U. S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. Use of the Argonne Leadership Computing Facility is supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-06CH11357. Use of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility, is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The authors also acknowledge the computing resources provided on Fusion and Blues, high performance computing clusters operated by the Laboratory Computing Resource Center (LCRC) at Argonne National Laboratory.

38


Page 38 of 47

Page 39 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References: 1. Berman, D.; Deshmukh, S. A.; Sankaranarayanan, S. K. R. S.; Erdemir, A.; Sumant, A. V. Macroscale Superlubricity Enabled by Graphene Nanoscroll Formation. Science 2015, 348, 1118-1122. 2. Jiang, Z.; He, J.; Deshmukh, S. A.; Kanjanaboos, P.; Kamath, G.; Wang, Y.; Sankaranarayanan, S. K. R. S.; Wang, J.; Jaeger, H. M.; Lin, X.-M. Subnanometer Ligand Shell Asymmetry Leads to JanusLike Nanoparticle Membranes. Nat. Mater. 2015, 14, 912-917. 3. Erdemir, A.; Ramirez, G.; Eryilmaz, O. L.; Narayanan, B.; Liao, Y.; Kamath, G.; Sankaranarayanan, S. K. R. S. Carbon-Based Tribofilms from Lubricating Oils. Nature 2016, 536, 67-71. 4. Zhang, Z.; Schwanz, D.; Narayanan, B.; Kotiuga, M.; Dura, J. A.; Cherukara, M.; Zhou, H.; Freeland, J. W.; Li, J.; Sutarto, R. et al. Perovskite Nickelates as Electric-Field Sensors in Salt Water. Nature 2017, 553, 68. 5. Zuo, F.; Panda, P.; Kotiuga, M.; Li, J.; Kang, M.; Mazzoli, C.; Zhou, H.; Barbour, A.; Wilkins, S.; Narayanan, B. et al. Habituation Based Synaptic Plasticity and Organismic Learning in a Quantum Perovskite. Nat. Commun. 2017, 8, 240. 6. Cherukara, M. J.; Narayanan, B.; Chan, H.; Sankaranarayanan, S. K. R. S. Silicene Growth through Island Migration and Coalescence. Nanoscale 2017, 9, 10186-10192. 7. Narayanan, B.; Chan, H.; Kinaci, A.; Sen, F. G.; Gray, S. K.; Chan, M. K. Y.; Sankaranarayanan, S. K. R. S. Machine Learnt Bond Order Potential to Model Metal–Organic (Co–C) Heterostructures. Nanoscale 2017, 9, 18229-18239. 8. Sun, Y.; Kotiuga, M.; Lim, D.; Narayanan, B.; Cherukara, M.; Zhang, Z.; Dong, Y.; Kou, R.; Sun, C.-J.; Lu, Q., et al. Strongly Correlated Perovskite Lithium Ion Shuttles. Proc. Natl. Acad. Sci. 2018. 9. Asadi, M.; Sayahpour, B.; Abbasi, P.; Ngo, A. T.; Karis, K.; Jokisaari, J. R.; Liu, C.; Narayanan, B.; Gerard, M.; Yasaei, P., et al. A Lithium–Oxygen Battery with a Long Cycle Life in an Air-Like Atmosphere. Nature 2018, 555, 502. 10. Pang, Q.; Shyamsunder, A.; Narayanan, B.; Kwok, C. Y.; Curtiss, L. A.; Nazar, L. F. Tuning the Electrolyte Network Structure to Invoke Quasi-Solid State Sulfur Conversion and Suppress Lithium Dendrite Formation in Li–S Batteries. Nat. Energy 2018, 3, 783-791. 11. Ulvestad, A.; Sasikumar, K.; Kim, J. W.; Harder, R.; Maxey, E.; Clark, J. N.; Narayanan, B.; Deshmukh, S. A.; Ferrier, N.; Mulvaney, P. et al. In Situ 3d Imaging of Catalysis Induced Strain in Gold Nanoparticles. J. Phys. Chem. Lett. 2016, 7, 3008-3013. 12. Bagchi, B.; Chakravarty, C. Interplay Between Multiple Length and Time Scales in Complex Chemical Systems. J. Chem. Sci. 2010, 122, 459-470. 13. Sen, F. G.; Kinaci, A.; Narayanan, B.; Gray, S. K.; Davis, M. J.; Sankaranarayanan, S. K.; Chan, M. Towards Accurate Prediction of Catalytic Activity in Iro 2 Nanoclusters Via First Principles-Based Variable Charge Force Field. J. Mater. Chem. A 2015, 3, 18970-18982. 14. Narayanan, B.; Kinaci, A.; Sen, F. G.; Davis, M. J.; Gray, S. K.; Chan, M. K. Y.; Sankaranarayanan, S. K. R. S. Describing the Diverse Geometries of Gold from Nanoclusters to Bulk− a First-Principles Based Hybrid Bond Order Potential. J. Phys. Chem. C 2016, 120, 13787-13800. 15. Narayanan, B.; Sasikumar, K.; Mei, Z.-G.; Kinaci, A.; Sen, F. G.; Davis, M. J.; Gray, S. K.; Chan, M. K.; Sankaranarayanan, S. K. Development of a Modified Embedded Atom Force Field for Zirconium Nitride Using Multi-Objective Evolutionary Optimization. J. Phys. Chem. C 2016, 120, 17475-17483. 16. Cherukara, M. J.; Narayanan, B.; Kinaci, A.; Sasikumar, K.; Gray, S. K.; Chan, M. K.; Sankaranarayanan, S. K. Ab Initio-Based Bond Order Potential to Investigate Low Thermal Conductivity of Stanene Nanostructures. J. Phys. Chem. Lett. 2016, 7, 3752-3759. 17. Kinaci, A.; Narayanan, B.; Sen, F. G.; Davis, M. J.; Gray, S. K.; Sankaranarayanan, S. K. R. S.; Chan, M. K. Y. Unraveling the Planar-Globular Transition in Gold Nanoclusters Through Evolutionary Search. Sci. Rep. 2016, 6, 34974. 18. Stillinger, F.; Weber, T. A. Computer Simulation of Local Order in Condensed Phases of Silicon. Phys. Rev. B 1985, 31, 5262-5272.

39



19. Vashishta, P.; Kalia, R. K.; Rino, J. P.; Ebbsjö, I. Interaction Potential for Sio2: A MolecularDynamics Study of Structural Correlations. Phys. Rev. B 1990, 41, 12197-12209. 20. Rafii-Tabar, H.; Sutton, A. P. Long-Range Finnis-Sinclair Potentials for Fcc Metallic Alloys. Philos. Mag. Lett. 1991, 63, 217-224. 21. Daw, M. S.; Baskes, M. I. Embedded-Atom Method: Derivation and Application to Impurities, Surfaces and Other Defects in Metals. Phys. Rev. B 1983, 29, 6443-6454. 22. Baskes, M. I. Modified Embedded-Atom Potentials for Cubic Materials and Impurities. Phys. Rev. B 1992, 46, 2727-2743. 23. Tersoff, J. New Empirical Approach for the Structure and Energy of Covalent Systems. Phys. Rev. B 1988, 37, 6991-7001. 24. Tersoff, J. Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon. Phys. Rev. Lett. 1988, 61, 2879. 25. Abell, G. C. Empirical Chemical Pseudopotential Theory of Molecular and Metallic Bonding. Phys. Rev. B 1985, 31, 6184. 26. Brenner, D. W. Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor Deposition of Diamond Films. Phys. Rev. B 1990, 42, 9458-9472. 27. Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. A SecondGeneration Reactive Empirical Bond Order (Rebo) Potential Energy Expression for Hydrocarbons. J. Phys.: Condens. Matter 2002, 14, 783–802. 28. Jelinek, B.; Groh, S.; Horstemeyer, M. F.; Houze, J.; Kim, S. G.; Wagner, G. J.; Moitra, A.; Baskes, M. I. Modified Embedded Atom Method Potential for Al, Si, Mg, Cu, and Fe Alloys. Phys. Rev. B 2012, 85, 245102. 29. Cornell, W.; Cieplak, P.; Bayly, C.; Gould, I.; Merz, K. J.; Ferguson, D.; Spellmeyer, D.; Fox, T.; Caldwell, J.; Kollman, P. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. 30. Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the Opls AllAtom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225–11236. 31. Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Karplus, S. S. M. Charmm: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187-217. 32. Sun, H.; Mumby, S. J.; Maple, J. R.; Hagler, A. T. An Ab Initio Cff93 All-Atom Force Field for Polycarbonates. J. Am. Chem. Soc. 1994, 116, 2978-2987. 33. Chenoweth, K.; Duin, A. C. T. v.; Goddard, W. A. Reaxff Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040-1053. 34. Duin, A. C. T. v.; Dasgupta, S.; Lorant, F.; Goddard, W. A. Reaxff: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2003, 105, 9396-9409. 35. Shan, T.-R.; Devine, B. D.; Hawkins, J. M.; Asthagiri, A.; Phillpot, S. R.; Sinnott, S. B. SecondGeneration Charge-Optimized Many-Body Potential for Si/Sio 2 and Amorphous Silica. Phys. Rev. B 2010, 82, 235302. 36. Devine, B.; Shan, T.-R.; Cheng, Y.-T.; McGaughey, A. J.; Lee, M.; Phillpot, S. R.; Sinnott, S. B. Atomistic Simulations of Copper Oxidation and Cu/Cu 2 O Interfaces Using Charge-Optimized ManyBody Potentials. Phys. Rev. B 2011, 84, 125308. 37. Rappe, A. K.; Goddard III, W. A. Charge Equilibration for Molecular Dynamics Simulations. J. Phys. Chem. 1991, 95, 3358-3363. 38. Monti, S.; Duin, A. C. T. v.; Kim, S. Y.; Barone, V. Exploration of the Conformational and Reactive Dynamics of Glycine and Diglycine on Tio2: Computational Investigations in the Gas Phase and in Solution. J. Phys. Chem. C 2012, 116, 5141-5150. 39. Fenter, P.; Kerisit, S.; Raiteri, P.; Gale., J. D. Is the Calcite–Water Interface Understood? Direct Comparisons of Molecular Dynamics Simulations with Specular X-Ray Reflectivity Data? J. Phys. Chem. C 2013, 117, 5028-5042.

40


Page 40 of 47

Page 41 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


40. Deshmukh, S. A.; Kamath, G.; Suthar, K. J.; Mancini, D. C.; Sankaranarayanan, S. K. NonEquilibrium Effects Evidenced by Vibrational Spectra During the Coil-to-Globule Transition in Poly (NIsopropylacrylamide) Subjected to an Ultrafast Heating–Cooling Cycle. Soft Matter 2014, 10, 1462-1480. 41. Kamath, G.; Deshmukh, S. A.; Baker, G. A.; Mancini, D. C.; Sankaranarayanan, S. K. Thermodynamic Considerations for Solubility and Conformational Transitions of Poly-N-IsopropylAcrylamide. Phys. Chem. Chem. Phys. 2013, 15, 12667-12673. 42. Deshmukh, S. A.; Sankaranarayanan, S. K.; Suthar, K.; Mancini, D. C. Role of Solvation Dynamics and Local Ordering of Water in Inducing Conformational Transitions in Poly (NIsopropylacrylamide) Oligomers Through the Lcst. J. Phys. Chem. B 2012, 116, 2651-2663. 43. ErkoÇ, A. Empirical Potential Energy Functions Used in the Simulations of Materials Properties. In Annual Reviews of Computational Physics Ix, WORLD SCIENTIFIC: 2001; Vol. Volume 9, pp 1-103. 44. Tersoff, J. Modelings Solid-State Chemistry- Interatomic Potentials for Multicomponent Systems. Phys. Rev. B 1989, 39, 5566-5568. 45. Lee, B. J.; Shim, J. H.; Baskes, M. I. Semiempirical Atomic Potentials for the Fcc Metals Cu, Ag, Au, Ni, Pd, Pt, Al, and Pb Based on First and Second Nearest-Neighbor Modified Embedded Atom Method. Phys. Rev. B 2003, 68. 46. McCall, J. Genetic Algorithms for Modelling and Optimisation. J. Comput. Appl. Math. 2005, 184, 205-222. 47. Y, G.; I, T.; DL, T.; AF, W.; M, M. Interpolating Moving Least-Squares Methods for Fitting Potential Energy Surfaces: Applications to Classical Dynamics Calculations. J. Chem. Phys. 2007, 127, 214106. 48. Maragakis, P.; van der Vaart, A.; Karplus, M. Gaussian-Mixture Umbrella Sampling. J. Phys. Chem. B 2009, 113, 4664-4673. 49. Ceriotti, M.; Tribello, G. A.; Parrinello, M. Simplifying the Representation of Complex FreeEnergy Landscapes Using Sketch-Map. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 13023-13028. 50. Coifman, R. R.; Lafon, S.; Lee, A. B.; Maggioni, M.; Nadler, B.; Warner, F.; Zucker, S. W. Geometric Diffusions as a Tool for Harmonic Analysis and Structure Definition of Data: Multiscale Methods. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 7432-7437. 51. Li, Z.; Kermode, J. R.; Vita, A. D. Molecular Dynamics with on-the-Fly Machine Learning of Quantum-Mechanical Forces. Phys. Rev. Lett. 2015, 114, 096405. 52. Bartók, A. P.; Payne, M. C.; Kondor, R.; Csányi, G. Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, Without the Electrons. Phys. Rev. Lett. 2010, 104, 136403. 53. Diwekar, U. Introduction to Applied Optimization; Springer Science & Business Media, LLC, 2008; Vol. 22. 54. Dantzig, G. B. Linear Programming and Extensions; Princeton University Press: Princeton, NJ, 1963. 55. Marquardt, D. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. SIAM J. Appl. Math. 1963, 11, 431-441. 56. Wild, S. M.; Sarich, J.; Schunck, N. Derivative-Free Optimization for Parameter Estimation in Computational Nuclear Physics. J. Phys. G: Nucl. Part. Phys. 2015, 42. 57. Wild, S. M.; Regis, R. G.; Shoemaker, C. A. Orbit: Optimization by Radial Basis Function Interpolation in Trust-Regions. SIAM J. Sci. Comput. 2008, 30, 3197-3219. 58. Sen, F. G.; Narayanan, B.; Larson, J.; Kinaci, A.; Sasikumar, K.; Davis, M. J.; Wild, S. M.; Gray, S. K.; Sankaranarayanan, S. K. R. S.; Chan, M. K. Y. Comparing Optimization Strategies for Force Field Parameterization. arXiv preprint arXiv:1812.00326 2018. 59. Fonseca, C. M.; Fleming, P. J. Genetic Algorithms for Multiobjective Optimization: Formulation, Discussion and Generalization. Genetic Algorithms: Proceedings of the Fifth International Conference (S. Forrest, ed.) 1993, San Mateo, CA: Morgan Kaufmann. 60. Wilde, M.; Hategan, M.; Wozniak, J. M.; Clifford, B.; Katz, D. S.; Foster, I. Swift: A Language for Distributed Parallel Scripting. Parallel Comput. 2011, 37, 633-652. 61. Pyykkö, P. Theoretical Chemistry of Gold. Angew. Chem. 2004, 43, 4412.

41



62. Schmidbaur, H.; Schier, A. A Briefing on Aurophilicity. Chem. Soc. Rev. 2008, 37, 1931-1951. 63. Li, X.; Tollberg, B.; Xiankui, H.; Lichang, W. Structural Study of Gold Clusters. J. Chem. Phys. 2006, 124, 114309-1-114309-10. 64. Serapian, S. A.; Bearpark, M. J.; Bresme, F. The Shape of Au-8: Gold Leaf or Gold Nugget? Nanoscale 2013, 5, 6445-6457. 65. Gruene, P.; Rayner, D. M.; Redlich, B.; van der Meer, A. F. G.; Lyon, J. T.; Meijer, G.; Fielicke, A. Structures of Neutral Au-7, Au-19, and Au-20 Clusters in the Gas Phase. Science 2008, 321, 674-676. 66. Bulusu, S.; Li, X.; Wang, L.-S.; Zeng, X. C. Evidence of Hollow Golden Cages. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 8326-8330. 67. Hakkinen, H.; Yoon, B.; Landman, U.; Li, X.; Zhai, H. J.; Wang, L. S. On the Electronic and Atomic Structures of Small Au-N(-) (N=4-14) Clusters: A Photoelectron Spectroscopy and DensityFunctional Study. J. Phys. Chem. A 2003, 107, 6168-6175. 68. Sutton, A. P.; Chen, J. Long-Range Finnis-Sinclair Potentials. Philos. Mag. Lett. 1990, 61, 139146. 69. Gupta, R. P. Lattice-Relaxation at a Metal-Surface. Phys. Rev. B 1981, 23, 6265-6270. 70. Masuda, K. Lattice-Relaxation at a Metal-Surface: Comments. Phys. Rev. B 1982, 26, 5968-5971. 71. Foiles, S. M.; Baskes, M. I.; Daw, M. S. Embedded-Atom-Method Functions for the Fcc Metals Cu, Ag, Au, Ni, Pd, Pt, and Their Alloys. Phys. Rev. B 1986, 33, 7983-7991. 72. Keith, J. A.; Fantauzzi, D.; Jacob, T.; van Duin, A. C. T. Reactive Forcefield for Simulating Gold Surfaces and Nanoparticles. Phys. Rev. B 2010, 81. 73. Sastry, K. Single and Multiobjective Genetic Algorithm Toolbox for Matlab in C++. 2007. 74. Zhu, F.-f.; Chen, W.-j.; Xu, Y.; Gao, C.-l.; Guan, D.-d.; Liu, C.-h.; Qian, D.; Zhang, S.-C.; Jia, J.f. Epitaxial Growth of Two-Dimensional Stanene. Nat. Mater. 2015, 14, 1020-1025. 75. Tachibana, Y.; Vayssieres, L.; Durrant, J. R. Artificial Photosynthesis for Solar Water-Splitting. Nat. Photonics 2012, 6, 511-518. 76. Blakemore, J. D.; Mara, M. W.; Kushner-Lenhoff, M. N.; Schley, N. D.; Konezny, S. J.; Rivalta, I.; Negre, C. F. A.; Snoeberger, R. C.; Kokhan, O.; Huang, J. et al. Characterization of an Amorphous Iridium Water-Oxidation Catalyst Electrodeposited from Organometallic Precursors. Inorg. Chem. 2013, 52, 1860-1871. 77. Huang, J.; Blakemore, J. D.; Fazi, D.; Kokhan, O.; Schley, N. D.; Crabtree, R. H.; Brudvig, G. W.; Tiede, D. M. Domain Structure for an Amorphous Iridium-Oxide Water-Oxidation Catalyst Characterized by X-Ray Pair Distribution Function Analysis. Phys. Chem. Chem. Phys. 2014, 16, 18141819. 78. Bolzan, A. A.; Fong, C.; Kennedy, B. J.; Howard, C. J. Structural Studies of Rutile-Type Metal Dioxides. Acta Crystallogr. B 1997, 53, 373-380. 79. Ono, S.; Kikegawa, T.; Ohishi, Y. High-Pressure and High-Temperature Synthesis of a Cubic Iro2 Polymorph. Physica B Condens. Matter 2005, 363, 140-145. 80. Panda, S. K.; Bhowal, S.; Delin, A.; Eriksson, O.; Dasgupta, I. Effect of Spin Orbit Coupling and Hubbard U on the Electronic Structure of Iro2. Phys. Rev. B 2014, 89, 155102. 81. Grindy, S.; Meredig, B.; Kirklin, S.; Saal, J. E.; Wolverton, C. Approaching Chemical Accuracy with Density Functional Calculations: Diatomic Energy Corrections. Phys. Rev. B 2013, 87, 075150. 82. Novell-Leruth, G.; Carchini, G.; Lopez, N. On the Properties of Binary Rutile Mo2 Compounds, M = Ir, Ru, Sn, and Ti: A Dft Study. J. Chem. Phys. 2013, 138, 194706. 83. Perron, H.; Domain, C.; Roques, J.; Drot, R.; Simoni, E.; Catalette, H. Optimization of Accurate Rutile Tio2 (110), (100), (101) and (001) Surface Models from Periodic Dft Calculations. Theor. Chem. Acc. 2007, 117, 565-574. 84. Zucker, R. V.; Chatain, D.; Dahmen, U.; Hagege, S.; Carter, W. C. New Software Tools for the Calculation and Display of Isolated and Attached Interfacial-Energy Minimizing Particle Shapes. J. Mater. Sci. 2012, 47, 8290--8302. 85. Rossmeisl, J.; Qu, Z. W.; Zhu, H.; Kroes, G. J.; Norskov, J. K. Electrolysis of Water on Oxide Surfaces. J. Electroanal. Chem. 2007, 607, 83-89.

42


Page 42 of 47

Page 43 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


86. Sasikumar, K.; Narayanan, B.; Cherukara, M.; Kinaci, A.; Sen, F. G.; Gray, S. K.; Chan, M. K. Y.; Sankaranarayanan, S. K. R. S. Evolutionary Optimization of a Charge Transfer Ionic Potential Model for Ta/Ta-Oxide Heterointerfaces. Chem. Mater. 2017, 29, 3603-3614. 87. Zhou, X. W.; Wadley, H. N. G. A Charge Transfer Ionic–Embedded Atom Method Potential for the O–Al–Ni–Co–Fe System. J. Phys.: Condens. Matter 2005, 17, 3619. 88. Ivanov, M. V.; Perevalov, T. V.; Aliev, V. S.; Gritsenko, V. A.; Kaichev, V. V. Electronic Structure of Δ-Ta2o5 with Oxygen Vacancy: Ab Initio Calculations and Comparison with Experiment. J. Appl. Phys. 2011, 110, 024115. 89. Nashed, R.; Hassan, W. M. I.; Ismail, Y.; Allam, N. K. Unravelling the Interplay of Crystal Structure and Electronic Band Structure of Tantalum Oxide (Ta2o5). Phys. Chem. Chem. Phys. 2013, 15, 1352-1357. 90. Li, Y.; Li, H.; Pickard, F. C.; Narayanan, B.; Sen, F. G.; Chan, M. K. Y.; Sankaranarayanan, S. K. R. S.; Brooks, B. R.; Roux, B. Machine Learning Force Field Parameters from Ab Initio Data. J. Chem. Theory Comput. 2017, 13, 4492-4503. 91. Artrith, N.; Urban, A. An Implementation of Artificial Neural-Network Potentials for Atomistic Materials Simulations: Performance for Tio2. Comput. Mater. Sci. 2016, 114, 135-150. 92. Germann, T.; Richards, D.; McPherson, A.; Belak, J. Exascale Co-Design Center for Materials in Extreme Environments (Exmatex) Annual Report - Year 2. 2013. 93. Hamelberg, D. Accelerated Molecular Dynamics: A Promising and Efficient Simulation Method for Biomolecules. J. Chem. Phys. 2004, 120, 11919. 94. Molinero, V.; Moore, E. B. Water Modeled as an Intermediate Element Between Carbon and Silicon. J. Phys. Chem. B 2009, 113, 4008-4016. 95. Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-Set Convergence of Correlated Calculations on Water. J. Chem. Phys. 1997, 106, 9639-9646. 96. Sedlak, R.; Riley, K. E.; Řezáč, J.; Pitoňák, M.; Hobza, P. Mp2.5 and Mp2.X: Approaching Ccsd(T) Quality Description of Noncovalent Interaction at the Cost of a Single Ccsd Iteration. ChemPhysChem 2013, 14, 698-707. 97. Bozkaya, U.; Sherrill, C. D. Orbital-Optimized Mp2.5 and Its Analytic Gradients: Approaching Ccsd(T) Quality for Noncovalent Interactions. J. Chem. Phys. 2014, 141, 204105. 98. Feller, D. Application of Systematic Sequences of Wave Functions to the Water Dimer. J. Chem. Phys. 1992, 96, 6104-6114. 99. Fang, H.; Demir, H.; Kamakoti, P.; Sholl, D. S. Recent Developments in First-Principles Force Fields for Molecules in Nanoporous Materials. J. Mater. Chem. A 2014, 2, 274-291. 100. Reeve, S. T.; Strachan, A. Error Correction in Multi-Fidelity Molecular Dynamics Simulations Using Functional Uncertainty Quantification. J. Comput. Phys. 2017, 334, 207-220. 101. Hunt, M.; Haley, B.; McLennan, M.; Koslowski, M.; Murthy, J.; Strachan, A. Puq: A Code for Non-Intrusive Uncertainty Propagation in Computer Simulations. Comput. Phys. Commun. 2015, 194, 97107.

43



TOC Graphic

44


Page 44 of 47

Page 45 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Biography: Henry Chan received his Ph.D. in Computational Chemistry from the University of Illinois at Chicago in 2015. He is currently a postdoctoral researcher at the Center for Nanoscale Materials, Argonne National Laboratory. He has extensive experience in simulations of inorganic and soft-matter material systems. His current research focus is on the application of machine learning principles in the development of interatomic force fields and techniques that enable accurate large-scale molecular simulations.

Badri Narayanan received his Ph.D. in Materials Science from Colorado School of Mines in 2013. He is currently an Assistant Professor of Mechanical Engineering at University of Louisville, Kentucky, USA. Prior to this, he was a postdoctoral researcher at the Center for Nanoscale Materials, Argonne National Laboratory from 2013-2016 and subsequently a staff scientist in Materials Science Division at Argonne. His research focuses on using machine learning methods to develop materials models at atomic and mesoscopic scales; and employing first-principles, atomistic, and meso-scale methods to gain molecularlevel understanding of structure-property-processing relationships in materials for a variety of energy technologies, including batteries, tribology, sensing, catalysis, and advanced manufacturing.

Mathew Cherukara received his Ph.D in computational materials engineering from Purdue University in 2015 and his bachelors in materials engineering from the Indian Institute of Technology (IIT) Madras. He is now an assistant scientist at the Center for Nanoscale Materials at Argonne National Laboratory. Dr. Cherukara's work revolves around the use of advanced optimization and machine learning techniques, in particular deep learning and generative AI for model development and data analysis of materials simulations and 4D X-ray imaging.

45



Fatih G. Sen received his Ph.D. in Engineering Materials from University of Windsor in 2013. Following that he joined the Center for Nanoscale Materials in Argonne National Laboratory as a postdoctoral researcher. His research interests are computational materials design, energy storage and conversion, nanotechnology, surface science, and artifical intelligence. He is currently a scientist at Novelis Global Research and Technology Center.

Kiran Sasikumar received his Ph.D. in Materials Science and Engineering from Rensselaer Polytechnic Institute in 2014. He is currently a research scientist and scientific software developer at Avant-garde Materials Simulation, Deutschland. Formerly, as a postdoctoral researcher at the Center for Nanoscale Materials in Argonne National Laboratory, he worked on investigating materials interaction at nanoscale interfaces and surfaces using a wide variety of tools, including coherent diffractive imaging, molecular dynamics, and continuum finite element simulations. As a computational materials scientist, his research interests span crystal structure prediction, thermodynamics, nucleation and phase change, heat transfer, and reaction kinetics in nanomaterials.

Stephen K. Gray obtained a BSc. degree in Chemistry from Carleton University, Ottawa, and a Ph.D. degree in Chemistry from the University of California at Berkeley. After post-doctoral studies at Oxford University and the University of Chicago he was an Assistant Professor of Chemistry at Northern Illinois University before joining the staff at Argonne National Laboratory. He is currently Senior Scientist and Group Leader of the Theory and Modeling Group of Argonne’s Center for Nanoscale Materials. His research interests include light-interactions with nanoscale materials and quantum information science.

46


Page 46 of 47

Page 47 of 47 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Maria Chan received her PhD in Physics from the Massachusetts Institute of Technology. Since 2012, Chan has been a staff scientist at the Center of Nanoscale Materials, part of Argonne National Laboratory. She is also a member of the University of Chicago Consortium for Advanced Science and Engineering, and a senior fellow at the Northwestern-Argonne Institute of Science and Engineering. Chan's research focuses on the computational prediction of materials properties, using first principles, atomistic, and machine learning methods, particularly in applications towards materials relevant to energy technologies, such as energy storage, photovoltaics, catalysis, and thermal management.

Subramanian Sankaranarayanan is a Scientist in the Nanoscale Science and Technology (NST) Division at Argonne and a Senior Fellow at the Institute of Molecular Engineering at University of Chicago. He obtained his Ph.D. in Chemical Engineering from University of South Florida in 2007 and was a postdoctoral fellow at the School of Engineering and Applied Sciences at Harvard University from 2007-2010. His research at Argonne focuses on the use of machine learning to develop predictive models by bridging the electronic, atomistic and mesoscopic scales. Other interests include inverse design of materials and machine learning for integrating atomistic simulations with X-ray imaging. His interests span a diverse range of applications from tribology and corrosion to neuromorphic computing and thermal management.

47


Machine Learning Classical Interatomic Potentials for Molecular

Recommend Documents