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))).

Background & Summary

Consistent and curated datasets have facilitated progress in the natural sciences. High-quality reference data sets were, for example, essential in the development of accurate computational methodology, in particular in quantum chemistry. With the rise of machine learning, datasets have increased in size and have transformed from reference status to a primary source of data for predictions[1-7] and discovery[8-12].

In this article we present a new dataset for molecular spectroscopy applications. Spectroscopy is ubiquitous in science as one of the primary ways of determining a material’s or molecule’s properties. However, publicly available spectroscopic datasets for technologically relevant molecules are rare. Examples include a dataset of chemical shifts for structures taken from the CSD[13,14], the Harvard Clean Energy Project[15] as well as the QM8[16,17] and QM9[18] datasets. The QM8 database offers optical spectra computed with time-dependent density functional theory (TDDFT) for 22 k organic molecules, while QM9, widely known as one of the standard benchmark sets for machine learning in chemistry, provides a variety of properties for 134 k organic molecules computed with density functional theory (DFT)[19,20], including energy levels for the highest occupied and the lowest unoccupied molecular orbitals (HOMO and LUMO, respectively). Although QM8 and QM9 are of unprecedented size compared to previous, common benchmark sets in quantum chemistry of several hundred to thousands of molecules, they still contain only small molecules with restricted elemental diversity (H, C, N, O and F) and with simple bonding patterns[6]. They lack larger, more complex molecules with, e.g., extended heteroaromatic backbones and attached functional groups, as commonly targeted in organic synthesis[21,22] and applied in (opto-)electronic[23-26] or pharmaceutical research[22,27,28].

We have based the spectroscopic dataset presented in this article on a diverse collection of 64,725 organic crystals that were extracted from the Cambridge Structural Database (CSD)[29] by Schober et al.[30,31]. This 64 k dataset of experimental crystal structures gathered from a variety of application areas was originally compiled to optimize the charge carrier mobility for applications in organic electronics. For our OE62 dataset, we used 61,489 unique organic molecular structures, extracted from the respective organic crystals. All extracted geometries were then relaxed in the gas phase with density-functional theory (DFT).

The molecules in OE62 cover a considerable part of chemical space, as illustrated in Fig. 1. The dataset contains molecules with up to 174 (or 92 non-hydrogen) atoms and a diverse composition of 16 different elements. A large number of different scaffolds and functional groups are included, representing a multifaceted sample of the design space available in organic chemistry[6,30,32,33].

Fig. 1

Chemical space spanned by OE62. (a) Molecular size distributions (including hydrogen atoms) for the OE62 dataset and its 31 k and 5 k subsets. (b) Distribution of the 16 different element types in the datasets. (c) Typical structures found in the 62 k dataset, with chemical diversity arising from a rich combinatorial space of scaffold-functional group pairings: The dataset contains aliphatic molecules, as well as molecules with conjugated and complex aromatic backbones and diverse functional groups of technological relevance. The refcode_csd identifiers of depicted molecules are (from left to right): ZZTVO01, VOCMIK, FATVEC, WASVAN, BIDLUW, KETZAL, EHORAU.

To go into more detail, all molecules in OE62 are fully relaxed at the Perdew-Burker-Ernzerhof (PBE)[34] level of DFT including Tkatchenko-Scheffler van der Waals (TS-vdW) corrections[35]. For these equilibrium structures, we then report molecular orbital energies at the PBE and PBE hybrid (PBE0)[36,37] level, in the following referring to this part as 62 k set. Partial charges and total energies for DFT-calculations are also included. In two subsets, randomly drawn to span more than half (31 k) and more than 5000 (5 k) of the molecular structures, we provide additional computational results: the influence of solvation – in this case implicit water – on the energy levels is addressed on the PBE0 level for a subset of 30,876 molecules. For the second subset of 5,239 molecules, we computed the quasi-particle energies with many-body perturbation theory in the G0W0 approximation[19,38,39] and extrapolated to the complete basis set (CBS) limit. Figure 2 gives a schematic overview of the dataset nesting in OE62 while Table 1 lists computational settings and computed properties. Figure 3(a,b) illustrate the HOMO level and solvation free energy distributions of the 5 k subset.

Fig. 2

Schematic overview of the three datasets and the applied computational methods. The 31 k set includes all structures from the 5 k set and the 62 k all structures from the 31 k and 5 k sets.

Table 1

Overview of the data (sub)sets in OE62: Applied computational method, resulting molecular properties and DOI-based references to the input and output files of corresponding calculations deposited in the NOMAD repository.

Set	Method	Computed properties	Access to data records on NOMAD
62 k	DFT PBE + vdW (vacuum) Tier2 basis set, tight settings	• relaxed geometry	[61–67]
		• occupied & unoccupied MO energies
		• total energy
		• Hirshfeld charges
62 k	DFT PBE0 (vacuum) Tier2 basis set, tight settings	• geometry fixed at the PBE + vdW level	[68]
		• occupied & unoccupied MO energies
		• total energy
		• Hirshfeld charges
31 k	DFT PBE0 (water) Tier2 basis set, tight settings, MPE implicit solvation	• geometry fixed at the PBE + vdW level	[69]
		• occupied & unoccupied MO energies
		• total energy
		• Hirshfeld charges
5 k	DFT PBE0 (vacuum) def2-TZVP & def2-QZVP basis sets (see text), tight settings	• geometry fixed at the PBE + vdW level	[70]
		• occupied & unoccupied MO energies
		• total energy
5 k	G0W0@PBE0 (vacuum) def2-TZVP & def2-QZVP basis sets (see text), tight settings	• geomety fixed at the PBE + vdW level	[71]
		• occupied & unoccupied MO energies
		• CBS energies of occupied & unoccu pied MOs

Fig. 3

The GW5000 subset compared to the other (sub)sets in OE62. Panel (a) shows distributions of HOMO energies from G0W0@PBE0 (vacuum), PBE + vdW, PBE0 (vacuum) and PBE0 (water) computations. Panel (b) shows the distribution of solvation free energies Δ. In (a and b), distribution medians are marked by dotted lines. Panel (c) depicts a correlation plot for the approximately linear relationship between the G0W0@PBE0 CBS quasiparticle energies and the DFT HOMO energies (PBE and PBE0 in vacuum).

We refer to the 5 k subset of G0W0 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 G0W0 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 G0W0@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 G0W0 quasiparticle energies from the computationally cheaper PBE0 method without having to perform G0W0 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 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 G0W0 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-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-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.

62 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 fmax < 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 fmax < 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).

62 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).

31 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 of the solute’s electronic density. Outside of this cavity the dielectric constant of water was applied[52]. The density isovalue 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 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 lmax,R and lmax,O to 14 and 8, respectively, and the degree of overdetermination dod 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: G0W0@PBE0 (vacuum)

For the 5 k subset, the relaxed PBE + vdW structures in vacuum were used as input for the G0W0[19,39,53] calculations, using the FHI-aims G0W0 implementation based on the analytic continuation[44]. The PBE0 hybrid functional was used for the underlying DFT calculation (G0W0@PBE0) in combination with the atomic ZORA approximation.

In these G0W0 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 G0W0 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 G0W0 self-energy elements were calculated for a set of imaginary frequencies 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 was performed using a modified Gauss-Legendre grid[44] with 200 grid points. The same grid was employed for the set of frequencies , 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:

The dataset and related files can be freely downloaded from the media repository of the Technical University of Munich (mediaTUM) under 10.14459/2019mp1507656[59]. The dataset is provided as JSON output data of Pandas[60] DataFrames. Within Python, these dataframes allow structured access to data in a tabular format, where each molecule is stored in a row of the dataframe, while the data is organized in columns. The content of the dataframe is summarized and explained in 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/).

Table 2

Dataframe structure of all three dataframes df_62 k, df_31 k and df_5 k.

No.	Column name	Unit	Method	Dataframes	Description
1	refcode_csd	—	—	62 k, 31 k, 5 k	CSD reference code, unique identifier for the crystal from which the molecule was extracted
2	canonical_smiles	—	Open Babel	62 k, 31 k, 5 k	Molecular string representations derived from DFT PBE + vdW relaxed geometries.
3	inchi	—	Open Babel	62 k, 31 k, 5 k
4	number_of_atoms	—	—	62 k, 31 k, 5 k	Number of atoms in the molecule
5	xyz_pbe_relaxed	Å	PBE + vdW (vacuum)	62 k, 31 k, 5 k	String in XYZ-file format of DFT PBE + vdW relaxed geometry. Line 1 contains the number of atoms. Line 2 is empty. The remaining lines contain atomic type and coordinate (x, y, z).
6	energies_occ_pbe	eV	PBE + vdW (vacuum)	62 k, 31 k, 5 k	List of eigenvalues of occupied molecular Kohn-Sham orbitals. Given in ascending order, the last value is the HOMO energy.
7	energies_occ_pbe0_vac_tier2	eV	PBE0 (vacuum)	62 k, 31 k, 5 k
8	energies_occ_pbe0_water	eV	PBE0 (water)	31 k, 5 k
9	energies_occ_pbe0_vac_tzvp	eV	PBE0 (vacuum)	5 k
10	energies_occ_pbe0_vac_qzvp	eV	PBE0 (vacuum)	5 k
11	energies_occ_gw_tzvp	eV	G0W0@PBE0 (vacuum)	5 k
12	energies_occ_gw_qzvp	eV	G0W0@PBE0 (vacuum)	5 k
13	cbs_occ_gw	eV	G0W0@PBE0 (vacuum)	5 k	List of CBS energies of occupied states computed from G0W0@PBE0 TZVP and QZVP energies from 10 and 11. Same order as lists described above.
14	energies_unocc_pbe	eV	PBE + vdW (vacuum)	62 k, 31 k, 5 k	List of eigenvalues of virtual (unoccupied) molecular Kohn-Sham orbitals. Given in ascending order, the first value is the LUMO energy. Only virtual states below the vacuum level (i.e. with negative eigenvalue) are listed. If the LUMO energy is positive, only the LUMO energy is listed. If 20 has more negative eigenvalues than 19, we also include positive eigenvalues in 19 so that both lists in 19 and 20 have equal length.
15	energies_unocc_pbe0_vac_tier2	eV	PBE0 (vacuum)	62 k, 31 k, 5 k
16	energies_unocc_pbe0_water	eV	PBE0 (water)	31 k, 5 k
17	energies_unocc_pbe0_vac_tzvp	eV	PBE0 (vacuum)	5 k
18	energies_unocc_pbe0_vac_qzvp	eV	PBE0 (vacuum)	5 k
19	energies_unocc_gw_tzvp	eV	G0W0@PBE0 (vacuum)	5 k
20	energies_unocc_gw_qzvp	eV	G0W0@PBE0 (vacuum)	5 k
21	cbs_unocc_gw	eV	G0W0@PBE0 (vacuum)	5 k	List of CBS energies of unoccupied states computed from G0W0@PBE0 TZVP and QZVP energies from 19 and 20. Same order as lists described above.
22	total_energy_pbe	eV	PBE + vdW (vacuum)	62 k, 31 k, 5 k	Total energy of the DFT calculations. Note, for consistency with 22, 23 and 24 also include the vdW contribution to the total energy. 25 and 26 do not include it.
23	total_energy_pbe0_vac_tier2	eV	PBE0 (vacuum)	62 k, 31 k, 5 k
24	total_energy_pbe0_water	eV	PBE0 (water)	31 k, 5 k
25	total_energy_pbe0_vac_tzvp	eV	PBE0 (vacuum)	5 k
26	total_energy_pbe0_vac_qzvp	eV	PBE0 (vacuum)	5 k
27	hirshfeld_pbe	qe	PBE + vdW (vacuum)	62 k, 31 k, 5 k	List of Hirshfeld partial charges on atoms. Same order as atoms in xyz_pbe_relaxed.
28	hirshfeld_pbe0_vac_tier2	qe	PBE0 (vacuum)	62 k, 31 k, 5 k
29	hirshfeld_pbe0_water	qe	PBE0 (water)	31 k, 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.

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, an empty comment line and n 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. G0W0@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 G0W0@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.

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.

Fig. 4

Coulomb matrix distances between initial crystal geometries and PBE + vdW relaxed geometries. Panel (a) shows the distribution of Coulomb matrix distances for all 62 k molecules and panel (b) shows the distribution of Coulomb matrix distances for the 284 cases that did not pass the consistency check. Two example molecules are shown in (a) for short and large distances between Coulomb matrices (the refcode_csd identifiers are CILWUP (1) and ODAHUW (2)). In (b), 2D structures of three example molecules that failed the consistency check are shown (DAZIND (3), YOMDUA (4) and FODBAC (5)).

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 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.

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).

Fig. 5

Accuracy assessment of HOMO- and atomization energies computed at the PBE0 (vacuum) DFT level of theory. (a) Four example molecules and their refcode_csd identifiers. (b) For the example molecules, the HOMO energy convergence of the Tier1 and Tier2 basis sets is compared against the Tier4 basis set provided with FHI-aims, always employing tight integration settings. (c) Difference in HOMO-energy between the Tier2 (T2) and QZVP basis sets for all molecules of the 5 k set. The distribution-median is given by a dotted line, located at −0.008 eV. (d) Same as (b), but for atomization energies Ef.

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 G0W0 quasiparticle energies

Figure 6(a) shows the convergence of the G0W0@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 G0W0 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 G0W0 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.

Fig. 6

Accuracy assessment of G0W0 quasiparticle energies. (a) Convergence of the HOMO G0W0 energies with respect to the inverse of the number of basis functions NBF for the four example molecules shown in Fig. 5. Dashed lines represent linear straight line fits using the def2-QZVP and def2-TZVP points. The intersection of the straight line with the ordinate gives an estimate for the complete basis set limit (CBS) as indicated for BMLTAA. (b) Deviation of the HOMO G0W0 energies from the CBS limit for the 5 k subset. Median values of the distributions are indicated by black dashed lines. (c) Percentage of states with negative slope of the CBS fit. (d) Average G0W0@PBE0 quasiparticle spectrum, where each energy state was artificially broadened by a Gaussian distribution.

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 G0W0 self-energy 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.

G0W0 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 G0W0 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 G0W0 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.

Measurement(s)	organic molecule
Technology Type(s)	digital curation • spectroscopy
Factor Type(s)	computational method