MOB-ML

Reading notes on T. Husch, J. Sun, L. Cheng, S. J. R. Lee, and T. F. Miller III, “Improved accuracy and transferability of molecular-orbital-based machine learning,” J. Chem. Phys. 154, 064108 (2021).

Abstract

Molecular-orbital-based machine learning (MOB-ML) predicts accurate correlation energies at the cost of obtaining molecular orbitals from a Hartree-Fock (HF) calculation. MOB-ML is a $\Delta$-ML method, meaning it learns the difference between a cheap HF calculation and an expensive post-HF method such as CCSD(T). The method builds on Nesbet’s theorem, which decomposes the total correlation energy into a sum of orbital-pair contributions. This decomposition reframes what would otherwise be an extrapolation task, training on small molecules and predicting for large ones, into an interpolation task, since the pair energies depend on local orbital features rather than on global molecular size.

Introduction

The overarching goal is the efficient calculation of accurate potential energies for molecular systems. Density functional theory (DFT) is widely used but suffers from a lack of systematic improvability and uncertain quality across applications. By contrast, wavefunction-based methods (MP2, CCSD(T), and so on) are systematically improvable, but their steep computational scaling limits routine use.

MOB-ML learns the sum of orbital-pair correlation energies, where each pair energy $\epsilon_{ij}$ is regressed against an individual feature vector $\mathbf{f}_{ij}$ representing the interacting molecular orbital (MO) pair $(i, j)$. The total correlation energy is recovered by summing over all pairs,

$E^{\text{corr}} = \sum_{ij} \epsilon_{ij}.$

This is Nesbet’s theorem. For single-reference methods the correlation energy decomposes exactly into these pair contributions, where $i$ and $j$ index occupied orbitals.

Because each pair energy depends only on local orbital information, the learning problem becomes size-independent once the feature vector is sufficiently descriptive, which converts the extrapolation problem into an interpolation one.

Theory

MOB-ML as a $\Delta$-learning model. MOB-ML learns the map

$\epsilon_{ij} \approx \epsilon^{\text{ML}}(\mathbf{f}_{ij}),$

where $\mathbf{f}_{ij}$ is a feature vector encoding information about the molecular orbital pair $(i, j)$. As a $\Delta$-ML model, MOB-ML predicts the difference between a cheap (HF) and an expensive (e.g. CCSD(T)) calculation. This formulation is generally more data-efficient than learning the total energy directly, because the correlation energy is a smoother, smaller quantity than the total energy. $\Delta$-ML models are common in computational chemistry, the idea being that the ML model only needs to learn the residual that the baseline method misses.

Gerade/ungerade rotation and canonical labeling. Given a pair of localized occupied orbitals $\phi_i$ and $\phi_j$, the assignment of labels $i$ and $j$ is arbitrary. To resolve this ambiguity, the orbitals are rotated into gerade (symmetric) and ungerade (antisymmetric) combinations (gerade and ungerade come from the German for “even” and “odd”),

$\tilde{\phi}_i = \frac{1}{\sqrt{2}}(\phi_i + \phi_j), \qquad \tilde{\phi}_j = \frac{1}{\sqrt{2}}(\phi_i - \phi_j).$

This defines a unique canonical labeling; $\tilde{\phi}_i$ is always the gerade combination and $\tilde{\phi}_j$ is always the ungerade one, regardless of how $i$ and $j$ were originally assigned. The rotation is purely a device to make the labeling unambiguous, unrelated to any energy ordering. The rotation has two side effects. Some feature vector elements acquire a sign that depends on the arbitrary original assignment of $i$ versus $j$, which is addressed below under orbital-index permutation invariance. And certain off-diagonal feature vector elements no longer decay to zero at large orbital separation, which is addressed under size consistency via a damping function.

Construction of the feature vector. The feature vector $\mathbf{f}_{ij}$ assembles information on the molecular orbital interactions between the pair. Its components include the following.

• Orbital energies of valence-occupied and valence-virtual orbitals ($F_{pp}$).
• Mean-field interaction energies of valence-occupied and valence-virtual orbitals ($F_{pq}$, the off-diagonal Fock matrix elements).
• Coulomb interaction integrals of valence-occupied and valence-virtual orbitals ($[\mathbf{K}^{pq}]_{qq}$, two-electron repulsion integrals).
• Exchange integrals of valence-occupied and valence-virtual orbitals ($[\mathbf{K}^{pq}]_{pq}$).

All of these quantities enter the MP2 or MP3 correlation energy expressions, and the model is thus fed the same information that perturbation theory uses to approximate correlation.

Defining importance of feature vector elements. Careful ordering of the elements within each feature vector block is necessary because Gaussian process regression (GPR), the regression method used here, is sensitive to permutations of its input. GPR also requires fixed-length feature vectors, so a consistent ordering and truncation scheme is essential. Additionally, since only a limited number of orbital-pair interactions contribute significantly to the pair correlation energy (due to the near-sighted nature of dynamic electron correlation), a cutoff criterion is introduced to keep the feature vector at a fixed, manageable length.

For feature vector elements involving a single auxiliary index (e.g. the Coulomb integrals $[\mathbf{K}^{pq}]_{aa}$ for orbital $a$), the importance is determined by the spatial distance between the orbital centroids. For elements involving two auxiliary indices (most prominently the exchange integrals $\{[\mathbf{K}^{ab}]_{ab}\}$ between valence-virtual orbitals $a$ and $b$), the importance is instead defined via the MP3 energy approximation, where the element that appears at lowest order in perturbation theory is deemed most important. The importance of a feature is thus tied to its role in the perturbative expansion of the correlation energy rather than to a generic statistical criterion.

Orbital-index permutation invariance. The Fock, Coulomb, and exchange matrix elements that comprise MOB features are naturally invariant to rotation and translation of the molecule. However, the gerade/ungerade rotation introduces a sign ambiguity for energetically degenerate orbital pairs $i, j$, since certain feature vector elements of the form

$M_{jp} = \frac{1}{\sqrt{2}}(M_{ip} - M_{jp})$

flip sign when $i$ and $j$ are swapped. Here $M_{jp}$ stands generically for any feature vector element that could be $F_{jp}$, $[\mathbf{K}^{ij}]_{pp}$, or $[\mathbf{K}^{ip}]_{jp}$. To restore permutation invariance, the absolute value is taken,

$M_{jp} = \frac{1}{\sqrt{2}}|M_{ip} - M_{jp}|.$

The corresponding gerade (sum) combination $\frac{1}{\sqrt{2}}(M_{ip} + M_{jp})$ is already permutation-invariant and requires no modification.

Size consistency. Size consistency requires that the energy of two isolated molecules equals the sum of their individual energies at infinite separation. In MOB-ML, this translates to two conditions, namely that (i) the diagonal feature vectors $\mathbf{f}_i$ must not be affected by distant, non-interacting molecules, and (ii) the pair correlation energy $\epsilon_{ij}$ for off-diagonal pairs with orbitals on different molecules must vanish, i.e. $\epsilon^{\text{ML}}(\mathbf{f}_{ij}) = 0$ when $r_{ij} = \infty$. No explicit training on dimers is required, because size consistency emerges from the feature vector construction alone.

Three modifications ensure these conditions are met. First, the diagonal feature vector $\mathbf{f}_i$ contains blocks ($\{F_{ab}\}$, $\{[\mathbf{K}^{aa}]_{bb}\}$, $\{[\mathbf{K}^{ab}]_{ab}\}$) whose elements depend only on virtual orbitals $a$ and $b$, not on orbital $i$, and therefore do not naturally decay with the distance between $i$ and the virtual orbitals. These elements are multiplied by their estimated importance from the ranking above, which is distance-based, causing them to decay smoothly to zero as the virtual orbitals move far from orbital $i$. To equalize decay rates across blocks, the cube root of $[\mathbf{K}^{pp}]_{qq}$-type integrals is used as a normalization factor.

Second, for off-diagonal feature vectors $\mathbf{f}_{ij}$, the same importance-weighting is applied, but this alone is insufficient because the gerade/ungerade rotation mixes contributions from orbitals $i$ and $j$, so even after importance-weighting some elements remain nonzero at $r_{ij} \to \infty$. A damping function inspired by the semi-classical limit of the MP2 expression is applied to each off-diagonal feature vector element,

$\frac{1}{1 + \frac{1}{6}\!\left(\frac{r_{ij}}{r_0}\right)^{\!6}},$

with $r_0 = 5.0$ bohrs chosen empirically. The damping radius must be large enough not to interfere with short-range machine learning, but small enough that off-diagonal features are fully damped at long range. The value could be optimized systematically, but $r_0 = 5.0$ bohrs works throughout this study. This suppresses all off-diagonal features at large $r_{ij}$, ensuring $\mathbf{f}_{ij} \to \mathbf{0}$ as $r_{ij} \to \infty$.

Third, the constraint $\epsilon^{\text{ML}}(\mathbf{0}) = 0$ is enforced, so that a zero feature vector yields zero pair correlation energy. Together, these three modifications guarantee size consistency. For molecules $A$ and $B$ at infinite separation, all cross-molecule features vanish, cross-molecule pair energies are zero, and $E^{\text{corr}}(AB) = E^{\text{corr}}(A) + E^{\text{corr}}(B)$.

Computational Details

Five datasets are used to benchmark MOB-ML, namely (i) a series of alkane molecules (ethane, propane, butane, isobutane with about 1000 geometries each), (ii) the potential energy surface of malonaldehyde (50 thermalized structures), (iii) a thermalized version of the QM7b dataset (QM7b-T, approximately 7000 organic molecules with $\leq 7$ heavy atoms) and its larger counterpart GDB13-T (13 heavy atoms), (iv) backbone-backbone interactions (BBI, 98 pairs), and (v) a subset of mononuclear octahedral transition-metal complexes (TM-T). All molecules require a cheap HF calculation for feature generation, while only training molecules require an expensive MP2 or CCSD(T) calculation for labels. Test-set reference energies are computed solely for validation purposes; in a production setting, these would not be needed.

Features are generated from HF calculations using a cc-pVTZ basis for organic molecules and a def2-TZVP basis for transition metals. Occupied orbitals are localized using either the Boys-Foster scheme or intrinsic bond orbital (IBO) localization; IBO is preferred for datasets involving large geometric distortions, as it yields smoother descriptions of localized orbitals across conformational changes. Boys localization tends to yield more efficient learning curves for rigid molecules, while IBO localization is more robust for flexible ones such as malonaldehyde.

ML pipeline. Feature selection is performed first using a random forest regression to identify the most informative features. Diagonal and off-diagonal pair correlation energies are then regressed separately using GPR with a Matérn 5/2 kernel and white noise regularization. Kernel hyperparameters are optimized by minimizing the negative log marginal likelihood. Both random sampling and active learning strategies are used to select training molecules. In active learning, the GPR variance identifies which unlabeled molecule would be most informative, and only that molecule’s expensive correlated calculation is performed.

Results

Transferability within a molecular family. Training on ethane and propane correlation energies, MOB-ML predicts butane and isobutane to within 0.11 and 0.10 kcal/mol peak MARE, respectively, improving on previous MOB-ML results. Each GPR prediction comes with an uncertainty estimate. The variance is small for ethane and propane (the interpolation regime) but larger for butane and isobutane (extrapolation), confirming that the model knows when it is uncertain. Including butane in the training set reduces errors further, illustrating the interpolation-versus-extrapolation effect directly.

MOB-ML is also tested on the malonaldehyde proton-transfer surface. Trained only on basin structures (reactants and products), the model predicts the transition-state barrier to within 0.35 kcal/mol, without ever seeing transition-state geometries in training. The highest errors coincide with the highest GPR variance, as expected for an extrapolative regime.

Transferability across organic chemistry space. Training on QM7b-T and predicting GDB13-T demonstrates size-transferability. Only approximately 70 actively sampled QM7b-T molecules (about 1% of the dataset) are needed for sub-kcal/mol accuracy on the remaining 99%. Active learning matches or beats random sampling with fewer training molecules, approximately 70 versus approximately 100 to reach 1 kcal/mol MAE. When transferred to GDB13-T, the error per electron pair remains comparable to within-QM7b-T predictions; the total error is larger simply because bigger molecules have more pairs. This constant per-pair error is the direct consequence of the interpolation framing via Nesbet’s theorem, as the model sees similar orbital pairs regardless of molecular size. MOB-ML outperforms DeepHF and FCHL18 at comparable or smaller training set sizes. Chemically accurate relative conformer energies are obtained for QM7b-T when trained on only 30 randomly chosen molecules (0.4% of the dataset).

Non-covalent interactions. Extending a QM7b-T-trained model to the backbone-backbone interaction (BBI) dataset initially yields approximately 1 kcal/mol MAE with large GPR uncertainty, indicating an extrapolative regime. Active learning identifies just two BBI data points whose inclusion reduces the MAE to 0.28 kcal/mol for the remaining 98 pairs. This demonstrates the active learning strategy for extending MOB-ML to new regions of chemical space with minimal additional computational effort.

Transition-metal complexes. MOB-ML is applied to closed-shell Fe, Co, and Ni complexes with various ligands spanning weak-field to strong-field. The learning curve per valence-occupied orbital is comparable to QM7b-T, and no modifications to the MOB-ML framework are needed. MP2 is not expected to be quantitative for transition metals, so this experiment is a proof of concept for the learning efficiency of MOB-ML in this chemical domain rather than a claim of accuracy.

Conclusions

MOB-ML provides a general framework for learning correlation energies at the cost of molecular orbital generation. The work makes three main points. First, enforcing physical constraints (size consistency, permutation invariance, and energy invariances) in the feature vector design leads to substantial improvements in both accuracy and data efficiency. Second, the orbital-pair decomposition via Nesbet’s theorem converts the extrapolation problem of going from small molecules to large molecules into an interpolation problem over similar orbital pairs across molecular sizes, enabling transferability across chemical space. Third, the GPR uncertainty estimates provide a natural active learning strategy for extending MOB-ML to new chemical domains with minimal additional reference calculations. Future directions include open-shell systems and electronically excited states, regimes where the single-reference decomposition underlying Nesbet’s theorem may require modification.