Anatomy of molecular properties evaluated with explicitly correlated electronic wave functions

ABSTRACT Static electric dipole and quadrupole moments were evaluated at the explicitly correlated second-order Møller–Plesset (MP2-F12) level for BH, CO, H2O, and HF molecules. The electron correlation contributions to the multipole moments were further decomposed into the direct (unrelaxed) and indirect (orbital response) components; we found that both components are equally important for the conventional (MP2) contribution, whereas the F12 correction to these properties originates primarily from the orbital response effects. Finally, the direct contribution dominates in the perturbative Hartree–Fock basis set incompleteness (CABS singles) correction. Two basis set families were employed: the standard aug-cc-pVXZ series and the cc-pVXZ-F12 series designed specifically for the F12 methods. The aug-cc-pVXZ MP2-F12 multipole moments usually have smaller basis set errors than the cc-pVXZ-F12 counterparts, albeit their differences are small at the triple (X = T) and quadruple (X = Q) zeta level. With the MP2-F12 calculations, the basis set errors of dipole and quadrupole moments can be reduced to ∼0.001 a.u., or roughly 0.1%, at the aug-cc-pVDZ and aug-cc-pVTZ levels, respectively.


Introduction
Rapid development of high-order many-body electronic structure models, such as coupled-cluster [1] including up to quadruple excitations, in combination with the basis set extrapolation and/or explicit correlation approaches, has recently allowed confident prediction of chemical differences of electronic energies in chemicals (1 kcal/mol) with a better accuracy [2]. These developments also allow predictive computation of other molecular properties, such as the electromagnetic response properties. For example, the best theoretical determination of the vibrational expectation value for the electric dipole moment CONTACT  of a water molecule by Lodi et al. [3] using high-order coupled-cluster methods combined with basis set extrapolation and incorporation of vibrational, relativistic, and post-Born-Oppenheimer [4] effects agrees with the best experimentally derived estimate [5,6] to 10 −3 a.u. This paper concerns with accurate computation of molecular response properties without resorting to the empirical basis set extrapolation, but rather via the use of the efficient explicitly correlated R12/F12 wave function ansatz [7]. Explicitly correlated wave functions, by including two-electron basis functions that depend explicitly on interelectronic distances, have smaller and more rapidly convergent basis set errors than the conventional wave function methods based on Slater determinants only. The most generally applicable explicitly correlated methods are the R12/F12 methods pioneered by Kutzelnigg [7,8] and developed to fruition by many others [9][10][11][12]. Although the majority of applications of F12 methods focus on accurate molecular energies, F12 methods have been also applied to analytic computation of molecular properties, such as dipole moments, polarisabilities, and first and second hyperpolarisabilities [13][14][15][16]. First, Kordel et al. [13] implemented the analytic calculation of the relaxed one-electron reduced density matrix (1-RDM) at the MP2-R12 level (with R12 indicating the use of the original linear correlation factor). Later, Klopper et al. [14] reported the analytic evaluation of 1-RDM at the MP2-F12 level of theory based on the modern F12 technology. Neiss and Hättig [15] were the first to report the analytic evaluation of response properties with an explicitly correlated coupled-cluster model, namely the CCSD(R12) method based on the older R12 technology. Most recently, Yang and Hättig extended these efforts to the CCSD(F12) method using the modern F12 formulation [16]. Within the CCSD(F12) model, Hanauer and Köhn [17] studied the response properties with a modified version of Ten-no's fixed amplitude ansatz (XSP approach [18]), which includes the explicitly correlated geminals generated from singly excited determinants.
In this paper, we report a detailed assessment of the MP2-F12 method for computing molecular one-electron properties such as the electric multipole moments. The motivation for this work is to develop efficient methods for evaluating molecular properties with explicitly correlated wave functions. Consider that the efficient use of the coupled-cluster F12 methods for computing molecular energies involves drastic (but benign) approximations [19][20][21][22] relative to the extremely complex and expensive CCSD-F12 method [23,24]; for molecular properties, one expects that by the judicious use of approximations, it will be possible to compute response properties with coupledcluster F12 wave functions in a cost-effective manner. The first step is to understand the anatomy of the prototypical molecular properties obtained with the simplest explicitly correlated method, namely MP2-F12. This approach will be extended to the approximate coupled-cluster F12 models elsewhere.
Our work is based on the state-of-the-art F12 technology using the diagonal orbital-invariant (SP) F12 ansatz [25] and CABS approximation [26]. The geminal Fock matrix (intermediate B) was evaluated using the approach of Kedžuch et al. (so-called approximation C) [27]; previous work used other formulations of this term [13][14][15][16]. We limited our study to static response properties (electric dipole and quadrupole moments) of relatively small systems (BH, HF, H 2 O, and CO) for which accurate basis set limits can be obtained by extrapolation using the standard aug-cc-pVXZ basis set family [28][29][30]. Explicitly correlated wave functions were computed with aug-cc-pVXZ series as well as the cc-pVXZ-F12 series of basis sets specifically designed for the F12 methods by Peterson et al. [31][32][33]. Although the absolute electronic energies obtained with a cc-pVXZ-F12 basis are closer to the CBS limit than the aug-cc-pVXZ counterpart, our previous work has demonstrated that neither one of these two basis set series is a clear winner for prediction of chemical energy differences [34].
This paper is organised as follows. First, in Section 2, we briefly review the MP2-F12 method, and then present the formalism for the evaluation of the MP2-F12 properties via the relaxed 1-RDM. Section 3 reports the analysis of various contributions to the MP2-F12 dipole and quadrupole moments. Finally, we summarise our findings in Section 4.

MP2-F12 one-electron reduced density matrix
The MP2-F12 correlation energy is obtained by minimising the second-order Hylleraas functional, where the zeroth-order wave function, |0 , is the Hartree-Fock Slater determinant,F N ≡ F μ νã ν μ is the corresponding normal-ordered Fock operator,Ŵ N ≡ (1/4)ḡ μν ρσã ρσ μν is the normal-ordered perturbation operator, and the firstorder wave function, |1 , is composed of standard and geminal double excitations. In this work, the diagonal fixed-coefficient ansatz (also known as SP ansatz) of Tenno [25] is adopted. By restricting ourselves to the closedshell reference states, the first-order wave function is written as where t i j ab are the amplitudes of conventional double excitations (note that these amplitudes differ from the conventional MP1 amplitudes due to the coupling to the geminal terms, see Equation (14)), andR i j αβ are the antisymmetrised integrals of the projected geminal correlation factor, f(r 12 ), with the bar denoting antisymmetrisation:R i j ProjectorQ 12 is introduced to ensure that the geminal functions are orthogonal to both the reference and standard double excitations [35]: whereÔ andV are the projectors on the occupied and virtual orbitals. The notation is explained in detail in Appendix.
In the SP approach, the geminal correlation factor is determined to satisfy the first-order (natural) singlet and triplet coalescence conditions [25]. Thus, the spinadapted geminal correlation factor can be written as [34] f (r 12 where C 0 = 1/2, C 1 = 1/4, andP 0 andP 1 are the singlet and triplet spin projectors, respectively, and γ (r 12 ) is the standard exponential correlation factor introduced by Ten-no [36]. In this work, it is approximated with an expansion in terms of Gaussian geminals [37]: where l c (∼1 a 0 ) is a correlation length scale and depends empirically on the orbital basis set. Static response of the molecular energy with respect to a perturbation is expressed as a Taylor series, with the n-th order derivatives identified as the corresponding response molecular properties [38]. The electric dipole moment (μ), for example, is obtained as the first-derivative of the energy with respect to an external electric field (F); the result can be expressed as a trace of a product of the corresponding operator matrixμ with the relaxed 1-RDM D: The appearance of the response of the Hartree-Fock orbitals to the perturbation is what distinguishes D from the regular (unrelaxed) 1-RDM of MP2-F12. The correlation contribution to the latter is obtained simply as the expectation value 1|a κ λ |1 . The detailed expressions for the unrelaxed correlation contributions to D are: where intermediate A, is related to the coupling intermediate C: The 1-RDM expressions are simplified significantly if we assume the extended Brillouin condition (EBC), F a α = 0. For energies, the effect of EBC is generally small: it is only noticeable for double-zeta basis sets, and is completely negligible for triple-zeta and larger basis sets [37].
The assumption of EBC eliminates intermediates C and A and turns t i j ab into the standard MP1 amplitudes, i.e. the EBC assumption eliminates most of the 'coupling' terms which are bilinear in conventional and F12 amplitudes. The only remaining coupling contribution that survives EBC isR αc i j t i j bc , the first term in Equation (12). The 1-RDM obtained by ignoring this term and assuming the EBC will be termed the coupling-free MP2-F12 1-RDM.
To obtain relaxed one-electron properties, the response of the Hartree-Fock orbitals to the perturbation has to be taken into account. This is achieved by including the Hartree-Fock stationarity conditions, and, if the frozen-core approximation is invoked, as constraints into the Hylleraas functional by the Lagrange method of undetermined multipliers; the multipliersκ m a and κ m i are determined by making the resulting Lagrangian, stationary with respect to variation of the non-redundant orbital rotation parameters. κ m a and κ m i entering the Lagrangian through the orthogonal orbital rotation operator is implicit in the bras/kets: The resulting stationarity conditions, ∂L (2) /∂κ = 0, are linear equations for the multipliers: these are also known as the Z-vector equations [39] and bear close similarity to the coupled-perturbed Hartree-Fock (CPHF) equations. The right-hand sides of Equations (22) and (23), Y m a and Y m i , are obtained by differentiating the MP2-F12 Hylleraas functional with respect to the orbital rotation parameters. In the canonical basis, the MP2-F12 Hylleraas functional is expressed as where the elements of intermediates V, B, and X familiar from the MP2-F12 formalism are defined as By differentiating H (2) in the above expression, we obtained two types of contributions to Y m a : from active (non-frozen) occupied orbitals only, (Y act ) i a , and from all occupied orbitals, (Y all ) m a : and one type of contribution to Y i m : wherer mn κλ represent the antisymmetrised integrals of the correlation factor, f(r 12 ), and the definition of V, B, and X intermediates here is the same as in Equations (25)- (27). D k l and D α β in (29) are the unrelaxed/direct oneelectron densities defined in Equations (10)-(13); we will denote them by D dir in the following sections. Using these expressions for Y, we solve Equations (22) and (23) to produce the multipliersκ m a andκ m i . The multipliers are the 'orbital response' contribution to the relaxed 1-RDM and will be denoted as D or . Similar to the (unrelaxed) one-electron density, the coupling-free 'orbital response' contribution is produced neglecting the terms that involve both geminal integrals and amplitudes t in Y.
The infinite-range summations in the intermediates V, C, and X were approximated with finite sums over the complementary auxiliary basis set (CABS) [26]. For the intermediate B, approximation C of Kedžuch et al. [27] and the CABS approach [26] were used. After integrating over the spin degrees of freedom, the spin-free expressions for these intermediates were obtained. For example, where the r integrals are the two-electron integrals of the bare correlation factor, γ (r 12 ). The procedure for integrating over the spin coordinates is the same as the one described in our previous paper [34], except we only considered the closed-shell reference state in this work.
As is usual (and necessary) in the F12 methods, the basis set incompleteness of the Hartree-Fock wave function was corrected by a first-order perturbative correction due to single excitations into CABS (called as 'CABS singles' [21]). The CABS singles correction contributes to the (unrelaxed) one-electron density: where t α i are CABS singles amplitudes. Following the same procedure described previously, the CABS singles contribution to the right-hand side of the orbital response equations can be obtained: with which we can solve the orbital response equations for the CABS singles correction and obtain its orbital response contribution. We implemented the F12 and CABS singles corrections to 1-RDM in the Massively Parallel Quantum Chemistry (MPQC) package [40]. The implementation was validated by comparing the electric dipole moments computed with the explicitly correlated 1-RDM against the values obtained by numerical differentiation of the MP2-F12 energy with respect to the electric field. While the agreement to seven to eight decimal places is possible for the HF, CABS singles, and MP2 contributions, for the F12 correction, the agreement is only 0.001 a.u. or better. Although this is sufficient for our target accuracy, it is important to discuss the reason for the poorer agreement between numerical and analytic results. The limited precision of the computed multipole moments is due to the use of 1-RDM to evaluate molecular properties. Consider the following contribution to the coupling-free MP2-F12 1-RDM from intermediate B: Expectation value of any operator with respect to this contribution to the density will be slowly convergent with respect to the basis set size. This can be easily seen by recognising that the expectation value of the identity operator with respect to this density is nothing but a matrix element of intermediate X: Direct evaluation of X i j i j by a double insertion of a finite CABS basis results in a slowly convergent sum and is the whole rationale of the R12/F12 method where such sums are evaluated with partial analytic integration. Similarly, for any other operatorÔ, the expectation value O β α (D F12 ) α β will be slowly convergent. For example, ifÔ = F, this expectation value is up to a constant equivalent to intermediate B.
The rigorous strategy for the evaluation of O β α (D F12 ) α β in the R12/F12 framework depends on whether opera-torÔ commutes with the correlation factor. For operators that do not commute, such as the magnetic multipole operators or relativistic corrections [41], one has to carry out the R12/F12-style manipulations similar to those involved in the evaluation of intermediate B [8,27]; this will lead to the need to compute integrals of the commutator ofÔ with the correlation factor. This is inconvenient since the property evaluation calls for evaluation of property-dependent two-electron integrals. For property operators that do commute with the correlation factor, such as the electric multipole moments, it is possible to avoid the need for such integrals and evaluate such expectation values robustly. In this work, we choose the density-based route for the sake of universality since it is appropriate for any property operator. We will consider the alternative analytic formulation of molecular properties, based on the strategies discussed above, in future work.

Computational details
The electric dipole operator is expressed as where α represents a Cartesian axis (x, y, or z), Q is the nuclear charge, r Aα and r iα are positions of nuclei and electrons, respectively. For the electric quadrupole moment, we used the traceless moment operator defined by Buckingham [42]: where β also represents a Cartesian axis. For the four molecules studied, there is only one non-zero component of the dipole moment, μ z , whereas all the diagonal components of the traceless quadrupole moment, xx , yy , and zz , are non-zero. Moreover, in the linear molecules, the diagonal components of the quadrupole moment have the following relationship: xx = yy = −1/2 zz . Hence, we are only concerned with μ z and zz for BH, CO, and HF. For the nonlinear H 2 O molecule, xx and yy are also studied in addition to μ z and zz . The molecular geometries of BH, HF, H 2 O, and CO were obtained from the literature [3,43,44]. We employed two orbital basis set series: the standard aug-cc-pVXZ family (X = D, T, Q, 5) designed by Dunning et al. [28][29][30] for computations with standard correlated methods, and the cc-pVXZ-F12 family (X = D, T, Q) designed by Peterson et al. [31] specifically for F12 methods. The matching complementary auxiliary basis sets (CABS) were constructed by using the CABS+ approach [26] from the augcc-pVXZ/OptRI [45] and cc-pVXZ-F12/OptRI basis sets [32]. In our F12 calculations, the recommended correlation length scales [31] were used, and the correlation factors were expanded in terms of six Gaussian geminals [46]. In addition, the evaluation of all post-Hartree-Fock quantities utilised the robust density fitting [47]. The aug-cc-pV(X + 1)Z − RIandcc − pV(X + 1)Z-RI basis sets were used for the density fitting in the aug-cc-pVXZ and cc-pVXZ-F12 calculations, respectively, where X is the cardinal number of the corresponding orbital basis. All calculations were performed with the developmental version (git version: 518d8ae) of MPQC (freely available under the GPL license at http://www.mpqc.org/).

Direct and orbital response contributions to dipole and quadrupole moments
We first investigated the direct (from D dir ) and the orbital response (from D or ) contributions to the dipole and quadrupole moment components, μ z and zz , within the CABS singles ((2) S ), MP2 correlation, and F12 correlation corrections. We also studied the effects of the coupling terms in the F12 correction. In these calculations, the aug-cc-pVXZ basis sets (X = D, T, Q) were used, and the corresponding results are presented in Figures 1  and 2.
We found that the CABS singles correction to the dipole moment is only significant at the double-ζ level for BH, HF, and H 2 O, whereas at the triple-and quadruple-ζ levels, it is very small. For CO, even the double-ζ value of the CABS singles correction is small (−0.00133 a.u.). Nevertheless, for all four molecules, the direct density contribution is the major component in the double-ζ CABS singles correction, whereas the orbital response contribution is much smaller and ranges from −0.002 to −0.001 a.u. On the contrary, the orbital response effects are much more important in the MP2 correlation correction. In fact, the direct and orbital response contributions are both significant in the MP2 dipole moment correction. Finally, we found that the direct density contributions from the F12 corrections are negligible at all basis set levels. Even at the double-ζ level, their absolute values are less than 0.001 a.u. Hence, the F12 contributions result mostly from the orbital response effects. The inclusion of the coupling appears to have a small but not always negligible effect on the F12 calculations. For CO, for example, the coupling contributes 0.00367 a.u. to the total F12 correction (−0.0076 a.u.).
Similar to the dipole calculations, we found that the orbital response contribution is a relatively small portion of the CABS singles correction to the quadrupole moment for the linear molecules (see Figure 2(a,b,d)). For HF and CO, we found that the double-ζ CABS singles correction appears to be similar to or larger compared to the MP2 correlation correction. This suggests that the double-ζ CABS singles correction may be more important for computing quadrupole than dipole moments. In the MP2 correlation correction to zz , the orbital response contribution increases with the size of the basis set, and it becomes much larger than the direct density contribution at the triple-and quadruple-ζ levels for BH and CO. For HF, orbital response effects at the double-ζ level are almost negligible (on the order of 0.0001 a.u.), but their contribution becomes comparable to the direct density contribution at the quadruple-ζ level. Again, for these linear molecules, most of the F12 correlation contributions come from the orbital response effects (>90% for F12 and >70% for F12 C ).
The observed pattern for zz of H 2 O is similar to that of the HF molecule (Figure 2(c)): (1) the double-ζ CABS singles correction is significant and much larger than the corresponding MP2 correlation correction (−0.01027 vs. −0.00165 a.u.), and the orbital response effects are negligible for the CABS singles correction; (2) the orbital response effects on the MP2 correlation correction are still very small at the double-ζ level (0.00012 a.u.), but become comparable to the direct density contribution at higher levels; (3) in the F12 contributions, the orbital response effects dominate. On the other hand, xx and yy for H 2 O resemble the behaviour of zz for BH. The exceptions here are their F12 correlation corrections (see the Supporting Information), which do not decrease monotonically with respect to the basis set size. Nevertheless, the F12 correction for xx and yy still results primarily from the orbital response effects.

Basis set convergence of Hartree-Fock and correlation components of the dipole moment.
To investigate the basis set convergence of various contributions to the dipole moment component, μ z , we used two basis set families: aug-cc-pVXZ and cc-pVXZ-F12. The HF/aug-cc-pV6Z values were used as the Hartree-Fock basis set limit estimate, and the basis set limits for the MP2 correlation corrections were estimated using the standard X −3 extrapolation formula [48] from the aug-cc-pVQZ and aug-cc-pV5Z values.
The basis set convergence of the HF and HF(2) S (HF with the CABS singles correction) contributions to the dipole moments of the four molecules are shown in Figure 3. For BH, HF, and CO, the CABS singles correction reduces the HF basis set errors considerably at the double-and triple-ζ levels. With the aug-cc-pVTZ basis set, the HF(2) S calculations are essentially converged. For H 2 O, even though we saw a relatively fast convergence of the HF contribution to the basis set limit, the HF basis set errors are still large at the double-ζ level (see Figure 3(c)). Yet, the basis set error is essentially zero already at the   HF(2) S /aug-cc-pVDZ level of theory. In general, the aug-cc-pVXZ family leads to better convergence to the basis set limit in the dipole calculations. The exceptions are the cc-pVDZ-F12 HF(2) S dipole moments for BH and CO, which are closer to the basis set limits than the corresponding aug-cc-pVDZ values. Figure 4 shows the basis set convergence of various correlation contributions to the z component of the dipole moment for the four molecules, where MP2-F12 and MP2-F12 C refer to the F12-corrected correlation contributions with and without the coupling between the conventional and F12 terms in the wave function (see Section 2 ). Such coupling is often neglected when computing energies, due to the often better results obtained without the coupling for small basis sets. We find that for properties, the cancellation is less reliable. As expected, the MP2 correlation contribution from the aug-cc-pVXZ calculations converges very slowly to the basis set limit, but the F12 correction significantly improves the basis set convergence. For BH, H 2 O, and CO, the MP2-F12/aug-cc-pVDZ values are already very close to the basis set limits, and are even better than the MP2/aug-cc-pV5Z values. Moreover, the coupling has a negligible effect on the aug-cc-pVXZ calculations of the BH dipole moment: the MP2-F12 C values are almost identical to the corresponding MP2-F12 values. On the other hand, for H 2 O and CO, the MP2-F12 C calculations yield larger errors than the MP2-F12 calculations at the double-ζ level, due to the fortuitous cancellation of errors in the latter. However, such cancellation is not always observed: the MP2-F12 C /aug-cc-pVXZ μ z for the HF molecule is more accurate than the corresponding MP2-F12 value. As expected, the difference between MP2-F12 and MP2-F12 C is essentially negligible with triple-ζ and larger basis sets.
We found that in all cases except the MP2-F12 C dipole moment of CO, the basis set errors of dipole moments are smaller with the aug-cc-pVDZ basis set than with the corresponding (larger) cc-pVDZ-F12 basis. The differences between the two series become smaller at the tripleand quadruple-ζ levels, but are not negligible. The better relative performance of aug-cc-pVXZ series, albeit for a small test set, is unexpected as the basis set errors of absolute correlation energies obtained with the cc-pVXZ-F12 series are usually much better than their aug-cc-pVXZ counterparts (however, we already observed that the better absolute energies do not necessarily mean better other properties, such as chemical energy differences [34]).
The difference between the two basis set families is especially pronounced for the dipole calculations on H 2 O, which contains two hydrogen atoms. At the doubleζ level, the differences between the F12-corrected correlation values with the two basis set series are around 0.01 a.u. for H 2 O, whereas the corresponding differences for the other molecules are around or less than 0.005 a.u. The cc-pVXZ-F12 basis sets for non-hydrogen atoms are considerably larger than the corresponding aug-cc-pVXZ basis sets, but for the hydrogen atom they are smaller when X = T, Q (cc-pVTZ-F12 has one less d function than aug-cc-pVTZ, and cc-pVQZ-F12 has one less d and f functions compared to aug-cc-pVQZ). Hence, it is possible that the relatively worse performance of the cc-pVXZ-F12 basis sets for dipole calculations could be attributed to the hydrogen basis functions. To test this hypothesis, we performed the cc-pVXZ-F12 calculations with the aug-cc-pVXZ basis sets for hydrogen atoms (labelled as (aug-)cc-pVXZ-F12 in Table 1). We found that the basis set errors of the calculations do become significantly smaller (see Table 1), and are very close to the aug-cc-pVXZ results, especially at the double-ζ level (see the Supporting Information for the comparison of the three). For H 2 O, the basis set errors of the correlation contributions are reduced by one order of magnitude when using the aug-cc-pVXZ basis sets for the hydrogen atoms. The results here suggest that the cc-pVXZ-F12 basis sets for hydrogen might need to be revisited.

Basis set convergence of Hartree-Fock and correlation components of the quadrupole moment.
Following the same approach, we studied the basis set convergence of the HF and correlation contributions to  zz for the four test molecules. The basis set limits were determined as before: the HF CBS was approximated by the HF/aug-cc-pV6Z values, while the CBS of MP2 correlation contribution was estimated by the standard X −3 extrapolation [48] from the aug-cc-pVQZ and aug-cc-pV5Z values. Figure 5 shows the basis set convergence of the HF and HF(2) S contributions to zz for the four molecules. For BH and H 2 O, the CABS singles correction significantly reduces the HF basis set errors in the aug-cc-pVDZ calculations, and the resultant HF(2) S values are near the basis set limits. Compared to the aug-cc-pVXZ series, the cc-pVXZ-F12 HF(2) S values have larger errors for these two molecules. However, the opposite is observed for the HF molecule, where cc-pVDZ-F12 outperforms aug-cc-pVDZ and gives a result very close to the basis set limit. Finally, for CO, both basis set families yield almost identical results after the CABS correction, and the resultant basis set errors are negligible at the double-ζ level. Despite the slightly worse performance of the double-ζ calculation on the HF molecule, the use of the aug-cc-pVXZ family in the HF(2) S calculations leads to a more consistent behaviour of convergence compared to that of the cc-pVXZ-F12 family; this is seemingly related to the erratic convergence of the uncorrected HF/cc-pVXZ-F12 values.
The basis set convergence of the MP2-and F12corrected correlation contributions to zz for the test molecules are presented in Figure 6. For BH and CO, the basis set convergence of various correlation contributions to zz resembles that of the corresponding μ z : (1) the basis set convergence of the MP2 correlation contributions is slow, but the F12 correction significantly accelerates the convergence; (2) the two basis set families generate very different results at the double-ζ level (the discrepancy is on the order of 0.01 a.u.), but the difference diminishes as we increase the basis set size. Nevertheless, the aug-cc-pVTZ basis set is needed in the MP2-F12 calculations to obtain values near the basis set limits. Similarly, the aug-cc-pVTZ basis set is also required for the MP2-F12 or MP2-F12 C calculations on HF and H 2 O to approach the basis set limits (within 0.001 a.u.). Although the cc-pVXZ-F12 series gives smaller errors than the aug-cc-pVXZ series at the double-ζ F12 level of theory on these two molecules, at the triple-and quadruple-ζ level, the cc-pVXZ-F12 values are significantly worse. For example, the cc-pVDZ-F12 values of the F12-corrected correlation contributions to zz for H 2 O are very close to the basis set limit, whereas the corresponding TZ and QZ values are ∼0.001 a.u. lower than the basis set limit. This is mostly due to the poor convergence of the underlying MP2/cc-pVXZ-F12 results. In addition, the MP2-F12 and The cc-pVXZ-F basis sets were used for the oxygen atom, while the aug-cc-pVXZ basis sets were used for the hydrogen atoms. b The value was obtained using the standard X − extrapolation with the aug-cc-pVQZ and aug-cc-pVZ values.
MP2-F12 C calculations of zz seem to give similar results for all four molecules, and the only obvious discrepancy is between the double-ζ calculations on CO using these two methods.
We repeated the comparative analyses of cc-pVXZ-F12 and aug-cc-pVXZ MP2-F12 dipole moments for zz by replacing the cc-pVXZ-F12 basis sets for H with the corresponding aug-cc-pVXZ sets. For the linear molecules (BH and HF), we found that the basis set errors of the correlation contributions become smaller as a result (0.001-0.003 a.u.) and match closely the aug-cc-pVXZ results (see the figures in the Supporting Information). For H 2 O, the convergence of the correlation corrections to xx and yy is also much improved. However, the correlation corrections to zz still converge to a value much lower than the basis set limit. Since there are three Cartesian contributions, x 2 , y 2 , and z 2 (also called second moment of charge), to the electronic parts of xx , yy , and zz (see Equation (40)), their different convergence behaviours may be due to the inconsistent convergence of those Cartesian contributions.
To clarify the issue above, we computed x 2 , y 2 , and z 2 of H 2 O with the aug-cc-pVXZ, cc-pVXZ-F12, and (aug-)cc-pVXZ-F12 basis sets. The results are presented in Table 2. We found that the convergences of the MP2 values are consistent: they become smaller as the basis set size increases, and slowly approach the basis set limits from below. On the other hand, the MP2-F12 values of x 2 and z 2 with aug-cc-pVXZ and cc-pVXZ-F12 decrease with the basis set size, and thus converge to the basis set limit from above, whereas the y 2 values approach the basis set limit from the opposite direction. Although the MP2-F12 values of each Cartesian contribution converge quickly and monotonically to the basis set limits, convergence of the traceless quadrupole moment can be erratic. This inconsistency of convergence among the x 2 , y 2 , and z 2 values is especially obvious for the results from the (aug-)cc-pVXZ-F12 calculations. Despite the fact that the MP2-F12/(aug-)cc-pVXZ-F12 values closely follow the corresponding aug-cc-pVXZ values, the former sometimes differ non-negligibly from the CBS estimate even with Table . Variation of the CABS singles and F corrections to μ z (a.u.) with the CABS basis, relative to the value obtained with the aug-cc-pVXZ and aug-cc-pVXZ-CABS OBS and CABS, respectively.
(uc)aug-cc-pVQZ cc-pVXZ-F-CABS  Variation of the CABS singles and F corrections to zz (a.u.) with the CABS basis, relative to the value obtained with the aug-cc-pVXZ and aug-cc-pVXZ-CABS OBS and CABS, respectively.
(uc)aug-cc-pVQZ cc-pVXZ-F-CABS  a F C refers to the F correction that neglects the coupling between F and conventional terms in the wave function.
the quadruple-zeta basis set. This is probably more obvious from the plots in the Supporting Information (the MP2-F12 C results are also included). In addition, we noticed that the x 2 , y 2 , and z 2 values from the MP2-F12/cc-pVXZ-F12 calculations are actually above the basis set limits at all basis set levels, even though the corresponding zz values converge to a value below the basis set limit. Since H 2 O is a nonlinear molecule, the results here indicate that the calculations of the quadrupole moments for nonlinear molecules may be more challenging; and more tests are needed.

CABS approximation and the (2) S and F12 contributions to the multipole moments.
Next, we investigated how the quality of the CABS influences the MP2-F12 components of dipole and quadrupole moments. We computed the MP2-F12 values with the aug-cc-pVXZ (X = D, T) orbital basis and two additional CABS basis sets: (1) the uncontracted augcc-pVQZ basis, and (2) cc-pVXZ-F12-CABS basis set (the matching CABS for cc-pVXZ-F12). In Table 3, we list the CABS singles and F12 contributions to μ z evaluated with these two CABS relative to the reference values obtained with the standard matching aug-cc-pVXZ-CABS. We found that the results from the calculations with the different CABSs are very similar. The variation of the CABS singles correction contribution is close to or smaller than 0.001 a.u. For the F12 corrections (F12 and F12 C ) to μ z , the effect of the choice of the CABS is even smaller: the difference among results using different CABSs is on the order of 0.0001 a.u. For the calculations of zz (Table 4), the choice of the CABS has a slightly larger influence on the CABS singles correction at the double-ζ level especially for the calculations using cc-pVDZ-F12-CABS, but the variation of the CABS singles correction with the different CABSs is again reduced to or smaller than 0.001 a.u. with the triple-ζ basis set. Similar observations can be made to the F12 corrections (F12 and F12 C ) to zz , and their variation with the different CABSs is very small (ࣘ0.0001 a.u.). Thus, we conclude that the default CABS is sufficient for the dipole moment calculations and triple-or higher ζ quadrupole calculations.

Variation of the multipole moments with the geminal exponent.
All computations in this paper, unless noted, used the recommended geminal length scales [31] for the F12 calculations. The recommended geminal exponents are obtained by minimising the average error for a test set of molecules, and thus do not guarantee optimality for any given system. For example, in our previous paper [34], we found that, as expected, the recommended values are not necessarily optimal for any given chemical energy difference. Thus, we investigated how the F12 contribution to the dipole and quadrupole moments varies with the geminal exponent. Table 5 lists the aug-cc-pVXZ (X = D, T) MP2-F12 basis set errors of the multipole moments obtained with a range of geminal exponents around the recommended values for these basis sets. As expected, the recommended geminal exponents do not necessarily yield the smallest errors. Also, as expected, unlike the absolute MP2-F12 correlation energy obtained from the Hylleraas functional, the molecular properties are not variational relative to the geminal exponent: their basis set errors can be both positive and negative. For HF and H 2 O, the variation of the multipole moments is negligible (ࣘ0.001 a.u.), but for BH and CO, variation with the geminal exponent is somewhat stronger. The use of the recommended geminal exponents seems warranted in most cases, since these values tend to minimise the errors. The exceptions are the double-ζ calculations of quadrupole moments, where the geminal exponent 0.8 seems to be optimal and the discrepancies between results using 0.8 and 1.1 as the geminal exponent are 0.005 and 0.003 a.u. for BH and a The values were computed as the difference between the MP/aug-cc-pCVXZ results with and without a frozen core.
CO, respectively. In general, we conclude that the recommended geminal exponents are a reasonable choice for the F12 calculations of dipole and quadrupole moments.

Effect of core correlation on the dipole and quadrupole moments
Until this point, all dipole and quadrupole moment computations correlated the valence electrons only. To investigate the influence of the core-core and core-valence correlation effects, we compared the results from the MP2/aug-cc-pCVXZ (X = D, T, Q) calculations with and without frozen core (Table 6). We found that the effect of core correlation on μ z lacks universal pattern. While for HF and H 2 O the core correlation correction grows with the basis set cardinal number, it becomes smaller for BH and oscillates for CO. The core correlation corrections to μ z are relatively small, around 0.001 a.u.; this is comparable to the residual basis set error of frozencore MP2-F12/aug-cc-pVDZ value of μ z . Therefore, it is reasonable to conclude that the influence of the corecore and core-valence electron correlations on the dipole moment calculations is not significant, and calculations with only valence electrons correlated are sufficient. The core correlation effect on zz is more systematic: the correction increases slightly with respect to the basis set size, but its magnitude varies for different molecules. While these effects are very small for the zz calculations on HF and H 2 O (on the order of 0.0001 a.u.), they are noticeably larger for the calculations on CO (∼0.003 a.u.). For the zz calculations on BH, however, we found that the core-core and core-valence correlation effects are much larger than those for other molecules. Hence, there is no The basis set limits were estimated as the sum of the HF/aug-cc-pVZ value and the extrapolated MP correlation contribution from the aug-cc-pVQZ and aug-cc-pVZ values. b The CABS singles correction was included in the calculations. c The corresponding value without the CABS singles correction is much smaller (. a.u.). d The value without the CABS singles correction is −. a.u. e The value without the CABS singles correction is −. a.u. The basis set limits were estimated as the sum of the HF/aug-cc-pVZ value and the extrapolated MP correlation contribution from the aug-cc-pVQZ and aug-cc-pVZ values. b The CABS singles correction was included in the calculations. c The corresponding value without the CABS singles correction is much smaller (. a.u.). d The value without the CABS singles correction is . a.u. e The corresponding value without the CABS singles correction is much smaller (. a.u.). f The value without the CABS singles correction is . a.u. consistent conclusion for the core-core and core-valence correlation effects on the calculations of zz , and for spectroscopic quality treatment, it will be necessary to correlate all electrons.

Summary and conclusions
In this work, we presented the formalism for computing static one-electron properties with the MP2-F12 method based on the modern F12 technology with the minimal number of approximations. We tested the performance of MP2-F12 in detail for accurate computation of the electric dipole and quadrupole moments for four molecules, BH, HF, H 2 O, and CO, with two basis set families, the standard aug-cc-pVXZ series as well as the F12-optimised cc-pVXZ-F12 series. We found that the contributions from the direct (unrelaxed) density and orbital response effects are very different in the CABS singles and MP2 correlation corrections to the dipole and quadrupole moments. In the CABS singles correction, we can simply neglect the orbital response effects, whereas both the direct density and orbital response contributions are important in the MP2 correlation contribution. On the other hand, orbital response effects are the dominant contribution in the F12 correlation, whereas the direct F12 density contribution is negligibly small.
We conclude that the CABS singles correction is generally necessary at the double-ζ level for the Hartree-Fock calculations of dipole and quadrupole moments. The inclusion of the F12 correction significantly reduces the basis set errors of the MP2 correlation contributions to both dipole and quadrupole moments. The F12 correction to the 1-RDM can be greatly simplified by neglecting the coupling between the conventional and F12 terms, with little to no increase of the basis set error in most cases. At the MP2-F12/aug-cc-pVDZ level of theory, we can obtain dipole moments within approximately 0.001 a.u. from the CBS limit (Table 7). For quadrupole moments, the MP2-F12/aug-cc-pVTZ calculations are needed to approach the MP2 basis set limits to the same accuracy ( Table 8). The performance of the aug-cc-pVXZ and cc-pVXZ-F12 basis sets usually differs significantly at the double-ζ level, but the difference becomes smaller at triple-and quadruple-ζ levels. In general, the F12 calculations with the aug-cc-pVXZ series give more accurate results compared to those with the cc-pVXZ-F12 family.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
The study was supported by the US National Science Foundation [award CHE-0847295], [award CHE-1362655]; and by Camille and Henry Dreyfus Foundation.