Annealing for prediction of grand canonical crystal structures: Efficient implementation of n-body atomic interactions

We propose an annealing scheme usable on modern Ising machines for crystal structures prediction (CSP) by taking into account the general n-body atomic interactions, and in particular three-body interactions which are necessary to simulate covalent bonds. The crystal structure is represented by discretizing a unit cell and placing binary variables which express the existence or non-existence of an atom on every grid point. The resulting quadratic unconstrained binary optimization (QUBO) or higher-order unconstrained binary optimization (HUBO) problems implement the CSP problem and is solved using simulated and quantum annealing. Using the example of Lennard-Jones clusters we show that it is not necessary to include the target atom number in the formulation allowing for simultaneous optimization of both the particle density and the configuration and argue that this is advantageous for use on annealing machines as it reduces the total amount of interactions. We further provide a scheme that allows for reduction of higher-order interaction terms that is inspired by the underlying physics. We show for a covalently bonded monolayer MoS2 crystal that we can simultaneously optimize for the particle density as well as the crystal structure using simulated annealing. We also show that we reproduce ground states of the interatomic potential with high probability that are not represented on the initial discretization of the unit cell.


INTRODUCTION
Crystal structure prediction (CSP) from chemical composition alone is still one of the most difficult problems in materials science, even for the simplest structures [1].The reason why this problem is still a challenge is that the variation of possible structures grows exponentially as the number of atoms increases, making an exhaustive search for the most stable structure, i.e. finding the global minimum on the Born-Oppenheimer surface, unfeasible even with today's supercomputers.For a small number of atoms, a brute force approach is possible, but reliably finding global optima for larger systems is out of reach of current computers.
In recent years, the use of quantum computers has attracted a great deal of attention as a means of searching for globally optimal solutions [25][26][27][28][29][30].Quantum computers are characterized by their ability to escape from locally stable solutions and accelerate the search for globally optimal solutions by utilizing the quantum tunneling effect [27,31].Quantum annealing (QA) machines [25,[32][33][34][35] and gate-based quantum computers are the two main current architectures in development.Exhaustive structure search using gate-based quantum computers has been reported recently [36,37].In the method described in Ref. [37], space is divided into meshes, and the presence or absence of atoms on each mesh is represented as a {0,1} digital number, allowing the crystal structure to be encoded onto qubits as a bit sequence.On the qubits, various atomic coordination structures can be prepared at once by using the quantum superposition states on the qubits.The idea is to perform exhaustive structural optimization by applying a probabilistic arXiv:2307.03123v2[quant-ph] 11 Sep 2023 imaginary-time evolution technique reported in [30].
In this paper, we report a method to perform exhaustive structural optimization using QA.In particular, we discuss how to reformulate structural optimization as a quadratic unconstrained binary optimization problem (QUBO) or higher-order unconstrained binary optimization (HUBO).We provide a scheme for implementing an empirical three-body interatomic potential on QA hardware, and we provide a detailed analysis of preliminary SA and QA results.In particular we argue that providing more physical information in the form of penalty terms does not necessarily speed up the computation.
The remainder of the paper is structured as followed.In Section II we present the HUBO formulation for the CSP.In Section III we introduce the methods and general parameters used for optimization.In Section IV we outline the parameters for a Lennard-Jones cluster of Krypton atoms for which we optimized both structure and density using SA and QA.In Section V we present a covalently bonded MoS 2 crystal modeled by a Stillinger-Weber potential for which we optimized again the structure and density using SA.We then close with the conclusions in Section VI.

II. HUBO FORMULATION
In this section we discuss the construction of our HUBO.In Section II A we discuss the notation of our unit cell discretization and the encoding into a HUBO of the CSP.In Section II B we disuss the penalty terms we use and finally in Section II C we discuss a physicallymotivated scheme to reduce the interaction terms of interaction terms of order higher than quadratic.

A. CSP problem encoding and Hamiltonian
Consider a unit cell that is spanned by a given basis {⃗ a i } with periodic boundary conditions along a chosen set of basis vectors and a set of atom species S. We look at a set of N lattice points X in this unit cell generated by partitioning each basis vectors into g +1 points and forming the corresponding lattice.The lattice points have the form i ki g ⃗ a i where k i ∈ {0, . . ., G i } with G i = g if we have no periodic boundary conditions along ⃗ a i and G i = g − 1 otherwise.Consider a set b s x of binary variables that we define such that if b s x = 1 there is an atom of species s ∈ S on x ∈ X. Assume that we have a set of potential functions V s1,...,sm m (x 1 , . . ., x m ) for a configuration of atoms of species s i on x i for m ∈ {1, . . ., M }.As is usual for interatomic potential functions we assume that it does not depend on the order in which the argument, species pairs are supplied, i.e.
for any permutation σ.Assuming that we have no periodic boundary conditions, we define our Hamiltonian as where the prime indicates that the x i ∈ X should be chosen such that x i ̸ = x j for any pair i, j (the species are chosen freely).Defined as such, finding the optimal nuclear structure on the lattice X corresponds to finding an optimal binary string that minimizes this Hamiltonian, as energy contributions only arise if all binary variables involved in an interaction are 1, i.e. all atoms involved in the interaction are present.Generalising this to the case with periodic boundary conditions requires a careful consideration of the selfinteractions of atoms with their periodic images and a fitting definition of the Hamiltonian.This is done in detail in Appendix A.

B. Penalty terms
Eq. (2) allows us to calculate the cohesive energy of a given configuration (see Appendix B).Thus, for wellconstructed interatomic potentials that accurately model a wide range of configurations of a material, Eq. ( 2) not only gives the optimal configuration, but by simultaneously finding the optimal amount of binary variables that should have the value 1 we optimize for the optimal density of atoms in the unit cell.
It is possible to a priori fix a target atom number in the unit cell by adding a penalty term such as to the Hamiltonian for an appropriately large positive P and all s ∈ S where C s is the target particle number for species s atoms.We call this an absolute penalty term.
Equivalently, knowing the chemical formula (e.g.Al 2 (SO 4 ) 3 ) but not the optimal density, a penalty term like ensures that the ratios of atoms are respected, where c s1,s2 is the target ratio (in the above example c S,O = 1/4).This penalty term allows for finding the optimal density in the range that the ratio is respected.We call this a relative penalty term.

C. Reduction of interaction terms
Interatomic potentials will usually include a cutoff distance.To reduce pairwise interaction terms, it is crucial to choose the right penalty terms as an absolute penalty term will introduce interactions between any pair of binary variables for the same species, even if their pairwise distance is higher than the cutoff distance.Similarly, relative penalty terms introduce pairwise interactions between any pair of binary variables of the two involved species.Choosing the wrong penalty terms can make the difference between having a sparse or fully connected graph of pairwise interactions.Ideally, no penalty terms would be introduced, but this is dependent on the quality of the chosen potential.
The number of interaction terms for the cubic or higher order terms in the HUBO will be orders of magnitudes higher than for the pair interactions.Often, alongside the total number of spins, the density of the interaction graph is the main bottleneck for modern annealing machines [35,38] and as such it is crucial to devise schemes that reduce the interaction number beyond just applying a cutoff.To this end we use the 'deduc-reduc' method from [39].In particular, we make the assumption that if the pairwise interaction between two binary variables is higher than a user-set threshold T , then any higher-order interaction containing this pair can safely be set to 0 without influencing the ground state.At the same time we replace any pairwise interaction J ij by min(J ij , T ).The intuition behind this is that for the interatomic potentials we use in this work, the pairwise interaction rapidly increases if the atoms are too close, and thus the ground state does not contain atoms on the two involved locations and we do not need to evaluate the higher-order terms.This is a simplification that does not loose any generality with respect to the ground state of the HUBO and which in particular also does not require any a priori knowledge like atomic radii of the involved species.

III. METHODS
We will find optimal binary strings for the HUBO problems using SA and QA.In this section we outline the notation, parameters and settings we used for the optimization.

A. Simulated Annealing
Simulated annealing is a classic algorithm for optimizing cost functions with several local minima [40].We assume some basic knowledge of the algorithm and will only discuss the specifics of our implementation.We use a geometric cooling schedule where T min and T max are the minimum and maximum temperature.The number of steps N steps is the number of Monte Carlo steps per spin to perform.Choosing the right neighbourhood for a configuration in SA (i.e.defining legal transitions of the Markov chain) is crucial and generally one aims to have a smooth energy landscape with not too rugged local minima [41][42][43].Traditionally, SA for HUBOs performs single bit flips.As this is equivalent to removing or adding an atom from the configuration, especially in the presence of penalty terms, this can be a costly operation.Thus, for each step in the schedule we loop over every binary variable and attempt to flip it and then we loop over every opposite valued pair in the current configuration and attempt to exchange their values.This latter flip moves an existing atom to a random location and does not break penalty terms such as the absolute penalty (3) or relative penalty (4), thus ensuring a smoother energy landscape.So when we speak of Monte Carlo steps per spin we mean that we attempt

B. Quantum Annealing
We also assume familiarity with the basic concepts of quantum annealing [25,27].We use the Advantage system available through the D-Wave leap cloud service [44].Our HUBO and QUBO problems are very densely connected and if the cutoff of the potential function is large enough or the system small enough, the problem might even be fully connected.Embedding these onto the Pegasus architecture of the Advantage system [35] requires us to calculate a minor embedding [45][46][47].Instead of manually calculating an embedding best fit for our problem, we use the standard implementation for clique embedding in the D-Wave Ocean SDK.This procedure can lead to results with broken chains which require a fitting unembedding.While there is evidence that designing a problem specific unembedding algorithm [48] can be advantageous we choose the simple majority vote which sets the binary value of a chain to the one that occurs most often on the chain.

C. Benchmarking
For benchmarking the various optimization schemes for the HUBO and QUBO formulation we use the time-to-FIG.1: The target FCC configuration of the Krypton system with Krypton atoms in pink (graphics due to Vesta).The solid atoms on the origin and the three incident face centers are the locations encoded in the HUBO while the remaining transparent ones are copies due to the periodic boundary conditions and not part of X.
solution [49,50] given by where τ is the running annealing time as measured on the local machine and P GS (τ ) is the probability of the corresponding algorithm to return the ground state with a running time of τ .The time-to-solution can be understood as the average time it takes to get the ground state with probability p r which we set to 0.99.

IV. KRYPTON SYSTEM
In this section we introduce an LJ cluster system consisting of Krypton atoms in Section IV A and the related SA and QA results in Section IV B.

A. Setup
For the calculation of the potential functions we rely on the Open Knowledgebase of Interatomic Models (OpenKIM) [51].In particular we will look at a three dimensional cubic unit cell of side length 5.653 Å with the Lennard-Jones potential parameters due to Bernades for Krypton [51][52][53][54][55] and periodic boundary conditions along all three basis vectors.We will look for the ground state configuration of Krypton atoms in this unit cell discretized into a equipartitioned lattice of size g 3 , which is equal to the face-centered cubic configuration and can be seen in Fig. 1.The energy of the FCC configuration is −0.431eV and for any interaction value J ij we take min(J ij , 1eV).While for the SA calculations this is not strictly necessary, it helps for the QA calculations as the energy range is normalised to be between 0 and 1 on D-Wave machines, thus upper bounding the energy ensures that the physically interesting energy range takes up a larger portion of the renormalised energy range.We will FIG.2: The SA time to solution results for the Krypton system with a penalty term in blue crosses and without in orange plus-symbols plotted against various grid granularities g.The solid line corresponds to a fit of the measured points to a(N + N 2 ) where a = 21.015 is the fitting parameter.
simply refer to this system as the Krypton system.We perform SA calculations without any penalty terms and with an absolute number penalty term setting C Kr = 4, we call the former grand canonical and the latter microcanonical.As the unit cell is smaller than the cutoff distance of the potential, even the grand canonical calculation QUBO is fully connected.We use a penalty strength of P = 1, and vary the temperature from 10 −2 to 10 −4 .The various probabilities correspond to the measured probability across 1000 annealing runs.
Since the systems are fully connected, for the QA calculations, we simplify the QUBO by fixing the binary variable for the origin to be 1 and removing any binary variable that had an interaction with the origin of more then 1eV.This corresponds in essence to removing the translational invariance of the problem.Further, we use pausing [56,57].We use a base length of the schedule of 20µs and we pause for 3µs.We consider the success probability, i.e. the ratio of obtained ground states over 40000 annealing runs, plotted against the pause location s p ∈ (0, 1) so that the dimensionless time in the annealing schedule goes from 0 to s p at (17 • s p )µs until (17 • s p )µs + 3µs and then goes to 1 linearly until 20µs.We use a chain strength of 1.28.These parameters were heuristically found to provide reasonable results.

B. Results and discussions
In Fig. 2 we plot the TTS against various grid spacings g for SA calculations for the grand and microcanonical system.We performed SA until we found the ground state FCC configuration with a probability of more than 90% and take the minimum TTS across the schedule FIG.3: The histogram for the residual energy of the Krypton system after running SA for 3 Monte Carlo steps per spin and the various g together with their average residual energy, i.e. energy above the ground state, ⟨H⟩.This is the full histogram, no results have been cut.
steps as the data point for g.This takes at most 30 schedule steps for both systems and it is apparent that both systems have comparable performance.In particular note the fit to the function N + N 2 , which is the scaling of the number of flips the SA algorithm attempts with the spin number N .There are two main mechanisms that increase the required TTS.The first is that, as we attempt more spin flips per schedule step with increasing g, SA requires more time per schedule step to perform the increasing amount of flips.The second is that with increasing g the atoms have more fine-grained displacement possibilities so that there are more local minima of the QUBO problem with energies closer to the actual ground state leading to an increased time to escape the local minima to find the ground state.
If the global minimum were harder to find due to increasing amounts of local minima, we would expect an increasing number of required schedule steps with increasing g.What we see is that the fit a(N + N 2 ) with a constant a = 21.015reconstructs the data well for g ≥ 12 for both systems.Thus there is no significant scaling ∼ T T S(τ )/(N + N 2 ) of the required scheduled steps with g for the microcanonical and the grand canonical system.
Further, in Fig. 3 we show a representative energy histogram for the grand canonical calculations with 3 Monte Carlo steps per spin for g ∈ {12, 14, 16, 18, 20}.Despite not putting any particle number restrictions the annealing process, even for this low amount of schedule steps, only returns solutions with the correct atom density and in fact all returned energies are lower than the first excited state energy corresponding to an FCC configuration with an atom taken out (see Fig. 8 in the appendix), a state we call FCC−1.Using the Broy-FIG.4: Ground state probabilities for the g = 4 Krypton system using the D-Wave Advantage 4.1 system with various pause locations s p ranging from 0.01 to 0.8.In blue the grand canonical calculation and in orange the microcanonical with a penalty strength of 0.05.The dashed lines correspond to the ground state probability after applying BFGS on the results from the solid lines and the dotted line to the probability of running annealing with no pauses and an annealing time of 18.9µs.
den-Fletcher-Goldfarb-Shanno (BFGS) algorithm [58][59][60][61] to converge to a local minimum off the grid X we confirmed that all states with 4 atoms converge to the ground state meaning that the TTS of the combination of annealing combined with BFGS is considerably lower than that of only annealing.
We also confirmed these tendencies on the D-Wave Advantage 4.1 system available on D-Wave Leap.We performed calculations only for the g = 4 system since the minor embedding for the g = 6 system had chain lengths of up to 20 spins which proved too hard to optimize.In Fig. 4 we plot the pause location s p against the success probability for the grand and microcanonical system for just QA with pausing and without pausing and a schedule length of 18.9µs and quantum annealing with pausing followed by BFGS.The penalty strength in the microcanonical calculations is 0.05 as it provided the best ground state probability.First we see that pausing improves the performance as for both systems the probabilities without pausing are around 0.001 and with pausing the maximum probabilities for the grand canonical system are 0.0067 at s p = 0.45 and 0.005425 at s p = 0.34 for the microcanonical one.Since there are no samedensity local minima, performing BFGS optimization on the results with pausing, is equivalent to looking at the results that have the correct density.We see that for QA+BFGS calculations both systems have success probabilities between 0.15 and 0.22 with the grand canonical consistently having a higher probability.
Without pausing QA has a TTS of around 0.9 * 10 5 µs comparable with the microcanonical system TTS for SA in the g = 4 case (see Fig. 2).With pausing we find a TTS of 13700µs and 16931µs respectively for the grand and microcanonical system providing comparable times to the grand canonical SA calculations albeit the QA calculations are a bit slower.Thus we find no indications of a quantum speedup.Possible reasons for this result may include the embedding of full connectivity on the sparse hardware graph and noise effects.We leave it for future research to analyse this problem with a wider set of parameters and using more intricate embedding techniques.Note though, that while the SA+BFGS algorithm did not provide any other minima than the global one, QA+BFGS returns the FCC-1 configuration with probabilities between 0.3 and 0.33 across all pause locations s p for the grand canonical system and 0.21 and 0.25 for the microcanonical system.Thus while we might not expect a quantum speedup there might be an advantage due to the higher breadth of results returned by QA compared to SA allowing a wider exploration of the potential energy landscape.
Summarizing, we see that also for QA, at least in this very simple system, there are no performance costs in leaving out the penalty and in fact we can expect performance increases confirming the tendencies found in SA.

V. MOS2 SYSTEM
In this section we introduce a MoS 2 system governed by the three-body Stillinger-Weber potential in Section V A and the related SA results in Section V B.

A. Setup
For the second system we consider the Stillinger Weber potential [62,63] which is a simple three-body po-tential that reflects covalent bond dynamics.We use the parametrization for hexagonal monolayer Molybden-Disulfide due to Wen et al. [64][65][66][67].We do this on the supercell consisting of a 2 × 2 lattice of hexagonal lattice unit cells with a single unit cell having a lattice constant of 3.20 Å and thickness of 3.19 Å.Thus the lattice vectors for our system are ⃗ a ), ⃗ a 3 = (0, 0, 3.19 Å) We build the lattice by partitioning both ⃗ a 1 and ⃗ a 2 into g = 6 equal parts each and applying periodic boundary conditions and partioining ⃗ a 3 into three equal parts without periodic boundary conditions.Thus the amount of required bits scales like 6g 2 , where the additional 2 comes from the amount of species.The target ground state is the 2H configuration (see Fig. 5) and has an energy of −55.5283eV.The first excited state that we expect to see is the 1T configuration, with the same amount of atoms and an energy that is 1.4755eV above the ground state (see Fig. 9 in the appendix).We will refer to this system as the MoS 2 system.
We use our deduc-reduc with a threshold of 10eV which in this particular case reduced the amount of non-zero three-body interaction terms by 18.8% (from 1573728 to 1277267) in the g = 6 system.Any lower threshold seemed to impact the ground state configuration on our SA calculations.There is no general-use scheme known to the authors, that would allow to quadratize this HUBO so as to make it runnable on any modern Ising machine [68][69][70] and so while our deduc-reduc step reduces the interactions it can only be a first step in conjunction with other approaches yet to be found and we perform no QA for this system.
We perform SA for the system with both absolute penalty terms (C Mo = 4, C S = 8) and relative penalty terms (C Mo,S = 1/2).For simplicity we call the former the absolute system and the latter the relative system.Grand canonical calculations as in the Krypton system without penalty terms do not work for this potential, as it is more favourable to produce configurations with a single atom species rather than a MoS 2 mix, so we limit our analysis to the relative and absolute system and recall that the former retains the function of simultaneously optimizing for the atom density.The number of pairwise interaction terms without interactions increases by 1.2% using the absolute penalty (from 21420 to 21708) and by 8.4% using the relative penalty (23220), underlining again the importance of finding potentials that can be used without penalties to reduce the number of pairwise interactions necessary.In fact, since this potential is parametrized for hexagonal MoS 2 we cannot expect it to yield accurate results for non-hexagonal configurations.This is a problem that does not pertain to the parametrization but the Stillinger-Weber potential in general.Since this one of the simplest three-body potentials we use it anyway for this proof-of-concept calculation.
We use a penalty strength of P = 10 and a temperature range of 10 to 0.1 for SA.The various probabilities FIG.6: Plot in solid lines of the ground state probability for SA for the MoS 2 system with schedule steps going from 2 to 500 for both relative penalties and absolute penalties (blue and orange respectively).The scale for the probability is to the left.In dashed lines the average residual energy with the corresponding scale to the right.
correspond to the measured probability across 1000 annealing runs.

B. Results and discussions
The MoS 2 system proves harder to optimize than the Krypton system.In Fig. 6 the ground state probabilities for schedule steps going from 2 to 500 are plotted.As opposed to the Krypton system where even for g = 20 we need only 30 schedule steps to reach a ground state probability of above 0.9 we see that it hovers around 0.4 for the absolute penalty and around 0.15 for the relative penalty at 500 schedule steps.In particular note that here the used penalty terms have an effect on the ground state probability, and that supplying more information (in form of the absolute penalty) leads to higher ground state probabilities.As expected the ground state probability increases with increasing amount of schedule steps but the slope does not offset the increase in calculation length and so the TTS turns out to be minimized for a number of schedule steps in the single digits for both system.In Fig. 6 the average residual energies are plotted and we see that both systems seem to converge to an average residual energy that is well above the target 0eV.To understand this, consider the energy histogram in Fig. 7 for the resulting states of only SA (top) and SA followed by BFGS with the same potential (bottom) after 500 schedule steps.First, note that despite not fixing an absolute number of atoms in the relative penalty, we find the correct density of Mo 4 S 8 in 42.8% of the configurations (in green in Fig. 7) and that the average residual energy for the states with the correct density is FIG.7: Histogram of the residual energy for the MoS 2 system with SA with 500 schedule steps (top) and SA + BFGS (bottom) applied to the MoS 2 system with the absolute number penalty Eq. ( 3) in blue and the relative number penalty Eq. ( 4) in orange for results with sub-optimal density and green for the optimal density.Found local minima are marked by a dotted line and the shaded area to the left (see Appendix C for the configurations).This is not the full histogram, i.e. there are configurations with energies higher than 10eV.
2.3826eV while it is 10.6117eV for the states with the wrong density (in orange) so that the relative penalty calculations allow for simultaneous optimization of the atom density and the optimal configuration.The probability to obtain either 2H or 1T configurations is 42% for the absolute penalty system and 18.9% for the relative penalty system.To understand the physical nature of the remaining local minima, which form the majority of found states, we performed BFGS on all the resulting states from SA.While the probability for 2H and 1T rose to 42.8% and 20.9% for the absolute penalty and relative penalty system respectively we see that most states converge to a local minimum that has an energy below that of 2H.First, for the relative penalty system we see that 57.2% of all observed configurations have 5 Molybden atoms and 10 Sulfur and form configurations that have an energy that is more than 2.5eV lower than that of 2H.In Fig. 7 we only shade the region as the BFGS algorithm does not converge well for these configurations so that we do not get well formed peaks but rather a distribution in the shaded area.The next lower state is a state we call orthorombic (see Appendix C for an image of both the orthorombic and an example Mo 5 S 10 configuration) and has an energy that is 0.9313eV lower than that of 2H.We find this configuration with a probability of 21.8% for the relative penalty system and 57.2% for the absolute penalty system.
Using the Vienna ab initio simulation package [71][72][73] with the projector augmented-wave method [74,75] we find that the energy of the 2H configuration is in fact the lowest of the four found local minima, followed by the 1T, the orthorombic and finally the Mo 5 S 10 configurations.The fact that this order is not represented is due to the fact that the Stillinger-Weber potential is parametrized to model hexagonally ordered MoS 2 configurations and thus does not correctly model other configurations.The potential is not fit to provide new physical insights in our application and these results should be taken merely as a proof of concept.
Noteworthy about these results is that, despite the orthorombic and locally optimal Mo 5 S 10 states not being representable on the discretization of the unit cell, the combination of SA and BFGS managed to find these states in a majority of attempts.This is a strong indication that if we are able to provide a fitting potential or directly a fitting HUBO we can find a wide array of globally and locally optimal configurations even if they are not part of the initial discretization.Thus, in particular it might suffice to have rougher discretizations with spin numbers that fit onto current quantum hardware instead of trying to be fine grained enough to represent all possible local minima.

VI. CONCLUSIONS
In this paper we have presented an annealing scheme for crystal structure prediction based on n-body atomic interactions.We discretized a given unit cell with a lattice and placed binary variables on the lattice points to express the existence or non-existence of an atom at every grid point.In particular this is done for 3-body atomic interactions which is the minimum order necessary for covalent crystals.We solved the resulting HUBOs using SA and QA giving insights into the crystal structure.We have shown that a grand canonical calculation without penalty terms allows for the simultaneous optimization of both the nuclear structure as well as the particle density inside the unit cell.Further, we have also shown evidence that the difficulty of solving the nuclear structure problem does not necessarily scale with the mesh size.These results show that it might not always be advantageous to put all the available information into the QUBO to speed up calculations in particular as this also increases the amount of total interaction terms the reduction of which is crucial for embedding problems into modern hardware with limited graphs.
We also considered a Molybden-Disulfide monolayer system modeled by a three-body Stillinger-Weber type potential.
Using our interaction number reduction scheme we reduced the amount of cubic interactions by 18.8% while maintaining physical accuracy to the extent of the used potential.We have shown that the potential contained unphysical ground states that are due to the limited transferability of the potential outside the context of hexagonal monolayer MoS 2 .While these results do not provide physical insights, we show that our algorithm reproduces the ground state of the system even if they are not representable on the chosen discretization of the unit cell in the annealing step of the algorithm.Thus, while we could only optimize the roughest discretization for the Krypton system on the D-Wave quantum annealer, this could be a hint that rougher discretizations, that are easier to embed onto quantum annealers, are enough for the local optimization algorithm to find a wide array of ground state and locally optimal configurations.
An immediate future research question is to choose a more fitting potential to construct a HUBO that accurately models a wide array of covalent crystal configurations to test the performance with rough unit cell meshes on larger unit cells.
Another research direction is to investigate the nature of returned local minima by QA and to confirm the tendency we found where QA provided a more varied insight into the energy than SA which tended to favour only ground states.
Note added.During the writing of this manuscript we have become aware of a similar proposal for the construction of the QUBO [76] for ionic crystals.That paper does not address higher-order optimization problems and thus does not address covalent bonds and did not consider the grand canonical case, their focus is on classical computation and providing guarantees that ground truths to the crystal structure prediction problem are found using their algorithm.They have similar findings with respect to the reproducibility of the ground state even if it is not contained in the initial discretization.We give an overview of the local minima indicated by dotted lines in the histograms Figs. 3 and 7.The local minima for the Krypton system are given in Fig. 8 and for the MoS 2 system in Fig. 9.

2
spin flips where N • |S| is the (unreduced) binary variable number.

FIG. 5 :
FIG. 5: The target 2H ground state configuration of the MoS 2 system with Sulfur in yellow and Molybden in violet (graphics due to Vesta).The bottom six sulfur atoms on the boundary, two at (⃗ a 1 + ⃗ a 2 )/2 and four Molybden atoms with z-coordinate given by ⃗ a 3 /2 are the locations encoded in the HUBO (non transparent atoms).The remaining ten Sulfur atoms (transparent) are copies due to the periodic boundary conditions and not part of X.

FIG. 8 :FIG. 9 :
FIG.8: Kr 3 configuration that corresponds to an FCC configuration with a single atom taken out and which has a residual energy of 0.2029eV