Counterpoise-corrected geometries and harmonic frequencies of N-body clusters : Application to „ HF ... n „ n Ä 3 , 4 ...

The differences between three previously defined counterpoise ~CP! schemes for removing the BSSE in molecular complexes formed by more than two subunits have been assessed by CP-corrected geometry optimizations and frequency calculations for the hydrogen fluoride trimer and tetramer. The types of the functional counterpoise ~FC! procedures included the site–site ~SSFC!, pairwise additive, and hierarchical Valiron–Mayer ~VMFC! schemes. The latter approach takes into account the basis set extension of the dimers in the trimer, dimers and trimers in the tetramer, etc. The number of different calculations required to apply this counterpoise scheme increases very rapidly with the cluster size. The symmetry of the chosen systems makes the test of this approach computationally feasible. All the optimizations and frequency calculations have been carried out automatically using a new program that generates the necessary input files and repeatedly calls a slightly modified version of a Gaussian link. The results show that geometrical parameters, zero-point vibrational energies, and redshifts computed on the CP-corrected potential energy surfaces differ considerably from those evaluated on the uncorrected surfaces. The structural and energetic properties obtained with the conventional SSFC procedure are almost identical to those predicted by the more costly and complex VMFC method. Hence, the former seems to be more appropriate in the present case. Furthermore, symmetry-adapted perturbation theory calculations show the importance of computing the interaction energies at the CP-corrected geometries. ©2003 American Institute of Physics. @DOI: 10.1063/1.1527011 #


I. INTRODUCTION
It has been recognized that the Basis Set Superposition Error ͑BSSE͒ ͑Ref.1͒ represents more than just an unbalance between the energies of the complex and its fragments in the computations the interaction energy. 2,3The BSSE is a phenomenon related to the LCAO approximation that affects the whole description of the complex, i.e., stationary points, 4 vibration frequencies, wave function, 5 etc.The appearance of the so-called a priori BSSE-correction methods 6,7 helped to understand this point of view and to recognize that the well known and widely used counterpoise ͑CP͒ correction 8 may be viewed as an energetic correction to the complex's energy.This interpretation has the advantage of permitting a straightforward definition of counterpoise-corrected derivatives of the energy 9 in order to obtain the stationary points on the counterpoise-corrected potential energy surface, vibration frequencies, dipole moments, spin-spin coupling constants, 10 etc.In the alternative, equally valid argument the CP correction is viewed as the energetic correction to the interaction energy.It is, in principle, possible to formulate a gradient-optimization procedure based on the minimization of the CP-corrected interaction energy in combination with the monomer relaxation effects.An optimization of such a quantity, henceforth referred as the stabilization energy, was suggested by van Lenthe et al. 3 This idea was not implemented analytically until Ref. 9.
The effect of BSSE on the structural properties has been demonstrated in the hydrogen fluoride dimer.A conventional MP2/6-31G** geometry optimization leads to a centrosymmetric C 2h stationary point.When correcting for BSSE by using either the Chemical Hamiltonian Approach ͑CHA͒, 6 or by optimizing with the counterpoise-corrected gradient, the correct quasilinear C s structure is obtained. 4he removal of the BSSE in molecular complexes composed of more than two fragments has not been extensively discussed in the literature.A few years ago, Turi and Dannenberg 11 pointed out the ambiguity of the counterpoise correction when studying growing chains of hydrogen fluoride.They showed that the BSSE computed for the addition of a new HF monomer to the (HF) n aggregate depends upon whether the incoming monomer is added to the H or to the F end of the aggregate.Hence, one can obtain different interaction energies for the same chemical process, which is unacceptable.They proposed the use of the counterpoise method by defining as many fragments as there are monomer subunits in the complex, with the BSSE defined as the difference between the energy of each monomer in its own basis set and that of the whole aggregate.
This method clearly solves the problem of the ambiguity of the CP correction but is unable to explain all the effects of the incoming monomer on the interaction ͑and BSSE͒ already present in the molecular aggregate.Valiron and Mayer 12 illustrated this deficiency with the example of three interacting H atoms described by the Slater 1s orbitals.In this particular case, the counterpoise scheme above will not predict any BSSE in the system whereas all the H-HϩH interactions bear some BSSE.Hence, there is a second-order BSSE due to the basis set extensions of all the H 2 descriptions within H 3 .Indeed, these diatomic basis set extensions are as natural as the atomic ones.The BSSE is due to the improvement of the description of the atoms ͑fragments͒ within the complex by using the other fragment basis sets to expand the genuine atomic ͑single fragment͒ contributions to the Hamiltonian.Analogously, the genuine diatomic ͑fragment pair͒ descriptions, including the respective interaction contribution within the atom ͑fragment͒ pair, are also artificially improved due to their expansion in the whole complex basis set ͑this is the particular case of the basis set extensions present on the H-HϩH interaction commented above͒.In this sense, the hierarchical partition of an aggregate into atomic ͑single fragment͒, diatomic ͑fragment pairs͒, etc., arises naturally.
One way to take into account those high-order BSSE effects within the counterpoise framework was first introduced by White and Davidson 13 and later generalized by Valiron and Mayer. 12They proposed a hierarchical counterpoise scheme for N-body clusters that treats the basis set extension effects of all the monomers, dimers, trimers, and so on, present in an aggregate.However, their proposed scheme was never tested in CP-corrected geometry optimizations.In a recent paper, Mierzwicki and Latajka 14 analyzed the behavior of these two counterpoise methods in the calculation of many-body interactions of Li͑NH 4 ) n and Li͑NH 4 ) n ϩ clusters at several levels of theory.They also used another, rather unusual scheme, introduced by Wells and Wilson, 15 where the counterpoise correction is carried out over pairs of fragments.
In the present paper we intend to go one step further.As commented above, in order to properly take BSSE into account, the counterpoise correction will henceforth be viewed as a correction to be added to the aggregate's description.This allows one to compute not only interaction energies, but also gradients and harmonic frequencies for the three different counterpoise schemes.Furthermore, the location of the stationary points on the BSSE-corrected PES is essential to obtaining the reliable counterpoise-corrected energies and to avoiding the artifacts which are sometimes referred to as an overcorrection. 16e wish to assess the differences between the various CP methods in terms of geometries, vibration frequencies, and interaction energies.For the first time, the full geometry optimizations using both the pairwise additive and the hierarchical counterpoise methods will be performed.The use of the hierarchical counterpoise scheme will help elucidate the effects of the high-order BSSE terms and will help to determine whether or not they can be neglected.The validity of the pairwise additive scheme will also be analyzed.
These methods will be applied to the hydrogen fluoride trimer and tetramer.The hydrogen fluoride clusters have received a great deal of attention lately.Recent experimental 17 and theoretical 16,18 -22 studies predict planar ring structures of C nh symmetry for the (HF) n 3рnϽ6 gas phase oligomers.In the case of the tetramer and pentamer, however, there is still some debate. 23,24X-ray and neutron diffraction experiments 25 have shown that solid HF tends to form infinite zig-zag chains with very large cooperative effects.Therefore, there must be an inversion of the relative stability of the cyclic and chain isomers as the aggregate grows.In this paper, both the cyclic and chainlike arrangements are considered in order to compare the BSSE effect for two structures where the importance of the cooperative effects is very different.Also, the high symmetry of the cyclic aggregates will allow us to perform hierarchical counterpoisecorrected geometry optimizations with a relatively large basis even for the tetramer.
Finally, we will perform symmetry-adapted perturbation theory ͑SAPT͒ analysis for the hydrogen fluoride trimer in order to gain insights into the nature of the interaction in this system.Even though SAPT is a genuinely BSSE-free methodology, its results depend on the choice of geometry that, in turn, is affected by the BSSE.One of the goals of this paper is to compare the SAPT results at the uncorrected and the counterpoise-corrected geometries and to assess the effects of BSSE-induced changes in the geometry upon the interaction energy components.
In the next section we briefly discuss the three different counterpoise methods used throughout the paper and derive the corresponding expressions for the counterpoise-corrected cluster energies by using a many-body partitioning of the energy of the aggregate.

II. METHOD
Let us consider first a dimer AB.The energy of this dimer at a given geometry with rigid monomers 26 can be expressed simply as where ⌬E AB represents the two-body interaction energy.According to the counterpoise philosophy, this value must be computed using the same basis set for all the terms involved, where the superscript AB means that the whole complex basis set is used ͑if no superscript is used, it is assumed that the energy is computed with the fragment's own basis set͒.In this way, the counterpoise-corrected dimer energy is recovered.Note that the one-body interaction energies, i.e., the fragment energies, are computed with the so-called monomer-centered basis sets ͑MCBS͒, whereas only the interaction energy term is computed with the dimer-centered basis set ͑DCBS͒.It is very important to point out that this is conceptually similar to the case of the a priori methods, such as CHA, where the diagonal ͑fragment-only͒ blocks of the Hamiltonian are maintained, and the BSSE-correction takes place only in the off-diagonal blocks ͑intermolecular interac-tion͒.
When a complex is composed of three interacting units, ABC, the energy of the system can be expressed as the sum of one-, two-, and three-body interaction energies, where the last term is due to the nonadditivity of the interaction.In order to obtain a counterpoise-corrected energy of the trimer, the three-body energy term must be computed following the standard counterpoise prescription, i.e., using the same, trimer, basis set for all the terms, The point now is to determine which two-body interaction energies must be used.If no counterpoise-correction is taken into account at all for those terms, the following expression is obtained by substituting Eqs.͑1͒ and ͑4͒ into Eq.͑3͒, There are three counterpoiselike terms related to basis set extension for all the dimers in the trimer, and also three terms corresponding to the basis set extensions of the monomers, which contribute to the counterpoise correction with opposite signs.The application of this scheme, however, yields meaningless results because the monomer basis set extensions are usually larger than those for the dimers, and hence the BSSE is negative.In other words, the energy of the supermolecule ͑and so the stabilization energy͒ is lowered upon counterpoise correction, which is unacceptable.Alternatively, one can consider using counterpoisecorrected two-body interaction terms in Eq. ͑3͒ but using the basis set of the whole trimer ABC ͑TCBS͒.

͑6͒
In this case the conventional counterpoise scheme is obtained, which includes only the basis set extensions of the monomers in the whole basis set.Wells and Wilson 15 called this approach site-site function counterpoise.However, the same considerations as those for the dimer case may suggest that the two-body interaction energy terms should be described with the respective DCBS basis set, in the same way as the monomer ͑one-body͒ contributions are expressed in the MCBS.According to these considerations, the counterpoise-corrected trimer energy will have the following expression: Rearranging the terms, the Valiron and Mayer's hierarchical counterpoise expression for the energy of a complex is obtained, The last three extra terms with respect to the conventional counterpoise scheme of Eq. ͑6͒, correspond to the differences, for each dimer in the aggregate, between the dimer interaction energy computed within the DCBS and TCBS.These effects will henceforth be dubbed second-order basisset extension effects.
Another counterpoise scheme previously proposed by Wells and Wilson, 15 the pairwise additive function counterpoise ͑PAFC͒ can also be obtained in a systematic manner like the other two schemes discussed above.In this case, the three-and higher-body interaction terms are not corrected according to the counterpoise scheme.Instead, only the twobody interaction energies are corrected by using DCBS.In the case of a trimer, the expression for the corrected energy can be easily obtained from Eq. ͑3͒, The main feature of this approach is that the whole complex basis set is never used for any subunit's calculation, except for the trivial case of a dimer.The counterpoise-correction is obtained by summing up, over all the subunits, the differences between the MCBS and all the different DCBS descriptions of a given fragment.For a given N-body cluster, the energetic difference between the MCBS and the wholecomplex basis set description of each fragment, as defined in the conventional counterpoise correction, is substituted by NϪ1 energy differences calculated using only the corresponding DCBS.Therefore, one might expect that this scheme may have problems in dealing with cyclic or highly packed clusters where the presence of many close-by DCBS representations for each fragment may lead to an overestimation of the BSSE.
The N-body cluster generalization of these three function counterpoise schemes is straightforward and the derivations can be found elsewhere. 12The final expressions for the socalled site-site, pairwise-additive and Valiron-Mayer ͑hier-archical͒ function counterpoise schemes, SSFC, PAFC, and VMFC, respectively, are as follows:

͑12͒
In Eq. ͑12͒ the third, fourth, and nth term on the right-hand side will be referred to as the second-, third-, and nth-order CP contributions.An important point is the scalability of these methods.Obviously, in the VMFC approach the number of needed calculations rapidly increases with the cluster size.The SSFC method needs 2N extra energy calculations.For the PAFC, N(NϪ1) DCBS calculations plus N MCBS calculations must be carried out, that is N 2 extra energy calculations.In case of the VMFC, it can be proved that the total number of the energy calculations is given by the relation, ).This means that the full hierarchical CP treatment of the nonsymmetric trimer through hexamer series would involve 19, 65, 211, and 665 energy evaluations, respectively.The treatment including only a second-order CPcorrection, ͓VMCP͑2͔͒ would involve N(Nϩ1) monomer plus 2( 2 N ) dimer calculations, that is a total of 2N 2 ϩ1 energy evaluations.In this case only 19, 33, 51, and 73 calculations are needed for the trimer up through the hexamer series.The use of the hierarchical scheme is clearly prohibitive even for relatively small oligomers, however, the high symmetry may enable such calculations.
Once the CP-corrected energy of an aggregate is obtained, the interaction and stabilization energies of the complex are obtained by subtracting the energies of the monomers computed at the CP-corrected complex geometry and isolated, respectively ͑note that since the BSSE is already taken into account in the complex energy, the monomers energies are computed with the MCBS͒, The next step is to characterize the stationary points on the CP-corrected potential energy surface of the complex, and to compute the vibrational frequencies.It has been shown that the gradient, Hessian, and in general any derivative of energy, can be obtained by a linear combination of all the terms properly differentiated. 9n all the geometry optimization and frequency calculations, we have used our code to automatically generate all the necessary input files and repeatedly call a slightly modified GAUSSIAN 98 package. 27We have rewritten the program 28 in order to accommodate the molecular symmetry, and to use the VMFC and PAFC methods.

A. HF n cyclic
The results of geometry optimizations for the cyclic trimer and tetramer are gathered in Table I.As shown in Fig. 1, both the trimer and tetramer structures are determined by three parameters: the intermolecular F-F (R F-F ) and intramolecular F-H (R F-H ) distances, and the angle HFF angle ͑␣͒.We have studied only the C 3h and C 4h configurations of the trimer and tetramer, respectively.
It is seen that in the trimer the intermolecular distance increases upon CP corrections, with the larger differences corresponding to the PAFC method.Upon CP-correction, the intramolecular F-H distance shortens by Ͻ0.01 Å.However, this difference is still larger than the variation of this distance with respect to basis sets.The cyclic nature of the complexes precludes large effect of BSSE correction on the angular parameter.In all cases, ␣ increases ͑by up to 3°͒, leading to a larger deviation from the triangular arrangement and hence to larger H-bond distances.
The addition of diffuse functions to the 6-31G(d,p) basis set dramatically decreases the effect of BSSE.The differences between the uncorrected and CP-corrected intermolecular distances decrease from more than 0.1 Å to Ϸ0.05 Å upon inclusion of diffuse functions.The medium polarized basis set, specifically designed to correctly describe intermolecular interactions, yields large BSSE.Indeed, both the uncorrected and the CP-corrected geometrical parameters are close to the values obtained with the 6-31G(d,p) basis set.
Adding another HF unit to the complex results in a shortening of the intermolecular distance by Ϸ0.06 Å.The intramolecular H-F distance decreases, whereas angle ␣ slightly increases.Cooperative effects are also evident in the energetics of the complex.The stabilization energy per hydrogen bond increases by more than 1 kcal/mol, thus providing the extra stabilization energy of Ϸ6 kcal/mol for the tetramer The similar effects are observed in the tetramer.The intermolecular F-F distance lengthens and the ␣ angle slightly increases.Even though the CP correction increases with respect to the trimer, the differences in geometrical parameters are comparable to those found in the trimer.
As for the energies, the CP correction to the trimer and tetramer energies is always overestimated at the uncorrected geometry.The CP-corrected stabilization energies computed at the uncorrected minima ͑single-point counterpoise calculation͒ are smaller than those evaluated at the corresponding CP-corrected stationary point.The differences in the case of the trimer range from Ϸ0.1 kcal/mol for the 6-31ϩϩG(d,p) basis set to more than 1 kcal/mol for both the 6-31G(d,p) and the Sadlej basis set.In the tetramer these differences are twice as large.It is important to note that after the CP-correction the basis set dependence of both the calculated stabilization energies and geometrical parameters decreases.The uncorrected stabilization energies obtained with the 6-31G(d,p) and the Sadlej basis sets are far too large.All the CP-corrected values are within 3 and 5 kcal/mol for the trimer and the tetramer, respectively.The same situation has been observed in previous studies of weakly bound complexes. 4,31he differences between the SSFC and the VMFC corrected values are still appreciable for the smallest basis set 6-31G(d,p).The intermolecular distance is again the most sensitive geometrical parameter.The inclusion of high-order terms in the CP method leads to larger intermolecular distances, the differences being 0.015 and 0.024 Å for the trimer and tetramer, respectively.The effects on the cluster energy are much more evident.The SSFC method leads to the stabilization energy higher in magnitude by Ϸ1 and 2 kcal/mol, for the trimer and tetramer, respectively.In the tetramer case, the inclusion of the third order CP-correction terms, VMFC͑3͒, feasible here due to the high symmetry, shows no significant effect on either the geometry or the energy of the complex, provided the basis set used is flexible enough.Only in the 6-31G(d, p) basis set, the intermolecular distance still increases by 0.007 Å and the stabilization energy decreases by 0.3 kcal/mol.Fortunately, in the remaining, more diffused basis sets, the high-order correction VMFC terms have practically no effect on both geometrical parameters and energies.It is worth to point out two facts.First, in the 6-31ϩϩG(d, p) tetramer calculation, the inclusion of high-order CP-correction terms induces a smaller BSSE.The intermolecular distance slightly shortens, and the stabilization energy increases by 0.16 kcal/mol upon correction.Even though it is a rather unexpected result, it should be emphasized that the high-order terms in the VMFC method can actually be of opposite sign.The fact that the dimer correction term is negative does not mean that its energetical description is better with the DCBS than with the TCBS.Instead, it is the dimer interaction energy, which is larger ͑more negative͒.The reason why this happens is that the lowering of the monomer energies is larger than the dimer counterparts as the basis set increases.Second, it is remarkable that, despite the large BSSE exhibited at the monomer level by the Sadlej basis set, the effects of high-order CP correction are rather insignificant.It is confirmed again that a basis set should not be considered bad or unbalanced just because it bears a large BSSE.Indeed, we will show that the Sadlej basis set provides very accurate results, provided the BSSE is properly taken into account.
Finally, the PAFC method leads to larger CP correction than SSCP and VMCP in all the cases except for the calculations involving the 6-31G(d, p) basis set.In this case, the PAFC results show that the method seems to mimic the effect of the high-order CP-correction terms.However, the differences observed in both the geometrical parameters and cluster energies when using more suitable basis sets make this method not advisable.
It has been shown how the CP methods affect the location of the stationary points at the PES.Obviously, the harmonic frequencies on a corrected and uncorrected PES are expected to be different as well because of the two main factors.First, the geometrical parameters of the stationary points are different, so the differences in frequencies are predicated on how large is the CP-correction on geometry.Second, the higher-order derivatives of the CP-correction term are non-zero, so the CP-corrected second derivatives are expected to differ from the derivatives evaluated for the uncorrected PES.By comparing the corrected and uncorrected frequencies at the CP-corrected stationary points one can determine whether the rather expensive calculations of the CP-corrected second order derivatives are needed.The differences between the frequencies properly computed on the corresponding uncorrected and CP-corrected PES contain both the ''geometrical'' and the ''differential'' factors.
The uncorrected and CP-corrected harmonic frequencies for the cyclic HF trimer calculated for the three basis sets and CP methods are shown in Table II.In all the cases, the uncorrected low frequencies are overestimated whereas the frequencies of the two stretching modes are underestimated with respect to the CP-corrected values.For the 6-31G(d,p) and Sadlej basis sets the differences between the uncorrected and the CP-corrected frequencies range from 60 cm Ϫ1 for the lowest frequency to more than 200 cm Ϫ1 for the frequencies labeled 1 and 2 .In general, the BSSE modifies the low frequencies by 10%-25%.The differences in the fundamental stretching frequency 4 are Ͼ100 cm Ϫ1 .As expected, the 6-31ϩϩG(d,p) frequencies are modified very little.The maximum differences are Ϸ50 cm Ϫ1 , even for the most sensitive frequencies 1 and 2 .
The inclusion of the second-order CP-correction seems to induce no appreciable changes in the frequencies.Only for the smallest basis set the frequencies are further shifted by up to 7% with respect to the SSFC values.The PAFC frequencies are very similar to both the SSFC or VMFC values.
Out of the two mentioned factors affecting frequencies, the ''geometrical'' one is clearly more important.In general, the ''geometrical'' and ''derivative'' factors act in opposite directions on the frequency shift.For instance, the uncorrected 6-31G(d,p) lowest frequencies decrease when computed on the CP-corrected PES, but tend to increase when computed using the CP-corrected second derivatives.However, the opposite trend is observed for the Sadlej basis set.
The effect on the zero-point vibrational energy ͑ZPVE͒ correction can be explained also on the basis of these opposite effects.For the small basis sets, the ZPVE decreases when computed at the CP-corrected PES but then increases when using CP-corrected second derivatives.The opposite occurs for the Sadlej basis set, which shows the largest effect in the ZPVE correction.In this case, the uncorrected ZPVE correction of 5.72 kcal/mol decreases to 4.84 kcal/mol when computed on the VMFC PES and further reduces to 4.54 kcal/mol upon correcting the second derivatives.
The results obtained in the Sadlej basis set ͑see Tables I and II͒ for the C 3h trimer may be compared to the results obtained by using an empirically refined SC-2.9ϩHF3BG potential of Quack, Stohner, and Suhm. 20Our CP-corrected binding energy and ZPVE correction are slightly smaller than the reference results.Also, the R F-F distance is too long reflecting the fact that this basis set ͑and the level of correlation treatment͒ is expected to underestimate the attraction in this system.The predicted redshift of the HF stretching frequency ͑249 cm Ϫ1 ͒ agrees very well with the harmonic value obtained ͑250 cm Ϫ1 ͒ by Quack, Stohner, and Suhm.It should be stressed that the CP correction appears to be essential for the calculation of this quantity.The CPuncorrected value of redshift is severely overestimated ͑371 cm Ϫ1 ͒ in this basis set.
Recently, Liedl 16 studied the concerted hydrogen exchange process of the HF trimer at the MP2/aug-cc-pVXZ, Xϭ2,4.He found that the uncorrected energies for the C 3h and D 3h structures were less basis set dependent than the counterpoise-corrected ones.On this basis, he claimed the counterpoise-corrected results were useless.However, this assertion was based on the misinterpretation of his own data.First, the uncorrected values for the minimum and the transition state indeed show a weaker basis set dependence, but they lack a monotonic trend.Therefore these values cannot be used to properly extrapolate to the basis set limit, and his complete basis set ͑CBS͒ limit is completely arbitrary.The CP-corrected values, on the other hand, vary monotonically and a CBS extrapolation can be carried out.
Second, the CP-corrected values, contrary to Leidl's assertions, provide a much better description of barrier height dependence on the basis set than the uncorrected values.For example, the reported barriers for the hydrogen exchange for the aug-cc-pVXZ, Xϭ2 -4 series are 20.17,17.89, and 18.61 kcal/mol, and 23.83, 20.48, and 20.21 kcal/mol for the uncorrected and the CP-corrected values, respectively.The best estimate of this barrier using explicitly correlated coupled cluster calculations is 20.33 kcal/mol.The CPcorrected CBS values are clearly closer to this value than the uncorrected ones Finally, his chief criticism of the usefulness of the CPcorrection was based on the results obtained for the transition-state D 3h structure of the trimer where the monomers are highly stretched ͑with a fragment relaxation of Ϸ60 kcal/mol͒.It should be emphasized that the transition state structures are much more sensitive to the basis set ͑and to the inclusion of correlation effects͒ than the equilibrium structures.Therefore this criticism is completely unwarranted.
Another effect that was not taken into account by Liedl was the use of CP-correction in the geometry optimization.Our calculations performed at the MP2/aug-cc-pVDZ indicate that the stabilization energy for the C 3h structure computed at the CP-corrected stationary point is Ϫ13.30kcal/ mol, i.e., 0.16 kcal/mol lower than the single-point CPcorrected value.It should be mentioned that the effects of BSSE on geometry are much lower in the aug-cc-pVDZ basis set than in other basis sets of comparable size. 4

B. "HF… n linear
The zig-zag linear structures of both the HF trimer and tetramer were also studied at the same level of theory.The results obtained for the uncorrected and the CP-corrected geometry optimizations are shown in Tables III and IV for the trimer and tetramer, respectively.The definitions of geometrical parameters are depicted in Fig. 2. In the this case, only the SSFC method was used for the corrected optimization, since including the second-order CP or full VMFC corrections would involve 33 and 65 gradient calculations, respectively.
The observed trends are similar to those obtained for the cyclic structures.Upon the CP correction, intermolecular distances increase while the intramolecular HF bonds shorten, leading to a weaker interaction.The differences between the uncorrected and CP-corrected parameters are of the same order than for the cyclic complexes, except for the intermolecular bond angles.In this case, the 6-31G(d,p) basis set poorly describes the directionality of the interaction.The intermolecular bond angles ␣ 1 and ␣ 2 are overestimated whereas ␤ 1 and ␤ 2 are clearly underestimated by up to 20°.This is not surprising since in (HF) 2 the uncorrected geometry optimization at this level of theory leads to a spurious cyclic structure. 4The corresponding CP-corrected optimizations, however, overcome this problem.Indeed, the CPcorrected angular parameters are in good agreement with the values obtained using more flexible basis sets.Again, the effect of BSSE is minimized by the addition of diffuse functions.The discrepancies between uncorrected and CPcorrected values are Ͻ1°for the angles and 0.06 Å for intermolecular distances.The Sadlej basis set bears the largest BSSE.However, whereas the intermolecular distances are underestimated by more than 0.12 Å in the absence of the CP-correction, the angular features of the complex are well described even at the uncorrected level.
Regarding the CP-corrected stabilization energies, the linear structures are about 4 and 8 kcal/mol less stable compared to the cyclic trimer and tetramer, respectively.The CP correction in the trimer ranges from 1.6 and 2.76 kcal/mol for the 6-31G(d,p) basis set to 7.95 and 12.31 kcal/mol for the Sadlej basis set for the trimer and the tetramer, respectively.The dependence of the CP correction on the geometry seems to be less important than for the cyclic case.However, the BSSE can still be overestimated by up to 1 kcal/mol ͑depending on the basis set͒ when computed at the uncorrected geometry.The cooperative effects are obviously less important than in the cyclic structures.However, the addition of another HF unit to the linear trimer enhances the stabilization energy per hydrogen bond by Ϸ0.5 kcal/mol The performance of the two first-order CP methods, i.e., SSFC and PAFC, is similar.Both methods modify the values of the geometrical parameters in the same direction, even though the PAFC method leads to larger corrections than the remaining treatments.
Our results show again that the differences between the SSFC and VMFC approaches are only appreciable in the context of small basis sets.The inclusion of second-order CP terms in the 6-31G(d,p) basis set increases the intermolecular distances and angles by up to 0.01 Å and 0.5°, whereas it   V. Description of SAPT corrections ⑀ (i j) ͑where i and j correspond to the interaction and the intramonomer correlation operators, respectively͒ which are implicitly present in the two-and three-body supermolecular Møller-Plesset interaction energy terms ⌬E int at the SCF level and in the second order.
The inclusion of high-order CP-correction terms was computationally feasible only for the trimer, therefore only the SSFC method was used to compute the CP-corrected geometry of the linear tetramer.Nevertheless, we performed a single-point second-order CP-correction at the SSFC corrected geometry with the 6-31G(d, p) basis set.The value of CP correction increased by 0.64 kcal/mol.
As pointed out recently by Rinco ´n et al., 32 the open chain structures for the HF trimer and tetramer are first-order saddle points connecting two equivalent cyclic configurations.Our results are consistent with their findings in all the cases.The CP-correction does not change the topology of the PES in any case, not even for the 6-31G(d,p) basis set, where the effect on the geometry is very large.
In the second-order of Møller-Plesset theory, one of the dominating SAPT terms is the second order dispersion term ⑀ disp (20) .This term is additive and thus only contribute the twobody components.The first nonadditive dispersion component appears in the third order of perturbation theory as the ⑀ disp (30) term.This term usually dominates the nonadditivity of the ⌬E int (3)  The SAPT contributions were calculated for the cyclic and linear trimer at three basis sets.The calculations were carried out at uncorrected and VMFC-corrected minima of the trimers and the results are shown in Tables VI-VIII.As a trimer of highly polar molecules, (HF) 3 is dominated, at the level of two-body interactions, by the electrostatic attraction.These effects are counterbalanced, to a certain degree, by the repulsive exchange effects.The two-body induction effects  (20) ͑2-body͒ Ϫ4.07 Ϫ3.47 0.60 4.5 ⑀ disp (30) ͑3-body͒ 0.00 0.00 0.00 0.0 ⑀ es (10)  are also important.The two-body dispersion interaction is the third in importance.At the level of three-body interactions the bulk of nonadditive interaction originates from the SCF-deformation term.The exchange nonadditivity is quite small while the three-body dispersion is nearly zero.The cyclic configuration is stabilized over the linear one at the level of two body interactions, because of more favorable electrostatic and induction effects in the cyclic arrangement.The three-body terms also favor the cyclic structure.Calculations of SAPT terms at two different geometries, one uncorrected and a counterpoise-corrected ͑at the VMFC level͒ are displayed in Tables VI-VIII.In the 6-31G(d,p) basis set the evaluation of SAPT terms at the uncorrected minimum geometry leads to large discrepancies in SAPT terms resulting in considerable percent errors in the electrostatic, exchange, and induction terms.For example, in the cyclic configuration, the errors in these terms range from 50% to over 100%.These errors are reduced in the 6-31ϩϩG(d,p) basis set to 15%-30%.The results obtained in Sadlej basis set also indicate large discrepancies in SAPT terms derived using these two geometries.Although, this basis set produces reliable values of SAPT terms, it also generates large values of BSSE that result in large distortions of geometrical parameters.It is worth noting that because of this difference in geometry it is possible to obtain a false picture of the interaction energy composition if the calculations are done for the uncorrected minimum.For example, in Sadlej basis set the Heitler-London interaction energy, ⌬E HL is negative ͑Ϫ3.11 kcal/mol͒ in the VMFC-corrected minimum whereas in the uncorrected minimum it has a repulsive ͑5.71 kcal/mol͒ value for the cyclic configuration.A similar sign reversal of ⌬E HL also takes place in the linear configuration.

IV. CONCLUSIONS
This paper has examined the methods of BSSE-free geometry optimization and frequency calculations in clusters larger than a dimer.Three different counterpoise schemes have been critically examined.It has been shown that the counterpoise-corrected supermolecule energy can be easily obtained in all the cases by using the many-body partitioning of energy.The expressions for the so-called site-site, pairwise-additive, and hierarchical function counterpoise are reproduced.
A computer program for such calculations using three counterpoise schemes has been coded and tested for gradient optimizations and harmonic frequency calculations of the HF trimer and tetramer.The high symmetry of the cyclic complex has made possible the study of the cyclic HF tetramer within the hierarchical CP approach.
Calculations performed in three different basis sets 6-31G(d,p), 6-31ϩϩG(d,p), and Sadlej basis sets indicate that only the latter two are suitable to judge the performance of the CP-procedures.The first basis set leads to results that are too erratic.This basis set performs very poorly even after the CP-correction, and should be avoided in the studies of intermolecular interactions.
Generally, both SSFC and VMFC lead to very similar values of the CP correction.The PAFC scheme leads to the larger, most likely overestimated, CP values.Therefore, in our opinion it does not represent a valid correction scheme.
A comparison of the results obtained with the conventional ͑SSFC͒ and the hierarchical ͑VMFC͒ CP methods indicates that, except for unsuitable basis sets such as the 6-31G(d,p), the high-order BSSE effects are not important.We conclude that the conventional CP scheme is clearly preferred in this case, mainly due to the extra computational cost required by the application of the VMFC ͑which is practically inapplicable in larger clusters with low symmetry͒.
The presented SAPT results show that already at the level of two-body interactions the cyclic configuration is stabilized over the linear one.The three-body terms also favor the cyclic structure.The large differences between the results obtained at the uncorrected and counterpoise-corrected geometries underscore the need for performing the analysis of the interaction energy at the counterpoise-corrected minimum geometries.

TABLE I .
Geometrical parameters ͑Å, deg͒, total ͑a.u.͒, and stabilization energies ͑kcal/mol͒, and CP corrections ͑kcal/mol͒ for the cyclic HF trimer and tetramer in several basis sets and counterpoise methods.See Fig.1for the definition of the geometrical parameters.
FIG. 1. Geometrical parameters of the C 3h and C 4h cyclic hydrogen fluoride trimer and tetramer.

TABLE II .
Harmonic frequencies ͑cm Ϫ1 ͒, ZPVE correction ͑kcal/mol͒ to the stabilization energy and frequency shifts ͑cm Ϫ1 ͒ for the cyclic HF trimer.Values in parentheses correspond to the uncorrected harmonic frequencies and ⌬ZPVE computed at the CP-corrected geometry.Redshifts calculated with respect to the monomer H-F stretching frequencies obtained at each level of theory ͓4193.4,4118.7, and 4082.4 cm Ϫ1 for the 6-31G(d,p), 6-31ϩϩG(d,p), and Sadlej basis sets, respectively͔.

TABLE IV .
Uncorrected and SSFC-corrected geometrical parameters ͑Å, deg͒, total ͑a.u.͒ and stabilization energies ͑kcal/mol͒, and CP corrections ͑kcal/mol͒ ͑single-point CP in parentheses͒ for the linear HF tetramer in several basis sets.See Fig.2for the definition of the geometrical parameters.

TABLE VI .
2-body and 3-body SAPT contributions ͑kcal/mol͒ to the interaction energy for the cyclic and linear (HF) 3 with the 6-31G(d,p) basis set at the uncorrected and VMFC-corrected geometries.The difference between the corrected and uncorrected values is also reported.͑The total 2ϩ3-body effect, where not indicated; see text for details.͒ supermolecular Møller-Plesset term.The physical sense of the SAPT corrections considered in this work, and their correspondence to the supermolecular Møller-Plesset terms are summarized in Table V ͑see also Ref. 34͒.The basis set dependence of two-body SAPT terms in (HF) 2 was analyzed previously ͑see Ref. 34͒.Monomer properties of HF in the Sadlej basis set can be found in Ref. 18.

TABLE VII .
Energetic SAPT contributions ͑kcal/mol͒ to the interaction energy for the cyclic and linear (HF) 3 with the 6-31ϩϩG(d,p) basis set at the uncorrected and VMFC-corrected geometries.

TABLE VIII .
Energetic SAPT contributions ͑kcal/mol͒ for the cyclic and linear HF trimer with the Sadlej basis set.