Molecular electric moments calculated by using natural orbital functional theory

The molecular electric dipole, quadrupole and octupole moments of a selected set of 21 spin-compensated molecules are determined employing the extended version of the Piris natural orbital functional 6 (PNOF6), using the triple-$\zeta$ Gaussian basis set with polarization functions developed by Sadlej, at the experimental geometries. The performance of the PNOF6 is established by carrying out a statistical analysis of the mean absolute errors with respect to the experiment. The calculated PNOF6 electric moments agree satisfactorily with the corresponding experimental data, and are in good agreement with the values obtained by accurate ab initio methods, namely, the coupled-cluster single and doubles (CCSD) and multi-reference single and double excitation configuration interaction (MRSD-CI) methods.


I. INTRODUCTION
The interpretation and understanding of intermolecular forces, particularly those relating to long-range electrostatic interactions, require knowledge of the electrostatic moments [1,2]. The electric moments are essential to provide simple ways to figure out the electric field behaviour of complex molecules. These electrical properties provide also information about the molecular symmetry since the electric moments depend on the geometry and charge distribution of the molecule. It has long recognized the role of electrostatic interactions in a wide range of biological phenomena [3]. The electrostatic energy is frequently the ruling contribution to molecular interactions in large biological systems, hence it is extremely important to describe properly the electrostatic potentials around these molecules. In order to improve the current treatment of the electrostatics for biomolecular simulations, which are traditionally modeled using a set of atom-centered point charges, the knowledge of higher multipole moments is required to include the effects of non-spherical charge distributions on intermolecular electrostatic interactions. In principle, one can experimentally find the components of the electric field at each point, but it turns into a formidable task for large molecular systems. There are several techniques to determine experimentally the dipole moments [4,5], but it is still very difficult to obtain precise experimental values of higher multipole moments such as quadrupole or octupole moments [1,6,7], independently of the experimental conditions. Theoretical calculations are therefore essential but challenging for quantum chemistry methods. The accurate calculation of these properties is highly dependent on the method employed [8], either regarding approximate density functionals [9] or methods based on wavefunctions [10,11]. Consequently, calculating the multipole moments is a way to assess any electronic structure method. The natural orbital functional (NOF) theory [12] has emerged in recent years [13,14] as an alternative method to conventional ab initio approaches and density func-tional theory (DFT). A series of functionals has been proposed by Piris and collaborators (PNOFi, i=1, 6) [15,16] using a reconstruction of the two-particle reduced-density matrix (2-RDM) in terms of the one-particle RDM (1-RDM) by ensuring necessary N-representability positivity conditions on the 2-RDM [17]. In this work, we employ the PNOF6 [16], which has proved a better treatment of both dynamic and non-dynamic electron correlations than its predecessors [18][19][20][21][22].
The aim of the present paper is to apply PNOF6, in its extended version, to the determination of molecular dipole and quadrupole moments of selected spincompensated molecules, namely, H 2 , HF, BH, HCl, H 2 O, H 2 CO, C 2 H 2 , C 2 H 4 , C 2 H 6 , C 6 H 6 , CH 3 CCH, CH 3 F, HCCF, ClF, CO, CO 2 , O 3 , N 2 , NH 3 , and PH 3 . Moreover, the octupole moment of CH 4 , a molecule without dipole and quadrupole moments is also studied. The Gaussian basis set of Sadlej [23,24], which has been specially developed to compute accurately molecular electric properties, is employed to perform all calculations.
We compare the obtained PNOF6 results with the experimental values reported in the literature [6,10,11,[25][26][27][28], as well as with the theoretically computed values of Bundgen et al. who used the multi-reference single and double excitation configuration interaction (MRSD-CI) method, and the coupled-cluster single and doubles (CCSD) values calculated by us. Recall that the CCSD values for one-electron properties differ from full-CI results in only 2% if no multiconfigurational character is observed [29], so they can be considered as benchmark calculations. To our knowledge, this is the first NOF study of higher multipole moments such as quadrupole and octupole moments.
This article is organized as follows. We start in section II with the basic concepts and notations related to PNOF6 and electric multipole moments. The section III is dedicated to present our results and those obtained by using CCSD and MRSD-CI methods. Here, we discuss the outcomes obtained for the dipole, quadrupole and octupole principal moments, in separate sections. The performance of PNOF6 is established by carrying out a statistical analysis of the mean absolute errors (MAE) with respect to the experimental marks.

A. The NOF Theory
We briefly describe here the theoretical framework of our approach. A more detailed description of PNOF6 can be found in reference [16]. We focus on the extended version of PNOF [30], which provides a more flexible description of the electron pairs in the NOF framework. Recall that PNOF6 is an orbital-pairing approach, which is reflected in the sum rule for the occupation numbers, namely, where p denotes a spatial natural orbital and n p its occupation number. This involves coupling each orbital g, below the Fermi level (F ), with N c orbitals above it (p > F ), so the orbital subspace Ω g ≡ {g, p 1 , p 2 , · · · , p Nc } . Taking into account the spin, each subspace contains an electron pair. Henceforth, we will denote PNOF6(N c ) the method we use, emphasizing the number N c of usually weakly-occupied orbitals employed in the description of each electron pair. The PNOF6(N c ) energy for a singlet state of an Nelectron molecule can be cast as The first term of the energy (2) draws the system as independent F = N/2 electron pairs described by the following NOF for two-electron systems, where H pp is the matrix element of the kinetic energy and nuclear attraction terms, whereas J pp =< pp|pp > is the Coulomb interaction between two electrons with opposite spins at the spatial orbital p. It is worth noting that the interaction energy, the last term of equations (2) and (3), is equal for electrons belonging to the same subspace Ω g or two different subspaces (Ω g = Ω f ), therefore, the intrapair and interpair electron correlations are equally balanced in PNOF6(N c ). The interaction energy E int pq is given by where J pq = pq|pq and K pq = pq|qp are the usual direct and exchange integrals, respectively. L pq = pp|qq is the exchange and time-inversion integral [31], which reduces to K pq for real orbitals. ∆ and Π are the auxiliary matrices proposed in reference [17] in order to reconstruct the 2-RDM in terms of the occupation numbers. The conservation of the total spin allows to determine the diagonal elements as ∆ pp = n 2 p and Π pp = n p [32], whereas known analytical necessary N -representability conditions provide bounds for the off-diagonal terms [33]. In the case of PNOF6(N c ), the off-diagonal terms of ∆ and Π matrices are where h p = (1 − n p ) is the hole in the spatial orbital p.
The other magnitudes are defined as It is noteworthy that the reconstruction of the 2-RDM, and therefore the functional (2), are independent of the orbital-pairing sum rules (1). These additional constraints are imposed to ensure that no fractional electron numbers appear when non-dynamic electron correlation effects become important [34][35][36]. Additionally, this allows the constraint-free minimization of the PNOF6(N c ) energy (2) with respect to the occupation numbers, which yields substantial savings of computational time.
At present, the procedure for the minimization of the energy (2) with respect to both the occupation numbers and the natural orbitals is carried out by the iterative diagonalization method developed by Piris and Ugalde [37] implemented in the DoNOF program package. The matrix element of the kinetic energy and nuclear attraction terms, as well as the electron repulsion integrals are inputs to our computational code. In the current implementation, we have used the GAMESS program [38,39] for this task.

B. Dipole, Quadrupole and Octupole Moments
The potential of the electric field at any point outside a distribution of charges is simply related to the electric multipole moments. As any distribution function, the essential features of the charge distribution can be characterized by its moments, thereby for an uncharged molecule the first (dipole), second (quadrupole) and third (octupole) electric moments are the most important terms in the multipole expansion, therefore, are usually sufficient to characterize its interaction with an external field. The components of the symmetric dipole, quadrupole, and octupole moments were defined by Buckingham [1,2] as where the Greek subscripts denote the Cartesian directions x , y and z . Note that the nuclear contribution is taken into account separately from the electronic contribution, which arises from the negative charge distribution over all the space. The formulas (8) and (9) define symmetric tensors in all subscripts. Moreover, equation (8) defines a traceless tensor for the quadrupole moment, namely, Θ xx + Θ yy + Θ zz = 0, similarly, equation (9) leads to Ω xxz + Ω yyz + Ω zzz = 0 for octupole tensor, and respective permutations between the subscripts x , y and z .

III. RESULTS AND DISCUSSION
In the following sections, we show the PNOF6(N c ) results obtained for the dipole, quadrupole, and octupole moments with respect to the center of mass. The chosen basis set is known to be an important factor in the calculation of molecular electric properties. We used the Gaussian basis set of Sadlej [23,24], which is a correlation-consistent valence triple-ζ basis set augmented with additional basis functions selected specifically for the correlated calculation of electric properties. Thus, it contains sufficient diffuse and polarization functions in order to give an accurate description of the outervalence region. It has been shown [40,41] that the Sadlej basis set has effectively the same accuracy as the aug-cc-pVTZ basis set.
Since the number N c of usually weakly-occupied orbitals is related to the description of the electron pairs, we begin studying the H 2 molecule, where there are not interpair correlation effects. This molecule has zero dipole moment, hence the calculated quadrupole moment values for different N c values are shown in table I. As expected, the best description of the electron pair is obtained when the number of usually weakly-occupied orbitals is maximum, in fact, the calculated quadrupole moment converges to the CCSD value, which is the full CI result for this molecule. In the present work, we will carry out all PNOF6(N c ) calculations by using the maximum N c value allowed by the Sadlej basis set for each molecule. In our calculations, we have set to one the occupancies of the core orbitals. Consequently, the maximum possible value of N c is given by the number of basis functions above the Fermi level, divided by the number of the considered strongly occupied orbitals. For comparison, we have included the available experimental data, and the calculated Hartree-Fock (HF) and CCSD values using the GAMESS program [38,39]. The experimental equilibrium geometries [10,11,42,44] have been used to carry out all calculations. The performance of theoretically obtained results is established by carrying out a statistical analysis of the mean absolute errors (MAE) with respect to the experimental data. Atomic units (a.u.) are used throughout.

A. Dipole Moment
In this work, the dipoles are aligned along the principal symmetry axis of the studied molecules, set on z direction. Table III A shows the independent component µ z of the dipole moments obtained at the HF, PNOF6(N c ), and CCSD levels of theory.
Overall, the inclusion of electron correlation effects through, both PNOF6(N c ) and CCSD, improves signifi-  cantly the performance of the HF method. PNOF6(N c ) and CCSD afford MAEs with respect to experimental data of 0.0309 a.u. and 0.0177 a.u., respectively. It is worth noting the agreement between PNOF6(N c ) and CCSD results, as well as with the experimental data. Note that the aug-cc-pVTZ basis set of Dunning [47] was used for the BH molecule since there is no Sadlej-pVTZ basis set available for Boron. In this case, the PNOF6(38) result is very close to the Full-CI/aug-cc-pVTZ value obtained by Halkier et al. [29], 0.5433 a.u., showing a result as good as the CCSD one.
The electronic structure and bonding situation of carbon monoxide is of special interest for modern electronic structure methods. The dipole moment of CO, extensively studied in Refs. [41,48,49], is very small (0.0481 a.u.) and ends at the carbon atom, although carbon is less electronegative than oxygen. The result shown in Table III A is representative, while HF gives the wrong direction for the CO dipole moment, PNOF6(9) corrects the sign, giving a result that is in excellent agreement with the experimental value. Remarkably, the result obtained at CCSD level is 34% away from the experimental value, so that it is necessary to include third order triplet excitations in the cluster theory in order to obtain a reasonable value, such as the one reported by Maroulis [49] at the CCSD(T) level, 0.0492 a.u.. Accordingly, the relevant electron correlation for CO is well accounted by the PNOF6(9) method.
Regarding the values obtained for HF, H 2 O, H 2 CO, HCl, NH 3 , and ClF, PNOF6(N c ) competes with Coupled Cluster, providing values that differ from experimental data in less than a 7%. In the case of HCCF and PH 3 , PNOF6(N c ) seems to lack relevant dynamic electron correlation and thereby the obtained dipole moments are not as accurate as the CCSD ones. Conversely, our values are in excellent agreement with experimental data in the case of the methyls CH 3 F and CH 3 CCH, often attached to large organic molecules, giving dipole moments with errors of 0.4% and 2% respectively, with respect to experimental values.
A special case is ozone, which is a molecule with strong multiconfigurational character. The PNOF6(7) dipole moves into the right direction from the HF value, but overestimates the effects of the electron correlation. Taking into account the good CCSD result for O 3 , which is not valid for higher electric moments, it seems that the dynamic electron correlation compensates for the lack of non-dynamical in this method, and could improve our numerical value of the dipole.

B. Quadrupole Moment
Tables III B and III B list the molecular quadrupole moments obtained at the HF, CCSD, MRSD-CI and PNOF6(N c ) levels of theory, along with the experimental values taken from Refs. [6, 10, 11, 25-28, 43, 45, 50, 51]. Inspection of these Tables shows that PNOF6(N c ) quadrupole moments agree satisfactorily with the experimental data, whereas the discrepancies are consistent with those observed using the CCSD and MRSD-CI methods in most cases.
In the case of linear molecules (H 2 , HF, BH, HCl, HCCF, ClF, CO, C 2 H 2 , CO 2 and N 2 ), NH 3 and PH 3 , belonging to the C 3v point symmetry group, the D 6h C 6 H 6 molecule, and the trigonal planar C 2 H 6 , which has D 3d symmetry, the relation Θ xx = Θ yy = − 1 2 Θ zz holds for quadrupole moment tensor, so Θ zz alone is sufficient to determine completely the quadrupole moment. Setting the main axis of symmetry in the z direction of the coordinate system, the results for these molecules are reported in Table III B. From the latter, one can observe that PNOF6(N c ) yields a MAE of 0.15 a.u., hence considering the added complexity of the quadrupole moment, the performance of PNOF6(N c ) is within a reasonable accuracy. Taking into account the experimental uncertainty, PNOF6(N c ) results agree with the experimental data for H 2 , HCl, CO, N 2 , PH 3 , ClF, CH 3 F, C 2 H 6 , and C 6 H 6 . The value obtained for H 2 reproduces the experimental one with high precision. It is also worth noting the excellent agreement with the experiment obtained for the quadrupole moment of Benzene, which is of great interest for many fields of chemistry and biology [28,53]. Indeed, the quadrupole moment of Benzene plays an important role in determining the crystal structures and molecular recognition in biological systems because it is the key to the intermolecular interactions between π-systems. For HCCF, NH 3 , C 2 H 2 , and CO 2 , the quadrupole moments fall out of the experimental error intervals, however, in the case of HCCF, C 2 H 2 , and CO 2 the mean relative percentage error is below 11%, whereas the results obtained for NH 3 is only 0.05 a.u. away from the higher limit of the experimental uncertainty. For CH 3 CCH the PNOF6(12) result deviates from the experimental value in a 13%, ergo more dynamic correlation is clearly necessary to improve this result, an effect not observed for the dipole moment of this molecule. For the Hydrogen fluoride, the HF result is the closest to the experimental value, however, the PNOF6(7) result is in outstanding agreement with the full-CI/aug-cc-pVTZ value of 1.6958 a.u. [29]. For the Boron monohydride, the experimental quadrupole moment is not available, so we use the full-CI/aug-cc-pVTZ calculation reported by Halkier et al. [29], 2.3293 a.u., in order to carry out the comparison. The agreement between PNOF6(38) and full-CI is good, according to the relative percentage error obtained below 1.7%. Table III B shows the Θ zz and Θ xx components obtained for H 2 O, H 2 CO, C 2 H 4 , and O 3 . In this work, we use the traceless quadrupole moment, hence two components are sufficient to determine completely this magnitude. On the other hand, MRSD-CI values are significantly better than CCSD calculations when many components of the quadrupole tensor are studied [10], thereby MRSD-CI is used as benchmark theoretical method in Table III B. According to the results reported in Table III B, PNOF6(N c ) performs better than the MRSD-CI method for this selected set of molecules. For H 2 O and H 2 CO, the PNOF6(N c ) values fall into the experimental error interval, which is specially broad for H 2 CO. In the case of the C 2 H 4 molecule, the longitudinal component Θ zz obtained with PNOF6(13) is near the limit of the experimental error interval, as well as the Θ xx component. Finally, we have the results obtained for O 3 , which is a stringent test for quadrupole calculations due to its twoconfigurational character [26,54]. One can observe that Ozone is well described by PNOF6 (7) comparing to the results obtained by using HF and MRSD-CI methods.

C. Octupole Moment
The octupole moment is particularly interesting in the case of methane. The octupole moment is the first nonzero term in the multipole expansion of the electrostatic interaction for methane molecule, so it is crucial in order to describe properly its interactions with external fields. Actually, the octupole-octupole interaction is the main long-range orientation dependent interaction in methane. Moreover, the complex charge distribution of methane, which has long been studied in the literature [44,55,56], is mainly dependent on its octupole moment, thus, the octupole moment is essential to characterize the charge distribution of tetrahedral molecules.
For tetrahedral molecules the octupole moment is simply given by one component, namely Ω = Ω xyz . Employing PNOF6 (14) with the Sadlej-pVTZ basis set at the experimental equilibrium geometry [44], the result obtained for CH 4 is Ω xyz = 2.1142 a.u., whereas the experimental mark reported in Ref. [7] is Ω xyz = 2.95 ± 0.17 a.u. Although the PNOF6(14) result falls out of the experimental interval error, this value is reasonable taking into account the discrepancies between experimental marks obtained by different experimental techniques [7]. Besides, comparing to theoretical calculations, the PNOF6 (14) value is very close to the result obtained by using CCSD, Ω xyz = 2.0595 a.u.. Consequently, we can conclude that PNOF6 (14) describes properly the octupole moment of methane.

IV. CONCLUSIONS
The PNOF6 method, in its extended version, has been assessed by comparing the molecular electric moments with  Our results show that PNOF6(N c ) is able to predict electric properties as accurate as high-level electronic structure methods such as CCSD or MRSD-CI, therefore the functional computes quite accurately the charge distribution of molecular systems. To our knowledge, this is the first NOF study of higher multipole moments such as quadrupole and octupole moments.
For PNOF6(N c ) dipole moments, the obtained MAE with respect to experimental data is 0.0309 a.u., being consistent with the theoretical benchmark calculations.
Remarkable is the result obtained by PNOF6(9) for Carbon monoxide, for which, HF gives a wrong direction of the dipole and CCSD overestimates it severely, whereas PNOF6(9) corrects the sign, giving a result that is in excellent agreement with the experimental mark.
The high performance of PNOF6(N c ) in computing electric quadrupole moments has been shown by most of the studied molecules, for which the computed values fall into the experimental interval error. It has been shown that the method is capable of providing the different components of the quadrupole moment tensor. The PNOF6(N c ) MAE with respect to the experiment is 0.1291 a.u., which is very close to the corresponding MAEs of 0.0902 a.u. and 0.1448 a.u. obtained by using the well-established CCSD and MRSD-CI methods, respectively. In particular, the results obtained for the ozone molecule with a marked multiconfigurational character, show that the method is able to treat properly non-dynamic and dynamic electron correlations. Finally, the study of the octupole moment was focused here on methane, due to its important role in the description of the long-range electrostatic interactions for this molecule. The PNOF6 (14) result is in excellent agreement with the value provided by the CCSD method.