How does basis set superposition error change the potential surfaces for hydrogen-bonded dimers?

We describe a simple method to automate the geometric optimization of molecular orbital calculations of supermolecules on potential surfaces that are corrected for basis set superposition error using the counterpoise ~CP! method. This method is applied to the H-bonding complexes HF/HCN, HF/H2O, and HCCH/H2O using the 6-31G~d,p! and D9511~d,p! basis sets at both the Hartree–Fock and second-order Mo ” ller–Plesset levels. We report the interaction energies, geometries, and vibrational frequencies of these complexes on the CP-optimized surfaces; and compare them with similar values calculated using traditional methods, including the ~more traditional ! single point CP correction. Upon optimization on the CP-corrected surface, the interaction energies become more negative ~b fore vibrational corrections ! and the H-bonding stretching vibrations decrease in all cases. The extent of the effects vary from extremely small to quite large depending on the complex and the calculational method. The relative magnitudes of the vibrational corrections cannot be predicted from the H-bond stretching frequencies alone. © 1996 American Institute of Physics. @S0021-9606 ~96!04047-0#


I. INTRODUCTION
The importance of the basis set superposition error ͑BSSE͒ to the calculation of intermolecular interactions using ab initio calculations with basis sets below the Hartree-Fock limit has been appreciated for some time.The origin of this error lies in the possibility that the unused basis functions of the second unit in the associated complex may augment the basis set of the first unit, thereby lowering its energy compared to a calculation of this unit alone.The first unit will cause a similar error on the second.Although several other approaches to correcting this error have been discussed in the literature, the counterpoise ͑CP͒ correction proposed by Boys and Bernardi 1 has been the most popular means of correcting for BSSE.The CP method calculates each of the units with just the basis functions of the other ͑without the nuclei or electrons͒, using so-called ''ghost orbitals.''This method has proven to be somewhat controversial. 2 A problem with the normal use of the CP correction in accurate calculations of intermolecular interactions arises from the fact that the CP correction is usually added to the previously optimized geometry of the complex.In principle, since the BSSE causes the intermolecular interactions to be artifactually too attractive, the CP correction should make the complexes less stable.Consequently, the intermolecular distance will be greater when the complex is optimized with the CP correction included in the energetic expression.Furthermore, the vibrational force constants are generally reported on the uncorrected surface.This tends to make the intermolecular vibrations appear too strong, resulting in zero-point vibration energies ͑ZPVE͒, and vibrational corrections to enthalpy calculations that are incorrect.A striking result of this problem is the interaction energy of acetylene with ozone, which has a well-defined minimum but becomes repulsive after both CP and ZPVE corrections. 3everal examples of molecular orbital ͑MO͒ calculations where CP has been included in the optimization have been performed. 4However, the optimizations were done point by point as there are no options for this procedure in the common ab initio programs.Very recently, several authors have addressed the importance of relocating stationary points in the CP-corrected potential energy surface.Also, they have suggested the convenience of having an automated optimization procedure which uses the CP-corrected energy.2͑h͒,5 In this paper, we outline a procedure that allows automatic calculation of the CP correction within a normal ab initio optimization calculation including analytic first and second derivatives.This allows us to optimize geometries, locate transition states, and perform vibrational analyses on the CP-corrected potential service.We describe a program that automates GAUSSIAN 94 to perform this procedure, then provide several examples where the CP-optimized complexes differ significantly both energetically and structurally from the analogous structures calculated the traditional way, a͒ Fundacio ´n Banco Bilbao Vizcaya Visiting Professor at the Universitat Auto `noma de Barcelona.
The basic problem can be stated as the need to optimize E super CP where E super CP is described in Eq. ͑1͒, and E super represents the total energy of the supermolecular aggregate containing n monomeric units.Using the notation employed previously, 3͑d͒ the CP correction is stated in Eq. ͑2͒, where, the E m 's represent the energies of the individual monomers with the subscripts ''opt,'' and ''f '' denoting the individually optimized and the monomers frozen in their supermolecular geometries, and the asterisk ͑*͒ represents monomers calculated with ''ghost'' orbitals.Equation ͑1͒ can be rewritten as Eq.͑3͒.The CP-corrected interaction energy calculated at the CP-optimized geometry, E interaction CP , is expressed by Eq. ͑4͒,

͑5͒
In order to find a stationary point with respect to geometrical variation of the supermolecule, we require that the derivatives of E super CP with respect to all internal coordinates of  the supermolecule be zero ͓Eq.͑5͔͒.The energies of the optimized monomers are not a function of the supermolecular calculation, so their derivatives with respect to the geometrical parameters of the supermolecule are always zero.Equation ͑5͒ illustrates that the first derivative of either E interaction CP , or E super CP with respect to any internal coordinate, p j , can be expressed as a simple sum of the first derivatives of E super and the energies of each monomer frozen in its supermolecular geometry with and without ghost orbitals.Thus 2nϩ1 derivatives must be evaluated for each internal parameter.Since each p j will be the same for the supermolecule and the monomers, the derivatives at each geometric point are readily available from GAUSSIAN 94 or any other program that provides these derivatives.
Force constants and vibrational frequencies can be derived from the matrix of second derivatives.Each element of the Hessian matrix can be calculated in a manner similar to the first derivatives as indicated by In principle, not only the geometric variables, p j , will differ from those normally obtained from optimizations that do not contain CP corrections, but other molecular properties such as one-electron density, electric field values at nuclei, electrostatic potentials, dipole moments, polarizabilities, IR frequencies and intensities, etc., will differ as well.In general, any property that can be defined as a derivative of the energy can be calculated by using a variant of Eq. ͑5͒.
The purpose of this paper is twofold: first, to devise a procedure to build CP-corrected potential energy surfaces; and second, to apply this procedure to systems of chemical interest.Thus, in this paper we utilize the procedure described below to examine the CP-corrected surfaces of three complexes: ͑a͒ HF/HCN; ͑b͒ HF/H 2 O; and ͑c͒ HCCH/H 2 O.One should note that the effects of CP correction on potential energy surfaces can be considered to be similar to those due to basis set changes, the inclusion of electron correlation, the application of electric fields, etc.All of these will change the energies, equilibrium geometry, and curvatures at stationary points, i.e., harmonic frequencies.These three aspects will be analyzed in this paper.We report the CP-corrected surfaces of three complexes: ͑a͒ HF/HCN; ͑b͒ HF/H 2 O; and ͑c͒ HCCH/H 2 O.The first of these systems was previously studied by Bouteiller. 4͑a͒,4͑c͒ We shall compare our results with his.The second and third systems have been studied using more common procedures.HF/H 2 O has been found to be a nonplanar complex, 4͑c͒,12 while HCCH/H 2 O has been reported to be either planar 3 or nonplanar 6 depending upon the calculational methods used.In some cases, CP, applied in the traditional way ͑as a single point correction͒ lowers the energy of the planar below that of the nonplanar system.

II. METHODS
We realized the procedure outlined in Sec.I by writing a short segment of FORTRAN code designed to drive the energy optimizer, and several UNIX shell command files.These drive GAUSSIAN 94, the program used to perform the quantum chemical calculations of energies and analytical derivatives of the energy.For that purpose, we set up 2nϩ1 Z-matrices for the five types of geometrical inputs ͑in this paper, since nϭ2, this amounts to 5 Z matrices͒.The normal supermolecule Z matrix, as well as similar Z matrices containing either dummy atoms or ghost atoms, as appropriate, were used for monomers (m f ) and monomers with dimer basis set (m f *).Consequently, each of the 2nϩ1 calcula- tions yields similarly structured output, facilitating data manipulation.We used the direct inversion in the iterative subspace ͑GDIIS͒ method of Pulay, 7 which converges rapidly for systems having smooth potential energy surfaces with flat regions around energy minima, to optimize geometrical variables.The derivatives of the CP-corrected energy were taken from the GAUSSIAN 94 results using Eq.͑5͒.Particularly rigorous convergence criteria were applied ͑gradients were minimized to 10 Ϫ4 hartrees/bohr or hartrees/rad͒ to ensure proper location of minima in the flat surface necessary for meaningful low frequency vibrational calculations.In each calculation, we started the search on the CP-corrected sur-face with the optimized geometry on the uncorrected surface.For the normal, uncorrected calculations, we had to calculate the initial exact Hessians to ensure rapid convergence.However, use of GDIIS did not necessitate calculation of second derivatives.The unit matrix was used as Hessian.
We used GAUSSIAN 94 to calculate the harmonic frequencies from the second derivatives of the CP-corrected surface derived from the five different force constant matrices by application of Eq. ͑6͒.Minima on the CP-corrected surface were characterized using these frequencies in the usual manner.
Details of the code employed to automate the procedure used in this paper can be furnished upon request from the authors.

A. HF/HCN
The results for HF/HCN are presented in Tables I-III and Fig. 1.Our results are analogous to those reported by Boutellier.However, there are some notable differences which are probably due to the different basis sets employed.In addition to calculating the geometry and frequencies, we  compare the normal and CP-corrected potential surfaces at the minimum geometries for each.Note that the corrected interaction energies on the normal surface include corrections for both CP ͑E d ϪE a in Fig. 2͒ and zero-point vibration energy ͑ZPVE͒, while the corrected interaction energy on the CP-corrected surface only includes the ZPVE correction.A CP correction is recorded in Table I only to indicate the energy difference between the two surfaces at the CPcorrected minimum ͑E c ϪE b in Fig. 2͒.As expected, the H...N H-bonding distance is always longer on the CP-corrected surface.The largest change ͑0.092 Å͒ occurs for the MP2/D95ϩϩ͑d,p͒ calculations.As further expected, the interactions become more attractive ͑by from 0.1 to 0.9 kcal/mol͒.The MP2/D95ϩϩ͑d,p͒ calculation shows the largest effect.Here again, the MP2 calculations show the largest differences between the normal and CPcorrected surfaces.The correction is not always lower on the CP-corrected surface, as might be anticipated from the changes in the intermolecular stretching frequencies.While the H-bond stretch was shifted to lower frequencies in all four cases; as expected, increases in other frequencies more than overcame these shifts in both HF calculations.
The H-bonding stretching frequency calculated on the CP-corrected surface agrees remarkably well with the experimental 10 value of 168 cm Ϫ1 , but not with the reported harmonic CP-corrected value previously reported. 11The MP2/D95ϩϩ͑d,p͒ value for the F...N distance is in good agreement with the reported experimental values. 12

B. HF/H 2 O
The results for HF/H 2 O 13 are presented in Tables IV-VI and Fig. 3. Unlike the previous example, here the HF surfaces are more affected by CP correction than the MP2 surfaces.This can be seen from the differences in the normal and CP-corrected interaction energies.The CP correction is significantly diminished on the CP-corrected surfaces for the HF calculations, but relatively unchanged for the MP2's.The changes in the H...O H-bonding distances are most significant for the D95ϩϩ͑d,p͒ basis set in both HF and MP2 optimizations.The complex is predicted to be nonplanar in FIG. 3. Geometrical parameters for HF/H 2 O. all calculations.However, each of the geometries becomes closer to planar upon CP correction.This can be seen by the decrease in the angles a 1 and d 1 , both of which should be 90°for a planar complex.
Bevan et al. have determined structural parameters of HF/H 2 O from analysis of the microwave spectra of various isotopically labeled species. 14They have concluded that the complex contains a single H...O hydrogen bond with a F...O separation of 2.662 Å.While they emphasize the difficulty in distinguishing between C 2v ͑planar͒ and rapidly interconverting C s ͑pyramidal͒ geometries, they prefer a C s geometry based upon an analysis of the intensities of the vibrational satellites due to the thermal population of the lowest vibrational modes of the complex.The experimental 15 enthalpy of interaction of 6.2 kcal/mol is slightly greater than our best value of 5.4.

C. HCCH/H 2 O
The results of HCCH/H 2 O are presented in Tables VII-IX and Fig. 4. The energies and geometries obtained on the normal surface confirm the results previously published 3 for the HF calculations.The slight differences between the present MP2 results with those previously published result 3 from the use of the frozen core approximation in the current calculations.Other calculations on this system have been reported by Miller et al. 5 As in the other examples, the O...Hbonding distances increase, and the H-bond stretching fre-quencies decrease upon optimization with CP.In our previous report, 3 we noted that application of CP to an optimized nonplanar geometrical minimum and a planar saddle point had the effect of lowering the saddle point below the minimum.The present calculations show the optimized CPcorrected surfaces to be planar ͑or almost planar͒ in all cases, while the normal optimized geometries are nonplanar in the cases of both HF and MP2/6-31G ͑d,p͒, as previously reported.The difference in energy between the normal and CP-corrected surfaces is always less at the CP-corrected minimum.However, this difference is particularly large in the MP2/6-31 ͑d,p͒ case ͑going from 2.04 to 1.18 kcal/mol͒.All of the calculations are in reasonable agreement with the experimental 16 O...H distance.

III. DISCUSSION
The optimization of CP-corrected potential surfaces provides several interesting insights.Clearly, the CP-optimized geometry must be of lower energy than the normally optimized geometry plus the ͑single-point͒ CP correction.However, the difference in energy between these species can vary greatly.Since the CP correction must go to zero as the basis set approaches the Hartree-Fock limit, the two surfaces must converge at this point.However, it does not follow that any particular augmentation to the basis set will reduce either the energetic or geometric CP correction.For HF calculations, the variational principle dictates that the CP-corrected sur-  face must lie above the normal surface at all points.Thus, the CP correction will always decrease the interaction stabilization for calculations that satisfy the variational principle.As MP2 calculations are not variational, this may not be true for MP2 calculations.In these cases, the CP-corrected and normal surfaces may cross.Nevertheless, the corrected and uncorrected surfaces must converge at large intermolecular separations.
As previously noted, the BSSE provides a nonphysical attractive interaction.One might expect a correction for this interaction should cause the interacting molecules to separate and the frequency of the stretching vibration that separates the entities to decrease.We have observed these trends in all the cases studied here.Since vibrations involve the normal modes of molecules, these simple expectations may become incorrect in very complex systems.The simple conclusion that the ZPVE should decrease upon going from the normal to the CP-corrected surface 3 has proven incorrect in several instances ͑as noted above͒.The increase in the others more than counteracts the decrease in the H-bond stretching frequency.Moreover, the vibration that most represents the H-bond stretch is really a delocalized normal mode.The fact that some other higher frequencies increase is due to mixing of the intermolecular modes with the stiffer vibrations normally associated with the intramolecular modes.If the systems become sufficiently complex, the unique identification of a primary H-bond stretch may become obscure.
Equations ͑5͒ and ͑6͒ demonstrate that the derivatives of the energy with respect to any parameter of the system can be calculated as a simple sum of individual derivatives.Thus, any molecular property that can be written as a derivative of the energy with respect to some parameter can be calculated at any point upon the CP-corrected surface.For The question of whether the CP correction is the best method for correcting for BSSE has been extensively discussed in the literature.2͑f͒ The objection that it overcorrects for BSSE has often been disputed.The ambiguity of how it is performed has been noted.For example, adding HF molecules to a growing chain of HF's gives different CP corrections depending upon how one defines the interacting species to which ghosts orbitals are assigned.2͑d͒ We now can recognize that this is due to the normal practice of applying CP as a single-point correction.On a CP-optimized surface, all three CP methods used in the study of HF aggregates would necessarily converge to the same energy.
Other methods of correcting for BSSE have been proposed in the literature. 17However, none of these have been programmed to obtain analytical derivatives of the BSSEfree surface ͑except at the Hartree-Fock limit which is by definition, BSSE-free͒.The methodology that we have outlined here provides a simple procedure to calculate CPcorrected potential energy surfaces.We have implemented the procedure to run the GAUSSIAN 94 program.Analogous procedures to run with other MO packages can easily be developed.

IV. CONCLUSIONS
Since the CP-corrected energy at any point on the potential surface can be expressed as a sum of energies at the same point on the surface, the first and second derivatives of this energy can be expressed as a sum of derivatives of the individual energies with respect to each parameter.If the individual derivatives in the summation can be expressed analytically, it follows that the corresponding derivatives of the CP-corrected energy can also be expressed analytically.Calculating the derivatives of the CP-correct energy with respect to geometric parameters easily leads to the optimization of the geometry of aggregates on the corresponding surface.
Determining other derivatives, such as with respect to electric fields, can lead to the determination of other properties of aggregates on the CP-corrected surface.
Since CP correction always leads to reduction of the apparent attraction between molecules ͑by correcting for the nonphysical attraction due to BSSE͒, the intermolecular separation of H-bonding dimers increases upon optimization with CP.For HF calculations that are far from the Hartree-Fock limit, the magnitude of these effects are not yet predictable.Similarly, while the H-bonding stretching frequencies generally decrease with CP-corrected optimization, the ZPVEs do not always do so, as other vibrational frequencies can increase sufficiently to counteract this effect.

Fig. 2 .b
Point b in Fig. 2. c Point c in Fig. 2. d Point d in Fig. 2.

TABLE I .
Selected energetic results for HF/HCN.Total energies are in hartrees, all others in kcal/mol.
b Point b in Fig. 2. c Point c in Fig. 2. d Point d in Fig. 2. e Reference 11͑b͒.

TABLE II
a Reference 11͑a͒.b Reference 11͑b͒.c Reference 4͑c͒ ͑calculated͒.

TABLE III .
Comparison of the calculated vibrational frequencies ͑cm Ϫ1 ͒ on the normal and CP-corrected surfaces for HF/HCN.
a Reference 9.

TABLE IV .
Selected energetic results for HF/H 2 O.Total energies are in hartrees, all others in kcal/mol.
a Point a in

TABLE V .
Comparison of the geometric parameters of HF/H 2 O on the normal and CP-optimized surfaces.Distances in Å and angles in degrees.
a Reference 13.

TABLE VI .
Comparison of the calculated vibrational frequencies on the normal and CP-corrected surfaces for HF/H 2 O.

TABLE VII .
Selected energetic results for HCCH/H 2 O.Total energies are in hartrees, all others in kcal/mol.
a Point a in

TABLE VIII .
Comparison of the geometric parameters of HCCH/H 2 O on the normal and CP-optimized surfaces.Distances in Å and angles in degrees.
example, polarizabilities and hyperpolarizabilities can be calculated this way.We shall present some examples of this type of calculation in a forthcoming paper.