The first aim was to identify the DFT method that yields the most accurate electron densities, and to employ it for computations with a larger basis set. Therefore, the iso-density surfaces computed using different DFT methods were compared with those from CCSD(T) computations for the def2-TZVPD basis set. By comparing all DFT iso-density surfaces for all studied cut-off densities (from 0.0008 to 0.0025 a.u., with 0.0001 a.u. intervals) with the respective CCSD(T) iso-density surfaces, we found the best agreement for the double-hybrid DSD-PBEP86 functional showing only 0.08% average absolute deviation (Table 1).Table 1 Agreement of DFT iso-density surfaces with CCSD(T) results calculated for all studied cut-off densities using the def2-TZVPD basis set. The mean unsigned percentage error (MUPE) and Pearson correlation coefficient (R) are givenAccordingly, we selected this double-hybrid functional for computations of iso-density surfaces with the larger quadruple-zeta basis set and refer to these DSD-PBEP86/def2-QZVPD calculations in the rest of this study for investigating iso-density surfaces. As a further verification, we also computed iso-density surfaces at the CCSD(T)/def2-QZVPD level of theory, which due to the computational demands could be completed for the lowest energy conformers of a subset of 27 compounds. For them, a comparison of iso-density surfaces between DSD-PBEP86 and CCSD(T) methods, both with the def2-QZVPD basis set, yielded a MUPE of 0.16% and a correlation coefficient of 0.99999, underlining the robustness of the double-hybrid DFT results.Comparing the theoretically calculated iso-density surfaces (at DSD-PBEP86/def2-QZVPD level of theory) with experimentally derived TE surfaces as reference values for the 104 studied molecules, the best agreement was observed for the cut-off density of 0.0016 a.u., which yielded MUPE and Pearson correlation coefficient of 1.59% and 0.995, respectively. Details of the computed surfaces are provided in the supplementary material.To further investigate the robustness of this suggested value of the cut-off density, we also evaluated the optimum cut-off densities for 10 datasets, each containing 50 randomly selected compounds from the benchmark set. For all cases, we found the value of 0.0016 a.u. to yield the best agreement between TE and iso-density surfaces.Furthermore, for a larger dataset containing 184 additional compounds with higher uncertainty in the experimental phase-change data, the best agreement between TE and iso-density surfaces was observed for electron density cut-off of 0.0015 a.u. with a MUPE of 6.5%. This result was close to the MUPE of 6.6% obtained using a cut-off density of 0.0016 a.u. (Pearson correlation coefficient is 0.94 in both cases). Despite the larger experimental uncertainties for this dataset, the agreement between the found optimum cut-off densities for the benchmark set and the test set is very comparable. We thus interpret these results as an additional support of the robustness of the suggested cut-off density of 0.0016 a.u.Figure 1 shows the agreement of the molecular iso-density surfaces with the TE surfaces as a function of the density cut-off. Density cut-offs in the range 0.0015 to 0.0019 a.u. provide the most accurate results, with an optimum at 0.0016 a.u. The cut-off densities of 0.002 a.u. proposed by Bader and co-workers36 and 0.001 a.u. by Boyd37 yield MUPEs of 3.17% and 6.98%, respectively, showing the impact of the cut-off density on the computed surfaces and their agreement with the experimental values.Fig. 1: Agreement of iso-density surfaces from quantum chemical calculations with experimental results.The mean unsigned percentage error (MUPE) is plotted as a function of the electron iso-density cut-off. Source data are provided as a Source Data file.Figure 2 shows correlation plots of the molecular surfaces obtained via the two different approaches, experimentally derived (TE) and iso-density surfaces from DSD-PBEP86/def2-QZVPD computations. Different iso-density surfaces are shown, with cut-off densities of 0.0016 a.u. suggested by us (upper panel), 0.002 a.u. as suggested by Bader (middle panel), and 0.001 a.u. as suggested by Boyd (lower panel). The cut-off density of 0.0016 a.u. provides the best agreement, while the two other cut-off densities applied yield slightly larger deviations. The close to perfect correlation (R = 0.995) for all density cut-offs applied underlines the close correspondence between the TE and iso-density surfaces. Considering that the employed experimental TE approach for quantifying molecular surfaces is based on thermodynamic phase-change data while the theoretical results are obtained from the computed electron densities and thus in a completely different way, this remarkable agreement is interpreted as a strong mutual validation of both approaches.Fig. 2: Comparison of iso-density surfaces (in Å2) from quantum chemical calculations with experimental surfaces.a, b Cut-off density 0.0016 a.u. c, d 0.002 a.u. e, f 0.001 a.u. The residuals are plotted at the right. Source data are provided as a Source Data file.Further analysis of the results demonstrates the importance of conformer sampling and Boltzmann averaging. When considering only the single lowest energy conformer of each molecule, the best results were obtained again for the cut-off density of 0.0016 a.u., but the agreement to the TE surfaces (MUPE of 1.87%) was lower than those obtained via conformer sampling (MUPE of 1.59%). The importance of conformer sampling becomes more evident when studying larger and more flexible molecules, which have a greater number of low-energy conformers. Accordingly, in the studied set of molecules, the three most pronounced improvements by conformer sampling were observed for 2-methyl-octane, heptane, and 1-hexanethiol, where the agreement with TE surfaces for single lowest energy conformers in terms of MUPE were reduced via conformer sampling (Boltzmann averaging) from of 2.30, 4.21, and 2.9% to 0.44, 2.12, and 0.19%, respectively. An illustration of the importance of Boltzmann averaging and the variation of iso-density surfaces for two low-energy conformers is depicted in Fig. 3 for 1-pentanethiol as an example.Fig. 3: Iso-density surfaces of two low-energy conformers of 1-pentanethiol contoured at cut-off density of 0.0016 a.u.The surface areas of 159.8 Å2 (left) and 149.2 Å2 (right) are significantly different from the experimentally determined surface (156.5 Å2). The Boltzmann averaged iso-density surface for multiple conformers is 157.2 Å2, in agreement with the experimental value.One of the important applications of the proposed method is the evaluation of atomic radii of elements. In Table 2, we present the atomic radii of a number of elements and their comparison with experimental estimations. Details of the computations, including employed phase-change data are provided in supplementary material.Table 2 Comparison of different estimations of atomic radii of selected elements (in Å). Experimental estimations of atomic radii for noble gas atoms are taken from Ref. 39. and for other elements from Ref. 28Table 2 shows that for noble gases (for which the uncertainty in experimental estimations are significantly lower, as discussed above), there is an excellent agreement between the atomic radii estimated via TE surfaces, iso-density surfaces with cut-off density of 0.0016 a.u., and reference values from Ref. 39. Also, for N and F elements, a close agreement between the radii determined via iso-density surfaces and the experimental estimation is observable. For the radii estimated with iso-density cut-off of 0.001 a.u., the results closely match the values reported by Rahm et al. calculated for the same cut-off density at PBE0 level of theory38. However, the computed radii with cut-off density of 0.001 a.u. are significantly larger than those predicted by the other studied methods, highlighting the importance of the applied density cut-off. Noteworthy, the estimations of atomic radii for N, O, and F elements via iso-density surfaces are obtained from the diatomic molecule. In another study40, we recently proposed a rigorous approach to estimate radii of open-shell atoms in their isolated states on the basis of the Tkatchenko-Scheffler method14 relating radii of atoms in molecules to the radii of free atoms via:$${R}_{A}={R}_{A,{free}}{\left(\frac{{V}_{A}^{{eff}}}{{V}_{A}^{{free}}}\right)}^{1/3},$$
(1)
$$\frac{{V}_{A}^{{eff}}}{{V}_{A}^{{free}}}=\frac{\int {r}^{3} {w}_{A}({\bf{r}}) n({\bf{r}}) \, {d}^{3}{\bf{r}}}{\int {r}^{3}{n}_{A}^{{free}}({\bf{r}}) \, {d}^{3}{\bf{r}}},$$
(2)
where \({w}_{A}\) is the atomic weight for partitioning the molecular space into atomic sub-spaces, and \({R}_{{A},{free}}\) and \({n}_{A}^{{free}}\) are the radius and the electron density of atom A in the free (isolated) state.Using iso-density surfaces with a cut-off density of 0.0016 a.u., we optimized the radii of free atoms in a way to obtain the best match between solvent-excluded surfaces (SES) constructed from the radii of atoms in molecules based on Eq. (1) and iso-density surfaces for a set of 1235 molecules40. The molecular surfaces predicted on the basis of these two different approaches closely match, with MUPE of 0.75% and Pearson correlation coefficient of 0.9996. Also, the optimized radii of N and O atoms, found to be 1.65 Å and 1.49 Å, respectively, perfectly agree with the experimental estimations (1.66 Å and 1.50 Å, Table 2). The consistency of the two alternative computational methods for evaluating molecular surfaces and the agreement of the optimized radii with the experimental estimations strongly reinforce the robustness of the proposed approaches.
Electron iso-density surfaces provide a thermodynamically consistent representation of atomic and molecular surfaces
