# Fileset

[lin-et-al-2025-determination-of-stable-proton-configurations-by-black-box-optimization-using-an-ising-machine.pdf](https://mdr.nims.go.jp/filesets/8169036a-e7e2-499b-878f-955045565992/download)

## Creator

Jianbo Lin, [Tomofumi Tada](https://orcid.org/0000-0003-3093-3779), Ai Koizumi, [Masato Sumita](https://orcid.org/0000-0002-3506-1028), [Koji Tsuda](https://orcid.org/0000-0002-4288-1606), [Ryo Tamura](https://orcid.org/0000-0002-0349-358X)

## Rights

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

## Other metadata

[Determination of Stable Proton Configurations by Black-Box Optimization Using an Ising Machine](https://mdr.nims.go.jp/datasets/5c36e926-074d-439f-b56f-2fcd8c20e4b7)

## Fulltext

Determination of Stable Proton Configurations by Black-Box Optimization Using an Ising MachineDetermination of Stable Proton Configurations by Black-BoxOptimization Using an Ising MachineJianbo Lin,* Tomofumi Tada,* Ai Koizumi, Masato Sumita, Koji Tsuda, and Ryo Tamura*Cite This: J. Phys. Chem. C 2025, 129, 2332−2340 Read OnlineACCESS Metrics & More Article Recommendations *sı Supporting InformationABSTRACT: Stable proton configurations in solid-state materials are aprerequisite for the theoretical microscopic investigation of solid-stateproton-conductive materials. However, a large number of initial atomisticconfigurations should be considered to find stable proton configurations, andrelaxation calculations using the density functional theory approach arerequired for each initial configuration. Consequently, the determination ofstable configurations is a difficult and time-consuming task. Furthermore,when the size of the simulation cells or the number of doped atoms increases,the number of initial configurations leads to a combinatorial explosion,rendering the computation infeasible. In this study, black-box optimizationwas combined with an Ising machine and density functional calculations toperform an efficient search for stable proton configurations. Scandium-dopedbarium zirconate, a typical high-proton conductive oxide, was selected as the model system. The Ising machine was able to rapidlyselect the initial atomistic configuration, ultimately leading to stable proton configurations after subsequent relaxation calculations.This optimization strategy should be able to solve various issues related to configuration optimization in solid-state materials, therebypromoting novel scientific discoveries.■ INTRODUCTIONProtonic ceramic fuel cells (PCFCs) are expected to exhibithigh power generation efficiencies. As previously reported,doping is a key process in achieving high proton conductivitiesin PCFC oxide electrolytes.1−6 More specifically, doping withlow-valence elements creates oxygen vacancies that maintainelectrical neutrality, and these vacancies can subsequently befilled with hydroxy groups to introduce protons into thesystem. Previously, the relationship between the conductivityand the doping level has been studied experimentally.7,8However, to understand the mechanism of proton conductivityfrom a microscopic viewpoint, atomistic simulations based onmolecular dynamics are essential,9,10 and stable atomisticconfigurations are required to commence these simulations. Atpresent, determination of the stable atomistic configurations ofdoped atoms and protons is a difficult and time-consumingtask since the initial configurations of the doping atoms andprotons can change depending on the size of the simulationcell and the number of doped atoms. In addition, relaxationcalculations using the density functional theory (DFT)approach are required for each initial configuration, therebyrendering it impossible to perform relaxation calculations forall possible initial configurations and to find stable structuresusing such approaches.Using machine learning to optimize atomistic configurations,Ju et al.11 demonstrated that the configurations of silicon andgermanium atoms in a lattice could be rapidly explored bymeans of Bayesian optimization to maximize or minimizephonon transport. Bayesian optimization is a method foroptimizing black-box functions,12 and is characterized by theuse of a Gaussian process as a surrogate model.13,14 Followingthe work of Ju et al., many structure search problems aimed atatomistic configurations have been solved using Bayesianoptimization.15−19 However, Bayesian optimization has aproblem in that when the number of candidates causes acombinatorial explosion, the time required for the selection ofpromising candidates increases exponentially because acquis-ition function calculations are required for all candidates. Toovercome this limitation, Kitai et al. proposed a black-boxoptimization (BBO) algorithm using a quantum annealerknown as the factorization machines with quantum annealing(FMQA) algorithm.20 In the FMQA algorithm, the factoriza-tion machine (FM)21 model is adopted as a surrogate model.Because the FM model is expressed by the Ising model, thelower-energy states of the FM can be effectively obtained usingquantum annealers and Ising machines.22 Indeed, thisalgorithm has been widely used to solve BBO problems inmaterials science and chemistry.23−30 For example, the FMQAReceived: October 18, 2024Revised: January 4, 2025Accepted: January 21, 2025Published: January 28, 2025Articlepubs.acs.org/JPCC© 2025 The Authors. Published byAmerican Chemical Society2332https://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−2340This article is licensed under CC-BY 4.0Downloaded via NATL INST FOR MATLS SCIENCE (NIMS) on May 22, 2025 at 00:43:55 (UTC).See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.https://pubs.acs.org/action/doSearch?field1=Contrib&text1="Jianbo+Lin"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Tomofumi+Tada"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Ai+Koizumi"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Masato+Sumita"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Koji+Tsuda"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Ryo+Tamura"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://pubs.acs.org/action/showCitFormats?doi=10.1021/acs.jpcc.4c07104&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?goto=articleMetrics&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?goto=recommendations&?ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?goto=supporting-info&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=tgr1&ref=pdfhttps://pubs.acs.org/toc/jpccck/129/5?ref=pdfhttps://pubs.acs.org/toc/jpccck/129/5?ref=pdfhttps://pubs.acs.org/toc/jpccck/129/5?ref=pdfhttps://pubs.acs.org/toc/jpccck/129/5?ref=pdfpubs.acs.org/JPCC?ref=pdfhttps://pubs.acs.org?ref=pdfhttps://pubs.acs.org?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://pubs.acs.org/JPCC?ref=pdfhttps://pubs.acs.org/JPCC?ref=pdfhttps://acsopenscience.org/researchers/open-access/https://creativecommons.org/licenses/by/4.0/https://creativecommons.org/licenses/by/4.0/https://creativecommons.org/licenses/by/4.0/algorithm was adopted to optimize the atomic configurationsof the magnesium and germanium ions in a magnetic tunneljunction structure,31 and its efficiency was confirmed for smallproblems. Several examples have been reported that theFMQA algorithm tends to achieve better optimizationperformance than other optimization methods such asBayesian optimization,20 genetic algorithm,26,32 and particleswarm optimization26 in previous studies. Furthermore, itsapplication in crystal structure prediction has also beenconsidered.33,34In this study, the FMQA algorithm combined with theVienna Ab initio Simulation Package (VASP) is employed toeffectively search for stable proton configurations in PCFCs(Figure 1). In each optimization cycle, the FM model is trainedto predict the total energy after the relaxation calculation fromthe initial proton configuration. By solving the FM model usingan Ising machine, an initial proton configuration that leads tostable proton configurations after the relaxation calculation issuggested. For the suggested configuration, the total energy iscalculated by relaxation calculations using VASP, and theamount of training data for the FM model is increased.Subsequently, in each cycle, a relaxation calculation must beperformed using VASP to acquire the total energy with therequired calculation accuracy. The optimization task istherefore to determine the stable proton configurations usingas few cycles as possible. For the purpose of this study,scandium-doped barium zirconate,35 which is a representativeelectrolyte candidate for PCFCs, is employed as a modelsystem to perform the stable proton configuration searches bycombining FMQA and VASP. In addition, a graphicsprocessing unit-based Ising machine is used to solve the FMmodel.■ METHODSTarget Material. As a model system for protonconfiguration optimization, scandium-doped barium zirconatewas employed according to the crystal structure shown inFigure 2a. The mother compound was the BaZrO3 perovskite,and fast proton conduction has previously been observed in itsdoped materials,4,5,37 wherein a high proton conductivity wasachieved by substituting zirconium with scandium.35 Underthese conditions, oxygen vacancies are generated, and hydroxygroups can be introduced into the system in the presence ofwater vapor to compensate for these vacancies. In the currentstudy, eight initial proton positions are considered at eachoxygen atom, as determined by DFT structural optimizationsfor several configurations, including those of the protons. Eachproton position was directed to the second nearest oxygenatom to form a hydrogen bond (see Figure 2b). Thus, theseproton positions can be used as the initial protonconfigurations for the relaxation calculations. A 2 × 2 × 2supercell model was adopted to search for a stableconfiguration. Based on the number of scandium ions (x) inthe supercell, the composition of the model was expressed asBa8Zr8−xScxO24Hx, corresponding to the condition that alloxygen vacancies are compensated for by hydroxy groups;therefore, the number of introduced protons is the same as thenumber of scandium ions. When the doped positions of thescandium ions are fixed, the number of proton configurations is24Cx × 8x, without considering symmetry. Thus, when thevalue of x increases, many proton configurations must beconsidered. Of course, when symmetry is considered, thenumber of independent initial proton configurations isreduced; however, a combinatorial explosion of the initialconfigurations occurs when larger systems are considered.Structural Relaxation Based on DFT Calculations.Based on one of the initial proton configurations, structuralrelaxation calculations were performed using VASP.38−40 Toeffectively obtain accurate stable configurations, multiple-stepstructure relaxation is often performed.41,42 Thus, for thisstudy, three total energies were used as metrics to determine astable configuration. The first is the total energy of the initialproton configuration without the relaxation calculation, whichis expressed as E0 and is known as the no-relaxation totalenergy. The second is the total energy with normal accuracy inthe structural relaxation (E1), and the third is the highlyaccurate total energy after structural relaxation (E2). Thedetailed setup of the VASP calculations employed for multiple-step structural relaxation to obtain the total energy issummarized in Supplementary Note A.Optimization of the Proton Configuration Using theFMQA Algorithm. In quantum annealers and Ising machines,only 0/1 variables can be used; thus, it is necessary to encodeatomistic configurations into a bit string composed of 0/1variables. Oxygen ions can take the form of protonated (i.e.,hydroxy) or nonprotonated species. For protonated oxygenions, a proton is placed at any one of eight possible positions inthe initial state for structural relaxation (Figure 2b). Thus, ninebits were used to express the position of a proton at eachoxygen atom site, and a one-hot constraint was considered toensure that no more than two protons were attached to thesame oxygen. In total, the proton configuration can beexpressed as 24 × 9 = 216 bits, while the bit string can beexpressed as q = {qni}n=1,...,24, i=1,...,9 where qni = 0 or 1. Theindices n and i represent the indices of the oxygen ions andFigure 1. Optimization cycle for the BBO based on use of the FMQAalgorithm and VASP to obtain a stable proton configuration. Theatomistic configurations were drawn using VESTA.36Figure 2. (a) Atomistic structure of the scandium-doped bariumzirconate. (b) Eight stable proton positions. The atomisticconfigurations were drawn using VESTA.36The Journal of Physical Chemistry C pubs.acs.org/JPCC Articlehttps://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−23402333https://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig1&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig1&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig1&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig1&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig2&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig2&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig2&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig2&ref=pdfpubs.acs.org/JPCC?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-astheir proton positions, respectively. For the proton position, i =1,...,8 refers to the positions shown in Figure 2b, while i = 9indicates the absence of a proton. The optimization cycleemployed to obtain a stable proton configuration (i.e., thelowest total energy) using the FMQA algorithm thereforeinvolved four key steps. First, in step 1, M protonconfigurations were generated as the initial training data setfrom 24Cx × 8x configurations, which were used as the initialconfigurations for structural relaxation. Relaxation calculationswere performed using VASP and the total energies wereobtained. These data were used as the training data for the FMmodel, as shown in Figure 1. In step 2, the FM model was usedas a surrogate model to predict the total energy from the initialproton configuration and was trained when bit string q wasinput. The values of the total energy predicted by the FM areas follows:= * + * *= = = = = = =E w q v v q qq( )n ini nin i m j kKnik mjk ni mj1241912419124191(1)where the hyperparameter of K is fixed at a default value of 8.21The trained parameters are denoted wni* and vnik* . In step 3, theFM is solved by the Ising machine under the followingconstraints:==q x24nn1249(2)==q n1,ini19(3)The former constraint determines the number of hydrogenatoms, whereas the latter ensures that at most, one hydrogenatom is attached to one oxygen atom. For the purpose of thisstudy, Fixstars Amplify AE43 was used as the Ising machine. Asingle solution for bit string q* that minimizes the FM modeldefined by eq 1 was selected using the Fixstars Amplify AE. Ifthe selected q* was already included in the training data, thenthe bit string was randomly generated under the constraints ofeqs 2 and 3, which has the effect of exploration in black-boxFigure 3. (a) Eight zirconium ion positions in BaZrO3 and the configuration of the two scandium ions in the Sc@(1,2) structure. (b) Number ofcycles required to obtain the most stable proton configurations from the FMQA algorithm and random sampling using the total energies E0 (norelaxation), E1 (relaxation calculation with normal accuracy), and E2 (relaxation calculation with high accuracy) for the Sc@(1,2) structure. Thehorizontal line in the box represents the median value (50th percentile data), while the top and bottom areas represent the 75th and 25thpercentiles, respectively. The points denote the outliers, while the error bars represent the maximum and minimum values without outliers. (c)Relationship between E0 and E2 for the same initial state. (d) Relationship between E1 and E2 for the same initial state. The atomistic configurationswere drawn using VESTA.36The Journal of Physical Chemistry C pubs.acs.org/JPCC Articlehttps://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−23402334https://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig3&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig3&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig3&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig3&ref=pdfpubs.acs.org/JPCC?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-asoptimization. Finally, in step 4, DFT relaxation calculationswere conducted using the initial proton configuration ex-pressed by the selected q*, and the total energy was obtained.The number of training data points was increased to M + 1.Steps 2−4 were subsequently repeated as desired.The optimization task was defined as finding a superiorinitial state for structural relaxation, denoted by bit string q,which will lead to the minimum total energy after therelaxation calculations, using as few cycles as possible. Notethat in the initial configurations of 24Cx × 8x, multipleconfigurations can be considered to have the same symmetry.Thus, from a symmetry consideration, in a case where theconfiguration suggested by q* in the third step had alreadybeen investigated in the previous cycles, the calculated totalenergy would have been used without the relaxationcalculation, thereby reducing the computation time. Thenumber of structures that are equivalent in terms of symmetrydepends on the doping types and the symmetry of the crystals,so the reduced computation time will be highly system-dependent. However, for the ease of implementation,symmetry was not considered in the current study.■ RESULTS AND DISCUSSIONOptimization Performance for Stable Proton Config-urations. Initially, the optimization performance of the systemwas investigated using a fixed scandium ion configuration. Acase was considered in which the number of scandium ions wastwo, i.e., x = 2. In this case, by considering the symmetry, threepossible configurations exist for the scandium ions for a 2 × 2× 2 supercell. By assigning each zirconium ion position inBaZrO3 with a number between 1 and 8 (see Figure 3a), thescandium ion configurations were denoted as Sc@(1,2)(Figure 2a), Sc@(1,6), and Sc@(1,8), respectively. For eachcase, the number of possible proton configurations as the initialstates for the relaxation calculations was defined as 24Cx × 8x =17,664, without considering symmetry. For all configurations,DFT relaxation calculations were performed in advance using aFugaku supercomputer, and the total energies of E0, E1, and E2were obtained with different accuracies. Thus, for x = 2, theoptimal solutions for the stable proton configurations wereknown for the various relaxation steps, and these solutionswere used to estimate the efficiency of the FMQA algorithm infinding stable proton configurations. In addition, during eachoptimization cycle, the total energy of the selected protonconfiguration can be accessed from the previously calculatedtotal energies. Consequently, it is not necessary to perform anadditional DFT calculation to estimate the efficiency.To verify the BBO performance using the FMQA algorithmfor exploring stable proton configurations, optimization wasperformed for cases where the total energies of E0, E1, and E2were targeted. The number of initial training data points wasfixed at M = 10, and optimization runs were performed untilthe optimal solution was obtained. Ten runs were performedindependently by changing the initial training data selections,and Figure 3b shows box plots of the number of cyclesrequired to obtain the most stable proton configuration forFigure 4. Minimum total energies depending on the number of iteration cycles when (a) E0, (b) E1, and (c) E2 were targeted as the total energiesfor Sc@(3,4,6,7,8). Three independent runs were performed using the FMQA and by random sampling. The means and deviations are indicated bylines and shaded areas. The inset of (a) shows the configuration the scandium ions in Sc@(3,4,6,7,8). (d) Histogram of the total energies with thehigh accuracies (E2), as determined by performing further relaxation calculations based on the configurations found by FMQA optimization usingE0 and E1. The total energy of the proton configuration found by FMQA optimization using E2 is also plotted. The atomistic configurations weredrawn using VESTA.36The Journal of Physical Chemistry C pubs.acs.org/JPCC Articlehttps://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−23402335https://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig4&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig4&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig4&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig4&ref=pdfpubs.acs.org/JPCC?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-asSc@(1,2). Here, the most stable configurations are defined asthose with total energies ranging from the lowest total energyup to 0.005 eV. For comparison, five random samplings wereconducted, starting from each initial training data set used inthe FMQA optimizations; thus, 50 runs were performed.Figure S1 summarizes the results for obtained Sc@(1,6) andSc@(1,8), wherein the most stable proton configurations wereidentified within 150 BBO cycles using the FMQA algorithm.Consequently, it was found that the FMQA outperformed therandom sampling approach. Notably, when random samplingis employed, the number of cycles required to find stableproton configurations for E0 is larger than those required for E1and E2. This can be attributed to the fact that the number ofmost stable configurations for E0 is smaller than those for E1and E2. Furthermore, Figure 3c,d shows the relationshipsamong E0, E1, and E2 for the same initial proton configurations.It can be seen that the order between E0 and E2 is largelyinterchanged when the values of E2 are small, indicating thatfind stable proton configurations using E0, i.e., without the useof a relaxation calculation, are not the true stable protonconfigurations using E2. However, the order is comparablebetween E1 and E2, indicating that optimization using E1 issufficient for finding stable proton configurations in the presentproblem setting.Subsequently, to examine the optimization performance formore complex problems, the case of x = 5 was considered, andthe results are shown in Figure 4. In this case, three possibleconfigurations of scandium ions exist, which were denoted asSc@(3,4,6,7,8) (Figure 4a), Sc@(1,3,4,6,8), andSc@(1,2,4,6,7). When the scandium ion configuration wasfixed, the number of initial candidate structures for therelaxation calculations was 24C5 × 85 = 1,392,771,072 withoutconsidering symmetry, and so it was not possible to calculatethe total energies of all configurations. Consequently, thenumber of iteration cycles of the BBO runs was fixed at 200,including the generation of ten initial training data sets. Threeindependent runs were performed by altering the initialtraining data selections. Figure 4a−c shows the minimum totalenergies over the various cycles for Sc@(3,4,6,7,8) when E0,E1, and E2 were targeted for optimization. In all cases, use ofthe FMQA approach generated more stable structures overfewer cycles compared to those required for random sampling.For the stable configurations determined by FMQAoptimizations using E0 and E1, configurations were changedby performing further relaxations. More specifically, startingfrom the top five configurations with small E0 and E1 values foreach independent FMQA optimization run, further relaxationcalculations were performed, which have the same accuracy asthose performed using E2. The total energies obtained by thehighly accurate relaxation calculations are summarized inFigure 4d. As a result, it was confirmed that stableconfigurations with a value equivalent to the total energyobtained from the FMQA optimizations performed using E2could be obtained from the optimizations using E1. Thisimplies that by performing FMQA optimization using E1,followed by further relaxation for a small number of stableconfigurations, the computation time can be reduced to obtainstable proton configurations. However, as observed for x = 2,the discovery of stable configurations was also meaninglessusing E0. Figure S2 summarizes the results obtained forSc@(1,3,4,6,8) and Sc@(1,2,4,6,7) when the FMQA opti-mization was performed using E1. Further relaxationcalculations were also performed on the identified protonconfigurations, and the total energies obtained with the sameaccuracy as that of E2 were compared. Consequently, theSc@(3,4,6,7,8) structure was defined as possessing the lowesttotal energy.We discuss the computational time required for config-uration optimization by FMQA. The computation timerequired to perform one cycle of FMQA was about 0.1 s forFM training, about 1 s for promising configuration selectionusing Fixstars Amplify, and 900−1800 s for total energycalculations by DFT for each task. Thus, most of thecomputational time was spent on the DFT calculations.Since FMQA algorithm and random sampling do not differsignificantly in overall computation time due to largecomputational time in DFT calculations, we consider FMQAto be a powerful method for finding optimal configurations.Recent Ising machines are capable of handling problems on theorder of 10,000 bits, and optimization calculations usingsimulation cells as large as 1,000 atoms are in principlepossible. However, when considering larger systems, thecomputational time required for DFT calculations becomes abottleneck. To overcome this problem, the use of the effectiveDFT package for large target systems such as OpenMX,44CP2K,45 and CONQUEST46,47 might be one of the work-arounds. In addition, if the accuracy of DFT calculations is nota concern, another solution is to use the neural networkpotentials,48,49 which can dramatically reduce the computa-tional time required for total energy calculations.Discussion Regarding the Stable Proton Configura-tions.When the number of Sc ions (i.e., x) is equal to seven oreight, only one configuration exists wherein the scandium ionscan be substituted for symmetry. However, the number ofinitial proton configurations for the relaxation calculationdrastically increases under these conditions, reaching 1011 and1012 for x = 7 and 8, respectively, without symmetry. In thesecases, owing to the large number of candidates, conventionalBBO methods such as Bayesian optimization cannot beperformed. Thus, stable proton configuration searches wereperformed by FMQA optimization using E1. More specifically,three independent FMQA optimizations were carried out, andfurther relaxation calculations were performed for the top 15configurations. The results of the optimization cycle using E1are shown in Figure S2c,d, and the most stable configurationsare summarized in Figure 5, wherein the numbers of scandiumions (x) are 2, 5, 7, and 8.From the obtained configurations, it was possible to discussthe characteristics of the stable proton configurations. First, inmany cases, hydrogen atoms attach to the oxygen sitesbetween the scandium ions, forming Sc−O−Sc species. Thisindicates that the position of hydrogen is more stable on theoxygen atom of Sc−O−Sc than on Sc−O−Zr or Zr−O−Zr.This stabilization of protons near the scandium ions wasconsistent with the results of previous ab initio calcula-tions50,51,37 and nuclear magnetic resonance experiments.52Subsequently, hydrogen bond pairs appear, as denoted by thered triangles in Figure 5. Figure 6 compares the most stableconfiguration containing a pair of hydrogen bonds and ametastable configuration containing two independent hydro-gen bonds when x = 2. The lengths of the hydrogen bondswere calculated to be 1.89 and 1.81 Å in the most stable case,which are longer than those calculated for the metastable case,i.e., 1.73 and 1.67 Å. These results indicate that the metastablestate has a lower energy (i.e., it is more stable) if the length ofthe hydrogen bonds is considered. In contrast, latticeThe Journal of Physical Chemistry C pubs.acs.org/JPCC Articlehttps://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−23402336https://pubs.acs.org/doi/suppl/10.1021/acs.jpcc.4c07104/suppl_file/jp4c07104_si_001.pdfhttps://pubs.acs.org/doi/suppl/10.1021/acs.jpcc.4c07104/suppl_file/jp4c07104_si_001.pdfhttps://pubs.acs.org/doi/suppl/10.1021/acs.jpcc.4c07104/suppl_file/jp4c07104_si_001.pdfpubs.acs.org/JPCC?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-asdeformation was found to increase the total energy, and in themost stable state, the lattice deformation was reduced by theformation of a pair of hydrogen bonds. This property offorming a hydrogen bond pair in a stable proton configurationwas observed even when the number of scandium ions wasincreased (see Figure 5). The contributions of latticedeformation and the formation of a pair of hydrogen bondsfor total energy depend on simulation conditions. Future workwill therefore examine whether the characteristics of the stablestructures observed in this study are universal for simulationconditions using the FMQA algorithm.Simultaneous Optimization of the Doping andProton Configurations. It is also possible to optimize thedoping configuration of scandium ions simultaneously with theproton configuration. The doping configuration should also betreated as a bit string, and so the bits for the scandium/zirconium sites were prepared as 1 for scandium ions and 0 forzirconium ions. In the 2 × 2 × 2 supercell model, eight bitswere used to represent the doping configuration, which wasdefined by {ql}l=1,...,8. In other words, when simultaneouslyconsidering the doping and proton configurations, in additionto the 216-bit, {qni}n=1,...,24, i=1,...,9 was also employed for doping.Therefore, a 224-bit FM was utilized as a surrogate model topredict the total energy. When solving the trained FM with theIsing machine, in addition to the constraints defined in eqs 2and 3, the following constraint is required:==q xll18(4)However, in the present study, symmetry was observed in thedoping positions. For example, for x = 2, one scandium ion canbe fixed at index 1, and the remaining scandium ions can beplaced at one of the three sites (i.e., 2, 6, or 8) indicated inFigure 3a. Therefore, instead of using eq 4, the followingconstraints can be used to narrow the search space, whereinonly three bits are considered (i.e., q2, q6, and q8):=q 11 (5)==q 1ll2,6,8 (6)In the case where x = 2, the doping and protonconfigurations were simultaneously optimized by BBO usingFMQA, and the results are presented in Figure 7 for searchingwith the total energies E0, E1, and E2. The results of randomsampling are also shown for comparison, and the optimizationperformance of FMQA was observed to outperform that ofrandom sampling.■ CONCLUSIONSIn this study, black-box optimization (BBO) calculations wereperformed to determine stable proton configurations using thefactorization machines with quantum annealing (FMQA)algorithm combined with the Vienna Ab initio SimulationPackage (VASP). A graphics processing unit-based Isingmachine was used to select a promising initial state for therelaxation calculations. As a model system for protonconductors, scandium-doped barium zirconate was used,wherein the number of scandium ions (x) was set as 2, 5, 7,or 8 in the 2 × 2 × 2 supercell model. For x = 2, the number ofpossible proton configurations was 52,992, and so it waspossible to calculate the total energies for all cases using VASP,and the most stable configuration was identified. In this case,the optimization performance of FMQA algorithm alwaysoutperformed that of random sampling for finding the moststable configurations. Subsequently, using the FMQA algo-rithm, the proton configurations were optimized for x = 5, 7,and 8, wherein the numbers of initial configurations were inFigure 5. Most stable proton configurations found for (a) 2, (b) 5,(c) 7, and (d) 8, scandium ions. When x = 2, the optimalconfiguration is shown. The red triangles indicate pairs of hydrogenbonds. Barium ions are omitted for clarity. The atomisticconfigurations were drawn using VESTA.36Figure 6. Hydrogen bonds in the (a) most stable and (b) metastableconfigurations when x = 2. The atomistic configurations were drawnusing VESTA.36Figure 7. Number of cycles required when the most stable protonconfiguration was identified by FMQA optimization and randomsampling using the total energies E0 (no relaxation), E1 (relaxationcalculation with normal accuracy), and E2 (relaxation calculation withhigh accuracy) for optimization of the doped scandium ion andproton configurations.The Journal of Physical Chemistry C pubs.acs.org/JPCC Articlehttps://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−23402337https://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig5&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig5&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig5&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig5&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig6&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig6&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig6&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig6&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig7&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig7&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig7&ref=pdfhttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?fig=fig7&ref=pdfpubs.acs.org/JPCC?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-asthe order of 109, 1011, and 1012, respectively. Notably, such vastnumber of possible configurations are difficult to treat usingconventional BBO methods, such as Bayesian optimization,because of the computation times required. Moreover, theoptimization results identified the formation of hydrogen bondpairs in the stable proton configurations, regardless of thenumber of scandium ions. Notably, simultaneous optimizationof the doped scandium ions and proton configurations wasachieved. It is therefore expected that this optimizationstrategy based on a combination of Ising machines and densityfunctional theory calculations will be able to solve variousissues related to configuration optimization in solid-statematerials, and promote novel scientific discoveries. Indeed, ifthe proton and dopant configurations are given discretely, ourmethod is useful for most cases without exception. In addition,the case of multiple proton configurations in a single oxygencan be handled by introducing inequality constraints forQUBO, and various types of configuration optimization areachieved. However, if the arrangement of protons and dopantatoms is continuous, our method cannot handle it directly, andanother development method is required.We show that FMQA can handle large search spaces, e.g., onthe order of 1012. However, it is difficult to determine whetherthe stable structure found by FMQA is the most stablestructure or not, because it is impossible to compute allconfigurations by DFT calculations. Here, we discuss thestrategies for determining whether the most stable structurehas been obtained. The first is to compare the results of thesimulations with the experimental results. Molecular dynamicscalculation for proton conductivities is possible by adoptingstable proton configurations found by FMQA optimizations asinitial structures. We believe that the comparison of protonconductivities by simulations and experiments will allow us toverify the correctness of the stable structure obtained from theFMQA optimizations. The second is the introduction of astopping criterion for the black-box optimization. The stoppingcriterion for Bayesian optimization leading to the effectivefinding of the optimal solution is controversial.53,54 Concerningthese studies, introducing the stopping criterion for the FMQAalgorithm is interesting, and we expect to be able to terminatethe optimization tasks when an optimal solution is found.■ ASSOCIATED CONTENT*sı Supporting InformationThe Supporting Information is available free of charge athttps://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104.Detailed setup for the VASP calculations; parametersemployed for the VASP input files for each relaxationstep; optimization results for Sc@(1,6) and Sc@(1,8)structures; and optimization results depending on thenumber of Sc atoms (PDF)■ AUTHOR INFORMATIONCorresponding AuthorsJianbo Lin − Center for Basic Research on Materials, NationalInstitute for Materials Science, Ibaraki 305-0044, Japan;Email: linjb86@gmail.comTomofumi Tada − Kyushu University Platform of Inter/Transdisciplinary Energy Research, Kyushu University,Fukuoka 819-0395, Japan; orcid.org/0000-0003-3093-3779; Email: tada.tomofumi.054@m.kyushu-u.ac.jpRyo Tamura − Center for Basic Research on Materials,National Institute for Materials Science, Ibaraki 305-0044,Japan; RIKEN Center for Advanced Intelligence Project,Tokyo 103-0027, Japan; Graduate School of FrontierSciences, The University of Tokyo, Chiba 277-8561, Japan;orcid.org/0000-0002-0349-358X; Email: tamura.ryo@nims.go.jpAuthorsAi Koizumi − Center for Basic Research on Materials,National Institute for Materials Science, Ibaraki 305-0044,JapanMasato Sumita − RIKEN Center for Advanced IntelligenceProject, Tokyo 103-0027, Japan; orcid.org/0000-0002-3506-1028Koji Tsuda − Graduate School of Frontier Sciences, TheUniversity of Tokyo, Chiba 277-8561, Japan; Center forBasic Research on Materials, National Institute for MaterialsScience, Ibaraki 305-0044, Japan; RIKEN Center forAdvanced Intelligence Project, Tokyo 103-0027, Japan;orcid.org/0000-0002-4288-1606Complete contact information is available at:https://pubs.acs.org/10.1021/acs.jpcc.4c07104Author ContributionsT.T. and R.T. conceived and designed the study. J.L. and A.K.conducted the investigations and calculations. R.T., M.S., andK.T. devised the original idea for the optimization method. Allmembers contributed to manuscript preparation.FundingThis study was supported by a project subsidized by the CoreResearch for Evolutional Science and Technology (CREST)(grant number JPMJCR2234) and KAKENHI 21H01008. Thecomputational DFT work was partly supported by the MEXTProgram: Data Creation and Utilization Type MaterialResearch and Development Project (grant numberJPMXP1122683430) and the MEXT Program for PromotingResearches on the Supercomputer Fugaku (Data-DrivenResearch Methods Development and Materials InnovationLed by Computational Materials Science; grant numberJPMXP1020230327).NotesThe authors declare no competing financial interest.■ ACKNOWLEDGMENTSThe authors would like to thank Yoshihiro Yamazaki andShusuke Kasamatsu for their fruitful discussions. Thecomputational resources of the Fugaku supercomputer wereprovided by the RIKEN Center for Computational Science(Project ID: hp230212 and hp240223). The computationsperformed in this study were performed on an Ising machine inFixstars Amplify. We would like to thank Editage (www.edi-tage.jp) for English language editing.■ REFERENCES(1) Iwahara, H.; Esaka, T.; Uchida, H.; Maeda, N. ProtonConduction in Sintered Oxides and Its Application to SteamElectrolysis for Hydrogen Production. Solid State Ionics 1981, 3−4,359−363.(2) Norby, T. Solid-State Protonic Conductors: Principles, Proper-ties, Progress and Prospects. Solid State Ionics 1999, 125 (1), 1−11.(3) Kreuer, K. D. Proton-Conducting Oxides. Annu. Rev. Mater. Res.2003, 33, 333−359.The Journal of Physical Chemistry C pubs.acs.org/JPCC Articlehttps://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−23402338https://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?goto=supporting-infohttps://pubs.acs.org/doi/suppl/10.1021/acs.jpcc.4c07104/suppl_file/jp4c07104_si_001.pdfhttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Jianbo+Lin"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfmailto:linjb86@gmail.comhttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Tomofumi+Tada"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://orcid.org/0000-0003-3093-3779https://orcid.org/0000-0003-3093-3779mailto:tada.tomofumi.054@m.kyushu-u.ac.jphttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Ryo+Tamura"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://orcid.org/0000-0002-0349-358Xhttps://orcid.org/0000-0002-0349-358Xmailto:tamura.ryo@nims.go.jpmailto:tamura.ryo@nims.go.jphttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Ai+Koizumi"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://pubs.acs.org/action/doSearch?field1=Contrib&text1="Masato+Sumita"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://orcid.org/0000-0002-3506-1028https://orcid.org/0000-0002-3506-1028https://pubs.acs.org/action/doSearch?field1=Contrib&text1="Koji+Tsuda"&field2=AllField&text2=&publication=&accessType=allContent&Earliest=&ref=pdfhttps://orcid.org/0000-0002-4288-1606https://orcid.org/0000-0002-4288-1606https://pubs.acs.org/doi/10.1021/acs.jpcc.4c07104?ref=pdfhttp://www.editage.jphttp://www.editage.jphttps://doi.org/10.1016/0167-2738(81)90113-2https://doi.org/10.1016/0167-2738(81)90113-2https://doi.org/10.1016/0167-2738(81)90113-2https://doi.org/10.1016/S0167-2738(99)00152-6https://doi.org/10.1016/S0167-2738(99)00152-6https://doi.org/10.1146/annurev.matsci.33.022802.091825pubs.acs.org/JPCC?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-as(4) Yamazaki, Y.; Hernandez-Sanchez, R.; Haile, S. M. High TotalProton Conductivity in Large-Grained Yttrium-Doped BariumZirconate. Chem. Mater. 2009, 21 (13), 2755−2762.(5) Hyodo, J.; Kitabayashi, K.; Hoshino, K.; Okuyama, Y.; Yamazaki,Y. Fast and Stable Proton Conduction in Heavily Scandium-DopedPolycrystalline Barium Zirconate at Intermediate Temperatures. Adv.Energy Mater. 2020, 10 (25), No. 2000213.(6) Cao, J.; Ji, Y.; Shao, Z. Perovskites for Protonic Ceramic FuelCells: A Review. Energy Environ. Sci. 2022, 15 (6), 2200−2232.(7) Yamazaki, Y.; Babilo, P.; Haile, S. M. Defect Chemistry ofYttrium-Doped Barium Zirconate: A Thermodynamic Analysis ofWater Uptake. Chem. Mater. 2008, 20 (20), 6352−6357.(8) Yamazaki, Y.; Yang, C.-K.; Haile, S. M. Unraveling the DefectChemistry and Proton Uptake of Yttrium-Doped Barium Zirconate.Scr. Mater. 2011, 65 (2), 102−107.(9) Kitamura, N.; Akola, J.; Kohara, S.; Fujimoto, K.; Idemoto, Y.Proton Distribution and Dynamics in Y- and Zn-Doped BaZrO3. J.Phys. Chem. C 2014, 118 (33), 18846−18852.(10) Hashizume, K.; Hossain, M.; Khan, M.; Akash, T.; Hashizume,K. Molecular Dynamics Simulation For Barium Zirconate ProtonConducting Oxide Energy Materials: A Mini Review. In InternationalExchange and Innovation Conference on Engineering & Science(IEICES), 2021; Vol. 7, p 1.(11) Ju, S.; Shiga, T.; Feng, L.; Hou, Z.; Tsuda, K.; Shiomi, J.Designing Nanostructures for Phonon Transport via BayesianOptimization. Phys. Rev. X 2017, 7 (2), No. 021024.(12) Terayama, K.; Sumita, M.; Tamura, R.; Tsuda, K. Black-BoxOptimization for Automated Discovery. Acc. Chem. Res. 2021, 54 (6),1334−1346.(13) Ueno, T.; Rhone, T. D.; Hou, Z.; Mizoguchi, T.; Tsuda, K.COMBO: An Efficient Bayesian Optimization Library for MaterialsScience. Mater. Discovery 2016, 4, 18−21.(14) Motoyama, Y.; Tamura, R.; Yoshimi, K.; Terayama, K.; Ueno,T.; Tsuda, K. Bayesian Optimization Package: PHYSBO. Comput.Phys. Commun. 2022, 278, No. 108405.(15) Garijo del Río, E.; Mortensen, J. J.; Jacobsen, K. W. LocalBayesian Optimizer for Atomic Structures. Phys. Rev. B 2019, 100(10), No. 104103.(16) Todorovic,́ M.; Gutmann, M. U.; Corander, J.; Rinke, P.Bayesian Inference of Atomistic Structure in Functional Materials. npjComput. Mater. 2019, 5 (1), 1−7.(17) Yan, J.; Wei, H.; Xie, H.; Gu, X.; Bao, H. Seeking for LowThermal Conductivity Atomic Configurations in SiGe Alloys withBayesian Optimization. ES Energy Environ. 2020, 8 (4), 56−64.(18) Ono, S. Optimization of Configurations of Atomic Species onTwo-Dimensional Hexagonal Lattices for Copper-Based Systems. AIPAdv. 2022, 12 (8), No. 085313.(19) Kusaba, A.; Kangawa, Y.; Kuboyama, T.; Oshiyama, A.Exploration of a Large-Scale Reconstructed Structure onGaN(0001) Surface by Bayesian Optimization. Appl. Phys. Lett.2022, 120 (2), No. 021602.(20) Kitai, K.; Guo, J.; Ju, S.; Tanaka, S.; Tsuda, K.; Shiomi, J.;Tamura, R. Designing Metamaterials with Quantum Annealing andFactorization Machines. Phys. Rev. Res. 2020, 2 (1), No. 013319.(21) Rendle, S. Factorization Machines with Libfm. ACM Trans.Intell. Syst. Technol. 2012, 3 (3), 57.(22) Mohseni, N.; McMahon, P. L.; Byrnes, T. Ising Machines asHardware Solvers of Combinatorial Optimization Problems. Nat. Rev.Phys. 2022, 4 (6), 363−379.(23) Wilson, B. A.; Kudyshev, Z. A.; Kildishev, A. V.; Kais, S.;Shalaev, V. M.; Boltasseva, A. Machine Learning Framework forQuantum Sampling of Highly Constrained, Continuous OptimizationProblems. Appl. Phys. Rev. 2021, 8 (4), No. 041418.(24) Kim, S.; Shang, W.; Moon, S.; Pastega, T.; Lee, E.; Luo, T.High-Performance Transparent Radiative Cooler Designed byQuantum Computing. ACS Energy Lett. 2022, 7, 4134−4141.(25) Matsumori, T.; Taki, M.; Kadowaki, T. Application of QUBOSolver Using Black-Box Optimization to Structural Design forResonance Avoidance. Sci. Rep. 2022, 12 (1), 12143.(26) Inoue, T.; Seki, Y.; Tanaka, S.; Togawa, N.; Ishizaki, K.; Noda,S. Towards Optimization of Photonic-Crystal Surface-Emitting Lasersvia Quantum Annealing. Opt. Express, OE 2022, 30 (24), 43503−43512.(27) Urushihara, M.; Karube, M.; Yamaguchi, K.; Tamura, R.Optimization of Core−Shell Nanoparticles Using a Combination ofMachine Learning and Ising Machine. Adv. Photon. Res. 2023, 4 (12),No. 2300226.(28) Mao, Z.; Matsuda, Y.; Tamura, R.; Tsuda, K. Chemical Designwith GPU-Based Ising Machines. Digit. Discov. 2023, 2 (4), 1098−1103.(29) Tucš, A.; Berenger, F.; Yumoto, A.; Tamura, R.; Uzawa, T.;Tsuda, K. Quantum Annealing Designs Nonhemolytic AntimicrobialPeptides in a Discrete Latent Space. ACS Med. Chem. Lett. 2023, 14(5), 577−582.(30) Kim, S.; Jung, S.; Bobbitt, A.; Lee, E.; Luo, T. Wide-AngleSpectral Filter for Energy-Saving Windows Designed by QuantumAnnealing-Enhanced Active Learning. Cell Rep. Phys. Sci. 2024, 5 (3),No. 101847.(31) Nawa, K.; Suzuki, T.; Masuda, K.; Tanaka, S.; Miura, Y.Quantum Annealing Optimization Method for the Design of BarrierMaterials in Magnetic Tunnel Junctions. Phys. Rev. Appl. 2023, 20 (2),No. 024044.(32) Suga, Y.; Maruo, A.; Jippo, H. A Feasibility Study for QuantumComputing Methodologies in Automotive Advanced MaterialInvestigation. Trans. Soc. Automot. Eng. Jpn. 2024, 55 (3), 621−627.(33) Couzinie, Y.; Seki, Y.; Nishiya, Y.; Nishi, H.; Kosugi, T.;Tanaka, S.; Matsushita, Y. Machine Learning Supported Annealing forPrediction of Grand Canonical Crystal Structures. arXiv:2408.035562024.(34) Luo, T.; Xu, Z.; Shang, W.; Kim, S.; Lee, E. QALO: QuantumAnnealing-Assisted Lattice Optimization. npj Comput. Mater. 2024,11, 4.(35) Hoshino, K.; Kasamatsu, S.; Hyodo, J.; Yamamoto, K.;Setoyama, H.; Okajima, T.; Yamazaki, Y. Probing Local Environmentsof Oxygen Vacancies Responsible for Hydration in Sc-Doped BariumZirconates at Elevated Temperatures: In Situ X-Ray AbsorptionSpectroscopy, Thermogravimetry, and Active Learning Ab InitioReplica Exchange Monte Carlo Simulations. Chem. Mater. 2023, 35(6), 2289−2301.(36) Momma, K.; Izumi, F. VESTA 3 for Three-DimensionalVisualization of Crystal, Volumetric and Morphology Data. J. Appl.Crystallogr. 2011, 44 (6), 1272−1276.(37) Vera, C. Y. R.; Ding, H.; Peterson, D.; Gibbons, W. T.; Zhou,M.; Ding, D. A Mini-Review on Proton Conduction of BaZrO3-BasedPerovskite Electrolytes. J. Phys. Energy 2021, 3 (3), No. 032019.(38) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for LiquidMetals. Phys. Rev. B 1993, 47 (1), 558−561.(39) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for AbInitio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys.Rev. B 1996, 54 (16), 11169−11186.(40) Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total EnergyCalculations for Metals and Semiconductors Using a Plane-WaveBasis Set. Comput. Mater. Sci. 1996, 6 (1), 15−50.(41) Peressi, M.; Baldereschi, A. Chapter 2 - Ab Initio Studies ofStructural and Electronic Properties. In Characterization of Semi-conductor Heterostructures and Nanostructures, 2nd ed.; Lamberti, C.;Agostini, G., Eds.; Elsevier: Oxford, 2013; pp 21−73.(42) Musielewicz, J.; Wang, X.; Tian, T.; Ulissi, Z. FINETUNA:Fine-Tuning Accelerated Molecular Simulations. Mach. Learn.: Sci.Technol. 2022, 3 (3), No. 03LT01.(43) Fixstars Amplify. https://amplify.fixstars.com/en/.(44) Ozaki, T. Variationally Optimized Atomic Orbitals for Large-Scale Electronic Structures. Phys. Rev. B 2003, 67 (15), No. 155108.(45) Kühne, T. D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V. V.;Seewald, P.; Stein, F.; Laino, T.; Khaliullin, R. Z.; Schütt, O.;Schiffmann, F.; Golze, D.; Wilhelm, J.; Chulkov, S.; Bani-Hashemian,M. H.; Weber, V.; Borsťnik, U.; Taillefumier, M.; Jakobovits, A. S.;Lazzaro, A.; Pabst, H.; Müller, T.; Schade, R.; Guidon, M.;The Journal of Physical Chemistry C pubs.acs.org/JPCC Articlehttps://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−23402339https://doi.org/10.1021/cm900208w?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/cm900208w?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/cm900208w?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1002/aenm.202000213https://doi.org/10.1002/aenm.202000213https://doi.org/10.1039/D2EE00132Bhttps://doi.org/10.1039/D2EE00132Bhttps://doi.org/10.1021/cm800843s?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/cm800843s?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/cm800843s?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1016/j.scriptamat.2010.12.034https://doi.org/10.1016/j.scriptamat.2010.12.034https://doi.org/10.1021/jp502455v?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1103/PhysRevX.7.021024https://doi.org/10.1103/PhysRevX.7.021024https://doi.org/10.1021/acs.accounts.0c00713?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/acs.accounts.0c00713?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1016/j.md.2016.04.001https://doi.org/10.1016/j.md.2016.04.001https://doi.org/10.1016/j.cpc.2022.108405https://doi.org/10.1103/PhysRevB.100.104103https://doi.org/10.1103/PhysRevB.100.104103https://doi.org/10.1038/s41524-019-0175-2https://doi.org/10.30919/esee8c356https://doi.org/10.30919/esee8c356https://doi.org/10.30919/esee8c356https://doi.org/10.1063/5.0098517https://doi.org/10.1063/5.0098517https://doi.org/10.1063/5.0078660https://doi.org/10.1063/5.0078660https://doi.org/10.1103/PhysRevResearch.2.013319https://doi.org/10.1103/PhysRevResearch.2.013319https://doi.org/10.1145/2168752.2168771https://doi.org/10.1038/s42254-022-00440-8https://doi.org/10.1038/s42254-022-00440-8https://doi.org/10.1063/5.0060481https://doi.org/10.1063/5.0060481https://doi.org/10.1063/5.0060481https://doi.org/10.1021/acsenergylett.2c01969?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/acsenergylett.2c01969?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1038/s41598-022-16149-8https://doi.org/10.1038/s41598-022-16149-8https://doi.org/10.1038/s41598-022-16149-8https://doi.org/10.1364/OE.476839https://doi.org/10.1364/OE.476839https://doi.org/10.1002/adpr.202300226https://doi.org/10.1002/adpr.202300226https://doi.org/10.1039/D3DD00047Hhttps://doi.org/10.1039/D3DD00047Hhttps://doi.org/10.1021/acsmedchemlett.2c00487?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/acsmedchemlett.2c00487?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1016/j.xcrp.2024.101847https://doi.org/10.1016/j.xcrp.2024.101847https://doi.org/10.1016/j.xcrp.2024.101847https://doi.org/10.1103/PhysRevApplied.20.024044https://doi.org/10.1103/PhysRevApplied.20.024044https://doi.org/10.11351/jsaeronbun.55.621https://doi.org/10.11351/jsaeronbun.55.621https://doi.org/10.11351/jsaeronbun.55.621https://doi.org/10.48550/arXiv.2408.03556https://doi.org/10.48550/arXiv.2408.03556https://doi.org/10.21203/rs.3.rs-4518513/v1https://doi.org/10.21203/rs.3.rs-4518513/v1https://doi.org/10.1021/acs.chemmater.2c02116?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/acs.chemmater.2c02116?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/acs.chemmater.2c02116?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/acs.chemmater.2c02116?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/acs.chemmater.2c02116?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1107/S0021889811038970https://doi.org/10.1107/S0021889811038970https://doi.org/10.1088/2515-7655/ac12abhttps://doi.org/10.1088/2515-7655/ac12abhttps://doi.org/10.1103/PhysRevB.47.558https://doi.org/10.1103/PhysRevB.47.558https://doi.org/10.1103/PhysRevB.54.11169https://doi.org/10.1103/PhysRevB.54.11169https://doi.org/10.1016/0927-0256(96)00008-0https://doi.org/10.1016/0927-0256(96)00008-0https://doi.org/10.1016/0927-0256(96)00008-0https://doi.org/10.1088/2632-2153/ac8fe0https://doi.org/10.1088/2632-2153/ac8fe0https://amplify.fixstars.com/en/https://doi.org/10.1103/PhysRevB.67.155108https://doi.org/10.1103/PhysRevB.67.155108pubs.acs.org/JPCC?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-asAndermatt, S.; Holmberg, N.; Schenter, G. K.; Hehn, A.; Bussy, A.;Belleflamme, F.; Tabacchi, G.; Glöß, A.; Lass, M.; Bethune, I.; Mundy,C. J.; Plessl, C.; Watkins, M.; VandeVondele, J.; Krack, M.; Hutter, J.CP2K: An Electronic Structure and Molecular Dynamics SoftwarePackage - Quickstep: Efficient and Accurate Electronic StructureCalculations. J. Chem. Phys. 2020, 152 (19), 194103.(46) Bowler, D. R.; Miyazaki, T.; Gillan, M. J. Recent Progress inLinear Scaling Ab Initio Electronic Structure Techniques. J. Phys.:Condens. Matter 2002, 14 (11), 2781.(47) Nakata, A.; Baker, J. S.; Mujahed, S. Y.; Poulton, J. T. L.;Arapan, S.; Lin, J.; Raza, Z.; Yadav, S.; Truflandier, L.; Miyazaki, T.;Bowler, D. R. Large Scale and Linear Scaling DFT with theCONQUEST Code. J. Chem. Phys. 2020, 152 (16), 164112.(48) Behler, J. Constructing High-Dimensional Neural NetworkPotentials: A Tutorial Review. Int. J. Quantum Chem. 2015, 115 (16),1032−1050.(49) Takamoto, S.; Shinagawa, C.; Motoki, D.; Nakago, K.; Li, W.;Kurata, I.; Watanabe, T.; Yayama, Y.; Iriguchi, H.; Asano, Y.;Onodera, T.; Ishii, T.; Kudo, T.; Ono, H.; Sawada, R.; Ishitani, R.;Ong, M.; Yamaguchi, T.; Kataoka, T.; Hayashi, A.; Charoenphakdee,N.; Ibuka, T. Towards Universal Neural Network Potential forMaterial Discovery Applicable to Arbitrary Combination of 45Elements. Nat. Commun. 2022, 13 (1), 2991.(50) Björketun, M. E.; Sundell, P. G.; Wahnström, G. Structure andThermodynamic Stability of Hydrogen Interstitials in BaZrO3Perovskite Oxide from Density Functional Calculations. FaradayDiscuss. 2007, 134 (0), 247−265.(51) Yamazaki, Y.; Kuwabara, A.; Hyodo, J.; Okuyama, Y.; Fisher, C.A. J.; Haile, S. M. Oxygen Affinity: The Missing Link EnablingPrediction of Proton Conductivities in Doped Barium Zirconates.Chem. Mater. 2020, 32 (17), 7292−7300.(52) Buannic, L.; Sperrin, L.; Derviso̧ğlu, R.; Blanc, F.; Grey, C. P.Proton Distribution in Sc-Doped BaZrO3: A Solid State NMR andFirst Principle Calculations Analysis. Phys. Chem. Chem. Phys. 2018,20 (6), 4317−4328.(53) Makarova, A.; Shen, H.; Perrone, V.; Klein, A.; Faddoul, J. B.;Krause, A.; Seeger, M.; Archambeau, C. Automatic Termination forHyperparameter Optimization. arXiv:2104.08166 2022.(54) Ishibashi, H.; Karasuyama, M.; Takeuchi, I.; Hino, H. AStopping Criterion for Bayesian Optimization by the Gap of ExpectedMinimum Simple Regrets. In Proceedings of The 26th InternationalConference on Artificial Intelligence and Statistics; PMLR, 2023; pp6463−6497.The Journal of Physical Chemistry C pubs.acs.org/JPCC Articlehttps://doi.org/10.1021/acs.jpcc.4c07104J. Phys. Chem. C 2025, 129, 2332−23402340https://doi.org/10.1063/5.0007045https://doi.org/10.1063/5.0007045https://doi.org/10.1063/5.0007045https://doi.org/10.1088/0953-8984/14/11/303https://doi.org/10.1088/0953-8984/14/11/303https://doi.org/10.1063/5.0005074https://doi.org/10.1063/5.0005074https://doi.org/10.1002/qua.24890https://doi.org/10.1002/qua.24890https://doi.org/10.1038/s41467-022-30687-9https://doi.org/10.1038/s41467-022-30687-9https://doi.org/10.1038/s41467-022-30687-9https://doi.org/10.1039/B602081Jhttps://doi.org/10.1039/B602081Jhttps://doi.org/10.1039/B602081Jhttps://doi.org/10.1021/acs.chemmater.0c01869?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1021/acs.chemmater.0c01869?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-ashttps://doi.org/10.1039/C7CP08523Khttps://doi.org/10.1039/C7CP08523Khttps://doi.org/10.48550/arXiv.2104.08166https://doi.org/10.48550/arXiv.2104.08166pubs.acs.org/JPCC?ref=pdfhttps://doi.org/10.1021/acs.jpcc.4c07104?urlappend=%3Fref%3DPDF&jav=VoR&rel=cite-as