Crystal Structure Prediction via Oblivious Local Search
Dmytro Antypov Argyrios Deligkas Vladimir Gusev Matthew J. Rosseinsky Paul G. Spirakis Michail Theofilatos
cs.CE


Quantum-corrected scattering and absorption of a Schwarzschild black hole with GUP
    E. Passos
    March 2020
==================================================================================













 



We study Crystal Structure Prediction, one of the major problems in computational 
chemistry. This is essentially a continuous optimization problem, 
where many different, simple and sophisticated, methods have been proposed and applied. 





The simple searching techniques are easy to understand, usually easy to implement, 
but they can be slow in practice. On the other hand, the 
more sophisticated approaches perform well in general, however almost all of them 
have a large number of parameters that require fine tuning and, in the majority of the
cases, chemical expertise is needed in order to properly set them up. In addition, due 
to the chemical expertise involved in the parameter-tuning, these approaches can be 
biased towards previously-known crystal structures. 

Our contribution is twofold. Firstly, we formalize the Crystal Structure Prediction 
problem, alongside several other intermediate problems, from a theoretical computer
science perspective. Secondly, we propose an oblivious algorithm for Crystal Structure
Prediction that is based on local search. Oblivious means that our algorithm requires 
minimal knowledge about the composition we are trying to compute a crystal structure for. 
In addition, our algorithm can be used as an intermediate step by any method. 
Our experiments show that our algorithms outperform the standard basin hopping, a well studied algorithm for the problem.



Keywords: crystal structure prediction; local search; combinatorial neighborhood. 



§ INTRODUCTION


The discovery of new materials has historically been made by experimental investigation guided by chemical understanding. 
This approach can be both time consuming and challenging because of the large space to be explored. 
For example, a “traditional” method for discovering 
inorganic solid structures relies on knowledge of crystal chemistry coupled with 
repeating synthesis experiments and systematically varying elemental ratios, each
of which can take lots of time <cit.>. As a result 
there is a very large unexplored space of chemical systems: only 72% of 
binary systems, 16% of ternary, and just 0.6% of quaternary systems have been
studied experimentally <cit.>. 

These inefficiencies forced physical scientists to develop computational approaches in order 
to tackle the problem of finding new materials. The first approach is based on 
data mining where only pre-existing knowledge is used <cit.>. 
Although this approach has proven to be successful, there is the underlying risk of 
missing best-in-class materials by being biased towards known crystal structures. 
Hence, the second approach tries to fill this gap and aims at finding new materials 
with little, or no, pre-existing knowledge, by predicting the crystal 
structure of the material. This approach has led to the discovery of several new, 
counterintuitive, materials whose existence could not be deduced by the structures of  
previously-known  materials <cit.>. 

Several heuristic methods have been suggested for crystal structure prediction. 
All these methods are based on the same fundamental principle. 
Every arrangement of ions in the 3-dimensional Euclidean space corresponds to an 
energy value and it defines a point on the potential energy surface.




Then, the crystal structure prediction problem is formulated as a mathematical optimization 
problem where the goal is to compute the structure that corresponds to the global minimum 
of the potential energy surface, since this is the most likely structure that corresponds to a stable
material. The difficulties in solving this optimization problem is that the potential energy
surface is highly non convex, with exponentially many, with respect to the 
number of ions, local minima <cit.>. For this reason, several different algorithmic techniques were 
proposed ranging from simple techniques, like quasi-random sampling
<cit.>, basin hopping 
<cit.>, and simulated annealing 
techniques <cit.>, to more sophisticated techniques, like  
evolutionary and genetic algorithms
<cit.>, and tiling approaches
<cit.>. A recent comprehensive review on these techniques can be found 
in <cit.>.

The simple searching techniques are easy to understand, usually easy to implement, 
and they are unbiased, but they can be slow in practice. On the other hand, the 
more sophisticated approaches perform well in general, however almost all of them 
have a large number of parameters that require fine tuning and, in the majority of the
cases, chemical expertise is needed in order to properly set them up. In addition, due 
to the chemical expertise involved in the parameter-tuning, these approaches can be 
biased towards previously-known structures. 

The majority of the abovementioned heuristic techniques work, at a very high level, in a 
similar way. Given a current solution x for the crystal structure prediction problem, 
i.e., a location for every ion in the 3-dimensional space, they iteratively perform the 
following three steps. 


    
  *  Choose a new potential solution x'. This can be done by taking into 
    account, or modifying, x.
    
  *  Perform gradient descent on the potential energy surface starting from 
    x', until a local minimum is found. This process is called relaxation of x'.

    
  *  Decide whether to keep x as the candidate solution or to update it 
    to the solution found after relaxing x'.


For example, basin hopping algorithms randomly choose x', they relax x' and if the energy
of the relaxed structure is lower than the x, or a Metropolis criterion is satisfied, 
they accept this as a current solution; else
they keep x and they randomly choose x”. The procedure usually stops when the algorithm 
fails to find a structure with lower energy within a predefined number of iterations. The more sophisticated 
algorithms take into account knowledge harvested from chemists and put constraints on the
way x' is selected. For example, the MC-EMMA <cit.> and the FUSE <cit.>
algorithms use a set of building blocks to construct x'. These building blocks are local 
configurations of ions that are present in, or similar to,  known crystal structures. 


These approaches restrict the search space ,which accelerate search, but reduce the 
number of possible solutions.


This general algorithm is easy to understand, however there are some hidden difficulties that make the 
problem more challenging. Firstly, it is not trivial even how to evaluate the potential energy
of a structure. There are several different methods for calculating the energy of a structure, 
ranging from quantum mechanical methods, like density functional theory [<https://en.wikipedia.org/wiki/Density_functional_theory>], 
to force fields methods [<https://en.wikipedia.org/wiki/Force_field_(chemistry)>], 
like the Buckingham-Coulomb potential function. All of which 
though, are hard to compute (see Section <ref>) from the point of view of (theoretical) computer science and thus only
numerical methods are known and used in practice for them <cit.>; still there are cases
where some methods need considerable time to calculate the energy of a structure.

This yields another, 
more important, difficulty, the relaxation of a structure. Since it is hard to compute the  
energy of a structure, it is even harder to apply gradient descent on the potential energy surface. For these 
reasons, the majority of the heuristic algorithms depend on external, well established, codes <cit.> 
for computing the aforementioned quantities.

Put differently, both energy computations and relaxations of structures are treated as oracles
or black boxes.




 §.§ Our contribution.

Our contribution is twofold. Firstly, we formalize the Crystal Structure Prediction problem 
from the theoretical computer science perspective; to the best of our knowledge, this is 
among the few papers that attempt to connect computational chemistry and computer science.
En route to this, we introduce several intermediate open problems from computational chemistry 
in CS terms. Any (partial) positive solution to these questions can significantly help 
computational chemists to identify new materials. On the other hand, any negative result can 
formally explain why the discovery of new materials is a notoriously difficult task.

Our second contribution is the partial answer for some of the questions we cast. 
In general, our goal is to create oblivious algorithms that are easy to implement, 
they are fast, and they work well in practice. With oblivious we mean that we are seeking 
for general procedures  that require minimal input and they have zero, or just 
a few, parameters chosen by the user. 

    
  * We propose a purely combinatorial method for estimating the energy of a  structure, 
    which we term depth energy computation. 
    We choose to compare our method against 
GULP <cit.>, which is considered to be the state of the art for computing the 
energy of a structure and for performing relaxations when the Buckingham-Coulomb energy is used. 
Our method requires only the charges of the atoms and their corresponding Buckingham coefficients 
to work; see Eq. <ref> in Section <ref>. 
In addition, it needs only one parameter, the depth k. 
We experimentally demonstrate that our method monotonically approximates with respect to k 
the energy computed by GULP and that it achieves an error of 0.0032 for k = 6. 

Our experiments show 
that the structure that achieves the minimum energy in depth 1 is likely to be 
the structure with the minimum energy overall. 
In fact we show something much stronger. If the energy of x is lower than the energy of x' 
when it is computed via the depth energy computation for k = 1, then, almost always, 
the energy of x will be lower than the energy of x' when it is computed via GULP.
 
  * We derive oblivious algorithms for choosing which structure to relax next. All
 of our algorithms are based on local search. More formally, starting with x and using only local 
 changes we select x'. We define several “combinatorial neighborhoods” and we evaluate their 
 efficiency. Our neighborhoods are oblivious since they only need access to an oracle that 
 calculates the energy of a structure. We show that our method outperforms basin hopping. 
 Moreover, we view our algorithms as an intermediate step before 
 relaxation that can be applied to any existing algorithm.
 























































§ PRELIMINARIES
 

A crystal is a solid material whose atoms are arranged in a highly ordered 
configuration, forming a crystal structure that extends in all directions.
A crystal structure is characterized by its unit cell; a parallelepiped 
that contains atoms in a specific arrangement. 
The unit cell is the period of the crystal; unit cells are stacked in the 
three dimensional space to form the crystal. 
In this paper we focus on ionically bonded crystals, which we describe next;
what follows is relevant only on crystals of this type.
In order to fully define the unit cell of a ionically bonded crystal structure, we 
have to specify a composition, unit cell parameters, and an arrangement 
of the ions.












  
Composition.
A composition is the chemical formula that describes the ratio of ions that belong 
to the unit cell. The chemical formula contains anions, negatively charged
ions, and cations, positively charged ions.


The chemical formula is a way of presenting information about the chemical
proportions of ions that constitute a particular chemical compound, and it 
does not provide any information about the exact number of atoms in the unit cell. 
More formally, the composition is defined by a set of distinct chemical elements 
{ e_1, e_2, …, e_m}, their multiplicity n_i, and a non-zero integer charge 
q_i for each element i. 
The number m denotes the total number of distinct chemical elements, and 
n_i/∑_j=1^mn_j is the proportion of the atoms of type e_i in the 
unit cell.


It is required that the sum of the charges adds up to zero, i.e. , 
∑_i=1^mq_in_i = 0, so that the unit cell is charge neutral.
For example, the composition for Strontium Titanate, , denotes that the 
following  hold. For every ion of strontium (Sr) in the unit cell, there exists 
one ion of titanium (Ti) and three ions of oxygen (O). Furthermore, the charge 
of  every ion of Sr is +2, of Ti is +4, and of O is -2. Hence, when the ratios 
of the ions are according to the composition, the charge of the unit cell is zero.

Another parameter for every atom is the atomic radius. This usually 
corresponds to the distance from the center of the nucleus to the boundary of 
the surrounding shells of electrons. Since the boundary is not a well-defined 
physical entity, there are various non-equivalent definitions of atomic radius. 
In crystal structures though, the ionic radius is used and usually is 
treated as a hard sphere. Thus, we will use ρ_i to denote the ionic radius 
of the element e_i.






  
Unit cell parameters.
Unit cell parameters provide a formal description of the parallelepiped that represents 
the unit cell. These include the lengths y_1, y_2, y_3 of the
parallelepiped in every dimension and the angles θ_12, θ_13, 
and θ_23 between the corresponding facets. 
For brevity, we denote y = (y_1, y_2, y_3) and  θ = (θ_12, 
θ_13, θ_23), and we use (y,θ) to denote the unit cell parameters.



  
Arrangement.
An arrangement describes the position of each atom of the composition in 
the unit cell. The position of ion i is specified by a point 
x_i = (x_i1, x_i2, x_i3) in the parallelepiped defined by the 
unit cell parameters; fractional coordinates x_i denote the location of 
the nucleus of the ion i in the unit cell. A unit cell parameters-arrangement combination 
(y, θ, x) in a unit cell with n ions is a point in the 
3n+6-dimensional space.





For any two points x_i and x_j we will use (x_i,x_j) to denote 
their Euclidean distance.  



As we have already said, a unit cell parameters-arrangement configuration (y,θ,x)
defines the period of an infinite structure that covers the whole 3d space.
To get some intuition, assume that we have an orthogonal unit cell, i.e., 
all the angles are 90 degrees. Then for every ion with position 
(x_i1,x_i2,x_i3) in the unit cell, there exist “copies” of the ion 
in the positions (k_1 · y_1+x_i1, k_2 · y_2+x_i2, k_3 · y_3+x_i3) 
for every possible combination of integers k_1,k_2, and k_3.
A unit cell parameters-arrangement configuration is feasible if the hard spheres of any two 
ions of the crystal structure do not overlap; formally, it is feasible if for every two
ions i and j it holds that (x_i,x_j) ≥ρ_i + ρ_j.





 §.§ Energy


Any unit cell parameters-arrangement configuration of a composition corresponds to a potential energy. When the number of ions in the unit cell is fixed,  the
set of configurations define the potential energy surface. 



Buckingham-Coulomb potential is among the most well adopted 

methods for computing energy <cit.> and it is
the sum of the Buckingham potential and the Coulomb potential.


The Coulomb potential is long-range and depends only on the charges and the
distance between the ions; for a pair of ions i and j, the Coulomb energy 
is defined by 

    CE(i,j) := q_iq_j/(x_i,x_j).

Note, ions i and j can be in different unit cells.

The Buckingham potential is short-range and depends on the species of the 
ions and their distance. More formally, it depends on  
positive composition-dependent constants A_e_i,e_j, B_e_i,e_j, and
C_e_i,e_j for every pair of species e_i and e_j; here i can be equal 
to j [The Buckingham constants are composition-depended since they can 
have small discrepancies in different compositions. For example the constants 
A_Ti,O, B_Ti,O, and C_Ti,O  for can be different than those for . There is a long 
line of research in computational chemistry that tries to learn/estimate the
Buckingham constants for various compositions. In addition, more than one set of
Buckingham constants can be available for a given composition.]. 
So, for the pair of ions i, of specie e_i, and j, of specie e_j, 

the Buckingham energy is 

    BE(i,j) := A_e_i,e_j·exp(-B_e_i,e_j·(x_i,x_j)) - C_e_i,e_j/(x_i,x_j)^6.

Again, ions i and j can be in  different unit cells.


Let S(x_i,ρ) denote the sphere with centre x_i and radius ρ.
The total energy of a crystal structure whose unit cell is characterized 
by n ions with arrangement x = (x_1, …, x_n) is then defined

    E(y,θ, x) = lim_ρ→∞∑_i=1^n ∑_j ≠ i, j ∈ S(x_i, ρ)( BE(i,j)+ CE(i,j) ).

E(y,θ, x) conditionally converges to a certain value <cit.> and 
usually numerical approaches are used to compute it. For this reason, and since 
we aim for an oblivious algorithm, we view the computation of the energy of a
structure as a black box. More specifically, we assume that we have an 
oracle that given any structure (y,θ,x), it returns its corresponding 
energy. 



Given a composition and Buckingham parameters for it, find a simple, purely
combinatorial way that approximates the energy for every crystal structure.




Given a composition {e_1, e_2, …, e_m} and an oracle that computes the 
energy of every structure for this composition, learn efficiently (with respect to 
the number of oracle calls) the Buckingham parameters A_e_i,e_j, B_e_i,e_j, 
and C_e_i,e_j for every i,j ∈ [m].





  
Relaxation. The relaxation of a crystal structure 
(y,θ,x) computes a stationary point on the potential energy surface 
by applying gradient descent starting from (y,θ,x). The relaxation of 
a structure can change both the arrangement x of the ions in the unit 
cell and the unit cell parameters (y, θ) of the unit cell.

 
We follow a similar approach as we did with the energy and we assume that 
there is an oracle that given a crystal structure (y,θ,x) it returns 
the relaxed structure. 


Find an alternative, quicker, way to compute an approximate local minimum when:



  a) the unit cell parameters (y,θ) of the unit cell are fixed;

  b) the arrangement x of the ions is fixed;

  c) both unit cell parameters and arrangement are free.










 §.§ Crystal Structure Prediction problems

In crystal structure prediction problems the general goal is to minimize the
energy in the unit cell. There are two kinds of problems we are concerned. 
The first cares only about the value of the energy and the second one cares
for the arrangement and the unit cell parameters that achieve the minimum energy.
From the computational chemistry point of view, both questions are interesting
in their own right. The existence of a unit cell parameters-arrangement that achieves 
lower-than-currently-known energy usually suffices for constructing a new material.
On the other hand, identifying the arrangement and the unit cell parameters of a crystal
structure that achieves the lowest possible energy can help physical scientists to predict
the properties of the material.

[title = ]
Input: A composition with its corresponding Buckingham constants, a positive 
integer n, and a rational Ê.

Question: Is there a crystal structure (y,θ,x) for the composition with n ions that is neutrally charged and achieves Buckingham-Coulomb 
energy E(y,θ,x) < Ê?


[title = ]
Input: A composition with its corresponding Buckingham constants, a positive 
integer n.

Task: Find a crystal structure (y, θ, x) for 
the composition with n ions that is neutrally charged and the Buckingham-Coulomb 
energy E(y, θ, x) is minimized.


The second class of problems, the ones that ultimately computational chemists would
like to solve, take as input only the composition and the goal is to construct a unit 
cell, with any number of atoms, such that the average energy per ion is 
minimized. 

[title = ]
Input: A composition with its corresponding Buckingham constants and a 
rational Ê.

Question: Is there a crystal structure for the composition
that is neutrally charged and E(y, θ, x)/n < Ê?


[title = ]
Input: A composition with its corresponding Buckingham constants.

Task: Find a crystal structure for the composition that is
neutrally charged and the average Buckingham-Coulomb energy per ion in the unit
cell, E(y, θ, x)/n, is minimized.


Although the problems are considered to be intractable <cit.>, 
only recently the first correct NP-hardness result was proven for a variant of CSP <cit.>.
However, for the problems presented, there are no correct NP-hardness results in the literature.


Provide provable lower bounds and upper bounds for the four problems defined above.




Construct a heuristic algorithm that works well in practice.





























































































































  §.§.§  Electrostatic energy of C.
 





  
Damping. 

In <cit.> an alternative method, termed damping, was proposed for 
calculating the electrostatic energy. In this method, another parameter, the 
damping parameter R_d, is introduced and furthermore the error 
function, , are used. Recall, 

    (z) = 1 - 2/√(π)∫_0^z e^-t^2dt.
 
Then, the electrostatic energy E_i for ion i is given by 

    E_i := 1/2∑_j ∈ s_ic_i · c_j ·(r_ij/R_d)/r_ij + Δ E_i,
 
where 

    Δ E_i := -π· c_i ·ρ· R_a^2 + π· c_i ·ρ·R_a^2 - R_d^2/2·(R_a/R_c) + √(π)· c_i ·ρ· R_a· R_d · e^-R^2_a/R^2_d - 1/√(π)· R_d· c_i^2.
 
When R_c ≫ R_d the above simplifies to 

    Δ E_i := - π· c_i ·ρ·R_d^2/2 - 1/√(π)· R_d· c_i^2.
 

The electrostatic energy of the coloring is then E = ∑_i E_i. 









§ CHEAP COMPUTATION OF ELECTROSTATIC ENERGY


A very time-consuming step in any CSP optimization algorithm is the calculation of the electrostatic energy.
In general, the electrostatic energy is defined as an infinite sum  of the energies of all pairs of atoms.
Here, we provide a parameterized algorithm that computes the energy up to a predefined depth.
More formally, ...



Even though the potential energy surface might be different for each depth, the problem of finding the global minimum in each depth is equivalent to each other.




Experiments.

Convergence and precision.

Minima (+ global minimum) in depth 1 are minima in depth e.g. 20.









§ LOCAL SEARCH
 


Local search algorithms start from a feasible solution and iteratively obtain better 
solutions. The key concept for the success of such algorithms, is given a feasible
solution, to be able to efficiently find an improved one. 

Put formally, a local search algorithm is defined by a neighbourhood function
N and a local rule r. In every iteration, the algorithm does the following.


  * Has the current best solution x.

  * Computes the neighbourhood N(x).

  * If there is an improved solution x' in N(x), then it updates x according
to the rule r, i.e. x' = r(N(x)); else it terminates and outputs x.



The neighborhood N(x) of a solution x consists of all feasible solutions that are 
“close” in some sense to x. The size of the neighborhood can be constant or a 
function of the input. In principle, the larger the size of the neighborhood, the better 
the quality of the locally optimal solutions. However, the downside of choosing large 
neighborhoods is that, in general, it makes each iteration computationally more expensive.
Running time and quality of solutions are competing considerations, and 
the trade off between them can be determined through experimentation. 

We study the following combinatorial neighborhoods for  Crystal Structure  
Prediction. All of them keep the unit cell parameters fixed and change only the arrangement
x of the n ions. Thus, for notation brevity, we  define the neighbourhoods only 
with respect to the arrangement x. 


    
  * k-ion swap. This neighborhood consists of all feasible arrangements that 
    are produced by swapping the locations of k ions. The size of this neighbourhood is O(n^k k!). 
    
  * k-swap. This neighborhood is parameterized by a discretization step
    δ. Using δ we discretize  the unit cell and then we perform
    swaps of k ions with the content of every point of the discretization. So, an ion can swap
    positions with another ion, or simply move to another vacant position. Again, we take into
    account only the feasible arrangements of the ions. The size of this neighbourhood is 
    O(n^kk!/δ^3k ).
    
  * Axes. This neighborhood has a parameter δ and computes the following for 
    every ion i. Firstly, for every dimension it computes a plane parallel to the  
    corresponding facet of the unit cell and contains the ion i. The intersection of any pair 
    of these planes defines an “axis”. Then, this axis is discretized according to δ.
    The neighborhood locates the ion to every point on the discretization on the three axes and
    we keep only the arrangements that are feasible. The size of this neighbourhood is O(n/δ).











In all of our neighborhoods, we are using a greedy rule to choose x'; x' is an 
arrangement that achieves the minimum energy in N(x). 













In the CSP problem, not all solutions are feasible, that is, there are constraints on the minimum distances between the ions as described in Section <ref>. For this reason, we will let the set of feasible solutions be represented as subsets of a finite set.
We formalize this as follows:


Let E = {1, 2, …, m} be a finite set and F ⊆ 2^E be the set of all feasible solutions. 2^E denotes the set of all the subsets of E.
Consider now an optimization problem P with instances (F, c), where c is the objective (or cost) function.
An objective function maps each feasible solution to a value r ∈^1 (c: F →^1).
In a minimization problem, the goal is to find an f ∈ F, such that c(f) ≤ c(y),  ∀  y ∈ F.


We call the neighborhood of a solution f ∈ F, the set N(f) which contains all the feasible solutions that are “close” in some sense to the solution f. Then,


	N: F → 2^E



defined for each instance f ∈ F. This function associates each feasible solution f to a subset N(f) of E.

We say that a solution f^* ∈ F is locally optimal with respect to a neighborhood function N if c(f^*) ≤ c(f),  ∀  f ∈ N(f^*). 




 §.§ Axis neighborhood


In Section <ref> we presented a combinatorial version of the CSP problem, that is, the ions are placed on a finite (two or three dimensional) grid, while the vertices of the graph are in the Euclidean metric space. For simplicity, in this section we consider two-dimensional grids, but our approaches can be easily applied to any number of dimensions.

Let D=(d_1, d_2) be the size of the grid along the x and y axis respectively. We refer to D as the discretization of our model. We also refer to a position of the grid as (x,y), where x is the x-th row and y is the y-th column of the grid.


We define the axis neighborhood of a feasible solution (or valid configuration) f as being the set of configurations that can be obtained by moving any atom along its current row or column.
More formally,


    N(f) = {g: g ∈ F  and    g can be obtained from f as follows: remove an ion
       α in position (x,y) from f; then place it in a position
       (x',y') where either (1)x'=x and y'=[0,d_2], or
       (2)y'=y and x'=[0,d_1]}


The above neighborhood can be seen graphically in Figure <ref>.
A more generous neighborhood is, for instance, the 2-swap neighborhood, consisting of all the swaps between two ions, or an ion and an empty position.
The downside of it is that the size of a neighborhood becomes O(...), making each iteration computationally more expensive. We experimentally show that the running time is significantly larger, while the quality of the solutions......
In general, a larger neighborhood does not necessarily produce a more effective heuristic unless one can search the larger neighborhood in a very efficient manner.







A very significant design issue is the way that the neighborhood is explored, i.e. the heuristic algorithm.
An exponential to the input size neighborhood does not mean that in each step a heuristic should search all g ∈ N(f) in order to proceed to the next step.
In general, the heuristic algorithm can only search part of the neighborhood, make a local move based on a specific strategy (e.g, greedy strategy), and then move to the next step.
However, one can find efficient ways that allow searching the whole neighborhood without sacrificing the running time.
In addition, randomization can be an invaluable ally in local search. It is typically used in two ways: to pick a random initial solution, for instance a random configuration of ions, and to choose a local move when several are available.


Our simple heuristic uses randomization in order to produce an initial valid configuration f on the given graph G, and an arbitrary ordering O of the atoms.
In each step it proceeds by exploring part of the neighborhood of the current solution, and selects the new solution f based on a greedy strategy, thus, allowing only improvement steps.
In each step, using a round-robin approach, a single atom α is chosen from the ordering O, and then the set N(f, α) ⊂ N(f), where only α does not maintain its position, is explored. We do this in order to reduce the overall running time.
The algorithm terminates once f is a locally optimal solution (with respect to the given neighborhood). In particular, if after searching all N(f, α),  ∀α∈ O, the solution cannot be improved, it terminates. Observe that ⋃_α∈ O N(α, f) = N(f).
In some cases, the fact that a heuristic explores part, instead of the whole neighbourhood in each step, only affects the running time and not the quality of the solution; this is true in our case, because we eventually reach a locally optimal solution, which cannot be improved by searching all g ∈ N(f).







§ CRYSTAL STRUCTURE PREDICTION ALGORITHM


In this section, we present a simple algorithm that uses the simple heuristic described in the previous section, and tries to predict the lowest energy configuration of atoms.
The local neighbourhood that we use has been chosen based on the results of the experiments shown in Section <ref>.
The structure of the algorithm is shown in Fig. <ref>.




 §.§ Combining continuous and combinatorial techniques


Combinatorial algorithms have been proved to solve various combinatorial optimization problems.
However, many global optimization problems, such as the Crystal Structure Prediction, depend on continuous variables.
Hybrid methods, combining discrete and continuous techniques have been extensively used in the literature <cit.> to solve problems in the continuous space. Such algorithms allow to cover widely the solution space using combinatorial optimization techniques, while avoiding the risk of trapping into local minima.
In general, the combinatorial techniques are used in order to “explore” the whole domain and localize “promising areas” likely to contain the global optimum, while continuous techniques are used to “exploit” the already obtained best areas in order to reach the local optimum.
The balance between global and local methods affects the accuracy and time complexity of the algorithm, thus, the tuning of its parameters must be taken into consideration.

Our algorithm follows this approach, employing a gradient descent algorithm, which given any configuration of atoms, finds the local optimum by moving the atoms to the directions determined by the first derivatives until convergence, and the axes neighborhood as described in Section <ref>.
In each iteration of the algorithm, the combinatorial neighborhood of an initially random configuration is explored until a promising area is found, and then, the gradient descent algorithm relaxes that configuration in order to reach the locally optimal solution (or globally optimal).
This process is more accurately described in Figure <ref>.

In general, we cannot know whether a solution is actually optimal or not. In our experiments, we measure the number of iterations needed in order to reach the globally optimal solution, with respect to the number of atoms in the unit cell.
The results (mean value and deviation) can be used as an indication this value.














 §.§ Initialization


In addition to the main stages of the algorithm, neighborhood search and relaxation, there is also an initialization phase.
It is a very important step, as it affects the efficiency of its subsequent steps, as we show in Section <ref>.
Given the chemistry C and the radii of the atoms, the size and volume of the unit cell are determined, based on the experiments in Section <ref>.
In addition, the number of iterations is given as input. In general, the more the number of iterations, the more likely it is that the best configuration of atoms will be found.
Finally, the atoms are placed arbitrarily on the grid, with the only constraint not to overlap with its other.

...




 §.§ Exploration and stopping criteria


The second phase of the algorithm is responsible to localize a promising area of the potential energy surface (i.e. the objective function). For this purpose, the heuristic approach of Section <ref> is being executed.

...

















§ ALGORITHMS



We propose two algorithms. The first one is a step towards answering Open Question <ref> while the second is a heuristic for problem.

For Open Question <ref>, we propose the depth energy computation for estimating the energy of a structure. Our algorithm has a single parameter, the depth parameter k, and works as follows. Given a crystal structure, it creates k layers around the unit cell with copies of the structure. So, for the unit cell parameters (y,θ) and the arrangement x of n atoms the energy is E(y,θ, x) = ∑_i=1^n ∑_j ≠ i, j ∈ D(k)( BE(i,j)+ CE(i,j) ), where D(k) denotes the set of ions in the k layers of unit cells, and BE(i,j) and CE(i,j) are computed as in Equations <ref> and <ref> respectively.

For problem, we slightly modify basin hopping.
In a step of basing hopping, a structure is randomly chosen and it is followed by a relaxation. Our
algorithm applies a combinatorial local search using the Axes neighborhood, since this turned out to
be the best among our heuristics, before the relaxation. 
So, we will perform a relaxation, only after combinatorial local search cannot further improve 
the solution. Our algorithm can be used as a standalone one and it can also be integrated into any other heuristic algorithm for the Crystal Structure Prediction problem since it is oblivious. 
In addition, it provides a very fast criterion that when it succeeds it guarantees finding a lower 
energy crystal structure.





§ EXPERIMENTS

In this section we evaluate our algorithms via experimental simulations. We first focus on which we use as a benchmark. We do this because it is a well studied composition for which the Crystal Structure Prediction problem is solved.

We have implemented the algorithms in Python 2.7 and we use the Atomic Simulation Environment (<https://wiki.fysik.dtu.dk/ase/>) package for setting up, manipulating, running, visualizing and analyzing atomistic simulations.
All experiments were performed on a 4-core Intel i7-4710MQ with 8GB of RAM.









 
        2c|15 atoms          2c20 atoms
 
    k       Time     Energy difference       Time     Energy difference

    1         41          0.0639          73          0.0670    

    2         115         0.0226          174         0.0238    
 
    3         231         0.0114          390         0.0120    

    4         437         0.0068          778         0.0072    

    5         757         0.0045          1380        0.0047    

    6         1218        0.0032          2244        0.0033    

    GULP              20          -                 22          -           


tableComparison between depth approach and GULP for . Time is in milliseconds (ms) 
and energy in electronvolts (eV). The energy difference shows the average difference in energy between 
the depth approach and GULP. Results averaged over 2000 random feasible structures.










 
        6cEnergy difference 
 
    k       1      2     3     4     5     6

    15 atoms       0.0639      0.0226      0.0114      0.0068     0.0045        0.0032    

    20 atoms        0.0670        0.0238       0.0120       0.0072       0.0047       0.0033  


tableComparison between depth approach and GULP for . Energy is in electronvolts (eV). The energy difference shows the average difference in energy between the depth approach and the energy calculated by GULP. Results averaged over 2000 random feasible structures.






We evaluate the depth energy computation in several different dimensions. For all the experiments we performed for energy computation, we fixed the unit cell to be cubic. Firstly, we evaluate how depth energy computation behaves with respect to k. We see that the method converges very fast and k=6 already achieves accuracy of three decimal points. Then, we compare our depth approach against GULP; see Table <ref>.

Our goal is to provide an intuitively simpler to interpret and work with method for computing the energy. Even though the energy calculated by the depth approach differs from the one calculated by GULP, we observe that the relative energies between two random arrangements remain usually the same even for k=1. To be more precise, let E_1(x) denote the energy of a feasible arrangement x when k=1 and let E_G(x) denote the energy of this arrangement as it is computed by GULP. Our experiments show that if for two random feasible arrangements x_1 and x_2 it holds that E_1(x_1) <  E_1(x_2), then E_G(x_1) <  E_G(x_2) for 99.8% of 1000 pairs of arrangements. This percentage reaches 100 % for k=6. For the “special” arrangement of ions x^* that minimizes the energy computed by GULP, that is x^* =  E_G(x), our experiments show that it is always true that E_k(x^*) < E_k(x), for every k=1, …, 6, where x is a random feasible structure over 10000 of them. So, this is a good indication that the arrangement that minimizes the energy for k=1, also minimizes the energy overall. We view this as a striking result; it significantly simplifies the problem thus new, analytical, methods can be derived for the problem.






 
    3|c|Ordering of energy between pairs of random structures 

        15 atoms          20 atoms
 
        Success rate         Success rate
 
    Depth 1         98.9 %         99.3 %    

    Depth 2         99.8 %         99.8 %    
 
    Depth 3         99.8 %         99.8 %    

    Depth 4         99.8 %         99.8 %    

    Depth 5         99.9 %         99.8 %    

    Depth 6         100 %          99.9 %    

    
table. Results averaged over 1000 pairs of random structures.









 
    2|c|Ordering of energy between random structures and global optimum 

        Success rate 
 
    Depth 1         99.99 %    

    Depth 2         100.0 %    
 
    Depth 3         100.0 %    

    
table(5 atoms). Results averaged over 10000 random arrangements.








The next set of experiments compares the three neighborhoods described in 
Section <ref> for [The values of the Buckingham parameters can be found in the appendix]. We compare them in several different dimensions: the average 
CPU time they need in order to find a local optimum with respect to their combinatorial 
neighborhood and the average drop in energy until they reach such a local optimum 
(Tables <ref> and <ref>); the average CPU time the relaxation needs 
starting from such local minimum and the average drop in energy from relaxation 
(Tables <ref> and <ref>). In addition, for the case of , we compare how often we can find 
the optimal arrangement from a single structure. 
We observe that the Axes neighborhood has the best tradeoff between energy drop and CPU time.
The 2-ion-swap neighbourhood outperforms the other two in terms of running time, however it seems to decrease the probability of finding the best arrangement when performing a relaxation on the resulting structures. This renders the use of 2-ion-swap neighbourhood inappropriate.
Axes neighborhood is significantly faster and performs smoother in terms of running time than the 2-swap neighbourhood. However the latter one performs better with respect to the energy drop, which is expected since axes is a subset of the 2-swap neighbourhood. In addition, the relaxation from the local minimum found by 2-swap significantly improves the probability of finding the best arrangement with only one relaxation.




We should highlight that there exist structures where the relaxation cannot improve 
their energy, but the neighborhoods do; hence using Axes neighborhood we can escape from some local minima of the continuous space.




 
    Neighbourhood         Running time          Time stdev        Energy drop       Energy drop stdev   
 
    Axes              5.36        1.54        13.46       10.60 
 
    2-ion swap      0.96        0.33        7.75        8.45 

    2-swap          34.66       14.06        16.21       10.94 


tableComparison of local neighbourhoods for reaching a combinatorial minimum for with 15 atoms per unit cell and δ = 1 (375 grid points). 
Time is in seconds and energy in electronvolts (eV). Results averaged over 1000 arrangements.








 
    Neighbourhood      Running time     Time 
 stdev     Energy drop     Energy drop 
 stdev     Global 
 minimum
 
    Random structures-GULP        8.80        6.35        18.60       10.26        6.6% 

    Axes-GULP              7.92        6.16        5.53       2.28       10.0% 
 
    2-ion swap-GULP        8.82        6.25        11.09        5.64        4.7% 

    2-swap-GULP          5.14       5.08        2.79       1.05        14.8%


tableEvaluation of relaxation procedure
after using a combinatorial neighborhood for with 15 atoms per unit cell. 
Time is in seconds and energy in electronvolts (eV). Results averaged over 1000 arragnements.










 
    Neighbourhood         Running time          Time stdev        Energy drop       Energy drop stdev   
 
    Axes              7.56        2.08        8.23        3.71
 
    2-ion swap      0.88        0.43        1.16        1.80 

    2-swap          27.93       9.98        10.26       4.05


tableComparison of local neighbourhoods for reaching a combinatorial minimum for and δ = 1 (343 grid points). 
Time is in seconds and energy in electronvolts (eV). Results averaged over 1000 arrangements.










 
    Neighbourhood      Running time     Time 
 stdev     Energy drop     Energy drop 
 stdev 
 
    Random structures-GULP      2.81     1.45     12.84     4.88 

    Axes-GULP                   2.30     1.20     5.09      1.81 
 
    2-ion swap-GULP           2.47     2.25     12.29     4.38 

    2-swap-GULP               1.99     1.09     3.11      0.93 


tableEvaluation of relaxation procedure
after using a combinatorial neighborhood for . 
Time is in seconds and energy in electronvolts (eV). Results averaged over 1000 arragnements.











Next, we compare our algorithm for against basin hopping where the next structure to relax is chosen at random. Based on the results of our previous experiments, we have chosen the Axes neighbourhood as an intermediate step before the relaxation. We have run these algorithms 200 times for with 15 atoms per unit cell, and 25 times for with 20 atoms per unit cell.
We report how the energy varies with respect to time until the best arrangement is found (Fig. <ref>) and we report other statistics that further validate our approach (Table <ref>).

As we can see, it is relatively easy to reach low levels of energy and the majority of time is needed to find the absolute minimum.
In addition, the overhead posed by the use of the neighbourhood search divided by the time needed by the basin hopping to find the global minimum decreases as the number of the atoms in the unit cell increases.





 
    Algorithm      Total time 
 mean     Total time 
 median     Total time 
 stdev     Relaxations     Time for 
 relaxations     Time for 
 local search     Time for 
 energy calculations
 
    Axes-GULP           83.64        57.25     85.80      9.12        65.89       17.75        10.17 

    Basin hopping       169.74        112.51      168.58      18.45       169.73       -     -


tableComparison of relaxed structures for with 15 atoms. Time is in seconds and energy in electronvolts (eV). Results averaged over 500 arragnements.











 
    Algorithm      Total time 
 mean     Total time 
 stdev     Relaxations     Time for 
 relaxations     Time for 
 local search 
 
    Axes-GULP           227.89      287.21      13.24        126.26       101.63 

    Basin hopping       167.89      114.89      18.14       160.79       - 


tableStatistics from the experiments depicted in Figure <ref> (with 15 atoms per unit cell).









 
    Algorithm      Total time 
 mean     Total time 
 stdev     Relaxations     Time for 
 relaxations     Time for 
 local search 
 
    Axes-GULP           2280.57      781.66      104.33        1016.72       1049.13 

    Basin hopping       5766.20      4748.33      450.66       4895.60       - 


tableStatistics from the experiments depicted in Figure <ref> (with 20 atoms per unit cell).












 
    Algorithm      Number of 
 atoms      Total time 
 mean     Total time 
 stdev     Relaxations     Time for 
 relaxations     Time for 
 local search 
 
    Axes-GULP       15      227.89      287.21      13.24        126.26       101.63 

                    20      2280.57      781.66      104.33        1016.72       1049.13 

    Basin hopping      15      167.89      114.89      18.14       160.79       - 

                       20      5766.20      4748.33      450.66       4895.60       - 


tableStatistics from the experiments depicted in Figures <ref> and <ref> (with 15 and 20 atoms per unit cell).














In our last set of experiments, we compare our algorithm against basin hopping algorithm for which contains 22 atoms in its primitive unit cell.
In this set of experiments, we produced 3500 random structures. We simulated basin hopping by sequentially relaxing the constructed structures. However, none of the relaxations managed to find the optimal configuration.  Our algorithm, using the same order of structures as before, first used the Axes neighborhood as an intermediate step followed by a relaxation; it managed to find the optimal configuration after visiting only 720 structures (Fig. <ref>).




























§ CONCLUSIONS

In this paper we have introduced and studied the Crystal Structure Prediction problem through the lens
of computer science. This is an important and very exciting problem in computational chemistry, 
which computer scientists are not actively studying  yet. We have identified several open questions 
whose solution  would have significant impact to the discovery of new materials. These problems are 
challenging and several different techniques and machineries from computer science could be applied 
for solving them. Our simple-to-understand algorithms are a first step towards their solution. We hope 
that our algorithms will be used as benchmarks in the future, since more sophisticated techniques for 
basin hopping can be invented. For the energy computation via the depth approach, we conjecture that 
the arrangement that minimizes the energy for k=1 or k=2, matches the arrangement that minimizes the 
energy when it is computed via GULP. Our numerical simulations provide significant evidence towards this. 
A formal result of this would greatly simplify the objective function of the optimization problem and it 
would give more hope to faster methods for relaxation. In addition, it could provide the foundations for
new techniques for crystal structure prediction. Our algorithm that utilizes the Axes neighborhood as an 
intermediate step before relaxation, seems to speed up the time the standard basin hopping needs to find 
the global minimum. Are there any other neighborhoods that outperform the Axes one? Can local search, or
Axes neighborhood in particular, improve existing methods for crystal structure prediction by a simple
integration as an intermediate step? We believe that this is indeed the case.

















abbrv






§ APPENDIX














 
    Interaction      A (eV)     ρ  ()     C (eV ^-6) 

    O^2- - O^2-           1388.77      0.36262      175 

    Y^3+ - O^2-       23000      0.24203      0

    Sr^2+ - O^2-       1952.39      0.33685      19.22

    Ti^4+ - O^2-       4590.7279      0.261      0


tableThe values of Buckingham potential parameters we used in our experiments as they were found in <cit.>. All the missing parameters are set to zero.




