Molecular Dipole Moment Learning via Rotationally Equivariant Gaussian Process Regression with Derivatives in Molecular-orbital-based Machine Learning

This study extends the accurate and transferable molecular-orbital-based machine learning (MOB-ML) approach to modeling the contribution of electron correlation to dipole moments at the cost of Hartree-Fock computations. A molecular-orbital-based (MOB) pairwise decomposition of the correlation part of the dipole moment is applied, and these pair dipole moments could be further regressed as a universal function of molecular orbitals (MOs). The dipole MOB features consist of the energy MOB features and their responses to electric fields. An interpretable and rotationally equivariant Gaussian process regression (GPR) with derivatives algorithm is introduced to learn the dipole moment more efficiently. The proposed problem setup, feature design, and ML algorithm are shown to provide highly-accurate models for both dipole moment and energies on water and fourteen small molecules. To demonstrate the ability of MOB-ML to function as generalized density-matrix functionals for molecular dipole moments and energies of organic molecules, we further apply the proposed MOB-ML approach to train and test the molecules from the QM9 dataset. The application of local scalable GPR with Gaussian mixture model unsupervised clustering (GMM/GPR) scales up MOB-ML to a large-data regime while retaining the prediction accuracy. In addition, compared with literature results, MOB-ML provides the best test MAEs of 4.21 mDebye and 0.045 kcal/mol for dipole moment and energy models, respectively, when training on 110000 QM9 molecules. The excellent transferability of the resulting QM9 models is also illustrated by the accurate predictions for four different series of peptides.


I. INTRODUCTION
Applications of machine learning (ML) to electronic structure theory have grown rapidly with an increasing number of studies in a variety of chemical systems and applications , 1,2 such as directly predicting the molecular properties, developing force fields and interatomic potentials, and designing novel and efficient catalysts, 3,4 drugs 5-8 and materials. 9,10 Important applications of machine learning include predicting chemical properties directly to reduce computational costs via supervised learning, 11,12 detecting the patterns of chemical spaces via unsupervised learning, 13,14 and proposing more suitable chemical systems via reinforcement learning 7,15 and generative models . 16,17 Numerous approaches have been presented in the field of machine learning for electronic structure during the last decades to aid in the learning of molecular energies and other molecular properties .  Molecular-orbital-based machine learning (MOB-ML) 31,37,38,[46][47][48]59 is one such method that uses molecular orbital (MO) information from Hartree-Fock (HF) computation to create a simpler and more direct mapping from the MO-based (MOB) features to the correlation energies. By utilizing quantum level information as features, MOB-ML has shown great efficiency and transferability for highly accurate predicted molecular energies with few training data. The introduction of a general regression with a clustering learning framework and a scalable Gaussian process algorithm further scales up MOB-ML to learn a large amount of data. 48 The responses of energy to some external variables in the time-independent framework are frequently studied molecular properties 60,61 in ML for quantum chemistry, such as atomic force, electric or magnetic dipole moments, and electric or magnetic polarizability. The electric dipole moment (referred as dipole), defined as the first-order response of the molecular energy to an external electric field, is an essential component in computational spectroscopy.. 62 However, it remains challenging to efficiently and accurately compute dipole moments with density functional theory. 61,63 Several previous efforts have been made to model the dipole moment and other response properties via ML methods. 29,36,[50][51][52][53][54]58,64,65 MOB-ML also shows a great potential to model any molecular properties determined by electronic wavefunctions, including these response properties. For instance, highly-accurate atomic force predictions could be directly provided by the MOB-ML models trained on energy labels only. 47 In this work, we consider the application of MOB-ML to learn general time-independent linear response properties using dipole moments as an example. An accurate and transferable MOB-ML framework is introduced to learn the contribution of electron correlation to the dipole moment as a summation of pairwise contributions. These pair dipole moments can be learnt as the functions of MOB features and their derivatives to electric fields via a rotationally equivariant Gaussian process regression (GPR) algorithm, i.e., GPR with derivatives. 66 We also adapt the local GPR with the Gaussian mixture model (GMM) unsupervised clustering (GMM/GPR) algorithm 59 to scale up the training and reduce the training costs without sacrificing accuracy and transferability. 67 The accuracy and efficiency of the proposed MOB-ML approach are tested on various benchmark systems, including water, fourteen small molecules, the QM9 benchmark dataset, 68 and four series of peptides. 50 The literature methods to predict dipole moments and energies for QM9 datasets are also included for comparison.

II. THEORY AND METHODS
A. Review of MOB-ML to learn molecular energies MOB-ML is motivated by the Nesbet theorem 69 to learn the pairwise contribution of the correlation energy E corr as a universal function of occupied molecular orbitals (MO) represented by features computed from Hartree-Fock (HF): where ε i j is the pair energy corresponding to the occupied MO i and j, f ε i j is the MOB-ML feature vector for the MO pair i and j, and φ k is the set of all the MOs. The symmetrization of i, j orbitals is also adapted to ensure f ε i j = f ε ji . The feature vector is consisted of Fock matrix elements F pq and part of the electron repulsion integrals [κ pq ] mn = pq|mn : 37,46 which is referred as "energy feature set" in this study. The energy feature generation and sorting approach is the same as in Ref. 46. The reference pair energies can be calculated in most of post-HF methods, such as second-order Møller-Plesset perturbation theory 70 (MP2) and various coupled cluster theories, including CCSD and CCSD(T). 71,72 B. MO Decomposition of dipole moments in MOB-ML For a given system, the electric dipole moments µ can be expressed as the linear response of the energy E with respect to external electric field E E E , i.e., which could be further expressed as the sum of HF and correlation components: The correlation part µ corr can then be decomposed on pairs of occupied orbitals similar to Eq. 1: µ i j is referred to as pair dipole moments and is regressed by ML similar to Eq. 2. Compared with Eq. 2, we add the feature derivatives information ∇ E E E f i j (φ k ) as part of the features motivated by Figure 1 displays an example of the dipole moment decomposition on a water molecule to facilitate an understanding of this decomposition.
Example decomposition of dipole moments as a sum of pairwise MO contributions for a water molecule. The four vertices of the tetrahedron represent the self-interactions of four MOs (indexing as ii), and the six edges connecting the vertices represent the interactions between two MOs (indexing as i j). The pair dipoles µ ii and µ i j are shown in red and blue arrows with correct direction and scaling, respectively. Four MOs are also shown next to the corresponding vertex i. The relative length scaling is also shown using a grey line.

C. Feature design of dipole learning in MOB-ML
The energy feature set f ε i j (Eq. 3) includes enough information to model molecular energy and satisfies different invariance properties, including translational, rotational, and orbital permutational invariances. 37,46 To efficiently model the dipole moments, we additionally include the responses of feature vector to electric field E E E , i.e., However, the direct definition of ∇ E E E f ε i j does not satisfy translational invariance due to the dependence of the Fock matrix elements ∇ E E E F pq on the positional operator matrix element r pq .
where h pq are the one-electron Hamiltonian matrix elements, J k pq and K k pq are Coloumb and exchange matrix elements of the kth orbital, and r pq are the position operator matrix elements. The r pq · E E E and ∇ E E E r pq · E E E terms vanish when E E E = 0. However, the existence of r pq term in Eq. 10 results in a non-translataional invariant dipole feature design. Therefore, a redefinition of the derivatives of Fock matrix elements is adapted by subtracting the r pq term to make dipole feature vector (Eq. 8) translataional invariant.

D. Review of scalable GPR algorithms in MOB-ML
In MOB-ML, Gaussian process regression (GPR) is used to model the molecular energies with high accuracy and great transferability. 31,37,38,46 Gaussian process (GP) 66 describes a prior distribution of random functions such that, for any finite number of possible inputs x = {x 1 , . . . , x n }, the function values f (x) has a multivariate Gaussian distribution where K(x, x) is the kernel matrix. Assuming the training data (X, y) has a Gaussian distributed noise with variance σ 2 n (also referred as Gaussian likelihood), i.e., y ∼ N( f , Σ n ), where Σ n = σ 2 n I, the GPR prediction f (X ) for the test points is a multivariate Gaussian distribution with the mean (prediction) and variance (uncertainty) as Although GPR only requires a few hundred molecules to achieve highly-accurate models, it is difficult to increase the training size due to the high computational complexity (O(N 3 )) of full GPR. 37,46 To scale up the MOB-ML training, a regression-with-clustering framework is introduced by constructing a local regressor within each grouping of points. The resulting clusters from GMM agree with the chemically intuitive groupings of MO types, 59 and thus could be directly applied to facilitate the construction of local regression models of pair dipole moments. In addition, a scalable exact GPR regressor, termed Alternative Black-box Matrix-Matrix Multiplication (AltBBMM), 48 has been developed to lower the computational cost of each local regression. Ref. 59 accesses different combinations of clustering methods and local regressors and concludes the most efficient protocol is local regression by AltBBMM with GMM (GMM/GPR), which achieves the state-of-the-art accuracy on drug-like benchmark organic datasets.

E. GPR with derivatives for ML response and rotational equivariance
Several previous studies have recognized the importance of rotational equvariance for ML framework to efficiently learn molecular dipole moments. 51,54 For any rotational operatorÛ, a function g : From a physics perspective, the rotational equivariance guarantees that the predicted property will rotate correspondingly with the rotation of the system. Therefore, the rotational equivariance property is required for any tensorial molecular properties, such as force, dipole moment, and polarizability. In addition, the molecular energy is rotationally invariant, i.e., remains constant with the rotation of the system. The conditions for energy and dipole models can be formulated as following: For the MOB features, when applying a rotation operator U, the energy and dipole feature sets satisfy the relationship where U is the matrix representation ofÛ. Therefore, the MOB-ML energy model is always rotationally invariant for any regressor. However, it remains challenging and requires a special ML algorithm design to make the MOB-ML dipole model rotationally equivariant for a greater learning efficiency.
Assuming the training energy and dipole sets are respectively, we introduce a rotationally equivariant GPR with derivatives formalism that could accurately learn dipole moment and energy separately or simultaneously.
A single-task energy model could be directly learnt using the naive GPR in Eq. 12 with the prior distribution Since the derivative of a GP is also a GP, 66,73 dipole moments can be regressed by GPR with the prior distribution of µ ML (x) and the corresponding kernel matrix K µ could be written as following: where the superscripts of K ε represent the derivatives to the arguments, e.g.
including these derivatives in the dipole feature set is necessary to model µ ML . This mathematical deduction agrees with the physical intuition discussed in Sec. II C.
By using the Gaussian likelihood y µ ∼ N(µ ML , Σ µ ) with Σ µ = σ 2 µ I, for a set of test points X µ , the prediction mean and variance of the GPR with derivatives can be evaluated using Eq. 12 as where GPR with derivatives can be generalized to the multi-task learning of ε ML and µ ML simultaneously. In such case, their joint distribution is also a GP with the predictive mean and variance as Although there might be other rotationally equivariant GPR frameworks that provide models with similar accuracy, they might not ensure the learnt dipole model is a derivative of the energy model. As an analogy to the response in physics, it is desirable for the ML model of h and its response property model g to satisfy Eq. 22, termed as "ML response relationship".
This relationship requires model g to be conservative (or curl-free), i.e., ∇ × g = 0. 74,75 In this study, we apply this physically driven GPR with derivatives formalism to satisfy this response relationship. The detailed proofs to show that both single-task and multi-task GPR with derivatives satisfy the rotational equivariance (Eq. 14). Without any specification, we adapt the single-task framework of GPR with derivatives for all the following training. The performance comparisons of single-task and multi-task models are demonstrated in Sec. IV A.

A. Data generation
For water and fourteen small molecules, we sample 200 total configurations at 50 fs intervals from ab initio molecular dynamics (AIMD) trajectories using the ENTOS QCORE 76 software. The structures of the QM9 dataset and the four different series of peptides are directly obtained from Ref. 68 and Ref. 50, respectively. All the HF calculations and MOB-ML features generations are performed by the ENTOS QCORE software. The HF calculations are using the cc-pVTZ basis 77 and the cc-pVTZ-JKFIT density fitting basis. 78 The valenceoccupied and valence-virtual MOs are then localized using the Boys-Foster localization scheme. 46,79 The derivatives of the orbitals with respect to electric fields are calculated through the coupled-perturbed HF calculations implemented in the ENTOS QCORE software. According to the chain rule, the feature derivatives with respect to electric fields are then calculated from the orbitals and orbital derivatives.
The HF orbitals are then imported into the Molpro 2018.0 80 package via the matrop functionality to generate the local MP2 70 pair energies. The frozen core approximation is used in the local MP2 calculations. The derivatives of the pair energies with respect to electric fields are implemented and calculated in the Molpro package following Ref. 81.

B. Machine learning protocols
For all the datasets, we separately learn the energies and dipole moments using the energy feature set and dipole feature set, respectively. For the water and the small molecule datasets, the results for multi-task models, i.e., dipole + energy models, that learn both tasks simultaneously are also included for comparison. Table I summarizes the usage of energy and dipole features in this work. We unsupervisedly cluster the MOs represented by energy features instead of separately clustering energy and dipole feature space, i.e., the clustering models are identical for the energy and dipole learning with the same training sets. Feature selection is performed before all the models on energy labels using the random forest regression implementation in the SCIKIT-LEARN 82 package following the protocol in Ref. 37. For the dipole learning, we use the selected energy features and their derivatives as the selected dipole features. We apply the AltBBMM algorithm as the default GP regressor, and reimplement the GMM with full covariance matrix following SCIKIT-LEARN using CUPY 83 to enable multi-GPU training. The implementations for both algorithms are available online at https://github.com/SUSYUSTC/ BBMM.git For all GPR and GMM/GPR, we employ the Matérn 5/2 kernel with white noise regularization. 66 The parameters used in training for GPR and GMM/GPR are further discussed in the Supporting information.
All the results for water and small molecule datasets are collected with GPR without clustering. The local GPR with GMM unsupervised clustering is applied to scale the MOB-ML training in the QM9 dataset. In this work, we follow the same clustering protocol introduced in Ref. 59 to generate the GMM models. The GMM model is initialized by Kmeans clustering, and its number of clusters is automatically determined by minimizing the Bayesian information criterion. GMM could not be performed to cluster 50000 and 110000 QM9 molecules due to limited memory, and we thus apply the GMM model trained on 20000 QM9 molecules to approximate the clustering results of these two models. To reduce the learning costs, we also apply the same capping strategy described in Ref. 38. For the clusters containing a large number of points, we randomly select training points with the capping size defined ahead to regress these local GPR. The capping size is 1000000 and 300000 pairs for dipole and energy learning, respectively.

A. Dipole moment learning for small molecules via MOB-ML
To demonstrate the ability to learn dipole moments using the MOB representations, we first test the prediction accuracies of MOB-ML on water and other small molecules. Figure 2 displays the mean absolute errors (MAEs) of dipole moments (|µ|) and molecular energies E for water learnt by single-task and multi-task models. The sizes of the training set are varied from 2 to 100 geometries, and the test set is composed of 100 geometries not included in any training sets. MAEs at the same order of reference data accuracy for molecular energies are achieved by training on 100 geometries using energy labels (5.54 × 10 −5 kcal/mol) and dipole combined with energy labels (1.84 × 10 −5 kcal/mol). Across all the training sizes, single-task and multi-task models provide similar accuracies for molecular energies, but dipole models provide much better dipole predictions than the dipole + energy model. Table II lists the MAEs and maximum errors (Max) of dipole moments and total energies from single-task and multitask models training on 50 geometries and testing on different 50 geometries for small molecules with different sizes. MOB-ML provides very accurate predictions for all the test molecules. Comparing the results of molecules with different molecular sizes sharing similar MO properties, such as CH 4 , C 2 H 6 , and C 3 H 8 , it is clear that the larger molecules have much bigger errors than the smaller ones for both dipole and energy owing to the increasing number of pairwise contributions to the final result. The total errors should scale linearly with the increase of molecular size by summing up predicted pairwise contribution with Gaussian distributed pairwise errors. For the systems that share similar numbers of MOs, such as C 2 H 4 and C 2 H 6 , the more rigid molecule (C 2 H 4 ) is easier to learn for both dipole and energy.
Learning dipole moments and molecular energies simultaneously does not always provide better prediction accuracies  2. Prediction accuracies of dipole moments and energies of water using a reference theory of MP2/cc-pVTZ. The predicted MAEs are plotted versus the number of training geometries on a log-log scale ("learning curves"). The MOB-ML models for dipole moments and energies are constructed by training on the labels individually (dipole model and energy model) or simultaneously using multi-task learning (dipole + energy model). The primary and secondary y-axises represent the prediction MAE of dipole moments |µ| err in milli-Debye (mDebye) and the prediction MAE of energy in kcal/mol, separately. The results are also summarized in Table S1 in the Supporting Information (SI) of both dipole and energy for most of the molecules (12 out of total 14 molecules) and need higher computational costs since it trains more points within a model. This observation indicates that dipole vectors and energies of each pair of MOs might vary independently as functions of MOB features, and therefore no mutual supervision could be provided by multitask learning.
From a theoretical perspective, the ML response relationship could improve the learning efficiency of forces, but is not expected to enhance the learning efficiency of the dipole model with zero electric field. Therefore, the dipole + energy model for these small molecules cannot provide better accuracies than only training on dipoles. Furthermore, the inclusion of non-correlated tasks in the multi-task models further complicates the learning process and even leads to deterioration in training. Figure 3 could help explain this observation by using the water molecule as an example. The relative water MP2/cc-pVTZ total energy is plotted as a function of one of the O-H bond lengths and the strength of the external electric field along the bond direction. We fix the bond angel and only change one of the O-H bond lengths (red curve) with E = 0 to facilitate the understanding. The training set can be treated as samples of this simplified potential. We note that all the bond lengths and angles vary in the actual water dataset. From these data, the energy information can then directly infer the force information (black arrows) since the derivative with respect to bond length d could be estimated as differences between two sampled training points. However, since there is no change along the axis of the electric field E E E , no estimated information of dipole moments (blue arrows) is available from the training set.

B. MOB-ML for dipole moments and energies of organic molecules in QM9
In our previous studies, we have illustrated the excellent accuracy and transferability of MOB-ML to learn molecular energies using two thermalized organic molecule datasets, i.e., QM7b-T and GDB-13-T. 37,38,46,59 In this study, we systematically examine the learning performance of MOB-ML for both the dipole and energy using the benchmark organic chemistry dataset QM9, 68 which contains optimized structures of 133885 molecules with up to nine heavy atoms (HAs) of C, O, N, and F. QM9 is a standard benchmark dataset that has been assessed in many different literature studies. 29,36,41,[49][50][51][52][53][54][55][56][57][58] Figure 4 displays the predicted MAEs for dipole moments (in mDebye) and energies (in kcal/mol) as functions of number of training geometries on a log-log scale (learning curves). Since our GPR regression, AltBBMM, can only train at most 1 million points, we collect the results of MOB-ML (GPR) up to training on 1000 and 2000 dipole moments and molecular energies, respectively. The application of GMM/GPR scales the training of dipole moments and energy to the same training size (at most 110,000 QM9 molecules) as the literature models. The test sets of MOB-ML approaches remain the same across the entire learning curve with a size of 11,843 molecules. Different literature approaches computed at the B3LYP/6-31G(2df,p) level of theory are included for comparison.
Similarly, the prediction errors of energies from MOB-ML approaches are compared with SchNet, 29 SLATM, 56 PhysNet, 36 SOAP, 57 FCHL18, 55 OrbNet-Equi, 51 and QML 49 in Fig. 4b. MOB-ML (GPR) and MOB-ML (GMM/GPR) are still provide the best sets of results across all the training sizes. MOB-ML (GMM/GPR) achieves accuracies of 0.99 kcal/mol and 0.045 kcal/mol with only 100 and 110000 training molecules, respectively. Both numbers are the current best in this field. QML with an orbital based features and (∆-learning) (QML (MO, ∆-learning)) and OrbNet-Equi (∆learning) are other two most accurate approaches.
The top approaches to predict dipole moments and energies, i.e., MOB-ML, OrbNet-Equi (∆-learning), and QML (MO, ∆-learning), are further compared here from a theoretical perspective. To achieve the best accuracy, all the three approaches apply the idea of "∆-learning" by predicting the differences between low-level and high-level theories instead of directly predicting. In addition, all three approaches are orbital-based ML approaches that adapt features related to energy matrix elements (FJK features), and these features are considered to contain high-quality quantum-level information to make the ML map easier. Although the three approaches share several similarities, their differences might explain their prediction accuracy differences. Firstly, MOB-ML predicts the results from wavefunction theories using HF computations with the same basis, while OrbNet-Equi and QML (MO) predict the results from DFT using GFN-xTB and minimal basis HF computations, respectively. Since HF computed with cc-pVTZ basis set is more expensive and contains more accurate information about the orbitals than GFN-xTB and minimal basis HF, MOB-ML is more expensive in evaluation than the other two methods. MOB-ML explicitly decomposes the differences between low-and high-level theory results onto MOs and learns these pairwise contributions using GPR, while OrbNet-Equi and QML (MO) directly learn these differences by implicitly decomposing them to each kernel in ker-nel ridge regression (KRR) or nodes in graph neural network (GNN) by carefully designing the ML frameworks. This explicit decomposition brings an accuracy gain to MOB-ML but limits the applications of MOB-ML to decomposable properties. On the other hand, QML (MO) and OrbNet-Equi are able to predict more different molecular properties. energies. These models are collected on 8 NVIDIA Tesla V100-SXM2-32GB GPUs. Since each local GPR could be regressed independently on different GPUs, GMM/GPR is a highly-parallelized approach with excellent multi-GPU speedups. For both dipole and energy learning, it is clear that GMM/GPR provides much lower training costs compared with learning without clustering. Across all the training sizes that we could train directly with GPR, GMM/GPR could achieve more than 68.5 and 21.4 times speedups for dipole and energy learning, respectively, without loss of accuracy and transferability. For example, the best GMM/GPR energy model training on 110,000 molecules only takes 20.7 hrs, but it is over 7 times more expensive than the best dipole model (146.5 hrs). This is because the training points of dipole moment are 3 times larger than ones of energy for each molecule.

D. Predictions of dipole moments of four challenge cases
To illustrate the accuracy of MOB-ML in the actual biochemical systems, we further assess the prediction accuracy of the best MOB-ML (GMM/GPR) model on four different sets of peptides, termed as "challenging dataset". This challenging dataset is firstly introduced in Veit et al., 50 and in this study, we also included the literature results predicted by MuML model. 50 All the true and predicted dipole moments from different theories and ML models are plotted as a function of the chain length in Fig. 6. Table S3 in SI summarizes the predicted and true dipole moments and energies using the best QM9 energy model (GMM/GPR trained on 110,000 molecules). We note that MP2/cc-pVTZ is still an affordable theory for the QM9 benchmark dataset, but it is nearly impossible to obtain MP2 energy for the molecules with more than 30 heavy atoms, for instance, α-helix molecules with n = 7, 10, 13, 16 in the challenging dataset. Therefore, no true MP2/cc-pVTZ results are provided for these large α-helix molecules; meanwhile, MOB-ML provides reasonable dipole moment predictions for these molecules.
Except the large polyenonic amino acids in Fig. 6(c), MOB-ML provides nearly identical predicted dipole moments as MP2 for all other molecules in panels (a), (b), and (d), which indicates that the MOB-ML model for dipole moments has an excellent transferability to large molecules. MuML (combined) model also has a deteriorated accuracy for all the molecules in polyenoic amino acids. In this case, the MOB-ML model over-corrects the results from HF calculation, and the size of the errors increases dramatically with the increasing of sizes. The major contributors to the total dipole moments for both cases are the local polarization of the end groups, and therefore the total dipole moments should be nearly constant predictions with the increase of the molecular length. This observation agrees with the ML predictions and the true dipole moments of n-amino carboxylic acids, but not the true values of polyenoic amino acids. Veit et al. 50 assess the partial charges for each functional group carefully to understand the puzzling results of polyenoic amino acids. The partial charge distribution suggests that end groups in polyenoic amino acids carry a net positive and negative charge, leading to an increase of dipole moments with increasing molecular length.

V. CONCLUSION
In this study, we extend the MOB-ML framework to learn pairwise contributions of electron correlation part of dipole moments accurately and transferablely using the information computed from HF calculations. The feature derivatives are appended to the original energy features to form the dipole feature set to include more information in the MOB representation for dipole prediction. The introduction of GPR with derivatives algorithm leads to efficient and physical modeling of dipole moments by satisfying the significant properties of equivariance and ML response. For water and other small molecules, MOB-ML could provide more accurate predictions for the dipole moments and energies by learning two tasks separately than simultaneously. To generate a universal dipole model and energy model for organic molecules, we apply MOB-ML to the QM9 dataset and train the two sets of labels separately using the corresponding GPR training protocols. The GMM/GPR framework that constructs local scalable GPRs for clusters detected by GMM clustering approach is also applied to reduce the learning costs of MOB-ML and scale up its training to 110,000 molecules without loss of accuracy. MOB-ML (GMM/GPR) model could achieve accuracies of 4.21 mDebye and 0.045 kcal/mol by learning 110000 QM9 molecules for dipole moments and energies, respectively, and an accuracy of 0.99 kcal/mol by only training on 100 QM9 molecules for energies. MOB-ML is superior to all other literature results for both dipole moments and molecular energies and provides accurate and transferable results to most of the tested peptides with different three-dimensional structures. As a future study direction, this MOB feature design and GPR learning framework could be extended to model other response properties, such as polarizability and excited state properties. Furthermore, by taking advantage of the existence of the currently best MP2/cc-pVTZ models, multifidelity learning or ∆-learning 23 approaches could also be applied to learn the MOB-ML model at CCSD(T)/cc-pVTZ level by learning the differences between CCSD(T) and MP2 values.

ACKNOWLEDGMENTS
We thank Vignesh Bhethanabotla for his help in improving the quality of this manuscript. TFM acknowledges support from the US Army Research Laboratory (W911NF-12-2-0023), the US Department of Energy (DE-SC0019390), the Caltech DeLogi Fund, and the Camille and Henry Dreyfus Foundation (Award ML- . 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 reference pairwise decomposed dipole moments and energies and MOB features of QM7b-T, QM9, and GDB-13-T are available at Caltech Data: https://data.caltech. edu/records/1177. The corresponding HF dipole moments and energies are also available. The total dipole moment and energy data of challenging datasets are also included in this online dataset. The implementation of the multi-GPU Al-tBBMM and GMM are available online at https://github. com/SUSYUSTC/BBMM.git. The parameters of training the models are included in the Supporting information. Table S1, S2, and S3 list the MOB-ML results plotted in Fig. 2, 4, and 6, respectively.