# A machine learning potential for simulating infrared spectra of nanosilicate
clusters

Zeyuan Tang  Center for Interstellar Catalysis, Department of Physics and
Astronomy, Aarhus University, Ny Munkegade 120, Aarhus C 8000, Denmark  Stefan
T. Bromley  Departament de Ciència de Materials i Química Física & Institut de
Química Teòrica i Computatcional (IQTCUB), Universitat de Barcelona, c/Martí i
Franquès 1-11, 08028 Barcelona, Spain  Institució Catalana de Recerca i
Estudis Avançats (ICREA), Passeig Lluis Companys 23, 08010 Barcelona, Spain
Bjørk Hammer<EMAIL_ADDRESS>Center for Interstellar Catalysis, Department
of Physics and Astronomy, Aarhus University, Ny Munkegade 120, Aarhus C 8000,
Denmark

###### Abstract

The use of machine learning (ML) in chemical physics has enabled the
construction of interatomic potentials having the accuracy of ab initio
methods and a computational cost comparable to that of classical force fields.
Training an ML model requires an efficient method for the generation of
training data. Here we apply an accurate and efficient protocol to collect
training data for constructing a neural network based ML interatomic potential
for nanosilicate clusters. Initial training data are taken from normal modes
and farthest point sampling. Later on, the set of training data is extended
via an active learning strategy in which new data are identified by the
disagreement between an ensemble of ML models. The whole process is further
accelerated by parallel sampling over structures. We use the ML model to run
molecular dynamics (MD) simulations of nanosilicate clusters with various
sizes, from which infrared spectra with anharmonicity included can be
extracted. Such spectroscopic data are needed for understanding the properties
of silicate dust grains in the interstellar medium (ISM) and in circumstellar
environments.

##  I Introduction

Silicates are the main constituent of dust grains in the interstellar medium
(ISM) Henning (2010) and in circumstellar environments Tielens  _et al._
(1997). They provide surfaces for chemical reactions and nucleation sites for
the formation of ice mantles Potapov and McCoustra (2021). Silicate dust with
different compositions, sizes and shapes exhibit varied properties Min  _et
al._ (2007) that can be probed by infrared (IR) observations Kemper, Vriend,
and Tielens (2004); Spoon  _et al._ (2022). For example, a broad band at 9.7
$\mu$m is normally assigned to Si-O stretching modes and a band around 18
$\mu$m can be attributed to O-Si-O bending modes. Using an astrophysical model
to describe IR emission from carbonaceous and silicate dust grains Draine and
Li (2001); Li and Draine (2001a), it has been estimated that around 10% of
interstellar Si could be locked up in ultrasmall silicate grains with radii
$<$ 1.5 nm Li and Draine (2001b). Therefore, having efficient and accurate
methods to model the structures and vibrational properties of nanosilicate
clusters is beneficial for understanding IR features and other properties of
these, presumably very abundant, species.

Most previously reported spectroscopic data of silicate clusters are from
classical force field modelings Zamirri  _et al._ (2019) and quantum chemical
calculations Escatllar  _et al._ (2019); Zamirri  _et al._ (2019); Guiu,
Escatllar, and Bromley (2021). Although producing silicate clusters in
experiments is difficult Sabri  _et al._ (2013), recent cluster beam studies
have produced small silicate clusters and verified their structures via their
theoretically calculated IR spectra Mariñoso Guiu  _et al._ (2022). Many
theoretical studies of silicates are based on density functional theory (DFT)
calculations, whereby the vibrational spectra assume the harmonic
approximation. However, at finite temperature, an anharmonic treatment becomes
necessary due to a number of possible associated effects (e.g. thermal peak
broadening, frequency redshifting, combination bands, conformational
fluxionality) Guiu, Escatllar, and Bromley (2021). In this case, ab initio
molecular dynamics (AIMD) is a more appropriate method to compute IR spectra
with anharmonic and temperature effects included Guiu, Escatllar, and Bromley
(2021). Hybrid functionals are often used in IR calculations of nanosilicate
clusters in order to provide sufficient accuracy Guiu, Escatllar, and Bromley
(2021). However, the high computational costs of AIMD simulations with hybrid
functionals hinder their application for large cluster sizes. Computationally
inexpensive interatomic potentials are needed for MD simulations with large
system sizes. A number of such potentials have been developed and
parameterized classical force fields thus exist, that do describe bulk
silicate Price, Parker, and Leslie (1987); Walker, Wright, and Slater (2003)
and silicate clusters Escatllar  _et al._ (2019) quite well. They are often
used to roughly predict structure and energetic properties. For example, a
force field optimized for silicate clusters has been used for global
optimization of silicate clusters with different stoichiometries and sizes
Escatllar  _et al._ (2019). However, classical force fields are sometimes
limited by their accuracy. The low-energy isomers predicted by force fields
are thus typically refined by a more accurate method like DFT to ensure an
accurate energetic ordering is obtained.

Machine learning (ML) techniques such as artificial neural network Behler and
Parrinello (2007); Smith, Isayev, and Roitberg (2017); Schütt  _et al._
(2018); Unke and Meuwly (2019); Zaverkin and Kästner (2020); Zaverkin  _et
al._ (2021) and gaussian process regression Bartók  _et al._ (2010); Koistinen
_et al._ (2017); Denzel and Kästner (2018) are powerful tools to fit
interatomic potentials Behler (2021); Deringer  _et al._ (2021a) while
balancing accuracy and computational costs. They have been successfully
applied to global optimization Ouyang, Xie, and Jiang (2015); Chiriki, Jindal,
and Bulusu (2017); Bisbo and Hammer (2020); Paleico and Behler (2020a, b);
Timmermann  _et al._ (2020); Kaappa, del Río, and Jacobsen (2021); Bisbo and
Hammer (2022); Christiansen, Rønne, and Hammer (2022); Rønne  _et al._ (2022),
molecular dynamics Li, Kermode, and De Vita (2015); Caro  _et al._ (2018,
2020); Lim  _et al._ (2020); Noé  _et al._ (2020); Deringer  _et al._ (2021b);
Westermayr  _et al._ (2021); Unke  _et al._ (2021); Musaelian  _et al._ (2023)
and vibrational spectroscopy Gastegger, Behler, and Marquetand (2017); Lam,
Abdul-Al, and Allouche (2020); Gastegger, T. Schütt, and Müller (2021); Käser
_et al._ (2021); Beckmann  _et al._ (2022). In this work, we will develop
neural network based potentials for silicate clusters. The desired
applications are MD simulations of silicate clusters for variable
compositions, sizes and temperatures, from which IR spectra can be extracted.
The paper is organized as follows. We will first summarize the essential
elements to develop an accurate ML potential for three target properties
(total energy, atomic forces and total dipole moment) of silicate clusters.
This includes our choice of neural network architectures and the generation of
training data. We move on to a discussion of how the ML potential is trained
and validated on the three target properties. Its spectroscopic predictions
are then examined with respect to both harmonic and MD-based spectra.
Afterwards, the transferability of the ML potential is tested on the silicate
clusters that are not included in the training data. Finally, the ML potential
is used in MD simulations of silicate clusters and the astrophysical
implication of the obtained IR spectra is discussed.

##  II Methods

In order to compute IR spectra via an MD-based approach, it is needed to have
an accurate description of total energy, atomic forces and total dipole
moment. Total energy and atomic forces are trained in one model since atomic
forces are the negative gradient of the total energy and the inclusion of
forces in training often improves model accuracy for both energy and forces
Le, Huynh, and Raff (2009); Bartók  _et al._ (2010); Artrith and Behler
(2012); Christensen and von Lilienfeld (2020). The total dipole moment is
considered as a property which is independent from the energy and the forces.
As such, it is therefore trained in a separate model following Gastegger et al
Gastegger, Behler, and Marquetand (2017). Our choice of neural network models
is described in details in Section II.1. One of the challenges in training
such a model is to make sure the trained model is accurate and robust during a
long time MD simulations (needed for spectroscopic accuracy). The trained
model should also be applicable for MD simulations at high temperature (e.g.
800K) which are needed for mimicking some circumstellar conditions Guiu,
Escatllar, and Bromley (2021); Heese  _et al._ (2017). To address this issue,
we will use an active learning strategy Kolsbjerg, Peterson, and Hammer
(2018); van der Oord  _et al._ (2022) to train the ML model. The model gets
iteratively improved. During each iteration, the model is used in an MD
sampling for exploring the configuration space and selecting training data
Behler (2015); Gastegger, Behler, and Marquetand (2017). Both low and high
temperatures are chosen in the MD sampling so that different regions of the
configuration space are covered. More details can be found in Section II.2.
Another challenge is the selection of training data which is expected to be
representative and unique. This could also be addressed by the active
learning. At each iteration, only uncertain data that the current model cannot
predict well will be added to the training data Kolsbjerg, Peterson, and
Hammer (2018); Gastegger, Behler, and Marquetand (2017). More discussions will
be given in Section II.2. Finally, MD details will be described in Section
II.3 and DFT settings will be summarized in Section II.4.

###  II.1 Neural network models

The Gaussian moment neural network (GM-NN) Zaverkin and Kästner (2020);
Zaverkin  _et al._ (2021) developed by Kästner’s group was used to construct
potentials for total energies and atomic forces. The total energy of a given
system $\hat{\text{E}}$ is calculated as a sum of atomic energies
$\epsilon_{i}$.

$\hat{\text{E}}=\sum_{i=1}^{\text{N}_{\text{atoms}}}\epsilon_{i}$ (1)

Those atomic energies are the outputs from GM-NN models which take the local
atomic environments (learnable Gaussian moments in this model) as inputs. The
prediction of atomic force $\hat{\textbf{F}}_{i}$ on an atom $i$ is achieved
by taking the partial derivative of the total energy with respect to the
atomic position $\textbf{r}_{i}$ of atom $i$.

$\hat{\textbf{F}}_{i}=-\dfrac{\partial\hat{\text{E}}}{\partial\textbf{r}_{i}}$
(2)

During the training, a combined loss function is used.

$\mathcal{L}_{\text{E,F}}=\sum_{n=1}^{N_{\text{structures}}}\biggl{[}\lambda_{\text{E}}||\text{E}^{\text{ref}}-\hat{\text{E}}||^{2}+\frac{\lambda_{\text{F}}}{3\text{N}_{\text{atoms}}}\sum_{i=1}^{\text{N}_{\text{atoms}}}||\textbf{F}^{\text{ref}}_{i}-\hat{\textbf{F}}_{i}||^{2}\biggr{]}$
(3)

where the trade-off between energy and force is set to $\lambda_{\text{E}}=1$
and $\lambda_{\text{F}}=$ 10 Å2. The GM-NN code Zaverkin, Holzmüller, and
Kästner (2021) was used to train total energies and atomic forces.

The dipole moment was trained separately since the GM-NN code Zaverkin,
Holzmüller, and Kästner (2021) cannot provide dipole moment predictions yet.
We used the SchNet neural network Schütt  _et al._ (2018) implemented in
SchNetPack v1.0.0 Schütt  _et al._ (2019) to model the dipole moment. The
total electric dipole moment is computed as a sum of atomic charges weighted
by their positions with respect to the center of mass.

$\bm{\hat{\mu}}=\sum_{i=1}^{\text{N}_{\text{atoms}}}\hat{q}_{i}(\textbf{r}_{i}-\textbf{r}_{0})$
(4)

where $\hat{q}_{i}$ is the charge of atom $i$ predicted by a SchNet model and
$\textbf{r}_{0}$ is the center of mass for a given structure. Note that atomic
charges are not fixed for each individual atom, but depend on their atomic
environments. Consequently, atomic charges are fluctuating in MD simulations.
The loss function for dipole moment is defined as

$\mathcal{L}_{\mu}=\sum_{n=1}^{N_{\text{structures}}}\biggl{[}\frac{1}{3}||\bm{\mu}^{\text{ref}}-\hat{\bm{\mu}}||^{2}\biggr{]}$
(5)

###  II.2 Training data

The silicate clusters investigated in this work are taken from global minimum
energy structures of pyroxene (MgSiO3)${}_{\text{N}}$ and olivine
(Mg2SiO4)${}_{\text{N}}$ with a range of cluster sizes (N=1-10) Escatllar  _et
al._ (2019). The name of those silicate clusters is abbreviated as P1, P2, …
and P10 for pyroxene (MgSiO3)${}_{\text{N}}$ (N=1-10). Similarly, O1, O2, …
and O10 stand for olivine (Mg2SiO4)${}_{\text{N}}$ (N=1-10). Their structures
are shown in Fig. S1. During the construction of ML interatomic potentials,
the generation of training data is an important factor next to choosing
appropriate ML models. The ideal training data should be representative and
not redundant in the configuration space. We have used an active learning
strategy to generate our training data Gastegger, Behler, and Marquetand
(2017); Kolsbjerg, Peterson, and Hammer (2018); Gastegger and Marquetand
(2020). The active learning has been widely used to generate training data
while retaining small amounts of training data Smith  _et al._ (2018);
Bernstein, Csányi, and Deringer (2019); Schran, Brezina, and Marsalek (2020);
Vandermause  _et al._ (2020); Miksch  _et al._ (2021); Xie  _et al._ (2021);
Young  _et al._ (2021).

The first step is to make an initial set of training data. Initial training
data are generated by normal mode sampling Smith, Isayev, and Roitberg (2017);
Braams and Bowman (2009); Qu, Yu, and Bowman (2018) with a modification. Each
normal mode represents a unique way to displace a structure from its local
minimum as seen in Fig. 1. All displaced structures are further filtered by
farthest point sampling (FPS) Bartók  _et al._ (2017); Imbalzano  _et al._
(2018) in order to keep the size of training data small. For each silicate
cluster, one local minimum structure is selected along with
$3\times\text{N}_{\text{atoms}}$ structures from our modified normal mode
sampling. In the end, 2000 structures collected from all silicate systems are
used to train an initial ML model. The whole training data set is split into a
training subset and a validation subset with a ratio of 9:1 during initial and
subsequent trainings.

Later on, the ML model is iteratively improved by adding training data from
the active learning following the scheme in Fig. 2. The sampling is driven by
MD using the ML model and is paused upon the appearance of uncertain data that
the current ML model cannot describe well. A straightforward approach is to
run a DFT calculation at each MD step and compare with the model prediction.
Uncertain data can be easily located with this method which is, however,
computationally expensive. We use the query by committee method Smith  _et
al._ (2018); Schran, Brezina, and Marsalek (2020) which is a common approach
to get an uncertain estimate without heavy DFT calculations. It uses an
ensemble of ML models and finds uncertain data where different ML models
disagree with each other regarding the energy prediction Schran, Behler, and
Marx (2019); Schran, Brezina, and Marsalek (2020). The uncertain data is
recomputed at the DFT level and added to the training data. A new ML model is
trained after the collection of uncertain data. The active learning then moves
to next iteration in which MD sampling is performed again. The MD trajectory
from previous iteration is used as the starting geometry for the following
iteration. For each new iteration, the distribution of atomic velocities is
also taken from the previous iteration, and the overall translation and
rotation are removed. The target simulation temperature remains unchanged from
one iteration to another.

In order to further speed up the sampling, multiple systems are sampled in
parallel at different temperatures and stoichiometries Gastegger and
Marquetand (2020). When samplings are finished, DFT calculations on uncertain
data can be parallelized as well. Therefore, multiple structures (from one to
twenty) are added to the training data in one iteration. Only one job of ML
training in one iteration is conducted using the updated training data.

In summary, our protocol for collecting training data was composed of
iterations with these steps as shown in Fig. 2:

  1. 1.

Parallelize over structures

    1. (a)

MD propagation until the appearance of an uncertain data or reaching the
maximum simulation time

    2. (b)

Store uncertain data if appeared

  2. 2.

Wait until all MD runs have ended

  3. 3.

Perform DFT calculations on uncertain data in a massively parallel manner

  4. 4.

Train the ML model

The collection of training data will stop when the accuracy of the ML model
cannot be improved by adding more training data. In practice, we assess the
difference between DFT and ML in terms of harmonic frequencies of all
investigated silicate clusters. We stop collecting training data when this
difference cannot be decreased within a number of iterations. Finally, 1458
instances of uncertain data are collected and the final training data consists
of 3458 structures.

Figure 1:  Farthest point sampling on displaced structures along normal modes.
(a) Principal component analysis (PCA) on the global features of various P2
structures. The filled circles in the same color are the structures along one
specific mode and different color indicates different modes. The red crosses
are the selected points from farthest point sampling. (b) - (g) Snapshots of
the selected structures. The normal vectors are plotted in blue arrows. One
Si-O bond length is labeled to quantify the difference between those six
structures.  Figure 2:  Scheme of the active learning in a parallel manner. In
each sampling panel, blue and orange lines are energy predictions from two
independent ML models trained on the same data. Dashed black lines are the
target DFT potential. Each vertical red bar indicates the appearance of an
uncertain data.

###  II.3 MD simulations

For all MD simulations, a time step of 0.5 fs was used. During the active
learning, the temperature was controlled by a Berendsen thermostat Berendsen
_et al._ (1984) with a time constant of 100 fs implemented in ASE 3.21.1
Larsen  _et al._ (2017). The maximum simulation time was set to 5 ps during
the active learning. During the production runs for the IR spectra, the Nosé-
Hoover chain thermostat Martyna, Klein, and Tuckerman (1992) in SchNetPack
1.0.0 Schütt  _et al._ (2019) was used to keep a constant temperature. The
time constant was 100 fs and a chain length of 3 was used. The total
simulation time was set to 40 ps. The first 5 ps was discarded when extracting
IR spectra from the autocorrelation function of the dipole time derivative
Thomas  _et al._ (2013) according to Eq. 6.

$I_{IR}\propto\int_{-\infty}^{+\infty}\langle\dot{\mu}(\tau)\dot{\mu}(\tau+t)\rangle_{\tau}e^{-i\omega
t}dt$ (6)

All autocorrelation functions were computed using the Wiener-Khinchin theorem
Wiener (1930). In order to obtain IR spectra with high quality, a Hann window
function Blackman and Tukey (1958) and zero-padding were applied to the
autocorrelation functions before the Fourier transform. A maximum correlation
depth of 2048 fs was used. All processing of IR spectra was done with
SchNetPack 1.0.0 Schütt  _et al._ (2019).

###  II.4 DFT simulations

All DFT calculations were performed with ORCA 5.0.0 Neese  _et al._ (2020). A
hybrid functional PBE0 Adamo and Barone (1999) was used due to its excellent
accuracy in modeling structural and spectroscopic properties of silicate
clusters Guiu, Escatllar, and Bromley (2021). A double-zeta basis set def2-SVP
Weigend and Ahlrichs (2005) and an auxiliary basis def2/J Weigend (2006) were
used. All hybrid DFT calculations were accelerated by the RIJCOSX
approximation Neese  _et al._ (2009). Here, the RI-J approximation is the
resolution of identity (RI) approximation for Coulomb integrals (J) Neese
(2003). The ”chain of spheres” COSX approximation is used to speed up HF
exchange integrals Neese  _et al._ (2009). For all DFT calculations, the SCF
convergence criteria were set to tight. Default DFT and COSX grids ”DEFGRID2”
were used, and normally show small errors for energies, geometries and
frequencies compared with a larger grid setting. Examples of ORCA inputs for
single point calculations, structure optimization and harmonic calculations
are given in the supplementary material.

##  III Results

###  III.1 Validation of the ML model

To assess the quality of the trained ML model during the active learning, the
model was stored at the end of each iteration and tested against harmonic
frequencies of all investigated silicate clusters. Harmonic frequencies are
chosen since the reference frequencies from DFT calculations are already
available before the active learning and it is computationally cheap to
compute harmonic frequencies with the ML potential. The mean absolute error
(MAE) and root mean squared errors (RMSE) are used to quantify the difference
between the results from DFT and ML. Fig. 3 (a) shows MAE/RMSE of harmonic
frequency at each sampling iteration. Before the active learning, when the
training data are only sampled from normal modes, the MAE of harmonic
frequencies is 3.2 cm-1. This error decreases with the use of active learning,
indicating the improvement of the ML model when more uncertain data are
included in the training. After about 40 iterations, the MAE becomes stagnated
and the ML model cannot be further improved by adding more data according to
our sampling strategy. Therefore, the active learning is stopped with a
1.6/2.2 cm-1 MAE/RMSE of harmonic frequency. We also need to see if the final
ML model can be reliably used in MD simulations. To assess this, another set
of test data other than harmonic frequencies is constructed based on the MD
simulations using the final ML model. The MD simulations are performed for 100
ps for each silicate cluster at three different temperatures (400K, 800K and
1200K). No obscure configurations (e.g. overlapping atoms and atoms being far
from each other) were observed in the trajectory of any above-mentioned MD
run. DFT calculations at the PBE0/def2-SVP level are performed on 100 randomly
selected structures from each MD run. Finally, a test data set consisting of
6000 structures was collected. To assess the degree of improvement in terms of
energy and force predictions, energies and forces are recomputed on the test
data using the ML model at each iteration. Their corresponding MAEs and RMSEs
are shown in Fig. 3 (b) and Fig. 3 (c), and decrease as the iteration
proceeds. The final ML model gives an MAE/RMSE of 2.3/3.5 meV for energy per
atom and 0.043/0.060 eV/Å for atomic force, based on the test data set. The
accuracy of the final ML model on the training and validation data set are
included in Table 1.

A comparison between the ML model and classical force fields is given in Fig.
S2 in terms of energy and force predictions. For the energy comparison, we use
relative energies with respect to the energy of global minimum structure,
divided by the number of atoms. The force field by Escatllar et al Escatllar
_et al._ (2019) is chosen since it is parameterized with respect to relative
energies and cluster geometries. Our ML model shows a one order of magnitude
improvement for energy prediction over the Escatllar force field Escatllar
_et al._ (2019). For the force prediction, we choose the force field by Walker
et al Walker, Wright, and Slater (2003), as it is designed to study oxygen
diffusion in olivine and more suitable for force prediction than the Escatllar
force field Escatllar  _et al._ (2019). Our ML model shows a two order of
magnitude improvement for force prediction over the Walker force field Walker,
Wright, and Slater (2003).

Figure 3:  Accuracy of ML models at each iteration of active learning, in terms of mean absolute errors (MAEs) and root mean squared errors (RMSEs) of (a) harmonic frequency, (b) energy per atom, and (c) atomic force. Harmonic frequency is based on 20 silicate clusters (P1-P10 and O1-O10). Energy and forces are based on the test data. Note that the test data are not included in the training and are used to assess the accuracy of ML models. The generation of test data is described in details in Section III.1.  Table 1: Accuracy of the final ML model on training, validation and test data sets Property | Unit | Model | Training set | Validation set | Test set
---|---|---|---|---|---
MAE | RMSE | MAE | RMSE | MAE | RMSE
Energy per atom | meV | GM-NN | 1.6 | 3.2 | 1.7 | 3.0 | 2.3 | 3.5
Atomic forces | ev/Å | GM-NN | 0.037 | 0.050 | 0.045 | 0.077 | 0.043 | 0.060
Dipole moment | Debye | SchNet | 0.049 | 0.064 | 0.090 | 0.128 | 0.074 | 0.102

In order to compute infrared intensity, the prediction of dipole moment is
needed. We have trained SchNet models only for dipole moment based on two set
of training data, namely the initial dataset and the final dataset after the
active learning is finished. Fig. 4 compares the reference DFT dipole moment
against the predicted dipole moment for the test data when different training
data are used. The ML model shows good accuracy in predicting dipole moment
with MAE of 0.142 Debye and RMSE of 0.200 Debye when only structures sampled
from normal modes are included in the training. The ML model from the active
learning results in a better accuracy in predicting dipole moment and the
MAE/RMSE is only 0.074/0.102 Debye. Adding training data by the active
learning indeed improve the ML model’s accuracy regarding dipole moment
prediction.

Figure 4: Accuracy of ML models in predicting dipole moment for the test data
when different training data are used. Note that the test data are not
included in the training and are used to assess the accuracy of ML models. The
generation of test data is described in details in Section III.1.

###  III.2 Spectroscopic behavior of the ML model

Our ML models have reproduced the potential energy surface (PES) and dipole
surface of DFT to a reasonable degree of accuracy. The next step is to
benchmark the spectroscopic behavior of ML models against DFT results. IR
spectra containing both band positions and infrared intensities at the
harmonic level are evaluated first. Infrared intensities, which are not
examined in Section III.1 yet, are calculated from dipole derivatives along
normal modes. The DFT-calculated band positions and intensities are reproduced
almost exactly by ML models for all investigated silicate clusters, see Fig.
S3. MAEs of harmonic frequencies for most clusters are below 2.5 cm-1. P1 and
O1 show slightly higher MAEs of 4.2 cm-1 and 3.6 cm-1. For harmonic
intensities, the averaged MAE for all clusters is 26.8 km/mol. P2 shows the
lowest MAE of 2.8 km/mol and P10 shows the highest MAE of 49.4 km/mol.

Selected harmonic spectra are given in the bottom panel of Fig. 5 along with
MD-based spectra at three selected temperatures (100K, 400K and 800K). For MD-
based spectra, band positions are well reproduced by ML models while
intensities have larger variations than band positions. The variations in
intensities could arise from differences in velocity initializations in the MD
runs or inaccuracy of the ML models. In order to check the influence of
variable velocity initializations, we performed twenty five independent DFT-MD
runs for P2 at 400K. The averaged IR spectra over independent DFT-MD runs are
compared with those of ML-MD runs in Fig. 6. When only one MD run is used,
intensity variations are observed for peaks with frequencies around 464, 748,
827, and 1114 cm-1. The variation becomes smaller for bands at 464, 827 and
1114 cm-1 when the number of independent MD runs increased from one to five,
and barely changes when twenty five independent runs are used. The
disagreement between DFT and ML results of band position and intensity at 748
cm-1 (Si-O stretching modes between the central O atoms and Si atoms) is not
sensitive to the number of independent MD runs. This disagreement is therefore
likely due to model fitting. Performing independent DFT-MD runs for each
silicate size and temperature is computationally expensive. Therefore only one
DFT-MD run and one ML-MD run are compared in Fig. 5 for each system. In Fig.
S4, one DFT-MD run is compared against independent ML-MD runs. In general, DFT
results are within the range of independent ML runs. We expect better
agreement between DFT and ML when more independent runs are performed. Since
it is computationally expensive to run many DFT-MD simulations, the resulting
IR spectra are expected to be biased by the velocity initializations. However,
the use of an ML potential allows for a fast evaluation of many independent MD
runs and thus the generation of converged IR spectra.

Figure 5:  Comparing IR spectra from DFT and the ML model. Both harmonic and
MD-based spectra are given for selected silicate clusters (a) O2, (b) P2, (c)
P5 at 100K, 400K and 800K. Figure 6:  MD-based IR spectra of P2 at 400K
averaged over independent runs. DFT results are shown in solid black lines and
ML results are shown in colored dashed lines. Number of independent DFT-MD and
ML-MD runs is labeled on the right of each spectrum. Bands at 464, 748, 827,
and 1114 cm-1 are highlighted by vertical gray lines.

###  III.3 Transferability of the ML model

After the assessment of the ML model, we want to see whether the ML model is
accurate and transferable for the structures that are different from the
global minimum based training and test data. Examples of such structures are
high-energy isomers that could be obtained during global optimization searches
via the GOFEE (global optimization with first-principles energy expressions)
method Bisbo and Hammer (2022). At first, P5 high-energy isomers are chosen to
test the transferability of the ML model. Only harmonic spectra are used to
verify the transferability of the ML model, since harmonic frequencies are
sensitive to the quality of the PES and harmonic calculations at the DFT level
are computationally much cheaper than DFT-MD simulations. In addition,
harmonic frequencies are the main reservoir of vibrational transitions for a
given nanosilicate cluster. The synthetic IR spectra of various nanosilicate
clusters based on the harmonic approximation are already used to estimate the
abundance of nanosilicate clusters in the diffuse ISM Zeegers  _et al._ (2023)
and are valuable for understanding the property of interstellar dust. We
compute harmonic IR spectra of these isomers with DFT (PBE0/def2-SVP) and the
ML model and show their spectra comparisons in Fig. 7 (a). All isomers have
MAEs of harmonic frequencies lower than 4 cm-1 and infrared intensities are
well fitted. The good transferability of the ML model is expected when the
target system is similar with training data in the configuration space even
though they may look quite different in the Cartesian coordinate space. The
configuration space for the nanosilicates investigated in this work has high
dimensions and is hard to directly visualize. Therefore, principal component
analysis (PCA) is used to reduce the dimension of the configuration space with
each configuration represented by a global Smooth Overlap of Atomic Positions
(SOAP) feature Bartók, Kondor, and Csányi (2013). Among the four selected
isomers, isomer-I has the lowest MAE of harmonic frequency since it lies
closer to training data in the feature space than other isomers in Fig. 7 (b).
Similarly, isomer-II shows the highest MAE of harmonic frequency because it
has the lowest similarity to training data in the feature space. Later on,
transferability analysis is also performed for other pyroxene (P3-P8) and
olivine (O2-O7) isomers (see their structure snapshots in Fig. S6 and Fig. S7)
which are not directly included in the training data. The results are shown in
Fig. S5. The MAE of harmonic frequency is 2.1 cm-1 for pyroxene and 1.7 cm-1
for olivine. They indicate that the ML model tends to show good prediction
accuracy on the test data which looks similar with training data in the
feature space. We do not expect that our ML model will be reliable and
accurate for bulk silicates and chemical processes involving dramatic changes
in bonding types (e.g. oxygen diffusion, silicate growth and non-
stoichiometric silicates). The transferability to larger cluster sizes will be
examined in a separate study.

Figure 7:  Transferability of the ML model on P5 isomers found via the GOFEE
method Bisbo and Hammer (2022). (a) Harmonic IR spectra of selected P5 isomers
with DFT results in black solid lines and ML results in colored dashed lines.
MAEs of harmonic frequencies are labeled for each isomer. (b) Principal
component analysis (PCA) on global Smooth Overlap of Atomic Positions (SOAP)
features Bartók, Kondor, and Csányi (2013) of all pyroxene clusters. Gray
circles are training data. Blue circles are P5 training data. Colored crosses
are high-energy P5 isomers with their structures and relative energies to P5
global minimum attached.

###  III.4 Application of the ML model

As our motivation for developing a ML potential of nanosilicate clusters is to
aid the interpretation of silicate band in infrared observations, we compute
IR spectra of 20 silicate clusters (P1-P10 and O1-O10) via MD simulations with
the ML model. Fig. 8 shows the sum of IR spectra for two types (pyroxene and
olivine) of silicate clusters while assuming the distribution of cluster size
is uniform. For both pyroxene and olivine clusters, the most intense IR peaks
are in the 9-10 $\mu$m range. They both have a main peak with a slightly
smaller wavelength than the signature 9.7 $\mu$m feature at low temperature.
Our MD-based IR spectra at low temperature are in a good agreement with a
previous IR study of pyroxene and olivine clusters under the harmonic
approximation Escatllar  _et al._ (2019). As temperature increases, this
blueshifted peak broadens and becomes closer to 9.7 $\mu$m. In addition, the
center position of IR peaks of olivine is closer to 9.7 $\mu$m than that of
pyroxene. More systematic IR investigations of nanosilicate clusters and their
astrophysical implications are planned in the future.

Figure 8:  Sum of IR spectra of (a) pyroxene (P1-P10) (b) olivine (O1-O10) at
100K, 400K and 800K using the ML model, assuming equal weights of different
cluster sizes. A vertical black line highlights the 9.7 $\mu$m silicate band.

##  IV Conclusions

We have demonstrated that active learning is an efficient method to generate
training data for constructing ML interatomic potentials. After initial
training, our ML potential is iteratively improved employing active learning,
which can be further accelerated in a parallel manner. Our comparison of ML-
based and DFT-based infrared spectra (Fig. 5) reveals that our ML potential is
highly accurate for infrared spectroscopic studies of nanosilicate clusters
with both harmonic approximation and an MD-based approach. Our ML potential
also exhibits certain transferability regarding harmonic frequencies and IR
intensities for high-energy isomers of nanosilicate clusters that are not
directly included in the training data, but have similar cluster sizes as the
training data. Since the ML potential is computationally much cheaper than
DFT, longer MD runs can be conducted and better statistics can be achieved. It
remains an open question if our ML potential is still accurate and reliable
for larger silicate cluster sizes that are not covered in this study. We will
address this issue in a future study. Our work will allow for systemic
investigations into how the IR spectra of silicate clusters depend on size,
structure and temperature. Our work will be particularly beneficial for
understanding the properties of silicate dust grains in space, especially in
connection with the highly sensitive IR capabilities of the James Webb Space
Telescope Zeegers  _et al._ (2021).

##  V Supplementary Material

See the supplementary material for examples of ORCA input files, silicate
structures that are used in this work and more validations of the ML model.

###### Acknowledgements.

We thank Andreas Møller Slavensky for providing the structure files of high-
energy silicate isomers. This work has been supported by the Danish National
Research Foundation through the Center of Excellence “InterCat” (Grant
agreement no.: DNRF150) and the VILLUM FONDEN (Investigator grant, Project No.
16562). S.T.B acknowledges support from the MICINN funded project grants
PID2021-127957NB-I00 and TED2021-132550B-C21, project grant 2021SGR00354
funded by the Generalitat de Catalunya and through the María de Maeztu program
for Spanish Structures of Excellence (CEX2021-001202-M).

##  VI Conflict of Interest

The authors have no conflicts to disclose.

##  VII DATA AVAILABILITY

The ML model developed in this work and source data for figures in this
article are available in https://github.com/zyt0y/MLP-IR-nanosilicate-
clusters. Other data that support the findings of this work are available
within the article and its supplementary material.

## References

## References

  * Henning (2010) T. Henning, Annu. Rev. Astron. Astrophys. 48, 21 (2010).
  * Tielens  _et al._ (1997) A. Tielens, L. Waters, F. Molster, and K. Justtanont, Astrophys. Space Sci. 255, 415 (1997).
  * Potapov and McCoustra (2021) A. Potapov and M. McCoustra, Int Rev Phys Chem 40, 299 (2021).
  * Min  _et al._ (2007) M. Min, L. B. F. M. Waters, A. de Koter, J. W. Hovenier, L. P. Keller, and F. Markwick-Kemper, A&A 462, 667 (2007).
  * Kemper, Vriend, and Tielens (2004) F. Kemper, W. J. Vriend, and A. G. G. M. Tielens, ApJ 609, 826 (2004).
  * Spoon  _et al._ (2022) H. W. W. Spoon, A. Hernán-Caballero, D. Rupke, L. B. F. M. Waters, V. Lebouteiller, A. G. G. M. Tielens, T. Loredo, Y. Su, and V. Viola, ApJS 259, 37 (2022).
  * Draine and Li (2001) B. T. Draine and A. Li, ApJ 551, 807 (2001).
  * Li and Draine (2001a) A. Li and B. T. Draine, ApJ 554, 778 (2001a).
  * Li and Draine (2001b) A. Li and B. T. Draine, ApJ 550, L213 (2001b).
  * Zamirri  _et al._ (2019) L. Zamirri, A. Macià Escatllar, J. Mariñoso Guiu, P. Ugliengo, and S. T. Bromley, ACS Earth Space Chem. 3, 2323 (2019).
  * Escatllar  _et al._ (2019) A. M. Escatllar, T. Lazaukas, S. M. Woodley, and S. T. Bromley, ACS Earth Space Chem. 3, 2390 (2019).
  * Guiu, Escatllar, and Bromley (2021) J. M. Guiu, A. M. Escatllar, and S. T. Bromley, ACS Earth Space Chem. 5, 812 (2021).
  * Sabri  _et al._ (2013) T. Sabri, L. Gavilan, C. Jäger, J. L. Lemaire, G. Vidali, H. Mutschke, and T. Henning, ApJ 780, 180 (2013).
  * Mariñoso Guiu  _et al._ (2022) J. Mariñoso Guiu, B.-A. Ghejan, T. M. Bernhardt, J. M. Bakker, S. M. Lang, and S. T. Bromley, ACS Earth Space Chem. 6, 2465 (2022).
  * Price, Parker, and Leslie (1987) G. D. Price, S. C. Parker, and M. Leslie, Phys Chem Minerals 15, 181 (1987).
  * Walker, Wright, and Slater (2003) A. M. Walker, K. Wright, and B. Slater, Phys Chem Minerals 30, 536 (2003).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).
  * Smith, Isayev, and Roitberg (2017) J. S. Smith, O. Isayev, and A. E. Roitberg, Chem. Sci. 8, 3192 (2017).
  * Schütt  _et al._ (2018) K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, and K.-R. Müller, J. Chem. Phys. 148, 241722 (2018).
  * Unke and Meuwly (2019) O. T. Unke and M. Meuwly, J. Chem. Theory Comput. 15, 3678 (2019).
  * Zaverkin and Kästner (2020) V. Zaverkin and J. Kästner, J. Chem. Theory Comput. 16, 5410 (2020).
  * Zaverkin  _et al._ (2021) V. Zaverkin, D. Holzmüller, I. Steinwart, and J. Kästner, J. Chem. Theory Comput. 17, 6658 (2021).
  * Bartók  _et al._ (2010) A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Phys. Rev. Lett. 104, 136403 (2010).
  * Koistinen  _et al._ (2017) O.-P. Koistinen, F. B. Dagbjartsdóttir, V. Ásgeirsson, A. Vehtari, and H. Jónsson, J. Chem. Phys. 147, 152720 (2017).
  * Denzel and Kästner (2018) A. Denzel and J. Kästner, J. Chem. Phys. 148, 094114 (2018).
  * Behler (2021) J. Behler, Chem. Rev. 121, 10037 (2021).
  * Deringer  _et al._ (2021a) V. L. Deringer, A. P. Bartók, N. Bernstein, D. M. Wilkins, M. Ceriotti, and G. Csányi, Chem. Rev. 121, 10073 (2021a).
  * Ouyang, Xie, and Jiang (2015) R. Ouyang, Y. Xie, and D.-e. Jiang, Nanoscale 7, 14817 (2015).
  * Chiriki, Jindal, and Bulusu (2017) S. Chiriki, S. Jindal, and S. S. Bulusu, J. Chem. Phys. 146, 084314 (2017).
  * Bisbo and Hammer (2020) M. K. Bisbo and B. Hammer, Phys. Rev. Lett. 124, 086102 (2020).
  * Paleico and Behler (2020a) M. L. Paleico and J. Behler, J. Chem. Phys. 152, 094109 (2020a).
  * Paleico and Behler (2020b) M. L. Paleico and J. Behler, J. Chem. Phys. 153, 054704 (2020b).
  * Timmermann  _et al._ (2020) J. Timmermann, F. Kraushofer, N. Resch, P. Li, Y. Wang, Z. Mao, M. Riva, Y. Lee, C. Staacke, M. Schmid, C. Scheurer, G. S. Parkinson, U. Diebold, and K. Reuter, Phys. Rev. Lett. 125, 206101 (2020).
  * Kaappa, del Río, and Jacobsen (2021) S. Kaappa, E. G. del Río, and K. W. Jacobsen, Phys. Rev. B 103, 174114 (2021).
  * Bisbo and Hammer (2022) M. K. Bisbo and B. Hammer, Phys. Rev. B 105, 245404 (2022).
  * Christiansen, Rønne, and Hammer (2022) M.-P. V. Christiansen, N. Rønne, and B. Hammer, J. Chem. Phys. 157, 054701 (2022).
  * Rønne  _et al._ (2022) N. Rønne, M.-P. V. Christiansen, A. M. Slavensky, Z. Tang, F. Brix, M. E. Pedersen, M. K. Bisbo, and B. Hammer, J. Chem. Phys. 157, 174115 (2022).
  * Li, Kermode, and De Vita (2015) Z. Li, J. R. Kermode, and A. De Vita, Phys. Rev. Lett. 114, 096405 (2015).
  * Caro  _et al._ (2018) M. A. Caro, V. L. Deringer, J. Koskinen, T. Laurila, and G. Csányi, Phys. Rev. Lett. 120, 166101 (2018).
  * Caro  _et al._ (2020) M. A. Caro, G. Csányi, T. Laurila, and V. L. Deringer, Phys. Rev. B 102, 174201 (2020).
  * Lim  _et al._ (2020) J. S. Lim, J. Vandermause, M. A. van Spronsen, A. Musaelian, Y. Xie, L. Sun, C. R. O’Connor, T. Egle, N. Molinari, J. Florian, K. Duanmu, R. J. Madix, P. Sautet, C. M. Friend, and B. Kozinsky, J. Am. Chem. Soc. 142, 15907 (2020).
  * Noé  _et al._ (2020) F. Noé, A. Tkatchenko, K.-R. Müller, and C. Clementi, Annu. Rev. Phys. Chem. 71, 361 (2020).
  * Deringer  _et al._ (2021b) V. L. Deringer, N. Bernstein, G. Csányi, C. Ben Mahmoud, M. Ceriotti, M. Wilson, D. A. Drabold, and S. R. Elliott, Nature 589, 59 (2021b).
  * Westermayr  _et al._ (2021) J. Westermayr, M. Gastegger, K. T. Schütt, and R. J. Maurer, J. Chem. Phys. 154, 230903 (2021).
  * Unke  _et al._ (2021) O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko, and K.-R. Müller, Chem. Rev. 121, 10142 (2021).
  * Musaelian  _et al._ (2023) A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, Nat Commun 14, 579 (2023).
  * Gastegger, Behler, and Marquetand (2017) M. Gastegger, J. Behler, and P. Marquetand, Chem. Sci. 8, 6924 (2017).
  * Lam, Abdul-Al, and Allouche (2020) J. Lam, S. Abdul-Al, and A.-R. Allouche, J. Chem. Theory Comput. 16, 1681 (2020).
  * Gastegger, T. Schütt, and Müller (2021) M. Gastegger, K. T. Schütt, and K.-R. Müller, Chem. Sci. 12, 11473 (2021).
  * Käser  _et al._ (2021) S. Käser, E. D. Boittier, M. Upadhyay, and M. Meuwly, J. Chem. Theory Comput. 17, 3687 (2021).
  * Beckmann  _et al._ (2022) R. Beckmann, F. Brieuc, C. Schran, and D. Marx, J. Chem. Theory Comput. 18, 5492 (2022).
  * Le, Huynh, and Raff (2009) H. M. Le, S. Huynh, and L. M. Raff, J. Chem. Phys. 131, 014107 (2009).
  * Artrith and Behler (2012) N. Artrith and J. Behler, Phys. Rev. B 85, 045439 (2012).
  * Christensen and von Lilienfeld (2020) A. S. Christensen and O. A. von Lilienfeld, Mach. Learn.: Sci. Technol. 1, 045018 (2020).
  * Heese  _et al._ (2017) S. Heese, S. Wolf, A. Dutrey, and S. Guilloteau, A&A 604, A5 (2017).
  * Kolsbjerg, Peterson, and Hammer (2018) E. L. Kolsbjerg, A. A. Peterson, and B. Hammer, Phys. Rev. B 97, 195424 (2018).
  * van der Oord  _et al._ (2022) C. van der Oord, M. Sachs, D. P. Kovács, C. Ortner, and G. Csányi, “Hyperactive Learning (HAL) for Data-Driven Interatomic Potentials,”  (2022), arXiv:2210.04225 [physics, stat] .
  * Behler (2015) J. Behler, Int. J. Quantum Chem. 115, 1032 (2015).
  * Zaverkin, Holzmüller, and Kästner (2021) V. Zaverkin, D. Holzmüller, and J. Kästner, “The Gaussian Moment Neural Network Package,” https://gitlab.com/zaverkin_v/gmnn (accessed 21 November 2021).
  * Schütt  _et al._ (2019) K. T. Schütt, P. Kessel, M. Gastegger, K. A. Nicoli, A. Tkatchenko, and K.-R. Müller, J. Chem. Theory Comput. 15, 448 (2019).
  * Gastegger and Marquetand (2020) M. Gastegger and P. Marquetand, in  _Machine Learning Meets Quantum Physics_, edited by K. T. Schütt, S. Chmiela, O. A. von Lilienfeld, A. Tkatchenko, K. Tsuda, and K.-R. Müller (Springer International Publishing, Cham, 2020) pp. 233–252.
  * Smith  _et al._ (2018) J. S. Smith, B. Nebgen, N. Lubbers, O. Isayev, and A. E. Roitberg, J. Chem. Phys. 148, 241733 (2018).
  * Bernstein, Csányi, and Deringer (2019) N. Bernstein, G. Csányi, and V. L. Deringer, npj Comput Mater 5, 1 (2019).
  * Schran, Brezina, and Marsalek (2020) C. Schran, K. Brezina, and O. Marsalek, J. Chem. Phys. 153, 104105 (2020).
  * Vandermause  _et al._ (2020) J. Vandermause, S. B. Torrisi, S. Batzner, Y. Xie, L. Sun, A. M. Kolpak, and B. Kozinsky, Npj Comput. Mater. 6, 1 (2020).
  * Miksch  _et al._ (2021) A. M. Miksch, T. Morawietz, J. Kästner, A. Urban, and N. Artrith, Mach. Learn.: Sci. Technol. 2, 031001 (2021).
  * Xie  _et al._ (2021) Y. Xie, J. Vandermause, L. Sun, A. Cepellotti, and B. Kozinsky, npj Comput Mater 7, 1 (2021).
  * Young  _et al._ (2021) T. A. Young, T. Johnston-Wood, V. L. Deringer, and F. Duarte, Chem. Sci. 12, 10944 (2021).
  * Braams and Bowman (2009) B. J. Braams and J. M. Bowman, Int Rev Phys Chem 28, 577 (2009).
  * Qu, Yu, and Bowman (2018) C. Qu, Q. Yu, and J. M. Bowman, Annu. Rev. Phys. Chem. 69, 151 (2018).
  * Bartók  _et al._ (2017) A. P. Bartók, S. De, C. Poelking, N. Bernstein, J. R. Kermode, G. Csányi, and M. Ceriotti, Sci Adv 3, e1701816 (2017).
  * Imbalzano  _et al._ (2018) G. Imbalzano, A. Anelli, D. Giofré, S. Klees, J. Behler, and M. Ceriotti, J. Chem. Phys. 148, 241730 (2018).
  * Schran, Behler, and Marx (2019) C. Schran, J. Behler, and D. Marx, J. Chem. Theory Comput. (2019), 10.1021/acs.jctc.9b00805.
  * Berendsen  _et al._ (1984) H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984).
  * Larsen  _et al._ (2017) A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, Marcin Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, Kristen Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, Andrew Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, Tejs Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, and K. W. Jacobsen, J. Phys.: Condens. Matter 29, 273002 (2017).
  * Martyna, Klein, and Tuckerman (1992) G. J. Martyna, M. L. Klein, and M. Tuckerman, J. Chem. Phys. 97, 2635 (1992).
  * Thomas  _et al._ (2013) M. Thomas, M. Brehm, R. Fligg, P. Vöhringer, and B. Kirchner, Phys. Chem. Chem. Phys. 15, 6608 (2013).
  * Wiener (1930) N. Wiener, Acta Math. 55, 117 (1930).
  * Blackman and Tukey (1958) R. B. Blackman and J. W. Tukey, Bell Labs Tech. J. 37, 185 (1958).
  * Neese  _et al._ (2020) F. Neese, F. Wennmohs, U. Becker, and C. Riplinger, J. Chem. Phys. 152, 224108 (2020).
  * Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  * Weigend and Ahlrichs (2005) F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005).
  * Weigend (2006) F. Weigend, Phys. Chem. Chem. Phys. 8, 1057 (2006).
  * Neese  _et al._ (2009) F. Neese, F. Wennmohs, A. Hansen, and U. Becker, Chem. Phys. 356, 98 (2009).
  * Neese (2003) F. Neese, J. Comput. Chem. 24, 1740 (2003).
  * Zeegers  _et al._ (2023) S. Zeegers, J. M. Guiu, F. Kemper, J. Marshall, and S. T. Bromley, Faraday Discuss. (2023), 10.1039/D3FD00055A.
  * Bartók, Kondor, and Csányi (2013) A. P. Bartók, R. Kondor, and G. Csányi, Phys. Rev. B 87, 184115 (2013).
  * Zeegers  _et al._ (2021) S. Zeegers, J. Bouwman, J. Chiar, E. Costantini, M. Decleir, T. Dharmawardena, T. R. Geballe, K. D. Gordon, J. D. Green, T. K. Henning, O. Jones, F. Kemper, A. Li, J. P. Marshall, M. McClure, K. Misselt, Y. J. Pendleton, G. Perotti, K. M. Pontoppidan, A. Potapov, P. Scicluna, A. Tielens, R. Waters, and E. Zari, JWST Proposal. Cycle 1 , 2183 (2021).