A basic electro-topological descriptor for the prediction of organic molecule geometries by simple machine learning

This paper proposes a machine learning (ML) method to predict stable molecular geometries from their chemical composition. The method is useful for generating molecular conformations which may serve as initial geometries for saving time during expensive structure optimizations by quantum mechanical calculations of large molecules. Conformations are found by predicting the local arrangement around each atom in the molecule after trained from a database of previously optimized small molecules. It works by dividing each molecule in the database into minimal building blocks of different type. The algorithm is then trained to predict bond lengths and angles for each type of building block using an electro-topological fingerprint as descriptor. A conformation is then generated by joining the predicted blocks. Our model is able to give promising results for optimized molecular geometries from the basic knowledge of the chemical formula and connectivity. The method trends to reproduce interatomic distances within test blocks with RMSD under $0.05$ \r{A}


Introduction
This work assesses the problem of generating reliable conformers of molecules from proposed chemical compositions. Realistic initial bond lengths and angles are essential for efficient geometry optimizations. They are normally the first step of the usual computational workflow of systematic variations of the atomic coordinates inside a molecule and the calculation of the potential energy and forces of the system in order to find a minimum value for the potential energy, which indicates a theoretically optimized conformation. Moreover, approximate and reliable molecular geometries serve for many modelling and process simulation purposes from docking to Molecular Dynamics.
The function combining all possible variations of the geometry and the potential energy forms a high-dimensional surface and it is the well known potential energy surface (PES). All possible conformers for a given compound are comprehended as minima in the appropriate PES. Two processes are required in order to obtain the model of the conformer with the lowest potential energy of a given molecule: (1) a procedure to attain its corresponding and plausible PES and (2) a method to navigate it to search for minima. Procedures to find PES's and then computing potential energies can vary in efficiency depending on several factors, and mostly the number of atoms of the compound. Options exist from the easily computed empirical force fields based on classical considerations of bodies in a molecule, 1,2 passing fast quantum mechanical semi-empirical calculations being parameterised for specific scenarios, 3 to more general and reliable but computationally expensive ab-initio and DFT calculations. 4 The method to navigate the PES can also be computationally expensive depending on how fast the global minimum can be found. 5 Obviously, if the initial choice of a guess conformation is near to the final optimized structure, the correspnding PES's minimum will be faster reached. This is crucial for accurate geometry optimizations of large molecules.
Machine learning (ML) has been recently used in a variety of topics inside the field of quantum mechanics.  For molecular geometry optimizations, progress has been made mostly by using neural networks to parameterise classical force fields. 38,39 The proposal made in this work is to use ML, and specifically the Kernel Ridge Regression algorithm, for predicting molecular conformations by producing the local arrangement of each atom belonging to a molecule. It can be achieved by using a large and confident database of optimised small molecules 40 as a source for both the training and testing sets. For this purpose, certain molecular structural blocks are characterised and defined by an electro-topological descriptor 18 45 with a new ab-initio torsion angle database, we join all predicted blocks to produce the desired molecular structures with reasonable reliability. Results promise fast and confident predictions of molecular geometries and conformations from their formulas taken as structural graphs.

Learning data
The database used in this work is part of a larger one 46,47 of quantum PES's minima geometries of small organic molecules containing up to 8 C, O, N and/or F atoms. The optimised geometries of this database are reported to be found using DFT/B3LYP 48 with the 6-31G(2df,p) basis set as a commonly accepted reliable PES. We will refer to this database as 8CONF. The size of this resulting database subset is 21k molecules. 40 To facilitate predictions based on this data, we seek a representation which minimizes the amount of redundant information in the learning set, and also could group together similar kinds of data.
First, each molecule is split into blocks. A block is the main building part of our model, and it is characterised by: (1) a central atom with more than one bond and (2) the first neighbors of such central atom. Figure 1 shows a block decomposition for an example molecule from the 8CONF database. Note that each atom can normally be included in multiple blocks: the block centered around itself as well as each of the blocks surrounding its neighboring atoms. Atoms with only one neighbor are not considered to define a block. Blocks therefore have from two to four neighboring atoms in the selected molecular sets where all atoms belong to the first and second rows of the periodic  Figure 1: Block decomposition for the molecule C 3 H 8 . Each C atom is bonded to multiple atoms and hence defines a block. The molecule can therefore be divided in 3 blocks, two of which (for atoms 1 and 3) belong to the same block-class. In each block we have redundant atoms to indicate where blocks join together.
A unique Cartesian representation is not well suited for predictions because coordinate values depend on the chosen reference center. Instead, we represent local coordinates within a block by 1) subtracting the molecular Cartesian coordinates of the central atom position in that block to define it as the local coordinate origin, and 2) computing the matrix B of scalar products a i · a j between each pair (i, j) of coordinates of local position vectors corresponding to the non-central atoms in the block. This matrix is the feature we use for training and predictions.
The matrix B of scalar products is symmetric and at most 4 × 4 in size, and therefore has up to ten unique degrees of freedom. B contains enough information to rebuild the set of molecular Cartesian coordinates (see Appendix 8.1) for a block except for translations, rotations, and chirality, with which the matrix is invariant.
The non-central atoms in a block have no natural ordering. Hence a way must be chosen to assign an index i to each of them with a minimum of ambiguity for ML training and testing. To this end we define an equivalence relation for the set of all blocks, i.e., each block belongs to a single, specific equivalence class or block-class. The ML algorithm is then independently trained for each block-class.
Two blocks belong to the same class if 1) the species of the central atom in each block is the same, 2) the species of each neighboring atom is the same (some ambiguity is solved by sorting by atomic numbers) and 3) the arrangement of the atoms in space is the same, i.e. either tetrahedral (TH), triangular (TR) or linear (L) .
The definition of block-classes can be applied with different levels of restrictions. As a result, the number of different block-classes can vary, as blocks with the same atoms can appear in very different environments.
We denote a block-class by a series of chemical symbols followed by certain indicators of spatial configurations when necessary. The first symbol is that of the central atom, followed by the symbols of the neighboring atoms ordered with higher atomic numbers first. Figure 2 shows the distribution of blocks inside the database.  However doing so also decreases the size of the learning sets, resulting in block-classes with few or no members. In general block-classes must be defined in such a way to get neither too many nor too few.

Descriptors
To complete our learning data, a property is required that can correlate with the desired feature. This property is named as a descriptor or fingerprint. In this case, desired features are the scalar products of a block being directly related to the pursued molecular geometry predictions.
Descriptors needs to fulfill a number of required characteristics such as: (1) easiness to establish or compute, (2) good representability and (3) low dimensionality. 49 We took the so-called electro-  training to grasp correlation between the descriptor and the property, and T was used for that. Let R be the scalar product tensor of all blocks inside a block-class in the T. In matrix notation the training process can be defined as: where C represents the coefficients to be computed, are the kernel terms as calculated by using L 1 norm distances | · | between the descriptors d and α is a regularization parameter to be commented below. From solving eq. 1, we obtain the coefficients C, used to get the prediction P by: where D is the matrix of distances between descriptors d of known j learning points and the new i The KRR method depends on two so-called hyper-parameters: λ to control the regularization factor of the kernel K, and α which controls the radius of inclusion/similarity for the kernel func- A new set of Cartesian coordinates can be predicted with our algorithm as outlined below.

Scheme 1 shows the flow diagram. A SMILES 51 representation of the molecule is taken as input.
It is then analyzed to extract connectivity, the blocks and their associated classes. The corresponding e-states are then computed. The learning data-set is built by using the block-classes to compute the scalar product and the e-states in order to train the KRR algorithm. The newly parameterized KRR is then used to predict each scalar product. Finally, the blocks are reconstructed again from A. However, methods must be also evaluated by their worst cases as some block-classes rose them to around 0.3Å. Among the highest RMSD's values, two different kinds of problematic blocks are found: those belonging to a group on which there are not enough data to properly train the KRR algorithm (e.g. N-N-N-H with only 7 occurrences) and those belonging to molecules where the blocks are located inside rings that tend to distort the bond lengths and angles (e.g. C-C-C-C-C which is a very abundant block-class with multiple combinations). The second experiment with the C testing set involved the complete rebuilding of the molecules.
After predicting the blocks and reconstructing the coordinates, TTKDG-ML procedure was fed with the distance matrices of each block. Some restrictions were put in place in the form of an arbitrary high number (the decimal value 10 was used) to be interpreted as a weight by the algorithm. So, the better conformers were left within those non modifying the predicted distances by the KRR method. The best is selected by the minimun RMSD value. Histograms comparing both RMSD results and potential energy calculations using a simple Effective Medium Theory (EMT) 54,55 for the 450 molecules of the test-set are shown in Figure 6.
This method consists of an expression derived for obtaining a value related with the total energy of a system of interacting atoms. It is based on an ansatz for the total electron density of the system as a superposition of atom densities taken from calculations for the atoms embedded in a homogeneous electron gas. We use these energy values here to illustrate molecular complexity compared with total RMSD predictions. As EMT was developed only for solids and not for molecules, the resulting values of energy must be only taken as just a reference. It can be realized that there are eight molecules from the test-set on which the TTKDG does not perform well. It must be originated in the fact that they contain one or more blocks with poor prediction results (See Figure 5) such as N-C-C-C.Tr and C-C-C-C-C.Td. The results show that the proposed TTKDG-ML performs better than the original ETKDG implemeted in RDkit (See Figure 6 B and D).

Conclusions
To predict accurate geometries of molecules is a complex problem able to be multiple approached. A possible solution is described here by using a machine learning tool based on a Kernel Ridge Regression routine. The described procedure departing from bonded atom blocks results in a very flexible and easy-to-expand way to describe the canonical geometry of an atom in a molecule and its environment. The descriptor used showed a very good correlation with bond distances and angles when the coordinates were transformed into scalar products. RMSD values obtained from the experiments performed on blocks support our decision on this kind of descriptors and validate the block structure. To solve most of the problematic block predictions will be required to add more related training data because the involved atoms in such blocks are rare in common molecules and therefore scarce in the training set. For those atoms belonging to rings, a better data treatment could be a solution. But this needs to be exercised with caution, because too much pruning of the data base could lead to over-fitting the KRR method results.
For approaching the final purpose of getting a complete machine learning method, the selected descriptor was tested for torsion angles. It completely failed. The symmetries involved when joining two blocks were very difficult to grasp by the e-state with the only use of the KRR algorithm.
Even the selection of the correct torsion angle for the learning process proved to be non-trivial. 8 Appendix

Scalar product formulations and reconstruction
For the sake of argument clarity, tetravalent blocks are taken here as an example for the scalar product formulation. Each block has five atoms, or a total of 15 coordinates. Rotations and translations account for six degrees of freedom (DOF's), leaving nine of them to consider.
Translations were automatically eliminated by measuring neighbouring atom coordinates as displacements from the central atom. The elimination of rotational dependencies were achieved by forming the matrix of all scalar products between the nonzero position vectors X ac , where a = 1 . . . 4 specifies direct neighbouring atoms and c is one of x, y, z: The bond lengths can be obtained directly from the four diagonal elements. The information about the angle is provided by the six distinct scalar products. The system is over determined by one DOF. However the matrix has only rank 3 because it was formed from a 4 × 3 matrix.
Therefore one eigenvalue is zero, leaving a three-dimensional eigenspace which spans the nine DOF's remaining as well determined.
The matrix can be rewritten by using its eigendecomposition as: This defines a matrixX = QΛ 1/2 which produces the same scalar products. The matrixX will contain one row of zeros, corresponding to the eigenvalue 0.
This row/column is discarded, and the remaining matrix R will be an eligible set of reconstructed positions. In weeding out numerical garbage, the lowest eigenvalue is always discarded even if it is not exactly zero. The remaining eigenvalues must be positive since the scalar product is positive definite.