Statistical Analysis of Semiclassical Dispersion Corrections
Thomas Weymuth Jonny Proppe Markus Reiher
physics.chem-ph cond-mat.mtrl-sci physics.data-an





Statistical Analysis of Semiclassical Dispersion Corrections





Thomas Weymuth^a,
Jonny Proppe^a,
and Markus Reiher^a,[Corresponding author; e-mail: markus.reiher@phys.chem.ethz.ch]



^a Laboratorium für Physikalische Chemie, ETH Zürich, 

Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland

January 25, 2018









Abstract




Semiclassical dispersion corrections developed by Grimme and coworkers have become indispensable in applications of Kohn–Sham density functional theory.
A deeper understanding of the underlying parameterization might be crucial for well-founded further improvements of this successful approach.
To this end, we present an in-depth assessment of the fit parameters present in semiclassical (D3-type) dispersion corrections by means of a statistically rigorous analysis. 
We find that the choice of the cost function generally has a small effect on the empirical parameters of D3-type dispersion corrections with respect to the reference 
set under consideration.
Only in a few cases, the choice of cost function has a surprisingly large effect on the total dispersion energies.
In particular, the weighting scheme in the cost function can significantly affect the reliability of predictions.
In order to obtain unbiased (data-independent) uncertainty estimates for both the empirical fit parameters and the corresponding predictions, we carried out a nonparametric 
bootstrap analysis. 
This analysis reveals that the standard deviation of the mean of the empirical D3 parameters is small. Moreover, the mean prediction uncertainty obtained by bootstrapping 
is not much larger than previously reported error measures. 
On the basis of a jackknife analysis, we find that the original reference set is slightly skewed, but our results also suggest that this feature 
hardly affects the prediction of dispersion energies.
Furthermore, we find that the introduction of small uncertainties to the reference data does not change the conclusions drawn in this paper.
However, a rigorous analysis of error accumulation arising from different parameterizations reveals that error cancellation does not necessarily occur, leading to a monotonically 
increasing deviation in the dispersion energy with increasing molecule size.
We discuss this issue in detail at the prominent example of the C_60 “buckycatcher”.
We find deviations between individual parameterizations of several tens of kcal mol^-1 in some cases.
Hence, in combination with any calculation of dispersion energies, we recommend to always determine the associated uncertainties
for which we will provide a software tool.




§ INTRODUCTION



Over the last decades, density functional theory (DFT) and in  particular Kohn–Sham DFT (KS-DFT)<cit.> 
has become the major workhorse for a broad range of applications in computational and theoretical chemistry, because it offers a satisfactory trade-off between the 
accuracy achievable and the computational resources required<cit.>. 
However, long-range electron correlation, and hence, London dispersion interaction, are not described by most exchange–correlation functionals 
(see, e.g., Refs. <cit.>). 
This problem has been addressed in a number of ways<cit.>. 
For example, Corminboeuf and coworkers developed a density-dependent dispersion correction for KS-DFT applicable to the entire periodic table <cit.>.
Roethlisberger and coworkers designed atom-centered, nonlocal potentials to take into account noncovalent interactions<cit.>.
Truhlar and coworkers designed exchange–correlation functionals which are highly parameterized to describe weak noncovalent 
interactions better<cit.>, whereas Dion et al. constructed a functional form which includes dispersion interactions in a 
nonempirical way<cit.>.
In this respect, the nonlocal correlation energy functional by Van Voorhis and coworkers, where the number of empirical parameters is kept to a 
minimum<cit.>, can be considered an intermediate between the former approaches.

To date, the most widely used approach is the introduction of semiclassical dispersion corrections, which are added to the electronic energy a posteriori.
They have been proposed as early as the 1970s<cit.>.
Significant advances of this direction were achieved in the 2000s by Becke and Johnson<cit.>, 
by Tkatchenko and Scheffler<cit.>, and by Grimme and coworkers<cit.> (Dx-type dispersion corrections).
Such corrections are designed to avoid double-counting of electron correlation effects.
Especially the so-called D3 approach by Grimme <cit.> has found widespread use as is evident from bibliometric tools such as Web of Science
and Google Scholar.
Therefore, D3 is the de facto standard for addressing nonlocal interactions in KS-DFT.

Like almost all parameterized physicochemical property models (e.g., approximations to the exact exchange–correlation functional), 
semiclassical dispersion corrections have not been subjected to a rigorous statistical analysis yet.
Here, we present such an analysis of D3-type dispersion corrections with a focus on parameter optimization, reference set selection, 
uncertainty propagation, and model transferability.
All Dx-type dispersion corrections (x ∈{ 2,3,4 }) comprise a small number (usually two or three) of global empirical parameters (apart from the C_n^AB 
coefficients that account for pair-specific contributions; see below) that are calibrated against a 
set of reference dispersion energies.
In this work, we consider the D3-zero and D3-BJ dispersion corrections (i.e., D3 models with a zero-damping function<cit.> and 
a Becke–Johnson (BJ) damping function<cit.>, respectively), the latter being the most widely applied D3 variant.
It is advisable to subject the corresponding reference sets (called fit sets in the original literature<cit.>) as well as the  
calibration procedure to a detailed statistical analysis to estimate the (co)variance of the empirical parameters and its effect on the prediction of dispersion energies.

Among other things, this strategy will allow us to estimate reliably the model prediction uncertainty (MPU), i.e., the expected (random) deviation 
of a prediction from a measurement or computational result — as of now, the MPU of semiclassical dispersion corrections has usually been assessed by calculating 
some error measure for a given validation set. 
Such MPU measures, however, have a strong dependence on the  validation set chosen.
The tools and methods required for an improved statistical analysis have only recently been embraced in quantum chemistry<cit.>. 
In this study, we revise the original calibration procedure of the D3 approach, and employ nonparametric bootstrapping<cit.> to obtain 
reliable estimates for both the uncertainty of the D3-BJ parameters and the corresponding MPU that is approximately independent 
of the validation set chosen.

We should note that the fourth generation of dispersion corrections by Grimme and coworkers (D4) has recently been published<cit.>. 
We deliberately chose the D3-BJ model as the object of our analysis because of its importance in current computational studies. We may also note that revised 
parameterizations for D3-zero and D3-BJ have been developed<cit.>, which, however, have never attracted as much attention as the 
original parameterization.
Furthermore, since the dispersion corrections by Grimme and coworkers feature the same functional form for a given damping function (e.g., D3-BJ and D4-BJ only 
differ in the values of the dispersion coefficients C_n^AB, see below), one can expect our results to be generalizable.

This paper is organized as follows: 
In Section <ref>, we detail the methodology adopted in this study. 
In Section <ref>, we provide a detailed analysis of the D3-zero and D3-BJ parameterizations.
We further examine the effect of different weighting schemes on the D3-BJ parameterization in Section <ref>.
In Section <ref>, we demonstrate how to infer the uncertainty of the empirical D3-BJ parameters on the basis of a bootstrap analysis, before we harness 
the related jackknife method in Section <ref> to identify reference data that bias the D3-BJ  parameterization.
Subsequently, in Section <ref>, we discuss the effect of uncertainty in the reference data on the parameterization; an issue that is 
disregarded in most quantum chemical calibration studies.
Finally, the insights obtained will serve to provide an in-depth assessment of the molecule-size dependence of the model prediction uncertainty of the dispersion 
corrections in Section <ref>.



§ COMPUTATIONAL METHODOLOGY



The generic D3 model by Grimme <cit.> is an additive a posteriori correction to the electronic energy,

    E_disp^D3≡  -∑_A,B>A ( s_6 C_6^AB/R_AB^6f^(6)_damp(R_AB) + s_8 C_8^AB/R_AB^8f^(8)_damp(R_AB) )  ,

where A and B identify a pair of atoms in the system under consideration.
The dispersion coefficients C_n^AB (n = 6, 8) are completely determined by the nuclear coordinates and atomic numbers and, therefore, not affected by the calibration procedure 
chosen for estimating optimal parameter values.
The damping function, f^(n)_damp(R_AB), is a function of  pairwise nuclear distances, R_AB, and requires further specification.
The original zero-damping function<cit.>,

    f^(n)_damp,zero(R_AB) ≡1/1 + 6 (R_AB / (s_r,nR_0^AB))^a_n ,

contains one empirical (or adjustable) parameter, s_r,6, with fixed values for s_r,8, a_6, and a_8.
The improved BJ damping function<cit.>,

    f^(n)_damp,BJ(R_AB) ≡R_AB^n/R_AB^n + (a_1R_0^AB + a_2)^n ,

contains two empirical parameters, a_1 and a_2.
In both cases, the cut-off radius R_0^AB is defined as the square-root ratio of the corresponding C_6^AB and C_8^AB dispersion coefficients,

    R_0^AB≡√(C_8^AB/C_6^AB) .

Independent of the damping function chosen, the global scaling factor s_8 constitutes an additional empirical parameter, whereas s_6 is set to one.
The parameter s_8 is adjustable in order to compensate for double-counting effects<cit.> and is functional-dependent (like the empirical 
parameters of the damping functions).
In line with the work by Grimme and coworkers<cit.>, we define the dispersion coefficients C_n^AB to be strictly positive, 
i.e., the calculated dispersion energy will be smaller than zero and 
asymptotically approach zero for increasing distances R_AB.

For all calculations of dispersion energies, E_disp^D3, Grimme's implementation of D3-zero and D3-BJ retrieved from the web page given in 
Ref. <cit.> was employed. 
The program was minimally modified such that the values of the empirical parameters (s_r,6 and s_8 in the case of D3-zero, a_1, s_8, and a_2 in the case 
of D3-BJ) could be passed as command-line arguments instead of parsing an input file,
thereby leading to both significant computational speedup and easy parallelizability; features that are mandatory for efficient bootstrapping.

All reference molecular structures and dispersion energies were obtained from the following data sets of the GMTKN24 database<cit.> 
for D3-zero, the GMTKN30 database <cit.> for D3-BJ, and the supporting information of Ref. <cit.> for both D3-zero and D3-BJ:
S22, S22+, RG6, ADIM6, ACONF, CCONF, PCONF, and SCONF.
Note that the S22 and ACONF reference subsets are different in GMTKN24 and GMTKN30.
The reference dispersion energy is defined as the difference in the intermolecular electronic interaction energy obtained with a (highly accurate) benchmark model 
and with a density functional (DF),

    Δ E_disp^ref≡	 E_inter^benchmark - E_inter^DF .

The electronic interaction energy calculated with model X, in turn, is defined as the electronic energy difference between the entire system and 
its subsystems (here, two; “sub1” and “sub2”),

    E_inter^X ≡ E_el,system^X - ( E_el,sub1^X + E_el,sub2^X )  .

Both E_inter^benchmark and E_inter^DF are taken from the databases mentioned above and were not recalculated by us. I.e.,
in this work, no electronic structure calculations had to be carried out.

All statistical protocols discussed in this paper were already introduced by us before<cit.> and made available in our open-source calibration 
software reBoot<cit.>.
Here, we converted the Octave scripts of reBoot to the Python language.
These custom-made Python scripts rely on the packages NumPy 1.13.2<cit.> and 
SciPy 0.19.1<cit.>. 
For D3-zero, the calibration was accomplished by a brute-force search as implemented in SciPy 0.19.1, whereas for D3-BJ, differential 
evolution<cit.> (also implemented in SciPy 0.19.1) was chosen since the increased parameter space would render a brute-force 
search combined with bootstrapping inefficient. 
For D3-zero, we considered the following parameter ranges: 1.0 ≤ s_r,6≤ 1.6, 0.5 ≤ s_8 ≤ 1.1, and for D3-BJ, it is 
0.0 ≤ a_1 ≤ 0.7, 0.0 ≤ s_8 ≤ 3.5, 2.5 ≤ a_2 ≤ 6.5. 
We carefully set all numerical convergence thresholds tight enough. 
Our preliminary tests showed that the default values yield satisfactory convergent behavior, except for the relative (dimensionless)  parameter 
of the differential evolution algorithm, which was set to 0.001. 
For all bootstrap runs, a sample size of 10'000 was chosen (see also Section <ref>). 
For the generation of bootstrap samples, we exploited the default random number generator of Python as implemented in the module .
The parameters obtained in the bootstrap runs will be made available as part of our program BootD3 (see also Section <ref>).

All plots were produced with Mathematica 11.0.1<cit.>.



§ RESULTS AND DISCUSSION





 §.§ Parameter Optimization for D3-zero and D3-BJ



For our statistical analysis, we implemented a calibration routine which yields optimal values for the empirical parameters (s_r,6 and s_8 in case of D3-zero, 
and a_1, s_8, and a_2 in case of D3-BJ, respectively). 
This routine minimizes a cost function and forms the basis for our bootstrap analysis.
The cost function employed in this work is the mean absolute error (MAE) between the reference intermolecular dispersion energies, Δ E_disp^ref 
(taken from Refs. <cit.> for D3-zero and Refs. <cit.> for D3-BJ), and 
the intermolecular dispersion energies calculated for a given set of values of the empirical parameters, Δ E_disp^D3,

    MAE≡1/N∑_i=1^N w_i |Δ E_disp,i^ref - Δ E_disp,i^D3| ,

where

    Δ E_disp,i^D3≡ E_disp,system,i^D3 - ( E_disp,sub1,i^D3 + E_disp,sub2,i^D3)  .


Here, N is the number of reference data points and w_i refers to the weight of the i-th data point that will be discussed later.
Note that in the original work<cit.> by Grimme and coworkers, the MAE is termed MAD (mean absolute deviation).
The MAD is a more general quantity which is defined as the average of the absolute deviations from a central point such as the mean of the absolute residuals (difference 
between reference and prediction) <cit.>.
If that central point is set to zero, both MAD and MAE are identical quantities and as such employed in Grimme's and our studies.
We will apply the more specific term MAE in this paper, even though it is interchangeable with Grimme's definition of the MAD. 

First, we reproduced the results reported by Grimme and coworkers for D3-zero<cit.>. 
In fact, when we adopted the parameter values reported by Grimme et al. to be optimal for D3-zero, we reproduced Grimme's MAEs for all density functionals studied
(PBE, revPBE, TPSS, BP86, B97-D, B3LYP, and B2PLYP). 
In a few cases, we observed negligible differences on the order of 0.01 kcal mol^-1.
These results clearly suggest that our MAE definition provided in Eq. (<ref>) is essentially identical to that used by Grimme and coworkers.
Two alternative MAE definitions,

    MAE2≡1/N∑_i=1^N |Δ E_disp,i^ref - Δ E_disp,i^D3|

and

    MAE3≡1/∑_i=1^N w_i∑_i=1^N w_i |Δ E_disp,i^ref - Δ E_disp,i^D3| ,

resulted in significant deviations from the MAEs reported by Grimme and coworkers.

Concerning the MAEs reported by Grimme et al. for D3-BJ<cit.>, we found that we could reproduce these values for the S22 and PCONF subsets for all 
density functionals studied. 
However, when comparing the MAEs over the entire reference set, we found significant deviations for the density functionals BP86 (0.16 kcal mol^-1), B97-D 
(0.03 kcal mol^-1), B3LYP (0.04 kcal mol^-1), and B2PLYP (0.05 kcal mol^-1).
In all these cases, our MAEs were higher than those reported in Ref. <cit.> (note that the MAEs in that reference had to be extracted from a 
figure).
We may assume that these deviations result from differences in the other subsets; especially ACONF, since this reference data set is, as S22, different 
in GMTKN24 and GMTKN30<cit.>.
However, subset-specific MAEs were only reported for S22 and PCONF in Ref. <cit.>.

It is instructive to visually inspect the cost function to be minimized. The cost function for the PBE density functional is shown in Fig. <ref>.
The MAE varies between zero and very large values of more than 10 kcal mol^-1 over the large parameter range considered. 
However, the optimization problem is not too complex since the parameter dependence of the cost function does not exhibit a rugged form, but is usually convex.
Still, there can be cases where the cost function exhibits more than a single minimum, as shown in Fig. <ref> c). In this case, the cost function has two minima, 
one around s_8 ≈ 0.4 and a second one around s_8 ≈ 1.1. 
A more detailed investigation shows that the minimum around s_8 ≈ 1.1 is the global minimum. 
As indicated by the crosses in Fig. <ref>, the differential evolution algorithm adopted in our study reliably finds the global minimum, also for nonconvex cost 
functions.

The optimal parameters obtained with our optimization procedure are reported in Tables <ref> and <ref> for D3-zero and D3-BJ, respectively.
At this point, we should note that we report parameters with a precision of four decimal places in this work (in some exceptional and justified cases we will 
divert from this choice). 
We emphasize that this choice does not indicate the number of statistically significant figures; we chose this information content to be consistent 
with the representation of the empirical D3-BJ parameters in Grimme's D3 implementation<cit.>.
In fact, as we will see, (statistical) parameter uncertainty usually already affects the third decimal place. Still, we will report more digits for the 
sake of consistency and in order to rule out any rounding issues. Of course, many more digits were taken into account in the calculations, so that numerical
(cut-off) errors have a negligible effect on our results.





We obtained slightly different parameter values for D3-zero compared to those reported by Grimme and coworkers.
However, all qualitative trends between the Grimme parameterization and ours are the same.
For example, in both parameterizations, the largest value for s_r,6 was found for B2PLYP, while the smallest value was found for B97-D. 
The differences between the MAEs obtained with the two parameterizations were also very small. 
Still, the two parameterizations sometimes led to significant MAE differences, especially for BP86, where we found an MAE of 0.56 kcal mol^-1 while the original 
D3-zero parameterization produced an MAE of 0.65 kcal mol^-1.
These differences for D3-zero are most likely a result of different cost functions employed for parameter optimization.
Grimme and coworkers briefly mention the use of a weighted least-squares calibration procedure <cit.>.
However, when we substituted the absolute value of residuals, | E_disp,i^ref - E_disp,i^D3|, in Eq. (<ref>) by its 
squared value, ( E_disp,i^ref - E_disp,i^D3)^2, the deviations from Grimme's parameters were still significant and even 
increased in some cases.
In the following, all results reported here refer to the MAE cost function introduced in Eq. (<ref>).
We emphasize that the representation of residuals (here, absolute versus squared values) in the cost function does not affect the conclusions drawn in 
this paper.



For D3-BJ, we find more pronounced deviations from the parameters reported by Grimme and coworkers<cit.> compared to D3-zero in some cases, e.g., for 
s_8 when considering PBE and B97-D (see Table <ref>).
Due to the additional empirical parameter to be optimized in D3-BJ compared to D3-zero, the solution space is of higher dimension and may, therefore, be more sensitive 
to differences in the calibration procedure (e.g., cost function, reference set, optimization algorithm).
Still, as in the case of D3-zero, all qualitative trends are the same.
We should note that even though the differences between the MAEs are only on the order of a few hundredths kcal mol^-1, their relative differences can be significant.
For example, for both BP86 and B97-D the MAE obtained with our parameterization is lowered by 10 % compared to the MAE evaluated with the original parameterization, 
whereas for B2PLYP, the MAE obtained with the original parameterization is lower by 15 %. This last result further corroborates our assumption
mentioned above in the discussion of different MAE definitions that the 
reference set employed in this study for D3-BJ is different from the one used in the original parameterization of 
D3-BJ<cit.>.
Since we optimized the D3-BJ parameters by minimizing the MAE, another parameterization cannot lead to a lower MAE.
In fact, when adopting the original D3-BJ parameters, we consistently find higher MAEs compared to our parameterization.
However, as can be seen in Table <ref>, Grimme et al. report lower MAEs for BP86, B3LYP, and B2PLYP compared to our data.

It is imperative to test the performance of a given set of empirical parameters on an independent validation set, i.e., a data set which has not been employed 
in the calibration routine. 
To this end, we provide the MAEs resulting from our D3-BJ parameterizations for three validation sets — all taken from the 
GMTKN30 database<cit.> — in Table <ref>.





The reliability of our new parameterization is comparable to that of the original parameterization.
The difference between the corresponding MAEs for the two parameterizations is usually very small (on the order of 0.01 kcal mol^-1). 
In a few cases (e.g., AL2X with BP86 and DARC with B97-D), the difference between the MAEs resulting from the two parameterizations is significantly larger, 
reaching almost 1 kcal mol^-1.
If the overall effect of the two slightly different parameterizations may become large in some cases, this will raise the question whether such effects based on 
a specific choice for a parameter set are actually negligible in routine quantum chemical calculations.





Note that in the case of AL2X (with BP86) and DARC (with B97-D), also the absolute errors are comparatively large. 
To investigate the reasons for this, we calculated the errors for the individual data points of these two sets.
On average, the MAE resulting from Grimme's original parameterization was 135.6 % of the MAE resulting from our parameter values in the case of the AL2X data set 
(Table <ref>). 
This deviation was mostly created by Al_2Cl_6, but also by Al_2Me_4, Al_2Me_5, and Al_2Me_6.
When closer analyzing the dispersion-energy difference determined for the Al_2Cl_6 data point, we found that for AlCl_3 (one of the two identical monomers of 
Al_2Cl_6), the total dispersion energy is -5.87 kcal mol^-1 for Grimme's original parameterization and -4.86 kcal mol^-1 for our new parameters.  
When analyzing the individual pairs, we found that for all atom pairs, the dispersion energy is smaller in Grimme's case. 
We made the same observation for Al_2Cl_6. 
For the DARC validation set, we found the same trend, but the relative deviations are less pronounced (Table <ref>).



We suggest two reasons for large absolute dispersion-energy errors of individual systems.
First, the prediction model (here, D3-BJ) is inadequate such that certain domains of the chemical space are systematically misdescribed, no matter what parameterization 
would be chosen (see also Section <ref> below).
Second, since the number of pairwise interactions depends on system size, small inaccuracies in the calibration procedure may accumulate to large errors for large 
molecular systems.
This latter error type will be discussed in detail in Section <ref>.



 §.§ Consequences of Different Weights on Model Performance



In the original calibration routine of Grimme and coworkers<cit.>, the RG6 set (containing six rare gas dimers) was given a weight of 20, while all 
other data points were assigned a weight of one. 
It is not immediately clear why such a weighting scheme should be  adopted, nor is this choice explained in Refs. <cit.>. 
Therefore, we decided to carry out a reparameterization where this increased weighting was removed, hence considering all data points on an equal footing. 
In the following, we refer to this variant as the “uniform weighting scheme”, for which we find MAE = MAE2 = MAE3.
The resulting optimal parameter values together with their corresponding MAEs with respect to the reference set are listed in Table <ref>.

By comparison with Table <ref>, we note that the parameter values sometimes change significantly. 
For example, in the case of PBE, a_1 decreases from 0.4309 Bohr to 0.0181 Bohr when adopting the uniform weighting scheme. 
Also for TPSS, B97-D, and BP86 we note rather strong differences between the parameter values obtained with the original and the uniform weighting scheme.
For the other density functionals, however, the changes are much smaller. 
It is not surprising that the MAE with respect to the reference set is generally smaller in the case of the uniform weighting scheme.
Increasing the weight of a given data point can be achieved by repeatedly adding the data point to the reference set until the desired weight has been reached;
cf. Eq. (<ref>). 
For example, in the present case, each data point of the RG6 set would be added twenty times to the reference set. 
In the uniform weighting scheme, each of these data points occurs only once, and so the reference set is effectively smaller.
It is natural that the (cumulative) error with respect to a reference set always depends on its size. 
Therefore, we need to analyze the MAE obtained on some validation sets in order to develop a clear picture of the reliability of this new parameterization.






The MAEs obtained for the validation sets HEAVY28, AL2X, and DARC are provided in Table <ref>.
The uniform weighting scheme sometimes leads to a significantly reduced reliability (cf. Table <ref>). 
For instance, for PBE, the MAE (Eq. (<ref>)) with respect to the AL2X set increases from 2.14 kcal mol^-1 to 5.54 kcal mol^-1. 
In some cases, the uniform weighting scheme also leads to improved results. 
Most interestingly, for the DARC set, a lower MAE is found for most density functionals studied in this work.



Therefore, neither the original weighting scheme adopted by Grimme and coworkers nor the uniform weighting scheme is optimal in all cases, but 
both weighting schemes lead to a good overall reliability.
From the results above it can be expected that the numerical effect that different weighting schemes have on dispersion corrections in actual quantum chemical 
calculations can become significant. 
We will show in Section <ref> that different weighting schemes may have a dramatic impact on the dispersion energy error for large molecular systems.



 §.§ Bootstrapping Model Prediction Uncertainty



The MPU is an important key characteristic of every model. 
However, it is often difficult to quantify reliably the MPU. 
A first estimate for the MPU is the residual error with respect to the reference set (such as the MAE), but such an estimate depends on the actual reference set chosen. 
To obtain an MPU measure independent of the reference set, one often resorts to validation sets (as we have done above in Sections <ref> and <ref>).

Another way of estimating the MPU without the need of additional validation sets is the so-called nonparametric bootstrapping approach<cit.>, which 
we recently implemented for the calibration of physicochemical property models and discussed in detail at the example of ^57Fe Mössbauer isomer shift 
prediction<cit.>.
In this approach, one generates a large collection of data sets, sampled from the original reference set, to reproduce approximately the population distribution 
underlying the original reference set.
For each of these data sets (bootstrap samples), an individual calibration is carried out, leading to an ensemble of parameter sets.
This ensemble allows us obtain an estimate m for the expected value μ, 

    μ(v) ≈	m(v) ≡1/B∑_b=1^B v_b  ,
		
as well as an estimate s^2 for the (co)variances σ^2,

    σ^2(v,w) ≈ s^2(v,w) ≡ 1/B-M∑_b=1^B (v_b - m(v)) (w_b - m(w)) ,
		
 of the generic parameters v and w.
 Here, B and M denote the numbers of bootstrap samples and empirical parameters, respectively, and v_b refers to the optimal v-value for the b-th bootstrap sample.
 For v = w, we additionally define s^2(v,v) ≡ s^2(v).
 Moreover, instead of studying the low-order moments (expected value, covariance) of the parameter distributions, their specific functional form can be approximated 
through histograms that, e.g., facilitate to judge whether the original reference set <cit.> is balanced or skewed (as will be discussed below).
Furthermore, the ensemble of bootstrap samples allows us to generate an ensemble of MAEs.
The mean over all MAEs is an estimate to the actual MPU,

    MPU≈ m(MAE) = 1/B∑_b=1^B MAE_b  ,

where MAE_b constitutes the MAE determined from the b-th parameter set with respect to the original reference set.
The difference between max(MAE) and min(MAE) defines the maximum MAE spread sampled.



Table <ref> contains the mean values for all three empirical D3-BJ parameters as obtained from our bootstrapping procedure together with the 
corresponding standard deviations. 
By comparison with Table <ref>, we see that the parameters obtained from bootstrapping are not significantly different from the original parameters, 
and all qualitative trends are the same.
Concerning the standard deviations given in Table <ref>, we should note that these are the direct standard deviations of the parameter distributions 
as obtained from the bootstrapping procedure. 
Usually, however, one is not so much interested in this standard deviation but more in the standard deviation of the mean of these distributions. 
Assuming these parameters to be distributed normally, the standard deviation s(m(v)) of the mean m(v) of a given parameter v can be determined on the basis 
of the well-known formula s(m(v)) = s(v)/√(B), where B, in our case, equals 10'000. 
Therefore, the standard deviation of the mean of the empirical parameters a_1, s_8, and a_2 is smaller by a factor of 100 compared to the direct standard deviations. 
Hence, based on the results shown in Table <ref> and assuming a Gaussian distribution, we found the following empirical parameters 
together with a two-sided 95 % confidence interval when employing the PBE density functional: a_1 = 0.4191 Bohr ± 0.0014 Bohr; s_8 = 1.237 ± 0.007; 
a_2 = 4.955 Bohr ± 0.006 Bohr
To assess the reliability of assuming a Gaussian distribution, we also bootstrapped the mean of the PBE-D3-BJ parameters.
We found the following mean values and two-sided 95 % confidence intervals: a_1 = 0.4191 Bohr ± 0.0014 Bohr; s_8 = 1.237 ± 0.007; 
a_2 = 4.955 Bohr ± 0.006 Bohr
These results are literally identical to the ones based on the assumption of a normal distribution and we may safely assume that this
assumption is also true for all other density functionals.
Strictly speaking, it is not reasonable to report more than three decimal places for s_8 and a_2, since any figure after the third decimal place is affected 
by parameter uncertainty.
However, as we have already stated above, we will continue to report four decimal places for the sake of consistency (also recall that many more digits were
taken into account during the calculations so that cut-off errors can be neglected).


We now analyze the mean prediction errors obtained from the individual bootstrap samples. 
A summary is given in Table <ref>. 
The average MAEs obtained from bootstrapping are only slightly larger than the MAEs obtained with the original parameterization (cf. Table <ref>), 
even though the largest MAEs are often more than twice as large as the mean.





Histograms of the distributions of the results obtained for a_1, s_8, and a_2 are presented in Fig. <ref>.
While most parameters are centered around the optimal value obtained by the standard optimization procedure as discussed in Section <ref>, we note some distinct 
outliers for all parameters, i.e., there are a few bootstrap samples which lead to very different parameters. 
This finding suggests that the original reference set is not balanced but skewed in the sense that it contains a few data points which significantly bias the optimization.





To examine this issue in more detail, we study the MAE with respect to the original reference set as a function of the empirical D3-BJ parameters a_1, s_8, and 
a_2 for all 10'000 bootstrap samples (Fig. <ref>).
We observe an interesting pattern in each MAE–parameter plot.
It appears that subsets of MAE–parameter pairs form parabolic clusters.
When analyzing all MAEs smaller than 0.49 kcal mol^-1, we cannot find a single molecular system which is either present or absent in all data points.
Moreover, we do not find that any reference subset is over- or underrepresented (on average) in the bootstrap samples yielding an MAE <0.49 kcal mol^-1.
On the contrary, we find that every bootstrap sample that yields an MAE larger than 0.70 kcal mol^-1 does not contain any of the following three RG6 systems: 
the Rn–Rn dimer, the Rn–Xe dimer, and the Xe–Xe dimer.
It is not straightforward to decide whether the absence of these dimers biases the parameterization or not.
Intuitively, one may argue that this is not the case since each of these dimers is assigned a weight of 20.
Consequently, any contribution of these dimers to the MAE is also given this weight, cf. Eq. <ref>.

We discuss the issue of biasing reference data by means of the jackknife method <cit.> in the next Section.



 §.§ Jackknifing the Bias of Reference Data



In the first step of the jackknife method, one generates N new data sets (jackknife samples) from the original reference set containing N data points, each 
containing all but one data point from the original reference set. 
In each of these jackknife samples, the missing data point is a different one; hence, in the first jackknife sample, the first data point is absent, whereas in the 
second sample, the second point is absent, and so forth. 
Then, the optimal parameter values are determined for each of these jackknife samples. 
If the original reference set is not biased, the distribution of the parameter values obtained by the jackknife method should indicate a smooth function. 
If, however, the original reference set is biased by certain data points, the jackknife samples with those data points removed will indicate nonsmooth parameter distributions.

The parameter distributions obtained with the jackknife method are shown in Fig. <ref>. 
As already expected on the basis of the bootstrapping analysis, we find distinct outliers in all three histograms of Fig. <ref>.
For instance, there are two outliers in the distributions of a_1 and a_2, respectively. 
A closer analysis revealed that for both distributions, these outliers correspond to the jackknife samples where the Rn–Xe dimer and the Xe–Xe from the RG6 set have 
been removed.
Hence, these two dimers heavily bias the original calibration procedure (recall that all compounds of the RG6 set are given an increased weight of 20 in the original 
calibration routine, see Section <ref>).
Interestingly, the absence of these dimers and the Rn–Rn dimer led to a significant increase of the MAE with respect to the original reference set as inferred from 
our bootstrap analysis presented above.
While Rn–Xe and Xe–Xe seemed to be important for the parameterization in our bootstrapping analysis (though we noted that this interpretation may have been caused by 
the large weight assigned to these systems), they now appear to bias the parameterization.
Removing the Rn–Xe and Xe–Xe dimers from the original reference set and repeating the optimization, we found the following values for the PBE density functional: 
a_1 = 0.3969 Bohr; s_8 = 1.1801; a_2 = 5.0241 Bohr 
The MAE with respect to the original reference set amounts to 0.48 kcal mol^-1, which is only 0.01 kcal mol^-1 lower compared to our standard calibration 
procedure.
We emphasize that the MAE is not as large as 0.70 kcal mol^-1 in this case as it is now given with respect to the new reference data set in which these two dimers 
have been removed.
For the validation sets, we found the following MAEs: HEAVY28: 0.37 kcal mol^-1; AL2X: 2.24 kcal mol^-1; DARC: 4.00 kcal mol^-1. 
These values are slightly larger than for the original parameters, except for DARC, where the MAE is reduced by 0.02 kcal mol^-1, which is, however, not significant. 
Therefore, even though the original reference set is clearly biased, the resulting parameters are not less reliable than parameters resulting from a parameterization 
where this bias was removed, at least as measured by the validation sets chosen in this study.
For different validation sets, especially those representing chemically different or large molecular systems, the situation could change unpredictably.
We will discuss this issue in detail in Section <ref>.



In general, the fact that a given reference set is skewed indicates that  the chemical space is not represented uniformly well. 
Therefore, for domains that are chemically rather dissimilar to the domain covered by the reference set, significant prediction errors are to be expected.





 §.§ Consideration of Reference Uncertainty



Irikura et al. were the first to discuss the inherent uncertainty of quantum chemical data in detail <cit.>.
While a quantum chemical result is effectively reproducible, it is biased due to the introduction of approximations concerning, e.g., the Hamiltonian 
or the construction of the wave function including basis-set truncation.
This bias is generally unknown and can only be accessed through statistical modeling, which provides this bias as an uncertainty.
It is important to examine how robust the D3-BJ parameters are with respect to uncertainties in the reference data. 
To this end, we have repeated the bootstrapping procedure (B = 10'000) for the PBE density functional, but during the generation of the bootstrap samples, we 
replaced the reference energy of each data point by an artificially generated value.
This value is a random number, taken from a Gaussian distribution with the mean of the original reference value, and a standard deviation of 5 % of this reference 
value.
Clearly, the choice of a 5 % standard deviation is rather arbitrary, but it is intended to reflect the expectation that the reference data employed 
are already of high quality.



The new parameter values as well as their standard deviations are reported in Table <ref> and compared to the corresponding values obtained from 
the original bootstrapping analysis (see Section <ref>). 
None of the empirical parameter values changed significantly upon adding a small uncertainty to the reference data. 
Also the resulting standard deviations are almost unaffected.
Furthermore, the resulting prediction errors (also reported in Table <ref>) are virtually unchanged. 
Our findings suggest that the D3-BJ parameterization is rather robust with respect to small uncertainties in the reference data employed.





 §.§ Effect of Molecular Size on Model Prediction Uncertainty



In this section, we investigate how deviations of pairwise dispersion interactions induced by different parameterizations propagate to total dispersion energies.
For this purpose, we introduce some reduced notation for the sake of clarity.
First, we define E to be the total D3-BJ dispersion energy for the original parameterization,

    E ≡ E_disp^D3-BJ = -∑_AB∑_n=6,8 s_n C_n^AB/R_AB^n + (a_1 R_0^AB + a_2 )^n .

The pairwise dispersion energies are abbreviated as

    E_n^AB≡ -C_n^AB/R_AB^n + (a_1 R_0^AB + a_2 )^n ,

which sum up to the total dipole–dipole (n = 6) and dipole–quadrupole (n = 8) interaction terms

    E_n ≡∑_AB E_n^AB .

Since s_6 = 1, the total dispersion energy for the original parameterization reads

    E = E_6 + s_8E_8  .

Note that E_6 and E_8 are system dependent and that we do not indicate this dependence to keep the notation uncluttered.

Now, we determine the implications of a modified parameterization on the pairwise and total dispersion energies.
We begin with the empirical parameters of the BJ damping function, a_1 and a_2.
We express the new parameter a_1,b obtained from another calibration (e.g., against a new reference set indicated by the index variable b) as the 
sum of the original parameter a_1 and a sample-dependent parameter d_b,

    a_1,b≡ a_1 + d_b  .

Similarly, we define the new parameter a_2,b as the sum of a_2 and the product of d_b and another sample-dependent parameter ε_b,

    a_2,b≡ a_2 + ε_b d_b  .

As we have seen in the preceding sections, both a_1 and a_2 are always larger than zero. 
Consequently, we obtain the boundary conditions d_b > -a_1 and ε_b > -a_2/d_b for the sample-dependent parameters.
The empirical damping term a_1R_0^AB + a_2 changes accordingly for a different parameterization,

    a_1,bR_0^AB + a_2,b = a_1R_0^AB + a_2_≡ p^AB,  original term + d_b (R_0^AB + ε_b )_≡ q_b^AB ,

and can be written as the sum of the original empirical damping term, p^AB, and a sample-dependent term, q_b^AB.
With this definition, we can determine some general properties that hold for all pairwise interactions under certain conditions; ∀ AB holds

    p^AB
    < p^AB + q_b^AB if ε_b > -min(R_0^AB)  ,  d_b > 0  ,  or ε_b < -max(R_0^AB)  ,  d_b < 0 
    
    > p^AB + q_b^AB if ε_b > -min(R_0^AB)  ,  d_b < 0  ,  or ε_b < -max(R_0^AB)  ,  d_b > 0 
     .

For the pairwise dispersion energies arising from the parameterization of the b-th reference set,

    E_n,b^AB≡ -C_n^AB/R_AB^n + (p^AB + q_b^AB)^n ,

an analogous analysis is possible: ∀ AB we have

    E_n^AB
    < E_n,b^AB if ε_b > -min(R_0^AB)  ,  d_b > 0  ,  or ε_b < -max(R_0^AB)  ,  d_b < 0 
    
    > E_n,b^AB if ε_b > -min(R_0^AB)  ,  d_b < 0  ,  or ε_b < -max(R_0^AB)  ,  d_b > 0 
     .

Note here that all pairwise dispersion energies are strictly negative quantities.
This analysis reveals that for a significant domain of ε_b, the deviation between E_n and

    E_n,b≡∑_AB E_n,b^AB

will increase with every pairwise dispersion interaction, i.e., the difference between these two quantities,

    Δ_n,b≡ E_n,b - E_n < -E_n  ,

will constantly increase with an increasing number of pairwise interaction terms and, hence, for increasing sizes of molecular systems.

Error cancellation with respect to the original parameterization can only occur in two ways.
First, the condition -max(R_0^AB) < ε_b < -min(R_0^AB) is fulfilled as illustrated in Fig. <ref> (we will refer to
this a type-I error cancellation below).
For a diatomic system, the white error-cancellation domain in Fig. <ref> reduces to a single value, R_0^AB, which depends on the atomic 
numbers of the two atoms A and B.
This example illustrates that the error-cancellation domain is system-dependent, i.e., its position and width will be generally different for different molecular systems.
Consequently, if error cancellation may occur for one molecular system, it may not be found for another system.
Moreover, the further the distance between ε_b and the error-cancellation domain and the larger the absolute value of d_b, the larger 
|Δ_n,b| will become.



The second possibility for error cancellation (type-II error cancellation) depends on the s_8 parameter in the b-th reference set,

    s_8,b≡ s_8 + δ_b  ,

where we can postulate δ_b > -s_8 without loss of generality since dispersion energies are strictly negative quantities.
When rearranging the expression for the total dispersion energy corresponding to the b-th parameterization,

    E_b ≡ E_6,b + s_8,bE_8,b =    E_6 + s_8E_8    +  Δ_6,b + s_8Δ_8,b + δ_bE_8,b
    
    =    E    +  Δ_6,b + s_8Δ_8,b + δ_bE_8,b ,

we find that E_b is the sum of the original total dispersion energy, E, the change in the total dipole–dipole interaction, Δ_6,b, the change in the 
total dipole–quadrupole interaction, s_8Δ_8,b, and the product of the change in s_8, δ_b, and the total dipole–quadrupole interaction of the 
b-reference set, E_8,b.
If ε_b does not lie in the error-cancellation domain, both Δ_6,b and s_8Δ_8,b will be either strictly positive or negative; and since 
E_8,b is always strictly negative, error cancellation can only occur if δ_b has the same sign as Δ_n,b (n = 6, 8).
On the contrary, if ε_b does lie in the error-cancellation domain, Δ_6,b and s_8Δ_8,b will have smaller absolute values or even 
different signs, and the importance of δ_b as an error-cancellation parameter decreases.
Still, only if δ_b has the same sign as Δ_6,b + s_8Δ_8,b, further error cancellation is possible.
At the same time, no matter what value ε_b takes, δ_b can also have an error-enhancing effect, i.e., if it shows the opposite sign to 
Δ_6,b + s_8Δ_8,b.

We emphasize that all terms in Eq. (<ref>) are system-dependent and, therefore, the accumulation of error needs to be studied case by case.
Clearly, for a set of small molecular systems, the error spread can be expected to be smaller than for a set of large systems.
We understand that the accumulation of error is a general property of models building on sums of pairwise interactions and is not a result of the specific 
formulation of the D3 approach or any other semiclassical dispersion correction.
Rather, the success of such models depends on the construction of the pairwise interaction terms themselves — assuming that a sum over pairwise interactions 
is the dominant contribution to the total energy.
The more accurately these terms describe the underlying physical effects, the smaller will be the expected error per term.

To develop an intuition of how strongly errors may accumulate with an increasing number of pairwise D3-BJ interaction terms, we will first study the correlations 
between the three empirical parameters a_1, s_8, and a_2, which are reported in Table <ref>.
The correlation between two generic parameters v and w is defined as the covariance of v and w divided by the product of the standard deviations of v and w,

    r(v,w) = s^2(v,w)/s(v)s(w) .




It follows from Eq. (<ref>) that the first type of error cancellation can only occur if ε_b < 0 (see also Fig. <ref>).
In these cases, cf. Eqs. (<ref>) and (<ref>),  a_1 increases (decreases) and a_2 decreases (increases).
Consequently, our bootstrapped negative correlation between a_1 and a_2, r(a_1,a_2) = -0.74, indicates a tendency for type-I error cancellation.
Note however, that an increase (decrease) in a_1 coinciding with a decrease (increase) in a_2 does not imply error cancellation.

Next, we find that the correlation between a_1 and s_8 is relatively large and positive, r(a_1,s_8) = 0.51, compared to the positive correlation between 
a_2 and s_8, r(a_2,s_8) = 0.17.
Given that Δ_6,b + s_8Δ_8,b is positive, which indicates a tendency for an increase in a_2 whereas no direct tendency can be derived for the 
change in a_1.
However, based on the negative correlation between a_1 and a_2 which we estimated through bootstrapping, we expect a decrease in a_1 on average if 
Δ_6,b + s_8Δ_8,b > 0.
In this case, type-II error cancellation can only occur if δ_b is also positive.
Since we expect a_1 to be decreasing, we also expect s_8 to be decreasing due to the positive correlation between a_1 and s_8, which would however, imply, 
δ_b < 0.
Hence, we do not find direct evidence for a tendency towards type-II error cancellation.

Clearly, error accumulation does not occur uniformly over all pairwise dispersion interactions.
If the distance between two atoms is rather large, we expect the possible error in the corresponding interaction term to be smaller than for two atoms that are closer.
Therefore, one should expect the sampled deviation in the total dispersion energy to be a function of some effective system size, Ω_eff, which we define here as

    Ω_eff = ∑_ABf^(6)_damp,BJ(R_AB)/R_AB^6 = ∑_AB1/R_AB^6 + (a_1R_0^AB+a_2)^6 ,

where a_1, s_8, and a_2 are the parameters obtained from calibration against the original reference set <cit.>.
With this definition, the effective system size is a monotonically decreasing function of R_AB that asymptotically approaches zero. Note that
with this definition, the resulting effective size values are comparatively small.



In Fig. <ref>, the bootstrapped standard deviation of the total dispersion energy, s(E_disp^D3-BJ), is shown as a 
function of the effective system size, Ω_eff.
This and all following results were produced with the PBE-D3-BJ model.
We find a significant increasing trend, which is clearly linear for most of the individual reference subsets, e.g., for the homologous series of alkanes 
(ADIM6), we find maximum correlation between s(E_disp^D3-BJ) and Ω_eff (perfect linearity).
The increase of s(E_disp^D3-BJ) per unit Ω_eff is 366.59 kcal mol^-1 for the ADIM6 set and is the smallest among all reference 
subsets.
For the S22 and S22+ sets, which comprise many unsaturated organic molecules, the increase of s(E_disp^D3-BJ) per unit Ω_eff amounts 
already to 927.50 kcal mol^-1.
Averaged over all reference subsets, we find a slope of Δ s(E_disp^D3-BJ) / ΔΩ_eff = 538.76 Bohr^6 kcal mol^-1
with a negligible intercept of 0.16 kcal mol^-1.
Hence, for a molecular system with an effective size of 0.01 Bohr^-6 (see below for an example molecule of approximately this size), the effect of parameter 
uncertainty on the MPU may already be on the order of 5 kcal mol^-1.



We will now extrapolate these findings to a specific example with an effective size that surmounts all 152 systems contained in the original 
reference set.
For this purpose, we selected the C_60 “buckycatcher”, that has been a prototypical example for the development of dispersion corrections<cit.>,
with an effective size of Ω_eff = 1.221 × 10^-2 Bohr^-6 (Fig. <ref>) and a total dispersion energy of E_disp^D3-BJ = -229.27 kcal mol^-1 
as determined with our parameterization provided in Table  <ref>.
Based on the slope and intercept obtained for the entire reference set, we would expect a standard deviation of the total dispersion energy of 
s(E_disp^D3-BJ) = 6.74 kcal mol^-1.
However, when calculating E_disp^D3-BJ for all 10'000 bootstrapped parameter sets, we found a much larger standard deviation of 
s(E_disp^D3-BJ) = 14.08 kcal mol^-1.
As the buckycatcher is a highly unsaturated system, we expect its uncertainty in the total dispersion energy to rather follow the trend exhibited by 
the S22 and S22+ subsets.
With this assumption, we obtained s(E_disp^D3-BJ) = 11.23 kcal mol^-1, which is more in line with the result of our bootstrap analysis.
Noteworthy, the maximum bootstrap spread in E_disp^D3-BJ, max(E_n,b) - min(E_n,b), even amounts to 176.90 kcal mol^-1.


In Section <ref>, we found that the original and uniform weighting schemes yield D3-BJ parameters for PBE that are significantly different from each other.
In the case of the buckycatcher, the deviations in these parameter sets amount to a difference in E_disp^D3-BJ of 
63.47 kcal mol^-1, which is certainly too large to provide a quantitative description of dispersion interactions in that molecular system based on the D3-BJ model.
Contrary to the previous scenario, we find a small difference of 4.95 kcal mol^-1 for the skewed and balanced reference data sets (cf. 
Section <ref>).
These findings indicate that the weighting scheme has a more significant effect on the MPU than small modifications of the reference set.
Furthermore, when comparing E_disp^D3-BJ for the mean parameter sets obtained from the original and modified bootstrapping procedures 
(cf. Section <ref>), we find a difference of only 1.91 kcal mol^-1.

Another factor determining the calibration procedure is the representation of residuals in the cost function.
Most commonly, the absolute value of the residual and the squared residual are employed (cf. Section <ref>). 
While the conclusions drawn in this paper are not affected by the choice of the specific residual representation, it certainly has a numerical impact on 
the dispersion energy since different representations (and, hence, cost functions) generally yield different optimal parameter values.



§ CONCLUSIONS



In this work, we have subjected D3-type semiclassical dispersion corrections to a rigorous statistical analysis, with a focus on D3-BJ.
This study is the first assessing Dx-type semiclassical dispersion interactions based on a rigorous statistical analysis.
We have recalibrated the empirical parameters s_r,6 and s_8 for D3-zero as well as a_1, s_8, and a_2 for D3-BJ.
Overall, our results agree well with the original parameters reported by Grimme and coworkers<cit.>; in some cases, our parameterization results in a 
slightly lower mean prediction error compared to the original one, and in other cases, we arrive at the opposite result.
We therefore do not intend to suggest our new D3-BJ parameterization as an improvement over the original D3-BJ parameterization.

For the molecular systems contained in the reference set, these slightly different parameterizations only have a small effect on the MAE of Eq. (<ref>) 
(less than 1 kcal mol^-1). However, we studied in detail the effect of different parameterizations on increasingly larger molecular systems 
motivated by the fact that physical models based on sums of constantly attractive pairwise interactions are particularly prone to error accumulation.
We examined this issue for a fairly large range of molecular sizes from two angles.
First, we studied the effect of the specific calibration procedure (choice of reference set and cost function, and consideration of uncertainty in the reference data) 
on the dispersion energy.
Second, we also took into account the effect of uncertainty in the empirical D3 parameters on the dispersion energy.
For both, we carried out a rigorous analysis of error accumulation for the three empirical D3-BJ parameters.
While the error introduced by two different parameterizations can partially cancel out, we found that this effect is system-dependent, suggesting that error cancellation 
may occur in some molecules, but be completely absent in others (i.e., every pairwise interaction resulting from one parameterization is constantly larger (or 
smaller) compared to another parameterization).
In that case, the error in the dispersion energy may become dramatically large as it would monotonically increase with the number of pairwise interactions and, hence, 
with molecular size.
Although the correlations of the empirical D3-BJ  parameters obtained from our statistical analysis suggest that partial error compensation is likely to occur,
we found that error cancellation of individual pairwise interaction terms is not possible over wide parameter ranges.

We find that the nonuniform weighting scheme originally employed by Grimme and coworkers<cit.> sometimes leads to larger errors (as assessed through a variety of different validation sets) than those that would result from a uniform weighting scheme that considers all data points on the same footing. 
When considering a relatively large molecular system — here, the prominent C_60 buckycatcher — different weighting schemes lead to deviations in the dispersion energy of more than 60 kcal mol^-1, which is clearly too large for a quantitative description of the total dispersion energy.
Furthermore, we examined the effect of biased reference data on the dispersion energy.
When assuming a small uncertainty in the reference data of 5 %, we find only slightly different parameters, the deviation of which accumulates to
about 5 kcal mol^-1 for the buckycatcher and less than or about 1 kcal mol^-1 for all other molecular systems studied.

Next, we investigated the role of parameter uncertainty for error accumulation.
Here, parameter uncertainty was estimated on the basis of a nonparametric bootstrap analysis.
It turns out that the standard deviations of the means of the empirical parameters are small in all cases.
Moreover, the bootstrap analysis confirms the MAE to be a reliable estimate of prediction uncertainty.
Most interestingly, based on the bootstrap analysis we found a linear dependence of the expected energy uncertainty on the effective molecular size.
For all molecules contained in the reference set, this uncertainty is smaller than 2 kcal mol^-1, while for the buckycatcher, the small estimated parameter 
uncertainty leads to a standard deviation in the dispersion energy of 14 kcal mol^-1 (however, the maximum deviation in the dispersion energy determined with bootstrapping amounts to a surprisingly large value of almost 180 kcal mol^-1
in this case).

We also applied the jackknife method, which is related to the bootstrap approach, to examine the consistency of the reference data.
Our analysis shows that the original reference set employed by Grimme et al. is slightly biased, even though our results also suggest that the effect of this bias on the dispersion energy is negligible in the cases studied.

To support routine DFT-D3 calculations in efficiently estimating uncertainties of dispersion energies, we will make the program BootD3 available on our web page.
BootD3 reads a user-specified file of Cartesian coordinates and runs Grimme's D3 implementation<cit.> for each of the 10'000 bootstrapped parameter sets generated for this study.
In this way, a virtual error bar can be assigned to a dispersion energy.
All density functionals examined in this work and the corresponding bootstrapped parameter sets are available in BootD3.



§ ACKNOWLEDGMENTS


This work was financially supported by the Schweizerischer Nationalfonds (Project No. 200020_169120).




10

phys_rev_1965_140_a1133
Kohn, W.;  Sham, L. J.  Self-Consistent Equations Including Exchange and
  Correlation Effects,  Phys. Rev. 1965, 140,
  A1133–A1138.

kohn1999
Kohn, W.  Nobel Lecture: Electronic structure of matter–wave functions and
  density functionals,  Rev. Mod. Phys. 1999, 71,
  1253–1266.

koch2001
Koch, W.;  Holthausen, M. C. A Chemist's Guide to Density Functional
  Theory; Wiley-VCH Verlag GmbH: Weinheim, Germany, 2nd ed.; 2001.

theory_and_appl_of_comp_chem
Dykstra, C. E.;  Frenking, G.;  Kim, K. S.;  Scuseria, G. E.,  Eds.;
  Theory and Applications of Computational Chemistry: The First Forty
  Years; Elsevier B. V.: Amsterdam, Netherlands, 2005.

j_chem_phys_2014_140_18a301
Becke, A. D.  Perspective: Fifty years of density-functional theory in
  chemical physics,  J. Chem. Phys. 2014, 140,
  18A301.

j_phys_chem_a_2005_109_11015
Ruzsinszky, A.;  Perdew, J. P.;  Csonka, G. I.  Binding Energy Curves from
  Nonempirical Density Functionals II. van der Waals Bonds in Rare-Gas and
  Alkaline-Earth Diatomics,  J. Phys. Chem. A 2005,
  109, 11015–11021.

chem_phys_lett_2006_419_333
Johnson, E. R.;  DiLabio, G. A.  Structure and binding energies in van der
  Waals dimers: Comparison between density functional theory and correlated ab
  initio methods,  Chem. Phys. Lett. 2006, 419,
  333–339.

chem_phys_lett_1994_229_175
Kristyá, S.;  Pulay, P.  Can (semi)local density functional theory
  account for the London dispersion forces?,  Chem. Phys. Lett.
  1994, 229, 175–180.

chem_phys_lett_1995_233_134
Pérez-Jordá, J. M.;  Becke, A. D.  A density-functional study of van
  der Waals forces: rare gas diatomics,  Chem. Phys. Lett.
  1995, 233, 134–137.

wodrich2006
Wodrich, M. D.;  Corminboeuf, C.;  Schleyer, P. v. R.  Systematic
  Errors in Computed Alkane Energies Using B3LYP and Other Popular
  DFT Functionals,  Org. Lett. 2006, 8,
  3631–3634.

j_phys_org_chem_2009_22_1127
Johnson, E. R.;  Mackie, I. D.;  DiLabio, G. A.  Dispersion interactions
  in density-functional theory,  J. Phys. Org. Chem. 2009,
  22, 1127–1135.

j_chem_phys_2012_137_120901
Klimeš, J.;  Michaelides, A.  Perspective: Advances and challenges in
  treating van der Waals dispersion forces in density functional theory,
  J. Chem. Phys. 2012, 137, 120901.

j_phys_condens_matter_2012_24_073201
Dobson, J. F.;  Gould, T.  Calculation of dispersion energies,  J.
  Phys.: Condens. Matter 2012, 24, 073201.

adv_funct_mater_2015_25_2054
Tkatchenko, A.  Current Understanding of Van der Waals Effects in Realistic
  Materials,  Adv. Funct. Mater. 2015, 25,
  2054–2061.

rev_mod_phys_2010_82_1887
French, R. H. et al.   Long range interactions in nanoscale
  science,  Rev. Mod. Phys. 2010, 82, 1887–1944.

chem_rev_2016_116_5105
Grimme, S.;  Hansen, A.;  Brandenburg, J. G.;  Bannwarth, C.
  Dispersion-Corrected Mean-Field Electronic Structure Methods,
  Chem. Rev. 2016, 116, 5105–5154.

steinmann2010
Steinmann, S. N.;  Corminboeuf, C.  A System-Dependent
  Density-Based Dispersion Correction,  J. Chem. Theory Comput.
  2010, 6, 1990–2001.

phys_rev_lett_2004_93_153004
von Lilienfeld, O. A.;  Tavernelli, I.;  Rothlisberger, U.  Optimization
  of Effective Atom Centered Potentials for London Dispersion Forces in Density
  Functional Theory,  Phys. Rev. Lett. 2004, 93,
  153004.

phil_trans_r_soc_a_2014_372_20120476
Peverati, R.;  Truhlar, D. G.  Quest for a universal density functional: the
  accuracy of density functionals across a broad spectrum of databases in
  chemistry and physics,  Phil. Trans. R. Soc. A 2014,
  372, 20120476.

phys_rev_lett_2004_92_246401
Dion, M.;  Rydberg, H.;  Schröder, E.;  Langreth, D. C.;   Lundqvist, B. I.  Van der Waals Density Functional for General Geometries,
  Phys. Rev. Lett. 2004, 92, 246401.

phys_rev_lett_2005_95_109902
Dion, M.;  Rydberg, H.;  Schröder, E.;  Langreth, D. C.;   Lundqvist, B. I.  Erratum: Van der Waals Density Functional for General
  Geometries [Phys. Rev. Lett. 92, 246401 (2004)],  Phys.
  Rev. Lett. 2005, 95, 109902.

vydrov2010
Vydrov, O. A.;  Van Voorhis, T.  Nonlocal van der Waals density functional:
  The simpler the better,  J. Chem. Phys. 2010,
  133, 244103.

j_chem_phys_1974_61_2372
Cohen, J. S.;  Pack, R. T.  Modified statistical method for intermolecular
  potentials. Combining rules for higher van der Waals coefficients,
  J. Chem. Phys. 1974, 61, 2372–2382.

ahlrichs1977
Ahlrichs, R.;  Penco, R.;  Scoles, G.  Intermolecular Forces in Simple
  Systems,  Chem. Phys. 1977, 19, 119–130.

j_chem_phys_2005_123_154101
Becke, A. D.;  Johnson, E. R.  A density-functional model of the dispersion
  interaction,  J. Chem. Phys. 2005, 123, 154101.

j_chem_phys_2005_123_024101
Johnson, E. R.;  Becke, A. D.  A post-Hartree–Fock model of intermolecular
  interactions,  J. Chem. Phys. 2005, 123, 024101.

j_chem_phys_2006_124_174104
Johnson, E. R.;  Becke, A. D.  A post-Hartree–Fock model of intermolecular
  interactions: Inclusion of higher-order corrections,  J. Chem.
  Phys. 2006, 124, 174104.

phys_rev_lett_2009_102_073005
Tkatchenko, A.;  Scheffler, M.  Accurate Molecular Van Der Waals
  Interactions from Ground-State Electron Density and Free-Atom Reference
  Data,  Phys. Rev. Lett. 2009, 102, 073005.

phys_rev_lett_2012_108_236402
Tkatchenko, A.;  Robert A. DiStasio, J.;  Car, R.;  Scheffler, M.
  Accurate and Efficient Method for Many-Body van der Waals Interactions,
  Phys. Rev. Lett. 2012, 108, 236402.

j_comp_chem_2004_25_1463
Grimme, S.  Accurate description of van der Waals complexes by density
  functional theory including empirical corrections,  J. Comput.
  Chem. 2004, 25, 1463–1473.

j_comp_chem_2006_27_1787
Grimme, S.  Semiempirical GGA-type density functional constructed with a
  long-range dispersion correction,  J. Comput. Chem. 2006,
  27, 1787–1799.

grimme2010
Grimme, S.;  Antony, J.;  Ehrlich, S.;  Krieg, H.  A consistent and
  accurate ab initio parametrization of density functional dispersion
  correction (DFT-D)for the 94 elements H-Pu,  J. Chem. Phys.
  2010, 132, 154104.

grimme2011
Grimme, S.;  Ehrlich, S.;  Goerigk, L.  Effect of the damping function in
  dispersion corrected density functional theory,  J. Comput. Chem.
  2011, 32, 1456–1465.

j_chem_phys_2017_147_034112
Caldeweyher, E.;  Bannwarth, C.;  Grimme, S.  Extension of the D3
  dispersion coefficient model,  J. Chem. Phys. 2017,
  147, 034112.

irikura2004
Irikura, K. K.;  Johnson, III, R. D.;  Kacker, R. N.  Uncertainty
  associated with virtual measurements from computational quantum chemistry
  models,  Metrologia 2004, 41, 369–375.

mortensen2005
Mortensen, J. J.;  Kaasbjerg, K.;  Frederiksen, S. L.;   Nørskov, J. K.;  Sethna, J. P.;  Jacobsen, K. W.  Bayesian Error
  Estimation in Density-Functional Theory,  Phys. Rev. Lett.
  2005, 95, 216401.

sutton2016
Sutton, J. E.;  Guo, W.;  Katsoulakis, M. A.;  Vlachos, D. G.  Effects
  of correlated parameters and uncertainty in electronic-structure-based
  chemical kinetic modelling,  Nat. Chem. 2016, 8,
  331–337.

simm2016
Simm, G. N.;  Reiher, M.  Systematic Error Estimation for Chemical Reaction
  Energies,  J. Chem. Theory Comput. 2016, 12,
  2762–2773.

proppe2016
Proppe, J.;  Husch, T.;  Simm, G. N.;  Reiher, M.  Uncertainty
  quantification for quantum chemical models of complex reaction networks,
  Faraday Discuss. 2016, 195, 497–520.

simm2017
Simm, G. N.;  Proppe, J.;  Reiher, M.  Error Assessment of Computational
  Models in Chemistry,  Chimia 2017, 71, 202–208.

proppe2017
Proppe, J.;  Reiher, M.  Reliable Estimation of Prediction Uncertainty for
  Physicochemical Property Models,  J. Chem. Theory Comput.
  2017, 13, 3297–3317.

pernot2015
Pernot, P.;  Civalleri, B.;  Presti, D.;  Savin, A.  Prediction
  Uncertainty of Density Functional Approximations for Properties
  of Crystals with Cubic Symmetry,  J. Phys. Chem. A
  2015, 119, 5288–5304.

pernot2017
Pernot, P.;  Cailliez, F.  A critical review of statistical
  calibration/prediction models handling data inconsistency and model
  inadequacy,  AIChE J. 2017, 63, 4642–4665.

pernot2017a
Pernot, P.  The parameter uncertainty inflation fallacy,  J. Chem.
  Phys. 2017, 147, 104102.

pernot2018
Pernot, P.;  Savin, A.  Probabilistic Performance Estimators for
  Computational Chemistry Methods: The Empirical Cumulative Distribution
  Function of Absolute Errors,   2018,  arXiv:1801.03305
  [physics.chem-ph].

oung2018
Oung, S. W.;  Rudolph, J.;  Jacob, C. R.  Uncertainty quantification in
  theoretical spectroscopy: The structural sensitivity of X-ray emission
  spectra,  Int. J. Quantum Chem. 2018, 118,
  e25458.

ann_statist_1979_7_1
Efron, B.  Bootstrap Methods: Another Look at the Jackknife,  Ann.
  Statist. 1979, 7, 1–26.

j_phys_chem_lett_2016_7_2197
Smith, D. G. A.;  Burns, L. A.;  Patkowski, K.;  Sherrill, C. D.
  Revised Damping Parameters for the D3 Dispersion Correction to Density
  Functional Theory,  J. Phys. Chem. Lett. 2016, 7,
  2197–2203.

dftd3_website
https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3.

j_chem_theory_comput_2010_6_107
Goerigk, L.;  Grimme, S.  A General Database for Main Group Thermochemistry,
  Kinetics, and Noncovalent Interactions – Assessment of Common and
  Reparameterized (meta-)GGA Density Functionals,  J. Chem.
  Theory Comput. 2010, 6, 107–126.

gmtkn_website
https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/GMTKN/gmtkn.

j_chem_theory_comput_2011_7_291
Goerigk, L.;  Grimme, S.  Efficient and Accurate Double-Hybrid-Meta-GGA
  Density Functionals–Evaluation with the Extended GMTKN30 Database for
  General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions,
  J. Chem. Theory Comput. 2011, 7, 291–309.

proppe2017a
Proppe, J.;  Reiher, M. “reBoot: A Program for Statistical
  Calibration of Property Models, www.reiher.ethz.ch/software/reboot.html”,
  2017.

comput_sci_eng_2011_13_22
van der Walt, S.;  Colbert, S. C.;  Varoquaux, G.  The NumPy Array: A
  Structure for Efficient Numerical Computation,  Comput. Sci. Eng.
  2011, 13, 22–30.

scipy
Jones, E. et al.  “SciPy: Open source scientific tools for
  Python”,  2001.

j_global_optim_197_11_341
Storn, R.;  Price, K.  Differential Evolution – A Simple and Efficient
  Heuristic for Global Optimization over Continuous Spaces,  J. Global
  Optim. 1997, 11, 341–359.

mathematica11_0_1_0
Wolfram Research, Inc., “Mathematica, Version 11.0.1.0”,  2016.

riu2003
Riu, J.;  Bro, R.  Jack-Knife Technique for Outlier Detection and Estimation
  of Standard Errors in PARAFAC Models,  Chemometr. Intell. Lab.
  2003, 65, 35–49.

phys_chem_chem_phys_2008_10_2813
Zhao, Y.;  Truhlar, D. G.  Computational characterization and modeling of
  buckyball tweezers: density functional study of concave–convex
  π⋯π interactions,  Phys. Chem. Chem. Phys.
  2008, 10, 2813–2818.
















