XtalOpt: Evolutionary Crystal Structure Prediction

• XtalOpt is an evolutionary algorithm for crystal structure prediction, using stochastic structure generation and physics-based relaxations to discover low-enthalpy phases.
• It employs genetic operators such as lattice mutation, atomic-position mutation, and crossover to maintain diversity and optimize candidate structures.
• Recent versions incorporate multi-objective and Pareto optimization strategies, enabling variable-composition searches and enhanced performance with ML-based relaxations.

Searching arXiv for the cited XtalOpt papers and closely related versions. XtalOpt is an evolutionary algorithm for crystal structure prediction that maintains a population of candidate crystal structures and iteratively improves that population through local relaxation, fitness evaluation, parent selection, and genetic variation. In the versioned literature summarized here, it is described as a platform for global optimization of crystalline materials under fixed composition, structurally constrained searches for metastable phases, multi-objective search for functional materials, and, in a later release, variable-composition Pareto optimization (Hajinazar et al., 2024). Across these formulations, the recurring goal is to discover low-enthalpy or otherwise targeted crystal structures by combining stochastic structure generation with physics-based local optimization and evolutionary selection.

1. Core evolutionary search formalism

In its population-based form, XtalOpt evolves candidate structures generation by generation. A generation consists of spawning new structures, locally relaxing every structure, evaluating fitness or filtering structures, selecting parents by roulette-wheel sampling weighted by fitness, applying genetic variation operators to produce offspring, and then replacing or augmenting the population with offspring while respecting constraints (Hajinazar et al., 2024). The search terminates when the user-specified maximum number of structures has been generated, when no fitness improvement is observed over a given number of generations, or when a hard or soft exit flag is triggered (Hajinazar et al., 2024).

Each candidate structure is represented by a triplet consisting of lattice vectors, fractional atomic coordinates, and atom types. Internally, the lattice is stored as a $3\times 3$ matrix $A$, and atomic positions are stored as fractional vectors $\{s_i\}$, so that Cartesian coordinates are given by $x_i = A \cdot s_i$ (Hajinazar et al., 2024). In the fixed-composition setting, composition is constant throughout a run, although an extended grand-canonical version is noted elsewhere in the version-13 description (Hajinazar et al., 2024).

The standard evolutionary workflow described for metastable-phase discovery follows the same structure. Initialization can be purely random or seed-assisted; every trial structure is geometry-optimized; duplicates are removed with XtalComp; low-enthalpy structures form a breeding pool; and standard operators such as lattice mutation, atomic-coordinate mutation, crossover, and chemical transmutation generate children for subsequent generations (Wang et al., 2023). In that formulation, the top $N_{\rm pool}$ lowest-enthalpy structures are retained as parents, with $N_{\rm pool}=50$ in practice (Wang et al., 2023).

This continuity across releases suggests that XtalOpt’s identity is defined less by a single fitness formalism than by a stable evolutionary search engine onto which different objective schemes and constraints can be grafted.

2. Structure generation, relaxation, and diversity control

Initial structures are generated randomly within user-defined bounds. In version 13, the initial population is random within user-defined cell-volume and composition limits (Hajinazar et al., 2024). In the carbon superhardness study, random structures were generated using the RandSpg algorithm, with the number of formula units per primitive cell chosen by the user; for carbon, searches were performed with 8, 12, 16, and 20 atoms per cell, and random cell volumes were restricted to lie near that of diamond in order to bias the search toward realistic densities (Avery et al., 2019).

Local optimization is delegated to external engines. Version 13 states that every structure is locally relaxed via an external optimizer such as VASP or an interatomic potential (Hajinazar et al., 2024). In the 2019 superhard-materials study, every new structure was sent to VASP for a full geometry and cell-shape relaxation using PBE-GGA, PAW, and a 500 eV cutoff; $\Gamma$-centered k-meshes were chosen so that lattice constant times k-point count was approximately $25\,\text{\AA}$ in the final evolutionary-algorithm step, with tighter settings used for phonons and final energies, and convergence was imposed when all forces were below $1\times 10^{-5}\,\text{eV/\AA}$ (Avery et al., 2019). In the structurally constrained study, DFT relaxations were performed in three steps—internal coordinates with fixed cell, then volume-only, then full relaxation—and precise re-optimizations used cutoffs between 600 and 1000 eV depending on the system (Wang et al., 2023).

Diversity control is handled explicitly. In the carbon study, after relaxation each child was compared to the existing pool via the XtalComp duplicate-checking routine, and only genuinely new minima were retained (Avery et al., 2019). The structurally constrained study likewise states that XtalComp identifies and removes duplicates from the population (Wang et al., 2023). Duplicate removal functions as a mechanism for preventing the evolutionary search from repeatedly sampling symmetry-equivalent or near-identical minima.

3. Fitness functions and multi-objective optimization

For single-objective search, XtalOpt uses enthalpy minimization. In version 13, the fitness of structure $s$ is

$A$0

with enthalpy defined as

$A$1

where $A$2 is the internal energy from the local optimizer, $A$3 is the external pressure, and $A$4 is the unit-cell volume computed from the lattice matrix (Hajinazar et al., 2024). Pressure is either taken from the DFT run for variable-cell relaxations or set to the user-specified external pressure (Hajinazar et al., 2024).

Version 13 generalizes this to a decomposition-based multi-objective evolutionary search. If $A$5 are objectives to be minimized and $A$6 are objectives to be maximized, and each objective is normalized between current population extrema, the generalized fitness is

$A$7

Enthalpy is always included as the first objective in $A$8 unless the user explicitly assigns it zero weight, so that the search remains anchored to low-enthalpy, (meta)stable structures while allowing properties such as band gaps, hardness, or elastic moduli to be optimized in parallel (Hajinazar et al., 2024). The same source explicitly distinguishes this weighted linear scalarization from Pareto-based schemes such as NSGA-II: version 13 does not perform non-dominated sorting or crowding-distance selection, and no explicit Pareto front is constructed (Hajinazar et al., 2024).

A concrete earlier realization of multi-objective XtalOpt appears in the carbon study, where each relaxed individual acquires both its $A$9 enthalpy $\{s_i\}$0 and an estimated Vickers hardness $\{s_i\}$1, and the combined fitness is

$\{s_i\}$2

with $\{s_i\}$3 and $\{s_i\}$4 chosen so that $\{s_i\}$5 (Avery et al., 2019). There, roulette-wheel selection is applied after sorting the pool by descending $\{s_i\}$6 and renormalizing over the retained population (Avery et al., 2019).

Version 14 marks a methodological transition. It introduces Pareto optimization for variable-composition search and uses distance above the convex hull, $\{s_i\}$7, rather than raw enthalpy as the energetic stability objective (Hajinazar et al., 7 Jun 2025). In Basic mode, objectives are still normalized and combined through a scalar fitness, but in Pareto mode the raw objectives, including $\{s_i\}$8, enter non-dominated sorting directly and weights are not used (Hajinazar et al., 7 Jun 2025). The version-14 description therefore positions Pareto optimization as a replacement for weighted-sum-only optimization in that release (Hajinazar et al., 7 Jun 2025).

4. Genetic operators and constrained search mechanisms

The principal variation operators retained in version 13 are lattice mutation, atomic-position mutation, cut-and-splice crossover, and optional symmetry-based moves (Hajinazar et al., 2024). Lattice mutation applies a small random symmetric strain $\{s_i\}$9 to the lattice via $x_i = A \cdot s_i$0, with $x_i = A \cdot s_i$1 (Hajinazar et al., 2024). Atomic-position mutation displaces each fractional coordinate by up to $x_i = A \cdot s_i$2 in each direction and folds positions back into $x_i = A \cdot s_i$3 (Hajinazar et al., 2024). Cut-and-splice crossover selects a random plane in fractional coordinates, splits two parent crystals along that plane, and reassembles halves from each parent before rescaling to avoid atomic overlaps (Hajinazar et al., 2024). Symmetry-based moves can raise or lower the space-group symmetry by small perturbations followed by re-symmetrization (Hajinazar et al., 2024).

The carbon superhardness implementation describes crossover and single-parent mutation in less formal but compatible terms, including shuffling atoms and small cell-shape perturbations (Avery et al., 2019). Version 14 broadens the operator set to crossover, stripple, permustrain, permutomic, permucomp, and optional random supercell mutation, with operator probabilities derived from adjustable integer weights (Hajinazar et al., 7 Jun 2025). In fixed-composition or multi-composition searches, operators that would change composition are skipped; in variable-composition mode, composition can change through crossover of parents with different stoichiometries, atom addition or removal, and random generation of new compositions (Hajinazar et al., 7 Jun 2025).

Constrained search was developed explicitly in the structurally constrained evolutionary algorithm. Rather than continuously mixing energetic and structural metrics, each structure is assigned

$x_i = A \cdot s_i$4

with $x_i = A \cdot s_i$5 per formula unit, ensuring that failing structures never enter the top-$x_i = A \cdot s_i$6 breeding pool (Wang et al., 2023). Constraint checks occur after DFT relaxation and before selection (Wang et al., 2023).

The structural constraints supported in that work include coordination-number filters, local chemical-environment requirements, Bravais-lattice filtering, and space-group filtering, all implemented through Pymatgen analysis (Wang et al., 2023). Coordination numbers are defined by a user-chosen cutoff $x_i = A \cdot s_i$7, with examples such as $x_i = A \cdot s_i$8 for N–N or C–N bonds (Wang et al., 2023). Bravais lattice is determined with symmetry-analysis tolerance 0.02, and the space-group number can be required to match a target value (Wang et al., 2023).

Version 13 incorporates a related but more general “filtration” objective. A user may specify objective = fil /path/filter.sh [out](https://www.emergentmind.com/topics/outer-automorphism-out).dat; the filter script reads output.POSCAR or other files and writes exactly one integer, 1 for keep or 0 for discard (Hajinazar et al., 2024). Structures for which the script returns 0 are removed from the parent pool entirely, though they remain in the overall archive (Hajinazar et al., 2024). Parents are drawn only from the filtered pool, and any child that fails the filter is immediately rejected or replaced according to the user’s policy, prior to any relaxation or objective evaluation (Hajinazar et al., 2024).

A plausible implication is that XtalOpt’s constrained-search capability evolved from hard structural screening tightly coupled to enthalpy ranking toward a more general framework in which arbitrary user scripts can enforce admissibility conditions at the parent-pool level.

5. Software architecture, interfaces, and execution model

Version 13 describes the software architecture as a C++ core with Qt for the GUI, supplemented by a Python/C shell for job scripting and for user-provided objective scripts (Hajinazar et al., 2024). Relaxations and objective evaluations are submitted asynchronously either to a local queue or to a remote scheduler such as PBS or SLURM; the code tracks job submission IDs, monitors completion, and parses outputs (Hajinazar et al., 2024). The main modules listed are the Search Engine, Scheduler Interface, Objective Manager, and IO Manager (Hajinazar et al., 2024).

The Objective Manager is responsible for launching user scripts, parsing single-value outputs, applying weights, and handling filtering (Hajinazar et al., 2024). The IO layer reads and writes input flags, POSCAR outputs, structure.state logs, and results.txt live ranking files (Hajinazar et al., 2024). Important command-line input entries include structure limits, composition, optimizer type, population size, maximum number of structures, queue settings, and a list of objectives with associated scripts and weights (Hajinazar et al., 2024). The output includes one column per objective plus fitness in results.txt, per-structure output.POSCAR and state files, and various global logs (Hajinazar et al., 2024).

Version 14 extends the execution model by supporting local optimizers that can read and write VASP POSCAR/CONTCAR plus OUTCAR, including VASP, CP2K, and Quantum ESPRESSO, and by providing a Python wrapper, vasp_uip.py, to connect ML potentials such as MACE or CHGNet through ASE for rapid relaxations (Hajinazar et al., 7 Jun 2025). It also introduces explicit support for fixed-composition, multi-composition, and fully variable-composition search modes (Hajinazar et al., 7 Jun 2025).

The asynchronous job-scheduling model in version 13 is described as enabling scaling to hundreds of concurrent relaxations on modern HPC clusters (Hajinazar et al., 2024). That statement is paired with the observation that lightweight objective scripts contribute to efficient and robust high-throughput crystal-structure prediction workflows (Hajinazar et al., 2024). This suggests that the software is intended not only as an evolutionary search engine in the abstract, but also as an HPC-oriented orchestration layer for large numbers of external relaxation and property-evaluation jobs.

6. Applications, benchmarks, and methodological trajectory

The 2019 superhard-carbon study provides an early demonstration of XtalOpt as a multi-objective search tool directed toward low-enthalpy, superhard materials. The workflow combined VASP relaxations, AFLOW-ML estimates of shear modulus, the Teter hardness model $x_i = A \cdot s_i$9, and a combined energy-hardness fitness (Avery et al., 2019). Four searches on carbon cells with 8, 12, 16, and 20 atoms produced 5,624 relaxed minima (Avery et al., 2019). Screening criteria of $N_{\rm pool}$0 and $N_{\rm pool}$1 yielded 827 superhard, low-energy structures out of 1,324 low-energy structures; removing duplicates left 89 distinct minima, of which phonon calculations found 79 dynamically stable, including 43 newly predicted superhard phases (Avery et al., 2019). The paper also reports a diamond reference hardness of $N_{\rm pool}$2 and two newly predicted polytypes, $N_{\rm pool}$3-16 and $N_{\rm pool}$4-16, with $N_{\rm pool}$5 (Avery et al., 2019).

The structurally constrained study benchmarks targeted metastable-phase discovery across XeN$N_{\rm pool}$6, TiO$N_{\rm pool}$7, BaH$N_{\rm pool}$8, and WC$N_{\rm pool}$9N$N_{\rm pool}=50$0 (Wang et al., 2023). In XeN$N_{\rm pool}=50$1 at 50 GPa, a coordination-constrained run discarded 57.4% of candidates and located the desired $N_{\rm pool}=50$2–II 2D-polymeric phase at approximately the 250th structure, whereas an unconstrained run of approximately 1500 relaxations did not recover that phase (Wang et al., 2023). In TiO$N_{\rm pool}=50$3, an orthorhombic Bravais-lattice constraint identified brookite as the 345th structure, compared with the 1141st structure in the unconstrained run (Wang et al., 2023). In BaH$N_{\rm pool}=50$4 at 50 GPa, a space-group-constrained run found the correct $N_{\rm pool}=50$5 phase after about 65 relaxations, compared with a ranking of 354th in the unconstrained search (Wang et al., 2023). In WC$N_{\rm pool}=50$6N$N_{\rm pool}=50$7, a chemical-environment-constrained run discovered a new salt-like melaminate $N_{\rm pool}=50$8 phase as the 515th structure, 53 meV/atom above the $N_{\rm pool}=50$9 ground state and 29 meV/atom below a previously proposed $\Gamma$0 phase in D3(BJ) results; phonon spectra confirmed dynamical stability (Wang et al., 2023).

Version 13 itself does not report a new performance benchmark in the manuscript, but it cites earlier applications of XtalOpt’s multi-objective evolutionary search and states that searches with tens to hundreds of structures per generation typically converge to low-enthalpy or metastable candidates within a few hundred to one thousand total structures (Hajinazar et al., 2024). Version 14, in turn, reports in-house tests on the XtalOpt website showing that variable-composition plus Pareto searches locate new ternary oxide and nitride phases in fewer generations than fixed-composition or multi-composition alone, and that an ML-potential wrapper can speed local relaxations by 10–100× versus DFT (Hajinazar et al., 7 Jun 2025).

Taken together, these studies trace a methodological trajectory: initial enthalpy-driven evolutionary search, addition of explicit property-driven scalarized fitness, introduction of hard structural filters for metastable targets, general multi-objective weighted linear scalarization in version 13, and finally variable-composition Pareto optimization around distance above the convex hull in version 14 (Hajinazar et al., 2024). A common misconception is that all XtalOpt multi-objective searches are Pareto-based; the version history shows instead that weighted-sum scalarization and Pareto optimization correspond to distinct stages of the codebase, with the former characterizing version 13 and earlier property-coupled studies, and the latter introduced explicitly in version 14 (Hajinazar et al., 2024).