Background & summary

The goal of computationally designing novel materials and molecules with desired physicochemical properties is yet to be achieved. High-throughput screening represents the most straightforward approach towards materials design1. However, it presupposes that all assumptions and approximations inherent to the employed modeling techniques are applicable to the entire chemical compound space, which is the space populated by all stable molecules or materials2. Furthermore, due to the combinatorial scaling of chemical space with molecular size, it is difficult to explore or even navigate. Conclusive insights about the domain of applicability (transferability) are lacking even for the most popular first principle quantum chemistry methods. For example, the reliability and accuracy of density functional theory is known to dramatically depend on chemical composition and atomistic configurations3, highlighting the importance of reliable experimental4 or high-level quantum chemistry state-of-the-art results5,6. Unfortunately, the systems reported are typically small, which implies the existence of severe selection bias. One can therefore question how representative they are. The problem of representative diversity has triggered the design of special purpose chemical space libraries for method validation or molecular design7–10.

Here, we report molecular structures and properties obtained from quantum chemistry calculations for the first 134k molecules of the chemical universe GDB-17 data base11, covering a molecular property set of unprecedented size and consistency. The data-set corresponds to the GDB-9 subset of all neutral molecules with up to nine atoms (CONF), not counting hydrogen. The molecular size distribution of all 134k molecules is shown in Fig. 1. This data set contains small amino acids, such as GLY, ALA, as well as nucleobases cytosine, uracil, and thymine. Also pharmaceutically relevant organic building blocks, such as pyruvic acid, piperazine, or hydroxy urea are included. Among the 134k molecules, there are 621 stoichiometries, among which C7H10O2 dominates with 6,095 constitutional isomers for which atomization energies and radii of gyration also are on display in Fig. 1.

Figure 1: Illustration of the scaling of chemical space with system size.
figure 1

For the smallest 134k molecules, with up to 9 heavy atoms CONF (not counting hydrogens) taken from the chemical universe GDB-1711, the distribution of molecular size is shown as a function of number of occupied electron orbitals, i.e. number of electron pairs, Nep=Ne/2. Each black box denotes the number of constitutional isomers for one out of the 621 stoichiometries present in the 134k molecules. The two left-hand side insets correspond to zoom-ins for smaller compounds. The right-hand side inset zooms in on the predominant stoichiometry, C7H10O2, and features a scatter plot of G4MP2 relative (w.r.t. global minimum) potential energies of atomization E versus molecular radius of gyration, Rg. Joined projected distributions are shown as well.

Full size image

For all 134k molecules, we have calculated equilibrium geometries, frontier orbital eigenvalues, dipole moments, harmonic frequencies, polarizabilities, and thermochemical energetics corresponding to atomization energies, enthalpies, and entropies at ambient temperature. These properties have been obtained at the B3LYP/6-31G(2df,p) level of theory which forms the basis for the more accurate state-of-the art Gn methods which are on par with experimental accuracy12. For the 6,095 constitutional isomers of the predominant stoichiometry, C7H10O2, we report the energetics at the significantly more accurate G4MP212 level of theory.

This report is structured as follows. We first describe the genesis of the results. Thereafter, we discuss the validation of our DFT results by comparison to (i) G4MP2, (ii) G4, and (iii) CBS-QB3 results for 100 molecules, randomly chosen out of the 134k set. This data can serve the development, training and evaluation of inductive statistical data analysis-based machine learning (ML) models13. It might also assist the search and discovery of hitherto unknown trends, structure-property relationships, and molecular materials design1,14,15.

Methods

Generation of atomic coordinates

Starting with ref. 11, we use all SMILES16 strings for molecules with up to nine heavy atoms. Cations, anions, and molecules containing S, Br, Cl, or I, have been excluded, resulting in 133,885 molecules. 1,705 zwitterions have been kept in the data due to their occurrence in small biomolecules, such as amino acids. Initial Cartesian coordinates for all molecules were generated by parsing the corresponding SMILES strings using Corina (Version 3.491 2013)17. We subsequently carried out geometry relaxations at the PM7 semi-empirical level of theory using MOPAC (Version 13.136L 2012)18. In the PM7 calculations, we invoked tight electronic and geometric convergence thresholds, using precise keyword. PM7 equilibrium geometries have subsequently been used as input for B3LYP geometry relaxations using Gaussian 0919. We iteratively refined the electronic and geometry thresholds. For the first iteration, Gaussian 09’s default electronic and geometry thresholds have been used for all molecules. For those molecules which failed to reach SCF convergence ultrafine grids have been invoked within a second iteration for evaluating the XC energy contributions. Within a third iteration on the remaining unconverged molecules, we identified those which had relaxed to saddle points, and further tightened the SCF criteria using the keyword scf(maxcycle=200, verytight). All those molecules which still featured imaginary frequencies entered the fourth iteration using keywords, opt(calcfc, maxstep=5, maxcycles=1000). calcfc constructs a Hessian in the first step of the geometry relaxation for eigenvector following. Within the fifth and final iteration, all molecules which still failed to reach convergence, have subsequently been converged using opt(calcall, maxstep=1, maxcycles=1000). calcall constructs a Hessian for all steps through the the geometry relaxation. After all these measures taken, eleven problematic molecules still failed to converge to a minimal geometry. Out of these eleven molecules, six can be converged with low threshold using the opt(loose)-keyword. In the remaining five there were two near-linear molecules which converged to saddle points with very low imaginary frequencies (ω0<i10 cm−1 for the lowest mode). In the readme.txt file of this report, all these 11 molecules are specified using their indices in the database.

In the case of the 6,095 constitutional isomers of C7H10O2, all molecules converged to local minima during the B3LYP geometry relaxation. To compute atomization energies, we have also performed spin-unrestricted calculations for all atoms with spin-multiplicities 2,3,4,3,2 for the atoms H, C, N, O, F, respectively.

Data Records

Molecular structures and properties are publicly available at Figshare (Data Citation 1) in a plain text XYZ-like format described below. Deposited files include the 133, 885 GDB-1 to GDB-9 molecules (dsgdb9nsd.xyz.tar.bz2), the 6,095 constitutional isomers of C7H10O2 (dsC7O2H10nsd.xyz.tar.bz2), the 100 validation molecules (see Table 1) enthalpies of atomization (validation.txt), and atomic reference data (atomref.txt).

Table 1 Validation of atomization enthalpies at B3LYP/6-31G(2df,p)-level.
Full size table

File format

For each molecule, atomic coordinates and calculated properties are stored in a file named dataset_index.xyz. The XYZ format (originally developed for the XMol program by the Minnesota Supercomputer Center) is a widespread plain text format for encoding Cartesian coordinates of molecules, with no formal specification. It contains a header line specifying the number of atoms na, a comment line, and na lines containing element type and atomic coordinates, one atom per line. We have extended this format as indicated in Table 2. Now, the comment line is used to store all scalar properties, Mulliken charges are added as a fifth column. Harmonic vibrational frequencies, SMILES and InChI are appended as respective additional lines.

Table 2 XYZ-like file format for molecular structure and properties.
Full size table

Properties

All molecular geometries were relaxed, and properties calculated, at the DFT/B3LYP/6-31G(2df,p) level of theory. The list of properties of the 134k molecules is summarized in Table 3. For a subset of 6,095 isomers of C7H10O2, energetics (properties 12–16) were additionally calculated at the G4MP2 level of theory. For a validation set of 100 randomly drawn molecules from the 133,885 GDB-9 set, enthalpies of atomization were calculated at the DFT/B3LYP/6-31G(2df,p), G4MP2, G4 and CBS-QB3 levels of theory.

Table 3 Calculated properties.
Full size table

Technical Validation

Validation of geometry consistency

To validate the consistency of the relaxed B3LYP geometries, we have used them to generate the corresponding InChI20 strings with Corina and Open Babel (Version 2.3.0 2011)21. InChI corresponds to ‘IUPAC International Chemical Identifier’. The resulting strings have been compared to the InChI strings obtained from the initial Cartesian coordinates (generated by Corina using the original GDB-17 SMILES strings). See Fig. 2 for a flow-chart of this consistency check. Out of the 134k molecules, 3,054 molecules did not pass this test. This is due to the fact that SMILES and InChI representations are not unique because transformation of Cartesian coordinates to string based chemical identifiers is prone to implementation specific artifacts. For molecules with same topology, small differences in interatomic distances, bond and dihedral angles can lead to different molecular graphs encoded by the string. To quantify this artifact, the distribution of Coulomb-matrix distances13, DIJ [Hartree], using the Manhattan or L1 metric, between Corina generated Cartesian coordinates (see Fig. 2) and B3LYP coordinates is on display in Fig. 3 for all the 3,054 molecules.

Figure 2
figure 2

Schematic flow chart used for geometry consistency check.

Full size image
Figure 3: Histogram of Coulomb-matrix distances.
figure 3

For 3,054 molecules which failed the consistency test shown in Fig. 2 Coulomb-matrix distances, DIJ in Ha, between B3LYP and Corina geometries are shown.

Full size image

Consider, for example, molecule indexed 58 in the 134k set, which is among the 3,054 molecules for which the consistency check failed. Its original GDB-17 SMILES corresponds to NC(=N)C#N, and Cartesian coordinates can be generated using Corina. When feeding back these coordinates to Open Babel to perform the inverse task of reproducing the initial SMILES string, [NH]C(=[NH2])C#N is obtained instead. By contrast, performing first a geometry relaxation of the Corina generated Cartesian coordinates using PM7 followed by B3LYP, and only then parsing through Open Babel, recovers the original SMILES string. In this case, DIJ using the Coulomb-matrices of Corina and B3LYP geometries is rather small (~27 Ha) when compared to distances between geometries of other failed molecules, see distance distribution plot in Fig. 3. By contrast, failed molecules with large DIJ between Corina and B3LYP geometries, see Fig. 3, correspond to molecules for which the B3LYP relaxation induces dramatic distortion with significant alteration and rearrangement of covalent bonding patterns. As a result, parsing back these geometries with Open Babel yields different SMILES strings. Note that all the 6,095 constitutional isomers of C7H10O2 for which G4MP2 calculations have been performed, passed this geometry consistency check, shown as a flow-chart in Fig. 2.

Validation of quantum chemistry results

All 134k molecules have been modeled using B3LYP/6-31G(2df,p) based DFT. Previously, B3LYP has been validated for several subsets, containing up to a few hundred small molecules. These benchmarks are of limited use since they are not necessarily sufficiently representative for gauging B3LYP’s performance in general. In the case of DFT’s systematic errors this issue is particularly pertinent22. Experimental data assembled in the NIST database is very sparse by comparison to our 134k organic molecules made up of CHONF atoms. Consequently, we have performed additional benchmark calculations for a subset of 100 randomly selected molecules using high level theories G4MP212, G423, and CBS-QB324,25.

The predictive power of the G4MP2 method is widely considered to be on par with experimental uncertainties. For example, comparison to the G3/05 test set26,27 with 454 experimental energies (including enthalpies of formation, ionization potentials, electron affinities, proton affinities, and hydrogen bond energies) of small molecules yields MAE, and RMSE of 1.0, and 1.5 kcal/mol, respectively12. For the same properties and molecules, the slightly more accurate, and considerably more expensive method G412,23 yields errors of MAE=0.8 kcal/mol, RMSE=1.2 kcal/mol. G4MP2 has been shown to deviate only by 1.4 kcal/mol from 261 bond dissociation enthalpies computed with the highly accurate W1w composite procedure28,29 for the BDE261 data set28. Consequently, we believe these calculations to be sufficiently suitable to validate the quality of the B3LYP energetics predictions. Various resulting deviations are summarized in Table 1. For the 100 molecules, the mean absolute error of B3LYP heats of atomization amounts to no more than 5 kcal/mol.

Additional information

How to cite this article: Ramakrishnan, R. et al. Quantum chemistry structures and properties of 134 kilo molecules. Sci. Data 1:140022 doi: 10.1038/sdata.2014.22 (2014).