Atomic structures and orbital energies of 61,489 crystal-forming organic molecules

Data science and machine learning in materials science require large datasets of technologically relevant molecules or materials. Currently, publicly available molecular datasets with realistic molecular geometries and spectral properties are rare. We here supply a diverse benchmark spectroscopy dataset of 61,489 molecules extracted from organic crystals in the Cambridge Structural Database (CSD), denoted OE62. Molecular equilibrium geometries are reported at the Perdew-Burke-Ernzerhof (PBE) level of density functional theory (DFT) including van der Waals corrections for all 62 k molecules. For these geometries, OE62 supplies total energies and orbital eigenvalues at the PBE and the PBE hybrid (PBE0) functional level of DFT for all 62 k molecules in vacuum as well as at the PBE0 level for a subset of 30,876 molecules in (implicit) water. For 5,239 molecules in vacuum, the dataset provides quasiparticle energies computed with many-body perturbation theory in the G0W0 approximation with a PBE0 starting point (denoted GW5000 in analogy to the GW100 benchmark set (M. van Setten et al. J. Chem. Theory Comput. 12, 5076 (2016))).

lists computational settings and computed properties. Figure 3(a,b) illustrate the HOMO level and solvation free energy distributions of the 5 k subset.
We refer to the 5 k subset of G 0 W 0 quasiparticle energies as GW5000 in analogy to the GW100 benchmark set 40 . GW100 was a landmark dataset of 100 atoms and molecules that for the first time demonstrated the high numerical accuracy of the computationally costly G 0 W 0 approach. GW100 quickly became the standard reference for GW code development and validation. The GW5000 subset in OE62 is of the same high numeric quality as GW100, but extends the set of reference molecules by a factor of 50. To illustrate the value of multi-level computational results we present a first, preliminary finding in Fig. 3. Panel c) shows the correlation between the G 0 W 0 @ PBE0 quasiparticle HOMO energies and the DFT HOMO eigenvalues for the GW5000 subset. The correlation is to first approximation linear with PBE0 having a lower variance than PBE. This linear relation (slope of 1.195 and intercept of −0.492 for PBE0) could now be used to predict G 0 W 0 quasiparticle energies from the computationally cheaper PBE0 method without having to perform G 0 W 0 calculations. Applying this linear correction to the PBE0 results yields quasiparticle energy predictions with a root mean square error (RMSE) of only 0.17 eV to the respective GW5000 values.
Given the high-quality computational results from different levels of theory, the (subs)sets included in OE62 can be used to develop, train and evaluate machine learning algorithms, facilitating the search and discovery of diverse molecular structures with improved properties. In the following, we first describe the procedure used to . In (a and b), distribution medians are marked by dotted lines. Panel (c) depicts a correlation plot for the approximately linear relationship between the G 0 W 0 @PBE0 CBS quasiparticle energies and the DFT HOMO energies (PBE and PBE0 in vacuum).
www.nature.com/scientificdata www.nature.com/scientificdata/ compute molecular structures and properties, followed by a full description of the dataset format and content as well as by a validation of our DFT and G 0 W 0 results. OE62 is freely available as a download from the Technical University of Munich. The input and output files of all calculations performed for OE62 can be downloaded from the Novel Materials Discovery (NOMAD) laboratory (https://nomad-repository.eu).

Methods
All crystal structures collected from the CSD for the 64 k dataset are mono-molecular, i.e. they contain only a single type of molecule per unit-cell. A single molecular structure (conformer) from each crystal was extracted by a custom Python code 30,31 . This 64 k dataset of molecular structures provides the starting point for the dataset published here. A fraction of the crystals contained in the CSD have polymorphic forms or were added multiple times, coming e.g. from different experimental sources. Although they occur in different crystalline entries in the 64 k dataset, the same molecular structure could enter our molecular database multiple times. First, the SMILES identifiers were computed for the 64 k dataset 30,31 from a combination of Open Babel 41 (www.openbabel.org) and RDKit (www.rdkit.org) 42 . We subsequently excluded all extracted molecules whose non-isomeric, canonical SMILES identifier occurred multiple times, keeping only one case each. Further, molecules with an odd number of electrons were removed. After these filtering steps 61,539 molecules remained.
We relaxed the geometries of all molecules at the PBE + vdW level of theory, as implemented in the FHI-aims all-electron code [43][44][45] . We chose the PBE + vdW functional for three reasons: 1) It is an all-purpose functional with a favorable accuracy/computational cost ratio that is implemented in all the major electronic structure codes. 2) We would like to stay consistent with previous work 6,46 , in which PBE + vdW was also used for molecular structures optimization of large molecular data sets. 3) While there might be more accuracte semi-local functionals than PBE 47 , the addition of vdW corrections makes PBE + vdW appropriate for organic compounds. For organic crystals, for which highly accurate, low-temperature experimental geometries are available, PBE + vdW yields excellent agreement with typical root-mean-squared deviations of only 0.005-0.01 Å per atom [48][49][50] .
Given that slightly differing bond assignments in the newly obtained low-energy geometries might change some of the molecular identifiers, we generated new InChI 51 ('IUPAC International Identifier') and canonical SMILES identifiers using Open Babel (Version 2.4.1 2016), and report these in our dataset. We then checked these representations for duplicates and concurrently removed them. In addition, 6 molecules were removed for which geometry optimization or single point calculations had failed. In total, 61,489 unique molecules remained, which form the basis of the OE62 set.
From the OE62 set we generated two subsets: For the 31 k subset we randomly picked 30,876 molecules. The same was done for the 5 k set by randomly picking 5,239 molecules from the 31 k subset with the additional constraint that the largest molecule should not exceed 100 atoms. The size and element distributions of all three sets are shown in Fig. 1.
In the following we explain the data and additional subsets we created and provide the computational settings. All settings are also listed in Table 1.

k set: DFT PBE + vdW (vacuum).
We pre-relaxed all molecular geometries at the PBE level of theory.
For structure relaxation, we used the trust radius enhanced variant of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm as implemented in FHI-aims with a maximum atomic residual force criterion of f max < 0.01 eV Å −1 . The electronic wave functions were expanded in a Tier1 basis set at light integration settings 43 . Since our database only contains closed-shell molecules, we performed spin-restricted DFT calculations. Dispersive forces were included in the geometry relaxations using the Tkatchenko-Scheffler (TS) 35 method, while relativistic effects were treated on the level of the atomic zero-order regular approximation (atomic ZORA) 43 . The DFT self-consistency cycle was treated as converged when changes of total energy, sum of eigenvalues and charge density were found below 10 −6 eV, 10 −3 eV and 10 −5 e Å −3 , respectively. Starting from these pre-relaxed structures, we obtained the final geometries by performing a new relaxation with Tier2 basis sets, tight integration settings and a convergence criterion of f max < 0.001 eV Å −1 . The eigenvalues of the molecular states are then stored in our dataset alongside the molecular geometries. We refer to this part of the dataset as PBE + vdW (vacuum).

k set: DFT PBE0 (vacuum).
Using the relaxed geometries obtained at the PBE + vdW (vacuum) level of theory, we further carried out single point calculations for all structures using the PBE0 hybrid functional. Computational settings as described before were used, employing again the Tier2 basis set with a tight integration grid. Note, that tabulated total energies obtained at this level also include the vdW contribution computed through the TS method, while "vdW" was dropped from the name to emphasize the single point character of these computations. We correspondingly refer to this set as PBE0 (vacuum).

k subset: DFT PBE0 (water).
To study the influence of solvation-here by water-on the PBE0 results, we performed calculations using the Multipole Expansion (MPE) implicit solvation method as implemented in FHI-aims 52 for the 31 k susbset. The MPE method facilitates an efficient treatment of the solvation effects on a solute, by using a continuum model of the solvent around it. In detail, the solute molecule is placed within a cavity with the dielectric permitivity of vacuum. The position of the cavity surface is determined by an iso-value iso ρ of the solute's electronic density. Outside of this cavity the dielectric constant of water ε = .
78 36 b was applied 52 . The density isovalue iso ρ as well as the α and β parameters for non-electrostatic contributions to the solvation free energy were taken from the published SPANC parameter-set 52 .
In the MPE method, the solvation cavity is discretized using a large number of points homogeneously distributed at the density iso-surface. Sampling of these points was achieved using an inexpensive pseudo-dynamical optimisation, allowing up to 1000 optimisation steps and removing the worst 0.1% of walkers at each neighbor-list (2020) 7:58 | https://doi.org/10.1038/s41597-020-0385-y www.nature.com/scientificdata www.nature.com/scientificdata/ update step 52 , to account for the more complex molecules included in the 62 k dataset. To obtain highly converged eigenvalues, we increased the reaction field-and polarization potential expansion orders l max,R and l max,O to 14 and 8, respectively, and the degree of overdetermination d od to 16, keeping all other parameters at their default values 52 . Note that the molecular geometries were not further relaxed in the presence of the water solvent. We kept the structures fixed at the PBE + vdW level. Tabulated total energies again include the vdW contribution obtained by the TS method. The resulting data is referred to as PBE0 (water). 5 k subset: G 0 W 0 @PBE0 (vacuum). For the 5 k subset, the relaxed PBE + vdW structures in vacuum were used as input for the G 0 W 0 19,39,53 calculations, using the FHI-aims G 0 W 0 implementation based on the analytic continuation 44 . The PBE0 hybrid functional was used for the underlying DFT calculation (G 0 W 0 @PBE0) in combination with the atomic ZORA approximation.
In these G 0 W 0 and PBE0 calculations, we employed the def2 triple-zeta valence plus polarization (def2-TZVP) and the def2 quadruple-zeta valence plus polarization (def2-QZVP) basis sets 54 . The def2-TZVP and def2-QZVP basis sets are contracted Gaussian orbitals, treated numerically to be compliant with the numeric atom-centered orbital (NAO) technology in FHI-aims 43 . They are fully all-electron for all elements and do not contain effective core potentials. The def2 basis sets are available from the EMSL database 55,56 , except for iodine (see Supplementary  Information). Note that a basis set of def2-TZVP quality is not available for I and all def2-TZVP calculations for iodine-containing molecules were correspondingly performed with def2-QZVP for I and with def2-TZVP for all other elements.
Since G 0 W 0 calculations converge slowly with respect to basis set size 39 , we extrapolated the quasiparticle energies to the complete basis set (CBS) limit. Following the procedure for the GW100 benchmark set 40 , the extrapolated values are calculated from the def2-TZVP and def2-QZVP results by a linear regression against the inverse of the total number of basis functions (see Technical Validation).
The G 0 W 0 self-energy elements were calculated for a set of imaginary frequencies ω i { } and then analytically continued to the real frequency axis using a Padé approximant 57 with 16 parameters. The numerical integration along the imaginary frequency axis ω′ i { } was performed using a modified Gauss-Legendre grid 44 with 200 grid points. The same grid was employed for the set of frequencies ω i { }, for which the self-energy is computed. The analytic continuation in combination with the Padé model yields accurate results for valence states 40 , but is not reliable for core and semi-core states 58 . Therefore, we included only occupied states with quasiparticle energies larger than −30 eV in the data set, see also Technical Validation for more details.

Data Records
The curated data for all 61,489 molecules is publicly available from two sources:  Table 2. We also provide a tutorial file, which explains loading, filtering and data extraction from dataframes within Python. On mediaTUM, the dataset is distributed under a Creative Commons licence (https://creativecommons.org/ licenses/by-sa/4.0/). 2. The input and output files of all performed calculations can be downloaded from NOMAD. Due to the size of OE62 we provide an individual DOI for each applied computational method 61-71 . Dataframe format. We provide three dataframes: df_62 k, df_31 k and df_5 k. For each molecule in these dataframes, we provide three identifiers (refcode_csd, canonical_smiles and inchi in columns 1 to 3). In column 5, atomic coordinates of PBE + vdW (vacuum) relaxed structures are stored as a string in a standard XYZ format (xyz_pbe_relaxed): The structure information contains a header line specifying the number of atoms n a , an empty comment line and n a lines containing element type and relaxed atomic coordinates, one atom per line. The structure of all three dataframes is summarized in Table 2.
The following list provides a brief overview over the three dataframes: • Dataframe df_5 k includes 5,239 structures with results for all molecular properties in columns 5 to 29.
• Dataframe df_31 k accommodates 30,876 structures, including all structures from df_5 k. G 0 W 0 @PBE0 results are only available for molecules from its 5 k subset, while respective columns are left blank for the remaining molecules in df_31 k. • Dataframe df_62 k contains all 61,489 structures, including all structures from df_31 k and df_5 k.
PBE0 (water) results are only available for molecules from its 31 k subset, while respective columns are left blank for the remaining molecules in df_62 k. The same applies for G 0 W 0 @PBE0 results for the structures from the 5 k subset. The dataframe is ordered, such that the molecules included in the 5 k subset are included first, while the remaining molecules of 31 k and 62 k subsets follow subsequently. This data structure facilitates the filtering of the dataframe by single lines of code, as shown in the tutorial.
In addition, a spreadsheet file is provided in the distributed archive which contains the total energies of all atomic species of the dataset. They are computed for the respective levels of theory using similar computational settings, so that atomization energies for all molecules can be computed from the available molecular total energies. (2020) 7:58 | https://doi.org/10.1038/s41597-020-0385-y www.nature.com/scientificdata www.nature.com/scientificdata/ Finally, future updated versions of the dataset on mediaTUM will be distributed through the versioned DOI given above. In such cases, updated descriptions will be provided in the distributed archive alongside the dataset.

Technical Validation
Validation of relaxed geometries. To quantify the degree to which relaxation in vacuum changes the geometry of the structures compared to their crystalline form, we computed the distance between the two Coulomb matrices 72,73 of the original crystal geometry and the PBE + vdW relaxed geometry for each of the 62 k molecules. The distribution of these Coulomb matrix distances is shown in Fig. 4(a). Small distances signify small changes and large distances signify large differences. Most molecules exhibit only little changes in geometry during relaxation, where bond lengths are shifted by a small amount, as illustrated for the example of molecule 1. In some rare cases we find significant shifts in geometry caused by the environmental change from intermolecular interactions in the crystal to intramolecular interactions in vacuum, as shown for molecule 2. The crystal-extracted structure is shaped according to intermolecular van der Waals interactions that were present in the crystal. After relaxation, the intramolecular interactions cause a contraction of the molecular structure.
To validate that the chemical integrity of the majority of the 62 k molecules is preserved during the PBE + vdW relaxation, we perform a consistency check similarly to ref. 18 . We generated InChI strings from the relaxed PBE + vdW geometries and compare them to those obtained from the initial crystal-extracted cartesian coordinates. For 284 pairs, the two InChI strings did not match. Such mismatches can, for example, be caused by specifics in the implementation, in which Openbabel assigns different InChI strings to molecules with the same topology, possibly caused by changes in bond lengths, bond angles or dihedral angles. Examples are shown in Fig. 4(b) with molecule 3 exhibiting a small Coulomb matrix distance or molecule 5, which exhibits a large  Table 2. Dataframe structure of all three dataframes df_62 k, df_31 k and df_5 k. Columns 1 to 3 contain molecular identifiers. Columns 5 to 29 contain molecular properties computed at respective level of theory. All mentioned energies are given in eV.
www.nature.com/scientificdata www.nature.com/scientificdata/ Coulomb matrix distance due to stronger relaxation. Here, stereoassignments change in the molecular structure, causing the different InChI-identifiers. Conversely, the mismatch can be also caused by changes in molecular topology during relaxation. This is the case for molecule 4, for which an intramolecular ring-closure takes place. Compared to 3,054 such inconsistencies found during the collection of the 134 k molecules for the QM9 database 18 , the number of 284 found here is considerably small. The reason is most likely that our molecular starting geometries were derived from experimentally observed, well-resolved solid-form conformers.  www.nature.com/scientificdata www.nature.com/scientificdata/ Validation of DFT atomization and orbital energies. For PBE and PBE0 calculations, the Tier2 basis set of FHI-aims typically provides converged results for both the atomization energy as well as for molecular orbital energies 74,75 . The Tier2 basis set has also been used in other large molecular datasets 46,72 . We here illustrate the convergence for four selected cases featured in Fig. 5(a) for PBE0 vacuum calculations at tight settings. As expected, HOMO energies at the Tier2 level are well-converged, here estimated within 0.01 eV around reference values obtained with the largest standard basis set included in FHI-aims (Tier4), see Fig. 5(b). The lower lying orbital energies exhibit a similar convergence behavior (not shown).
A further quality assessment of predicted HOMO-energies comes from the comparison of Tier2 and QZVP basis set results, as contained in the 5 k subset, see Fig. 5(c). We find only a small RMSE of 0.009 eV between the Tier2 and the much larger QZVP basis sets. Figure 5 also shows the convergence of the atomization energy of the four molecules in panel d). Again, at the Tier2 level, convergence to better than 0.1 eV with respect to Tier4 is observed. This is consistent with results found in a previous benchmark study 74 . Validation of G 0 W 0 quasiparticle energies. Figure 6(a) shows the convergence of the G 0 W 0 @PBE0 quasi-particle energies with respect to basis set size and their extrapolation to the CBS limit for the four molecules displayed in Fig. 5(a). In all four cases, the G 0 W 0 energies are not converged even with the largest basis set and CBS extrapolation is required. The slow convergence is typical for the whole 5 k set, as demonstrated in Fig. 6(b), which reports the deviation of the HOMO G 0 W 0 energies computed at the TZVP and QZVP level from the CBS limit for all molecules of the 5 k subset. The distributions displayed in Fig. 6(b) are centered around −0.38 eV (TZVP) and −0.17 eV (QZVP) with a standard deviation of 0.02 eV (TZVP) and 0.01 eV (QZVP) from the median values. Similar results are obtained by including all occupied states above −30 eV in the analysis. In this case, the median value amounts to −0.35 eV for TZVP and −0.15 eV for QZVP. Respective distributions for the deviations of all occupied states from the CBS limit can be found in the Supporting Information.
The quasiparticle energies at the QZVP level are typically lower in energy than the TZVP values, i.e., the straight line determined from the linear extrapolation to the CBS limit has a positive slope, see Fig. 6(a). This empirical observation was already made in the GW100 benchmark study 40 for the HOMO level and we also observed it here in our GW5000 study for the valence states. There is no proof that for a given basis set the slope has to be positive. In fact, for ~4% of the energies level above −30 eV we find negative slopes, as shown in Fig. 6(c). This percentage increases considerably in the semi-core energy region between −50 and −30 eV. Such an increase is indicative of either 1) a failure of the analytic continuation used to continue the G 0 W 0 self-energy www.nature.com/scientificdata www.nature.com/scientificdata/ from imaginary-to the real-frequency axis or 2) the insufficiency of the def2-TZVP basis set to converge the deeper occupied states at the DFT level. Based on our analysis in Fig. 6(c), we therefore include only states with energies larger than −30 eV in the 5 k set. Figure 6(d) confirms that the spectral weight averaged over the whole 5 k subset is located mostly between −30 to −5 eV and thus, not much spectral information is lost by setting the cutoff threshold to −30 eV.
G 0 W 0 calculations were initially run for 5,500 structures randomly drawn from the 31 k set. From these 5,500 molecules, we filtered out molecules for which the analytic continuation of the G 0 W 0 self-energy is inaccurate or breaks down completely. In FHI-aims the quasiparticle equation is solved iteratively to determine the quasiparticle energies. For some molecules, the pole structure of the self-energy gives rise to multiple solutions and the iterative solution does not converge. We excluded all molecules from the dataset for which at least one TZVP or QZVP level did not converge. Moreover, large differences between the TZVP and QZVP quasiparticle energies are an indication of further problems in the G 0 W 0 calculation, since the median difference between TZVP and QZVP is only 0.21 eV (see Fig. 6(b)). We thus excluded molecules for which at least one level exceeded QZVP/ TZVP difference of 0.8 eV. This leaves 5,239 molecules in the 5 k set.

Code availability
All electronic structure data contained in this work was generated with the FHI-aims code [43][44][45] . The code is available for a license fee from https://aimsclub.fhi-berlin.mpg.de/aims_obtaining_simple.php. Parsing of outputs and data collection were performed with custom-made Python scripts, which will be available upon request. Finally, the published archive contains a tutorial detailing how to access the dataset.