# Atomistic Descriptor Optimization Using Complementary Euclidean and Geodesic
Distance Information

Gopal R. Iyer<EMAIL_ADDRESS>Brenda M. Rubenstein  Department of
Chemistry, Brown University, Providence, Rhode Island 02912, USA
<EMAIL_ADDRESS>
###### Abstract

Descriptors are physically-inspired, symmetry-preserving schemes for
representing atomistic systems that play a central role in the construction of
models of potential energy surfaces. Although physical intuition can be
flexibly encoded into descriptor schemes, they are generally ultimately guided
only by the spatial or topological arrangement of atoms in the system.
However, since interatomic potential models aim to capture the variation of
the potential energy with respect to atomic configurations, it is conceivable
that they would benefit from descriptor schemes that implicitly encode both
structural and energetic information rather than structural information alone.
Therefore, we propose a novel approach for the optimization of descriptors
based on encoding information about geodesic distances along potential energy
manifolds into the hyperparameters of commonly used descriptor schemes. To
accomplish this, we combine two ideas: (1) a differential-geometric approach
for the fast estimation of approximate geodesic distances [Zhu et al. J. Chem.
Phys. 150, 164103 (2019)]; and (2) an information-theoretic evaluation metric
- information imbalance \- for measuring the shared information between two
distance measures [Glielmo et al. PNAS Nexus, 1, 1 (2022)]. Using three
example molecules - ethanol, malonaldehyde, and aspirin - from the MD22
dataset, we first show that Euclidean (in Cartesian coordinates) and geodesic
distances are inequivalent distance measures, indicating the need for updated
ground-truth distance measures that go beyond the Euclidean (or, more broadly,
spatial) distance. We then utilize a Bayesian optimization framework to show
that descriptors (in this case, atom-centered symmetry functions) can be
optimized to maximally express a certain type of distance information, such as
Euclidean or geodesic information. We also show that modifying the Bayesian
optimization algorithm to minimize a combined objective function - the sum of
the descriptor$\leftrightarrow$Euclidean and
descriptor$\leftrightarrow$geodesic information imbalances - can yield
descriptors that not only optimally express both Euclidean and geodesic
distance information simultaneously, but in fact resolve substantial
disagreements between descriptors optimized to encode only one type of
distance measure. We discuss the relevance of our approach to the design of
more physically rich and informative descriptors that can encode useful,
alternative information about molecular systems.

supporting

##  1 Introduction

Machine learning-based interatomic potentials (MLPs) have emerged as
computationally tractable solutions to the problem of modeling atomistic
systems with high accuracy at larger length- and longer time-scales than are
typically accessible to direct ab initio calculations. By training flexible
functional forms, such as neural networks on large pre-existing collections of
high-accuracy ab initio data, MLPs have successfully modeled a very large
number of atomistic systems with high fidelity.1, 2, 3, 4, 5

Many MLP models are structured like standard machine learning models, which
may be summarized as: Model = Representation + Architecture, where the
architecture refers to the functional form (e.g., a neural network) that
processes the representation of the system to predict the output. In the
context of MLPs, the representation is usually a physics-inspired
featurization of the atomic positions referred to as a descriptor. Descriptors
are typically chosen to be functions of atomic positions that are invariant to
a global rotation, translation, and permutation of identical atoms.
Additionally, they may be chosen to encode some aspect of the physics of the
specific atomistic system being studied, such as local structural
arrangements. An important variation on this theme are kernel methods that
estimate the relative energies of atomic configurations using featurized
kernels of the structural similarity between them.1, 6, 7, 8

The field of descriptor design is vast and an exhaustive summary of the
various strategies merits an independent review in its own right; see, for
example, Ref. 9. In this work, we limit our discussion to atom-centered
symmetry function (ACSF) descriptors, although, as we note later, our
approaches readily generalize to other descriptor formulations. ACSF
descriptors are among the first representation schemes to be have been
demonstrated successfully for atomistic neural network potentials.10, 11 As
the name suggests, ACSF descriptors involve constructing atom-wise functions
that encode the local neighborhood of each atom (labeled $i$) in terms of
different local structures, such as its pair-wise separation from its
neighboring atoms ($r_{ij}$), the angles it forms with its neighboring
triplets ($\theta_{ijk}$), etc.10 Typical implementations of atomistic neural
networks, termed Behler–Parrinello neural networks, then involve passing these
representations to several individual neural networks, each corresponding to a
single chemical element contained in the system, the outputs of which are then
added to give a prediction for the energy of the atomic configuration. Force
calculations are performed by taking the derivative of the energy with respect
to the atom positions indirectly through derivatives with respect to the
descriptors.11

The construction of atomistic descriptors helps to determine whether or not
two distinct atomic configurations encode similar physics, either because they
are related by symmetry, or because the interactions between their atomic
components lead to similar expressions for the energy - the most commonly
modeled quantity - or other observables. In this work, we ask the following
question: Is it possible to encode information beyond structural features -
such as the approximate intrinsic geometry of the potential energy surface -
within a given descriptor scheme? If a descriptor scheme can be made to
express the rich geometric structure of the potential energy manifold, the
atomic fingerprints generated by such a descriptor will then, in principle, be
distributed according to structural and energetic similarity, making such
representations more amenable to accurate modeling via an interatomic
potential. Including energetic information content in descriptors is of
particular relevance for the efficient sampling of structures for the training
of data-driven models of force fields1, 3 and other quantities. When training
such models with the aim of representing reactive dynamics, it is often
crucial that higher energy structures surrounding transition states and
intermediates along reaction pathways are represented in the training data,
yet Euclidean distances alone often do not directly provide information about
the energetics or spacing of configurations along reaction coordinates.
Several schemes for structure sampling along reaction coordinates have been
proposed, both in the context of materials modeling12, 13 and biomolecular
simulation,14 with the latter strongly focusing on the sampling of
energetically favorable regions (configurations) separated by low-probability
structural transformations, nonetheless, gaps remain. The need to better
sample along reaction coordinates also raises the important question of which
physical quantity is the most informative for describing geometric information
about the potential energy manifold given a limited number of known atomic
configurations of interest.

Given the fundamental importance of distance measures in describing
geometries, we begin by noting that distance may serve as a meaningful,
readily interpretable candidate for such a physical quantity. Several distance
metrics have been proposed for atomistic systems, including direct Euclidean
distance between Cartesian coordinate representations (assuming prior
rotational alignment, e.g., via the Kabsch algorithm) and kernel-based
similarity,1, 6 including combined structural/alchemical dissimilarity
kernels.8 Despite significant developments in this regard, most popular
similarity metrics are fundamentally grounded in the spatial structure of
atomic arrangements, largely due to the apparent lack of any other type of
information that can be gleaned from the structure a priori. For instance,
when measuring the expressivity of ACSF and SOAP (smooth overlap of atomic
positions) descriptors, Glielmo et al. (Ref. 15) used Euclidean distance as
the ground-truth distance measure to assess whether ACSF descriptors were able
to encode distance information to a sufficient degree.

Here, we argue that the approximate geodesic distance between atomic
configurations along a potential energy manifold - which, as we will note, can
be computed relatively easily - is an underutilized, yet physically
informative, distance measure that contains information that may be
fundamentally inaccessible to purely structure-based distance measures.
Geodesic distances have been used in the context of dimensionality reduction
for unsupervised atomistic machine learning, for example, in the ISOMAP
algorithm.16 For the computation of approximate geodesic distances, we rely on
the differential geometric approach developed by Zhu et al. in Ref. 17 in the
context of reaction pathway interpolation. We then show that it is possible to
design and optimize descriptors to be simultaneously capable of encoding the
information contained in both Euclidean and geodesic distances, thereby
resolving the information imbalance15 between these contradictory metrics via
a single optimal representation. For this purpose, we utilize a Bayesian
optimization algorithm to query and optimize the hyperparameters of three
different families of ACSF descriptors (chosen to have different levels of
expressivity). We demonstrate our approach on three molecular sub-datasets of
the MD22 dataset - ethanol, malonaldehyde, and aspirin.

The rest of this work is organized as follows. We describe the details of our
approach in Section 2, including a description of the geodesic distance
interpolation method (Section 2.2), the details of the Bayesian optimization
algorithm (Section 2.4), and the information-theoretic objective function for
the Bayesian optimization originally proposed in Ref. 15 (Section 2.3).
Further computational details are provided in Section 3. In Section 4, we
first demonstrate the inequivalence of Euclidean and geodesic distances, and
then discuss the application of the Bayesian optimization approach to encode
Euclidean and geodesic distance information into the ACSF descriptors
constructed in Section 3.1. We conclude and present an outlook for this
approach in Section 5.

##  2 Theory

###  2.1 The Need for Alternative Ground-Truth Distance Measures

As discussed in Section 1, Euclidean distance (measured in Cartesian
coordinates)111Hereinafter, unless specified otherwise, the term Euclidean
distance is assumed to mean that it is measured in Cartesian coordinates. as
well as other, more intricate similarity measures have been devised and
utilized to interpret the relative arrangement of distinct atomistic
configurations in an energy landscape in order to enable more accurate and
generalizable predictive models. However, given that we are attempting to
model the energy, which is a unique manifestation of the underlying physics of
the atoms and electrons in the system, it seems limiting to rely only on
spatially-informed similarity/distance measures without any concrete (even
approximate) energetic information. This is especially significant when
attempting to sample atomic configurations that feature sufficient structural
and energetic diversity for potential models to be able to model both low- and
high-energy configurations accurately.

Consider the heuristic example of the Morse potential that defines the
interaction between two particles separated by a distance $r$ as

$V(r)=D_{e}\left(1-e^{-a(r-r_{e})}\right)^{2}-D_{e},$ (1)

where the parameter $D_{e}$ governs the range of the interaction, $r_{e}$
defines the equilibrium length, and $a$ determines the decay rate of the
interaction. Fig. 1 shows the Morse interaction potential for two values of
$D_{e}$.

Figure 1: The Morse potential,
$V(r)=D_{e}\left(1-e^{-a\left(r-r_{e}\right)}\right)^{2}-D_{e}$, for well-
depth parameter $D_{e}=1,2$ for a simple two-atom model system. $r_{e}$ is the
equilibrium distance between the two atoms. Despite the spatial similarity
between structures at $r=1.0$ and $r=2.0$ being independent of the parameters
of the potential, the geodesic distance between them changes as a consequence
of changing the interaction, i.e. $L_{1}\rightarrow L_{2}$.

As can be seen from the figure, two configurations, e.g., at $r=1.0$ and
$r=2.0$, will be separated by the same Euclidean distance - or indeed any
spatially-informed distance measure - regardless of the potential that
operates on the structures. However, the geodesic distance along the potential
energy manifold between the two configurations is a unique consequence of the
particular physical interaction (in this case, Morse potentials with
$D_{e}=1,2$) that operates on the structures. Referring to the “Model =
Representation + Architecture” mnemonic in Section 1, it can be seen that
using a purely spatial, potential-invariant representation in a model that may
be required to capture distinct types of physical interactions - e.g.,
attractive versus repulsive interactions - shifts much of the computational
burden of adapting to the energy landscape on to the architecture, since the
representation is not adapted to the physics of the system. (Note that the
short-/long-range interaction example is somewhat simplistic given the variety
of solutions that have been proposed to address this problem. See, for
instance, Ref. 4.)

A principal argument in this work is that the geodesic distance measure has
the capacity to efficiently encode information about the potential energy
manifold that is absent in conventional spatial similarity measures, and that
this information can be embedded into well-known descriptor schemes via
appropriate parameterization. This argument is based on the following
observations:

  1. 1.

Geodesic distances are information-rich: Given data about a discrete
collection of structures and their energies (or forces), it is the task of
interatomic potential models to make the best possible estimations of the
energies of the continuum of structures that lie between - or beyond - those
structures within a certain representation (descriptor) scheme. For this
reason, one can also expect that any a priori structure$\rightarrow$energy
data about the continuum of states between these discrete structures is likely
to benefit from a statistically optimized potential energy model. While this
information is not directly available in any limited, discrete set of training
data, knowledge of geodesic distances (between a large number of pairs of
structures) necessarily contains this information implicitly. Consider, for
instance, two cases of the geodesic distance, $l_{AB}$, between structures $A$
and $B$, in relation to the direct, Euclidean distance, $l^{0}_{AB}$, between
the structures, calculated as the length of the linear path in Cartesian
coordinates between structures $A$ and $B$. Case 1: If $l_{AB}\approx
l^{0}_{AB}$, this indicates the absence of a significant energy barrier or
basin between the structures and, by extension, that any structures lying
directly between them (in Cartesian coordinates) are unlikely to have very
high potential energy compared to $A$ and $B$. Case 2: Conversely, if
$l_{AB}\gg l^{0}_{AB}$, there are likely several high- or low-energy (compared
to A and B) structures along the direct path between $A$ and $B$, due to which
the minimum-energy path takes a more non-linear (thereby, longer) route
between the structures. (Of course, $l_{AB}$ and $l^{0}_{AB}$ are intimately
connected to the coordinate system that we utilize. We address this connection
in more detail in Section 2.2.) If one could devise a way to capture this
information directly within the descriptor scheme used to represent an
atomistic system, this would - in principle - enable a fitting algorithm to
access the rich geometric structure of the energy manifold, even if this is
only approximately known.

  2. 2.

Fast approximations exist for computing geodesic distances: At first glance,
it seems causally incorrect to expect that geodesic distances along the
potential manifold will be available for parameterizing descriptors when the
very purpose of the descriptors is to be able to approximate the potential
that one does not know a priori via some (learnable) function. However, our
goal is not to compute the exact geodesic distances between structures using
the (theoretically) true potential energy function, but only to obtain
reasonable approximations of geodesic paths that capture the general geometric
features of the energy landscape. As we will discuss in Section 2.2, borrowing
ideas from chemical reaction theory, it is possible to rapidly generate
approximate geodesic paths and distances between a large number of pairs of
structures at only a fraction of the cost of computing them with full ab
initio accuracy and using a pathway-search algorithm such as nudged elastic
band.

  3. 3.

Descriptor schemes contain several optimizable hyperparameters: Popular
schemes for generating atomic fingerprints often contain several
hyperparameters that determine how an atomistic system is featurized. Beyond a
fairly limited set of specifications regarding the range and scale of these
values, they may be largely considered a set of optimizable hyperparameters
that can be adjusted to best capture the relevant structural features of the
system. Furthermore, descriptor schemes can utilize multiple copies of the
same hyperparameter with different values to increase the structural
“completeness” or uniqueness of the atomic fingerprint, introducing additional
degrees of freedom that may be optimized or fine-tuned. Hyperparameter tuning
generally follows a train/validation pipeline, wherein a subset of the overall
dataset - termed the validation set - is used to evaluate, and possibly
improve, a model’s hyperparameters without being directly used to train the
model. Instead, our thesis here is that the hyperparameters of a descriptor
scheme can be optimized to maximize the descriptor’s informativeness about the
energy landscape without making reference to any particular model
architecture, purely through the use of information from pre-computable
distance measures, such as geodesic distances. We discuss this in greater
detail in the sections that follow, and demonstrate it in the context of ACSF
descriptors.

It is important to clarify here that geodesic distances are not necessarily
always “better” or “more informative” than purely spatially-informed distance
measures, such as Euclidean distance. Rather, our aim is to illustrate that,
because geodesic distance incorporates a fundamentally different category of
physical information - namely, energetic - than spatial distance, it is likely
to contain at least some information that is beyond the capacity of spatial
similarity measures to express. The goal then becomes twofold: First, to
estimate the relative information content of two distance measures, e.g.,
Euclidean and geodesic. Second, to design a descriptor so as to maximize its
information content given that different ground-truth distance measures may
contain partially exclusive physical information.

We now discuss a potential strategy for the fast approximation of geodesic
distances and its relevance for the design of more informative descriptors.

###  2.2 Fast Approximation of Geodesic Distances and Application to
Descriptor Design

In the context of determining the minimum-energy (geodesic) pathways between
equilibrium structures in reaction chemistry, it is of interest to generate
initial guess pathways that are as accurate and close to the true minimum-
energy pathways as possible. While early attempts at generating these guess
pathways relied on naive approximations such as linear interpolation between
the Cartesian coordinates of the initial and final structures,18 more advanced
interpolation techniques have since been developed.

By assessing various interpolation schemes using a unifying differential-
geometric approach, Zhu et al. (Ref. 17) noted that the relative success of
certain interpolation schemes - particularly those involving inverse bond
length-based internal coordinate representations - over others, can be
explained by the implicit choice of metric tensor (mapping length elements of
the Riemannian manifold of the potential energy to the coordinate system) that
accompanies these representations. More intuitively, they highlight that
linear interpolation between two points (atomic configurations) in some system
of coordinates is likely to resemble the geodesic path along the energy
manifold if and only if the potential energy surface is “almost flat” when
represented in those coordinates. Mathematically, given a potential energy
function, $U$, and mass-weighted Cartesian coordinates, $\mathbf{x}$, the
construction of an internal coordinate representation $\mathbf{q}$ such that

$\frac{\partial q^{k}}{\partial x_{i}}\frac{\partial q^{k}}{\partial
x_{j}}\approx\frac{\partial U}{\partial x_{i}}\frac{\partial U}{\partial
x_{j}}$ (2)

will necessarily result in straight-line paths in $\mathbf{q}$ approximating
geodesic paths along the potential energy manifold. (Note that the Einstein
summation convention is assumed in Eq. 2.)

An alternative way to state this would be to note that, the more closely an
internal coordinate representation reflects the true geometry of the potential
energy surface, the more accurately its Euclidean paths will correspond to
geodesic paths along the manifold. They showed that it is possible to achieve
this approximation using a simple system of internal coordinates given by

$q_{i}\equiv\exp\left[-\alpha\frac{r_{kl}-r_{e}^{kl}}{r_{e}^{kl}}\right]+\beta\frac{r_{e}^{kl}}{r_{kl}},$
(3)

where $r_{kl}$ is the distance between atoms $k$ and $l$, and $r_{e}^{kl}$ is
a parameter representing the equilibrium distance between $k$ and $l$. By
adjusting the parameters $\alpha$ and $\beta$, it is possible to control the
extent to which short- or long-range interactions operate between pairs of
atoms given their separation; it is then straightforward to modify this scheme
to include different parameters for chemically distinct pairs of atoms. Using
this approach, they demonstrated that good approximate geodesic pathways
between vastly different configurations may be constructed rapidly, even for
biomolecular systems involving $\mathcal{O}(10^{3})$ atoms. Furthermore,
numerical approximations enable the geodesic distances between two
configurations to be computed to arbitrary precision, which we do not discuss
here for brevity.

In the context of fitting physics-inspired potential energy forms to ab initio
data, descriptors are of particular significance due to the flexibility with
which they can be designed to capture and emphasize various physical and
chemical phenomena that one anticipates in a given class of atomistic systems.
Descriptors are generally constructed in a way that captures the translation-,
rotation-, and identical-atom permutation-invariance of the potential energy.
As a consequence, many descriptors naturally fall under the definition of the
internal coordinates referenced in the foregoing discussion.

This raises the question of why linear interpolation in the space of
descriptors does not immediately result in geodesic interpolation along the
manifold given that descriptors are designed to capture the features of
potential energy surfaces. First, descriptors are not always constructed to
satisfy the equivalence of their metric tensor with the Fukui metric for mass-
weighted Cartesian coordinates, as depicted in Eq. 2. If they were, this would
largely obviate the need for highly non-linear fitting functions, such as
neural networks, to model the energy landscape. Second, even if a potential
manifold were overwhelmingly flat in the coordinate space of some descriptor,
the repeated evaluation of analytical or numerical derivatives of these
descriptors with respect to Cartesian coordinates - a requirement for the
procedure of Ref. 17 \- can be a highly non-trivial computational task that
may be hard to generalize to new descriptor forms, except those explicitly
designed to be used with gradient-based optimization algorithms.

Instead of attempting to construct descriptors that achieve the metric
equivalence of Eq. 2, we take an alternative approach. As described below, we
use a data-driven procedure to encode the information content of a large
collection of easily computable pair-wise geodesic distances into the
hyperparameters of descriptor forms without making any reference to the metric
that governs their correspondence to the energy manifold. This allows us to
retain the flexibility afforded by various descriptor schemes while utilizing
additional information about the geometry of the potential energy landscape
that can be determined prior to fitting the parameters of a statistical
learning model to the training data, or indeed without even having defined
such a model.

By using pre-computable data on the approximate geodesic distances between
structures along the energy manifold we are able to endow the descriptor
scheme with a more meaningful notion of how atomic configurations should be
distributed in descriptor space. In making the underlying assumption that
geodesic distances between structures - a quantity defined on the
$(3N-6)$-dimensional energy manifold of an $N$-atom system - can by captured
within the $\mathcal{O}$(1–10)-hyperparameter set of a descriptor, we are
effectively invoking the manifold hypothesis19 before even passing the atomic
fingerprints computed using the descriptor to a regression algorithm that
attempts to fit its parameters to the energies. We note here that our
descriptor optimization approach is not restricted to one particular method of
computing geodesic distances, such as in Ref. 17. It can be modified to
accommodate other distance measures, enabling different types of information
to be expressed.

We demonstrate this procedure for atom-centered symmetry function (ACSF)
descriptors (Section 3.1), which we have chosen due to their scalable local
construction, physical interpretability in featurizing atomic environments,
and widespread use.9 We encode the information content within the descriptor
hyperparameters using a Bayesian optimization algorithm (Section 2.4) whose
objective function is set to be the information imbalance between descriptor
distances and a ground-truth distance measure, which we set to be the
Euclidean and geodesic measures for comparison (Section 2.3).

###  2.3 Information Imbalance as an Objective Function

We now address the question of the choice of objective function to use for
capturing the information content of a distance measure, e.g., the geodesic
measure, in the descriptor’s hyperparameters. First, we note that it is
impractical to optimize the Euclidean distance between descriptor fingerprints
to be identical to the geodesic distances between the structures themselves
because the latter are unitless scalar quantities and there is no physical
motivation for equalizing them. Instead, we propose using a statistical,
information-theoretic quantity termed the symmetric information imbalance
between two distance measures.15 This may be described briefly as follows.
Consider two distance measures, $\mathcal{A}$ and $\mathcal{B}$. Given a
collection of points distributed according to some unknown distribution, if we
select any single point as a reference and rank each of the other points
according to their distance from the reference point, it is possible that the
rank ordering will not be identical between the two distance measures. This is
equivalent to the observation that distances between points in two distinct
distance measures need not vary monotonically with each other. The information
imbalance in going from $\mathcal{A}$ to $\mathcal{B}$ is defined as

$\Delta(\mathcal{A}\rightarrow\mathcal{B})=\frac{2\langle
r_{ij}^{\mathcal{B}}|r_{ij}^{\mathcal{A}}=1\rangle_{ij}}{N},$ (4)

where $r_{ij}^{\mathcal{X}}$ is the rank of point $j$ with respect to point
$i$ in distance measure $\mathcal{X}$, and $N$ is the total number of points
in the dataset. In essence, $\Delta(\mathcal{A}\rightarrow\mathcal{B})$
answers the question, on average, for two points $i$ and $j$ that are nearest
neighbors (rank $r_{ij}^{A}=1$) in distance measure $A$, what is their rank in
distance measure B ($r_{ij}^{B}$)? The symmetric information imbalance is then
defined as

$\bar{\Delta}(\mathcal{A}\leftrightarrow\mathcal{B})=\frac{\Delta(\mathcal{A}\rightarrow\mathcal{B})+\Delta(\mathcal{B}\rightarrow\mathcal{A})}{\sqrt{2}},$
(5)

where $\bar{\Delta}$ simply symmetrizes the imbalance in going from one
distance measure to the other. In general, if
$\Delta(\mathcal{A}\rightarrow\mathcal{B})<\Delta(\mathcal{B}\rightarrow\mathcal{A})$,
distance measure $\mathcal{A}$ may be said to be more informative than
$\mathcal{B}$, while a low value of
$\bar{\Delta}(\mathcal{A}\leftrightarrow\mathcal{B})$ indicates greater
agreement between $\mathcal{A}$ and $\mathcal{B}$ than a higher (closer to 1)
value.

Information imbalance has previously been applied to guide descriptor
compression and has been used to assess the correspondence between distances
in various descriptor representations and Euclidean distances in a dataset of
Si atoms slightly perturbed from their bulk lattice positions.15
Interestingly, in Ref. 15, the Euclidean distance is taken to be the ground-
truth distance measure between structures by showing that increasingly
expressive descriptors show lower information imbalance with respect to the
Euclidean measure. As we will illustrate in Section 4.1, however, the
Euclidean and geodesic distance measures can fundamentally disagree. This,
combined with the discussion in Section 2.1, suggests that there is no
underlying reason for one of these contradictory distance measures to be
chosen as the ground-truth over the other.

###  2.4 Bayesian Optimization

Having defined an objective function, we formulate the optimization problem as
follows. For a set of structures, we pre-compute all pair-wise distances,
$\mathbf{r}^{\mathcal{T}}$, in the ground-truth (Euclidean or geodesic)
distance measure, $\mathcal{T}$. Then, for a given set of descriptor
hyperparameters $\mathbf{p}$, we evaluate the atomic fingerprints of each
configuration in the structure set and compute the pair-wise distances,
$\mathbf{r}^{\mathcal{D}(\mathbf{p})}$ in the space of this descriptor,
$\mathcal{D}(\mathbf{p})$. We then evaluate the objective function, i.e., the
symmetric information imbalance between the descriptor and the ground-truth
measure given by

$\bar{\Delta}(\mathcal{D}(\mathbf{p})\leftrightarrow\mathcal{T})=\frac{\sqrt{2}}{N}\left(\langle
r^{\mathcal{D}(\mathbf{p})}_{ij}|r^{\mathcal{T}}_{ij}=1\rangle_{ij}+\langle
r^{\mathcal{T}}_{ij}|r^{\mathcal{D}(\mathbf{p})}_{ij}=1\rangle_{ij}\right).$
(6)

The task of the optimization algorithm then becomes to iteratively update the
set of hyperparameters $\mathbf{p}$ in order to identify the set that
minimizes the information imbalance of the descriptor distance measure with
respect to the ground-truth measure.

The question now arises of which optimization algorithm to use. Since the
objective function is a statistical quantity, it is difficult to rely on
gradient-based algorithms, such as gradient descent. Moreover, since new
hyperparameter sets must be generated and passed to Eq. 6 on-the-fly, it is
desirable to use an algorithm that operates in data-sparse regimes and
intelligently queries new data points for inexpensive, data-efficient
optimization. Computing the information imbalance from the pair-wise distances
of a large set of structures is a moderately expensive computational
operation, which further encourages the use of such an efficient algorithm. As
we discuss below, Bayesian optimization presents an elegant solution to these
problems.

Bayesian optimization has emerged as a powerful strategy for the gradient-
free, stochastic optimization of black-box function parameters. Its data
efficiency, coupled with the ability to systematically explore the landscape
of the objective function through a probabilistic framework, makes it
particularly appealing for tasks in which function evaluations are
computationally expensive. This optimization technique leverages a
probabilistic model - typically a Gaussian process - to model the objective
function and utilizes an acquisition function to guide the search for the
optimum by balancing exploration and exploitation.20 This enables the
algorithm to incorporate prior knowledge and iteratively update it with new
observations, thereby refining its model of the objective function with each
iteration. Such an approach is in stark contrast to grid or random search
methods, which lack the ability to learn from prior evaluations and can be
prohibitively expensive in high-dimensional spaces or when evaluations are
costly.21

In mainstream machine learning, Bayesian optimization has been successfully
applied to the task of hyperparameter tuning. Traditional methods for
hyperparameter optimization, such as exhaustive search, are often impractical
due to the exponential growth of the search space with the number of
hyperparameters. Bayesian optimization addresses this challenge by
constructing a probabilistic model of the performance function and iteratively
selecting the most promising hyperparameters to evaluate based on the current
model. This approach has been shown to significantly reduce the number of
evaluations needed to identify optimal or near-optimal hyperparameters.22

In the present context, we use Bayesian optimization to intelligently query
new descriptor hyperparameter sets and compute the information imbalance
between that descriptor scheme and the ground-truth distance measure, e.g.,
the Euclidean or geodesic measure.

##  3 Computational Details

###  3.1 ACSF Descriptors

Atom-centered symmetry function (ACSF) descriptors rely on the construction of
functions centered around each atom ($i$) of a system by featurizing only the
local neighborhood of that atom.10, 11 Here, we begin by focusing on
$G^{\textrm{II}}$\- and $G^{\textrm{IV}}$-type ACSF descriptors, which take
the form

$G_{i}^{\textrm{II}}=\sum_{j\neq
i}^{R_{ij}<R_{c}}e^{-\eta\left(R_{ij}-R_{s}\right)^{2}/R_{c}^{2}}f_{c}\left(R_{ij}\right),$
(7)

and

$G_{i}^{\textrm{IV}}=2^{1-\zeta}\sum_{j,k\neq
i}^{R_{ij},R_{ik}<R_{c}}\left(1+\lambda\cos\theta_{ijk}\right)^{\zeta}e^{-\eta\left(R_{ij}^{2}+R_{jk}^{2}+R_{ki}^{2}\right)}f_{c}\left(R_{ij}\right)f_{c}\left(R_{jk}\right)f_{c}\left(R_{ki}\right).$
(8)

In Eqs. 7 and 8, $R_{ij}$ is the separation between atoms $i$ and $j$, $R_{c}$
is a cutoff radius parameter, $R_{s}$ is a small displacement to the center of
the spherical shell around which the descriptor is constructed, $\eta$,
$\zeta$, and $\lambda$ are tunable parameters, and $f_{c}(\cdot)$ is a smooth
cosine-type tapering function that goes to 0 at $R_{c}$.

We now describe the three schemes of ACSF descriptors that we optimize using
our approach.

  1. 1.

Scheme I: $G^{\textrm{II}}$-type descriptors with 2 $\eta$ parameters,

  2. 2.

Scheme II: $G^{\textrm{II}}$-type descriptors with 4 $\eta$ parameters, and

  3. 3.

Scheme III: $G^{\textrm{II}}$\- and $G^{\textrm{IV}}$-type descriptors with 4
$\eta$ parameters in the $G^{\textrm{II}}$ descriptor, 1 (different) $\eta$
parameter in the $G^{\textrm{IV}}$ descriptor, 1 $\zeta$ parameter in the
$G^{\textrm{IV}}$ desciptor, and 2 $\lambda$ parameters in the
$G^{\textrm{IV}}$ descriptor.

###  3.2 Optimization Settings

The hyperparameters of each descriptor defined in Section 3.1 are then treated
as the variables that we attempt to optimize with the Bayesian optimization
algorithm. We perform all Bayesian optimization calculations using the scikit-
optimize package.23 For each loop of the algorithm, we use 100 function calls
and 50 random initializations, which we found to be sufficient for minimizing
the information imbalance given the dimensionality of the hyperparameter
space. The information imbalance with respect to a given ground-truth distance
measure - whether Euclidean or geodesic - was computed in each case using in-
house Python implementations of Eqs. 4–6.15 At the end of an optimization
cycle, the hyperparameter set that achieves the lowest value of the objective
function is chosen as the optimized set. If more than one hyperparameter set
achieves an optimal value of the objective function, then one of these sets is
chosen at random by the algorithm, which we report directly.

Further details on the constrained optimization of the hyperparameters are
provided in Section 4.2.

###  3.3 Dataset and Ground-Truth Distance Computation

We perform our analysis using molecular configurations of the ethanol,
malonaldehyde, and aspirin datasets of the MD22 database.24 These datasets
include forces and energies computed using coupled-cluster theory with single,
double, and perturbative triple excitations (CCSD(T)) on atomic configurations
taken from ab initio molecular dynamics (AIMD) trajectories. These AIMD
simulations themselves utilize DFT-PBE calculations with van der Waals
corrections and are performed at 400–500 K.24 The use of high-temperature MD
simulations is significant because it indicates a higher likelihood of
sufficiently broad sampling of the energy landscape including the traversal
and crossing of potential energy barriers. We choose these molecules as test
cases due to the conformational flexibility observed in the molecular dynamics
trajectories of these molecules in the MD22 database. Furthermore, the correct
measurement of Euclidean distances and geodesic distances between different
configurations first requires them to be maximally aligned, both
translationally and rotationally, prior to evaluating the displacement of each
of the individual atoms and the geodesic pathways between configurations. The
relatively small size of these molecules makes such alignment fairly
unambiguous - in this case, by first aligning their center-of-mass coordinates
and then using the Kabsch algorithm25.

Euclidean distances are computed in the usual way, i.e., as the $L^{2}$ norm
of the vector difference between the Cartesian representation of two
structures. Geodesic distances are computed as in Ref. 17 with the values of
the parameters $\alpha$ and $\beta$ in Eq. 3 being set to 1.7 and 0.01,
respectively, in all cases. We note, however, that it might even be useful to
calibrate these parameters using energies from the training dataset prior to
computing geodesic distances. We compute Euclidean and geodesic distances
between all pairs of structures from a subset of 100 structures from the
dataset of each molecule, chosen using a farthest-point sampling algorithm.12

##  4 Results and Discussion

###  4.1 Inequivalence of Euclidean and Geodesic Distance Measures

First, we provide a more concrete demonstration of the inequivalence of the
Euclidean and geodesic distances via the example of pair-wise distances in a
100-structure subset of the ethanol training dataset of the MD22 database. In
Fig. 2(a), we clearly observe that the geodesic distances vary highly non-
monotonically with the Euclidean distance, with only a slight positive
correlation emerging when many - i.e., ${100\choose 2}$ in this case - pairs
are considered. It can also be observed that there is a greater spread of the
Euclidean distances than the geodesic distances, suggesting that pairs of
configurations may seem more different from each other based on the Euclidean
distance despite being relatively similar based on the geodesic distance. This
shows that the Euclidean distance alone can be somewhat misleading as a metric
for the separation between configurations if one wants to optimally sample
configurations along a potential energy landscape. Now, plotting the ranks
(indexed closeness) of every point with respect to every other point in the
data subset using the geodesic distance measure versus the Euclidean distance
measure yields Fig. 2(b). Given that distance measures in good agreement
should display an almost linear relationship in such a parity plot, the highly
noisy nature of Fig. 2(b) indicates generally poor agreement between the
Euclidean and geodesic distance measures. Equivalent plots for malonaldehyde
and aspirin are shown in Fig.
LABEL:fig:malonaldehyde_and_aspirin_demonstration of the Supporting
Information.

((a))

((b))

Figure 2: Non-equivalence of Euclidean and geodesic distance measures
demonstrated on all pairwise distances in a set of 100 structures from the
MD22 ethanol training dataset. (a) Non-monotonicity between Euclidean and
geodesic distances. Despite the general collective increase of geodesic
distance with increasing Euclidean distance as seen by the positive slope of
the black regression line, the relationship is much more ambiguous at the
level of individual configurations. (b) Scatter plot of ranks between pairs of
points in both distance measures.

Next, we demonstrate this more quantitatively by plotting the information
imbalance in going from the geodesic distance measure ($\mathcal{G}$) to the
Euclidean distance measure ($\mathcal{R}$), i.e.,
$\Delta(\mathcal{G}\rightarrow\mathcal{R})$, versus the reverse,
$\Delta(\mathcal{R}\rightarrow\mathcal{G})$. This is depicted in Fig. 3 for a
100-structure subset of the ethanol, malonaldehyde, and aspirin training
datasets of the MD22 database. Recall that, in this representation, points
above the $y=x$ line indicate that the Euclidean distance measure is “more
informative” than the geodesic distance measure and vice versa, while a value
that does not lie close to the origin indicates a significant, finite
disagreement between the two distance measures.

Figure 3: Information imbalance between Euclidean and geodesic distances
measured on all pairwise distances in a set of 100 structures from the MD22
training datasets of ethanol, malonaldehyde, and aspirin.
$\Delta(\mathcal{R}\rightarrow\mathcal{G})$ refers to the information
imbalance in going from the Euclidean measure ($\mathcal{R}$) to the geodesic
measure ($\mathcal{G}$), and vice versa. The dashed line ($y=x$) broadly
divides the plot into regions of relative informativeness. For example,
regions above the dashed line indicate that $\mathcal{R}$ contains the
information in $\mathcal{G}$ and vice versa.

It can be seen from Fig. 3 that the Euclidean and geodesic distance measures
are not necessarily in good agreement. Conceptually, Fig. 3 may be interpreted
as follows. If a point lies above (below) the dashed line $y=x$, then any two
atomic configurations in the dataset that are separated by a small Euclidean
(geodesic) distance are also statistically likely to also be separated by a
small geodesic (Euclidean) distance. However, the converse is less true, i.e.,
two configurations separated by a small geodesic (Euclidean) distance are
statistically less likely to also be separated by a small Euclidean (geodesic)
distance. An alternative way to state this would be to say that low Euclidean
(geodesic) distance is a better indicator of low geodesic (Euclidean) distance
than the reverse. In more concrete, physical terms, this interpretation may be
stated in the context of potential energy manifolds following the discussion
in Section 2.1. Consider two atomic configurations $A$ and $B$ separated by a
moderate geodesic distance. Of several possibilities for the geodesic path
connecting $A$ and $B$, consider the following two cases: (1) $A$ and $B$ are
spatially distant but are separated by a smooth, flat region on the potential
energy surface due to which the geodesic distance is directly correlated with
the linear (Euclidean) path between the structures; and (2) $A$ and $B$ are
spatially close but are separated by significant energy barriers or basins due
to which the direct Euclidean path does not accurately reflect the geodesic
path. This entails a possible non-uniqueness in the mapping from geodesic
distance to Euclidean distance and vice-versa. Consider the case where the
mapping from geodesic distance to Euclidean distance is less unique (more
degenerate) than the inverse map, i.e.,
$\Delta(\mathcal{R}\rightarrow\mathcal{G})<\Delta(\mathcal{G}\rightarrow\mathcal{R})$.
This suggests that a small spatial (Euclidean) separation between
configurations $A$ and $B$ is a better indicator of geodesic closeness than a
small geodesic distance between $A$ and $B$ is an indicator of Euclidean
closeness.

Of course, as seen in Fig. 3, the extent to which this (or its opposite) is
true depends on the chemistry of the specific system being studied. For
example, geodesic distances appear more informative for the aspirin sub-
dataset, while Euclidean distances appear more informative (to varying
degrees) for the ethanol and malonaldehyde sub-datasets. We note that
providing mechanistic explanations for the observed values of information
imbalance is a non-trivial task. This is because information imbalance is a
statistical quantity computed over many pairs of structures that may display a
wide range of behaviors in the relationship between Euclidean and geodesic
distances. Instead, our goal with these examples is primarily to show that
Euclidean and geodesic distances between atomic configurations tend to be
inequivalent, as evidenced by the non-uniqueness in the mapping from one
distance measure to the other.

We now discuss the optimization of the ACSF descriptors with the objective of
minimizing the symmetric information imbalance between distances in descriptor
space and one or both of the Euclidean and geodesic distance measures.

###  4.2 Combined-Objective Descriptor Optimization

The descriptor schemes defined in Section 3.1 contain few to several
hyperparameters which can be treated as the optimizable parameters for the
Bayesian optimization framework. However, indiscriminately optimizing them not
only increases the burden on the optimizer, it may also wastefully iterate
through many meaningless or unphysical functions. To avoid this, we optimize
the (hyper)parameters of the descriptors with the following constraints:

  1. 1.

Scheme I: First, we optimize $R_{c}$ and $R_{s}$ restricting them to the
continuous range of $(1,7)$ Å. The two $\eta$ parameters, $\eta_{1}$ and
$\eta_{2}$, are optimized separately. We first optimize two auxiliary
parameters, $x_{1}\in[\epsilon,1-\epsilon]$ and $x_{2}\in[1+\epsilon,20]$,
where $\epsilon$ is some small value $\sim 10^{-6}$, and subsequently set
$\eta_{1}=\log_{10}(x_{1})$ and $\eta_{2}=\log_{10}(x_{2})$.

  2. 2.

Scheme II: Similar to Scheme I, in addition to $R_{c}$ and $R_{s}$, two
parameters $x_{1}$ and $x_{2}$ are first optimized, following which four
$\eta$ parameters are defined to lie evenly spaced on the $\log_{10}$ scale,
with the endpoints being set to $\eta_{1}=\log_{10}(x_{1})$ and
$\eta_{4}=\log_{10}(x_{4})$.

  3. 3.

Scheme III: In addition to $R_{c}$, $R_{s}$, and the four $\eta$ parameters
optimized according to Scheme II, we optimize a single $\eta$ parameter for
the $G^{\textrm{IV}}$ descriptors. Further, we optimize $\zeta\in[2,64]$,
requiring that it is an integer and explored on the log scale. $\lambda$ is
also set to $-1$ and $1$ alternately.

Using each of these schemes, we now utilize the Bayesian optimization
algorithm described in Section 2.4 to optimize the descriptor hyperparameters
so as to minimize the symmetric information imbalance between the Euclidean
distance measure in descriptor space, $\mathcal{D}$, and either the Euclidean
distance measure (in the original Cartesian coordinates), $\mathcal{R}$, or
the geodesic distance measure, $\mathcal{G}$. This results in two different
objective functions, $\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})$ and
$\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})$, with which the
descriptor can be designed. However, to demonstrate that it is possible for a
single descriptor to reconcile the unshared information contained in both the
Euclidean and geodesic distance measures, we also optimize the descriptors to
minimize a third objective function, which we define to be the sum of the two
previous information imbalances, i.e.,
$\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})+\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})$.

In Fig. 4, we show the results of Bayesian optimization of each of these
objectives as plots of $\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})$
versus $\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})$ for each of three
systems - ethanol, malonaldehyde, and aspirin - for descriptor Scheme 3.
Similar plots for Schemes 1 and 2 are shown in Figs.
LABEL:fig:supp_bayesian_optimization_1 and
LABEL:fig:supp_bayesian_optimization_2 of the Supporting Information. We note
that, in each of these plots, the closer a point (i.e., an optimized
descriptor) lies to the origin the better, with the best descriptors lining
the Pareto frontier of the plot. In Bayesian optimization, if multiple
explored points optimize the objective equally well, they are considered
equivalent and any single one of them may be chosen and reported by the
algorithm as the optimal choice. (We do, in fact, observe this phenomenon
frequently when the optimizer is tasked with minimizing either
$\bar{\Delta}(\mathcal{D}\rightarrow\mathcal{G})$ or
$\bar{\Delta}(\mathcal{D}\rightarrow\mathcal{R})$.) Fig. 4 shows that -
notably in the case of ethanol and malonaldehyde - given the objective of
minimizing only $\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})$,
Bayesian optimization identifies descriptors that capture the information in
the Euclidean distance measure reasonably well, albeit without much shared
information with the geodesic measure. The reverse is also found to be true
when optimizing only $\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})$.
However, the use of the combined objective function,
$\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})+\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})$,
generates descriptors that optimally balance the information content contained
in both the Euclidean and geodesic distance measures, as seen by the bold blue
cross points, which are consistently the closest to the origin. Comparing
Figs. 4, LABEL:fig:supp_bayesian_optimization_1, and
LABEL:fig:supp_bayesian_optimization_2, it is also interesting to note that
this theme is observed regardless of the number of hyperparameters in the
descriptor scheme.

((a))

((b))

((c))

Figure 4: Bayesian optimization samples of descriptor parameters (Scheme 3)
for (a) ethanol, (b) malonaldehyde, and (c) aspirin. The objective function is
chosen to minimize the symmetric information imbalance between the descriptor
distance measure ($\mathcal{D}$) and the geodesic distance ($\mathcal{G}$)
[i.e. $\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})$, green circles],
the Euclidean distance ($\mathcal{R}$) [i.e.
$\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})$, red triangles], and the
sum of both imbalances [i.e.
$\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})+\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})$,
blue squares]. Crosses represent the optimal descriptor parameters suggested
at the end of each Bayesian optimization loop. Combined optimization leads to
descriptors closest to the origin (blue crosses), i.e., encoding information
about both ground-truth distance measures. Black dashed lines mark the
theoretical limit for 100 samples.

Given these newly constructed descriptors, it is now desirable to check
whether descriptors constructed by minimizing one objective, say
$\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})$, are more, less, or
equally informative compared to descriptors constructed by minimizing the
other objectives. We present this in Fig. 5. In these heatmaps, $D^{c}$,
$D^{\mathcal{R}}$, and $D^{\mathcal{G}}$ refer to descriptors that are
optimized by minimizing (a) the combination of Euclidean and geodesic
information imbalances,
$\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})+\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})$;
(b) only the Euclidean imbalance,
$\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{R})$; and (c) only the
geodesic imbalance, $\bar{\Delta}(\mathcal{D}\leftrightarrow\mathcal{G})$. The
entries of each matrix in Fig. 5 are arranged as follows:
$(i,j)=\Delta(D^{i}\rightarrow D^{j})$, where we now note the asymmetric
information imbalance.

Figure 5: Information imbalance between three types of descriptors.
$D^{\mathcal{R}}$ is optimized to minimize the symmetric information imbalance
($\bar{\Delta}$) between descriptor distances and Euclidean distances.
$D^{\mathcal{G}}$ is optimized to minimize $\bar{\Delta}$ between descriptor
distances and geodesic distances. $D^{c}$ is optimized to minimize the sum of
symmetric information imbalances between descriptor distances and both
Euclidean and geodesic distances. In each matrix, entry
$(i,j)=\Delta(D^{i}\rightarrow D^{j})$. The confinement of large imbalance
values ($>$ 0.5) to the $D^{\mathcal{R}}\times D^{\mathcal{G}}$ sub-matrices
indicates that the joint Euclidean–geodesic optimization of descriptors leads
them to contain more comprehensive distance information.

In the case of ethanol and malonaldehyde, we note that high ($>0.5$) values of
information imbalance - in some cases, significantly high, $\sim 1$ \- are
restricted to the $D^{\mathcal{R}}\times D^{\mathcal{G}}$ sub-block of each of
the $3\times 3$ matrices. We do not include the results for aspirin since all
descriptors are found to be in complete agreement with each other in this
case. In the case of ethanol, the highest of these correspond to the imbalance
term $\Delta(D^{\mathcal{G}}\rightarrow D^{\mathcal{R}})$, whereas for
malonaldehyde, these correspond to $\Delta(D^{\mathcal{R}}\rightarrow
D^{\mathcal{G}})$. This highlights the unintuitive fact that substantial
disagreements in distance information between two similarly constructed
descriptors can arise from the seemingly inconsequential fact that their
hyperparameters better encode one type of ground-truth distance information
over the other. Conversely, optimizing the descriptors to minimize the
combined objective leads to better agreement with descriptors constructed
using only one objective, indicating that combined-objective optimization is
capable of resolving the disagreement between contradictory structural
representations.

##  5 Conclusions and Outlook

In this work, we have proposed a novel paradigm for optimizing and evaluating
atomistic descriptors by taking into account the partially exclusive
information content of Euclidean and geodesic distance measures in molecular
configuration space. We have presented arguments in favor of utilizing updated
ground-truth distance measures that account for physical information beyond
spatial similarity, demonstrating this via the fundamental disagreements
between Euclidean and geodesic distances. We have then introduced the
relevance of this approach in atomistic descriptor design through a Bayesian
optimization algorithm that utilizes information imbalance as an objective
function, thereby showing that it is possible to construct descriptors that
maximize their shared information with a given distance measure, either
Euclidean or geodesic. We have also shown that disagreements between
descriptors that encode distances in only one distance measure can be resolved
by a combined-objective Bayesian optimization strategy. This is noteworthy
because, by minimizing the information imbalance between a descriptor scheme
and conflicting distance measures, we maximize the uniqueness of atomic
fingerprints generated using that descriptor. This, in turn, makes them more
amenable to unambiguous sampling - structural or energetic - and statistical
modeling.26, 27

Atomistic representations in which direct (linear) distances correspond to
geodesic distances along potential energy manifolds may have the capacity to
disentangle the often highly complex, nonlinear structure of these manifolds,
rendering them tractable to easier interpretation, partitioning (into similar
and dissimilar regions), and fitting by statistical learning models. Given
that geodesic distances represent a fundamentally different category of
physical information in addition to the structural information that is
available from an atomic configuration, this approach enables us to encode
some approximate knowledge of the potential energy surface a priori into the
hyperparameters of the descriptor, thereby maximizing the physical
informativeness of an otherwise limited, closed-form structural representation
of the atomistic system.

Given the utility of machine learning-based interatomic potentials, the
application of our approach to an atomistic machine learning task would
constitute a compelling demonstration of its utility. However, given the
extensive scope for design, testing, and analysis this entails, we leave this
avenue of exploration for subsequent work. Our goal with this work has been to
introduce a unique way of thinking about atomistic representations and the
modeling of potential energy landscapes. One may extend this concept to the
inclusion of yet other types of information, besides structural and energetic,
for instance, by treating state likelihood or prediction uncertainty as a
physical quantity that can be encoded into descriptor representations. This
may be considered similar in spirit to metadynamics28 or uncertainty-driven
methods that have been previously suggested in the context of active learning
for atomistic machine learning potentials.29

##  6 Data Availability Statement

The code necessary to reproduce the results of this study is openly available
online at https://github.com/gopal-
iyer/geodesic_descriptor_optimization/tree/main.

##  7 Supporting Information

The Supporting Information document contains plots demonstrating (1) the non-
equivalence of Euclidean and geodesic distance measures in the examples of
malonaldehyde and aspirin; and (2) plots depicting Bayesian optimization
samples for optimizing the hyperparameters of descriptor Schemes 2 and 3.

##  8 Disclosure Statement

The authors have no competing interests to declare.

G.R.I. (concept, derivation, implementation, drafting of the manuscript) was
funded by AFOSR Award Number FA9550-19-1-9999 and the Brown University
Chemistry Department Vince Wernig Fellowship. B.R. (concept, drafting of the
manuscript) was funded by NSF CTMC CAREER Award 2046744. This research was
conducted using computational resources and services at the Center for
Computation and Visualization, Brown University.

## References

  * Deringer et al. 2019 Deringer, V. L.; Caro, M. A.; Csányi, G. Machine Learning Interatomic Potentials as Emerging Tools for Materials Science. _Advanced Materials_ 2019, _31_ , 1902765
  * Zuo et al. 2020 Zuo, Y.; Chen, C.; Li, X.; Deng, Z.; Chen, Y.; Behler, J.; Csányi, G.; Shapeev, A. V.; Thompson, A. P.; Wood, M. A.; Ong, S. P. Performance and Cost Assessment of Machine Learning Interatomic Potentials. _The Journal of Physical Chemistry A_ 2020, _124_ , 731–745, PMID: 31916773
  * Unke et al. 2021 Unke, O. T.; Chmiela, S.; Sauceda, H. E.; Gastegger, M.; Poltavsky, I.; Schütt, K. T.; Tkatchenko, A.; Müller, K.-R. Machine Learning Force Fields. _Chemical Reviews_ 2021, _121_ , 10142–10186, PMID: 33705118
  * Anstine and Isayev 2023 Anstine, D. M.; Isayev, O. Machine Learning Interatomic Potentials and Long-Range Physics. _The Journal of Physical Chemistry A_ 2023, _127_ , 2417–2431, PMID: 36802360
  * Batzner et al. 2022 Batzner, S.; Musaelian, A.; Sun, L.; Geiger, M.; Mailoa, J. P.; Kornbluth, M.; Molinari, N.; Smidt, T. E.; Kozinsky, B. E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. _Nature Communications_ 2022, _13_ , 2453
  * Glielmo et al. 2017 Glielmo, A.; Sollich, P.; De Vita, A. Accurate interatomic force fields via machine learning with covariant kernels. _Phys. Rev. B_ 2017, _95_ , 214302
  * Dhaliwal et al. 2022 Dhaliwal, G.; Nair, P. B.; Singh, C. V. Machine learned interatomic potentials using random features. _npj Computational Materials_ 2022, _8_ , 7
  * De et al. 2016 De, S.; Bartók, A. P.; Csányi, G.; Ceriotti, M. Comparing molecules and solids across structural and alchemical space. _Phys. Chem. Chem. Phys._ 2016, _18_ , 13754–13769
  * Musil et al. 2021 Musil, F.; Grisafi, A.; Bartók, A. P.; Ortner, C.; Csányi, G.; Ceriotti, M. Physics-Inspired Structural Representations for Molecules and Materials. _Chemical Reviews_ 2021, _121_ , 9759–9815, PMID: 34310133
  * Behler 2011 Behler, J. Atom-centered symmetry functions for constructing high-dimensional neural network potentials. _The Journal of Chemical Physics_ 2011, _134_ , 074106
  * Behler and Parrinello 2007 Behler, J.; Parrinello, M. Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces. _Phys. Rev. Lett._ 2007, _98_ , 146401
  * Imbalzano et al. 2018 Imbalzano, G.; Anelli, A.; Giofré, D.; Klees, S.; Behler, J.; Ceriotti, M. Automatic selection of atomic fingerprints and reference configurations for machine-learning potentials. _The Journal of Chemical Physics_ 2018, _148_
  * Qi et al. 2024 Qi, J.; Ko, T. W.; Wood, B. C.; Pham, T. A.; Ong, S. P. Robust training of machine learning interatomic potentials with dimensionality reduction and stratified sampling. _npj Computational Materials_ 2024, _10_ , 43
  * Mehdi et al. 2024 Mehdi, S.; Smith, Z.; Herron, L.; Zou, Z.; Tiwary, P. Enhanced Sampling with Machine Learning. _Annual Review of Physical Chemistry_ 2024, _75_ , null, PMID: 38382572
  * Glielmo et al. 2022 Glielmo, A.; Zeni, C.; Cheng, B.; Csányi, G.; Laio, A. Ranking the information content of distance measures. _PNAS Nexus_ 2022, _1_ , pgac039
  * Ceriotti 2019 Ceriotti, M. Unsupervised machine learning in atomistic simulations, between predictions and understanding. _The Journal of Chemical Physics_ 2019, _150_ , 150901
  * Zhu et al. 2019 Zhu, X.; Thompson, K. C.; Martínez, T. J. Geodesic interpolation for reaction pathways. _The Journal of Chemical Physics_ 2019, _150_ , 164103
  * 18 Jónsson, H.; Mills, G.; Jacobsen, K. W. _Classical and Quantum Dynamics in Condensed Phase Simulations_ ; pp 385–404
  * Gorban and Tyukin 2018 Gorban, A. N.; Tyukin, I. Y. Blessing of dimensionality: mathematical foundations of the statistical physics of data. _Philos Trans A Math Phys Eng Sci_ 2018, _376_
  * Shahriari et al. 2016 Shahriari, B.; Swersky, K.; Wang, Z.; Adams, R. P.; de Freitas, N. Taking the Human Out of the Loop: A Review of Bayesian Optimization. _Proceedings of the IEEE_ 2016, _104_ , 148–175
  * Snoek et al. 2012 Snoek, J.; Larochelle, H.; Adams, R. P. Practical Bayesian Optimization of Machine Learning Algorithms. Advances in Neural Information Processing Systems. 2012
  * Bergstra et al. 2011 Bergstra, J.; Bardenet, R.; Bengio, Y.; Kégl, B. Algorithms for Hyper-Parameter Optimization. Advances in Neural Information Processing Systems. 2011
  * Head et al. 2018 Head, T. et al. scikit-optimize/scikit-optimize: v0.5.2. 2018; https://doi.org/10.5281/zenodo.1207017
  * Chmiela et al. 2023 Chmiela, S.; Vassilev-Galindo, V.; Unke, O. T.; Kabylda, A.; Sauceda, H. E.; Tkatchenko, A.; Müller, K.-R. Accurate global machine learning force fields for molecules with hundreds of atoms. _Science Advances_ 2023, _9_ , eadf0873
  * Kabsch 1976 Kabsch, W. A solution for the best rotation to relate two sets of vectors. _Acta Crystallographica Section A_ 1976, _32_ , 922–923
  * Pozdnyakov et al. 2020 Pozdnyakov, S. N.; Willatt, M. J.; Bartók, A. P.; Ortner, C.; Csányi, G.; Ceriotti, M. Incompleteness of Atomic Structure Representations. _Phys. Rev. Lett._ 2020, _125_ , 166001
  * Nigam et al. 2024 Nigam, J.; Pozdnyakov, S. N.; Huguenin-Dumittan, K. K.; Ceriotti, M. Completeness of atomic structure representations. _APL Machine Learning_ 2024, _2_ , 016110
  * Bussi and Laio 2020 Bussi, G.; Laio, A. Using metadynamics to explore complex free-energy landscapes. _Nature Reviews Physics_ 2020, _2_ , 200–212
  * Kulichenko et al. 2023 Kulichenko, M.; Barros, K.; Lubbers, N.; Li, Y. W.; Messerly, R.; Tretiak, S.; Smith, J. S.; Nebgen, B. Uncertainty-driven dynamics for active learning of interatomic potentials. _Nature Computational Science_ 2023, _3_ , 230–239