Annealing for prediction of grand canonical crystal structures: 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 nonexistence 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 ${\mathrm{MoS}}_2$ 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.

Figure 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 $\mathrm{X}$.

Figure 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+\left(\genfrac{}{}{0pt}{}{N}{2}\right)]$, where $a=21.015$ is the fitting parameter.

Figure 3

The histogram for the residual energy of the krypton system after running SA for three Monte Carlo steps per spin $g\in \{12,14,16,18,20\}$ [colored from light red (gray) to dark purple (gray)] together with their average residual energy, i.e., energy above the ground state, $\langle H\rangle$. This is the full histogram, no results have been cut.

Figure 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 (dark gray) the grand canonical calculation and in orange (light gray) 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\mu \mathrm{s}$.

Figure 5

The target $2\mathrm{H}$ ground-state configuration of the ${\mathrm{MoS}}_2$ system with sulfur in yellow and molybdenum in violet (graphics due to Vesta). The bottom six sulfur atoms on the boundary, two at $({\vec{a}}_1+{\vec{a}}_2)/2$ and four molybdenum atoms with $z$ coordinate given by ${\vec{a}}_3/2$ are the locations encoded in the HUBO (nontransparent atoms). The remaining ten sulfur atoms (transparent) are copies due to the periodic boundary conditions and not part of $X$.

Figure 6

Plot in solid lines of the ground-state probability for SA for the ${\mathrm{MoS}}_2$ system with schedule steps going from 2 to 500 for both relative penalties and absolute penalties [blue (dark gray) and orange (light gray), respectively]. The scale for the probability is to the left. In dashed lines the average residual energy with the corresponding scale to the right.

Figure 7

Histogram of the residual energy for the ${\mathrm{MoS}}_2$ system with SA with 500 schedule steps (top) and SA $+$ BFGS (bottom) applied to the ${\mathrm{MoS}}_2$ system with the absolute number penalty Eq. (3) in light red (light gray) and the relative number penalty Eq. (4) in dark purple (dark gray) for results with suboptimal density and red (gray) for the optimal density. Found local minima are marked by a dotted line and the shaded area to the left (see Appendix pp3 for the configurations). This is not the full histogram, i.e., there are configurations with energies higher than $10\mathrm{eV}$.

Figure 8

${\mathrm{Kr}}_3$ configuration that corresponds to an fcc configuration with a single atom taken out and which has a residual energy of $0.2029\mathrm{eV}$.

Figure 9

Local minima of the ${\mathrm{MoS}}_2$ system marked with a dotted line in Fig. 7 from the main text. From left to right, (a) an example ${\mathrm{Mo}}_5{\mathrm{S}}_{10}$ configuration with a residual energy of $-6.2161\mathrm{eV}$, (b) the orthorhombic state with a residual energy of $-0.9313\mathrm{eV}$, and (c) the $1\mathrm{T}$ configuration of ${\mathrm{MoS}}_2,1.4755\mathrm{eV}$.