A Universal Density Matrix Functional from Molecular Orbital-Based Machine Learning: Transferability across Organic Molecules

TL;DR

This paper demonstrates that molecular-orbital-based machine learning (MOB-ML) can accurately predict correlation energies across diverse organic molecules with significantly fewer training data, showing high transferability and efficiency.

Contribution

The study introduces refined feature strategies for MOB-ML and proves its transferability and data efficiency in predicting correlation energies for a wide range of organic molecules.

Findings

MOB-ML achieves chemical accuracy with three times fewer training geometries than Δ-ML.

MOB-ML predicts energies for larger molecules with 36 times fewer training calculations.

High transferability of MOB-ML models across different molecular sizes is demonstrated.

Abstract

We address the degree to which machine learning can be used to accurately and transferably predict post-Hartree-Fock correlation energies. Refined strategies for feature design and selection are presented, and the molecular-orbital-based machine learning (MOB-ML) method is applied to several test systems. Strikingly, for the MP2, CCSD, and CCSD(T) levels of theory, it is shown that the thermally accessible (350 K) potential energy surface for a single water molecule can be described to within 1 millihartree using a model that is trained from only a single reference calculation at a randomized geometry. To explore the breadth of chemical diversity that can be described, MOB-ML is also applied to a new dataset of thermalized (350 K) geometries of 7211 organic models with up to seven heavy atoms. In comparison with the previously reported $\Delta$-ML method, MOB-ML is shown to reach chemical accuracy with three-fold fewer training geometries. Finally, a transferability test in which models trained for seven-heavy-atom systems are used to predict energies for thirteen-heavy-atom systems reveals that MOB-ML reaches chemical accuracy with 36-fold fewer training calculations than $\Delta$-ML (140 versus 5000 training calculations).

Tables

Table 1. Table 1: Number of features selected as a function of the number randomly chosen training molecules for the QM7b-T dataset at the CCSD(T)/cc-pVDZ level. The number of features that exceed an importance threshold as well as the inverse participation ratio (IPR) of the feature importance scores are reported (see text).

	# of important features		feature weight IPR
Training size	${ϵ}_d^{\mathrm{ML}}$	${ϵ}_o^{\mathrm{ML}}$	${ϵ}_d^{\mathrm{ML}}$	${ϵ}_o^{\mathrm{ML}}$
20	50	28	4.720	1.116
50	46	28	3.718	1.097
100	46	26	3.450	1.115
200	42	24	3.430	1.120

Equations

$$E_{\text{c}}=\sum _{ij}^{\text{occ}}{ϵ}_{ij},$$

$${ϵ}_{ij}=ϵ\left[\{{\varphi }_p{\}}^{ij}\right].$$

$${ϵ}_{ij}^{\text{MP2}}=\frac{1}{4}\sum _{ab}^{\mathrm{virt}}\frac{{\left|⟨ij||ab⟩\right|}^2}{e_a+e_b-e_i-e_j}$$

$${ϵ}_{ij}\approx {ϵ}^{\mathrm{ML}}\left[{\mathbf{f}}_{ij}\right]$$

$${\mathbf{f}}_{ij}=$$

$$J_{ii},J_{ij},J_{jj},{\mathbf{J}}_i^{\text{o}},{\mathbf{J}}_j^{\text{o}},{\mathbf{J}}_i^{\text{v}},{\mathbf{J}}_j^{\text{v}},{\mathbf{J}}_{ij}^{\text{vv}},$$

$$K_{ij},{\mathbf{K}}_i^{\text{o}},{\mathbf{K}}_j^{\text{o}},{\mathbf{K}}_i^{\text{v}},{\mathbf{K}}_j^{\text{v}},{\mathbf{K}}_{ij}^{\text{vv}}).$$

$${ϵ}_{ij}\approx \{\begin{array}{ll}{ϵ}_{\mathrm{d}}^{\mathrm{ML}}\left[{\mathbf{f}}_i\right]& \text{if i=j}\\ {ϵ}_{\mathrm{o}}^{\mathrm{ML}}\left[{\mathbf{f}}_{ij}\right]& \text{if i=j}\end{array}$$

$$M_{ii}$$

$$M_{jj}$$

$$M_{ij}$$

$$M_{ip}$$

$$M_{jp}$$

$$R_a^{ij}=\Vert ⟨{\varphi }_i|\hat{R}|{\varphi }_i⟩-⟨{\varphi }_a|\hat{R}|{\varphi }_a⟩\Vert +\Vert ⟨{\varphi }_j|\hat{R}|{\varphi }_j⟩-⟨{\varphi }_a|\hat{R}|{\varphi }_a⟩\Vert ,$$

Full text

Machine learning (ML) has recently seen wide application in chemistry, including in the fields of drug discovery,Lavecchia (2015); Gawehn, Hiss, and Schneider (2016); Popova, Isayev, and Tropsha (2018) materials design,Kim et al. (2017); Ren et al. (2018); Butler et al. (2018); Sanchez-Lengeling and Aspuru-Guzik (2018) and reaction prediction.Wei, Duvenaud, and Aspuru-Guzik (2016); Raccuglia et al. (2016); Ulissi et al. (2017); Segler and Waller (2017); Segler, Preuss, and Waller (2018) In the context of quantum chemistry, much work has focused on predicting electronic energies or densities based on atom- or geometry-specific features,Smith, Isayev, and Roitberg (2017); S Smith et al. (2018); Lubbers, Smith, and Barros (2018); Bartók et al. (2010); Rupp et al. (2012); Montavon et al. (2013a); Hansen et al. (2013); Gasparotto and Ceriotti (2014); Ramakrishnan et al. (2015); Behler (2016); Kearnes et al. (2016); Paesani (2016); Schütt et al. (2017); Brockherde et al. (2017); Smith, Isayev, and Roitberg (2017); Wu et al. (2018); Nguyen et al. (2018); Yao et al. (2018); Fujikake et al. (2018); Li et al. (2018); Grisafi et al. (2019); Zhang et al. (2018) although other strategies have also been employed.McGibbon et al. (2017) Recently, we reported an accurate and transferable molecular-orbital-based machine learning (MOB-ML) approach to the prediction of correlated wavefunction energies based on input features from a self-consistent field calculation, such as the Hartree-Fock (HF) method.Welborn, Cheng, and Miller III (2018) In this communication, we present refinements to the MOB-ML method with comparisons to test cases from our previous work. We then demonstrate the performance of MOB-ML across a broad swath of chemical space, as represented by the QM7bMontavon et al. (2013a) and GDB-13Blum and Reymond (2009) test sets of organic molecules.

I Theory

The current work aims to predict post-Hartree-Fock correlated wavefunction energies using features of the Hartree-Fock molecular orbitals (MOs). The starting point for the MOB-ML methodWelborn, Cheng, and Miller III (2018) is that the correlation energy can be decomposed into pairwise occupied MO contributionsNesbet (1958); Szabo and Ostlund (1996)

[TABLE]

where the pair correlation energy $\epsilon_{ij}$ can be written as a functional of the full set of MOs, $\{\phi_{p}\}$, appropriately indexed by $i$ and $j$

[TABLE]

The functional $\epsilon$ is universal across all chemical systems; for a given level of correlated wavefunction theory, there is a corresponding $\epsilon$ that maps the HF MOs to the pair correlation energy, regardless of the molecular composition or geometry. Furthermore, $\epsilon$ simultaneously describes the pair correlation energy for all pairs of occupied MOs (i.e., the functional form of $\epsilon$ does not depend on $i$ and $j$). For example, in second-order Møller-Plessett perturbation theory (MP2)Møller and Plesset (1934), the pair correlation energies are

[TABLE]

where $a$ and $b$ index virtual MOs, $e_{p}$ is the Hartree-Fock orbital energy corresponding to MO $\phi_{p}$, and $\left\langle ij\right|\left|ab\right\rangle$ are antisymmetrized electron repulsion integrals.Szabo and Ostlund (1996) A corresponding expression for the pair correlation energy exists for any post-Hartree-Fock method, but it is typically costly to evaluate in closed form.

In MOB-ML, a machine learning model is constructed for the pair energy functional

[TABLE]

where $\mathbf{f}_{ij}$ denotes a vector of features associated with MOs $i$ and $j$. Eq. 4 thus presents the opportunity for the machine learning of a universal density matrix functional for correlated wavefunction energies, which can be evaluated at the cost of the MO calculation.

Following our previous work,Welborn, Cheng, and Miller III (2018) the features $\mathbf{f}_{ij}$ correspond to unique elements of the Fock ($\mathbf{F}$), Coulomb ($\mathbf{J}$), and exchange ($\mathbf{K}$) matrices between $\phi_{i}$, $\phi_{j}$, and the set of virtual orbitals. In the current work, we additionally include features associated with matrix elements between pairs of occupied orbitals for which one member of the pair differs from $\phi_{i}$ or $\phi_{j}$ (i.e., non-$i$,$j$ occupied MO pairs). The feature vector takes the form

[TABLE]

where for a given matrix ($\mathbf{F}$, $\mathbf{J}$, or $\mathbf{K}$) the superscript o denotes a row of its occupied–occupied block, the superscript v denotes a row of its occupied–virtual block, and the superscript vv denotes its virtual–virtual block. Redundant elements are removed, such that the virtual–virtual block is represented by its upper triangle and the diagonal elements of $\mathbf{K}$ (which are identical to those of $\mathbf{J}$) are omitted. To increase transferability and accuracy, we choose $\phi_{i}$ and $\phi_{j}$ to be localized molecular orbitals (LMOs) rather than canonical MOs and employ valence virtual LMOsKnizia (2013) in place of the set of all virtual MOs (as detailed in Ref. 35). We separate Eq. 4 to independently machine learn the cases of $i=j$ and $i\neq j$,

[TABLE]

where $\mathbf{f}_{i}$ denotes $\mathbf{f}_{ii}$ (Eq. 5) with redundant elements removed; by separating the pair energies in this way, we avoid the situation where a single ML model is required to distinguish between the cases of $i=j$ and $\phi_{i}$ being nearly degenerate to $\phi_{j}$, a distinction which can represent a sharp variation in the function to be learned.

In the current work, several technical refinements are introduced to improve training efficiency (i.e., the accuracy and transferability of the model as a function of the number of training examples). These are now described.

Occupied LMO symmetrization. The feature vector is preprocessed to specify a canonical ordering of the occupied and virtual LMO pairs. This reduces permutation of elements in the feature vector, resulting in greater ML training efficiency. Matrix elements $M_{ij}$ ($\mathbf{M}=\mathbf{F}$, $\mathbf{J}$, $\mathbf{K}$) associated with $\phi_{i}$ and $\phi_{j}$ are rotated into gerade and ungerade combinations

[TABLE]

with the sign convention that $F_{ij}$ is negative. Here, $p$ indexes any LMO other than $i$ or $j$ (i.e. an occupied LMO $k$, such that $i\neq k\neq j$, or a valence virtual LMO).

LMO sorting. The virtual LMO pairs are sorted by increasing distance from occupied orbitals $\phi_{i}$ and $\phi_{j}$. Sorting in this way ensures that features corresponding to valence virtual LMOs are listed in decreasing order of heuristic importance, and that the mapping between valence virtual LMOs and their associated features is roughly preserved. We recognize this issue could also potentially be addressed through the use of symmetry functions,Behler and Parrinello (2007) but these are not employed in the current work.

For purposes of sorting, distance is defined as

[TABLE]

where $\phi_{a}$ is a virtual LMO, $\hat{R}$ is the Cartesian position operator, and $\|.\|$ denotes the 2-norm. $\left\|\left\langle\phi_{i}\right|\hat{R}\left|\phi_{i}\right\rangle-\left\langle\phi_{a}\right|\hat{R}\left|\phi_{a}\right\rangle\right\|$ represents the Euclidean distance between the centroids of orbital $i$ and orbital $a$. Previously,Welborn, Cheng, and Miller III (2018) distances were defined based on Coulomb repulsion, which was found to sometimes lead to inconsistent sorting in systems with strongly polarized bonds. The non-$i$,$j$ occupied LMO pairs are sorted in the same manner as the virtual LMO pairs.

Orbital localization. We employ Boys localizationBoys (1960) to obtain the occupied LMOs, rather than IBO localizationKnizia (2013) which was employed in our previous work.Welborn, Cheng, and Miller III (2018) Particularly for molecules that include triple bonds or multiple lone pairs, it is found that Boys localization provides more consistent localization as a function of small geometry changes than IBO localization; and the chemically unintuitive mixing of $\sigma$ and $\pi$ bonds in Boys localization (“banana bonds”)Kaldor (1967) does not present a problem for the MOB-ML method.

Feature selection. Prior to training, automatic feature selection is performed using random forest regression Breiman (2001a) with the mean decrease of accuracy criterion (sometimes called permutation importance).Breiman (2001b) This technique was found to be more effective than our previous useWelborn, Cheng, and Miller III (2018) of the Gini importance scoreBreiman (2001a) which led to worse accuracy and failed to select any features for the case of methane.

The reason for using feature selection in this way is twofold. First, GPR performance is known to degrade for high-dimensional datasets (in practice 50-100 features);Tripathy, Bilionis, and Gonzalez (2016) and second, the use of the full feature set with small molecules can lead to overfitting as features can become correlated.

A reference Python implementation for generating MOB-ML features is provided online. 111See https://github.com/thomasfmiller/MOB-ML for the available code.

II Computational details

Results are presented for a single water molecule; a series of alkane molecules; a thermalized version of the QM7b set of 7211 molecules with up to seven C, O, N, S, and Cl heavy atoms; and a thermalized version of the GDB-13 set of molecules with thirteen C, O, N, S, and Cl heavy atoms. All datasets employed in this work are provided in the Supporting Information.

Training and test geometries are sampled at 50 fs intervals from ab initio molecular dynamics trajectories performed with the Q-Chem 5.0 software package,Shao et al. (2015) using the B3LYPVosko, Wilk, and Nusair (1980); Lee, Yang, and Parr (1988); Becke (1993); Stephens et al. (1994)/6-31g*Hariharan and Pople (1973) level of theory and a Langevin thermostatBussi and Parrinello (2007) at 350 K.

The features and training pair energies associated with these geometries are computed using the Molpro 2018.0 software packageWerner et al. (2018) in a cc-pVTZ basis set unless otherwise noted.Dunning (1989) Valence virtual orbitals used in feature construction are determined with the Intrinsic Bond Orbital method.Knizia (2013) Reference pair correlation energies are computed with second-order Møller-Plessett perturbation theory (MP2)Møller and Plesset (1934); Saebo and Pulay (1993) and coupled cluster with singles and doubles (CCSD)Čížek (1966); Hampel and Werner (1996) as well as with perturbative triples (CCSD(T)).Bartlett et al. (1990); Schütz (2000) Density fitting for both Coulomb and exchange integrals Polly et al. (2004) is employed for all results below except those corresponding to the water molecule. The frozen core approximation is used in all cases.

Gaussian process regression (GPR)Rasmussen and Williams (2006) is employed to machine learn $\epsilon^{\mathrm{ML}}_{\mathrm{d}}$ and $\epsilon^{\mathrm{ML}}_{\mathrm{o}}$ (Eq. 6) using the GPy 1.9.6 software package. GPy (2012) The GPR kernel is Matérn 5/2 with white noise regularizationRasmussen and Williams (2006). Kernel hyperparameters are optimized with respect to the log marginal likelihood objective for the water and alkane series results, as well as for $\epsilon^{\mathrm{ML}}_{\mathrm{d}}$ of the QM7b results. We use the Matérn 3/2 kernel instead of the Matérn 5/2 kernel for the case of $\epsilon^{\mathrm{ML}}_{\mathrm{o}}$ for QM7b results, as it was empirically found to yield slightly better accuracy. 222In principle, the smoothness of the Matérn kernel could be taken as a kernel hyperparameter; however, this possibility was not explored in this work.

Feature selection is performed using the random forest regression implementation in the Scikit-learn v0.20.0 package. Pedregosa et al. (2011)

III Results

The ML model of Eq. 6 is a universal functional for any molecular Hamiltonian. In principle, with an adequate feature list and unlimited training data (and time), it should accurately and simultaneously describe all molecular systems. In practice, we must train the ML model using a truncated feature list and finite data. These choices determine the accuracy of the model.

Below, we examine the performance of the MOB-ML method in three increasingly broad regions of chemical space: (i) training on randomized water molecule geometries and predicting the energies of other water molecule geometries; (ii) training on geometries of short alkanes and predicting the energies of longer alkanes, and (iii) training on a small set of organic molecules and predicting the energies of a broader set of organic molecules. The first two test cases were introduced in our previous work,Welborn, Cheng, and Miller III (2018) and we explore how the refined methodology reported herein leads to improvements in accuracy and transferability. The last case represents a demanding new test of transferability across chemical space. In all cases, we report the ML prediction accuracy as a function of the number of training examples.

As a first example, we consider the performance of MOB-ML for a single water molecule. A separate model is trained to predict the correlation energy at the MP2, CCSD, and CCSD(T) levels of theory, using reference calculations on a subset of 1000 randomized water geometries to predict the correlation energy for the remainder. Feature selection with an importance threshold of $1\text{\times}{10}^{-3}$ results in 12, 11 and 10 features for $\epsilon^{\mathrm{ML}}_{\mathrm{o}}$ for MP2, CCSD and CCSD(T), respectively; ten features are selected for $\epsilon^{\mathrm{ML}}_{\mathrm{d}}$ for all three post-Hartree-Fock methods.

Figure 1 presents the test set prediction accuracy of each MOB-ML model as a function of the number of training geometries (i.e., the “learning curve"). MOB-ML predictions are shown for MP2, CCSD, and CCSD(T), and the model shows the same level of accuracy for all three methods. Remarkably, all three models achieve a prediction mean absolute error (MAE) of 1 mH when trained on only a single water geometry, indicating that only a single reference calculation is needed to provide chemical accuracy for the remaining 999 geometries at each level of theory. Since it contains 10 distinct LMO pairs, this single geometry provides enough information to yield a chemically accurate MOB-ML model for the global thermally accessible potential energy surface.

For all three methods (Fig. 1), the learning curve exhibits the expectedCortes et al. (1994) power-law behavior as a function of training data, and the total error reaches microhartree accuracy with tens of water training geometries. As compared to our previous results, where training on 200 geometries resulted in a prediction MAE of 0.027 mH for the case of CCSD,Welborn, Cheng, and Miller III (2018) the current implementation of the MOB-ML model is substantially improved; the improvement for this case stems primarily from the use of Boys localization,Boys (1960) which specifies unique and consistent LMOs corresponding to the oxygen lone pairs.

Next, we explore the transferability of MOB-ML predictions for a model that is trained on thermalized geometries of short alkanes and then used for predictions on thermalized geometries of larger and more branched alkanes (n-butane and isobutane). For these predictions, the absolute zero of energy is shifted for each molecule to compare relative energies on its potential energy surface (i.e., parallelity errors are removed). These shifts are reported in the figure caption; for no other results reported in the paper are parallelity errors removed.

In our previous work,Welborn, Cheng, and Miller III (2018) this test was performed using training data that combined of 100 geometries of methane, 300 of ethane, and 50 or propane; the resulting predictions are reproduced here in Fig. 2a. This earlier implementation of MOB-ML led to predictions for n-butane and isobutane with substantial errors (0.59 mH for n-butane and 0.93 mH for isobutane) and noticable skew with respect to the true correlation energy.

The predictions of MOB-ML in the current work (Fig. 2b) are markedly improved. First, the overall prediction accuracy is improved for all four summary statistics (inset in Fig. 2) despite substantial reduction in the number of training examples used. (The current work uses only 50 geometries of ethane, 20 geometries of propane, and no methane data.) Second, n-butane and isobutane are predicted with nearly identical accuracy. Finally, the prediction errors are no longer skewed as a function of true correlation energy. The primary methodological sources of these improvements are found to be symmetrization of occupied orbitals (Eq. 7) and the improved feature selection methodology. The MOB-ML features in the current work are selected with an importance threshold of $1\text{\times}{10}^{-4}$, resulting in 27 features for $\epsilon^{\mathrm{ML}}_{\mathrm{d}}$ and 12 features for $\epsilon^{\mathrm{ML}}_{\mathrm{o}}$; results presented in Fig 2b for CCSD(T) are qualitatively identical to those obtained for CCSD (not shown).

We now examine the transferability of the MOB-ML method across a broad swath of chemical space. Specifically, we consider the QM7b dataset,Montavon et al. (2013b) which is comprised of 7,211 plausible organic molecules with 7 or fewer heavy atoms. The chemical elements in QM7b are limited to those likely to be found in drug-like compounds: C, H, O, N, S, and Cl. We refer to the dataset used herein as QM7b-T to reflect the fact that it contains geometries sampled at a temperature of 350 K (as described in Sec. II), as opposed to DFT-optimized geometries. The MOB-ML model is trained on a randomly chosen subset of QM7b-T molecules and used to predict the correlation energy of the remainder. Active learning was also tested as a training data selection strategy, but was not found to improve the predictions in the regime of chemical accuracy, and in fact led to slightly worse transferability.

For comparison, a $\Delta$-ML modelRamakrishnan et al. (2015) was trained on the same molecules using kernel-ridge regression using the FCHL representationFaber et al. (2018) with a Gaussian kernel function (FCHL/$\Delta$-ML), as implemented in the QML package.Christensen, Faber, and von Lilienfeld (2019) All hyperparameters of the model were set to those obtained in Ref. 70, which have previously been demonstrated to work well for datasets containing structures similar to those in QM7b-T.Christensen, Faber, and von Lilienfeld (2019)

A possible source of concern for MOB-ML is that the number of selected features would grow with the chemical complexity of the training data. For example, 27 features for $\epsilon^{\mathrm{ML}}_{\mathrm{d}}$ and 12 features for $\epsilon^{\mathrm{ML}}_{\mathrm{o}}$ were selected in the alkane test case using ethane + propane training data (Fig. 2b), whereas only 10 features for $\epsilon^{\mathrm{ML}}_{\mathrm{d}}$ and 10 features for $\epsilon^{\mathrm{ML}}_{\mathrm{o}}$ were selected for the water test case at the CCSD(T) level of theory (Fig. 1). To examine this, we perform feature selection on increasing numbers of randomly selected molecules from the QM7b-T dataset. Table 1 presents two statistics on the feature importance as a function of the number of training molecules: (i) the number of “important features" (i.e., those whose permutation importanceBreiman (2001b) exceeds a set threshold of $2\text{\times}{10}^{-4}$ and $5\text{\times}{10}^{-5}$ for $\epsilon^{\mathrm{ML}}_{\mathrm{d}}$ and $\epsilon^{\mathrm{ML}}_{\mathrm{o}}$, respectively) and (ii) the inverse participation ratioKramer and MacKinnon (1993) of the feature importance scores. The latter is a threshold-less measure of the number of important features; it takes a value of 1 when only 1 feature has nonzero importance and $N$ when all $N$ features have equal importance. Although the QM7b-T dataset contains many different chemical elements and bonding motifs, Table 1 reveals that the selected features remain compact and do not grow with the number of training molecules. Indeed, for a large number of training molecules, the number of selected features slightly decreases, reaching 42 and 24 selected features for $\epsilon^{\mathrm{ML}}_{\mathrm{d}}$ and $\epsilon^{\mathrm{ML}}_{\mathrm{o}}$, respectively, for the largest training sizes considered.

The learning curves for MOB-ML models trained at the MP2/cc-pVTZ and CCSD(T)/cc-pVDZ levels of theory are shown in Fig. 3a, as well as the FCHL/$\Delta$-ML learning curve for MP2/cc-pVTZ. At the MP2 level of theory, the MOB-ML model achieves an accuracy of 2 mH with 110 training calculations (representing 1.5% of the molecules in the QM7b-T dataset), whereas the FCHL/$\Delta$-ML requires over 300 training geometries to reach the same accuracy threshold. Fig. 3a also illustrates the relative insensitivity of MOB-ML to the level of electronic structure theory, with the learning curve for CCSD(T)/cc-pVDZ reaching 2 mH accuracy with 140 training calculations. An analysis of the sensitivity of the MOB-ML predictions to the number of selected features is presented in supporting information Fig. S1, which indicates that the reported results are robust with respect to the number of selected features.

As a final test of transferability of the MOB-ML and FCHL/$\Delta$-ML methods across chemical space, Figs. 3b and 3c show results in which the ML methods are trained on QM7b-T molecules and then used to predict results for a dataset of 13-heavy-atom organic molecules at thermalized geometries, GDB-13-T, which includes six thermally sampled geometries each of 1,000 13-heavy-atom organic molecules chosen randomly from the GDB-13 dataset.Blum and Reymond (2009) Like QM7b, the members of GDB-13 contain C, H, N, O, S, and Cl. The size of these molecules precludes the use of coupled cluster theory to generate reference data; we therefore make comparison at the MP2/cc-pVTZ level of theory, noting that MOB-ML has consistently been shown to be insensitive to the employed post-Hartree-Fock method (as in Fig. 3a). Transfer learning results as a function of the number of training molecules are presented in Figs. 3b (on a linear-linear scale) and 3c (on a log-log scale).

Using the MOB-ML model that is trained on 110 seven-heavy-atom molecules (corresponding to a prediction MAE of 1.89 mH for QM7b-T), we observe a prediction MAE of 3.88 mH for GDB-13-T. Expressed in terms of size-intensive quantities, the prediction MAE per heavy atom is 0.277 mH and 0.298 mH for QM7b-T and GDB-13-T, respectively, indicating that the accuracy of the MOB-ML results are only slightly worse when the model is transferred to the dataset of larger molecules. On a per-heavy-atom basis, MOB-ML reaches chemical accuracy with the same number of QM7b-T training calculations (approximately 100), regardless of whether it is tested on QM7b-T or GDB-13-T.

In contrast with MOB-ML, the FCHL/$\Delta$-ML method is found to be significantly less transferable from QM7b-T to GDB-13-T. For models trained using 100 seven-heavy-atom molecules, the MAE per heavy atom of FCHL/$\Delta$-ML is over twice that of MOB-ML (Fig. 3b). Moreover, whereas MOB-ML reaches the per-heavy-atom chemical accuracy threshold with 140 training calculations, the FCHL/$\Delta$-ML method only reaches that threshold with 5000 training calculations.

IV Conclusions

Molecular-orbital-based machine learning (MOB-ML) has been shown to be a simple and strikingly accurate strategy for predicting correlated wavefunction energies at the cost of a Hartree-Fock calculation, benefiting from the intrinsic transferability of the localized molecular orbital representation. The starting point for the MOB-ML method is a rigorous mapping from the Hartree-Fock molecular orbitals to the total correlation energy, which ensures that the use of sufficient training data and molecular orbital features will produce a model that matches the corresponding correlated wavefunction method across the entirety of chemical space. The current work explores this possibility within the subspace of organic molecules. It is shown that MOB-ML predicts energies of the QM7b-T dataset to within a 2 millihartree accuracy using only 110 training calculations at the MP2/cc-pVTZ level of theory and using 140 training calculations at the CCSD(T)/cc-pVDZ level of theory. Direct comparison with FCHL/$\Delta$-ML reveals that MOB-ML is threefold more efficient in reaching chemical accuracy for describing QM7b-T. Further, a transferability test of a MOB-ML model trained on QM7b-T to GDB-13-T reveals that MOB-ML exhibits negligible degradation in accuracy; as a result, chemical accuracy is achieved with 36-times fewer training calculations using MOB-ML versus FCHL/$\Delta$-ML. These results suggest that MOB-ML provides a promising approach toward the development of density matrix functionals that are applicable across broad swathes of chemical space.

Acknowledgements.

We thank Daniel Smith (Molecular Sciences Software Institute) and Alberto Gobbi (Genentech) for a helpful discussion about available training datasets. T.F.M. acknowledges support from AFOSR award no. FA9550-17-1-0102. A.S.C. acknowledges support from the National Centre of Competence in Research (NCCR) Materials Revolution: Computational Design and Discovery of Novel Materials (MARVEL) of the Swiss National Science Foundation (SNSF). We additionally acknowledge support from the Resnick Sustainability Institute postdoctoral fellowship (M.W.) and the Camille Dreyfus Teacher-Scholar Award (T.F.M.). Computational resources were provided by the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the DOE Office of Science under contract DE-AC02-05CH11231.

Supporting Information

The datasets used in this work are available for download;Cheng et al. (2019) they include MOB-ML features, HF energies, pair correlation energies, and geometries. MOB-ML and FCHL/$\Delta$-ML predictions corresponding to Fig. 3 and an analysis of the sensitivity of the results of Fig. 3 to the number of selected features are available in the attached supporting information. A reference Python implementation for generating MOB-ML features is available at https://github.com/thomasfmiller/MOB-ML.