Boosting the OEP2-sc method with spin-component scaling

The semi-canonical second-order optimised effective potential (OEP2-sc) method [R.J. Bartlett, I. Grabowski, S. Hirata and S. Ivanov, J. Chem. Phys. 122, 034104 (2005). doi:10.1063/1.1809605] is one of the most stable and efficient correlated OEP approaches. In this work, we introduce its scaled-opposite-spin (SOS) variant, named OEP2-SOS-sc. This new method further improves over the OEP2-sc for various properties, including correlation potentials, relaxed densities, and energies, fully retaining the good stability of the original approach. Thus, the OEP2-SOS-sc is an efficient and promising tool for high-level density functional theory simulations. GRAPHICAL ABSTRACT

Another approach that has been used to improve the OEP-GL2 performance is the spin-resolved (SR) version of the OEP-GL2 method [47,58,59]. In such approach, the correlation energy functional and, correspondingly, the correlation potential are decomposed into samespin (SS) and opposite-spin (OS) contributions as it is done in MP2-based approaches [60][61][62][63][64]. Those terms are then scaled by appropriately selected coefficients to reduce the large overestimation of the OEP-GL2 functional and improve the performance for different properties [47,58,59,65]. Moreover, the observation that there exists a well-defined proportionality between the OS and SS contributions allowed a further simplification of the method, whereby only the OS contribution is considered together with an appropriate c OS scaling coefficient. In this way, the OEP2-scaled-opposite-spin (SOS) method [47,58] is obtained. This has proved to have good accuracy and a superior efficiency, thanks to the higher efficiency in the computation of the OS part to the SS one. Despite its success, however, the OEP2-SOS method still suffers some limitations, because of the poor quality of the OEP-GL2 that is used as a starting reference. Nevertheless, one interesting feature of the SOS approach is that it can be applied to any MBPT formula where the spin resolution is available.
In this work, we consider applying the SOS method to the ab initio OEP2-sc energy expression, which is a better starting point than the OEP-GL2 to obtain the OEP-SOSsc approach. This new method can improve, with no additional computational effort, both the OEP2-sc and the OEP2-SOS, taking advantage of the good features of both methods. As a result, we obtain a new workhorse computational tool in the framework of ab initio DFT.
The article is organised as follows: in Section 2, we present the theory of the OEP-SOS-sc method; computational details are reported in Section 3; in Section 4, we present results for the correlation potential, electronic density, IP and correlation energy of atoms and molecules. Finally, conclusion are drawn in Section 5.

Theory
For a given orbital-and eigenvalue-dependent XC energy functional (E xc ), the OEP equation for the XC potential reads [1,3,4,6,20,21,30] which is an integral equation (Fredholm of the first kind) with the inhomogeneity given by where the index σ labels the spin (throughout we will denote the spin indexes with σ , τ , . . . ), is the static KS linear response function, and φ pσ and ε pσ denote the KS orbitals and eigenvalues, respectively, determined by the standard KS equations with v ext (r) and v J (r) the external (nuclear) and the Coulomb potentials, respectively. In all equations, we use the convention that i, j, k label occupied KS orbitals, a, b, c label virtual ones, while the indexes p, q, r, s are used otherwise.
As usual, the XC energy expression is separated into exchange and correlation contributions E xc = E x + E c . The exchange energy functional has the form of the usual Hartree-Fock (HF) exchange energy with (p σ q σ |r σ s σ ) being two-electron integrals in the Mulliken notation computed from KS orbitals. In the OEP2-sc method, the KS orbitals undergo a semicanonical transformation that makes the second-order energy expression invariant with respect to orbital rotations (mixing of occupied or virtual orbitals among themselves). The corresponding correlation energy functional is therefore with f σ pq = ε pσ δ pq − p σ |K + v xc |q σ being the Fock matrix elements defined in terms of the KS spin orbitals. For more details on orbital-dependent equations corresponding to the OEP2-sc method, we refer the reader to Refs. [8,44,46].
In the SR OEP2-sc variant, the simplest energy expression which can be considered is the OS part of the OEP2-sc correlation functional given by Equation (6), i.e.
where c OS is a single scaling factor that will be discussed later. We will name this method as OEP2-SOS-sc in order to underline that the KS orbitals undergo the same transformation as in the OEP2-sc method. This specific choice for the correlation energy functional is certainly favourable from the computational point of view. Additionally, the proper choice of the OS parameter in the energy expression can compensate for the lack of SS energy contributions (second term in Equation (6)) and of the singly excited energy terms (the last term in Equation (6)) [47,58,59,65]. This approximation is naturally justified by the proportionality between SS and OS parts of the MP2 energy expression [61] and, in addition, in our KS-DFT case this is transferred to proportionality of other quantities like electron densities [47,58,65], correlation potentials [47,58], IPs [59,66,67], etc. A similar proportionality is also observed in the case of SR OEP2sc correlation energy expression, correlation potentials, and densities which are reported in Figures S1-S3 in the supplementary file for few representative systems. This is essential in our orbital-dependent KS-DFT method, because these quantities depend directly on each other in the KS OEP scheme.
To make the method feasible, one needs to estimate the value of the c OS parameter used in the OEP2-SOS-sc method. One of the obvious choices is to fix the OEP2-SOS-sc coefficient by the direct utilisation in this context of the ansatz proposed in Ref. [47], namely where E c OS, (2) @HF and E c OS, (2) @OEPx(sc) denote the OS part of the second-order correlation energy expression computed with HF and semi-canonical transformed OEPx orbitals, respectively. However, in the present case, in contrast to what happens for the OEP2-SOS method, the E c OS, (2) @OEPx(sc) term is always almost identical in value to the E c OS, (2) @HF term, because of the use of semicanocal orbitals in the OEPx expression. This leads in all systems to a ratio E c OS, (2) @HF/E c OS, (2) @OEPx(sc) ≈ 1, as shown numerically in Figure S4 in the supplementary file (see Supporting Information). For this reason, for OEP2-SOS-sc we neglect Equation (8) and we adopt in Equation (7) simply c OS = 1.3, which is the same coefficient used in the SOS-MP2 [61] approach. Note that this choice is not only slightly simplifying the whole method but makes the scaling parameter completely system independent, thus resolving some size-dependency problems of other OEP2-SOS approaches. The choice c OS = 1.3 is further analyzed in Section 4.5 where we discuss the impact of variations of the c OS parameter on the quality of predictions.

Computational details
All calculations were performed with a locally modified ACESII [68] software package. A uncontracted cc-pVTZ basis set [69] was employed for all systems except the Be and Ne atoms for which a fully uncontracted ROOS-ATZP basis set [70] was used. We remark that the choice of basis set was mostly dictated by the need to ensure the best possible expansion for both the wave functions and the OEP potential in the linear combination of atomic orbitals (LCAO) OEP procedure. In all calculations, tight convergence criteria (SCF: 10 −8 ) was imposed to guarantee stable results. Moreover, to solve the algebraic OEP equations, the truncated singular-value decomposition (TSVD) of the response matrix was employed. The cutoff criteria in the TSVD procedure was set to 10 −6 . For technical details on this type of calculations, we refer the reader to Refs. [44,58].
To assess the quality of the new method, the results have been benchmarked against CCSD(T) [71] data obtained in the same basis set. Moreover, a comparison with CCSD, MP2, SOS-MP2 [61] and the OEP2-sc methods has been considered. In the assessment, we considered several properties: • Correlation potentials and densities as in our previous studies [44,45,58,72]. We investigated the quality of correlation potentials and densities [44,73,74] looking at their spatial behaviour. Both quantities were obtained from fully self-consistent OEP calculations. The densities were analyzed in terms of correlated densities defined as where ρ X is the density obtained from the exact exchange only (X = OEPx) [21] or HF (X = HF) calculations, for DFT and wave-function theory (WFT) methods, respectively. To provide a quantitative measure of the quality of the computed densities, we additionally computed the integrated density differences ( IDDs) [58,65,75] defined as • Total and correlation energies were calculated for the 36 systems (9 atoms and 27 molecules) listed in Table S1 in the Supporting Information. For molecular systems, geometries were taken from Ref. [76]. • Vertical ionisation potentials (VIPs) were computed as in Refs. [48,59,77,78] for 32 molecular systems listed in Table S1 in Ref. [48] using two different approaches: (i) as the energy difference between the neutral and the ionic species [66] (VIP = E(N) − E(N − 1)) and (ii) as the opposite of the energy of the highest occupied molecular orbital (HOMO) (VIP = −ε HOMO ). We note that for the exact XC functional, these two approaches should yield the same value of VIP for any system. Thus, the inspection of both is essential to assess the quality of the XC functional and potential. • Atomisation and reaction energies using the same setup as in Ref. [48]. The atomisation energies have been calculated for 27 molecular systems listed in Table S3 in the Supporting Information, whereas the open-shell and close-shell reactions are listed in Tables S4 and S5, respectively.

Correlation potential and densities
We start our analysis by inspecting correlation potential and correlated densities (Equation (9)). These are reported in Figures 1 and 2. One can note that the OEP2-SOS-sc method, in general, improves the OEP2-sc results. This is especially visible for Ne atom, where the OEP2-SOS-sc correlation potential is almost identical to the one obtained from KS[CCSD(T)] method [44,79]. For Be and the CO molecule, the improvement is less pronounced, being in some regions identical with OEP2sc results. Similar behaviour is observed for correlated densities. Also, for this property, the OEP2-SOS-sc outperforms the OEP2-sc results, being closer to CCSD(T) reference. Comparing current data with the one reported in Ref. [58] for the OEP2-SOSa method, which we recall was constructed using the GL2 correlation energy expression, shows that the OEP2-SOS-sc results are less overestimated in this context. This is properly because the OEP2-sc method provides a better starting point for constructing SR second-order functional.
To show this more quantitatively, in Table 1 we report the IDD, defined in Equation (10), corresponding to both OEP2 methods. The reported IDDs confirm that the best densities are found in the case of OEP2-SOS-sc, which improves over OEP2-sc results of about 40%. The OEP2-SOS-sc method gives mean absolute error (MAE) of 2.0 (MAE/N e = 0.19) whereas the OEP2-sc gives on average slightly worse performance (MAE = 3.25; MAE/N e = 0.31). The only exception is the N 2 dimer, where the opposite trend in densities is observed. This is a peculiar case, since the OEP2-SOS-sc correlation potential is slightly closer to the reference one than the OEP2-sc potential but nevertheless the OEP2-sc correlated density is a bit better than the OEP2-SOS-sc one (see S5 in the Supporting Information). This may depend on fine details of the correlation potential that will be investigated in the future. Thus, we see that in almost all cases the one-parameter OEP2-SOS-sc method inherits all good features of the OEP2-sc method (stability and convergence), giving in the same time, improved results.

Correlation energies
In Figure 3, we report the mean absolute relative errors (MAREs) obtained for correlation energies computed with respect to CCSD(T) data. The OEP2-SOS-sc method improves almost twice with respect to the OEP2sc method results (MARE = 7.9%, MAE = 14.44 mHa), giving for correlation energies MARE = 4.46% (MAE = 11.22 mHa), which is much closer to the CCSD MARE of 2.71% (MAE = 9.78 mHa) than for other reported methods. Almost the same improvement is observed for SOS-MP2 (MARE = 6.30%, MAE = 19.49 mHa) method with respect to MP2 (MARE = 10.42%, MAE = 24.18 mHa).
We observe an almost identical qualitative behaviour for total energies (see Table S2 in the Supporting Information). The best results in terms of MARE are obtained for the OEP2-SOS-sc method (MARE = 0.029%, In this context, the performance of MP2 and OEP2-sc methods are also quite similar. All second-order-based methods give total energy errors distinctly larger than the CCSD method, which yields the total energies MARE of 0.008% (MAE = 8.71 mHa). This might indicate that the former suffers from a slight relaxation error issue coming from the density errors discussed above.

Ionisation potentials
The MARE of the VIPs computed both as energy differences and as opposite of the HOMO energy is reported in Figure 4. One can immediately note a rather large improvement of the quality of predictions of OEP2-SOSsc with respect to its unscaled counterpart for VIPs  The VIPs obtained from the energy difference method exhibit much smaller errors [48]. This is due to two facts. First, the VIP energy of the HOMO is quite sensitive to the quality of XC potential [48]. Second, the energy difference approach fully considers the relaxation correction related to lack of electron in ionised species, which is only partially included in HOMO energy [66,80]. The OEP2-SOS-sc method outperforms here also the ionization potential equation of motion CCSD (IP-EOM-CCSD) (MARE = 1.07%, MAE = 0.13 eV) results giving a MARE of 0.63% (MAE = 0.08 eV). These are very  Table S1 in Ref. [48]. The MP2 and OEP2-sc results are taken from Ref. [48].
interesting results since the cost of a single calculation with the IP-EOM-CCSD method is O(N 6 ), whereas the cost of OEP2-SOS-sc can be reduced to O(N 3 ) [61]. The OEP2-sc and MP2 methods give here larger MAREs of 1.27% and 2.14%, respectively.
Finally, in Figure 5, we report the relative difference between VIP obtained from both approaches for OEP2-sc and OEP2-SOS-sc methods. One can note that the utilisation of the latter method gives in this case smaller differences that indicate the improved quality of the XC potentials with respect to the OEP2-sc method.

Atomisation and reaction energies
In Figure 6, we report MAEs and MAREs on atomisation, close-and open-shell reaction (OSR) energies calculated with respect to CCSD(T) data. In the case of atomisation energies, both the MAE and the MARE indicate that OEP2-SOS-sc (MAE = 3.92 kcal/mol; MARE = 3.14%) outperforms the unscaled OEP2-sc variant (MAE = 8.07 kcal/mol; MARE = 7.35%). A similar but smaller improvement is observed for the SOS-MP2 method with respect to MP2 with the former outperforming the later by about 2.1 kcal/mol. Because here the improvement is only ascribable to the better description of correlation energies due to the scaling procedure (see Section 4.2), we can infer that in the case of OEP2-SOS-sc the improvement with respect to OEP2-sc traces instead back to both a better description of the energy and a better electron density (relaxation effect on total energy). In any case, it is worth noting that both SOS-MP2 and OEP2-SOS-sc methods give better predictions than CCSD which yields here MAE of 4.08 kcal/mol (MARE = 3.15%), being at the same time computationally less demanding.  Concerning the reaction energies, the improvement of the OEP2-SOS-sc with respect to OEP2-sc is less evident: for OSRs it performs slightly better (MAEs are 8.4 kcal/mol and 7.1 kcal/mol for OEP2-sc and OEP2-SOS-sc, respectively) while for closed-shell reactions (CSRs) the performance is substantially identical (MAEs are 3.0 kcal/mol for both methods).

Impact of the c OS value on the results
As discussed so far, the results indicate that the utilisation of c OS = 1.3 scaling parameter value (the same as in the SOS-MP2 method) leads to an improved, general behaviour of the OEP2-SOS-sc method with respect to the original OEP2-sc. However, many studies have shown that spin-component scaling parameters are system [58,62,63,81] and property [64][65][66][67]   The results of this study are summarised in Figure 7. As one can note in the case of correlation energies, the value c OS = 1.3 seems to be the optimal choice. Only a slightly different situation is observed for atomisation energies. In fact, in this case, the MARE reaches its minimum around c OS ≈ 1.2 (MARE = 2.84%). For c OS = 1.3, we obtain MARE = 3.14%, which is in line with CCSD results (MARE = 3.15%). Therefore for both the energy-related properties, the choice c OS = 1.3 provides very accurate results.
For IPs (calculated from HOMO energies as IP = −ε H ), on the contrary, the minimum MARE is obtained for c OS ≈ 0.9. Thus, for best performance, we need a scaling value that is rather different from our choice. We have observed similar behaviour in the OEP2-SOS approach (OEP2-SOS(opt) method), where the scaling parameter needs to be optimised for IP only [59]. Nonetheless, for all investigated cases, the OEP2-sc (indicated by dashed horizontal line) performance is always worse than for our OEP2-SOS-sc method.
All these results indicate that, in line with all other SR methods, there exists no single, universal parametrization of the OEP2-SOS-sc method yielding the best possible results for all quantities. On the contrary, the optimal parametrization is property-dependent and systemdependent. Nevertheless, the OEP2-SOS-sc method with c OS = 1.3 is a simple and efficient choice that gives an improved performance with respect to the initial OEP2-sc variant (indicated in Figure 7 by dashed horizontal line) and other OEP2 methods based on SR second-order energy expression such as the OEP2-SOSa method.

Conclusions
We have proposed a SOS version of the OEP2-sc method (OEP2-SOS-sc) as an efficient approach to supplement OEPx with an ab initio correlation functional. The performance of the OEP2-SOS-sc method was validated for benchmark atomic and molecular systems, closed-and open-shell ones. We have shown that the OEP2-SOS-sc results outperform OEP2-sc for the calculation of correlation potentials, correlated densities, IPs and atomisation energies. In addition, OEP-SOS-sc is computationally less costly than OEP-sc as no SS and singly excited terms in the correlation functional and potential need to be computed. By fixing a single c OS parameter, we were also able to simplify the SR OEP method, making the scaling parameter completely system independent. This also resolves some size-dependency problems of the other OEP2-SOS approaches.
The new OEP2-SOS-sc method complements the family of ab initio DFT approaches, and thanks to favourable computational cost concerning similar correlated OEP approaches, good stability and accuracy, should allow for applications to the larger and more complex systems and is an alternative to standard DFT and MP2 type calculations.

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

Note about Prof. Wolniewicz
We dedicate our research to Professor Lutosław Wolniewicz, our friend, colleague, teacher and mentor. But also, and perhaps most of all, a world-class and outstanding scientist who was and will be an authority for many generations of physicists and computational and quantum chemists. Professor Wolniewicz was a theoretician physicist, incredibly talented, with great possibilities and broad knowledge, and most importantly, he was a person who was able to use his talent. His research and discoveries were of fundamental importance for developing quantum physics, chemistry and astrophysics. The super-accurate calculations of the dissociation energy of the H 2 particles carried out by Kołos and Wolniewicz [82,83] showed a contradiction between the theoretical value and the experimental value known at that time and led to the verification of the experiment. Because they were one of the first theoretical calculations carried out based on quantum mechanics that correct the experiment, they bring the collaboration between theoreticians and experimentalists to a new level.