E(3)-equivariant models cannot learn chirality: Field-based molecular generation

## Abstract

Abstract This work introduces FMG, a field-based model for drug-like molecule generation. We show how the flexibility of this method provides crucial advantages over the prevalent, point-cloud based methods, and achieves competitive molecular stability generation. We tackle optical isomerism (enantiomers), a previously omitted molecular property that is crucial for drug safety and effectiveness, and thus account for all molecular geometry aspects. We demonstrate how previous methods are invariant to a group of transformations that includes enantiomer pairs, leading them invariant to the molecular R and S configurations, while our field-based generative model captures this property.

## 1 Introduction

Developing generative models for a specific data modality requires carefully choosing the generative method, feature representation, and neural network (NN) architecture. Field-based representations offer a unifying perspective since modeling continuous domain functions is needed in many domains. In this work, we show how field-based molecule representations can be modeled successfully, achieving competitive performance, compared to the point-cloud molecule representation counterpart.

Organic molecules generally contain a restricted number of atom types, bound mostly through covalent bonds. Most atoms in organic molecules are non-ionized, and distances between two bonded atoms should follow empirically attested values. In an attempt to achieve realistic bond distances, rotation equivariant NN architectures have been devised (Peng et al., ; Guan et al., ). The developed methods only receive pairwise atom distances, which are rotation invariant.

While the aforementioned properties clearly incentivize the development of rotational equivariant architectures, this may come at the cost of potentially incorrect enantiomer generation. Enantiomers are mirror-image pairs of atoms, where one of the molecules (chiral configurations) cannot be superimposed to the other through any rotation and translation (see Fig 1). They are present in any molecule with a tetrahedral carbon configuration with four unique groups. In the popular drug-like dataset GEOM, about 27% of the molecules contain such chiral configurations, and they can be crucial in the potency and safety of drugs. In ibuprofen, albuterol, and naproxen, one of the two enantiomers was shown to be more effective, while in thalidomide and propoxyphene, the incorrect enantiomer is toxic and can lead to severe side effects (Smith, ).

Most prior work that generates molecular conformations focuses on methods that are invariant to the special orthogonal group of transformations in three dimensions, denoted SO(3). This means that any rotation $R$ multiplied with a vector of atom coordinates $\mathbf{m}$ will leave their corresponding NN activations unchanged. In this work, we discuss and show why these methods are, in fact, invariant to any orthogonal group transformation, of which SO(3) is a subset, and how this leads to unwanted invariance wrt different chiral configurations of molecules.

Mirror-imaging a molecule may lead to different chiral molecules for certain tetrahedral configurations and can be achieved by multiplying its atom positions with a reflection matrix $R$, which is orthogonal and satisfies $\operatorname{det}R=-1$. Expressing pairwise differences as vectors $\mathbf{d}=\mathbf{x}_{i}-\mathbf{x}_{j}$, we note that Euclidean distance will therefore be invariant to reflections, since $||\mathbf{d}||_{2}=||R\mathbf{d}||_{2}$, for any orthogonal matrix $R$. Neural networks that only receive pairwise atom Euclidean distances will be unable to distinguish chiral configurations. We conduct experiments that also empirically attest to this fact, and illustrate how the problem cannot be trivially solved in the context of point-cloud formulations, even when removing the SO(3) invariance of the NN. Conversely, our field-based method perfectly captures the chiral distribution in the simplified chiral experiment we design.

We devise a three-dimensional multi field-based molecule representation, where atoms and bonds are represented as isotropic radial basis functions. Previous work using similar input representations (Pinheiro et al., ; Ragoza et al., ) has shown that working in this input space is challenging. Namely, the methods rely on evaluating the RBF fields at specific grid locations, creating very large tensors, compared to point-cloud representations. As such, GNN or transformer-based architectures are computationally prohibitive, and CNNs are the only viable option. While CNNs only have a translation invariance inductive bias, we show that defining simple, yet effective, reference rotations for molecules enables powerful U-Net architecture to achieve very good results, while preserving correct chirality.

Figure: Figure 1: Relative bond angles and distances are not sufficient to represent chirality. Top: A chiral molecule can not be rotated to its mirror image. Bottom: Chiral S-ibuprofen is an efficient COX-inhibitor resulting in pain and inflammation relief, while the mirror R-ibuprofen is not (Evans, ).
Refer to caption: x1.png

## 2 Related works

A large number of generative models to represent molecules have been proposed, which we distinguish primarily by their molecular representation instead of generative method families.

#### SMILES

SMILES is a string-based representation in which molecular graphs can be sufficiently encoded (Weininger, ), and has been employed in autoregressive models (Segler et al., ) and variational autoencoder (VAE) based generative models (Gómez-Bombarelli et al., ; Kusner et al., ). There are a few limitations when using smiles representations directly. First, similar molecular graphs may determine very different smiles encodings (e.g., measured by Levenshtein distance). Second, there can be multiple different smile representations for the same molecule, and the atom ordering in the string does not reflect the true, physical distance between the atoms.

#### Graphs

An alternative molecule representation consists of molecular graphs, where vertices encode atom type information, and their covalent bonds are represented through adjacency matrices. VAEs (Liu et al., ; Jin et al., ), flow-based models (Shi et al., ), and reinforcement learning (Shi et al., ; You et al., ; Zhou et al., ) methods were developed. The main drawback of only generating molecular graphs is that the stereochemical properties of molecules cannot be captured, and the functionality of various drug-like molecules depends on their three-dimensional conformation.

#### Point clouds

The three-dimensional configuration of a molecule can be encoded as $\mathbb{R}^{3}$ point clouds of atoms. Molecules reside in the $\mathbb{R}^{3}$ space, and their conformations are symmetric w.r.t. global translations and rotations. Inspired by this, methods that adhere to the symmetries of the Euclidean group in three dimensions (E(3) equivariant) have been developed. Autoregressive models (Gebauer et al., , ) have been developed, with the notable downside that they introduce an artificial atom order. GNN parametrized denoising diffusion models (Guan et al., ; Peng et al., ; Vignac et al., ; Xu et al., ; Wu et al., ) have also been employed, which, compared to their autoregressive counterparts, remove the arbitrary order of atom generation. Although rotational invariance is unquestionably well motivated, these methods also encode the undesired reflection invariance, which is part of the E(3) symmetries.

#### Fields

The $\mathbb{R}^{3}$ coordinates of atoms in molecules can be used as the mean positions of Gaussian or other densities, which fill the three-dimensional space. These densities can be interpreted as fields, where each molecule is represented by a function on $\mathbb{R}^{3}$, with values given by the densities of the atoms. Parameterizing the generative models with CNNs, we crucially distinguish this approach from the previously developed methods that use point cloud formulations by the E(3) symmetries they adhere to. Namely, our approach only respects translation invariance, while the point cloud approaches respect all E(3) symmetries. At the expense of losing the rotational invariance, this allows us to generate chiral conformations according to data statistics, which cannot be captured by reflection invariant models. This input representation has been previously modeled using VAE (Ragoza et al., ) and score-based models (Pinheiro et al., ), with some degree of success.

## 3 Molecular field diffusion

In this Section, we propose a diffusion model to generate molecular fields.

### 3.1 Field representation

We consider a molecule as a set of atoms $\{a_{i},\mathbf{m}_{i}\}$ and a set of bonds $\{b_{ij},\mathbf{m}_{ij}\}$, with atom types $a_{i}\in\mathcal{A}=\{\texttt{H},\texttt{C},\texttt{N},\texttt{O},\texttt{P},
\texttt{S},\ldots\}$, bond types $b_{ij}\in\mathcal{B}=\{b_{1},b_{2},b_{3}\}$, and we assume atom locations $\mathbf{m}_{i}\in\mathbb{R}^{3}$ and mid-point bond locations $\mathbf{m}_{ij}=\frac{\mathbf{m}_{i}+\mathbf{m}_{j}}{2}$ in Ångstroms. We represent a molecule using a collection of continuous 3D fields, one for each atom and bond type. We define the field for atom type $a$, $u_{a}$, as a mixture of RBF functions positioned on the atoms’ centers

$$ $\displaystyle u_{a}(\mathbf{x})=\gamma_{a}\sum_{i=1}^{N}\mathds{1}[a_{i}=a] \exp\left(-\frac{1}{2\sigma_{a}^{2}}||\mathbf{x}-\mathbf{m}_{i}||^{2}\right),$ (1) $$

where $N$ is the number of atoms, $\mathds{1}[\cdot]$ is the indicator functions, the variance $\sigma_{a}^{2}$ represents the size of atom type $a$, and $\gamma_{a}$ normalises the channels to unit interval $u_{a}(\mathbf{x})\in[0,1]$ (see Appendix B for field parameter values).
We create bond fields similarly

$$ $\displaystyle u_{b}(\mathbf{x})$ $\displaystyle=\gamma_{b}\sum_{(i,j)}\mathds{1}[b_{ij}=b]\exp\left(-\frac{1}{2 \sigma_{b}}||\mathbf{x}-\mathbf{m}_{ij}||^{2}\right).$ (2) $$

We evaluate a molecule’s atom and bond fields $u_{a},u_{b}$ at a fixed set $\mathcal{G}$ of $H\times W\times D$ locations, and create a vector-valued field based on which multi-channel 4D tensors are derived

$$ $\displaystyle\mathbf{u}$ $\displaystyle=(\mathbf{u}_{\texttt{H}},\mathbf{u}_{\texttt{C}},\ldots;\mathbf{ u}_{b_{1}},\mathbf{u}_{b_{2}},\mathbf{u}_{b_{3}})\in\mathbb{R}^{K\times H \times W\times D},$ (3) $$

where $K=|\mathcal{A}|+|\mathcal{B}|$ is the number of channels. We use tensor $\mathbf{u}$ as the input molecule representation for the diffusion model (See Figure 2).

Figure: Figure 2: The hydrogen channel voxel grid values $\mathbf{u}_{\text{H}}$ s.t., $\mathbf{u}_{\text{H}}>0.1$, are visualized. All atom positions and bonds are depicted through circles and line plots, respectively, for reference.
Refer to caption: x2.png

Figure: Figure 3: The faster atom schedule in the forward process results in the model first generating the molecule structure through its bonds, followed by atoms.
Refer to caption: x3.png

### 3.2 Diffusion model

We consider a diffusion generative model to generate molecular fields $\mathbf{u}$. In diffusion models, the data is incrementally noised by a forward process, and our goal is to learn a denoising reverse process parameterised by $\theta$. We consider the field tensors $\mathbf{u}_{t}$ over time $t\in[0,1]$, where we assume the observed or generated tensors are $\mathbf{u}_{0}$ at time $t=0$.

We follow the denoising diffusion probabilistic model (DDPM) proposed in (Sohl-Dickstein et al., ; Ho et al., ; Kingma et al., ), and refer the reader to the original publications for full derivations. We note that our tensor fields are effectively multi-channel 3D images, and thus image-based diffusion models translate to molecular voxel grid fields in a straightforward manner.

#### Forward process

We begin by assuming a Gaussian forward process

$$ $q(\mathbf{u}_{t}|\mathbf{u}_{0})=\mathcal{N}(\mathbf{u}_{t}|\alpha_{t}\mathbf{ u}_{0},\sigma_{t}^{2}\mathbf{I}),$ (4) $$

where $\alpha_{t}\in[0,1]$ monotonically reduces the original signal $\mathbf{u}_{0}$, while $\sigma_{t}\in[0,1]$ monotonically increases noise. We choose a variance-preserving noise specification, where

$$ $\displaystyle\alpha_{t}=\sqrt{1-\sigma_{t}^{2}}$ (5) $$

At $t=1$ we obtain pure noise $q(\mathbf{u}_{1}|\mathbf{u}_{0})=\mathcal{N}(\mathbf{0},\mathbf{I})$.

#### Forward posterior

The forward process admits a tractable conditional posterior distribution for an intermediate state $\mathbf{u}_{s}$ if we know both the original state $\mathbf{u}_{0}$ and further noised state $\mathbf{u}_{t}$ for $s<t$,

$$ $\displaystyle q\big{(}\mathbf{u}_{s}|\mathbf{u}_{t},\mathbf{u}_{0})$ $\displaystyle=\mathcal{N}\left(\frac{\alpha_{ts}\sigma^{2}_{s}}{\sigma^{2}_{t} }\mathbf{u}_{t}+\frac{\alpha_{s}\sigma_{ts}^{2}}{\sigma_{t}^{2}}\mathbf{u}_{0} ,\frac{\sigma_{ts}^{2}\sigma_{s}^{2}}{\sigma_{t}^{2}}\mathbf{I}\right),$ (6) $$

where $\alpha_{ts}=\alpha_{t}/\alpha_{s}$ and $\sigma_{ts}^{2}=\sigma_{t}^{2}-\alpha^{2}_{ts}\sigma^{2}_{s}$ (See cf. Kingma et al. ()). Due to Gaussianity of the forward, the posterior interpolates linearly between the endpoints while adding variance.

#### Reverse process

Our goal is to learn a generative process that reverses the forward by matching their marginals. The reverse process can be framed as the posterior of the forward

$$ $\displaystyle p_{\theta}(\mathbf{u}_{s}|\mathbf{u}_{t})$ $\displaystyle=q\Big{(}\mathbf{u}_{s}\big{|}\mathbf{u}_{t},\mathbf{u}_{0}:=\hat {\mathbf{u}}_{\theta}(\mathbf{u}_{t},t)\Big{)},$ (7) $$

where we replace the observation conditioning $\mathbf{u}_{0}$ with a fully denoised prediction $\hat{\mathbf{u}}_{\theta}(\mathbf{u}_{t},t)$ from a neural network parameterised by $\theta$.

#### Parameterisation

In practice, it is often more efficient to train a neural network to predict the added noise $\boldsymbol{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ of

$$ $\displaystyle\mathbf{u}_{t}=\alpha_{t}\mathbf{u}_{0}+\sigma_{t}\boldsymbol{ \epsilon},$ (8) $$

than to predict the signal $\mathbf{u}_{0}$. We then train a neural network $\boldsymbol{\epsilon}_{\theta}(\mathbf{u}_{t})$ to predict the noise $\boldsymbol{\epsilon}$, and follow the above equations with a substitution

$$ $\displaystyle\hat{\mathbf{u}}_{\theta}=\mathbf{u}_{t}/\alpha_{t}-\sigma_{t} \hat{\boldsymbol{\epsilon}}_{\theta}(\mathbf{u}_{t},t)/\alpha_{t}.$ (9) $$

#### Loss

The diffusion model admits a variational lower bound (Kingma et al., ), but earlier works have argued for a simplified loss (Ho et al., )

$$ $\displaystyle\mathcal{L}_{\mathrm{simple}}$ $\displaystyle=\mathbb{E}_{\mathbf{u},t,\boldsymbol{\epsilon}}\Big{[}|| \boldsymbol{\epsilon}-\hat{\boldsymbol{\epsilon}}_{\theta}(\mathbf{u}_{t},t)|| ^{2}\Big{]},$ (10) $$

where we learn to denoise arbitrary combinations of observations $\mathbf{u}$, noises $\boldsymbol{\epsilon}$ and noise scales $t\sim\mathrm{U}(0,1)$. We also opt for the simple loss $\mathcal{L}_{\mathrm{simple}}$.

### 3.3 Noise schedules

The cosine schedule for the forward process variance $\sigma_{t}^{2}$ has been shown to increase the negative log-likelihood and image quality (Nichol & Dhariwal, ). The effect is that data points change at a slower pace between consecutive time steps compared to the linear scheduler proposed by Ho et al. (). To further emphasize the stability of generated molecules, Vignac et al. () and Peng et al. () employed separate noise schedulers for the bonds and atom types, using a cosine scheduler with an additional parameter, and a new, sigmoid scheduler, respectively. Interestingly, these methods reported similar improvements while using reverse noise process order, where the bond type noising process is faster than atoms’ in Vignac et al. (), and vice-versa in Peng et al. ().

We achieved better results using the approach of Vignac et al. (), and therefore employ a model that first generates the bonds (molecular structure), and then fills in the atoms. Their proposed noise scheduler is:

$$ $\alpha_{t}=\cos\left(\frac{\pi}{2}\frac{(t+s)^{\nu}}{1+s}\right)^{2},$ (11) $$

where $\alpha_{t}$ determines the signal amplitude for a given time step $t$ in Equation 4, and we set $\nu$ to $1$ and $1.5$ for atom and bond channels, respectively. This induces a tendency to generate bonds earlier than atoms. In Figure 3, $\alpha_{t}$ values are plotted for atom and bond fields, respectively, along with an illustration of how bond RBF components $\mathbf{u}_{b}$ are generated before atoms $\mathbf{u}_{a}$, over time steps $t\in[0,1]$, for one sample generated by FMG (please refer to Appendix G for details on how this figure is created).

### 3.4 Classifier free guidance

Classification guidance in denoising diffusion models was introduced as a technique to control the tradeoff between the sample diversity and quality (Dhariwal & Nichol, ). This method requires training a separate classifier for a feature or label $\mathbf{y}$, $p_{\phi}(\mathbf{y}|\mathbf{u}_{t})$, on corrupted data points $\mathbf{u}_{t}$ at arbitrary noise scales $t$. The gradients of the log probability $\nabla_{\mathbf{u}}\log p_{\phi}(\mathbf{y}|\mathbf{u}_{t})$ are then used to guide the data points $\mathbf{u}_{t}$ in the direction of label $\mathbf{y}$ region. Following this work, Ho & Salimans () empirically showed that using an approximation of the implicit classifier $p(\mathbf{y}|\mathbf{u}_{t})\propto p(\mathbf{u}_{t}|\mathbf{y})/p(\mathbf{u}_
{t})$, instead of a separately trained classifier $p_{\phi}(\mathbf{y}|\mathbf{u}_{t})$ achieves similar results. Concretely, given the conditional and unconditional noise estimators $\boldsymbol{\epsilon}_{\theta}(\mathbf{u}_{t};\mathbf{y})$ and $\boldsymbol{\epsilon}_{\theta}(\mathbf{u}_{t})$, they use the following linear interpolation during sampling:

$$ $\tilde{\boldsymbol{\epsilon}}_{\theta}(\mathbf{u}_{t};\mathbf{y})=(1+\beta) \boldsymbol{\epsilon}_{\theta}(\mathbf{u}_{t};\mathbf{y})-\beta\boldsymbol{ \epsilon}_{\theta}(\mathbf{u}_{t}),$ (12) $$

where we note that $\boldsymbol{\epsilon}_{\theta}(\mathbf{u}_{t};\mathbf{y})$ and $\boldsymbol{\epsilon}_{\theta}(\mathbf{u}_{t})$ are up to a constant equal to the conditional and marginal score estimates $\nabla_{\mathbf{u}}\log p_{\theta}(\mathbf{u}_{t}|\mathbf{y})$ and $\nabla_{\mathbf{u}}\log p_{\theta}(\mathbf{u}_{t})$.

We use classification-free guidance to control the generation on the number of atoms. As expected, we observe improvements in all metrics and we use this, as it provides a better comparison to point-cloud molecule generation methods, which always require pre-specifying the number of atoms (e.g. sampled according to training data statistics).

### 3.5 Architecture choice and rotation invariance

We parametrize the denoiser $\boldsymbol{\epsilon}_{\theta}$ with a three-dimensional U-Net architecture, which consists of ResNet blocks, scaling layers, and attention layers. The choice is well-suited for translational invariance, as convolutions are in principle invariant to translations. However, convolutions are not rotationally invariant.

We propose a simple technique that side-steps the issue by orienting each molecule wrt its principal components $\mathbf{V}=(\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{3})$ of their atom positions $\mathbf{M}=(\mathbf{m}_{1},\ldots,\mathbf{m}_{N})$. Oriented molecules $\tilde{\mathbf{M}}=\mathbf{V}^{T}\mathbf{M}$ have a canonical frame of reference $\mathbf{V}$, where the axis of the largest variation is placed on the horizontal axis, the second-largest on the vertical axis, and the smallest on the depth axis. This ensures a high overlap between the large values of fields from different molecules.

### 3.6 Extracting molecular graphs

To extract a molecular graph from the generated fields $\mathbf{u}=(\mathbf{u}_{\texttt{H}},\mathbf{u}_{\texttt{C}},\ldots,\mathbf{u}_
{b_{1}},\mathbf{u}_{b_{2}},\mathbf{u}_{b_{3}})$, we first initialize the number of atoms $N$, their positions $\mathbf{m}_{i}$, and their types $a$, using a computationally efficient peak picking algorithm which retrieves sets $A_{a}$ of atom positions $\mathbf{m}_{i}$, for each atom type. We then optimize the mean-squared error loss between the values of the discretized generated fields $\mathbf{u}_{a}$, and their corresponding RBF function values determined by $\mathbf{m}_{i}\in A_{a}$

$$ $\begin{split}\operatorname*{arg\,min}_{\{\mathbf{m}_{i}\}}:\sum_{a\in\mathcal{ A}}&\sum_{\mathbf{x}\in\mathcal{G}}\Big{[}\mathbf{u}_{a}(\mathbf{x})-\gamma_{a }\sum_{i=1}^{N}\mathds{1}[a_{i}=a]\\ &\quad\exp\Big{(}-\frac{1}{2\sigma_{a}^{2}}||\mathbf{x}-\mathbf{m}_{i}||^{2} \Big{)}\Big{]}^{2},\end{split}$ (13) $$

by matching the generated atom fields $\mathbf{u}_{a}$ against the RBF mixtures induced by $\{\mathbf{m}_{i}\}$ across the grid points $\mathcal{G}$.

Similarly, given the optimized atom positions $\{\mathbf{m}_{i}\}$, we initialize the bonds (see Appendix D for details).
Briefly, considering the atom positions $\mathbf{m}_{i}$ fixed, we determine the bonds by first considering atoms within lenient distances (determined based on data statistics, with large margins), and filtering for those that contain large field values $\mathbf{u}_{b}$ on any bond channel $b\in\mathcal{B}$, around the mean points $\mathbf{m}_{ij}$ determined by the (potentially) bonded atom positions $\mathbf{m}_{i}$ and $\mathbf{m}_{j}$. Denoting the initial bond sets for each channel as $B_{ij,b}$, we optimize the RBF component $\gamma_{ij}$ values, which have $B_{ij,b}\neq\varnothing$:

$$ $\begin{split}\operatorname*{arg\,min}_{\gamma_{ij}}\sum_{b\in\mathcal{B}}&\sum _{\mathbf{x}\in\mathcal{G}}\Big{(}\mathbf{u}_{b}(\mathbf{x})-\sum_{(i,j)} \gamma_{ij}\mathds{1}[b_{ij}=b]\\ &\exp\Big{(}-\frac{1}{2\sigma_{b}^{2}}||\mathbf{x}-\mathbf{m}_{ij}||^{2}\Big{) }\Big{)}^{2},\end{split}$ (14) $$

where the $\gamma_{ij}$ values corresponding to incorrect bonds should optimally be close to 0 0.

Figure 4 illustrates a bond that was initially extracted based on $B_{ij,b}\neq\varnothing$, but removed after optimizing and finding the low $\gamma_{ij}$ value. Please refer to Appendix D for a detailed explanation of the molecular graph extraction procedure.

Figure: Figure 4: The bond marked with “x” on the left is removed since its optimized $\gamma_{ij}$ value (right bar plot) is significantly lower than the rest. The illustration shows the bond type 1 values of a noisy field $\tilde{\mathbf{u}}_{b_{1}}=\mathbf{u}_{b_{1}}+\boldsymbol{\epsilon}$, where $\boldsymbol{\epsilon}$ is a tensor of the same dimension as $\mathbf{u}_{b_{1}}$, containing samples of independent uniform noise $\mathcal{U}(0,0.1)$. The bond is initially considered due to our very inclusive initial bond extraction procedure applied to a noisy field. The bond field values $\mathbf{u}_{b_{1}}$ above a threshold of $0.2$ are shown and encoded according to the color bar.
Refer to caption: x4.png

## 4 Evaluation

#### Basic molecule properties

We use atom and molecule neutrality, validity, and novelty, as defined in (Guan et al., ). We call the stability metric introduced in (Guan et al., ) neutrality, to clarify and distinguish from other work that considered stable atoms even those without a formal charge of 0 (Vignac et al., ). Atom neutrality is the percentage of atoms that have a formal charge of 0, while molecule neutrality is the percentage of molecules having all atoms with 0 charge. Although in practice generating H atoms explicitly can be argued to be redundant, we emphasize that this still serves as one of the best benchmarking tools. Unlike validity, this metric cannot be trivially increased by e.g., extracting only covalent bond type 1 between atoms, or other such trivial solutions that may increase one metric at the expense of other distribution properties (Guan et al., ).

Note how the neutrality metric is strictly appropriate only when generating bond types $\mathcal{B}=\{b_{1},b_{2},b_{3}\}$. A few methods, such as the model from (Vignac et al., ) and (Peng et al., ) opted for explicitly generating aromatic rings, which creates ambiguity in the neutrality and validity metrics we have used. Please refer to Appendix I for details about this problem.

#### Molecule graph quality

We compute atom TVa and bond TVb total variation between the datasets and generated distributions, and report the sum of the distances over all atom and bond types, respectively.

Next, we compute the average maximum similarity between all generated molecules and the test dataset. We compute the Morgan fingerprints (Zhong & Guan, ) of all the generated and test dataset molecular graphs and determine the maximum Tanimoto Similarity between each generated molecule and the test set, which we denote by MST.

#### Molecular conformation metrics

To assess the quality of the generated 3D molecular conformations, we employ the Wasserstein metrics introduced in (Vignac et al., ), and compute bond angle (BA) and bond length (BL) W1 distance between generated and data distributions.

Please refer to Appendix H for more details about the performance metrics.

## 5 Experiments and results

### 5.1 QM9

QM9 (Ramakrishnan et al., ) is a small molecule dataset, containing 134k molecules of at most 29 atoms, including hydrogen. All molecules in the data are described by their atom positions and computed geometric, energetic, electronic, and thermodynamic properties. We use the same splits defined in (Guan et al., ), where we split the data into 100k, 18k, and 10k molecules for training, validation, and testing, respectively.

Out of the 100k training molecules, we use 98522 molecules. The molecules were selected such that the 95% percentile of the RBF components corresponding to each atom $\mathbf{m}_{i}\pm 2\sigma$ would fit within the selected field evaluated positions $\mathcal{G}$.

In Table 1, we compare our model to a wide variety of previous methods. We achieve state-of-the-art performance in terms of molecule and atom neutrality, indicating that our method captures the basic chemical properties very well.

In Table 2, the maximum test similarity and the bond distribution are competitive to the point-cloud approaches, while MiDi performs very well on TVa. Since we lack the strong enforcement of the atom number for each generated molecule, and only use classification-free guidance to stochastically guide the model, we further analyze how well the generated molecules adhere to the number of atoms we condition the model on in Appendix F. Our generative method likely trades this property for higher molecular stability, by removing or replacing the RBF components in the atom fields $a$, such that the atoms’ formal charges are 0. We note how this is not possible for current point-cloud methods, as the number of atoms in each molecule is strictly enforced, and these methods can only change the atom types and locations that are initially randomly assigned, but not the number of atoms in the generated molecules.

Remarkably, even though our NN architecture is not invariant to SO(3) transformations, our method achieves satisfactory conformational accuracy, with a very low distance between generated and training data molecules. We report the W1 distance for bond distances and bond angles distributions of the data and generated molecules in Appendix A.1, Table 6, and the corresponding CDFs based on which W${}_{1}$ is computed in Figure 6.

**Table 1: Basic molecule properties for molecular graphs generated by models trained on the QM9 dataset. Results marked with (*) use deterministic, post-processing techniques to determine the bonds or assign formal charges to their atoms when computing the metrics. Results marked with ($\dagger$) have been re-tested by us, and those marked with $(\ddagger)$ generate an additional, explicit aromatic bond channel. We compute the mean and standard deviation of the metrics for three different runs, each generating 10k samples. All properties are expressed in percentages.**
|  | Neutrality |  |  |  |
| --- | --- | --- | --- | --- |
| Model | Atom | Mol | Validity | Nov (↑) |
| EDM($\dagger$) | 98.1 | 81.9 | 91.6 | 65.7 |
| GDM-AUG | 97.6 | 71.6 | 90.4 | 74.6 |
| EDM-Bridge | 98.8 | 84.6 | 92.0 | 65.7 |
| MolDiff(*, $\ddagger$) | - | - | 97.0 | - |
| MiDi($\dagger$) | 98.3 | 84.5 | 96.7 | 54.1 |
| G-Schnet | 95.7 | 68.1 | 85.5 | - |
| GeoLDM | 98.9 | 89.4 | 93.8 | - |
| VoxMol(*) | 99.2 | 89.3 | 98.7 | - |
| FMG | 99.4$\pm$ 0.1 | 91.5$\pm$ 1.4 | 98.8 $\pm$ 0.3 | 68.8 $\pm$ 0.3 |
| QM9 | 99.9 | 99.6 | 99.6 | - |

**Table 2: Molecule graph quality of generated samples from models trained on the QM9 dataset. We report mean and standard deviation over three runs, with 10k samples per model.**
| Model | MST (↑) | TVa(↓) | TVb(↓) |
| --- | --- | --- | --- |
| EDM($\dagger$) | 0.57 | 0.49 | 0.98 |
| MiDi($\dagger$) | 0.55 | 0.19 | 0.81 |
| FMG | 0.54$\pm$ 0.00 | 0.52 $\pm$ 0.02 | 0.95 $\pm$ 0.01 |

### 5.2 QM9 conditional generation

We employ the same setting as in previous work (Guan et al., ), where we split molecular properties in one thousand equal bins, and form, based on training data statistics, the conditional and marginal categorical distributions $p(N_{a}|\mathbf{c}_{i})$, and $p(\mathbf{c}_{i})$, where $\mathbf{c}_{i}$ is the $i^{\text{th}}$ property bin and $N_{a}$ is the number of atoms. We test our conditional model by sampling two thousand pairs of $[\mathbf{c}_{i},N_{a}]\sim p(\mathbf{c}_{i})p(N_{a}|\mathbf{c}_{i})$, and further sample the actual conditional property $c\sim\mathcal{U}(\mathbf{c}_{i,0},\mathbf{c}_{i,1})$, with which we condition the generative model on $\boldsymbol{\epsilon}_{\theta}(\mathbf{f}_{t},t,N_{a},c)$.

The results and further details of the experiment are shown in Appendix A.2, where we compare the difference between the highest and lowest occupied molecular energies $\Delta\epsilon$ and the polarizability $\alpha$ of the generated molecules, compared to the true, conditional values. FMG outperforms EDM in terms of molecular stability, while EDM is better in its fidelity to the conditioning variables $c$.

Additionally, our method generates similar fidelity to the conditioning variable $c$, even without classification-free guidance on $N_{a}\sim p(N_{a}|c)$. This could prove useful for future situations where the conditional $p(N_{a}|c)$, cannot be accurately determined based on data statistics (e.g., when conditioning on multiple properties $c_{1},c_{2}$, where data may be insufficient to determine $p(N_{a}|c_{1},c_{2})$ accurately).

### 5.3 GEOM

The Geometric Ensemble Of Molecules (GEOM) dataset (Axelrod & Gómez-Bombarelli, ) contains over 450 thousand molecules of up to 181 atoms and 37 million conformations along with their corresponding energies. Following previous work (Guan et al., ; Peng et al., ), we select the 30 lowest energy conformers for each molecule graph. Similar to (Peng et al., ), we also restrict the data points to only those molecules that contain atoms from {H, C, N, O, F, P, S, Cl}. Finally, we extract those molecules that can be contained in the pre-defined evaluated positions that determine the tensors $\mathbf{u}_{a}$ and $\mathbf{u}_{b}$, similar to our QM9 experiments. After filtering, the data contains 96.7% of the training molecular conformations and 99.7% unique molecular graph coverage.

The stability results can be found in Table 3, where our method achieves highly competitive neutral atoms and molecule percentages (please refer to our explanations on the challenges of benchmarking against methods that generate explicit aromatic bonds from Appendix I). In Table 4, we see the same effect mentioned in the QM9 experiments, where the method has a worse TV${}_{a}$ metric, likely as a cause of our method favoring stability instead. On the contrary, when we do not generate H atoms explicitly, our generated atom distribution becomes significantly closer to the data.

Finally, we analyze the conformational properties of our method in Appendix A.1. Even on this significantly more complex dataset, our method generates most molecular geometries very well.

**Table 3: Basic molecule properties for molecular graphs generated by models trained on the GEOM dataset. We report the mean result of 10k generated samples, for both explicit and implicit H models. An H checkmark means the method explicitly generates the H atoms. All properties are expressed in percentages.**
|  |  | Neutrality (↑) |  |  |  |
| --- | --- | --- | --- | --- | --- |
| Model | H | Atom | Molecule | Val | Nov(↑) |
| EDM($\dagger$) | ✓ | 81.0 | 0.5 | 3.0 | 100 |
| MiDi($\dagger,\ddagger$) | ✓ | 96.8 | 31.3 | 62.7 | 100 |
| GeoLDM | ✓ | 84.4 | - | 99.3 | - |
| VoxMol(*) | ✓ | 98.1 | 75.0 | 93.4 | - |
| FMG | ✓ | 94.7 | 16.5 | 64.9 | 100 |
| GEOM | ✓ | 99.4 | 86.2 | 89.6 | - |
| MolDiff($\dagger$, $\ddagger$) |  | - | - | 85.0 | 100 |
| FMG |  | - | - | 80.8 | 100 |

**Table 4: Molecule graph quality of generated samples from models trained on the GEOM dataset. We report the mean result of 10k generated samples, for both explicit and implicit H models.**
| Model | H | MST(↑) | TVa(↓) | TVb(↓) |
| --- | --- | --- | --- | --- |
| EDM($\dagger$) | ✓ | 0.46 | 1.51 | 5.72 |
| MiDi($\dagger$,$\ddagger$) | ✓ | 0.53 | 0.76 | - |
| FMG | ✓ | 0.49 | 2.3 | 4.78 |
| MolDiff($\dagger$, $\ddagger$) |  | 0.62 | 2.03 | - |
| FMG |  | 0.50 | 0.64 | 4.58 |

### 5.4 Chirality analysis

We aim to also empirically attest that our proposed method can capture all molecule chiral properties, while current implementations of SO(3) invariant NN lead to invariance w.r.t. enantiomers. We create a simplified experiment, where we attempt to overfit on 48 unique molecules having a carbon tetrahedral configuration connected to O, N, H, and C. We generated the molecule atom positions using RDKit, and analyzed the distribution of the determinants of the atom positions $\mathbf{m}$ of O, N, and H.

We consider matrices defined by $R=[\mathbf{m}_{\texttt{O}}-\mathbf{m}_{\texttt{C}},\mathbf{m}_{\texttt{N}}-
\mathbf{m}_{\texttt{C}},\mathbf{m}_{\texttt{H}}-\mathbf{m}_{\texttt{C}}]^{T}
\in\mathbb{R}^{3\times 3}$, where $\mathbf{m}_{C}$ is the central carbon atom of the tetrahedral configuration, and the triplets $(\mathbf{m}_{O},\mathbf{m}_{C},\mathbf{m}_{H})$ are the positions of three of the four atoms it connects to. The enantiomer configuration of the generated and the data molecules can be uniquely identified by the sign of $\operatorname{det}R$, while the volume between the four atoms is defined by its absolute value $|\operatorname{det}R|$.
In Figure 5, the determinant values of molecules generated by our method match those found in the data, while the equivariant method (EDM) has any of the two enantiomers equally likely. We further assess the SO(3) variant method defined by (Guan et al., ), which concatenates atom positions to its GNN node states, that should, in theory, be able to capture enantiomer states. Surprisingly, although the method now clearly starts recognizing the two different enantiomers, it fails to produce realistic conformations and generates many molecules with approximately 0 volume, i.e., the central tetrahedral carbon C and the three connecting atoms, N, H, and O, are approximately coplanar.

Please refer to Appendix A.3 for details of this experiment and theoretical analysis of the current SO(3) implementations’ behavior related to enantiomer configurations.

Figure: Figure 5: Kernel density estimation of the triplet vector $[\mathbf{m}_{\texttt{O}}$ - $\mathbf{m}_{\texttt{C}}$, $\mathbf{m}_{\texttt{N}}$ - $\mathbf{m}_{\texttt{C}}$, $\mathbf{m}_{\texttt{H}}$-$\mathbf{m}_{\texttt{C}}]^{T}$ determinants, resulted from the data and generated molecules. EDM is unable to distinguish and correctly generate the enantiomer distribution, while the SO(3) variant method (GNN) generates improper, 0 volume conformations. FMG matches the enantiomer distribution.
Refer to caption: x5.png

### 5.5 Ablation

We aim to determine which aspects of the method contribute the most to our competitive molecule generation performance. In Table 5, FMGrot receives fields based on randomly rotated sets of atom coordinates $\mathbf{m}_{i}$ during training, FMG${}_{\beta}$, $\beta\in\{0,2,5\}$ denotes the original model with various degrees of classification free guidance on the number of atoms (see Section 3.4), and FMGdist is the performance of the molecular graphs, where bonds are extracted according to bond pairwise atom distances. The reference alignment is the main factor that increases the neutrality metric to state-of-the-art performance, as FMGrot performance drops significantly. Classification guidance on the number of atoms also has a clear effect on TVa, which supports our hypothesis that having a fixed, pre-specified number of atoms in point-cloud based methods is likely the leading cause of better atom type distribution, at the expense of molecule stability.

For FMGdist, we use the same lookup table from (Guan et al., ) to determine bonds based on atom distances. In Table 5, the performance of FMGdist is remarkably good. This further confirms that our method generates very good conformations (having bond distances per type within acceptable margins) since it would otherwise be impossible for this method to have such high atom and molecule stability.

**Table 5: FMG performance with various design choices. All basic molecular properties (Val and Neutrality Atom and Mol) are expressed in percentages.**
|  |  | Neutrality |  |  |  |
| --- | --- | --- | --- | --- | --- |
| Model | Val | Atom | Mol | TVa(↓) | TVb(↓) |
| EDM | 91.6 | 98.1 | 81.9 | 0.49 | 0.98 |
| FMGrot | 95.8 | 98.6 | 82.6 | 0.69 | 0.94 |
| FMG0 | 98.6 | 99.3 | 92.2 | 1.02 | 1.16 |
| FMG2 | 98.8 | 99.4 | 91.5 | 0.52 | 0.95 |
| FMG5 | 98.6 | 99.4 | 92.2 | 0.48 | 1.02 |
| FMGdist | 97.7 | 99.1 | 90.9 | 0.54 | 1.00 |

## 6 Conclusion

In this work, we show how continuous field-based molecular representations can be successfully employed for drug-like molecule generation. We achieve highly competitive molecule stability and analyze enantiomers, a previously omitted molecular property. We show how our method can capture this molecular property, while point-cloud methods cannot, and fail even in the SO(3) variant context.

## 7 Acknowledgements

We acknowledge the computational resources provided by the Aalto Science-IT Project.

## Appendix A Additional experiments and details

### A.1 Conformational analysis

We report the Wasserstein distances between the generated and data bond angle distributions in Table 6. We calculated the Wasserstein W1 bond angle (BA) and bond length (BL) distances for each atom and bond type and respectively averaged over atom and bond types (please refer to Appndix H for details). In Figure 6, we also plot the CDF of the distributions, based on which the W1 distances were computed.

Note that for the GEOM dataset, MiDi uses explicitly generated aromatic bond types, and thus has a different bond type set $\mathcal{B}$ compared to ours. This also affects the distributions of the bond types 1 and 2 and we therefore omit MiDi from the GEOM bond distance metric BL.

**Table 6: Conformation results on GEOM and QM9 datasets, measured by the W${}_{1}$ distance between the generated and data molecules. For QM9, we report the mean and standard deviation over three models, with 10k samples from each model, and 10k samples from a single model for GEOM.**
|  | GEOM | QM9 |  |  |
| --- | --- | --- | --- | --- |
| Model | BL (↓W${}_{1}\cdot 10^{2}$ ) | BA (↓W${}_{1}$) | BL (↓W${}_{1}\cdot 10^{2}$ ) | BA (↓W${}_{1}$) |
| EDM($\dagger$) | 9.09 | 5.48 | 0.47 | 0.68 |
| MiDi($\dagger,\ddagger$) | - | 1.55 | 1.96 | 0.95 |
| FMG | 2.99 | 1.99 | 0.52 $\pm$ 0.02 | 0.95 $\pm$ 0.01 |

Figure: Figure 6: Cumulative distribution function plot of generated and training data for MiDi, EDM, and FMG for the GEOM and QM9 datasets. All models capture these conformational properties very well on QM9, and reasonably well on GEOM.
Refer to caption: x6.png

### A.2 Conditional generation

We train classifiers $\phi_{c}$ for the two considered properties, $\alpha$ and $\Delta\epsilon$. The property classifiers $\phi_{c}$ are E(3) invariant graph neural networks that receive pairwise distances between atoms $||\mathbf{m}_{i}-\mathbf{m}_{j}||$ and atom types as categorical features, very similar to the architecture used for parametrizing EDM. Please refer to (Guan et al., ) for more details.

We test our generative models by predicting the generated molecules’ properties, and computing the L1 loss between those predictions and the variable we used to condition the method during generation, $c$.

In Table 7, FMG is the trained conditional method for the two properties, respectively, and FMG${}_{\beta=0}$ is the same model, where we do not use classification-free guidances on the number of atoms according to the property $p(N|c)$.

We also include additional baselines: FMGuc is an unconditional model, where we sample molecules according to an appropriate number of atoms $p(N|c)$, and FMGrnd are completely random samples, from an unconditional model.

**Table 7: Conditional generation results. Atom and Mol Neutrality are reported in percentages. The mean results of 2k samples generated by our methods are shown.**
|  | $\alpha$ | $\Delta\epsilon$ |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | Neutrality |  | Neutrality |  |  |  |
| Model | L1(↓) | Atom | Mol | L1 (↓) | Atom | Mol | Conditional |
| EDM | 2.76 | - | 80.4 | 655 | - | 81.7 | ✓ |
| FMG | 3.55 | 99.0 | 86.8 | 746 | 98.9 | 84.5 | ✓ |
| FMG${}_{\beta=0}$ | 3.92 | 99.0 | 86.7 | 734 | 98.6 | 83.0 | ✓ |
| FMG${}_{\text{uc}}$ | 6.92 | 99.4 | 91.5 | 1259 | 99.4 | 91.5 |  |
| FMG${}_{\text{rnd}}$ | 9.03 | 99.4 | 91.5 | 1451 | 99.4 | 91.5 |  |

We implement property conditioning by shift-scaling the channels, as presented in (Perez et al., ), for all U-Net layers (similar to time and number of atoms conditioning). Note that we did not use classifier-free guidance for molecular properties $c$.

### A.3 Chirality Analysis

#### Experiment

We create a dataset of 48 molecules, each one having a central tetrahedral carbon that is always connected to the set of four atoms $\mathcal{C}=\{\text{{C}, {H}, {N}, {O}}\}$. We create multiple unique molecules, by connecting the three atoms $\text{{C},{N},{O}}\in\mathcal{C}$ to other arbitrary residues. In Figure 7, we show half of the molecular graphs within our data. For each of them, we create molecular conformations using RDKit. Note how in this experiment we are not necessarily interested in creating realistic molecules, but simply any molecules that contain chiral configurations.

Figure: Figure 7: Molecular graphs of half of the chiral dataset.
Refer to caption: x7.png

We train FMG for 10k iterations, and EDM for 300k iterations, and test it for multiple checkpoints. Since the dataset is small, we do not use mini-batches and set our batch size to the data length.

#### Theoretical analysis

Let $\mathbf{v}_{1},\mathbf{v}_{2}\in\mathbb{R}^{3}$ be two, three-dimensional vectors, and let $\mathbf{d}=\mathbf{v}_{1}-\mathbf{v}_{2}$. By definition, orthogonal matrices $Q$ have an inverse equal to $Q^{-1}=Q^{T}$. Therefore, it is easy to see that $||Q\mathbf{d}||=||\mathbf{d}||$ since

$$ $||Q\mathbf{d}||^{2}=(Q\mathbf{d})^{T}(Q\mathbf{d})=\mathbf{d}^{T}Q^{T}Q\mathbf {d}=\mathbf{d}^{T}\mathbf{I}\mathbf{d}=\mathbf{d}^{T}\mathbf{d}=||\mathbf{d}|| ^{2},$ (15) $$

where $\mathbf{I}$ is the identity matrix, and $||\cdot||$ is the L2 norm. Since reflection and rotation matrices are orthogonal, it follows, that the Euclidean distances between any two vectors $\mathbf{v}_{1},\mathbf{v}_{2}$ are invariant to orthogonal matrix multiplication, and, therefore, to any rotation or reflection matrix.

Consider now the cosine of the angle between any two vectors $\mathbf{v}_{1},\mathbf{v}_{2}\in\mathbb{R}^{3}$, transformed by an orthogonal matrix $Q$. We have that

$$ $\cos(Q\mathbf{v}_{1},Q\mathbf{v}_{2})=\frac{Q\mathbf{v}_{1}\cdot Q\mathbf{v}_{ 2}}{||Q\mathbf{v}_{1}||||Q\mathbf{v}_{2}||}=\frac{(Q\mathbf{v}_{1})^{T}(Q \mathbf{v}_{2})}{||\mathbf{v}_{1}||||\mathbf{v}_{2}||}=\frac{\mathbf{v}_{1} \mathbf{I}\mathbf{v}_{2}}{||\mathbf{v}_{1}||||\mathbf{v}_{2}||}=\frac{\mathbf{ v}_{1}\cdot\mathbf{v}_{2}}{||\mathbf{v}_{1}||||\mathbf{v}_{2}||}=\cos(\mathbf{ v}_{1},\mathbf{v}_{2}),$ (16) $$

where we have used equation 15 for $||Q\mathbf{v}_{1}||=||\mathbf{v}_{1}||$. All distances and angles between all vectors in the $\mathbb{R}^{3}$ space are therefore preserved.

Any function $\phi$ (e.g., a neural network) that receives pairwise atom distances or angles will therefore be invariant to orthogonal matrix multiplication and implicitly, to reflections.

Figure: Figure 8: R and S configurations of a general tetrahedral configuration with three, distinct residues.
Refer to caption: x8.png

Consider some tetrahedral configuration, where the central carbon atom is positioned at $\mathbf{m}_{C}$, and is connected to four unique residues $R1,R2,R3$, and $R4$, initially positioned at $\mathbf{m}_{R1},\mathbf{m}_{R2},\mathbf{m}_{R3}$, and $\mathbf{m}_{R4}$, respectively (see Figure 8). Observing the molecule along the $R_{1}$-C-axis, a way of defining the chirality of the molecule is to specify the ordering of the remaining residues in clock-wise or counter-clockwise orientation. Assuming we start with $R2$ position fixed at $\mathbf{m}_{R2}$, swapping $R3$ and $R4$ positions has the effect of reversing the rotation direction of $R3$ and $R4$, and creating a distinct chiral configuration. Now, considering the plane defined by $m_{R1}$, the central carbon position $\mathbf{m}_{C}$, and $\mathbf{m}_{R2}$, the same effect can be achieved by a reflection about this plane.

In general, a reflection about any plane can be achieved by a specific translation, a rotation, and a reflection, applied to the whole $\mathbb{R}^{3}$ space (after reflection, the molecule may be moved back with an opposite translation and rotation). As Euclidean distances and angles between any two vectors are preserved under each of these operations, it follows that Euclidean distances and angles become invariant to the chirality.

To our knowledge, all current point-cloud-based methods use Euclidean distances or relative angles between atom position vectors $\mathbf{m}_{i},\mathbf{m}_{j}$, e.g., for message passing in their GNN architectures. These methods cannot distinguish enantiomers, and, as such, will always produce them completely at random, and the molecular graphs will not change based on whether the S or R enantiomer configuration is present in that molecule.

## Appendix B Field parameters

For the QM9 dataset, we determined that a resolution of 0.33Å and $32^{3}$ voxels for each molecule is a good compromise between the sample size and the percentage of the dataset’s samples (the total size being $\mathbf{u}\in\mathbb{R}^{8\times 32\times 32\times 32}$, for 5 and 3 atom and bond channels, respectively). We set the RBF component’s $\sigma=0.224$, for all fields, and $\gamma=0.176$. The component values $\gamma=0.176$ were determined as the maximum value of an RBF component having its positions $\mathbf{m}$ in one of the evaluated positions $\mathcal{G}$. Although the field parameters do not reflect the actual, physical size of the atoms, we note how true physical fields may be recovered by appropriately scaling the variance of each atom RBF component, based on the channels $a$ they belong to.

We used the same resolution for GEOM. After rotating all molecules using the principle components of the atom position vectors, the resulting tensors $\mathbf{u}\in\mathbb{R}^{11\times 64\times 40\times 32}$ were enough to cover 96.66% of the total data points and 99.71% unique molecular graphs.

## Appendix C Model and optimization parameters

#### Diffusion Model

We use 100 timesteps $t\in[0,1]$ for our diffusion processes, the $\boldsymbol{\epsilon}$ parametrization, and the modified cosine scheduler for $\alpha_{t}$, as presented in main Section 3. The conditional scale for the classifier-free guidance was kept to $\beta=2$ during the generation phase, for our main results (unless explicitly mentioned otherwise).

#### Optimization

We optimize the $\mathcal{L}_{\mathrm{simple}}$ objective using Adam optimizer with $8\cdot 10^{-5}$ learning rate and (0.9, 0.99) for the Adam’s $\beta$ parameters. The models are trained using batch sizes of 32 and 64, for GEOM and QM9 datasets, respectively.

#### U-Net

We use a three-dimensional adaptation of the U-Net employed by (Ho et al., ). For both GEOM and QM9 datasets, we use the same hyperparameters, except for the batch size. The architecture consists of three U-Net layers, with linear attention for each resolution level, except the last one, which uses full attention. Its initial convolution converts the 8 and 11 channels (for QM9 and GEOM datasets, respectively), into a hidden representation of $\mathbb{R}^{128\times H\times W\times D}$, which is halved for each of the $H,W,D$ dimensions, while its channel size is multiplied by 1, 2, and 3.

#### Resources

For the QM9 experiments, we use four A100 GPUs (40GB memory version) for 180 hours, and 1.2 million iterations (or about 780 epochs). On the GEOM dataset, we used the same four A100 GPUs for 330 hours and trained our explicit H GEOM model for 1.32 million iterations (6 epochs), and our implicit one for 1.2 million iterations (5.5 epochs).

## Appendix D Molecule graph extraction parameters and algorithms

The molecular graph extraction has four components, used in the following order: deterministic peak extraction, atom positions optimization, deterministic bond extraction, and $\gamma_{ij}$ parameter optimization that removes false bonds.

#### Deterministic peak extraction

The algorithm works by defining a simple, recursive flood-fill algorithm for retrieving the neighbors. Assuming the RBF functions have a small $\sigma$ parameter and reasonably small amounts of noise produced by the diffusion model, it will never fail to extract the correct number of atoms and their approximate positions from the fields.

Figure: Algorithm 1 Peak extraction

#### Optimizing the atom positions

The initial atom positions retrieved in the $A_{\mathbf{m}}$ sets, are further optimized using Equation 13. Note how both atom numbers and their types are obtained using the peak extraction algorithm, and the optimization only fine-tunes their positions.

#### Extracting bonds

The distances between any two bonds $b_{ij}$ and $b_{jk}$ are half as large as the distances between the atoms $a_{i}$ and $a_{k}$ ($||\mathbf{m}_{ij}-\mathbf{m}_{jk}||_{2}=||\mathbf{m}_{i}-\mathbf{m}_{k}||_{2}/2$). This may create a significantly larger overlap between the bond RBF components, compared to atoms. Therefore, we treat the bond assignment differently and do not reuse the same peak extraction algorithm. Using the optimized atom positions $\mathbf{M}_{a}$, we check for all the atom pairs $\mathbf{m}_{i},\mathbf{m}_{j}\in\mathbf{M}_{a}$ that are within some predetermined threshold distance $||\mathbf{m}_{i}-\mathbf{m}_{j}||_{2}<d_{1}$ and determine $B_{ij,b}=\{u_{b}(\mathbf{x});||\mathbf{x}-\mathbf{m}_{ij}||_{2}<d_{2},u_{b}(
\mathbf{x})>r_{2}\}$. A covalent bond type $b$ is added between atoms $i$ and $j$ if $B_{ij,b}\neq\varnothing$.

To ensure we don’t miss bonds due to the low resolution, we select a large sphere (determined by its radius $d_{2}$) around the approximate bond positions $\mathbf{m}_{ij}$, within which we check for values greater than $r$. Additionally, since the true $\mathbf{m}_{ij}$ positions may land at the edge of a voxel, the maximum value corresponding to bonds $b_{ij}$ in the tensors $\mathbf{u}_{b}$ can be small, and therefore we chose a low threshold $r$. This may then recover bonds at positions where the model did not generate a true RBF bond component.

#### Optimizing γ i ⁢ j subscript 𝛾 𝑖 𝑗 \gamma_{ij} italic_γ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT

To remove erroneous bonds, extracted due to our lenient bond-extraction parameters, we optimize the $\gamma_{ij}$ assigned to each bond $b_{ij}$, assuming the mean bond positions $\mathbf{m}_{ij}$ fixed. Let $\gamma_{b}=\{\gamma_{ij};b_{ij}=b\}$ (where $b_{ij}$ is initially assigned based on the condition $B_{ij,b}\neq\varnothing$), then optimize the RBF components $\gamma_{b}$ using equation 14. In Figure 9, we analyze how robust the method is to various degrees of uniform noise. We choose 300 molecules from the QM9 dataset and gradually add more and more noise to the voxelized fields $\mathbf{u}$.

Figure: Figure 9: Performance of the graph extraction method for various degrees of noise. The bond optimization step is crucial and makes the method highly robust to noise.
Refer to caption: x9.png

#### Parameters

We choose threshold parameter $t=0.3$ in Algorithm 1. For bond detection, we choose a large distance $d_{1}$ of $0.35$Å (for reference, the margin used to detect bond type 1 based on distance in the ablation Table 5 is $0.1$Å). The distance $d_{2}=0.45$ Å will result in the method checking for the closest midpoint voxel to $\mathbf{m}_{ij}$ and a few others directly in contact with it, while the bond threshold is set to $r=0.3$. Both the mean positions $\mathbf{m}$ and bond RBF component factors $\gamma_{ij}$ are optimized using gradient descent for 500 iterations and a learning rate of $5$. All parameters were selected according to robustness analyses, such as the one illustrated in Figure 9, based on the training data.

## Appendix E Generated samples

Figure: Figure 10: Random valid samples from FMG trained on QM9.
Refer to caption: x10.png

Figure: Figure 11: Random valid samples from FMG trained on GEOM.
Refer to caption: x11.png

## Appendix F Classificaiton guidance analysis

We determine how accurately the classification free guidance results in molecules of the actual desired number of atoms. We compute the accuracy as the percentage of generated molecules in the correct conditioning bin $N$, the L1 distance between the binned number of atoms generated by the model $\tilde{N}$, and the conditioning bin $N$ as $\min(|\tilde{N}_{\text{high}}-N_{\text{low}}|,|N_{\text{high}}-\tilde{N}_{
\text{low}}|)$ (or 0, when $N=\tilde{N}$), and “L1 lower” is the fraction of the total L1 distance attributed to incorrectly generated bins because of a strictly lower number of atoms generated $\tilde{N}<N$ (as a measure of a left or right shift of the distribution).

In Figure 12, the generated bin accuracy is relatively low but is generally off by approximately 1 atom (measured by the L1 distance) for the QM9 dataset. In terms of the distribution shift, L1 lower suggests that the QM9 FMG model tends to generate more atoms than the conditioning variable.

Similar results in terms of accuracy and L1 distance can be noted in the GEOM dataset (relative to the much larger possible number of atoms), but now the shift is in the opposite direction. This could be an artifact of the number of training iterations: early on in the QM9 dataset training experiments, we observed that the model tends to generate a lower number of atoms and bonds, which likely explains this opposite shift of the GEOM method. For reference, our GEOM model has been trained for about 6 epochs, compared to its QM9 counterpart (about 780).

Figure: Figure 12: Histogram of the generated number of atom classes $\tilde{N}$ and the conditioning used $N$. Small and more significant atom number distribution shifts occur for QM9 and GEOM datasets, respectively. Together with the result of Table 5 which shows how critical it is to enforce the correct number of atoms for TVa, this result suggests that our method may remove atoms for stability, inducing an atom distribution shift.
Refer to caption: x12.png

## Appendix G Visualization

Figure 3 is a simplified illustration meant to conceptually show the noise scheduler’s effect, and it is perhaps not accurately representing an actual visualization of FMG generating a molecule (please refer to Figure 13 for such an illustration).

Figure: Figure 13: Three-dimensional visualization of a generated molecule by FMG on QM9. The upper and lower images depict the summed values (normalized to $[0,1]$ at each step) over all bond $\mathbf{u}_{b,1}+\mathbf{u}_{b,2}+\mathbf{u}_{b,3}$ and atom $\mathbf{u}_{\texttt{H}}+\mathbf{u}_{\texttt{C}}+\mathbf{u}_{\texttt{O}}+ \mathbf{u}_{\texttt{N}}+\mathbf{u}_{\texttt{F}}$ fields, respectively. The same effect as in the main paper Figure 3 can be seen, where bonds become visible through the noise quicker than the atoms. Circles and lines are added for reference, denoting atoms and bonds extracted at the last time step.
Refer to caption: extracted/5429350/appendix_figures/generation_mols.png

Figure 3 of the main text is done using the maximum $L_{2}$ loss between noisy $\mathbf{u}_{t}$ (for various $t\in[0,1]$) and the determined RBF function at step $t=0$. Using a similar loss function as in equation 13, we determine the fidelity of each timestep atom and bond fields ($\mathbf{u}_{b,t}$, $\mathbf{u}_{b,a}$) for atom $a$ and bond $b$ channels separately, to the final extracted RBF components based on $\mathbf{u}_{0}$. Using the maximum loss values over each channel separately, we vary the illustration’s transparency parameter for various bond and atom types over time $t$.

## Appendix H Performance metric details

#### Validity:

Percentage of valid molecules, computed using RDKit’s SanitizeMol function.

#### Neutrality:

Percentage of generated atoms that have a formal charge of 0, and percentage of molecules that contain all atoms with a formal charge of 0.

#### Novelty:

percentage of generated molecular graphs that are not present in the training data (after canonicalization).

Among our molecular quality metrics, the atom and bond total variation distances TVa and TVb are a good measure of how well most of the data modes are captured, whereas MST should better measure the quality of the molecules. For instance, if the method captures a single mode of distribution well, and generates very high-quality molecules, it may still have relatively bad TVa and TVb.

#### TV a , TV b

We define TVa and TVb as:

$$ $\text{TV}_{\text{a}}=\sum_{a\in\mathcal{A}}\sum_{N_{a}=i}^{N}|p_{d}(N_{a})-p_{ g}(N_{a})|$ $$

$$ $\text{TV}_{\text{b}}=\sum_{b\in\mathcal{B}}\sum_{N_{b}=i}^{N}|p_{d}(N_{b})-p_{ g}(N_{b})|,$ $$

where $N_{a},N_{b}$ are the numbers of atoms and bonds of type $b$, $a$ and $b$ are the atom and bond types, $\mathcal{A}$ and $\mathcal{B}$ are the sets of all unique atoms and bonds, $p_{d}$ and $p_{g}$ are probability distributions determined on data and generated statistics, and $|\cdot|$ is the absolute value.

#### MST:

Denoting the Tanimoto Similarity as TS$=(\mathbf{v}_{1}\cap\mathbf{v}_{2})/(\mathbf{v}_{1}\cup\mathbf{v}_{2})$, where $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ are Morgan Fingerprints based on molecular graphs, we compute the average maximum similarity to test (MST) as:

$$ $\frac{1}{N}\sum_{i=1}^{N}\max_{\mathbf{v}_{t}\in\mathcal{V}_{\text{test}}} \text{TS}(\mathbf{v}_{i},\mathbf{v}_{t}),$ $$

where $\mathcal{V}$ is the set of all test set Morgan Fingerprints.

#### BL and BA:

We determine the Wasserstein W1 distance between the generated and data bond angles and bond lengths.

The angle distributions are defined, for each atom type, by all atoms that have two or more bonds, and we consider all pair-wise angles between position vectors $\mathbf{m}_{i}-\mathbf{m}_{k}$, and $\mathbf{m}_{j}-\mathbf{m}_{k}$, for any two atoms $i$ and $j$ connected to the same central atom $k$. We consider the average W1 distance over all atom types, and we denote it by BA.

Similarly, we create bond Euclidean distance distributions for all bond types $b\in\mathcal{B}$ and denote our averaged W1 distance over all bond types as BL.

We use discretized bins of length $l_{\theta}=0.1$ degrees for angles, and $l_{d}=0.01$Å for bond Euclidean distances, and approximate the generated and data cumulative distribution functions $F_{g}$ and $F_{d}$ and compute:

$$ $\text{BA}=l_{\theta}\frac{1}{N_{\mathcal{A}}}\sum_{a\in\mathcal{A}}\sum_{x\in B _{a}}|F_{d,a}(x)-F_{g,a}(x)|$ $$

$$ $\text{BL}=l_{b}\frac{1}{N_{\mathcal{B}}}\sum_{b\in\mathcal{B}}\sum_{x\in B_{b} }|F_{d,b}(x)-F_{g,b}(x)|,$ $$

where the bin upper limits $B_{a}$ and $B_{b}$ are chosen to contain all maximum and minimum values of the generated and training data samples, $N_{\mathcal{A}}$ and $N_{\mathcal{B}}$ are the number of unique atoms and bonds considered (e.g., we excluded H and F for BA computations, as these atoms never have two bonds).

## Appendix I Explicit aromatic bonds and neutrality

#### Problem explanation

In our settings, we have used the Kekulé molecule representations, which represent aromatic bonds as alternating covalent bond types 1 and 2. In Figure 14, we show three different such aromatic molecules, using the Kekulé representation. Noting that the N and C valencies are 3 and 4, respectively, all formal charges are 0 in these cases (all carbons are assumed to be implicitly connected to a H).

Figure: Figure 14: Three aromatic molecules, with increasingly harder to determine neutrality. In the first case
Refer to caption: x13.png

Instead of the alternating bond types 1 and 2, an equivalent representation denotes these bonds as a separate “aromatic” bond. By considering an electron occupancy of 1.5 for aromatic bonds, Pyrrole has the problem that the N in its molecular graph is apparently (but not in actuality) ionized, having a formal charge of -1. Furthermore, in Pyrrazole, only one of the N atoms should have the (apparent) formal charge -1. Lastly, the C atoms connecting to two rings in Phenanthroline would have a formal charge of -0.5 (as it connects with three 1.5 bonds).

This induces ambiguity. One way to solve this would be to use RDKit and convert the explicit aromatic representation to Kekulé representations. However, in the case of invalid molecules (about 38% of the generated molecules by MiDi), it becomes hard, if not impossible, to determine exactly how many atoms connecting to a (potentially incorrect) aromatic ring have the formal charge 0.

#### Explaining the reported results for MiDi

We have used the molecules reported in the MiDi github repository. We note that on the QM9 dataset, both our method and MiDi generate molecules in Kekulé form, and they are directly comparable in Table 1. For the GEOM neutrality results in Table 3, we have reported an upper-bound for MiDi, by considering bond type “aromatic” as 1.5 electro occupancy and allowing C to have both 0 and -0.5 formal charge. Since many N atoms in GEOM are (in actuality) ionized, it is not ”allowed” to have a formal charge of -1. Regardless, the comparison is very hard to do, and it is difficult to say which model generates more neutral atoms or molecules.

#### MTS bias

We compute the maximum test similarity using valid molecules. This results in a bias in the methods that generate aromatic bonds explicitly, as RDKit will require these molecular graphs to be converted to Kekulé representations. In this case, this metric will not consider molecules containing incorrect aromatic rings (as determined by RDKit’s SanitizeMol). However, in our case, since we directly generate Kekulé molecule representations, this quality filter does not exist anymore, and samples that have incorrect aromatic rings (but are not identified as aromatic because of e.g., one circle atom being incorrectly hybridized) will be considered in the final metric. Note how this bias is unrelated to the performance of the generative model, as those unkekulizable molecules are still generated by explicit-aromatic generating methods, but are simply filtered out by this metric.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: Please see weakness.
soundness: 3
strengths: Chiral awareness in E(3)-invariant architectures is essential for accurately modeling molecular properties. I appreciate the authors’ approach of representing molecules as collections of fields corresponding to atoms and bonds. These fields, as probability densities, closely mirror the quantum mechanics formalism than point clouds, where molecules are wavefunctions. This aligns with methods like Density Functional Theory (DFT), which models charge density and provides a more expressive representation than point clouds. However, the impact on memory overhead from this approach warrants further investigation. Overall, I found this paper to be well-written, clearly delineating previous work, the contributions of this study, and its limitations.
summary: This paper tackles the limitations of E(3)-invariant features in modeling chirality, as mirror symmetry is inherently implied in such features. The authors first provide a theoretical motivation, illustrating why E(3)-invariant features are inadequate for representing chirality. They then propose using fields to model molecular graphs, learning these fields through a reverse diffusion process. The effectiveness of their approach is benchmarked on the QM9 and Geom datasets.
weaknesses: The impact of this approach on memory overhead requires further investigation. How does the memory overhead compare to point cloud representations? Could you provide insights on how memory usage scales with increasing molecule size and greater diversity of atom and bond types?

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - What's the trade-off between losing rotation invariance and gaining the ability of learning chirality?
soundness: 2
strengths: - Since most deep learning methods use GNNs for molecules, using CNNs is novel.
- As stated by the authors, capturing chirality is important for representing molecules. And GNN-based methods may have difficulty in breaking  the reflections symmetry.
summary: This paper proposes to generate molecules by representing them as functions in 3D space. Each atom is represented as a gaussian distribution centered at the atom position. The bonds are represented similarly with the center of the gaussian being placed at the center of the bonds. Compared to  such a representation can learn the chirality since it bears the reflection invariance. Experiments are conducted on QM9 and GEOM.
weaknesses: - Although alignment methods are employed, the proposed method cannot guarantee the rotation invariance.
- As shown in table 2, the molecule graph quality metric is not good on the QM9 dataset.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. Could adaptive grid resolutions or sparse representations help reduce memory requirements while maintaining accuracy?
2. How sensitive is the model's performance to the choice of field parameters (σ, γ)? What's the theoretical basis for choosing these values?
3. Why did the authors choose 0.33Å resolution for the grid, and how would different resolutions affect the quality-efficiency trade-off?
4. How might this approach be extended to other molecular properties beyond chirality? 
5. Since chirality is a crucial property of molecules that are involved in drug discovery, is there any specific case that the author can provide on the FMG? Other than that, some properties that are not related to the chirality of the molecule (i.e. orbital energies), why FMG outperforms?
6. Can this model be utilized in total synthesis and retrosynthesis?
soundness: 3
strengths: The work identifies and addresses a fundamental limitation in current E(3)-equivariant models regarding chirality, which has been largely overlooked in previous research. Besides, it provides original theoretical proofs about the inherent limitations of E(3)-equivariant models in handling chirality.

The overall quality is high, reflected by the rigorous theoretical foundation with formal proofs about E(3) invariance and chirality (Propositions 1-3) and comprehensive empirical validation on standard benchmarks (QM9 and GEOM datasets). Well-structured presentation with clear motivation and problem statement. Ablation studies and
summary: The authors identify a critical limitation in current E(3)-equivariant models for 3D molecular generation: they cannot properly account for molecular chirality, which is crucial for drug efficacy and safety. To address this, they introduce a field-based representation approach that replaces rotational symmetry constraints with reference rotations, allowing their model to correctly capture chiral properties while maintaining competitive performance with existing E(3)-based methods. They provide theoretical proofs showing why E(3)-equivariant models must necessarily ignore chirality, and demonstrate empirically that their method (FMG) achieves state-of-the-art performance in molecular neutrality while successfully capturing chiral properties, though at the cost of increased computational complexity for larger molecules.
weaknesses: The first weakness is the method's computational complexity scales with the volume of the 3D grid (H×W×D), making it potentially intractable for larger molecules. As the author mentioned chirality in amino acids, will be a bottleneck for implementing the model in bio- and pharma-tasks. The current implementation requires high-resolution grids to maintain accuracy, leading to significant memory requirements. Another limitation is the analysis of field parameter sensitivity (σ, γ) is not very clear.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: Questions to the Authors:

1. Your model underperforms compared to JODO, EQGAT-Diff, and Semla-Flow but claims SOTA performance. Can you add a comparison to one of these models?

2. Since non-equivariant models like MCF achieve high performance, can you check whether it still cannot distinguish chiral conformers?

3. Your neutrality metric seems inaccurate (or at least description) (e.g., nitromethane is neutral despite charged atoms). Can you address this issue?

4. Is the strong focus on chirality awareness justified when environmental factors often determine enantiomer generation in practice?
soundness: 2
strengths: The authors proposed a novel approach to molecular generation using field representations and developed a required theoretical background for the approach. The authors prove challenges to incorporate the same restriction in current E(3) equivariant models.
summary: The paper proposes a field-based diffusion method for molecule generation that is not invariant with respect to the chirality of the generated molecules. Recognizing that chirality is a critical property influencing the bioactivity of molecules—with several famous examples underscoring its impact on molecular toxicity and drug efficacy—the approach addresses an important aspect of molecular design.

The results demonstrate comparable performance to EDM and MiDi in terms of molecule validity and stability. The authors also highlight molecular neutrality as an important metric for generated molecules. Furthermore, they assess MMFF energies for the 3D molecules and compute the Wasserstein distance between the bond length distributions and bond angle distributions of the generated molecules and those of the dataset molecules.
weaknesses: 1. Overstated Claims of State-of-the-Art Performance: Although the paper claims to achieve state-of-the-art (SOTA) performance, it falls significantly behind previously published models such as JODO, EQGAT-Diff, and the more recent Semla-Flow—all of which were published before ICLR 2025. The performance gap is quite substantial. For larger and more realistic molecules in the GEOM Drugs dataset, only 16% have all valid atoms according to the paper's model, whereas other models are approaching perfect scores on these metrics.

2. Equivariance Requirement Not Essential for High Performance: Some models that operate directly in E(3) space are eliminating the need for equivariance. For instance, Simplified Scalable Conformer Generation (MCF) demonstrates significant performance improvements in conformer generation without the requirement for equivariance. While MCF tackles a slightly different task, it illustrates that we can develop equivariance-unaware models for molecular modeling that achieve genuine SOTA performance.

3. Incorrect Basis for Neutrality Metrics: The neutrality metrics in the paper are based on inaccurate assertions. Specifically, the statement: "Atom neutrality is the percentage of atoms with a formal charge of 0, while molecules are neutral if all their atoms are neutral." is incorrect. For example, nitromethane has charged oxygen and nitrogen atoms but is an overall neutral molecule. Atoms within a molecule can carry formal charges that balance out, resulting in a neutral molecule despite individual atoms being charged.

4. Inadequate and Incomplete 3D Structure Evaluation: The evaluation of the 3D structures is poor and incomplete. The paper only compares distribution differences for bond lengths, angles, and energies. These benchmarks are difficult to interpret and can be significantly affected by a single outlier result. Depending on the implementation—for example, in MiDi's implementation—one problematic molecule can drastically alter the results. Even with these limitations, the model does not outperform MiDi on these benchmarks.

5. Overstated Importance of Chirality Awareness in the Model: The emphasis on chirality awareness in the model may be overstated. 3D molecule generation serves as a foundational model for structure-based drug design or molecular docking. In these applications, the provided environment (such as a protein or part of a protein) creates asymmetric restrictions on the generated molecule, naturally leading to the production of the appropriate enantiomer. Additionally, in conditional generation, the classifier itself can learn chirality in a simpler way.

Meta-review: This paper proposes a CNN-based method for molecular generation, focusing on effectively capturing chirality which usually fails to be 
 incorporated in E(3) invariant features. The method represents molecules as functions in 3D space. Each atom is represented as a Gaussian distribution centered at the atom position. The method employs fields to generate model molecular graphs, learning these fields through a reverse diffusion process. Given that the chiral awareness in E(3)-invariant architectures is essential for molecular generation and the proposed method is technically novel and empirically strong, I recommend the acceptance of the paper.

Decision: Accept (Poster)