Evolutionary algorithm (EA)

2025 April 2

Background

Evolutionary algorithms (EAs) are metaheuristic methods inspired by the theory of evolution. EA can effectively generate new structures (offspring) by inheriting the local environments of the stable structures (parents) explored so far. Oganov group’s USPEX is a well-known software, and there are others such as XtalOpt. For example, there are previous studies such as the following papers.

  • T. S. Bush, C. R. A. Catlow, and P. D. Battle, J. Mater. Chem. 5, 1269 (1995).
  • A. R. Oganov and C. W. Glass, J. Chem. Phys. 124, 244704 (2006).
  • A. R. Oganov, A. O. Lyakhov, and M. Valle, Acc. Chem. Res. 44, 227 (2011).
  • A. O. Lyakhov, A. R. Oganov, H. T. Stokes, and Q. Zhu, Comput. Phys. Commun. 184, 1172 (2013).

fig_EA_pes fig_EA_pes

Procedure

  1. Initialize population
  2. Evaluate fitness
  3. Natural selection
  4. Select parents
  5. Create next generation
  6. Repeat from step 2: Evaluate fitness

Initialize population

In the first generation, a set of random structures is generated according to the number specified by n_pop. tot_struc is not used in EA or EA-vc.

Evaluate fitness

Currently, energy is the only property that can be used as fitness in CrySPY. By setting fit_reverse = False, the algorithm is configured to search for the minimum value. The fit_reverse setting is designed for future cases where fitness may be based on properties other than energy.

Natural selection

DFT calculations occasionally fail and produce extremely unreasonable energy values. emax_ea and emin_ea can be used to filter out structures with unreasonably low (or high) energy values:    $$ \mathrm{emin\_ea} \le E \ (\mathrm{eV/atom}) \le \mathrm{emax\_ea} $$ For example, if emin_ea is set, any structure with an energy lower than that value will be ignored.

In natural selection, the current population and elite individuals preserved from previous generations are first ranked based on fitness. The number of elite individuals used here is specified by n_elite. Only the top n_fittest individuals among the current population and elite individuals are selected, while all others are eliminated. During the natural selection process, duplicates are removed using the StructureMatcher class provided by pymatgen, and then the top n_fittest individuals are selected. n_fittest is often set to about half of n_pop (the population size). Note that in the current implementation, if n_fittest = 0, all individuals are retained. The figure below shows an example of the natural selection process when n_fittest = 5.

fig_EA_natural_selection fig_EA_natural_selection

Select parents

Two parent selection methods are implemented in CrySPY to select a single parent individual from the candidate parents. Both methods are designed so that individuals with higher fitness have a higher probability of being selected. Setting slct_func = TNM enables tournament selection, while slct_func = RLT enables roulette selection. Tournament selection requires fewer parameters and is easier to use.

Create next generation

The next generation consists of offspring produced by evolutionary operations on candidate parents, along with some randomly generated structures. Random structures are added in each generation to maintain diversity and to help escape from local minima.

Evolutionary operations

Here, we introduce the operations of the fixed-composition EA implemented in CrySPY.

Population size

The sum of structures from crossover, permutation, strain, and random generation must be equal to n_pop.

  • n_pop = n_crsov + n_perm + n_strain + n_rand

Subsections of Evolutionary algorithm (EA)

Crossover

Overview

Crossover is an evolutionary operation that creates a new structure (offspring) by exchanging sliced regions between two parent structures. This promotes structural diversity and enables the inheritance of locally stable features. It is one of the main operators used to explore low-energy configurations in structure search.

How it works

  1. Select two distinct individuals from the candidate parents
  2. Perform a random translation
  3. Randomly select a lattice vector
  4. Slice the parents near the center
  5. Swap the sliced halves
  6. Select the offspring with more atoms
  7. Adjust the number of atoms near the border
  8. Perform a minimum interatomic distance check

fig_EA_crossover fig_EA_crossover

4. Slice the parents near the center

The slice point is placed near the center and slightly varied each time.

slice_point = np.clip(np.random.normal(loc=0.5, scale=0.1), 0.3, 0.7)

If any of the subsequent steps fail, the process may be restarted from step 4. However, the number of retries is limited to maxcnt_ea, and if this limit is exceeded, the parent selection step is repeated.

5. Swap the sliced halves

When crs_lat is set to random, the lattice vectors of one of the two parent structures are randomly selected. When crs_lat is set to equal, the average of the lattice vectors of the two parent structures is used. The default is random.

6. Select the offspring with more atoms

Swapping the sliced parts of the parent structures results in two structures with different numbers of atoms. Temporarily, we select the structure with more atoms.

fig_EA_crossover_natoms fig_EA_crossover_natoms

However, if the composition differs too much from the target, the process restarts from step 4 (Slice the parents near the center). The tolerance for the difference in the number of atoms is set by nat_diff_tole. The default value of nat_diff_tole is 4, which allows a tolerance of ±4 atoms per element. In the figure above, the number of blue atoms is -1 and the number of green atoms is +2 relative to the original composition.

7. Adjust the number of atoms near the border

Deletion

When adjusting the number of atoms, the process starts with atom deletion. The number of green atoms is excessive and needs to be reduced. As illustrated in the figure below, atoms that do not satisfy the minimum interatomic distance defined by mindist are preferentially removed.

fig_EA_crossover_rm_mindist fig_EA_crossover_rm_mindist

As shown below, if atoms that violate the minimum interatomic distance remain after the deletion process, the procedure is restarted from step 4 (Slice the parents near the center).

fig_EA_crossover_rm_mindist2 fig_EA_crossover_rm_mindist2

If there are no atoms violating the minimum interatomic distance but atoms still need to be deleted, atoms are removed in order of increasing distance from the border, as shown in the figure below. Note that, in addition to the central slicing point, positions with internal coordinates of 0 are also considered borders.

fig_EA_crossover_rm_border fig_EA_crossover_rm_border

Addition

When atoms are lacking, the missing atoms are added near the border. The internal coordinate along the selected axis is determined as shown below. Here, mean refers to either the slice point or 0.

coords[axis] = np.random.normal(loc=mean, scale=0.08)

The remaining two components of the coordinate are determined randomly. Atoms are added until the target number is reached, while checking for violations of the minimum interatomic distance.

fig_EA_crossover_addition fig_EA_crossover_addition

Permutation

Overview

Permutation is an evolutionary operation that generates new structures (offspring) by modifying the atomic arrangement within a single structure. It enables the exploration of alternative configurations without changing the lattice or the overall composition.

How it works

The positions of atoms of different elements are swapped. The number of swaps can be specified by ntimes, which is set to 1 by default. After the swap, a minimum interatomic distance check is performed.

fig_EA_permutation fig_EA_permutation

Strain

Overview

Strain is an evolutionary operation that generates a new structure (offspring) by applying a small random distortion to the lattice of a parent structure. It helps to explore nearby regions of the configuration space while preserving atomic connectivity and composition. This operator is useful for fine-tuning structural candidates and escaping local minima.

How it works

The lattice vectors are $ \mathbf{a} $ transformed to $ \mathbf{a}' $ by applying a strain matrix, as follows:

$$ \mathbf{a}' = \begin{pmatrix} 1 + \eta_1 & \frac{1}{2} \eta_6 & \frac{1}{2} \eta_5 \\ \frac{1}{2} \eta_6 & 1 + \eta_2 & \frac{1}{2} \eta_4 \\ \frac{1}{2} \eta_5 & \frac{1}{2} \eta_4 & 1 + \eta_3 \end{pmatrix} \mathbf{a}. $$

Here, $ \eta_i $ are given by a Gaussinan distribution $ \mathcal{N}\left( 0, \ \sigma_{\mathrm{st}}^2 \right) $. $ \sigma_{\mathrm{st}} $ is specified by the input parameter sigma_st (by default, sigma_st = 0.5). As shown in the figure below, the lattice is deformed and then rescaled to restore the original volume. Finally, the minimum interatomic distance constraint is checked.

fig_EA_strain fig_EA_strain

Tournament selection

Overview

Tournament selection is a method used to choose parent individuals from the candidate parents based on their fitness. It is designed to balance selection pressure and diversity in the population. The figure below shows an example with n_fittest = 10 and t_size = 3.

fig_EA_tournament fig_EA_tournament

How it works

  1. A fixed number of individuals (t_size) are randomly selected from the candidate parents.
  2. Among them, the individual with the highest fitness (i.e., lowest energy) is chosen as the parent.
  3. This process is repeated until the required number of parents is selected.

Advantages

  • Simple and efficient
  • Requires only one parameter (t_size)
  • Can control selection pressure by adjusting t_size

Notes

  • The default value of t_size is 3.
  • If t_size is small, diversity is promoted.
  • If t_size is large, selection pressure increases, favoring the fittest individuals.
  • Unlike roulette selection, tournament selection never chooses the bottom (t_size - 1 ) individuals from the candidate parents.

Roulette selection

Overview

Roulette selection is a probabilistic method used to select parent individuals from the candidate parents based on their fitness. In roulette selection, each individual’s chance of being selected is proportional to its fitness.

How it works

  1. When fit_reverse is set to False (default), corresponding to minimization mode where energy is used as fitness, the fitness values of the candidate parents are multiplied by –1.
  2. The fitness values $ f_i $ are linearly scaled into $ f'_i $ using the following equation, where $ a $ and $ b $ are parameters specified by a_rlt and b_rlt, respectively (with the condition that $ a > b $). $$ f_i' = \frac{a - b}{f_{\mathrm{max}} - f_{\mathrm{min}}} f_i + \frac{b f_{\mathrm{max}} - a f_{\mathrm{min}}}{f_{\mathrm{max}} - f_{\mathrm{min}}} $$
  3. The scaled fitness values $ f_i’ $ are converted into selection probabilities using the following equation:    $$ p_i = \frac{f_i’}{\sum_k f_k’} $$ Each probability $ p_i $ represents the likelihood of selecting the $ i $-th individual.
  4. Parent individuals are then selected one by one according to the probabilities $ p_i $ using roulette wheel sampling, until the required number of parents is obtained.

Advantages

  • All individuals have a non-zero chance of being selected
  • Selection pressure can be adjusted by scaling the fitness values

Notes

  • By default, a_rlt = 10.0 and b_rlt = 1.0
  • Proper scaling of fitness values is important to ensure meaningful selection pressure.The figure below shows examples of $ p_i $ when $ a $ is relatively small (left) and relatively large (right). If $ a $ is too small, the selection pressure becomes weak, making it more difficult to favor individuals with higher fitness.

fig_EA_roulette fig_EA_roulette