Gaussian Process-Based Refinement of Dispersion Corrections

TL;DR

This paper introduces a Gaussian process regression method to refine dispersion corrections in molecular systems, improving accuracy through active learning and variance-based data selection.

Contribution

The study develops a novel GP-based refinement model for dispersion corrections, utilizing active learning and batch-wise variance sampling for efficient training.

Findings

The D3-GP model improves dispersion correction accuracy.

Variance-based sampling reduces the number of reference calculations needed.

The approach is applicable to diverse molecular conformations.

Abstract

We employ Gaussian process (GP) regression to adjust for systematic errors in D3-type dispersion corrections introducing the associated, statistically improved model D3-GP. We generated a data set containing interaction energies for 1,248 molecular dimers, which resemble the dispersion-dominated systems contained in the S66 data set. Our systems do not only represent equilibrium structures, but also dimers with various relative orientations and conformations at both shorter and longer distances. To train our D3-GP model, we engineered two different vectorial representations of (supra-)molecular systems, both derived from the matrix of atom-pairwise D3(BJ) interaction terms: (a) a distance-resolved interaction energy histogram and (b) eigenvalues of the interaction matrix ordered according to their decreasing absolute value. We harness the prediction variance obtained from GP regression to select optimal training sets in an automated fashion. The larger the variance, the more information the corresponding data point may add to the training set. For a given set of molecular systems, variance-based sampling can approximately determine the smallest subset being subjected to reference calculations such that all dispersion corrections for the remaining systems fall below a predefined accuracy threshold. Our refined learning algorithm selects multiple (instead of single) systems which can be subjected to reference calculations simultaneously. We refer to the underlying selection strategy as batch-wise variance-based sampling (BVS). BVS-guided active learning is an essential component of our D3-GP workflow, which is implemented in a black-box fashion. This approach leads overall to a self-improving model (D3-GP) that predicts system-focused and GP-refined D3-type dispersion corrections for any given system of reference data.

Tables

Table 1. Table 1: Overview of the 1,248 molecular reference systems (dimers) studied in this work.

Data set	#Dimers	Description
S13x8	104	dispersion-dominated subset (dimers #34 to #46) of the S66x8 data set[57]
ROTA	1,100	ethyne–pentane dimers; varying relative orientations and centroid distances (3.5–10.0 Å), fixed pentane conformation
CONF	44	ethyne–pentane dimers; varying relative orientations and pentane conformations, fixed centroid distance (5.2 Å)

Table 2. Table 2: Performances (as measured by MARE and MAE) of four different PBE-specific D3(BJ) parameter sets with respect to three validation sets ( S13x8-V , ROTA , and CONF ). The highest accuracies are highlighted in bold. The D3(BJ) parameters a 1 subscript 𝑎 1 a_{1} and s 8 subscript 𝑠 8 s_{8} are dimensionless. Values for the D3(BJ) parameter a 2 subscript 𝑎 2 a_{2} and the MAE are reported in bohr and kcal mol -1 , respectively.

					S13x8-V			ROTA			CONF
Reference	$a_1$	$a_2$	$s_8$		MARE	MAE		MARE	MAE		MARE	MAE
Grimme et al.[25]	0.4289	4.4407	0.7875		13.4 %	0.27		27.8 %	0.10		34.7 %	0.21
Smith et al.[28]	0.0121	5.9390	0.3589		9.5 %	0.19		23.8 %	0.08		29.7 %	0.18
Weymuth et al.[30]	0.4309	4.8327	1.0892		12.0 %	0.26		25.6 %	0.09		30.8 %	0.19
This work	0.0000	5.6841	0.0000		8.2 %	0.19		20.1 %	0.07		23.8 %	0.15

Equations

$$\mathbb{E}\big{[}\hat{y}(\mathbf{x}_{\ast})\big{]}\equiv\hat{\mu}_{y}(\mathbf{x}_{\ast})=\mathbf{k}(\mathbf{x}_{\ast},\mathbf{X})\Big{(}\mathbf{K}(\mathbf{X},\mathbf{X})+\alpha_{0}\mathbf{I}\Big{)}^{-1}\mathbf{y}\$$

$$\mathbb{V}\big{[}\hat{y}(\mathbf{x}_{\ast})\big{]}\equiv\hat{\sigma}_{y}^{2}(\mathbf{x}_{\ast})=k(\mathbf{x}_{\ast},\mathbf{x}_{\ast})-\mathbf{k}(\mathbf{x}_{\ast},\mathbf{X})\Big{(}\mathbf{K}(\mathbf{X},\mathbf{X})+\alpha_{0}\mathbf{I}\Big{)}^{-1}\mathbf{k}(\mathbf{X},\mathbf{x}_{\ast})\ ,$$

$$E_{\text{elec},\mathcal{S}}^{\text{DLPNO-CCSD(T)/CBS}}\approx \frac{E_{\text{elec},\mathcal{S}}^{\text{HF/QZ}}\exp \{1.63X\}-E_{\text{elec},\mathcal{S}}^{\text{HF/TZ}}\exp \{1.63(X-1)\}}{\exp \{1.63X\}-\exp \{1.63(X-1)\}}$$

$$+\frac{E_{\text{elec},\mathcal{S}}^{\text{corr/QZ}}{X}^3-E_{\text{elec},\mathcal{S}}^{\text{corr/TZ}}(X-1{)}^3}{X^3-(X-1{)}^3}.$$

$$E_{\text{int}}^{\mathcal{M}}\equiv\Delta E_{\text{elec}}^{\mathcal{M}}=E_{\text{elec},\mathcal{AB}}^{\mathcal{M}}-\Big{(}E_{\text{elec},\mathcal{A}}^{\mathcal{M}}+E_{\text{elec},\mathcal{B}}^{\mathcal{M}}\Big{)}\ ,$$

$$E_{\text{disp},\mathcal{S}}^{\text{D3(BJ)}}\equiv E^{\text{D3(BJ)}}=\sum _{I>J}\sum _{n=6,8}{s}_n\frac{C_n^{IJ}}{R_{IJ}^n+{\left(a_1\sqrt{C_8^{IJ}/C_6^{IJ}}+a_2\right)}^n}.$$

$$\Delta E_{\text{disp}}^{\text{D3(BJ)}}=E_{\text{disp},\mathcal{AB}}^{\text{D3(BJ)}}-\Big{(}E_{\text{disp},\mathcal{A}}^{\text{D3(BJ)}}+E_{\text{disp},\mathcal{B}}^{\text{D3(BJ)}}\Big{)}\ .$$

$$\Delta E_{\text{int}}=E_{\text{int}}^{\text{DLPNO-CCSD(T)/CBS}}-\left(E_{\text{int}}^{\text{PBE/ma-def2-QZVPP}}+\Delta E_{\text{disp}}^{\text{D3(BJ)}}\right),$$

$$E_{IJ}^{\text{D3(BJ)}}=\{\begin{array}{ll}{\sum }_{n=6,8}{s}_n\frac{C_n^{IJ}}{R_{IJ}^n+{\left(a_1\sqrt{C_8^{IJ}/C_6^{IJ}}+a_2\right)}^n}& \forall I\ne J\\ 0& \forall I=J\end{array}.$$

$$e_m^{histD3(\mathsf{BJ})}=\sum _{\begin{array}{c}I>J\\ r_m^{\text{min}}<R_{IJ}\le r_m^{\text{max}}\end{array}}\sum _{n=6,8}{s}_n\frac{C_n^{IJ}}{R_{IJ}^n+{\left(a_1\sqrt{C_8^{IJ}/C_6^{IJ}}+a_2\right)}^n}.$$

$$\mathbf{x}^{\mathsf{eig/histD3(BJ)}}\equiv\Delta\mathbf{e}^{\mathsf{eig/histD3(BJ)}}=\mathbf{e}_{\mathcal{AB}}^{\mathsf{eig/histD3(BJ)}}-\Big{(}\mathbf{e}_{\mathcal{A}}^{\mathsf{eig/histD3(BJ)}}+\mathbf{e}_{\mathcal{B}}^{\mathsf{eig/histD3(BJ)}}\Big{)}\ .$$

$$k_{\text{iso}}^{\text{M12}}({\mathbf{x}}_i,{\mathbf{x}}_j)={\alpha }_1\exp \left\{-d_{\text{iso}}({\mathbf{x}}_i,{\mathbf{x}}_j)\right\},$$

$$d_{\text{iso}}({\mathbf{x}}_i,{\mathbf{x}}_j)=\sqrt{\frac{({\mathbf{x}}_i-{\mathbf{x}}_j)({\mathbf{x}}_i-{\mathbf{x}}_j{)}^{\top }}{{\alpha }_2}},$$

$$k_{\text{iso}}^{\text{M32}}({\mathbf{x}}_i,{\mathbf{x}}_j)={\alpha }_1\left(1+\sqrt{3}{d}_{\text{iso}}({\mathbf{x}}_i,{\mathbf{x}}_j)\right)\exp \left\{-\sqrt{3}{d}_{\text{iso}}({\mathbf{x}}_i,{\mathbf{x}}_j)\right\}.$$

$$k_{\text{aniso}}^{\text{M32}}({\mathbf{x}}_i,{\mathbf{x}}_j)={\alpha }_1\left(1+\sqrt{3}{d}_{\text{aniso}}({\mathbf{x}}_i,{\mathbf{x}}_j)\right)\exp \left\{-\sqrt{3}{d}_{\text{aniso}}({\mathbf{x}}_i,{\mathbf{x}}_j)\right\},$$

$$d_{\text{aniso}}({\mathbf{x}}_i,{\mathbf{x}}_j)=\sqrt{({\mathbf{x}}_i-{\mathbf{x}}_j)\left(\begin{array}{ccc}{\alpha }_{2,a_1}& 0& 0\\ 0& {\alpha }_{2,a_2}& 0\\ 0& 0& {\alpha }_{2,s_8}\end{array}\right)({\mathbf{x}}_i-{\mathbf{x}}_j{)}^{\top }}.$$

$$u_{\text{LCB}}(\mathbf{x})={\hat{\mu }}_y(\mathbf{x})-{\hat{\sigma }}_y(\mathbf{x}).$$

$$\hat{\mu }\equiv \text{ME}=\frac{1}{N}\sum _{n=1}^N{\Delta }_n.$$

$$\hat{\sigma }\equiv \text{RMSD}=\sqrt{\frac{1}{N-1}\sum _{n=1}^N{\left({\Delta }_n-\hat{\mu }\right)}^2}.$$

$$\text{RMSE}=\sqrt{\frac{1}{N-1}\sum _{n=1}^N{\Delta }_n^2}=\sqrt{\frac{N}{N-1}{\hat{\mu }}^2+{\hat{\sigma }}^2}.$$

$$\text{MAE}=\frac{1}{N}\sum _{n=1}^N\left|{\Delta }_n\right|,$$

$${\delta }_n=\frac{{\Delta }_n}{\left|y_n\right|}.$$

$$\text{MRE}=\frac{1}{N}\sum _{n=1}^N{\delta }_n,$$

$$\text{MARE}=\frac{1}{N}\sum _{n=1}^N\left|{\delta }_n\right|.$$

$${\delta }_n^c=\frac{{\Delta }_n}{\max \left\{\left|y_n\right|,c\right\}},$$

$$E^{\text{D3(BJ)}}=\sum _{I>J}\frac{C_6^{IJ}}{R_{IJ}^6+a_2^6}.$$

$${\rho }_{\text{inv}}=\sum _{I\in \mathcal{A}}\sum _{J\in \mathcal{B}}{R}_{IJ}^{-1}.$$

Full text

Keywords: noncovalent interactions, dispersion corrections, Gaussian process regression, active learning, uncertainty quantification.

aDepartments of Chemistry and Computer Science, University of Toronto, Toronto, Ontario, Canada.

bLaboratory of Physical Chemistry, ETH Zurich, Vladimir-Prelog-Weg 2, 8093 Zurich, Switzerland.

† These authors contributed equally.

1 Introduction

The inability of (semi)local exchange–correlation functionals to capture long-range electron correlation and, hence, London dispersion interactions, was demonstrated several times in the past two decades.[1, 2, 3, 4, 5, 6] Three major groups of approaches were developed to account for dispersion interactions in density functional theory (DFT): effective one-electron potentials, nonlocal density functionals, and semiclassical dispersion corrections. Important examples of the first two groups are the atom-centered nonlocal potentials by Röthlisberger and co-workers,[7] and the Vydrov–Van Voorhis nonlocal density functional.[8] A discussion of these nonadditive approaches is beyond the scope of this paper, and we refer the reader to Refs.  9 and  10 for general reviews addressing dispersion corrections for electronic-structure methods. Here, we consider semiclassical dispersion corrections, which are simply added to the electronic energy. The computational efficiency resulting from this additive character indicates why the semiclassical approach is the most widely used approach to take dispersion interactions into account. It has been proposed as early as the 1970s.[11, 12] Significant advances were achieved in the 2000s and 2010s by Becke and Johnson,[13, 14, 15] Tkatchenko and Scheffler,[16, 17] Steinmann and Corminboeuf,[18, 19] Sato and Nakai,[20, 21] and Grimme and co-workers[22, 23, 24, 25, 26, 27] (D$x$-type dispersion corrections with $x=1,2,3,4$).

The popularity of D$x$-type dispersion correction models[22, 23, 24, 25, 26, 27] may be attributed to (i) their good performance across the periodic table, (ii) their availability for many mean-field electronic-structure approximations, (iii) their incorporation in various quantum-chemical computer programs, and (iv) their very low computational complexity. Originally optimized for equilibrium structures, Smith et al.[28] showed that a mere reparametrization of the most popular D$x$ models, D3(0)[24] (with the original, zero-damping function) and D3(BJ)[25] (with the Becke–Johnson damping function), also allows for an accurate description of noncovalent interactions across configuration space. Further improvements over the original D3 approach have recently been released with D4.[26, 27] Compared to D3-type dispersion corrections, the estimation of dispersion coefficients is not only based on the molecular structure, but also involves partial charges obtained from a classical model (default) or, alternatively, from electronic-structure calculations. In the latter case, the dependence on the electron density renders the D4 approach more similar to other semiclassical dispersion correction schemes that also incorporate electronic-structure information.[13, 14, 15, 16, 17, 18, 19, 20, 21]

Predictions derived from truncated physical models (such as those discussed in this work) can be expected to be biased in a nontrivial way with respect to the underlying model variables.[29] To identify such systematic errors, one could, in principle, carry out reference measurements or reference calculations in a case-by-case fashion. However, this approach is unfeasible and would render any prediction unnecessary. One would rather prefer a measure that informs about the reliability of a prediction, here, an a posteriori correction of the approximate interaction energy of a given (supra-)molecular system. The variance of a system-specific dispersion correction obtained from a statistical analysis is one possible measure. Recently, we showed how to estimate such variances by sampling the uncertainty in the D3(BJ) parameters with bootstrapping.[30] However, parameter uncertainty may explain only a small fraction of the total prediction error, which is usually dominated by the inadequacy of the model itself.[31, 32]

Machine learning[33] offers a remedy for alleviating model inadequacy as the performance of machine-learned models improves with data. Only a few studies have been published to date which probe the possibilities of machine learning for noncovalent-interaction research. Mezei et al.[34] applied kernel ridge regression[35] in combination with many-body descriptors[36] to learn corrections for interaction energies obtained from DFT calculations. Instead of learning corrections to the interaction energy directly, they learned differences to CCSD(T) atomization energies to also improve on intramolecular interactions. McGibbon et al.[37] learned corrections for interaction energies obtained from second-order Møller–Plesset perturbation theory calculations. They applied neural-network regression[33] on descriptors that were derived from various electronic-structure approximations, e.g., density matrix overlaps and components of noncovalent interactions (dispersion, electrostatics, exchange, and induction). CCSD(T) interaction energies for 200,000 molecules across the periodic table were employed to train their neural network. They estimated system-specific prediction uncertainties by generating 1,000 instances of their neural network via dropout.[38] Gao et al.[39] used partial least-squares regression[40] to select descriptors such as electronic energies or HOMO–LUMO gaps, which were obtained from DFT calculations. A neural network was trained on these descriptors to learn corrections to DFT interaction energies with respect to a CCSD(T) reference. An alternative approach to learning corrections for interaction energies was proposed by Bereau et al.[41] They applied kernel ridge regression to estimate local environment-dependent properties such as electrostatic multipole coefficients and polarizabilities. These machine-learned quantities were then fed into classical intermolecular potentials to improve the prediction of many-body dispersion and other types of noncovalent interactions. Machine-learned models for electron densities[42, 43] have been adapted for highly scalable models[44] and to gain insights into noncovalent interaction[45] by Corminboeuf, Ceriotti, and co-workers. Machine learning usually is a data-extensive endeavor as the corresponding prediction models potentially improve with an increasing number of benchmark data. However, the calculation of benchmark interaction energies constitutes a computational bottleneck. This limitation offers the incentive to generate informative training sets, which contain a minimum number of entries to achieve a predefined accuracy over a given domain of application, i.e., the chemical and configuration space covered by the set of systems at hand. It is not sufficient to sample the existing set of systems uniformly, as the number of systems that are required to achieve that accuracy is not known a priori. Therefore, a dynamic extension of the training set is necessary, e.g., by means of active learning.[46] In active-learning approaches, predictions learned from a given training set are harnessed to determine a batch of new systems for which benchmarks are required. The training set will be extended by these benchmarks to optimally improve the underlying prediction model for those systems that are not yet contained in the training set.

Regarding research in computational chemistry, controlling the number of benchmarks is particularly important in data-extensive scenarios such as molecular-dynamics simulations and chemical-space exploration. In such scenarios, not even the set of systems is known a priori. Botu and Ramprasad[47] compared similarities between molecular fingerprints[48] (vectorial representations of molecules) to actively learn energies and forces for molecular-dynamics simulations. If a new configuration along the trajectory is too dissimilar to the training set, electronic-structure calculations are performed to update the machine-learned force field. Another approach was described by Behler[49] to select training data for learning high-dimensional neural-network potentials. Several neural networks with different architectures but similar performances are employed to predict potential energies for a given trajectory obtained from one of these networks. Due to the choice of multiple neural networks, a prediction variance can be estimated for each configuration along the trajectory. High variances indicate a disagreement between the different neural-network models. Therefore, high-variance configurations are subjected to electronic-structure calculations to update these models. This approach is sometimes denoted “Query by Committee”.[50] Active-learning approaches have also been applied to accelerate chemical-space exploration (based on Thompson sampling[51]) and for the prediction of molecular properties (based on the D-optimality criterion[52]). The seminal work by Csányi et al. can be considered an early predecessor to the active-learning approach developed here.[53]

In this work, we adjust for systematic errors in D3(BJ) dispersion corrections by exploiting the capabilities of Gaussian process (GP) regression,[54] the Bayesian analog to kernel ridge regression. We refer to the resulting prediction model as D3(BJ)-GP. The only input it requires is the structure of a (supra-)molecular system, i.e., its atomic positions and corresponding proton numbers, which are transformed to linear combinations of atom-pairwise D3(BJ) interaction terms. This way, the computational efficiency of the original D3(BJ) approach can be preserved to a high degree. Due to the Bayesian nature of GP regression, D3(BJ)-GP dispersion corrections are being sampled from a posterior probability distribution. The variance of this distribution can be considered a reliability measure. If the variance is small with respect to a user-defined tolerance, the mean of the posterior distribution can be considered a reliable dispersion correction. If the variance is too large, however, the user is advised to perform reference electronic-structure calculations in order to generate a new benchmark, extend the training set, and update the GP posterior distribution. Here, we harness the posterior variance as a simple selection criterion to actively learn dispersion corrections. Hence, the training set employed to optimize the D3(BJ)-GP model is dynamic, which constitutes an essential novelty to the active research field of dispersion-corrected electronic-structure methods.[9, 10] As a result, a system-focused model is obtained that by virtue of automated reference calculations can improve itself in a rolling fashion. Naturally, a meta machine-learning model may be applied to extract transferable knowledge from sets of such system-focused self-improving models. Furthermore, we present an extension to our original, sequential formulation of variance-based active learning,[55] which allows us to generate multiple benchmarks in parallel. We refer to our strategy of selecting multiple molecular systems prior to benchmark generation as batch-wise variance-based sampling (BVS). BVS is particularly interesting for data-extensive scenarios as mentioned above.

This paper is structured as follows: A concise introduction to GP regression and related topics (Bayesian optimization and BVS) is provided in Section 2. In Section 3, we introduce the data sets studied in this work, provide details on electronic-structure calculations, and describe the settings employed for GP regression and subsequent analyses. The results of this study are presented and discussed in Section 4. We report on a problem-specific D3(BJ) reparametrization, the performance of the D3(BJ)-GP model based on a fixed training set, and the improvement of the D3(BJ)-GP performance resulting from BVS-guided active learning.

2 Gaussian Process Regression

The GP posterior distribution of the target variable $y$, which may be a physical property, is Gaussian by definition and can, therefore, be fully specified by both its mean and its variance. As with any other regression method, GP regression is applied to learn a mapping from a domain $\{x_{1},\cdots,x_{M}\}$ of $M$ continuous input variables (features) to $y$. The features may be descriptors of a physical system. Given $N$ observations of the target variable, $\mathbf{y}=(y_{1},\cdots,y_{N})^{\top}$, mean and variance mapped from any feature vector $\mathbf{x}_{\ast}=(x_{\ast,1},\cdots,x_{\ast,M})$ read

[TABLE]

and

[TABLE]

respectively. The rows of the design matrix $\mathbf{X}$ represent the $N$ feature vectors $\{\mathbf{x}_{n}=(x_{n,1},\cdots,x_{n,M})\}$ of the observed (training) set. The kernel $k(\mathbf{x}_{i},\mathbf{x}_{j})$ is a measure of similarity between feature vectors $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$, and the row vector $\mathbf{k}(\mathbf{x}_{i},\mathbf{X})$ as well as the matrix $\mathbf{K}(\mathbf{X},\mathbf{X})$ are the similarity between one feature vector, $\mathbf{x}_{i}$, or a matrix of features, $\mathbf{X}$, with the design matrix $\mathbf{X}$. The hyperparameter $\alpha_{0}$, multiplied with the $(N\times N)$-dimensional identity matrix $\mathbf{I}$, represents the constant (homoscedastic) noise of each observation. We refer the reader to Ref. 54 for a comprehensive introduction to the topic.

2.1 Bayesian Optimization

The mean of the GP posterior distribution, $\hat{\mu}_{y}(\mathbf{x})$, represents the best approximation to the target function $y(\mathbf{x})$. As evaluations (observations) of $\hat{\mu}_{y}(\mathbf{x})$ are much faster than evaluations of $y(\mathbf{x})$, the former can also be used to accelerate searches on $y(\mathbf{x})$ if one is interested in optimization tasks. In such cases, $y(\mathbf{x})$ serves as an objective function. To ensure finding the global minimum (or other extrema) of $y(\mathbf{x})$, we target an optimal balance between probing its already known, promising regions and exploring its yet unknown regions. Bayesian optimization[56] achieves this goal by combining mean and variance of the GP posterior distribution, $\hat{\mu}_{y}(\mathbf{x})$ and $\hat{\sigma}^{2}_{y}(\mathbf{x})$, to construct an acquisition function ${u(\mathbf{x})=f\left(\hat{\mu}_{y}(\mathbf{x}),\hat{\sigma}^{2}_{y}(\mathbf{x})\right)}$, which can be interpreted as a purposely biased objective function. The true objective $y(\mathbf{x}_{\ast})$ is evaluated at the current minimum $\mathbf{x}_{\ast}$ of $u(\mathbf{x})$, which, in turn, is used to update the posterior distribution.

2.2 Variance-Based Active Learning

Consider a training set $\mathcal{T}^{N_{y}}_{N_{\mathbf{x}}}$ with $N_{\mathbf{x}}$ feature vectors and $N_{y}$ observations, where $N_{\mathbf{x}}=N_{y}$, and a query set $\mathcal{Q}_{M_{\mathbf{x}}}$ with $M_{\mathbf{x}}$ feature vectors. We seek to add the smallest subset of $\mathcal{Q}_{M_{\mathbf{x}}}$ to $\mathcal{T}^{N_{y}}_{N_{\mathbf{x}}}$ such that every feature vector of the remaining query set is assigned a posterior standard deviation $\hat{\sigma}_{y}(\mathbf{x}_{\ast})$ that is below a predefined accuracy threshold $t$. Promoting the query vector with the highest posterior variance to a training vector is a straightforward strategy in this respect. Note that the hyperparameters are a function of the training set defined by $\mathbf{X}$ and $\mathbf{y}$ and, hence, the variance $\hat{\sigma}_{y}^{2}$ as defined in Eq. (2) depends implicitly on $\mathbf{y}$. As a consequence, updating $\hat{\sigma}_{y}^{2}$ requires an observation after each promotion of a query vector. This sequential approach[55] to making observations constitutes a critical computational bottleneck.

To render our approach more efficient, we determine a batch size $L$ that reflects the number of observations we can actually make in parallel. We determine the query vector with the highest posterior variance, $\hat{\sigma}_{y}^{2}$, and add its feature vector to $\mathbf{X}$, resulting in $\mathcal{Q}_{M_{\mathbf{x}}-1}$ and $\mathcal{T}_{N_{\mathbf{x}}+1}^{N_{y}}$. Subsequently, we update $\hat{\sigma}_{y}^{2}$ without a preceding hyperparameter optimization. As we temporarily fix the hyperparameters, $\hat{\sigma}_{y}^{2}$ does no longer depend on $\mathbf{y}$ and, therefore, we can update the variance without making observations. We perform this subalgorithm (query vector promotion followed by a fixed-hyperparameter variance update) at most $L$ times to obtain $\mathcal{Q}_{M_{\mathbf{x}}-l}$ and $\mathcal{T}_{N_{\mathbf{x}}+l}^{N_{y}}$, where $l\in\{1,\cdots,L\}$. The subalgorithm will stop if either $\hat{\sigma}_{y}(\mathbf{x}_{m})<t~{}\forall\mathbf{x}_{m}\in\mathcal{Q}_{M_{\mathbf{x}}-l}$ or $l=L$. After adding a batch of $l$ feature vectors from the query set to the training set, we make $l$ observations in parallel to obtain $\mathcal{T}_{N_{\mathbf{x}}+l}^{N_{y}+l}$. If either $\hat{\sigma}_{y}(\mathbf{x}_{m})<t~{}\forall\mathbf{x}_{m}\in\mathcal{Q}_{M_{\mathbf{x}}-l}$ or the query set is empty, the overall algorithm stops. Otherwise, we sample another batch of $l\leq L$ query vectors as described above.

3 Methodology

3.1 Molecular Reference Systems

We considered three different data sets (Table 1). The first set, S13x8, consists of 13 chemically distinct dimers at 8 different intermolecular distances each (104 dimers in total). It is the dispersion-dominated subset (dimers #34 to #46) of the S66x8 data set,[57] which can be retrieved from the BEGDB database.[58] The second set, ROTA, contains 1,100 ethyne–pentane dimers with 16 different centroid distances (3.5–10.0 Å). The relative orientation between ethyne and pentane was sampled from a uniform spherical distribution (see the Supporting Information for further details). All of these dimers contain the same pentane conformer (the one of the ethyne–pentane dimer at equilibrium). The third set, CONF, contains 44 ethyne–pentane dimers with a fixed centroid distance of 5.2 Å. Here, random conformations of pentane were sampled with RDKit[59] and the relative orientation between ethyne and pentane was, again, sampled from a uniform spherical distribution. The latter two sets, ROTA and CONF, were created for this work. Both molecular structures and interaction energies of all 1,248 dimers considered are provided in the Supporting Information.

3.2 Calculation of Interaction Energies

3.2.1 Electronic-Structure Calculations

All quantum chemical calculations were carried out with Orca 4.0.1.[60, 61] For all chemical systems studied, we performed DLPNO-CCSD(T)[62, 63] and PBE calculations[64] on both the interacting and the isolated monomers. In the DLPNO-CCDS(T) calculations, we employed the aug-cc-pV$n$Z[65, 66] basis sets ($n=$ T, Q) and aug-cc-pV$m$Z[67, 68] auxiliary basis sets ($m=n+1$, i.e., $m=$ Q, 5). We chose the TightPNO keyword, which is recommended for studying noncovalent interactions.[69] We extrapolated the triple-$\zeta$ and quadruple-$\zeta$ DLPNO-CCSD(T) energies to the complete-basis-set (CBS) limit with the extrapolation scheme developed by Halkier et al.[70, 71] that has already been applied for the calculation of DLPNO-CCSD(T)/CBS energies.[72, 73, 74] The scheme comprises two separate extrapolations for the Hartree–Fock (HF) and correlation (corr) fractions of the electronic energy,

[TABLE]

Here, the subscript $\mathcal{S}$ denotes the molecular system and $X$ is the cardinal number of the larger basis set, i.e., $X=n_{\text{max}}$ (here, $n_{\text{max}}=4$). For our PBE calculations, we chose the ma-def-QZVPP[75, 76] basis set and the def2-QZVP[68] auxiliary basis set. We applied counterpoise (CP) corrections[77] and TightSCF criteria to both DLPNO-CCSD(T) and PBE calculations. The intermolecular interaction energy was obtained as the difference between the electronic energies of the isolated monomers $\mathcal{A}$ and $\mathcal{B}$, and the corresponding dimer $\mathcal{AB}$,

[TABLE]

where $\mathcal{M}:=\{$DLPNO-CCSD(T)/CBS, PBE/ma-def2-QZVPP$\}$.

3.2.2 D3-Type Dispersion Corrections

All D3-type dispersion corrections to PBE interaction energies were obtained from Grimme’s stand-alone program dftD3[24, 25] and incorporate the Becke–Johnson (BJ) damping scheme,[25]

[TABLE]

Here, $R_{IJ}$ is the distance between atoms $I$ and $J$ in system $\mathcal{S}$, and $C_{6}^{IJ}$ and $C_{8}^{IJ}$ are the corresponding dispersion coefficients taken from Refs. 24 and 25, which are negative by definition. Spatial position and proton number of the $I$-th atom are denoted $\mathbf{R}_{I}$ and $Z_{I}$, respectively. The scaling factor $s_{6}$ is set to unity by definition,[25] whereas $a_{1}$, $a_{2}$, and $s_{8}$ are functional-dependent, empirical parameters. We employed several parametrizations of the D3(BJ) model;[25, 28, 30] see also Section 3.3.2. D3(BJ) dispersion corrections were calculated for both interacting and isolated monomers to yield relative D3(BJ) dispersion corrections,

[TABLE]

3.3 Gaussian Process Regression

We employed the gpml[78] code written in Matlab[79] for all GP regression and prediction tasks discussed in this work. Optimal hyperparameters for the GP models considered were obtained with Rasmussen’s implementation of the conjugate gradient algorithm[78] in gpml. For estimating the GP posterior distribution, we employed a zero-mean prior and a Gaussian likelihood.

3.3.1 D3-GP-Type Dispersion Corrections

The residual (intermolecular) interaction energy,

[TABLE]

served as the target variable. We considered two different featurizations, eigD3(BJ) and histD3(BJ), both being derived from the D3(BJ) matrix $\mathbf{E}^{\text{D3(BJ)}}$ with elements

[TABLE]

To obtain the eigD3(BJ) feature vector, the eigenvalues of $\mathbf{E}^{\text{D3(BJ)}}$ were sorted according to decreasing absolute value. As the dimensionality of the eigD3(BJ) vector represents the number of atoms in the corresponding molecular system, eigD3(BJ) vectors of smaller systems were padded with zeros to match the dimensionality of eigD3(BJ) vectors corresponding to larger systems. This procedure has been adopted from Rupp et al. who generated feature vectors from Coulomb matrices.[80]

By contrast, the histD3(BJ) feature vector is strictly 16-dimensional where each element represents a sum over $E_{IJ}^{\text{D3(BJ)}}$ corresponding to a predefined range of interatomic distances (Table S1),

[TABLE]

Both eigD3(BJ) and histD3(BJ) vectors were determined for the interacting and isolated monomers. The final feature vectors employed in GP regression read

[TABLE]

We considered both the isotropic Matérn-$\nicefrac{{1}}{{2}}$ kernel (also known as exponential kernel or Laplacian kernel),[54]

[TABLE]

[TABLE]

and the isotropic Matérn-$\nicefrac{{3}}{{2}}$ kernel,[54]

[TABLE]

As the numerical uncertainty in the interaction energies we consider can be controlled by the SCF convergence criterion, we have prior knowledge of the noise hyperparameter, $\alpha_{0}$. Therefore, we exclude $\alpha_{0}$ from hyperparameter optimization and set it to a constant value of $1.0\times 10^{-5}$ kcal mol*-1*. This way, hyperparameter optimization can be accelerated and is less prone to overfitting. We optimized the remaining hyperparameters $\alpha_{1}$ and $\alpha_{2}$ by minimizing the leave-one-out negative log marginal lihood (maximum of 1,000 optimization steps).

3.3.2 D3(BJ) Reparametrization

We applied Bayesian optimization[56] to reparametrize the D3(BJ) model for PBE. We chose the mean absolute relative error (MARE, Eq. (23)), of $\Delta E_{\text{int}}$ over a subset of S13x8 (S13x8-T, see the Supporting Information), as objective function, $y(\mathbf{x})$. Here, $\mathbf{x}$ is a 3-dimensional vector composed of the empirical D3(BJ) parameters $a_{1}$, $a_{2}$, and $s_{8}$. 25 parameter vectors (initial training set) were generated by Latin hypercube sampling[81] over the domain $\mathcal{X}=\mathcal{X}_{a_{1}}\times\mathcal{X}_{a_{2}}\times\mathcal{X}_{s_{8}}$, where ${\mathcal{X}_{a_{1}}:=[0.0,0.7]}$, ${\mathcal{X}_{a_{2}}:=[4.0,6.5]}$ bohr, and ${\mathcal{X}_{s_{8}}:=[0.0,3.5]}$. We rounded each value to the fourth decimal place according to previous D3 parametrizations. We evaluated $y(\mathbf{x})$ for each of these 25 parameter triples to approximate the objective function over the same domain by GP regression with the anisotropic Matérn-$\nicefrac{{3}}{{2}}$ kernel,[54]

[TABLE]

[TABLE]

Here, $\alpha_{1}$, $\alpha_{2,a_{1}}$, $\alpha_{2,a_{2}}$, and $\alpha_{2,s_{8}}$ represent four hyperparameters in addition to the noise hyperparameter $\alpha_{0}$ introduced in Eqs. (1) and (2). We optimized the five hyperparameters by minimizing the log marginal likelihood (maximum of 10,000 optimization steps). We constructed a lower-confidence-bound acquisition function,

[TABLE]

At each iteration, we defined an equidistant grid with a resolution of 0.1 units over the domain $\mathcal{X}$. The initial and terminal grid points of each dimension were located at the bounds of the domain. The grid point that minimized $u_{\text{LCB}}(\mathbf{x})$ defined the center of a cube with an edge length of 0.2 units, being aligned with the three axes of the coordinate system. In cases where the cube exceeded the bounds of the domain, we only considered its volume (cuboid) being located inside the domain. The faces of this cube or cuboid redefined the bounds of the domain over which we built an equidistant grid with a ten-times-higher resolution of 0.01 units. Again, the initial and terminal grid points of each dimension were located at the bounds of the domain. We repeated this procedure until we determined the minimum of $u_{\text{LCB}}(\mathbf{x})$ with a resolution of 0.0001 units. Our algorithm stopped after determining the same parameter set twice; here, after 44 iterations. The optimal parameter set was used to calculate both D3(BJ) energies and eigD3(BJ)/histD3(BJ) feature vectors.

3.4 Performance Measures

All performance measures considered in this work are defined with respect to the two vectors $\mathbf{y}$ (observations) and $\hat{\mathbf{y}}$ (predictions). The mean error (ME) represents the centroid of the error vector $\boldsymbol{\Delta}=\mathbf{y}-\hat{\mathbf{y}}$,

[TABLE]

The root-mean-square deviation (RMSD) is a measure of the scatter of $\mathbf{\Delta}$,

[TABLE]

The root-mean-square error (RMSE) can be expressed as a function of both ME and RMSD,

[TABLE]

The mean absolute error (MAE),

[TABLE]

is the $L_{1}$-norm analog to the RMSE. All performance measures presented can be transformed into dimensionless quantities by dividing each element of the error vector by the absolute value of the corresponding element in the observation vector,

[TABLE]

All $\{\Delta_{n}\}$ in Eqs. (17)–(20) need to be replaced by $\{\delta_{n}\}$ to obtain relative-performance measures. Examples are the mean relative error,

[TABLE]

and the mean absolute relative error,

[TABLE]

Note that relative-performance measures yield unreasonably large values in cases where the absolute value of an observation is smaller than or comparable to the underlying measurement/numerical uncertainty, i.e., when it is close to zero. In such cases, a capped version of the relative error,

[TABLE]

avoids running into singularity issues. The cap $c$ is a strictly positive real number.

4 Results and Discussion

4.1 Assessment of Interaction Energies

Smith et al.[28] included the S66x8 data set developed by Řezáč et al.[57] among other data sets for their reparametrization of PBE-D3(BJ) (other density functional were also considered). They also adopted the corresponding CP-corrected CCSD(T)/CBS∗ interaction energies. We introduce the asterik to distinguish the extrapolation scheme employed by Řezáč et al. from the one considered in this work. While Smith et al. employed the same basis set (def2-QZVP) as Grimme and co-workers[25] for calculating PBE interaction energies, they applied CP corrections instead of neglecting basis set superposition errors. Considering the S13x8 subset, we compared the results presented in this work against those published by Smith et al. In the following, we report all relative and signed errors with respect to their study.

We find a squared correlation coefficient of $r^{2}=0.9999$ between PBE/def2-QZVP (Smith et al.[28]) and PBE/ma-def2-QZVPP (this work) interaction energies. The MAE and the corresponding RMSD both amount to 0.0087 kcal mol*-1*. Given that the CP-corrected CCSD(T)/CBS∗ interaction energies are reported with an accuracy of 0.01 kcal mol*-1*, we conclude that the effect of basis set alteration on PBE interaction energies is insignificant in this case.

Furthermore, we find a squared correlation coefficient of $r^{2}=0.9954$ between CCSD(T)/CBS∗ (Řezáč et al.[57]) and DLPNO-CCSD(T)/CBS (this work) interaction energies. MAE and RMSD amount to 0.11 kcal mol*-1* and 0.12 kcal mol*-1*, respectively. We consider these deviations significant given the limiting resolution of 0.01 kcal mol*-1* mentioned above. For the five smallest relative distances (0.90, 0.95, 1.00, 1.05, and 1.10 with respect to the equilibrium distance, $r_{\text{eq}}$), we find that the unsigned ME ($|$ME$|$) equals the MAE. As the ME is negative in all of these cases, the corresponding DLPNO-CCSD(T)/CBS dispersion interactions are, without exception, determined to be less attractive than the ones obtained from CCSD(T)/CBS∗ calculations. For the three largest relative distances (1.25$r_{\text{eq}}$, 1.50$r_{\text{eq}}$, and 2.00$r_{\text{eq}}$), the unsigned ME represents at least 70 % of the MAE. In these cases, the ME is positive, suggesting that DLPNO-CCSD(T)/CBS tends to predict more attractive dispersion corrections at long distances compared to CCSD(T)/CBS∗. However, the ME over the entire S13x8 set amounts to $-$0.09 kcal mol*-1*.

The CBS extrapolation scheme employed by Řezáč et al.[57] involves MP2 energies obtained with aug-cc-pVTZ and aug-cc-pVQZ basis sets, and approximates the CBS limit of the correlation energy by adding to the MP2/CBS energy a CCSD(T)/aug-cc-pVDZ correction. In this work, however, we extrapolate directly from DLPNO-CCSD(T)/aug-cc-pV$\nicefrac{{T}}{{Q}}$Z energies as defined in Eq. (3.2.1). Therefore, it is not only the underlying electronic-structure approximation but also the extrapolation scheme that is different from the methodology reported by Řezáč et al. While CCSD(T) should by construction capture more long-range interactions than DLPNO-CCSD(T), our results are based on triple/quadruple-$\zeta$ basis sets (aug-cc-pV$\nicefrac{{T}}{{Q}}$Z) compared to a double-$\zeta$ basis set (aug-cc-pVDZ) employed by Řezáč et al. It is debatable whether diffuse basis functions can compensate for the inaccuracies introduced by a double-$\zeta$ basis set. The positive ME determined for the three largest intermolecular distances indicates that DLPNO-CCSD(T)/aug-cc-pV$\nicefrac{{T}}{{Q}}$Z calculations may capture more long-range interactions than CCSD(T)/aug-cc-pVDZ corrections to MP2/aug-cc-pV$\nicefrac{{T}}{{Q}}$Z energies.

To gain more insight, we examined the effect of the CBS extrapolation scheme introduced in Eq. (3.2.1) on predicting interaction energies (as measured by the MAE over S13x8) by varying its two parameters. We find that the interaction energies are insensitive with respect to smaller (more negative) values of the Hartree–Fock parameter (default: $-$1.63, dimensionless, cf. Eq. (3.2.1)). For larger (less negative or positive) values, the MAE increases. Increasing the correlation parameter from its default value of $-$3.00 (dimensionless, cf. Eq. (3.2.1)) to $-$1.17 minimizes the MAE to 0.03 kcal mol*-1* (RMSD $=$ 0.04 kcal mol*-1*), which constitutes a significant reduction of the deviation between DLPNO-CCSD(T)/CBS and CCSD(T)/CBS∗ interaction energies. Moreover, the overall ME is now slightly positive (0.01 kcal mol*-1*). For comparative purposes, we note that ME, MAE, and RMSD of DLNPO-CCSD(T)/aug-cc-pVQZ interaction energies (with respect to CCSD(T)/CBS∗) are larger in magnitude than those of their DLNPO-CCSD(T)/CBS analogs. These quantities amount to $-$0.13 kcal mol*-1*, 0.15 kcal mol*-1*, and 0.16 kcal mol*-1*, respectively. We do not adopt the tuned correlation parameter ($-$1.17) for this study as a reduced deviation from the CCSD(T)/CBS∗ results over S13x8 does not necessarily imply an improvement in accuracy. Therefore, we consider a reparametrization of the D3(BJ) model for PBE as it has previously been trained on CCSD(T)/CBS∗ interaction energies.[25, 28, 30]

4.2 System-focused D3(BJ) Reparametrization

A reparametrization of the D3(BJ) model for PBE is not only important for an unbiased assessment of its performance with respect to DLPNO-CCSD(T)/CBS reference interaction energies. As the eigD3(BJ) and histD3(BJ) feature vectors (cf. Section 3.3.1) are functions of the empirical D3(BJ) parameters $a_{1}$, $a_{2}$, and $s_{8}$, mean and variance obtained from GP regression will also depend on the specific D3(BJ) parameters. A system-focused reparametrization would, therefore, ensure an unbiased performance test of the D3(BJ)-GP model as we expect to obtain close-to-optimal PBE-D3(BJ) interaction energies and, hence, maximally reduced residual interaction energies, Eq. (7), as input for GP regression. Following the line of argument by Smith et al.,[28] we decided to optimize the D3(BJ) parameters for PBE with respect to a relative rather than an absolute error, as the latter tends to bias the parameters toward systems with large absolute interaction energies. In this work, for instance, the largest unsigned DLPNO-CCSD(T)/CBS interaction energy is 320 times larger than the smallest one.

The objective function (MARE) of the optimization procedure is shown in Figure 1. At the global minimum, $a_{1}$ and $s_{8}$ are zero, and $a_{2}=5.6841$ bohr. This remarkable reduction of the number of parameters may be explained with the underlying training set (S13x8-T), which solely consists of dispersion-dominated dimers. The training sets employed in previous D3(BJ) parametrizations[25, 28, 30] were more heterogeneous and also included, e.g., hydrogen-bonded aggregates. The smaller the number of effects that need to be captured, the more likely it is that a simple model captures these effects as accurately as a more complex model. To support this hypothesis, we repeated the optimization with the identical settings and training systems, but this time, we incorporated the CCSD(T)/CBS∗ interaction energies published by Řezáč et al.[57] Again, $a_{1}$ and $s_{8}$ are zero at the global minimum, and $a_{2}$ changes to a slightly lower value of 5.5530 bohr. We do not consider this parameter set in the following. These findings suggest that a simplified expression of the D3(BJ) dispersion energy, defined in Eq. (5), can be applied when studying dispersion-dominated systems (at least in the case of PBE),

[TABLE]

This expression does no longer involve dipole–quadrupole interactions ($n=8$), neither through the dipole–quadrupole term itself nor through the term $\sqrt{C_{8}^{IJ}/C_{6}^{IJ}}$ of the damping function, which is multiplied by $a_{1}$. For $R_{IJ}\rightarrow 0$, the corresponding D3(BJ) term converges to the finite negative value of the simple expression $C_{6}^{IJ}/a_{2}^{6}$.

As the D3(BJ) parameter set for PBE is a function of the training set (here, S13x8-T), the reference interaction energies (here, $E_{\text{int}}^{\text{DLPNO-CCSD(T)/CBS}}$), and the objective (here, MARE), it is not surprising that optimal parameter sets obtained in other studies may differ significantly from the one determined in this study. For instance, Grimme et al.[25] and our group[30] chose a weighted MAE as objective, CCSD(T)/CBS∗ energies, and a training set comprising equilibrium structures only. Whereas Smith et al.[28] also incorporated CCSD(T)/CBS∗ energies, they employed a similar objective (a capped version of the MARE, cf. Eq. (24)) and a training set that also included nonequilibrium structures. These methodological differences may explain why the PBE-specific D3(BJ) parameters by Smith et al. are much more similar to the ones obtained in this work compared to the parameters published by Grimme et al.[25] and previously by our group[30] (note that these parameter sets are, in turn, similar to each other). The similarity to the parameters by Smith et al. underpins once more that the change from CCSD(T)/CBS∗ to DLPNO-CCSD(T)/CBS has only a minor effect on the D3(BJ) model.

The accuracy of the D3(BJ) parameter sets mentioned above with respect to three different validation sets (S13x8-V, ROTA, and CONF) are reported in Table 2. The validation sets considered have no data points in common with the training set (S13x8-T). Our new D3(BJ) parameters for PBE result in the lowest MAREs and MAEs over all validation sets (in the case of S13x8-V, the parameters by Smith et al. yield the same MAE). The new parameter set also leads to the smallest unsigned MRE ($|$MRE$|$) independent of the validation set and, hence, reveals the strongest reduction of systematic errors in D3(BJ)-type dispersion corrections. Based on this problem-specific validation, we conclude that our new PBE-specific D3(BJ) parameter set maximizes the statistical validity of D3(BJ)-type dispersion corrections for the given domain of application. If not otherwise mentioned, the term “D3(BJ)” and all other terms including it refer to that parameter set, i.e., $a_{1}=s_{8}=0$ and $a_{2}=5.6841$ bohr.

4.3 Learning Curves of the D3(BJ)-GP Model

To examine GP regression for interpolation and extrapolation tasks, we present three prediction scenarios with the same underlying training set on three different test sets (Figure 2). Recall that the target variable is the residual interaction energy, $\Delta E_{\text{int}}$, defined in Eq. (7). For the training set, 1,000 dimers were randomly drawn (10 repetitions) from ROTA, which we refer to as ROTA-T. Its complementary holdout set consisting of the 100 remaining ROTA dimers, denoted ROTA-H, is used as the first test set. The other two test sets are CONF and S13x8. Starting from a 1 % fraction of ROTA-T data, we successively added data points in 20 equidistant intervals until the entire ROTA-T set was used for training. After each addition, we performed a hyperparameter optimization. We determined the MARE of both PBE-D3(BJ) interaction energies (independent of the training set and, hence, constant) and the corresponding PBE-D3(BJ)-GP interaction energies. For the latter, we employed the histD3(BJ) featurization. Furthermore, we considered two different kernels: The Matérn-$\nicefrac{{1}}{{2}}$ (Laplacian/exponential) kernel (Figure 2a) and the Matérn-$\nicefrac{{3}}{{2}}$ kernel (Figure 2b) introduced in Eqs. (11) and (13), respectively. During preliminary testing, we found that these kernels perform best across the Matérn kernel family (to which also the Gaussian/squared exponential kernel belongs[54]) for predicting residual interaction energies.

We find that the D3(BJ)-GP model for PBE (solid lines) almost immediately — after training on less than five data points — predicts improved dispersion corrections compared to D3(BJ) (dashed lines) for both the ROTA-H and CONF test sets. The predictions improve at a high rate until about 30 data points of the training set (ROTA-T) are considered. Afterwards, the predictions continue to improve only marginally. The small associated standard deviations (gray bands) indicate that the learning curves are rather insensitive toward the actual composition of the training set. While the Matérn-$\nicefrac{{1}}{{2}}$ kernel leads to slightly smaller prediction errors (3 %) and standard deviations ($\leq$1 %) than the Matérn-$\nicefrac{{3}}{{2}}$ kernel (5–6 % and $\leq$2 %, respectively), both kernels yield qualitatively identical results. Note that, since ROTA-H is complementary to ROTA-T, the PBE-D3(BJ) error for ROTA-H is also associated with a nonzero standard deviation. While each dimer included in the training set (ROTA-T) contains the same pentane conformer, the GP model can, without any delay, generalize to 44 other pentane conformers contained in the CONF dimers. This finding suggests that the CONF set has a large overlap with the ROTA set in feature space (cf. Section 4.4).

Regarding the S13x8 test set, the D3(BJ)-GP model based on the Matérn-$\nicefrac{{1}}{{2}}$ kernel results in consistently larger prediction errors than the D3(BJ) model. By contrast, the D3(BJ)-GP model based on the Matérn-$\nicefrac{{3}}{{2}}$ kernel quickly improves over the D3(BJ) model, even though only slightly. Beyond a 5 % fraction of training data, adding more data to ROTA-T does not improve predictions in both cases. Given that the S13x8 set contains chemically different systems compared to ROTA and CONF (who are solely composed of ethyne–pentane dimers), we consider this test an extrapolation task. To understand this result, we analyzed the hyperparameter $\alpha_{2}$ introduced in Eq. (12) for both kernels, which describes on what length scale in feature space information decays. The smaller the length scale, the higher is the decay rate of the kernel $k(\mathbf{x}_{i},\mathbf{x}_{j})$ when the distance between $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ increases. As we specified a zero-mean prior, the mean of the posterior distribution converges to zero for large distances between $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$. In such cases, D3(BJ)-GP- and D3(BJ)-type dispersion corrections become identical. This property of the posterior distribution is highly appealing as it preserves the accuracy of PBE-D3(BJ) interaction energies when the distance (in feature space) to the training set increases.

For the Matérn-$\nicefrac{{3}}{{2}}$ kernel, $\alpha_{2}$ is on the order of 0.1 kcal mol*-1*. Given that mean (and standard deviation) of the ROTA and S13x8 reference interaction energies are 0.28 kcal mol*-1* (0.25 kcal mol*-1*) and 2.03 kcal mol*-1* (1.19 kcal mol*-1*), respectively, it seems that most of the S13x8 systems have a distance to the ROTA set that is significantly larger than $\alpha_{2}$. To test this hypothesis, we examined the kernel matrix $\mathbf{K}\left(\mathbf{X}^{\textsf{ROTA}},\mathbf{X}^{\textsf{S13x8}}\right)$, where $\mathbf{X}^{\textsf{ROTA}}$ and $\mathbf{X}^{\textsf{S13x8}}$ represent the feature vectors of training and test sets, respectively. The kernel matrix was normalized such that $k(\mathbf{x}_{i},\mathbf{x}_{i})=1$ (arbitrary units), which constitutes the maximum possible value of $k(\mathbf{x}_{i},\mathbf{x}_{j})$. Of the 104 dimers contained in the S13x8 set, only 24 exhibit a kernel for which $k\left(\mathbf{x}_{i}^{\textsf{ROTA}},\mathbf{x}_{j}^{\textsf{S13x8}}\right)>0.1$, which indicates that the ROTA and S13x8 sets have only a small overlap in feature space (confirming the extrapolation assumption stated above). Moreover, for 50 % of the S13x8 dimers, we find that $k\left(\mathbf{x}_{i}^{\textsf{ROTA}},\mathbf{x}_{j}^{\textsf{S13x8}}\right)<0.0005$. Hence, the slight improvement of PBE-D3(BJ)-GP over PBE-D3(BJ) (Figure 2b, right panel) originates from the few S13x8 dimers that are close to ROTA in the histD3(BJ) feature space. For the remaining majority of the S13x8 test set, the D3(BJ)-GP model yields essentially identical dispersion corrections as its D3(BJ) analog.

By contrast, $\alpha_{2}$ is on the order of 104 kcal mol*-1* in the Matérn-$\nicefrac{{1}}{{2}}$ case. This length scale is by several orders of magnitude larger than the range of reference interaction energies considered in this work. The largest element of the normalized kernel matrix $\mathbf{K}\left(\mathbf{X}^{\textsf{ROTA}},\mathbf{X}^{\textsf{S13x8}}\right)$ amounts to 0.9997. Hence, the D3(BJ)-GP model can be considered an almost global model just as its D3(BJ) analog. That the former performs worse than the latter (Figure 2a, right panel) can be attributed to this extensive length scale and the remaining kernel-specific hyperparameter $\alpha_{1}$. They were learned from samples of the ROTA set, which spans a much smaller energy range than the S13x8 set. While these hyperparameters perform well for interpolation tasks (cf. Figure 2a, left and central panels), they perform poorly for extrapolation tasks. Note that the deviation of the S13x8-related learning curve to the PBE-D3(BJ) error (dashed line) is significantly larger than the standard deviation of the learning curve (gray band). This behavior is an indication of underfitting; the prediction model is not yet flexible enough to generalize to data outside the training domain.

The results discussed above and shown in Figure 2 were derived from the histD3(BJ) featurization. The complementary eigD3(BJ) results are reported in the Supporting Information (Figure S1). They are in qualitative agreement with the histD3(BJ) results, but mean and standard deviation of the eigD3(BJ) learning curves are systematically larger than those of the histD3(BJ) learning curves. This relatively poor performance can be explained in part with the properties of the D3(BJ) matrix, introduced in Eq. (8), of the interacting system. As its diagonal-block elements are effectively zero, each nonzero element in the eigD3(BJ) vector with an odd index is nearly identical in magnitude to the next (evenly indexed) element. They only differ by their signs, i.e., there is always a pair of one positive eigenvalue and one negative eigenvalue which share almost identical absolute values. The order of signs seems to be randomly flipped, which introduces some degree of noise to the eigD3(BJ) vector. Applying a strict order of signs led to systematically smaller learning curves compared to a random order, but they are still systematically larger than the histD3(BJ) learning curves.

We also analyzed the learning curves of a linear model with respect to the same training and test sets (Figure S2). Both histD3(BJ) and eigD3(BJ) feature vectors served (separately) as input variables. In all cases, the learning curves obtained with the linear models exhibit significantly larger mean values and standard deviations than those obtained with our GP models. We also explored other polynomial models (constant, quadratic, cubic, and quartic); the results are of the same quality. Please see the Supporting Information for further details. We conclude that linear regression based on polynomial models is not well suited for the prediction tasks considered in this study.

4.4 Data Diversity Analysis

To interpret the results of the previous section from a data perspective, we put the different data sets (ROTA, CONF, and S13x8) into a mutual context by analyzing the distributions of their molecular systems over different variables. As the intermolecular distance and the relative spatial orientation both have an essential effect on the magnitude of noncovalent interactions, an intuitive way to structure the data is to unravel these two properties (Figure 3a). For this purpose, we introduce two measures of the intermolecular distance: the minimum interatomic distance, $R_{\text{min}}$, and the cumulative reciprocal interatomic distance,

[TABLE]

While $R_{\text{min}}$ correlates with the largest interatomic noncovalent interaction between two molecules, $\rho_{\text{inv}}$ encodes the overall noncovalent interactions between two molecules. This data decomposition confirms that (i) CONF has a large overlap with ROTA and (ii) S13x8 is well-separated from ROTA. The uniform sampling of the ROTA dimers with respect to both distance and relative orientation is reflected by its uniformly distributed data cloud. The relatively large $\rho_{\text{inv}}$ values of the S13x8 set can be mainly attributed to larger monomer sizes compared to the ROTA and CONF set. It is not surprising that the ethyne–pentane and (subsequently) ethene–pentane series of the S13x8 set are closest to the ROTA–CONF data cloud in this representation. Note that the ethyne–pentane series represents an upper bound of $\rho_{\text{inv}}$ to the ROTA–CONF cloud, which originates from the parallel orientation of its dimers.

The majority of systems over all data sets lies in the high-energy regime as measured by the intermolecular D3(BJ) dispersion energy. Hence, the majority of systems is located in a dense cluster in the eigD3(BJ) and histD3(BJ) feature spaces, respectively. The most informative linear representation of a feature space in two dimensions can be obtained by selecting its first two principal components. Here, we applied a principal component analysis[40] to the histD3(BJ) feature space (Figure 3b). Its first principle component reveals an almost perfect correlation with the intermolecular D3(BJ) dispersion energy (${r^{2}=0.9917}$). To minimize the number of training data (and hence, the number of electronic-structure calculations) needed to achieve a predefined accuracy over the entire feature domain considered, we would choose a sampling strategy that samples this domain uniformly. The black crosses in Figure 3 represent such a training set, which could be obtained from, e.g., $k$-means clustering.[40] By contrast, random sampling would mainly select data located in the ROTA–CONF cluster due to its large number of instances, which would lead to poor generalizations into the low-energy feature subspace. However, sampling the feature space uniformly is not sufficient as we would not know how many data points to pick in order to achieve a predefined accuracy. For this purpose, we harness the variance of the GP posterior distribution.

4.5 Variance-Based Selection of Training Data

We considered an initial training set $\mathcal{T}_{8}^{8}$ comprising eight randomly drawn target values and the corresponding feature vectors. This set was randomly drawn from the overall data set of 1,248 dimers (S13x8 + ROTA + CONF). The set of all 1,240 remaining dimers, the query set $\mathcal{Q}_{1240}$, incorporates only the associated feature vectors. In a practical setup, electronic-structure calculations would not yet be available for the query set. After GP regression on $\mathcal{T}_{8}^{8}$ including hyperparameter optimization, the GP posterior distribution was determined for all feature vectors in $\mathcal{Q}_{1240}$ (query vectors). We chose three different batch sizes (${L=}$ 1, 40, 1,240), to demonstrate the advantages and limitations of batch learning (cf. Section 2.2) compared to sequential learning[55] (${L=1}$). In Figure 4, we compare BVS against random sampling by examining the maximum posterior standard deviation, $\hat{\sigma}_{y,\text{max}}$, of the query set as a function of the number of training data. The maximum standard deviation can be interpreted as the maximum prediction uncertainty. To cover the entire range of training data, we set the accuracy threshold $t$ to zero. The results shown correspond to the Matérn-$\nicefrac{{1}}{{2}}$ kernel and the histD3(BJ) featurization.

Initially, $\hat{\sigma}_{y,\text{max}}$ fluctuates strongly in the case of variance-based sampling with sequential learning (${L=1}$, blue line). As the effect of a new data point on the GP posterior distribution is larger for small training set sizes, the hyperparameters are more sensitive to new data points in the initial phase. This hyperparameter sensitivity explains the observed fluctuations. After adding about 50 training data points, $\hat{\sigma}_{y,\text{max}}$ decays rapidly, indicating a rapid improvement of the prediction accuracy with respect to the query set. The profile of $\hat{\sigma}_{y,\text{max}}$ for BVS with a batch size of ${L=40}$ (black solid line) is very similar to the ${L=1}$ case, but does not exhibit fluctuations in the initial phase. This apparently more stable behavior can be explained by the initially fixed hyperparameters, which are only reoptimized after the addition of 40 query vectors to the training set. As $t=0$, all following hyperparameter optimizations (indicated by the black circles) are also performed after selecting batches of 40 query vectors. In a practical setup, this approach would be 40 times more efficient than the ${L=1}$ approach, without significantly altering the number of training data points for a given threshold beyond 50 training data points. For the hypothetical thresholds ${t_{1}=0.05}$ kcal mol*-1* and ${t_{2}=0.01}$ kcal mol*-1* (gray dashed lines), the algorithm would stop at about 200 and 770 training data points, respectively, in both cases (${L=1}$ and ${L=40}$).

The profile of $\hat{\sigma}_{y,\text{max}}$ will change significantly if the hyperparameters are not updated while adding training data points as in the case of BVS with a batch size of ${L=}$ 1,240 (red line). Here, the maximum posterior standard deviation is significantly overestimated throughout. For instance, the first hypothetical threshold ${t_{1}=0.05}$ kcal mol*-1* is only reached after adding 821 data points to the training set, which would amount to an overhead of about 3,700 DLPNO-CCSD(T) single-point calculations compared to the ${L=1}$ and ${L=40}$ cases. Hence, the actual batch size employed in BVS is an important parameter that requires careful adjustment.

Compared to BVS, the decay rate of $\hat{\sigma}_{y,\text{max}}$ with respect to the number of training data is much lower in the case of random sampling (black dashed line and gray band resulting from 5 repetitions; $L=40$). For instace, the first hypothetical threshold ${t_{1}=0.05}$ kcal mol*-1* is only reached after adding almost all query vectors (1,238) to the training set. As the posterior variance is not taken into account, potentially many uninformative (low-variance) query vectors are drawn, which renders random sampling highly inefficient for the generation of small informative training sets.

It seems that BVS with a batch size of ${L=40}$ is the most efficient option among the ones considered. However, the posterior variance is only an estimate of the scatter of the prediction error and not the prediction error itself. As the GP posterior distribution is Gaussian, 95.5% of all D3(BJ)-GP errors should be smaller than $\hat{\sigma}_{y,\text{max}}$. We find that the Matérn-$\nicefrac{{1}}{{2}}$ kernel fulfills this condition after the addition of the first and all following batches. While the Matérn-$\nicefrac{{3}}{{2}}$ yields slightly higher decay rates of $\hat{\sigma}_{y,\text{max}}$ (Figure S3), it violates this condition in an unpredictable manner. As the statistical validity of the posterior variance is of utmost importance for the efficient and reliable generation of informative training sets, we choose to incorporate the Matérn-$\nicefrac{{1}}{{2}}$ kernel into the D3(BJ)-GP model.

The actual prediction error of the D3(BJ)-GP model for PBE obtained from both BVS and random sampling is shown in Figure 5. It is based on the same initial training and query sets as those employed for the production of Figure 4. In all cases, a batch size of ${L=40}$ was chosen. We employed the Matérn-$\nicefrac{{1}}{{2}}$ kernel and both featurizations introduced in this work, eigD3(BJ) (black lines) and histD3(BJ) (red lines). Mean (dashed lines) and standard deviation (gray bands) of the prediction errors obtained via random sampling are the result of 5 repeated draws. Additionally, we report the prediction error of the D3(BJ) model for PBE (green line). Note that we consider the MAE instead of the MARE as the posterior variance is learned from absolute residual interaction energies. Learning their relative analog is inefficient as the division by the reference interaction energy would annihilate the smoothness of the target function gained by applying the $\Delta$-learning approach[82] introduced by von Lilienfeld and co-workers.

BVS leads to significantly smaller prediction errors than random sampling, independent of the number of training data considered. This finding confirms the hypothesis that query data associated with a high posterior variance likely contribute to the overall error in a significant way. For both sampling strategies, we find that the histD3(BJ) featurization leads to systematically smaller prediction errors than its eigD3(BJ) analog, which is in agreement with the findings discussed in Section 4.3. Note that the query sets obtained with random sampling and BVS are different and, therefore, a comparison of the corresponding prediction errors is somewhat biased. For this reason, we only refer to the prediction errors associated with BVS in the following.

As the query set decreases constantly, two effects are intermingled in the D3(BJ)-GP prediction error: The degree of improvement by adding more training data (learning rate) and the improvement by removing instances with large absolute errors from the query set. The change of the D3(BJ) prediction error supports unfolding these effects as it is not affected by the training data (it would be constant for a fixed test set). Hence, the change in the difference between the D3(BJ)-GP and D3(BJ) prediction errors can be interpreted as the learning rate of the D3(BJ)-GP model. It only increases until the addition of about 200 training data points (and again for $>$1,000 training data points). However, the learning rate is positive throughout as the D3(BJ)-GP prediction error is constantly smaller than its D3(BJ) complement. This finding is independent of the featurization.

5 Conclusions

We introduced the D3(BJ)-GP model, which adjusts for systematic errors in D3(BJ)-type dispersion corrections by harnessing the capabilities of Gaussian process (GP) regression. The model learns mappings from linear combinations of atom-pairwise D3(BJ) terms (input) to errors in PBE-D3(BJ)/ma-def2-QZVPP interaction energies (target) with respect to DLPNO-CCSD(T)/CBS. We generated input–target pairs for 1,248 molecular dimers, which resemble the dispersion-dominated systems contained in the S66x8 data set.[58, 57] Our systems represent various out-of-equilibrium distances, conformations, and relative orientations. As the empirical D3(BJ) parameters ($a_{1}$, $a_{2}$, $s_{8}$) were previously optimized with respect to CCSD(T)/CBS interaction energies, we performed a reparametrization of the D3(BJ) model for PBE with respect to 66 dimers selected from our data set. At the global minimum, $a_{1}$ and $s_{8}$ equal zero and $a_{2}=5.6841$ bohr, effectively eliminating all dipole–quadrupole terms from the D3(BJ) energy expression. We found that these three parameters perform best for the remaining 1,182 dispersion-dominated dimers compared to previously published PBE-specific D3(BJ) parameter sets.[25, 28, 30] When repeating the optimization after replacing the DLPNO-CCSD(T)/CBS interaction energies by their CCSD(T)/CBS analogs, the global minimum was still located at $a_{1}$ and $s_{8}$ equal to zero. This finding suggests that a single empirical D3(BJ) parameter, $a_{2}$, is sufficient to yield reliable dispersion corrections for dispersion-dominated systems, at least for the PBE functional.

Compared to the D3(BJ) model, which has a fixed functional form and is global with respect to its parameters, the D3(BJ)-GP model has a high degree of functional flexibility as it is based on a local, nonlinear kernel. Therefore, D3(BJ)-GP-type dispersion corrections can be significantly improved for a given domain of application by adding new input–target pairs to the training set representing that domain. Regarding interpolation, we showed that the D3(BJ)-GP model is associated with a steep learning curve and starts to improve over D3(BJ) after training it on less than five data points. After the addition of about 30 data points, the predictions continued to improve only marginally. Regarding extrapolation, the accuracy of dispersion corrections improved only slightly (by about 1 %) or even declined with a close-to-zero or even negative learning rate. These results show that the performance improvement of the D3(BJ)-GP model critically depends on the composition of the training set. It is not enough to sample the domain of application uniformly as we do not know the number of data points required to achieve a certain accuracy beforehand. Furthermore, uniform sampling does not distinguish between relevant and irrelevant input dimensions, which are also unknown beforehand. Hence, selecting the training set a priori, as done in the case of the original D3(BJ) model, is a poor strategy for learning D3(BJ)-GP-type dispersion corrections.

To optimally increase the performance of the D3(BJ)-GP model, it is built on dynamic training sets. We aimed at selecting training sets that are maximally informative, i.e., that contain the least number of data points to achieve a predefined accuracy over a given domain of application. For this purpose, we exploited the variance of the GP posterior distribution. Of all potential training systems considered, the system with the highest posterior variance is selected (referred to as variance-based sampling). If interaction energies are not yet available for the new training system, electronic-structure calculations must be carried out. Subsequently, the GP posterior distribution is updated and the next training system is selected. Through this active learning process, we obtain a system-focused, self-improving model for dispersion interactions equipped with confidence intervals. We demonstrated that the posterior variance can be approximately updated from only the input variables of the new training system, which are obtained efficiently from D3(BJ) calculations. This way, we collect a batch of new training systems before the corresponding electronic-structure calculations are being carried out at the same time. We refer to this selection approach as batch-wise variance-based sampling (BVS). BVS-guided active learning renders our D3(BJ)-GP workflow efficient and controllable. We showed that dynamic training sets significantly improve the learning rates of the D3(BJ)-GP model. For the data set studied, it systematically yields improved dispersion corrections, which can be calculated almost as efficiently as their semiclassical D3(BJ) complements. Like its predecessor D3, the D3-GP approach is implemented in a black-box fashion, the difference being that the latter requires a specification of the training set (representing the domain of application). Once benchmark data are provided for the associated systems (we recommend an automated selection by our BVS-guided active learning scheme), the D3(BJ)-GP model will update itself through an optimization of the underlying GP hyperparameters. Note that the D3 parameters (here, of the BJ damping function) are not affected by this procedure. They are only determined once as in the original D3 approach.

Although we have examined the D3-GP approach with respect to the PBE functional and the BJ damping function only, the D3-GP workflow can be straightforwardly generalized without further effort as it is independent of the electronic-structure approximation and the damping scheme. As it is also independent of the dispersion coefficient model, the D3-GP workflow would essentially be the same for the D4 approach. Moreover, modifications of the input variables (molecular representations) will render any D$x$-GP approach also applicable for other, non-D$x$-type dispersion corrections. One could also learn reference interaction energies directly — instead of indirectly via the residual interaction energy defined in Eq. (7). In that case, however, the physical information already encoded in the quantity to be corrected (here, PBE-D3(BJ)/ma-def2-QZVPP interaction energies) would need to be encoded into the input variables, which is an intricate and ambiguous task. Eventually, the workflow is not even problem-dependent. Our BVS-guided active learning scheme combined with GP regression could be helpful for a diverse range of prediction tasks.

To make our D$x$-GP approach generally accessible, we will release a version-controlled database in future work. This database can then be used and updated to collectively improve on the D$x$-GP model. Version control is crucial as the database is dynamic, which affects the composition of the training set generated by BVS-guided active learning. With version control, previous results can be permanently reproduced even though the database has been altered meanwhile.

Acknowledgments

J.P. acknowledges funding through an Early Postdoc.Mobility fellowship by the Swiss National Science Foundation (project no. 178463). S.G. and M.R. are grateful for financial support by the Swiss National Science Foundation (project no. 200021_182400).