PaperlayerView exact revision on arXiv

Source-linked AI summary

MolCrystalFlow: Molecular Crystal Structure Prediction via Flow Matching

Cheng Zeng, Harry W. Sullivan, Thomas Egg, Maya M. Martirossyan, Philipp Höllmer, Jirui Jin, Richard G. Hennig, Adrian Roitberg, Stefano Martiniani, Ellad B. Tadmor, Mingjie Liu

arXiv:2602.16020v3cs.LGcond-mat.mtrl-sci

TL;DR

Molecular crystal structure prediction remains difficult because fully periodic crystals combine complex intra- and intermolecular interactions with large configuration spaces. MolCrystalFlow models molecules as rigid bodies and learns lattice, orientation, and centroid variables on geometric manifolds. It outperforms MOFFlow, performs competitively with Genarris, and supports an integrated pipeline that identifies experimentally similar polymorphs for two target compounds.

  • Problem

    Fully periodic molecular crystal prediction remains elusive because large molecules, complex interactions, and periodic constraints make structure generation difficult.

    Source: Introduction · pp. 1–2

  • Method

    MolCrystalFlow uses a periodic E(3)-invariant flow model that represents molecular conformers as rigid bodies and generates lattice matrices, centroids, and orientations on their native manifolds.

    Source: Introduction · p. 2, Overview of MolCrystalFlow · p. 3, Neural Network Architecture · p. 6

  • Results

    MolCrystalFlow outperforms MOFFlow, achieves competitive performance against Genarris, and identifies polymorphs for two target compounds with energies and packing geometries close to experiment.

    Source: Introduction · pp. 2–3, Discussion · p. 10

  • Takeaways & Limitations

    The model and integrated generative-evaluation pipeline provide a basis for more efficient exploration and discovery of molecular crystal polymorphs.

    Source: Introduction · p. 3, Discussion · p. 10

    The data-driven model uses structural information without explicit energetic information, making intrinsically low-energy packing generation challenging.

    Source: Discussion · p. 10

Abstract

from arXiv · show

Molecular crystal structure prediction represents a grand challenge in computational chemistry due to large sizes of constituent molecules and complex intra- and intermolecular interactions. While generative modeling has revolutionized structure discovery for molecules, inorganic solids, and metal-organic frameworks, extending such approaches to fully periodic molecular crystals is still elusive. Here, we present MolCrystalFlow, a flow-based generative model for molecular crystal structure prediction. The framework disentangles intramolecular complexity from intermolecular packing by embedding molecules as rigid bodies and jointly learning the lattice matrix, molecular orientations, and centroid positions. Centroids and orientations are represented on their native Riemannian manifolds, allowing geodesic flow construction and graph neural network operations that respects geometric symmetries. We benchmark our model against a state-of-the-art generative model (MOFFlow) for large-size periodic crystals and a rule-based structure generation method (Genarris) on two open-source molecular crystal datasets. MolCrystalFlow outperforms MOFFlow while achieving competitive performance against Genarris. We also demonstrate an integration of MolCrystalFlow model with universal machine learning potential to accelerate molecular crystal structure prediction, paving the way for data-driven generative discovery of molecular crystals.

Introduction

Molecular crystal structure prediction is difficult because polymorphism, competing interactions, and periodic constraints create vast search spaces with major pharmaceutical and materials consequences. MolCrystalFlow addresses this gap with a periodic generative framework that models molecular crystals as rigid-body packings and integrates structure generation with energy evaluation.

Motivation

  • Polymorphism can substantially change drug solubility, bioavailability, stability, and manufacturability, as illustrated by the emergence of Ritonavir Form B after commercialization.The unexpected form caused precipitation, reduced bioavailability, product withdrawals, and reformulation costs.Source: Introduction · pp. 1–2
  • Conventional CSP workflows generate and rank enormous candidate pools, often requiring millions of CPU hours while remaining system-specific and computationally intensive.Candidate structures must also be deduplicated, clustered, filtered, and evaluated before stability rankings are computed.Source: Introduction · p. 2

Generative modeling

  • Generative models offer a route beyond brute-force search, but larger periodic crystals remain difficult because they require long-range interactions, periodic symmetry, and scalable configuration-space exploration.An all-atom inorganic model matched about 70% of structures below 20 atoms but only 27.4% near 50 atoms.Source: Introduction · p. 2

Related work

  • Existing rigid-body approaches for frameworks and clusters do not directly enforce the periodic translational structure required for molecular crystals, leaving stable and diverse periodic polymorph generation unresolved.Related methods either lack natural periodic lattices, do not enforce periodic translational invariance, or require post hoc lattice inference.Source: Introduction · p. 2

Contribution

  • MolCrystalFlow is a periodic E(3)-invariant flow-based model that conditions on molecular conformers, treats them as rigid bodies, and generates the crystal packing with explicit lattice periodicity.The model is benchmarked on two public datasets against MOFFlow and is integrated with universal machine-learning potentials and DFT calculations.Source: Introduction · p. 2
  • The integrated pipeline identifies polymorphs for two of three open CSP competition targets with energies and packing geometries close to experimental structures.This demonstrates an end-to-end workflow combining generative sampling with energetic and structural evaluation.Source: Introduction · p. 3

Results

MolCrystalFlow represents molecular crystals as rigid-body building blocks and jointly generates lattice, centroid, and orientation variables on symmetry-aware manifolds. Across open datasets and CSP targets, it improves structure-generation metrics over MOFFlow, remains competitive with Genarris-3, and supports a u-MLIP/DFT ranking workflow.

Overview of MolCrystalFlow

  • MolCrystalFlow reduces crystal-packing search by treating molecular conformers as rigid bodies and generating lattice matrices, centroids, and orientations for physically plausible polymorphs.The rigid-body approximation fully specifies packing through the lattice matrix, molecular centroids, and rotational orientations.Source: Overview of MolCrystalFlow · p. 3
  • An EGNN creates invariant building-block embeddings augmented with 18 molecular descriptors before joint flow matching generates crystal-level structure variables.The second stage evolves lattice parameters, fractional centroid coordinates, and orientations jointly, reconstructing all-atom structures from fixed internal coordinates.Source: Overview of MolCrystalFlow · p. 3

Neural Network Architecture

  • MolCrystalNet enforces periodic E(3) symmetry by combining invariant building-block features with time and geometric features while predicting coordinate velocities and denoised lattice and orientation targets.Its message passing includes periodic fractional-coordinate differences, relative orientations, lattice interactions, and permutation invariance over molecules in the unit cell.Source: Neural Network Architecture · p. 6
  • Axis-flip states χ are modeled as intrinsic building-block degrees of freedom, with χ embeddings and grouped optimal transport reducing cross-links between flow paths.The design accounts for molecules sharing or opposing local-frame axis-flip states that cannot be changed through rotational evolution.Source: Neural Network Architecture · p. 6

Accurate and fast generation of densely packed molecular crystals

  • Direct generation substantially outperforms MOFFlow and Genarris-3 across site tolerances, while optimized MolCrystalFlow and Genarris-3 outperform optimized MOFFlow.Genarris-3 is stronger at stricter tolerances after rigid-press optimization, whereas MolCrystalFlow is better at looser tolerances with rigid-press and UMA-OMC optimization.Source: Accurate and fast generation of… · pp. 7–8, Figure 3 · p. 7
  • 3.86±0.07% RMAD measures MolCrystalFlow lattice-volume deviation across 200–1600 Å3, versus 18.8±0.6% for MOFFlow and 59.0±0.35% for directly generated Genarris-3.Rigid-press optimization reduces Genarris-3 deviation to 10.7±0.20%. MolCrystalFlow generates structures in 22 ms each, faster than Genarris-3 at 43 ms but slower than MOFFlow at 6 ms under the stated settings.Source: Accurate and fast generation of… · p. 8

Integrating MolCrystalFlow with u-MLIP for Crystal Structure Prediction

  • On three CSP competition targets, the MolCrystalFlow–u-MLIP–DFT workflow places generated structures in low-energy basins, with DFT identifying near-experimental-energy Z = 4 polymorphs for targets VIII and XI.For target X, DFT places a Z = 2 polymorph suggested by UMA-OMC approximately 3 kJ/mol above the experimental Z = 4 structure.Source: Integrating MolCrystalFlow with u-MLIP… · p. 8, Figure 4 · p. 9
  • Packing similarity is close for targets VIII and XI but not target X, with RMSD7 = 0.397 Å, RMSD10 = 1.107 Å, and RMSD5 = 2.703 Å, respectively.The target-X result highlights difficulty recovering the exact experimental polymorph, while the overall comparison with Genarris-3 is comparable across the three targets.Source: Integrating MolCrystalFlow with u-MLIP… · p. 10

Discussion

MolCrystalFlow combines hierarchical representations with Riemannian flow matching to generate densely packed molecular crystals, outperforming flow-based models and remaining competitive with rule-based methods. Its current scope is limited by the absence of energetic information and by rigid-body assumptions that exclude conformational polymorphism.

Discussion

  • MolCrystalFlow achieves state-of-the-art performance among flow-based models for large periodic systems and competitive performance against multi-shot, rule-based crystal structure prediction methods.The approach uses hierarchical representations and Riemannian flow matching to efficiently sample densely packed molecular crystals.Source: Discussion · p. 10

Limitations

  • The model is trained only on structural information without explicit energetic information, making intrinsically low-energy crystal packing generation challenging.The paper identifies energy-based formulations and inference-time scaling or reweighting as possible remedies.Source: Discussion · p. 10
  • The rigid-body assumption limits the model on crystals where conformational polymorphism is important, motivating extensions that include torsional degrees of freedom.The study’s open datasets satisfy the rigid-body assumption, but many experimentally observed crystals do not.Source: Discussion · p. 10

Conclusion

  • MolCrystalFlow provides an important step toward end-to-end generative design of molecular crystals and lays groundwork for future generative prediction advances.Source: Discussion · p. 10

Methods

MolCrystalFlow represents molecular crystals hierarchically, separating invariant molecular information from global packing variables. Rigid molecules are embedded with an equivariant graph network and represented by fractional centroids, orientations, and a lattice matrix for flow-based generation.

Hierarchical Representation of Molecular Crystals

  • An EGNN embeds each rigid molecular building block through distance-based message passing, producing an E(3)-invariant representation.The distance-only message function ensures invariance under translations, rotations, and reflections.Source: Hierarchical Representation of… · pp. 10–11
  • The building-block embedding is a learnable weighted mean of final atomic features, with weights predicted from node features and atom-to-centroid distances.Source: Hierarchical Representation of… · p. 11
  • Eighteen auxiliary chemical, basic, and geometric descriptors are concatenated with the EGNN embedding as input to the subsequent packing-flow stage.The descriptors include molecular size, chemical properties, and principal-component geometry.Source: Hierarchical Representation of… · p. 11

Molecules in the crystal

  • Each molecule is treated as a rigid body with a fractional centroid and rotational orientation, while fixed PCA-frame atomic coordinates are used to reconstruct the final all-atom crystal.Fractional centroids enforce lattice periodicity, and PCA orientations provide the molecular frame.Source: Hierarchical Representation of… · p. 11

Hierarchical Representation of Molecular Crystals

  • The crystal is parameterized by fractional building-block centroids, rotation matrices, and a lattice matrix, separating local molecular geometry from global crystal packing.F contains centroid positions, R stores orientations, and L defines the lattice.Source: Hierarchical Representation of… · p. 11

Extracting the axis-flip state χ

The axis-flip state χ resolves ambiguities in PCA-based molecular orientations. It is constructed from an equivariant displacement vector projected onto PCA axes, then used to normalize local coordinates and rotations consistently across equivalent building blocks.

Extracting the axis-flip state χ

  • A molecular building block’s PCA axes are obtained by eigendecomposing the covariance matrix of centered atomic coordinates and ordering eigenvectors by descending eigenvalues.Source: Extracting the axis-flip state χ · p. 11
  • The equivariant displacement vector is defined from the center of mass and geometric centroid, with a nearest-neighbor fallback for symmetric molecules.Source: Extracting the axis-flip state χ · p. 12
  • Projecting the displacement onto PCA axes yields signed components whose signs form the binary axis-flip vector η.The method also tracks the number of flipped axes to address symmetry-induced degeneracies.Source: Extracting the axis-flip state χ · p. 12
  • The unique axis-flip state χ is used to renormalize local coordinates and rotations, removing spurious inversions while preserving SE(3) equivariance.This makes local frames consistent across symmetry-equivalent building blocks.Source: Extracting the axis-flip state χ · pp. 11–12

Geometric Symmetry of Molecular Crystal Structures

MolCrystalFlow imposes permutation, rotational, and periodic translational invariance so its learned molecular-crystal distribution remains consistent with crystal geometry.

Geometric Symmetry of Molecular Crystal Structures

  • Permutation invariance ensures that reordering identical building blocks does not change the crystal distribution, while graph neural networks handle this naturally.Source: Geometric Symmetry of Molecular… · p. 12
  • SO(3) invariance requires global crystal rotations to leave the distribution unchanged, with lattice and molecular orientations transformed consistently.Fractional coordinates remain unchanged because they are defined relative to the lattice basis.Source: Geometric Symmetry of Molecular… · p. 12
  • Periodic translational invariance identifies configurations that differ by lattice-vector translations of all building-block centroids.Coordinates are wrapped back into the unit cell after translation.Source: Geometric Symmetry of Molecular… · p. 12

Joint Flow Matching on Riemannian manifolds

MolCrystalFlow learns a joint distribution over lattice matrices, fractional coordinates, and molecular orientations using flow matching on their corresponding Riemannian manifolds. Geodesic interpolation and manifold-aware objectives preserve valid structural states during generation.

Joint manifold flow matching

  • MolCrystalFlow learns the joint distribution p(L, F, R) with a conditional velocity field defined directly on the relevant Riemannian manifolds.Numerical integration therefore produces states that remain on the manifolds.Source: Joint Flow Matching on Riemannian… · p. 12

Structural modalities

  • Centroids use wrapped fractional coordinates on a three-dimensional torus, orientations use axis-angle coordinates on SO(3), and lattices use normalized lower-triangular Euclidean matrices.Lattice lengths are ordered as a ≤ b ≤ c to remove redundant representations.Source: Joint Flow Matching on Riemannian… · p. 12

Geodesic interpolation

  • Flow matching jointly interpolates lattice matrices, fractional coordinates, and orientations along geodesic paths between t = 0 and t = 1.Fractional-coordinate interpolation unwraps the endpoint to choose the shortest torus path, then wraps the result back onto the torus.Source: Joint Flow Matching on Riemannian… · pp. 12–13

Training objective

  • The loss combines denoising endpoint predictions for lattices and rotations with a velocity-field objective for fractional coordinates.The reported weights are λL = 0.1, λR = 1, and λF = 2.Source: Joint Flow Matching on Riemannian… · p. 13

Neural Network Architecture for Periodic E(3)-Invariant Flow Matching

The architecture uses a periodic E(3)-invariant graph neural network to parameterize velocity fields for lattice, coordinate, and orientation updates. Periodic embeddings, relative geometric features, and numerical integration together enforce the required symmetries during generation.

Network architecture

  • A periodic E(3)-invariant graph neural network predicts velocity fields for fractional coordinates, molecular orientations, and lattice matrices.Node embeddings combine building-block representations with a time embedding, while the lattice is shared across the crystal.Source: Neural Network Architecture for… · p. 13

Geometric features

  • Message passing uses relative fractional coordinates, relative orientations, and the current lattice matrix to encode periodic crystal geometry.Fourier-based coordinate embeddings enforce periodicity on the torus, while rotational features use axis-angle quantities and trigonometric embeddings for periodic angles.Source: Neural Network Architecture for… · p. 13

Symmetry constraints

  • Normalized lattice-direction and displacement vectors make message-passing inputs invariant to global rotations and translations.Additional compositional operations preserve lattice-update invariance and rotational equivariance.Source: Neural Network Architecture for… · p. 13

Inference integration

  • During inference, fractional-coordinate velocities are integrated directly, whereas lattice and orientation velocities are computed from denoised and interpolant states.All ordinary differential equations are solved with Euler integration using 50 discrete time steps.Source: Neural Network Architecture for… · pp. 13–14

Base Distributions and Optimal Transport

MolCrystalFlow combines structured base distributions with data-dependent optimal transport to initialize valid crystal modalities and reduce trajectory crossings. The base priors respect periodic, lattice, rotational, and symmetry constraints.

Base distributions

  • The base distributions assign uniform fractional coordinates on [0,1), structured lattice parameters, and a symmetric-rotation prior for molecular orientations.Lattice lengths use Gaussian distributions, angles use a uniform [60°, 120°] distribution, and orientations are generated from a reference SO(3) rotation through symmetry relations.Source: Base Distributions and Optimal Transport · p. 14

Base-distribution symmetry

  • The symmetric-rotation prior preserves full SE(3) invariance when propagated through the equivariant flow-matching process.This keeps the generated crystal distribution p(L, F, R) invariant under the stated rigid transformations.Source: Base Distributions and Optimal Transport · p. 14

Optimal transport

  • Data-dependent optimal transport aligns source and target samples to mitigate trajectory crossings and improve sampling efficiency.Building-block indices are permuted only within the same χ group.Source: Base Distributions and Optimal Transport · p. 14

Geodesic alignment

  • For rotations, the transport assignment minimizes geodesic distances on SO(3), while fractional coordinates use wrapped distances.The optimization acts on paired samples through a permutation operator restricted to each χ group.Source: Base Distributions and Optimal Transport · p. 14

Open-Source Molecular Crystal Dataset

The study preprocesses molecular crystals into rigid molecular building blocks and evaluates MolCrystalFlow on two open-source datasets with filtered, split subsets. The datasets differ substantially in source, elemental coverage, scale, and selection criteria.

Data preprocessing

  • Molecular crystals are decomposed into connected building blocks using element-dependent cutoffs, minimum-image periodic distances, and Cartesian unwrapping.Lattice vectors are sorted by length and transformed into normalized lower-triangular form for learning.Source: Open-Source Molecular Crystal Dataset · p. 14

Thurlemann dataset

  • The Thurlemann dataset contains 11,489 unstrained candidate structures originally, with 11,488 retained after filtering and split into 10,000 training, 738 validation, and 750 test structures.Identical chemical formulas are excluded from different splits.Source: Open-Source Molecular Crystal Dataset · p. 14

OMC25-MCF dataset

  • OMC25-MCF retains 46,801 single-component molecular crystals from OMC25 by selecting one lowest-lattice-energy packing per conformer under a rigid-body approximation.The selected subset spans 12 elements and applies a maximum atomic-force threshold of 0.05 eV/Å.Source: Open-Source Molecular Crystal Dataset · p. 14
  • OMC25-MCF is divided into 42,120 training, 3,681 validation, and 1,000 test structures.The split corresponds to 90% training data, with the remaining structures assigned to validation and testing.Source: Open-Source Molecular Crystal Dataset · pp. 14–15

Evaluation Metrics

Evaluation combines structure matching on raw generated crystals with targeted CCDC CSP benchmarks, including sampling across unit-cell multiplicities and molecular conformational states.

Evaluation Metrics

  • Raw predictions are counted correct when at least one generated structure matches the ground truth under StructureMatcher tolerances.Volume scaling is disabled to preserve intramolecular bond lengths, and multiple candidates are generated because molecular conformers can have numerous polymorphs.Source: Evaluation Metrics · p. 15
  • CCDC CSP evaluation covers targets VIII, XI, and X, whose experimental unit cells contain Z = 8, 8, and 4 molecules, respectively.Sampling is restricted to Z ∈{2, 4} because Z = 1 structures were less stable in preliminary tests and the dataset had limited Z = 8 examples.Source: Evaluation Metrics · p. 15
  • Each CCDC conformer receives 800 sampled structures across two Z values and two χ-state scenarios, with sampling repeated using two random seeds.Structures are optimized and ranked with UMA-OMC before DFT calculations, while 15-molecule assembly fidelity is assessed using COMPACK RMSDs.Source: Evaluation Metrics · p. 15

Computational settings

The study uses matched neural-network implementations, rigid-body optimization, DFT validation, and comparable Genarris and press-optimization workflows to evaluate generated structures.

Computational settings

  • MolCrystalFlow and MOFFlow are implemented in PyTorch and trained and inferred on one Nvidia B200 GPU with 192 GB memory.MOFFlow uses the original work's hyperparameter settings, while MolCrystalFlow hyperparameters are reported in Table S1.Source: Computational settings · p. 15
  • UMA-OMC optimization first relaxes generated structures with rigid bodies and fixed lattices, then relaxes all atomic and lattice degrees of freedom.The stages allow up to 100 and 1000 steps, respectively, with force thresholds of 0.05 and 0.01 eV/Å; the 10 lowest-energy structures per round proceed to DFT.Source: Computational settings · p. 15
  • DFT validation uses VASP with PBE, PAW, a 520 eV cutoff, specified force and energy convergence thresholds, and D3(BJ) or MBD dispersion corrections.Relaxations use conjugate gradients for up to 1500 ionic steps with automatically generated Monkhorst-Pack meshes.Source: Computational settings · p. 15
  • Genarris-3.0 generates trial structures across compatible space groups and compacts them with Symmetric Rigid Press while preserving space-group symmetry.The workflow generates 800 accepted structures per space group for the CCDC targets and 200 for the Thurlemann validation set.Source: Computational settings · pp. 15–16
  • MOFFlow- and MolCrystalFlow-generated structures undergo rigid-body press optimization after 1.5× lattice scaling, with the space group fixed to 1 and volume restored afterward.Nonconverged structures are retained, and approximately 86% of structures converge for both models.Source: Computational settings · p. 16

Related papers

Loading 2602.16020v3…