# Peptide Isomerization is Suppressed at the Air-Water Interface

Aditya N. Singh  Department of Chemistry, University of California, Berkeley,
CA 94720, USA  Chemical Sciences Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  David T. Limmer<EMAIL_ADDRESS>Department of Chemistry, University of California, Berkeley, CA 94720, USA
Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley,
CA 94720, USA  Materials Science Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  Kavli Energy Nanoscience Institute at
Berkeley, Berkeley, CA 94720, USA

(August 27, 2024)

###### Abstract

We use molecular dynamics simulations to study the thermodynamics and kinetics
of alanine dipeptide isomerization at the air-water interface.
Thermodynamically, we find an affinity of the dipeptide to the interface. This
affinity arises from stabilizing intramolecular interactions that become
unshielded as the dipeptide is desolvated. Kinetically, we consider the rate
of transitions between the $\alpha_{L}$ and $\beta$ conformations of alanine
dipeptide and evaluate it as a continuous function of the distance from the
interface using a recent extension of transition path sampling, TPS+U. The
rate of isomerization at the Gibbs dividing surface is suppressed relative to
the bulk by a factor of 3. Examination of the ensemble of transition states
elucidates the role of solvent degrees of freedom in mediating favorable
intramolecular interactions along the reaction pathway of isomerization. Near
the air-water interface, water is less effective at mediating these
intramolecular interactions.

A number of chemical and biological processes have been purported to be
catalyzed by air-water interfaces.1; 2; 3; 4; 5; 6; 7; 8; 9; 10 Evidence of
peptide-bond formation at the air-water interface has been reported using
infrared reflection spectroscopy1 and phosphorylation reactions have been
claimed to become thermodynamically spontaneous in microdroplets generated
with electrospray ionization.3; 10 The interface of heated gas microbubbles
has been shown to accumulate a range of biomolecules and act as a site for
their condensation, enrichment and assembly.4 The intrinsic physical and
chemical properties of the liquid-vapor interface have been implicated in
these studies, however understanding the relative roles of reduced
dimensionality and altered dielectric and steric environments on these
processes is difficult. Motivated by these observations, we study peptide
isomerization to investigate the effects of the air-water interface on a
simple unimolecular reaction.11; 12; 13; 14 In this way we isolate the role of
the altered solvation environment of the air-water interface on a biologically
relevant reaction, free of complications from diffusion limitations and
covalent bond breaking.

Computational and theoretical studies of reactions at interfaces are
challenging. Simulations of interfaces require large systems, and often
reactions of biological relevance are rare relative to the time-scales
accessible to ab initio methods15; 16; 17; 18; 19. As a means of understanding
some general effects of air-water interfaces on biochemical processes, we
consider the effect of the air-water interface on the conformational changes
of alanine dipeptide. Alanine dipeptide has long served as the model to
explore peptide dynamics, as it possesses many of the essential features of
larger proteins including flexible torsional angles and the presence of CO, NH
and CH3 groups that mediate both hydrophilic and hydrophobic interactions.13;
14 Previous studies have demonstrated the importance of solvent in the
reaction coordinate of isomerism20; 21; 22; 23; 24. Moreover, prior works on
the mechanism of peptide isomerization at water-hexane25, water-gold26 and
water-bilayer27 interfaces have found significant discrepancies in both
kinetic and thermodynamic observables near the interface. These
inconsistencies along with the aforementioned characteristics emphasize the
need of a systematic study of peptide isomerization at the air-water
interface.

In this work, we use molecular dynamics simulation to investigate the
thermodynamics and kinetics of alanine dipeptide at the air-water interface.
We observe an affinity of the peptide towards the interface, a phenomenon that
has been observed in previous studies of other interfaces26; 25. We find that
not all conformational states are equally stabilized, and particularly the
relative likelihood of observing conformations intermediate between long lived
dihedral states is decreased, suggesting a suppression of rate at the
interface. An extension of transition path sampling, TPS+U28 is used to probe
the relative rate of isomerization between the metastable conformations
$\alpha_{L}$ and $\beta$ as a function of the distance from the air-water
interface, without prior mechanistic assumptions or the conservation of the
mechanism near and away from the interface. This calculation confirms a 3-fold
slow-down of isomerization near the air-water interface relative to its bulk
value. We find an increase in favorable intramolecular peptide interactions
that stabilize particular dihedral states at interface. This favorable
intramolecular interaction is not observed for intermediate dihedral angles,
from which we conclude that both the increased stability as well as the
suppression of rate can be attributed to the availability of electrostatically
favorable conformations that are inaccessible in bulk due to shielding from
the solvent.

Figure 1: Stability of alanine dipeptide at the air-water interface. (A)
Representative snapshot of the molecular dynamics simulation. (B) The free
energy surface as a function of the dihedral angle $\phi$ and the distance
from the Gibbs dividing surface $z$. Positive values of $z$ correspond to the
bulk whereas negative values correspond to the vapor phase. The contours lines
are spaced by 2 $k_{\mathrm{B}}T$. The red dashed line in both figures denotes
the Gibbs dividing surface.

Simulation details. The AMBER ff14SB forcefield29 parameters and TIP3P water
model30 are used to describe the molecular interactions of a solvated alanine
dipeptide. The SETTLE algorithm31 is used to keep the water molecules rigid.
All simulations are performed in OpenMM32. The system is comprised of 3000
water molecules and 1 alanine dipeptide and are placed in a periodic box of
size 3x3x18 nm3. A cutoff of 1.4 nm is used for computing short range
interactions, and the Particle Mesh Ewald method33 is employed to calculate
long-range electrostatics. A Langevin integrator with a characteristic time-
scale of 1 ps-1 is used to maintain the temperature at 300 K for all
simulations. A velocity Verlet discretization34 is used to ensure time-
reversibility.

For all calculations, the system is first equilibrated for 1 ns. After
equilibration, the water slab is approximately 10.4 nm thick along the
coordinate orthogonal to the air-water interface. The Gibbs dividing surface,
defined as the plane where the density of water is half of that in bulk on
average, is located 5.2 nm away from the center of mass of all water
molecules. We employ a coordinate system where $z$ is defined as the distance
from the Gibbs dividing surface, with positive values corresponding to liquid
phase. A representative snapshot of the simulation is shown in Fig. 1(A).

We investigate the transition between two lived conformations of alanine
dipeptide, denoted the $\alpha_{L}$ to $\beta$ states. The C-N-Cα-C torsional
angle, denoted $\phi$, is used as the order parameter to separate these two
conformations. Specifically, we employ the indicator functions

$h_{\alpha}(t)=\Theta[\phi(t)]\;\;\;\;\;h_{\beta}(t)=\Theta[-\phi(t)]$ (1)

where $\phi$ is evaluated at time $t$ and $\Theta$ is the Heaviside function.
Within each basin defined by the above indicator functions in the bulk liquid,
both the $\alpha_{L}$ to $\beta$ states exhibit additional long lived
conformations distinguished by other intramolecular coordinates.35; 36
Thermodynamically, we must take care to ensure that sampling is ergodic within
each basin. Kinetically, all rate calculations are started in the $\alpha_{L}$
state, for which there is a direct path to the $\beta$, free of long-lived
intermediates.23

Thermodynamics near the air-water interface. Umbrella sampling is employed to
investigate the thermodynamics of the two metastable isomers of alanine
dipeptide near the air-water interface.37 Simulations are performed with a
harmonic bias on $\phi$ and the $z$ component of the center of mass of alanine
dipeptide of the form
$U(z,\phi)=K_{\phi}(\phi-\phi_{0})^{2}/2+K_{z}(z-z_{0})^{2}/2$ with $K_{\phi}$
= 500 kJ mol-1 rad-2 and $K_{z}=200$ kJ mol-1 nm-2. A total of 20 windows are
spaced 0.1 nm apart along $z$ using $z_{0}$ between -0.2 nm and 1.2 nm, and a
total of 36 windows spaced 0.18 rad apart along $\phi$ using $\phi_{0}$
between $-\pi$ and $\pi$, equating to 720 windows in total. Each window is run
for 5 ns and 1 ns of the simulation is discarded for equilibration. The
Weighted Histogram Analysis Method 38 is used for estimating the free energies
from the joint histograms in $\phi$ and $z$.

The free energy is shown in Fig. 1(B), defined from the joint probability,

$F(z,\phi)=-k_{\mathrm{B}}T\ln\langle\delta(z-z(0))\delta(\phi-\phi(0))\rangle$
(2)

expressed as an average, $\langle\dots\rangle$, over Dirac delta functions
$\delta$ within the canonical ensemble at fixed temperature $T$ and
$k_{\mathrm{B}}$ is Boltzmann’s constant. For all values of $z$, two stable
conformations are observed around $\phi=\pi/3$ and $\phi=-2\pi/3$ with a
barrier between them. These two regions correspond to the $\alpha_{L}$ and the
$\beta$ phase, respectively. The free energy plateaus for $z>$ 1 nm away from
the air-water interface, suggesting that the effect of the interface on the
peptide stability decays within a few molecular layers. We define
$z_{\textrm{B}}\,=\,1.15\,\textrm{nm}$ as the reference $z$ value for
estimating bulk properties.

The $\beta$ conformation is more stable than $\alpha_{L}$ by about 2.5
$k_{\mathrm{B}}T$, in accordance with previous work with a similar
forcefield.36 A considerable decrease in free energy is observed for moving
either conformation to the water-vapor interface. This can be quantified by
computing the free energy as a function of just $z$,

$F(z)=-k_{\mathrm{B}}T\ln\left[\int_{-\pi}^{\pi}e^{-\beta
F(z,\phi)}d\phi\right]$ (3)

by integrating over $\phi$. This relative free energy is shown in Fig. 2 and
reveals a negligible change from $z_{\mathrm{B}}$ to $z=$0.6 nm, below which
the free energy starts to decrease. A minima of -1.4 $k_{\mathrm{B}}T$
relative to the bulk is observed around 0.2 nm away from the Gibbs dividing
surface. This interfacial absorption is similar to studies of alanine
dipeptide at gold-water26 and water-hexane25 interfaces. While the
destabilization of the dipeptide in the vapor can be explained by the
unfavorable process of desolvation, the increased in stability near the air-
water interface reflects an interplay between energetic and entropic changes
to the full system, explored below.

The decrease in free energy near the interface is different for different
$\phi$ values of the peptide. Both the $\alpha_{L}$ and the $\beta$ phases
exhibit a 1.4 $k_{\mathrm{B}}T$ decrease in free energy near the interface, as
computed by restricting Eq. 3 over the domains associated with the indicator
functions in Eq. 1. We will refer to the $z$ dependent free energies for each
state as $F_{\alpha}(z)$ and $F_{\beta}(z)$. As the change is the same for
both conformers, the free energy difference between the two states,
$F_{\alpha}(z)-F_{\beta}(z)$ is constant, with a value of
$2.5\,k_{\mathrm{B}}T$ that is independent of $z$ for the range of $z$
considered. However, the free energy of the dividing surface between the two
phases, identified as $\phi=0$, only decreases by 0.4 $k_{\mathrm{B}}T$. This
inhomogeneous change in stability of different conformations results in an
increase in the barrier height by about 1 $k_{\mathrm{B}}T$ near the
interface, suggesting a suppression of fluctuations correlated with
transitions between $\beta$ and $\alpha_{L}$ at the interface.

Kinetics near the air-water interface. The presence of the air-water interface
breaks translational symmetry, and as a consequence the rate of isomerization
between the $\alpha_{L}$ and $\beta$ conformations can in principle depend on
the distance from the interface. In the limit that the diffusion of the
dipeptide along $z$ is minimal during the characteristic transition path
time,39; 40 we can associate an isomerization rate to each value of $z$. This
separation of timescales exists for alanine dipeptide as the translation
diffusion is small with a typical time to diffuse its diameter of nearly 1/2
ns41, permitting us to define and compute $k(z)$, the rate of isomerization as
a continuous function of $z$.

We define the spatially dependent rate constant using the general relationship
between rate constants and path partition functions.42 Specifically, for
observations times $t$ where the phenomenological rate constant is defined,
the relative rate at $z$ and $z^{\prime}$ is equal to

$\frac{k(z)}{k(z^{\prime})}=\frac{\Omega_{\alpha,\beta}(z)}{\Omega_{\alpha,\beta}(z^{\prime})}e^{\beta[F_{\alpha}(z)-F_{\alpha}(z^{\prime})]}$
(4)

where $\Omega_{\alpha,\beta}(z)$ is the partition function associated with
trajectories that start in $\alpha_{L}$ and end in $\beta$ at time $t$.This
path partition function takes the form

$\Omega_{\alpha,\beta}(z)=\int\mathcal{D}[\mathbf{X}(t)]P[\mathbf{X}(t)]h_{\alpha}(0)h_{\beta}(t)\delta[z-{z}(0)]$
(5)

where $P[\mathbf{X}(t)]$ is the probability to observe trajectory
$\mathbf{X}(t)$, and $\delta[z-z(0)]$ constrains the transition to occur at
$z$. This identification of the rate constant follows from the definition of a
transition probability, or the standard time correlation function
formalism.43; 44 TPS+U 28 provides a method to compute relative rates without
prior mechanistic assumptions using this relation to path partition functions.
The algorithm combines biasing potentials as in umbrella sampling with
Transition Path Sampling45; 42 to probe the relative rates, analogous to
calculation of relative free energies.

To perform TPS+U, we generated initial reactive trajectories by propagating
configurations at $\phi=0$ backward and forward in time for each value of $z$
studied. The length of the trajectories are kept fixed at 4 ps, which is
sufficient for the dipeptide to diffuse over the barrier. The Shooting from
the Top algorithm46 is used to generate new decorrelated trajectories from old
ones. This is done by randomly selecting configurations with $\phi$ in the
range of $-0.3\leq\phi\leq 0.3$ from the the old trajectory, drawing new
velocities from the Maxwell-Boltzmann distribution, and propagating it forward
and backward in time.

The number of shooting moves is varied according to $z$, as convergence is
slow near the interface. Approximately 6000 shooting moves are used near the
bulk, and around 12,000 moves are used for $z\leq 0.6$ nm. The acceptance rate
varies between 0.2 to 0.3, resulting in 2000-4000 independent reactive
trajectories. Harmonic biases of the same form used for the umbrella sampling
but with $K_{z}$ between 100 and 200 nm-2 used to constrain the trajectories
at a given $z$. In total, 32 windows with minima spaced by 0.05 nm between
$-0.2\leq z\leq 1.4$ nm are used together with histogram reweighting to
compute the relative path partition function as a function $z$.

Figure 2: Kinetic and thermodynamic observables as a function of $z$. The red
dashed line denotes the free energy relative to the bulk, $\Delta
F(z)=F(z)-F(z_{\mathrm{B}})$. The black line and blue markers corresponds to
the log ratio of the rate of $\alpha_{L}$ to $\beta$ transition as a function
of $z$, computed using TPS+U and Transition State Theory respectively. The
errorbars and the shading denote the standard error. Inset depicts the $z$
dependent transmission coefficient, $\kappa(z)$.

The results from TPS+U are shown in Fig. 2. The relative rate plateaus for
$z>$ 0.8 nm to the bulk value, analogously as $F(z)$. Around $z=$0.6 nm, the
rate decreases, and reaches a minima around $z=0.1$ nm before increasing
towards its eventual value in the vapor. The magnitude of the suppression of
the rate is approximately 2.5 times at $z=0$. The suppression of the
isomerization rate at the air-water interface is distinct from previous
observations of smaller polar molecules which exhibit slight isomerization
rate enhancements,47 and more consistent with observed suppression of ion pair
dissociation.48; 16; 15

The observed rate suppression can be anticipated from the destabilization of
conformations with $\phi=0$ observed in $F(z,\phi)$. We can relate the rate
suppression at the interface to this destabilization by comparing $k(z)$ with
that approximated by transition state theory.49; 43 If we assume that the
isomerization mechanism is invariant along $z$, and employ a diving surface as
$\phi^{{\ddagger}}=0$, transition state theory gives us an estimate of the
rate of $\alpha_{L}$ to $\beta$ isomerization. Specifically, using
$F(z,\phi)$, the transition state theory rate can be computed from

$k_{\mathrm{TST}}(z)=\frac{\langle|\dot{\phi}|\rangle}{2}\frac{e^{-\beta
F(z,0)}}{\int_{0}^{\pi}e^{-\beta F(z,\phi)}d\phi}$ (6)

where $\langle|\dot{\phi}|\rangle$ corresponds to the average angular
velocity, and the limits of the integral correspond to the $\alpha_{L}$ state.
We find that the denominator in Eq. 6 does not change with $z$. As a
consequence the log ratio of rates $\ln[k(z)/k(z_{\mathrm{B}})]$ is simply
equal to the difference in barrier heights. This estimate is shown in Fig. 2,
and is consistent with TPS+U calculations, illustrating that the change in
stability of conformations around $\phi=0$ is correlated with the fluctuations
that result in isomerization.

In order to determine the absolute rate, we have evaluated the rate in the
bulk using the Bennett-Chandler method.43 The Bennett-Chandler method corrects
the transition state theory approximation by defining a correlation function
called the transmission coefficient, $\kappa(z,t)$,

$\kappa(z,t)=\frac{\langle\dot{\phi}(0)\;\Theta[-\dot{\phi}(0)]\left[h_{\beta}(t)-h_{\beta}(-t)\right]\rangle_{\phi^{\ddagger},z}}{\langle\dot{\phi}(0)\Theta[-\dot{\phi}(0)]\rangle_{\phi^{\ddagger},z}}$
(7)

where $\langle\dots\rangle_{\phi^{{\ddagger}},z}$ denotes an average over
configurations on top of the dividing surface at a given value of $z$. The
transmission coefficient decays quickly and reaches a plateau denoted
$\kappa(z)$, and the value at the plateau determines the phenomenological
rate. Using the transmission coefficient, the Bennett-Chandler rate is given
by $k_{\mathrm{BC}}(z)=k_{\mathrm{TST}}(z)\kappa(z)$. To compute
$\kappa(z_{\mathrm{B}})$, 30 configurations on top of the dividing surface are
selected. Maxwell-Boltzmann velocities are drawn randomly, and the system is
propagated forward and backwards in time for 2 ps. A total of 2000 forward and
backward trajectories are run for each configuration. The transmission
coefficient plateaus within 0.3 ps seconds and is found to be small, 0.15
$\pm$ 0.03, yielding a rate in the bulk of 0.119 ns-1. Additional calculations
of $\kappa$ at different values of $z$ are unchanged from this bulk value as
shown in Fig. 2. Together with the observation of minimal changes in the
torsional diffusion constant near the interface, these findings suggest that
dynamical correlations between the solvent and $\phi$ are unaffected by the
interface.

Enhanced stability and rate suppression. As the relative suppression of rate
near the air-water interface can be quantitatively captured by the
$z$-dependence in the barrier height along the $\phi$ direction, the
thermodynamics of the dipeptide is further investigated. We have analyzed the
various components of the potential energy as a function of $z$ leveraging the
pairwise additivity of the potential. To understand how these different
interactions change near the interface, 10 ns simulations are performed with a
harmonic potential along a range of $z$ values and two $\phi$ values denoting
the $\alpha_{L}$ conformation and the dividing surface.

A prominent difference between the two states is observed in the average
intramolecular potential of the dipeptide, $\langle\Delta
U_{\mathrm{P-P}}(z)\rangle$ = $\langle U_{\mathrm{P-P}}(z)\rangle$ \- $\langle
U_{\mathrm{P-P}}(z_{\mathrm{B}})\rangle$, shown in Fig. 3(A). While the
average intramolecular interactions are invariant along $z$ for conformations
with $\phi=0$, these interactions decrease by around 2 $k_{\mathrm{B}}T$ for
the $\alpha_{L}$ conformation at the interface relative to the bulk. A further
decomposition of these interactions into bonded and nonbonded terms reveals
that this enhanced stability arises exclusively from electrostatic
interactions. These changes are not well approximated by a point dipole, due
to their short ranged nature and the complexity of the charge distribution of
alanine dipeptide. These results suggest that the decrease in water density at
the interface allows the dipeptide to access certain electrostatically
favorable conformations that are unfavorable in bulk. Analogous conformations
are not accessible to the barrier state. In this way, as the dipeptide is
desolvated, intramolecular interactions become unshielded by the water,
stabilizing the metastable states $\alpha_{L}$ and $\beta$ and destabilizing
configurations between them.

These observations are in line with multiple studies that have underscored the
role of water in the isomerization of alanine dipeptide.20; 21; 22; 23; 24
Bolhuis et al.21 analyzed reactive pathways of $\alpha_{R}$ to $\beta$
isomerization of alanine dipeptide in both vacuum and solvated phases and
found the dihedral angles to be a poor descriptor of the reaction coordinate.
Ma and Dinner22 suggested that information on the electrostatic interactions
between the water and peptide was essential to capture the isomerization in
explicit solvent. This was confirmed in a subsequent study by Velez-Vega et
al23. More recently, Jung et al.24 used neural networks with symmetry
functions to encapsulate the effect of water required to describe the reaction
coordinate.

Figure 3: Importance of intramolecular peptide degress of freedom. (A) Change
in the peptide-peptide potential energy as a function of $z$ for the barrier
and the $\alpha_{L}$ conformation. Errorbars denote the standard error
obtained from block averaging. (B) Commitor probability distribution for the
barrier along $\phi$ ($-0.15\leq\phi\leq 0.15$). (C) Visualization of the
dipeptide along with the water in its first solvation shell at the GDS (top)
and the bulk (bottom) is shown as reference.

In order to reconcile the previously observed importance of solvent degrees of
freedom with the success of transition state theory to predict the rate
profile, we have analyzed the commitor probability distribution conditioned
along $\phi=0$.44 The commitor, $p_{\beta}$, is a direct indicator of the
progression of a reaction. It is defined as the fraction of trajectories that
go to $\beta$ before $\alpha_{L}$. By computing the commitor probability
conditioned on a particular value of an order parameter, we can infer its
correlation with the transition state ensemble.44; 50 The commitor conditioned
on $\phi=0$ at $z$ is given by,

$p_{\beta}=\frac{\int\mathcal{D}\left[\mathbf{X}(t)\right]P\left[\mathbf{X}(t)\right]\delta\left[\phi(0)\right]\delta\left[z-z(0)\right]h_{\beta}(t)}{\int\mathcal{D}\left[\mathbf{X}(t)\right]P\left[\mathbf{X}(t)\right]\delta\left[\phi(0)\right]\delta\left[z-z(0)\right]}$
(8)

and calculated for 1000 configurations within a small window around $\phi=0$
that are selected from the TPS+U calculations. For each configuration we
compute an estimate of $p_{\beta}$ by propagating 10 trajectories with random
Maxwell-Boltzmann velocities, until they reach the $\alpha_{L}$ or $\beta$
state. The resultant distributions of $p_{\beta}$ from configurations drawn
from the bulk and the interface are shown in Fig. 3 (B).

A moderately flat distribution of commitment probabilities is observed in the
bulk, suggesting that the transition occurs diffusively over the the
thermodynamic barrier defined by $\phi=0$. This finding is consistent with
previous studies that highlight the importance of solvent degrees of freedom
in the reaction pathway of isomerization. By contrast, the distribution of
$p_{\beta}$ at the interface resembles a Gaussian distribution peaked at 0.5.
The probability of $p_{\beta}=0.5$ is nearly 3 times larger than that of
$p_{\beta}=0$ and $p_{\beta}=1$, suggesting that that $\phi$ becomes a better
reaction descriptor near the interface or that the effect of water becomes
less significant to the isomerization pathway near the interface.

The commitor calculations along with the potential decomposition suggest that
isomerization is modulated by the reorganization of intramolecular dipeptide
interactions. In bulk, these interactions are strongly dependent on the water
structure around the peptide, making the water coordinates important in
describing the reaction. Near the interface however, the role of the water
structure diminishes, making $\phi$ an accurate descriptor of the reaction
coordinate. However, due to the increased stability of the peptide near the
interface, the energy fluctuations required to induce this reorganization
become rarer, in turn suppressing the rate of peptide isomerization.

The physical and chemical properties of air-water interfaces can substantially
alter the kinetics of biologically relevant reactions. Using isomerization of
alanine dipeptide as a model, we have illuminated the interplay of
interactions that arise at the air-water interface. Electrostatic interactions
between the peptide dominate at the interface leading to a 4 times increase in
stability. These interactions are strongly dependent on the conformation of
the peptide, and we find the transition state of the $\alpha_{L}$ to $\beta$
isomerization to be unaffected by this stabilization. As a result, a three-
fold suppression in the rate is observed at the air-water interface.

Our observations of selective affinity of peptides to the air-water interface
as well as sluggish conformational dynamics they exhibit there both have their
origins in the altered ability of the air-water interface to solvate polar
solutes. These observations are consistent with similar mechanisms underlying
the suppression of rates of other charge organization processes like
hydrolysis15 and ion pair dissociation.16 Direct comparison to the
thermodynamics of adsorption and the role of energetic driving forces could be
accomplished with spectroscopic titration measurements like those done with
simple ions.51 Understanding the implication of these effects for more complex
chemistry like peptide bond formation1 and acid dissociation52; 53 is an
important direction for future work.

## Acknowledgements

ANS and DTL were supported by NSF Grant CHE2102314.

## References

  * Griffith and Vaida (2012) E. C. Griffith and V. Vaida, Proc. Natl. Acad. Sci. U. S. A. 109, 15697 (2012).
  * Vaida (2017) V. Vaida, Proc. Natl. Acad. Sci. U. S. A. 114, 12359 (2017).
  * Nam  _et al._ (2017) I. Nam, J. K. Lee, H. G. Nam, and R. N. Zare, Proc. Natl. Acad. Sci. U. S. A. 114, 12396 (2017).
  * Morasch  _et al._ (2019) M. Morasch, J. Liu, C. F. Dirscherl, A. Ianeselli, A. Kühnlein, K. L. Vay, P. Schwintek, S. Islam, M. K. Corpinot, B. Scheu, D. B. Dingwell, P. Schwille, H. Mutschler, M. W. Powner, C. B. Mast, and D. Braun, Nat. Chem. 11, 779 (2019), 2104.05537 .
  * Yadav  _et al._ (2020) M. Yadav, R. Kumar, and R. Krishnamurthy, Chem. Rev. 120, 4766 (2020).
  * Ruiz-Lopez  _et al._ (2020) M. F. Ruiz-Lopez, J. S. Francisco, M. T. Martins-Costa, and J. M. Anglada, Nat. Rev. Chem. 4, 459 (2020).
  * Guo  _et al._ (2021) W. Guo, A. B. Kinghorn, Y. Zhang, Q. Li, A. D. Poonam, J. A. Tanner, and H. C. Shum, Nat. Commun. 12, 1 (2021).
  * Lee  _et al._ (2019) J. K. Lee, D. Samanta, H. G. Nam, and R. N. Zare, J. Am. Chem. Soc. 141, 10585 (2019).
  * Stroberg and Schnell (2018) W. Stroberg and S. Schnell, Biophys. J. 115, 3 (2018).
  * Nam  _et al._ (2018) I. Nam, H. G. Nam, and R. N. Zare, Proc. Natl. Acad. Sci. U. S. A. 115, 36 (2018).
  * Srisailam  _et al._ (2002) S. Srisailam, T. K. S. Kumar, T. Srimathi, and C. Yu, J. Am. Chem. Soc. 124, 1884 (2002).
  * Dentino  _et al._ (1991) A. R. Dentino, P. A. Raj, K. K. Bhandary, M. E. Wilson, and M. Levine, J. Biol. Chem. 266, 18460 (1991).
  * Smith (1999) P. E. Smith, J. Chem. Phys. 111, 5568 (1999).
  * Brooks and Case (1993) C. Brooks and D. A. Case, Chem. Rev. 93, 2487 (1993).
  * Galib and Limmer (2021) M. Galib and D. T. Limmer, Science 371, 921 (2021).
  * Niblett and Limmer (2021) S. P. Niblett and D. T. Limmer, J. Phys. Chem. B 125, 2174 (2021).
  * Kattirtzi  _et al._ (2017) J. A. Kattirtzi, D. T. Limmer, and A. P. Willard, Proc. Natl. Acad. Sci. U. S. A. 114, 13374 (2017).
  * Serva  _et al._ (2018) A. Serva, S. Pezzotti, S. Bougueroua, D. R. Galimberti, and M.-P. Gaigeot, J. Mol. Struct. 1165, 71 (2018).
  * Huggins  _et al._ (2019) D. J. Huggins, P. C. Biggin, M. A. Dämgen, J. W. Essex, S. A. Harris, R. H. Henchman, S. Khalid, A. Kuzmanic, C. A. Laughton, J. Michel, _et al._ , WIREs Comput. Mol. Sci. 9, e1393 (2019).
  * Drozdov  _et al._ (2004) A. N. Drozdov, A. Grossfield, and R. V. Pappu, J. Am. Chem. Soc. 126, 2574 (2004).
  * Bolhuis  _et al._ (2000) P. G. Bolhuis, C. Dellago, and D. Chandler, Proc. Natl. Acad. Sci. U. S. A. 97, 5877 (2000).
  * Ma and Dinner (2005) A. Ma and A. R. Dinner, J. Phys. Chem. B 109, 6769 (2005).
  * Velez-Vega  _et al._ (2009) C. Velez-Vega, E. E. Borrero, and F. A. Escobedo, J. Chem. Phys. 130, 06B611 (2009).
  * Jung  _et al._ (2019) H. Jung, R. Covino, and G. Hummer, arXiv (2019), 1901.04595 .
  * Pohorille and Wilson (1994a) A. Pohorille and M. A. Wilson, Jerus. Sym. Q. , 207 (1994a).
  * Bellucci and Corni (2014) L. Bellucci and S. Corni, J. Phys. Chem. C 118, 11357 (2014).
  * Pohorille and Wilson (1994b) A. Pohorille and M. Wilson, in  _ACS Symposium Series_ , Vol. 568 (AMERICAN CHEMICAL SOCIETY, 1994) pp. 395–395.
  * Schile and Limmer (2019) A. J. Schile and D. T. Limmer, J. Chem. Phys. 150, 191102 (2019), 1902.05227 .
  * Maier  _et al._ (2015) J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser, and C. Simmerling, J. Chem. Theory Comput. 11, 3696 (2015).
  * Jorgensen  _et al._ (1983) W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983).
  * Miyamoto and Kollman (1992) S. Miyamoto and P. A. Kollman, J. Comput. Chem. 13, 952 (1992).
  * Eastman  _et al._ (2017) P. Eastman, J. Swails, J. D. Chodera, R. T. McGibbon, Y. Zhao, K. A. Beauchamp, L.-P. Wang, A. C. Simmonett, M. P. Harrigan, C. D. Stern, _et al._ , PLoS Comput. Biol. 13, e1005659 (2017).
  * Essmann  _et al._ (1995) U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995).
  * Sivak  _et al._ (2013) D. A. Sivak, J. D. Chodera, and G. E. Crooks, Phys. Rev. X 3, 011007 (2013).
  * Ferguson  _et al._ (2011) A. L. Ferguson, A. Z. Panagiotopoulos, P. G. Debenedetti, and I. G. Kevrekidis, J. Chem. Phys. 134, 04B606 (2011).
  * Vymětal and Vondrášek (2010) J. Vymětal and J. Vondrášek, J. Phys. Chem. B 114, 5632 (2010).
  * Frenkel and Smit (2001) D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ , Vol. 1 (Elsevier, 2001).
  * Kumar  _et al._ (1992) S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Kollman, J. Comput. Chem. 13, 1011 (1992).
  * Zwanzig (1990) R. Zwanzig, Accounts Chem. Res. 23, 148 (1990).
  * Chekmarev  _et al._ (2004) D. S. Chekmarev, T. Ishida, and R. M. Levy, J. Phys. Chem. B 108, 19487 (2004).
  * Rossky and Karplus (1979) P. J. Rossky and M. Karplus, J. Am. Chem. Soc. 101, 1913 (1979).
  * Dellago  _et al._ (1998) C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, J. Chem. Phys. 108, 1964 (1998).
  * Chandler (1978) D. Chandler, J. Chem. Phys. 68, 2959 (1978).
  * Dellago  _et al._ (2006) C. Dellago, P. G. Bolhuis, and P. L. Geissler, in  _Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 1_ (Springer, 2006) pp. 349–391.
  * Bolhuis  _et al._ (2002) P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Annu. Rev. Phys. Chem. 53, 291 (2002).
  * Jung  _et al._ (2017) H. Jung, K.-i. Okazaki, and G. Hummer, J. Chem. Phys. 147, 152716 (2017).
  * Benjamin and Pohorille (1993) I. Benjamin and A. Pohorille, J. Chem. Phys. 98, 236 (1993).
  * Dang and Schenter (2018) L. X. Dang and G. K. Schenter, J. Chem. Phys. 148, 222820 (2018).
  * Nitzan (2006) A. Nitzan,  _Chemical dynamics in condensed phases: relaxation, transfer and reactions in condensed molecular systems_ (Oxford university press, 2006).
  * Geissler  _et al._ (1999) P. L. Geissler, C. Dellago, and D. Chandler, J. Phys. Chem. B 103, 3706 (1999).
  * Otten  _et al._ (2012) D. E. Otten, P. R. Shaffer, P. L. Geissler, and R. J. Saykally, Proc. Natl. Acad. Sci. U. S. A. 109, 701 (2012).
  * Shamay  _et al._ (2007) E. S. Shamay, V. Buch, M. Parrinello, and G. L. Richmond, J. Am. Chem. Soc. 129, 12910 (2007).
  * Wang  _et al._ (2009) S. Wang, R. Bianco, and J. T. Hynes, J. Phys. Chem. A 113, 1295 (2009).


# Accelerating Materials Discovery with Bayesian Optimization and Graph Deep
Learning

Yunxing Zuo  Mingde Qin  Chi Chen  Weike Ye  Xiangguo Li  Jian Luo  Shyue Ping
Ong  Department of NanoEngineering, University of California San Diego, 9500
Gilman Dr, Mail Code 0448, La Jolla, CA 92093-0448, United States

###### Abstract

Machine learning (ML) models utilizing structure-based features provide an
efficient means for accurate property predictions across diverse chemical
spaces. However, obtaining equilibrium crystal structures typically requires
expensive density functional theory (DFT) calculations, which limits ML-based
exploration to either known crystals or a small number of hypothetical
crystals. Here, we demonstrate that the application of Bayesian optimization
with symmetry constraints using a graph deep learning energy model can be used
to perform “DFT-free” relaxations of crystal structures. Using this approach
to significantly improve the accuracy of ML-predicted formation energies and
elastic moduli of hypothetical crystals, two novel ultra-incompressible hard
materials MoWC2 ($P6_{3}/mmc$) and ReWB ($Pca2_{1}$) were identified and
successfully synthesized via in-situ reactive spark plasma sintering from a
screening of 399,960 transition metal borides and carbides. This work
addresses a critical bottleneck to accurate property predictions for
hypothetical materials, paving the way to ML-accelerated discovery of new
materials with exceptional properties.

###### keywords:

materials discovery, bayesian optimization, graph neural network, deep
learning

††journal: Materials Today

##  1 Introduction

The accurate prediction of novel stable crystals and their properties is a
fundamental goal in computation-guided materials discovery. While ab initio
approaches such as density functional theory (DFT) [1, 2] have been
phenomenally successful in this regard [3, 4, 5, 6], their high computational
cost and poor scalability have limited the broad application across vast
chemical and structural spaces. As a result, high-throughput DFT screening has
been mostly performed on $\sim O(100-1000)$ crystals with relatively small
unit cells.

To circumvent this limitation, machine learning (ML) has emerged as a new
paradigm for developing efficient surrogate models for predicting materials
properties at scale [7, 8, 9, 10, 11, 12]. Such ML models are usually trained
on large databases of materials properties [13, 14, 15] to learn the
relationship between input chemical and/or structural features and target
properties (e.g., formation energies, band gaps, elastic moduli, etc.). Only
ML models utilizing structural as well as chemical features can distinguish
between polymorphs and be universally applied in materials discovery across
diverse crystal structures. In particular, graph neural networks, where atoms
and bonds in crystals are represented as nodes and edges in a mathematical
graph, have emerged as a particularly promising approach with state-of-the-art
accuracy in predicting a broad range of energetic, electronic and mechanical
properties [16, 17, 18, 19, 20, 21, 22, 23].

Ironically, a critical bottleneck in the application of structure-based ML
models for materials discovery is the requirement for equilibrium crystal
structures as the inputs. These are obtained by “relaxing” initial input
structures along their potential energy surfaces, which are typically computed
via expensive DFT calculations. While there have been recent efforts [24, 25]
at deriving accurate interatomic forces from graph representations, the
application has been limited to a few molecular systems or constrained
chemical spaces.

Here, we propose a Bayesian Optimization With Symmetry Relaxation (BOWSR)
algorithm to obtain equilibrium crystal structures for accurate ML property
predictions without DFT. Utilizing a highly efficient MatErials Graph Network
(MEGNet) formation energy model [16], we demonstrate that BOWSR-relaxed
structures can serve as accurate inputs to ML property models, yielding far
higher accuracy in the prediction of various materials properties. Finally, we
demonstrate the power of this approach by screening $\sim$400,000 transition
metal borides and carbides for ultra-incompressible hard materials. Two new
materials with relatively low predicted energy above hull [26] were
successfully synthesized and demonstrated to have exceptional mechanical
properties, in line with the ML prediction.

##  2 Results

###  2.1 Bayesian Optimization With Symmetry Relaxation Algorithm

Bayesian optimization (BO) is an adaptive strategy for the global optimization
of a function. In crystal structure relaxation, the function to be optimized
is the potential energy surface, which expresses the energy of the crystal as
a function of the lattice parameters and atomic coordinates. In the BOWSR
algorithm, the symmetry (space group) of the lattice and the Wyckoff positions
of the atoms are constrained during the relaxation process, i.e., only the
independent lattice parameters and atomic coordinates are allowed to vary. The
BO goal is then the minimization of the following mapping:

$x:=\\{a,b,c,\alpha,\beta,\gamma,\vec{c_{1}},\vec{c_{2}},...\\}$ (1) $x_{\rm
opt}=\mathop{\mathrm{argmin}}_{x}U(x),U:R^{n}\xrightarrow{}R$ (2)

where {$a,b,c,\alpha,\beta,\gamma$} and {$\vec{c_{1}},\vec{c_{2}},...$}
represent the independent lattice parameters and the atomic positions for a
$P1$ crystal, respectively, and $U(.)$ is the energy of the system. The
schematic of the BOWSR algorithm as well as two examples of the geometry
parameterization for a high-symmetry cubic crystal and a low-symmetry
triclinic crystal are shown in Figure 1.

Figure 1: Bayesian Optimization With Symmetry Relaxation (BOWSR) algorithm. a,
The BOWSR algorithm parameterizes each crystal based on the independent
lattice parameters and atomic coordinates based on its space group. The
potential energy surface is then approximated by initializing a set of
training observations and energies from the ML energy model. Bayesian
optimization is then used to iteratively propose lower energy geometries based
on prior observations. b,  Two examples of the geometry parameterization for
cubic perovskite SrTiO3 and triclinic PdN2Cl2. For the high-symmetry cubic
perovskite, the lattice parameter $a$ is the only independent parameter, and
all atoms are in special Wyckoff positions with no degrees of freedom in the
fractional coordinates. For the triclinic crystal, all six lattice parameters
and all atomic coordinates are independent parameters.

The convergence accuracy and speed of BOWSR are set by the energy evaluator
$U(.)$, which can be any computational model, including ab initio methods,
empirical potentials, and surrogate ML models. In this work, we have elected
to use a MatErials Graph Network (MEGNet) formation energy model previously
trained on the DFT-computed formation energies of 133,420 Materials Project
crystals [13]. This MEGNet model has a cross-validated mean absolute error
(MAE) of 26 meV atom-1, which is among the best accuracy among general ML
models thus far [22, 18]. Examples of the convergence of the BOWSR algorithm
using the MEGNet energy model are shown in Figure S1.

###  2.2 Properties Prediction

Elemental substitution is a common, chemically intuitive approach to deriving
potential new compounds. For instance, the rock salt LiCl can be derived from
rock salt NaCl by substituting Na for the chemically similar Li. Here, we
demonstrate the potential for BOWSR to substantially enhance ML property
predictions of the formation energies and elastic moduli (bulk and shear
moduli) of substituted crystals. The dataset comprises a total of 12,277 and
1,672 unique crystals with pre-computed DFT formation energies and elastic
moduli, respectively, from the Materials Project [13]. These crystals belong
to 144 (35 binary, 91 ternary, and 18 quaternary) common structure prototypes
in the Inorganic Crystal Structure Database (ICSD) [27, 28]. Each prototype
comprises at least 30 compositions (statistical distribution shown in Figure
S2). For each crystal in the dataset (e.g., rock salt GeTe), another crystal
with the same prototype but a different composition (e.g., rock salt NaCl) was
selected at random and multi-element substitutions (Na$\rightarrow$Ge,
Cl$\rightarrow$Te) were performed to arrive at an “unrelaxed” structure. The
BOWSR algorithm was then applied to obtain a BOWSR-relaxed structure. The
unrelaxed, BOWSR-relaxed, and the original DFT-relaxed structures were then
used as inputs for property predictions using MEGNet models. These MEGNet
models were trained on the DFT-computed formation energies and elastic moduli
of 133,420 and 12,179 crystals, respectively, from the Materials Project.

Figure 2: Performance of MEGNet materials properties predictions for three
levels of relaxation. Three levels of relaxation are served as inputs for
MEGNet materials properties predictions — unrelaxed, BOWSR-relaxed, and DFT-
relaxed structures. a, b, and c  show the parity plot of MEGNet prediction on
formation energies, bulk moduli, and shear moduli, respectively. d, e, and f
show the errors distribution on the corresponding materials property.

Figure 2 compares the MEGNet model predictions using the unrelaxed, BOWSR-
relaxed, and DFT-relaxed structures as inputs with respect to DFT-computed
values. The mean absolute errors (MAEs) of the MEGNet models using the DFT-
relaxed structures provide a best-case performance baseline. It should be
noted that the MEGNet models were trained using a superset of data from the
Materials Project that includes the DFT-relaxed structures [13]. Hence, the
reported MAEs of MEGNet predictions using DFT-relaxed structures in this work
are much smaller than the previously reported MAEs of these MEGNet models from
cross-validation and should not be considered as a metric for MEGNet
performance. Using the unrelaxed structures as inputs results in much higher,
positively skewed MAEs in the MEGNet formation energy prediction compared to
using DFT-relaxed structures. This is because the unrelaxed structures have
lattice parameters and atomic positions that can deviate substantially from
the ground state DFT-relaxed structures, resulting in higher energies. Using
the BOWSR-relaxed structures as inputs reduces the MAEs by a factor of four,
from 363 meV atom-1 to 88 meV atom-1. The $R^{2}$ also substantially increases
from 0.67 to 0.96, and the error distribution is Gaussian-like with a mean
close to 0. Similarly, large improvements in the MEGNet predictions of the
elastic moduli are also observed using the BOWSR-relaxed structures compared
to using unrelaxed structures, with MAEs in the $\log_{10}K_{\rm VRH}$ and
$\log_{10}G_{\rm VRH}$ reducing by half. We tested the sensitivity of the
BOWSR algorithm to the initial structures used to perform elemental
substitution. Using four randomly chosen parent structures with different
lattice parameters for each prototype, the above procedures were repeated and
the results are shown in Figure S3. The BOWSR-relaxed structures exhibit
consistently low errors regardless of initial structure selection.

We also tested the sensitivity of the BOWSR algorithm to the accuracy of the
energy evaluator by artificially introducing Gaussian noise into the MEGNet
formation energy prediction. The energy errors from the BOWSR-relaxed
structures are linearly correlated with the errors of the surrogate ML model
with the root mean square error (RMSE) ranging from 27 to 1000 meV $\rm
atom^{-1}$ (Figure S4), which indicates the robustness of the BO propagation
and the broad applicability of the BOWSR algorithm to any surrogate ML models.
The same linear correlations are also observed between the elastic moduli
errors and the errors of the surrogate ML model, as shown in Figure S4c and
S4d.

###  2.3 Discovery of Ultra-incompressible Hard Materials

Figure 3: Flowchart of ultra-incompressible materials discovery leveraged by
the BOWSR algorithm and MEGNet models. The materials candidates were generated
by elemental substitution to structure prototypes. These candidates were
relaxed by the BOWSR algorithm and subsequently screened by their predicted
thermodynamic stability and mechanical properties. The screened candidates
were further verified by DFT calculations, and the high-ranked candidates were
directed to synthesis attempts. Quantity in the parenthesis shows the number
of candidates at each stage.

We used the BOWSR algorithm with the MEGNet models to rapidly screen hundreds
of thousands of candidates for novel ultra-incompressible hard materials, as
shown in Figure 3. Given that binary compounds have already been extensively
explored in the literature [29, 30], we targeted our search in 12 ternary
M${}^{\prime}_{x}$M${}^{\prime\prime}_{y}$Xz chemical spaces, where M′, M′′ =
Mo, W, Os, or Re and X = B or C. These elements were selected based on their
common occurrences in ultra-incompressible hard binary compounds. By
combinatorially applying elemental substitutions to 5,555 ternary structures
prototypes in the ICSD [27, 28], 399,960 candidates were generated and relaxed
using the BOWSR algorithm with the MEGNet energy model. The BOWSR-relaxed
candidates were then screened for stability and mechanical properties using
MEGNet property models. The stability metric used was the energy above hull
$E_{\rm hull}^{\rm MEGNet}$, which was computed using the predicted formation
energy $E_{\rm f}^{\rm MEGNet}$ with the 0 K phase diagram in the Materials
Project database [26, 13, 31]. At this intermediate stage, a relatively
generous threshold of $E_{\rm hull}^{\rm MEGNet}<100$ meV atom-1 was used to
obtain candidates that are likely to be thermodynamic stable [32]. Of these,
candidates with relatively high MEGNet-predicted bulk and shear moduli
($K_{\rm VRH}^{\rm MEGNet}>250$ GPa and $G_{\rm VRH}^{\rm MEGNet}>100$ GPa)
were identified. Similar to the stability criterion, the mechanical criteria
used are slightly lower than the conventional threshold for ultra-
incompressibility to account for the higher MAE of the MEGNet elastic moduli
predictions [30]. DFT relaxations and energy calculations were then carried
out on the 1,603 candidates that passed all three ML-based screening criteria.
Subsequently, expensive DFT elastic tensor calculations [33] were performed on
the 143 candidates that have DFT $E_{\rm hull}^{\rm DFT}<100$ meV atom-1.

Table 1 summarizes the computed elastic properties of the top ten candidates
with the highest computed bulk modulus together with other well-known ultra-
incompressible materials. Attempts were then made to synthesize all ten
candidates with eight unique compositions via in-situ reactive spark plasma
sintering (SPS) using elemental precursors in the appropriate ratios (see
Methods). Two crystals, MoWC2 ($P6_{3}/mmc$) and ReWB ($Pca2_{1}$), were
successfully synthesized and confirmed via X-ray diffraction (XRD, Figure 4a)
as single phase, while the synthesis of the other six compositions yielded
multiple phases (see Figure S5 - S10). Henceforth, we will refer to the two
novel phases of MoWC2 ($P6_{3}/mmc$) and ReWB ($Pca2_{1}$) simply as MoWC2 and
ReWB, respectively.

Table 1: DFT-computed bulk modulus ($\bm{K}_{\rm\textbf{VRH}}$), shear modulus
($\bm{G}_{\rm\textbf{VRH}}$), Young’s modulus ($\bm{E}_{\rm\textbf{VRH}}$),
Poisson’s ratio ($\bm{\nu}$) and energy above hull
($\bm{E}_{\rm\textbf{hull}}$) for the top 10 candidates with regard to
$\bm{K}_{\rm\textbf{VRH}}$ in descending order. MoWC2 and ReWB are bolded as
they are successfully synthesized by experiments. Some of the known ultra-
incompressible materials are used as references.

$K_{\rm VRH}$ (GPa) $G_{\rm VRH}$ (GPa) $E_{\rm VRH}$ (GPa) $\nu$ $E_{\rm
hull}$ (meV atom-1) Candidates ReOsB ($P\bar{6}m2$) 370.7 241.3 594.7 0.233
31.7 ReOsB2 ($P6_{3}/mmc$) 367.3 220.9 552.0 0.250 87.8 MoWC2
($\textbf{{P}}\textbf{6}_{\textbf{3}}/\textbf{{mmc}}$) 357.9 260.5 628.8 0.207
96.3 ReWB ($Fddd$) 356.8 176.9 455.5 0.287 20.6 Re13WB9 ($P\bar{6}m2$) 353.1
177.0 455.1 0.285 88.4 ReWB ($\textbf{{Pca}}\textbf{2}_{\textbf{1}}$) 352.6
144.1 380.4 0.320 33.1 OsWB ($Pbam$) 351.1 183.1 467.9 0.278 43.3 ReWB
($Cmce$) 350.9 161.5 420.1 0.301 32.6 Re6W7B8 ($P6/m$) 348.4 182.8 466.8 0.277
22.2 ReW2B2 ($P4/mbm$) 345.8 156.0 406.8 0.304 72.1 Known materials C
($Fd\bar{3}m$) 430.3 503.6 1086.9 0.079 136.4 WC ($P\bar{6}m2$) 389.8 280.0
677.8 0.210 1.1 BN ($F\bar{4}3m$) 370.1 382.8 852.4 0.116 77.3 ReB2
($P6_{3}/mmc$) 334.9 272.3 642.7 0.180 4.7

Figure 4: Two new materials proposed by the BOWSR algorithm and MEGNet
prediction confirmed by XRD characterizations. a, Measured and calculated XRD
patterns of two materials (ReWB and MoWC2). The major peaks are indexed for
reference. The pymatgen library [34] was used to calculate the XRD patterns of
the DFT-relaxed crystal structures. Minor shifts in peak positions between the
measured and calculated XRD patterns can be attributed to the systematic
errors between DFT and experimentally-measured lattice parameters. The
comparison in lattice parameters between BOWSR-relaxed, DFT-relaxed, and
experimental structures is shown in Table S1. b, Crystal structures and space
group of these two materials. The predicted energy above hull for ReWB and
MoWC2 are 66 and 7 meV atom-1, respectively. Figure 5: Experimental
measurements and theoretical prediction of mechanical properties for the new
materials. Comparisons of a,  bulk moduli and b,  shear moduli between MEGNet,
DFT and pulse-echo measurements. c, Measured Vickers microhardness of ReWB and
MoWC2 under loads ranging from 0.49 N to 9.8 N. d, Comparisons between the
hardness derived from MEGNet predicted bulk modulus and shear modulus [35] and
the hardness measured by Vickers indentation at the low load (0.49 N). The DFT
calculated results are referenced. The error bar represents the standard
deviation of multiple experimental measurements.

The mechanical properties of MoWC2 and ReWB were measured using the pulse-echo
method [36, 37]. As shown in Figure 5a and 5b, the experimentally-measured
bulk and shear moduli are in excellent agreement with both the MEGNet and DFT
predictions. Both new materials exhibit ultra-incompressibility, with bulk
modulus close to or larger than 300 GPa [30]. MoWC2 also exhibits high
estimated Vickers hardnesses $H_{v}$ of 36.6 at 0.49 N indentation load and
20.9 GPa at 9.8 N load (Figure 5c). ReWB has a comparatively lower measured
hardness of 29.5 at 0.49 N load and 17.6 GPa at 9.8 N load. The $H_{v}$ values
at 0.49 N load are within 20-25% of those derived from the MEGNet and DFT
predicted shear moduli via the empirical relation $H_{v}=0.151G$ [35], as
shown in Figure 5d.

##  3 Discussion

Many materials properties, such as formation energies, mechanical properties,
etc., exhibit a strong dependence on the crystal structure. However, obtaining
equilibrium crystal structures as inputs to accurate ML models still requires
expensive ab initio computations. By coupling an accurate MEGNet energy model
with Bayesian optimization of symmetry-constrained parameters, we demonstrate
that the new algorithm can reasonably approximate equilibrium structures. The
resulting substantial improvements in ML property predictions enable the rapid
screening of $\sim$400,000 candidate crystals for stability and exceptional
mechanical properties, $10^{3}$ \- $10^{4}$ orders of magnitude larger than
that accessible by high-throughput DFT calculations.

The effectiveness of the BOWSR algorithm is limited by the accuracy of the
energy evaluator. While the MEGNet formation energy model has been selected in
this work, we foresee that the future development of more accurate ML energy
models may improve the quality of the BOWSR-relaxed structures, and subsequent
ML property predictions, even further. Here, the search for ultra-
incompressible materials has been chosen as a model problem due to the high
cost of acquisition of elastic moduli via standard DFT approaches, but the
approach outlined can be readily extended to any property for which a reliable
ML model can be developed. It should be noted, however, that there is an
inverse relationship between the cost of acquisition and the training data
size; hence, datasets on expensive properties (e.g., elastic moduli, optical
properties, etc.) tend to be much smaller in size compared to cheaper
properties, making it more difficult to build reliable ML models for high-cost
properties. While approaches such as transfer learning or multi-fidelity
models have been shown to mitigate this trade-off to some extent [16, 17, 38,
39], the generally higher errors in ML models for high-cost properties should
be factored into the screening process in the thresholds.

##  4 Acknowledgements

This research was primarily supported by the National Science Foundation
Materials Research Science and Engineering Center program through the UC
Irvine Center for Complex and Active Materials (DMR-2011967). The authors also
acknowledge the use of data and software resources from the Materials Project,
funded by the U.S. Department of Energy, Office of Science, Office of Basic
Energy Sciences, Materials Sciences and Engineering Division under contract
no. DE-AC02-05-CH11231: Materials Project program KC23MP, and computing
resources provided by the Extreme Science and Engineering Discovery
Environment (XSEDE) under grant ACI-1548562.

##  5 Material and Methods

###  5.1 Bayesian Optimization With Symmetry Relaxation Algorithm

Geometry relaxation of a crystal structure requires the optimization of up to
$3N+6$ variables - six lattice parameters and three fractional coordinates for
each of the $N$ atoms. By constraining the symmetry to remain unchanged during
relaxation can reduce the number of independent variables considerably [40,
41, 42]. The open-source spglib [43] package was used for symmetry
determination. The search for optimized symmetry-constrained lattice
parameters and atomic coordinates that minimize the total energy was then
carried out via Bayesian optimization (BO). The changes in the variables were
used as the optimization inputs to reduce the tendency of the BO process being
dominated by parameters with large magnitudes. This approach has been
previously used for geometry optimization along the imaginary phonon modes
[42].

Using a Latin hypercube sampling, a set of training observations
$D\sim\\{(x_{i},U(x_{i}))\ i=1:m\\}$ were initialized, where the $x_{i}$ are
the $m$ independent lattice parameters and atomic coordinates and $U(.)$ is
the energy of the corresponding structure evaluated by the surrogate model
(see Equations 1 and 2). The BO strategy comprises two steps [44]:

  1. 1.

A Gaussian process (GP) model is trained on the initialized training
observations $D$ to approximate the $U(x)$. The Rational Quadratic kernel [45]
is adopted as the covariance function of GP. The noise level of GP model is
set to the root mean square error (RMSE) of the energy model.

  2. 2.

The acquisition function that balances the exploitation and exploration is
calculated for samples in the search space apart from the training
observations and the candidate with optimal acquisition function is proposed
to be evaluated (formation energy prediction by surrogate ML model).
Exploitation represents the samples with high predicted mean from the GP, and
exploration accounts the samples with high predictive uncertainty [46, 47,
48]. Here, we chose the expected improvement as the acquisition function,
which can be analytically expressed as [49, 50]:

$E[I(x)]=(\mu(x)-U(x^{+})-\xi)\cdot\Phi(Z)+\sigma(x)\cdot\phi(Z)$ (3)

and

$x^{+}=\mathop{\mathrm{argmax}}_{i=1,\dots,n}U(x_{i})$ (4)
$Z=\frac{\mu(x)-U(x^{+})-\xi}{\sigma(x)}$ (5)

where $\mu(x)$ and $\sigma(x)$ are the mean and standard deviation of the
posterior distribution on $x$ from the GP, respectively, and $\Phi(x)$ and
$\phi(x)$ are the cumulative distribution function (CDF) and probability
density function (PDF), respectively. The $\xi$ parameter can be tuned to
balance the trade-off between the first term (exploitation) and the second
term (exploration) in Equation (3). Until the maximum number of iteration
steps is reached, the sample with optimal acquisition function was iteratively
augmented to the training observations and used to update the GP surrogate
model in the next loop.

###  5.2 Materials Graph Network Models

The Materials Graph Network (MEGNet) models used in this work are based on the
same architecture as our previous work [16]. Briefly, three graph
convolutional layers with [64, 64, 32] neurons were used in each update
function, and the shifted softplus function was used as the non-linear
activation function. A set2set readout function with two passes was used after
the graph convolution steps. The cutoff radius for constructing the neighbor
graph was 5 Å. The MEGNet formation energy ($E_{f}$) and elasticity
($K_{\mathrm{VRH}}$ and $G_{\mathrm{VRH}}$) models were trained using the
2019.4.1 version of Materials Project database [13] containing 133,420
structure-formation energy and 12,179 structure-bulk/shear modulus data pairs.
Each dataset was split into 80%:10%:10% train:validation:test ratios. During
the model training, we used a batch size of 128 structures, and set the
initial learning rate to 0.001 in the Adam optimizer. All models were trained
for a maximum of 1500 epochs with an early stopping callback, which terminates
the model training if the validation error does not reduce for 500 consecutive
steps. The mean absolute errors (MAEs) of $E_{f}$,
$\log_{10}(K_{\mathrm{VRH}})$ and $\log_{10}(G_{\mathrm{VRH}})$ models in test
data are 26 meV atom-1, 0.07, and 0.12, respectively.

###  5.3 DFT calculations

The DFT relaxations, energy and elastic tensor calculations for the small
number of candidates that passed the ML screening were carried out using
Vienna ab initio simulation package (VASP) [51] within the projector augmented
wave approach [52]. The exchange-correlation interaction was described using
the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) [53]
functional for structural relaxations and energy calculations. The plane wave
energy cutoff was set to 520 eV, and the k-point density of at least 1,000 per
number of atoms was used. All structures were relaxed with energies and forces
converged to $10^{-5}$ eV and 0.01 eV/Å, respectively, consistent with the
calculation setting used in the Materials Project [13]. The elastic tensor
calculations were performed using the procedure described in previous work
[33]. A tighter energy convergence criterion of $10^{-7}$ eV was used, and
strains with magnitude of (-1%, -0.5%, 0.5%, 1%) were applied to each of the 6
independent components of strain tensor.

###  5.4 Synthesis

Bulk specimens of candidates ReOsB, ReOsB2, MoWC2, ReWB, Re13WB9, OsWB,
Re6W7B8, and ReW2B2 were synthesized via in-situ reactive spark plasma
sintering (SPS). Elemental powders of Mo, W ($>$99.5% purity, $\sim$325 mesh,
Alfa Aesar), Re ($\sim$99.99% purity, $\sim$325 mesh, Strem Chemicals), Os
($\sim$99.8% purity, $\sim$200 mesh, Alfa Aesar), boron ($\sim$99% purity, 1-2
$\mu$m, US Research Nanomaterials), and graphite ($\sim$99.9% purity, 0.4-1.2
$\mu$m, US Research Nanomaterials) were utilized as precursors. For each
composition, stoichiometric amounts of elemental powders were weighted out in
batches of 5 g. The powders were first mixed by a vortex mixer, and then high
energy ball milled (HEBM) in a Spex 8000D mill (SpexCertPrep) by tungsten
carbide lined stainless steel jars as well as 11.2 mm tungsten carbide milling
media (ball-to-powder ratio $\approx$ 4.5:1) for 50 min. 0.05 g or $\sim$1 wt%
of stearic acid was used as lubricant in the milling process. After HEMB, the
as-milled powder mixtures were loaded into 10 mm graphite dies lined with
graphite foils in batches of 2.5 g, and subsequently consolidated into dense
pellets via SPS in vacuum ($<10^{-2}$ Torr) by a Thermal Technologies 3000
series SPS machine. The HEBM and powder handing were conducted in an argon
atmosphere (with O2 level $<$ 10 ppm) to prevent oxidation.

During the SPS, specimens were initially heated to 1400 ∘C at a rate of 100
∘C/min under constant pressure of 10 MPa. For the final densification, the
temperature was subsequently raised at a constant rate of 30 ∘C/min to a final
isothermal sintering temperature, which was set at different levels for
different target compositions — 1800 ∘C (ReWB), 1700 ∘C (MoWC2 and Re6W7B8),
1600 ∘C (Re13WB9 and ReW2B2), or 1500 ∘C (ReOsB, ReOsB2, and OsWB), and
maintained isothermally for 10 min. Meanwhile, the pressure was increased to
50 MPa at a ramp rate of 5 MPa/min. The final densification temperature was
optimized for each specimen to achieve a high relative density while prevent
specimen melting due to overheating. The in-situ reactions between elemental
precursors likely took place during the initial temperature ramping. After
sintering, the specimens were cooled down naturally inside the SPS machine
(with power off).

###  5.5 Experimental Characterization

Sintered specimens were first ground to remove the carbon-contaminated surface
layer from the graphite tooling, and polished for further characterizations.
X-ray diffraction (XRD) experiments were conducted using a Rigaku Miniflex
diffractometer with the Cu K$\alpha$ radiation at 30 kV and 15 mA. The Vickers
microhardness tests were carried out on a LECO diamond microindentor with
loading force varying from 0.49 N (50 gf) to 9.8 N (1 kgf) and constant
holding time of 15 s, abiding by the ASTM Standard C1327. Over 20 measurements
at different locations were conducted for each specimen at each indentation
load to ensure statistical validity and minimize the microstructural and grain
boundary effects. In particular, over 30 measurements were conducted for each
specimen at 9.8 N indentation load.

The Young’s and shear moduli of the specimens were calculated from the
ultrasonic velocities measured with a Tektronix TDS 420A digital oscilloscope,
following the ASTM standard A494-15. Multiple measurements were conducted at
different locations.

##  6 Data Availability

The raw/processed data required to reproduce these findings cannot be shared
at this time due to legal or ethical reasons.

##  7 References

## References

  * [1] W. Kohn, L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys. Rev. 140 (4A) (1965) A1133–A1138.  doi:10.1103/PhysRev.140.A1133.
  * [2] P. Hohenberg, W. Kohn, Inhomogeneous Electron Gas, Phys. Rev. 136 (3B) (1964) B864–B871.  doi:10.1103/PhysRev.136.B864.
  * [3] J. Greeley, T. F. Jaramillo, J. Bonde, I. Chorkendorff, J. K. Nørskov, Computational high-throughput screening of electrocatalytic materials for hydrogen evolution, Nature Mater 5 (11) (2006) 909–913.  doi:10.1038/nmat1752.
  * [4] R. Gautier, X. Zhang, L. Hu, L. Yu, Y. Lin, T. O. L. Sunde, D. Chon, K. R. Poeppelmeier, A. Zunger, Prediction and accelerated laboratory discovery of previously unknown 18-electron ABX compounds, Nature Chem 7 (4) (2015) 308–316.  doi:10.1038/nchem.2207.
  * [5] Q. Yan, J. Yu, S. K. Suram, L. Zhou, A. Shinde, P. F. Newhouse, W. Chen, G. Li, K. A. Persson, J. M. Gregoire, J. B. Neaton, Solar fuels photoanode materials discovery by integrating high-throughput theory and experiment, Proc Natl Acad Sci USA 114 (12) (2017) 3040–3043.  doi:10.1073/pnas.1619940114.
  * [6] Z. Wang, J. Ha, Y. H. Kim, W. B. Im, J. McKittrick, S. P. Ong, Mining Unexplored Chemistries for Phosphors for High-Color-Quality White-Light-Emitting Diodes, Joule 2 (5) (2018) 914–926.  doi:10.1016/j.joule.2018.01.015.
  * [7] A. Seko, A. Togo, H. Hayashi, K. Tsuda, L. Chaput, I. Tanaka, Prediction of Low-Thermal-Conductivity Compounds with First-Principles Anharmonic Lattice-Dynamics Calculations and Bayesian Optimization, Phys. Rev. Lett. 115 (20) (2015) 205901.  doi:10.1103/PhysRevLett.115.205901.
  * [8] D. Xue, P. V. Balachandran, J. Hogden, J. Theiler, D. Xue, T. Lookman, Accelerated search for materials with targeted properties by adaptive design, Nat Commun 7 (1) (2016) 11241.  doi:10.1038/ncomms11241.
  * [9] F. A. Faber, A. Lindmaa, O. A. von Lilienfeld, R. Armiento, Machine Learning Energies of 2 Million Elpasolite ( A B C 2 D 6 ) Crystals, Phys. Rev. Lett. 117 (13) (2016) 135502.  doi:10.1103/PhysRevLett.117.135502.
  * [10] B. Sanchez-Lengeling, A. Aspuru-Guzik, Inverse molecular design using machine learning: Generative models for matter engineering, Science 361 (6400) (2018) 360–365.  doi:10.1126/science.aat2663.
  * [11] J. Schmidt, M. R. G. Marques, S. Botti, M. A. L. Marques, Recent advances and applications of machine learning in solid-state materials science, npj Comput Mater 5 (1) (2019) 83.  doi:10.1038/s41524-019-0221-0.
  * [12] B. Kim, S. Lee, J. Kim, Inverse design of porous materials using artificial neural networks, Sci. Adv. 6 (1) (2020) eaax9324.  doi:10.1126/sciadv.aax9324.
  * [13] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, K. A. Persson, Commentary: The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 1 (1) (2013) 011002.  doi:10.1063/1.4812323.
  * [14] J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, C. Wolverton, Materials Design and Discovery with High-Throughput Density Functional Theory: The Open Quantum Materials Database (OQMD), JOM 65 (11) (2013) 1501–1509.  doi:10.1007/s11837-013-0755-4.
  * [15] S. Curtarolo, W. Setyawan, S. Wang, J. Xue, K. Yang, R. H. Taylor, L. J. Nelson, G. L. Hart, S. Sanvito, M. Buongiorno-Nardelli, N. Mingo, O. Levy, AFLOWLIB.ORG: A distributed materials properties repository from high-throughput ab initio calculations, Computational Materials Science 58 (2012) 227–235.  doi:10.1016/j.commatsci.2012.02.002.
  * [16] C. Chen, W. Ye, Y. Zuo, C. Zheng, S. P. Ong, Graph Networks as a Universal Machine Learning Framework for Molecules and Crystals, Chem. Mater. 31 (9) (2019) 3564–3572.  doi:10.1021/acs.chemmater.9b01294.
  * [17] C. Chen, Y. Zuo, W. Ye, X. Li, S. P. Ong, Learning properties of ordered and disordered materials from multi-fidelity data, Nat Comput Sci 1 (1) (2021) 46–53.  doi:10.1038/s43588-020-00002-x.
  * [18] T. Xie, J. C. Grossman, Crystal Graph Convolutional Neural Networks for an Accurate and Interpretable Prediction of Material Properties, Phys. Rev. Lett. 120 (14) (2018) 145301.  doi:10.1103/PhysRevLett.120.145301.
  * [19] D. Duvenaud, D. Maclaurin, J. Aguilera-Iparraguirre, R. Gomez-Bombarelli, T. Hirzel, A. Aspuru-Guzik, R. P. Adams, Convolutional Networks on Graphs for Learning Molecular Fingerprints 9.
  * [20] S. Kearnes, K. McCloskey, M. Berndl, V. Pande, P. Riley, Molecular graph convolutions: Moving beyond fingerprints, J Comput Aided Mol Des 30 (8) (2016) 595–608.  doi:10.1007/s10822-016-9938-8.
  * [21] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, G. E. Dahl, Neural Message Passing for Quantum Chemistry, ArXiv170401212 CsarXiv:1704.01212.
  * [22] K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, K.-R. Müller, SchNet – A deep learning architecture for molecules and materials, The Journal of Chemical Physics 148 (24) (2018) 241722\.  doi:10.1063/1.5019779.
  * [23] P. B. Jørgensen, K. W. Jacobsen, M. N. Schmidt, Neural Message Passing with Edge Updates for Predicting Properties of Molecules and Materials, ArXiv180603146 Cs StatarXiv:1806.03146.
  * [24] K. T. Schütt, P. Kessel, M. Gastegger, K. A. Nicoli, A. Tkatchenko, K.-R. Müller, SchNetPack: A Deep Learning Toolbox For Atomistic Systems, J. Chem. Theory Comput. 15 (1) (2019) 448–455.  doi:10.1021/acs.jctc.8b00908.
  * [25] G. Cheon, L. Yang, K. McCloskey, E. J. Reed, E. D. Cubuk, Crystal Structure Search with Random Relaxations Using Graph Networks, ArXiv201202920 Cond-Mat PhysicsphysicsarXiv:2012.02920.
  * [26] S. P. Ong, L. Wang, B. Kang, G. Ceder, Li-Fe-P-O 2 Phase Diagram from First Principles Calculations, Chem. Mater. 20 (5) (2008) 1798–1807.  doi:10.1021/cm702327g.
  * [27] G. Bergerhoff, R. Hundt, R. Sievers, I. D. Brown, The inorganic crystal structure data base, J. Chem. Inf. Model. 23 (2) (1983) 66–69.  doi:10.1021/ci00038a003.
  * [28] A. Belsky, M. Hellenbrandt, V. L. Karen, P. Luksch, New developments in the Inorganic Crystal Structure Database (ICSD): Accessibility in support of materials research and design, Acta Cryst Sect A Found Cryst 58 (3) (2002) 364–369.  doi:10.1107/S0108768102006948.
  * [29] G. Akopov, M. T. Yeung, R. B. Kaner, Rediscovering the Crystal Chemistry of Borides, Adv. Mater. 29 (21) (2017) 1604506.  doi:10.1002/adma.201604506.
  * [30] M. T. Yeung, R. Mohammadi, R. B. Kaner, Ultraincompressible, Superhard Materials, Annu. Rev. Mater. Res. 46 (1) (2016) 465–485.  doi:10.1146/annurev-matsci-070115-032148.
  * [31] W. Ye, C. Chen, Z. Wang, I.-H. Chu, S. P. Ong, Deep neural networks for accurate predictions of crystal stability, Nat Commun 9 (1) (2018) 3800.  doi:10.1038/s41467-018-06322-x.
  * [32] W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain, W. D. Richards, A. C. Gamst, K. A. Persson, G. Ceder, The thermodynamic scale of inorganic crystalline metastability, Sci. Adv. 2 (11) (2016) e1600225.  doi:10.1126/sciadv.1600225.
  * [33] M. de Jong, W. Chen, T. Angsten, A. Jain, R. Notestine, A. Gamst, M. Sluiter, C. Krishna Ande, S. van der Zwaag, J. J. Plata, C. Toher, S. Curtarolo, G. Ceder, K. A. Persson, M. Asta, Charting the complete elastic properties of inorganic crystalline compounds, Sci Data 2 (1) (2015) 150009.  doi:10.1038/sdata.2015.9.
  * [34] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, G. Ceder, Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis, Computational Materials Science 68 (2013) 314–319.  doi:10.1016/j.commatsci.2012.10.028.
  * [35] X.-Q. Chen, H. Niu, D. Li, Y. Li, Intrinsic Correlation between Hardness and Elasticity in Polycrystalline Materials and Bulk Metallic Glasses, Intermetallics 19 (9) (2011) 1275–1281.  arXiv:1102.4063, doi:10.1016/j.intermet.2011.03.026.
  * [36] H. B. Huntington, Ultrasonic Measurements on Single Crystals, Phys. Rev. 72 (4) (1947) 321–331.  doi:10.1103/PhysRev.72.321.
  * [37] D. Lazarus, The Variation of the Adiabatic Elastic Constants of KCl, NaCl, CuZn, Cu, and Al with Pressure to 10,000 Bars, Phys. Rev. 76 (4) (1949) 545–553.  doi:10.1103/PhysRev.76.545.
  * [38] G. Pilania, J. Gubernatis, T. Lookman, Multi-fidelity machine learning models for accurate bandgap predictions of solids, Computational Materials Science 129 (2017) 156–163.  doi:10.1016/j.commatsci.2016.12.004.
  * [39] D. Jha, K. Choudhary, F. Tavazza, W.-k. Liao, A. Choudhary, C. Campbell, A. Agrawal, Enhancing materials property prediction by leveraging computational and experimental data using deep transfer learning, Nat Commun 10 (1) (2019) 5316.  doi:10.1038/s41467-019-13297-w.
  * [40] Y. Wang, J. Lv, L. Zhu, Y. Ma, Crystal structure prediction via particle-swarm optimization, Phys. Rev. B 82 (9) (2010) 094116.  doi:10.1103/PhysRevB.82.094116.
  * [41] C. W. Glass, A. R. Oganov, N. Hansen, USPEX—Evolutionary crystal structure prediction, Computer Physics Communications 175 (11-12) (2006) 713–720.  doi:10.1016/j.cpc.2006.07.020.
  * [42] M.-O. Lenz, T. A. R. Purcell, D. Hicks, S. Curtarolo, M. Scheffler, C. Carbogno, Parametrically constrained geometry relaxations for high-throughput materials science, npj Comput Mater 5 (1) (2019) 123.  doi:10.1038/s41524-019-0254-4.
  * [43] A. Togo, I. Tanaka, $\texttt{}Spglib{}$: A software library for crystal symmetry search, ArXiv180801590 Cond-MatarXiv:1808.01590.
  * [44] T. Lookman, P. V. Balachandran, D. Xue, R. Yuan, Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design, npj Comput Mater 5 (1) (2019) 21.  doi:10.1038/s41524-019-0153-8.
  * [45] C. E. Rasmussen, C. K. I. Williams, Gaussian Processes for Machine Learning, Adaptive Computation and Machine Learning, MIT Press, Cambridge, Mass, 2006\.
  * [46] R. Martinez-Cantin, N. de Freitas, E. Brochu, J. Castellanos, A. Doucet, A Bayesian exploration-exploitation approach for optimal online sensing and planning with a visually guided mobile robot, Auton Robot 27 (2) (2009) 93–103.  doi:10.1007/s10514-009-9130-2.
  * [47] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, N. de Freitas, Taking the Human Out of the Loop: A Review of Bayesian Optimization, Proc. IEEE 104 (1) (2016) 148–175.  doi:10.1109/JPROC.2015.2494218.
  * [48] P. V. Balachandran, D. Xue, J. Theiler, J. Hogden, T. Lookman, Adaptive Strategies for Materials Design using Uncertainties, Sci Rep 6 (1) (2016) 19660.  doi:10.1038/srep19660.
  * [49] D. R. Jones, M. Schonlau, Efficient Global Optimization of Expensive Black-Box Functions 38.
  * [50] T. Lookman, F. J. Alexander, K. Rajan (Eds.), Information Science for Materials Discovery and Design, Vol. 225 of Springer Series in Materials Science, Springer International Publishing, Cham, 2016.  doi:10.1007/978-3-319-23871-5.
  * [51] G. Kresse, J. Furthmüller, Efficient iterative schemes for _ab initio_ total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (16) (1996) 11169–11186.  doi:10.1103/PhysRevB.54.11169.
  * [52] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50 (24) (1994) 17953–17979.  doi:10.1103/PhysRevB.50.17953.
  * [53] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77 (18) (1996) 3865–3868.  doi:10.1103/PhysRevLett.77.3865.