* Li [1997] T. Li.  Numerical solution of multivariate polynomial systems by homotopy continuation methods.  _Acta numerica_ , 6:399–436, 1997.  doi: 10.1017/s0962492900002749.
  * Wigner [1931 pp 251-254] E. P. Wigner.  _Gruppentheorie_.  Vieweg, 1931 pp 251-254.  doi: doi.org/10.1007/978-3-663-02555-9.
  * Wigner [1959 pp 233-236] E. P. Wigner.  _Group Theory_.  Academic Press Inc., 1959 pp 233-236.  doi: 10.1016/b978-0-12-750550-3.x5001-0.
  * Schumacher and Westmoreland [2010] B. Schumacher and M. Westmoreland.  _Quantum processes systems, and information_.  Cambridge University Press, 2010.  doi: 10.1017/CBO9780511814006.
  * Blum et al. [1996] L. Blum, F. Cucker, M. Shub, and S. Smale.  Complexity and real computation: A manifesto.  _Int. J. Bifurc. Chaos_ , 6(01):3–26, 1996.  doi: 10.1142/s0218127496001818.
  * Faugere [1999] J. Faugere.  A new efficient algorithm for computing Gröbner bases (F4).  _J. Pure Appl. Algebr._ , 139(1):61–88, 1999\.  doi: 10.1016/s0022-4049(99)00005-5.
  * Buchberger [1976] B. Buchberger.  A theoretical basis for the reduction of polynomials to canonical forms.  _ACM SIGSAM Bulletin_ , 10(3):19–29, 1976.  doi: 10.1145/1088216.1088219.
  * [33] Private communication with Dr. Allan Steel.
  * Cox et al. [1992] D. Cox, J. Little, and D. O’Shea.  _Ideals, varieties, and algorithms_ , volume 3.  Springer, 1992.  doi: 10.1007/978-1-4757-2181-2.
  * [35] Note that this ‘tracking’, a numerical technique of polynomial homotopy continuation, is unrelated to the efficient ‘tracking’ of open quantum systems, which is the physics problem under investigation.
  * Hill et al. [2018] C. Hill, K. Lee, A. Leykin, T. Duff, A. Jensen, and J. Sommars.  Solving polynomial systems via homotopy continuation and monodromy.  _IMA J. Numer. Anal._ , 04 2018.  doi: 10.1093/imanum/dry017.
  * Gorini et al. [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan.  Completely positive dynamical semigroups of N-level systems.  _J. Math. Phys._ , 17(5):821–825, 1976.  doi: 10.1063/1.522979.
  * [38] To prove that $L=1$ PREs are not possible for $K=4,5$ for the example MEs chosen, we rely upon ruling out various graphs using the considerations of Sec. 3.1. For the graphs that are not ruled out in this way, we then obtain computational Nullstellensätze showing that no PRE can exist. In contrast to the case where a Nullstellensatz is obtained for a fully connected graph (for the $K=3$), for these $L=1$, $K=4,5$ examples there are many isomorphically different graphs that can be drawn when the number of non-zero rates is less than the maximal number. That is, ruling out the graph with the optimal (largest) number of transitions is not sufficient to prove that no PRE can exist. Instead, we must consider all the isomorphically different graphs with fewer rates also and then either rule them out via Sec. 3.1 or by Nullstellensatz. This difficulty does not arise if a Nullstellensatz can be performed on the fully connected graph as all graphs with fewer transitions can be obtained from this graph just by setting some of the transition rates to zero. Consequently, one must first identify all the relevant graphs (there are about 50) and then obtain Nullstellensätze for them in order to have an overall proof of $K\leq 5$ PRE non-existence for a single example ME. Some graph Nullstellensätze take longer than others, with a total runtime for each ME of around 2 days.
  * Li et al. [1989] T. Y. Li, T. Sauer, and J. A. Yorke.  The cheater’s homotopy: an efficient procedure for solving systems of polynomial equations.  _SIAM J. Numer. Anal._ , 26(5):1241–1251, 1989\.  doi: 10.1137/0726069.
  * [40] There is no guarantee that a numerical search for solutions to a polynomial system will actually find all existing solutions. Algorithms have various defined tolerances that may not be tight enough to faithfully explore the search space for a given ME. This is in contrast to the case where a Hilbert Nullstellensatz is found, which is a guarantee that no solutions can exist. However, in the case of square polynomial systems, it has always been the case that there will exist non-PRE solutions to the constraints (that is, complex solutions, whereas PREs require real valued solutions as the real and imaginary equations have been separated out). These complex solutions will mean that a Nullstellensatz will not exist, even in the absence of PREs, and cannot be utilized as an exact proof.
  * Hardy [2001] L. Hardy.  Quantum theory from five reasonable axioms.  _arXiv:quant-ph/0101012_ , 2001.
  * Monras et al. [2011] A. Monras, A. Beige, and K. Wiesner.  Hidden quantum Markov models and non-adaptive read-out of many-body states.  _Appl. Math. and Comp. Sciences_ , 3:93, 2011.
  * Gu et al. [2012] M. Gu, K. Wiesner, E. Rieper, and V. Vedral.  Quantum mechanics can reduce the complexity of classical models.  _Nat. Commun._ , 3, Mar 2012.  doi: 10.1038/ncomms1761.
  * Mehta [2011a] D. Mehta.  Numerical polynomial homotopy continuation method and string vacua.  _Adv. High Energy Phys._ , 2011, 2011a.  doi: 10.1155/2011/263937.
  * Mehta [2011b] D. Mehta.  Finding all the stationary points of a potential-energy landscape via numerical polynomial-homotopy-continuation method.  _Phys. Rev. E_ , 84(2):025702, 2011b.  doi: 10.1103/physreve.84.025702.
  * Mehta et al. [2016] D. Mehta, H. D. Nguyen, and K. Turitsyn.  Numerical polynomial homotopy continuation method to locate all the power flow solutions.  _IET Gener. Transm. Distrib._ , 10(12):2972–2980, 2016.  doi: 10.1049/iet-gtd.2015.1546.
  * Sommese and Wampler II [2005] A. J. Sommese and C. W. Wampler II.  _The Numerical solution of systems of polynomials arising in engineering and science_.  World Scientific, 2005.  doi: 10.1142/5763.
  * Allgower and Georg [2012] E. L. Allgower and K. Georg.  _Numerical continuation methods: an introduction_ , volume 13.  Springer Science & Business Media, 2012.  doi: 10.1137/1.9780898719154.
  * Leykin et al. [2006] A. Leykin, J. Verschelde, and Y. Zhuang.  Parallel homotopy algorithms to solve polynomial systems.  In _ICMS_ , pages 225–234. Springer, 2006.  doi: 10.1007/11832225˙22.
  * Malajovich and Meer [2007] G. Malajovich and K. Meer.  Computing minimal multi-homogeneous Bézout numbers is hard.  _Theory Comput. Syst._ , 40(4):553–570, 2007\.  doi: 10.1007/s00224-006-1322-y.
  * [51] T. Duff, C. Hill, A. Jensen, K. Lee, A. Leykin, and J. Sommars. MonodromySolver: a Macaulay2 package for solving polynomial systems via homotopy continuation and monodromy. Available at http://people.math.gatech.edu/∼aleykin3/MonodromySolver.


# Efficient Modelling of Anharmonicity and Quantum Effects in PdCuH2 with
Machine Learning Potentials

Francesco Belli  Department of Chemistry, State University of New York at
Buffalo, Buffalo, NY 14260-3000, USA  Eva Zurek<EMAIL_ADDRESS>Department
of Chemistry, State University of New York at Buffalo, Buffalo, NY 14260-3000,
USA

###### Abstract

Quantum nuclear effects and anharmonicity impact a wide range of functional
materials and their properties. One of the most powerful techniques to model
these effects is the Stochastic Self-Consistent Harmonic Approximation
(SSCHA). Unfortunately, the SSCHA is extremely computationally expensive,
prohibiting its routine use. We propose a protocol that pairs machine learning
interatomic potentials, which can be tailored for the system at hand via
active learning, with the SSCHA. Our method leverages an upscaling procedure
that allows for the treatment of supercells of up to thousands of atoms with
practically minimal computational effort. The protocol is applied to PdCuHx
($x=0-2$) compounds, chosen because previous experimental studies have
reported superconducting critical temperatures, $T_{\text{c}}$s, as high as 17
K at ambient pressures in an unknown hydrogenated PdCu phase. We identify a
$P4/mmm$ PdCuH2 structure, which is shown to be dynamically stable only upon
the inclusion of quantum fluctuations, as being a key contributor to the
measured superconductivity. For this system, our methodology is able to reduce
the computational expense for the SSCHA calculations by $\sim$96%. The
proposed protocol opens the door towards the routine inclusion of quantum
nuclear motion and anharmonicity in materials discovery.

##  I Introduction

In the framework of standard Density Functional Theory (DFT), two fundamental
approximations are commonly made. The first is concerned with the classical
treatment of the nuclei, assumed to be fixed points in space generating an
effective potential within which the electrons move. The second is the
harmonic treatment of this potential, which is approximated through a second-
order expansion in the vicinity of its local minima. These approximations are
key for enabling the calculation of thermodynamic properties such as the
harmonic free energy, the harmonic vibrational spectra, and the heat capacity
of chemical compounds and materials [1].

In some circumstances however, these fundamental approximations are no longer
valid. Two examples are systems at high temperatures whose atoms undergo large
amplitude atomic vibrations, or those containing light nuclei whose quantum
oscillations cannot be disregarded. In both cases the atoms explore features
of the potential extending far beyond the limits of a second-order expansion
around the local minima [2], defining the presence of anharmonicity, driven by
temperature, or quantum nuclear effects. Anharmonicity is key for a broad
range of materials properties including the emergence of charge density waves
(CDW) in materials [3, 4], the enhancement or suppression of ferroelectricity
in perovskites [5], and thermodynamic properties of materials at planetary
temperatures [6]. However, the effects of anharmonicity are most prominent in
hydrogen-rich systems, where they can alter structural phase stability,
renormalize phonon spectra and dramatically modify superconducting properties
[7, 8, 9, 10, 11].

The importance of anharmonic effects has led to the development of a number of
computational techniques that can be used to treat them. One class of
approaches is based upon molecular dynamics simulations carried out above the
Debye temperature (assuming classical nuclei) [12, 13, 14, 15, 16, 17, 18,
19], or those of the path integral variant [20] (both of these methods are
able to produce accurate results in the low temperature quantum limit).
Another approach to describe anharmonicity to a non-perturbative level is
based upon the self-consistent harmonic approximation (SCHA) [21], whose
latest incarnation is termed the Stochastic SCHA (SSCHA) [22, 23, 24, 25, 26].
The SSCHA is able to comprehensively treat anharmonicity arising from finite
temperature and from quantum nuclear effects, and it has been applied to a
variety of systems including the high pressure hydrides [7, 8, 9, 10, 11].
Unfortunately, the predictive capacity of all of these methods has been
limited by the computational cost associated with their requisite molecular
dynamics or Monte-Carlo simulations, thereby imposing a limitation on the size
of systems that can be modelled and the accuracy that can be attained.

In addressing this challenge, Machine Learning Potentials (MLPs) emerge as a
viable solution [27, 28, 29, 30, 31, 32]. MLPs rely on adaptable mathematical
expressions used to derive analytical potential energy surfaces. The
parameters of these expressions are trained on electronic structure reference
data that usually includes geometric coordinates, as well as their associated
energies, forces and stresses (EFS). As a result, MLPs are able to preserve
the high accuracy of electronic structure methods, but at a much lower
computational cost. Some of the most popular MLP variants include neural
network potentials [33, 34], Gaussian approximation potentials (GAPs) [35],
and moment tensor potentials (MTPs) [36, 37]. These potentials have proven to
be powerful tools in atomistic modelling including accelerating materials
discovery [38, 39], reducing computational load, upscaling simulations to
bigger sizes with minimal computational effort [40], and in simulating
superstructural ordering and mixed chemical occupancy [41]. Additionally, MLPs
have already been used to capture anharmonicity and quantum fluctuations,
setting the stage for the work described herein [5, 42, 43].

In what follows, an efficient protocol to combine MTPs with the SSCHA to model
the effect of anharmonicity and quantum nuclear effects is proposed. Although
MLPs have already been applied in conjunction with the SSCHA [5, 42] there are
still open questions related to the proper efficient training of the potential
and its accuracy when extrapolating to larger supercells. Here these questions
are addressed by creating an active learning procedure combining the SSCHA
with MTPs through the MLIP (Machine Learning Interatomic Potential) package
[37]. The outlined methodology allows for a drastic reduction of the
computational cost while providing a solid framework to assess the accuracy
and scalability of the developed MLPs.

Further, the computational framework is applied to PdCuH2, a cousin of PdH – a
highly anharmonic conventional superconductor famous for its inverse isotope
effect wherein the superconducting critical temperature, $T_{\text{c}}$, of
the deuterated form is higher than that of its lighter brethren [44]. In 1974
a PdCuHx compound with $x\sim 1.4$ was synthesized and its $T_{\text{c}}$ was
reported to be 17 K [45], nearly double the value measured for PdH. A recent
computational study concluded that the most likely candidate for the
superconducting phase was a $P4/mmm$ PdCuH2 stoichiometry structure, but an
accurate calculation of its $T_{\text{c}}$ could not be performed due to the
presence of imaginary phonon modes, which were hypothesized to be a
manifestation of the strong anharmonicity [46]. In what follows, we show that
quantum nuclear effects render this PdCuH2 phase dynamically stable, yielding
a $T_{\text{c}}$ that is in-line with experiment. Additionally, we searched
for a possible inverse isotope effect by calculating the $T_{\text{c}}$ of
PdCuD2, but none was found.

##  II Results

### SSCHA + MLIP Workflow

The first step of the SSCHA involves the generation of a set of displaced
atomic configurations, sampled from a Gaussian distribution with a width
linked to the magnitude of the computed force constants [25]. Consequently, an
initial guess for the dynamical matrices, typically obtained within the
harmonic approximation for the chosen supercell, is required before
initialising the SSCHA. For accurate results, the chosen supercell should be
large enough to capture the necessary atomic interactions. In the procedure
implemented here the initial guess is made for a supercell smaller then the
one needed for convergence – foreshadowing the upscaling procedure that
distinguishes this workflow from those proposed by others.

Figure 1:  SSCHA + MLIP workflow. Flowchart for the workflow that combines the
Stochastic Self-Consistent Harmonic Approximation (SSCHA) with the Machine
Learning Interatomic Potential (MLIP) package for the generation of moment
tensor potentials (MTPs). In this flowchart, parallelograms represent input
and output files, rectangles denote individual workflow steps, and rhombi
indicate logic operations. Large rounded rectangles encompass sections of the
program related to the SSCHA and MLIP packages. This workflow fully integrates
both programs and incorporates an on-the-fly procedure for MTP training,
significantly reducing the number of required DFT calculations and automating
dataset training.

It is at this step that the SSCHA+MLIP integration, based upon an active
learning protocol, begins: the workflow is depicted schematically in Figure 1.
Displaced configurations for the supercell are created, and single point DFT
calculations are performed to obtain the EFS of all of the members of this
first population, to be used for training the MTP. In subsequent populations,
the extrapolation grade, $\gamma$, of each structure is determined via the
generalised D-optimality criterion, which is able to detect the most diverse
configurations [37]. To ensure accurate predictions we use a very conservative
value of $\gamma_{\text{select}}=2$, chosen so that every geometry with
$\gamma_{\text{select}}>2$, which is not well described by the MTP, has their
structures and EFS added to the training set for retraining the potential.

The retrained MTP is then employed to predict the EFS of the members of this
population not chosen for training. After the EFS of all of the configurations
in the population have been obtained (via MTP or DFT) they are combined and
fed to the SSCHA to perform the optimisation of the force constants and
centroid positions. This procedure is repeated for subsequent populations
until the SSCHA is converged. Though the MTP parameters are updated for each
population, they are not used to predict the EFS of members of prior
populations because of the cautious value of $\gamma_{\text{select}}$ and
because in the SSCHA each geometry is independent of the others, unlike in a
molecular dynamics run or in a relaxation trajectory where the same set of MTP
parameters must be employed to avoid the introduction of discontinuities in
the potential energy surface.

Figure 2:  Effective process time. The CPU time, in hours, required for
executing the full active learning and upscaling procedure starting with the
(a) 1$\times$1$\times$1 cell and extending to the (b) 2$\times$2$\times$2 and
the (c) 3$\times$3$\times$3 supercells. The orange lines represent the total
CPU time that would be necessary without the machine learning protocol; the
dark green lines represent the actual CPU time utilised. In (a) these lines
overlap. The evolution of the difference between the anharmonic free energy,
$\Delta F$, relative to its minimum calculated value as a function of the
population during training for (d) 1$\times$1$\times$, (e) 2$\times$2$\times$2
and (f) 3$\times$3$\times$3 supercells. The potential was trained at a
temperature of 300 K and can be used to reproduce supercells larger than
3$\times$3$\times$3.

The cycle consisting of generating displaced atomic configurations, passing
them through the active learning procedure for EFS predictions, and the SSCHA
optimization is repeated until the SSCHA reaches convergence. Subsequently, an
upscaling procedure is performed, in which the SSCHA is reinitialized with an
increased supercell size. This step is crucial for properly modelling
phenomena such as CDWs, which might not be captured in smaller cells. A tight-
binding model is employed to extrapolate the starting guess for the new
upscaled dynamical matrices, from the converged dynamical matrix obtained in
the previous smaller cell.

The SSCHA+MLIP cycle is repeated again, except now the first iteration of the
SSCHA employs the final MTP developed from the previous run on the smaller
cell, instead of DFT. As before, when extrapolative configurations are
encountered they are stored and may be used to retrain the MTP. This upscaling
process is then iterated until reaching the desired size for the supercell.
The rationale underlying the upscaling procedure is rooted in the idea that
the majority of atomic interactions required for training the potential are
already present in the small cells [47] and that progressively fewer
configurations are required to retrain the potential when considering bigger
cells.

The results of the training procedure for the specific case of $P4/mmm$
PdCuH2, whose phonon spectra and superconducting properties are investigated
in the following section, is reported in Figure 2. The SSCHA procedure began
from calculations performed on the four atom unit cell. In the first three
populations all of the 140 configurations were evaluated by DFT. Consequently,
the yellow points (CPU time required for a traditional SSCHA run) overlap with
the green points (CPU time required for the MLIP accelerated SSCHA) in the top
panel. After these three populations, the SSCHA was converged for the unit
cell and the calculations were upscaled to a 2$\times$2$\times$2 supercell
containing 32 atoms. In this supercell, the SSCHA was converged after the
6${}^{\text{th}}$ population and only 342 out of the 1500 configurations
(22.8%) were deemed extrapolative requiring further DFT calculations, whose
results were employed to retrain the potential. The system was further
upscaled to a 3$\times$3$\times$3 supercell containing 108 atoms; from the
2180 configurations required for SSCHA convergence only 56 (2.57%) were sent
for DFT calculations. After the training on the 3$\times$3$\times$3 supercell
was complete, the MTP did not require any additional DFT data while upscaling
to even larger cells. Though the potential was explicitly trained on a cells
containing up to 108 atoms, it could be used to perform the SSCHA on a
6$\times$6$\times$5 supercell containing 720 atoms. In this specific case
12,000 configurations and four populations were necessary to achieve
convergence for the 6$\times$6$\times$5 supercell. Because of the formal
$N^{3}$ scaling in DFT and the increasingly large number of configurations
required by the SSCHA for larger unit cells, such calculations would have been
impossible without the MLP acceleration. When taking into account the amount
of time required to converge the 3$\times$3$\times$3 supercell, the SSCHA +
MLIP workflow results in a reduction in computational resources by an
impressive $\sim$96%. A full DFT SSCHA run was not performed for the
6$\times$6$\times$5 supercell, since this would have been impossible during a
typical postdoc lifetime, and particular effort was invested in safeguarding
the sanity and mental well-being of the one involved in this investigation.

The accuracy and the predictive capabilities of the potential were analysed by
inspecting the root mean square (RMS) prediction errors for the EFS as a
function of the size of the supercell (Supplementary Figure 2). The prediction
errors were obtained by comparing the DFT EFS of an auxilliary SSCHA-generated
set of structures for the given supercell size with those predicted by the MTP
trained on a particular supercell with the SSCHA+MLIP procedure, as described
more in depth in Supplementary note 1. Although the MTP was relatively
accurate after training on a 2$\times$2$\times$2 supercell, the accuracy we
desire was only achieved after training on the 3$\times$3$\times$3\. The
resulting MTP gives an RMS error of less than: 0.5 meV/atom for the energies,
0.05 meV/(atom Å ) for the forces, and 0.5 GPa for the stresses. The necessity
for training on the 3$\times$3$\times$3 supercell can also be noticed by
inspecting the behavior of the free energy at 300 K as a function of the
supercell and population, reported in Figure 2. Although, one might suspect
convergence is reached at the end of the 2$\times$2$\times$2 run, a small
decrease in the free energy is observed when the training is completed on the
3$\times$3$\times$3 cell.

The SSCHA + MLIP workflow, including the active learning protocol and the
upscaling procedure, is more efficient and complete than recently proposed
workflows [5, 42]. A key advantage of our approach is that it enables the
automation of the training and facilitates the evaluation of its efficacy,
ensuring accuracy in the predictions, via the active learning phase.
Additionally, this methodology not only provides the possibility of
drastically decreasing the computational cost required for calculating the EFS
during the upscaling procedure, but it also enables the identification of
emerging phenomena such as CDWs, which potentially require retraining of the
potential. Furthermore, a MLP developed at a particular temperature and
pressure can be employed as a starting guess for a simulation performed for a
different set of thermodynamic variables. For example, an MTP trained at high
temperatures could be used for a calculation performed at low temperatures,
requiring only a few additional DFT calculations for retraining via the active
learning protocol. Similarly, an MTP trained at a few pressures in a narrow
pressure range could be employed to perform the SSCHA dynamics allowing for a
full cell relaxation with minimal computational expense.

Figure 3: Crystal structures. The most stable PdCuH2 phases proposed in (a)
Ref. [46], and (b) that found here. Both studies identified the same PdCuH
phase as the ground state. It can be obtained from $P4/mmm$ PdCuH2 by removing
the hydrogen atoms that are at the same sites as the hydrogen in the center of
the blue transparent octahedron.

### PdCuHx ($\mathbf{x=0-2}$) Structures through Crystal Structure Prediction

As the hydrogen content and positions within the synthesized PdCuHx compound
are unknown [45], crystal structure prediction (CSP) calculations with the
XtalOpt evolutionary algorithm were performed for $x=0-2$ at zero temperature
and pressure with 1-4 formula units in the cell. The equimolar PdCu alloy
adopts a disordered fcc structure and an ordered CsCl (B2) phase, and at high
temperatures experiments find a mixture of the two [48, 49]. Our CSP searches
identified the B2 ($Pm\bar{3}m$ symmetry) ordered phase as being the ground
state at 0 K, verifying the CSP framework and choice of computational
parameters.

The small hydrogen atoms are often adsorbed into the interstitial sites of
transition metal solids. Therefore, it was not a surprise that the most stable
PdCuH phase found in our CSP runs could be obtained from $Pm\bar{3}m$ PdCu by
filling half of its octahedral sites, specifically those that are coordinated
equatorially by four palladium atoms and axially by two copper atoms (cf.
Figure 3 (a)) with hydrogen. This structure was also predicted to be the
ground state by Vocaturo et al. who found it via random placement of hydrogen
within the parent PdCu unit cell followed by local relaxation [46]. Filling
the second octahedral site, this time with copper on the equatorial and
palladium on the axial positions leads to the most stable PdCuH2 structure
identified by Vocaturo and co-workers with $P4/mmm$ symmetry [46], as shown in
Figure 3 (a). Another PdCuH2 phase they identified, which was $\sim$50
meV/atom higher in energy than $P4/mmm$ PdCuH2, was characterized by hydrogen
atoms on both of the tetrahedral holes, highlighting the preference of
hydrogen occupying the octahedral rather than the tetrahedral sites in
hydrogenated PdCu alloys.

The evolutionary searches performed here, on the other hand, predicted an
$I4_{1}/amd$ symmetry PdCuH2 structure (Figure 3(b)) as being the ground
state. It was calculated to be 18 meV/atom _more stable_ than the previously
proposed $P4/mmm$ PdCuH2. This newfound $I4_{1}/amd$ symmetry phase is also
characterized by hydrogen-filling of the octahedral sites, however it differs
from $P4/mmm$ in the arrangement of the metal atoms. Whereas the latter
contains alternating layers of Pd and Cu stacked in an ABA… fashion, in
$I4_{1}/amd$ each layer is constituted by an equal proportion of Cu and Pd
atoms with an ABCDA… stacking along the $c$-axis. In $I4_{1}/amd$ each
hydrogen atom occupies the center of an octahedron where two copper atoms and
two palladium atoms, diagonal to each other, constitute the equatorial sites.
Similarly, one of each type of metal atom lies on the axial positions. As a
result, the Pd-H and Cu-H distances are the same in the equatorial plane, but
they differ along the $c$ lattice vector with slightly shorter Pd-H (1.97 Å)
and slightly longer Cu-H (2.25 Å) distances so that the hydrogen atoms do not
all lie in the same $ab$ plane, unlike in $P4/mmm$ PdCuH2.

### Electronic Structure and Superconductivity

Now that we have described the structural peculiarities of the most stable
PdCuH phase identified via both random [46] and evolutionary searches, as well
as two PdCuH2 phases – a metastable one that can be obtained by filling the
octahedral holes of B2 PdCu with hydrogen, and the ground state, which
possesses a different arrangement of the transition metal atoms, let us
investigate their electronic structure and propensity for superconductivity.
Specifically, we aim to determine which of these compounds, if any, could be a
potential contributor to the measured $T_{\text{c}}$ of 17 K [45] in a
hydrogenated PdCu sample. Since previous computations classified PdCuH as a
sub-Kelvin superconductor [46], we turned our focus to the PdCuH2
stoichiometry.

Figure 4:  Electronic band structures. Electronic band structure and atom
projected densities of states (DOS) of (a) $P4/mmm$ PdCuH2 and (b)
$I4_{1}/amd$ PdCuH2. The total DOS is given by the black line.

Figure 4 reports the electronic band structure and the density of states (DOS)
for both $P4/mmm$ and $I4_{1}/amd$ PdCuH2. Whereas the former contains one
formula unit (FU) in the primitive cell, the unit cell of the latter is
doubled, and therefore so is the total number of bands. Importantly, the DOS
at the Fermi level, $E_{\text{F}}$, in the lower-energy $I4_{1}/amd$ structure
is only half the value as in the $P4/mmm$ phase, 0.424 vs. 0.784 states eV-1
FU-1, respectively. A computational experiment where the DOS was calculated
for two hypothetical PdCuH2 structures: (i) with the same PdCu lattice as in
the $I4_{1}/amd$ phase but with the hydrogen atoms situated in the center of
the octahedral holes, and (ii) with the same distorted hydrogen lattice as in
the $I4_{1}/amd$ phase, but with the metal atoms arranged as in the $P4/mmm$
geometry, revealed that the metal atom placement was the main driver of the
decreased DOS in the ground state.

In both structures the contributions to the DOS at $E_{\text{F}}$ were the
lowest for hydrogen, followed by copper then palladium (16%, 31% and 53% for
$P4/mmm$ vs. 16%, 35% and 49% for $I4_{1}/amd$). The profile of the DOS for
both PdCuH2 structures resembles the one calculated for pure PdH [50].
However, the extra half electron per MH unit in PdCuH2 shifts the Fermi energy
further away from the bulk of the d states of both copper and palladium, with
a concomitant reduction of the DOS at $E_{\text{F}}$ in the ternary hydride by
almost a factor of two.

Figure 5: Phonon band structures. The phonon band structure, Eliashberg
spectral function, $\alpha^{2}F(\omega)$, and the integral of the electron
phonon coupling, $\lambda(\omega)$, for (a) $P4/mmm$ PdCuH2 and (b)
$I4_{1}/amd$ PdCuH2 calculated with classical nuclei (blue dashed lines and
shading) and incorporating the effects of anharmonicity and quantum nuclei
(red solid lines and shading).

The phonon spectra for both PdCuH2 phases is shown in Figure 5. Since the
$P4/mmm$ unit cell contains four atoms while eight atoms comprise the
$I4_{1}/amd$ unit cell, a total of 12 and 24 phonon modes, respectively, is
obtained. In both phases, the phonons can be divided into three distinct
branches: the six (or 12) lowest frequency modes are related to the acoustic
vibrations of both the Pd and Cu atoms, while the intermediate and upper
branches, having two (or four) and four (or eight) modes respectively,
correspond to optical vibrations linked to the hydrogen atoms.

The classical $P4/mmm$ PdCuH2 phonons (blue dashed curve) reveal an
instability at the $Z$ high-symmetry point and a softening at $X$, and they
are in good agreement with those previously calculated by Vocaturo and co-
workers [46]. Additionally, the intermediate optical branch associated with
the hydrogen atom motions shows signs of strong anharmonicity through soft
modes at $Z$ and $R$. The introduction of quantum nuclear effects (red solid
line) has a profound effect on the phonon band structure. Most notably, they
are able to heal the instability at $Z$ and the soft mode at $X$ associated
with the acoustic branches. Both the intermediate and upper optical branches
are subject to a strong red shift, which has a greater impact on the
intermediate branch that is pushed up by about 400 cm-1.

The phonon band structure for the lower symmetry $I4_{1}/amd$ PdCuH2 phase is
markedly different. Importantly, no instabilities are present in the harmonic
phonons and the introduction of quantum nuclear effects does not significantly
affect the vibrational frequencies associated with a particular mode,
suggesting the motions of atoms in this phase are mostly harmonic. The most
notable difference between the harmonic and the SSCHA-computed phonon modes
are the minor renormalizations of the intermediate and upper optical branches,
which undergo a blue shift separating them in energy.

Finally, the $T_{\text{c}}$ of these two phases was estimated via the Allen-
Dynes equation [51] (Table 1). For classical nuclei, and ignoring the
instability around the $Z$-point, the $T_{\text{c}}$ of $P4/mmm$ PdCuH2 was
estimated to be 37-44 K, in good agreement with previous estimates of 40-45 K,
also obtained by neglecting the imaginary phonon modes [46]. The large value
of $\lambda$ is related to the soft hydrogen optical modes that enhance the
Eliashberg spectral function at low frequencies, as seen in Figure 5(a).

When quantum nuclear effects are introduced, the integral of the spectral
function up to $\sim$260 cm-1 significantly decreases, and the region between
260-460 cm-1 falls to zero, shifting upwards in frequency, with a concomitant
decrease of $\lambda$ and increase of $\omega_{\textrm{log}}$. The resulting
$T_{\text{c}}$, which incorporates the quantum nuclear fluctuations, is
significantly reduced to 13-20 K, in-line with the maximum value of 17 K
experimentally reported for a presumed Pd1.1Cu0.9Hx, $x\sim 1.4$ phase [46].
About 75% of $\lambda$ derives from the part of the Eliashberg spectral
function associated with the motions of the hydrogen atoms, suggesting that
they provide the key contribution to superconductivity in this phase, similar
to previous results obtained for PdH [52]. However, in $P4/mmm$ PdCuH2 the
acoustic branches furnish a greater contribution towards $\lambda$ increasing
its value by $\sim$33% compared to PdH [52], while $\omega_{\textrm{log}}$
remains almost unchanged with just a small reduction of 8%. This suggests that
the more-than-doubled quantum $T_{\text{c}}$ of $P4/mmm$ PdCuH2 as compared to
PdH is likely due to the reduction in mass related to the substitution of half
of the palladium atoms with the lighter copper atoms.

Considering the similarity of PdCuH2 with PdH, we wondered if it might be
subject to the inverse isotope effect? Towards this end, we calculated the
phonon spectrum, Eliashberg spectral function (Supplementary Figure 2) and
resulting $T_{\text{c}}$ of the deuterated compound with and without the
SSCHA. While the $\lambda$ of $P4/mmm$ PdCuH2 and PdCuD2 were nearly
identical, the latter had a reduction of $\sim$15% in $\omega_{\textrm{log}}$,
as expected, resulting in a decreased $T_{\text{c}}$ of 11-17 K. Therefore, no
inverse isotopic effect was found.

The computed trends in the superconducting properties differ markedly for
$I4_{1}/amd$ PdCuH2, whose phonon band structure and Eliashberg spectral
function was not strongly perturbed by the introduction of quantum nuclear
effects. Assuming classical nuclei yields a $T_{\text{c}}$ of 5-10 K, whose
value is reduced slightly to 2-5 K via the introduction of quantum nuclear
effects. As compared to the $P4/mmm$ phase, the reduction, which results from
the slight blue shift of the optical phonon branches that moves
$\alpha^{2}F(\omega)$ to higher frequencies (Figure 5(b)), is not significant.
The lower $T_{\text{c}}$ of this phase as compared to $P4/mmm$ PdCuH2 is
likely due to the decreased value of the DOS at the Fermi energy.

Table 1:  Superconductivity. The values of the electron phonon coupling
parameter, $\lambda$, the logarithmic average frequency, $\omega_{\text{log}}$
(cm-1), and the superconducting critical temperature, $T_{\text{c}}$ (K), for
the specified systems calculated for classical (harmonic) and quantum (SSCHA)
protons. These values can be compared with experimental reports of a maximum
$T_{\text{c}}=$ 17 K in a PdCuHx compound whose stoichiometry and structure
were not determined [45].

| Harmonica |  | SSCHA
---|---|---|---
| $\lambda$ | $\omega_{\text{log}}$ | $T_{\text{c}}^{b}$ |  | $\lambda$ | $\omega_{\text{log}}$ | $T_{\text{c}}^{b}$
$P4/mmm$ PdCuH2 | 1.64 | 208 | 37-44 |  | 0.73 | 344 | 13-20
$I4_{1}/amd$ PdCuH2 | 0.57 | 351 | 5-10 |  | 0.48 | 349 | 2-5
PdHc | 1.55 | 205 | 34 |  | 0.40 | 405 | 5
$P4/mmm$ PdCuD2 | 1.65 | 163 | 29-34 |  | 0.74 | 291 | 11-17

a The imaginary frequencies in the harmonic phonons of $P4/mmm$ PdCuH2 and
PdCuD2 were neglected.
b $T_{\text{c}}$ calculated using the Allen-Dynes modified McMillan equation
with $\mu^{*}=$ 0.15 and 0.1.
c Reference [52].

In the experiments PdCu alloys were made by arc-melting in ultrapure argon,
but their atomistic arrangement was not reported [45]. Hydrogen was added to
the resulting PdCu metal foils by precharging under an H2 gas pressure of 4
atm at 300∘C followed by H${}_{2}^{+}$-ion implantation. It is therefore
likely that the hydrogen atoms were incorporated into octahedral holes already
present in the PdCu alloys, and the position of the metal atoms did not
rearrange during the hydrogenation.

The PdCu lattice of the $I4_{1}/amd$ structure is 11 meV/atom less stable than
the CSP-found lowest energy B2 phase, and also differs markedly from the
disordered fcc structure that can be found in the alloy. Therefore, it is very
unlikely that the hydrogenated alloy made in experiment contains regions that,
on a local scale, resembled $I4_{1}/amd$ PdCuH2. Thus, though this phase is
thermodynamically preferred, it is kinetically disfavored, explaining why the
measured maximum $T_{\text{c}}$ was significantly higher than the one
estimated for this phase. Similarly, the good agreement between the
experiments and the calculated $T_{\text{c}}$ of the $P4/mmm$ phase, which
could be derived from the B2 PdCu phase by placing hydrogen in its octahedral
holes, suggests that local regions with this motif were present in the
superconducting phase, whose stoichiometry was estimated as being
Pd1.1Cu0.9H1.4.

Calculations on supercell models derived from $P4/mmm$ PdCuH2 where some of
the hydrogen atoms are removed from the octahedral positions resulting in a
$\sim$Pd1.1Cu0.9H1.4 composition are directions for future work, as are
studies of metastable PdCu lattices that are subsequently hydrogenated. The
results presented here reveal a $\sim$4-fold enhancement of $T_{\text{c}}$,
depending upon the atomistic arrangement of the Pd and Cu atoms in the alloy,
which in turn affect the placement of the H atoms thereby influencing the DOS
at $E_{\text{F}}$, the vibrational properties and the electron phonon
coupling. It is conceivable that $T_{\text{c}}$ could be raised further, and
the atomistic arrangement of the alloy be engineered via synthesis routes far
from equilibrium. In addition, these avenues may prove fruitful for designing
hydrogenated AgPd or AuPd alloys with superconducting properties surpassing
those that have been previously reported in this class of compounds [45].

##  III Discussion

This study outlines an effective protocol that combines the Stochastic Self-
Consistent Harmonic Approximation (SSCHA) with the Machine Learning
Interatomic Potential (MLIP) package, enabling the automated and rapid
calculation of anharmonic and quantum nuclear effects for a wide range of
materials properties. In the proposed workflow DFT calculations are performed
on geometries (generated by the Monte Carlo algorithm of the SSCHA) that are
deemed to be extrapolative, and a subset of these is employed to train
quantum-accurate moment tensor potentials (MTPs) that are used to predict the
energies, forces and stresses of the remaining geometries required for the
SSCHA. Furthermore, we introduce an upscaling protocol, where the MTPs are
first trained on small cells, and retrained on larger ones, allowing for the
simulation of systems with sizes that were previously out-of-computational-
reach. The upscaling protocol is based on the assumption that most of the
interatomic interactions are already captured in the small cells, and that
only minor corrections to the MTPs are required when larger cells are
considered. This procedure opens the door for the routine modelling of quantum
and anharmonic effects, on cells of up to thousands of atoms, while at the
same time dramatically enhancing the computational efficiency.

The protocol is put to the test on the PdCuH2 system, which was hypothesized
[46] to be an important contributor to the superconductivity reported for
hydrogenated PdCu alloys with a maximum superconducting critical temperature,
$T_{\text{c}}$, of 17 K at ambient pressures [45]. Hydrogenation of the most
stable $Pm\bar{3}m$ symmetry PdCu phase would likely yield the kinetically
favored $P4/mmm$ PdCuH2 structure, metastable by 18 meV/atom from the newly
predicted $I4_{1}/amd$ PdCuH2 ground state. Similar to the highly anharmonic
PdH phase, $P4/mmm$ PdCuH2 is found to be dynamically stable only upon the
inclusion of quantum and anharmonic effects, which heal the imaginary modes in
its phonon band structure, and renormalize the frequencies of others. These
effects decrease the predicted $T_{\text{c}}$ of $P4/mmm$ PdCuH2 from 37-44 K
to 13-20 K, mirroring previous trends computed for PdH [52], whose anharmonic
$T_{\text{c}}$ is notably smaller ($\sim$5 K), likely because of the lighter
mass of the copper, as compared to the palladium, atoms. An inverse isotope
effect on superconductivity was not observed when hydrogen was replaced by
deuterium, with an estimated $T_{\text{c}}$ of 11-17 K for $P4/mmm$ PdCuD2.
The thermodynamically preferred $I4_{1}/amd$ PdCuH2 phase was shown to be
relatively well treated within the harmonic approximation, with the predicted
$T_{\text{c}}$ decreasing from 5-10 K to 2-5 K upon inclusion of
anharmonicity.

The good agreement between the computed $T_{\text{c}}$ of $P4/mmm$ PdCuH2 and
that measured for the hydrogenated alloy, with composition likely being closer
to Pd1.1Cu0.9H1.4, suggests structural similarities, on a local level.
Further, the results presented herein may have broader implications for the
discovery of ambient-condition superconducting hydrides, as many of the
recently predicted compounds involve ternaries featuring Pd (e.g. Mg2PdH6
[53]), Cu (e.g. Li2CuH6 [54]) and other combinations of alkaline earth and
transition metal atoms [55]. In addition This procedure opens the door towards
studies of anharmonicity using more costly electronic structure methods, such
as hybrid or meta-GGA functionals, which were previously out of reach.
Furthermore, it enables a determination of the relative stability of
enthalpically close phases, which could not be separated well by the PBE
functional as for example in the high pressure phases of lithium [56].

##  IV Methods

###  IV.1 Computational Details

The CSP searches were performed using the XtalOpt open-source evolutionary
algorithm [57, 58, 59, 60]. The initial generation consisted of random
symmetric structures created by randspg [61], and duplicate structures were
identified via XtalComp [61] and discarded from the breeding pool. The DFT
calculations for the CSP were performed with the Vienna _ab-initio_ Simulation
Package (VASP) version 5.4.1 [62, 63], while the calculations of the
electronic band structures, relative enthalpies, and the SSCHA were performed
with the Quantum ESPRESSO (QE) [64, 65] package. PAW and ultrasoft
pseudopotentials were used for VASP and QE, respectively, both within the
Perdew-Burke-Ernzerhof (PBE) [66] parametrization of the exchange-correlation
functional. The Pd 4d95s15p04f0, Cu 3d104s14p04f0, H 1s12p0, and Pd 5s15p04d9,
Cu 4s1.54p03d9.5, H 1s1 electronic configurations were treated explicitly for
VASP and QE, respectively. The CSP relaxations for each structure were
performed using a 4 step procedure where (i) the atomic positions, (ii) cell
parameters and (iii+iv) cell parameters and ionic positions were relaxed. In
these steps the energy cutoffs for the wavefunctions and a $\mathbf{k}$-point
grid density were set to 120, 240, 300 and 400 eV, and $2\pi\times 0.1$,
$2\pi\times 0.07$, $2\pi\times 0.05$, and $2\pi\times 0.025$ Å -1,
respectively. The single point QE calculations were performed imposing energy
cutoffs for the wavefunctions and density as 1224 eV and 12,240 eV,
respectively. The integration over the Brillouin zone in the self-consistent
calculations was performed with a first-order Methfessel-Paxton [67] smearing
of 0.272 eV coupled with a $\mathbf{k}$-point grid with density $2\pi\times
0.0175$ Å -1.

The electron-phonon coupling parameters were calculated via QE, on a k-point
grid with density $2\pi\times 0.00875$ Å -1, with a Gaussian smearing of
0.1088 eV, and a q-point grid with density $2\pi\times 0.0574$ Å -1. The
superconducting critical temperature was calculated using the Allen-Dynes
equation [51] and values of the Coulomb repulsion parameter, $\mu^{*}$, of 0.1
and 0.15.

In the SSCHA [25], the displaced configurations for the various supercells
used for the MTP training were generated at 300 K. The relaxation was
performed by fixing the size and shape of the cell, but allowing the atoms to
relax to their equilibrium positions. The MLIP package [37] was employed to
generate an MTP of level 10, with cutoff distances of $0.9<x<9$ Å . Every
structure with extrapolation grade, $\gamma_{\text{select}}>2$, became a
candidate for retraining the potential.

With this approach, potential training was considered complete within the
first 3 supercells, involving a total of 532 DFT configurations (140 in the
unit cell, 342 in the 2$\times$2$\times$2 and 50 in the 3$\times$3$\times$3
cells). The training performed at 300 K appears to be a key factor for the
fast training, as temperature increases the magnitude of the atomic
oscillations, thereby sampling areas of the potential that are further away
from the equilibrium position. This procedure, therefore, allows for a drastic
reduction in the need for further DFT calculations in SSCHA runs performed at
lower temperatures, since the smaller atomic vibrations are already contained
in the area already sampled at higher temperatures. The final phonon spectra
was computed at 0 K within the bubble approximation (Eq. 59 of Ref. [25]).

Code Availability. The code will be released in a subsequent publication, and
in the meanwhile, may be obtained from the corresponding authors upon
reasonable request.

Data Availability The datasets generated during and/or analyzed during the
current study are summarized in the supplementary information, and are
available from the corresponding author on reasonable request.

Acknowledgements Funding for this research is provided by the National Science
Foundation, under award DMR-2136038. Calculations were performed at the Center
for Computational Research at SUNY Buffalo
(http://hdl.handle.net/10477/79221).

Author Contributions The project was conceived and supervised E.Z. F.B.
developed the method and performed the calculations and analysis. F.B. wrote
the manuscript, which was edited by E.Z.

Competing Interests: The authors declare no competing interests.

Additional Information: Supplementary information The online version contains
supplementary material available at

## References

  * [1] Baroni, S., De Gironcoli, S., Dal Corso, A. & Giannozzi, P.  Phonons and related crystal properties from density-functional perturbation theory.  _Rev. Mod. Phys._ 73, 515 (2001).
  * [2] Wei, B., Sun, Q., Li, C. & Hong, J.  Phonon anharmonicity: a pertinent review of recent progress and perspective.  _Sci China Phys. Mech. Astron._ 64, 117001 (2021).
  * [3] Gutierrez-Amigo, M. _et al._ Purely anharmonic charge density wave in the two-dimensional dirac semimetal snp.  _Phys. Rev. B_ 109, 174112 (2024).
  * [4] Diego, J. _et al._ van der waals driven anharmonic melting of the 3d charge density wave in vse2.  _Nat. Commun._ 12, 598 (2021).
  * [5] Ranalli, L. _et al._ Temperature-dependent anharmonic phonons in quantum paraelectric ktao3 by first principles and machine-learned force fields.  _Adv. Quantum Technol._ 6, 2200131 (2023).
  * [6] Blancas, E. _et al._ Thermodynamics of solids including anharmonicity through quasiparticle theory.  _Chem. Rxiv._ (2024).
  * [7] Belli, F. & Errea, I.  Impact of ionic quantum fluctuations on the thermodynamic stability and superconductivity of labh 8.  _Phys. Rev. B_ 106, 134509 (2022).
  * [8] Hou, P., Belli, F., Bianco, R. & Errea, I.  Quantum anharmonic enhancement of superconductivity in p63/mmc sch6 at high pressures: A first-principles study.  _J. Appl. Phys._ 130 (2021).
  * [9] Hou, P., Belli, F., Bianco, R. & Errea, I.  Strong anharmonic and quantum effects in p m 3¯ n alh 3 under high pressure: A first-principles study.  _Phys. Rev. B_ 103, 134305 (2021).
  * [10] Errea, I. _et al._ Quantum crystal structure in the 250-kelvin superconducting lanthanum hydride.  _Nature_ 578, 66–69 (2020).
  * [11] Errea, I. _et al._ Quantum hydrogen-bond symmetrization in the superconducting hydrogen sulfide system.  _Nature_ 532, 81–84 (2016).
  * [12] Hellman, O., Abrikosov, I. & Simak, S.  Lattice dynamics of anharmonic solids from first principles.  _Phys. Rev. B_ 84, 180301 (2011).
  * [13] Hellman, O. & Abrikosov, I. A.  Temperature-dependent effective third-order interatomic force constants from first principles.  _Phys. Rev. B_ 88, 144301 (2013).
  * [14] Hellman, O., Steneteg, P., Abrikosov, I. A. & Simak, S. I.  Temperature dependent effective potential method for accurate free energy calculations of solids.  _Phys. Rev. B_ 87, 104111 (2013).
  * [15] Ljungberg, M. P. & Íñiguez, J.  Temperature-dependent classical phonons from efficient nondynamical simulations.  _Phys. Rev. Lett._ 110, 105503 (2013).
  * [16] Magdău, I. B. & Ackland, G. J.  Identification of high-pressure phases iii and iv in hydrogen: Simulating raman spectra using molecular dynamics.  _Phys. Rev. B_ 87, 174110 (2013).
  * [17] Wang, C., Chan, C. T. & Ho, K.  Tight-binding molecular-dynamics study of phonon anharmonic effects in silicon and diamond.  _Phys. Rev. B_ 42, 11276 (1990).
  * [18] Zhang, D.-B., Sun, T. & Wentzcovitch, R. M.  Phonon quasiparticles and anharmonic free energy in complex systems.  _Phys. Rev. Lett._ 112, 058501 (2014).
  * [19] Magdau, I. B. & Ackland, G. J.  Identification of high pressure phases iii and iv in hydrogen: Simulating raman spectra using molecular dynamics.  _Rev. Rev. B_ 87 (2013).
  * [20] Ceperley, O. M.  Path integrals in the theory of condensed helium.  _Rev. Mod. Phys._ 867, 279–355 (1995).
  * [21] Hooton, D.  Li. a new treatment of anharmonicity in lattice thermodynamics: I.  _Lond. Edinb. Phil. Mag_ 46, 422–432 (1955).
  * [22] Errea, I., Calandra, M. & Mauri, F.  Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides.  _Phys. Rev. B_ 89, 064302 (2014).
  * [23] Bianco, R., Errea, I., Paulatto, L., Calandra, M. & Mauri, F.  Second-order structural phase transition, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: theory and stochastic implementation.  _Phys. Rev. V_ 96 (2017).
  * [24] Bianco, R., Errea, I., Calandra, M. & Mauri, F.  High-pressure phase diagram of hydrogen and deuterium sulfides from first principles: Structural and vibrational properties including quantum and anharmonic effects.  _Phys. Rev. B_ 97, 214101 (2018).
  * [25] Monacelli, L. _et al._ The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects.  _J. Phys. Cond. Matt._ 33 (2021).
  * [26] Monacelli, L., Errea, I., Calandra, M. & Mauri, F.  Pressure and stress tensor of complex anharmonic crystals within the stochastic self-consistent harmonic approximation.  _Phys. Rev. B_ 98, 024106 (2018).
  * [27] Behler, J.  Perspective: Machine learning potentials for atomistic simulations.  _J. Chem. Phys_ 145 (2016).
  * [28] Bartók, A. P. _et al._ Machine learning unifies the modeling of materials and molecules.  _Sci. Adv._ 3, e1701816 (2017).
  * [29] Deringer, V. L., Caro, M. A. & Csányi, G.  Machine learning interatomic potentials as emerging tools for materials science.  _Adv. Mater._ 31, 1902765 (2019).
  * [30] Noé, F., Tkatchenko, A., Müller, K.-R. & Clementi, C.  Machine learning for molecular simulation.  _Annu. Rev. Phys. Chem._ 71, 361–390 (2020).
  * [31] Westermayr, J., Gastegger, M., Schütt, K. T. & Maurer, R. J.  Perspective on integrating machine learning into computational chemistry and materials science.  _J. Chem. Phys_ 154 (2021).
  * [32] Blank, T. B., Brown, S. D., Calhoun, A. W. & Doren, D. J.  Neural network models of potential energy surfaces.  _J. Chem. Phys_ 103, 4129–4137 (1995).
  * [33] Behler, J. & Parrinello, M.  Generalized neural-network representation of high-dimensional potential-energy surfaces.  _Phys. Rev. Lett._ 98, 146401 (2007).
  * [34] Unke, O. T. _et al._ Spookynet: Learning force fields with electronic degrees of freedom and nonlocal effects.  _Nat. Commun._ 12, 7273 (2021).
  * [35] Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  _Phys. Rev. Lett._ 104, 136403 (2010).
  * [36] Shapeev, A. V.  Moment tensor potentials: A class of systematically improvable interatomic potentials.  _Multiscale Model. Simul._ 14, 1153–1173 (2016).
  * [37] Novikov, I. S., Gubaev, K., Podryabinkin, E. V. & Shapeev, A. V.  The mlip package: moment tensor potentials with mpi and active learning.  _Mach. Learn. Technol._ 2, 025002 (2020).
  * [38] Merchant, A. _et al._ Scaling deep learning for materials discovery.  _Nature_ 1–6 (2023).
  * [39] Roberts, J. _et al._ Machine learned interatomic potentials for ternary carbides trained on the aflow database.  _arXiv preprint arXiv:2401.01852_ (2024).
  * [40] Leverant, C. J., Harvey, J. A. & Alam, T. M.  Machine learning-based upscaling of finite-size molecular dynamics diffusion simulations for binary fluids.  _J. Phys. Chem. Lett._ 11, 10375–10381 (2020).
  * [41] Cerasoli, D., F. Donadio.  Effective optimization of atomic decoration in giant and superstructurally ordered crystals with machine learning.  _Res. Sq_ (2024).
  * [42] Lucrezi, R. _et al._ Quantum lattice dynamics and their importance in ternary superhydride clathrates.  _Commun. Phys._ 6, 298 (2023).
  * [43] Lam, J., Abdul-Al, S. & Allouche, .  Combining quantum mechanics and machine-learning calculations for anharmonic corrections to vibrational frequencies.  _J. Chem. Theory Comput._ 16, 1681–1689 (2020).
  * [44] Skośkiewicz, T.  Superconductivity in the palladium-hydrogen system.  _Phys. Status Solidi B_ 59, 329–334 (1973).
  * [45] Stritzker, B.  High superconducting transition temperatures in the palladium-noble metal-hydrogen system.  _Z. Phys._ 268, 261–264 (1974).
  * [46] Vocaturo, R., Tresca, C., Ghiringhelli, G. & Profeta, G.  Prediction of ambient-pressure superconductivity in ternary hydride pdcuhx.  _J. Appl. Phys._ 131 (2022).
  * [47] Meziere, J. A. _et al._ Accelerating training of mlips through small-cell training.  _J. Mater. Res._ 38, 5095–5105 (2023).
  * [48] Subramanian, P. & Laughlin, D.  Cu-pd (copper-palladium).  _J. Phase Equilib. Diff._ 12, 231–243 (1991).
  * [49] Zhang, K. & Way, J. D.  Palladium-copper membranes for hydrogen separation.  _Sep. Purif. Technol._ 186, 39–44 (2017).
  * [50] Yang, X., Li, H., Ahuja, R., Kang, T. & Luo, W.  Formation and electronic properties of palladium hydrides and palladium-rhodium dihydride alloys under pressure.  _Sci. Rep._ 7, 3520 (2017).
  * [51] Allen, P. & Dynes, R.  Transition temperature of strong-coupled superconductors reanalyzed.  _Phys. Rev. B_ 12, 905 (1975).
  * [52] Errea, I., Calandra, M. & Mauri, F.  First-principles theory of anharmonicity and the inverse isotope effect in superconducting palladium-hydride compounds.  _Phys. Rev. Lett._ 111, 177002 (2013).
  * [53] Sanna, A. _et al._ Prediction of ambient pressure conventional superconductivity above 80 k in hydride compounds.  _Npj Comput. Mater._ 10, 44 (2024).
  * [54] Cerqueira, T. F. T., Fang, Y. W., Errea, I., Sanna, A. & Marques, M. A. L.  Searching materials space for hydride superconductors at ambient pressure.  _arXiv preprint arXiv:2403.13496_ (2024).
  * [55] Dolui, K. _et al._ Feasible route to high-temperature ambient-pressure hydride superconductivity.  _Phys. Rev. Lett._ 132, 166001 (2024).
  * [56] Racioppi, S. _et al._ Phase boundaries, isotope effect and superconductivity of lithium under hydrostatic conditions.  _arXiv preprint arXiv:2312.17498_ (2023).
  * [57] Lonie, D. C. & Zurek, E.  Xtalopt: An open-source evolutionary algorithm for crystal structure prediction.  _Comput. Phys. Commun._ 182, 372–387 (2011).
  * [58] Avery, P., Toher, C., Curtarolo, S. & Zurek, E.  Xtalopt version r12: An open-source evolutionary algorithm for crystal structure prediction.  _Comput. Phys. Commun._ 237, 274–275 (2019).
  * [59] Xtalopt evolutionary crystal structure prediction.  http://xtalopt.github.io/ (2024).
  * [60] Falls, Z., Avery, P., Wang, X., Hilleke, K. P. & Zurek, E.  The xtalopt evolutionary algorithm for crystal structure prediction.  _J. Phys. Chem. C_ 125, 1601–1620 (2020).
  * [61] Avery, P. & Zurek, E.  Randspg: an open-source program for generating atomistic crystal structures with specific spacegroups.  _Comput. Phys. Commun._ 213, 208–216 (2017).
  * [62] Kresse, G. & Hafner, J.  Ab initio molecular dynamics for liquid metals.  _Phys. Rev. B._ 47, 558–561 (1993).
  * [63] Kresse, G. & Joubert, D.  From ultrasoft pseudopotentials to the projector augmented–wave method.  _Phys. Rev. B._ 59, 1758–1775 (1999).
  * [64] Giannozzi, P. _et al._ Quantum espresso: a modular and open-source software project for quantum simulations of materials.  _J. Condens. Matter Phys._ 21, 395502 (2009).
  * [65] Giannozzi, P. _et al._ Advanced capabilities for materials modelling with quantum espresso.  _J. Condens. Matter Phys._ 29, 465901 (2017).
  * [66] Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Phys. Rev. Lett._ 77, 3865 (1996).
  * [67] Methfessel, M. & Paxton, A.  High-precision sampling for brillouin-zone integration in metals.  _Phys. Rev. B_ 40, 3616 (1989).