GAtor: A First Principles Genetic Algorithm for Molecular Crystal Structure Prediction

1 Introduction

• Polymorphs can differ in stability, solubility, compressibility, and optoelectronic properties, affecting pharmaceutical and device performance.Molecular crystals are bound by weak dispersion interactions, and the same molecule may crystallize in several solid forms.Source: §1 · pp. 1–2
• Crystal structure prediction seeks putative crystal structures from a molecule’s two-dimensional chemical diagram, but requires accurate electronic structure methods and efficient configuration-space exploration.CSP can reveal molecular behavior, predict new polymorphs, and complement experiments.Source: §1 · p. 2
• Polymorph energy differences are typically within a few kJ/mol, motivating quantum-mechanical methods and dispersion-inclusive DFT for stability ranking.Vibrational and finite-temperature contributions can also affect some relative rankings.Source: §1 · p. 2
• GAtor is a modular, massively parallel first-principles genetic algorithm for (semi-)rigid molecular crystals, with customizable fitness, selection, crossover, and mutation schemes.It interfaces with FHI-aims for dispersion-inclusive DFT evaluations and local optimizations; parallel replicas share a population without waiting for a full generation.Source: §1 · p. 3

2 DFT Settings

• The GA search locally optimizes structures with PBE+TS using light numerical settings and tier 1 basis sets to accelerate geometry relaxations.Space-group symmetry may vary; the calculations use a 2 × 2 × 2 k-point grid and reduced angular grids, with SCF force and stress evaluations omitted.Source: §2 · pp. 3–4
• The best 5–10% of final GA structures are re-relaxed and re-ranked with higher-level PBE+TS settings, then evaluated by PBE+MBD and finally PBE0+MBD.The authors consider the PBE0+MBD ranking the most reliable among the methods used.Source: §2 · p. 4
• Thermal contributions are excluded from the study’s final energy ranking, although they have changed rankings for approximately 9% of organic compounds.The search settings are intended to save computational time, while postprocessing provides more precise rankings.Source: §2 · pp. 3–4

3 Code Description

• GAtor’s core tasks—including selection, similarity checks, crossover, and mutation—can be interchanged or modified, while energy evaluations and local optimizations currently use FHI-aims.The implementation uses spglib, scikit-learn, and pymatgen.Source: §3 · p. 4
• Independent GA replicas run across computing nodes and read from and write to a dynamically updated common structure pool.Additional multiprocessing may be used within each replica for child generation.Source: §3 · p. 4, Figure 2 · p. 4
• After selection, crossover or mutation generates offspring that must pass geometry, uniqueness, and energy checks before entering the common population.The population fitness is then updated, allowing another selection to occur immediately.Source: §3 · p. 5

3.1 GA Initialization

• Genarris supplies the initial pool with random symmetric structures, using fragment-based DFT and machine-learning clustering to produce a diverse starting population.The structures are pre-relaxed with PBE+TS and their total energies are stored before GAtor updates their starting fitness values.Source: §3.1 · p. 5

3.2.1 Energy-based Fitness

• Energy-based fitness assigns higher values to structures with lower relative energies, using the population’s dynamically updated highest and lowest total energies.The energies are evaluated with dispersion-inclusive DFT.Source: §3.2.1 · p. 5

3.2.2 Cluster-Based Fitness

• Traditional energy-based selection may repeatedly explore the same landscape regions, as population, packing-motif, or energy-method biases can contribute to genetic drift.A desirable packing motif may lie in a narrow, rarely visited well.Source: §3.2.2 · pp. 5–6
• GAtor dynamically clusters structures by similarity using RDF, RCD, or lattice-parameter descriptors and offers K-Means and Affinity Propagation clustering.Affinity Propagation determines the number of clusters from a similarity matrix rather than requiring it in advance.Source: §3.2.2 · p. 6
• Fitness sharing scales structures according to the size of their shared cluster, increasing fitness for under-sampled low-energy motifs and suppressing densely populated regions.The authors identify further investigation of descriptor and fitness-function effects as future work.Source: §3.2.2 · p. 6

3.3 Selection

• GAtor offers roulette-wheel and tournament selection, giving higher-fitness structures a greater probability of mating.These are the two selection strategies currently available in GAtor.Source: §3.3 · p. 7

3.3.1 Roulette wheel selection

• Roulette-wheel selection chooses a parent using a uniform random value and the sorted, normalized fitness values.A structure is selected when its fitness value exceeds the random value.Source: §3.3.1 · p. 7

3.3.2 Tournament Selection

• Tournament selection chooses the winner and runner-up from a user-sized random group, controlling selection pressure without sorting the population.The user-defined tournament size controls selection pressure.Source: §3.3.2 · p. 7

3.4 Crossover

• Molecular-crystal crossover operators avoid relying only on cut-and-splice, which can break space-group symmetries associated with packing and lower energies.GAtor’s operators can preserve or break selected parent symmetries.Source: §3.4 · p. 7

3.4.1 Standard Crossover

• Standard crossover combines parent lattice parameters, molecular orientations, geometries, and center-of-mass positions to generate offspring.Molecular geometries can relax during optimization, so a relaxed parent geometry is inherited for semirigid molecules.Source: §3.4.1 · pp. 7–8
• In the standard-crossover example, mixed parent genes produce a Z′=2 child with P¯1 symmetry before optimization.The child has lower symmetry than either parent but retains inversion symmetry.Source: §3.4.1 · p. 8

3.4.2 Symmetric Crossover

• Symmetric crossover combines asymmetric-unit genes and applies one parent’s space-group operations, so the child inherits that parent’s space group before optimization.The procedure accounts for differences in asymmetric-unit molecule counts and lattice types when selecting compatible symmetry operations.Source: §3.4.2 · pp. 8–9

3.5 Mutation

• Mutation operators modify one parent to explore the energy surface or exploit promising regions, with user control over mutation rates and operator choice.GAtor can also apply two mutation operations in combination, an approach that encourages exploration and may reduce duplicate offspring.Source: §3.5 · p. 9

3.5.1 Strains

• Strain mutations alter lattice vectors using selected strain components, with options including random, directional, volume-preserving, and single-angle changes.Random strain strength is controlled by a user-defined standard deviation; molecular centers of mass move according to their fractional coordinates.Source: §3.5.1 · p. 9

3.5.2 Molecular Rotations

• Molecular rotations reorient selected molecules; applying the same rotation about each molecule’s principal axes preserves P21/c symmetry in the illustrated Z′=1 case.The illustrated rotation changes the asymmetric-unit orientation.Source: §3.5.2 · pp. 9–10, Figure 4 · p. 10

3.5.3 Translations

• Translational mutations move selected molecular centers of mass; equal-and-opposite translations of paired enantiomers can break glide symmetry and form a tightly packed dimer.The example produces an asymmetric unit containing two molecules.Source: §3.5.3 · p. 10

3.5.4 Permutations

• Permutation mutations swap molecular center-of-mass positions and may preserve, add, or break space-group symmetries depending on molecular, lattice, and permutation symmetry.The example swaps symmetry-related positions while retaining the parent’s P21/c symmetry.Source: §3.5.4 · p. 10

3.5.5 Permutation-Rotations and Permutation-Reflections

• Permutation-rotation and permutation-reflection mutations swap molecules and then reorient them, producing examples that break parent symmetries and yield same-chirality asymmetric units.For crystals with chiral molecules, permutation-reflection can swap the relative positions of enantiomers.Source: §3.5.5 · pp. 10–11, Figure 5 · p. 11

3.6 Rejection Criteria

• Pre-optimization rejection criteria screen offspring to preserve population diversity and avoid expensive optimization of unreasonable or redundant structures.The checks address geometry, similarity, and energy.Source: §3.6 · p. 11

3.6.1 Geometry Checks

• Geometry checks can reject structures with intermolecular contacts below a threshold based on the atoms’ van der Waals radii and a user-defined scale.The scale is typically set between 0.6 and 0.9; users may also constrain center-of-mass distances or unit-cell volume.Source: §3.6.1 · p. 11

3.6.2 Similarity Checks

• Similarity checks compare offspring before and after optimization to avoid redundant calculations and identify duplicates despite different crystal representations.Pre-relaxation checks use loose site and lattice tolerances; post-relaxation checks use stricter tolerances and normalized RMS distances between equivalent atomic sites.Source: §3.6.2 · pp. 11–12

3.6.3 Single Point Energy (SPE) Cutoff

• Single-point PBE+TS calculations reject unrelaxed offspring above a user-defined energy cutoff before local optimization.The relative cutoff is typically 70–120 kJ/mol per molecule, but may depend on the system.Source: §3.6.3 · p. 12

3.7 Termination

• A run may stop when the top structures in the common population remain unchanged for a user-defined number of iterations.For example, termination can follow 50 iterations without change in the top 20 structures; an iteration is when a structure passes rejection criteria and enters the common pool.Source: §3.7 · p. 12
• Users can instead stop when the total population reaches a specified size or terminate the code manually.Source: §3.7 · p. 12
• After interruptions such as wall-time limits or hardware failures, a restart can finish leftover calculations before new selection and can change GA settings without affecting the common population.Settings examples include tightening energy cutoffs or changing mutation schemes.Source: §3.7 · p. 12

4.1 Target XXII

• All Target XXII runs generated the experimental structure; tournament selection with energy-based fitness also found the second-lowest-energy structure, while symmetric crossover found the experimental structure sooner on average.The experimental structure was approximately 3.3 kJ/mol per molecule above the global minimum. Energy-based runs brought the top-20 average to within approximately 5 kJ/mol per molecule of the minimum; the cluster-based run reached within 7 kJ/mol per molecule.Source: Figure 6 · p. 13, §4.1 · pp. 13–14
• For Target XXII, cluster-based fitness sampled the landscape more uniformly, reducing over-selection of planar structures and reaching regions outside the initial pool, but at increased computational cost.Energy-based fitness showed initial-pool bias and genetic drift; lattice-parameter clustering promoted under-represented structures, including bent-conformation and elongated-cell structures.Source: Figure 8 · p. 15, §4.1 · pp. 15–16

4.2 Target II

• All Target II runs generated the experimental structure and at least one lower-energy structure; the 100% mutation run with tournament selection had the lowest top-20 average energy.Strain mutations were particularly effective at generating new low-energy structures. The top 100 structures fell within 5.3, 5.5, and 6.1 kJ/mol per molecule using PBE+TS, PBE+MBD, and PBE0+MBD, respectively.Source: Figure 10 · p. 17, §4.2 · p. 17, Figure 11 · p. 18
• The #1 PBE0+MBD Target II structure has a previously unreported scaffold motif and Z′>1; GAtor’s crossover and mutation operators enabled its generation.The experimental structure ranked #10 with PBE0+MBD, while the scaffold structure was ranked as the global minimum by all three DFT methods.Source: §4.2 · pp. 17, 19, Figure 11 · p. 18

4.3 Target XIII

• For Target XIII, only one of six GA runs found the experimental structure, but higher-level PBE+TS re-ranking placed it first; symmetric crossover was essential for producing low-energy structures.The experimental structure was nearly 6 kJ/mol per molecule above the global minimum with lower-level settings, illustrating how the hierarchical search settings could bias structure selection.Source: §4.3 · pp. 19–21, Figure 13 · p. 20

4.4 Target I

• For Target I, all Z=4 runs found structures below the experimental polymorph at lower-level settings, while some Z=8 symmetric-crossover runs were slower and missed the Z=8 polymorph.The Z=4 runs had similar minimum-energy convergence across tested settings; the Z=8 runs using 25% or 50% symmetric crossover with roulette-wheel selection did not locate the Z=8 polymorph before termination.Source: Figure 14 · p. 21, §4.4 · p. 21
• Four of the six top PBE+MBD structures re-ranked with PBE0+MBD have packing motifs similar to metastable Pbca, with co-facial dimers oriented oppositely.The top 100 structures were re-relaxed with higher-level PBE+TS settings, re-ranked with PBE+MBD, and the top six were re-ranked with PBE0+MBD.Source: §4.4 · p. 22
• The experimental Z=4 P21/c polymorph ranks #11 with PBE+TS and #1 with PBE+MBD and PBE0+MBD; metastable Z=8 Pbca ranks #1, #4, and #3, respectively.The methods find the Pbca and putative Z=8 P21/c structures practically energetically degenerate, but the latter has Z′=2 and a slightly different lattice, making it a unique lattice-energy minimum. The study identifies many-body dispersion and long-range screening in MBD as important to the ranking.Source: §4.4 · pp. 22–23
• The six highlighted structures have very similar densities, but the most stable experimental P21/c structure has the lowest density; relative energy differences increase with PBE0+MBD versus PBE+MBD.There is no significant re-ranking between PBE+MBD and PBE0+MBD for the structures considered.Source: §4.4 · p. 22
• Earlier computational studies after the second blind test ranked Z=8 Pbca as most stable, although recrystallization attempts yielded stable Z=4 P21/c; a cited study proposed Pbca may be a saddle point.The proposed symmetry breaking could produce the Z′=2 P21/c structure, but whether structures are minima or saddle points may depend strongly on the energy method used (Refs. 156–158; Ref. 68; Refs. 159,160).Source: §4.4 · pp. 22–23

Conclusion and Best Practices

• GAtor combines energy- or cluster-based fitness with multiple selection, crossover, and mutation options; its molecular-crystal operators support exploration of high-Z′ structures.The authors recommend diverse initial pools, several shorter runs with different settings, and evolutionary niching in at least one run.Source: §5 · pp. 23–24

Results

• The experimental structures of all four targets were predicted, alongside additional low-energy structures, but the most effective GA settings differed by target.For Target XIII, all three ranking methods placed the experimental form first; for Target XXII, only PBE0+MBD did so.Source: §5 · p. 23
• Target-specific settings mattered: symmetric crossover and tournament selection helped for Target XXII, standard crossover for Target II, and symmetric crossover for low-energy Target XIII structures.For Target I, tested settings were equally effective for low-energy Z=4 structures, while some symmetric-crossover and roulette-wheel combinations were less effective for Z=8.Source: §5 · p. 23
• PBE+MBD and PBE0+MBD ranked Target I’s Z=4 polymorph first and its Z=8 polymorph fourth and third, respectively; MBD gave the correct lattice-energy ordering.Target XXII’s experimental structure ranked first only with PBE0+MBD.Source: §5 · p. 23

Conclusion and Best Practices

• Low-energy structures from different runs were re-relaxed and re-ranked with PBE+TS, PBE+MBD, and PBE0+MBD.The authors identify cluster-based fitness as useful for uniform sampling and for overcoming initial-pool and selection biases.Source: §5 · pp. 23–24

Results

• For Target II, the previously unreported Z′=2, P¯1 scaffold structure ranked first with all three DFT methods, 1.8 kJ/mol per molecule below the experimental form.Several low-lying putative structures have higher computed densities than the observed structure, motivating possible high-pressure experiments.Source: §5 · p. 24

• Supporting information includes experimental-versus-predicted structure comparisons with RMS differences, initial-pool space-group distributions, and CIFs and total energies for structures used in re-ranking.Source: Supporting Information Available · p. 24