The multicomponent orbital-optimized second-order Møller–Plesset perturbation theory (OOMP2) method is the first multicomponent MP2 method that is able to calculate qualitatively accurate protonic densities, protonic affinities, and geometrical changes due to nuclear quantum effects in multicomponent systems. In this study, two approximations of the multicomponent OOMP2 method are introduced in an effort to demonstrate that, in orbital-optimized multicomponent methods, performing the orbital-optimization process with only electron–proton correlation is sufficient to obtain accurate protonic properties. Additionally, these approximations should reduce the computational expense of the multicomponent OOMP2 method. In the first approximation, the first-order wave function is written as a linear combination of one-electron one-proton excitations rather than as a linear combination of one-electron one-proton and two-electron excitations as in the original multicomponent OOMP2 method. Electron–electron correlation is included perturbatively after the orbital-optimization procedure has converged. In the second approach, the first approximation is further modified to neglect all terms in the orbital-rotation gradients that depend on the two-electron molecular-orbital integrals, which, assuming a fixed-sized protonic basis set, reduces the computational scaling for the orbital-optimization iterations to $Ne3$, where *N*_{e} is a measure of the electronic system size, compared to the $Ne5$ scaling of the original multicomponent OOMP2 method. The second approximation requires one $Ne5$ step after orbital convergence to compute the electron–electron correlation energy. The accuracy of the calculated protonic densities, protonic affinities, and optimized geometries of these approximations is similar or improved relative to the original multicomponent OOMP2 method.

## I. INTRODUCTION

Recent studies^{1–7} have consistently demonstrated the utility of electronic and protonic orbitals different from those of multicomponent Hartree–Fock (HF) in multicomponent methods^{8–12} culminating with the multicomponent orbital-optimized second-order Møller–Plesset perturbation theory (OOMP2) method.^{6} Single-component MP2^{13–15} is a workhorse method in quantum chemistry, but without orbital optimization, multicomponent MP2-based methods are qualitatively inaccurate^{16,17} for all properties related to the quantum proton. Even with the extra computational expense incurred by orbital optimization, multicomponent OOMP2 has the potential to be transformative in the field of multicomponent methods due to it having the same balance of accuracy and computational expense as single-component MP2.

The essential need for orbital optimization in multicomponent methods is because of the qualitative failure of multicomponent HF,^{8,10} which leads to a bad set of reference orbitals for post-HF multicomponent calculations. The multicomponent HF orbitals also introduce significant multireference character into the wave function,^{7} which would seem to indicate the need for multireference methods. However, by using orbitals different from those of multicomponent HF, as has been done in multicomponent density functional theory (DFT),^{1–3} coupled-cluster theory (CC),^{5,6} perturbation theory,^{6} and configuration interaction (CI),^{7} it appears that such multireference character can be significantly reduced.

In our opinion, an underappreciated aspect of multicomponent methods is that while they possess the same formal scaling as their corresponding single-component methods, in practice, their computational expense can be much greater due to either an increase in the number of wave function parameters or the need to perform many iterations in orbital-optimized methods to converge both the protonic and electronic orbitals.^{11,18} As an example, in our tests, a multicomponent DFT calculation on the FHF^{−} molecule using the B3LYP^{19–23} electron–electron and epc17-1 electron–proton correlation functionals^{2} and the aug-cc-pVTZ electronic^{24,25} and 8*s*8*p*8*d* protonic basis sets took 281 electronic and protonic SCF cycles to converge, while the corresponding single-component DFT calculation took only 19 SCF cycles. This is also seen in the multicomponent OOMP2 method. In our code, the multicomponent OOMP2 method takes 71 combined electronic and protonic iterations to converge the energy to 10^{−8} Ha for a calculation on the FHF^{−} molecule with the aug-cc-pVTZ electronic and 8*s*8*p*8*d*8*f* protonic basis sets, while the corresponding single-component OOMP2 calculation takes only 10 electronic iterations to converge.

Because of the large number of SCF iterations required in multicomponent orbital-optimized methods relative to their single-component counterparts, it is essential that the iterations be performed as efficiently as possible. In this study, we explore two approximations to the original multicomponent OOMP2 method. In the first method, the first-order wave function is modified to be a linear combination of only one-electron one-proton excitations. Electron–electron correlation can then be included after the calculation has converged using the single-component MP2 method. In the second method, we further approximate the first method by neglecting all terms in the electronic orbital-rotation gradient that contain two-electron molecular-orbital (MO) integrals. The two-electron MO integrals then do not need to be calculated during the orbital-optimization steps, which should be the computational bottleneck of most multicomponent OOMP2 calculations.

We show that by neglecting electron–electron correlation during orbital optimization, the accuracy of the calculated proton densities, proton affinities, and geometric perturbations due to nuclear quantum effects is similar or improved relative to the multicomponent OOMP2 method. This is similar to how the NEO-SOS′-OOMP2 method gives better results than full multicomponent OOMP2^{6} due to the former method biasing the orbital-optimization procedure to focus on electron–proton correlation rather than treating electron–proton and electron–electron correlation equally, as in multicomponent OOMP2. Based on the results in this study, we conclude by hypothesizing that for all orbital-optimized multicomponent methods, it is sufficient to perform the orbital-optimization process with only electron–proton correlation included.

## II. THEORY

Throughout this article, the letters *i*, *j*, *k*, … refer to occupied electronic spin orbitals, the letters *a*, *b*, *c*, … refer to unoccupied electronic spin orbitals, and the letters *p*, *q*, *r*, … refer to arbitrary electronic spin orbitals. Protonic spin orbitals are defined analogously but with uppercase letters. We present the multicomponent OOMP2 equations for a single proton treated quantum mechanically. Generalization to a non-protonic particle or multiple protons is straightforward.

The multicomponent second-order Møller–Plesset correction to the energy can be obtained by minimizing the Hylleraas functional,^{26}

with respect to the multicomponent MP2 amplitudes **t**, where |Ψ_{0}⟩ is the multicomponent HF wave function; $\u01240$ is the multicomponent Fock operator; $V^$ is the perturbation operator, $V^=\u0124\u2212\u01240$; $\u0124$ is the multicomponent Hamiltonian; and *E*_{0} is the zeroth-order energy. In multicomponent MP2, the first-order wave function is written as a linear combination of two-electron excitations and one-electron one-proton excitations as

Here, $tabij$ and $taAiI$ are multicomponent MP2 electron–electron and electron–proton amplitudes,

where *ε*_{p} and *ε*_{P} are the energies of spin orbitals *p* and *P*, respectively.

In the multicomponent orbital-optimized MP2 method (OOMP2), the extended multicomponent Hylleraas functional,

is optimized with respect to the multicomponent MP2 amplitudes and the orbital-rotation coefficients **R**^{e} and **R**^{p} of the electronic and protonic orbitals, respectively. *E*_{HF} is the multicomponent HF energy. The orbital-rotation coefficients are defined as

where **c**^{e} and **c**^{p} are the MO coefficients of the electronic and protonic orbitals, respectively. The orbital-rotation coefficient matrices are anti-Hermitian matrices. The change in an electronic or protonic orbital upon rotation by exp(**R**^{e}) or exp(**R**^{p}), respectively, is

The multicomponent OOMP2 extended Hylleraas function is explicitly written as

where *h*^{e} and *h*^{p} are the one-particle electronic and protonic Hamiltonians, respectively,

and $Fpqe$ and $FPQp$ are elements of the electronic and protonic Fock matrices, respectively,

The spin orbitals satisfy the pseudo-canonical conditions,

and the MP2-like density blocks are

The orbital-rotation gradient is obtained by expanding the change in the multicomponent OOMP2 extended Hylleraas functional to first order and taking its derivative with respect to the orbital-rotation coefficient matrix element $Rpqe$ or $RPQp$. In the multicomponent OOMP2 method, only the occupied-virtual orbital-rotation elements are non-redundant and need to be optimized. The electronic orbital-rotation gradient is

The protonic orbital-rotation gradient is

A sufficient condition for convergence is that the orbital-rotation gradients are equal to zero.

Similar to the single-component OOMP2 method,^{27–31} in multicomponent OOMP2, Brillouin’s theorem no longer holds and single excitations should be included in the first-order wave function. Previous single-component OOMP2 studies have neglected these terms,^{30} as orbital optimization should minimize the importance of single excitations. Additionally, the inclusion of single excitations leads to poor convergence of the orbitals during the orbital-optimization procedure.^{30} Therefore, we do not include any single excitations in any of our single-component or multicomponent OOMP2 calculations.

We converge the orbitals using a DIIS scheme^{32,33} similar to that of Neese’s^{29,30} implementation of single-component OOMP2.^{27–31} In this scheme, a matrix *B*^{e} is constructed for the electronic system as

with a protonic *B*^{p} matrix defined analogously. Δ is a level-shift parameter, which for all calculations in this study was set to 1.0 hartree for the electronic *B*^{e} matrix and 0.005 hartree for the protonic *B*^{p} matrix. When the electronic orbitals have converged, the *B*^{e} matrix is diagonal, so a new set of orbitals can be obtained by diagonalizing this matrix. To accelerate convergence, the standard DIIS procedure was used to construct these ** B** matrices. In practice, for all multicomponent OOMP2 calculations, we alternated construction of the

*B*^{e}and

*B*^{p}matrices each iteration. Therefore, only one of either the electronic or protonic orbitals is updated each iteration.

Orbital optimization in the multicomponent OOMP2 method simplifies the calculation of first-order properties, as the first-order protonic density matrix is a sum of the multicomponent HF MO protonic density matrix and the protonic *D*′^{p} density matrix in Eq. (11). This differs from the case of standard multicomponent MP2, where the multicomponent perturbed MP2 equations for the orbital response need to be solved to obtain the relaxed protonic density matrix.

In the multicomponent OOMP2 method, the electronic and protonic orbitals are optimized with electron–electron and electron–proton correlation included. However, for all of the systems that have been previously studied with the multicomponent OOMP2 method, orbital-optimization is likely to be superfluous in a single-component framework, as none of the systems are open shell or contain transition metals. For multicomponent systems, it is the poor multicomponent HF description of the electron–proton interactions that necessitates the need for orbital-optimized methods, which leads us to hypothesize that in the multicomponent OOMP2 method, the electron–electron and electron–proton correlation can be separated and the orbital optimization performed only with electron–proton correlation included. The electron–electron correlation can then be included perturbatively after the electronic and protonic orbitals have converged during the modified multicomponent OOMP2 procedure in a manner similar to how the energy from single excitations can be included perturbatively in single-component OOMP2.^{29} Further evidence for this hypothesis is seen in the excellent performance of the NEO-SOS′-OOMP2 method^{6} in which the contributions to the energy due to electron–proton correlation are scaled by 1.2 to increase their importance relative to those due to electron–electron correlation.

Our modification of the multicomponent OOMP2 method defines the first-order wave function to be a linear combination of one-electron one-proton terms,

With such a first-order wave function, the electronic orbital-rotation gradients are

with protonic orbital-rotation gradients identical to Eq. (13). After convergence of the orbitals, the electronic MP2 energy is then included perturbatively using the single-component expression for the MP2 energy ignoring any single-excitation contributions. We call this modification the multicomponent OOMP2(ee) method.

As this modification is variational with respect to all proton-related parameters, the proton density is calculated identically to the multicomponent OOMP2 method. Calculating the exact electronic density requires solving the perturbed MP2 equations, but for multicomponent systems, the protonic density is normally of much more interest than the electronic density, so this is not an issue at present.

The multicomponent OOMP2(ee) method should be computationally more efficient than the multicomponent OOMP2 method because in the multicomponent OOMP2(ee) method, the electronic MP2-like density matrix, *D*^{e}, and the electronic MP2 energy are not computed each iteration. The electronic orbital-rotation gradient is also simplified with the elimination of all terms that formally scale $Ne5$, where *N*_{e} is a measure of the size of the electronic system. However, for large systems, the computational bottleneck of both the multicomponent OOMP2 and OOMP2(ee) methods is likely to be identical: the calculation of the two-electron MO integrals, as the electronic orbital-rotation gradient in the OOMP2(ee) method has terms that depend on these integrals. This is especially true considering that large electronic basis sets are necessary in multicomponent methods to calculate accurate quantum protonic properties. In both methods, these two-electron MO integrals must be recomputed each iteration.

Because of the large computational expense of constructing the two-electron MO integrals, we investigate an approximation to the multicomponent OOMP2(ee) method where all terms depending on the two-electron MO integrals are neglected in the electronic orbital-rotation gradient. The electronic orbital-rotation gradient is then

We call this approximation the multicomponent OOMP2(ee′) method. Except for the differences between Eqs. (16) and (17), all steps in the multicomponent OOMP2(ee) and OOMP2(ee′) methods are identical.

The multicomponent OOMP2(ee′) method only requires the calculation of the two-electron MO integrals a single time during the perturbative calculation of the electronic MP2 energy. The multicomponent OOMP2(ee′) method still has a formal scaling of *N*^{5} with respect to the total system size as it has terms that scale $Ne2Np3$ and $Ne3Np2$, where *N*_{p} is a measure of the size of the protonic system. However, for multicomponent calculations with a single quantum proton, which is typical for multicomponent methods, the number of proton MOs remains fixed as the size of the system increases. Assuming this, the orbital-optimization iterations in the multicomponent OOMP2(ee′) method scale $Ne3$ rather than $Ne5$ as in the multicomponent OOMP2 and OOMP2(ee) methods. The multicomponent OOMP2(ee′) method then requires one perturbative step after convergence that scales $Ne5$.

The multicomponent OOMP2(ee′) method is no longer variational with respect to all proton-related parameters. Therefore, calculation of the exact first-order properties requires the solution of the multicomponent perturbed MP2 equations. In our multicomponent OOMP2(ee′) calculations, we will compute the first-order properties using the unrelaxed protonic density matrix. As will be demonstrated, the multicomponent OOMP2(ee) and OOMP2(ee′) methods give nearly identical results even when the multicomponent OOMP2(ee′) method uses an unrelaxed protonic density matrix.

## III. RESULTS AND DISCUSSION

The multicomponent OOMP2, OOMP2(ee), and OOMP2(ee′) methods were implemented in a locally modified version of PySCF.^{34} Presently, this code is not optimized, so we focus on the relative accuracy of the methods rather than their timings with the assumption that, for large systems, the calculation of the two-electron MO integrals will be the bottleneck in any multicomponent OOMP2 calculation.

We test the multicomponent OOMP2(ee) and OOMP2(ee′) methods by computing the protonic density of the FHF^{−} molecule, the proton affinities of 12 molecules, and the minimum-energy F–F distance of the FHF^{−} molecule when nuclear quantum effects of the proton are included. These three tests have emerged as the standard benchmark calculations for new multicomponent methods.^{1,2,4–6}

All single-component calculations in this study were performed with the ORCA program.^{35,36} All single-component calculations used the aug-cc-pVTZ electronic basis set for all atoms. All multicomponent calculations that computed protonic densities and affinities used the aug-cc-pVTZ basis set for the non-quantum proton atoms and the aug-cc-pVQZ basis set for the quantum proton. For the multicomponent calculations that determined the equilibrium F–F distance of FHF^{−}, the aug-cc-pVTZ electronic basis set was used for all atoms as this electronic basis set gives results that agree with the previous multicomponent OOMP2 study.^{6}

The multicomponent calculations used an 8*s*8*p*8*d*8*f* even-tempered^{37} protonic basis set with the exponents ranging from $22$ to 32. The quantum proton electronic and protonic basis sets were centered at the location of the hydrogen atom in the single-component calculations.

For all multicomponent OOMP2 results, the single-component method used for comparison is the OOMP2 method, while for the multicomponent OOMP2(ee) and OOMP2(ee′) results, the single-component method used for comparison for the density error is OOMP2 as this allows comparisons with the multicomponent OOMP2 method and standard single-component MP2 for the proton affinities and the minimum-energy F–F distance of the FHF^{−} molecule as the multicomponent OOMP2(ee) and OOMP2(ee′) methods do not include electron–electron correlation in the orbital-optimization procedure.

For all multicomponent FHF^{−} protonic density calculations, the coordinates of the molecule were obtained from a single-component OOMP2 optimization with the molecule constrained on the z-axis and with the classical hydrogen atom at the origin. The error in the protonic density for a method, *α*, is computed as the root mean-squared deviation relative to the protonic density calculated using the Fourier-grid Hamiltonian (FGH) method^{38,39} on a set of grid points. Mathematically, this is defined as

where $\rho i\alpha $ is the protonic density at point *i* with method *α* and *N*_{grid} is the number of points in the grid.

The FGH method solves the nuclear Schrödinger equation on a grid and is numerically exact for the adiabatic systems in this study. The grid in this study is a direct product of 32 uniformly spaced points in the x-, y-, and z-directions with the grid ranging from −0.5610 Å to 0.5984 Å in each dimension for a total of 32 768 points. For the FGH calculation, the fluorine atoms are fixed at −1.145 Å and 1.145 Å on the z-axis.

The minimum-energy F–F distance of FHF^{−} is found by performing a series of calculations at fixed F–F distances with the FHF^{−} molecule constrained to be on the z-axis with the fluorine atoms equidistant from the origin. For all of these multicomponent calculations, the quantum proton basis functions were centered at the origin.

The proton affinity of a molecule *A* is calculated as the difference between the energy of molecule *A* calculated with a single-component method and the energy of the molecule *AH*^{+} calculated with the corresponding multicomponent method. 5/2RT is added to the energy difference to account for the change in the total translational energy upon protonation and to define the proton affinities in terms of enthalpy rather than energy. By treating the proton quantum mechanically in the *AH*^{+} molecule, the zero-point vibrational energy of the proton including any anharmonic effects is explicitly included. The vibrational energy due to all other atoms is assumed not to change upon protonation, which has been shown to be a reasonable assumption in the previous multicomponent calculations.^{1,4–6} For the multicomponent proton affinity calculations, the geometry of the system was obtained from a single-component OOMP2 geometry optimization for the multicomponent OOMP2 calculations and from a single-component MP2 geometry optimization for the multicomponent MP2, OOMP2(ee), and OOMP2(ee′) calculations. For all of these multicomponent calculations, the quantum proton electronic and protonic basis functions were centered at the location of the classical proton in the single-component geometry optimization. For *AH*^{+} molecules with nonequivalent protons, the most acidic proton was treated quantum mechanically as the addition of this proton is what would be measured experimentally.

The absolute energies and density error for multicomponent OOMP2, OOMP2(ee), and OOMP2(ee′) methods are shown in Table I. The absolute energies of the corresponding single-component methods are also included in Table I. The energy changes by 11.6 mHa for both the multicomponent OOMP2 and OOMP2(ee) methods compared to their reference single-component method, while for the OOMP2(ee′) method, it changes by 11.5 mHa, indicating that all three methods treat the change in energy due to nuclear quantum effects similarly. The density error in the multicomponent OOMP2 method is slightly larger than in the previous multicomponent OOMP2 study,^{6} but this is attributed to the reference density FGH for the multicomponent OOMP2 methods being obtained using single-component orbital-optimized coupled-cluster doubles in that study rather than single-component OOMP2 as in this study.

. | Equilibrium F–F distance . | Energy . | ||||
---|---|---|---|---|---|---|

Multicomponent method . | Density error . | Multi . | FGH^{a}
. | Single . | Multi . | Single . |

OOMP2 | 0.397 | 2.313 | 2.312 | 2.290 | −200.182 31 | −200.193 88 |

OOMP2(ee) | 0.367 | 2.305 | 2.303 | 2.282 | −200.174 41 | −200.185 97 |

OOMP2(ee′) | 0.368 | 2.305 | 2.303 | 2.282 | −200.174 46 | −200.185 97 |

. | Equilibrium F–F distance . | Energy . | ||||
---|---|---|---|---|---|---|

Multicomponent method . | Density error . | Multi . | FGH^{a}
. | Single . | Multi . | Single . |

OOMP2 | 0.397 | 2.313 | 2.312 | 2.290 | −200.182 31 | −200.193 88 |

OOMP2(ee) | 0.367 | 2.305 | 2.303 | 2.282 | −200.174 41 | −200.185 97 |

OOMP2(ee′) | 0.368 | 2.305 | 2.303 | 2.282 | −200.174 46 | −200.185 97 |

^{a}

Reference 6.

Table I shows that the multicomponent OOMP2(ee) and OOMP2(ee′) methods have a slightly lower density error than the multicomponent OOMP2 method. This is hypothesized to be because the orbital-optimization procedure in the multicomponent OOMP2 method includes both electron–electron and electron–proton correlation, while the orbital-optimization procedure in the multicomponent OOMP2(ee) and OOMP2(ee′) methods only includes electron–proton correlation, allowing these methods to focus better on optimizing the protonic orbital. Overall, the density error is similar for all three methods, indicating that the multicomponent OOMP2(ee) and OOMP2(ee′) methods accurately reproduce the multicomponent OOMP2 method.

The results for the proton affinity calculations are presented in Table II. All multicomponent methods that include orbital optimization of the protonic orbitals perform significantly better than the multicomponent MP2 method without any orbital optimization. Inclusion of the electron–electron correlation energy as a perturbative correction after convergence of the orbital optimization is shown to be valid with a mean absolute difference (MAD) for the 12 molecule test set relative to the experimental proton affinity for the multicomponent OOMP2(ee) and OOMP2(ee′) methods of 0.22 eV and 0.22 eV, respectively, compared to an MAD of 0.25 eV for the original multicomponent OOMP2 method. We note that for the NO_{2}^{−}, HCOO^{−}, and N_{2} molecules, we obtain a slightly different proton affinity for the multicomponent OOMP2 method than a previous study.^{6} At present, the reason for this difference is unknown.

Molecule . | Experiment^{a}
. | OOMP2 . | OOMP2(ee) . | OOMP2(ee′) . | MP2^{b}
. |
---|---|---|---|---|---|

CN^{−} | 15.31 | 0.29 | 0.27 | 0.27 | 0.37 |

NO_{2}^{−} | 14.75 | 0.31 | 0.29 | 0.29 | 0.37 |

NH_{3} | 8.85 | 0.21 | 0.19 | 0.19 | 0.28 |

HCOO^{−} | 14.97 | 0.34 | 0.26 | 0.26 | 0.35 |

HO^{−} | 16.95 | 0.42 | 0.32 | 0.32 | 0.41 |

HS^{−} | 15.31 | 0.31 | 0.31 | 0.31 | 0.41 |

H_{2}O | 7.16 | 0.22 | 0.20 | 0.20 | 0.28 |

H_{2}S | 7.31 | 0.14 | 0.15 | 0.15 | 0.24 |

CO | 6.16 | 0.04 | 0.05 | 0.05 | 0.14 |

N_{2} | 5.12 | 0.19 | 0.20 | 0.20 | 0.28 |

CO_{2} | 5.60 | 0.26 | 0.22 | 0.22 | 0.30 |

CH_{2}O | 7.39 | 0.21 | 0.22 | 0.22 | 0.29 |

MAD | 0.25 | 0.22 | 0.22 | 0.31 |

Molecule . | Experiment^{a}
. | OOMP2 . | OOMP2(ee) . | OOMP2(ee′) . | MP2^{b}
. |
---|---|---|---|---|---|

CN^{−} | 15.31 | 0.29 | 0.27 | 0.27 | 0.37 |

NO_{2}^{−} | 14.75 | 0.31 | 0.29 | 0.29 | 0.37 |

NH_{3} | 8.85 | 0.21 | 0.19 | 0.19 | 0.28 |

HCOO^{−} | 14.97 | 0.34 | 0.26 | 0.26 | 0.35 |

HO^{−} | 16.95 | 0.42 | 0.32 | 0.32 | 0.41 |

HS^{−} | 15.31 | 0.31 | 0.31 | 0.31 | 0.41 |

H_{2}O | 7.16 | 0.22 | 0.20 | 0.20 | 0.28 |

H_{2}S | 7.31 | 0.14 | 0.15 | 0.15 | 0.24 |

CO | 6.16 | 0.04 | 0.05 | 0.05 | 0.14 |

N_{2} | 5.12 | 0.19 | 0.20 | 0.20 | 0.28 |

CO_{2} | 5.60 | 0.26 | 0.22 | 0.22 | 0.30 |

CH_{2}O | 7.39 | 0.21 | 0.22 | 0.22 | 0.29 |

MAD | 0.25 | 0.22 | 0.22 | 0.31 |

^{a}

References 40–43.

^{b}

Single-component method for comparison is MP2.

The equilibrium F–F bond distances for the three multicomponent methods in this study are shown in Table I. All three methods predict a distance within 0.02 Å of the appropriate FGH value, which further demonstrates that all three methods are able to include electron–proton correlation successfully.

The OOMP2(ee′) results are nearly identical to those of the OOMP2(ee) method with a difference in the density error of only 0.001, and no difference in either the MAD of the proton affinities or the equilibrium F–F bond distances of the FHF^{−} molecule, which validates the approximation of ignoring all two-electron terms in the electronic orbital-rotation gradient in the multicomponent OOMP2(ee) method.

## IV. CONCLUSIONS

We have introduced two approximations to the multicomponent OOMP2 method. The first approximation, OOMP2(ee), defines the first-order wave function as a linear combination of one-electron one-proton excitations. The electron–electron correlation energy is then included perturbatively after convergence of the orbital-optimization. The second approximation, OOMP2(ee′), is identical to the first approximation, except that all terms in the orbital-rotation gradient depending on the two-electron MO integrals are neglected, which leads to orbital-optimization iterations that formally scale $Ne3$ for a fixed protonic basis set compared to the $Ne5$ scaling of the multicomponent OOMP2 and OOMP2(ee) methods. This second approximation introduces very little error as the multicomponent OOMP2(ee) and OOMP2(ee′) methods give similar results. The accuracy of the multicomponent OOMP2(ee) and OOMP2(ee′) methods for the properties calculated in this study is similar or better than the multicomponent OOMP2 method due to the orbital-optimization procedure focusing solely on including electron–proton correlation. Based on these results, we offer a general hypothesis that when using orbital-optimized multicomponent methods, the orbital-optimization process can be performed with only electron–proton correlation included, assuming that a separation of the electron–electron and electron–proton correlation can be realized. We additionally recommend that future multicomponent MP2 studies use the multicomponent OOMP2(ee) and OOMP2(ee′) methods for calculating protonic densities and proton affinities due to the increased computational efficiency of these methods relative to the multicomponent OOMP2 method.

## DATA AVAILABILITY

The data that support the findings of this study are available within this article.

## ACKNOWLEDGMENTS

The computations for this work were performed on the high-performance computing infrastructure provided by Research Computing Support Services and, in part, by the National Science Foundation under Grant No. CNS-1429294 at the University of Missouri, Columbia, MO.

## REFERENCES

^{−}, He, Li

^{+}, Be

^{++}usw