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

John A. Keith<EMAIL_ADDRESS>[  Valentin Vassilev-Galindo  [  Bingqing
Cheng  [  Stefan Chmiela  [  Michael Gastegger  [  Klaus-Robert Müller  klaus-
<EMAIL_ADDRESS>[  Alexandre Tkatchenko
<EMAIL_ADDRESS>[

###### Abstract

Machine learning models are poised to make a transformative impact on chemical
sciences by dramatically accelerating computational algorithms and amplifying
insights available from computational chemistry methods. However, achieving
this requires a confluence and coaction of expertise in computer science and
physical sciences. This review is written for new and experienced researchers
working at the intersection of both fields. We first provide concise tutorials
of computational chemistry and machine learning methods, showing how insights
involving both can be achieved. We then follow with a critical review of
noteworthy applications that demonstrate how computational chemistry and
machine learning can be used together to provide insightful (and useful)
predictions in molecular and materials modeling, retrosyntheses, catalysis,
and drug design.

Pitt]Department of Chemical and Petroleum Engineering Swanson School of
Engineering, University of Pittsburgh, Pittsburgh, USA Lux]Department of
Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg
City, Luxembourg Camb]Accelerate Programme for Scientific Discovery,
Department of Computer Science and Technology, 15 J. J. Thomson Avenue,
Cambridge CB3 0FD, United Kingdom, and Cavendish Laboratory, University of
Cambridge, J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom
TUBerlin]Department of Software Engineering and Theoretical Computer Science,
Technische Universität Berlin, 10587, Berlin, Germany TUBerlin]Department of
Software Engineering and Theoretical Computer Science, Technische Universität
Berlin, 10587, Berlin, Germany TUBerlin]Machine Learning Group, Technische
Universität Berlin, 10587, Berlin, Germany; Department of Artificial
Intelligence, Korea University, Anam-dong, Seongbuk-gu, Seoul, 02841, Korea;
Max-Planck-Institut für Informatik, Saarbrücken, Germany; Google Research,
Brain team, Berlin, Germany Lux]Department of Physics and Materials Science,
University of Luxembourg, L-1511 Luxembourg City, Luxembourg

###### Contents

  1. 1 Introduction
    1. 1.1 Background
    2. 1.2 Motivation for this review
  2. 2 CompChem and Notable Intersections with ML
    1. 2.1 Computational modeling, data, and information across many scales
      1. 2.1.1 Models and levels of abstraction
      2. 2.1.2 CompChem representations
      3. 2.1.3 Method accuracy
      4. 2.1.4 Precision and reproducibility
    2. 2.2 Hierarchies of methods
      1. 2.2.1 Wavefunction theory methods
      2. 2.2.2 Correlated wavefunction methods
      3. 2.2.3 Density-functional theory
      4. 2.2.4 Semi-empirical methods
      5. 2.2.5 Nuclear quantum effects
      6. 2.2.6 Interatomic Potentials
    3. 2.3 Response properties
    4. 2.4 Solvation models
    5. 2.5 Insightful predictions for molecular and material properties
  3. 3 Machine Learning Tutorial and Intersections with Chemistry
    1. 3.1 What is ML?
      1. 3.1.1 What is ML good at?
      2. 3.1.2 What is ML not good at?
    2. 3.2 Types of learning
      1. 3.2.1 Supervised learning
      2. 3.2.2 Unsupervised learning
      3. 3.2.3 Reinforcement learning
    3. 3.3 Universal approximators
    4. 3.4 The ML workflow
      1. 3.4.1 Datasets
      2. 3.4.2 Descriptors
      3. 3.4.3 Training
  4. 4 Applications of Machine Learning to Chemical Systems
    1. 4.1 Representing chemical systems
      1. 4.1.1 Descriptors
      2. 4.1.2 Representing local environments
      3. 4.1.3 The locality approximation
      4. 4.1.4 Advantages of built-in symmetries
      5. 4.1.5 End-to-end NN representations
    2. 4.2 From descriptors to predictions
    3. 4.3 CompChem data
      1. 4.3.1 Benchmark data sets
      2. 4.3.2 Visualization of data sets
      3. 4.3.3 Text and data mining for chemistry
    4. 4.4 Transforming atomistic modelling
      1. 4.4.1 Predicting thermodynamic properties
      2. 4.4.2 Nuclear quantum effects
    5. 4.5 ML for structure search, sampling, and generation
    6. 4.6 Multiscale modeling
  5. 5 Selected applications and paths toward insights
    1. 5.1 Molecular and material design
    2. 5.2 Retrosynthetic technologies
    3. 5.3 Catalysis
    4. 5.4 Drug design
  6. 6 Conclusions and Outlook
  7. 7 Author Bios

##  1 Introduction

###  1.1 Background

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:

  1. 1.

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

  2. 2.

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

  3. 3.

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.

  4. 4.

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.

Figure 1: Heatmaps illustrating the extent that ML terms appear in scientific
papers aligned by American Chemical Society (ACS) technical divisions from
2000-2010 a) and from 2010-present b). c) A line graph showing the number of
occurrences of any ML term being found in papers attributed to the ACS PHYS
division, from 2000-present. Figures were made by Charles D. Griego. Python
scripts used to generate these figures and corresponding Table 1 are freely
available with a creative commons attribution license. Readers are welcome to
use, adapt, and share these scripts with appropriate attribution:
https://github.com/keithgroup/scopus_searching_ML_in_chem_literature)

Table 1: List of top ranked keywords (per ACS division) with corresponding percentage of matches for any ML term. Division | Rank 1 | Rank 2 | Rank 3 | Rank 4
---|---|---|---|---
PHYS | electro* (56.3%) | spectroscopy (9.7%) | ion* (6.0%) | nano* (5.5%)
ANYL | sensor* (55.4%) | spectroscopy (13.1%) | characterization* (11.6%) | spectrometry (4.3%)
ENVR | *sensor* (60.9%) | soil* (14.5%) | water quality (4.2%) | environmental monitor* (3.0%)
AGFD | protein* (31.9%) | agricultur* (18.2%) | food (10.8%) | fruit* (5.7%)
ENFL | fuel* (19.2%) | petroleum (11.6%) | energy efficiency (11.1%) | batter* (10.7%)
AGRO | soil (43.3%) | crop* (25.4%) | groundwater (11.5%) | developing countr* (4.4%)
ORGN | protein* (64.6%) | amino acid* (19.7%) | peptide* (8.4%) | aromatic* (3.2%)
POLY | plastic* (51.1%) | polymer* (37.7%) | polymeriz* (5.0%) | polymeric (2.8%)
PMSE | *polymer* (50.4%) | *peptide* (30.4%) | thin film* (8.9%) | tissue engineering (4.3%)
BIOT | biochemi* (37.9%) | biophysic* (18.3%) | systems biology (10.3%) | biotechnology (9.9%)
GEOC | groundwater (33.4%) | mining (31.6%) | *geochem* (12.4%) | anthropogenic (10.9%)
MEDI | protein interaction* (25.7%) | drug discovery (19.1%) | drug design (19.1%) | antibiotic* (11.3%)
COMP | drug discovery (18.7%) | drug design (18.6%) | molecular model* (14.3%) | protein database* (13.1%)
COLL | nanoparticle* (21.2%) | adsorption (19.6%) | thin film* (14.9%) | tribolog* (9.6%)
BIOL | drug discovery (41.9%) | protein folding (19.3%) | biosynthesis (12.2%) | cytochrome* (12.2%)
TOXI | toxi* (99.2%) | chemical exposure* (0.6%) | antibody drug conjugate* (0.1%) |
CATL | cataly* (64.0%) | metal oxides (20.2%) | photocataly* (5.3%) | surface chemistry (2.8%)
CINF | drug discovery (51.7%) | computational chemistry (17.5%) | bio* modeling (7.8%) | chem* database* (7.6%)
INOR | electrochem* (67.1%) | nanomaterial* (14.9%) | organometallic* (5.8%) | metal organic framework* (4.0%)
NUCL | nuclear fuel* (28.6%) | isotope* (27.3%) | radioisotope* (15.0%) | nuclear medicine* (9.0%)
CARB | carbohydrate* (43.0%) | glycoprotein* (42.4%) | glycan* (6.6%) | oligosaccharide* (5.6%)
RUBB | rubber* (100.0%) |  |  |
CELL | cellulose (41.7%) | polysaccharide* (26.3%) | lignin (16.3%) | lignocellulos* (9.7%)
I&EC | water purification (38.1%) | industrial chem* (23.7%) | rare earth element* (11.9%) | industrial and engineering chemistry (8.3%)
FLUO | fluorine* (99.8%) | radiopharmaceutical chem* (0.2%) |  |
CHED | chem* class* (76.4%) | chem* communication* (8.5%) | chem* educat* (5.5%) | lab* safety (5.5%)
CHAS | chem* safety (51.5%) | lab* safety (16.7%) | environmental health and safety (16.7%) | chem* regulations (15.2%)
BMGT | chem* compan* (65.4%) | chem* enterprise* (28.8%) | chem* business* (3.8%) | chem* research and development (1.9%)
SCHB | commercial chem* (50.0%) | chem* sector* (28.6%) | academic entrepreneur* (14.3%) | science advoca* (7.1%)
HIST | chem* histor* (53.8%) | evolution of chem* (30.8%) | history of chem* (15.4%) |
PROF | chem* education (100.0%) |  |  |
CHAL | pharmaceutical patent* (60.0%) | chem* in commerce (30.0%) | chem* patent* (10.0%) |

###  1.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.

Figure 2: A data-information-knowledge-insight-wisdom-impact hierarchy
courtesy of the Gapingvoid Culture Design Group,
https://www.gapingvoidart.com.

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.

##  2 CompChem and Notable Intersections with ML

###  2.1 Computational modeling, data, and information across many scales

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.

####  2.1.1 Models and levels of abstraction

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.

Figure 3: Hierarchy of computational methods and corresponding time and length
scales

####  2.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.

####  2.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

####  2.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.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.

Figure 4: Potential energy surface (PES) of a fictional system with the two
coordinates $R_{1}$ and $R_{2}$. The minima of the PES correspond to stable
states of a system, such as equilibrium configurations and reactants or
products. Minima can be connected by paths (red line), along which
rearrangements and reactions can occur. The maximum along such a path is
called a transition state. Transition states are first order saddle points, a
maximum in one coordinate and minima in all others. They correspond to the
minimum energy required to transition between two PES minima and play a
crucial role in the description of chemical transformations.

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.

Figure 5: a) A ‘magic cube’58 depiction of hierarchies of correlated
wavefunction approaches. b) A ‘Jacob’s Ladder’59 depiction of hierarchies of
Kohn-Sham DFT approaches. c) An illustration of expected hierarchies of
atomistic potentials adapted from Ref. 60. d) An illustration of overall
hierarchies in predictive atomic scale modeling methods.

####  2.2.1 Wavefunction theory methods

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:

$\hat{H}\Psi=(\hat{T}+\hat{V})\Psi=E\Psi,$ (1)

Here, the Hamiltonian operator ($\hat{H}$) is the sum of the kinetic
($\hat{T}$) and potential ($\hat{V}$) operators, $\Psi$ is the wavefunction
(i.e. an eigenfunction) that represents particles in the system, and $E$ is
the energy (i.e. an eigenvalue). In this way, nuclei can be treated as fixed
point charges, and then Eq. 1 can be transformed into the so-called electronic
Schrödinger equation, where the Hamiltonian $\hat{H}_{\mathrm{el}}$ and
wavefunction $\Psi_{\mathrm{el}}(\mathbf{r};\mathbf{R})$ now only depend on
the nuclear coordinates $\mathbf{R}$ in a parametric fashion:

$\displaystyle\hat{H}_{\mathrm{el}}\Psi_{\mathrm{el}}(\mathbf{r;R})$
$\displaystyle=\left[\hat{T}_{\mathrm{e}}(\mathbf{r})+\hat{V}_{\mathrm{eN}}(\mathbf{r;R})+\hat{V}_{\mathrm{NN}}(\mathbf{R})+\hat{V}_{\mathrm{ee}}(\mathbf{r})\right]\Psi_{\mathrm{el}}(\mathbf{r;R})$
$\displaystyle=E_{\mathrm{el}}\Psi_{\mathrm{el}}(\mathbf{r;R})$ (2)

The above expression has $\hat{H}_{\mathrm{el}}$ composed of single electron
($\mathrm{e}$), electron-nuclear ($\mathrm{eN}$), nuclear-nuclear
($\mathrm{NN}$), and electron-electron ($\mathrm{ee}$) terms. Here, we will
now implicitly assume the Born-Oppenheimer approximation throughout and leave
off the subscript indicating the electronic problem. However, we note that the
Born–Oppenheimer approximation is not always sufficient and the
computationally intensive nonadiabatic quantum dynamics may be required.65 In
certain cases, semi-classical treatments are appropriate, for example,
nonadiabatic effects between electrons and nuclei can be considered using
nuclear-electronic orbital methods.66

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:

$\Psi(\mathbf{r}_{1},\ldots,\mathbf{r}_{n})=\frac{1}{\sqrt{n!}}\begin{vmatrix}\phi_{1}(\mathbf{r}_{1})&\ldots&\phi_{n}(\mathbf{r}_{1})\\\
\vdots&\ddots&\vdots\\\
\phi_{1}(\mathbf{r}_{n})&\ldots&\phi_{n}(\mathbf{r}_{n})\\\ \end{vmatrix}$ (3)

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

$\phi_{i}=\sum_{\mu}c_{\mu i}\chi_{\mu}$ (4)

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:

$\displaystyle E_{\text{HF}}=$ $\displaystyle-\sum_{i}\int
d\mathbf{r}_{1}\phi^{*}_{i}(\mathbf{r}_{1})\underbrace{\frac{1}{2}\nabla^{2}}_{\hat{T}_{\mathrm{e}}}\phi_{i}(\mathbf{r}_{1})-\sum_{i}\int
d\mathbf{r}_{1}\phi^{*}_{i}(\mathbf{r}_{1})\underbrace{\sum^{N}_{\alpha}\frac{Z_{\alpha}}{|\mathbf{r}_{1}-\mathbf{R}_{\alpha}|}}_{\hat{V}_{\mathrm{eN}}}\phi_{i}(\mathbf{r}_{1})$
$\displaystyle+\sum_{i,j>i}\int\int
d\mathbf{r}_{1}d\mathbf{r}_{2}\phi^{*}_{i}(\mathbf{r}_{1})\phi_{i}(\mathbf{r}_{1})\underbrace{\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}}_{\hat{V}_{\mathrm{ee}}(\mathrm{Coulomb})}\phi^{*}_{j}(\mathbf{r}_{2})\phi_{j}(\mathbf{r}_{2})$
$\displaystyle-\sum_{i,j>i}\int\int
d\mathbf{r}_{1}d\mathbf{r}_{2}\phi^{*}_{i}(\mathbf{r}_{1})\phi_{j}(\mathbf{r}_{1})\underbrace{\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}}_{\hat{V}_{\mathrm{ee}}(\mathrm{exchange})}\phi^{*}_{j}(\mathbf{r}_{2})\phi_{i}(\mathbf{r}_{2})$
(5)

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:

$\mathbf{F}\mathbf{C}=\boldsymbol{\epsilon}\mathbf{S}\mathbf{C}$ (6)

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.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:

$\Psi_{\mathrm{FCI}}=a_{0}\Psi_{HF}+\sum_{\alpha,\beta}a^{\alpha}_{\beta}\Psi^{\alpha}_{\beta}+\sum_{\alpha,\beta,\gamma,\delta}a^{\alpha\gamma}_{\beta\delta}\Psi^{\alpha\gamma}_{\beta\delta}+\ldots$
(7)

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.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$:

$E[\rho]=T[\rho]+V[\rho]$ (8)

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:

$E[\rho]=T_{\mathrm{ni}}[\rho]+V_{\mathrm{eN}}[\rho]+V_{\mathrm{ee}}[\rho]+\Delta
T_{\mathrm{ee}}[\rho]+\Delta V_{\mathrm{ee}}[\rho]$ (9)

Here, $T_{\mathrm{ni}}[\rho]$ is the kinetic energy of the non-interacting
electrons, $V_{\mathrm{eN}}[\rho]$ is the exact nuclear-electron potential and
$V_{\mathrm{ee}}[\rho]$ is the Coulombic (classical) energy of the non-
interacting electrons. The last two terms are corrections due to the
interacting nature of electrons and non-classical electron-electron repulsion.
KS-DFT also expands the three dimensional electron density into a spin
orbital-basis $\phi$ similar to HF theory to define the one-electron kinetic
energy in a straightforward manner. This allows the $T_{\mathrm{ni}}$,
$V_{\mathrm{eN}}$ and $V_{\mathrm{ee}}$ expressions to be evaluated exactly
and one arrives at the KS energy:

$\displaystyle E_{\text{KS}}[\rho]=$ $\displaystyle-\underbrace{\sum_{i}\int
d\mathbf{r}_{1}\phi^{*}_{i}(\mathbf{r}_{1})\frac{1}{2}\nabla^{2}\phi_{i}(\mathbf{r}_{1})}_{{T}_{\mathrm{ni}}}-\underbrace{\sum_{i}\int
d\mathbf{r}_{1}\phi^{*}_{i}(\mathbf{r}_{1})\sum^{N}_{\alpha}\frac{Z_{\alpha}}{|\mathbf{r}_{1}-\mathbf{R}_{\alpha}|}\phi_{i}(\mathbf{r}_{1})}_{V_{\mathrm{eN}}}$
$\displaystyle+\underbrace{\sum_{i}\int
d\mathbf{r}_{1}\phi^{*}_{i}(\mathbf{r}_{1})\left[\frac{1}{2}\int
d\mathbf{r}^{{}^{\prime}}\frac{\rho(\mathbf{r}^{{}^{\prime}})}{|\mathbf{r}_{1}-\mathbf{r}^{{}^{\prime}}|}\right]\phi_{i}(\mathbf{r}_{1})}_{\hat{V}_{\mathrm{ee}}(\mathrm{Coulomb})}+E_{\text{xc}}[\rho]$
(10)

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.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:

$E[\delta\rho_{0}]=E_{\mathrm{rep}}+E_{\mathrm{Coul}}[\delta\rho_{0}]+E_{\mathrm{BS}}[\delta\rho_{0}].$
(11)

$E_{\mathrm{rep}}$ is a repulsive potential containing interactions between
the nuclei and contributions from the exchange correlation functional (these
are typically approximated via pairwise potentials). The charge fluctuation
term $E_{\mathrm{Coul}}$ is modeled as a Coulomb potential of Gaussian charge
distributions computed from the approximate density. Finally,
$E_{\mathrm{BS}}$ refers to the ‘band structure’ term which considers the
electronic structure and contains contributions from $T_{\mathrm{ni}}$,
$V_{\mathrm{eN}}$, and the exchange correlation functional (see Eq. 2.2.3). To
compute $E_{\mathrm{BS}}$, the density is expressed in a minimal basis of
atomic orbitals, similar as in NDDO. The necessary Hamiltonian and overlap
integrals are then evaluated via an approximate scheme based on Slater–Koster
transformations. In addition to the energy, atomic partial charges are also
computed in this step, which are then used in $E_{\mathrm{Coul}}$. As a
consequence, DFTB equations can also be solved self-consistently. DFTB methods
are parameterized by finding suitable forms for the repulsive potential and
adjusting the parameters used in the Slater–Koster integrals. Non-self
consistent and self-consistent tight-binding DFT methods193, 194 have been
developed for simulating large scale systems. Semi-empirical methods have also
been a target of different ML schemes, yielding improved parametrization
schemes and more accurate functional approximations.195, 196, 197, 198

####  2.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.2.6 Interatomic Potentials

Table 2: Types of interatomic potentials and their areas of application. Potential |  Reactive |  Typical applications |  Examples
---|---|---|---
Pairwise-distance-based |  sometimes |  materials, liquids |  Lennard-Jones 208, 209, Morse 210, Coulomb interactions, Buckingham 211
Distance and angle-based |  usually no |  materials, liquids |  many water potentials (e.g. SPC, TIP4P, mW)212, Stillinger-Weber 213
Class I (non-polarizable) force fields |  no |  Proteins, lipids, polymers, nucleic acids, carbohydrates, organic molecules, liquids |  AMBER,214, 215 GAFF,216 CHARMM,217 GROMOS,218, 219, 220 OPLS,221, 222 DREIDING,223 MMFF94,224 UFF,225 COMPASS,226 INTERFACE227, interatomic potentials for ionic systems228
Class II (polarizable) |  no |  Proteins, lipids, polymers, nucleic acids, carbohydrates, organic molecules, liquids |  AMOEBA,229 classical Drude oscillator models,230 fluctuating charge (FQ) models,231 MB-Pol,232 distributed point polarizable models (DPP2),233 and many more.234
Embedded atom method (EAM)-like |  yes |  Reactions within solid materials |  EAM,235 MEAM,236 Finnis-Sinclair,237 Sutton-Chen238
Bond-order potentials (BOPs) |  yes |  Reactions within solids, liquids, gases |  Brenner,239 Tersoff,240, 241 REBO,239, 242 COMB,243, 244 ReaxFF,245, 246 APT247
Other quantum mechanics-derived forcefields |  yes |  Reactions within liquids and gases |  EVB248 and related models249, 250

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

$E_{\mathrm{LJ}}=\sum_{i,j>i}4\varepsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right].$
(12)

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:

$\displaystyle E_{\mathrm{tot}}=$
$\displaystyle\sum_{\mathrm{bonds}}k_{ij}\left(r_{ij}-\bar{r}_{ij}\right)^{2}+$
$\displaystyle\sum_{\mathrm{angles}}k_{ijk}\left(\theta_{ijk}-\bar{\theta}_{ijk}\right)^{2}+$
$\displaystyle\sum_{\mathrm{dihedral}}\sum_{\gamma}k_{ijkl}^{(\gamma)}\left[1-\cos(\gamma\phi_{ijkl}-\bar{\phi}_{ijkl}^{(\gamma)})\right]+$
$\displaystyle\sum_{i,j>i}\frac{q_{i}q_{j}}{r_{ij}}+E_{\mathrm{LJ}}$ (13)

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

$E_{\mathrm{tot}}=\sum_{i}^{N}\left[F_{i}(\tilde{\rho}_{i})+\frac{1}{2}\sum_{j\neq
i}V_{ij}(r_{ij})\right].$ (14)

$F_{i}$ is an embedding function and $\tilde{\rho}_{i}$ an approximation to
the local electron density based on the environment of atom $i$.
$F_{i}(\tilde{\rho}_{i})$ can be seen as a contribution due to non-localized
electrons in a metal. $V_{ij}$ is a term describing to the core-core repulsion
between atoms. An EAM potential is determined by the functional forms used for
$F_{i}$ and $V_{ij}$ as well as how the density is expressed. Its dependence
on the local environment without the need for predefined bonds make EAM well
suited for modeling material properties of metals. An extension of EAM is
modified EAM (MEAM)236 which includes directional dependence in the
description of the local density $\tilde{\rho}_{i}$, but this brings greater
computational cost. EAMs also form the conceptual basis of the embedded atom
neural network (EANN) MLPs.269

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:

$E_{\mathrm{tot}}=\sum_{i,j>i}\left[V_{\mathrm{rep}}(r_{ij})-b_{ij(k)}V_{\mathrm{att}}(r_{ij})\right]f_{\mathrm{cut}}(r_{ij}).$
(15)

$V_{\mathrm{rep}}$ and $V_{\mathrm{att}}$ are repulsive and attractive
potentials depending on the interatomic distance $r_{ij}$. A cutoff function
$f_{\mathrm{cut}}$ restricts all interactions to the local atomic environment.
$b_{ij(k)}$ is the bond order term, from which the potential takes its name.
This term measures the bond order between atoms $i$ and $j$ (i.e., ‘1’ for a
single bond, ‘2’ for a double bond, ‘0.6’ for partially dissociated bond).
Bond orders can also depend on neighboring atoms $k$ in some implementations.
BOPs are typically used for covalently bound systems, such as bulk solids and
liquids containing hydrogen, carbon or silicon (e.g. carbon nanotubes and
graphene). Depending on the exact form of the expressions in Eq. 15, different
types of BOPs are obtained, such as Tersoff241, 240 and REBO239, 242
potentials. BOPs can also be extended to incorporate dynamically assigned
charges, yielding potentials like COMB243, 270 or ReaxFF.245, 246 As with
EAMs, BOPs have also been used as a starting point for constructing more
elaborate machine learning potentials (MLPs)271, 272, 273 that will also be
discussed in more detail in Section 4.

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.

###  2.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

$\boldsymbol{\Pi}(n_{\mathbf{R}},n_{\boldsymbol{\epsilon}},n_{\mathbf{B}},n_{\mathbf{I}_{i}})=\frac{\partial^{n_{\mathbf{R}}+n_{\boldsymbol{\epsilon}}+n_{\mathbf{B}}+n_{\mathbf{I}_{i}}}E(\mathbf{R},\boldsymbol{\epsilon},\mathbf{B},\mathbf{I}_{i})}{\partial\mathbf{R}^{n_{\mathbf{R}}}\partial\boldsymbol{\epsilon}^{n_{\boldsymbol{\epsilon}}}\partial\mathbf{B}^{n_{\mathbf{B}}}\partial\mathbf{I}_{i}^{n_{\mathbf{I}_{i}}}},$
(16)

where the $n$s indicate the $n$-th order partial derivative with respect to
the quantity in the subscript.102

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).

Table 3: Response properties of the potential energy.

$n_{\mathbf{R}}$ | $n_{\boldsymbol{\epsilon}}$ | $n_{\mathbf{B}}$ | $n_{\mathbf{I}}$ | Property
---|---|---|---|---
0 | 0 | 0 | 0 | Energy
1 | 0 | 0 | 0 | Forces
2 | 0 | 0 | 0 | Hessian (harmonic frequencies)
0 | 1 | 0 | 0 | Dipole moment (IR)
1 | 1 | 0 | 0 | Infrared absorption intensities (IR)
0 | 2 | 0 | 0 | Polarizability (Raman)
0 | 0 | 1 | 1 | Nuclear magnetic shielding (NMR)
0 | 0 | 0 | 2 | Nuclear spin-spin coupling (NMR)
0 | 1 | 1 | 0 | Optical rotation (circular dichroism)

$\vdots$

Infrared spectra can be simulated based on dipole moments
$\boldsymbol{\mu}=-\boldsymbol{\Pi}(0,1,0,0)$, while molecular
polariziabilities $\boldsymbol{\alpha}=-\boldsymbol{\Pi}(0,2,0,0)$ offer
access to polarized and depolarized Raman spectra. A central response property
of the magnetic field are nuclear magnetic shielding tensors
$\boldsymbol{\sigma}_{i}=\boldsymbol{\Pi}(0,0,1,1)$. These allow the
computation of chemical shifts recorded in nuclear magnetic resonance (NMR)
spectroscopy via their trace
$\sigma_{i}=\frac{1}{3}\mathrm{Tr}[\boldsymbol{\sigma}_{i}]$. The beauty of
this formalism lies in the fact, that a single energy calculation method
provides access to a wide range of quantum chemical properties in a highly
systematic manner. A large number of modern MLPs use the response of the
potential energy with respect to nuclear positions to obtain energy conserving
forces. However, far fewer applications model perturbations with respect to
electric and magnetic fields. Ref. 299 extends the descriptor used in the FCHL
Kernel by an explicit field dependent term making it possible to predict
dipole moments across chemical compound space. Ref. 300 introduces a general
NN framework to model interactions of a system with vector fields, which was
then used to predict dipole moments, polarizabilities and nuclear magnetic
shielding tensors as response properties.

###  2.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:

$-\nabla\left[\epsilon(\mathbf{r})\nabla
V(\mathbf{r})\right]=4\pi\rho_{\mathrm{m}}(\mathbf{r})$ (17)

Here, $\rho_{\mathrm{m}}(\mathbf{r})$ is the charge distribution of the solute
and $\epsilon(\mathbf{r})$ is the position dependent permittivity, which
usually is set to one within the cavity and the $\varepsilon$ of the solvent
on the outside. $V(\mathbf{r})$ is the electrostatic potential composed of the
two terms

$V(\mathbf{r})=V_{\mathrm{m}}(\mathbf{r})+V_{\mathrm{s}}(\mathbf{r}),$ (18)

where $V_{\mathrm{m}}(\mathbf{r})$ is the solute potential and
$V_{\mathrm{s}}(\mathbf{r})$ is the apparent potential due to the surface
charge distribution $\sigma(\mathbf{s})$,

$V_{\mathrm{s}}(\mathbf{r})=\int_{\Gamma}\frac{\sigma(\mathbf{s})}{|\mathbf{r}-\mathbf{s}|}d\mathbf{s}.$
(19)

$\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

$\hat{H}_{\mathrm{eff}}=\hat{H}+V_{\mathrm{s}}(\mathbf{r}),$ (20)

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.

###  2.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.

##  3 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.

###  3.1 What is ML?

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

$\hat{f}=\underset{f\in\mathcal{F}}{\arg\min}\left[\sum_{i}^{M}\mathcal{L}\left(f\left(\mathbf{x}_{i}\right),y_{i}\right)\right],$
(21)

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

$\hat{f}=\underset{f\in\mathcal{F}}{\arg\min}\left[\sum_{i}^{M}\mathcal{L}\left(f\left(\mathbf{x}_{i}\right),y_{i}\right)+\|\Gamma\Theta\|^{2}\right].$
(22)

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

$\hat{f}=\underset{f\in\mathcal{F}}{\arg\min}\left[\int\mathcal{L}\left(f\left(\mathbf{x}_{i}\right),y_{i}\right)dp(\mathbf{x},\mathbf{y})\right],$
(23)

that we would actually like to solve. In an ideal world, we would minimize the
loss function over the _complete_ distribution of inputs and labels
$p(\mathbf{x},\mathbf{y})$. However, this is obviously impossible in practice,
so we apply the principle of Occam’s razor that presumes that simpler
(parsimonious) hypotheses are more likely to be correct. With this additional
consideration we hope to be able to recover a reasonably general model,
despite only having seen a finite training set. A common way to favor simpler
models is via an additional term in the cost function, which is what
$\|\Gamma\Theta\|^{2}$ in Eq. 22 expresses. Here, $\Gamma$ is a matrix that
defines “simplicity” with regard to the model parameters $\Theta$. Usually,
$\Gamma=\lambda\mathbb{I}$ (with $\lambda>0$) is chosen to simply favor a
small L2-norm on the parameters, such that the solution does not rely on
individual input features too strongly. This particular approach is called
Tikhonov regularization, 346, 347, 348 but other regularization techniques
also exist. 349, 350

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).

####  3.1.1 What is ML good at?

##### Implicit knowledge from data

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

##### Redundancy in first principles calculations

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.

####  3.1.2 What is ML not good at?

##### Lack of generality and precision

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.

##### Reliance on high-quality data

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.

##### Inability to derive high-level concepts

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

##### Prone to artifacts

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.

###  3.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.

####  3.2.1 Supervised learning

_Supervised learning_ addresses learning problems where the ML model
$\hat{f}:\mathcal{X}\stackrel{{\scriptstyle\mathrm{ML}}}{{\longrightarrow}}\mathcal{Y}$
connects a set of known inputs $\mathcal{X}$ and outputs $\mathcal{Y}$, either
to perform a regression or classification task. While the former maps onto a
continuous space (e.g. energy, polarizability), the latter outputs a
categorical value (e.g. acid or base; metal or insulator) for each data point.

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

####  3.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.

####  3.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.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

$f({\bf x})={\rm sign}(\sum_{i=1}^{N}w_{i}x_{i}-b)$ (24)

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

$f_{k}({\bf x})={\rm g}\left(\sum_{j=1}^{H}w_{kj}{\rm
g}\left[\sum_{i=1}^{N}w_{ji}x_{i}-b\right]\right),$ (25)

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:

$f({\bf x})=\sum_{j=1}^{N}y_{j}\alpha_{j}K({\bf x}_{j},{\bf x})-b,$ (26)

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

###  3.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:

  * •

Gathering and preparing the data

  * •

Choosing a representation

  * •

Training the model

    * –

Train model condidates

    * –

Evaluate model accuracy

    * –

Tune hyperparameters

  * •

Testing the model out of sample

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.

####  3.4.1 Datasets

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.

####  3.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.

####  3.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.

##### Kernel-based learning algorithms

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

$\nabla_{\vec{\alpha}}\mathcal{L}(\hat{f}(\mathbf{x}),\mathbf{y})=(\mathbf{K}+\lambda\mathbb{I})\vec{\alpha}-\mathbf{y}=0$
(27)

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:

$\vec{\alpha}^{t}=\vec{\alpha}^{t-1}-\gamma\underbrace{\nabla_{\vec{\alpha}}\mathcal{L}(\hat{f}(\mathbf{x}),\mathbf{y})}_{\text{e.g.
}(\mathbf{K}+\lambda\mathbb{I})\vec{\alpha}^{t-1}-\mathbf{y}},$ (28)

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).

##### Neural networks

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

##### Hyper-parameters

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.

##### Model selection

Figure 6: Supervised learning algorithms have to balance two sources of error
during training: the bias and variance of the model. A highly biased model is
based on flawed assumptions about the problem at hand (under-fitting).
Conversely, a high variance causes a model to follow small variations in the
data too closely, therefore making it susceptible to picking up random noise
(over-fitting). The optimal bias-variance trade-off minimizes the
generalization error of the model, e.g. how well it performs on unknown data.
It can be estimated with cross-validation techniques.

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

##  4 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.

###  4.1 Representing chemical systems

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

$\psi(\\{\mathcal{R}^{3N},\mathcal{Z}^{N}\\})=\psi(S(\\{\mathcal{R}^{3N},\mathcal{Z}^{N}\\})),$
(29)

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

####  4.1.1 Descriptors

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.

Table 4 summarizes representations for chemical systems. 411, 418, 276, 412,
417, 452, 453, 415 All respect the aforementioned physical symmetries and
invariances needed for chemical systems. Many have similar theoretical
foundations that can be understood as the basis onto which the atomic density
is projected,454 and the connection between them has been summarised in a
recent review. 455 The descriptors in Table 4 can be classified into two
categories: global and atomic, as indicated under the column “globality”.
Traditional descriptors used in cheminformatics are global descriptors based
on the covalent connectivity of atoms. These include simple valence counting