Matrix of orthogonalized atomic orbital coefficients representation for radicals and ions Open Access

Molecular (also often called chemical or quantum) machine learning (ML) has rapidly entered the domain of chemical science, facilitating both the development of simulation techniques and the design of improved reagents, catalysts, functional materials, etc.^1–8 Molecular representations are the prerequisite to machine-learn chemical properties, as they uniquely encode information about molecular composition and structure into a numerical format.^9–11 These representations can be based on simple molecular properties (such as molecular weight and van der Waals volume,¹² number of heteroatoms,¹³ partial charges,¹⁴ lipophilicity,¹⁵ etc.; also known as fingerprints or descriptors), molecular coordinates, graphs,¹⁶ topologies,¹⁷ and images¹⁸ and can take a form of a string (e.g., SMILES¹⁹ and SELFIES²⁰), a vector, or a matrix. Among them, the physics-inspired representations⁹ arguably aim to encode molecular geometries and compositions most comprehensively. Some of these representations are built from atom-centered continuous basis functions, such as the smooth overlap of atomic positions (SOAP)²¹ and the atomic cluster expansion (ACE).²² Others use potentials, such as the Coulomb Matrix (CM),⁹ the spectrum of London and Axilrod-Teller-Muto (SLATM) potential,²³ and the Faber-Christensen-Huang-von Lilienfeld representations (FCHL18,19),^24,25 or transform input structural data into an internal coordinate system, for example, in the many-body tensor representation (MBTR).²⁶ However, these representations only consider the nuclear charge and position of the nuclei in space, assuming charge neutrality. As a result, they violate the injectivity requirement of machine learning representations, i.e., are unable to distinguish between compounds with distinct electronic configurations but identical atomic compositions and geometries (Fig. 1). To satisfy the injectivity condition, several so-called quantum-inspired representations were developed by including additional electronic structure information. Molecular orbital basis machine learning (MOB-ML)²⁷ and the F (Fock), J (Coulomb), and K (exchange) matrices (FJK) representation²⁸ use post-Hartree-Fock molecular orbital properties, i.e., they require an a priori ab initio computation and thus operate in a ∆-ML fashion. In contrast, the spectrum of approximated Hamiltonian matrices (SPA^HM)²⁹ from Corminboeuf and colleagues employs a guess (i.e., pre-Hartree-Fock) electronic Hamiltonian, making it a simpler and quicker way to encode not only the composition and geometry but also the charge, spin multiplicity, and electronic state of a molecule.

In this work, we present matrix of orthogonalized atomic orbital coefficients (MAOC), a quantum-inspired representation that uses meta-Löwdin orthogonalized atomic orbitals to generate localized molecular orbital coefficients. MAOC is an atomic and molecular representation that can describe charged and open-shell compounds and distinguish molecules with identical structures but distinct electronic states. We tested the performance of MAOC in comparison with other coordinate- and Hamiltonian-based representations using kernel ridge regression machine learning model for predicting orbital and single-point energies for a broad range of conventional and redox-active species from known and newly constructed datasets.

Because the size of MAOC is related to the number of orbitals used to represent an atom, representation size grows quickly with compound size. To address this issue, a dimensionality reduction technique—principal component analysis (PCA)—can be used to reduce the size of MAOC, making it more compact, as well as easier to generate and use in machine learning tasks. This version of MAOC is denoted as PCX-MAOC, where X is the number of principal components.

The code for generating MAOC and PCX-MAOC is freely available from https://github.com/hits-ccc/MAOC/tree/main/Codes/MAOC_mol_rep.

Three sets of non-overlapping chemical data were used to test MAOC in supervised machine learning. The QM7b dataset^38,39 was selected since it is among the most widely used datasets for machine learning in chemistry. It consists of 7211 relatively small neutral molecules with no conjugated bonds within cyclic moieties. To evaluate the performance of MAOC for radicals and ions, geometries of the anionic, cationic, and dicationic forms of the compounds in the dataset were optimized and their various vertical and adiabatic properties were computed in this work; strongly spin-contaminated species were eliminated from the dataset (see Computational procedures below). Therefore, in total 7197 geometry-optimized anion radicals, 6999 geometry-optimized cation radicals, 7198 geometry-optimized dications, as well as 7208 anion radicals and 7208 cation radicals in the geometry of the parent neutral molecule were added to the original QM7b dataset of neutral molecules. We denote this extended dataset as QM7b^X.

The N-HPC-1 dataset, constructed in this work,⁴⁰ consists of 3735 nitrogen-doped polycyclic compounds that were generated by introducing up to 11 nitrogen atoms into prototypical polyaromatic skeletons (Fig. 2). For each generated compound, canonical SMILES were produced, and any duplicates were removed from the dataset. To generate respective open-shell compounds, one, two, or three electrons were removed, or one electron was added to the neutral molecules to produce eight groups within the dataset: neutral singlets, neutral triplets, anionic doublets, cationic doublets, dicationic singlets, dicationic triplets, tricationic doublets, and tricationic quartets. Stable open-shell N-heteropolycycles display a range of practically useful electronic,⁴¹ magnetic,⁴² and optical⁴³ properties, whose facile prediction with a tailor-made representation could provide valuable molecular design guidelines.

The REDOX dataset, also constructed for this study,⁴⁴ contains 4146 neutral redox-active molecules⁴⁵ with 1–4 unpaired electrons (Fig. 2). The dataset consists of several classes of compounds carrying diverse functional groups: organic radicals (nitroxyl, phenoxyl, and galvinoxyl), carbonyl compounds (quinones, carboxylates, and phenazine-derived radicals), and cyanides. These systems are widely utilized across chemical sciences in synthesis, characterization, radical trapping, organic batteries, etc. Nitroxyls undergo one electron oxidation to form oxoammonium cations, while all other species in the dataset instead favor one electron reduction toward the corresponding anions and anion radicals. With the inclusion of these adducts, the entire REDOX dataset thus comprises 4146 neutral radicals, as well as 4018 anionic and 687 cationic open- and closed-shell molecules.

The QM7b^X, N-HPC-1, and REDOX datasets can be freely obtained from GitHub (https://github.com/hits-ccc/MAOC/tree/main/Datasets).

Geometries of all systems in the N-HPC-1 and REDOX datasets and of all charged systems in the QM7b^X dataset were optimized at the PBE0-D3/def2-TZVP level of theory using the ORCA 5.0 package.⁴⁶ Geometries of neutral compounds in the QM7b dataset were used to compute vertical properties at the same level of theory. In all computations, wavefunction stability checks were performed. For all open-shell compounds, the expectation value of the spin-squared operator was assessed, and species with $Ŝ2$ above 10% of the expectation value were excluded (see Sec. II in the supplementary material). Single-point energies (SPEs), the energies of the highest, lowest, and singly occupied molecular orbitals (HOMOs, LUMOs, and SOMOs, respectively), as well as HOMO-LUMO gaps for the closed-shell and SOMO-LUMO gaps for the open-shell species, were also computed at the PBE0-D3/def2-TZVP level of theory.

Global SLATM representations²³ were generated using QML code⁴⁷ with the variables set to sigmas = [0.05, 0.05], dgrids = [0.03, 0.03], rcut = 4.8, and rpower = 6. Coulomb Matrix representations⁹ were generated using the DScribe package⁴⁸ and sorted using the L2 norm to ensure permutation invariance. The SOAP representations²¹ were generated using the DScribe⁴⁸ with the cutoff for the local regions (rcut) set to 6 Å, the number of radial basis functions set to 8, the maximum degree of spherical harmonics set to 100, and the spherical Gaussian-type orbitals chosen as the type of radial basis functions. The SPA^HM representations were generated using the code provided in Ref. 29with the LB guess and the MINAO basis (see also Sec. VII in the supplementary material). The MAOC representations were generated using the localized orbital (lo) package of the Python-based Simulations of Chemistry Framework (PySCF) package.⁴⁹ Unless specified otherwise, MAOC representations were generated using the pcseg-0 basis set and ANO⁵⁰ as the reference basis. In this work, a sorted flattened version of MAOC is used. The PCX-MAOC representation is created by reducing the size of an unflattened MAOC array from M × M to M × 1, M × 2, …, M × X (where X < M) using the PCA dimensionality reduction technique. In this work, the number of principal components, X, is set to 3 (see Sec. III in the supplementary material), and the representation is denoted as PC3-MAOC.

A DELL XPS 15 9510 laptop with a processor 11th Gen Intel(R) Core(TM) i7-11800H @ 2.30 GHz, 64-bits operating system (x64-based processor), and 32.0 GB of installed RAM was used to generate all representations. The CPU timings were measured using the datetime built-in Python library in a Jupyter notebook. The storage space requirements were evaluated as the amount of space needed to store the NumPy ndarrays.

A kernel ridge regression (KRR) model was used to test MAOC and other representations in predicting the properties of compounds in the QM7b^X, N-HPC-1, and REDOX datasets. All learning curves were generated using KRR with a Laplacian kernel, and the sigma and lambda hyperparameters were optimized for each representation and training set size using a two-fold cross-validation splitting. The sigma values examined ranged from 1 to 8000. The grid step size was set to be between 100 and 500, depending on the size of the dataset under consideration. The lambda values were chosen to be between 1 and 10⁻¹⁰ (i.e., 1, 10⁻³, 10⁻⁵, 10⁻⁶, 10⁻⁷, 10⁻⁸, and 10⁻¹⁰). The Laplacian kernel was chosen after evaluating its performance vs the Gaussian kernel (see Sec. IV in the supplementary material). Unless specified otherwise, the train-to-test ratio for all properties under investigation is 80/20 and all the results in this study are the mean of five measurements using random shuffling of the data. The Cholesky decomposition (in QML code⁴⁷) is used to solve the equation $a=K+λI−1y$ of the KRR model. The step size (training set size) for the learning curves is chosen based on the dataset, ranging from 400 to 1500. All machine learning operations were performed on a MacBook Pro (2019) with a 1.5 GHz Quad-Core Intel Core i5 processor and 8 GB of 2133 MHz LPDDR3 memory.

Principal component analysis (PCA) was conducted utilizing the sklearn Python library and an “auto” SVD solver. In this work, the number of components was set to either two (when mapping the chemical space) or three (when constructing PC3-MAOC), and the explained variance ratio was recorded close to the axis of the PCA plots. t-distributed Stochastic Neighbor Embedding (t-SNE) was performed using the sklearn Python library. The number of components was set to two, and all the other values were left as defaults.

Since MAOC uses orthogonalized atomic orbitals to construct the guess localized molecular orbitals, the basis set used to generate these atomic orbitals affects the matrix generation time, the training time for machine learning models, the amount of space needed to store the representation, and the quality of the property predictions. To identify the optimal choice, we tested several basis sets of varying sizes and flexibilities [Fig. 3(a)]. Generating MAOC for the N-HPC-1 dataset with the Karlsruhe (def2-SVP) or Dunning (cc-pVDZ) basis sets requires significantly more time and storage space than with the Pople, Jensen (pc-0 and pcseg-0), and LANL2DZ basis sets. STO-3G produces the smallest representations file in a relatively short time. In comparison to other representations (Table I), generation of MAOC requires the longest time; however, even for the largest studied dataset, QM7b^X, this is only ∼two CPU-hours. Moreover, both MAOC and SPA^HM are well parallelized, i.e., MAOC for QM7b^X requires just over 15 min on an eight-CPU laptop. Executing principal component analysis on the generated MAOC is virtually instant. In terms of storage space, PC3-MAOC and SPA^HM are by far the most compact representations, while full MAOC is roughly in the same order of magnitude as SLATM and SOAP.

The learning curves in Fig. 3(b) illustrate the accuracy of the KRR machine learning model for the prediction of HOMO- and SOMO-LUMO gaps for a random sample of 4000 compounds from the N-HPC-1 dataset using MAOC built from various basis sets in terms of mean absolute error (MAE). The train-to-test ratio is 70/30, and the sigma values for each set of training data are optimized. Surprisingly, the highest accuracy is achieved with relatively small basis sets (3-21G, 6-31G, and pcseg-0). However, results with all tested basis sets ultimately fall within a narrow range of 0.31–0.35 eV in MAE. We note that adding diffuse and polarization functions to the basis sets increases the computational costs but does not affect (polarization functions) or even reduces (diffuse functions) the accuracy of ML predictions and, therefore, is not advised [see Fig. 3(a) and Secs. V A and V B in the supplementary material]. The reason behind the counterintuitively superior performance of smaller basis sets is that MAOC constructed from larger basis sets includes many more virtual orbitals, which, as shown in Section E, do not carry significant weight in machine learning of the studied properties. While this is the case for organic molecules in this work, it might not always be: polarization functions are likely important (and must thus be included when constructing MAOC) for molecules containing heavy metals, etc. Finally, using principal component analysis, we evaluated how well the “atomic” space of chemical elements across the Periodic Table is mapped with MAOC depending on the basis set (see Sec. V C in the supplementary material). Ultimately, pcseg-0 was chosen as an optimal basis set for constructing MAOC representations in this work, but a similarly good ratio of computational cost to accuracy can be achieved with 3-21G and 6-31G.

By construction from orthogonalized atomic orbitals, MAOC can serve not only as a molecular but also an atomic representation, differentiating single-atom systems and enabling prediction of their properties. This contrasts with many conventional molecular representations, which either rely on many-atom environments or, for a single atom, produce a representation unsuitably small for machine learning. Instead, prediction of monatomic (monoatomic) properties employs such traditional descriptors as electronegativity and atomic radii,⁵¹ which, however, are not uniquely defined.^52,53 Instead, MAOC offers a rigorously defined method to encode single atoms across the entire Periodic Table for machine learning tasks. With MAOC, the “atomic” space is well resolved and the proximity between two atoms reflects the similarity in their electronic (orbital) configurations (Fig. 4).

Coordinate-based representations, for example, Coulomb Matrix, SOAP, and SLATM, distinguish compounds based on their compositions and geometries, while the quantum-inspired representations SPA^HM and MAOC also account for electronic structure features. As a result, only the latter can resolve the chemical space of molecules with similar geometries but nonidentical charges and spins [Fig. 5(a)]. This is further illustrated in the t-SNE plot of the chemical space of anionic and cationic compounds in the geometries of the parent neutral species from the QM7b dataset [Fig. 5(b)].

In this work, we focus on the fundamental properties of charged (and) open-shell compounds: energies of the frontier molecular orbitals, energy gaps between SOMO or HOMO and LUMO, and ground state single-point electronic energies. We first discuss the results for the commonly available QM7b dataset,^38,39 which consists of small neutral organic compounds composed of elements from the first to third rows of the Periodic Table and mostly lacking conjugated π-bonds in cyclic moieties. To create an extended QM7b^X dataset, anion radicals, cation radicals, and dications were considered for each parent molecule in QM7b; for the open-shell compounds, both vertical and adiabatic properties were evaluated. Overall, the best predictive performance of kernel ridge regression in terms of mean absolute errors (MAEs) and root-mean-square errors (RMSEs) for the orbital energies is achieved with SLATM, SOAP, and SPA^HM representations, followed by MAOC, PC3-MAOC, and Coulomb Matrix [see Table II and Figs. 6(a) and 6(c)]. Compared with geometry-based representations, their quantum-inspired counterparts perform worse in predicting the single-point energies [Fig. 6(b)] but are generally better in predicting the orbital energies of the charged open-shell compounds. For example, SPA^HM and PC3-MAOC are the best representations for machine-learning the SOMO-LUMO gap of anion radicals [Fig. 6(d)] and cation radicals [Fig. 6(e)], respectively. We have also combined the data for non-optimized (vertical) anion and cation radicals into a single subset to test the representations on systems with identical geometries but distinct electronic states. In this case, only the quantum-inspired representations, SPA^HM and MAOC, are able to learn from the data, with SPA^HM delivering the lowest MAEs [Fig. 6(f)].

For the N-HPC-1 dataset⁴⁰ (Fig. 2) comprising neutral N-heteropolycycles, their open-shell (doublet) anions and cations, closed- and open-shell (triplet) dications, and open-shell (doublet and quartet) trications, CM generally performs worse than other representations (Table III and Fig. 7). Performance of the latter varies: quantum-inspired representations offer the best accuracy in predicting the single-point energies, while SOAP and SLATM generally produce better orbital energies. However, in several cases, such as singlet dications and doublet trications, PC3-MAOC outperforms other representations, including the full MAOC (see Sec. III E). Finally, the performance of MAOC for the N-HPC-1 dataset containing π-conjugated compounds is improved (i.e., is associated with lower MAEs) relative to QM7b^X, which is comprised mostly of saturated systems. To explain this, we turn to the way MAOC is constructed. The localized orbitals of conjugated π-bonds are no longer atom-centered, but rather are distributed over multiple atoms and account for the long-range interactions between them. This captures the global character of MAOC.

Next, we assess the performance of the kernel ridge regression with various representations for the REDOX dataset⁴⁴ (Table IV and Fig. 8), which includes redox-active molecules from diverse chemical classes and with a range of functional groups (Fig. 2). The errors in property predictions for this dataset are markedly higher than for QM7b^X and N-HPC-1 due to the greater chemical diversity within and non-combinatorial generation of the REDOX dataset. Each class of redox-active species is represented by ∼1000 compounds formed from a variety of structural motifs and functional group combinations, resulting in a much broader range of values for a given property compared with the QM7b^X and N-HPC-1 datasets (see Sec. VIII in the supplementary material) and limiting the model’s ability to “learn”. Nonetheless, the quantum-inspired representations—SPA^HM, MAOC, and PC3-MAOC—as well as the geometry-based SLATM, outperform CM and SOAP in all cases. SPA^HM predicts significantly better single-point energies, while MAOC and PC3-MAOC afford the best accuracy for frontier orbital energies. Functional groups, present in REDOX but absent in the N-HPC-1 dataset, enhance the global character of MAOC through resonance, which results in molecular orbitals that are centered not on a single bond but distributed over at least three atoms.

Finally, to further test the generalizability of MAOC, we applied it to predict the HOMO-LUMO gaps in a dataset containing in total 5000 compounds, randomly selected from the QM7b, N-HPC-1, and REDOX datasets. We find that the performance of MAOC and PC3-MAOC for this “mixed” dataset average between the corresponding errors for the individual datasets (see Sec. VI in the supplementary material). This further demonstrates that MAOC is a suitable representation for compositionally, structurally, and electronically diverse sets of chemical compounds.

For the N-HPC-1 dataset of N-heteropolycycles, the PCA-reduced version of MAOC, PC3-MAOC, outperforms the full representation in predicting all considered properties. In search for a reason behind this puzzling behavior, we have analyzed the orbital contributions to the three principal components of PC3-MAOC. Specifically, for all species in a given subset, localized molecular orbitals contributing the most to each principal component (above 0.2 variance) are listed. We then count how many times an orbital of interest, e.g., HOMO-2, HOMO-1, etc., occurs in this list. Finally, this occurrence count is normalized to the total number of L-MOs in the list. For example, six frontier and near-frontier localized MOs—from HOMO-2 to LUMO+2—constitute ∼11% of all orbitals with variance above 0.2 in PC1-3 of MAOC for the parent neutral compounds from the QM7b dataset (Fig. 9). Among the three studied datasets, N-HPC-1 stands out in that its MAOC principal component lists feature HOMO and LUMO much less often than the respective higher- and lower-lying orbitals (HOMO-2, HOMO-1, LUMO+1, and LUMO+2). In contrast, all six frontier and near-frontier MOs occur with similar frequency in the >0.2 variance lists of the QM7b^X and REDOX datasets. An in-depth quantitative analysis of these contributions to the molecular representations and their effects on the quality of machine learning predictions is planned for a future study.

Principal component analysis of the MAOC representation allows for an elaboration of how the electronic structure governs the chemical properties. For example, in line with the chemical intuition, predicting orbital energies of the anionic compounds requires information on the occupied molecular orbitals, while properties of cations are instead dominated by the unoccupied MOs (Fig. 9). While these relationships are straightforward when predicting the energy levels of molecular orbitals themselves, as is done in this study, less obvious effects can potentially be uncovered through such an analysis of more complex properties, such as reaction energies, EPR spectra characteristics, and excited state features.

In this work, we present a quantum-inspired molecular representation for chemical machine learning tasks—the matrix of orthogonalized atomic orbital coefficients (MAOC). It utilizes not only the structures and compositions of the systems it represents but also their electronic information, i.e., charge and spin multiplicity. The matrix of orthogonalized molecular orbitals is constructed from the coefficients of the non-optimized (guess) localized molecular orbitals, which in turn are linear combinations of predefined orthogonalized atomic orbitals. In generating the latter, even very small atom-centered basis sets, such as pcseg-0 and STO-3G, afford the accuracy of larger basis sets at a fraction of the computational cost. MAOC is a representation with dual character: local when the core and lone pair orbitals of atoms are constructed, and global when the orthogonalized atomic orbitals are combined to form molecular orbitals, which can either be localized or distributed over a part or even an entire molecule. MAOC is uniquely suited to represent not only molecules but also monatomic (by construction) and periodic (using meta-Löwdin orthogonalization scheme) systems. Furthermore, MAOC allows learning and predicting the properties of systems with identical compositions and geometries but distinct electronic configurations, such as vertical properties or redox properties of rigid compounds. Finally, we also present a reduced version of MAOC based on the principal component analysis dimensionality reduction technique, PCX-MAOC, which possesses all of the above-mentioned features of full MAOC while being significantly more compact, enabling further analysis of the connection between electronic structure and chemical properties, and, in certain cases, leading to higher prediction accuracy.

The performance of the kernel ridge regression ML model, in conjunction with MAOC, PC3-MAOC, and several other coordinate- and Hamiltonian-based representations, was tested in predicting the frontier orbital energy levels and ground state single-point energies for an extended QM7b dataset and two new datasets, N-HPC-1 and REDOX. We observed that MAOC and PC3-MAOC outperform geometry-based representations for sets of compounds with identical geometries but distinct electronic configurations and for a chemically diverse set of redox-active compounds, while affording accuracy similar to that of SOAP, SLATM, and SPA^HM for a set of N-heteropolycycles. Further improvements in MAOC’s performance are likely achievable with larger datasets (in cases where property values are spread over a broad range) and with other machine learning architectures, such as neural networks.

The supplementary material contains additional computational and data analysis details, additional data, and a full set of learning curves.

The authors gratefully acknowledge support from the Klaus Tschira Foundation, funding from the German Research Foundation (DFG) through the collaborative research center No. SFB1249 (Project No. 281029004-SFB 1249, Project C10), and funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 101042290 PATTERNCHEM), as well as access to computational resources from the state of Baden-Württemberg through bwHPC (bwForCluster JUSTUS2). The authors also thank Ksenia Briling and Professor Clemence Corminboeuf for their help with the SPA^HM representation, Professor Olexandr Isayev and Professor Michele Ceriotti for fruitful discussions, and Dr. John Lindner for proofreading the manuscript.

The authors declare no conflicts to disclose.

Stiv Llenga: Conceptualization (equal); Formal analysis (lead); Investigation (lead); Methodology (lead); Visualization (equal); Writing – original draft (lead). Ganna Gryn’ova: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Supervision (lead); Visualization (equal); Writing – review & editing (lead).

The data that support the findings of this study are available within the article and its supplementary material, as well as from GitHub at https://github.com/hits-ccc/MAOC/tree/main/Datasets.