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 Ne 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.

Recent studies1–7 have consistently demonstrated the utility of electronic and protonic orbitals different from those of multicomponent Hartree–Fock (HF) in multicomponent methods8–12 culminating with the multicomponent orbital-optimized second-order Møller–Plesset perturbation theory (OOMP2) method.6 Single-component MP213–15 is a workhorse method in quantum chemistry, but without orbital optimization, multicomponent MP2-based methods are qualitatively inaccurate16,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 B3LYP19–23 electron–electron and epc17-1 electron–proton correlation functionals2 and the aug-cc-pVTZ electronic24,25 and 8s8p8d 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 8s8p8d8f 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 OOMP26 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.

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 

EMP2=mint2Ψ1V^Ψ0+Ψ1Ĥ0E0Ψ1,
(1)

with respect to the multicomponent MP2 amplitudes t, where |Ψ0⟩ is the multicomponent HF wave function; Ĥ0 is the multicomponent Fock operator; V^ is the perturbation operator, V^=ĤĤ0; Ĥ is the multicomponent Hamiltonian; and E0 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

Ψ1=14ijabtabijΨijab+iIaAtaAiIΨiIaA.
(2)

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

tabij=ij||abεa+εbεiεj,taAiI=iI|aAεa+εAεiεI,
(3)

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,

Lt,Re,Rp=EHF+2Ψ1V^Ψ0+Ψ1Ĥ0E0Ψ1,
(4)

is optimized with respect to the multicomponent MP2 amplitudes and the orbital-rotation coefficients Re and Rp of the electronic and protonic orbitals, respectively. EHF is the multicomponent HF energy. The orbital-rotation coefficients are defined as

ce,new=expRece,old,cp,new=expRpcp,old,
(5)

where ce and cp 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(Re) or exp(Rp), respectively, is

expRei=i+aRaiea+,expRea=aiRaiei+,expRpI=I+ARAIpA+,expRpA=AIRAIpI+.
(6)

The multicomponent OOMP2 extended Hylleraas function is explicitly written as

Lt,Re,Rp=iihei+12ijij||ij+IIhpIiIiI|iI+ijabtabijijab2iIaAtaAiIiIaA+ijDijFije+abDabFabe+ijDijeFij+abDabeFab+IJDIJpFIJ+ABDABpFAB,
(7)

where he and hp are the one-particle electronic and protonic Hamiltonians, respectively,

he=12iNeleci2iNelecANclassZArierAc,hp=12mpINprotI2+INprotANclassZArIprAc,
(8)

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

Fpqe=pheq+ipi||qiIpI|qI,FPQp=PhpQiiP|iQ.
(9)

The spin orbitals satisfy the pseudo-canonical conditions,

Fije=δijFiie=δijεi,Fabe=δabFaae=δabεa,FIJp=δIJFIIp=δIJεI,FABp=δABFAAp=δABεA,
(10)

and the MP2-like density blocks are

Dije=12kabtabiktabjk,Dabe=12ijctacijtbcij,Dije=IaAtaAiItaAjI,Dabe=iIAtaAiItbAiI,DIJp=iaAtaAiItaAiJ,DABp=iIataAiItaBiI.
(11)

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

Lt,Re,RpRcke=gcke=2Fcke+2iabtabkiciba2ijatcaijijka2IaAtaAkIcIaA+2iIAtcAiIiIkA+2iDkieFcie2aDcaeFkae+2iDkieFcie2aDcaeFkae+2ijDijeic||jk+ik||jc+2abDabeac||bk+ak||bc+2ijDijeic||jk+ik||jc+2abDabeac||bk+ak||bc2IJDIJpcI|kJ+kI|cJ2ABDABpcA|kB+kA|cB.
(12)

The protonic orbital-rotation gradient is

Lt,Re,RpRCKp=gCKp=2FCKp2iaAtaAiKiCaA+2iaItaCiIiIaK+2IDKIpFCIp2ADCApFKAp2ijDijeiC|jK2abDabeaC|bK.
(13)

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 scheme32,33 similar to that of Neese’s29,30 implementation of single-component OOMP2.27–31 In this scheme, a matrix Be is constructed for the electronic system as

Bije=δijFiie,Babe=δabFaae+Δ,Baie=Biae=gaie,
(14)

with a protonic Bp matrix defined analogously. Δ is a level-shift parameter, which for all calculations in this study was set to 1.0 hartree for the electronic Be matrix and 0.005 hartree for the protonic Bp matrix. When the electronic orbitals have converged, the Be 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 Be and Bp 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 Dp 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 method6 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,

Ψ1=iIaAtaAiIΨiIaA.
(15)

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

Lt,Re,RpRcke=gcke=2Fcke2IaAtaAkIcIaA+2iIAtcAiIiIkA+2iDkieFcie2aDcaeFkae+2ijDijeic||jk+ik||jc+2abDabeac||bk+ak||bc2IJDIJpcI|kJ+kI|cJ2ABDABpcA|kB+kA|cB
(16)

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, De, 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 Ne 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

gcke=2Fcke2IaAtaAkIcIaA+2iIAtcAiIiIkA+2iDkieFcie2aDcaeFkae2IJDIJpcI|kJ+kI|cJ2ABDABpcA|kB+kA|cB.
(17)

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 N5 with respect to the total system size as it has terms that scale Ne2Np3 and Ne3Np2, where Np 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.

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 8s8p8d8f even-tempered37 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) method38,39 on a set of grid points. Mathematically, this is defined as

Density Error=iNgridρiαρiFGH2Ngrid,
(18)

where ρiα is the protonic density at point i with method α and Ngrid 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.

TABLE I.

Density error in atomic units, equilibrium F–F bond distances in Å for the respective multicomponent methods, FGH methods, and corresponding single-component methods, and absolute energies in hartree of the FHF molecule for the multicomponent OOMP2, OOMP2(ee), and OOMP2(ee′). For the multicomponent OOMP2 method, the corresponding single-component method is OOMP2. For the multicomponent OOMP2(ee) and OOMP2(ee′) methods, the corresponding single-component method for the density error is OOMP2, while for the equilibrium F–F bond distances and absolute energies, it is MP2.

Equilibrium F–F distanceEnergy
Multicomponent methodDensity errorMultiFGHaSingleMultiSingle
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 distanceEnergy
Multicomponent methodDensity errorMultiFGHaSingleMultiSingle
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 NO2, HCOO, and N2 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.

TABLE II.

Experimental proton affinities and absolute deviations relative to the experimental value of the calculated proton affinities of the multicomponent OOMP2, OOMP2(ee), OOMP2(ee′), and MP2 methods in eV. MAD is the mean absolute deviation average over all 12 molecules.

MoleculeExperimentaOOMP2OOMP2(ee)OOMP2(ee′)MP2b
CN 15.31 0.29 0.27 0.27 0.37 
NO2 14.75 0.31 0.29 0.29 0.37 
NH3 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 
H27.16 0.22 0.20 0.20 0.28 
H27.31 0.14 0.15 0.15 0.24 
CO 6.16 0.04 0.05 0.05 0.14 
N2 5.12 0.19 0.20 0.20 0.28 
CO2 5.60 0.26 0.22 0.22 0.30 
CH27.39 0.21 0.22 0.22 0.29 
MAD  0.25 0.22 0.22 0.31 
MoleculeExperimentaOOMP2OOMP2(ee)OOMP2(ee′)MP2b
CN 15.31 0.29 0.27 0.27 0.37 
NO2 14.75 0.31 0.29 0.29 0.37 
NH3 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 
H27.16 0.22 0.20 0.20 0.28 
H27.31 0.14 0.15 0.15 0.24 
CO 6.16 0.04 0.05 0.05 0.14 
N2 5.12 0.19 0.20 0.20 0.28 
CO2 5.60 0.26 0.22 0.22 0.30 
CH27.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.

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.

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

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.

1.
K. R.
Brorsen
,
Y.
Yang
, and
S.
Hammes-Schiffer
, “
Multicomponent density functional theory: Impact of nuclear quantum effects on proton affinities and geometries
,”
J. Phys. Chem. Lett.
8
,
3488
3493
(
2017
).
2.
Y.
Yang
,
K. R.
Brorsen
,
T.
Culpitt
,
M. V.
Pak
, and
S.
Hammes-Schiffer
, “
Development of a practical multicomponent density functional for electron-proton correlation to produce accurate proton densities
,”
J. Chem. Phys.
147
,
114113
(
2017
).
3.
K. R.
Brorsen
,
P. E.
Schneider
, and
S.
Hammes-Schiffer
, “
Alternative forms and transferability of electron-proton correlation functionals in nuclear-electronic orbital density functional theory
,”
J. Chem. Phys.
149
,
044110
(
2018
).
4.
F.
Pavošević
,
T.
Culpitt
, and
S.
Hammes-Schiffer
, “
Multicomponent coupled cluster singles and doubles theory within the nuclear-electronic orbital framework
,”
J. Chem. Theory Comput.
15
,
338
347
(
2018
).
5.
F.
Pavošević
and
S.
Hammes-Schiffer
, “
Multicomponent coupled cluster singles and doubles and Brueckner doubles methods: Proton densities and energies
,”
J. Chem. Phys.
151
,
074104
(
2019
).
6.
F.
Pavosevic
,
B. J. G.
Rousseau
, and
S.
Hammes-Schiffer
, “
Multicomponent orbital-optimized perturbation theory methods: Approaching coupled cluster Accuracy at lower cost
,”
J. Phys. Chem. Lett.
11
,
1578
1583
(
2020
).
7.
K. R.
Brorsen
, “
Quantifying multireference character in multicomponent systems with heat-bath configuration interaction
,”
J. Chem. Theory Comput.
16
,
2379
2388
(
2020
).
8.
M.
Tachikawa
,
K.
Mori
,
H.
Nakai
, and
K.
Iguchi
, “
An extension of ab initio molecular orbital theory to nuclear motion
,”
Chem. Phys. Lett.
290
,
437
442
(
1998
).
9.
M.
Tachikawa
, “
Multi-component molecular orbital theory for electrons and nuclei including many-body effect with full configuration interaction treatment: Isotope effects on hydrogen molecules
,”
Chem. Phys. Lett.
360
,
494
500
(
2002
).
10.
S. P.
Webb
,
T.
Iordanov
, and
S.
Hammes-Schiffer
, “
Multiconfigurational nuclear-electronic orbital approach: Incorporation of nuclear quantum effects in electronic structure calculations
,”
J. Chem. Phys.
117
,
4106
4118
(
2002
).
11.
R.
Flores-Moreno
,
E.
Posada
,
F.
Moncada
,
J.
Romero
,
J.
Charry
,
M.
Díaz-Tinoco
,
S. A.
González
,
N. F.
Aguirre
, and
A.
Reyes
, “
LOWDIN: The any particle molecular orbital code
,”
Int. J. Quantum Chem.
114
,
50
56
(
2014
).
12.
F.
Pavošević
,
T.
Culpitt
, and
S.
Hammes-Schiffer
, “
Multicomponent quantum chemistry: Integrating electornic and nuclear quantum effects via the nuclear-electronic orbital method
,”
Chem. Rev.
120
,
4222
4253
(
2020
).
13.
A.
Szabo
and
N. S.
Ostlund
,
Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory
(
Courier Corporation
,
2012
).
14.
T.
Helgaker
,
P.
Jorgensen
, and
J.
Olsen
,
Molecular Electronic-Structure Theory
(
John Wiley & Sons
,
2014
).
15.
S.
Hirata
,
X.
He
,
M. R.
Hermes
, and
S. Y.
Willow
, “
Second-order many-body perturbation theory: An eternal frontier
,”
J. Phys. Chem. A
118
,
655
672
(
2014
).
16.
H.
Nakai
and
K.
Sodeyama
, “
Many-body effects in nonadiabatic molecular theory for simultaneous determination of nuclear and electronic wave functions: Ab initio NOMO/MBPT and CC methods
,”
J. Chem. Phys.
118
,
1119
1127
(
2003
).
17.
C.
Swalina
,
M. V.
Pak
, and
S.
Hammes-Schiffer
, “
Alternative formulation of many-body perturbation theory for electron–proton correlation
,”
Chem. Phys. Lett.
404
,
394
399
(
2005
).
18.
D.
Mejía-Rodríguez
and
A.
de la Lande
, “
Multicomponent density functional theory with density fitting
,”
J. Chem. Phys.
150
,
174115
(
2019
).
19.
S. H.
Vosko
,
L.
Wilk
, and
M.
Nusair
, “
Accurate spin-dependent electron liquid correlation energies for local spin density calculations: A critical analysis
,”
Can. J. Phys.
58
,
1200
1211
(
1980
).
20.
A. D.
Becke
, “
Density-functional exchange-energy approximation with correct asymptotic behavior
,”
Phys. Rev. A
38
,
3098
(
1988
).
21.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
, “
Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density
,”
Phys. Rev. B
37
,
785
(
1988
).
22.
A. D.
Becke
, “
Density-functional thermochemistry. I. The effect of the exchange-only gradient correction
,”
J. Chem. Phys.
96
,
2155
2160
(
1992
).
23.
P. J.
Stephens
,
F. J.
Devlin
,
C. F.
Chabalowski
, and
M. J.
Frisch
, “
Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields
,”
J. Phys. Chem.
98
,
11623
11627
(
1994
).
24.
T. H.
Dunning
, Jr.
, “
Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen
,”
J. Chem. Phys.
90
,
1007
1023
(
1989
).
25.
R. A.
Kendall
,
T. H.
Dunning
, Jr.
, and
R. J.
Harrison
, “
Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions
,”
J. Chem. Phys.
96
,
6796
6806
(
1992
).
26.
E. A.
Hylleraas
, “
Über den grundterm der zweielektronenprobleme von H, He, Li+, Be++ usw
,”
Z. Phys.
65
,
209
225
(
1930
).
27.
R. C.
Lochan
and
M.
Head-Gordon
, “
Orbital-optimized opposite-spin scaled second-order correlation: An economical method to improve the description of open-shell molecules
,”
J. Chem. Phys.
126
,
164101
(
2007
).
28.
W.
Kurlancheek
and
M.
Head-Gordon
, “
Violations of N-representability from spin-unrestricted orbitals in Møller–Plesset perturbation theory and related double-hybrid density functional theory
,”
Mol. Phys.
107
,
1223
1232
(
2009
).
29.
F.
Neese
,
T.
Schwabe
,
S.
Kossmann
,
B.
Schirmer
, and
S.
Grimme
, “
Assessment of orbital-optimized, spin-component scaled second-order many-body perturbation theory for thermochemistry and kinetics
,”
J. Chem. Theory Comput.
5
,
3060
3073
(
2009
).
30.
S.
Kossmann
and
F.
Neese
, “
Efficient structure optimization with second-order many-body perturbation theory: The RIJCOSX-MP2 method
,”
J. Chem. Theory Comput.
6
,
2325
2338
(
2010
).
31.
U.
Bozkaya
,
J. M.
Turney
,
Y.
Yamaguchi
,
H. F.
Schaefer
 III
, and
C. D.
Sherrill
, “
Quadratically convergent algorithm for orbital optimization in the orbital-optimized coupled-cluster doubles method and in orbital-optimized second-order Møller-Plesset perturbation theory
,”
J. Chem. Phys.
135
,
104103
(
2011
).
32.
P.
Pulay
, “
Convergence acceleration of iterative sequences. The case of SCF iteration
,”
Chem. Phys. Lett.
73
,
393
398
(
1980
).
33.
P.
Pulay
, “
Improved SCF convergence acceleration
,”
J. Comput. Chem.
3
,
556
560
(
1982
).
34.
Q.
Sun
,
T. C.
Berkelbach
,
N. S.
Blunt
,
G. H.
Booth
,
S.
Guo
,
Z.
Li
,
J.
Liu
,
J. D.
McClain
,
E. R.
Sayfutyarova
,
S.
Sharma
,
S.
Wouters
, and
G. K. L.
Chan
, “
PySCF: The Python-based simulations of chemistry framework
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1340
(
2018
).
35.
F.
Neese
, “
The ORCA program system
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
73
78
(
2012
).
36.
F.
Neese
, “
Software update: The ORCA program system, version 4.0
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1327
(
2018
).
37.
R. D.
Bardo
and
K.
Ruedenberg
, “
Even-tempered atomic orbitals. VI. Optimal orbital exponents and optimal contractions of Gaussian primitives for hydrogen, carbon, and oxygen in molecules
,”
J. Chem. Phys.
60
,
918
931
(
1974
).
38.
C. C.
Marston
and
G. G.
Balint-Kurti
, “
The Fourier grid Hamiltonian method for bound state eigenvalues and eigenfunctions
,”
J. Chem. Phys.
91
,
3571
3576
(
1989
).
39.
S. P.
Webb
and
S.
Hammes-Schiffer
, “
Fourier grid Hamiltonian multiconfigurational self-consistent-field: A method to calculate multidimensional hydrogen vibrational wavefunctions
,”
J. Chem. Phys.
113
,
5214
5227
(
2000
).
40.
J. B.
Cumming
and
P.
Kebarle
, “
Summary of gas phase acidity measurements involving acids AH. Entropy changes in proton transfer reactions involving negative ions. Bond dissociation energies D(A—H) and electron affinities EA(A)
,”
Can. J. Chem.
56
,
1
9
(
1978
).
41.
W. L.
Jolly
,
Modern Inorganic Chemistry
(
McGraw-Hill College
,
1984
).
42.
S. T.
Graul
,
M. E.
Schnute
, and
R. R.
Squires
, “
Gas-phase acidities of carboxylic acids and alcohols from collision-induced dissociation of dimer cluster ions
,”
Int. J. Mass Spectrom. Ion Processes
96
,
181
198
(
1990
).
43.
E. P. L.
Hunter
and
S. G.
Lias
, “
Evaluated gas phase basicities and proton affinities of molecules: An update
,”
J. Phys. Chem. Ref. Data
27
,
413
656
(
1998
).