# Field-based Molecule Generation

Alexandru Dumitrescu  Dani Korpela  Markus Heinonen  Yogesh Verma  Valerii
Iakovlev  Vikas Garg  Harri Lähdesmäki

###### 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., 2023b; Guan et
al., 2023). 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, 2009).

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., 2023; Ragoza et
al., 2020) 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 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, 2001).

##  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, 1988), and has been employed in
autoregressive models (Segler et al., 2018) and variational autoencoder (VAE)
based generative models (Gómez-Bombarelli et al., 2018; Kusner et al., 2017).
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., 2018; Jin et al.,
2018), flow-based models (Shi et al., 2020), and reinforcement learning (Shi
et al., 2020; You et al., 2018; Zhou et al., 2019) 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., 2018, 2019) have been
developed, with the notable downside that they introduce an artificial atom
order. GNN parametrized denoising diffusion models (Guan et al., 2023; Peng et
al., 2023b; Vignac et al., 2023; Xu et al., 2023; Wu et al., 2022) 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., 2020) and score-based models (Pinheiro et al., 2023), 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 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.

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.

###  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., 2015; Ho et al., 2020b; Kingma et al., 2021), 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. (2021)). 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., 2021),
but earlier works have argued for a simplified loss (Ho et al., 2020b)

$\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, 2021). 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. (2020a). To further emphasize the stability of generated molecules,
Vignac et al. (2023) and Peng et al. (2023b) 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. (2023), and vice-versa in Peng et al. (2023b).

We achieved better results using the approach of Vignac et al. (2023), 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, 2021). 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 (2022) 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$.

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 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.

##  4 Evaluation

#### Basic molecule properties

We use atom and molecule neutrality, validity, and novelty, as defined in
(Guan et al., 2023). We call the stability metric introduced in (Guan et al.,
2023) neutrality, to clarify and distinguish from other work that considered
stable atoms even those without a formal charge of 0 (Vignac et al., 2023).
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., 2023).

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., 2023) and (Peng et al., 2023a) 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, 2023) 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., 2023), 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., 2014) 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., 2023), 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 W1 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., 2023), 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, 2022) 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., 2023; Peng et al., 2023a), we select the 30 lowest
energy conformers for each molecule graph. Similar to (Peng et al., 2023a), 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 TVa 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.,
2023), 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 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.

###  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\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., 2023) 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.

## References

  * Axelrod & Gómez-Bombarelli (2022) Axelrod, S. and Gómez-Bombarelli, R.  Geom, energy-annotated molecular conformations for property prediction and molecular generation.  _Scientific Data_ , 9(1):185, 2022.
  * Dhariwal & Nichol (2021) Dhariwal, P. and Nichol, A.  Diffusion models beat gans on image synthesis.  In _NeurIPS_ , 2021.
  * Evans (2001) Evans, A.  Comparative pharmacology of s(+)-ibuprofen and (rs)-ibuprofen.  _Clin Rheumatol_ , 20:9–14, 2001.
  * Gebauer et al. (2018) Gebauer, N., Gastegger, M., and Schütt, K.  Generating equilibrium molecules with deep neural networks.  _arXiv_ , 2018.
  * Gebauer et al. (2019) Gebauer, N., Gastegger, M., and Schütt, K.  Symmetry-adapted generation of 3d point sets for the targeted discovery of molecules.  In _NeurIPS_ , 2019.
  * Guan et al. (2023) Guan, J., Qian, W. W., Peng, X., Su, Y., Peng, J., and Ma, J.  3d equivariant diffusion for target-aware molecule generation and affinity prediction.  In _ICLR_ , 2023.
  * Gómez-Bombarelli et al. (2018) Gómez-Bombarelli, R., Wei, J., Duvenaud, D., Hernández-Lobato, J. M., Sánchez-Lengeling, B., Sheberla, D., Aguilera-Iparraguirre, J., Hirzel, T., Adams, R., and Aspuru-Guzik, A.  Automatic chemical design using a data-driven continuous representation of molecules.  _ACS Central Science_ , 4(2):268–276, 2018.
  * Ho & Salimans (2022) Ho, J. and Salimans, T.  Classifier-free diffusion guidance.  In _NeurIPS workshop DGM_ , 2022.
  * Ho et al. (2020a) Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  In _NeurIPS_ , 2020a.
  * Ho et al. (2020b) Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  In _NeurIPS_ , 2020b.
  * Jin et al. (2018) Jin, W., Barzilay, R., and Jaakkola, T.  Junction tree variational autoencoder for molecular graph generation.  In _ICML_ , 2018.
  * Kingma et al. (2021) Kingma, D. P., Salimans, T., Poole, B., and Ho, J.  Variational diffusion models.  In _NeurIPS_ , 2021.
  * Kusner et al. (2017) Kusner, M. J., Paige, B., and Hernández-Lobato, J. M.  Grammar variational autoencoder.  In _ICML_ , 2017.
  * Liu et al. (2018) Liu, Q., Allamanis, M., Brockschmidt, M., and Gaunt, A.  Constrained graph variational autoencoders for molecule design.  In _NeurIPS_ , 2018.
  * Nichol & Dhariwal (2021) Nichol, A. and Dhariwal, P.  Improved denoising diffusion probabilistic models.  In _ICML_ , 2021.
  * Peng et al. (2023a) Peng, X., Guan, J., Liu, Q., and Ma, J.  Moldiff: Addressing the atom-bond inconsistency problem in 3d molecule diffusion generation.  In _ICML_ , 2023a.
  * Peng et al. (2023b) Peng, X., Guan, J., Liu, Q., and Ma, J.  MolDiff: Addressing the atom-bond inconsistency problem in 3D molecule diffusion generation.  In _ICML_ , 2023b.
  * Perez et al. (2018) Perez, E., Strub, F., de Vries, H., Dumoulin, V., and Courville, A. C.  Film: Visual reasoning with a general conditioning layer.  In _AAAI_ , 2018.
  * Pinheiro et al. (2023) Pinheiro, P., Rackers, J., Kleinhenz, J., Maser, M., Mahmood, O., Watkins, A. M., Ra, S., Sresht, V., and Saremi, S.  3d molecule generation by denoising voxel grids.  In _NeurIPS_ , 2023.
  * Ragoza et al. (2020) Ragoza, M., Masuda, T., and Koes, D. R.  Learning a continuous representation of 3d molecular structures with deep generative models, 2020.
  * Ramakrishnan et al. (2014) Ramakrishnan, R., Dral, P. O., Rupp, M., and von Lilienfeld, O. A.  Quantum chemistry structures and properties of 134 kilo molecules.  _Scientific Data_ , 1(1):140022, 2014.
  * Segler et al. (2018) Segler, M., Kogej, T., Tyrchan, C., and Waller, M.  Generating focused molecule libraries for drug discovery with recurrent neural networks.  _ACS Central Science_ , 2018.
  * Shi et al. (2020) Shi, C., Xu, M., Zhu, Z., Zhang, W., Zhang, M., and Tang, J.  Graphaf: a flow-based autoregressive model for molecular graph generation.  In _ICLR_ , 2020.
  * Smith (2009) Smith, S. W.  Chiral Toxicology: It’s the Same Thing…Only Different.  _Toxicological Sciences_ , 110(1):4–30, 2009.
  * Sohl-Dickstein et al. (2015) Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N., and Ganguli, S.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _ICML_ , 2015.
  * Vignac et al. (2023) Vignac, C., Osman, N., Toni, L., and Frossard, P.  Midi: Mixed graph and 3d denoising diffusion for molecule generation.  In _ICLR workshop MLDD_ , 2023.
  * Weininger (1988) Weininger, D.  Smiles, a chemical language and information system. 1. introduction to methodology and encoding rules.  _Journal of Chemical Information and Computer Sciences_ , 28(1):31–36, 1988.
  * Wu et al. (2022) Wu, L., Gong, C., Liu, X., Ye, M., and qiang liu.  Diffusion-based molecule generation with informative prior bridges.  In _NeurIPS_ , 2022.
  * Xu et al. (2023) Xu, M., Powers, A. S., Dror, R. O., Ermon, S., and Leskovec, J.  Geometric latent diffusion models for 3D molecule generation.  In _ICML_ , 2023.
  * You et al. (2018) You, J., Liu, B., Ying, Z., Pande, V., and Leskovec, J.  Graph convolutional policy network for goal-directed molecular graph generation.  In _NeurIPS_ , 2018.
  * Zhong & Guan (2023) Zhong, S. and Guan, X.  Count-based morgan fingerprint: A more efficient and interpretable molecular representation in developing machine learning-based predictive regression models for water contaminants’ activities and properties.  _Environmental Science & Technology_, 57(46):18193–18202, 2023.
  * Zhou et al. (2019) Zhou, Z., Kearnes, S., Li, L., Zare, R. N., and Riley, P.  Optimization of molecules via deep reinforcement learning.  _Scientific Reports_ , 2019.

##  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 W1
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 (↓W1) |  BL (↓W${}_{1}\cdot 10^{2}$ ) |  BA (↓W1)
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 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.

###  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., 2023) 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β=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β=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., 2018), 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 7: Molecular graphs of half of the chiral dataset.

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 8: R and S configurations of a general tetrahedral configuration with
three, distinct residues.

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.,
2020b). 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.

Algorithm 1 Peak extraction

Input: fields $\mathbf{u}\in\mathbb{R}^{N_{x}\times N_{y}\times N_{z}}$,
threshold $t$, $q\leftarrow\\{\mathbf{x};\mathbf{u}(\mathbf{x})>t\\}$;
$A_{\mathbf{m}}\leftarrow\\{\\}$

repeat

$p=q.pop()$

$neigh\leftarrow\text{get\\_neigh}(p,t,\mathbf{u})$

$A_{\mathbf{m}}.\text{insert}(\text{mean}(neigh,p))$

$q.pop(neigh)$

until q is empty

return $A_{\mathbf{m}}$

#### 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 $\gamma_{ij}$

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 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.

#### 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 10: Random valid samples from FMG trained on QM9.

Figure 11: Random valid samples from FMG trained on GEOM.

##  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 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.

##  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 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.

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.

#### TVa, TVb

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 14: Three aromatic molecules, with increasingly harder to determine
neutrality. In the first case

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.