CGWGAN: crystal generative framework based on Wyckoff generative adversarial network

INTRODUCTION

Crystals play extremely important roles in various scientific disciplines, including physics, chemistry, biology, etc., especially in materials science and engineering. They embody distinct and precise structures of the atomic arrangement. A perfect crystal is characterized by a repeating basic unit, known as the crystal lattice, which, when translated in parallel, completely fills the whole three-dimensional (3D) space without any gaps. Crystals are categorized into 14 Bravais lattices, belonging to seven crystal systems, according to point group theory^[1] and are further divided into 230 space groups^[2]. In addition to crystal structure, chemical composition is also essential information of crystals. For a given molecule with a specific chemical formula, its atoms can arrange into entirely different crystal structures, a phenomenon known as isomerism. To fully define a crystal, one must know its structure, the types and arrangement of atoms, and the number of atoms, which leads to an extremely large variety of possible crystals. The arrangement of atoms within the crystal lattice is usually expressed by fractional coordinates and satisfies the crystal symmetry operations. In 1922, Wyckoff introduced the concept of Wyckoff positions, which summarize the equivalent positions within space groups. A Wyckoff position in a space group includes all sites whose site-symmetry groups are conjugate subgroups. Each Wyckoff position in a space group is labeled by a letter, known as the Wyckoff letter. The number of distinct Wyckoff positions in a space group is finite, with a maximum of 27 in space group Pmmm. In total, there are 1,731 Wyckoff positions across all space groups^[3]. Based on Wyckoff positions, one can define an irreducible unit in a lattice cell, termed the asymmetry unit (ASU)^[3]. An asymmetric unit of a space group is the smallest, indivisible part from which the entire lattice cell can be constructed by applying symmetry operations to it. Constructing a lattice cell from an ASU necessitates the precise determination of Wyckoff positions. However, in crystal generative deep learning models, atom sites are often generated with noise, making it challenging to achieve exact matches between atom sites across subgroups. The present work studies this issue deeply.

In addition to the commendable efforts in high-throughput computations and experiments^[4], the data-driven exploration of novel crystals is rapidly advancing in tandem with the significant progress in artificial intelligence (AI) and machine learning (ML)^[5]. Various ML models have been proposed to design and discover new crystals, demonstrating remarkable success. For example, the Universal Structure Predictor, Evolutionary Xtallography (USPEX), is an evolutionary algorithm-based method^[6] and enables crystal structure prediction under arbitrary conditions of pressure and temperature, relying solely on the chemical composition of the material. Wang et al. developed a method for crystal structure prediction “from scratch” using the particle swarm optimization (PSO) algorithm^[7]. This method efficiently minimizes free-energy surfaces, combining total-energy calculations with the PSO technique, and requires only chemical compositions to predict stable or metastable structures under specified external conditions, such as pressure. Shao et al. adopted a symmetry tree graph and an AI-based symmetry selection strategy to decompose the search space into Wyckoff position, thereby simplifying the prediction process^[8]. Their method was validated by successful predictions in binary Lennard-Jones mixtures and high-pressure ice phases with over 100 atoms. Ryan et al. employed atomic fingerprint (AFP)^[9], which represents purely topological relations between atoms in the crystal structure and reduces the 3D crystal structures to one-dimensional (1D) representations. They used a variational autoencoder (VAE) to discern the atoms of different elements based on the topology of their crystallographic environment^[9]. The trained model is able to generate crystal templates in order to predict the likelihood of forming new compounds by placing elements into these structural templates in a combinatorial fashion. Zhao et al. developed the physics guided crystal generative model (PGCGM) for efficient crystal material design, incorporating physical constraints related to interatomic distances^[10]. They conducted density functional theory (DFT) calculations to validate the generated structures in terms of thermodynamic stability. Xie et al. proposed the crystal diffusion variational autoencoder (CDVAE) model^[11]. CDVAE incorporates SE(3)-equivariance^[12] to learn from the data distribution of stable crystals. During the diffusion process, it adjusts atomic positions toward lower energy states and updates atom types to satisfy bonding preferences between neighboring atoms. Zhu et al. proposed a framework that combines Wyckoff position encoding with DFT to generate symmetry-compliant, stable inorganic crystals, facilitating the discovery of novel materials and the accurate prediction of ground-state and polymorphic structures^[13].

In AI crystal generation, a lattice cell is commonly characterized by a triplet of matrixes (or vectors), containing two types of fundamental components. The first type of components is topology- and position-related, such as the three lattice constants, fractional coordinates, and space groups. The second type of component concerns the arrangement of various types of atoms, a critical factor influencing both the energy of crystal formation and atomic bonding. The two steps can be integrated into a single model using either a conditional generative framework or sequence-based methods. In a conditional generative model, a well-pretrained model can capture the interaction between symmetry and energy. However, developing such a state-of-the-art (SOTA) model is quite challenging. The sequence strategy helps alleviate the model’s training burden. Our study adopts a similar approach: we generate a crystal site topology that precisely adheres to the space group symmetry, referred to as a crystal template. Subsequently, atoms are positioned based on the potential energy surface.

Feedback from the potential energy surface is crucial during crystal generation. Relying solely on DFT calculations is highly resource-intensive and time-consuming. Therefore, the crystal generation framework also incorporates SOTA surrogate models for energy calculations. Chen et al. developed the universal MatErials graph network (MEGNet) for accurate property prediction of molecules and crystals, achieving an accuracy of 86.2% in crystal formation enthalpy prediction^[14]. The universal graph deep learning interatomic potential (IAP), named M3GNet^[15], uses graph neural networks (GNNs) with three-body interactions. After being trained on the comprehensive DFT crystal dataset from the Materials Project (MP), M3GNet has demonstrated its efficiency in predicting stable materials by identifying 1.8 million potentially stable materials from a screening of 31 million hypothetical crystal structures. Merchant et al. explored scaling deep learning for materials discovery and developed graph networks for materials exploration (GNoME)^[16]. GNoME focused on large-scale active learning and structural prediction, achieving a prediction accuracy with a mean absolute error of 11 meV/atom for crystal formation energy and a stable prediction rate above 80%. This led to the discovery of over 2.2 million stable crystal structures. Furthermore, the application of transformers and attention-based algorithms has shown great promise in the study of crystals. For example, CrysalGPT^[17] demonstrates exceptional sequence-to-sequence transferability while maintaining accurate time-series predictions. Additionally, the molecular set transformer^[18], a computational screening tool, efficiently prioritizes molecular pairs capable of forming stable multicomponent structures, such as co-crystals.

This work develops the crystal generative framework based on Wyckoff generative adversarial network (CGWGAN). Figure 1 illustrates the architecture of CGWGAN, which consists of three modules: crystal template generation, atom filling, and energy/phonon calculations. In the crystal template generation module, a triplet of the asymmetric unit ($$ \mathcal{A} $$), space group ($$ \mathcal{G} $$), and lattice constant vector ($$ \mathcal{L} $$) ensures the structural identity of the crystals. The ASU, with its Wyckoff positions, seamlessly integrates with the group theory-guaranteed symmetry and periodic translation conditions, establishing a novel Wyckoff self-fulfillment crystal generation framework. After training, the generator can produce crystal templates to enhance diversity in crystal generation. These templates are then filled with atoms of different and/or the same elements to create atoms-loaded crystals. The formation energy and phonon stability of these crystals are screened using M3GNet. Finally, the surviving crystals are evaluated through DFT and phonon calculations to confirm their synthesis feasibility.

Figure 1. The architecture of the CGWGAN, including three modules of the crystal template generation, the atoms-filling in of different or the same atomic elements, and the DFT and phonon calculations. In the crystal template generation, a novel triplet of ASU ($$ \mathcal{A} $$), space group ($$ \mathcal{G} $$), and lattice constant vector ($$ \mathcal{L} $$) ensure the structural identity of crystals. CGWGAN: Crystal generative framework based on Wyckoff generative adversarial network; DFT: density functional theory; ASU: asymmetry unit.

MATERIALS AND METHODS

A crystal can be decomposed into two key components: (1) the arrangement of sites, which reflects the lattice geometry constrained by the symmetry group; and (2) the atomic attributes, which determine the energy distribution and bonding between atoms. The generation of a CGWGAN model follows this sequence, beginning with the template (site) generation, followed by the atom infill steps. A crystal can be described by three fundamental elements: the asymmetric unit $$ \boldsymbol{\mathcal{A}} $$; the space group $$ \boldsymbol{\mathcal{G}} $$; and the lattice parameters $$ \boldsymbol{\mathcal{L}} $$. These elements enable the derivation of all symmetry operations - including translations, rotations, reflections, and inversions - from any element g ∈ $$ \boldsymbol{\mathcal{G}} $$. Using a finite set of group operations, a perfect crystal can be constructed in 3D space.

For a crystal sample $$ \boldsymbol{\mathcal{C}}^i $$, the local site encoder_L configuration is embedded as:

Here, $$ \bigcup_{g \in \mathcal{G}} g\left(\boldsymbol{r}_{\boldsymbol{\mathcal{A}}}\right) $$ represents all the sites in the $$ \boldsymbol{\mathcal{A}} $$ and g is a given space group from $$ \boldsymbol{\mathcal{G}} $$.

where $$ \mathcal{H} $$ represents the letter in Wyckoff set for a given space group, and x, y, z represent the corresponding fractional coordinates. Key_j signifies the embedding of ASU sites j = 1, 2, 3, …, k in the $$ \boldsymbol{\mathcal{A}} $$, where k is a hyperparameter related to the number of total atoms in the crystal cell. In fact, once the specific space group g and ASU sites are known, x, y, z can all be obtained through operation expansion and are not necessarily used as input when training the model. In addition to the site configuration, crystals are born with physical constraints on global site encoder_G:

where E_j denotes atoms of different and/or same elements in $$ \boldsymbol{\mathcal{A}} $$, (a, b, c, α, β, γ) represents lattice constant vector $$ \boldsymbol{\mathcal{L}} $$. Therefore, a crystal can be fully represented as a vector $$ \hat{C} $$:

The atoms are masked in crystal templates,

The asymmetric unit $$ \boldsymbol{\mathcal{A}} $$ is the minimal unit to represent a crystal geometry and atom sites. CGWGAN generates a crystal template based on the asymmetric unit $$ \boldsymbol{\mathcal{A}} $$ with Wyckoff coordinates, space groups $$ \boldsymbol{\mathcal{G}} $$ and the lattice vectors $$ \boldsymbol{\mathcal{L}} $$.

Generative adversarial networks (GANs)^[19,20] are founded on the principles of game theory and implemented using neural networks. The architecture of GAN is displaced in Figure 2, which has two components: generator and discriminator. The generator begins by sampling data from a latent space and, through a network, transforms the sampled data back into the original feature space, thereby generating fake data. Both the fake data and the real data are then fed into the discriminator. The discriminator’s task is to distinguish between these two types of data. When the discriminator fails to differentiate fake data from real data, the GAN generates potential new data. CGWGAN uses Wasserstein Generative Adversarial Network^[21] (WGAN) with a gradient penalty (GP) training strategy. CGWGAN loss $$ \mathcal{L}_g $$ focuses solely on the overall coherence of the crystal templates,

Figure 2. (A) Generator architecture configuration; (B) A 1D convolutional self-attention block with residual connections in the discriminator. 1D: One-dimensional.

where c is the real sample with distribution p(c), z is a random noise with the p(z) distribution, and G_θ and f_w represent the generator and discriminator, respectively. A self-attention module is integrated into the discriminator. Figure 2A illustrates the architecture of the generator, which consists of N = 3 fully connected layers with ReLU activation functions, followed by an output layer comprising a fully connected layer with a Tanh activation function. Figure 2B illustrates the configuration of the discriminator. The input crystal geometry undergoes, three times, fully connected layers with LeakyReLu and is subsequently fed into two sub-discriminators. The first sub-discriminator employs three 1D convolutional layers to transform the input vector into three distinct output vectors: K (Key), Q (Query), and V (Value)^[22], thereby facilitating the computation of self-attention mechanisms. The second sub-discriminator is constructed similarly to a Residual Network, which serves to mitigate the issue of vanishing gradients. The information from the two sub-discriminators is aggregated by a fully connected layer, culminating in the final decision-making process.

All atom-filled crystals are screened first with those in the database (MP) to ensure that they are new crystals. Then, the pre-trained M3GNet model is adopted to examine their formation energy and kinetic stability, quickly. The crystal candidates that pass M3GNet are labeled with the formation energy and the lowest phonon frequency and stored in an ASE database. The corresponding automation scripts are provided in the data available section. After that, the Vienna ab initio simulation package (VASP), which employs DFT is employed to calculate the formation energy, with the exchange-correlation functional^[23] of generalized gradient approximation (GGA). Detailed specifications of the cutoff energy, K-spacing, and convergence criteria for first-principles calculations are provided in the Supplementary Materials. The formation energy E_form is defined by

where E_total is the total energy of the crystal, n_i and E_i are, respectively, the number and the energy of the i-th atom in its standard state. The structural reliability assessment is conducted through the calculation of phonon spectra utilizing the density functional perturbation theory (DFPT) method, in conjunction with the phonopy library^[24]:

where ω(q, s) is the phonon frequency at wave vector q and phonon mode s, M is the atomic mass, and Φ (q, s) is the force constant matrix. The absence or weakness of imaginary frequencies in the phonon spectrum indicates that there is no tendency for spontaneous alteration in the atomic arrangement within the crystal structure, signifying a state of structural stability. Structures that pass both formation energy and phonon tests of DFT are novel crystals with great potential to be experimentally synthesized.

RESULTS AND DISCUSSION

This study uses 154,714 crystal structures from the MP database, encompassing 228 space groups, excluding space groups No168 and No207, across the seven crystal systems. Among these, 36,368 structures, exhibiting three and four asymmetric unit sites and possessing symmetry beyond P1, underwent screening. The distributions of 154,714 crystal structures in the 228 space groups and seven crystal systems are presented in Supplementary Figures 1 and 2.

The CGWGAN generates complete structures by first creating site templates and then filling them with suitable elements, making the generator inherently efficient. During the generation process, templates with inter-site distances less than 1 Å are automatically filtered out. A minimum distance of 1 Å is set as the lower limit because the smallest diameter of H atom is 1.06 Å^[25], slightly larger than 1 Å, thereby structures with inter-site distances below this threshold are highly unlikely to form viable crystals after atomic infill. This tolerance significantly reduces computational cost without compromising the quality of the generated structures. Notably, 1 Å is a hyperparameter that can be adjusted based on the research system. This restriction effectively removes abnormal high-density structure prototypes. The atomic number density (ρ_A), expressed as $$ \frac{N_{atoms}}{V_{cell}} $$, serves as a common indicator reflecting the distribution pattern of atoms in real space^[26,27]. For instance, the CGWGAN generator produces a total of 140,000 crystal templates, each with a maximum of four Wyckoff positions, allowing for the accommodation of up to four elemental types. Supplementary Figure 3 categorizes these 140,000 crystal templates across all crystal systems. The distribution of synthetic crystal templates as a function of atomic number density (ρ_A), as shown in Figure 3A, overlaps with that of the training dataset by more than 93%, demonstrating a high consistency between the templates and the inherent distribution of the MP database. It is noteworthy that ρ_A was not explicitly provided to CGWGAN during the encoding and training process; however, CGWGAN successfully captured the underlying distribution pattern.

Figure 3. (A) The distribution of the synthetic database and the MP database; (B) The latent space distribution after the 1D convolutional self-attention block of the discriminator. MP: Materials Project; 1D: one-dimensional.

The alignment of features between synthetic templates and real crystals can be examined by analyzing the feature distributions in the latent space of the discriminator in CGWGAN. Figure 3B presents a t-distributed stochastic neighbor embedding (t-SNE) plot of the latent space features, demonstrating that the manifold features of the synthetic template structures closely follow those of the MP real crystals. All generated crystal templates are available as open-source on the Hugging Face platform.

The convergence percentage of generated crystals is used as a metric to evaluate our crystal generators. We use the popular physical method PyXtal^[26] as a baseline. PyXtal automatically finds suitable combinations of Wyckoff positions using only the chemical composition and symmetry group information, employing a stepwise merging scheme. Using the open PyXtal package, we generated crystals ranging from just two atoms (Ni1Se1) to 32 atoms (P32) and, correspondingly, generated these crystals with CGWGAN. Figure 4A shows the convergence percentages of crystals, generated by CGWGAN and PyXtal, against the number of atoms in a unit cell. The convergence percentage of crystals generated by PyXtal decreases approximately and linearly with the increase in the number of atoms in a unit cell, as shown by the dashed blue line. The convergence percentage of crystals generated by CGWGAN is almost independent of the number of atoms in a unit cell. In ternary systems with more than 20 atoms per unit cell (e.g., Mg₄Si₄O₁₂), the convergence percentage of crystals generated by PyXtal is about 50%, while the crystals generated by CGWGAN maintain a convergence percentage of around 95%. Thus, the performance of crystal formation is almost independent of the number of atoms and the type of elements. Furthermore, CGWGAN achieves convergence percentages of 88% and 74%, respectively, for binary and mono-elemental crystal generation. The results show that only if the number of atoms per unit cell is less than 10 for mono-elemental crystal generation, PyXtal is better than CGWGAN, whereas for most cases, CGWGAN performs better than PyXtal. A detailed illustration of the outcomes for CGWGAN and PyXtal is provided in the Supplementary Materials.

Figure 4. (A) The convergence percentage versus the number of atoms in a unit cell of crystals generated by CGWGAN and Pyxtal, covering unary, binary, and ternary systems; (B) Filtering funnel diagram for the selection of Ba-Ru-O crystal candidates; (C) The distribution of formation energy calculated by the DFT relaxation on the Ba-Ru-O 1,156 crystal candidates. CGWGAN: Crystal generative framework based on Wyckoff generative adversarial network; DFT: density functional theory.

Atoms of different elements are infilled into the crystal templates under the supervision of M3GNet to evaluate the formation energy and phonon spectrum of the candidate crystals. Candidates that pass the M3GNet screening have negative formation energies and no imaginary phonon modes in the entire Brillouin zone. These candidates then undergo final validation with CGWGAN through DFT calculations for formation energy and dynamic stability. DFT relaxation calculations identify energetically favorable atomic positions by allowing lattice constants and all atoms to shift slightly from their original locations. The convergence percentage is defined as the proportion of structures that converge under the calculation conditions of an energy convergence criterion of 1 × 10^-7 eV, a maximum of 300 relaxation steps, and a maximum shift magnitude of 0.1 Å. Before phonon calculations, it is crucial to perform a higher-precision optimization of the candidate structures. This includes refining the k-point mesh in reciprocal space with a spacing of 0.2 Å^-1 and applying Van der Waals force corrections. All electronic minimization algorithms update all bands for simultaneous orbital updates with plane-wave cutoff energy based on the pseudopotential. Notably, due to the absence of specific symmetry information in space group P1, we further calculated six eigenvalues related to mechanical stability based on elastic stability criteria^[28].

Figure 4B summarizes the crystal generation process for the Ba-Ru-O system through a sequential selection approach within the CGWGAN workflow. Initially, 120,000 crystal templates, each with a maximum of four Wyckoff positions, were generated. After inserting Ba-Ru-O atoms into the ASU sites, 1,795,032 crystal candidates were generated, as shown in Figure 4B. Subsequently, the generated crystal candidates were deduplicated against the Ba-Ru-O chemical compositions recorded in the MP dataset using the crystal matching function of Pymatgen^[29]. This process removed two crystals already present in the MP database, leaving 1,795,030 newly discovered candidates. These Ba-Ru-O crystal candidates were then screened by M3GNet and ranked according to their predicted formation energy. Due to the computational complexity of phonon calculations, only the top 5,000 candidates with the lowest formation energies were selected for phonon spectrum calculations. However, only 1,156 of these 5,000 candidates passed the M3GNet phonon spectrum stability criteria. The phonon spectrum stability screened by M3GNet is based on the criterion that the minimum phonon frequency should be greater than -0.5 THz. DFT relaxation was then performed on the 1,156 crystal candidates. Figure 4C shows the distribution of their formation energies, all of which are negative. Phonon frequencies were calculated using DFPT^[30] by establishing a 2 × 2 supercell. The stability criterion in the DFT-based phonon spectrum is greater than 0 THz. Finally, the top 30 candidates with the lowest energies underwent rigorous phonon spectrum calculations, resulting in seven newly discovered crystals.

Table 1 lists the seven generated novel computational crystals. Using Pymatgen symmetry checking, we found that these novel structures correspond to four crystal systems and seven space groups: Triclinic (#1), Orthorhombic (#25, #35, #71), Tetragonal (#99, #123) and Trigonal (#166). Additionally, the structures were verified using AFLOW’s XtalFinder, and the results are provided in the Supplementary Materials.

Figure 5 presents the phonon spectra for these candidate crystals, corresponding to space group numbers: (a) 1, (b) 25, (c) 35, (d) 71, (e) 99, (f) 123, (g) 166. For the crystal structures shown in Figure 5E and F, only minor discrepancies are observed in the phonon spectra. Additional X-ray diffraction (XRD) simulations using WPEM software^[32], as depicted in Figure 5H, confirm that these structures correspond to different configurations with varying asymmetric unit compositions in distinct space groups. This demonstrates that CGWGAN algorithm effectively learns the deeper rules of crystal geometry distribution. Apart from the structures shown in Figure 5, numerous metastable structures were also generated and most of them have negative formation energies but exhibit imaginary frequencies at some high-symmetry points in phonon calculations.

Figure 5. Stable phonon spectra of generated samples in Ba-Ru-O system: this calculation is conducted on a supercell two times along each axis. Visualization is conducted using VESTA^[31].

Some of the newly discovered crystals might belong to the superoxides^[33,34], whose valences might be higher than the typical oxidation states observed in common oxides. The Bader method^[35,36] is widely used to partition charges to each ion. For example, the charges in the newly discovered crystal BaRuO₆ are partitioned by the Bader charge analysis. The valence electrons are 11 e for Ba, 13 e for Ru, and 6 e for O. Figure 6A visualizes the isosurface of charge density in the irreducible lattice cell of BaRuO₆, including one Ba atom, one Ru atom and six O atoms. Figure 6B shows the corresponding electronic band structure, indicating that the O atoms between O₂-O₆ and O₁-O₅ exhibit significant electron cloud overlap, whereas O₃ and O₄ do not significantly exchange electrons with other O atoms and acquire more effective charge. Furthermore, the bands in the electronic band structure near the Fermi level are mainly contributed by O atoms with a minor contribution from Ru atoms, while Ba remains at a deep level; consequently, no charge associated with Ba is visualized on the isosurface of charge density. Table 2 lists fractional charge to each ion inside the crystal from the Bader charge analysis.

Figure 6. (A) The isosurface charge density in the irreducible lattice cell of BaRuO₆; and (B) the corresponding electronic band structure.

First-principles calculations select generated crystals with a negative formation energy, ensuring that the generated crystals are thermodynamically stable at zero K. Phonon spectrum calculations guarantee that the generated crystals are mechanically stable at zero K. These two properties are essential and necessary conditions for experimental synthesis of these generated crystals. In general, experimental synthesis of crystals is analyzed by two aspects. The first is the thermodynamic aspect that the Gibbs free energy of products (synthesized crystals) must be lower than the Gibbs free energy of reactants (raw materials). The second is the kinetic aspect that the resistance of paths from raw materials to products must be low, which is usually represented by energy barriers. Both synthesis thermodynamics and kinetics can be varied by adjusting the synthesis temperature, pressure, and other applied electric, magnetic, or microwave fields, and by utilizing catalysts. All of these are beyond the scope of the present study.

Regarding the dataset distribution used in this study, the generated distributions, statistical details on generated crystals and relaxed lattice parameters, a step-by-step breakdown of CGWGAN’s generative process, comparisons with other models or algorithms, and information on computational efficiencies, please refer to the Supplementary Materials.

CONCLUSIONS

In summary, CGWGAN has demonstrated its capability to efficiently generate high-quality crystal structures of target crystal systems and space groups. As an example, the present work discovers seven novel crystals within the Ba-Ru-O system, boasting high thermodynamic and kinetic stability. Obtaining the ASU through known operations in group theory, this physical encoding enables a fully reversible generative process. Beyond the architecture applied in the paper based on generative adversarial neural networks, this framework is versatile and applicable to any generative model framework, including VAEs and diffusion models. Apart from crystal generation, CGWGAN can be easily integrated with AI XRD analyzer^[37] service as out-library phase identification methods, thereby considerably fostering the development of materials informatics. The current CGWGAN has not paid sufficient attention to material properties, while the dream of material scientists and engineers is to design and discover novel crystals based on the need of material properties in scientific and engineering practice. To materialize this dream, CGWGAN will be improved by adding the conditions of certain designed material properties. In addition to M3GNet which is trained on the MP database, experimental crystal databases will be used to enhance the synthesizability of generated crystals from CGWGAN.