Improvement of look ahead based on quadratic approximation for crystal structure prediction

ABSTRACT Crystal structure prediction based on first-principles calculations is usually time-consuming since a lot of candidate structures have to be locally optimized. Look Ahead based on Quadratic Approximation, which is one of the selection-type algorithms we previously developed, can control the optimization priority of the candidates. It can efficiently reduce the computational costs in many cases, however, it has been found to be inefficient in some data sets. In the present study, we proposed an improved score of Look Ahead based on Quadratic Approximation, where the stress term is added to overcome the drawbacks of the previous score. Crystal structure prediction simulations by this improved algorithm are performed to investigate the searching efficiency for typical materials such as Si, Al2O3, NaCl, and SrCO3. These results demonstrate that this improved algorithm reduces the searching cost to less than 40% with respect to random search in most cases. The introduction of the stress term makes this algorithm more robust and versatile. GRAPHICAL ABSTRACT


Introduction
Crystal structure prediction (CSP) becomes a key to the advancement of new materials design, with the development of computation capabilities. It has become possible to predict the stable crystal structure for given chemical composition using searching algorithms in conjunction with first-principles calculations. However, CSP is basically a quite difficult global optimization problem due to the exponential increase of the number of potential energy minima with respect to the system size. A great deal of effort has been put into overcoming this problem. So far, several searching algorithms in CSP have been successfully developed such as random search (RS) [1][2][3][4], simulated annealing [5,6], minima hopping [7,8], evolutionary algorithm (EA) [9][10][11][12], and particle swarm optimization (PSO) [13,14]. Although significant advances of the searching algorithms have been made over the years, CSP simulations are computationally expensive even now. Firstprinciples calculations for local structure optimization and energy evaluation are time-consuming in CSP. Previously, we proposed a unique method using selection-type algorithms [15][16][17][18] to reduce computational costs. All generated structures have to be locally optimized one by one in order in the conventional methods, whereas we can efficiently select good candidates for priority optimization using machine learning in the selection-type algorithms. We have proposed two selection-type algorithms: Bayesian optimization (BO) and Look Ahead based on Quadratic Approximation (LAQA). Our CSP simulations with BO demonstrated a reduction of the number of searching trials required to find the global minimum structure by 30-40% in comparison with RS [15]. The searching efficiency strongly depends on a structure descriptor, which is required for machine learning. For non-experts, there are difficulties in determining the parameters for descriptors.
To simplify and improve the efficiency, we have also proposed another selection-type algorithm, LAQA [17]. LAQA is classified as reinforcement learning such as multi-armed bandit algorithms. Unlike supervised learning, no training data is required in LAQA. Usually, in CSP, first-principles calculations are performed to evaluate the force acting on each atom and the stress tensor for a given structure. Local structure optimization is carried out using gradient methods such as the quasi-Newton or conjugate gradient methods. However, we can skip unpromising calculations by stopping their optimizations if the final relaxed energy of a structure can be estimated in the middle of its calculation. By controlling the priority of the optimizations, LAQA allows us to find the most stable structure with a small number of total optimization steps. In LAQA, local structure optimization is not performed all at once. Candidates are selected according to the score calculated from the current energy and force, and only a few optimization steps for the selected candidates proceed. It should be noted that LAQA does not require a structure descriptor.
LAQA can efficiently reduce the computational costs in many cases, however, it has been found that LAQA is inefficient in some data sets. In the present study, we propose improved LAQA to overcome the drawback, where the stress term is added. We perform CSP simulations by RS and LAQA using the CrySPY code [18], testing several parameters. Specifically, we compare the number of total optimization steps required to find the most stable structure in each searching algorithm for the systems of Si of which the unit cell is composed of 16 atoms (hereinafter notated as Si 16 ), Al 4 O 6 , Na 8 Cl 8 , and Sr 4 C 4 O 12 . Moreover, we demonstrate that the introduction of the stress term makes LAQA more robust and versatile.

Random search
In this study, we first perform CSP simulations with RS to prepare the data sets for the typical systems and to compare the searching efficiency with LAQA. In our RS method with CrySPY, initial structures with specific chemical composition are randomly generated. CrySPY employs PyXtal [19] internally to generate random structures with specific space groups. In the random structure generation, a space group and a combination of the Wyckoff positions are randomly selected. The lattice parameters are automatically set to appropriate values depending on the number of atoms. The generated initial structures are locally optimized one by one using first-principles calculations.
For the system of Sr 4 C 4 O 12 , we used the VASP code [20][21][22][23] to evaluate the energy and to carry out the local structure optimization. For the other systems, the QUANTUM ESPRESSO (QE) [24][25][26] with pseudopotentials from PSlibrary [27] was also used. The totalenergy calculations were performed using the density functional theory (DFT) with the projectoraugmented wave (PAW) method [28,29]. The generalized gradient approximation (GGA) by Perdew, Burke, and Ernzerhof [30] was used for exchange-correlation functional. The internal atomic coordinates, as well as the cell parameters, were fully optimized using the conjugate gradient algorithm and BFGS quasinewton algorithm for VASP and QE, respectively. Cutoff energies of 789, 830, 1170, and 500 eV for the plane-wave expansion of the wave function were adopted for Si (both Si 8 and Si 16 ), Al 4 O 6 , Na 8 Cl 8 , and Sr 4 C 4 O 12 , respectively. k-point meshes with the density of 80 Å −3 were automatically generated in CrySPY with pymatgen [31]. In RS, the energy, force, and stress data were saved for each optimization step for LAQA.

Laqa
LAQA allows us to find the most stable structure with a small number of total optimization steps by controlling local optimization. Figure 1 illustrates the basic idea of acceleration of CSP by LAQA. In LAQA, local structure optimization of the selected candidates is not performed all at once. Only a few optimization steps for the selected candidates proceed. One of the simplest ways to choose a candidate is to decide on the order of energies, which is referred to as the greedy method (Greedy). However, Greedy does not always work since a candidate that reaches high energy even after the local optimization is possibly selected. In the example of Figure 1, candidate A is preferentially selected, resulting in wasteful computational costs. Therefore, we have introduced the score, which is based on a quadratic approximation from the current energy and force acting on each atom, to roughly estimate the final energy [17,18]. The score L is expressed by the following equation: where E is the energy per atom in eV, F is the averaged norm of the atomic force in eV/Å, ΔF is the absolute difference of F from the previous step, and w F is the weight parameter for the force term. We fix ΔF ¼ 1 for the first step and ΔF ¼ 10 À 6 if ΔF is less than 10 À 6 . We preferentially proceed with local optimization of candidates with high values of the score L. According to the force term, the score increases as the force applied to each atom increases. Additionally, if the force does not change significantly even if the step changes, that is, if it can be thought that large structural change is continuing, the score increases. Note that the score of fully optimized structures is set to À 1 so that they are not selected. In the example of Figure 1, candidates B and C are preferentially selected by LAQA because of the force term in Equation (1) , thereby enabling an efficient search for the global minimum without unnecessary optimization calculations.
LAQA efficiently works in many cases, however, we have found that it is occasionally not so efficient. As explained in more detail in section 3.1.1, not only force but also stress can be important for the score during the local optimization. In this study, we propose an improved score L to overcome this drawback, expressed as: where S is the average of the absolute values of the components of the stress tensor in eV/Å 3 , and w S is the weight parameter for the stress term. We perform the CSP simulations by LAQA with the improved score using several weight parameters and compare the searching efficiency with that of RS and Greedy (w F ¼ w S ¼ 0). In LAQA, We calculate the first step data of the energy, force, and stress for all the structures. Then, we evaluate the score L to select a candidate for proceeding with the local optimization. The candidates are simply chosen in order of the score L. In this study, the data sets obtained from the CSP simulations by RS were used. Scoring was done for each step, and a candidate with the best score was selected.

Simple example for improvement of LAQA
First, we discuss the importance of the stress term in (2) using a simple example of Si 8 . The CSP simulation by RS was performed, in which 10 initial structures were generated. Here, isotropically shrunk diamond structure (the lattice constant a ¼ 3:0 Å) as shown in Figure 2 drawn by VESTA [32] was intentionally included in the initial structures. This shrunk structure is a good candidate because it is easy to reach the most stable structure by local optimization. However, Greedy and LAQA without the stress term in Equation (2) would not select this good candidate. Figure 3 shows the energy E with respect to that of the most stable structure, force F, and stress S for each optimization step (step data of the E, F, and S). The step data of the shrunk diamond structure are represented in red bold lines in Figure 3. It can be seen that the values of E and S in the initial structure of the shrunk one are quite large because of its small interatomic distance, while the value of F keeps zero owing to its symmetry. This structure is certainly a good candidate, reaching the stable diamond structure in only 10 optimization steps.
We performed the CSP simulation by Greedy and LAQA using the same data set as RS. Here, we tested the values of 0.1, 0.5, and 1.0 for w F , and 0, 10, 100 for w S in LAQA. The reason for the large values of w S is that the values of E and F are roughly comparable, and the values of S are one order of magnitude smaller than those values. Figure 4 demonstrates the number of total optimization steps required to find the most stable structure, obtained by the CSP simulations using RS, Greedy, and LAQA. In RS (blue bar in Figure 4), the required number of steps was evaluated by taking the average of 10,000 times. Greedy (orange bar) and LQAQ without the stress term (green bars) did not work at all. The good candidate was not selected until the end for Greedy and LAQA without the stress term because the score was the  Step data of (a) energy E, (b) force F, and (c) stress S during local optimization. Red bold lines represent the step data of the shrunk diamond structure. lowest due to the high energy and zero force of the initial structure. The improved LAQA overcomes this problem and significantly reduces the total number of steps required to find a stable structure (red and purple bars in Figure 4). The stress term results in a high score for the good candidate, leading to early selections of the good one. It seems that a value of w S around 10 is appropriate for this system, and a value of around 100 tends to overestimate the stress term. Although this system is an extreme case where the force of the good candidate is zero, systems where the force term is small and the stress term becomes important can occur in practice.

Searching efficiency of improved LAQA
We performed test CSP simulations for four typical materials such as Si 16 , Al 4 O 6 , Na 8 Cl 8 , and Sr 4 C 4 O 12 , and used two data sets generated by the CSP simulations by RS for each material. We denote these two data sets as Si 16 -A and Si 16 -B. The most stable structures of Si, Al 2 O 3 , NaCl, and SrCO 3 are the diamond, corundum (α-Al 2 O 3 ), rocksalt, and aragonite-type structure, respectively. In the RS, we randomly prepared 400, 100, 100, 600 candidate structures for Si 16 , Al 4 O 6 , Na 8 Cl 8 , and Sr 4 C 4 O 12 , respectively, as listed in Table 1. The number of most stable structures reached by local optimization using first-principles calculations and the average number of optimization steps are also summarized in Table 1.
To show the effectiveness of the proposed method, we performed CSP simulations on the calculated data sets of Si 16 , Al 4 O 6 , Na 8 Cl 8 , and Sr 4 C 4 O 12 . In the CSP simulations, local optimization was controlled by Greedy and LAQA. As in the case of Si 8 , we tested the values of 0.1, 0.5, and 1.0 for w F , and 0, 10, 100 for w S in LAQA. Figure 5 illustrates the number of total optimization steps required to find the most stable structure, obtained by the CSP simulations using RS, Greedy, and LAQA. Greedy (orange bars in Figure 5) is very efficient on some data sets, however, it is worse than RS for Si 16 -B, Sr 4 C 4 O 12 -B. For LAQA without the stress term (green bars in Figure 5), e.g. in Sr 4 C 4 O 12 -A, its searching efficiency significantly improves Greedy, indicating that the force term is effective. On the other hand, the searching efficiencies are worse than those of RS in some cases for Si 16 -B and Sr 4 C 4 O 12 -B. Therefore, LAQA needs to be improved and the stress term has been introduced. Comparing the results of LAQA with w S of 10 and w S of 100 (red and purple bars, respectively), the searching efficiencies of LAQA with w S of 10 are superior except for the data sets of Sr 4 C 4 O 12 -A and Sr 4 C 4 O 12 -B. w S of 100 tends to overestimate the stress term in many cases. Using a value of 100 for w S often makes the relative contribution of the energy and force terms too small. The value of w S should be around 10, and the results using LAQA with w S of 10 are efficient enough even for the data sets of Sr 4 C 4 O 12 -A and Sr 4 C 4 O 12 -B. For example, in the data set of Al 4 O 6 -B, the introduction of the stress term makes it less efficient, however, the difference is insignificant. Although LAQA with the stress term using w S of 10 is not always the best, it demonstrates high searching efficiency on various data sets. In other words, the introduction of the stress term makes LAQA more robust and versatile.
Finally, we summarized the computational reduction of the improved LAQA in Figure 6. Figure 6 shows the number of total optimization steps required to find the stable structures by improved LQAQ with respect to RS, expressed as a percentage. Here, the values of weight parameters were set to be 0.1 and 10 Figure 4. Number of total optimization steps required to find the most stable structure for Si 8 . The weight parameters of w F ¼ 0:1; 0:5; 1:0, and w S ¼ 0; 10; 100 were used. In RS (blue bar), the required number of steps was evaluated by taking the average of 10,000 times. The required number of steps for Greedy (orange bar) and LQAQ without the stress term (green bars) become large. LAQA with the stress term (red and purple bars) is the most efficient. for w F and w S , respectively, since LAQA with these values was able to search for stable structures efficiently for all of the data sets. The improved LAQA reduced the searching cost to less than 40% in most cases, while the cost was reduced to about 60% even for the least efficient data set (Sr 4 C 4 O 12 -B). The characteristic of the improved LAQA is that it can reliably reduce the computational cost for any data set.

Conclusion
CSP based on first-principles calculations is usually computationally expensive since a lot of candidate structures have to be locally optimized. LAQA, which is one of the selection-type algorithms we previously developed, can control the optimization priority of the candidates. LAQA allows us to reduce the computational costs in many cases, however, it has been found to be inefficient in some data sets. We proposed an improved score of LAQA, where the stress term is added to overcome the drawbacks of the previous algorithm. The CSP simulations by RS, Greedy, and LAQA were performed to investigate the searching efficiency for typical materials such as Si 16 , Al 4 O 6 , Na 8 Cl 8 , and Sr 4 C 4 O 12 . Greedy is very efficient on some data sets, however, it is worse than RS for the other data sets. Although LAQA without the stress term improves Greedy for most cases, it is not always efficient. LAQA with the stress term demonstrates high searching efficiency for all of the data sets. The improved LAQA significantly reduced the searching cost to less than 40% with respect to RS in most cases. The introduction of the stress term makes LAQA more Figure 5. Number of total optimization steps required to find the most stable structure. The weight parameters of w F ¼ 0:1; 0:5; 1:0, and w S ¼ 0; 10; 100 were used. In RS (blue bar), the required number of steps was evaluated by taking the average of 10,000 times. Although LAQA with the stress term using w S of 10 (red bars) is not always the best, it demonstrates high searching efficiency on various data sets and is the most versatile. robust and versatile. Although the candidate structures were randomly generated in this study, LAQA can be applied to candidate structures generated by EA or other algorithms. Developing a combination method of EA and LAQA would be the next step.