Combining Machine Learning and Computational Chemistry for Predictive Insights Into Chemical Systems

Introduction

A lasting challenge in applied physical and chemical sciences has been to answer the question: how can one identify and make chemical compounds or materials that have optimal properties for a given purpose? A substantial part of research in physics, chemistry, and materials science concerns the discovery and characterization of novel compounds that can benefit society, but most advances still are generally attributed to trial-and-error experimentation, and this requires significant time and cost. Current global challenges create greater urgency for faster, better, and less expensive research and development efforts. Computational chemistry (CompChem) methods have significantly improved over time, and they promise paradigm shifts in how compounds are fundamentally understood and designed for specific applications.

Machine learning (ML) methods have in the last decades witnessed an unprecedented technological evolution enabling a plethora of applications, some of which have become daily companions in our lives. 1, 2, 3. Applications of ML include technological fields such as web search, translation, natural language processing, self-driving vehicles, control architectures, and in the sciences, e.g. medical diagnostics 4, 5, 6, 7, particle physics 8, nano sciences,9 bioinformatics, 10, 11 brain-computer interfaces, 12 social media analysis, 13 robotics, 14, 15 and team, social, or board games. 16, 17, 18 These methods have also become popular for accelerating the discovery and design of new materials, chemicals, and chemical processes. 19 At the same time, we have witnessed hype, criticism, and misunderstanding about how ML tools are to be used in chemical research. From this, we see a need for researchers working at the intersection of CompChem+ML to more critically recognize the true strengths and weaknesses of each component in any given study. Specifically, we wanted to review why and how CompChem+ML can provide useful insights into the study of molecules and materials.

While developing this review, we polled the scientific community with an anonymous online survey that asked for questions and concerns regarding the use of ML models with chemistry applications. Respondents raised excellent points including:

ML methods are becoming less understood while they are also more regularly used as black box tools.

Many publications show inadequate technical expertise in ML (e.g. inappropriate splitting of training, testing, and validation sets).

It can be difficult to compare different ML methods and know which is the best for a particular application, or whether ML should even be used at all.

Data quality and context are often missing from ML modeling, and data sets need to be made freely available and clearly explained.

Additionally, when asked about the most exciting active and emerging areas of ML in the next five years, respondents mentioned a wide range of topics from catalysis discovery, drug and peptide design, “above the arrow” reaction predictions, and generative models that promise to fundamentally transform chemical discovery. When asked about challenges that ML will not surmount in the next five years, respondents mentioned modeling complex photochemical and electrochemical environments, discovering exact exchange-correlation functionals, and completely autonomous reaction discovery. This review will give our perspective on many of these topics.

As context for this review, Figure 1 shows a heatmap depicting the frequency of ML keywords found in scientific articles that also have keywords associated with different American Chemical Society (ACS) technical divisions. Preparing this figure required several steps. First, lists of ML keywords were chosen. Second, lists of keywords were created by perusing ACS division symposia titles from over the past five years. Third, Python scripts used Scopus Application Programming Interfaces (APIs) to identify the number of scientific publications that matched sets of ML and division symposia keywords. Figure 1 elucidates several interesting points. First, the most popular ML approaches across all divisions are clearly neural networks, followed by genetic algorithms and support vector machines/kernel methods. Second, divisions such as physical (PHYS), analytical (ANYL), and environmental (ENVR) are already using diverse sets of ML approaches, while divisions such as inorganic (INOR), nuclear (NUCL), and carbohydrate (CARB) are primarily employing more distinct subsets of approaches, while other divisions such as educational (CHED), history (HIST), law (CHAL), and business-oriented divisions (BMGT and SCHB), i.e. divisions that produce much fewer scholarly journal articles, are not linking to publications that mention ML. Third, ML has more prevalence across practically all divisions over time. For further insight, Table 1 lists the top four keywords obtained from recent ACS symposium titles as well as their respective contribution percentage reflected in Figure 1. There, one sees that a handful of keywords can dominate matches in some of the bins, e.g.: ‘electro⁢’, ‘sensor’, ‘protein’, and ‘plastic’. With any ML application, there will be a risk of imperfect data and/or user bias, but this is a useful launch point to appreciate how and where ML is being used in chemical sciences.

2 Motivation for this review

The survey results and literature analysis above showed an opportunity for a tutorial reference to help readers address future research challenges that will require joint applications of CompChem, ML, and chemical and physical intuition (CPI). We will classify concepts in this review using a popular rendition of a “data to wisdom” hierarchy, Figure 2. Scholars have noted shortcomings with similar constructs,20 but we mean for this figure to reflect scientific progress, starting from data and ending in impact. CompChem, ML, and CPI each bring something to the table since all have different strengths and weaknesses. CPI can lead to knowledge, insight, and wisdom from data and information, but the applicability of CPI may be limited when one is faced with large data sets. Alternatively, CompChem is extraordinarily well-suited for generating high quality data that contain useful information (vide infra). Subsequently, non-linear ML models can provide numerical representations of partial or complete sets of data. The application of CPI to these ML models can be used to recognize relationships with respect to chemical and physical concepts to produce knowledge and then ideally insightful predictions. At the time of writing, none of the authors of this review would say that a combination of CompChem, ML, and CPI has reached its full potential of consistently providing wisdom (nor impact), since computational predictions must still be tested by experiment to ensure their validity.

ML brings many reasons to be optimistic about increased knowledge and insights. ML models are extremely well-suited for recognizing and accurately quantifying non-linear relationships (vide infra), a task that is especially difficult for even the most expert-level CPI alone. However, useful ML requires robust datasets that can be provided by CompChem as long as the CPI component is selecting and correctly interpreting appropriate methods for the task at hand. We furthermore note that the knowledge generation process shown in Fig. 2 is by no means a linear one — on the contrary, it contains many loops and dead ends. As we show later, within the troika of CompChem+ML+CPI, ML acts as a catalyst since it can be used in place of explorative data-driven hypotheses generation. Automatically generated hypotheses are then validated and calibrated with CompChem and CPI to yield further improved ML modeling (enriched by more physical prior knowledge), which then loops back with improved hypotheses. This feedback loop is the key to the modern knowledge discovery leading to insight, wisdom and hopefully a positive impact to society.

CompChem and Notable Intersections with ML

We consider quantum mechanics as described by the non-relativistic time-independent Schrödinger equation as our “standard model” because it accurately represents the physics of charged particles (electrons and nuclei) that make up almost all molecules and materials. Indeed, this opinion has been held by some for almost a century:

The fundamental laws necessary for the mathematical treatment of a large part of physics and the whole of chemistry are thus completely known, and the difficulty lies only in the fact that application of these laws leads to equations that are too complex to be solved. – P. A. M. Dirac, 1929

Any theoretical method for predicting molecular and/or material phenomena must first be rooted in quantum mechanics theory and then suitably coarse-grained and approximated so that it can be applied in a practical setting. CompChem, or more precisely, computational quantum chemistry defines computationally-driven numerical analyses based on quantum mechanics. In this section we will explain how and why different CompChem methods capture different aspects of underlying physics. Specifically, this section provides a concise overview of the broad range of CompChem methods that are available for generating datasets that would be useful for ML-assisted studies of molecules and/or materials.

A traditional example of a “good model” is the ideal gas equation: $PV=nRT$, which can be considered ‘simple’, ‘useful’, and ‘insightful’.21 The ideal gas equation relates macroscopic pressure ($P$), volume ($V$), number of molecules ($n$), and temperature ($T$) of gases under idealized conditions, without requiring explicit knowledge of the processes occurring on an atomic scale. Its simple functional form needs just one parameter, the ideal gas constant $R$, and this makes it possible to formulate useful insights, such as how at constant pressure a gas expands with rising temperature. On the other hand, this elegant equation only holds for conditions where the gas behaves as an ideal gas. The derivation of more accurate models of gases requires more mathematically complicated equations of state that rely on more free parameters22 that in turn obfuscate physical insights, require more computational effort to solve, and thus make the model less “good”. This example also offers a convenient connection to ML models that will be discussed later in Section 3. As mathematical models for complex phenomena become more complicated and less intuitive to derive, ML models that can infer non-linear relationships from data become more applicable when increasing amounts of empirical data become available.

Alternatively, the conventional CompChem treatment entails first determining the system’s relevant geometry and its total ground state energy, and from that physical properties of interest (e.g. pressure, volume, band gap, polarizability, etc) can be obtained using quantum and/or statistical mechanics. In this Section, we discuss the relevant CompChem methods for these. While the mathematical physics for these methods might occasionally be too complicated for a user to fully understand, many algorithms exist so that they can still be easily run in a ‘black-box’ way with modern computational chemistry software and accompanying tutorials.23, 24, 25, 26 CompChem thus serves as an invaluable tool to generate data and information for knowledge and insights across many length and time scales. Fig. 3 is an adaption of a multiscale hierarchy24 of different classes of CompChem methods, showing their applicability for modeling different length and time scales, and a depiction of how large scale models may be developed based on smaller scale theories.

1.2 CompChem representations

Integral to every CompChem study is the user’s representation for the system, i.e. how the user chooses to describe the system. CompChem representations can range from simple and lucid (e.g. a precise chemical system such as a water molecule isolated in a vacuum) to complex and ambiguous (e.g. a putative but unproved depiction of a solid-liquid interface under electrochemical conditions). Approximate wavefunctions (expressed on a basis set of mathematical functions) or approximate Hamiltonians (referred to as levels of theory) as described below in this section can also be considered representations. One might then say that many representations for different components of a system will constitute an overall representation, and this is true. The point we make is that the validity of any computational result depends on the overall representation, and sometimes an incorrect representation may provide a correct result due to “fortuitous error cancellation”. In CompChem studies, a valid representation is one that captures the nature of the physical phenomena of a system. For a molecular example, if one is determining the bond energy of a large biodiesel molecule using CompChem methods,27 it may or may not be justified to approximate a nearby long-chain alkyl group (-CnH(2n+1)) simply as a methyl (-CH3) or even a hydrogen atom. Indeed, choosing such a representation can sometimes be a useful example of CPI since alkyl bonds usually exhibit relatively short-ranged interactions (a feature that will be discussed in the context of ML in more detail in Section 4.1.3). An atomic scale geometry with fewer atoms would reduce the computational cost of the study or allow a more accurate but more computationally expensive calculation to be run. On the other hand, it might also be a poor choice if the chemical group, e.g. a substituted alkyl group participated in physical organic interactions such as subtle steric, induction, or resonance effects.28 For a solid-state example, a user might exercise good CPI by assuming that a relatively small unit cell under periodic boundary conditions would capture salient features of a bulk material or a material surface (as is often the case for many metals). On the other hand, subtle symmetry-breaking effects in materials (e.g. distortions arising from tilting octahedra groups in perovskites29, or surface reconstruction phenomena that occur on single crystals)30 might only be observed when considering larger and more computationally expensive unit cells. Relevant to both examples, it may also be that the CompChem method itself brings errors that obfuscate phenomena that the user intends to model. In general, CompChem errors may be due to errors in the initial set up the CompChem application, or if they were due to errors in how the CompChem method is treating the physics of the system, and both factors reflect the representation used in the study. In Section 3, we will discuss how the choice of ML representation also plays similarly critical roles in determining whether and to what extent an ML model is useful.

1.3 Method accuracy

The quantitative accuracy of a CompChem model stems from its suitability in describing the system. As explained above, an observed accuracy will depend on the representation being used. High quality CompChem calculations have traditionally been benchmarked against datasets that consist of well-controlled and relatively precise thermochemistry experiments on small, isolated molecules.31, 32 The error bars for standard calorimetry experiments are approximately 4 kJ/mol (or 1 kcal/mol or 0.04 eV), and computational methods that can provide greater accuracy than this are stated as achieving ‘chemical accuracy’. Note that this term should be used when describing the accuracy of the method compared to the most accurate data possible; for example, if one CompChem method was found to reproduce another CompChem method within 1 kJ/mol, but both methods reproduce experimental data with errors of 20 kJ/mol, then neither method should be called chemically accurate. There are many well-established reasons why CompChem models can bring errors. For example, errors may be due to size consistency33 or size extensivity34 problems that are intrinsic within the CompChem method, larger systems sometimes embody significant medium and long range interactions (e.g. van der Waals forces)35 or self-interaction errors36 that might not be noticeable in small test cases. The recommended path forward is to consider which fundamental interactions are in play in the system of interest, and then to use a CompChem model that is adequate at describing those interactions. Besides this, users should make use of existing tutorial references that provide practical knowledge about which parameters in a CompChem calculation should be carefully noted, for example Ref. 37. Historically the most popular CompChem methods for molecular and materials modeling (the B3LYP38 and PBE39 exchange correlation functionals, see section 2.2.3) are often said to have an expected accuracy of about 10-15 kJ/mol (or 2-4 kcal/mol or 0.1-0.2 eV) when modeling differences between the total energies of two similar systems, and errors are expected to be somewhat larger when considering transition state energies. Though this is used as a simple rule, it is obviously an oversimplification and actual accuracy is only assessed by thoughtful benchmarking of the case being considered.40, 41, 42, 43, 44

1.4 Precision and reproducibility

In CompChem, one normally assumes that that any two users using the same representation for the system with the same code on the same computing architecture will obtain the exact same result within the numerical precision of the computers being used. This is not always the case, especially for molecular dynamics (MD) simulations that often rely on stochastic methods.45 Computational precision also becomes more concerning when there are different versions of codes in circulation, errors that might arise from different compilers and libraries, and a lack of consensus in the community about which computational methods and which default settings should be used for specific application systems, e.g. grid density selections,46 or standard keywords for molecular dynamics simulations.45, 47 There have been efforts to confirm that different codes can reproduce energies for the same system representation,47, 48 but some commercial codes hold proprietary licenses that restrict publications that critically benchmark calculation accuracy and timings across different codes. A path forward to benefit the advancement of insight is the development of (open) source codes49 that perform as well if not better than commercial codes. While increased access to computational algorithms is beneficial, it also raises the need for enforcing high standards of quality and reproducibility.50, 51 We are also glad to see active developments to more lucidly show how any set of computational data is generated, precisely with which codes, keywords, and auxiliary scripts and routines.52, 53, 54, 55 We are now in an era where truly massive amounts of data and information can be generated for CompChem+ML efforts. To go forward, one needs to know what constitutes good and useful data, and the next section provides an overview of how to do this using CompChem.

2 Hierarchies of methods

Earlier we mentioned that a usual task in CompChem is to calculate the ground state energy of an atomic scale system. Indeed, CompChem methods can determine the energy for a hypothetical configuration of atoms, and this constitutes the potential energy surface (PES) of the system (Fig. 4). The PES is a hypersurface spanning $3N$ dimensions, where $N$ is the number of atoms in the system. Since the PES is used to analyze chemical bonding between atoms within the system, the PES can also be simplified by ignoring translational and rotational degrees of freedom for the entire system. This reduces the dimensionality of the PES from $3N$ to $3N-5$ for linear systems (e.g. diatomic molecules or perfectly linear molecules such as acetylene) or $3N-6$ for all other non-linear systems. Furthermore, since visualization is difficult beyond three dimensions, PES drawings will show a 1-D or 2-D projection of this hypersurface where the $z$-axis is conventionally used to represent the scale for system energy.

Any arbitrary PES will contain several interesting features. Minima on the PES correspond to mechanically stable configurations of a molecule or material, for example reactant and product states of a chemical reaction or different conformational isomers of a molecule. Because they are minima, the second derivative of the energy given by the PES with respect to any dimension will be positive. Minima can also be connected by pathways, which indicate chemical transformations (Fig. 4, red line). Along such pathways, the second derivative can be positive, zero, or negative, but all other second derivatives must be positive. Transition states are first order saddlepoints, and thus represent a maximum in one coordinate and a minimum along all others. They correspond to the lowest energy barriers connecting two minima on the PES and are hence important for characterizing transitions between PES minima (e.g. chemical reactions). Second order saddle points56 and bifurcating pathways57 can also exist, but these are not discussed further here.

A wide range of higher-level properties of the system can be predicted and/or derived using the PES, including predicted thermodynamic binding constants, kinetic rate constants for reactions, or properties based on dynamics of the system. The task then to choose an appropriate CompChem method that can carry out energy and gradient calculations on the system’s PES. Fig. 5 shows several different hierarchies for CompChem methods capable of doing this. Note that all of these methods mentioned in this figure fall in the categories of the bottom two regions in the multiscale hierarchy Fig. 3. All of these methods in principle could be used to develop coarse-grained or continuum models as well. Also note that methods in Fig. 5 will bring very different computational costs and opportunities for methods involving ML.

In standard computational quantum chemistry, a system’s energy can be computed in terms of the Schrödinger equation.61, 62, 63 The wavefunction that will be used to represent the positions of electrons and nuclei in the system ($\Psi(\mathbf{r},\mathbf{R})$) is hard to intuit since it can be complex valued. However, its square describes the real probability density of the nuclear ($\mathbf{R}$) and electronic positions ($\mathbf{r}$). In a real system, the position and interactions of a single particle in the system with respect to all other particles will be correlated, and this makes exactly solving the Schrödinger equation impossible for almost all systems of practical interest. To make the problem more tractable, one may exploit the Born-Oppenheimer approximation;64 since nuclei are expected to move much slower than the electrons they can be approximated as stationary at any point along the PES. This allows the energy to be calculated using the time-independent Schrödinger equation and solving the eigenvalue problem:

A second common approximation is to expand the total electronic wavefunction in terms one-electron wave functions (i.e. spin orbitals): $\phi(\mathbf{r}_{i})$. Electrons are fermions and therefore exhibit antisymmetry, which in turn results in the Pauli exclusion principle. Antisymmetry means that the interchange of any two within the system should bring an overall sign change to the wavefunction (i.e. from $+$ to $-$, or vice versa). This property is conveniently captured mathematically by combining one electron spin orbitals into the form of a Slater determinant:

Note that a determinant’s sign changes whenever two columns or rows are interchanged, and in a Slater determinant this corresponds to interchanging electrons and thus the physically appropriate sign change for the overall wavefunction. Additionally, $\frac{1}{\sqrt{n!}}$ is a normalizing factor to ensure the wavefunction is unitary.

The spin orbitals can be treated as a mathematical expansion using a basis set of $\mu$ functions $\chi_{\mu}$, each having coefficients $c_{\mu i}$, which are generally Gaussian basis functions,67, 68, 69 Slater-type hydrogenic orbitals,70 or plane waves under periodic boundary conditions: 71, 72, 73

The different types of mathematical functions bring different strengths and weaknesses, but these will not be discussed further here. A universal point is that larger basis sets will have more basis functions and thus give a more flexible and physical representation of electrons within the system. On one hand this can be crucial for capturing subtle electronic structure effects due to electron correlation. On the other hand, larger basis sets also necessitate significantly higher computational effort. A standard technique to avoid high computational effort in electronic structure calculations is to replace non-reacting core electrons with analytic functions using effective core potentials (ECPs, i.e. pseudopotentials).74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89 This requires reformulating the basis sets that describe the valence space of the atoms, for example see Refs. 90, 91. Larger nuclei that bring higher atomic numbers and larger numbers of electrons will also exhibit relativistic effects,92 and relativistic Hamiltonians are based on the Dirac equation93, 94 or quantum electrodynamics.95 These methods can range from reasonably cost-effective methods96, 97 to those bringing extremely high computational cost.98 Practical applications have traditionally used standard non-relativistic Hamiltonian methods along with ECPs (or pseudopotentials) that have been explicitly developed to account for compressed core orbitals that result from relativistic effects.

Using the Born-Oppenheimer approximation (Eq. 2.2.1) together with a Slater determinant wavefunction (Eq. 3) expressed in a finite basis set (Eq. 4) brings about the simplest wavefunction based method, the Hartree-Fock (HF) approach (for historical context see Refs. 99, 100, 101). The HF method is a mean field approach, where each electron is treated as if it moves within the average field generated by all other electrons. It is generally considered inaccurate when describing many chemical systems, but it continues to serve as a critical pillar for CompChem electronic structure calculations since it either establishes the foundation for all other accurate methods or provides energy contributions (i.e. exact exchange) that is not provided in some CompChem methods. CompChem methods that achieve accuracy higher than HF theory are said to contain electron correlation, a critical component for understanding molecules and materials (as described in more detail in section 2.2.2). Expressing $\Psi$ as a Slater determinant and rearranging Eq. 2.2.1 while temporarily neglecting nuclear-nuclear interactions allows one to define the HF energy in terms of integrals of the electronic spin orbitals:

where the first two terms are referred to as one-electron integrals and represent the kinetic energy due to the electrons and the potential energy contributions from electron-nuclei interactions. The remaining terms are two-electron integrals that describe the potential energy arising from electron-electron interactions and are called Coulomb and exchange integrals. Using Lagrange multipliers, one can express the Hartree–Fock equation in a compact matrix form, the so-called Roothan–Hall equations,102, 103, 104 which allow for an efficient solution:

Each matrix has a size of $\mu\times\mu$, where $\mu$ the number basis functions used to express the orbitals of the system. $\mathbf{C}$ is a coefficient matrix collecting the basis coefficients $c_{\mu i}$ (see Eq. 4), while $\mathbf{S}$ is the overlap matrix measuring the degree of overlap between individual basis functions and $\boldsymbol{\epsilon}$ is a diagonal matrix of the spin orbital energies. Finally, $\mathbf{F}$ is the Fock matrix, with elements of a similar form as in Eq. 5, but expressed in terms of basis functions $\chi_{\mu}$. One important detail not readily apparent in Eq. 6 is that the Fock matrix depends on the orbital coefficients that must be provided before Eq. 6 can be solved. As such, Eq. 6 cannot be solved in closed form, but instead requires a so-called self-consistent field approach. Starting from an arbitrary set of trial (i.e. initial guess) functions, one iteratively solves for optimal molecular orbital coefficients which are then used to construct a new Fock matrix, until a minimum energy is reached in accordance with the variational principle of quantum mechanics. Evaluating and transforming the two-electron integrals in Eq. 5 are a significant bottleneck for these calculations and thus the computational effort of the HF methods formally scales as $\mathcal{O}(\mu^{4})$ with the number of basis functions. This means, that a calculation on a system twice as large will require at least $2^{4}=16$ times as much computing time. The electronic exchange interaction resulting from the antisymmetry of the wavefunction imposes a strong constraint on the mathematical form of ML models for electronic wavefunctions. Construction of efficient and reliable antisymmetric ML models for the many-body wavefunction is an important area of current research. 105, 106

2.2 Correlated wavefunction methods

The system’s correlation energy is defined as sum of electron-electron interactions that originate beyond the mean-field approximation for electron-electron interactions that is provided by HF theory. While correlation energy makes up a rather small contribution to the overall energy of a system (usually about 1% of the total energy), because internal energies in molecular and material systems are so enormous, this contribution becomes rather significant. As an example, most molecular crystals would be unstable as solids if calculated using the HF level of theory. The missing component is attractive forces that are obtained from levels of theory that account for correlation energy. Correlation energies are obtained by calculating additional electron-electron interaction energies that arise from different arrangements of electron configurations (i.e. effectively, different possible excited states) that are not treated with the mean field approach of HF theory.

The most complete correlation treatment is the full configuration interaction (FCI) method, which is the exact numerical solution of the electronic Schrödinger equation (in the complete basis limit) that considers interactions arising from all possible excited configurations of electrons. The FCI wavefunction takes the form of a linear combination of all possible excited Slater determinants which can be generated from a single HF reference wavefunction by electron excitations:

where $\Psi^{\alpha}_{\beta}$ represents the Slater determinant obtained by exciting an electron from orbital $\alpha$ into an unoccupied orbital $\beta$ and the $a$s are expansion coefficients determining the weight of the different contributing configurations. Expectedly, FCI calculations scale extremely poorly with the number of electrons in the system ($\mathcal{O}(n!)$), as the number of possible configurations grows rapidly, making them feasible only for small molecules. For an example of the state of the art, FCI calculations have been used to benchmark highly accurate methods on calculations on a benzene molecule.107

Most correlated wavefunction methods use a subset of the possible configurations in Eq. 7 to be computationally tractable. The configuration interaction (CI)108 method for example only includes determinants up to a certain permutation level (e.g. ‘s’ingle and ‘d’ouble excitations in CISD). Alternatively, MP$n$34 (e.g. MP2) recovers the correlation energy by applying different orders of perturbation theory. Coupled cluster theory, another widely used post-HF method, includes additional electron configurations via cluster operators.109 One coupled cluster method that involves single, double, and perturbative triples excitations, CCSD(T), is referred to as the “gold-standard” approach for CompChem electronic structure methods since it brings high accuracy for molecular energies. However, there are many newer advances that improve upon CCSD(T).110, 107 Note that just because a method has a reputation for being accurate does not mean that it will be for all systems. For example, consider again the benzene molecule, which is best illustrated having dotted resonance bond depicting a planar molecule with equal C–C bond lengths. Such a geometry will not be found to be stable with many different CompChem methods, in part because of subtle chemical bonding interactions and/or errors that arise from specific choices of basis sets used with different levels of theory.111, 112

A key point to reiterate is that correlated wavefunction methods are founded on the HF theory, and so they are even more computationally demanding than HF calculations, e.g. $\mathcal{O}(n^{5})$ for MP2, $\mathcal{O}(n^{6})$ for CCSD and CISD and $\mathcal{O}(n^{7})$ for CCSD(T). However, this computational expense is alleviated by continually improving computing resources (e.g. the usability of graphics processing units (GPUs)) 113, 114, 115, 116 and the development of efficiency enhancing algorithms such as pseudospectral methods,117, 118, 119 resolution of the identity (RI),120 domain-based local pair natural orbital methods (DLPNO),121 and explicitly correlated R12/F12 methods.122 There are also ongoing efforts to develop other CompChem methods based on quantum Monte Carlo123 and density matrix renormalization group theory (DMRG)124 to provide high accuracy with competitive scaling with other computational methods. Efforts are beginning to become implemented that use ML to accelerate these types of calculations.125, 126, 127, 128, 129, 105, 106

Schemes have also been developed to exploit systematic errors between different levels of theory with different basis sets so that approximations can be extrapolated toward an exact result. Examples include the complete basis set (CBS),130 Gaussian G$n$,131 Weizmann (W-)$n$132 methods, and high accuracy extrapolated ab initio thermochemistry (HEAT)133 methods. For a recent review on these and other methods see Ref. 134. These schemes are also becoming a target of recent work using ML methods.135

HF determinants provide good baseline approximations of the ground state electronic structure of many molecules, but they may describe poorly more complicated bonding that arises during bond dissociation events, excited states, and conical intersections.136, 137, 138, 139 Some many-body wavefunctions are best described as a superposition of two or more configurations, e.g. when other configurations in Eq. 7 can have similar or higher expansion coefficients $a$ than the HF determinant. For this reason, high quality single reference methods like CCSD(T) fail because the theory assumes that salient electronic effects are captured by the initial single HF configuration. (In fact, methods such as CCSD(T) have been implemented with diagnostic approaches available that let users know when there may be cause for concern).140, 141, 142 In these cases, it may no longer be trivial to find reliable black-box or automated procedures (e.g. in situations involving resonance states, chemical reactions, molecular excited states, transition metal complexes, and metallic materials, etc.).136 So-called multiconfiguration approaches,136 such as the generalized valence bond (GVB) method143 or the complete active space self consistent field (CASSCF),144 the multireference CI (MRCI) methods,145 complete active space perturbation theory (CASPT2),146 or multireference coupled cluster (MRCC),147, 148 can more physically model these systems since they employ several suitable reference configurations with different degrees of correlation treatments. These methods are not black-box and should be expected to require an experienced practitioner with CPI to choose the reference states that can substantially influence the quality of results.149 This is an area though where ML can bring progress in automating the selections of physically justified active spaces.129

In closing, there are a large number of available correlated wavefunction methods but many are even more costly than HF theory by virtue of requiring an HF reference energy expression shown in Eq. 5. Fig. 5a depicts a so-called ‘magic cube’ (that is an extension beyond a traditional ‘Pople diagram’150, 135) that concisely shows a full hierarchy of computational approaches across different Hamiltonians, basis sets, and correlation treatment methods. This makes it easy to identify different wavefunction methods that should be more accurate and more likely to provide useful atomic scale insights (as well as those that would be more computationally intensive). Another important aspect highlighted in the ‘magic cube’ is that higher level wavefunction methods require larger basis sets to successfully model electron correlation effects. A CCSD(T) computation carried out with a small basis set for example might only offer the same accuracy as MP2 while being two orders of magnitude more expensive to evaluate.108 As was mentioned earlier with the benzene system, spurious errors with different basis sets might still be found that indicate problems with specific combinations of levels of theory and basis sets. The deep complexity of correlated wavefunction methods makes this a promising area for continued efforts in CompChem+ML research.

2.3 Density-functional theory

Density-functional theory (DFT)151 is another method to calculate the quantum mechanical internal energy of a system using an energy expression that relies on functionals (i.e., a function of a function) of electronic density $\rho$:

Compared to wavefunction theory, DFT should be far more efficient since the dimensionality of a density representation for electrons will always be in three rather than the 3$n$ dimensions for any $n$-electron system described by a many-body wavefunction method. DFT brings an important drawback that the exact expression for the energy functional is currently unknown, all approximations bring some degree of uncontrollable error, and this has precipitated disagreeable opinions from purists in chemical physics, especially those who are developing correlated wavefunction methods. However, there is also substantial evidence that DFT approximations are reasonably reliable and accurate for many practical applications that bring information, knowledge, and sometimes insight. We now provide a bird’s-eye view of DFT-based methods.

One thrust of DFT developments since its inception has focused on designing accurate expressions strictly in terms of a density representation, and these approaches are referred to as ‘kinetic energy (KE-)’ or ‘orbital-free (OF-)’ DFT.152 Some energy contributions (e.g. nuclear-electron energy and classical electron-electron energy terms) can be expressed exactly, but other terms such as the kinetic energy as a function of the density are not known and must be approximated. OF-DFT is very computationally efficient (these methods should scale linearly with system size153, 154) but these formulations have not yet been developed to rival the accuracy or transferability of wavefunction methods, though they have been used for studying different classes of chemical and materials systems.155, 156, 157 OF-DFT methods are also used in exciting applications modeling chemistry and materials under extreme conditions.158, 159, 160 One should expect that once highly accurate forms are developed and matured, accurate CompChem calculations on electronic structures on systems having more than a million of atoms might become commonplace. Indeed, there are efforts to use ML to develop more physical OFDFT methods.161, 162

The most commonly used form of DFT (which is also one of the most widely used CompChem methods in use today) is called Kohn–Sham (KS-)DFT.163 In KS-DFT, one assumes a fictitious system of non-interacting electrons with the same ground state density as the real system of interest. This makes it possible to split the energy functional in Eq. 8 into a new form that involves an exact expression of the kinetic energy for non-interacting electrons:

The last two correction terms in Eq. 9 arise from electron interactions, and these are combined into the so-called ’exchange-correlation’ term ($E_{\text{xc}}$), which uniquely defines which scheme of KS-DFT is being used. In theory, an exact $E_{\text{xc}}$ term would capture all differences between the exact FCI energy and the system of non-interacting electrons for a ground state.

The KS-DFT equations can be cast in a similar form as the Roothan–Hall equations (Eq. 6), which allows for a computationally efficient solution. Moreover, the elements of the Kohn-Sham matrix (which replaces the Fock matrix $\mathbf{F}$) are easier to evaluate due to the fact that several of computationally intensive integrals are now accounted for via $E_{\text{xc}}$. Hence, the formal scaling for KS-DFT is $\mathcal{O}(n^{3})$ with respect to the number of electrons. Even though this is much poorer scaling than ideally linear scaling OF-DFT, the exact treatment of non-interacting electrons makes KS-DFT more accurate. Furthermore, there are several modern exchange-correlation functionals that routinely achieve much higher accuracy than HF theory with less computational cost, and thus KS-DFT is a competitive alternative with many correlated wavefunction methods in many modern applications.

A remaining problem is constructing a practical expression for the exchange-correlation functional, as its exact functional form remains unknown. This has spawned a wealth of approximations that have been founded with different degrees of first principles and/or empirical schemes. Classes of KS-DFT functionals are defined by whether the exchange-correlation functional is based on just the homogeneous electron gas (i.e. the ‘local density approximation’, LDA), that and its derivative (i.e. the ‘generalized gradient approximation’, GGA), as well as other additional terms that should result in physically improved descriptions and/or error cancellations. The resulting hierarchy of KS-DFT functionals is often referred to as a ‘Jacob’s Ladder’ of DFT (Fig. 5b). Generally, the higher up the ladder one goes, the more accurate but more computationally demanding the calculation.164 While the intrinsic inexactness in DFT makes it difficult to assess which functionals are physically better than others.165, 166 The Jacob’s Ladder hierarchy is useful for clearly designating how and why newer methods should perform in specific applications (for perspective see Refs. 167, 168, 169), though there remains substantial benefit to those that bring CPI to the development efforts.

Indeed, by being based on a ground-state representation for homogeneous electron gas, DFT calculations can sometimes bring more easily physical insight into some systems that are very challenging for wavefunction theory to examine (e.g. metals170, 171 where HF theory provides divergent exchange energy behaviors). On the other hand, DFT is also generally not well-suited for studying physical phenomena involving localized orbitals or band structures such as those found in semiconducting materials with small band gaps, molecular or material excited charge transfer states, or interaction forces that can arise due to excited states, e.g. dispersion (or London) forces. The former features can normally be treated using Hubbard-corrected DFT+U models that require a system-specific $U-J$ parameter172, 173 or more generalizable but much more computationally expensive hybrid DFT approaches. Dispersion forces (i.e. van der Waals interactions) are non-existent in semilocal DFT approximations, and it is now commonplace to introduce them into DFT calculations using a variety of different methods.35

There is also growing interest in using embedded QC calculation schemes that can partition systems into discrete regions that could be treated with highly accurate correlated wavefunction theory and computationally efficient KS-DFT schemes separately.174, 175, 176, 177, 178 DFT has also been extended to the modeling of excited states in the form of time-dependent (TD-)DFT.179 Similar to ground state DFT, TDDFT is a computationally inexpensive alternative to excited state wavefunction-based methods. The approach yields reasonable results where excitations induce only small changes in the ground state density, e.g. low lying excited states.179, 180 However, due to its single reference nature, TDDFT tends to break down in situations where more than one electronic configuration contribute significantly to the excited state. Just as with correlated wavefunction methods, there are already signs of CompChem+ML efforts to improve the applicability of DFT-based methods.181, 182, 183, 184, 185

2.4 Semi-empirical methods

Correlated wavefunctions and, to a lesser degree, KS-DFT are still very computationally demanding and only of limited use for large scale simulations. Further approximations based on wavefunctions and DFT methods have been developed to simplify and accelerate energy calculations. These so-called semi-empirical methods still explicitly consider the electronic structure of a molecule, but in a more approximate way than methods described above.

Semi-empirical approaches based on wavefunction theory include methods like extended Hückel theory and neglect of diatomic differential overlap (NDDO).186 Both approaches are simplifications of the Hartree–Fock equations (Eq. 5) by introducing approximations to the different integrals. In the NDDO approach,187 only the two-electron integrals in Eq. 5 are considered, where the two orbitals on the right and left hand side of the $\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}$ operator are located on the same atom. The remaining two-center (and one-center) integrals are then approximated by introducing a set of empirical functions, one for each unique type of integral. Moreover, the overlap matrix in Eq. 6 is assumed to be diagonal, which greatly simplifies the energy evaluation. This reduces the required computational effort tremendously and allows the scaling of these approaches to be reduced to $\mathcal{O}(N^{2})$. NDDO serves as a basis for more sophisticated semi-empirical schemes, such as AM1,188 PM7189 and MNDO,190 where the energy is usually determined self-consistently using a minimally sized basis set. Inadequacies in theory can be compensated by different empirical parametrization schemes that can allow these calculations to rival the accuracy of higher level theory for some systems. For example Dral et al.191 provided a recent ‘big-data’ analysis of the performance of several semi-empirical methods with large datasets.

Semi-empirical schemes are also carried over to approximate KS-DFT with so-called density functional tight binding (DFTB).192 DFTB simplifies the Kohn–Sham equations (Eq. 2.2.3) by decomposing the total electron density $\rho$ into a density of free and neutral atoms $\rho_{0}$ and a small perturbation term $\delta\rho_{0}$ ($\rho=\rho_{0}+\delta\rho_{0}$). Expanding Eq. 2.2.3 in the perturbation $\delta\rho_{0}$ makes it possible to partition the total energy into three terms amendable to different approximation schemes:

2.5 Nuclear quantum effects

The quantum nature of lighter elements such as H–Li and even heavier elements that form strong chemical bonds (C–C bond in graphene for example 199) gives rise to significant nuclear quantum effects (NQEs). Such effects are responsible for large differences from the Dulong-Petit limit of the heat capacity of solids, isotope effects, and the deviations of the particle momentum distribution from the Maxwell-Boltzmann equation. 200 To capture NQEs, path-integral molecular dynamics (PIMD) 201, 202 or centroid molecular dynamics (CMD) 203, 204 can be used, but these methods are associated with much higher computational costs (usually about 30 times higher) compared with classical MD simulations using point nuclei. Moreover, because systems may be influenced by competing NQEs, the extent of NQEs is sensitive to the potential energy surface assumed. (Semi-)local DFT approaches may not even qualitatively predict isotope fractionation ratios, and usually hybrid DFT is needed to reach quantitative accuracy. 205 However, employing hybrid DFT calculations or beyond in PIMD/CMD simulations can accrue extremely high computational costs. For this reason, ML force fields have been proposed as efficient means to carry out PIMD simulations, enabling essentially exact quantum-mechanical treatment of both electronic and nuclear degrees of freedom, at least for small molecules with dozens of atoms. 206, 207

2.6 Interatomic Potentials

Interatomic potentials introduce an additional level of abstraction compared to methods described above. Instead of using exact quantum mechanical expressions to create the PES for the system, analytic functions are used to model a pre-supposed PES that contains explicit interactions between atoms, while electrons are treated in an implicit manner (sometimes using partial charge schemes).251, 252, 253, 254, 255, 256 Interatomic potentials thus are (often times dramatically) more computationally efficient than correlated wavefunction, DFT, and semi-empirical approaches. This efficiency makes it possible to study even larger systems of atoms (e.g. biomolecules, surfaces, and materials) than is possible with other computational methods. Note that different empirical potentials bring substantially different computational efficiencies; for example LJ potentials are more efficient than classical forcefields like AMBER and CHARMM, while those are more efficient than most bond-order potentials such as ReaxFF.245, 246 The degree of efficiency arises from the balance of using accurate and/or physically justified functional forms, approximations, and model parameterizations. There are many different formulations (see Fig. 5c), and we will discuss the most general classes. An overview of the different types of potentials and their features is provided in Tab. 2. For extensive discussions on these methods including semi-empirical approaches, we refer to the extensive review by Akimov and Prezhdo (Ref. 257). An excellent review for interatomic potentials is provided by Harrison et. al. (Ref. 258), and an excellent overview of modern methods can be found in a special issue of J. Chem. Phys.259

The distinctions between different types of forcefields can be blurry sometimes, and we will differentiate categories in ascending complexity. One of the simplest interatomic potentials is the Lennard–Jones potential:260

It models the total energy as the sum of all pairwise interaction between atoms $i$ and $j$ using an attractive and repulsive term depending on the interatomic distance $r_{ij}$. $\varepsilon_{ij}$ modulates the strength of the interaction function, while $\sigma_{ij}$ defines where it reaches its minimum. The Lennard–Jones potential is a prototypical “good model” of interatomic potentials, as it has a sufficiently simple physical form with only two parameters while still yielding useful results.

For covalent systems such as bulk carbon or silicon, just pairwise distances are not sufficient to capture the local coordination of the atoms, and many empirical potentials 212, 213, 261 for these systems were expressed as a function of the pairwise distances and three-body terms within a certain cutoff distance. The pairwise term can take the form of LJ-type, electrostatic, or harmonic potentials, and the three-body term is usually a function of the angles formed by sets of three atoms.

So-called Class I classical force fields introduce a more complicated energy expression:

The first three terms are the energy contributions of the distances ($r_{ij}$), angles ($\theta_{ijk}$) and dihedral angles ($\phi_{ijkl}$) between bonded atoms. Because of this, they are also referred to as bonded contributions. Bond and angle energies are modeled via harmonic potentials, with the $k_{ij}$ and $k_{ijk}$ parameters modulating the potential strength and $\bar{r}_{ij}$ and $\bar{\theta}_{ijk}$ being the equilibrium distances and angles. The dihedral term is modeled with a Fourier series to capture the periodicity of dihedral angles, with $k_{ijkl}$ and $\phi_{ijkl}$ as free parameters. The last two terms account for non-bonded interactions. The long range electrostatics are modeled as the Coulomb energy between charges $q_{i}$ and $q_{j}$ and the van der Waals energy is treated via a Lennard–Jones potential (Eq. 12). In Class I/II force fields, empirical parameters are tabulated for a variety of elements in wide ranges of chemical environments (for example Ref. 262). Parameters for any one system should not necessarily be assumed to transfer well to other systems, and reparametrizations may be needed depending on the application. Different sets of parametrization schemes give rise to different types of classical FFs, with CHARMM,217 Amber,214, 215 GROMOS218, 219, 220 and OPLS221, 222 being a few of many examples.

An extension beyond these FFs are Class II (i.e. “polarizable”) force fields, where the static charges are replaced by environment dependent functions (e.g. AMOEBA263). A significant advantage to Class I and II types of forcefields is that they are computationally efficient, which makes them will suited for molecular dynamics simulations of complex and extended (bio)molecules, such as proteins, lipids or polymers. Implementations of forcefield calculations on GPUs makes these simulations extremely productive.264, 265, 266, 267, 268 A disadvantage of Class I and II types of interatomic potentials is that they rely on predefined bonding patterns to compute the total energy, and this limits their transferability. In general, bonds between atoms are defined at the beginning of the simulation run and cannot change. Furthermore, bonding terms make use of harmonic potentials that are not suitable for modeling bond dissociation.

Reactive potentials, which eschew harmonic potential dependencies and thus can describe the formation and breaking of chemical bonds, include the embedded atom method (EAM, Fig. 5c), which is used widely in materials science.235 EAM is a type of many-body potential primarily used for metals, where each atom is embedded in the environment of all others. The total energy is given by

Another common type of reactive potentials are bond order potentials (BOPs). In general, BOPs model the total energy of a system as interactions between the neighboring atoms:

While efficient and versatile, all interatomic potentials described above are inherently constrained by their functional forms. A different approach is pursued by machine learned potentials (MLPs), such as Behler-Parinello Neural Networks,274 q-SNAP,275 and GAP potentials276 (Fig. 5c). In MLPs, suitable functional expressions for interactions and energy are determined in a fully data-driven manner and ultimately only limited by the amount and quality of available reference data. One can then use substantially more data to generate a much more accurate MLP than would be possible when using, for instance, a ReaxFF potential trained on similar data sets.277

For the sake of completeness, we note that all approaches described here are fully atomistic – each atom is modeled as an individual entity. It is also possible to combine groups of atoms into pseudo-particles giving rise to so-called coarse grained methods. On an even higher level of abstraction, whole environments can be modeled as a single continuum. As such approaches are not subject of the present review, we refer the interested reader e.g. to Refs. 278 and 279.

3 Response properties

Once an energy calculation is completed by one of the CompChem methods above, many other interesting molecular properties can be calculated. Most of these properties can be obtained as the response of the energy to a perturbation, e.g. changes in nuclear coordinates $\mathbf{R}$, external electric ($\boldsymbol{\epsilon}$) or magnetic ($\mathbf{B}$) fields or the nuclear magnetic moments $\{\mathbf{I}_{i}\}$. Given an expression for the energy, which depends on the above quantities, so-called response properties can be computed via the corresponding partial derivatives of the energy. A general response property $\boldsymbol{\Pi}$ then takes the form

A common response property are nuclear forces $\mathbf{F}=-\boldsymbol{\Pi}(1,0,0,0)$ that are the negative first derivative of the energy with respect to the nuclear positions. Such calculations allow a plethora of different geometry optimization schemes for chemical structures on the PES. Hessian calculations corresponding to the second derivative of energy with respect to nuclear positions are necessary to confirm the location of first order saddle points on the PES and identify normal modes and their frequencies for vibrational partition functions that are useful for modeling temperature dependencies based on statistical thermodynamics. Hessian calculations are computationally costly, since they normally involve calculations based on finite differences methods involving many nuclear force calculations. Many methods have been developed to allow CompChem algorithms to sample minimum energy regions of the PES280, 281, 282, 283, 284 or precisely locate points of interest.285, 286 Historically, many of these techniques have relied on approximate or full Hessian calculations,287 but other approaches such as the nudged-elastic band288, 289 and string290, 291, 292 methods are popular alternatives that do not require a Hessian calculation. There have also been efforts using different forms of ML to accelerate procedures or overcome long-standing challenges in efficient sampling of and optimization on the PES.293, 294, 295, 296, 297, 298

The general expression above can provide a wealth of other quantities, some of which are relevant for molecular spectroscopy and/or provide a direct connection to experiment (see Tab. 3).

4 Solvation models

An important aspect of CompChem is molecular descriptions from with a solution environment. Simulating a dynamical environment composed of many surrounding molecules is usually not feasible with electronic-structure methods. To circumvent this problem, solvation modeling schemes have been devised (see Refs. 301, 302, 303, 304, 305, 306 for discussions on this topic).

The most popular approach are so-called polarizable continuum solvent models (PCM).279 They model the electrostatic interaction of a solute molecule with its environment by representing the charge distribution of the solvent molecules as a continuous electric field, the reaction field. This dielectric continuum can be interpreted as a thermally averaged representation of the environment and is typically assigned a constant permittivity depending on the particular solvent to be modeled ($\varepsilon=80.4$ for water). The solute is placed inside a cavity embedded in this continuum. The charge distribution of the molecule then polarizes the continuous medium, which in turn acts back on the molecule. To compute the electrostatic interactions arising from this mutual polarization with electronic structure theory, a self consistent scheme is employed. After constructing a suitable molecular cavity, a Poisson problem of the following form is solved:

$\Gamma$ indicates the surface of the cavity. Eq. 17 is solved numerically to obtain the surface charge distribution $\sigma(\mathbf{s})$. Once $\sigma(\mathbf{s})$ has been determined in this fashion, the potential is computed according to Eq. 19 and used to construct an effective Hamiltonian of the form

where $\hat{H}$ is the vacuum Hamiltonian. These equations are then solved self consistently in a Roothan–Hall or Kohn–Sham approach, yielding the electrostatic solvent-solute interaction energy. This scheme is also called the self-consistent reaction field approach (SCRF).

Continuum models differ in how the cavities are constructed and how Eq. 17 is solved to obtain the surface charge distribution. Variants include the original PCM model, also refered to as dielectric PCM (D-PCM),307 the integral equation formulation of PCM (IEFPCM),308 SMD,309 conductor PCM (C-PCM)310 or the conductor-like screening model (COSMO).311 The latter two approaches replace the dielectric medium by a perfect conductor to allow for a particularly efficient computation of $\sigma(\mathbf{s})$. PCMs can be further extended with statistical thermodynamics treatments to account for solutes having different size and concentration effects, and this leads to models such as COSMO-RS.312

A drawback of most PCM-like approaches is that they neglect local solvent structures. Thus, they cannot reliably account for situations where explicit solvent interactions are important, e.g. when for stabilizing specific sites for a transition state through hydrogen bonding.301 Furthermore, while implicit models might be parametrized to fit bulk-like properties of mixed or ionic solvents (e.g. Ref. 313), the complex local solvent environment presented by these systems are treatable by other means. For mixed solvent systems a range of hybrid schemes such as COSMO-RS,305 reference interaction site models (RISMs)314, 315 or QM/MM316, 317, 318 approaches have been developed. As an in-depth discussion of these alternative schemes exceeds the scope of this review, we instead refer to other references.319, 320

ML models are becoming used to describe solvent effects. Ref. 300 introduces a continuum ML model based on a reaction field that can predict energies and response properties for continuum solvents, it can extrapolate to solvents not seen during training, and it can be extended to operate in a QM/MM fashion to account for explicit solvents effects in a Claisen rearrangement reaction. Ref. 321 implemented automatable calculation schemes and unsupervised ML to allow predictions of single ion solvation energies for monovalent and divalent cations and anions based on physically rigorous quasi-chemical theory.322, 323 Ref. 324 used convolutional neural networks and molecular dynamics simulations to carry out high-throughput screening of mixed solvent systems. Ref. 325 implemented efficient ways to carry out ML-based QM/MM molecular dynamics simulations.

5 Insightful predictions for molecular and material properties

By solving for electronic structures, by whatever means is appropriate, one obtains molecular energies and energy spectrum (typically corresponding to quasiparticles given by Kohn-Sham or Hartree-Fock orbitals). From these, one can then compute molecular or material properties that arise from quantum mechanical and statistical operators, e.g. thermodynamic energies, response properties, highest and lowest occupied molecular orbital energies, band gaps, among other properties. Many properties are defined by the characters of the orbitals, and having knowledge of these should always be helpful and aid in deriving useful insight into designing molecules and materials for a particular function. Furthermore, one is often interested in how these molecules behave over time (i.e. the dynamics given some statistical ensemble that depends on temperature, pressure, etc) over all possible degrees of freedom. By understanding how energies and forces change over time, one can predict thermal and pressure dependencies as well as spectroscopic properties for advanced knowledge that builds toward insightful predictions.

Molecular and materials chemistry is vastly complex and variable, and one often faces a question of whether to span wider chemical spaces versus take deeper explorations of a specific phenomenon. A key problem is that even after the effort of either approach, it is also not as clear how information for one system might be related to another to provide more knowledge. For instance, one may decide to calculate all possible properties of ethanol with a CompChem method, but understanding how any calculated property would be correlated to an analogous property of isopropanol is still usually difficult to do. There is great interest in understanding chemical and materials space through applications of quantitative structure activity/property relationships,326, 327 cheminformatics,328 conceptual DFT,329 and alchemical perturbation DFT.330 All these applications benefit from greater access to CompChem data, and all have promise as being interfaced with ML for transformative applications to catalyze wisdom and impact.

Machine Learning Tutorial and Intersections with Chemistry

Machine learning (ML) has had a dramatic impact on many aspects of our daily lives and has arguably become one of the most far-reaching technologies of our era. It is hard to overstate its importance in solving long standing computer science challenges such as image classification 331, 332, 333, 334 or natural language processing 335, 336, 337, 338, 339 – tasks that require knowledge that is hard to capture in a traditional computer program 340, 341, 342. Previous classical artificial intelligence (AI) approaches relied on very large sets of rules and heuristics, but these were unable to cover the full scope of these complex problems. Over the past decade, advances in ML algorithms and computer technology made it possible to learn underlying regularities and relevant patterns from massive datasets that enable automatic constructions of powerful models that can sometimes even outperform humans at those tasks.

This development inspired researchers to approach challenges in science with the same tools, driven by the hope that ML would revolutionize their respective fields in a similar way. Here, we give an overview of these developments in chemistry and physics to serve as an orientation for newcomers to ML. We will first explain what tasks ML is good at and when it might not be the best solution to a problem. We will start by introducing the field of ML in general terms and dissect its strengths and weaknesses.

In the most general sense, ML algorithms estimate functional relationships without being given any explicit instructions of how to analyze or draw conclusions from the data. Learning algorithms can recover mappings between a set of inputs and corresponding outputs or just from the inputs alone. Without output labels, the algorithm is left on its own to discover structure in the data.

Universal approximators 343, 344 are commonly used for that purpose. These reconstruct any function that fulfills a few basic properties, such as continuity and smoothness, as long as enough data is available. Smoothness is a crucial ingredient that makes a function learnable, because it implies that neighboring points are correlated to $\mathcal{Y}$ in similar ways. That property means that one can draw successful conclusions about unknown points as long as they are close to the training data (coming from the same underlying probability distribution).341 In contrast, completely random processes in the above sense allow no predictions.

An association that immediately springs to mind is traditional regression analysis, but ML goes a step further. Regression analyses aim to reconstruct the function that goes through a set of known data points with the lowest error, but ML techniques aim to identify functions to predict interpolations between data points and thus minimize the prediction error for new data points that might later appear.345 Those contrasting objectives are mirrored in the different optimization targets: In traditional regression the optimization task

only measures the fit to the data, but learning algorithms typically aim to find models $\hat{f}$ that satisfy

Both optimization targets reward a close fit, often using the squared loss $\mathcal{L}\left(\hat{f}(\mathbf{x}),\mathbf{y}\right)=\left(\hat{f}(\mathbf{x})-\mathbf{y}\right)^{2}$. However, the key difference is an additional regularization term in Eq. 22, which influences the selection of candidate models by introducing additional properties that promote generalization. To understand why this is necessary, it is helpful to consider that Eq. 22 is only a proxy for the optimization problem

A model that is heavily regularized (i.e. using a large $\lambda$) will eventually become biased in that it is too simplistic to fit the data well. In contrast, a lack of regularization might yield an overly complex model with high variance. Such an “overly fit” model will follow the data exactly to the point that it also models the noise components and consequently fail to generalize (see Fig. 6). Finding the appropriate amount of regularization $\lambda$ to manage under- and over-fitting is known as attaining a good bias-variance trade-off.351 We will introduce a process called cross-validation to address this challenge further below (see Section 3.4.3).

ML algorithms can infer functional relationships from data in a statistically rigorous way without detailed knowledge about the problem at hand. ML thus captures implicit knowledge from a dataset – even aspects where CPI might not be available. Traditional modeling approaches, like classical forcefields discussed in Section 2.2.6, rely on preconceived notions about the PES it is modeling and thus the way the physical system behaves. In contrast, ML algorithms start from a loss function and a much more general model class. Within the limits permitted by the noise inherent to the data, generalization can be improved to arbitrary accuracy given increasingly larger informative training data sets. This process allows us to explore a problem even before there is a reasonably full understanding. An ML predictor can serve as starting point for theory building and be regarded as a versatile tool in the modelling loop: building predictive models, improving them, enriching them by formal insight, improving further and ultimately extracting a formal understanding. More and more research efforts start to combine data-driven learning algorithms with rigorous scientific or engineering theory to yield novel insights and applications.352, 15, 9

For a quantum chemical property for compounds in a dataset, first principles calculations need to be repeated independently for each input, even if they are very similar. No formally rigorous method exists to exploit redundancies in the calculations in such a scenario. The empiricism of learning algorithms however does provide a pathway to extract information based on compound structure similarity. A data-driven angle allows one to ask questions in new ways and can give rise to new perspectives on established problems. For example, unsupervised algorithms like clustering or projection methods group objects according to latent structural patterns and provide insights that would remain hidden when only looking at individual compounds.

1.2 What is ML not good at?

Some difficult problems in chemistry and physics can be solved accurately with CompChem, but doing so would require significant resources. For example, enumerating all pair-wise interactions in a many-body system will inevitably scale quadratically, and there is no obvious path around this. One might ask if empirical approaches can address such fundamental problems more efficiently, but this is unfortunately not possible since ML is more suited for finding solutions in general function spaces rather than in deterministic algorithms where constraints guide the solution process. However, if we were not as interested in finding a full solution but rather some aspect of it, the stochastic nature of ML can be beneficial. For instance, a traditional ML approach might not be the best tool for explicitly calculating the Schrödinger equation, but it might be a far more useful tool for developing a forcefield that returns the energy of a system without the need for a cumbersome wavefunction and a self-consistent algorithm. As an example, Hermann and Noé 105 used deep neural networks to show how ML methods may be suitable for overcoming challenges faced by traditional CompChem approaches.

ML algorithms require a large amount of high quality data, and it is hard to decide a priori when a data set is sufficient. Sometimes, a data set may be large, but it does not adequately sample all the relevant systems one intends to model, e.g. a molecular dynamics simulation might generate many thousands of molecular confirmations used to train an ML forcefield, but perhaps that sampling only occurred in a local region of the PES. In this case, the ML forcefield would be effective at modeling regions of the PES it was trained to but useless in other regions until more data and broader sampling occurred. This feature is general to all empirical models that are generally limited in their extrapolation abilities.

Standard ML algorithms cannot conceptualize knowledge from a data set. Two of the main reasons are the non-linearity and excessive parametric complexity of most models that allow many equally viable solutions for the same problem.353, 354 It can be hard to gain insight into the modeled relationship because it is not based on a small set of simple rules. Techniques have emerged to make ML models interpretable (explainable AI - XAI355). While helpful, drawing scientific insight clearly still requires human expertise.356, 357, 352, 358, 359, 360, 361, 355 Furthermore, the path from an ML model back to a physical set of equations is being explored, but it is far from being fully established automatically.362, 363, 364, 365, 366, 367, 368

Despite following the rules of best practice, ML algorithms can give unexpected and undesired results. Instead of extracting meaningful relationships, they may occasionally exploit nuisance patterns within the underlying experimental design, like the model architecture, the loss function or artifacts in the dataset. This results in a “clever Hans” predictor,360 which technically manages the learning problem but uses a trivial solution that is only applicable within the narrow scope of the particular experimental setup at hand. The predictor will appear to be performing well, while actually harvesting the wrong information, and therefore not allowing any generalization or transferable insights.

For example, a recently proposed random forest predictor for the success of Buchwald-Hartwig coupling reactions 369 was later revealed to give almost the same performance when the original inputs were replaced by Gaussian noise. 370, 371 This finding strongly suggested that the ML algorithm exploited some hidden underlying structure in the input data, irrespective of the chemical knowledge that was provided through the descriptor. Even though the model might appear quite useful, any conclusions that rely on the importance of the chemical features used in the model were thus rendered questionable at best. This example demonstrates that out-of-sample validation alone is often not sufficient to establish that a proposed model has indeed learned something meaningful. Therefore, the hypothesis described by the model must be challenged in extensive testing in practically relevant scenarios like actual physical simulations. In other words the ML model needs to lead to a better understanding of the modeling itself and the underlying chemistry.

2 Types of learning

ML models are classified by the type of learning problem they solve. Consider for instance a data scientist who develops an ML model that can predict acidity constants (p$K_{\text{a}}$’s) for any molecule. A researcher with knowledge of physical organic chemistry might be aware of the empirical Taft equation28 that provides a linear free energy relationship between molecules on the basis of empirical parameters that account for a molecules fundamental field, inductive, resonance, and steric effects (e.g. values related to Hammett $\rho$ and $\sigma$ values). There are several ways the data scientist might develop an ML model to do this. Examples mentioned here include supervised, unsupervised, and reinforcement learning.

Using the p$K_{\text{a}}$ predictor example, a supervised learning algorithm could be trained to correlate recognizable chemical patterns or structures to experimentally known p$K_{\text{a}}$s. The goal would be to deduce the relationship between these inputs and outputs, such that the model is able to generalize beyond the known training set. A standard universal approximator has to accomplish this learning task without any preconceived notion about the problem at hand and will therefore likely require many examples before it can make accurate predictions. Recently, a lot of research is being carried out that investigates ways to incorporate high-level concepts into the learning algorithm in the form of prior knowledge.372, 207 In this vein, one could take into account chemically relevant parameters such as Hammett constants so that the parameterized ML model incorporates the modified Hammett or Taft equation. An example of a classification problem in materials science is the categorization of materials, where identifying characteristics of the electronic structure can be used to distinguish between insulators and metals.373

2.2 Unsupervised learning

Unsupervised learning describes problems were only the inputs $\mathcal{X}$ are known, with no corresponding labels. In this setting, the goal is to recover some of the underlying structure of the data to gain a higher-level understanding. Unsupervised learning problems are not as rigorously defined as supervised problems in the sense that there can be multiple correct answers, depending on the model and objective function that is applied.

For example, one might be interested in separating conformers of a molecule from a molecular dynamics trajectory, given exclusively the positions of the atoms. A clustering algorithm (like the k-means algorithm) could identify those conformers by grouping the data based on common patterns.374, 375 Alternatively, a projection technique could reveal a low-dimensional representations of the dataset.376 Often data is represented in high dimension, despite being intrinsically low-dimensional. With the right projection technique, it is possible to retain the meaningful properties in a representation with less degrees of freedom. A conceptually simple embedding method is principal component analysis (PCA) in which the relationship that is sought to be preserved is the scalar product between the data points.340 There are many other linear and non-linear projection methods, such as multi-dimensional scaling,377 kernel PCA,378, 379 t-distributed stochastic neighbor embedding (t-SNE),380 sketch-map,381 and the uniform manifold approximation and projection (UMAP).382 Finally, anomaly detection is a further variant of unsupervised learning, where ’outliers’ to the available data can be discovered.383 However, without knowing the labels (in this example, the potential energy associated with each geometry), there is no way to conclusively verify that the result is correct. The literature is gradually seeing more instances of unsupervised learning, particular to reveal important chemical properties to efficiently explore chemical/materials spaces.

2.3 Reinforcement learning

Reinforcement learning describes problems that combine aspects of supervised and unsupervised learning. Reinforcement learning problems often involve defining an agent within an environment that learns by receiving feedback in the form of punishments and rewards. The progress of the agent is characterized by a combination of explorative activity and exploitation of already gathered knowledge.384 For chemistry applications, reinforcement learning techniques are being increasingly used for finding molecules with desired properties in large chemical spaces.9

3 Universal approximators

Universal approximators have their origins in the 1960s, where the hope was to construct “learning machines” that have similar capabilities as the human brain. An early mathematical model of a single simplified neuron emerged that was called a perceptron (Eq. 24). 385, 386

Here, ${\bf x}$ denotes the $N$-dimensional input to the perceptron. It has $N+1$ parameters consisting of $w_{i}$ (so-called weights) and a single $b$ (a so-called threshold) that are adapted to the data. This adaption process is typically called “learning” (vide infra) and it amounts to minimizing a predefined loss function.

In the 1960s, this simple neural network had very limited use, as it was only able to model a linear separating hyperplane. Even simple non-linear functions like the XOR were out of reach.387 Thus, excitement waned but then reappeared two decades later with the emergence of novel models consisting of more neurons and their arrangement in multi-layer neural network structures 388 (see Eq. 25). Recent algorithmic and hardware advances now allow deep and increasingly complex architectures.1, 2

In Eq.(25), $g(\cdot)$ denotes an activation function that is a non-linear transformation that allows complex mappings between input and output. As with the perceptron, the parameters of multi-layer NNs can be learned efficiently using iterative algorithms that compute the gradient of the loss-function using the so-called back-propagation (BP) algorithm.388, 389, 390 In the late 1980s, artificial neural networks were then proven to be universal approximators of smooth nonlinear functions, 343, 391, 392 and so they gained broad interest even outside the ML community that then was still relatively small.

In 1995, a novel technique called Support Vector Machines (SVM) 393, 345 and kernel-based learning were then proposed, 379, 394, 395, 396 which came with some useful theoretical guarantees. SVMs implement a nonlinear predictor:

where $K$ is the so-called kernel. The kernel implicitly defines an inner product in some feature space and thus avoids an explicit mapping of the inputs. This “kernel trick”397 makes it possible to introduce non-linearity into any learning algorithm that can be expressed in terms of inner products of the input.379 It has since been applied to many other algorithms beyond SVMs,394 such as Gaussian Processes (GP),348 PCA,378, 379 independent component analysis (ICA).398

The most effective kernels are tailored to the specific learning task at hand, but there are many generic choices, such as the polynomial kernel $K({\bf x}_{j},{\bf x})=({\bf x}_{j}\cdot{\bf x}-b)^{d}$, which describes inner products between degree $d$ polynomials. Another popular choice is the Gaussian kernel $K({\bf x}_{j},{\bf x})=\exp({-{1/2\sigma^{2}}({\bf x}_{j}-{\bf x})^{2}})$. It is one of the most versatile kernels because it only imposes smoothness assumptions on the solution depending on the width parameter $\sigma$. 347, 395

As seen in Eq. (26), an SVM can also be understood as a shallow neural network with a fixed set of non-linearities. In other words, the kernel explicitly defines a similarity metric to compare data points, whereas neural networks have some freedom to shape this transformation during training because they nest parameterizable non-linear transformations on multiple scales. This difference gives both techniques some unique strengths and drawbacks. Despite that, there exists a duality between both approaches that allows neural networks to be translated into kernel machines and analyzed more formally (see Refs. 399, 400, 401).

In the context of computational chemistry, both NNs and kernel-based methods are the most used ML approaches. Simpler learners, such as nearest neighbors models or decision trees can still be surprisingly effective. Those have also been successfully used to solve a wide spectrum of problems including drug design, chemical synthesis planning, and crystal structure classification.402, 403, 404, 405, 406, 407

4 The ML workflow

In the following, we summarize the overall ML process, starting from a dataset all the way to trained and tested model. The ML workflow typically includes the following stages:

Note, that the progression to a good ML model is not necessarily linear and some steps (except the out of sample test) may require reiteration as we learn about the problem at hand.

On a fundamental level, ML models could be simply regarded as sophisticated parametrizations of datasets. While the architectural details of the model matter, the reference dataset forms the backbone that ultimately determines its effectiveness. If the dataset is not representative of the problem at hand, the model will be incomplete and behave unpredictably in situations that have been improperly captured. The same applies to any other shortcomings of the dataset, such as biases or noise artifacts that will also be reflected in the model. Some of these dataset issues are likely to remain unnoticed when following the standard model selection protocol since training and test datasets are usually sampled from the same distribution. If the sampling method is too narrow, errors seen during the cross-validation procedure may appear to be encouragingly small, but the ML model will fail catastrophically when applied to a real problem. If the training and test sets come from different distributions, then techniques to compensate this covariate shift can be used.408, 409

Robust models can generally only be constructed from comprehensive datasets, but it is possible to incorporate certain patterns into models to make them more data-efficient. Prior scientific knowledge or intuition about specific problems can be used to reduce the function space from which a ML algorithm has to select a solution. If some of the unphysical solutions are removed a priori, less data are necessary to identify a good model. This is why NNs and kernel methods, despite both being broad universal function classes, bring different scaling behaviors. The choice of the kernel function provides a direct way to include prior knowledge such as invariances, symmetries or conservation laws, whereas NNs are typically used if the learning problem cannot be characterized as specifically.410, 372, 207 In general, without prior knowledge NNs often require larger datasets to produce the same accuracy as well-constrained kernel methods that embody problem knowledge. This consideration is particularly important if the data is expensive, e.g. if it comes from high quality experiments and/or expensive computations.

4.2 Descriptors

In order to apply ML, the dataset needs to be encoded into a numerical representation (i.e. features/descriptors) that allow the learning algorithm to extract meaningful patterns and regularities.411, 412, 413, 414, 415, 416, 417, 418, 419 This is particularly challenging for unstructured data like molecular graphs that have well-defined invariable or equivariable characteristics that are hard to capture in a vectorial representation. For example, atoms of the same type are indistinguishable from each other, but it is hard to represent them without imposing some kind of order (which inevitably assigns an identity to each atom). Furthermore, physical systems can be translated and rotated in space without affecting many attributes. Only a representation that is adapted to those transformations can solve the learning problem efficiently.

It turned out to be a major challenge to reconcile all invariances of molecular systems in a descriptor without sacrificing its uniqueness or computability. Some representations cannot avoid collisions, where multiple geometries map onto the same representation. Others are unique, but prohibitively expensive to generate. Many solutions to this problem have been proposed, based on general strategies such as invariant integration,207 parameter sharing,420, 352, 421, 422 density representations,276 or finger printing techniques.423, 424, 425, 426, 427, 428, 429, 430, 431, 432 Alternatively, an NN model infers the representation from data.352, 433, 434, 423 To date, none of the proposed approaches are without compromise, which is why the optimal choice of descriptor depends on the learning task at hand.

4.3 Training

The training process is the key step that ties together the dataset and model architecture. Through the choice of the model architecture, we implicitly define a function space of possible solutions, which is then conditioned on the training dataset by selecting suitable parameters. This optimization task is guided by a loss function that encodes our two somewhat opposing objectives: (1) achieving a good fit to the data, while (2) keeping the parametrization general enough such that the trained model becomes applicable to data that is not covered in the training set (see the two terms in Eq. 22). Satisfying the latter objectives involves a process called model selection in which a suitable model is chosen from a set of variants that have been trained with exclusive focus on the first objective. Depending on the model architecture, more or less sophisticated optimization algorithms can be applied to train the set of model candidates.

are typically linear in their parameters $\vec{\alpha}$ (see Eq. 26). Coupled with a quadratic loss function, $\mathcal{L}(\hat{f}(\mathbf{x}),\mathbf{y})=(\hat{f}(\mathbf{x})-\mathbf{y})^{2}$, they yield a convex optimization problem. Convex problems can be solved quickly and reliably due to only having a single solution that is guaranteed to be globally optimal. This solution can be found algebraically by taking the derivative of the loss function and setting it to zero. For example, KRR and GPs then yield a linear system of the form

which is typically solved in a numerically robust way by factorizing the kernel matrix $\mathbf{K}$. There exist a broad spectrum of matrix factorization algorithms such as the Cholesky decomposition that exploit the symmetry and positive definiteness properties of kernel matrices.435, 436, 437, 438, 439 Factorization approaches are however only feasible if enough memory is available to store the matrix factors, and this can be a limitation for large-scale problems. In that case, numerical optimization algorithms provide an alternative: they take a multi-step approach to solve the optimization problem iteratively by following the gradient:

where $\gamma$ is the step size (or learning rate). Iterative solvers follow the gradient of the loss function until it vanishes at a minimum, which is much less computationally demanding per step, because it only requires the evaluation of the model $\hat{f}$. In particular, kernel models can be evaluated without storing $\mathbf{K}$ (see Eq. 28).

are constructed by nesting non-linear functions in multiple layers, which yields non-convex optimization problems. Closed-form solutions similar to Eq. 27 do not exist, which means that NNs can only be trained iteratively, i.e. analogous to Eq. 28. Several variants of this standard gradient descent algorithm exist including stochastic or mini-batch gradient descent, where only an $n$-sized portion of the training data $(\mathbf{x},\mathbf{y})_{i:i+n}$ is considered in every step. Due to multiple local minima and saddle points on the loss surface, the global minimum is exponentially hard to obtain such that these algorithms usually converge to a local minimum. However, thanks to the strong modelling power of NNs, local solutions are usually good enough.440

In addition to the parameters that are determined when fitting a ML model to the dataset (i.e. the node weights/biases or regression coefficients), many models contain so-called hyper-parameters that need to be fixed before training. Two types of hyper-parameters can be distinguished: ones that influence the model, such as the type of kernel or the NN architecture, and ones that affect the optimization algorithm, e.g. the choice of regularization scheme or the aforementioned learning rate. Both tune a given model to the prior beliefs about the dataset and thus play a significant role in model effectiveness. Hyper-parameters can be used to gauge the generalization behavior of a model.

Hyper-parameter spaces are often rather complex: certain parameters might need to be selected from unbounded value spaces, others could be restricted to integers or have interdependencies. This is why they are usually optimized using primitive exhaustive search schemes like grid or random search in combination with educated guesses for suitable search ranges. Common gradient-based optimization methods typically cannot be applied for this task. Instead, the performance of a given set of hyper-parameters is measured by evaluating the respective model on another training dataset called the validation dataset (see Fig. 6). This process is also referred to as model selection.

Cross-validation or out-of-sample testing is a technique to assess how a trained ML model will generalize to previously unseen data.395, 441 For a reasonably complex model, it is typically not challenging to generate the right responses for the data known from the training set. This is why the training error is not indicative of how the model will fulfill its ultimate purpose of predicting responses for new inputs. Alas, since the probability distribution of the data is typically unknown, it is not possible to determine this so-called generalization error exactly. Instead, this error is often estimated using an independent test subset that is held back and later passed through the trained model to compare its responses to the known test labels. If the model suffers from over-fitting on the training data, this test will yield large errors. It is important to remember not to tweak any parameters in response to these test results, as this will skew this assessment of the model performance and will lead to overfitting on the test set.442

Besides cross-validation, there are alternative ways to estimate the generalization error, for example via maximization of the marginal likelihood in Bayesian inference.443, 444, 445 Some well-defined learning scenarios even allow the computation of rigorous upper bounds for the generalization error.345, 446, 447, 448

Applications of Machine Learning to Chemical Systems

We now discuss ways that CompChem methods described in Section 2 and ML methods in Section 3 can be implemented as CompChem+ML approaches for insights into chemical systems. We often notice the lack of details about why an ML model is used and how it actually contributes to worthwhile and scientific insights. Thus, we will summarize the underlying attributes of conventional CompChem+ML efforts and then explain why these attributes are important for specific applications.

To begin, consider molecules or materials in a dataset, and any entry will be related to another based on an abstract concept of “similarity”. While similarity is an application-dependent concept, it should go hand in hand with CPI. For instance, physical properties of chemical systems can be attributed to the structure and/or composition of the chemical fragments within those systems. Thus, if chemical structures and/or compositions of two entries in the database were similar, then their physical properties would also likely be similar.

For CompChem+ML using a supervised algorithm, a CompChem prediction might be made on a hypothetical system, pinpointed by an ML model that was trained to identify chemical fragments that correlate with labeled physical properties. This would be a direct exploitation of chemical similarity. Alternatively, for CompChem+ML using an unsupervised algorithm, the ML model would identify an underlying distribution or key features based on the similarity between pairs of entries in the dataset without labeled data. This would be a more nuanced leveraging of chemical similarity. In both cases the accuracy, efficiency and reliability of the ML models depends strongly on the how similarity is defined and measured.

In this section we will first describe state-of-the-art descriptors and kernels for atomic systems that can be used to quantify the similarity between chemical systems. We will then explain the essential attributes of good atomic descriptors. Lastly for this section, we will elucidate why and how specific combinations of these descriptors and ML algorithms are beginning to revolutionize the field of CompChem.

In CompChem, molecules and materials are usually represented by the Cartesian coordinates and the chemical elements of all the atoms. Thus, the size of the vector representation containing the coordinates and charges will be $\{\mathcal{R}^{3N}\text{ and }\mathcal{Z}^{N}\}$, respectively, for a system of size $N$. Even though these atomic coordinates provide a complete description of the system, they are hardly ever used as the input of a ML model because this vector would introduce substantial superfluous redundancy. For instance, an ML model might treat two identical molecules that are rotated or translated as different molecules, and that in turn might cause the ML model to predict different physical properties for the two otherwise indistinguishable molecules. There are further difficulties when comparing molecules having different numbers of atoms. To work around these problems, atomic coordinates are usually converted into an appropriate representation $\psi$ that is suitable for a particular task. Such conversions are useful because they allow the incorporation of physical invariances. Mathematically speaking, the representation fulfills

where $S$ indicates a symmetry operation, e.g. a rigid rotation about an axis $C_{i}$, an exchange of two identical atoms, or a translation of the whole system in the Cartesian space, etc. It can also be advantageous to adopt a coarse-grained representation of the system. 449, 450 For example, dihedral angles of a peptide might be accounted for without the positions of the side-chains, positions of ions in a solution might be accounted for without the explicit coordinates of solvents, or just the center of mass for a water molecule might be accounted for in place of the full three-centered atomistic representation. The choice of these coarse-grained representations provides a way to incorporate prior knowledge of the data, or such representations can be learned from an unsupervised learning step. 451

Atomistic systems can be represented in a myriad of ways. Some descriptions are designed to emphasise particular aspects of a system, while others aim to disambiguate similar chemical or physical principles across a wide range of molecules or materials. The set of desirable properties in a representation thus depends on the task at hand. The following overview gives a coarse characterization of the most popular representations.

where the sum runs over all $N_{A}$ atoms $i$ in structure $A$ and $\mathcal{X}_{i}$ is the environment around atom $i$. When there are multiple chemical species, the descriptors for the local environments of different species can either be included in the single sum, or the averaging can be performed for the environments of each species separately and the species-specific averaged local descriptors can be concatenated. This can be done by considering the Root Mean Square Displacement (RMSD), 452 the best match between the environments of the two structures (best-match), 458 or by combining local descriptors using a regularized entropy match (REMatch). 458

1.2 Representing local environments

We will now describe the Smooth Overlap of Atomic Positions (SOAP) descriptors 412 since many other descriptors based on the atomic density are similar and differ mainly by how the density if projected onto basis functions.454, 419 To construct SOAP descriptors, one first considers an atomic environment $\mathcal{X}$ that contains only one atomic species, and a Gaussian function of width $\sigma$ is then placed on each atom $i$ in $\mathcal{X}$ to make an atomic density function:

to construct the power spectrum of the density using the expansion coefficients:

One then obtains a vector of descriptors $\mathbf{\psi}=\{\psi_{nn^{\prime}l}\}$ by considering all components $l\leq l_{\text{max}}$ and $n,n^{\prime}\leq n_{\text{max}}$ that act as band limits that control the spatial resolution of the atomic density. The generalization to more than one chemical species is straightforward:458 one constructs separate densities for each species $\alpha$ and then computes the power spectra $\psi_{nn^{\prime}l}^{\alpha\alpha^{\prime}}(\mathcal{X})$ for each pair of elements $\alpha$ and $\alpha^{\prime}$ where the two species indices correspond to the $c^{*}$ and $c$ coefficients, respectively. The resulting vectors corresponding to each of the $\alpha$ and $\alpha^{\prime}$ pairs are then concatenated to obtain the descriptor vector of the complete environment.

Atom-centered Symmetry Functions (ACSFs) or sometimes called Behler-Parrinello Symmetry Function411 descriptors differ from SOAP descriptors in that it projects the atomic densities over selected 2-body or 3-body symmetry functions. Faber–Christensen–Huang–Lilienfeld (FCHL) 417 descriptors follow similar principles while also considering the correlations between the atomic densities coming from different chemical species. The Many-Body Tensor Representation (MBTR)418 approach involves taking the histograms of atom counts, inverse pair-wise distances, and angles. Atomic cluster expansion (ACE) descriptors 419 first express atomic densities using spherical harmonics and then generate invariant products by contracting the spherical harmonics with the Clebsch-Gordan coefficients.

Most atomic descriptors use length-scale hyper parameters specifically chosen for a given problem and system.411, 418, 276, 412, 417, 452, 453, 415 There are several ways to automate hyper parameter selections. Ref. 374 introduced general heuristics for choosing the SOAP hyper-parameters for a system with arbitrary chemical composition based on characteristic bond lengths. Ref. 464 adopts the strategy to first generate a comprehensive set of ACSFs and then select a subset using the sparsification methods such as farthest point sampling (FPS) 465 and CUR matrix decomposition.466

A structural descriptor is complete when there is no pair of configurations that produces the same descriptor.467 For atomic descriptors, this means that different atomic environments — after considering the invariances of rotation, translation and permutation of identical atoms — should adopt distinct descriptors. Without completeness, any ML model using the descriptors as input will give identical predictions of physically different systems. Ensuring completeness while preserving the invariances is non-trivial, however. One of the simplest descriptors is based on permutationally-invariant pairwise atomic distances (2-body descriptors), and Ref. 412 demonstrated that these are generally not complete since one can construct two distinct tetrahedra using the same set of distances. Many have assumed that permutationally-invariant 3-body atomic descriptors uniquely specify atomic environments due to the tremendous success of ML models for chemical systems and particularly MLPs. However, Refs. 468 and 467 exemplify that structural degeneracies can be found even when using 3- or 4-body descriptors. This underscores an important shortcoming of state-of-the-art 3-body descriptors such as ACSF, 411 SOAP, 412 FCHL, 417 and MBTR.418 Atomic cluster expansion 419 should be a complete descriptor of local environments, but its reliance on spherical harmonic expansion and the subsequent contraction makes their evaluations expensive. Hence, there are still opportunities to develop improved atomic descriptors.

1.3 The locality approximation

Representing a many-body chemical system in terms atomic environments brings physical significance since certain extensive physical properties (e.g. the total energy, total electrostatic charge, and polarizability of a system) can be approximated by the sum of the atomic contributions coming from each atomic environment, e.g. $\Theta=\sum_{i}\theta(\mathcal{X}_{i})$. This approximation is valid because the atomic contribution associated with a central atom is largely determined by its neighbors, and long-range interactions can be approximated in a mean-field manner without explicitly considering distant atoms. Such “locality” is tacitly assumed in many ML models for CompChem, and it is a crucial necessity for most common atomistic potentials and MLPs (Section 2.2.6). Most MLPs (e.g. BPNN,274 GAP,276 and DeepMD461) approximate the total energy of a system as sums of local atomic energies.

Figure 7 illustrates locality by showing a KPCA map of the atom environments of carbon in the QM9 set (see Section 4.3 for more detailed descriptions of the data set). By color-coding the KPCA plot with the local energies from a SOAP-based GAP model trained on QM9 energies,469 one observes a systematic and smooth trend in energies across clusters. The total molecular energy can then be accurately predicted by the sum of local energies, which means the total energy can be approximated on the basis of all the local environments contained in the molecule. For example, an NN potential trained on liquid water simulations can predict the densities, lattice energies, and vibrational properties of diverse ice phases because the local atomic environments found in liquid water span the similar environments as those observed in ice phases.470 Another GAP potential of carbon trained on amorphous structures and other crystalline phases predicted novel carbon structures in random structure searches as well as approximate reaction barriers.471, 472

The locality approximation is typically rationalized based on the multiscale nature of interatomic interactions in chemical systems. It is generally expected that shorter interatomic distances correspond to stronger interactions, such that a cutoff may be imposed after a certain radial distance $d$ given a certain energy accuracy threshold $\epsilon$. The multiscale nature of interactions underlies the usual classification of chemical interactions, from strong covalent bonds and ionic interactions to weaker non-covalent hydrogen bonds and van der Waals interactions. However, our understanding of non-covalent interactions in large molecules and materials is still emerging 35 and no general rules-of-thumb exist to define the cutoff distance $d$ corresponding to a defined $\epsilon$. Hence, for systems having long-range interactions (which includes most chemical systems), the locality assumption needs revision. Models such as (s-)GDML 372, 207 that learn global interactions can be better suited in these cases. Global models tend to be more data-efficient because they are more specialized toward learning a full PES for a particular molecule or material, but this restricts the use of a single ML model to only the system it was trained upon.

1.4 Advantages of built-in symmetries

Built-in symmetry in ML models substantially compresses the dimensionality of atomic representations and ensures that physically equivalent systems are predicted to have identical properties. One of the most rigorous ways of imposing symmetry onto a model $f$ is via the invariant integration over the relevant group $\mathcal{S}$

where $\mathbf{P}_{\pi}\mathbf{x}$ is a permutation of the input. However, the cardinality of even basic symmetry groups is exceedingly high, which makes this operation prohibitively expensive. This combinatorial challenge can be solved by limiting the invariant integral to the physical point group and fluxional symmetries that actually occur in the training dataset, as done in sGDML.207 Alternative approaches such as parameter sharing 420, 352, 421, 422 or density representations 276 have also proven effective. For example, the DeepMD potential has two versions, the Smooth Edition (DeepPot-SE) explicitly preserves all the natural symmetries of the molecular system, and the other version that does not.461 The DeepPot-SE offers much improved stability and accuracy.207, 461

For ML predictions of scalar properties, the rotationally-invariant atomic descriptor framework described earlier is appropriate. One may wish to predict vectorial or tensorial properties including dipole moments, polarizability, and elasticity. A co-variant version of descriptors may be advantageous, and this can be expressed as

where $S$ indicates a symmetry operation such as a rigid rotation about an axis. Ref. 473 proposed a general method for transforming a standard kernel for fitting scalar properties into a co-variant one. Ref. 474 derived a rotational-symmetry-adapted SOAP kernel, which can be understood as using the angular-dependent SOAP vectors based on spherical harmonics expansions as the descriptors. Note that the SOAP kernels for learning scalar properties introduced in Ref. 412 remove angular dependencies by summing up the SOAP vectors in separate spherical harmonics channels.

Symmetry can be further exploited into “alchemical” representations that incorporate similarity between chemical species that are relatable by changing one atom into another. The FCHL 417 representation considers the similarity between elements in the same row and columns of the periodic table and performs very well on chemical compounds across chemical space. Ref. 475 compiled a data-driven periodic table of the elements by fitting to an elpasolite data set using an alchemical representation.

1.5 End-to-end NN representations

All descriptors introduced above rely on a suitable set of hyperparameters (e.g. length scales, radial and angular resolution). Determining an optimal set of hyperparameters can be a tedious process, especially when heuristics are unavailable or fail due to the structural and compositional complexity of the system. A poor choice of descriptors can limit the accuracy of the final ML model, e.g. when certain interatomic distances can not be resolved.

End-to-end NN representations follow a different strategy to learn a representation directly from reference data. Using atom types and positions of a system as inputs, end-to-end NNs construct a set of atom-wise features $\mathbf{x}_{i}$. These features are then used to predict the property of interest, e.g. the energy as a sum of atom-wise contributions. Unlike static descriptors, the representation is also optimized as part of the overall training process. This way end-to-end NNs can adapt to structural features in the data and the target properties in a fully automatic fashion to eliminate the need for extensive feature engineering from the practitioner.

The deep tensor NN framework (DTNN)352 introduced a procedure to iteratively refine a set of atom-wise features $\{\mathbf{x}_{i}\}$ based on interactions with neighboring atoms. Higher order interactions can then be captured in an hierarchical fashion. For example, a first information pass would only capture radial information, but further interactions would recover angular relations and so on. In DTNN, a learnable representation depending only on atom types $\mathbf{x}^{0}_{i}=\mathbf{e}_{Z_{i}}$ serves as an initial set of features. These are then refined by successive applications of an update function depending on the atomic environment that takes the general form:

Other ML models introduce additional physical information. The hierarchical interacting particle NN (HIP-NN)476 enforces a physically motivated partitioning of the overall energy between the different refinement steps, while the PhysNet architecture477 introduces explicit terms for long range electrostatic and dispersion interactions. In Ref. 420, Gilmer et al. categorize graph networks of this general type as message-passing NNs (MPNNs) and introduce the concept of edge updates. These make it possible to use interatomic information beside the radial distance metric in the refinement procedure, and they have since been adapted for other architectures.478 Another interesting extension are end-to-end NNs incorporating higher-order features beside the scalar $\mathbf{x}_{i}$ used in the original DTNN framework. These are equivariant features that encode rotational symmetry and can be based on angles, dipole moment vectors, or features that can be expressed as spherical harmonics with $l>0$. This enables the exchange of only radial information between atoms in each interaction pass and instead include higher structural information, such as dipole-dipole interactions or angular information. In addition, equivariant end-to-end NNs can also be used to predict vectorial or tensorial properties in a manner similar to the rotational-symmetry-adapted SOAP kernel. Examples include TensorField networks,463 Cormorant,462 DimeNet,479 PiNet,480 and FieldSchNet300.

2 From descriptors to predictions

After a descriptor vector for each chemical structure is defined, one can then construct the design matrix and the kernel matrix for a set of structures. These matrices can then be used as the input of ML models. As described in Section 3, supervised ML methods such as NNs and GPs can be used to approximate non-linear and high-dimensional functions, particularly when massive amounts of training data become available. Thus, one should expect that using CompChem would be very useful for generating a large amount of almost noise-free training data of specific systems or atomic configurations, as long as a physically accurate method is being applied in the right way with appropriate computational resources. In contrast, experimental observations can be difficult to measure and reproduce precisely. Note that the aim of most CompChem+ML efforts have a similar scope as decades-old quantitative structure activity/property relationship (QSAR/QSPR) models that are often based on experiments or CompChem modeling. 481, 327, 326 Thus, researchers in CompChem+ML should be aware of potentially relatable work done by the QSAR/QSPR communities, and to what extent questions being posed have been sufficiently answered. On the other hand, ML usually provides higher accuracy than non-ML statistical models, and so QSAR/QSPR efforts have been turning toward ML models as well.482

Depending on the types of chemical problems being solved, different CompChem methods suitable for treating different chemical interactions will be used for training ML models having different approximations for underlying high-dimensional functions. For example, research efforts are towards learning electron densities,483 density functionals,162 and molecular polarizabilities.484

Besides these direct learning strategies, ML has been used to enhance the performance and suitability of CompChem models. As mentioned in Section 2, the $\Delta$-ML 485 approach is now a common technique for adapting an ML model that improves the quality of a theoretically insufficient but computationally affordable method. This approach has been used to learn many body corrections for water molecules to allow a relatively inexpensive KS-DFT approach like BLYP to more accurately reproduce CCSD(T) data.486 Along similar lines, Shaw and coworkers used CompChem features along with a neural network to reweight terms from an MP2 interaction energy to provide ML-enhanced methods with increased performance.126 Miller and coworkers have developed ML-models where molecular orbitals themselves are learned to generate a density matrix functional that provides CCSD(T)-quality potential energy surfaces with a single reference calculation.487 von Lilienfeld and coworkers have investigated how the choice of regressors and molecular representations for ML models impacts accuracy, and their finding suggest ways that ML models may be trained to be more accurate and less computationally expensive than hybrid DFT methods.488 Burke and co-workers have studied how ML methods can result in improved understanding and more physical exact KS-DFT181, 489, 490, 491 and OFDFT functionals.161 Brockherde et al. have presented an approach where ML models can directly learn the Hohenberg-Kohn map from the one-body potential efficiently to find the functional and its derivative.162, 184 Akashi and coworkers have also reported the out-of-training transferability of neural networks that capture total energies, which shows a path forward to generalizable methods.492

Besides the above-mentioned applications of supervised learning, one can exploit the “universal approximator” nature of ML architectures to find a function that gives the best solution in a variational setting. For instance, using restricted Boltzmann machines 493 or deep neural networks as a basis representation of wavefunctions 494, 105, 106 in Quantum Monte Carlo calculations.

3 CompChem data

We have laid the general framework for CompChem+ML studies, but this direction would not be complete without more details about training data (i.e., garbage in, garbage out). We now review the landscape of data sets in CompChem and how they will likely evolve over time. The past decade has seen continually increasing usefulness and availabilty of “big data” from CompChem that include community-wide data repositories comprised of millions of atomistic structures along with diverse physical and chemical properties.495, 496, 497, 498 Such repositories are becoming the norm, and it is more customary for different users to deposit raw or processed simulation data there for the benefit of the research community. This brings the possibility of robust validation tests for ML models, but it also necessitates approaches that are well-equipped to handle large and complex data sets. Typical data sets may come from diverse origins such as MD trajectories from ab initio simulations, data sets of small molecules and molecular conformers, or other training sets used for developing ML and non-ML forcefields for specific applications. As the data sets grow, so do the scope of publications that involve ML as shown in Fig. 1.

ML model must be validated before they can be trusted for predictions. Validations of descriptors or model trainings are performed on benchmark data sets, and the most popular ones are summarized in Table 5. These allow ML models to be compared on the same ground and provide large amounts of data for robust training. Their availability to the public also ensures that the data sets can evolve with time and be extended as a part of community efforts.517

Amongst the entries in Table 5, the most often used one is the QM9 set, which consists of approximately 134,000 of the smallest organic molecules that contain up to 9 heavy atoms (C, O, N, or F; excluding H) along with their CompChem-computed molecular properties such as total energies, dipole moments, HOMO-LUMO gaps, etc. Several ML studies have already been published using this dataset (see Fig. 8, Ref. 488). A popular challenge is to develop a next-generation ML model that learns the electronic energies of random assortments of organic molecules with higher accuracy and less required training data than other existing models. Doing so tests next generation molecular representations and training algorithms. The next significant advance will potentially be due to a combination of supervised and unsupervised learning models.

3.2 Visualization of data sets

As the structural data sets grow it becomes infeasible to manually identify hidden patterns or curate the data. Data-driven and automated frameworks for visualizing these data sets become increasingly popular.518, 519, 520 Dimensionality reduction effectively translates the high dimensional data (i.e. the xyz-coordinates for molecules or materials in different atomic configurations) into a low-dimensional space easily visualized on paper or a computer screen. In this way, entries such as those in the QM9 set can be shown (see Fig. 9). This map thus is a useful tool to help navigate the QM9 set by showing how molecules of different compositions (represented as colored dots) are similar or dissimilar based on molecular properties (i.e. atomization energies) along the principal axes. Similarly, Ref. 321 used SOAP-sketchmaps in conjunction with quasi-chemical theory to show an unsupervised learning procedure for identifying local solvent environments that significantly impact solvation energies of small ions.

These data-driven maps are generated by processing the design matrix (or kernel matrix) associated with a data set using dimensionality reduction techniques introduced in Section 3.2. A simple option is to use the ASAP code, 374 a Python-based command line tool, that automates analysis and mapping. Fig. 7 and 9 were generated using ASAP using only two commands that are displayed in the figure. Data sets can also be explored in an intuitive manner using interactive visualizers 521 that run in a web browser and display 3D-structures corresponding to each atomistic structure in the data set.

3.3 Text and data mining for chemistry

Conventional publications are an essential part of CompChem knowledge base, and ML is becoming useful at accelerating information extraction from the scientific literature via text mining.522, 523, 524 This topic was previously comprehensively reviewed in the context of cheminformatics.525, 526 Natural language processing has already driven text-mining efforts materials science discovery525 and experimental synthesis conditions of oxides.516, 527 CompChem+ML can also amplify existing efforts in chemometrics,528 the science of data-driven extraction of chemical information.529 This area has also branched into related disciplines of data mining for specific classes of materials530 and catalysis informatics.531 These approaches have great promise, especially for deriving information and knowledge from data, but it remains challenging to implement these in ways that achieve insight (and true impact).

Some have shown paths forward for doing so. For example, ML models can obtain knowledge from failed experimental data more reliably than humans who are more susceptible to survivor bias,532 and it can also be used to distill physical laws and fundamental equations using experimental363 and computational data.533 ML models can also be used to reliably predict SMILES representations that allow encoded information to be derived from low-resolution images found in the literature.534 ML models can interpret experimental X-ray absorption near edge structure (XANES) data and predict real space information about coordination environments.535 Likewise, scanning tunneling microscopy (STM) data can be used to classify structural and rotational states on surfaces,536 and name indicators can be used to predict in tandem mass spectrometry (MS/MS) properties.537 In closing, we see exciting opportunities for future applications that complement data and text mining to chemometrics through chemical space.

4 Transforming atomistic modelling

Training an MLP to reproduce a system’s PES usually requires generating diverse and high quality CompChem data points that cover the relevant temperature and pressure conditions, reaction pathways, polymorphs, defects, compositions, etc.542, 543, 544, 545, 546, 547, 548, 549 After data points comprised of atomic configurations, system energies, and forces are obtained, different methods for constructing MLPs employ either different descriptors (see a list of examples in Table 4) or different ML architectures to perform interpolations of the full PES. Again, smoothness is an essential feature for any PES, so special considerations are needed to avoid numerical noise that would result in discontinuities. 550, 551 Kernel method-based MLPs such as GAP276, 552 and sGDML372, 207, 553 ensure smoothness by relying on smoothly varying basis functions, but the scaling of kernel-based methods with respect to the number of training points is challenged without reduction mechanisms. 396, 554 As a much more efficient but somewhat less accurate alternative to GAP, Spectral Neighbor Analysis Potential (SNAP) 555 uses the coefficients of the SOAP descriptors and assumes a linear or quadratic relation between energies and the SOAP bispectrum components.556 The most popular MLPs are currently NN-based due to their flexibility and capacity to train based on large amounts data. Amongst these, ANI 500, 502 and BPNN432, 274, 557 potentials use ACSF descriptors as inputs, while Deep neural networks such as SchNet421, 433, 558 and DeepMD559 use the coordinates and nuclear charges of atoms. We now focus on a few example applications.

Many CompChem efforts focus on predicting thermodynamic properties at finite temperatures, such as heat capacity, density, and chemical potential. Although many physical properties are already accessible from MD simulations, doing estimations of free energies that establish the relative stability of different states using electronic structure methods remains difficult. The configurational part of the Gibbs free energy of a bulk system that has $N$ distinguishable particles with atomic coordinates $\textbf{r}=\{\textbf{r}_{1\ldots N}\}$, and the associated potential energy $U(\textbf{r})$ can be expressed as:

integrated over all possible coordinates r, where $k_{B}$ is the Boltzmann constant. In order to rigorously determine $G$, one must exhaustively sample the configuration space that has relatively high weight $\exp\left[-\dfrac{U(\textbf{r})+PV}{k_{B}T}\right]$. This normally requires thermodynamic integration or enhanced sampling methods (e.g. umbrella sampling, 560 metadynamics, 561 transition path sampling562), that require simulation times and scales far beyond what is accessible with MD simulations based on KS-DFT or correlated wavefunction methods.

However, MLPs have unleashed both limits on the time scale and system size. An early example, 563 used an MLP with umbrella sampling 560 and the free energy perturbation method 564 to reveal the influence of van der Waals corrections on the thermodynamic properties of liquid water. Later, the combination of an MLP trained from hybrid DFT data and free energy methods reproduced several thermodynamic properties of water from quantum mechanics, including the density of ice and water, the difference in melting temperature for normal and heavy water, and the stability of different forms of ice.565, 566 Ref. 567 employed the DeepMD approach to study the relatively long time-scale nucleation of gallium. MLPs for high-pressure hydrogen provided evidence on how hydrogen gradually turns into a metal in giant planets. 568 In all these examples, high accuracy and long timescales were required to model the specific phenomena and reveal physical insights, and it is precisely the combination of CompChem+ML that enables both.

4.2 Nuclear quantum effects

As mentioned in Sec. 2.2.5, NQEs of chemical systems having light elements bring challenges for atomistic modelling because the added mobility of lighter atoms in dynamics simulations requires higher computational cost to treat. To make the matter even more complicated, many atomistic potentials (see Sec. 2.2.6), particularly the ones for water or organic molecules, cannot be used to model NQEs, because they often describe colavent bonds as rigid, and thus cannot describe the fluctuations of the bond lengths and angles. As a remedy, several studies have been performed by training an MLP using higher rungs of KS-DFT (e.g. hybrid-DFT or meta-GGA) and this use this potential in PIMD simulations.569, 570, 565, 571 The study of water that was mentioned in the previous section and used MLPs trained from hybrid DFT revealed that NQEs were critical for promoting the hexagonal packing of molecules inside ice that ultimately lead to the six-fold symmetry of snowflakes.565 Highly data efficient ML potentials can even trained on reference data at the computationally very expensive quantum-chemical CCSD(T) level of accuracy. For example, the sGDML207, 572, 206 approach has been shown to faithfully reproduce such force fields for small molecules, which were then used to perform simulations with effectively fully quantized electrons and nuclei.

5 ML for structure search, sampling, and generation

Locating stationary points on the potential-energy surface (PES) is a frequent task in CompChem, since minima represent thermodynamically stable states and first order saddle points for transition states that define barrier heights that explaining reaction kinetics. Explorations for stationary points normally require many energy and force evaluations. ML approaches are being implemented to dramatically accelerate minimum energy as well as saddle-point optimizations.294, 295, 293, 573, 574, 575, 552 Bernstein et al. proposed an automated protocol that iteratively explores structural space using a GAP potential.552 Bisbo and Hammer employed an actively-learned surrogate model of the potential energy surface to performs local relaxations while only performing single-point quantum-mechanical calculations for selected structures with high values of acquisition.573 Refs. 293, 295, 296, 297 accelerate nudged elastic band (NEB) calculations by incorporating a surrogate ML models.

ML can also dramatically accelerate the challenge of efficiently sampling equilibrium or transition states by accelerating enhanced sampling methods such as umbrella sampling560 and metadynamics.561 These procedures make use of collective variables (CVs) that define a reaction coordinate, and computing the associated free energy surface (FES) amounts to generating the marginal probability distribution in these CVs. Unfortunately, the choice of the CVs is not always clear for specific systems, and ML have shown some promise in guiding their determination.576, 577, 578 Another direction is to exploit that ML models can be considered as universal approximators of free energy surfaces.579 For example, there are reports of adaptive enhanced sampling methods using a Gaussian Mixture model,580 using an NN architecture to represent the FES581 or the bias function in variational sampling simulations. 582

ML methods also offer fundamentally new ways to explore chemical compound and configuration space. Generative models can learn the structural and elemental distribution underlying chemical systems, and once trained, these models can then be used to directly sample from this distribution. It is furthermore possible to bias the generated structures towards exhibiting desired properties, e.g. drug activity or thermal conductivity. As a consequence, generative models offer exciting new avenues in drug and materials design.583, 584 Generative methods in CompChem include recurrent neural networks (RNNs), which can be used for the sequential generation of molecules encoded as SMILES strings (a string based representation of molecular graphs).585, 586, 587 Segler et al. demonstrated how such a recurrent model can first learn general molecular motifs and then be fine-tuned to sample molecules exhibiting activity against a variety of medical targets.586 Autoencoders (AE) are another frequently used ML method for molecular generation. AEs learn to transform of molecular graphs or SMILES into a low-dimensional feature space and backwards. The resulting feature vector represents a smooth encoding of the molecular distribution and can be used to effectively sample chemical space.588, 589, 590, 591, 592, 593 By applying a variational AE to the QM9 and ZINC databases, Gomez-Bombarelli et al. could generate several optimized functional compounds.594 An interesting extension to AEs are conditional AEs, which not only capture the distribution of molecular structures but also its dependence on various properties.425, 595 This makes it possible to directly generate structures exhibiting certain property ranges or combinations without the need for biasing or additional optimization steps. AEs can also form the basis of another approach for exploring chemical space called generative adversarial networks (GANs).596, 597 In a GAN, a generator model (often an AE) attempts to create samples that closely match the underlying data, while a discriminator tries to distinguish true from generated samples. These architectures can be enhanced by using reinforcement learning (RL) objectives. RL learns an optimal sequence of actions (e.g. placement of atoms) leading to a desired outcome (e.g. molecule with certain property). This makes it possible to drive generative process towards certain objectives, allowing for the targeted generation of molecules with particular properties.598, 599, 600, 601 RL in general is a promising alternative strategy for generative models,602, 603 and they offer the possibility for tightly integrating them into drug design cycles.604 Alternative approaches combine autoregressive models with graph convolution networks.605, 606

While these methods use SMILES or graphs to encode molecular structures, generative models have recently been extended to operate on 3D coordinates of molecules and materials.607, 608 Gebauer et. al proposed a auto-regressive generative model based on the SchNet architecture, called g-SchNet.609 Once trained on the QM9 dataset, g-SchNet was able to generate equilibrium structures without the need for optimization procedures. It was further found, that the model could be biased towards certain properties. In another promising approach, Noé et al. used an invertible neural network based on normalizing flows to learn the distribution of atomic positions (e.g. sampled from a molecular dynamics trajectory). This network can then be used to directly sample molecular configurations by sampling from this distribution without performing costly simulations.298

6 Multiscale modeling

Multiscale modeling is a term for including simulation or information from different scales (see Fig. 3). ML has been introduced into QM/MM-like schemes that enable improved multiscale simulations,610, 325, 300 and on the side of coarse-graining.611 Coarse-graining have been developed, 612 but the inherent functional form for these potentials relies on CPI as well as trial-and-error procedures. Several works used ML for constructing coarse-grained potentials by matching mean forces.613, 614, 449, 450 In closing, we see promise for experimental priors into ML models, for instance, using experimental measurements to improving a ML potential energy surface by complementing with experimental data. We are not aware of such efforts for developing highly accurate MLPs beyond the atomic scale, although much work has been done along this line to refine force fields of RNAs and proteins, often incorporating methods from ML, including the maximum entropy approach.615

Selected applications and paths toward insights

The central challenge posed at the beginning of this review was how to identify and make chemical compounds or materials having optimal properties for a given purpose. To do so would help address critical and broad issues from pollution, global warming to human diseases. Traditional developments are often slow, expensive, and restricted by non-transferable empirical optimizations, and so efforts have turned to CompChem+ML to alleviate this.616, 617, 504, 513

CompChem+ML are enabling searches through larger areas of chemical space much faster than before.618, 619, 620, 621, 19 This section is not to extensively review the large amount of work using CompChem+ML in these different areas, but rather to highlight examples of applications that have resulted in notable insights so that others might use these notable works as templates for future efforts.

Molecules and materials design is usually considered to be an optimization problem.270, 594, 589, 622, 425 Thus, a comprehensive understanding of chemical space is needed to identify compounds with desired properties that are subject to certain required constraints (e.g. a specific thermal stability or a suitable optical gap for absorbing sunlight). Those properties will also depend on many key variables (e.g. constitutive elements, crystal forms, geometrical and electronic characteristics, among others), which make the property prediction complex.519 CompChem calculations as explained in Section 2 should provide a continuous description of properties across a continuous representation (i.e. a descriptor or fingerprint) of molecules that is used to map molecular configurations to target properties, and vice-versa. ML methods then can be implemented to search large databases to extract structure-property relationships for designing compounds with specific characteristics.519, 623, 622, 624 Optimizations would then be performed on the structure-based function learned from training configurations, and the composition of the chemical compound would then be recovered back from the continuous representation.

As a protoypical example of molecular design via high-throughput screening, Gomez-Bombarelli et al.619 showed a computation-driven search for novel thermally activated delayed fluorescence organic light-emitting diode (OLED) emitters. That work first filtered a search space of 1.6 million molecules down to approximately 400,000 candidates using ML to anticipate criteria for desirable OLEDs. For the purpose of evaluating candidates, they estimated an upper bound on the delayed fluorescence rate constant (kTADF). Time-dependent DFT calculations were then used to provide refined predictions of specific properties of thousands of promising novel OLED molecules across the visible spectrum so that synthetic chemists, device scientists, and industry partners would be able to choose the most promising molecules for experimental validation and implementation. Notably, this example of CompChem+ML resulted in new devices that exhibited an external quantum efficiency of over 22$\%$. Fig. 10 shows the high accuracy of ML in predicting useful properties for high-throughput screening of molecules and materials based on kTADF calculations. This work exemplifies how ML can accelerate the design of novel compounds in such a way that could not be possible using traditional CompChem methods alone.

One can also exploit the flexibility of ML models to increase the quality of the chemical insights. Integrations of features relevant to learning tasks allow one to improve the accuracy of ML predictions for a given target property. Park and Wolverton625 improved the performance of the crystal graph convolution neural network (CGCNN)626 by adding to the original framework information about Voronoi tessellated crystal structure, which are explicit 3-body correlations of neighboring constituent atoms, and an optimized representation of interatomic bonds. The new approach that was labeled as iCGCNN achieved a predictive accuracy 20$\%$ higher than that of the original CGCNN when determining thermodynamic stabilities of compounds (i.e. predictions of hull distances). When used for high-throughput searches, iCGCNN exhibited a success rate higher than an undirected high-throughput search and higher than that of CGCNN. Fig. 11 shows the improvement in predictions of nearly stable compounds after using more appropriate descriptors. This study showcases how descriptors can be tailored to further enhance the success of ML-aided high-throughput screening.

2 Retrosynthetic technologies

A grand challenge in chemistry is to understand synthetic pathways to desired molecules.627, 628 Retrosynthesis involves the design of chemical steps to produce molecules and materials that would be crucial to drug discovery, medicinal chemistry, and materials science. As a different kind of optimization problem, the general tactic is to analyze atomic scale compounds recursively, map them onto synthetically achievable building blocks, and then assemble those blocks into the desired compound.629, 630, 631

Three main issues make retrosynthesis a formidable intellectual challenge:632 First, simple combinatorics make the space of possible reactions greater than the space of possible molecules. Second, reactants seldom contain only one reactive functional group, and thus require predictions of multiple functional groups. Third, one failed step in the route can invalidate the entire synthesis because organic synthesis is a multistep process.

Given these challenges, ML is becoming more established in determining reaction rules from CompChem data.628 Computer-aided synthesis planning was actually first attempted in the 1960s.633 Many have since attempted to formalize chemical perception and synthetic thinking using computer programs.634, 635, 636 These programs are typically based on one of three possible algorithms:636

Algorithms that use reaction rules (manually encoded or automatically derived from databases).

Algorithms that use principles of physical chemistry based on ab initio calculations to predict energy barriers.

ML approaches are used to try to overcome the generalization issues of rule-based algorithms (that normally suffering from incompleteness, infeasible suggestions, and human bias) while also avoiding the high cost of CompChem calculations. It is now possible to obtain purely data-driven approaches for synthesis planning, which are promoting a rapid advancement in the field. For example, Coley and co-workers637 designed a data-driven metric, SCScore, for describing a real synthesis modeled after the idea that products are, on average, more synthetically complex than each of their reactants. The definition of a metric for selecting the most promising disconnections that produce easily synthesizable compounds is crucial for avoiding combinatorial explosion. Fig. 12 shows that a data-driven metric, the SCScore, is more suitable than other heuristic metrics to perceive the complexity of each step in a given synthesis. This work offered a valuable contribution to the retrosynthesis working pipeline by providing a method that implicitly learns what structures and motifs are more prevalent as reactants.

Apart from isolated approaches or algorithms to deal with specific tasks within retrosynthesis, there is already software available to advance this field. One example is the Chematica program,638 which has implemented a new module that combines network theory, modern high-power computing, artificial intelligence, and expert chemical knowledge to design synthetic pathways. A scoring function is used to promote synthetic brevity and penalize any reactivity conflicts or non-selectivities, thus allowing it to find solutions that might be hard for a human to identify. Fig. 13A shows the decision tree for one of the almost 50,000 reaction rules used in Chematica. Reaction rules can be considered as the allowed moves from which the synthetic pathways are built, and such moves lead to an enormous synthetic space (the number of possibilities within n steps scales as 100n) as the one shown by the graph in Fig. 13B. Chematica explores this large synthetic space by truncating and reverting from unpromising connections and drives its searches to the most efficient sequences of steps. Moreover, in the pathways presented to the user, each substance can be further analyzed with molecular mechanics tools. This software was used to obtain insights into the synthetic pathways to eight targets (seven bioactive substances and one natural product). All of the computer-planned routes were not only successfully carried out in the laboratory, but they also resulted in improved yields and cost savings over previous known paths. This work opened an avenue for chemists to finally obtain reliable pathways from in silico retrosynthesis. For further reading we recommend the two-part reviews of Coley and co-workers.639, 640

3 Catalysis

Catalysis research requires multiscale approaches to determine chemical compounds that can influence barrier heights of reaction mechanisms to impact product yields and selectivities without otherwise being generated or consumed by the reaction.641 Traditional catalysis is normally discussed in textbooks in terms of homogeneous (i.e. within a solution phase), heterogeneous (occurring at a solid/liquid interface), and biological (occurring within enzymes and riboenzymes), but it is best not to use these terms too strictly because actual reaction mechanisms can be quite complex and overall processes may sometimes exhibit characteristics (by design) of two or more of these classical processes.642, 643, 644 Modern research in catalysis has been interested in studying chemical reactivity and reaction selectivity arising from stimuli from solar thermal energy,645, 646 electrochemical potentials,647 photons,648, 649, 650, 651 plasmas,652, 653 or other external resonances.654 Catalysis makes up roughly 35% of the world’s gross domestic product,655 and it is important to guide toward the end goal of achieving greater sustainability with catalytic processes.656, 657, 658

These reasons help make catalysis a fertile training ground for applying and developing theoretical models (e.g. Refs. 659, 660, 661) that can be used along with CompChem or CompChem+ML. The research field is also burgeoning with many reports and review articles662, 663, 664, 531, 665, 666 that discuss perspectives and progress using ML methods for catalysis science; here, we will mention notable examples that present a broad range of ways that CompChem+ML can be used for insights. For example, CompChem+ML methods are enabling more data generation by allowing costly QC calculations to be run more efficiently, and more information means more comprehensive predictions of chemical and materials phase diagrams for catalysis667, 668 as well as stability and reactivity descriptors identified on the fly.669, 670, 671, 672, 673 Figure 14 shows examples of the palettes of insight available using state-of-the-art CompChem+ML modeling for identifying activity and selectivity maps, as well as visualizations of data using t-Distributed Stochastic Neighbor Embedding (t-SNE).674

Regarding modeling of deeply complex chemical environments, Artrith and Kolpak developed MLPs for investigating the relationships between solvent, surface composition and morphology, surface electronic structure, and catalytic activity in systems composed of thousands of atoms interfaces.676 We expect such simulations for electro- and photo-catalysis elucidation will continue to improve in size, scale, and accuracy. For other physical insights, new approaches by Kulik and Getman and co-workers have also focused on developing ML models appropriate for elucidating complex $d$-orbital participation in homogeneous catalysis.677 Rappe and co-workers have used regularized random forests to analyze how local chemical pressure effects adsorbate states on surface sites for the hydrogen evolution reaction.678 Almost trivially simple ML approaches can be used in catalysis studies to deduce insights into interaction trends between single metal atoms and oxide supports,679 to identify the significance of features (e.g. adsorbate type or coverage) where CompChem theories break down,680 or they can be used to identify trends that result in optimal catalysis across multiple objectives such as activity and cost (Fig. 15).681

ML is also opening opportunities for CompChem+ML studies on highly detailed and complex networks of reactions.682, 683, 684, 685, 686, 687 Such models in principle can then significantly extend the range of utility of microkinetics modeling for predictions of products from catalysis.688, 689 ML also enables studies of complicated reaction networks that can allow predictions of regioselective products based on CompChem data,690 asymmetric catalysis important for natural product synthesis,691, 692 and biochemical reactions.693 Efforts to better understand ‘above-the-arrow’ optimizations of reaction conditions relate back to the challenge of retrosynthetic challenges.694, 695 Ideally, these efforts will continue while making use of rapid advances in CompChem+ML that enable predictive atomistic simulations to be run faster and more accurately. We see reason for excitement for different approaches, but we again stress the importance of ensuring that models will provide unique and physical results (see Section 3 where we discuss the risk of “clever Hans” predictors360).

4 Drug design

The central objective for drug discovery is to find structurally novel molecules with precise selectivity for a medicinal function. This involves identifying new chemical entities and obtaining structures with different physicochemical and polypharmacological properties (i.e., combinations of beneficial pharmacological effects and/or adverse side-effects).696, 697 Drug discovery involves the identification of targets (a property optimization task, as in material design) and the determination of compounds with good on-target effects and minimal off-target effects.698 Traditionally, a drug discovery program may take around six years before a drug candidate can be used in clinical trials and additional six or seven years are required for three clinical phases. Thus, it is important to identify adverse effects as soon as possible to minimize time and monetary costs.699 Accelerating drug discovery relies on predicting how and where a certain drug binds to more than one protein, a phenomenon that sometimes results in polypharmacology. Researchers are developing ready-to-use tools aimed to facilitate research for drug discovery,700 but CompChem+ML is expected to continue providing even more benefits to the drug development pipeline.701

In a recent study, Zhavoronkov et al.604 developed a deep generative model for de novo small-molecule design: the generative tensorial reinforcement learning (GENTRL) model that was used to discover potent inhibitors of discoidin domain receptor 1 (DDR1), a kinase target implicated in fibrosis and other diseases. The drug discovery process was carried out in only 46 days, beginning with the recollection of appropriate data for training and finishing with the synthesis and experimental test of some compounds (Fig. 16A). GENTRL was used to screen a total of 30,000 structures (some examples compared to the parent DDR1 kinase inhibitor are shown in Fig. 16B) down to only 40 structures that were randomly selected ensuring a coverage of the resulting chemical space and distribution of root-mean squared deviation values. Six of these molecules were then selected for experimental validation (see Fig. 16C), with one of them demonstrating favorable pharmacokinetics in mice. The predicted conformation of the successful compound according to pharmacophore modelling was very similar to the one predicted to be preferred and stable by CompChem methods. This work illustrates the utility of CompChem+ML approaches to give insights into drug design by rapidly giving compound candidates that are synthetically feasible and active against a desired target.

Besides generating new chemical structures with favorable pharmacokinetics, ML methods are also used in pharmaceutical research and development for peptide design, compound activity prediction and for assisting scoring protein-ligand interaction (docking).702, 703, 696, 704 An example of the latter was proposed by Batra et al.705 for efficiently identifying ligands that can potentially limit the host-virus interactions of SARS-CoV-2. Those authors designed a high-throughput strategy based on CompChem+ML that involved high-fidelity docking studies to find candidates displaying high-binding affinities. The ML model was used to search through thousands of approved ligands by the Food and Drug Administration (FDA) and a million biomolecules in the BindingDB database.503 From these, insights were obtained for more than 19,000 molecules satisfying the Vina score (i.e. an important physicochemical measure of the therapeutic process of a molecule that is used to rank molecular conformations and predict free energy of binding). Fig. 17 shows the Vina score predictions that led to the selection of the best candidates, some of which are also illustrated in the figure. The Vina scores for the top ligands were further confirmed using expensive docking approaches, resulting in the identification of 75 FDA-approved and 100 other ligands potentially useful to treat SARS-CoV-2. This study highlights a reasonable CompChem+ML strategy for making useful suggestions to aid expert biologists and medical professionals to focus in fewer candidates when performing either robust CompChem efforts or synthesis and trial experiments.

Conclusions and Outlook

Recent CompChem methods, algorithms, and codes have empowered new studies for a wealth of physical and chemical insights into molecules and materials. Today, the combination of CompChem+ML can be equipped to address new and more challenging questions in different domains of physics, materials science, chemistry, biology, and medicine. Productive research efforts in this direction necessitate interdisciplinary teams and increasing availability of high-quality data across appropriate regions of chemical compound space. Discovering new chemicals and materials requires thorough investigations. One needs to predict reaction pathways and interactions between molecules, optimize environmental conditions for catalytic reactions, enhance selectivities that eliminate undesired side reactions and/or side effects, and navigate other system-specific degrees of freedom. Addressing this complexity calls for a statistical view on chemical design and discovery, and CompChem+ML provides a natural synergy for obtaining predictive insights to lead to wisdom and impact.

This review provided an essential background of CompChem and ML and how they can be used together to make transformative impacts in the chemical sciences by amplifying insights available from CompChem methods. The successes of CompChem+ML are particularly visible in physical chemistry and include drastic acceleration of molecular and materials modeling, discovery and prediction of chemicals with desired properties, prediction of reaction pathways, and design of new catalysts and drug candidates. Nevertheless, we have only begun to scratch the surface of how successful applications of ML in chemistry can bring impact. There are many conceptual, theoretical, and practical challenges waiting to be solved to enable further synergies within the troika of CompChem, ML, and CPI. Here we enumerate some of the challenges that we consider to be the most pressing and interesting at this moment:

Reliance on ML in CompChem algorithms must be increased: ML algorithms can be integrated into CompChem algorithms at almost any simulation level (Fig. 3). ML algorithms are already available to accelerate calculations of CompChem energies, navigations along reaction pathways, and sampling of larger regions of the PES, but the reluctance of their use impedes progress. In general, these algorithms must be made more effective, efficient, accessible, user-friendly, and reproducible to benefit fundamental and applied research (see for example, Ref. 706).

More general ML approaches are needed: ML methods must continue to evolve beyond now-common applications of learning a narrow region of a PES or identifying straightforward structure/property relationships. New ML methods should have the capacity to predict energetic and electronic properties and their more convoluted relationships across chemical space. Such approaches should grow toward uniformly describing compositional (chemical arrangement of atoms in a molecule) and configurational (physical arrangement of atoms in space) degrees of freedom on equal footing. Further progress in this field requires developing new universal ML models suitable for insights across diverse systems and physicochemical properties.

ML representations must to include the right physics: ML methods that are claimed to be accurate but incorrectly describe the true physics of a system will eventually fail to achieve meaningful insights while lowering the reputation of other work in the field. Current ML representations (descriptors) can successfully describe local chemical bonding, but few if any are treating long-range electrostatics, polarization, and van der Waals dispersion interactions that are critical for rationalizing physical systems, both large and small. Combining intermolecular interaction theory (a key focus of advanced CompChem methods) with ML is an important direction for future progress towards studying complex molecular systems.

CompChem + ML applications need to strive toward achieving realistic complexity: Investigations using highly accurate CompChem methods normally require overly simplified model systems while more realistic model systems necessitate less accurate but computationally efficient CompChem methods. This compromise should no longer be necessary. We are due for a paradigm shift in how thermodynamics, kinetics, and dynamics of systems in complex chemical environments (e.g. for multiscale biological processes like drug design and/or catalytic processes at solid liquid interfaces under photochemical excitations, etc.) can be treated more faithfully with less corner-cutting. An emerging idea is to dispatch ML approaches into computationally efficient model Hamiltonians for electronic interactions based on correlated wavefunction, KS-DFT, tight-binding, molecular orbital techniques, and/or the many-body dispersion method. ML can predict Hamiltonian parameters and the quantum-mechanical observables would be calculated via diagonalization of the corresponding Hamiltonian. The challenge is to find an appropriate balance between prediction accuracy and computational efficiency to dramatically enhance larger scale simulations.

(Much) more experimental data is needed: Validations of ML predictions require extensive comparisons with experimental observables such as reaction rates, spectroscopic observations, solvation energies, melting temperatures. Such experiments may have previously been considered too routine, too mundane, or not insightful enough alone, but all high quality brings great value for future CompChem+ML efforts that tightly integrate quantum mechanics, statistical simulations, and fast ML predictions, all within a comprehensive molecular simulation framework. 707

(Much) more comprehensive data sets need to be assembled and curated: Current CompChem+ML efforts have profited heavily by the availability of benchmark data sets for relatively small molecules that allow a comparison of existing models.413, 515 While efforts fixated on boosting prediction accuracies and shrinking down requisite training set sizes for ML models have had their merits, it is time to move on as further improvements are meaningless if the ML models are not making useful and insightful predictions themselves. More useful predictions will require knowledge from larger datasets, and these will inevitably contain heterogeneous combinations of different levels of theory and/or experiments that must be analyzed, ‘cleaned’, and uncertainties adequately quantified in order for models to productively learn. Such hybrid data sets may be the key to arrive at novel hypotheses in chemistry that could then be experimentally tested.

Bolder and deeper explorations of chemical space are needed: So far most efforts to generate chemical data have focused on exploring parts of chemical space for new compounds for a targeted purpose. This should change. Combining ML model uncertainty estimates across broader swaths of chemical space could open pathways for fruitful statistical explorations, say, in an active learning framework. This could lead to discovering new synergies between data that otherwise would not have been possible to enable advances in scientific understanding and improve ML models. Generative models can bridge the gap between sampling and targeted structure generation imposing optimal compound properties, e.g. for inverse chemical design.608, 609, 125

This and other reviews708, 621, 709, 19, 707, 541, 710, 711 have stated how ML has become instrumental for recent progress in CompChem. We would like to also mention inspirations that ML has drawn from being applied to physical and chemical problems.

ML methods generally assume that data is subject to measurement noise while CompChem data is generally approximate but also noise-free from a statistical perspective. ML modeling still requires regularization, but regularizers should here particularly reflect the underlying physics of molecular and materials systems. ML models used in applications of vision contain discrete convolution filters that are suboptimal for chemical modeling, but recognition of this shortcoming has led to novel continuous convolution filters that are well suited for chemistry and have also become a popular novel architecture for core ML methods.433

Furthermore, invariances, symmetries, and conservation laws are key ingredients to physical and chemical systems. Incorporating them into ML has led to novel and useful models for chemistry since they can learn from significantly less data, which then makes it possible to build force fields at unprecedentedly high levels of theory.372, 207, 206 Using these powerful ML techniques for computer vision, natural language processing, and other applications is currently being explored. Structural information from molecular graphs provide the basis for novel tensor neural networks or message passing architectures352, 420 as well as graph explanation methods.712

Many further challenges exist that have led or will lead to mutual bidirectional cross-fertilization between ML and chemistry. These interdisciplinary efforts also initiate progress in respective application domains. The power of this path is that solving a burning problem in chemistry with a novel crafted ML model may also result in unforeseen insights in how to better design core ML methods. Interestingly, the exploratory usage of ML for knowledge discovery in chemistry typically requires novel ML models and unforeseen scientific innovations, and this can lead to interesting insight that is not necessary limited to chemistry alone, rather it is likely to go beyond.

To conclude, the past decade has shown that it has not been enough to just apply existing ML algorithms, but breakthroughs are happening by a handshaking of innovations resulting in novel ML algorithms and architectures driven by the pursuit of novel insights in chemistry while retaining a deep understanding about the underlying physical and chemical principles. Research programs that foster interdisciplinary exchange such as IPAM (www.ipam.ucla.edu) have seeded this progress, and these should be continued. Mixed teams with members educated in different aspects of physics, chemistry and ML have been instrumental. This also brings the need to solve the new educational challenge of developing new generations of researchers with an academic curriculum that interweaves chemistry, physics and computer science to enable a meaningful (multilingual) research contribution to this exciting emerging field.

JAK was supported by the Luxembourg National Research Fund (INTER/MOBILITY/19/13511646) and the U.S. National Science Foundation (CBET-1653392 and CBET-1705592).

VVG acknowledges financial support from the Luxembourg National Research Fund (FNR) under the program DTU PRIDE MASSENA (PRIDE/15/10935404)

BC acknowledges funding from the Swiss National Science Foundation (Project P2ELP2-184408).

KRM was supported in part by Institute of Information & Communications Technology Planning & Evaluation (IITP) grants funded by the Korea Government (No. 2017-0-00451, Development of BCI based Brain and Cognitive Computing Technology for Recognizing User’s Intentions using Deep Learning) and funded by the Korea Government (No. 2019-0-00079, Artificial Intelligence Graduate School Program, Korea University), and was partly supported by the German Ministry for Education and Research (BMBF) under Grants 01IS14013A-E, 01GQ1115, 01GQ0850, 01IS18025A, 031L0207D and 01IS18037A; the German Research Foundation (DFG) under Grant Math+, EXC 2046/1, Project ID 390685689.

AT acknowledges financial support from the European Research Council (ERC Consolidator Grant BeStMo and ERC-POC Grant DISCOVERER).

We gratefully acknowledge helpful comments on the manuscript by Hartmut Maennel.

Address correspondence to: jakeith@pitt.edu, klaus-robert.mueller@tu-berlin.de, alexandre.tkatchenko@uni.lu

Author Bios

John A. Keith is an associate professor and R.K. Mellon Faculty Fellow in Energy at the University of Pittsburgh in the department of chemical and petroleum engineering. He obtained his bachelors’ in chemistry at Wesleyan University and a Ph.D. degree in computational chemistry at Caltech in 2007. After an Alexander von Humboldt postdoctoral fellowship at the Universität Ulm, he was an Associate Research Scholar at Princeton University. He was a recipient of an NSF-CAREER award in 2017. His research interests lie in the applications and development of computational chemistry for engineering chemical reactions and materials for electrocatalysis, anticorrosion coatings, and the development of chemicals having less of an environmental footprint. He was a recipient of a Luxembourg Science Foundation INTER Mobility award in 2019-2020 to do a research sabbatical in Prof. Alexandre Tkatchenko’s group at the University of Luxembourg. This review is a primary product of that visit.

Valentin Vassilev-Galindo graduated with honors from University of Veracruz (Mexico) with a Bachelor’s degree in Chemical Engineering in 2014. Then, he enrolled to the Master program in Physical Chemistry at Cinvestav-Mérida (Mexico) where he worked under the supervision of Professor Gabriel Merino until receiving the MSc. degree in 2017. He is currently pursuing a PhD degree at the University of Luxembourg in the research group of Professor Alexandre Tkatchenko. His research is mainly related to machine learning potentials.

Bingqing Cheng is a Departmental Early Career Fellow at the Computer Laboratory, University of Cambridge, and a Junior Research Fellow at Trinity College. She received her Ph.D. from the École polytechnique fédérale de Lausanne (EPFL) in 2019. Her work focuses on theoretical predictions of material properties.

Stefan Chmiela is a senior researcher at the Berlin Institute for the Foundations of Learning and Data (BIFOLD). He received his Ph.D. from Technische Universität Berlin in 2019. His research interests include Hilbert space learning methods for applications in quantum chemistry, with particular focus on data efficiency and robustness.

Michael Gastegger is a postdoctoral researcher in the BASLEARN project of the Machine Learning Group at Technische Universität Berlin. He received his Ph.D. in Chemistry from the University of Vienna in Austria in 2017. His research interests include the development of machine learning methods for quantum chemistry and their application in simulations.

Klaus-Robert Müller has been a professor of computer science at Technische Universität Berlin since 2006; at the same time he is directing and co-directing the Berlin Machine Learning Center and the Berlin Big Data Center, respectively. He studied physics in Karlsruhe from 1984 to 1989 and obtained his Ph.D. degree in computer science at Technische Universität Karlsruhe in 1992. After completing a postdoctoral position at GMD FIRST in Berlin, he was a research fellow at the University of Tokyo from 1994 to 1995. In 1995, he founded the Intelligent Data Analysis group at GMD-FIRST (later Fraunhofer FIRST) and directed it until 2008. From 1999 to 2006, he was a professor at the University of Potsdam. He was awarded the Olympus Prize for Pattern Recognition (1999), the SEL Alcatel Communication Award (2006), the Science Prize of Berlin by the Governing Mayor of Berlin (2014), the Vodafone Innovations Award (2017). In 2012, he was elected member of the German National Academy of Sciences-Leopoldina, in 2017 of the Berlin Brandenburg Academy of Sciences and also in 2017 external scientific member of the Max Planck Society. In 2019 and 2020 he became a Highly Cited researcher in the cross-disciplinary area. His research interests are intelligent data analysis and Machine Learning in the sciences (Neuroscience (specifically Brain-Computer Interfaces), Physics, Chemistry) and in industry.

Alexandre Tkatchenko is a Professor of Theoretical Chemical Physics at the University of Luxembourg and Visiting Professor at Technische Universität Berlin. He obtained his bachelor degree in Computer Science and a Ph.D. in Physical Chemistry at the Universidad Autonoma Metropolitana in Mexico City. Between 2008 and 2010, he was an Alexander von Humboldt Fellow at the Fritz Haber Institute of the Max Planck Society in Berlin. Between 2011 and 2016, he led an independent research group at the same institute. Tkatchenko serves on editorial boards of two society journals: Physical Review Letters (APS) and Science Advances (AAAS). He received a number of awards, including elected Fellow of the American Physical Society, the 2020 Dirac Medal from WATOC, the Gerhard Ertl Young Investigator Award of the German Physical Society, and two flagship grants from the European Research Council: a Starting Grant in 2011 and a Consolidator Grant in 2017. His group pushes the boundaries of quantum mechanics, statistical mechanics, and machine learning to develop efficient methods to enable accurate modeling and obtain new insights into complex materials.

References