# Fileset

[Surface-structure search with variable composition and periodicity via machine learning and evolutionary algorithms  applications to Pt Ge oxidation a.pdf](https://mdr.nims.go.jp/filesets/a095dd7f-083a-4c9b-93b8-0e02f6a8686a/download)

## Creator

F. Kuroda, M. Otani

## Rights

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

## Other metadata

[Surface-Structure Search with Variable Composition and Periodicity via Machine Learning and Evolutionary Algorithms: Applications to Pt/Ge Oxidation and Au--Sn Alloying](https://mdr.nims.go.jp/datasets/c8c9d667-d954-4296-a6f3-583c66cabb8a)

## Fulltext

Microsoft Word - AMO_TSTA_A_2688059.docxScience and Technology of Advanced MaterialsISSN: 1468-6996 (Print) 1878-5514 (Online) Journal homepage: www.tandfonline.com/journals/tsta20Surface-structure search with variablecomposition and periodicity via machine learningand evolutionary algorithms: applications to Pt/Ge oxidation and Au–sn alloyingF. Kuroda & M. OtaniTo cite this article: F. Kuroda & M. Otani (07 Jul 2026): Surface-structure search with variablecomposition and periodicity via machine learning and evolutionary algorithms: applicationsto Pt/Ge oxidation and Au–sn alloying, Science and Technology of Advanced Materials, DOI:10.1080/14686996.2026.2688059To link to this article:  https://doi.org/10.1080/14686996.2026.2688059© 2026 The Author(s). Published by NationalInstitute for Materials Science in partnershipwith Taylor & Francis Group.Accepted author version posted online: 07Jul 2026.Submit your article to this journal View related articles View Crossmark dataFull Terms & Conditions of access and use can be found athttps://www.tandfonline.com/action/journalInformation?journalCode=tsta20https://www.tandfonline.com/journals/tsta20?src=pdfhttps://www.tandfonline.com/action/showCitFormats?doi=10.1080/14686996.2026.2688059https://doi.org/10.1080/14686996.2026.2688059https://www.tandfonline.com/action/authorSubmission?journalCode=tsta20&show=instructions&src=pdfhttps://www.tandfonline.com/action/authorSubmission?journalCode=tsta20&show=instructions&src=pdfhttps://www.tandfonline.com/doi/mlt/10.1080/14686996.2026.2688059?src=pdfhttps://www.tandfonline.com/doi/mlt/10.1080/14686996.2026.2688059?src=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1080/14686996.2026.2688059&domain=pdf&date_stamp=07%20Jul%202026http://crossmark.crossref.org/dialog/?doi=10.1080/14686996.2026.2688059&domain=pdf&date_stamp=07%20Jul%202026https://www.tandfonline.com/action/journalInformation?journalCode=tsta20ACCEPTED MANUSCRIPTSurface-Structure Search with Variable Composition and Periodicity via Machine Learning and Evolutionary Algorithms: Applications to Pt/Ge Oxidation and Au–Sn Alloying  F. Kurodaa and M. Otanib  a Materials DX Research Center, National Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono, Tsukuba, Ibaraki 305-8560, Japan,; bCenter for Computational Sciences (CCS), University of Tsukuba, 1-1-1, Tenno-dai, Tsukuba, Ibaraki 305-8577, Japan  ARTICLE HISTORY Compiled May 23, 2026 ABSTRACT First-principles structure prediction is essential for discovering functional materials; however, surface structure searches remain challenging because most search algorithms assume fixed in-plane periodicity and composition. Here we develop a global search framework that treats both two-dimensional superlattice periodicity and stoichiometry as dynamic variables, enabling the direct identification of the most stable surface structures across competing supercell shapes and compositions. The proposed method integrates an evolutionary algorithm with surface-specific variation operators and symmetry-enriched initialization, and accelerates screening via Bayesian optimization using atomic cluster expansion descriptors. Case studies on FCC Pt(111) and diamond Ge(100) surfaces yield oxygen-induced surface structures consistent with experimental observations, and the same framework identifies Sn alloying motifs on FCC Au(111) that agree with reported surface-structure trends. Overall, the framework delivers accurate structure prediction with substantially fewer high-cost first-principles evaluations and provides a general route to exploring complex materials landscapes–such as heterogeneous catalysis, electronics, and spintronics–in which coupled structural and compositional degrees of freedom govern functionality.  KEYWORDS Sections; lists; figures; tables; mathematics; fonts; references; appendices  CONTACT . F. Kuroda. Email: fumiaki.kuroda@aist.go.jp Publisher: Taylor & Francis & The Author(s). Published by National Institute forMaterials Science in partnership with Taylor & Francis Group.  Journal: Science and Technology of Advanced Materials  DOI: 10.1080/14686996.2026.2688059https://crossmark.crossref.org/dialog/?doi=10.1080/14686996.2026.2688059&domain=pdfACCEPTED MANUSCRIPT1.  Introduction  Atomic structure is the most fundamental information in crystalline materials, determining their properties and functionalities. Once a structural topology is specified, a precise structural model and a wide range of physical properties can be computed using modern quantum-mechanical methods. In particular, the atomic structure of surfaces and interfaces plays a central role in functional materials, governing key processes in heterogeneous catalysis [4,5] and enabling emergent phenomena in spintronics [1–3]. However, identifying the atomic structures realized experimentally is often challenging, especially for surfaces and interfaces, where structural motifs can be diverse and complex. This difficulty stems from the vastness of the configurational space [6]: even for relatively small systems, the number of distinct local minima increases rapidly with the number of atoms, leading to a rugged energy landscape populated by numerous metastable structures. As a result, exhaustive enumeration is intractable, and local relaxation can easily become trapped far from the global minimum. The problem is further complicated in practical surface searches (Fig. 1), where the in-plane periodicity and stoichiometry are frequently fixed a priori. While such constraints simplify the search, they can introduce a strong bias and may preclude the true ground state when the most stable structure occurs at a different superlattice periodicity, coverage, or composition. Consequently, a reliable framework should explore not only atomic coordinates but also the degrees of freedom associated with periodicity and stoichiometry, enabling unbiased competition among candidate reconstructions. Crystal structure prediction (CSP) aims to identify low-energy atomic configurations on a highly rugged potential-energy surface [7,8]. To explore such complex landscapes efficiently, many CSP frameworks adopt population-based global search strategies for structure generation. Among them, particle-swarm optimization (PSO) updates a population of “particles” using collective search dynamics; CALYPSO (Crystal structure AnaLYsis by Particle Swarm Optimization) [9–11] is a CSP package explicitly built on PSO, whereas evolutionary algorithms (EA) generate new candidates through variation operators such as crossover, mutation. A representative EA-based framework is USPEX (Universal Structure Predictor: Evolutionary Xtallography) [6,12–15], which is characterized by its rich set of variation EA operators for generating new candidate structures. Furthurmore, programs that combine EA with Bayesian optimization (BO) include CrySPY (Crystal Structure prediction in Python) [16,17] and GOFEE (Global Optimization with First-principles Energy Expressions) [18,19]. These frameworks accelerate the search for low-energy structures by generating candidate configurations via EA and then selectively prioritizing expensive evaluations (such as first-principles calculations) using BO guided by a surrogate model, thereby reducing the total number of energy/relaxation calculations required. USPEX, CALYPSO, and GOFEE can also be applied to slab models, enabling CSP on surfaces. In typical workflows, however, the degrees of freedom treated as variables differ across codes. GOFEE and CALYPSO usually explore structures within a fixed simulation cell and fixed stoichiometry, optimizing mainly the atomic coordinates in the prescribed slab supercell. In contrast, USPEX can perform variable-composition searches for slab/surface problems [15], allowing the stoichiometry to change during the search, while the in-plane surface cell (periodicity) is generally kept fixed to the user-defined slab supercell.  Although user-defined slab supercells reduce computational cost, they require manual preparation of multiple candidate cells to examine different surface periodicities. Because the ACCEPTED MANUSCRIPTcombinations of surface superlattices and compositions increase combinatorially, this fixed-cell treatment can substantially limit the accessible structure-search space. Further details are provided in Appendix A.  In this study, we propose a surface-structure prediction method that simultaneously explores atomic coordinates, surface superlattice periodicity, and chemical composition within a unified search space. The method couples an EA equipped with surface-specific variation operators and symmetry-enriched initialization to a Bayesian-optimization layer, thereby accelerating the discovery of low-energy reconstructions under less restrictive constraints. The paper is outlined as follows. First, we describe the details of the proposed method in Secs. 2.1–2.4, including random structure generation, parent selection, evolutionary operations, and candidate evaluation via BO. Then, we present three case studies to which the proposed framework is applied. In Sec. 2.5, we briefly describe the structural relaxation procedure and the energy evaluation settings used for these three case studies. As the first of the three case studies, we examine oxidation on FCC Pt(111), which is of central importance for heterogeneous catalysis and electrocatalysis. In this case study, we also discuss whether the BO-driven CSP scheme provides an effective optimization strategy for identifying low-energy surface structures. In addition, we discuss surface reconstruction and oxidation on diamond Ge(100), a post-Si semiconductor, and finally examine Au-Sn alloys on the FCC Au(111) surface, which have attracted attention as catalytic materials and as Rashba-splitting systems.   2.  Method  First of this section provides an computational workflow (Fig. 2) of the proposed surface-structure search framework, named BACCHUS (Bayesian Atomic Cluster Construction High-efficiently Unveiling of Surface-structures). In BACCHUS, (i) we first generate an initial pool of random structures by varying the supercell (in-plane periodicity) of slabs and adsorption configurations on slabs. (ii) Each structure is locally relaxed, and its total energy is evaluated using a chosen calculation engine, such as first-principles (FP) methods or neural network potentials (NNPs). (iii) From the resulting energies, we construct a convex hull and quantify the thermodynamic stability of each structure by its distance to the hull (energy above hull). (iv) Promising structures are then selected through a combination of fit-based screening, elite selection, and roulette-wheel selection, and are used as parents in an evolutionary algorithm. Offspring structures are generated via variation operators such as mutation and crossover, yielding on the order of 310  new candidates per generation. (v)–(vi) Using the accumulated dataset, we next build an Gaussian-process regression model based on atomic cluster expansion (ACE) surrogate model and employ BO to prioritize the next candidates for exact energy evaluations. By iterating this loop from (ii) to (vi), BACCHUS efficiently explores the surface-structure space while substantially reducing the number of expensive FP calculations required. Then, we describe each component of the workflow in detail from Sec. 2.1 to Sec. 2.5.  2.1.  Random structure generation  Initial structures are generated by random sampling. In BACCHUS, we treat both the in-plane supercell periodicity and the stoichiometry as dynamic variables. Therefore, we first generate random supercells by varying symmetrically distinct lattice vectors within user-defined supercell size (see Sec. 2.1.1). Next, we place adsorbates on the slab surface by identifying adsorption sites ACCEPTED MANUSCRIPTbased on the symmetry of the chosen superlattice (see Sec. 2.1.2).   2.1.1.  Supercell construction  We generate derivative superlattices from a given parent slab structure using the method proposed by A. Santoro [20,21] and G. L. W. Hart [22]. Starting from a parent cell of arbitrary lattice type, the procedure first enumerates all symmetry-distinct derivative superlattices, which are then used to construct the corresponding derivative structures. Consider the transformation   = ,B HA  (1)  where T1 2= ( , )a aA  is a 2 2×  matrix whose row vectors are the parent in-plane lattice vectors ,  = 1,2i ia . H  is a 2 2×  integer matrix, and T1 2= ( , )b bB  is the matrix of the superlattice vectors ,  = 1, 2i ib . When the determinant of H  is | |= 1±H , H  simplly rotates the parent cell. This transformation H  can be represented by a Hermite normal form (HNF) matrix:   1121 220= ,hh h   H  (2)  where 11 22, > 0h h  and 21 220 <h h≤ . The | |H  is given by 11 22h h , which corresponds to the supercell size Sn  (the number of parent cells contained in the supercell). We can systematically enumerate all possible HNF matrices for a given Sn  by finding all integer pairs 11 22 21( , , )h h h  that satisfy the above conditions. Then, we further reduce the number of superlattices by considering the symmetry of the parent lattice. In BACCHUS, the in-plane point group symmetry of the parent slab is identified using the spglib library [23], and we eliminate rotational operations { }R . If one obtains two supercell matrices with the same Sn , 1 1=B H A  and 2 2=B H A . These two supercells are considered equivalent, if they satisfy the relation T T 11 2= ( ) −C RB B  , where C  is an unimodular matrix and R  is a rotational operation of the cartesionian coordinates. For example (also see Appendix A), we consider the point group system of the FCC (111) surface. For S = 3n , we can generate four HNF matrices from the above procedure. However, two diagonal matrices in these four are found to be equivalent by the rotational symmetry (Fig. 3).  2.1.2.  Symmetry-enriched initialization  Initial structures are generated by placing symmetric sites on the slab surface of each supercell. We fist identify point group symmetry of the supercell surface using following relation: 1S = −R H RH , where SR  should be an integer matrix if R  is a symmetry operation of the given supercell. From this relation, we can construct the point symmetry group of the supercell surface. Moreover, the symmetry subgroup of the supercell surface can be obtained. Next, we discretize the supercell surface with a uniform two-dimensional mesh and compile a list of candidate positions. Using the symmetry group (including a selected subgroup) obtained above, we classify these mesh points into symmetry-equivalent sets. We then choose point sets that remain invariant under the symmetry operations and place adsorbates on them to generate symmetry-preserving initial configurations. This symmetry-enriched initialization allows us to sample adsorption ACCEPTED MANUSCRIPTstructures that respect the supercell surface symmetry by construction. Figure 4 shows examples of symmetry-enriched site generation on a supercell surface.   2.2.  Parent structures selection  After evaluating the energies of all structures in the current generation, we select parent structures for the evolutionary operations. We first construct the convex hull of formation energies from the evaluated structures. The formation energy formE  of an adsorbate structure containing AN  atoms of element A and BN  atoms of element B is defined as   form slab ads S slab A A B B= ,E E n E N Nμ μ+ − − −  (3)  where slab adsE +  is the total energy of the slab with adsorbates, slabE  is the energy per parent cell of the clean slab, Sn  is the supercell size (i.e., the number of parent cells in the supercell), and Aμ  ( Bμ ) is the chemical potential of adsorbates of element A (B). Unless otherwise noted, Aμ  and Bμ  are set to the energies per atom in their most stable phases. From the convex hull, we compute the energy above hull, referred to as the hull distance hdE , for each structure as a measure of its thermodynamic stability (Fig. 5) [24]. After computing hdE  for all structures, we select parent structures using a combination of three schemes: elite , fittest , and roulette . In elite  selection, we select the structures with the lowest formation energies. This scheme retains the most stable structures across all generations. In contrast, fittest  and roulette  selection are applied to structures in the current generation, where stN  structures are available after local relaxation and energy evaluation. In fittest  selection, we introduce two threshold parameters, minhdE  and fitr . Here, minhdE  is a minimum hull-distance threshold used to exclude structures that are too similar to the candidates already retained by elite  selection. fit fit(0 < < 1)r r  is the fittest ratio: we sort the stN  structures in ascending order of hdE  and select the lowest- hdE  top st fitN r  structures. In roulette  selection, rouN  structures are selected from st fitN r  structures probabilistically based on their hdE  values. The selection probability ip  of structure i  is defined as   = ,iiiifpf′′ (4)  max minmax min max min= ,i ibf afa bf ff f f f−−′ +− − (5)  hd= 1* ( ),if E i−  (6)  where maxf  and minf  are the maximum and minimum values of if  in the population, and a  and b  are user-defined parameters.  2.3.  Offspring structure generation  The purpose of the evolutionary operations is to generate more global structural diversity while retaining locally stable atomic motifs inherited from low-energy parents. Using the parent ACCEPTED MANUSCRIPTstructures selected in Sec. 2.2, we generate offspring structures via evolutionary operations. In BACCHUS, we implement crossover and eleven mutation operators (Fig. 6): single-permutation, multi-permutation, substitution, moving, building, reflection, vertical-reflection, modulation, chair-switching, rattle, and phonon-rattle. For crossover , two parent structures are each split approximately in half, and the resulting fragments are combined (stitched together) to create offspring structures. For single -permutation  and multi - permutation , we randomly select one or more adsorbates and change their element types to other species. In substitution , we randomly select an adsorbate and replace it with a different adsorbate species. In moving , we randomly select an adsorbate and significantly change its lateral position while keeping its z  coordinate fixed. In building , we randomly select adsorbates and stack them on top of another adsorbate to create a cluster-like structure. In reflection , we reflect the positions of adsorbates with respect to a vertical plane passing through a surface atom. In vertical - reflection , we reflect the adsorbate positions with respect to a horizontal plane passing through a surface atom. In modulation , we apply a sinusoidal modulation to the adsorbate positions along a chosen direction. In chair - switching , we switch the adsorbate configuration by using a different supercell lattice. In rattle , we randomly displace the adsorbates by a distance within a prescribed range. In phonon - rattle , phonon modes are first obtained for a given structure using the finite-displacement method. The adsorbates are then displaced along the phonon eigenvectors with amplitudes determined from the corresponding eigenvalues, following the procedure proposed by West [25] and implemented in the hiPhive package [26]. In addition, we assign an element-specific cutoff radius to each species. During structure generation, if any pair of atoms approaches closer than the corresponding cutoff distance (defined from the element-specific radii), we perform a short structural relaxation using a simple classical repulsive potential with a characteristic length on the order of the cutoff, thereby removing unphysical atomic overlaps and stabilizing the generated structures.  2.4.  Offspring structure evaluation  After generating offspring structures via evolutionary operations, we select candidates for expensive energy evaluations using BO. In BACCHUS, we construct ACE-based descriptors for BO and use them to evaluate the auction function for all candidates. Finally, we select the top evalN  candidates with the lowest auction function values for exact energy evaluations. In the following, we briefly describe ACE-based descriptors and the BO method used in BACCHUS.  2.4.1.  Atomic cluster expansion (ACE)  To build an accurate surrogate model, it is crucial to encode local atomic environments around each atom with high fidelity. Various descriptors for atomic environments have been proposed, including the Oganov–Valle fingerprint [27], the smooth overlap of atomic positions (SOAP) [28], Behler–Parrinello symmetry functions [29], and moment tensor potentials (MTP) [30]. In CrySPY, the Oganov–Valle fingerprint is used, whereas GOFEE employs an Oganov–Valle-type descriptor augmented with angular terms. Among these approaches, the atomic cluster expansion (ACE) [31] provides a systematic framework for constructing a (in principle) complete basis for ACCEPTED MANUSCRIPTrepresenting local atomic environments. In BACCHUS, we use ACE descriptors generated with ACEpotentials.jl [32]. ACE descriptors are constructed as follows. We define atom-centered basis coefficients by projecting one-particle basis functions onto the neighbor density around atom i :   , = ( , , ),i znlm znlm ij i jj iZ Zψ∈ r  (7)  where ijr  is the relative position vector from atom i  to atom j , and i  denotes the set of neighbors within a cutoff radius. Here, z  and Z  denote the elemental species. The one-particle basis functions ψ  are given by   ˆ( , , ) = ( ) ( ) ,znlm ij i j nl ij lm ij zZ jZ Z R r Yψ δr r  (8)  where =ij ijr r   and ˆ = /ij ij ijrr r , ( )nlR r  is a radial basis function (with a cutoff), and ˆ( )lmY r  is a spherical harmonic. The explicit form of ( )nlR r  is taken from the implementation in ACEpotentials.jl, and this implementation assume ( ) = ( )nl nR r R r , i.e., the radial function is independent of l . Using these one-particle coefficients, product basis functions up to v -body order are constructed as   , ,=1= ,vi i z n l mt t t tt∏znlm   (9)  where =1( ) = ( , , , )vt t t t tz n l mz,n,l,m  denotes a multi-index. From these products, rotationally invariant ACE features are obtained as   = ,ii   (10)  where  i  is the vector whose components are  ,i znlm , and   is a sparse linear operator that performs the contraction (equivalently, integration over the rotation group (3)O ) to ensure rotational invariance. In practice, the matrix elements of   can be expressed in terms of Wigner 3 j  symbols. The resulting vector i  constitutes the ACE descriptor for atom i  and encodes v -body geometric correlations in a symmetry-consistent manner.  2.4.2.  Dimensionality reduction and descriptors  The ACE descriptor i  for each atom i  in a structure is computed as described in Sec. 2.4.1. However, the dimensionality of i  can be quite high, which increases the computational cost of GPR. To mitigate this issue, we first average the ACE descriptors over atoms of the same element type ( z ) within each predefined domain (Fig. 7). Here, a domain refers to one of the following groups of atoms: (i) adsorbate atoms, (ii) surface (substrate) atoms, or (iii) bulk atoms. Moreover, we add the composition ratios of each element type within the adsorbate to the descriptor vector. Then, we apply principal component analysis (PCA), which is implemented in the Scikit-learn package [33], to reduce the dimensionality of the resulting descriptor vector. The length of the final descriptor vector can be determined by specifying the desired cumulative contribution ratio of PCA.  2.4.3.  Bayesian optimization (BO) ACCEPTED MANUSCRIPT In BACCHUS, we employ BO using PYSBO library [34]. The PYSBO library provides a random feature map-based GPR implementation, which is enabled to avoid the computationally expensive training process. This method can be described as bellow. Given a training set =1= {( , )}Mi i iD yx , where ix  denotes the descriptor of structure i  (e.g., an ACE-based representation) and iy  is its hull distance, GPR assumes that the training targets T1= ( , , )My yy   follow a multivariate Gaussian distribution,   ( | ) = ( ).p D Σy ,μ  (11)  When the Bayesian linear regression model is used, y  can be expressed as   T( ) = ( ) ,y φ ε+x w x  (12)  where w is the coefficient vector, ( )φ x  is the feature vector of x , and ε  is noise generated according to a Gaussian distribution (0, )σ . Let M ×Φ ∈   be the design matrix whose i -th row is T( )iφ x , and let M∈y  be the vector of training targets. Then, under the Bayesian linear regression model with Gaussian noise variance 2σ , the mean vector μ  and covariance matrix Σ  of the joint Gaussian ( | ) = ( , )p D y μ Σ  are given by   ( ) 1T 2= ,σ−+ I yμ Φ Φ Φ  (13)  ( ) 12 T 2= ,σ σ−+ IΣ Φ Φ  (14)  where I  is the identity matrix and 1= ( ( ), , ( ))Mφ φx xΦ . From this method, the Gaussian kernel ( , ')k x x  can be approximated as   22|| ' ||( , ') = exp2kη −−  x xx x   T( ) ( '),φ φ x x  (15)  T, ,1 1( ) = ( ( / ), , ( / )) ,z zβ βφ η ηw w x x x  (16)  where T, ( ) = 2 cos( )z β β+w x w x , and d -dimensional random vectors iw  is generated from the probability /2 22 exp( || || /2)dπ − − w . β  is uniformly sampled from [0,2 )π . Using the above GPR model, BO selects the next 1M + -th data vaule 1Mx +  by the aquistion function. In this study, we use expected improvement (EI) as the acquisition function. The EI expected value ( ) of how much maxy  updates when x  is observed:   bestEI( ) = [max(0, ( ))],y y−x x   c= ( )[ ( ) ( ( )) ( ( ))],x t x F t x f t xσ +  (17)  where c best c( ) = ( ( ) ) /t x x yμ σ− . Here, f  is the probability density of (0,1) , and cμ  and cσ  can be computed from the mean and standard deviation of the posterior distribution:   ( | ) = ( | )p D p Dy w   1 121= ( = , = ).A Aμσ− −Φ Σ y  (18)  Here, the matrix A  is defined as T 2= /A IσΦΦ + .   ACCEPTED MANUSCRIPT2.5.  Exact energy evaluation  Following the procedures described in Sec. 2.4, we obtain candidate structures that require high-fidelity energy evaluations. In principle, these exact evaluations are performed using FP calculations based on density functional theory (DFT). However, although BO optimizes the number of structures that should be evaluated, FP calculations remain computationally expensive when applied to a large number of candidates. To reduce the number of costly FP evaluations, we employ a universal neural network potential (UNN) to pre-relax and screen structures. In this work, we use the pretrained MACE [35] foundation potential medium-mpa-0 for this purpose. By combining UNN-based pre-evaluation with selective FP refinement, we substantially lower the overall computational cost while retaining FP-level accuracy for the final reported energetics.  3.  Results  In this section, we present three case studies to which the proposed BACCHUS framework is applied.  The parameters used for the BACCHUS framework in each case study are summarized in Appendix B.  The details of the computational settings for structural relaxation and energy evaluation are provided in Appendix C. Moreover, in Appendix D, we validate the accuracy of the MACE potential used for pre-evaluation by comparing its predictions with FP. First, we examine oxidation on FCC Pt(111) surface in Sec. 3.1 and discuss the effectiveness of the BO-driven CSP scheme. Next, in Sec. 3.2, we investigate surface reconstruction and oxidation on diamond Ge(100) surface. Finally, in Sec. 3.3, we explore Au–Sn alloys on FCC Au(111) surface.  3.1.  Pt–O surface structures on FCC Pt(111)  We first apply the proposed BACCHUS framework to the oxidation of the FCC Pt(111) surface, which is of central importance for heterogeneous catalysis and electrocatalysis. To validate the effectiveness of the BO-driven CSP scheme, we also perform a reference EA-based search without BO acceleration, in which candidates for exact energy evaluations are selected randomly. For this benchmark, we use a supercell size of S = 16n  and a total coverage of = 1.0θ  monolayer (ML), while varying the oxygen composition within the range of 0.25–0.90. In other words, we search surface structures spanning compositions from Pt 12 O 4  to Pt 2 O 14 . At each EA generation, the algorithm produces approximately 3000 offspring structures, from which we select 55  structures for exact (high-fidelity) energy evaluations. Each structural search is repeated three times starting from different randomly generated initial populations. Figure 8(a) shows the evolution of the formation-energy convex hull obtained by the BO-driven EA search. Here, the convex-hull area is defined as the area enclosed by the convex hull and the two end-member points (FCC Pt and O 2 ). The BO-driven EA search progressively discovers lower-energy structures across compositions compared with the initial random pool through iterative generations, and the convex hull obtained after the 12th generation exhibits a substantially deeper hull than the initial one. Figure 8(b) compares the growth of the convex-hull area between the EA search without BO (random selection) and the BO-driven search. From this figure, the BO-driven search increases the convex-hull area more rapidly than the random-ACCEPTED MANUSCRIPTselection baseline, indicating that incorporating BO based on ACE descriptors (as mentioned in Sec. 2.4) can substantially accelerate the EA-based search. From Fig. 8(a), we further confirm the efficiency of the BO-driven search and observe that Pt–O structures with higher oxygen concentrations tend to be more stable. Motivated by this trend, we next focus on searching stable structures in an oxygen-rich regime (composition range of 0.50–0.75). Here, we set S = 8n , which is smaller than in the benchmark above. To compensate for the smaller supercell, we increase the total coverage from = 1.0θ  to = 1.5θ . That is, we span compositions from Pt 6 O 6  to Pt 3 O 9 . Under these conditions, we search for low-energy surface structures and compare them with structures reported in experimental studies. Figures 9(a)–(d) show the lowest-energy structure at each composition obtained using the BACCHUS framework. For Pt 6 O 6 , Pt 5 O 7 , and Pt 4 O 8 , the most stable configurations feature square-planar PtO 4  units, in which four oxygen atoms are coordinated around a Pt atom, and neighboring units share edges. Similar square-planar PtO 4  motifs also appear in bulk Pt 3 O 4  [36]. In particular, the PtO 4  units in the Pt 6 O 6  and Pt 5 O 7  surfaces exhibit a row-ordered arrangement on top of the underlying Pt adatoms, whereas those in Pt 4 O 8  form a checkerboard-like pattern. At higher oxygen contents, the structure evolves into an edge-sharing network of PtO 6  octahedra, as found for Pt 3 O 9 . Although the nearest-neighbor and second nearest-neighbor Pt–Pt distances in our PtO 6 -octahedral motifs are about 2.8 and 3.2 Å, respectively, similar PtO 6-octahedral motifs are also present in bulk α -PtO 2  (CdI 2 -type structure), where the nearest-neighbor Pt–Pt distance is 3.16 Å [36].  Although Pt 3 O 9  is slightly less stable than the structures on the convex hull when the bulk-oxide chemical-potential limit is considered (see Fig. 14 and Appendix D), its hull distance is small. Therefore, Pt 3 O 9  may still be regarded as a possible metastable structure in the early stage of oxide-film growth.  In experimental studies using X-ray photoelectron spectroscopy (XPS) and scanning tunneling microscopy (STM) [37–40], Pt(111) oxidation has been reported to yield various ordered surface-oxide motifs, including row-ordered and other commensurate superstructures, depending on the preparation conditions. The row-ordered arrangements of PtO 4  units identified here provide atomistic models for such ordered phases. Moreover, the formation of α -PtO 2 -like trilayers has also been reported experimentally. Therefore, the structures obtained in our calculations appear to be broadly consistent with the experimentally reported motifs, suggesting that the proposed framework captures key structural features of Pt(111) oxidation.  3.2.  Ge and Ge–O surface structures on diamond Ge(100)  We focus on the surface reconstruction of diamond Ge(100), including the effects of oxidation. In this study, we search structures spanning the composition range from Ge 8  to Ge 4 O 4  with supercell size S = 4n . Figures 10(a) and (b) show the lowest and second-lowest-energy reconstruction structures of diamond Ge(100) obtained in this work. Both structures contain Ge dimers. The difference between them lies in the arrangement of the Ge dimers and in the buckling direction of the dimer atoms. The lowest-energy structure (a) exhibits anti-parallel dimer buckling. These degrees of ACCEPTED MANUSCRIPTfreedom of the Ge dimers lead to the energy difference between the two reconstructions, which is estimated to be 110 meV per dimer from our DFT calculations. Previous DFT studies [45,48–50] have proposed two stable dimer arrangements, namely the (4 2)c ×  and (2 2)p ×  reconstructions, which are nearly degenerate in energy. The (2 2)p ×  reconstruction corresponds to our lowest-energy structure. In experimental STM studies [45], both phases have been observed to coexist on the surface. Then, we discuss the effect of oxidation on the Ge(100) surface reconstruction. Previous DFT studies [51–54] have proposed several adsorption configurations of an O atom on a Ge dimer (Fig. 10(c)). In particular, the dimer-bridge site and the back-bond-down site have been reported to be energetically stable. In our study, the back-bond-down site is found to be preferred in the lowest-energy structures at each low-oxidation composition (Ge 7 O 1  and Ge 6 O 2  in Fig. 11(a)–(b)). On the other hand, oxygen adsorption at the back-bond-down site induces a lateral shift of the Ge dimer rows. At higher oxidation levels (Ge 5 O 3  and Ge 4 O 4  in Fig. 11(c)–(d)), more substantial structural changes occur. In Ge 5 O 3 , two O adatoms induce the protrusion of a Ge adatom. In Ge4 O 4 , the Ge dimers vanish upon O adsorption, and the reconstructed Ge 4 O 4  structure becomes condensed into a (1 1)p ×  unit cell. In STM experiments [45,48–50], oxidation has been reported to progressively remove the Ge dimers. In addition, bright protrusions are observed in STM experiments, which have been theoretically attributed to Ge adatom, and a (1 1)p ×  reconstruction emerges at higher oxidation levels. These experimental observations are consistent with our obtained structures: the disappearance of Ge dimers and the appearance of protruding Ge adatoms are captured in the Ge 5 O 3  and Ge 4 O 4  structures (Fig. 11(c)–(d)), and the condensation into a (1 1)p ×  unit cell in Ge 4 O 4  is consistent with the experimentally reported (1 1)p ×  phase [55].  3.3.  Au–Sn alloying on FCC Au (111)  To examine the performance of our framework for surface-alloy systems, we finally investigate Au–Sn alloying on FCC Au(111) using a supercell size of S = 9n  and a total coverage of = 3θ  ML (Fig. 12). At the Au 24 Sn 3  composition, Sn atoms segregate to the surface and form a hexagonal overlayer; the resulting top-layer stoichiometry corresponds to Au 2 Sn. This Au 2 Sn surface-alloy motif is consistent with experimental STM studies [51–54] on Sn/Au(111), which report ordered Au 2 Sn surface-alloy phases and commensurate superstructures.  In addition, as a new insight obtained from the present theoretical approach, we found that, at the Sn-rich composition Au14 Sn 13 , the system no longer preserves the underlying FCC-Au lattice; instead, the optimized structures exhibit a strongly disordered, “degradeda arrangement.  4.  Conclusion  We have developed a global surface-structure search framework, BACCHUS, that simultaneously explores atomic coordinates, in-plane superlattice periodicity, and ACCEPTED MANUSCRIPTstoichiometry/coverage within a unified search space. By combining an evolutionary algorithm equipped with surface-specific variation operators and symmetry-enriched initialization with Bayesian optimization based on ACE descriptors, the method efficiently navigates highly rugged energy landscapes while reducing the number of expensive first-principles evaluations. In the Pt(111) oxidation case study, the BO-driven search was shown to accelerate the growth of the formation-energy convex hull compared with a random-selection baseline, and it revealed oxygen-rich low-energy motifs consistent with experimentally discussed surface-oxide trends. For Ge(100), the framework reproduced the characteristic dimer reconstructions and captured oxidation-induced restructuring, including dimer-row shifts at low oxidation and the disappearance of dimers toward a (1 1)p × -like phase at higher oxidation, in line with reported STM observations. Finally, for Au–Sn surface alloying on Au(111), the BACCHUS framework identified surface segregation and hexagonal ordering of Sn at Au 24 Sn 3  with an Au 2 Sn-like top-layer stoichiometry consistent with experimental reports, while predicting strongly distorted, bulk-alloy-like arrangements at Sn-rich compositions such as Au14 Sn13 .  Because the present framework treats the surface periodicity as a search variable, it can, in principle, systematically explore long-period surface reconstructions. Nevertheless, complex long-period reconstructions, such as the Dimer–Adatom–Stacking-fault structure of the Si(111) 7×7 surface, remain challenging in practice because evolutionary structure searches generally become less efficient as the number of atoms and the system size increase [55].  The proposed framework provides a practical and general route to predicting complex surface and interface structures relevant to heterogeneous catalysis, electronics, and spintronics, where coupled structural and compositional degrees of freedom govern functionality. Future work will extend the approach to broader chemistries, environmental conditions, long-period reconstructions, and fully automated integration with high-throughput first-principles workflows.   Appendix A.  Variation of surface superlattice choices  Here, we illustrate the restriction imposed by fixed user-defined slab supercells using the case of S = 3n  (Fig. 13), where Sn  denotes the superlattice size, i.e., the area ratio of the surface supercell to the primitive surface unit cell. In the HNF representation (Eq. 2), four candidate supercells are generated for S = 3n  before symmetry reduction. After considering the symmetry of the hexagonal surface lattice, these four candidates are reduced to two symmetry-inequivalent supercell choices. For a larger superlattice size, the number of candidates increases further (the number of HNF candidates is given by the sum of the positive divisors of Sn ); for example, 13 candidate supercells are generated for S = 9n  before symmetry reduction. This example shows that many possible surface periodicities must be considered even for moderate superlattice sizes. In a fixed-cell approach, these candidate slab supercells have to be manually prepared and examined, and this burden increases further when different surface compositions are included. Therefore, fixed-cell treatments can substantially restrict the accessible structure-search space and may miss reconstructions requiring non-diagonal or differently shaped surface cells.  ACCEPTED MANUSCRIPT Appendix B.  Parameterization of the BACCHUS framework  In this appendix, we summarize the key parameters used in the BACCHUS framework for the case studies presented in the main text. In Eq. 5, the parameters were set to = 10a  and = 1b . We assigned an element-specific cutoff radius to each species. The cutoff radii were set to 1.0 Å for Ge, 0.7 Å for O, 1.0 Å for Pt, 1.0 Å for Au, and 1.0 Å for Sn. ACE descriptors were generated with a cutoff radius of 5.0 Å, and the maximum body order was set to 4. For the reference chemical potentials used in the formation-energy calculations, Aμ  and Bμ  were evaluated as the total energy per atom of the corresponding reference phase calculated using the same computational settings. Specifically, diamond Ge, diamond Sn, an isolated O 2  molecule, FCC Au, and FCC Pt were used as the reference phases for Ge, Sn, O, Au, and Pt, respectively. For oxygen, the reference chemical potential was taken as one half of the total energy of the isolated O 2  molecule.   Appendix C.  Computational settings  In this appendix, we provide the computational settings used for structural relaxation and energy evaluation. The global structure search and the associated energy/force evaluations were performed using the MACE pretrained model mace-mpa-0. During the search stage, structural relaxations were carried out only for the explored surface region (i.e., the adsorbates and the topmost surface layer), while the remaining slab atoms were kept fixed to reduce the computational cost. After the search, the most stable structures for each composition were re-optimized using DFT as implemented in Quantum ESPRESSO (QE) [56,57]. In these QE calculations, all atoms including the slab layers were fully relaxed. We employed the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) for the exchange–correlation functional [58]. Ultrasoft pseudopotentials (USPPs) [59] from the GBRV (Garrity–Bennett–Rabe–Vanderbilt) library [60] were used. For structural relaxations, the plane-wave cutoff and charge-density cutoff were set to 40 Ry and 400 Ry, respectively, with a 6 6 1× ×  Monkhorst–Pack k -point mesh. For final total-energy evaluations, we used the higher cutoffs (50 Ry and 500 Ry) with a denser 12 12 1× ×  k -point mesh. To obtain reliable formation energies involving oxygen, we applied an overbinding correction to the DFT total energy of O 2  within the GGA-PBE framework [61].  Appendix D.  Validation of MACE potential and discussion on thermodynamic stability.   In this appendix, we validate the accuracy of the MACE potential used for pre-screening by comparing its predictions with DFT calculations. Figure 14 presents the comparison. For the Pt–O and Ge–O systems, we computed the formation energies of the lowest-energy structures at each composition. For the Au–Sn system, we computed the formation energies of the structures ACCEPTED MANUSCRIPTon the formation-energy convex hull. Based on this validation, we confirmed that the MACE potential provides sufficiently accurate relative energies for the pre-screening stage.  In addition, we estimate the oxygen chemical-potential limits from bulk oxide phases. Among the considered bulk oxide structures from the Materials Project database [62] , PtO 2  (hydrophilite-type structure) and GeO 2  (rutile-type structure) are found to have the most negative formation energies for the Pt–O and Ge–O systems, respectively. The corresponding oxide-formation limits are evaluated as   b bP P2O/P 2= ,2ulk ulktO ttOE Eμ− (19)  b bG G2O/G 2= ,2ulk ulkeO eeOE Eμ− (20)  where bP 2ulktOE  and bG 2ulkeOE  are the total energies of bulk PtO 2  and GeO 2  of formula unit, respectively, and bPulktE  and bGulkeE  are the total energies of bulk Pt and Ge, respectively. These chemical potentials represent the upper bounds of the oxygen chemical potential above which bulk oxide formation becomes thermodynamically favorable. On the other hand, these oxide-derived chemical potentials correspond to the bulk-oxide limit, where the oxide layer grown on the substrate is sufficiently thick and can be regarded as a bulk crystalline phase. In the initial stage of oxide growth, the oxide layer is not fully developed, and surface, interface, and lattice-mismatch effects can contribute significantly to the total energy. Thus, the effective oxide energy may be written as f b s /o o=ilm ulk urf interfacexide xideE E E+ Δ . If the excess contribution s /urf interfaceEΔ  is positive, as is commonly expected for thin or incipient oxide layers, the corresponding oxide-formation chemical potential becomes higher than that estimated from the bulk oxide. Therefore, the chemical-potential limits derived from bulk PtO 2  and GeO 2  should be interpreted as bulk-oxide limits, while the actual formation threshold in the early stage of oxide growth may be shifted to higher oxygen chemical potentials. In Fig. 14, we also show the chemical-potential limits estimated from the corresponding bulk oxides. The results indicate that the oxide structures obtained by the present framework are stable with respect to these bulk-oxide limits, except for Pt 3 O 9 . Even for Pt 3 O 9 , however, the hull distance is small. Considering that surface, interface, and lattice-mismatch contributions can be significant in the initial stage of oxide-film growth, Pt 3 O 9  may still exist as a metastable structure during the early oxidation process.   Acknowledgments  This work was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (Grant No. JP23K13537). MO acknowledge financial supports from MEXT based on the Promotion of Development of a Joint Usage/Research System Project: Coalition of Universities for Research Excellence Program (CURE) (Grant No. JPMXP1323015474). The computation in this work has been done using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo (ISSPkyodo-SC-2025-Ca-002X, 2025-D-000X)   ACCEPTED MANUSCRIPTDisclosure statement No potential conflict of interest was reported by the author(s).     Reference  [1] Wenpei Gao, Zachary D. Hood, and Miaofang Chi. Interfaces in heterogeneous catalysts: Advancing mechanistic understanding through atomic-scale measurements. Accounts of Chemical Research, 50(4):787–795, 2017.  [2] Chenlu Xie, Zhiqiang Niu, Dohyung Kim, Mufan Li, and Peidong Yang.  Surface and interface control in nanoparticle catalysis. Chemical Reviews, 120(2):1184–1249, 2020.  [3] Fumiaki Kuroda, Masaya Ibe, and Minoru Otani. Screening multicomponent alloys using evidence theory and high-throughput dft calculations for c2-selective eco2rr. The Journal of Physical Chemistry C, 128:17302–17312, 2024.  [4] Igor Zutic, Jaroslav Fabian, and S. Das Sarma. Spintronics: Fundamentals and applications. Reviews of Modern Physics, 76:323–410, 2004.  [5] Xiong Zhou, Qian Shen, Yongfeng Wang, Yafei Dai, Yongjun Chen, and Kai Wu. Surface and interfacial sciences for future technologies. National Science Review, 11(9):nwae272, 2024.  [6] Artem R. Oganov and Colin W. Glass. Crystal structure prediction using ab initio evolutionary techniques: Principles and applications. The Journal of Chemical Physics, 124(24):244704, 2006.  [7] Scott M. Woodley and Richard Catlow. Crystal structure prediction from first principles.Nature Materials, 7(12):937–946, December 2008.  [8] Artem R. Oganov, Chris J. Pickard, Qiang Zhu, and Richard J. Needs. Structure prediction drives materials discovery. Nature Reviews Materials, 4(5):331–348, May 2019.  [9]  Yanchao Wang, Jian Lv, Li Zhu, and Yanming Ma.   Crystal structure prediction via particle-swarm optimization. Physical Review B, 82(9):094116, 2010.  [10] Yanchao Wang, Jian Lv, Li Zhu, and Yanming Ma. Calypso: A method for crystal structure prediction.  Computer Physics Communications, 183(10):2063–2070, 2012.  [11] Shaohua Lu, Yanchao Wang, Hanyu Liu, Mao-sheng Miao, and Yanming Ma.   Self-assembled ultrathin nanotubes on diamond (100) surface. Nature Communications, 5:3666, 2014.  [12] Colin W. Glass, Artem R. Oganov, and Nikolaus Hansen. Uspex—evolutionary crystal structure prediction. Computer Physics Communications, 175(11-12):713–720, 2006.  [13] Andriy O. Lyakhov, Artem R. Oganov, Harold T. Stokes, and Qiang Zhu. New develop-ments ACCEPTED MANUSCRIPTin evolutionary structure prediction algorithm uspex. Computer Physics Communications, 184(4):1172–1182, 2013.  [14] Artem R. Oganov, Andriy O. Lyakhov, and Mario Valle. How evolutionary crystal structure prediction works—and why. Accounts of Chemical Research, 44(3):227–237, 2011.  [15] Qiang Zhu, Li Li, Artem R. Oganov, and Philip B. Allen. Evolutionary method for predict-ing surface reconstructions with variable stoichiometry. Physical Review B, 87(19):195317, 2013.  [16] Tomoki Yamashita, Nobuya Sato, Hiori Kino, Takashi Miyake, Koji Tsuda, and Tamio Oguchi.   Crystal structure prediction accelerated by bayesian optimization.   PhysicalReview Materials, 2(1):013803, 2018.  [17] Tomoki Yamashita, Shinichi Kanehira, Nobuya Sato, Hiori Kino, Kei Terayama, Hikaru Sawahata, Takumi Sato, Futoshi Utsuno, Koji Tsuda, Takashi Miyake, and Tamio Oguchi. Cryspy: a crystal structure prediction tool accelerated by machine learning. Science and Technology of Advanced Materials: Methods, 1(1):87–97, 2021.  [18] Malthe K. Bisbo and Bjørk Hammer. Efficient global structure optimization with a machine-learned surrogate model. Physical Review Letters, 124:086102, 2020.  [19] Malthe K. Bisbo and Bjørk Hammer. Global optimization of atomic structure enhancedby machine learning. Physical Review B, 105:245404, 2022.  [20] A. Santoro and A. D. Mighell. Properties of crystal lattices: The derivative lattices and their determination. Acta Crystallographica Section A, 28:284–287, 1972.  [21] A. Santoro and A. D. Mighell. Coincidence-site lattices. Acta Crystallographica SectionA, 29:169–175, 1973.  [22] Gus L. W. Hart and Rodney W. Forcade. Algorithm for generating derivative structures.Physical Review B, 77:224115, 2008.  [23] Atsushi Togo, Kohei Shinohara, and Isao Tanaka. Spglib: a software library for crystal symmetry search. Science and Technology of Advanced Materials: Methods, 4(1):2384822– 2384836, 2024.  [24] Fumiaki Kuroda, Satoshi Hagiwara, and Minoru Otani. Structural changes in the lithium cobalt dioxide electrode: A combined approach with cluster expansion and bayesian op-timization. Physical Review Materials, 7(11):115402, 2023.  [25] D. West and S. K. Estreicher. First-principles calculations of vibrational lifetimes and decay channels: Hydrogen-related modes in si. Physical Review Letters, 96(11):115504, 2006.  [26] Fredrik Eriksson, Erik Fransson, and Paul Erhart. The hiphive package for the extraction of ACCEPTED MANUSCRIPThigh-order force constants by machine learning. Advanced Theory and Simulations, 2(5):1800184,  2019.  [27] Artem R. Oganov and Mario Valle. How to quantify energy landscapes of solids. The Journal of Chemical Physics, 130(10):104504, 2009.  [28] Albert P. Bartók, Risi Kondor, and Gábor Csányi. On representing chemical environments. Physical Review B, 87(18):184115, 2013.  [29] Jörg Behler and Michele Parrinello. Generalized neural-network representation of high-dimensional potential-energy surfaces. Physical Review Letters, 98(14):146401, 2007.     [30] Alexander V. Shapeev. Moment tensor potentials: A class of systematically improvable interatomic potentials. Multiscale Modeling & Simulation, 14(3):1153–1173, 2016.  [31] Ralf Drautz. Atomic cluster expansion for accurate and transferable interatomic potentials. Physical Review B, 99(1):014104, January 2019.  [32] William C. Witt, Cas van der Oord, Elena Gelžinytė, Teemu Järvinen, Andres Ross,James  P.  Darby,  Cheuk  Hin  Ho,  William  J.  Baldwin,  Matthias  Sachs,  James  Kermode, Noam Bernstein, Gábor Csányi, and Christoph Ortner. Acepotentials.jl: A julia implementation of the atomic cluster expansion. The Journal of Chemical Physics, 159(16):164101, 2023.  [33] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher,Matthieu Perrot, and Édouard Duchesnay.  Scikit-learn: Machine learning in python.Journal of Machine Learning Research, 12:2825–2830, 2011.  [34] Yuichi Motoyama, Ryo Tamura, Kazuyoshi Yoshimi, Kei Terayama, Tsuyoshi Ueno, and Koji Tsuda. Bayesian optimization package: Physbo. Computer Physics Communications, 278:108405, 2022.  [35] Ilyes Batatia, David P Kovacs, Gregor Simm, Christoph Ortner, and Gabor Csanyi. Mace: Higher order equivariant message passing neural networks for fast and accurate force fields.In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Ad-vances in Neural Information Processing Systems, volume 35, pages 11423–11436. Curran Associates, Inc., 2022.  [36] Ricardo K. Nomiyama, Maurício J. Piotrowski, and Juarez L. F. da Silva. Bulk structures of pto and pto2 from density functional calculations. Physical Review B, 84(10):100101, 2011.  [37] Sunil P. Devarajan, Jose A. Jr. Hinojosa, and Jason F. Weaver. Stm study of high-coverage structures of atomic oxygen on pt(1 1 1): p(2 × 1) and pt oxide chain structures. Surface Science, 602(19):3116–3124, 2008.  ACCEPTED MANUSCRIPT[38] C. Ellinger, A. Stierle, I. K. Robinson, A. Nefedov, and H. Dosch. Atmospheric pressure oxidation of pt(111). Journal of Physics: Condensed Matter, 20(18):184013, 2008.  [39] Matthijs A. van Spronsen, Joost W. M. Frenken, and Irene M. N. Groot. Observing theoxidation of platinum. Nature Communications, 8(1):429, 2017.  [40] Daniel Miller, Hernan Sanchez Casalongue, Hendrik Bluhm, Hirohito Ogasawara, An-ders Nilsson, and Sarp Kaya. Different reactivity of the various platinum oxides and chemisorbed oxygen in co oxidation on pt(111). Journal of the American Chemical Society, 136(17):6340–6347, 2014.  [41] M. Needels, M. C. Payne, and J. D. Joannopoulos. High-order reconstructions of the ge(100) surface. Physical Review B, 38(8):5543–5546, 1988.  [42] L. Spiess, A. J. Freeman, and P. Soukiassian.  Ge(100) 2×1 and c(4×2) surface recon-structions studied by ab initio total-energy molecular-force calculations. Physical Review B, 50(4):2249, 1994.  [43] K. Noatschk, E. V. S. Hofmann, J. Dabrowski, N. J. Curson, T. Schroeder, W. M. Klesse, and G. Seibold.  Ge(001) surface reconstruction with sn impurities.  Surface Science, 713:121912, 2021.  [44] Claudia Fleischmann, Koen Schouteden, Clement Merckling, Sonja Sioncke, Marc Meuris,C. Van Haesendonck, K. Temst, and A. Vantomme. Adsorption of O2 on Ge(100): Atomic geometry and site-specific electronic structure. The Journal of Physical Chemistry C, 116(18):9925–9929, April 2012.  [45] Tyler J. Grassman, Sarah R. Bishop, and Andrew C. Kummel. An atomic view of fermi level pinning of ge(100) by o2. Surface Science, 602(14):2373–2381, 2008.  [46] Ghazanfar  Shah,  Marian  W.  Radny,  and  Phillip  V.  Smith.    Initial  stages  of  oxy-gen chemisorption on the Ge(001) surface. The Journal of Physical Chemistry C, 118(29):15795–15803, July 2014.  [47] Kai Liu, In Won Yeu, Cheol Seong Hwang, and Jung Hae Choi. Initial oxidation and surface stability diagram of ge(100) as a function of the temperature and oxygen partial pressure through ab initio thermodynamics. Physica Scripta, 95(2):025701, 2020.  [48] T. Fukuda and T. Ogino.  Initial oxygen reaction on Ge(100) 2 × 1 surfaces.  Physical Review B, 56(20):13190–13193, 1997.  [49] T. Fukuda. Oxygen-induced missing dimer row formation on the Ge(100) surface. Surface Science, 417(2-3):L149–L153, 1998. Surface Science Letters.  [50] T. Fukuda and T. Ogino.  Initial oxidation stage of the Ge(100) 2 × 1 surface studied by scanning tunneling microscopy and ultra-violet photoelectron spectroscopy. Applied ACCEPTED MANUSCRIPTSurface Science, 130-132:165–169, 1998.  [51]  M. Maniraj, D. Jungkenn, W. Shi, S. Emmerich, L. Lyu, J. Kollamana, Z. Wei, B. Yan, M. Cinchetti, S. Mathias, B. Stadtmüller, and M. Aeschlimann. Structure and electronic properties of the (√3 × √3)r30◦  SnAu2/Au(111) surface alloy. Phys. Rev. B, 98:205419,Nov 2018.  [52] Pampa Sadhukhan, Dhanshree Pandey, Vipin Kumar Singh, Shuvam Sarkar, Abhishek Rai, Kuntala Bhattacharya, Aparna Chakrabarti, and Sudipta Roy Barman. Electronic structure and morphology of thin surface alloy layers formed by deposition of sn on au(1 1 1). Applied Surface Science, 506:144606, 2020.  [53] Jalil Shah, W. Wang, Hafiz M. Sohail, and R. I. G. Uhrberg. Atomic and electronic structures of the au2sn surface alloy on au(111). Physical Review B, 104(12):125408, 2021.  [54] Julian A. Hochhaus, Stefanie Hilgers, Marie Schmitz, Lukas Kesper, Ulf Berges, and Carsten Westphal.  Structural analysis of sn on au(111) at low coverages: Towards theau2sn surface alloy with alternating fcc and hcp domains.  Scientific Reports, 15:7953,2025.  [55] Andriy O. Lyakhov, Artem R. Oganov, and Mario Valle. How to predict very large and complex crystal structures. Computer Physics Communications, 181(9):1623–1632, 2010.  [56] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, et al. QUANTUM ESPRESSO: a modular and open-source  software project for quantum simulations of materials. Journal of Physics: Condensed Matter, 21(39):395502, 2009.  [57] Paolo Giannozzi, Oliviero Andreussi, Tobias Brumme, Olivier Bunau, Marco Buon-giorno Nardelli, Matteo Calandra, et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. Journal of Physics: Condensed Matter, 29(46):465901, 2017.  [58] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approxi-mation made simple. Physical Review Letters, 77(18):3865–3868, 1996.  [59] David Vanderbilt. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Physical Review B, 41(11):7892–7895, 1990.  [60]  Kevin F. Garrity, Joseph W. Bennett, Karin M. Rabe, and David Vanderbilt.   Pseu-dopotentials for high-throughput DFT calculations. Computational Materials Science, 81:446–452, 2014.  [61] Lei Wang, Thomas Maxisch, and Gerbrand Ceder. Oxidation energies of transition metal oxides within the GGA+U framework. Physical Review B, 73(19):195107, 2006.  [62] Matthew K. Horton, Patrick Huck, Ruo Xi Yang, Jason M. Munro, Shyam Dwaraknath,Alex M. Ganose, Ryan S. Kingsbury, Mingjian Wen, Jimmy X. Shen, Tyler S. Mathis, Aaron D. ACCEPTED MANUSCRIPTKaplan, Karlo Berket, Janosh Riebesell, Janine George, Andrew S. Rosen, EvanC. Spotte-Smith, Matthew J. McDermott, Orion A. Cohen, Alex Dunn, Matthew C. Kuner, Gian-Marco Rignanese, Guido Petretto, David Waroquiers, Sinead M. Griffin, Jeffrey B. Neaton, Daryl C. Chrzan, Mark Asta, Geoffroy Hautier, Shreyas Cholia, Gerbrand Ceder, Shyue Ping Ong, Anubhav Jain, and Kristin A. Persson. Accelerated data-driven materials science with the materials project. Nature Materials, 24(10):1522–1532, 2025. ACCEPTED MANUSCRIPTFigure  1. Schematic illustration of surface-structure search with variable periodicity and composition.   Figure  2. Schematic BACCHUS workflow  Figure  3. Two symmetry-equivalent supercells on FCC (111) surface  Figure  4. Symmetry-enriched site generation on a supercell surface. (a) three fold rotational symmetry, (b) mirror symmetry.  Figure  5. Schematic view of the formation-energy convex hull and the parent-selection procedure.  Figure  6. Evolutionary operations implemented in BACCHUS. Different colors represent different element types.  Figure  7. Schematic illustration of the ACE descriptor matrices. Light blue, yellow, and gray spheres represent adsorbate, surface, and bulk atoms, respectively.  Figure  8.  (a) Evolution of the formation-energy convex hull obtained by the BO-driven EA search. (b) Comparison of the growth of the convex-hull area between the EA search without BO (random selection) and the BO-driven search. Solid lines indicate the median, while dashed lines represent the minimum and maximum values.   Figure  9. The lowest-energy structure at each composition for Pt oxidation on FCC Pt (111). Blue spheres indicate Pt adatoms on the surface, light-blue spheres indicate the topmost Pt layer in the slab, white spheres indicate the remaining Pt atoms in the slab, and red spheres indicate O atoms.  Figure  10. The lowest (a) and second-lowest (b) reconstruction structures of diamond Ge(100) obtained by this work. Black spheres indicate Ge adatoms on the surface, light-blue spheres indicate the topmost Ge layer in the slab, and white spheres indicate the remaining Ge atoms in the slab. Panel (c) schematically illustrates the adsorption site of an O atom on the widely discussed Ge dimer   Figure  11. The lowest-energy structure at each composition for Ge oxidation on diamond Ge (100). Black spheres indicate Ge adatoms on the surface, light-blue spheres indicate the topmost Ge layer in the slab, white spheres indicate the remaining Ge atoms in the slab, and red spheres indicate O atoms.   Figure  12. The lowest-energy structure at Au 24 Sn 3  and Au 14 Sn 13  on FCC Au (111). Dark gold spheres indicate topmost Au adatoms on the surface, gold spheres indicate the Au surface layer, white spheres indicate the remaining Au atoms in the slab, and silver spheres indicate Sn atoms.  Figure  13.  Schematic illustration of superlattice choices for a hexagonal lattice with S = 3n . ACCEPTED MANUSCRIPTBlue solid lines indicate the supercell lattices generated by the HNF representation with symmetry reduction.   Figure  14. Formaiton energy comparison with DFT (QE) and MACE (medium-mpa-0).   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT   ACCEPTED MANUSCRIPT