A Distributional Framework for Generative Modeling of Molecular Crystals

Source: arXiv:2607.05266 · Published 2026-07-06 · By Michael Kilgour, Alex Dong, Mark E. Tuckerman, Jutta Rogal

TL;DR

This paper addresses the critical challenge of modeling molecular crystal polymorphism by proposing a principled framework, MXtalGFlow, to generate and analyze thermodynamically meaningful distributions of molecular crystal structures. Existing crystal structure prediction (CSP) approaches often yield disconnected sets of low-energy candidate structures without a clear distributional interpretation or guarantee of mode coverage. MXtalGFlow leverages a canonical, physics-grounded parameterization of molecular crystals together with energy-based generative flow networks (GFlowNets) trained to approximate the Boltzmann distribution over crystal parameters at a target temperature. By combining these elements, MXtalGFlow produces a coherent thermodynamic distribution over crystal packings, enabling rigorous distributional analysis beyond simple sampling of candidates. The authors demonstrate MXtalGFlow on two molecules and two intermolecular potentials (Lennard-Jones and the Universal Model for Atoms (UMA)), recovering known experimental polymorphs and discovering new packing motifs with competitive probabilities. The work pioneers treating molecular CSP as modeling physical distributions rather than isolated structures, enabling richer thermodynamic insights and quantification of sample diversity, mode coverage, and basin geometry.

Key findings

• A canonical 12-dimensional crystal parameterization was developed that is applicable to all molecules and space groups with Z' = 1 rigid molecules, removing unit cell degeneracy.
• The energy-based GFlowNet diffusion sampler trained on the reward function R(x) = exp(−βE(x)) with β = (kBT)−1 at room temperature successfully approximates the Boltzmann distribution of crystal structures.
• Compared to Lennard-Jones, UMA potentials produce substantially different, mostly non-overlapping crystal distributions for MIPCAS and NEHZOR molecules, influenced by hydrogen bonding.
• Distributional analysis using a radial distribution fingerprint Earth Mover's Distance (RDF EMD) metric identifies distinct structural basins and local maxima in the learned crystal landscape tied to known polymorphs (COMPACK matches RMSD < 0.3Å).
• For MIPCAS under UMA, one basin (basin 1) was predicted as slightly thermodynamically preferred over the known experimental polymorph (basin 2) based on integrated basin probability density.
• For NEHZOR under UMA, two basins (1 and 2) have higher configurational density and probability density than basins containing known polymorphs, suggesting additional plausible packing modes.
• The three-stage training protocol (expand support, thermalize on prior, global equilibration) stabilizes GFlowNet training on this highly multimodal, high-dimensional space without mode collapse.
• Sampling 10,000 crystals from trained models yields diverse, physically meaningful distributions that capture polymorphism and enable explicit computation of thermodynamic observables like basin probabilities.

Threat model

The adversary is hypothetical and not directly applicable here; the work assumes a fixed, known energy function defines a Boltzmann distribution over crystal structures at a given temperature. The model aims to learn and sample from this distribution, with no adversarial manipulation or attacks considered.

Methodology — deep read

The threat model assumes an energy function E(x) defines a Boltzmann distribution over crystal packings at fixed temperature, modeling the thermodynamics; the generative model seeks to approximate this distribution.

The data are generated during training as a support dataset of molecular crystal structures parameterized via a canonical 12D representation (including unit cell vectors, angles, and molecular coordinates of the asymmetric unit). The study focuses on rigid molecules with one asymmetric molecule (Z'=1), in fixed space groups P¯1 for MIPCAS and P21/c for NEHZOR.

The core generative architecture is a generative flow network (GFlowNet) diffusion sampler, implemented as a neural stochastic differential equation (SDE) that models forward and backward policies over trajectories through latent space. The loss is the trajectory balance (TB) loss ensuring consistent forward-backward probabilities relative to the Boltzmann reward R(x) = exp(−βE(x)) where β = 1/kBT (set to ~2.5 kJ/mol representing room temperature). The energy functions tested include classical Lennard-Jones (LJ) and the higher fidelity Universal Model for Atoms (UMA).

Training proceeds in three stages to ensure convergence and mode coverage: (1) expand model support via dataset pre-construction to define the sample space, (2) thermalize model on prior to encourage exploration, and (3) global equilibration to balance forward and backward distributions and approximate the Boltzmann distribution. Careful regularization addresses mode collapse, variance explosion, and overfitting.

Evaluation uses quantitative thermodynamic and distributional metrics: sample energy and density distributions, UMAP embeddings on RDF Earth Mover's Distance (EMD) over atomic radial distributions to visualize structural clusters, and a custom local probability density estimator P(x) based on RDF EMD neighborhoods. Basins are defined via discrete graph hill climbing to local probability maxima, yielding sets of polymorph basins. COMPACK RMSD analysis validates physical similarity of sampled structures to known experimental polymorphs (≤0.3 Å).

10,000 samples per model/molecule were generated, filtering out high-energy (>15 kBT) outliers. Analysis compared the resulting Boltzmann-like distributions under different potentials. Convergence quality was analyzed in supplement but remains imperfect, with known overweighting of high-energy tail samples.

No end-to-end code or weights are currently released, but all parameterization details and training protocols are thoroughly documented. The datasets are synthetically generated on the fly from the energy functions. The approach is thus reproducible in principle but computationally intensive due to the high dimensionality and complexity of the energy landscape.

A concrete example: for MIPCAS in P¯1, the GFlowNet was trained on UMA energies to sample crystal parameters, producing a 12D distribution over asymmetric unit and unit cell parameters. From 10k samples, using RDF EMD distances, probability density maxima identified a major basin containing the known experimental polymorph with match RMSD <0.3 Å. Other basins with slightly lower energy minima but higher configurational entropy were also found. This demonstrates the model captures both energetic and entropic features of the thermodynamic distribution.

Technical innovations

• A canonical, physics-grounded 12-dimensional parameterization of molecular crystals accounting for unit cell uniqueness and crystallographic symmetry for use in generative modeling.
• The application and adaptation of energy-based generative flow networks (GFlowNets) with a novel three-stage training protocol to sample Boltzmann distributions over complex, high-dimensional, multimodal crystal structure spaces.
• The introduction of a distributional analysis framework based on a radial distribution function Earth Mover's Distance metric to identify local probability maxima and structural basins in generated crystal distributions.
• Demonstration that the choice of energy function (Lennard-Jones vs UMA) drastically alters the sampled crystal distribution landscapes, underscoring the framework's utility to explore physical effects on polymorphism.

Datasets

• Synthetic molecular crystal samples generated on-the-fly via MXtalGFlow training for two molecules: MIPCAS (~rigid, P¯1), NEHZOR (~rigid, P21/c).

Baselines vs proposed

• Lennard-Jones potential: energy and density distributions show broader coverage but less physical fidelity vs UMA: models yield sharper, more chemically realistic polymorph basins consistent with experiments.
• Classical CSP tools: Mode discovery more efficient but yield disconnected structure sets vs MXtalGFlow: yields full approximate Boltzmann distributions with explicit basin probabilities.
• No direct quantitative metric comparison reported against prior generative molecular CSP since this is a novel distributional modeling approach.

Figures from the paper

Figures are reproduced from the source paper for academic discussion. Original copyright: the paper authors. See arXiv:2607.05266.

Fig 1: Illustration of our crystal parameterization scheme and crystal building workflow. In energy-

Fig 2: Diagrams of MIPCAS and NEHZOR conformers in (a) and (b).

Fig 3: Summary energy and density statistics for distributions of crystals sampled from trained

Fig 4: UMAP embeddings of RDF EMD matrices for 5k samples each from LJ and UMA GFlowNet

Fig 5: Landscape analysis for the learned P¯1 MIPCAS distribution under UMA. (a) A UMAP

Fig 6 (page 6).

Fig 7 (page 6).

Fig 8 (page 6).

Limitations

• Models trained only on rigid molecules with a fixed Z' = 1 and single space groups, limiting generalizability to flexible molecules or mixed symmetry.
• Sampling convergence is imperfect: the GFlowNet models overweight higher-energy tail structures, reflecting incomplete equilibration.
• Energy functions (LJ and UMA) are approximate: true ab initio and quantum effects are not captured, so thermodynamic predictions depend on chosen potential validity.
• Distributional analysis relies on heuristic RDF EMD metrics and cutoff parameters, which may affect basin assignment and probability estimates.
• The latent transform from crystal parameters to model input is not volume preserving, requiring Jacobian corrections that may introduce modeling complexity.
• No released code or pre-trained model weights yet, hindering immediate reproduction.

Open questions / follow-ons

• How to extend MXtalGFlow to flexible molecules, multiple asymmetric units, and varied space groups in a unified framework?
• Can more accurate, learned ab initio-based energy potentials be integrated to improve quantitative thermodynamic predictions?
• How to further improve convergence and mode coverage in extremely rugged, high-dimensional crystal energy landscapes?
• What are the best metrics and methods for rigorous validation of learned thermodynamic distributions against experimental crystallization outcomes?

Why it matters for bot defense

Although this work focuses on molecular crystal generation rather than bot detection or CAPTCHA defense, its core idea of sampling from a complex, multimodal distribution using an energy-based generative flow network is broadly relevant. For bot-defense, analogous challenges exist in modeling distributions of user behavior or request patterns, where ensuring full mode coverage and sampling realistic diverse behaviors is crucial. The demonstrated training strategies to stabilize GFlowNets and avoid mode collapse on rugged distributions could inspire robust modeling of legitimate vs. malicious traffic distributions. Moreover, the formal grounding in Boltzmann distributions highlights how introducing physical or domain-specific energies as training rewards enables interpretable probability models beyond naive generation. Distributional analysis methods like local density basins and graph hill climbing from this work could be adapted to detect clusters of suspicious bot activity or subtle polymorphic attack patterns. In sum, the paper provides an instructive example of principled generative modeling over complex combinatorial spaces, with lessons applicable to CAPTCHA and bot-defense researchers aiming to model rich, multimodal adversarial distributions.

Cite

@article{arxiv2607_05266,
  title={ A Distributional Framework for Generative Modeling of Molecular Crystals },
  author={ Michael Kilgour and Alex Dong and Mark E. Tuckerman and Jutta Rogal },
  journal={arXiv preprint arXiv:2607.05266},
  year={ 2026 },
  url={https://arxiv.org/abs/2607.05266}
}

Read the full paper

• arXiv:2607.05266 (abstract)
• arXiv:2607.05266 (PDF)