Year: 2021

Abstract:
Automating molecular design using deep reinforcement learning (RL) has the potential to greatly accelerate the search for novel materials. Despite recent progress on leveraging graph representations to design molecules, such methods are fundamentally limited by the lack of three-dimensional (3D) information. In light of this, we propose a novel actor-critic architecture for 3D molecular design that can generate molecular structures unattainable with previous approaches. This is achieved by exploiting the symmetries of the design process through a rotationally covariant state-action representation based on a spherical harmonics series expansion. We demonstrate the benefits of our approach on several 3D molecular design tasks, where we find that building in such symmetries significantly improves generalization and the quality of generated molecules.

text: Automating molecular design using deep reinforcement learning (RL) has the potential to greatly accelerate the search for novel materials. Despite recent progress on leveraging graph representations to design molecules, such methods are fundamentally limited by the lack of three-dimensional (3D) information. In light of this, we propose a novel actor-critic architecture for 3D molecular design that can generate molecular structures unattainable with previous approaches. This is achieved by exploiting the symmetries of the design process through a rotationally covariant state-action representation based on a spherical harmonics series expansion. We demonstrate the benefits of our approach on several 3D molecular design tasks, where we find that building in such symmetries significantly improves generalization and the quality of generated molecules.

heading: 1 INTRODUCTION
text: The search for molecular structures with desirable properties is a challenging task with important applications in de novo drug design and materials discovery (Schneider et al., 2019). There exist a plethora of machine learning approaches to accelerate this search, including generative models based on variational autoencoders (VAEs) (Gómez-Bombarelli et al., 2018), recurrent neural networks (RNNs) (Segler et al., 2018), and generative adversarial networks (GANs) (De Cao & Kipf, 2018). However, the reliance on a sufficiently large dataset for exploring unknown regions of chemical space is a severe limitation of such supervised models. Recent RL-based methods (e.g., Olivecrona et al. (2017), Jørgensen et al. (2019), Simm et al. (2020)) mitigate the need for an existing dataset of molecules as they only require access to a reward function.
Most approaches rely on graph representations of molecules, where atoms and bonds are represented by nodes and edges, respectively. This is a strongly simplified model designed for the description of single organic molecules. It is unsuitable for encoding metals and molecular clusters as it lacks information about the relative position of atoms in 3D space. Further, geometric constraints on the design process cannot be included, e.g. those given by the active site of an enzyme. A more general representation closer to the physical system is one in which a molecule is described by its atoms’ positions in Cartesian coordinates. However, it would be very inefficient to naively learn a model based on this representation. That is because molecular properties such as the energy are invariant (i.e. unchanged) under symmetry operations like translation or rotation of all atomic positions. A model without the right inductive bias would thus have to learn those symmetries from scratch.
In this work, we develop a novel RL approach for designing molecules in Cartesian coordinates that explicitly encodes these symmetry operations. The agent builds molecules by consecutively placing atoms such that if the generated structure is rotated or translated, the agent’s action is rotated and translated accordingly; this way, the reward remains the same (see Fig. 1 (a)). We achieve this through a rotationally covariant state representation based on spherical harmonics, which we integrate into a novel actor-critic network architecture with an auto-regressive policy that maintains the desired covariance. Building in this inductive bias enables us to generate molecular structures with more complex coordination geometry than the class of molecules that were attainable with previous approaches. Finally, we perform experiments on several 3D molecular design tasks, where
we find that our approach significantly improves the generalization capabilities of the RL agent and the quality of the generated molecules.
In summary, our contributions are as follows:
• we propose the first approach for 3D molecular design that exploits symmetries of the design process by leveraging a rotationally covariant state representation;
• we integrate this state representation into an actor-critic neural network architecture with a rotationally covariant auto-regressive policy, where the orientation of the atoms to be placed is modeled through a flexible distribution based on spherical harmonics;
• we demonstrate the benefits of our approach on several 3D molecular design tasks, including a newly proposed task that showcases the generalization capabilities of our agent.

heading: 2 BACKGROUND

heading: 2.1 REINFORCEMENT LEARNING FOR MOLECULAR DESIGN
text: In the standard RL setting (Sutton & Barto, 2018), an agent interacts with the environment to maximize its reward. Formally, such an environment is described by a Markov decision process (MDP) M = (S,A, T , µ0, γ, r) with states st ∈ S , actions at ∈ A, transition dynamics T : S × A 7→ S, initial state distribution µ0, discount factor γ ∈ (0, 1], and reward function r : S × A 7→ R. The goal is to learn a stochastic policy π(at|st) that maximizes the expected discounted return J(θ) = Es0∼µ0 [V π(s0)], where the value function V π(st) = Eπ[ ∑T t′=t γ
t′r(st′ , at′)|st] is defined as the expected discounted return when starting from state st and following policy π.
Following Simm et al. (2020), we design molecules by iteratively picking atoms from a bag and positioning them on a 3D canvas. Such a sequential decision-making problem is described by an MDP where the state st = (Ct,Bt) comprises both the canvas Ct and the bag Bt. The canvas Ct = C0∪{(ei, xi)}t−1i=0 is a set of atoms with chemical element ei ∈ {H,C,N,O, . . . } and position xi ∈ R3 placed up to time t − 1, where C0 can either be empty or contain a set of initially placed atoms. The number of atoms on the canvas is denoted by |Ct|. The bag Bt = {(e,m(e))} is a multi-set of atoms yet to be placed, where m(e) is the multiplicity of the element e. Each action at = (et, xt) consists of the element et ∈ Bt and position xt ∈ R3 of the next atom to be added to the canvas. Placing an atom through action at in state st is modeled by a deterministic transition function T (st, at) that yields the next state st+1 = (Ct+1,Bt+1) with Bt+1 = Bt\et. The reward function r(st, at) = −∆E(st, at) is given by the negative energy difference between the resulting structure described by Ct+1, and the sum of energies of the current structure Ct and a new atom of element et placed at the origin, i.e. ∆E(st, at) = E(Ct+1) − [E(Ct) + E({(e,0)})]. Intuitively, the reward encourages the agent to build stable, low-energy structures. We evaluate the energy using the fast semi-empirical Parametrized Method 6 (PM6) (Stewart, 2007) as implemented in SPARROW (Husch et al., 2018; Bosia et al., 2020); see Appendix A for details.
An example of a rollout is shown in Fig. 1 (b). At the beginning of the episode, the agent observes the initial state (C0,B0) ∼ µ0(s0), e.g. C0 = ∅ and B0 = SOF41. The agent then iteratively constructs a molecule by placing atoms from the bag onto the canvas until the bag is empty.2

heading: 2.2 ROTATIONALLY COVARIANT NEURAL NETWORKS
text: A function f : X 7→ Y is invariant under a transformation operator Tg : X 7→ X if f(Tg[x]) = f(x) for all x ∈ X , g ∈ G, where G is a mathematical group. In contrast, f is covariant with respect to Tg if there exists an operator T ′g : Y 7→ Y such that f(Tg[x]) = T ′g[f(x)]. To achieve rotational covariance, it is natural to work with spherical harmonics. They are a set of complexvalued functions Y m` : S2 7→ C with ` = 0, 1, 2, . . . and m = −`,−` + 1, . . . , ` − 1, ` on the unit sphere S2 in R3. The first few spherical harmonics are in Appendix B. They are defined by
Y m` (ϑ, ϕ) = (−1)m √ 2`+ 1
4π (`−m)! (`+m)! Pml (cos(ϑ))e imϕ, ϕ ∈ [0, 2π], ϑ ∈ [0, π], (1)
where Pm` denotes the associated normalized Legendre polynomials of the first kind (Bateman, 1953), and each Y m` is normalized such that ∫∫ |Y m` (ϑ, ϕ)|2 sinϑdϑdϕ = 1. Any square-integrable function f : S2 7→ C can be written as a series expansion in terms of the spherical harmonics,
f(x̃) = ∞∑ `=0 ∑̀ m=−` f̂m` Y m ` (x̃), (2)
where x̃ = (ϑ, ϕ) ∈ S2. The complex-valued coefficients {f̂m` } are the analogs of Fourier coefficients and are given by f̂m` = ∫ f(x̃)Y m∗` (x̃)Ω(dx̃). Such a function f can be modeled by learning the coefficients {f̂m` } using CORMORANT (Anderson et al., 2019), a neural network architecture for predicting properties of chemical systems that works entirely in Fourier space. A key feature is that each neuron is covariant to rotation but invariant to translation; further, each neuron explicitly corresponds to a subset of atoms in the molecule. The input of CORMORANT is a spherical function f0 : S2 7→ Cd and the output is a collection of vectors f̂ = {f̂0, f̂1, . . . , f̂L}, where each f̂` ∈ τ × (2`+ 1) is a rotationally covariant vector with τ channels. That is, if the input is rotated by R ∈ SO(3), then each f̂` transforms as f̂` 7→ D`(R)f̂`, where D`(R) : SO(3) 7→ C(2`+1)×(2`+1) are the irreducible representations of SO(3), also called the Wigner D-matrices.

heading: 3 COVARIANT POLICY FOR MOLECULAR DESIGN
text: An efficient RL agent needs to exploit the symmetries of the molecular design process. Therefore, we require a policy π(a|s) with actions a = (e, x) that is covariant under translation and rotation with respect to the position x, i.e., x should rotate (or translate) accordingly if the atoms on the canvas C are rotated (or translated). In contrast, the policy needs to be invariant to the element e, i.e. the chosen element remains unchanged under such transformations (see Fig. 1 (a)). Since learning such a policy is difficult when working directly in global Cartesian coordinates, we instead follow Simm et al. (2020) and use an action representation that is local with respect to an already placed focal atom. If the next atom is placed relative to the focal atom, covariance under translation of x is automatically achieved and only the rotational covariance remains to be dealt with.
As shown in Fig. 2, we model the action a through a sequence of sub-actions: (1) the index f ∈ {1, . . . , |C|} of the focal atom around which the next atom is placed,3 (2) the element e ∈ {1, . . . , Ne} of the next atom from the set of available elements, (3) a distance d ∈ R+ between the focal atom and the next atom, and (4) the orientation x̃ = (ϑ, ϕ) ∈ S2 of the atom on a unit sphere around the focal atom. Denoting xf as the position of the focal atom, we obtain action a = (e, x) by mapping the local coordinates (x̃, d, f) to global coordinates x = xf + d · x̃, where
1Shorthand for {(S, 1), (O, 1), (F, 4)}. 2Hereafter, we drop the time index when it is clear from the context. 3If the canvas C0 is empty, the agent selects an element e0 ∈ B0 and places it at the origin, i.e. a0 = (e0,0).
x is now covariant under translation and rotation. We choose these sub-actions using the following auto-regressive policy:
π(a|s) = π(x̃, d, e, f |s) = p(x̃|d, e, f, s) p(d|e, f, s) p(e|f, s) p(f |s). (3)
A novel actor-critic neural network architecture that implements this policy is illustrated in Fig. 3. In the following, we discuss its state embedding, actor, and critic networks in more detail.

heading: 3.1 STATE EMBEDDING
text: The state embedding network transforms canvas C and bag B to obtain a rotationally covariant and translationally invariant representation. For that, we concatenate a vectorized representation of the bag with each atom on the canvas and feed it into CORMORANT, i.e. scov ← CORMORANT(C,B), where scov = {scov` } Lmax `=0 , s cov ` ∈ C|C|×τ×(2`+1), and τ is the number of channels. For the sake of exposition, we assume a single channel for each element in the bag, i.e. τ = Ne (cf. Fig 3); in practice, we use up to four channels per element.
However, not every sub-action in Eq. (3) should transform equally under rotation and translation. While the orientation x̃ needs to be covariant under rotation, the choice of focal atom f , element e, and distance d have to be invariant to rotation and translation. For these sub-actions, we additionally require an invariant state representation. To obtain such a representation sinv ∈ R|C|×k, we employ a combination of transformations from Anderson et al. (2019) as listed in Appendix C (e.g. for ` = 0, one can simply select the scov`=0 component), which we collectively denote as Tinv.

heading: 3.2 ACTOR
text: Focal Atom and Element The distribution p(f |s) over the focal atom f is modeled as categorical, f ∼ Cat(f ;hf ), where hf are the logits for each atom in C predicted by a multi-layer perceptron (MLP). Likewise, the distribution over the element e is given by p(e|f, s) = Cat(e;he) with he = MLPe(sinvf ), where s inv f is the invariant representation for the focal atom. Since the number of possible focal atoms f increases and the set of available elements e decreases during a rollout, we mask out invalid focal atoms f /∈ {1, . . . , |Ct|} and elements e /∈ Bt by setting their probabilities to zero and re-normalizing the categorical distributions. The agent does not make use of chemical concepts like bond connectivity to aid the choice of the focal atom.
Distance We select the channel τe corresponding to element e from scovf to obtain s cov f,e := {scov`f,e}`=0,...,Lmax and s inv f,e ← Tinv(scovf,e). Then we model the distribution over the distance d between the focal atom and the next atom to be placed as a mixture of M Gaussians, p(d|e, f, s) =∑M m=1 πm N (µm, σ2m), where πm is the mixing coefficient of the m-th GaussianN (µm, σ2m). The mixing coefficients and the means are predicted by a mixture density network (MDN) (Bishop, 1994), i.e. {πm, µm}Mm=1 = MDN(sinvf,e). The standard deviations {σm}Mm=1 are global parameters. We guarantee that the sampled distances are positive by clipping values below zero.
Combining Invariant and Covariant Features The choice of distance d can significantly affect the orientation x̃ of the atom. For example, if d has the length of a triple bond, then x̃ will be very different from if it was a single bond. Thus, we condition scovf,e on distance d through a non-linear and learnable transformation that preserves rotational covariance. We then use this representation
to model a spherical distribution over x̃. Kondor & Trivedi (2018) showed that a linear transformation with learnable parameters is only covariant if the operation combines fragments with the same `. Further, the Clebsch-Gordan (CG) non-linearity allows one to combine two covariant features such that the result is still covariant. Thus, we obtain a rotationally covariant representation r̂ := {r̂`}`=0,...,Lmax ← Tcov(d, scovf,e) conditioned on all previous sub-actions as follows:
r̂` = [ scov`f,e ⊕ d · s cov `f,e ⊕ (d · scovf,e ⊗cg d · scovf,e)` ] ·W` ∀`, (4)
where ⊕ denotes the appropriate concatenation of matrices, and W` is a learnable complex-valued matrix. As in Anderson et al. (2019), we perform the CG product ⊗cg only channel-wise to reduce computational complexity.
Orientation Next, we utilize r̂ to obtain a rotationally covariant spherical distribution for the orientation x̃ based on the series expansion in Eq. (2). Taking inspiration from commonly used distributions (Jammalamadaka & Terdik, 2019), we propose to use the following expression:
p(x̃|d, e, f, s) = 1 Z exp −β ∣∣∣∣∣ 1√k Lmax∑ `=0 ∑̀ m=−` r̂m` Y m ` (x̃) ∣∣∣∣∣ 2  , (5)
where β ∈ R is a scaling parameter, and the term 1/ √ k with k = ∑Lmax `=0 ∑` m=−`|r̂m` |2 regularizes the distribution so that it does not approach a delta function. The normalization constant Z is estimated via Lebedev quadrature (Lebedev, 1975; 1977). We sample from the distribution in Eq. (5) using rejection sampling (Bishop, 2009) with a uniform proposal distribution q(x̃) = (4π)−1. Note that in contrast to more commonly used parametric distributions (e.g. von Mises-Fisher), this formulation allows to model multi-modalities. We discuss alternatives to Eq. (5) in Appendix D.

heading: 3.3 CRITIC
text: The critic needs to compute a value V for the state s that is invariant under translation and rotation. Given sinv, we apply a permutation-invariant set encoding (Zaheer et al., 2017) of the atoms, i.e.
V (s) = MLPρ  |C|∑ i=1 MLPφ(sinvi )  . (6) Finally, we use PPO (Schulman et al., 2017) to learn the parameters of the actor-critic architecture. To encourage sufficient exploration, we add an entropy regularization term over the categorical subaction distributions of the auto-regressive policy. For offline evaluation, we evaluate the policy without any exploration by choosing the most probable action. While the mode of the distributions for f and e is available in closed form, we approximate the global mode of the distributions over d and x̃ by evaluating the density at S samples and picking the one with the highest density.

heading: 4 RELATED WORK
text: Reinforcement Learning for Molecular Design There exists a large variety of RL-based approaches for molecular design using either string- or graph-based representations of molecules (Olivecrona et al., 2017; Guimaraes et al., 2018; Putin et al., 2018; Neil et al., 2018; Popova et al., 2018; You et al., 2018; Zhou et al., 2019). However, the choice of representation limits the molecules that can be generated to a (small) region of chemical space for which the representation is applicable, i.e., single organic molecules. Such representations also prohibit the use of reward functions based on quantum-mechanical properties; instead, heuristics are often used. Lastly, geometric constraints on the design process cannot be imposed as the representation does not include any 3D information.
Molecular Design in Cartesian Coordinates Another downside of string- and graph-based approaches is their neglect of information encoded in the interatomic distances. In light of this, Gebauer et al. (2018; 2019) proposed a supervised generative neural network for sequentially placing atoms in Cartesian coordinates. While the model respects local symmetries by construction, atoms are placed on a 3D grid. Similar to other supervised approaches, one further requires a dataset that covers the particular class of molecules to be generated. Hammer and coworkers (Jørgensen et al., 2019; Meldgaard et al., 2020) employed a Deep Q-Network (Mnih et al., 2015) to build planar compounds and crystalline surfaces by placing atoms on a grid. Recently, Simm et al. (2020) presented an RL formulation for molecular design in continuous 3D space. The agent models the position of the next atom to be placed in internal coordinates— i.e. the distance, angle, and dihedral angle with respect to already existing atoms—which are invariant under translation and rotation. By mapping from internal to Cartesian coordinates, they then obtain a policy that is covariant under these symmetry operations. However, as shown in Fig 4, the angle and dihedral angle are only defined with respect to two reference points, which are
chosen to be the two closest points to a focal atom. In highly symmetric states, e.g. as commonly encountered in materials, this representation fails to distinguish different configurations as one cannot uniquely select the two closest atoms as reference points anymore. In contrast, we do not rely on such reference points as the agent directly samples the orientation from a spherical distribution.
Covariant Neural Networks in Chemical Science Prior work employed rotationally covariant neural networks to predict translation- and rotation-invariant physical properties (Thomas et al., 2018; Kondor et al., 2018; Weiler et al., 2018; Anderson et al., 2019; Miller et al., 2020; Finzi et al., 2020; Fuchs et al., 2020), e.g. scalars such as the electronic energy. In contrast, we propose a translation-invariant and rotation-covariant neural network architecture for generating molecules. For a more general treatment of covariance (or equivariance) in RL, see van der Pol et al. (2020).

heading: 5 EXPERIMENTS
text: We perform experiments to answer the following questions: (1) is the agent able to learn how to build highly symmetric molecules in Cartesian coordinates from scratch, (2) can we increase the validity, diversity, and stability of generated molecules, and (3) does our approach lead to improved generalization? We address (1) and (2) by evaluating the agent on a diverse range of tasks from the MOLGYM benchmark suite (Simm et al., 2020), and (3) on a newly proposed stochastic-bag task (see Section 5.1) where bags are sampled from a distribution over bags. In Appendix G, we show with an additional experiment that the agent can learn to place water molecules around a given solute to form a solvation shell.4
We compare our approach (COVARIANT) against the RL agent proposed by Simm et al. (2020), which iteratively builds molecules on a 3D canvas by working in internal coordinates (INTERNAL). As an additional baseline, we consider a classical, optimization-based agent (OPT) with access to a black-box function that yields the energy E(C) and the atomic forces F (C) for a given canvas.5 The agent constructs molecules by alternating between randomly placing an atom and optimizing the structure. Moreover, the agent applies several heuristics inspired by fundamental chemical concepts to guide the placement of atoms. To make the comparisons fair, we grant OPT a comparable computational budget in terms of the total number of energy computations. Finally, for some experiments, the best possible performance based on quantum-chemical calculations can be reported. See Appendices E and F for more details on the baselines and experimental settings, and Appendix H for an additional runtime comparison between the agents.

heading: 5.1 STOCHASTIC-BAG TASK
text: In Simm et al. (2020), a set of molecular design tasks was introduced: the single-bag task assesses an agent’s ability to build single stable molecules, whereas the multi-bag task focuses on building several molecules of different composition and size at the same time. A limitation of these tasks is that the initial bags were selected such that they correspond to known formulas, which in practice might not be known a priori. In the stochastic-bag task, we relax this assumption by sampling from a more general distribution over bags. Before each episode, we construct a bag B = {(e,m(e))} by sampling counts (m(e1), ...,m(emax)) ∼ Mult(ζ, pe), where the bag size ζ is sampled uniformly from the interval [ζmin, ζmax]. Here, we obtain an empirical distribution pe from the multiplicities m(e) of a given bag B∗. For example, with B∗ = {(H, 2), (O, 1)} we obtain pH = 23 and pO = 1 3 . Since sampled bags might no longer correspond to valid molecules when placed completely, we discard bags where the sum of valence electrons over all atoms contained in the bag is odd. This ensures that the agent can build a closed-shell system.

heading: 5.2 RESULTS
text: Building Highly Symmetric Molecules First, we evaluate the ability to build stable molecules featuring high symmetry and coordination numbers (e.g. trigonal bipyramidal, square pyramidal, and octahedral) on the single-bag task with bags SOF4, IF5, SOF6, and SF6. As shown in Fig. 5 (a), COVARIANT can solve the task for SOF4 and IF5 within 30 000 to 40 000 steps, whereas INTERNAL fails to build low-energy configurations as it cannot distinguish highly symmetric intermediates (cf. Fig. 4). Further results in Fig. 5 (b) for SOF6 and SF6 show that COVARIANT is capable of building such structures. Likewise, OPT found the optimal structures for all four bags. While the constructed molecules are small in size, they would be unattainable with graph- or string-based methods as such representations lack important 3D information. For example, RDKIT (Landrum, 2019), a state-of-the-art library for 3D structure generation of organic molecules, failed at this task.
Validity, Diversity, and Stability of Generated Molecules Since string and graph representations are not well-suited for designing molecules with complex 3D structure, it is difficult to directly compare to most prior work. To still enable comparisons, we follow the GuacaMol benchmark (Brown et al., 2019) and report the chemical validity, diversity, and stability of the molecules generated by the agents for different experiments. A generated structure is considered valid if it can be successfully converted to a molecular graph by the tool XYZ2MOL (Jensen, 2019; Kim & Kim,
4Source code of the agent and environment is available at https://github.com/gncs/molgym. 5For the calculation of E(C) and F (C) we employ PM6; the same method as in the reward function.
2015). The validity reported in Table 1 is the ratio of valid molecules generated during offline evaluation at the end of training over 10 seeds. Two molecules are considered identical if their molecular graphs yield the same SMILES strings under RDKIT. The diversity shown in Table 1 is the total number of unique and valid structures generated during offline evaluation during training over 10 seeds.6 In the two stochastic-bag experiments, the agents are trained on bags of sizes from the interval [16, 22] sampled with B∗ = C7H8N2O2 and C7H10O2, respectively. Finally, to assess the stability of the generated molecules, valid structures generated in the last iteration underwent a structure optimization using the PM6 method (see Appendix A for details). Then, the root-mean-square deviation of atomic positions (RMSD, in Å) between the original and the optimized structure was computed. In Table 1, the median RMSD is given per experiment.
Results are listed in Table 1. We observe that COVARIANT significantly outperforms the other agents on most experiments both in terms of validity and diversity. The difference is particularly large for the more challenging stochastic-bag tasks, where COVARIANT does similarly well as on the singlebag experiments. This finding is confirmed in Fig. 6 (a), showing that the exact stoichiometry does not need to be known a priori for the agent to build valid molecules. Moreover, the structures generated by COVARIANT are overall slightly more stable compared to INTERNAL. In contrast, OPT often fails to build valid structures. Inspection of the generated structures reveals that for larger bags the agent tends to build multi-molecular clusters, which are considered invalid in this experiment. The stability for OPT is omitted as all of its valid structures are stable by definition.
Compared to graph-based approaches (e.g., Jin et al. (2017); Bradshaw et al. (2019a); Li et al. (2018b;a); Liu et al. (2018); De Cao & Kipf (2018); Bradshaw et al. (2019b)), the average validity
6For a fairer comparison in the stochastic-bag task, we use the structures generated during offline evaluation, instead of those generated during training as in Simm et al. (2020).
and diversity achieved by COVARIANT are still relatively low. This can partly be explained by the fact that state-of-the-art graph-based approaches have the strict rules of chemical bonding in organic molecules encoded into their models. But as a result, they are limited to generating single organic molecules and cannot build molecules for which these rules do not apply (e.g., hypervalent iodine compounds such as IF5). In terms of stability, the supervised generative model by Gebauer et al. (2019) reported an average RMSD of approximately 0.25 Å. While their approach and the considered molecules are significantly different from ours, this suggests that the generated structures are more stable compared to COVARIANT. Nonetheless, the RL approach presented in this work remains particularly attractive if no dataset exists on which such a supervised model can be trained.
Generalization To evaluate the generalization capabilities of all agents to unseen bags, we train on a distribution over bags with B∗ = C7H10O2 and C7H8N2O2, and test on sets of larger, outof-distribution bags {C6H14O3,C7H16O,C7H16O2,C8H18O}, and {C8H12N2O,C6H12N2O3, C7H14N2O,C7H14N2O2} respectively. As shown in Fig. 6 (b), COVARIANT obtains higher average returns with lower variance compared to INTERNAL on C7H10O2, while performing only slightly worse than the agent trained directly on the test bags (purple). Results for C7H8N2O2 are in Appendix G. Although the difference in performance seems to be marginal, we stress that chemical validity is often determined by the last 10% of the returns. Indeed, Table 1 and Fig. 9 in Appendix G highlight the higher quality of the structures generated by COVARIANT, indicating better generalization to unseen bags of larger size. OPT fails at this task.

heading: 6 CONCLUSION
text: We proposed a novel covariant actor-critic architecture based on spherical harmonics for designing highly symmetric molecules in 3D. We showed empirically that exploiting symmetries of the molecular design process improves the quality of the generated molecules and leads to better generalization. In future work, we aim to employ more accurate quantum-chemical methods (e.g., density functional theory) required for building transition metal complexes or structures in which weak intermolecular interactions are important. For that, however, the sample-efficiency of our agent needs to be improved. Finally, we aim to explore reward functions specifically tailored towards drug design.

heading: ACKNOWLEDGEMENTS
text: We thank A. J. Tripp and K. T. Jensen for useful discussions and feedback. RP receives funding from iCASE grant #1950384 with support from Nokia. JMHL acknowledges support from a Turing AI Fellowship under grant EP/V023756/1. This work has been performed using resources operated by the University of Cambridge Research Computing Service (funded by grant EP/P020259/1).

heading: A REWARD CALCULATION
text: In the reward function, the energyE has to be computed using quantum-chemical methods. For that, we use the fast semi-empirical Parametrized Method 6 (PM6) (Stewart, 2007). In particular, we use the implementation in the software package SPARROW (Husch et al., 2018; Bosia et al., 2020). For each calculation, a molecular charge of zero and the lowest possible spin multiplicity are chosen. All calculations are spin-unrestricted.
Limitations of semi-empirical methods are highlighted in, for example, recent work by Husch & Reiher (2018). More accurate methods such as approximate density functionals need to be employed especially for systems containing transition metals.
Further, we enforce that atoms are not placed too close (< 0.6 Å) nor too far away from each other (> 2.0 Å). If the agent places an atom outside these boundaries, the minimum reward of −0.6 is awarded and the episode terminates. Further, the environment encourages the agents to build single molecular structures by terminating the episode and return a reward of −0.6 if elements forming stable bimolecular compounds (e.g, H2) are placed too far away from other atoms on the canvas.

heading: B SPHERICAL HARMONICS
text: The spherical harmonics form an orthonormal basis of the Hilbert space of square-integrable functions L2C(S2). The first few spherical harmonics are given by:
Y 00 (ϑ, ϕ) = 1
2 √ π , (7)
Y −11 (ϑ, ϕ) =
√ 3
8π sinϑe−iϕ, Y 01 (ϑ, ϕ) =
√ 3
4π cosϑ, Y 11 (ϑ, ϕ) = −
√ 3
8π sinϑeiϕ, (8)
Y −22 (ϑ, ϕ) =
√ 15
32π sin2 ϑe−i2ϕ, Y −12 (ϑ, ϕ) =
√ 15
8π cosϑ sinϑe−iϕ,
Y 02 (ϑ, ϕ) =
√ 5
16π
( 3 cos2 ϑ− 1 ) ,
Y 12 (ϑ, ϕ) = − √ 15
8π cosϑ sinϑeiϕ, Y 22 (ϑ, ϕ) =
√ 15
32π sin2 ϑei2ϕ.
(9)
The spherical harmonics are normalized such that:∫ 2π 0 ∫ π 0 |Y m` (ϑ, ϕ)|2 sinϑdϑdϕ = 1 ∀`,m. (10)

heading: C CALCULATION OF INVARIANT FEATURES
text: One can obtain scalar invariats from the covariant features f̂ (Anderson et al., 2019):
• Take the component ` = 0: ξ1(f̂) = f̂`=0. • Calculate the scalar product with itself: ξ2(f̂`) = Re[ξ̃2(f̂)] + Im[ξ̃2(f̂)], where ξ̃2(f̂) =∑` m=−`(−1)mf̂m` f̂ −m ` .
• Calculate the SO(3)-invariant norm: ξ3(f̂`) = ∑` m=−` f̂ m ` ( f̂m` )∗ , where ∗ denotes the
complex conjugate.
The invariant components are then concatenated f inv ← Tinv = ξ1(f̂)⊕ (⊕L `=0 ξ2(f̂`)⊕ ξ3(f̂`) ) .

heading: D PROBABILITY DISTRIBUTION FOR ORIENTATION
text: In the main paper, we propose the expression in Eq. (5) for the distribution p(x̃|d, e, f, s). An alternative, equally valid expression is
p(x̃|d, e, f, s) = ∣∣∣∣∣ Lmax∑ `=0 ∑̀ m=−` 1√ k r̂m` Y m ` (x̃) ∣∣∣∣∣ 2 , (11)
where the term 1/ √ k with k = ∑Lmax `=0 ∑` m=−`|r̂m` |2 normalizes the distribution. We found experimentally that an agent using this expression performs worse when generating molecular structures featuring complex geometries. We hypothesize that this is because the distribution cannot get peaked enough for Lmax ≤ 5. As larger Lmax would result in a significant increase in computational complexity, we chose the expression in Eq. (5) over that in Eq. (11).
The normalization constant Z in Eq. (5) is estimated via Lebedev quadrature with 1730 angular grid points (Lebedev, 1975; 1977). We sample from the distribution in Eq. (5) using rejection sampling (Bishop, 2009) with a uniform proposal distribution q(x̃) = 14π . In rejection sampling, one first draws a sample from x̃0 from q(x̃). Then, one generates a random number u0 from the uniform distribution over [0,Mq(x̃0)], where M is such that Mq(x̃) ≥ p(x̃|d, e, f, s). We determine M by evaluating p(x̃|d, e, f, s) on a uniform grid on S2 employing a Fibonacci ‘sunflower’ grid. Finally, if u0 > p(x̃|d, e, f, s) then the sample is accepted. We ran an experiment to compare the COVARIANT agent as described in the main paper with an agent employing Eq. 11 for the distribution p(x̃|d, e, f, s). In Fig. 7, the hypothesis that the distribution in Eq. 11 cannot get narrow enough is confirmed. After 40 000 steps, the online performance of the alternative agent COVARIANT (ALT.) converges to around 0.5, which is significantly lower compared to COVARIANT. Note that the difference is smaller when considering the offline return because the estimated global mode for each distribution could still be similar.

heading: E BASELINES
text: E.1 OPT AGENT
Below, we detail the algorithm of the OPT agent. At the beginning of each experiment, the agent is given a canvas C0, a bag B0, and a black-box function that can compute the energy E(C) and the atomic forces F (C) for a given canvas. We assume a total charge of zero and a low-spin configuration. At the end of each experiment, we compute the total reward obtained for the final structure on canvas CT and report the total number of energy and gradient computations.
1. If the canvas Ct is not empty, randomly choose a focal atom f from the list of available atoms on the canvas. An atom is considered available if its number of neighbors is less than a predefined number that depends on its element (e.g., one for hydrogen and four for carbon). Two atoms on the canvas are neighbors if their Euclidean distance is below 1.5 Å. If there are no available atoms on the canvas, a focal atom is randomly chosen from the list of atoms on the canvas.
2. Randomly choose an element et from the bag Bt. 3. Randomly place the atom at = (et, xt) on a sphere with radial distance d = 1.1 Å around xf to obtain Ct+1,raw. If the canvas is empty, place the atom at the origin.
4. Optimize only the position of at using F to obtain Ct+1,opt. 5. Compute the energy difference ∆E(t) = E(Ct+1,opt)− [E(Ct) + E({et,0})]. 6. If ∆E(t) > 0, return et to the bag and go back to step 1.
7. Optimize canvas Ct+1,opt using F to obtain Ct+1. 8. Increment t by 1.
9. If the bag is not empty, go back to step 1.
In the experiments, the different agents need to be given a comparable computational budget to ensure a meaningful comparison of their performance. This is difficult as they use different computational resources: OPT runs on a CPU whereas INTERNAL and COVARIANT perform many of their computations on a GPU. However, we found experimentally that the quantum-chemical calculations are the most computationally expensive ones. These calculations are performed in the same way for all approaches. Therefore, we believe that by granting each approach the same number of PM6 calculations we achieve a fair comparison.
E.2 OPTIMAL RETURN
The optimal return for the single-bag tasks was derived in the following way. First, we obtained molecular structures for the complexes SOF4, IF5, SF6, and SOF6. Subsequently, we performed a structure optimization using the PM6 method. Since the undiscounted return is path-independent, we determined the return R(s) by computing the total interaction energy in the canvas C, i.e.
R(s) =  |C|∑ i=1 E({ei,0}) − E(C). (12)

heading: F EXPERIMENTAL DETAILS
text: F.1 COMPUTING INFRASTRUCTURE
Experiments were run on an Intel Xeon E5-2650 v4 2.2GHz 12-core processor (96GiB RAM) and an Nvidia P100 GPU (16GiB). Our agent is implemented in the deep learning framework PyTorch (Paszke et al., 2019). Data analysis was performed with the Python libraries matplotlib (Hunter, 2007) and pandas (McKinney, 2010).
F.2 IMPLEMENTATION DETAILS
The model architecture is summarized in Table 2, where the dimensions of sinv and sinvf,e are dinv = (Lmax + 2) · τ · 2 and dinvf,e = (Lmax + 2) · τe · 2, respectively. If possible, we made similar architectural choices as Simm et al. (2020), e.g. regarding the number of hidden units/layers, activation functions, and initialization schemes. We initialize the biases of each network with 0 and each weight matrix as a (semi-)orthogonal matrix. After each hidden layer, a ReLU non-linearity is employed. As explained in the main text, both MLPf and MLPe use a masked softmax activation function to guarantee that only valid actions are chosen. To model the distance d, we employ a Gaussian mixture model consisting of M = 3 Gaussians. As we treat the standard deviations {σm}3m=1
as global parameters, the MDN has 6 outputs. Further, we rescale the means µm ∈ [−1, 1] to µm ∈ [dmin, dmax]. If the sampled distance is negative, we clip the value at 0.001. Hyperparameters for CORMORANT are listed in Table 3. In our experiments, we found it important to use multiple filters τe per element (e.g. 4) and to set Lmax = 4. This gives the model enough flexibility to represent complex spherical distributions while remaining computationally tractable. For more details on CORMORANT, see the original work (Anderson et al., 2019). Further hyperparameters used in our experiments are in Table 4. PPO is known to be relatively robust with respect to the choice of hyperparameters, and we found the default values to be sufficient in most cases. Within the actor, the scaling parameter β is important to avoid that the spherical distribution approaches a delta distribution. Note that values of β can vary significantly across experiments and might require some tuning. Lastly, the right number of samples S for the global mode estimation of the spherical distribution generally depends on the shape of the distribution. In particular, we would expect that more samples are required as the distribution becomes more peaked. Since we avoid pathological behaviors by scaling the distribution with β, we found S = 1024 to be sufficient for all our experiments.
G ADDITIONAL RESULTS
In Fig. 9, a selection of molecular structures generated by the three agents during offline evaluation is shown. The agents are trained on a distribution over bags with B∗ = C7H10O2 and tested on sets of larger, out-of-distribution bags {C6H14O3,C7H16O,C7H16O2,C8H18O} (cf. Fig. 6 in the main text). From visual inspection of the structures, it can be seen why OPT fails at this task: it tends to generate molecular clusters, often containing H2. Similarly, INTERNAL commonly builds H2O molecules instead of constructing a single organic molecule out of the atoms in the bag. By contrast, COVARIANT often builds valid molecules. Further, its generated structures are more often “branched” than those of INTERNAL, indicating a higher degree of complexity.
Next, we assess the ability of COVARIANT to generate solvation clusters—a type of molecular structure that cannot be built with graph-based approaches. Following Simm et al. (2020), we task the agent to place 5 water molecules around a formaldehyde molecule that is already on the canvas at the beginning of each episode. In addition, the reward function is augmented with a penalty term for placing atoms far away from the center, i.e. r(st, at) = −∆E−ρ‖x‖2, where ρ is a hyper-parameter that is set to 0.01 (see Simm et al. (2020) for details). Therefore, the agent needs to place the water molecules such that hydrogen bonds can be formed between water molecules and between water molecules and the solute.
From Fig. 10, it can be seen that COVARIANT can solve this task by constructing stable H2O molecules and placing them in the vicinity of the solute. From visual inspection of the generated structures, it can be observed that in many cases COVARIANT arranges the molecules such that intermolecular bonds can be formed. However, it should be noted that the quantum-chemical method used in the reward function is not very well suited for modeling these interactions. Finally, Fig. 10 shows that while INTERNAL learns faster at the beginning of training, COVARIANT is slightly outperforming INTERNAL towards the end.

heading: H RUNTIME EVALUATION
text: We compared the runtimes between OPT, INTERNAL, and COVARIANT. For instance, for the singlebag task with the bag C3H5NO3, T = 240 steps of the last rollout took COVARIANT and INTERNAL on average 12 and 11 seconds (s), respectively. The final offline evaluation took the agents on average 4 and 1s, respectively. This speed difference is mainly due to the relatively slow rejection sampling procedure in COVARIANT. Each iteration, policy optimization took on average 2 and 6s for the agents COVARIANT and INTERNAL, respectively. In this case, INTERNAL is slower than COVARIANT as it performed around twice as many epochs during optimization due to early stopping. Since there is no training for OPT, this agent was overall faster than the others. Further, we note that the largest fraction of time was spent on the quantum-chemical calculations which are the same for all agents. The time the quantum-chemical calculation takes to converge depends not only on the size but also on the geometry of the input structure. The entire experiment took COVARIANT approximately 4 hours, INTERNAL 5 hours, and OPT 3 hours.

Review summaries:
The paper proposes an actor-critic neural network architecture for autoregressive generation of 3D molecular structures with reinforcement learning (RL). It builds upon the RL approach by Simm et al. (2020) which makes use of internal coordinates in order to deal with the symmetries that occur when placing atoms in the molecular design process.

This work presents an approach for 3D molecular design using reinforcement learning that exploits rotational symmetries in molecular conformations. The formulation involves an MDP that selects atoms from a “bag” and positions them in 3D space. The reward function is based on PM6 energies to encourage the generation of low-energy structures. The manuscript is well-written and provides ample context and appropriate references.