# Fileset

[s42005-025-02380-y.pdf](https://mdr.nims.go.jp/filesets/12d30c4c-c74c-4d1c-bc66-843b1f683e67/download)

## Creator

Chen Liang, Diptesh Das, Jiang Guo, [Ryo Tamura](https://orcid.org/0000-0002-0349-358X), [Zetian Mao](https://orcid.org/0000-0002-4853-7574), [Koji Tsuda](https://orcid.org/0000-0002-4288-1606)

## Rights

[Creative Commons BY-NC-ND Attribution-NonCommercial-NoDerivs 4.0 International](https://creativecommons.org/licenses/by-nc-nd/4.0/)

## Other metadata

[Predicting symmetric structures of large crystals with GPU-based Ising machines](https://mdr.nims.go.jp/datasets/b7bab30b-3174-409d-9c1b-bcea2704e509)

## Fulltext

Predicting symmetric structures of large crystals with GPU-based Ising machinescommunications physics ArticleA Nature Portfolio journalhttps://doi.org/10.1038/s42005-025-02380-yPredicting symmetric structures of largecrystals with GPU-based Ising machinesCheck for updatesChen Liang1, Diptesh Das1, Jiang Guo1, Ryo Tamura 1,2, Zetian Mao 1 & Koji Tsuda 1,2,3Solving black-box optimization problems with Ising machines is increasingly common in materialsscience. However, their application to crystal structure prediction (CSP) is still ineffective due tosymmetry agnostic encoding of atomic coordinates. We introduce CRYSIM, an algorithm thatencodes the spacegroup, theWyckoff positions combination, andcoordinatesof independent atomicsites as separate variables. This encoding reduces the search space substantially by exploiting thesymmetry in space groups. When CRYSIM is interfaced to Fixstars Amplify, a GPU-based Isingmachine, its prediction performance is competitive with CALYPSO and Bayesian optimization forcrystals containingmore than 150 atoms in a unit cell. Although it is not realistic to interfaceCRYSIM tocurrent small-scale quantumdevices, it has the potential to become the standardCSP algorithm in thecoming quantum age.The advancement of various significant technology fields relies ondiscovery of innovative materials with desired chemical or physicalproperties1, and obtaining correct structures of materials, arrange-ment of atoms in the unit cell, is the prerequisite. To achieve the goal,crystal structure prediction (CSP)2,3, in which the most stable crystalstructure is inferred only from its chemical composition, has beenwidely adopted. The vast configuration space and the richness oflocal minima on potential energy surfaces (PESs) renders CSP achallenging task4. Optimization algorithms, such as geneticalgorithms5–9, particle-swarm optimization (PSO)10,11, Bayesian opti-mization (BO)12–14, are proposed and successfully applied in practice.Typically, they create roughly-shaped initial structures and the finaloptimization is done by a geometric relaxation software either basedon first-principles calculation or pretrained universal neural networkpotentials (NNPs). Nevertheless, these methods generally require agreat number of iterations. In recent years, deep learning-basedcrystal generative models15–21 are developing fast, but they might findproblems in extrapolation outside their training datasets. Therefore,as an example, due to the scarcity of corresponding data, CSP on 2Dmaterials22,23 and nanoclusters24–26 generally relies on optimizationmethods. Besides, in both categories, most of the methods work wellfor crystals containing less than 60 atoms18 in a unit cell, but are stillnot ideal for larger crystals. For example, the training data ofCDVAE15 and MatterGen21 does not include large crystals with morethan 20 atoms in a unit cell. GNoME16 successfully explores largerones approximating 100 atoms, but considerable computational costis required.Isingmachines27,28 are hardware-assisted discrete optimizers that solvea quadratic unconstrained binary optimization (QUBO) problem,x* ¼ argminx2f0;1gMXMi¼1hixi þXMi¼1XMj¼iþ1Jijxixj !; ð1Þwhere x is an M-dimensional bit vector and hi and Jij are real-valuedparameters. CSP can be represented as a QUBO problem, either bysimplifying the energy function29–31 or the use of a surrogate machinelearning model32,33. Among the existing studies, Gusev et al.29, Ichikawaet al.30, and Xu et al.32 provided QUBO formulations of CSP under a fixedspace group. Couzinié et al.31 and Couzinié et al.33 disregarded the spacegroup and employed a grid-based representation of atomic coordinates intheir QUBO formulation. Notably, they lack a feature of dynamicallyadjusting the space group, which is common in state-of-the-art CSPalgorithms such as CALYPSO10,11, USPEX7,9,34, and CRYSPY13,14. Further-more, these algorithms have not been tested on complex problems, mainlydue to the scale restriction of D-Wave35 quantum annealer.In thiswork, we develop amethodnamedCRYSIM (CRYstal structureprediction with Symmetry-encoded Ising Machine). Our bit vector repre-sents the lattice parameters, the symmetry information including the crystalsystem (CS), the space group (SG) and theWyckoff positions combination(WPC), and coordinates of independent sites. This bit vector is translated toa crystal structure, and M3GNet36 provides its potential energy. Our goal isto find the optimal bit vector that gives the lowest potential energy. Toenable the search with an Ising machine, a 2-order factorization machine1Department of Computational Biology and Medical Sciences, Graduate School of Frontier Sciences, The University of Tokyo, Kashiwa, Japan. 2Center for BasicResearch on Materials, National Institute for Materials Science, Ibaraki, Japan. 3RIKEN Center for Advanced Intelligence Project, Tokyo, Japan.e-mail: zt.mao97@gmail.com; tsuda@k.u-tokyo.ac.jpCommunications Physics |           (2025) 8:477 11234567890():,;1234567890():,;http://crossmark.crossref.org/dialog/?doi=10.1038/s42005-025-02380-y&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s42005-025-02380-y&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s42005-025-02380-y&domain=pdfhttp://orcid.org/0000-0002-0349-358Xhttp://orcid.org/0000-0002-0349-358Xhttp://orcid.org/0000-0002-0349-358Xhttp://orcid.org/0000-0002-0349-358Xhttp://orcid.org/0000-0002-0349-358Xhttp://orcid.org/0000-0002-4853-7574http://orcid.org/0000-0002-4853-7574http://orcid.org/0000-0002-4853-7574http://orcid.org/0000-0002-4853-7574http://orcid.org/0000-0002-4853-7574http://orcid.org/0000-0002-4288-1606http://orcid.org/0000-0002-4288-1606http://orcid.org/0000-0002-4288-1606http://orcid.org/0000-0002-4288-1606http://orcid.org/0000-0002-4288-1606mailto:zt.mao97@gmail.commailto:tsuda@k.u-tokyo.ac.jpwww.nature.com/commsphys(FM)37–40 is trained with available pairs of bit vectors and correspondingenergies with an active learning workflow41. Since the prediction function isquadratic, the optimal bit vector that minimizes the FM-approximatedpotential energy can be found with an Ising machine. It does not alwayscoincide with the real optimal solution, but one can expect that the errordecreases as the amount of training data increases during the searchprocess.Our information-rich bit vector inevitably inhibits the use of D-Wavequantum annealers. Instead, we employ a GPU-based Ising machine, Fix-starsAmplify42, to solve aproblemwithover several thousandbits. It is basedon simulated annealing and usesmulti-level parallel processing onmultipleGPUs to find optimal solutions. Fixstars Amplify relies on conventionalsemiconductor technologies, but can handle large-scale problems up to130,000 bits with full connectivity. It has been employed in moleculargeneration43, materials design44–46, and various engineering fields47,48.Ourmethod outperforms BO12 andCALYPSO10,11 on three small crystalsas well as large ones containingmore than 150 atoms in unit cells. Specifically,CRYSIM successfully generates the ground truth Ca24Al16(SiO4)24 structure,containing 160 atoms in the unit cell, within 300 relaxations in 4 out of thetotal 5 trials. In this work, GPUs are adopted, but CRYSIM can leverage anyIsing machines, including rapidly developing quantum devices.ResultsBit vector encodingThe binary representation in CRYSIM consists of the following three parts:lattice parameters, symmetry information, and 3D coordinates of inde-pendent sites. In the first part, the six-dimensional lattice parameters areindividually discretized and summarized into a bit vector with one-hotencoding. The second part includes a crystal system (CS), a space group(SG), a group ofWyckoff position combinations (WPCs). The sizes of eachvector segment depend on the set of all possible space groups compatiblewith the given chemical composition, which is determined bywhether thereexists at least one WPC for achieving symmetry of the SG. Similarly, onlycompatible CSs are included in the embeddings. Accordingly, if m crystalsystems are involved, each of which has s1,…, sm compatible space groups,the CS part hasm bits to represent the CS and the SG part has maxi¼1;...;msibits to represent the SG. If the crystal structure has the i-th CS and j-th SG,the corresponding bits are set as 1, and the remaining are 0s.Given an SG, we engineered theWPC generator in GN-OA package49to compute the list of plausible WPCs according to the input chemicalcomposition50,51. In the process, WPCs for each element species are calcu-lated recursively based on its frequency in the unit cell. They are sorted in adescending order according to the maximummultiplicity of involvedWPs,so thatmoreplausible combinations areprioritized, considering the fact thatmost crystals in nature tend to occupy WPCs with general Wyckoffpositions34,52. Then, the Cartesian product among the WPCs of elements isconducted to deriveWPCs for the whole system, which are appended in thelist until themaximumnumber, 30,000, is reached. Finally, they are dividedinto 300 groups of size 100, which is encoded in the WPC segment with a300-dimensional one-hot vector for specifying the group.The third part consists of k copies of a g3-dimensional bit vector inorder to represent a crystal containing k element species, in which g denoteslattice discretization resolution (LDR). A 3D g× g× g grid is assumedwithinthe unit cell. If an independent site of the atom species exists near a gridpoint, the corresponding bit is set to 1. In decoding, 100 structures aregenerated corresponding to all WPCs in the specifiedWPC group. Amongthem, the one with the largest minimum interatomic distance (MID) isselected to increase the possibility of deriving stable states4. Details of theencoding and decoding procedures are provided in the Method. Besides,Supplementary Note 1 and 2 present a detailed explanation about WPCsgeneration and application.CRYSIMWorkflowThe workflow of CRYSIM is depicted in Fig. 1. First, 1000 initial structuresare obtained by random generation (RG) developed in this work (seeMethod for details) with the given chemical composition, and converted tobit vectors xl. Their potential energies yl are estimated using M3GNetwithout structure relaxation, leading to a training dataset described as thepairs of bit vectors and energies, i.e., D = {(xl, yl)∣l = 1, 2, …, 1000}. Thedataset is then used to train an FM model 37,38, simulating the M3GNet-estimated potential energy surface of the materials family. The functionalform of FM is described asy ¼ bþXMi¼1hixi þXMi¼1XMj¼iþ1XKk¼1wkiwkjxixj; ð2Þwhere b, hi andwki are real-valued parameters, and xi is the i-th element of avector x. It is similar toQUBO, but the weightmatrix of quadratic terms is alow-rank matrix parameterized by wki and wkj. Then, Fixstars Amplify42 isused tooptimize the bit vector tominimize the energy approximatedbyFM.PyTorch53 is leveraged to implement FM and the training workflow, withdetailed training settings and hyperparameters of Amplify provided inSupplementary Note 3.Based on the solution, 100 structures with the same SG but differentWPCs from the solvedWPCsgroup are translated back to crystal structures,and the one with the largest MID is selected and relaxed using M3GNet.After relaxation, we sample 30 structure frames from the relaxation tra-jectory. Among the samples, if a structure has anMID lower than 0.5Å butis still assigned a negative energy, the energy is adjusted to a high positivevalue to mitigate negative impact due to inaccuracy of NNP. Then the datapoints are added toD and FM is retrained. The above procedure is repeatedT = 300 times and the most stable structure is recorded as the final result.Recovering benchmark materialsWe begin with relatively simple benchmark crystal tasks to demonstrateCRYSIM’s ability to address CSP. Wei et al.54 developed a series ofquantity measurements for evaluating CSP algorithms, and selected fivecrystals, including ScBe5, Ca4S4, Ba3Na3Bi3, Li4Zr4O8 and Li3Ti3Se6O3, asexamples to conduct tests. Ground states of all compounds are deter-mined based on the Materials Project (MP) database55, i.e., mp-11277,mp-1672, mp-31235, mp-755253 and mp-1211008, respectively. Classi-cal algorithms considered for comparison include CALYPSO10,11 andsimple BO that directly optimizes lattice parameters, fractional coordi-nates, the SG number and theWPC index of crystals, implemented basedon the hyperopt package56 following the scheme of GN-OA49, denoted as“Crystal param. + Hyperopt BO”. All methods are limited to perform300 times of structure relaxation during one run to make a fair com-parison. Accordingly, CALYPSO is leveraged for 30 generations, with thepopulation size per iteration set as 10. Since the number of atoms are notoptimized in CRYSIM, we treat the structure indicated by the inputchemical formula as conventional cells, i.e., symmetry is calculatedwithin the given lattice, instead of a combination of several cells.Accordingly, the NumberOfFormula parameter of CALYPSO,representing the smallest and largest times of the input chemical formulaof considered structures in one optimization process, is set as “1 1” in alltests. Besides, in all experiments in this article, structures containinginteratomic distances smaller than 1.0 Å are excluded from the statistics,to ensure that all crystals remain physically valid. Tests of each methodrepeat three times with different seeds. Hyperparameters of classicalalgorithms are reported in Supplementary Note 4.We use StructureMatcher function from the pymatgenpackage57 to determine structural similarity between predicted and knownground truth materials. The function can compute minimum average pair-wise displacement between two corresponding atoms in two configurationsamong all permutations. The predicted structure successfully matches theground truth as long as the displacement is computable, which suggests thatStructureMatcher is able todistinguish corresponding atomsbetweenthem. Details of criterion of matching is provided in Method section.When assessing results, we only include structures reaching the lowestenergy (Emin) among all obtained ones in 300 iterations, i.e., fEjE ¼ Eming,https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 2www.nature.com/commsphysand compare them against the ground truth. Several major metrics aredefined for this task: (1) IM,0 denotes the first iteration at which the groundtruth is identified; (2) NE denotes the number of iterations reaching thelowest energy, i.e., NE ¼ jfEjE ¼ Emingj; (3) NM denotes the number ofsuccessfully matched ones among all considered structures. We provide afurther illustration on evaluation of CSP algorithms inDiscussion section.Figure 2 summarizes results of two representative crystals, and compre-hensive information is presented in Supplementary Tables 1–5. Othermetrics, such as displacement calculated by StructureMatcher of thestructure in iteration IM,0, denoted as DM,0, as well as the minimumdisplacement DM;min and corresponding iteration IM;min, are also reported.Predicted configurations with the lowest estimated relaxed energies in thethree trials are shown in Fig. 3.Ground states of ScBe5 andCa4S4 can be readily discovered by all threemethods, but CRYSIM generates significantly more stable states than thetwo classical methods. Besides, smaller IM,0 of CRYSIM optimizers indicatethat FM can quickly and effectively characterize the PES by learning frominitial datasets. For themore complicated Ba3Na3Bi3 system,CRYSIM is theonly method that successfully discovers the stable state with correct esti-mated energies.Fig. 1 | The workflow of CRYSIM that contains T iterations, using Si4O8 as anillustration. Thin arrows denote the workflow at the t-th iteration, and thick arrowsdenote entering and exiting iterations. a Given the considered material system, adataset is obtained by random generation (RG) to provide training samples anddetermine the upper bound of lattice parameters for binary representation. Potentialenergy of each material is also estimated by pretrained neural network potential(NNP) without structure relaxation. b Structures in the dataset {S1, S2, …, S1000} areencoded into binary vectors {x1, x2, …, x1000} using symmetry-informed integerencoding. c Factorization Machine (FM) is used to perform regression from thebinary vectors to their corresponding estimated energies, obtaining the objectivefunction to be optimized. d An Ising solver is employed to solve the learned objectivefunction to minimize y in t-th iteration, resulting in x*,t. Amplify is used in this work.e The solved binary embeddings x*,t is decoded into crystal structures. Since one bitin the Wyckoff position combination (WPC) segment represents a group of 100WPCs, 100 structures are derived. The one with the largest minimum interatomicdistance (MID) is selected as S*,t. We note that the Si4O8 structures drawn in the figure(e) are indicative, which have different space groups (SGs). f The solved structure S*,tis relaxed by NNP, leading to a structure-energy pair ðS�;tr ; Etmin;rÞ. If iterations havenot finished, frames in the relaxation trajectory are sampled. g Among the sampledstructures, if one contains an MID smaller than 0.5Å but still is estimated to have anegative energy, the energy is reassigned with a high positive one before adding thepoints into the training dataset for the next iteration. After finishing all iterations, thefinal structure S�r , the one with the lowest relaxed energy among all crystals in all Titerations, will be regarded as the discovered stable structure of this system.https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 3www.nature.com/commsphysAll methods fail on Li4Zr4O8 and Li3Ti3Se6O3 structure prediction ifonly the crystals of Emin are counted, which, however, reach an even lowerestimated relaxed energy than the ground states. This may be attributed totwo main reasons. First, the selected benchmark structures may notrepresent the ground states of corresponding chemical compositions, sug-gested by positive energies above hull. Especially, Li3Ti3Se6O3 (mp-1211008) exhibits a substantial energy above hull, 0.617 eV/atom accordingto MP, indicating a potential to transform into alternative phases predictedby the optimizationmethods. Second, the pretrained NNP leveraged in thisstudy is not sufficiently accurate, and a 0.01 eV variation in potential energycan affect the behavior of optimization algorithms. As an instance, wefurther introduce predicting results on Li8Zr4O12 (mp-4156), a stablestructure of Li-Zr-O family that has been observed in experiments, inSupplementary Fig. 1 and Supplementary Table 6. In all configurationswiththe lowest energies, hexa-atomic rings absent in mp-4156 can be found,which may suggest a systematic energy underestimation on near-equilibrium structures, as is discussed in other related works58. On theother hand, CRYSIM is a flexible tool, which can be incorporated with anyenergy estimators. For small systems, first-principles calculationsoftware59,60 is available.We also expect that the continuous improvement ofNNP technology can mitigate the problem.Large crystal structures predictionThis section introduces experiment results on three large material systems,including Y6Co51, Ca24Al16(SiO4)24 and (SiO2)96, to demonstrate capabilityofCRYSIM.The crystal Y2Co17 has been a classical benchmark for assessingCSP algorithms12,61,62, but the two stable structures in this material familyrecorded in MP can only be achieved with unit cells Y4Co34 (mp-570718)and Y6Co51 (mp-1106140). Since the number of atoms in unit cells are notoptimized inCRYSIM,we start directlywithY6Co51.Ca24Al16(SiO4)24 (mp-6008)29,63 and (SiO2)3263 have also been discussed in previous works asexamples of CSP on large crystals. Here, (SiO2)96 (mp-1200292) is chosensince it is the largest SiO2 crystal in MP that has been observed inexperiments.Apart fromCALYPSO andBO introduced earlier, simple RG, which isemployed to generate initial training set for CRYSIM, and PyXtal52-basedRG in CRYSPY13, denoted as “CRYSPY RG”, are additionally included asbaseline CSP methods. 300 times of structure relaxation, i.e., 300 iterations,are conducted in one run, and tests of each method repeat five times withdifferent seeds. Fig. 4a-c presents averaged accumulated lowest energies ofcrystals during generation in the 300 cycles. Optimal materials found byeach algorithm are visualized in Fig. 4d. Corresponding data is summarizedin SupplementaryTable 7, aswell as SupplementaryTables 8–10 formetricsdefined in the last section.OnY6Co51, BO is the onlymethod that discovers the ground state witha−7.127 eV/atomenergy, the same as corresponding relaxed energy ofmp-1106140. However, computational complexity of BO scales with the totalnumber of atoms64–66, leading to a significantly reduced performance onCa24Al16(SiO4)24 and (SiO2)96, even falling below CRYSPY RG. CALYPSOexhibits a higher stability thanBO for large systems, and the implementationof pair-wise distance consideration in input renders it unaffected whenscreening out configurations with small MIDs, as is shown in Supplemen-tary Table 11. However, this feature accelerates the optimization of energiesonly in the first tens of iterations in Fig. 4b-c, and then the algorithm issurpassed by CRYSIM methods.On the other hand, length of CRYSIM embeddings is determinedsolely by the number of elements given the LDR, making it advantageousover other algorithms especially on large crystals. On the Ca24Al16(SiO4)24system, which contains 160 atoms in the unit cell, CRYSIM successfullyfinds the ground truth structure in four out of five trials. For (SiO2)96,CRYSIM identifies a configurationwith a relaxed energy (−7.890 eV/atom)close to the stable one (−7.891 eV/atom), significantly lower than the onesfound by other methods. CSP for large crystals has been a long-standingchallenging task. Energy distribution of the configuration space tends toconcentrate on unstable states as the system size grows, which means thatthedifficulty offinding thegroundstate viaRGexponentially increases63,67,68.Though CRYSIM does not outperform in all systems, the superiority onlarge crystals establishes it as a promising approach for CSP.Fig. 2 | Optimization performance comparison between CRYSIM and two clas-sical crystal structure prediction (CSP) algorithms, CALYPSO and Bayesianoptimization (BO), on two benchmark materials. The first iteration when thegenerated structure matches the ground truth (IM,0), and the number of successfullymatched structures among the generated ones with the lowest energy (NM) arerecorded for (a, b) Ca4S4 and (c, d) Ba3Na3Bi3 for the three optimization methods in300 iterations with 3 seeds. Shadowed bars in c indicate that the correspondingmethods fail to find the ground truth structure with these seeds.https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 4www.nature.com/commsphysEffects of Processing TechniquesLattice Discretization Resolution. When converting 3D structures intobinary vectors, a higher LDR can reduce information loss, nevertheless,leading to exponentially increasing solving difficulty. We investigate theinfluence of LDR on CRYSIM optimization performance by testingg∈ {5, 7, 9, 12, 15} across the three large crystals considered in this study,namely Y6Co51, Ca24Al16(SiO4)24 and (SiO2)96. Table 1 summarizeslowest energies with corresponding average accumulated energy curvesrecorded in Supplementary Fig. 2, where CRYSIM with different LDRsare denoted as CRYSIM-g, such as “CRYSIM-5” for g = 5. Each value isthe mean of three trials with different random seeds.Y6Co51 and Ca24Al16(SiO4)24 achieve the best CSP results at LDR 12,while (SiO2)96 performs the best at 15. This phenomenon stems from theimplementation of CRYSIM, where each atom species is encoded using ag × g × g grid of bits, with 0 or 1 indicating the presence of an atom at eachdiscretized unit. Accordingly, a higher LDR is advantageous when thenumber of atoms per element increases, rather than the total number ofatoms. An LDR of 12 appears to strike a balance between representabilityand optimization difficulty for Y6Co51 and Ca24Al16(SiO4)24, where themaximum number of atoms of a specific element is 51 and 96, respectively.In contrast, (SiO2)96, which contains 192 oxygen atoms, may require ahigher LDR to better capture interatomic spatial relationship in CRYSIM.According to the results, as a guidance on selecting LDRs for other tests, it isrecommended touse a 15LDRwhen themaximumfrequencyof an elementspecies in theunit cell apparently exceeds 100, otherwise 12.Numbers of bitsfor representing each parameter for the three systems are further reported inSupplementary Table 12. We also provide time costs in the Ising solvingprocess byAmplify and spatial resolutions of each block for different LDR inSupplementary Table 13 and Supplementary Table 14, respectively. Werefer the readers to the Method section for details of calculation.Factorizationmachine. In this work, FM is employed as the regressor inCRYSIM to build Ising objective functions, and here the fitting accuracyis investigated. Supplementary Fig. 3a–e shows distributions of predictedversus calculated energies of one of the initial Ca24Al16(SiO4)24 datasetsderived by RG, the system requiring the largest number of bits torepresent due to its chemical composition. Learnable parameters of FMare decided upon metrics on the validation set, comprised of 10% of thedataset (see Supplementary Note 3 for details). Changes of Pearsoncorrelation coefficients (PCCs) andmean absolute errors (MAEs) duringtraining are further provided in Supplementary Fig. 3f. The consistentlyhigh PCC values indicate effective optimization toward the global opti-mum, despite of fluctuations due to out-of-distribution energies. Similartrends are observed across all other systems and random seeds.Additionally, a comparison between FM and full-rank quadraticregression (QR), in which quadratic terms in regression functions areFig. 3 | Comparison between the most stable benchmark configurations dis-covered by three crystal structure prediction (CSP) algorithms, CALYPSO,Bayesian optimization (BO) and CRYSIM, and the ground truth structurerecorded in Materials Project (MP). The side views (left column for each method)and top views (right column) of materials are visualized by VESTA software84,including mp-11277 (Sc: purple, Be: green), mp-1672 (Ca: blue, S: yellow), mp-31235 (Ba: green, Na: yellow, Bi: pink), mp-755253 (Li: green, Zr: blue, O: red), andmp-1211008 (Li: green, Ti: blue, Se: orange, O: red), respectively, with M3GNet-estimated relaxed energies labeled above. We note that the optimization algorithmsare guided by the total energies of the systems, which are normalized afterwards inthis figure by system sizes. Most configurations are expanded into superlattices todisplay the patterns. Crystals with the lowest energies are selected. If there are morethan one crystals having the same energy, the one obtained in the earliest iterationis shown.https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 5www.nature.com/commsphysindependently learned instead of multiplications between linear terms, isexhibited in Supplementary Fig. 4. These experiments are conducted withCRYSIM-5 representation of Y6Co51 system, containing 801 bits in theembeddings. Under these conditions, QR involves more than 600,000trainable parameters, whereas FM requires only 13,617 ones. For reference,optimization results of BO and CALYPSO, previously shown in Supple-mentary Table 7 are also included as baselines, with three trials performedfor each method. On this system, CRYSIM-QR achieves a lower averageaccumulated energy than CRYSIM-FM, indicating superior optimizationperformance. However, for systems represented with more than 1,000 bits,QR requires millions of trainable parameters, making FM a more practicaloption for these tasks in terms of computational efficiency.Processing with minimum interatomic distance. Inaccuracy of NNPson configurations with extremely smallMIDs renders negative impact onregression models. Tomitigate the effect, procedures related toMIDs aredesigned and integrated in theworkflow, including selecting the structurewith the largest MID from one solution vector in Fig. 1f, and adjustingunphysical energy estimations in Fig. 1g. Effectiveness of the proceduresis demonstrated inTable 1 by comparingCRYSIMof all considered LDRswith andwithout including these steps during optimization. Tests of eachmethod repeat three times with different seeds. The correspondingaccumulated energy curves are provided in Supplementary Fig. 2. Opti-mizers equipped with the modules achieve a widespread performanceenhancement, particularly for larger systems with higher LDRs. Besides,inclusion of MID processing enables structures exploration with highLDRs. For Y6Co51 and Ca24Al16(SiO4)24, CRYSIM-12 performs the bestwith MID-related procedures, but it cannot realize the full potential andis outperformed byCRYSIM-9without them. The advantage is attributedto improved efficiency in obtaining valid crystals, indicated by a reduc-tion of filtered-out configurations reported in Supplementary Table 15.DiscussionIn this section, we provide a detailed interpretation covering: (i) symmetryrepresentation, (ii) atomic position representation, (iii) MID-related pro-cedures in theworkflow, (iv) evaluation on end-to-endCSP algorithms, and(v) the use of Ising solvers.First of all, incorporating symmetry representations is a fundamentalrequirement in developing algorithms for crystal generation. Most stablecrystal structures in the nature have symmetry69,70, making it necessary toFig. 4 | Optimization performance comparison between CRYSIM and classicalcrystal structure prediction (CSP) algorithms on large crystals. Normalizedaccumulated lowest M3GNet-estimated relaxed energies in the course of iterationsof Y6Co51, Ca24Al16(SiO4)24 and (SiO2)96 structures derived from various CSPalgorithms are shown in (a–c), respectively. Each curve is averaged on five tests withdifferent random seeds, and colored shaded areas cover the maximum and mini-mum in the five trials. Dash lines are relaxed energies of ground truth materials inMaterials Project (MP). Representative predicted configurations after relaxation, aswell as corresponding ground states in MP, are displayed in (d), including side view(left column for each method) and top view (right column) of mp-1106140 (Y: grey,Co: blue), mp-6008 (Ca: grey, Al: light blue, Si: deep blue, O: red), and mp-1200292(Si: blue, O: red), respectively, visualized by VESTA software84, with M3GNet-estimated relaxed energies labeled above. We note that the optimization algorithmsare guided by the total energies of the systems, which are normalized afterwards inthis figure by system sizes.https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 6www.nature.com/commsphysconsider symmetry constraints in the CSP task. Classical optimizationmethods, including evolutionary algorithms34,71, PSO10,11, and BO12–14,49, relyon RG algorithms50–52 to sample symmetric configurations as startingpoints. Following discussions proposed in RandSpg51, two strategies oforganizing and utilizing WPCs lists are summarized in SupplementaryNote 1. CRYSIM differs from the above strategies since the search is con-ducted within a space of symmetric structures, in which the SG number isoptimized directly. Other methods, including USPEX34 and CALYPSO10,11,only guarantee the symmetry requirement in initialization.Recently, symmetry of crystals is also emphasized in generativemodelsto increase the possibility of deriving reasonable configurations. As one ofthe earliest trials in this field, DiffCSP++17 does not represent symmetry ofcrystals explicitly, in which networks generate configurations given a spe-cific pre-set WPC constraints. WyCryst20, SymmCD72 and WyckoffDiff 73successfully explore symmetry representations with one-hot labelling,operation matrices and graphs, respectively. However, current generativemodel-based solutions are still not perfect from the perspective of end-to-end CSP. The related works, including the ones trained with an auto-regression strategy, such as CrystalFormer74 andWyFormer75, and diffusionmodels, such as SymmCD and WyckoffDiff, rely on a prior distribution ofSG from the training dataset for conditional generation.WyCryst is amongthemethods that predict both SG andWPC for constructing structures, butit requires PyXtal as an external tool to provide atomic positions and fill inWPCs, similar to WyFormer and WyckoffDiff. In comparison, CRYSIMtreats the CSP problem in an exact end-to-endmanner, in which SG,WPCindex and atomic positions are all optimized and do not condition on anygiven information, except for the input chemical formula.Second, an accurate description on the CSP problem relies on aneffective representation of atomic coordinates in the lattices. There are twomain strategies for representing fractional coordinates as variables to beoptimized. The first is to treat each coordinate (xi, yi and zi for the i-th atom)as independent variables11,49, and the second is to split the whole crystallattice and use the derived discrete blocks in the 3D space to encodepositions29–32. Most optimization methods based on Ising models adopt thesecond approach, as it aligns with the goal of achieving guaranteed optimalsolutions through quantum annealing by fitting the system’s PES. Byrepresenting atomic positions via lattice splitting, an Isingmodel can encodephysical interactions: first-order terms capture the energy contribution of asingle atom due to external fields, while second-order terms describe pair-wise atomic interactions. This allows the Ising model to approximate theinteratomic potential in a quadratic form.Nevertheless, in practical implementations that account for symmetry,such as CRYSIM and other works29, the solved atomic positions are notexternal coordinates for building structures, but internal or independentsites to insert in WPs. As a result, the learned Ising model does not fullyreflect an actual interatomic potential. One possible solution is to firstestimate or sample an SG and WPC, derive corresponding constraints onlattice parameters and coordinates, and then optimize the two parts.Accordingly, by adding penalty terms, Ising solver can optimize directly onexternal coordinates and preserving the symmetry simultaneously. How-ever, the order of constraints on coordinates would be the same as themultiplicity of correspondingWPs, making it challenging to implement forcurrent solvers.From a practical perspective, another advantage of the second strategyover the first one is that it requires less bits to encode coordinates in largesystems, especially those with few atomic species. This is because the secondstrategy scales linearly with the number of atomic species and remainsconstant with respect to the number of atoms, whereas the first strategyscales linearlywith the number of atomswithin the cell. The scalingwith thenumber of atoms is generally more computationally intensive in large sys-tems. Taking the (SiO2)96 system in Results as an example, suppose thelattice is split into 15 * 15 * 15 blocks. Following the first strategy, eachcoordinates require 15 bits to represent, leading to 15 * 3* (96 + 192) = 12960 bits in total, but only 15 * 15 * 15 * 2 = 6750 areneeded based on the second one.Third, the MID-related procedures employed in this study arenecessary for preserving effectiveness of the algorithm and enhancingoptimization efficiency practically. Ideally, a pretrained NNP shouldassign high energies to unstable structures, allowing the objectivefunctions to reflect an accurate structure-energy relationship throughactive learning for diversion structural configurations. Accordingly,CSP optimizers, designed to identify low-energy solutions, can cor-rectly discover the ground states. However, training sets for state-of-the-art NNPs generally lack out-of-distribution structures in theconfiguration space, especially for crystals containing extremely closeatom pairs, rendering their estimated energies unreasonably low. Aspresented in Supplementary Table 15, many CRYSIM optimizerswithout MID processing are encouraged to generate abnormalstructures, due to their low estimated relaxed energies, which are,however, ineffective from a practical perspective. Directly replacingpretrained NNP with first-principles calculation software59,60 cancircumvent the problem, but it is still unrealistic for large crystalsconsidering current computational power.According to an observation that most abnormal relaxed structuresoriginate from abnormal decoded unrelaxed ones, we design MID proces-sing techniques on generated configurations (Fig. 1f, g) to improve theefficiency of CRYSIMoptimizers. Besides, previous works4 also suggest thatatoms in stable structures tend to uniformly distribute, instead of clusteringin a small space. However, for some material systems, strategies aimed atcontrollingMIDs of generatedmaterials (e.g., CRYSPYRGandCALYPSO)or attempting to obtain materials with larger MIDs (e.g., CRYSIM) mayhinder discovering of ground states, as evidenced by the Y6Co51 system inSupplementary Table 7. Although many structures derived by RG and BOimplemented in this work are screened out due to very small MIDs, asreported in Supplementary Table 11, CSP of Y6Co51 is finally accomplishedafter structure relaxation on crystals that may not have high MIDs. Weexpect that when a more accurate pretrained NNP is proposed, the MID-related procedures can be discarded.Table 1 | Influence of lattice discretization resolutions (LDRs) and minimum interatomic distance-related (MID-related)procedures on optimization performanceSystem MID proc. Lattice Discretization Resolution5 * 5 * 5 7 * 7 * 7 9 * 9 * 9 12 * 12 * 12 15 * 15 * 15Y6Co51 N −6.822 ± 0.067 −6.852 ± 0.054 −6.878 ± 0.038 −6.862 ± 0.088 −6.796 ± 0.010Y 6:878 ± 0:011 �6:888± 0:070 �6:910± 0:032 �6:913 ± 0:038 �6:855 ± 0:013Ca24Al16(SiO4)24 N −6.670 ± 0.038 −6.542 ± 0.074 −6.858 ± 0.077 −6.766 ± 0.082 −6.654 ± 0.259Y �6:984 ± 0:068 �7:266± 0:103 �7:139± 0:289 �7:456 ± 0:029 �7:070 ± 0:244(SiO2)96 N / �6:888± 0:109 −6.559 ± 0.111 −6.507 ± 0.189 −7.120 ± 0.070Y / −6.732 ± 0.211 �6:992± 0:281 �6:737 ± 0:140 �7:366 ± 0:380Lowest energiesof structuresdiscoveredbyCRYSIMoptimizers of different LDRsare shown in the figure, aswell as results of integratingMID-relatedprocedures (Y) or absenceof it (N) as an ablation study.Underlined values are the lowest average energies for each material system achieved by each LDR, and bold values are the lowest ones of each MID processing strategy among all LDRs. Each value isaveraged on three seeds. (unit: eV/atom).https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 7www.nature.com/commsphysForth, this work adopts a simple but rigorous criterion to evaluate end-to-end CSP optimizers. The primary objective of our optimization-basedCSP algorithm, CRYSIM, is precisely to locate the global minimum on thePES, representing themost thermodynamically stable structure according tothe chosen energymodel. To rigorously assess this specific capability duringbenchmarking,we consider the number of “successfulmatch” (NM) only forstructures with the lowest relaxed energies out of the 300 total structuresinstead of all of them. This prevents crediting success to fortuitous samplinginto higher-energy local minima, thereby isolating the evaluation of opti-mization performance.More critically, this methodology also reflects the practicability in CSPtasks. When the target structure is unknown, researchers inevitably rely onthe calculated energy ranking, treating the lowest-energy prediction as themost likely candidate for experimental synthesis or validation. Higher-energy predictions, even if potentially correct, cannot be identified as suchwithout prior knowledge of the ground truth. Thus, by focusing our eva-luation on the lowest-energy structure, our metric is not only relevant toCRYSIM’s optimization objective but also aligned with the practical inter-pretation and utility of CSP results.Finally, as a solver for addressing combinatorial optimization pro-blems, quantum annealing (QA)76–81 has gained attention due to its theo-retical ability to escape from local minima, and D-Wave system35 is amongthemostwidely used implementations ofQA29,31,33.However, themaximumnumber of variables for the D-Wave system is limited to 124 bits43, severelyrestricting its application. In thiswork,we integrateAmplify42, aGPU-basedIsing solver, into CRYSIM as a substitute of quantum annealers. Never-theless, we claim that the present implementation can be applied directly onquantum annealers without adjustment as quantum computers continue todevelop.ConclusionIn conclusion,CRYSIM, aCSPoptimizerbasedon symmetry-encoded Isingmodels, is proposed and tested across various CSP tasks. To the best of ourknowledge, it is the first Ising machine-based optimizer for CSP thatdynamically optimizes symmetry. CRYSIM outperforms CRYSPYRG, BO,and CALYPSO on most systems, showcasing its strong optimization cap-abilities not only for small benchmark crystals but also for larger ones,including Ca24Al16(SiO4)24 and (SiO2)96. The predicting accuracy of FM inCRYSIM is also discussed, highlighting its expressivity in CSP tasks.CRYSIMoffers a promising Isingmachine-based optimization tool for CSPthat could potentially be applied to quantum annealers in the future.MethodsRandom generation of crystal structuresRandom generation (RG) constitutes the foundation of CSP11,13,68, in whichatomic positions are randomly sampled according to the given chemicalcomposition. In this work, a simple crystal RG tool is implemented as theCSP baseline, as well as preparing training data for Ising models.In RG, six lattice parameters (lattice lengths a, b, c, and lattice angles α,β, γ) and fractional coordinates are sampled independently from uniformdistributions. To determine the default lower and upper bounds of thedistributions for lattice lengths, we perform statistical analysis on materialsinMP. LetMbe the set of allmaterials in theMPdatabase.Define a functionjmj : M ! N thatmaps eachmaterialm to the number of atoms in its unitcell. We then partitionM into five categories, {M1,M2,M3,M4,M5}, suchthat for each materialm ∈M, it belongs to categoryMi if it satisfiesM1 ¼ fmjjmj≤ 20g;M2 ¼ fmj20 < jmj≤ 50g;M3 ¼ fmj50 < jmj≤ 80g;M4 ¼ fmj80 < jmj≤ 100g;M5 ¼ fmjjmj > 100g:8>>>>>><>>>>>>:ð3ÞNext, the averages of a, b, and c for materials in each category Mi arecomputed, denoted as aMi, bMi, and cMi, respectively. For a specificmaterialsystemm0 to be generated, if thenumberof atoms ∣m0∣ fallswithin oneof theranges, the lower and upper bounds are determined as follows:lm0¼ 0:8 × ðaMiþ bMiþ cMiÞ;um0¼ 2 × ðaMiþ bMiþ cMiÞ:(ð4ÞThese bounds are the same for the three lattice lengths. The lower andupperbounds for the lattice angles are set to 50∘ and 130∘, respectively.Then, space group (SG) and corresponding Wyckoff positions com-bination (WPC) are derived for building symmetry. Given the input che-mical composition, let Sþ be the set of SGs that are compatible with thestoichiometry. For each S 2 Sþ, letWS be the set of correspondingWPCs(see Supplementary Note 1). Let jWSj denote the number of distinctcompatibleWPCs for SG S. The process involves sampling an SGand then aWPC. An SG number is sampled from all compatible ones (Sþ) uniformly,i.e., an SG Sl is selected by sampling its identifier iSl uniformly from the set ofidentifiers for SGs in Sþ:iSl � UðfidðSÞjS 2 SþgÞ; Sl ¼ SGðiSl Þ: ð5ÞBased on Sl, a WPC is subsequently sampled. Define the maximum WPCcount Wmax ¼ maxS2SþjWSj. Sample an integer iWluniformly fromf0; 1; . . . ;Wmax � 1g:iWl� Uðf0; 1; . . . ;Wmax � 1gÞ: ð6ÞTheWPCWl is derived from theWPCs set corresponding to the chosen SGSl based on the sampled identifier:Wl ¼ WPCSliWl� jWSljWmax� �� �; ð7Þwhere ⌊x⌋ is the floor function, which returns the greatest integer less thanor equal to x.Finally, fractional coordinates are uniformly sampled from theinterval [0, 1). These coordinates are treated as independent sites andplaced into the Wyckoff positionsWl, where they are transformed intoexternal coordinates to satisfy symmetry constraints. Similarly, thegenerated lattice parameters are assigned to variables defined by thecrystal system (CS) Cl associated with the sampled SG Sl. Since WPCsimpose dependencies among coordinates, reducing the degrees offreedom, only the earliest generated coordinates are used. The sameapproach applies to lattice parameters constrained by a CS.When generating datasets for training, structures containing atompairs with distances smaller than 1.5Å are removed to ensure that mostgenerated structures have a reasonable estimated energy, which is essentialfor training an accurate objective function. Distance filtering is not involvedwhen evaluating the performance of the RG baseline. A much refined RGprocess is implemented by PyXtal52, which has been tested as the CRYSPYRG baseline in this work.Details of symmetry-informed integer encoding in CRYSIMInteger encoding can be interpreted as follows. Suppose a binary vectorsegment containingNv bits is leveraged for representing parameter v. For acontinuous parameter v 2 ½vmin; vmaxÞ, the iv-th bit of the vector segmentwill be assigned as 1 and other elements are 0s withiv ¼ v�vminuvj k; ð8Þhttps://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 8www.nature.com/commsphyswhere uv ¼ ðvmax � vminÞ=Nv is the unit or interval of the representation.For a discrete parameter, iv = 1 if the parameter v of the system correspondsto the i-th category.Lattice parameters encoding. In the workflow of CRYSIM, integerencoding is initially performed on training sets generated by RG. Theupper and lower bound of lattice length encoding, simultaneously thehighest and lowest lattice length of decoded materials, are calculatedbased on data points in the sets. Let ll0;max and ll0;min represent themaximum and minimum lattice lengths (for a, b, and c) among allstructures in the initial training set. Then, the upper and lower bounds forlattice length encoding are defined asllmax ¼ 1:1 � ll0;max;llmin ¼ ll0;min:(ð9ÞThe number of bits for representing lattice lengths Na = Nb = Nc = Nll isdependent on divisions of the lattice when encoding atomic coordinates,which is calculated byNll ¼ Cll � g �llmax � llminllmax; ð10Þinwhich gdenotes the current LDR, andCll=10bydefault. For lattice angles(α, β and γ), the lower and upper bound are 50∘ and 130∘, the same as RG.The unit for encoding angles is 2°, so that one lattice angle is encoded usingNlg ¼ ð130� 50Þ=2 ¼ 40 ð11Þbits. The number of bits for encoding lattice parameters would be3*Nll + 3*Nlg.Fractional coordinates encoding. Atomic configurations are repre-sented by discretizing the unit cell into a g × g × g voxel grid, where g is theLDR. A 3D binary matrix X ∈ {0, 1}g×g×g is constructed, where elementXl,m,n corresponds to the voxel region Rl,m,n defined by fractional coor-dinates y = (y1, y2, y3):Rl;m;n ¼ yj lg< y1 <l þ 1g;mg< y2 <mþ 1g;ng< y3 <nþ 1g� �ð12Þfor l, m, n ∈ {0, 1, …, g − 1}. Xl,m,n is set to 1 if an atom’s fractionalcoordinates fall within Rl,m,n, and 0 otherwise. For crystals containingmultiple element species, a separate flattenedmatrix is constructed for eachof them.After concatenation, encoded informationof eachelement is storedin separate regions of the final embedding.Wenote that the optimized coordinates are internal ones, whichwill beplaced into WPCs to satisfy symmetry constraints. Besides, similar to RG,the derived bits might be redundant. In implementation, 1-bits in the left-most positions in the vector segment for each element are used, until allvariables in the solvedWPC are decided. Special positions are not involvedin optimization, therefore not encoded in the discretized space. Once aspecific WPC is selected/decoded, the special positions participate in con-structing structures automatically. For instance, when filling the optimizedposition (0.2, 0.4, 0.8) into a WP:A1 : ð1=2; 1=2; 1=2Þ;A2 : ðx; 1=2; 1=2Þ;A3 : ð�x; 1=2; 1=2Þ;A4 : ðx; y; zÞ;A5 : ð�x;�y;�zÞ;8>>>>>><>>>>>>:ð13Þthe result is:A1 : ð1=2; 1=2; 1=2Þ;A2 : ð0:2; 1=2; 1=2Þ;A3 : ð�0:2; 1=2; 1=2Þ;A4 : ð0:2; 0:4; 0:8Þ;A5 : ð�0:2;�0:4;�0:8Þ:8>>>>>><>>>>>>:ð14ÞFor experiments on benchmark crystals, the Y6Co51 andCa24Al16(SiO4)24 system in this study, LDRs of CRYSIM are set to 12, with(SiO2)96 being 15.Symmetry information representation. Symmetry informationinvolves the CS, SG and WPC, which define the crystal’s symmetry.Numbers of bits for the three parts, i.e.,NC,NS andNW, are dependent onthe WPCs list calculated based on stoichiometry of the system. To bespecific, only compatible SGs and CSs are encoded, and whether an SGand CS is compatible or not is determined by the existence ofWPCs thatcan be used to build the corresponding symmetry, as is illustrated inSupplementary Note 1.Details of calculating NC,NS andNW are presented as follows. LetCþandSþ denote all compatibleCSs and SGs, respectively, andSC denotes theset of SGs associated with a CS C. We then define the set of compatible SGsfor C as SC;þ ¼ SC \ Sþ. The numbers of bits for representing CSs andSGs can be decided asNC ¼ jCþj;NS ¼ maxC2CþjSC;þj:(ð15ÞNW can be independently set. In many cases, the number of distinct com-patible WPCs for an SG S, denoted as jWSj, is so large that encoding allWPCs within a binary segment becomes impractical. To address this, onlyWPCs with indices within the set fbi�jWSjNWcji ¼ 0; 1; . . . ;NW � 1g areincluded. In this work, NW is set to 300 by default. However, duringdecoding, each bit i corresponds to a group of 100 WPCs with indi-ces fbi�jWSjNWc þ jjj ¼ 0; 1; . . . ; 99g.When encoding symmetry information of a crystal having CSCl, SG Sland WPC Wl, first of all, the corresponding bit representing the CS isassigned as 1. The bit for Sl is derived byiSl ¼ bidCl ;þðSlÞ �NSjSCl ;þjc; ð16Þin which idCl ;þðSlÞ is the identifier of Sl in the ascending ordered sequencewith respect to the set SCl ;þ. For instance, in the SG set of the cubic CS, i.e.,fP23; F23; :::; Ia3dg, the identifier of SG P23 is 0. Similarly, the bit forWl iscalculated byiWl¼ bidSl ðWlÞ �NWjWSljc: ð17Þwhere idSl ðWlÞ represents the identifier ofWl inWSl.During decoding processes, the following relationships are adoptedidCl ;þðSlÞ ¼ biSl �jSC;þjNSþ 0:5c; Sl ¼ SGCl ;þðidCl ;þðSlÞÞ;idSl ðWlÞ ¼ biWl� jWSjNWþ 0:5c; Wl ¼ WPCSlðidSl ðWlÞÞ;8<: ð18Þto make the encoding-decoding procedures stable and reversible.We note that according to the integer encoding strategy implementedin CRYSIM, only one bit should be assigned as 1 in a vector segment for oneparameter, so that decoding from the segment to real values can behttps://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 9www.nature.com/commsphysperformed directly. Amplify provides options to add constraints as penaltyterms to theobjective function to encourage generationof solutions fulfillingspecific requirements. In summary, the number of bits for symmetryencoding would be NC + NS + NW.An example on symmetry information representation. We furtherprovide an illustrative example on encoding and decoding symmetryinformation. For the A4B4 system, there are 2 SGs available for the tri-clinic CS, 12 for monoclinic, 56 for orthorhombic, 66 for tetragonal, 12for trigonal, 19 for hexagonal, and 17 for cubic. Since all CSs are com-patible with the system, NC = 7. The maximum number of compatibleSGs across all CSs is 66, and therefore NS = 66.NW can be independentlyset as 300 by default. Accordingly, the total number of bits for encodingsymmetry information is 7 + 66 + 300. When encoding SGs, as anexample, P41, the No.76 SG, is the second compatible SG belonging totetragonal CS, thus the 4th bit for CS and the 2nd bit for SG are set to 1. ForP23 (SG No.195) and F23 (SG No. 196), the first and second compatibleSG in the cubic category, the corresponding bit for SG is calculatedfollowingiP23 ¼ 0 � 6617� �¼ 0; ð19ÞandiF23 ¼ 1 � 6617� �¼ 3; ð20Þrespectively. On the other hand, in the decoding process, the bit 2, 3 and 4will correspond tob2 � 1766 þ 0:5c ¼ 0 ! P23;b3 � 1766 þ 0:5c ¼ 1 ! F23;b4 � 1766 þ 0:5c ¼ 1 ! F23:8><>: ð21ÞIn actual implementation, since it is impossible to determine WPCindex from an already generated structure, we only curate training setsobtained fromourRGalgorithm, inwhich crystals are constructedbased onan already chosen WPC.Criterion for matching structuresStructureMatcher function in pymatgen package57 is used tocompare configurations, in which parameters are set as stol=0.5,ltol=0.3,angle_tol=10.0, consistent with other relatedworks82,83.This function will calculate the minimum normalized average root meansquare pair-wise displacement between two input structures amongall atompermutations. But if corresponding atoms in the two structures are notdetected, which means that the function cannot identify any similaritybetween them, the calculation will not be proceeded. Accordingly, a struc-ture is recognized to be accordant with the ground truth if having a com-putable displacement with it, and a model successfully finds the groundtruth in one run if there is at least one such structure being generated.Factorization machine for quadratic regressionFactorization machine (FM) is a type of regression model proposed as asubstitute of Support VectorMachine to address its failure on sparse data37.FM creates a mapping between a vector x 2 RM and real value y byequation (2), where b, hi, wk,i are coefficients for bias, linear and quadraticinteractions, respectively. In principle, FM can be extended to model n-interaction terms, but we restrict it to quadratic terms since current com-binatorial optimizers are efficient only for solving quadratic objectivefunctions. In that case, FM can be reformulated asy ¼ bþXMi¼1hixi þ12XKk¼1XMi¼1wkixi !2�XMi¼1w2kix2i !; ð22Þreducing computational complexity from O(KM) to O(2K)37. One ofadvantages of FM in this work is that it requires less fitting parameters,enabling a quadratic regression on binary vectors containing thousands ofbits. Taking a vector of 2000 bits as an example, a full-rank quadraticregressor requires 2000 × 1999 terms for interactions, while FM only needs2000 × K terms, in which K is usually smaller than 30. In this work, weimplement FM with PyTorch53 based on equation (22) to the acceleratelearning process.Data availabilityGround state configurations considered in this study can be downloadedfrom theMP database55. Initial datasets for training FM are generated usingRG implemented in CRYSIM, and no external data is included. All relevantdata is available from the authors upon request. The original data for thefigures in this paper and Supplementary Information is stored in Supple-mentary Data.Code availabilityImplementation of CRYSIM is available at https://github.com/tsudalab/CRYSIM. As of March 2025, Fixstars Amplify is available via Python APIfree of charge, which can be accessed at https://amplify.fixstars.com/en/.Received: 13 April 2025; Accepted: 19 October 2025;References1. de Pablo, J. J. et al. New frontiers for the materials genome initiative.npj Comput. Mater. 5, 41 (2019).2. Woodley, S. M. & Catlow, R. Crystal structure prediction from firstprinciples. Nat. Mater. 7, 937–946 (2008).3. Oganov, A. R., Pickard, C. J., Zhu, Q. & Needs, R. J. Structureprediction drives materials discovery. Nat. Rev. Mater. 4, 331–348(2019).4. Wang, Y., Lv, J., Gao, P. & Ma, Y. Crystal structure prediction viaefficient samplingof thepotential energy surface.Acc.Chem.Res.55,2068–2076 (2022).5. Bush, T. S., Catlow, C. R. A. & Battle, P. D. Evolutionary programmingtechniques for predicting inorganic crystal structures.J.Mater.Chem.5, 1269–1272 (1995).6. M. Woodley, S., D. Battle, P., D. Gale, J. & Richard A. Catlow, C. Theprediction of inorganic crystal structures using a genetic algorithmand energy minimisation. Phys. Chem. Chem. Phys. 1, 2535–2542(1999).7. Glass,C.W., Oganov, A. R. &Hansen,N. Uspex—evolutionary crystalstructure prediction.Computer Phys. Commun. 175, 713–720 (2006).8. Lonie, D. C. & Zurek, E. Xtalopt: An open-source evolutionaryalgorithm for crystal structure prediction. Comput. Phys. Commun.182, 372–387 (2011).9. Oganov, A. R., Lyakhov, A. O. & Valle, M. How evolutionary crystalstructure prediction works—and why. Acc. Chem. Res. 44, 227–237(2011).10. Wang, Y., Lv, J., Zhu, L. & Ma, Y. Crystal structure prediction viaparticle-swarm optimization. Phys. Rev. B 82, 094116 (2010).11. Wang, Y., Lv, J., Zhu, L. & Ma, Y. Calypso: A method for crystalstructure prediction. Comput. Phys. Commun. 183, 2063–2070(2012).12. Yamashita, T. et al. Crystal structure prediction accelerated bybayesian optimization. Phys. Rev. Mater. 2, 013803 (2018).https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 10https://github.com/tsudalab/CRYSIMhttps://github.com/tsudalab/CRYSIMhttps://amplify.fixstars.com/en/www.nature.com/commsphys13. Yamashita, T. et al. Cryspy: a crystal structure prediction toolaccelerated by machine learning. Sci. Technol. Adv. Mater. 1, 87–97(2021).14. Yamashita, T., Kino, H., Tsuda, K., Miyake, T. & Oguchi, T. Hybridalgorithm of bayesian optimization and evolutionary algorithm incrystal structure prediction. Sci. Technol. Adv. Mater.: Methods 2,67–74 (2022).15. Xie, T., Fu, X., Ganea, O.-E., Barzilay, R. & Jaakkola, T. S. Crystaldiffusion variational autoencoder for periodic material generation. InInternational Conference on Learning Representations (2022).16. Merchant, A. et al. Scaling deep learning for materials discovery.Nature 624, 80–85 (2023).17. Jiao, R., Huang, W., Liu, Y., Zhao, D. & Liu, Y. Space groupconstrained crystal generation. In The Twelfth InternationalConference on Learning Representations (2024).18. Luo, X. et al. Deep learning generative model for crystal structureprediction. npj Comput. Mater. 10, 254 (2024).19. Chang, L. et al. Shotgun crystal structure prediction using machine-learned formation energies. npj Comput. Mater. 10, 298 (2024).20. Zhu, R., Nong, W., Yamazaki, S. & Hippalgaonkar, K. Wycryst: Wyckoffinorganic crystal generator framework. Matter 7, 3469–3488 (2024).21. Zeni, C. et al. A generative model for inorganic materials design.Nature, 639, 624–632. https://doi.org/10.1038/s41586-025-08628-5(2025).22. Singh, A. K. et al. Genetic algorithm prediction of two-dimensionalgroup-iv dioxides for dielectrics. Phys. Rev. B 95, 155426 (2017).23. M. Dieb, T., Hou, Z. & Tsuda, K. Structure prediction of boron-dopedgraphene by machine learning. J. Chem. Phys. 148, 241716 (2018).24. Walsh, A. & Woodley, S. M. Evolutionary structure prediction andelectronic properties of indium oxide nanoclusters. Phys. Chem.Chem. Phys. 12, 8446–8453 (2010).25. Wang, Y. et al. Accelerated prediction of atomically precise clusterstructures using on-the-fly machine learning. npj Comput. Mater. 8,173 (2022).26. Wang, H. et al. Seeded growth of single-crystal black phosphorusnanoribbons. Nat. Mater. 23, 470–478 (2024).27. Tanaka, S., Tamura, R. & Chakrabarti, B. K.Quantum spin glasses,annealing and computation (Cambridge University Press, 2017).28. Mohseni, N.,McMahon, P. L. &Byrnes, T. Isingmachines as hardwaresolvers of combinatorial optimization problems. Nat. Rev. Phys. 4,363–379 (2022).29. Gusev, V. V. et al. Optimality guarantees for crystal structureprediction. Nature 619, 68–72 (2023).30. Ichikawa, K., Ohuchi, S., Ueno, K. & Yokoyama, T. Acceleratingoptimal elemental configuration search in crystal using Isingmachine.Phys. Rev. Res. 6, 033321 (2024).31. Couzinié, Y. et al. Annealing for prediction of grand canonical crystalstructures: Implementation of n-body atomic interactions. Phys. Rev.A 109, 032416 (2024).32. Xu, Z., Shang, W., Kim, S., Lee, E. & Luo, T. Quantum annealing-assisted lattice optimization. npj Comput. Mater. 11, 4 (2025).33. Couzinié, Y. et al. Machine learning supported annealing for theprediction of grand canonical crystal structures. J. Phys. Soc. Jpn. 94,044802 (2025).34. Lyakhov, A. O., Oganov, A. R., Stokes, H. T. & Zhu, Q. Newdevelopments in evolutionary structure prediction algorithm uspex.Computer Phys. Commun. 184, 1172–1182 (2013).35. McGeoch, C. C., Harris, R., Reinhardt, S. P. & Bunyk, P. I. Practicalannealing-based quantum computing. Computer 52, 38–46 (2019).36. Chen, C. & Ong, S. P. A universal graph deep learning interatomicpotential for the periodic table. Nat. Comput. Sci. 2, 718–728 (2022).37. Rendle, S. Factorization machines. In 2010 IEEE InternationalConference on Data Mining, 995–1000, https://doi.org/10.1109/ICDM.2010.127 (2010).38. Kitai, K. et al. Designing metamaterials with quantum annealing andfactorization machines. Phys. Rev. Res. 2, 013319 (2020).39. Guo, J., Kitai, K., Jippo, H. & Shiomi, J. Boosting the quality factor oftamm structures to millions by quantum-inspired classical annealerwith factorization machine 2408.05799 (2024).40. Xu, Z. et al. Quantum-inspired genetic algorithm for designing planarmultilayer photonic structure. npj Comput. Mater. 10, 257 (2024).41. Podryabinkin, E. V., Tikhonov, E. V., Shapeev, A. V. & Oganov, A. R.Accelerating crystal structure prediction by machine-learninginteratomic potentials with active learning. Phys. Rev. B 99, 064114(2019).42. Fixstars amplify. https://amplify.fixstars.com/en/. Accessed: 2024-11-26.43. Mao, Z., Matsuda, Y., Tamura, R. & Tsuda, K. Chemical design withGPU-based Ising machines. Digital Discov. 2, 1098–1103 (2023).44. Endo, K., Matsuda, Y., Tanaka, S. & Muramatsu, M. A phase-fieldmodel by an Isingmachineand its application to thephase-separationstructure of a diblock polymer. Sci. Rep. 12, 10794 (2022).45. Urushihara,M., Karube,M., Yamaguchi, K. & Tamura, R. Optimizationof core-shell nanoparticles using a combination of machine learningand Ising machine. Adv. Photonics Res. 4, 2300226 (2023).46. Tamura, R. et al. Machine learning prediction of the mechanicalproperties of injection-molded polypropylene through X-raydiffraction analysis. Sci. Technol. Adv. Mater. 25, 2388016 (2024).47. Kagemoto, H. Possible application of quantum computing in the fieldof ocean engineering: optimization of an offshore wind farm layoutwith the Ising model. J. Ocean Eng. Mar. Energy 10, 773–782 (2024).48. Sukulthanasorn, N. et al. A novel design update framework fortopology optimization with quantum annealing: Application to trussand continuum structures.ComputerMethods Appl. Mech. Eng. 437,117746 (2025).49. Cheng, G., Gong, X.-G. & Yin, W.-J. Crystal structure prediction bycombining graph network and optimization algorithm. Nat. Commun.13, 1492 (2022).50. Deng,X. &Dong,C.SMEPOC– acomputer program for the automaticgeneration of trial structural models for inorganic compounds withsymmetry restriction. J. Appl. Crystallogr. 42, 953–958 (2009).51. Avery, P. & Zurek, E. Randspg: An open-source program forgenerating atomistic crystal structures with specific space groups.Computer Phys. Commun. 213, 208–216 (2017).52. Fredericks, S., Parrish, K., Sayre, D. & Zhu, Q. Pyxtal: A Python libraryfor crystal structure generation and symmetry analysis. Comput.Phys. Commun. 261, 107810 (2021).53. Paszke, A. et al. Pytorch: An imperative style, high-performance deeplearning library. In Wallach, H. et al. (eds.) Advances in NeuralInformation Processing Systems, vol. 32 (Curran Associates, Inc.,2019).54. Wei, L., Li, Q., Omee, S. S. & Hu, J. Towards quantitative evaluation ofcrystal structure prediction performance. Comput. Materials Science235, https://doi.org/10.1016/j.commatsci.2024.112802 (2024).55. Jain, A. et al. Commentary: Thematerials project: Amaterials genomeapproach to accelerating materials innovation. APL Mater. 1, 011002(2013).56. Bergstra, J., Komer, B., Eliasmith, C., Yamins, D. & Cox, D. D.Hyperopt: a Python library for model selection and hyperparameteroptimization. Comput. Sci. Discov. 8, 014008 (2015).57. Ong, S. P. et al. Python materials genomics (pymatgen): A robust,open-source Python library for materials analysis. Comput. Mater.Sci. 68, 314–319 (2013).58. Deng, B. et al. Systematic softening in universal machine learninginteratomic potentials. npj Comput. Mater. 11, 9 (2025).59. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energycalculations formetals and semiconductors using a plane-wavebasisset. Comput. Mater. Sci. 6, 15–50 (1996).https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 11https://doi.org/10.1038/s41586-025-08628-5https://doi.org/10.1038/s41586-025-08628-5https://doi.org/10.1109/ICDM.2010.127https://doi.org/10.1109/ICDM.2010.127https://doi.org/10.1109/ICDM.2010.127https://amplify.fixstars.com/en/https://amplify.fixstars.com/en/https://doi.org/10.1016/j.commatsci.2024.112802https://doi.org/10.1016/j.commatsci.2024.112802www.nature.com/commsphys60. Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to theprojector augmented-wave method. Phys. Rev. B 59, 1758–1775(1999).61. Terayama, K., Yamashita, T., Oguchi, T. & Tsuda, K. Fine-grainedoptimization method for crystal structure prediction. npj Comput.Mater. 4, 32 (2018).62. Ishikawa, T., Fukazawa, T., Xing, G., Tadano, T. & Miyake, T.Evolutionary search for cobalt-rich compounds in the yttrium-cobalt-boron system. Phys. Rev. Mater. 5, 054408 (2021).63. Lyakhov,A.O.,Oganov, A. R. &Valle,M.How topredict very large andcomplex crystal structures. Computer Phys. Commun. 181,1623–1632 (2010).64. Frazier, P. I. A tutorial on Bayesian optimization 1807.02811 (2018).65. Moriconi, R., Deisenroth,M. P. & SeshKumar, K. S. High-dimensionalBayesian optimization using low-dimensional feature spaces.Mach.Learn. 109, 1925–1943 (2020).66. Eriksson, D. & Jankowiak, M. High-dimensional Bayesianoptimization with sparse axis-aligned subspaces. In de Campos, C. &Maathuis, M. H. (eds.) Proceedings of the Thirty-Seventh Conferenceon Uncertainty in Artificial Intelligence, vol. 161 of Proceedings ofMachine Learning Research, 493–503 (PMLR, 2021).67. Oganov, A.R. &Valle,M.How toquantify energy landscapesof solids.J. Chem. Phys. 130, 104504 (2009).68. Pickard, C. J. & Needs, R. J. Ab initio random structure searching. J.Phys. Condens. Matter 23, 053201 (2011).69. Wales, D. J. Symmetry, near-symmetry and energetics. Chem. Phys.Lett. 285, 330–336 (1998).70. Urusov, V. S. & Nadezhina, T. N. Frequency distribution and selectionof space groups in inorganic crystal chemistry. J. Struct. Chem. 50,22–37 (2009).71. Avery, P., Toher, C., Curtarolo, S. & Zurek, E. Xtalopt version r12: Anopen-source evolutionary algorithm for crystal structure prediction.Computer Phys. Commun. 237, 274–275 (2019).72. Levy,D. et al. SymmCD:Symmetry-preservingcrystal generationwithdiffusion models. In The Thirteenth International Conference onLearning Representations (2025).73. Kelvinius, F. E. et al. Wyckoffdiff – a generative diffusion model forcrystal symmetry. In Forty-second International Conference onMachine Learning (2025).74. Cao, Z., Luo, X., Lv, J. & Wang, L. Space group informed transformerfor crystalline materials generation. Science Bulletin https://doi.org/10.1016/j.scib.2025.09.035 (2025).75. Kazeev, N. et al. Wyckoff transformer: Generation of symmetriccrystals. In Forty-second International Conference on MachineLearning (2025).76. Finnila, A. B., Gomez, M. A., Sebenik, C., Stenson, C. & Doll, J. D.Quantum annealing: A new method for minimizing multidimensionalfunctions. Chem. Phys. Lett. 219, 343–348 (1994).77. Kadowaki, T. & Nishimori, H. Quantum annealing in the transverseIsing model. Phys. Rev. E 58, 5355–5363 (1998).78. Johnson, M. W. et al. Quantum annealing with manufactured spins.Nature 473, 194–198 (2011).79. Albash,T.&Lidar,D.A.Adiabaticquantumcomputation.Rev.ModernPhys. 90, https://doi.org/10.1103/RevModPhys.90.015002 (2018).80. Hauke,P., Katzgraber, H.G., Lechner,W.,Nishimori, H. &Oliver,W.D.Perspectives of quantum annealing: methods and implementations.Rep. Prog. Phys. 83, 054401 (2020).81. King,A.D. et al. Coherentquantumannealing in aprogrammable2000qubit Ising chain. Nat. Phys. 18, 1324–1328 (2022).82. Jiao, R. et al. Crystal structure predictionby joint equivariant diffusion.In Thirty-seventh Conference on Neural Information ProcessingSystems (2023).83. Sriram, A., Miller, B. K., Chen, R. T. Q. & Wood, B. M. FlowLLM: Flowmatching for material generation with large language models as basedistributions. In The Thirty-eighth Annual Conference on NeuralInformation Processing Systems (2024).84. Momma, K. & Izumi, F. VESTA3 for three-dimensional visualization ofcrystal, volumetric, and morphology data. J. Appl. Crystallogr. 44,1272–1276 (2011).AcknowledgementsThis work is supported by JST ERATO JPMJER1903, JST CRESTJPMJCR21O2, and JSPS KAKENHI Young Scientist (23K16942). C.L.would like to gratefully acknowledge the financial support from the ChinaScholarship Council (CSC No. 202306210120). The authors thank YaotangZhang for discussions.Author contributionsC.L. implemented the CRYSIM package and conducted all experiments.D.D. contributed to insights into the design of Ising models and constraints.Z.M.contributed to implementationof the training frameworkofFM. J.G.andR.T. contributed to analysis about the usage of Amplify and other Isingsolvers.K.T. proposed the ideaof thework.Z.M.andK.T.providedguidanceon experiments design and results analysis. All authors reviewed andcontributed to the writing of the manuscript.Competing interestsThe authors declare no competing interests.Additional informationSupplementary information The online version containssupplementary material available athttps://doi.org/10.1038/s42005-025-02380-y.Correspondence and requests for materials should be addressed toZetian Mao or Koji Tsuda.Peer review information Communications Physics thanks Wei Nong,Kedar Hippalgaonkar and the other, anonymous, reviewer(s) for theircontribution to the peer review of this work. A peer review file is available.Reprints and permissions information is available athttp://www.nature.com/reprintsPublisher’s note Springer Nature remains neutral with regard tojurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative CommonsAttribution-NonCommercial-NoDerivatives 4.0 International License,which permits any non-commercial use, sharing, distribution andreproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the CreativeCommons licence, and indicate if you modified the licensed material. Youdo not have permission under this licence to share adapted materialderived from this article or parts of it. The images or other third partymaterial in this article are included in the article’s Creative Commonslicence, unless indicated otherwise in a credit line to thematerial. If materialis not included in thearticle’sCreativeCommons licenceandyour intendeduse is not permitted by statutory regulation or exceeds the permitted use,you will need to obtain permission directly from the copyright holder. Toview a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.© The Author(s) 2025https://doi.org/10.1038/s42005-025-02380-y ArticleCommunications Physics |           (2025) 8:477 12https://doi.org/10.1016/j.scib.2025.09.035https://doi.org/10.1016/j.scib.2025.09.035https://doi.org/10.1016/j.scib.2025.09.035https://doi.org/10.1103/RevModPhys.90.015002https://doi.org/10.1103/RevModPhys.90.015002https://doi.org/10.1038/s42005-025-02380-yhttp://www.nature.com/reprintshttp://creativecommons.org/licenses/by-nc-nd/4.0/http://creativecommons.org/licenses/by-nc-nd/4.0/www.nature.com/commsphys Predicting symmetric structures of large crystals with GPU-based Ising machines Results Bit vector encoding CRYSIM Workflow Recovering benchmark materials Large crystal structures prediction Effects of Processing Techniques Lattice Discretization Resolution Factorization machine Processing with minimum interatomic distance Discussion Conclusion Methods Random generation of crystal structures Details of symmetry-informed integer encoding in CRYSIM Lattice parameters encoding Fractional coordinates encoding Symmetry information representation An example on symmetry information representation Criterion for matching structures Factorization machine for quadratic regression Data availability Code availability References Acknowledgements Author contributions Competing interests Additional information