[d]Aditya Pathak

# Precision Top Mass Measurement Using Energy Correlators

Jack Holguin  Ian Moult  Massimiliano Procura  Robert Schöfbeck  Dennis
Schwarz

###### Abstract

Precision top mass measurements at hadron colliders have been notoriously
difficult. The fundamental challenge in the current approaches lies in
achieving simultaneously high top mass sensitivity and good theoretical
control. Inspired by the use of standard candles in cosmology, we overcome
this problem by showing that a single energy correlator-based observable can
be constructed that reflects the characteristic angular scales associated with
both the $W$-boson and top quark. This gives direct access to the
dimensionless quantity $m_{t}/m_{W}$, from which $m_{t}$ can be extracted in a
well-defined short-distance mass scheme as a function of the well-known
$m_{W}$. A Monte-Carlo-based study is performed to demonstrate the properties
of our observable and the statistical feasibility of its extraction from the
Run 2 and 3 and High-Luminosity LHC data sets. The resulting $m_{t}$ has
remarkably small uncertainties from hadronization effects and is insensitive
to the underlying event and parton distribution functions. Our proposed
observable provides a road map for a rich program to achieve a top mass
determination at the LHC with record precision.

##  1 Introduction

Figure 1: Challenges in the current paradigm for top mass measurements.
Adapted from Ref. [1].

Precision measurement of the top mass is crucial for the stability analysis of
the electroweak vacuum, consistency tests of the SM, and indirect searches for
physics beyond the SM. The top quark mass $m_{t}$ is not a direct physical
observable but a Lagrangian parameter that receives quantum corrections that
must be renormalized and defined within a specific mass scheme. Figure 1
represents a generic top mass sensitive observable currently employed for top
mass measurement. Here, the highest mass sensitivity arises through a
threshold feature such as a peaked structure or a shoulder. The physics of the
threshold region is, in general, extremely challenging to describe as it is
dominated by several widely separated scales, multiple soft and collinear
emissions, and large nonperturbative corrections, leading to the breakdown of
fixed-order perturbation theory (as shown by the steeply rising green curve
labeled ‘NLO’ as we move towards the threshold structure). On the other hand,
fixed-order corrections predominantly impact the tail where the
$m_{t}$-sensitivity is relatively poorer.

The approach so far has been to describe the effects in the threshold region
using the Monte Carlo (MC) parton shower event generators interfaced with a
hadronization model. As shown in figure 1, the region of highest top mass
sensitivity is essentially described entirely by the parton shower, and NLO
corrections have little impact in this region (as shown by the closeness of
the ‘LO+Parton Shower’ and ‘NLO+ Parton Shower’ curves in the peak region).
This has led to MC-based extractions called the direct measurements with the
smallest quoted uncertainties, with the current world average being
$m_{t}^{\mathrm{MC}}=172.76\pm 0.3$ GeV [2]. The superscript ‘MC’ signifies
that these direct measurements are measuring the “Monte Carlo top mass
parameter” that does not have a straightforward correspondence to a field-
theoretic mass scheme definition. While MCs are indispensable tools for LHC
physics capable of describing arbitrarily complex observables, they are
limited in accuracy and typically only leading color, leading logarithmic
accurate. Consequently, limitations of parton showers and their complicated
interface with hadronization models have been argued to induce an additional
theory/conceptual uncertainty of $\mathcal{O}(1~{}{\,\mathrm{GeV}})$ in the
top mass [1], which is larger than the experimental uncertainty.

Next, measurements of multi-differential $t{\bar{t}}$ and $t{\bar{t}}+1$ jet
cross sections have yielded top mass measurements in the pole mass scheme, and
the current world average stands at $m_{t}^{\rm pole}=172.5\pm
0.7\,\mathrm{GeV}$ [2]. However, in accordance with the general features
described in figure 1, the highest sensitivity arises in the region of low
$t{\bar{t}}$ masses where large threshold corrections from Coulomb and soft
gluon resummation are expected to dominate. At present, they are only known
with large uncertainties and are estimated to induce additional sizable $\sim
1\,\mathrm{GeV}$ uncertainty in the $m_{t}^{\rm pole}$ measurements [3, 4] on
top of the quoted measurement. Finally, the use of MCs also induces the top
mass interpretation problem in various alternative top mass measurements that
directly rely on MCs. In conclusion, the precision measurement of the top mass
with existing approaches is notoriously challenging.

##  2 Energy-energy correlators on boosted top quarks

Figure 2: The three-point correlator measured in an equilateral configuration
shows the imprint of three-body decay of the top quark (left) at angles
$\zeta_{t}\sim m_{t}^{2}/Q^{2}$ corresponding to the top quark boost. The
peaked structure is highly sensitive to the top mass (center) and extremely
robust against hadronization (right) [5].

Correlation functions are the most fundamental tools for specifying the
ensemble properties of probabilistic systems and have a long history in QCD
collider physics. Energy-Energy Correlators (EECs) quantify the angular
correlations of hadronic energy flux within jets. While EECs were already
proposed several decades ago [6, 7, 8, 9], their unique potential for
precision collider physics has only just started to emerge [10]. Measured on a
jet, EECs are given by inclusive cross sections where only the relative angles
are fixed, and the energies are integrated over:

$\displaystyle\text{E}^{N}\text{C}(\zeta_{12},\zeta_{23},\ldots)\equiv\sum_{i,j,k\in
J}\int\textrm{d}\zeta_{ijk\ldots}\delta\big{(}\zeta_{12}-\zeta_{ij}\big{)}\delta(\zeta_{23}-\zeta_{kl})\ldots\frac{E_{i}E_{j}E_{k}\ldots}{Q^{N}}\frac{\textrm{d}\sigma_{i,j,k,\ldots}}{\textrm{d}\zeta_{ijk\ldots}}\,,$
(1)

where ENC refers to the $N$-point correlation function, and $Q$ is the large
momentum (such as jet $p_{T}$) associated with the jet and $E_{i}$’s (or
$p_{T,i}$’s) the energy of hadron $i$. As eq. (1) shows, they are inherently
distinct from the typical observables mentioned above as they are not sampled
event-by-event. Just as the moments of the cosmic microwave background power
spectrum offer insights into the universe’s structure at various length scales
in the case of cosmology, EECs map out the structure of the underlying field
theory across the “sky” of the detector. EECs are especially robust for LHC
experiments: the energy weighting naturally suppresses soft extraneous
contamination without requiring jet grooming. Unlike traditional jet
substructure observables, they can be used directly on detector tracks with
excellent angular precision without posing any limitations for theoretical
calculations. Crucially, EECs retain strong links with QFT as they can be
directly mapped to correlation functions, and, therefore, a lot of fruitful
connections between EECs and highly developed theoretical techniques in
conformal field theories and light ray operator expansion have just begun to
be explored. On the other hand, due to complicated clustering and grooming
algorithms in the typical event-by-event observables, a direct link with QFT
correlation functions is obscured.

In Ref. [5], EECs measured on boosted top quarks were proposed as a promising
mass-sensitive probe. The unique three-body (hadronic) decay of the top quark
imprints itself at a characteristic angle $\zeta_{t}\sim m_{t}^{2}/Q^{2}$ in
the three-point correlator (EEEC) spectrum. Specifically, in the equilateral
configuration, the top quark feature in the spectrum can be isolated by
imposing an asymmetry cut. This is shown in figure 2, left for ${e^{+}e^{-}}$
collisions. Zooming in, we see that it is highly sensitive to the top quark
mass (figure 2, center). Crucially, this top-resonance feature appears in the
region of the spectrum where the perturbation theory provides the dominant
description. This is because, unlike the jet mass, the soft radiation enters
the correlator only via recoil on the energetic hadrons. Consequently, EECs
measured on the top quark are not sensitive to scales below $\Gamma_{t}$. This
makes them naturally robust against hadronization and the underlying event
(UE). In contrast, the other jet substructure-based proposal of employing the
jet mass [11, 12, 13] as a $m_{t}$-sensitive probe, including with application
of jet grooming [14], is sensitive to inter-jet soft radiation at scales
parametrically lower than the top width, resulting in significantly
$\mathcal{O}(1\,\mathrm{GeV})$ large hadronization effects. The right panel of
figure 2 shows that the shift due to hadronization in the top mass is $\Delta
m_{t}^{\text{Had.}}\sim 250$ MeV. Hence, this idea overcomes the problems
shown in figure 1, achieving high top mass sensitivity in the perturbative
region for the first time. However, a direct application of this approach to
the LHC poses challenges. Unlike the jet mass and the transverse momentum
$p_{T}$, EECs are sensitive to a dimensionless angular scale. At the LHC, $Q$
must be replaced by jet $p_{T}$. Despite the theoretical elegance of this
approach, as mentioned above, the jet $p_{T}$ has large experimental
uncertainties and sensitivity to PDFs, making a direct extraction of the
$m_{t}$ from a measurement of the $\zeta_{t}$ angle challenging. In
conclusion, identifying a top mass-sensitive observable that is simultaneously
experimentally feasible at the LHC, completely robust to hadronization and the
UE, and calculable to high perturbative orders has so far remained an open
problem.

Figure 3: The shape of the three-point correlator on boosted top quark jets.
By lowering the small side cut $\zeta_{S}$, the $W$ peak emerges [15]. On the
right, the projections on the front and back faces are shown, as well as the
huge statistical gain compared to the previous equilateral projection of Ref.
[5] shown in green.

##  3 The Standard Candle Approach

In cosmology, dimensionful quantities of interest, such as distances, are not
directly obtained but indirectly through precisely measured quantities, such
as luminosity and standard candles. In our case, the dimensionless quantities
of interest are the angles that can be measured extremely precisely, thanks to
the excellent tracking system of the LHC. Recently, Ref. [15] showed that a
single EEC-based observable can be constructed that contains imprints from
angular scales of the top quark and $W$ boson decays. It is a specific
integrated EEC cross-section $T(\zeta,\zeta_{S},\zeta_{A})$ where the angle
$\zeta$ is an average of the medium and long sides of the triangle on the
celestial sphere probed by the correlator, and its shape is shown in figure 3.
Here, the medium and long sides are required not to differ more than the
asymmetry cut $\zeta_{A}$, and the short side is required to be larger than
the $\zeta_{S}$. The plot of this distribution shown in figure 3 reveals an
interesting shape: for large $\zeta_{S}$, there is a single peak at the
angular scale associated with the top decay $\zeta_{t}\sim
m_{t}^{2}/p_{T}^{2}$. Upon lowering $\zeta_{S}$, while keeping $\zeta_{A}$
fixed, we notice an emergence of a peak at the angular scale of the $W$,
$\zeta_{W}\sim m_{W}^{2}/p_{T}^{2}$. This results from the imprint of two-body
resonant decay of the $W$ boson in the “squeezed configuration” when one of
the angles is much smaller than the other two. In the right plot in figure 3,
we show the projections on the front and back faces.

Figure 4: The ratio of three-point to two-point EEC provides a robust standard
candle for identifying the angular scale associated with the $W$.

The presence of the $W$ in the distribution introduces dependence on the
$m_{t}/m_{W}$ ratio, which can be extracted from the measurements around the
two peaks in the three-point EEC. This shows that despite the complexity of
the three-body decay of the top quark, the intermediate $W$ boson is, in fact,
an important gift, enabling an elegant and complete calibration mechanism. The
standard candle observable robust to hadronization and underlying event
effects and calculable to high perturbative orders is shown in figure 4. Here,
instead of directly using the projection shown in figure 3, we consider a
ratio of the three-point to two-point EEC measured on top quark jets (removing
any restriction of the short side and asymmetry cuts), $W(\zeta)$. This is
because, unlike for the top, where $\Gamma_{t}$ provides a natural cutoff, the
$W$ boson decay is exposed to nonperturbative scales. The dominant
nonperturbative effects drop out when considering such a ratio $W(\zeta)$. In
light quark and gluon jets, these ratios have been computed to high
perturbative orders and led to the most precise jet substructure $\alpha_{s}$
determination at the LHC [16]. This is the first time such a standard candle
approach has been used in the precision jet substructure measurements. This is
entirely distinct from the direct $m_{t}$ measurements, where the $W$ decay is
reconstructed solely to achieve a fine-grained calibration of energy scales to
reduce experimental uncertainties.

##  4 Robustness against hadronization, UE, and initial state effects

Figure 5: Simulations showing remarkable properties of the ratio
$\zeta_{t}/\zeta_{W}$ in terms of independence of the $p_{T,\,\text{jet}}$
(left, top), robustness against hadronization, the UE and variations in PDF
sets (left, bottom), and promising prospects of a top mass extraction at the
HL-LHC with record precision (right).

The $W$-resonance feature captured through the standard candle observable
$W(\zeta)$ is highly correlated with the top quark, such that in the ratio
$\zeta_{t}/\zeta_{W}$, the dependence on the $p_{T}$, as well as the
systematic effects that impact $p_{T}$, entirely cancels. This is because a
single shared $p_{T}$-dependent boost determines the peak locations. This was
thoroughly verified in Ref. [15] through extensive simulation studies. The
results are shown in figure 5. The top left panel shows that the ratio is
essentially independent of the $p_{T}$. Perturbative corrections modify the
ratio’s precise value and are directly sensitive to the ratio $m_{t}/m_{W}$.

Crucially, while the jet $p_{T}$ itself is significantly impacted by
hadronization and the UE and sensitively depends on the choice of PDFs, these
effects cancel to a remarkable degree in the ratio. We see this in the left
panel, bottom. Here, the relative shifts in the peak locations are represented
as equivalent shifts in the extracted top mass $\Delta m_{t}$. Variations in
PDFs impact the $p_{T}$ spectrum, leading to shifts equivalent to a few GeVs
in each peak individually. Still, in the combination, we see a remarkable
cancellation down to almost 200 MeV (consistent with zero shift within the
error bars). Likewise, the impact of the UE, which can be as large as $3\%$ in
the individual $W$ and top peaks (also shown in figure 4 in the inset for the
$W$ peak), or equivalently $2.5\,\mathrm{GeV}$ shift in the top mass, is
completely eliminated in this standard candle approach (the second entry in
the middle panel with $\Delta m_{t}=70\pm 260$ MeV). Finally, we also see a
robust behavior against hadronization. While the hadronization shifts here can
be somewhat larger, $\sim 500$ MeV, the approach developed in the second part
of this proposal, will yield powerful methods to account for this adequately.
Nevertheless, they are still significantly smaller than $\sim 1\,\mathrm{GeV}$
nonperturbative shifts in the peak of the groomed jet mass distribution.
Finally, this observable must be precisely measurable with the data collected
during HL-LHC. The right panel of figure 5 shows an estimate of statistical
uncertainty at Runs 2 and 3 of the LHC and the HL-LHC demonstrate this is
possible. These results forecast statistical precision on $m_{t}$ at a few
hundred MeV level for the HL-LHC and with $1\,\mathrm{GeV}$ or better
precision with the LHC luminosities.

In conclusion, the promising results in the above simulation studies provide a
strong motivation for investing effort into developing this standard candle
approach. This will necessitate the development of a new formalism to describe
exclusive top decay processes. In the simulations discussed above, the ratio
$\zeta_{t}/\zeta_{W}$ was accessed by fitting the peak locations. A theory
prediction for the full shape near the top and $W$ resonance peaks, on the
other hand, will enable a much finer control over the observable, resulting in
an increased precision on the extracted $m_{t}$ and a consistent treatment of
perturbative uncertainties. While this will be challenging, the required
effort will be well worth the reward.

## Acknowledgments

I.M. is supported by start-up funds from Yale University. J.H. is supported by
the Leverhulme Trust as an Early Career Fellow. This work is supported in part
by the GLUODYNAMICS project funded by the “P2IO LabEx (ANR-10-LABX-0038)” in
the framework “Investissements d’Avenir” (ANR-11-IDEX-0003-01) managed by the
Agence Nationale de la Recherche (ANR), France.

## References

  * [1] A. H. Hoang, _What is the Top Quark Mass?_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 225 [2004.12915].
  * [2] Particle Data Group collaboration, P. Zyla et al., _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * [3] J. Piclum and C. Schwinn, _Soft-gluon and Coulomb corrections to hadronic top-quark pair production beyond NNLO_ , _JHEP_ 03 (2018) 164 [1801.05788].
  * [4] CMS collaboration, A. M. Sirunyan et al., _Measurement of $\mathrm{t\bar{t}}$ normalised multi-differential cross sections in pp collisions at $\sqrt{s}=13$ TeV, and simultaneous determination of the strong coupling strength, top quark pole mass, and parton distribution functions_, _Submitted to: Eur. Phys. J._ (2019) [1904.05237].
  * [5] J. Holguin, I. Moult, A. Pathak and M. Procura, _New paradigm for precision top physics: Weighing the top with energy correlators_ , _Phys. Rev. D_ 107 (2023) 114002 [2201.08393].
  * [6] C. Basham, L. S. Brown, S. D. Ellis and S. T. Love, _Energy Correlations in electron - Positron Annihilation: Testing QCD_ , _Phys. Rev. Lett._ 41 (1978) 1585.
  * [7] C. L. Basham, L. S. Brown, S. D. Ellis and S. T. Love, _Electron - Positron Annihilation Energy Pattern in Quantum Chromodynamics: Asymptotically Free Perturbation Theory_ , _Phys. Rev. D_ 17 (1978) 2298.
  * [8] C. L. Basham, L. S. Brown, S. D. Ellis and S. T. Love, _Energy Correlations in Perturbative Quantum Chromodynamics: A Conjecture for All Orders_ , _Phys. Lett. B_ 85 (1979) 297.
  * [9] C. Basham, L. Brown, S. Ellis and S. Love, _Energy Correlations in electron-Positron Annihilation in Quantum Chromodynamics: Asymptotically Free Perturbation Theory_ , _Phys. Rev. D_ 19 (1979) 2018.
  * [10] H. Chen, I. Moult, X. Zhang and H. X. Zhu, _Rethinking jets with energy correlators: Tracks, resummation, and analytic continuation_ , _Phys. Rev. D_ 102 (2020) 054012 [2004.11381].
  * [11] S. Fleming, A. H. Hoang, S. Mantry and I. W. Stewart, _Top Jets in the Peak Region: Factorization Analysis with NLL Resummation_ , _Phys. Rev._ D77 (2008) 114003.
  * [12] S. Fleming, A. H. Hoang, S. Mantry and I. W. Stewart, _Jets from massive unstable particles: Top-mass determination_ , _Phys. Rev._ D77 (2008) 074010 [hep-ph/0703207].
  * [13] B. Bachu, A. H. Hoang, V. Mateu, A. Pathak and I. W. Stewart, _Boosted top quarks in the peak region with N3LL resummation_ , _Phys. Rev. D_ 104 (2021) 014026 [2012.12304].
  * [14] A. H. Hoang, S. Mantry, A. Pathak and I. W. Stewart, _Extracting a Short Distance Top Mass with Light Grooming_ , _Phys. Rev._ D100 (2019) 074021 [1708.02586].
  * [15] J. Holguin, I. Moult, A. Pathak, M. Procura, R. Schöfbeck and D. Schwarz, _Using the $W$ as a Standard Candle to Reach the Top: Calibrating Energy Correlator Based Top Mass Measurements_, 2311.02157.
  * [16] CMS collaboration, _Measurement of energy correlators inside jets and determination of the strong coupling constant_ , .


# Benchmarking phonon anharmonicity in machine learning interatomic potentials

Sasaank Bandi  Department of Applied Physics and Applied Mathematics, Columbia
University, New York, NY 10027, USA  Chao Jiang  Computational Mechanics and
Materials Department, Idaho National Laboratory, Idaho Falls, ID 83415, USA.
Chris A. Marianetti  Department of Applied Physics and Applied Mathematics,
Columbia University, New York, NY 10027, USA

###### Abstract

Machine learning approaches have recently emerged as powerful tools to probe
structure-property relationships in crystals and molecules. Specifically,
Machine learning interatomic potentials (MLIP) can accurately reproduce first-
principles data at a cost similar to that of conventional interatomic
potential approaches. While MLIP have been extensively tested across various
classes of materials and molecules, a clear characterization of the anharmonic
terms encoded in the MLIP is lacking. Here, we benchmark popular MLIP using
the anharmonic vibrational Hamiltonian of ThO2 in the fluorite crystal
structure. This anharmonic Hamiltonian was constructed from density functional
theory (DFT) using our highly accurate and efficient irreducible derivative
methods, and then used to generate molecular dynamics trajectories. This data
set was used to train three classes of MLIP: Gaussian Approximation
Potentials, Artificial Neural Networks (ANN), and Graph Neural Networks (GNN).
The results were assessed by directly comparing phonons and their
interactions, as well as phonon linewidths, phonon lineshifts, and thermal
conductivity. The models were also trained on a DFT molecular dynamics
dataset, demonstrating good agreement up to fifth-order for the ANN and GNN.
Our analysis demonstrates that MLIP have great potential for accurately
characterizing anharmonicity in materials systems at a fraction of the cost of
conventional first principles-based approaches.

##  I Introduction

Phonon anharmonicity plays an essential role in the thermodynamics of
materials systems. Thus, constructing accurate vibrational Hamiltonians is
critical for materials property prediction at finite temperatures. In
particular, faithfully resolving phonon interactions up to fourth-order is
essential for capturing leading order behavior of the phonon linewidths,
phonon lineshifts, and thermal conductivity. In order to systematically assess
the phonons and their interactions, it is critical to use the minimum set of
derivatives allowed by group theory, which may be computed using irreducible
derivative approaches [1, 2, 3, 4]. Though these irreducible derivative
methods are significantly more efficient than competing finite difference
approaches, higher-order derivatives in materials with large supercells may
still be prohibitively expensive. Recently, machine learning interatomic
potentials (MLIP) have emerged as powerful tools to represent the Born-
Oppenheimer potential of materials. MLIP methods can achieve near quantum
chemical accuracy relative to the ab initio method they were trained on while
avoiding the polynomial scaling [5]. Here, we assess the fidelity of
anharmonic derivatives computed from MLIP.

MLIP aim to represent the total energy of a given atomic configuration in
terms of individual contributions from each atom. Since the complete space of
atomic displacements is high dimensional for very large molecules and
crystals, several descriptors have been developed to encode local atomic
environments [6] such as coulomb matrices [7], atom centered symmetry
functions [8], smooth overlap of atomic positions [9], and histograms of
distances, angles, and dihedrals (HDAD) [10]. While understanding how well
certain descriptor sets can represent anharmonic energy surfaces is important,
it will not be directly discussed in this work and will be the subject of
further study. Our present goal is to assess the most prevalent MLIP in the
literature, using the commonly associated descriptors. Therefore, our
benchmarks will be limited to three popular classes of potentials with
publicly available software.

The first MLIP considered in this work is the Gaussian approximation potential
(GAP) [11]. Here, each atomic energy is expanded in terms of basis functions
of the descriptors, and the distribution of basis weights is assumed to be
normally distributed about zero. Therefore, the energy of a given
configuration is also normally distributed, where the mean provides a
prediction from the model, and the variance can be a measure of uncertainity.
The advantages of this type of potential are the characterization of model
uncertainty and that the analytical form of the prediction is a sum of basis
functions. In practice, the size of basis functions may be prohibitively
large, and the “kernel trick” is employed to reduce memory and cost
requirements. It should be noted that the energy prediction provided by GAP is
equivalent to the prediction of the kernel ridge regression on the same basis
set. Recently, GAP has been used to predict second and third-order derivatives
in two dimensional materials where good agreement with DFT is achieved with
integrated properties such as phonon lifetimes and thermal conductivity [12,
13, 14, 15, 16]. Online trained GAP models have also been implemented in the
Vienna ab initio software package (VASP) [17, 18] and have also been used
along with approaches such as the stochastic self-consistant harmonic
approximation (SSCHA) [19] to predict temperature dependent phonon properties
using monte carlo sampling [20, 21].

While not directly formulated as a MLIP, several potential energy surface
fitting approaches, such as compressive sensing lattice dynamics (CSLD) [22,
23], have a similar mathematical foundation as GAP. In both models, a Bayesian
polynomial regression framework is used to fit model parameters to sampled DFT
data. While GAP generally fits descriptors of the local atomic environment,
CSLD uses the complete Cartesian coordinate basis which ensures that fitting
coefficients are more analogous to force constants. Due to the combinatorial
explosion of high-order derivatives when using the Cartesian basis, least
squares solutions tend to overfit parameters. CSLD overcomes this impediment
by fitting force constants using the least absolute shrinkage and selection
operator (LASSO), which has the same objective function as ordinary least
squares but with an additional penalty on the $L_{1}$ norm of the model
parameters [24]. The sparse solutions provided by LASSO have yielded accurate
anharmonic properties in various systems [22, 25, 23].

The second class of machine learning approaches we consider is the high
dimensional artificial neural network first described by Behler and Parrinello
(BPNN) [26]. In this model, each atomic species is represented by a
convolutional neural network, and the total energy is computed as a sum of the
output of each neural network. Historically, atom-centered symmetry functions
[8] have been used as descriptors to encode rotational invariance and have
been systematically improved over time. Recently, a set of descriptors
constructed as irreducible representations of the Euclidean group has been
shown to encode rotational equivariance, providing better predictions of
tensor properties [27], such as the forces. BPNNs have been used to predict
phonon dispersions and thermal conductivity in semiconductors [28, 29, 30,
31], new thermoelectric candidate materials [32], and superlattices [33].

Finally, the most recent developments in machine learning interatomic
potentials have focused on graph neural networks [5, 34, 35, 36, 37, 38]. In
this architecture, graphs are constructed for all neighboring atoms where
nodes encode atomic information, and bonding information is incorporated
through edge connections. Of the three models described in this work, graph
neural networks have been shown to achieve the highest accuracy in reproducing
forces and energies while requiring fewer training samples than the other two
approaches. Recently, several developments have been aimed at developing a
universal interatomic potential trained on large DFT databases, such as the
materials project [39], and phonon dispersions computed using these universal
models have shown good agreement with DFT [34, 37, 38]. While graph neural
networks have been used to predict thermal conductivity directly by training
on experimental and first principles data from materials property databases
[40, 41], we are not aware of any other work using graph neural network
interatomic potentials to compute anharmonic observables via phonon
interactions.

In this work, we benchmark how well MLIP reconstructs anharmonic potentials.
We trained the GAP, BPNN, and GNN using two datasets: one containing
configurations evaluated by an anharmonic Taylor series, which does not
contain noise, and one generated using DFT calculations. In our analysis, we
quantify errors at the level of the irreducible derivatives of the potential
energy surface, making our benchmark more robust than benchmarks using
integrated observables. In all three models, third-order derivatives are in
good agreement with the reference. Up to fifth-order derivatives are mostly
reproduced by the ANN and GNN, with the GNN demonstrating particularly
promising accuracy. We also provide comparisons of observables such as phonon
dispersions, phonon linewidths, phonon lineshifts, and thermal conductivity to
quantify the effects of the errors in the irreducible derivatives.

##  II Computational Methods and Details

DFT calculations were performed using the projector augmented wave (PAW)
method implemented in the Vienna ab initio software package (VASP) using the
local density approximation [17, 18]. The plane-wave basis cutoff energy was
set to 500 eV. A $\Gamma$ centered 10$\times$10$\times$10 $k$-point mesh was
used for primitive cell calculations, and measurements in other supercells
were conducted with corresponding $k-$mesh denisties. Gaussian smearing with
$\sigma=0.1$ was used to avoid numerical errors during $k$-point summation.
DFT energies were converged to within $10^{-6}$ eV, and the unit cell was
relaxed until all forces were within $0.02$ eV/Å.

Second (phonons), third, fourth, fifth, and sixth-order derivatives of the
Born-Oppenheimer potential of ThO2 were computed from density functional
theory and the MLIP via using irreducible derivative approaches [1]. The
finite difference calculations were extrapolated to a discretization size of
zero using quadratic errortails with ten $\Delta$s (for more information see
[1, 3]). Phonons and third-order derivatives commensurate with the $2\times
2\times 2$ supercell, which contains 24 atoms, were computed using the lone
irreducible derivative approach with forces (LID1) and energy derivatives
(LID0) for DFT and MLIP, respectively. Fourth-order phonon interactions
commensurate with the conventional cubic supercell, which contains 12 atoms,
were computed using the bundled irreducible derivative (BID) approach for DFT
due to the large computational cost. While BID generally evaluates derivatives
using the smallest set of measurements allowed by group theory, in this case
the number of measurements was tripled to ensure robust derivatives. The
fourth-order derivatives were evaluated with LID1 for all three machine
learning models. Fifth-order derivatives associated with the following $Q$
points were selectively computed using LID1 for both DFT and MLIP:
$Q=\left(\left(\frac{1}{2},0,0\right),\left(\frac{1}{2},0,0\right),\left(0,0,0\right),\left(0,0,0\right),\left(0,0,0\right)\right)$
and
$Q=\left(\left(0,0,0\right),\left(0,0,0\right),\left(0,0,0\right),\left(0,0,0\right),\left(0,0,0\right)\right)$
Finally, only sixth-order interactions commensurate with the primitive cell
were measured with LID1.

Two datasets were prepared to train the machine learning models. The first
dataset derives from a vibrational Hamiltonian containing second, third, and
fourth-order irreducible derivatives computed from DFT. Datapoints are
generated by performing molecular dynamics on the irreducible derivative
Hamiltonian using the microcanonical ensemble, which we refer to as
irreducible derivative molecular dynamics (IDMD) [4]. As a result, all
energies and forces in the dataset are in perfect agreement with the
irreducible derivatives, yielding a noiseless dataset, which we refer to as
the IDMD dataset. The IDMD was conducted in the $4\times 4\times 4$ supercell,
containing 192 atoms, for $4000$ steps with a timestep of $4$fs. MD trajectory
velocities were initialized by sampling from the Maxwell-Boltzmann
distribution at 2000K and a temperature range of 1522K-2238K was observed
throughout the simulation. The second dataset, referred to as DFTMD, was
generated using ab initio molecular dynamics in the conventional cubic cell
tripled in all three dimensions ($3S_{C}$), which contains 324 atoms. MD was
conducted using the canonical ensemble at 2000K and 3000K, and 1346 snapshots
were taken along the two trajectories. While this dataset contains noise
inherent to the numerics of the DFT calculation, it probes all derivatives
commensurate with the $3S_{C}$ supercell up to infinite-order. Each dataset
was constructed to provide a unique test of the MLIP. Since IDMD is a
noiseless dataset, this benchmark purely probes each MLIP ability to learn
anharmonic interactions while avoiding any influence of how well the model
accuracy scales with noise. On the other hand, DFTMD is a dataset that has
been constructed in accordance with the current state-of-the-art and
represents a typical dataset used for MLIP training. Since DFTMD dataset
contains infinite-order interactions, it provides the opportunity to test the
limits of each MLIP and determine the order at which the MLIP fails to
replicate DFT results.

All three MLIP were trained on IDMD, and only the BPNN and GNN were trained on
DFTMD. The GAP model was generated using the quantum mechanics and interatomic
potentials (QUIP) code [42, 43, 11]. Descriptors were constructed using two-
body and three-body terms with a cutoff of 6 Å, along with the smooth overlap
of positions (SOAP) descriptor with a cutoff of 4 Å. The neural network
potential was trained using the n2p2 code [44, 26]. Atom-centered radial and
angular symmetry functions with a cutoff radius of 10Å and 6Å were used,
respectively. For each atom, a four-layer neural network was constructed which
contained two hidden layers with 21 neurons each. The size of the input layer
for thorium atoms was 41, while the size was 46 for oxygen atoms. Finally, the
graph neural network was trained using nequIP [5]. Graph edges were
constructed with a cutoff of 6 Å with four interaction blocks. The E(3)
equivariant spherical harmonics were used with radial basis functions to form
descriptor sets [27]. More detailed model information for all three machine
learning interatomic potentials is provided in supplementary information.

##  III results

Figure 1: Results of the MLIP trained on IDMD compared to the reference, where green, blue, and red dots indicate GAP, BPNN, and GNN results, respectively. Panels $a$ and $b$ show the energies and forces for a distinct testing set, respectively. Panels $c$ and $d$ compare the third and fourth-order irreducible derivatives computed from MLIP to the reference irreducible derivatives, respectively. Figure 2: Phonon dispersion computed in the $2\hat{\bm{I}}$ FTG, where green, blue, red, and black points represent explicit measurements of the irreducible derivatives using GAP, BPNN, GNN, and the reference, respectively. Solid lines are a Fourier interpolation. Figure 3: Phonon linewidth contribution from the bubble diagram, which uses the third-order irreducible derivatives. Panels $a$ and $b$ show the linewidth of an acoustic and optical branch along a high symmetry path, respectively. The GAP, BPNN, GNN, and the reference are shown in green, blue, red, and black, respectively. Figure 4: Phonon lineshift contribution from the loop diagram, which uses the fourth-order irreducible derivatives. The BPNN, GNN, and the reference results are shown in blue, red, and black, respectively. Figure 5: Thermal conductivity evaluated using lifetimes from the bubble diagram in panel $a$ and lifetimes exclusively from the sunset diagram in panel $b$. GAP, BPNN, and GNN results are shown with green, blue, and red solid dots and the reference results are shown with black diamonds. Solid lines are drawn to direct the eyes. Figure 6: Anharmonic derivatives of the MLIP trained on DFTMD compared to the DFT reference where BPNN and GNN results are shown in blue and red, respectively. | IDMD |  | DFTMD
---|---|---|---
| 3rd Derivative | 4th Derivative |  | 3rd Derivative | 4th Derivative | 5th Derivative | 6th Derivative
GAP | 0.2207 | 1.652 |  | - | - | - | -
BPNN | 0.3100 | 0.480 |  | 0.1853 | 0.1890 | 0.4707 | 0.7965
GNN | 0.0747 | 0.212 |  | 0.0587 | 0.1346 | 0.2991 | 0.1359
Table 1: Table of root mean squared error divided by the average magnitude of
the irreducible derivatives for each model trained on the IDMD and DFTMD
datasets.

We begin by analyzing the results of the machine learning models trained on
IDMD. In Fig. 1, we show the results of all three models where GAP, BPNN, and
GNN results are shown in green, blue, and red, respectively. In figure 1
($a$), we compare energies computed on a distinct testing set generated using
IDMD, showing excellent agreement. Figure 1 ($b$) compares the forces from all
three models on the testing set; the forces from GAP are in reasonable
agreement with the reference, while the BPNN and GNN are in much better
agreement. In panels 1($c$) and 1($d$) and table 1, the third and fourth-order
irreducible derivatives are compared to the reference. In all three models,
the third-order derivatives are in good agreement with the reference, with the
GNN yielding the best fidelity of the three. At fourth-order, the BPNN and GNN
are in reasonable and good agreement with the reference, respectively, while
the GAP model fails to provide accurate irreducible derivatives.

To demonstrate the practical impact of the errors in the computed irreducible
derivatives, we compare several observables computed using irreducible
derivatives obtained from each model. We begin by presenting the phonon
dispersion (see Fig. 2). Solid dots represent grid points where irreducible
derivatives are computed, and the solid lines are a Fourier interpolation. We
demonstrate that for most phonon branches, all three models are in good
agreement with the reference, though the BPNN and the GAP both have non-
trivial discrepancies at the $X$ and $L$ points. We proceed by presenting the
phonon linewidths from the bubble diagram [4], which uses the third-order
irreducible derivatives, in a typical acoustic and optical branch (see Fig.
3). All other branches are included in supplementary material. In the acoustic
branch (Fig 3($a$)), all three models are in relatively good agreement with
the reference, though GAP overestimates the linewidths near the $X$ point by
about $16\%$. The optical branch linewidths (Fig. 3($b$)) are more sensitive
to errors in the irreducible derivatives, and thus, both GAP and BPNN yield
discrepancies up to $30\%$ between $\gamma$ and $X$. We note that the GNN
yields near perfect agreement with the reference for both branches. Next, we
present the phonon lineshifts computed by evaluating the loop diagram [4]
using the fourth-order irreducible derivatives (see Fig. 4). Only the BPNN and
GNN lineshifts have been plotted due to the massive discrepancies in the GAP
results (see supplementary material). While both models are in good agreement
with the reference for the acoustic branches, the BPNN significantly
overestimates the shifts of the optical branches, while the GNN is still
within $18\%$ of the reference. Finally, we compare the thermal conductivity
computed using the phonon lifetimes in the Boltzmann transport equation within
the relaxation time approximation. In panel 5($a$), the phonon lifetimes have
been evaluated using the bubble diagram, which only uses the third-order
derivatives, while in panel 5$b$ phonon lifetimes have been evaluated using
only the sunset diagram [4], which uses the fourth-order irreducible
derivatives. In both cases, it is evident that the thermal conductivity is
less sensitive to errors in the irreducible derivatives, and all three models
are in good agreement with the reference.

We proceed by discussing the results of the BPNN and GNN models trained on
DFTMD. Since we have already demonstrated that the GAP model fails to capture
interactions beyond third-order, it was not used for this benchmark aimed at
determining the highest order of derivatives that can be faithfully reproduced
using MLIP. We present a comparison of the third, fourth, fifth, and sixth-
order irreducible derivatives computed from the MLIP and DFT (see Fig. 6 and
Table 1). For third 6($a$), fourth 6($b$), and fifth-order 6($c$) derivatives,
both models are in excellent agreement with DFT, with the GNN outperforming
the BPNN by a factor of $1.5$ to $3$ depending on the derivative order. At
sixth-order 6($d$), both models fail to reproduce the irreducible derivatives
computed from DFT, though the graph neural network generates reasonable
results, especially on low-magnitude derivatives.

##  IV Conclusions

In summary, we have benchmarked three popular machine learning interatomic
potentials on the anharmonic irreducible derivatives of ThO2. We have
developed two training datasets: one that is computed from the anharmonic
vibrational Hamiltonian of ThO2, containing up to quartic terms, and one that
is representative of training conducted via ab initio methods in the
literature. The Behler Parrinello artificial neural network and the graph
neural network yield robust irreducible derivatives up to fifth-order, while
the Gaussian approximation potential is only able to accurately capture
anharmonic interactions up to third-order. Our work demonstrates the promising
potential of machine learning methods in characterizing anharmonicity in
materials systems. While we obtained accurate results using MLIP in the
present study, it would be interesting to verify that our findings hold in
systems with sensitivities, such as soft phonon modes. Future work will
include extending this analysis to strongly correlated electronic materials,
where the generation of accurate training data comes at a significant premium.
Additionally, to study materials with strong anharmonicity, further work will
be conducted to develop training algorithms to more accurately capture sixth-
order and beyond interactions.

##  V Acknowledgements

The development of MLIP for ThO2 was supported by the INL Laboratory Directed
Research & Development (LDRD) Program under DOE Idaho Operations Office
Contract DE-AC07-05ID14517. The construction of anharmonic Hamiltonians was
supported by the grant DE-SC0016507 funded by the U.S. Department of Energy,
Office of Science. This research used resources of the National Energy
Research Scientific Computing Center, a DOE Office of Science User Facility
supported by the Office of Science of the U.S. Department of Energy under
Contract No. DE-AC02-05CH11231.

## References

  * Fu  _et al._ [2019] L. Fu, M. Kornbluth, Z. Cheng, and C. A. Marianetti, Group theoretical approach to computing phonons and their interactions, Physical Review B 100, 014303 (2019).
  * Mathis  _et al._ [2022] M. A. Mathis, A. Khanolkar, L. Fu, M. S. Bryan, C. A. Dennett, K. Rickert, J. M. Mann, B. Winn, D. L. Abernathy, M. E. Manley, D. H. Hurley, and C. A. Marianetti, Generalized quasiharmonic approximation via space group irreducible derivatives, Physical Review B 106, 014314 (2022).
  * Bandi and Marianetti [2023] S. Bandi and C. A. Marianetti, Precisely computing phonons via irreducible derivatives, Physical Review B 107, 174302 (2023).
  * Xiao and Marianetti [2023] E. Xiao and C. A. Marianetti, Anharmonic phonon behavior via irreducible derivatives: Self-consistent perturbation theory and molecular dynamics, Physical Review B 107, 094303 (2023).
  * Batzner  _et al._ [2022] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nature Communications 13, 10.1038/s41467-022-29939-5 (2022).
  * Langer  _et al._ [2022] M. F. Langer, A. Goeßmann, and M. Rupp, Representations of molecules and materials for interpolation of quantum-mechanical simulations via machine learning, npj Computational Materials 8, 10.1038/s41524-022-00721-x (2022).
  * Rupp  _et al._ [2012] M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, Fast and accurate modeling of molecular atomization energies with machine learning, Physical Review Letters 108, 058301 (2012).
  * Behler [2011] J. Behler, Atom-centered symmetry functions for constructing high-dimensional neural network potentials, The Journal of Chemical Physics 134, 10.1063/1.3553717 (2011).
  * Bartók  _et al._ [2013] A. P. Bartók, R. Kondor, and G. Csányi, On representing chemical environments, Physical Review B 87, 184115 (2013).
  * Faber  _et al._ [2017] F. A. Faber, L. Hutchison, B. Huang, J. Gilmer, S. S. Schoenholz, G. E. Dahl, O. Vinyals, S. Kearnes, P. F. Riley, and O. A. von Lilienfeld, Prediction errors of molecular machine learning models lower than hybrid dft error, Journal of Chemical Theory and Computation 13, 5255 (2017).
  * Bartók  _et al._ [2010] A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons, Physical Review Letters 104, 136403 (2010).
  * Kocabaş  _et al._ [2023] T. Kocabaş, M. Keçeli, . Vázquez-Mayagoitia, and C. Sevik, Gaussian approximation potentials for accurate thermal properties of two-dimensional materials, Nanoscale 15, 8772 (2023).
  * Ouyang  _et al._ [2020] Y. Ouyang, Z. Zhang, C. Yu, J. He, G. Yan, and J. Chen, Accuracy of machine learning potential for predictions of multiple-target physical properties*, Chinese Physics Letters 37, 126301 (2020).
  * Zhang  _et al._ [2020] Y. Zhang, C. Shen, T. Long, and H. Zhang, Thermal conductivity of h-bn monolayers using machine learning interatomic potential, Journal of Physics: Condensed Matter 33, 105903 (2020).
  * Wei  _et al._ [2022] Z. Wei, C. Zhang, Y. Kan, Y. Zhang, and Y. Chen, Developing machine learning potential for classical molecular dynamics simulation with superior phonon properties, Computational Materials Science 202, 111012 (2022).
  * Arabha  _et al._ [2021] S. Arabha, Z. S. Aghbolagh, K. Ghorbani, S. M. Hatam-Lee, and A. Rajabpour, Recent advances in lattice thermal conductivity calculation using machine-learning interatomic potentials, Journal of Applied Physics 130, 10.1063/5.0069443 (2021).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initiomolecular dynamics for liquid metals, Physical Review B 47, 558 (1993).
  * Monacelli  _et al._ [2021] L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea, and F. Mauri, The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects, Journal of Physics: Condensed Matter 33, 363001 (2021).
  * Ranalli  _et al._ [2023] L. Ranalli, C. Verdi, L. Monacelli, G. Kresse, M. Calandra, and C. Franchini, Temperature‐dependent anharmonic phonons in quantum paraelectric ktao3 by first principles and machine‐learned force fields, Advanced Quantum Technologies 6, 10.1002/qute.202200131 (2023).
  * Verdi  _et al._ [2023] C. Verdi, L. Ranalli, C. Franchini, and G. Kresse, Quantum paraelectricity and structural phase transitions in strontium titanate beyond density functional theory, Physical Review Materials 7, l030801 (2023).
  * Zhou  _et al._ [2014] F. Zhou, W. Nielson, Y. Xia, and V. Ozoliņš, Lattice anharmonicity and thermal conductivity from compressive sensing of first-principles calculations, Physical Review Letters 113, 185501 (2014).
  * Zhou  _et al._ [2019] F. Zhou, W. Nielson, Y. Xia, and V. Ozoliņš, Compressive sensing lattice dynamics. i. general formalism, Physical Review B 100, 184308 (2019).
  * Tibshirani [1996] R. Tibshirani, Regression shrinkage and selection via the lasso: A retrospective, Journal of the Royal Statistical Society Series B: Statistical Methodology 73, 273 (1996).
  * Tadano and Tsuneyuki [2015] T. Tadano and S. Tsuneyuki, Self-consistent phonon calculations of lattice dynamical properties in cubic srtio3 with first-principles anharmonic force constants, Physical Review B 92, 054301 (2015).
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, Generalized neural-network representation of high-dimensional potential-energy surfaces, Physical Review Letters 98, 146401 (2007).
  * Geiger and Smidt [2022] M. Geiger and T. Smidt, e3nn: Euclidean neural networks (2022).
  * Li  _et al._ [2020a] R. Li, E. Lee, and T. Luo, A unified deep neural network potential capable of predicting thermal conductivity of silicon in different phases, Materials Today Physics 12, 100181 (2020a).
  * Li  _et al._ [2020b] R. Li, Z. Liu, A. Rohskopf, K. Gordiz, A. Henry, E. Lee, and T. Luo, A deep neural network interatomic potential for studying thermal conductivity of β -ga2o3, Applied Physics Letters 117, 10.1063/5.0025051 (2020b).
  * Mangold  _et al._ [2020] C. Mangold, S. Chen, G. Barbalinardo, J. Behler, P. Pochet, K. Termentzidis, Y. Han, L. Chaput, D. Lacroix, and D. Donadio, Transferability of neural network potentials for varying stoichiometry: Phonons and thermal conductivity of mnxgey compounds, Journal of Applied Physics 127, 10.1063/5.0009550 (2020).
  * Minamitani  _et al._ [2019] E. Minamitani, M. Ogura, and S. Watanabe, Simulating lattice thermal conductivity in semiconducting materials using high-dimensional neural network potential, Applied Physics Express 12, 095001 (2019).
  * Choi  _et al._ [2022] J. M. Choi, K. Lee, S. Kim, M. Moon, W. Jeong, and S. Han, Accelerated computation of lattice thermal conductivity using neural network interatomic potentials, Computational Materials Science 211, 111472 (2022).
  * Qu  _et al._ [2023] R. Qu, Y. Lv, and Z. Lu, A deep neural network potential to study the thermal conductivity of mnbi2te4 and bi2te3/mnbi2te4 superlattice, Journal of Electronic Materials 52, 4475 (2023).
  * Takamoto  _et al._ [2022] S. Takamoto, S. Izumi, and J. Li, Teanet: Universal neural network interatomic potential inspired by iterative electronic relaxations, Computational Materials Science 207, 111280 (2022).
  * Batatia  _et al._ [2022] I. Batatia, D. P. Kovács, G. N. C. Simm, C. Ortner, and G. Csányi, Mace: Higher order equivariant message passing neural networks for fast and accurate force fields (2022).
  * Gasteiger  _et al._ [2020] J. Gasteiger, J. Groß, and S. Günnemann, Directional message passing for molecular graphs (2020).
  * Chen and Ong [2022] C. Chen and S. P. Ong, A universal graph deep learning interatomic potential for the periodic table, Nature Computational Science 2, 718 (2022).
  * Deng  _et al._ [2023] B. Deng, P. Zhong, K. Jun, J. Riebesell, K. Han, C. J. Bartel, and G. Ceder, Chgnet as a pretrained universal neural network potential for charge-informed atomistic modelling, Nature Machine Intelligence 5, 1031 (2023).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, Commentary: The materials project: A materials genome approach to accelerating materials innovation, APL Materials 1, 10.1063/1.4812323 (2013).
  * Ojih  _et al._ [2023] J. Ojih, A. Rodriguez, J. Hu, and M. Hu, Screening outstanding mechanical properties and low lattice thermal conductivity using global attention graph neural network, Energy and AI 14, 100286 (2023).
  * Sun and Hu [2021] Y. Sun and W. Hu, Novel machine learning framework for thermal conductivity prediction by crystal graph convolution embedded ensemble, SmartMat 3, 474 (2021).
  * Csányi  _et al._ [2007] G. Csányi, S. Winfield, J. R. Kermode, A. De Vita, A. Comisso, N. Bernstein, and M. C. Payne, Expressive programming for computational physics in fortran 95+, IoP Comput. Phys. Newsletter (2007).
  * Kermode [2020] J. R. Kermode, f90wrap: an automated tool for constructing deep python interfaces to modern fortran codes, Journal of Physics: Condensed Matter 32, 305901 (2020).
  * Singraber  _et al._ [2019] A. Singraber, J. Behler, and C. Dellago, Library-based lammps implementation of high-dimensional neural network potentials, Journal of Chemical Theory and Computation 15, 1827 (2019).