Subject: SUMMARY: BSSE in calc. of activation energies

From: "Kalju Kahn" <[email-masked]>
Date: Fri, 19 Mar 1999 21:21:43 -0800

Hello, everyone. Is about month now from the day I posted my
initial question. Thanks for good responses, I think I
learned a lot about this subject. Sorry to be so slow in responding,
had to wait to get some literature and then wait even more until I
understood some problems with BSSE. Of course, as many of
you noticed, my original equation was incorrect because
I assumed that there is no geometrical changes in reactant
geometries during the TS formation. Kalju.

Original question was:
#########################################################
From: Kalju Kahn
[email-masked]
##########################
I realize there has been some controversy about basis set
superposition error (BSSE), and that the current consensus
favors the view that full counterpoise procedure (CP) of
Boys and Bernardini is a good way to correct for BSSE.
(e.g. van Duijneveldt et. al. in Chem. Rev, 1994, 94, 1873-1885)

My question concerns application of BSSE correction via CP
to bimolecular reactions. Let's take gas phase SN2
displacement as an example. Here we have a two reactants, A and B,
forming first a bimolecular "reactants complex" (RC), followed by
transition state (TS):
          A + B -> (A...B)[RC] -> A-B[TS]
The activation barrier can be defined as an energy difference
between TS and RC. What is the correct way to calculate this
energy difference? (Assuming Born-Oppenheimer approximation,
and dealing with size-consistent methods, e.g. HF level of theory)

My own answer would be:
 If one accepts a view that energies of both RC and TS
 (relative to isolated reactants) are affected by BSSE,
 one should correct both according to the CP scheme, and now:
 E_act = E(TS,TS_basis_in_TS_geometry)-E(RC,RC_basis_in_RC_geometry)-
        E(A,TS_basis_in_TS_geometry) -E(B,TS_basis_in_TS_geometry)+
        E(A,RC_basis_in_RC_geometry) +E(B,RC_basis_in_RC_geometry)

 Obviously, the number of basis functions in TS_basis and RC_basis
 is the same, but the geometrical position of ghost orbitals relative
 to partner molecule are different in TS_geometry and RC_geometry.
 So, the last four terms do not cancel out.

This seems a right way to account for BSSE, but the above method
was quoted as "misinterpretation of the BSSE concept" by Dr. J. A. Sordo
in "Comment on "Ab initio investigation of internal rotation in the
ethylene-sulfur dioxide dimer"" (J. Chem. Phys., 106(14), 6204).

Could somebody kindly explain me how do I misinterpret the BSSE
concept here? Or do we need to be concerned about BSSE when
calculating ab initio activation barriers?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Summary of answers, stripped off from headers and tails:

##########################################################################
From: Frank Jensen
<[email-masked]>
#####################
	I am not familiar with that paper, but my opinion is
similar to yours, there is always a BSSE present when comparing
calculations using different geometries. For non-bonded interactions,
as your RC, the CP procedure is well-defined. For bonded species, however,
the CP procedure is much less clear. You may be interested in
the results in Chem. Phys. Lett. 261, 633.
	Frank

##########################################################################
From: Pr. Xavier Assfeld
[email-masked]
##############################
To be rigorous, you forget to take into account the energy of geometrical
deformation, i.e. :

 Your_formula + E(A,A_basis_in_TS_geometry) - E(A,A_basis_in_RC_geometry)
              + E(B,B_basis_in_TS_geometry) - E(B,B_basis_in_RC_geometry)

Moreover, Pr. Sordo et al. publish some papers about BSSE for this kind of
systems, showing that Delta_E_act is almost BSSE free :

1] J. A. Sordo, J. Chem. Phys. 1997, 106, 6204
2] V. M. Rayon, and J. A. Sordo, J. Phys. Chem. 1997, 107, 7591
3] D. Suarez, and J. A. Sordo, Chemical Communications, 1998, 3, 385

I also would like to point out that there are other (better) methods to
evaluate BSSE than couterpoise. For exemple have a look at the work of
Istvan MAYER about the Chemical Hamiltonian Approach. There is plenty of
papers dealing with that but here is a "recent" reference :

I. Mayer, and P. Valiron, J. Chem. Phys. 1998, 109, 3360.

#########################################################################
From: Huub van Dam
[email-masked]
########################
COUNTERPOISE recipe for complexes of monomers with internal geometry

Ra : geometry of  monomer A in complex
Rae : geometry of monomer A in its equilibrium
Rb : geometry of  monomer B in complex
Rbe : geometry of monomer B in its equilibrium
R, distance (geometrical parameters) between A and B
e.g. Ea(R,Ra,Rb)(AB) is the energy of A at geometry defined by
     (R,Ra,Rb) and the basis of AB

Generally a complete counterpoise corrected interaction energy is :

Delta E = Eab(R,Ra,Rb)(AB) - Ea(R,Ra,Rb)(AB) - Eb(R,Ra,Rb)(AB) +
          Ea(Ra)(A) - Ea(Rae)(A) + Eb(Rb)(B) -Eb(Rbe)(B)

The first line gives the counterpoise corrected interaction energy as
it would be for rigid monomers; The second gives the correction for
the deformation of the monomers

Applying this to a transition state (TS) and reaction complex (RC)
gives for a interaction energy with respect to isolated monomers for
both :

Delta E-TS= Eab(R-TS,Ra-TS,Rb-TS)(AB) - Ea(R-TS,Ra-TS,Rb-TS)(AB)
                                      - Eb(R-TS,Ra-TS,Rb-TS)(AB) +
            Ea(Ra-TS)(A) - Ea(Rae)(A) + Eb(Rb-TS)(B) -Eb(Rbe)(B)

Delta E-RC= Eab(R-RC,Ra-RC,Rb-RC)(AB) - Ea(R-RC,Ra-RC,Rb-RC)(AB)
                                      - Eb(R-RC,Ra-RC,Rb-RC)(AB) +
            Ea(Ra-RC)(A) - Ea(Rae)(A) + Eb(Rb-RC)(B) -Eb(Rbe)(B)

To get the barrier one have to substract the two :

Delta E-TS - Delta E-RC =
             Eab(R-TS,Ra-TS,Rb-TS)(AB) -  Eab(R-RC,Ra-RC,Rb-RC)(AB)
           - Ea(R-TS,Ra-TS,Rb-TS)(AB) - Eb(R-TS,Ra-TS,Rb-TS)(AB)
           + Ea(R-RC,Ra-RC,Rb-RC)(AB) + Eb(R-RC,Ra-RC,Rb-RC)(AB)
           +  Ea(Ra-TS)(A) +  Eb(Rb-TS)(B)
           -  Ea(Ra-RC)(A) -  Eb(Rb-RC)(B)

Where the last two lines (the extra part) is due to a different
deformation in transition state and reaction complex

See : F.B. van Duijneveldt
      "Basis Superposition Error"
      in : S. Scheiner (ed.) , "Molecular Interactions"
           Wiley (1997)

##########################################################################
From: Tanja van Mourik
[email-masked]
#####################

I think that your way to correct for BSSE is correct, assuming at least
that BSSE is only important *between* (and not within) A and B, as it is
in my ArHF example below.  If you would calculate your energy difference
simply as:

E_act = E(TS,TS_basis_in_TS_geometry)-E(RC,RC_basis_in_RC_geometry)

then you would not correct for BSSE at all!

Assuming the BSSE within the monomers can be neglected, you can compare with
a diatomic AB (like He2, for example). In this case, you also have to correct
for BSSE at each individual point of the potential energy curve.  Simply
substracting two AB dimer energies does not produce a BSSE-free energy
difference, since the uncorrected potential energy curve will be distorted.
Also in this case, for each different interatomic distance the ghost orbitals
are at different geometrical positions.

As another example:
For ArHF, which has 2 linear minima and a T-shaped TS, the correct way to
compute the stationary points would be, in my opinion, to substract from
the dimer energy the Ar and HF energies computed in the ArHF basis set, and
you would do this for the 2 different minima (Ar...HF and Ar...FH), as well
as for the TS. The barrier height is then the difference between the minimum
and TS interaction energy. This of course assumes that the BSSE within HF
can be neglected. For this example, I can reference 2 papers (I must confess
I am on the author list of both): JCP 107, 2451 (1997) and Adv. Quant. Chem.
31, 105, 1999. You already mentioned the review article by van Duijneveldt
et al.

##########################################################################
From: Qadir K. Timerghazin
[email-masked]
##########################
See good article on this subj:
Some difficulties in computing BSSE-corrected potential surfaces
of chemical reactions.
G. Lendvay 1 , I. Mayer Chemical Physics Letters 297 1998 365–373

May be you'll find there some answers on questions you put.

##########################################################################
From: [email-masked]
#######################
The AE is calculated between the TS1 and P1 structures, using the SAME
basis set (Sordo's eq1), therefore no CP is needed. CP is needed when
calculating, say, the P1(A..B) energy relative to the fragments (A and
B) where P1 has the functions on BOTH fragments.

##########################################################################
From: Steven.Creve
[email-masked]
##########################
Hello,
Maybe it would be interesting to come to a generalized view of BSSE
procedures, and summarize that to the list.

##########################################################################

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This is almost all, neglecting two private responses. What follows is my
attept to come up with a generalized view on BSSE and CP correction in
particular. I am not competent to discuss methods which exclude BSSE
a priori (Chemical Hamiltonian Approach and Symmetry Adapted Perturbation
Theory), and these apparently are not applicable for activation energies.

0. BSSE is not specific to molecular complexes. It affects geometries
   and energies of individual molecules when nuclei-centered basis
   sets are used. The effects on bond lengths are small but BSSE could
   affect the conformation of flexible molecule (Extra stabilization
   occurs when ends of the chain are next to each other)

1. The energy difference between the reactant complex and the
   transition state is affected by BSSE. I am not convinced that
   Delta_E_act is almost BSSE free.

2. The interaction energy can be partitioned into two components:
   (i) intermolecular interaction energy between two deformed reactants
   (ii) energy associated with deformation of isolated reactants.
   My initial equation neglected deformatiion energy, this is equal to
   assuption that reactants do not change geometry upon formation of
   reactant complex or transition state. While deformation of monomers
   may be negligable in reactant complex, it is very significant in the
   transition state. Assuming that there is no BSSE associated with
   deformation of monomers, the CP equation is:

  E_act = E(TS,TS_basis_in_TS_geometry)- E(RC,RC_basis_in_RC_geometry)-
          E(A,TS_basis_in_TS_geometry) - E(B,TS_basis_in_TS_geometry) +
          E(A,RC_basis_in_RC_geometry) + E(B,RC_basis_in_RC_geometry) +
          E(A, A_basis_in_TS_geometry) - E(A, A_basis_in_RC_geometry) +
          E(B, B_basis_in_TS_geometry) - E(B, B_basis_in_RC_geometry)

  or, rearranged:

  E_act = E(TS,TS_basis_in_TS_geometry)- E(RC,RC_basis_in_RC_geometry)+
  	  E(A, A_basis_in_TS_geometry) - E(A,TS_basis_in_TS_geometry) +
          E(B, B_basis_in_TS_geometry) - E(B,TS_basis_in_TS_geometry) +
          E(A,RC_basis_in_RC_geometry) - E(A, A_basis_in_RC_geometry) +
          E(B,RC_basis_in_RC_geometry) - E(B, B_basis_in_RC_geometry)

   It is not clear how important is BSSE in the deformation process
   E_def(A) = E(A,A_basis_in_TS_geometry) - E(A,A_basis_in_RC_geometry)
   Neither it is clear how to calculate BSSE for a deformation process.
   And, Rayon and Sordo (Theor. Chem. Acc. 99, 68-70(1998) have argued
   that when geometrical changes in fragments are large (as they are
   in formation of TS), then one should be careful when performing CP.

3. There appear to be serious conceptual difficulties when applying
   two-body CP correction to bimolecular reaction where group or
   atom is being transfered. The paper by Lendvay and Mayer (Chem.
   Phys.Lett. 297, 365-373) points out that assignement of transferred
   group to reactants becames ambiguous in TS, and different assignements
   will generally give different results.

4. From a practical point of view it could be better to ignore BSSE
   completery in calculation of activation energies. The arguments in
   favour of this neglect are:
   1) There are conceptual probles with performing BSSE
      correction (BSSE in deformation energies and p.3)
   2) The BSSE may not be too big compared to the current accuracy
      of describing transition states. This is probably fine when
      E_act is large.
   3) Correcting for BSSE by CP scheme involves four additional
      calculations using the large basis set and four additional
      calculations using smaller basis sets. No correction requires
      only two calculations with large basis sets (RC and TS).

5. Ignoring BSSE in TS calculations is expected to yield too tight
   transition states (in analogy with too short intermolecular
   distances on uncorrected potential energy surfaces between two
   reactants). Does anybody know of BSSE free (Chemical Hamiltonian
   Approach, or very large basis set LCAO) determinations of TS
   structures to evaluate this statement?

6. Ignoring BSSE in TS calculations is expected to yield too low
   activation energies.  BSSE in TS is larger than BSSE in RC since
   the two fragments are closer. I have very crude estimate for the
   MP2/AUG-cc-pVDZ basis set and RCl-HOH system, according to this
   BSSE near RC is about 0.6 kcal/mol and near TS about 2.6 kcal/mol.
   So the activation energy is underestimated by about 2 kcal/mol.
   (Note: these numbers have nothing to do with Sn2 displacement
   reaction, I havn't studied it). This underestimation is not
   terrible at HF level, since RHF activation energies are generally
   overestimated.

7. Some recent papers about BSSE:
   *Introduction to Computational Chemistry, pg. 172, Frank Jensen,
   Wiley, 1999
   *Basis set superposition error, F.B. van Duijneveldt, in "Molecular
   Interactions", ed. S. Scheiner, Wiley, 1997
   (These two still disagree whether CP is an approximate way
   or the method to deal with BSSE in weakly interacting molecular
   complexes)
   *Ab initio Calculations on Uracil-Water, van Mourik et al,
   J. Phys. Chem. A 103, 1611-1618 (1999). (Neat example of manual
   optimization to get a BSSE-free potential energy surface).
   *Comparison of basis set superposition error corrected perturbation
   theories for calculating intermolecular interaction energies.
   J. Comp. Chem. 20(2):274-283, 1999 (BSSE is significant at
   correlated levels)
   *On the validity of the counterpoise correction for the basis
   set superposition error including fragment relaxation terms
   Rayon & Sordo, Theor. Chem. Acc, 99, 68-70 (1998) (CP not valid)
   *Some difficulties in computing BSSE-corrected potential surfaces
   of chemical reactions. Lendvay & Mayer, Chem. Phys. Lett. 297,
   365-373 (1998) (CP correction results in discontinuous PES)

******************************************************************
Kalju Kahn                              [email-masked]
Chemistry Department                    tel: (805)-893-7158
University of California Santa Barbara
Santa Barbara, CA 93106
##################################################################