{'fraction_non_alphanumeric': 0.044333947633784204, 'fraction_numerical': 0.014047776348273583, 'mean_word_length': 4.796049979846836, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simple non-empirical procedure for spin-component-scaled MP2 methods applied to the calculation of dissociation energy curve of noncovalently-interacting systems

We present a simple and non-empirical method to determine optimal scaling coefficients, within the (spin-component)-scaled MP2 approach, for calculating intermolecular potential energies of noncovalently-interacting systems. The method is based on an observed proportionality between (spin-component) MP2 and CCSD(T) energies for a wide range of intermolecular distances and allows to compute with high accuracy a large portion of the dissociation curve at the cost of a single CCSD(T) calculation. The accuracy of the present procedure is assessed for a series of noncovalently-interacting test systems: the obtained results reproduce CCSD(T) quality in all cases and definitely outperform conventional MP2, CCSD and SCS-MP2 results. The difficult case of the Beryllium dimer is also considered.


I. INTRODUCTION
Noncovalent interactions of organic molecules play a fundamental role in biochemistry, solvation, surface science, and supermolecular chemistry. An accurate description of the noncovalent interactions is thus very important for many applications [1]. Actually, the "golden standard" for the simulation of noncovalently interacting complexes is the coupled cluster single and double with perturbative triple (CCSD(T)) approach [2]. However, the CCSD(T) method has a high computational cost, scaling as O(N 7 ), and it is not easily applicable in general, especially when numerous single-point calculations are required as in the determination of a potential energy surface (PES) or for geometry optimizations.
In various applications is thus necessary to recover to lower-level computational methods. However, most approaches (including density functional theory (DFT) and standard Møller-Plesset second-order perturbation theory (MP2), are not fully satisfactory for an accurate description of noncovalent interactions and the corresponding PESs [3][4][5][6][7]. Therefore, in the last years different computational schemes, mainly based on variants of MP2, have been proposed to treat noncovalent interactions with sufficient accuracy and relative small computational effort [3,[6][7][8].
One of such methods is the so called spin-componentscaled MP2 (SCS-MP2) method [9,10] which is based on the spin resolved MP2 formula for the correlation energy with where i, j denote occupied orbitals, a, b denote virtual orbitals, (ia|jb) denotes a two-electron integral in the Mulliken notation, ǫ i is the energy of the i-th orbital, and c OS and c SS are the scaling coefficients for the opposite spin (OS) and same spin (SS) correlation, respectively. This approach was shown to yield quite improved results with respect to the standard MP2, when a proper choice of the scaling parameters is performed [7,11,12]. Alternatively, a good performance was also obtained by limiting the method to consider only the opposite spin part (i.e. setting c SS = 0 and optimizing only c OS ), resulting in the scaled-opposite-spin MP2 (SOS-MP2) method [13], which allows, when properly implemented, to reduce the computational scaling to O(N 4 ). However, the choice of the optimal scaling parameters in SCS-and SOS-MP2 is not trivial and different proposals have been made [10], either based on theoretical considerations or on empirical fittings, showing that the "best" scaling parameters are somehow system and basis-set dependent. This issue may be not particularly relevant for covalent bonds, where the considered binding energies and improvements over MP2 are much larger that the inaccuracies due to a non-optimal choice of the parameters. It becomes anyway relevant for noncovalently bonded complexes where the energies to be computed are much smaller. For these latter cases, very good SCS parameters were proposed by Distasio and Head-Gordon (SCS(MI)-MP2 method) [11] by optimizing the scaling coefficients against benchmark noncovalent interaction energy data. The same authors, as well as Grant Hill and Platts in a separate paper, also showed that accurate results can even be achieved by just same-spin scaling (SSS-MP2; c OS = 0, c SS ≈1.75) [11,12] methods. Additional investigations concerned the optimization of the SOS-MP2 for long-range interaction by the introduction in the SOS-MP2 method of a distance-dependent scaling coefficient and a modified distance-dependent two-electron operator [14].
The effort spent to optimize the scaling parameters and the fact that the results of these works are not in an agreement with "standard" SCS-MP2 (c OS = 6/5, c SS = 1/3) or SOS-MP2 (c OS = 1.3) values nor with the values suggested by theoretical arguments within wave function theory [10,15,16], indicates the difficulty to fix optimal values of the scaling for the proper description of noncovalent interactions. Moreover, all these approaches pay the prize of introducing a high level of empiricism. Therefore, a simple non-empirical procedure to fix the value of the various scaling parameters appears highly desirable.
In this paper we address this issue and we show a simple non-empirical scheme to fix the scaling factor in calculations of the dissociation curve of noncovalent bonded complexes. Here, non-empirical denotes the fact that scaling factors will be fixed by direct use of information from high-level ab initio calculations and not from some empirical fitting procedure. Of course, some empiricism is still implied in the use of spin-resolved MP2 formulas. Our approach focuses on scaled MP2 calculations with a single parameter. To this end, we introduce the SOS(R)-MP2, SSS(R)-MP2, and scaled MP2 (S(R)-MP2, i.e. with c SS = c OS ) methods. In acronym (R) stands for "calculated from Reference non-empirical values", to avoid confusion and distinguish from empirically scaled versions of standard SCS-MP2 methods. Our procedure is based on the observation that there exist a well defined proportionality between scaled-(same/opposite)-MP2 energies and the correlation energy computed by high-level methods (e.g. CCSD(T)), which is almost independent on the intermolecular distance. Thus, the scaling parameter can be fixed once by using the information from only one expensive highlevel calculations and successively the whole PES can be computed with high accuracy by performing only relatively cheap scaled-(same/opposite)-MP2 calculations. In particular, in this way the efficiency of the SOS-MP2 (SOS(R)-MP2) method (O(N 4 ) scaling [13]) can be fully exploited for large scale explorations of dissociation potential energy surfaces (PES) without introducing empirical parameters and achieving almost the CCSD(T) accuracy.
We acknowledge that a similar approach was already recently used in Ref. [17], to compute accurate interaction energies of an Ar monolayer adsorbed on an MgO substrate. However, in Ref. [17] only local MP2 calculations were performed while spin-resolved MP2 calculations, which are the main topic of the present work, were not considered. In addition, in the present work, a systematic assessment of the non-empirical procedure for fixing the scaling parameters is carried out for different systems and interactions.

II. METHOD
The main quantity of interest in the present paper is the inter-fragment correlation (IFC) energy, defined as: where in our calculations X=CCSD(T), CCSD, MP2, SOS(R)-MP2, SSS(R)-MP2 or S(R)-MP2. The IFC is the difference between the correlation energy (E c,AB ) of a complex (AB) and that of its constituting fragments (A and B), computed with a method X. The IFC for two exemplary cases (H 2 O dimer and Ne 2 ) is reported in Fig. 1 as a function of the intermolecular distance (R) for different methods. The IFC is always negative and its absolute value rapidly decreases with R. It can be observed that for noncovalent interactions, with good approximation, the IFC energies computed at the MP2 (or equivalently MP2(SS) or MP2(OS)) level and at the CCSD(T) level are proportional to each other over a wide range of inter-fragment distances R. This proportionality can be verified by direct comparison of the IFC energies calculated at different levels of theory for various noncovalent interacting complexes (see Fig.  2). For MP2 the proportionality is very well satisfied for all systems and distances. MP2(SS) and MP2(OS) show instead some larger deviations at higher distances. However such deviations are not very relevant, because the IFC energy is rapidly approaching zero for large distances. In fact, as shown in Fig. 1, for R/R 0 1.5 the IFC is already close or below subchemical accuracy.
The proportionality of the (spin-resolved) MP2 and coupled cluster energies can be rationalized considering that in noncovalent complexes the variation of IFC with distance is mainly due to relaxation of intramolecular This proportionality can be used to define singleparameter scaled (spin-resolved) MP2 approaches suitable for the description of the dissociation of noncovalent complexes. To do that, we define with whereR is some reference distance. With this definition we fix (at the cost of a single expensive CCSD(T) calculation) the proportionality between (spin-resolved) MP2 and CCSD(T) results, by simply imposing that at some point, the two IFC energies, or analogously, the respective total binding energies (E b ) are the same.
The choice of the reference pointR in Eqs. (8)-(10) is not a major concern, since the value of c is almost independent on it (see Fig. 2 and subsection III A). In this work we have chosen for the the reference point the equilibrium distance R 0 , computed at the CCSD(T) level. This choice is not mandatory, nevertheless it appears the best compromise between the need to avoid too short distances (where the proportionality relation may be less accurate) and the necessity to avoid the use of too small energies (as would result for large R values) in the ratio of Eqs. (8)-(10) to minimize numerical noise(see subsection III A for a further discussion).
Each of the equations (8), (9), (10) provides a simple non-empirical procedure for fixing the (systemdependent) scaling factor in one-parameter (SCS-)MP2 calculations when dissociation energies of noncovalent systems are of interest. Thus, a whole dissociation curve, comprising many single-point results, can be simulated by performing only a single CCSD(T) calculation and without the introduction of any empirical parameter.
The main question is if the accuracy of the such procedure is high enough. In section III we will show that deviations from the CCSD(T) reference results are indeed well below subchemical accuracy (i.e. 0.1 kcal/mol).

A. Computational details
To test our scaling approach we considered a representative set of small noncovalently interacting systems: He 2 , Ne 2 , He-Ne, Ar 2 (dispersion interaction), (H 2 O) 2 , (HF) 2 (hydrogen bond), (H 2 S) 2 , (HCl) 2 (dipole-dipole interaction). Additionally, benzene-HCN, the stacked benzene dimer, and Be 2 have been also considered, as particular cases. The former two are in fact relatively large systems with the first one also displaying a mixed electrostatic-dispersion character, which changes with the bond distance [18]. The latter is a system where MP2 and CCSD completely fail even qualitatively [19]. The sizes of the molecules included in our test set were limited by the need to compute in all cases the full CCSD(T) dissociation curves for reference purposes, except for benzene-HCN and the stacked benzene dimer in which case the reference data were taken from Ref. [18].
To assess the performance of the different methods to reproduce the full dissociation curves we considered the mean absolute error: where X=CCSD, MP2, SOS(R)-MP2, ..., and Note that because all the methods considered here use the Hartree-Fock exchange, the difference between the IFC energies of different methods in Eq. (12) are exactly the same as the differences in the total interaction energies. The MAE(X) defined in Eq. (11) measures the average deviation of the results from the reference CCSD(T) data over a given interval. In our work we decided to remove from this evaluation the smallest values of R, because our approach might be not fully justified at such small inter-fragment distances. Hence, the lower bound for the integration was fixed at R 0 . The upper bound of the integral was fixed instead to the value of R beyond which the CCSD(T) IFC energy is lower than 10 −4 Hartree, in order to keep the normalization factor finite in Eq. (11). Finally, we note that in all calculations no correction for the basis set superposition error (BSSE) was included, for computational simplicity and because it is readsorbed in the scaling coefficient. For the water dimer, e.g., we check that, including a BSSE correction, the MAE of Eq. (11) is 36.7, 36.6 , 5.5 10 −3 kcal/mol, for SOS(R)-MP2, SSS(R)-MP2 and S(R)-MP2 respectively. These values are in good agreement with the ones (without BSSE) reported in Tab. II.

III. RESULTS
The values of the c S , c OS , and c SS coefficients, obtained by applying Eqs. (8), (9), and (10) to the systems considered in this work, are reported in Table I, together with the reference inter-fragment separation R 0 . Inspection of the table confirms that the optimal scaling coefficients for (spin-resolved) MP2 calculations are significantly system dependent. This dependence is also more pronounced for the SSS(R)-MP2 method, in which MP2(SS) includes only a small part of the total correlation contribution, while it is much weaker for MP2, which can describe better the whole correlation effects. It is interesting to note, in addition, that the variation of the scaling coefficients among different systems is not the same for the various methods. In fact, for example the c OS is maximum for HF-HF and almost twice as big as the c OS for Ar 2 , while the value of c SS for HF-HF is about 20% smaller than the Ar 2 one and it is maximum for He-Ne. Therefore, no clear trend can be identified for the values of the scaling coefficients when different systems (and methods) are considered.
The results of Tab. I also indicate that, due to the variability of the optimal scaling coefficients, global scaling factors can hardly be expected to achieve high accuracy for a broad range of applications. In fact, even if most of the coefficients displayed in the table agree reasonably well with the SOS(MI)-MP2 (c OS(MI) ≈ 1.8) and SSS(MI)-MP2 (c SS(MI) ≈ 1.75) values [11,12], which were especially optimized for intermolecular interaction energies, some remarkable differences appear. These reflect the peculiarities of the correlation in some systems (e.g. same-spin correlation in He 2 or opposite-spin correlation in HF-HF) which cannot be captured by "average" scaling procedures. Note finally, that all the c OS values are rather different from the conventional scaling coefficient of 1.3 proposed for quantum chemical applications [13].
To demonstrate the performance of the method we report in Fig. 3 12)) at various distances. By definition all our approaches coincide with CCSD(T) at the equilibrium (reference) distance R 0 . Remarkably, however, the accuracy of the method is found to be very good at any distance, in particular for R > R 0 , where deviations from CCSD(T) values are vanishing small (of the order of 10 −3 kcal/mol). Such accuracy is not only well below subchemical accuracy but also largely sufficient to yield a highly accurate description of these dispersion dimers, as demonstrated by the comparison of the binding energy curves.
As further examples Fig. 4 reports the ∆E values also for one hydrogen-bond complex (H 2 O dimer) and one dipole-dipole complex (H 2 S dimer). In this case the binding curves are instead not reported because due to the relatively large value of E b the differences between the different methods are not easily visible on that scale. Inspection of the figure shows again that the present scaled approaches provide very accurate IFC energies over a broad range of intermolecular distances, with discrepancies from the reference well below the subchemical accuracy. The analysis of these results, together with those of Figs. 2 and 3, indicates that the S(R)-MP2 approach is overall slightly more accurate than SOS(R)-MP2 and SSS(R)-MP2, but all yield highly accurate results. For short distances, however, the methods necessarily face some limitations. In fact, as we mentioned above, at short inter-system separations the proportionality relation between (spin-resolved) MP2 and CCSD(T) may deteriorate due to the raising importance of intermolecular correlation (that is, the direct interaction of electrons and pairs "localized" on different systems, starts to be relevant). Moreover, even in absence of this phenomenon, i.e. for systems and distances R < R 0 where Fig. 2 suggests that the proportionality holds, the values of ∆E may be expected to rapidly increase (in absolute values) as the IFC energy is very large for R < R 0 , so that small deviations are largely magnified.
To have a more quantitative and global assessment of the quality of the scaling coefficients calculated within our simple non-empirical scaling scheme we report in Tab. II the mean absolute error, computed via Eq. (11), for several test systems and various methods, including the spin-component scaled methods optimized for intermolecular interactions [11], the conventional SOS-MP2 with c OS = 1.3, the unscaled CCSD and empirical 1/R 6 , 1/R 5 , and Lennard-Jones (LJ) potentials (fixed by fitting to CCSD(T) at R 0 ). Moreover, for comparison, we report the MAEs obtained by the popular PBE-D3 density functional (PBE exchange-correlation functional [25] supplemented with the D3 empirical dispersion correction [26]; note that the PBE functional is one of the best generalized gradient approximations for non-covalent interactions [27,28]). For DFT calculation the MAEs were computed by considering E b in place of E IF C into Eq. (12).
The values in the table are all very small, and are reported in unit of 10 −3 kcal/mol. The first column reports the binding energy computed at the CCSD(T) level: the considered systems span a wide range of intermolecular forces. In the last line also the relative mean absolute error (RMAE) is reported, which is computed as: where the sum runs over all the N systems and E b [i] is the binding energy of the i-th complex. The RMAE allows a fair global assessment of all the results. The results of Tab. II clearly show that the MP2-based methods using the here proposed non-empirical scaling can reproduce the reference CCSD(T) results with an ac-  (13)) is also reported in the last line. curacy below 1.5% in all systems (with a maximum MAE of 0.078 kcal/mol), with an accuracy two-three times better than the CCSD method (RMAE=3.68%). In particular, the present S(R)-MP2 has the lowest RMAE (0.53%) confirming its very high accuracy (the S(R)-MP2 MAE is always extremely small, lower than 0.006 kcal/mol). On the other hand, SCS(MI)-, SOS(MI)-, and SSS(MI)-MP2, despite performing very well in absolute terms (the deviation from CCSD(T) ranges from 0.003 to 0.1 kcal/mol) yield significantly worse results, with RMAE larger than 5%. A relatively small RMAE is obtained by the empirical 1/R 6 fit. This result however depends on the fact that in our test set simple dispersion dimers are predominant. The empirical potential works however much worse for other cases. In any case the RMAE of the 1/R 6 fit is five times worse than that of the S(R)-MP2 approach.

System E b SOS(R) SSS(R) S(R) SCS(MI) SOS(MI) SSS(MI) SOS
Finally, the PBE-D3 approach works fairly well, yielding a RMAE comparable to that of the SOS(MI)-and SCS(MI)-MP2 methods, but in any case more than five times larger than all the methods using the here proposed non-empirical scaling procedure. In addition, we note that the PBE-D3 method yields very good results for dispersion dimers, but not so good accuracy for hydrogenbond and dipole-dipole complexes.

A. Role of the reference distance
In this subsection we shortly analyze the role of the reference distance used to compute the scaling coefficients in Eqs. (8), (9), and (10). To this end we report in Tab. III the relative variation of the scaling coefficients and the MAEs (Eq. (11)) when the reference distance is changed from the equilibrium value R 0 toR > R 0 .
The tabulated results show that scaling coefficients change very little, as can be already inferred from Fig. 2. This provides a validation of the working hypothesis and indicates the robustness of our simple non-empirical scaling procedure. By inspecting the relative MAEs reported in Tab. III we see that the accuracy of the methods is not only preserved but in many cases is also improved. The improvement is anyway rather small and can be considered to lay within the numerical uncertainty related to the definition of the MAE (Eq. (11)). This result is important because it remarks the possibility to utilize the present non-empirical scaling procedure also in all those cases where the equilibrium distance is not known, by simply performing a single CCSD(T) calculations at any "reasonable" bond distance.

B. Benzene complexes
In this subsection we consider the application of our new scaling method to larger systems, namely the benzene-HCN T-shaped complex and the stacked benzene dimer.
The first of these systems is particularly interesting because the interaction has a mixed character involving both electrostatic and dispersion contributions in a variable proportion along the dissociation curve [18]. Its dissociation curves computed at the (spin-resolved) MP2 level as well as the corresponding scaled results are reported in Fig. 5. Reference CCSD(T) data from Ref. [18] are also reported. The scaling coefficients have been computed at the equilibrium distance of 2.34Å and are c OS = 1.27, c SS = 1.35, and c S = 0.66.
Consistently with the results discussed in the main body of section III, the non-empirical scaling methods all perform remarkably well, substantially improving over "bare" MP2, MP2(OS), and MP2(SS) calcu-  Similar results are found also by using scaling coefficients obtained considering a different reference dis-tanceR. ForR = 2.81Å, indeed, we found c OS =1.15, c SS =1.14, c S =0.57, and essentially the same performance for the binding energy (or equivalently the IFC energy). In fact, the scaling coefficients computed at the new distance are in good agreement, although slightly smaller, with the ones computed at the equilibrium distance. We remark that the small difference between the coefficients computed at different reference distances shall, in this special case, be partially traced back to the computational noise originating from the fact that the reference CCSD(T) values were not computed with the same set up as the MP2 ones, but rather extracted from the total binding energies (HF+CCSD(T) at the CBS limit) reported in Ref. [18]. Despite this small issue, the non-empirically scaled methods display an impressive robustness.
As a further test we consider a typical π-stacking complex. In Fig. 6 we thus report the dissociation curve of the stacked benzene dimer as computed at the (spinresolved) MP2 level and with the various scaled approaches. Reference CCSD(T) data from Ref. [18] are also reported. The scaling coefficients have been computed at the equilibrium distance of 3.765Å and are c OS = 1.39, c SS = 1.66, and c S = 0.76. Inspection of the figure shows that, similarly with the case of the benzene-HCN complex, all the scaled methods perform remarkably well yielding deviations within 0.2 kcal/mol from the reference CCSD(T) data. In particular, the non-empirical scaling procedure proposed here demonstrates to be able to correct well both the overbinding of the MP2 method and the strong underbinding of the MP2(OS) and MP2(SS) approaches.
Overall, the dissociation curves computed with our scaling approach agree very well with reference data and show that the method is very effective even for such "large" and "difficult" systems as the present ones. We note that this performance is remarkable especially for SOS(R)-and SSS(R)-MP2, because the MP2(OS) and MP2(SS) binding energies are almost degenerate, so that the corresponding scaling coefficient are almost identical. Thus, the same scaling is obtained for the two methods. However, this has only a minor effect on the SSS(R)-MP2 results at relatively small distances, while a good performance is observed for bond distances larger than the equilibrium one. This shows once more that the here presented scaling procedure is rather robust in many different situations.

C. Beryllium dimer
Finally, we consider in Fig. 7 the Beryllium dimer, which is a very difficult system where MP2 and CCSD fail to give even an approximate quantitative description of the dissociation curve [19]. For this system even CCSD(T) shows some limitations [19]. However, this high-level method can at least capture most of the features of the dissociation curve. Therefore, it will be used anyway as reference in the present example, although caution must be paid to the quantitative analysis of the results.
In this case, the computed scaling coefficients are c OS = 1.79, c SS = 3.03, and c S = 1.13. However, an accurate description of the dissociation curve is found only using the SOS(R)-MP2 method. On the contrary, methods based on the same-spin correlation fail quite evidently. This may be due to the fact that in this system, where static correlation is important, the stretching of the bond is not only promoting polarization and induction phenomena in each atom, but is also changing the actual "multireference" description of the system, influencing the same-spin correlation which includes antisymmetrized integrals. This fact indicates that the use of MP2(SS) correlation energy in scaling procedures has strong limitations for this peculiar system. Due to this limitations also the S(R)-MP2 method, which includes important contributions of same-spin correlation, cannot be accurate. Nevertheless, we remark that S(R)-MP2 and even SSS(R)-MP2 can compensate the limitations of the underlying MP2(SS) energy relatively well, thanks to appropriate values of the scaling factor granted by Eqs. (8), (9), (10). However, the independence of the scaling factor from the reference separation is partially lost and only a limited portion of the dissociation curve can be properly described when same-spin correlation is considered.
The SOS(R)-MP2 method with coefficients calculated in our scaled procedure provides for the Be 2 the best dissociation curve (as compared to CCSD(T)) significantly outperforming MP2, CCSD or SOS(MI)-MP2 methods.

IV. CONCLUSIONS
We have proposed a simple non-empirical procedure which can be used to calculate optimal scaling coefficients for (spin-resolved) MP2 calculations of the dissociation of noncovalent complexes. The applicability of the proposed method was demonstrated for a series of test systems, considering equilibrium and non-equilibrium reference distances, as well as for a notoriously difficult case as Be 2 . The obtained results show that the presented method works well in most cases, confirming the robustness of our hypothesis.
The proposed method is especially attractive because it allows to obtain a full dissociation curve of almost CCSD(T) quality at the cost of just a single CCSD(T) calculation. Moreover, it is conceivable to replace the costly CCSD(T) calculation with a focal-point analysis [29,30] (the ∆CCSD(T) procedure) to obtain an esti-mate of the CCSD(T) correlation energy from a cheaper calculation. In this way the computational cost is much reduced, losing only little on accuracy. Thus, future applications on large systems can be foreseen. Moreover, it might be also possible to replace the CCSD(T) reference with an alternative accurate treatment of the correlation, e.g. the CCSD[T] [31,32] or the FNO CCSD(T) [33] method, to optimize the accuracy or the computational effort.
Finally we remark that, although in this paper we focused only on single-parameter scaled MP2 methods, the methodology presented here can be extended in the future to treat the more general case of spin-componentscaled (SCS) MP2 calculations with two parameters (c OS and c SS ). To this end Eqs. (5)-(7) will be generalized to and an additional constraint will be needed for the scaling parameters (e.g. fixing the sum or the ratio of the parameters from theoretical considerations). We also have to stress that this approach can be useful for the verification of the quality of SCS coefficients used in different SCS-MP2 methods.