The CALYPSO methodology for structure prediction*

Qunchao Tong(童群超),Jian Lv(呂健),Pengyue Gao(高朋越),and Yanchao Wang(王彥超)

Innovation Center of Computational Physics Methods and Software,State Key Laboratory of Superhard Materials,College of Physics,Jilin University,Changchun 130012,China

Keywords:structure prediction,CALYPSO method,crystal structure,potential energy surface

1.Introduction

The atomic structure is the basis for a deep understanding of the properties or functionalities of materials and is closely relevant to many areas of science as diverse as physics,chemistry,biology,pharmaceutics,etc.Predicting structures utilizing modern theory and supercomputers is highly desirable for many reasons.On the one hand,experimental measurements of structures are usually indirect or under extreme conditions,leading to poor or incomplete data from which the complete structure determination is impossible. For example,photoelectron spectroscopy or collision cross section measurements are usually used to study the structures of atomic clusters at the nanoscale,but only indirect data,such as the electron binding energy and cross section,can be extracted.[1,2]Powder x-ray diffraction combined with a diamond anvil cell apparatus is a major tool for solving the crystal structures of materials at high pressures,but the small sample size generally induces poor diffraction data.[3]In these cases,theoretical structure predictions can provide putative structure models to supplement experiments for a complete structure determination.On the other hand,traditional methods of functional material discovery require expensive and laborious experimental attempts involving long-period trial-and-error processes. In this context,structure prediction combined with quantum mechanical simulation can be much faster and cheaper than experiments to search a wide range of systems for discovering promising new materials,which will provide guidance for purposive experimental synthesis.

Structure prediction with given information of chemical composition can be mathematically formulated as a global optimization problem,since the most stable structure generally corresponds to the global free energy(reduced to enthalpy at 0 K)minimum of the potential energy surface(PES).In principle,the exact solution can be obtained by visiting all local minima on the PES to determine the global one.However,the number of local minima on the PES is expected to grow exponentially with the number of atoms in the system.This growth renders the exhaustively ergodic searching strategy unfeasible,even for systems with a few atoms,and the issues can be infinitely exacerbated with an increasing system size.Therefore,structure prediction has been long thought to be a formidable problem.[4]

Over the years,much effort has been devoted to solving the problem of structure prediction.[5,6]Many algorithms have been proposed either to reduce the searching space or to enhance the sampling efficiency of the PES;these approaches have led to steady progress in the structure prediction community. A number of powerful structure prediction methods have been proposed,where typical examples are simulated annealing,[7]basin hopping,[8]minima hopping,[9]metadynamics,[10]random sampling,[11]genetic algorithm,[12–17]and our developed CALYPSO(Crystal structure AnaLYsis by Particle Swarm Optimization)method.[18,19]These methods significantly push the upper limit of the system size that is tractable for structure prediction and endow computational scientists with the predictive power to guide experiments.More and increasingly more theory-driven discoveries of new materials have been made,[20–24]such as the sulfur hydrides[25,26]and lanthanum superhydrides[27,28]with near room-temperature superconductivity,[29–31]the sodium chlorides with unexpected stoichiometries[32]and the Cu2Si monolayer with planar hexacoordinate bonding.[33]The structure prediction method has become a standard tool in essentially all fields concerning the arrangement of atoms.There are considerable studies in the literature that review the methodology of different structure prediction methods and their applications in different areas.[5,6,11,34,35]In this work,we only focus on our CALYPSO method,paying particular attention to its basic theory and design principles.We refer the readers to Refs.[36]and[37]and review papers in the same issue for the method’s wide applications.

This review is organized as follows. In Section 2,we provide a brief introduction on the general feature of the PES and the problem of structure prediction.In Section 3,we review the basic theory and design principles of the CALYPSO method,as well as the various functional modules of the CALYPSO software. Finally,the current challenges and future developments of the structure prediction method are discussed.

2.Potential energy surfaces and theoretical structure prediction

In the context of theoretical structure prediction,the PES represents the energy as a function of geometric structures of assemblies of atoms(Fig.1). The dimension of PES increases linearly with the number of atoms N within the system,which are 3N+6 and 3N ?3 for periodic crystal and isolated molecules,respectively. Theoretical and numerical analyses have demonstrated that the number of local minima on the PES grows exponentially as the number of atoms increases.[38,39]For example,the number of local minima is approximately 103for the 13-atom Lennard–Jones(LJ)cluster,and it increases to at least 1012for the 55-atom cluster.[40]There is some consensus on the general features of the PES,as described in Refs.[5],[9],[11],and[41],which we briefly restate as follows:

(i)The PES is composed of basins of attraction on which local optimization leads to an energy minimum that corresponds to a dynamically stable structure.

(ii)The substantial fraction of the PES corresponds to structures in which some atoms are very close to each other.These structures generally have energies even much higher than that of a set of isolated atoms,and the corresponding high-energy regions of the PES contain almost no minima.

Fig.1.Schematic illustration of a one-dimensional potential energy surface.

(iii) According to the Bell–Evans–Polanyi principle,[42,43]it is more likely to find a low energy local minimum if one moves from the current basin over a low barrier into a new basin than if one is required to overcome a high barrier.As a consequence,the low-energy basins on the PES tend to gather together,forming a superbasin,which is normally called a funnel.The structures that reside in the same funnel are generally similar to each other and can be regarded as corresponding to the same structural motif.

(iv)The probability distribution of the energies of the local minima of a PES is close to Gaussian when the system is sufficiently large.Symmetric structures tend to correspond to very low or very high-energy minima.[44]

The goal of structure prediction is to find the global minimum of a PES as well as the corresponding structure.The exponential increase in the number of local minima indicates that it belongs to a notorious NP-hard(non-deterministic polynomial-time hard)problem. Therefore,any algorithm cannot guarantee working without failing in polynomial time.All effective algorithms would be bound to have an upper limit in the system size that can be tractable,and they may eventually converge when averaged over all the systems.[45]Although the prospect for completely solving the problem of structure prediction seems bleak,significant progress has been made in the community during the last few decades,which make structure prediction a popular tool routinely used in modern computational physics,chemistry,and material science.This is partially because knowledge about the general features of the PES allows one to disregard the most irrelevant part or bias the desired part of the PES when developing a structure prediction method.This significantly simplifies the problem and makes structure prediction methods applicable to systems with moderate sizes. For example,the CALYPSO method can reliably address systems with several tens of atoms when the energy calculation is taken at the quantum mechanical level,and it can address more than one hundred atoms when empirical potentials are used.A large number of structures have been predicted and experimentally confirmed,validating the method.Some of them possess a large system size,e.g. the 64-atom BC3[46]and 76-atom Bi4Si3O12.[47]Note that,the majority of inorganic compounds tend to have small system sizes.Approximately 90%of crystals in the Inorganic Crystal Structure Database have less than 60 atoms in the primitive cell.[48]This number indicates that state-of-theart structure prediction methods are already applicable to most ordered inorganic compounds.

3.The CALYPSO method

CALYPSO is a numerical method for st ructure prediction with given information with respect to the chemical composition alone for a material.[18,19]The essence of the method is a heuristic global minimization strategy for the PES,which is composed of several structure dealing methods,such as,structure generation with constraints of symmetry,structure characterization through the bond characterization matrix(BCM),and structure evolution via the swarm-intelligence based particle swarm optimization(PSO)algorithm.

Structure predictions through CALYPSO consist of mainly four steps,as depicted in the flow chart in Fig.2.First,the initial structures are randomly generated with the constraints of symmetry to allow a diverse sampling of the PES.Then,BCM will be used to characterize the newly generated structures and examine their distances/similarities with all the previous ones.Structures with distances less than a threshold will be eliminated.After a user-specified number of structures(a population or generation)have been generated,local structural optimizations are performed to eliminate the noise of the energy surface and drive the structures to the local minima.Eventually,the swarm-intelligence PSO algorithm is applied to produce new structures for the next generation.This process continues iteratively until convergence criteria are reached.

The various structure dealing methods within CALYPSO were designed by taking into account the following three aspects for solving the problem of structure prediction:(i)simplifying the problem using physical constraints,(ii)extracting subproblems through structure characterization,and(iii)solving subproblems via a heuristic global optimization algorithm.Below we introduce the various structure dealing methods within CALYPSO from the above three aspects.

Fig.2.The flow char of CALYPSO.Various structure dealing methods are colored according to their role in solving the problem of structure prediction.

3.1.Simplifying the problem using physical constraints

Knowledge regarding the PES,as stated in Section 2,allows us to impose various physical constraints on structures to simplify the problem of structure prediction. The most straightforward one is the constraint of the minimum interatomic distance,which has been adopted by most structure prediction methods.From feature ii of the PES,we known that a substantial fraction of the PES possesses high energy and contains almost no minima,which corresponds to structures with some atoms very close proximity to each other.Therefore,imposing the constraint of the minimum interatomic distance on structures during the structure prediction can disregard a large fraction of high-energy regions of the PES.This constraint will significantly reduce the search space and simplify the problem. In the CALYPSO method,we used the constraint of the minimum interatomic distance when generating new structures.Note that there is no risk at losing the global minimum for imposing this constraint,as we know that all real structures show interatomic distances within a reasonable range(even distances between H atoms under high compression are not less than ～0.74 ?A).

The second constraint on structures is symmetry. The CALYPSO structure search involves the generation of a large number of random structures to sample the PES.Feature iv of the PES tells us that fully random structures correspond to the local minima with a Gaussian-like energetic distribution.This has been numerically demonstrated using isolated LJ systems in our previous work,[49]where we randomly generated 10000 structures for LJ38and LJ100clusters,receptively,and then locally optimized them.As shown in Fig.3,both clusters demonstrate a Gaussian-like distribution with a peak centered at a high energy.As the number of atoms increases from 38 to 100,the peaks shift to higher energies and become apparently narrower. A careful inspection of the structures around the peak revealed a large number of disordered or liquid-like structures without symmetry.Therefore,although intuitively a fully random generation of structures is expected to give a diverse sampling of the PES,it in fact frequently produces nonsymmetric liquid-like structures with high energies,which is unfavorable for structure prediction.Given that symmetric structures tend to correspond to a very low-or a very high-energy minimum(Feature iv),it is beneficial to impose symmetry during the generation of random structures,which would achieve a more diverse sampling of the PES in terms of energies and symmetries.For example,when we randomly impose C1to C6point group symmetry for the generation of the LJ38and LJ100clusters(Fig.3),the energetic distributions are more spread,and more low-energy local minima have been detected.In this respect,biasing symmetric structures during random structure generation is more sensible and“random”.

Fig.3. Energetic distributions of randomly generated structures with and without symmetry constraints for the(a)LJ38 and(b)LJ100 clusters.Energies are shown relative to the global minimum.The LJ potential for a pair of atoms is given by Uij=ε[(r/rij)12?2(r/rij)6],where ε and r are the pair equilibrium well depth and separation,respectively.The reduced units,i.e.,ε=r=1 are employed throughout.Reprinted from Ref.[49],with the permission of AIP Publishing.

In CALYPSO,different symmetry groups are adopted for generating structures with different dimensions. In particular,230 space group symmetries are used to constrain bulk crystals,and 17 in-plane space group symmetries are used to constrain two-dimensional(2D)systems,such as layered structures and surfaces. In principle,an infinite number of point groups exist for zero-dimensional(0D)isolated systems(molecules or nanoclusters).Here,we adopt 48 point group symmetries that frequently appear in such systems. For example,when we generate a crystal structure,one of the 230 space groups will be randomly selected,then the lattice parameters are generated within the chosen symmetry according to a confined volume,and the atomic positions are obtained by combining several sets of symmetry-equivalent coordinates(Wyckoff Positions)in accordance with the symmetry and the number of atoms in the simulation cell.We have demonstrated that the inclusion of symmetry constraints during structure generation significantly increases the efficiency and success rate of structure prediction.[19,49]

The third constraint on the structures is the dynamical stability,which is a prerequisite for the real structures and can be readily achieved by local optimizing a structure to its corresponding local minimum(Feature i of the PES).Each newly generated structure is subject to local optimization during the CALYPSO structure search. Although this step is the most time-consuming part(since it involves a large number of energy calculations)during structure prediction,it plays a role as transforming the searching space from an entire continuous PES to a series of discrete local minima,which significantly simplifies the problem.Local structural optimization has been an indispensable step in modern structure prediction methods and is largely responsible for their success.

3.2.Extracting subproblems through structure characterization

If a structure prediction problem can be split into several subproblems,then each subproblem should be easier to solve than the whole problem.Feature iii of the PES tell us that the basins of attraction arrange themselves into funnels,each of which represents a class of structures with a certain structural motif.Therefore,it is reasonable to split the problem on the basis of the funnels.This approach requires a method to characterize structures and provide a metric for measuring the distance(similarity/dissimilarity)between structures in the configuration space.Structure characterization is widely involved in many areas,such as structure prediction,[18,50–52]machine learning for constructing interatomic potentials or materials discovery,[53,54]and chemical informatics.[55]There are different requirements for structure characterization in different applications.In the context of structure prediction,a real valued function is required to uniquely characterize a structure solely from the geometry and chemical identities of the constituent atoms.The function should be invariant to the translation and rotation of a structure,as well as the permutation of atoms of the same type within the structure. Moreover,a reasonable metric for measuring distances between structures should be defined based on the function,i.e.,the distances should be zero between the identical structures and gradually increase as structures depart from each other.

Based on the above requirements,we designed a BCM to characterize structures in the CALYPSO method.[19,49]This matrix is realized by extending the bond-orientational order parameters,which was previously introduced by Steinhardt et al.[56]The BCM is constructed based on all bond information of a structure.This method utilizes spherical harmonics and exponential functions to describe the orientations and lengths of bonds.A bond vectorbetween atoms i and j is defined if the interatomic distance between them is less than a cutoff distance.The vectoris associated with the spherical harmonicswhereandare the polar angles.A weighted average over all bonds formed by types A and B atoms is then performed by the equation:

where δABand NδABdenote the type and number of bonds,respectively. Only even-l spherical harmonics are used in the above equation to guarantee the invariant bond information with respect to the direction of the bonds. To avoid dependence on the choice of the reference frame,it is important to consider the rotationally invariant combinations,

where u and v denote two individual structures and Ntypeis the number of bond types.Since all Qlexcept for Q are zero for isotropic systems,e.g.,infinitely large liquid,it is therefore plausible to define a distance from a structure to the isotropic system as

which can be used to estimate the degree of order for a given structure. A numerical demonstration of the validity of the BCM for structure characterization has been shown for isolated clusters in Ref.[49]and for crystals in Ref.[19].

Although we are able to uniquely quantify structures through the BCM,the unambiguous definition of the funnels of a PES is still difficult since it requires a detailed understanding of the PES.Therefore,we designed a scheme to dynamically decide the center of funnels during the CALYPSO structure search.Given a set of structures that have been detected during a structure search,the lowest-energy structure is always selected as the center of a funnel.Then,the second one is defined as the lowest-energy structure among the structures that have BCM distances(relative to previous defined centers)larger than the threshold value.This procedure repeats until a predefined number of centers have been chosen.As the structure search proceeds,new structures are continually included at each generation,and the centers of the funnels will be updated when a new structure with lower energy is found.The threshold value of distance will also be adjusted to allow a sufficient number of centers to be defined.Based on the defined centers,CALYPSO can readily classify a structure into a funnel which it is nearest to(without the need of any prior information regarding the PES).Therefore,the problem of structure prediction can be divided into several funnels,each of which can be solved through a heuristic global optimization.

3.3.Solving subproblems via a heuristic global optimization algorithm

The CALYPSO method is within the evolutionary scheme where structures evolve according to the PSO algorithm.The PSO algorithm is a typical swarm-intelligence scheme for global optimization inspired by natural biological systems(e.g.,ants,bees,or birds),[57,58]and has been applied to a variety of fields in engineering and chemical science.In practice,a candidate structure in the configurational space is regarded as a particle,and a set of individual particles is called a population or a generation.Each particle explores the search space via a velocity vector,which is affected by both its personal best experience and the best position found by the population so far(Fig.4(a)).

Within the PSO scheme,the position of ith particle at the jth dimension is updated according to the following equation:

where t denotes the generation index and j ∈{1,2,3}.The new velocityis calculated on the basis of its previous location,previous velocitythe personal best locationwith an achieved best fitness of this individual and the population best locationwith the best fitness for the entire population,

where ω(in the range of 0.9–0.4)denotes the inertia weight,c1=2,c2=2,r1and r2are two random numbers that are uniformly distributed in the range[0,1].

We have implemented two versions of the PSO algorithm(global and local version)in CALYPSO,which are illustrated in Figs.4(b)and 4(c),respectively. In the global version,the problem is solved as a whole though one PSO evolution.All particles seek new positions according to one best location(gbest)for the entire population and their personal best locations(pbest). It has been demonstrated to be efficient and powerful for global structural convergence,especially for small systems.However,it is noteworthy that for large systems with much more complex PES,finer structural searches are desirable.Therefore,we introduced a local version of PSO in which the whole problem is solved though several PSO evolutions.[49,59]As stated in the above section,we assume that each particle belongs to a funnel based on the BCM distances.Then,its velocity is adjusted according to both its personal best position and the best position(lbest)achieved so far within the funnel,

Definitely,lbest denotes the center of the funnel to which a particle belongs.Therefore,the entire swam diffuses into several funnels,in which several PSO evolutions are performed simultaneously. The local PSO algorithm can be seen as a combination of several information-shared global PSOs.This method apparently avoids premature structures by maintaining multiple attractors and allows the exploration of the larger space of PESs,although it is computationally more demanding.

Fig.4.Schematics of(a)the velocity and position updates in PSO,(b)global PSO,and(c)local PSO.Reprinted from Ref.[37],with the permission of IOP Publishing.

Fig.5.(a)Energy as a function of the evolution step of a CALYPSO structure search performed on the LJ100 cluster.(b)Energy–distance distribution for the minima searched during the CALYPSO run.The energies are shown relative to the global minimum and the BCM distances were calculated with respect to the isotropic system.Symbols denote the four lbests at the 200th generation,which are shown in the right panel.The numbers below each structure are the energies in ε relative to the global minimum and BCM distances with respect to the isotropic system.See text for detailed structural descriptions.Reprinted from Ref.[49],with the permission of AIP Publishing.

We illustrate the local PSO through a CALYPSO structure search of the LJ100cluster. Figure 5(a)depicts the energies as a function of the evolution step of the structure search,where four lbests were used.In the first 50 evolution steps,the overall energies decreased rapidly.Then,structural searches were mainly conducted on the low-energy regions,and the global stable structure was successfully produced at the 163rd generation. We plot the energy–distance distribution for the minima detected during this structure search in Fig.5(b).Structures possessing similar structural motifs show similar BCM distances to the isotropic system,and these structures are expected to reside in the same energy funnel.Different funnels can be clearly identified by the calculated BCM distances in Fig.5(b).The symbols denote the four lbest structures determined at the 200th generation(the corresponding structures are shown in the right panel in Fig.5).These structures reside in different funnels and possess structural motifs that are very different from each other.The first lbest structure is the global stable structure with an icosahedral structural motif and a BCM distance of 0.16137.The second lbest structure contains a Marks decahedral motif and has a BCM distance of 0.36857.The third and fourth lbests consist of Marks decahedral and Mackay icosahedral cores with incomplete anti-Mackay overlayers,respectively.This example clearly shows that local PSO enables a simultaneous search in different energy funnels of a potential energy surface.

3.4.Other strategies for biasing low-energy structures and increasing the structure diversity

We have introduced the basic idea of the CALYPSO structure search.However,several other efficient techniques have also been implemented in CALYPSO to further increase the search efficiency. To maintain the structural diversity as well as further biasing the structure search toward low-energy regions,we have included a penalty function during the structure search.In each generation,a certain percentage of lowenergy structures are used to produce the next generation by PSO,while the remaining high-energy structures are rejected and replead by newly generated random structures.This technique has been demonstrated to be important for enriching the structure diversity and is crucial for increasing the search efficiency.Moreover,inspired by the Monte Carlo simulation,the Metropolis criterion can also be included in CALYPSO for acceptance or rejection of new structures.A structure will be accepted if it has a lower energy than its parent structure.Otherwise,a selection probability is imposed based on the relative energies according to the Boltzmann distribution.This implementation further improves the possibility of generating low-energy structures during structural evolution.

3.5.Features of the CALYPSO software

The CALYPSO method as stated above has been implemented into the same-name software package. Currently,it contains ten functional modules,which allows to perform unbiased search for the energetically stable/metastable structures of isolated clusters/nanoparticles,[49]2D layered materials,[60]surfaces,[61]interfaces,[62]adsorbate systems,[63]proteins,[64]and three-dimensional(3D)crystals.[18]It can also be used to design novel functional materials with desirable functionalities when we change the global fitness function from the total energy to a specific property(e.g.,hardness,bandgaps,and even experimental x-ray diffraction patterns).[65–67]Moreover,we have recently proposed a simple and unambiguous definition for crystal structure prototypes based on hierarchical clustering and constructed a crystal structure prototype database(CSPD)by filtering existing crystallographic structure databases.[68]The CSPD also accumulates a low-energy structure during CALYPSO structure prediction to construct a separate theoretical structure database.This database can be used for generating initial structures for structure prediction or determining the prototype of a structure,as well as highthroughput calculation.

4.Conclusion and perspectives

We have presented a short review on the basic theory of our developed CALYPSO method,paying particular attention to the design principles of the various structure dealing methods.The validity of the CALYPSO method has been extensively demonstrated through its application to various materials ranging from 0D to 3D systems.We refer the readers to review papers in the same issue for their latest applications on 0D nanoclusters,2D layered materials,and 3D bulk crystals for superconductors and superhard materials,as well as exotic materials at high pressure.

Long thought to be impossible,modern theoretical structure prediction has achieved substantial progress and can routinely deal with systems with several tens of atoms. However,large systems containing hundreds of atoms is still a great challenge for structure prediction,where a large number of structure-related problems exist concerning realistic materials(e.g.,alloys,minerals,and drugs).Both searching and ranking problems should be further addressed. Quantum mechanical methods for energy calculations can provide a reasonable ranking of structures,but the computational cost is very demanding for large systems,while empirical potential suffers from poor transferability. One promising solution is the utilization of machine learning potentials,[69]which show accuracy comparable to those of quantum mechanical methods,but with less computational cost by several orders of magnitude. The utilization of machine learning potentials to accelerate structure predictions is now in its early stages but has already generated several encouraging results.[70–72]Recently,we have developed an acceleration scheme for the CALYPSO structure prediction of large systems,in which a machine-learning potential(Gaussian approximation potential[73])is trained in an on-the-fly manner from scratch during the structure searches and is used to accelerate the structure searches.[74]The scheme substantially reduced the computational cost by at least ～1–2 orders of magnitude compared with full DFT-based structure searches.For the search problem,it is still difficult to give a definitive answer as to whether we can cross the“exponential wall”of the number of local minima.Further simplifying the PES for a specific problem is a practical way.For example,the partial or complete fixation of the geometry of molecules is a wise choice for predicting molecular crystals.[75]In addition to large systems,several other important issues are also faced by structure prediction,such as predicting materials at finite temperature and the synthesizability of a material.With these challenges addressed,we imagine more wide applications of the structure prediction method in all fields relevant to atomic structures.