\documentstyle[12pt,rotate]{article}
\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\def\rf#1{${}^{#1}$}
\def\frm#1#2{\raisebox{0.2ex}{#1}\hspace{-0.1em}/\hspace{-0.1em}%
\raisebox{-0.6ex}{#2}}%
\begin{document}
\begin{singlespace}
\centerline{\Large Proton Affinities Calculated}
\vskip 0.1in
\centerline{\Large by Traditional {\sl ab initio} Approaches}
\vskip 0.1in
\centerline{\Large and by Density Functional Methods.}
\vskip 0.3in
\noindent
\addtocounter{footnote}{1}
{\large \bf Jan K. Labanowski{}$^\ast$\footnote{Ohio Supercomputer Center},
Ronald A. Hill\footnote{Ohio State University}\footnote{Present
address: Northeast Louisiana University, School of Pharmacy, Monroe,
Louisiana 71209-0470},
David J. Heisterberg\footnotemark[2],
Duane D. Miller\footnotemark[3]\footnote{Present address: College of Pharmacy,
University of Tennessee, 847 Monroe Ave, Memphis, TN 38163},
Charles F. Bender\footnotemark[2],
Jan W. Andzelm\footnote{Cray Research}\footnote{Present address:
Biosym Technologies,
9685 Scranton Road, San Diego, CA 92121-2777} }
\vskip 0.2in
\noindent
{\sl Contribution from the
Ohio State University, College of Pharmacy, 500 West 12th Ave.,
Columbus, Ohio 43210-1291,
Ohio Supercomputer Center, 1224 Kinnear Rd., Columbus, Ohio 43212-1163,
and
Cray Research, Inc., Industry Science and
Technology Department, Eagan, MN 55121}
\\
\\
{\bf Abstract:} Proton affinities for ethanol, formate, formic acid, methanol
and methylamine calculated by {\em ab initio} Hartree-Fock and MP2 approaches,
and various Density Functional Theory (DFT) methods are compared. Local Spin
Density (LSD)
calculations did not produce satisfactory results. However, Becke/Perdew
nonlocal gradient correction applied as a perturbation to the LSD energy,
or in a self-fconsistent
fashion dramatically improve calculated affinities yielding values of the same
quality as {\em ab initio} MP2 approach.
\end{singlespace}
\newpage
\section*{Introduction}
\renewcommand{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}
Protonation reactions, i.e., A $+$ H$^+$ $\rightarrow$ AH$^+$, are
among the most important in chemistry and biology.
Protonation/deprotonation is the first
step in many fundamental chemical rearrangements and in most
enzymatic reactions.
Two quantities are used to characterize the ability of a molecule
in the gas phase to accept a proton. The gas phase basicity is the
negative of the free energy change associated with the reaction.
The more frequently used index,
the proton affinity, is the negative of the enthalpy change at standard
conditions. Experimental
determination of these parameters is not easy (for an excellent review
on this topic see Dixon and Lias\rf{1}), and with the phenomenal growth in
computer power in recent years,
much attention has been given to the possibility of
calculating these parameters by quantum methods. {\sl Ab initio}
approaches are very successful in providing reliable values of proton
affinities and gas phase basicities for small molecules even at lower
levels of theory\rf{2}. However, due to computational
expense, application of {\sl ab initio} methods to the estimation of
proton affinities is still impractical for larger molecules.
Semiempirical methods such as AM1, MNDO and PM3, are not consistently
reliable
in calculations of proton affinities as shown by Ozment and
Schmiedekamp\rf{3}.
The recent progress in the Density Functional Theory (DFT)
approaches (for review see refs 4-6)
make this method another candidate for reliable calculation of proton
affinities, however, the performance of the method in this field
is still mostly untested. This prompted
us to analyze its performance on a few representative
molecules spanning a wide range of proton affinity values.
DFT methods are computationally less demanding than correlated
{\em ab initio} approaches and formally scale with the size of the molecule
as $N^3$ or $N^4$, depending on implementation.
For the the simplest Hartree-Fock {\sl ab initio} approach the scaling
is $N^4$. Moreover, the advantage of DFT methods is that they
should in principle include electron correlation energy via
the correlation/exchange potential, while the Hartree-Fock
approach by definition does not include this component of energy.
The simplest of the routinely-used {\sl ab initio} correlated approaches,
based on the second order many body perturbation theory,
MBPT(2) (frequently called second order M\o{}ller-Plesset
theory -- MP2), recovers only a portion of the correlation energy and
scales as $n\times N^4$ ($n$ is the number of
occupied molecular orbitals).
The major weakness of DFT approaches is that the exact mathematical
form of the exchange-correlation potential is not known. For that reason
approximations are used.
The most popular is the
Local Spin Density (LSD) approximation, which simply assumes that the
exchange-correlation potential dependence upon charge density
is represented by the functional form
found for the homogenous electron gas. This approximation works well in
many cases; however, it suffers from severly underestimating electron
exchange\rf{29} and overestimating the correlation energy\rf{30}. Appropriate
corrections to the local density approximation are therefore sought,
and several of such schemes already exist\rf{42,43}.
\section*{Computational Methods}
Proton affinity $P({\rm A})$ is defined as the negative of the molar enthalpy
change at 298.15 K ($-\Delta H_{\rm Rn}$) for the reaction
A $+$ H$^+$ $\rightarrow$ AH$^+$, or in the case of anions,
A$^-$ $+$ H$^+$ $\rightarrow$ AH.
To calculate these values from theory for gas phase reactions we may,
in most cases, obtain adequate results assuming ideal gas behaviour:
$\Delta H_{\rm Rn} = \Delta E_{\rm Rn} - RT$. Assuming that there is only
one conformer present, the energy $E(T)$ of a mole
of gas consisting of nonlinear polyatomic molecules can be approximated as\rf{33}:
\begin{equation}
E(T) = \underbrace{\frac{3}{2}RT}_{E^{trans}} +
\underbrace{\frac{3}{2}RT}_{E^{rot}} +
\underbrace{ZPE + E^{\prime vib}(T)}_{E^{vib}} + E^{elec}
\end{equation}
\begin{equation}
ZPE = \sum_{i=1}^{3n-6}{\frac{N h \nu_i}{2}}
\end{equation}
\begin{equation}
E^{\prime vib}(T) = \sum_{i=1}^{3n-6}{
\frac{N h \nu_i}{e^{N h \nu_i/R T} - 1} }
\end{equation}
where $n$ denotes the number of atoms in the molecule,
$ZPE$ is the {\sl zero point energy}, $E^{\prime vib}(T)$ represents the
temperature dependent portion of vibrational energy,
$\nu_i$ are the calculated vibrational frequencies and
$E^{elec}$ is the electronic energy\rf{7}.
The change in the energy occuring during protonation
of 1 mole of gas at 298.15{}$^\circ$K,
$\Delta E_{\rm Rn} = E_{{\rm AH}^+}(298.15^\circ{\rm K})
- E_{{\rm A}}(298.15^\circ{\rm K}) - E_{{\rm H}^+}(298.15^\circ{\rm K})$,
will therefore consist
of the following components:
\begin{description}
\item[$\Delta E^{rot}_{\rm Rn}$] -- the change in rotational energy upon
conversion of the reactants to the product.
Since a proton does not possess rotational kinetic energy,
this term is nonzero
only if the protonated molecule has a different fundamental shape than
the parent one
(e.g., when A is linear and AH$^+$ is nonlinear).
In our case, all parent and protonated molecules are nonlinear, hence,
$\Delta E^{rot}_{\rm Rn} \equiv 0$.
\item[$\Delta E^{trans}_{\rm Rn}$] -- the change in energy associated with
translational degrees of freedom. Since the proton on the left hand side
of the equation brings $\frac{3}{2} RT$,
the net $\Delta E^{trans}_{\rm Rn} = -\frac{3}{2} RT$; i.e.,
for 298.15 K $\Delta E^{trans}_{\rm Rn} \approx $ 0.88873249 kcal/mol.
\item[$\Delta E^{vib}_{\rm Rn}$] -- the change in energy associated with
internal vibrations of reactants and products. Only the change in the
{\sl zero point energy}, $\Delta ZPE$ is significant.
The $\Delta E^{\prime vib}$ is usually much less then 1 kcal/mol
at room temperatures (i.e., much smaller than the experimental error)
and is included here only for completeness.
\item[$\Delta E^{elec}_{\rm Rn}$] -- the change in the electronic energy
upon reaction. In our case it is the difference between ground state energies
(electronic $+$ nuclear) taken from quantum calculations with full geometry
optimization for the protonated and parent molecules.
The ground state energy of a proton is zero in this formulation.
\end{description}
Combining the above, the
following expression for proton affinity was used in actual computations:
\begin{equation}
P({\rm A}) = - \Delta E^{elec}_{\rm Rn} - \Delta ZPE - \Delta E^{\prime vib}
+ \frac{5}{2}RT
\end{equation}
The values of $E^{elec}$ and the vibrational frequences $\nu_i$ were obtained from
quantum calculations. Programs Gaussian90\rf{8} and ACES2\rf{9}
were used for traditional {\sl ab initio} calculations.
The original double-zeta basis set of Dunning/Huzinaga\rf{10}
(9s,6p/4s) $\longrightarrow$ (6111,411/31) was augmented with
polarization functions:
d-type on C, N, O with exponents
equal to: 0.8, 0.8, and 0.9 respectively, and p-type with an
exponent equal to 1.0 for H. Two sets of {\it ab initio} calculations
were performed. DH6D set used cartesian gaussians (i.e., six d-type functions)
and the DH5D\rf{11} set used five d-type functions with the angular part
represented by spherical harmonics.
These DZP basis sets had similiar characteristics
to the basis sets we used
in our DFT calculations. The RHF and MP2 calculations were also performed
for formic acid and its anion with a DH6D basis set augmented
with diffuse s-type and p-type functions (with exponent 0.0845\rf{12})
on oxygen atoms. It is well known that
adding diffuse functions on nonhydrogen atoms
dramatically improves results for proton affinities\rf{13} though
it is less important for molecular geometry.
This basis set is referred to later as DH6D${}^{(+)}$.
The {\sl ab initio} results were compared with those obtained from DFT
calculations performed with the programs: DGauss\rf{14},
DMol\rf{15} and an academic version of deMon\rf{16}.
The DGauss and deMon programs use the LSD potential developed by
Vosco {\sl et al.}\rf{17}
while DMol incorporates the Barth and Hedin LSD
potential\rf{18}.
Basis sets of double-zeta quality with polarization functions were used in
all DFT programs. In DMol, the DNP (Double Numerical with Polarization)
basis set included with the program was selected. The basis sets in
this program are given numerically
as cubic spline functions. The 300 radial points spaning a range
from nucleus to an outer distance of 10 bohrs is a default.
The angular portion of each basis function
corresponds to the appropriate spherical harmonic.
The DGauss and deMon basis functions
are analogous to those widely used in traditional {\sl ab initio}
calculations, i.e., they are represented analytically by a combination
of primitive gaussian functions; however, they have been reoptimized for the
DFT calculations\rf{20}. For DGauss and deMon DZVP atomic basis set
and A2 auxiliary fitting set was used for 1st row atoms
and DZVPP/A1 set for hydrogens. The contraction pattern for the atomic
basis set was
(9s,5p,1d/5s,1p) $\longrightarrow$ (621,41,1/41,1). In other words, for 1st
row atoms there were three s-type basis functions combining six, two, and
one gaussian primitives, respectively; two sets of p-type basis functions
containing
four and one primitives, and a single primitive for each of the six d-type
polarization functions\rf{19,20}. For hydrogens two s-type
functions combining four and one gaussian primitives, and p-type polarization
functions consisting of one gaussian were used as a basis set.
Beside atomic basis set, in DGauss and deMon, auxiliary sets of functions
are used to fit charge density and exchange-correlation potential\rf{20}.
For each atom there is one set for fitting density and another set for fitting
exchange-correlation potential. The sets consist of
uncontracted s-type gaussians and a few groups of s, p, and d-type cartesian
gaussians sharing the same exponent. The composition of these sets is denoted
as (A,B;C,D), where A,B pair specifies that A s-type gaussians, and B groups
of spd gaussians were used for fitting density, and C s-type gaussians and
D groups of spd gaussians were used to fit exchange-correlation potential.
The A2 auxiliary set was denoted (4,4;4,4), while the smaller A1 set used
for hydrogens was (3,1;3,1).
Starting geometries for molecules were obtained from model building with
Sybyl\rf{21} using previous computational results\rf{22,23} and experimental
data\rf{24-26}.
Geometries of protonated molecules were obtained by placing
a proton at position of the lone electron pair of the accepting atom.
The X--H bond lengths from the corresponding unprotonated molecules
were used for these new bonds as starting values.
These geometries were then fully optimized within given method.
It was absolutely necessary since these final geometries were
used for normal modes and vibrational frequency calculations.
Atom numbering used throughout this paper is shown in Figure 1 which
was produced with the MindTool software\rf{27}.
\\
\\
\centerline{\large\bf Figure 1}
\\
\\
Ab initio geometry optimizations were
performed in DIRECT SCF mode with Gaussian90,
and the resulting geometries were used with ACES2 for
vibrational frequency and normal mode calculations.
ACES2 computes normal modes and frequencies from the analytical
Hessian for RHF and MP2 calculations.
We also performed geometry optimizations
at the MP2 level (full core) with Gaussian90 in an analogous manner,
followed by frequency
calculations with ACES2. In the MP2 case we used the RHF-optimized geometries
as starting data.
The DFT calculations with both DGauss and DMol were carried out with
full geometry optimization at the Local Spin Density (LSD) level.
Geometry optimizations with nonlocal gradient corrections were performed with
deMon. Geometry optimization was followed by frequency calculations.
For all DFT codes the Hessian
needed for frequency and normal modes evaluation was calculated
from analytical energy gradients by
finite differences with a step of 0.01 bohr.
Integration grid in DGauss and deMon was of similar quality (MEDIUM and
FINE respectively). In DMol the FINE integration grid density was used
and the the maximum angular momentum number (LMAX) of the multipolar functions
used to analitically fit electron density and exchange-correlation potential
was 2 for hydrogens and 3 for other atoms.
Nonlocal density gradient corrections to the
LSD energies were not available in this version of DMol.
The version of DGauss used by us provided for single-point nonlocal
gradient corrections to LSD energies, while the deMon program
offered gradients of nonlocal gradient corrections,
i.e., provided for geometry optimization at this level.
Based on our previous experience\rf{28},
we chose the Becke\rf{29} functional for the exchange
and Perdew\rf{30} functional for the correlation potential,
which we label Becke-Perdew corrections.
The Becke-Perdew corrections were shown to improve results
compared to LSD approximation\rf{31}.
Since gradients of gradient-corrected energies could not be calculated
with this version of DGauss, single-point
corrected energies were calculated for the LSD-optimized geometries
and vibrational contributions to proton affinities were assumed to be the
same as for the LSD. For deMon, the fully-consistent LSD and NLSD results
are reported.
\section*{Results and Discussion}
\subsection*{Geometries}
Structural parameters obtained by different methods
are compared in Tables I--IV.
Since geometries resulting from DH6D and DH5D basis sets are practically
identical, only results of RHF and MP2 calculations for DH6D set
are reported.
It is well known that bond lengths calculated by the Hartree-Fock method are
usually too short compared to experimental values, while the corresponding
MP2 values\rf{2} are much closer to reality. It was also
supported by these calculations and clearly expressed for the polarized
C--O and O--H bonds. This effect was smaller for C--N and N--H bonds and quite
small for C--C and C--H bonds. Bond lengths from DFT bracket
{\em ab initio} values on both sides. Bonds involving
hydrogen (C--H, N--H or O--H) which
resulted from LSD calculations are the longest of all, and the effect is
most visible for the O--H bonds and the smallest for the C--H bonds.
The situation is opposite for bonds between the carbon atom and another
1st row element. These bonds calculated with LSD are the shortest of all
methods and the magnitude of this effect follows bond polarity.
Becke-Perdew corrections substantially improve bond lengths.
The H--X bonds are now shorter than for the case
of uncorrected LSD. The C--X bonds with Becke-Perdew corrections
are longer than their LSD counterparts. This brings
bond lengths calculated with Becke-Perdew corrections to a much closer
agreement with experimental and MP2 values.
Valence angles usually follow the known trend in which
the angle between shorter bonds tends to be larger, while
angle between overestimated bonds is frequently
smaller than it should be. This trend is evident if one
compares HF and MP2 results. It is also pronounced
for H--X--Y angles, which are usually too small for
LSD calculations due to the underestimation of X--Y bond
length.
Generally, the structural parameters calculated by
different methods are similar. Larger discrepancies are
visible for torsional angles in ethanol and protonated ethanol.
It is not surprising, however, since the potential energy surface
for torsional angles pivoted on the C--O bond of ethanol is
very flat, and it is still being debated if the minimum energy
of ethanol corresponds to a gauche or trans conformation.
\subsection*{Proton Affinities}
Total energies calculated by different quantum approaches are listed
in Table V.
\\
\\
\centerline{\large\bf Table V}
\\
\\
Direct comparison of absolute total energies calculated by different methods
used in this work is only possible among the results obtained by the same
program. For {\sl ab initio}, as expected,
the inclusion of six d-type polarization functions
(in the DH6D basis set), compared with five pure d-type functions (in the
DH5D basis set)
lowers the calculated energy. Using six cartesian d-type gaussians
is equivalent to using the set of five pure functions of d-type augmented with
a 3s-type (i.e., $Nr^2exp[-\alpha r^2]$) function.
This energy change is negligible for our RHF calculations
(on average 0.00060 hartree, i.e., $\approx$0.38 kcal/mol), however,
this change is roughly 10 times larger
(on average 0.0059 hartree, i.e., $\approx$3.7 kcal/mol)
for MP2 calculations. This reflects the
fact that polarization functions contribute substantially
to low-lying virtual orbitals, and the contamination with
the 3s orbital brought by combination of six cartesian d-type functions
is expressed much more strongly in correlated methods than at the
RHF level. On the other hand, the difference between corresponding MP2(5d)
and MP2(6d) energies is essentially identical for the parent
and protonated molecule
thus, this basis set modification
does not significantly affect the electronic energy contribution,
$\Delta E^{elec}_{\rm Rn}$, to the calculated proton affinity.
This should be expected since
the parent and the protonated molecule have the same number of ``heavy''
atoms whose basis sets contain d-type functions and they differ only
in the number of hydrogen atoms.
The effect of adding Becke-Perdew corrections to the LSD energy can be
examined for DGauss and deMon calculations. It is known that LSD approximation
tends to grossly overestimate correlation energy (frequently as much as 100\%)
while it underestimates electron exchange energy by as much as 10\%\rf{6}.
Both, the exchange and the correlation energy, are negative but typically
the magnitude of electron exchange energy is of several orders of magnitude
larger than the correlation energy for the same system. Therefore, the
magnitude of the 10\% correction to the exchange energy is typically
larger than the magnitude of the correlation energy. For this reason,
NLSD corrected energies are lower than the corresponding LSD energies.
DMol uses five d-type functions) while
DGauss and deMon use six cartesian gaussians as d-type functions.
Judging from {\sl ab initio} results,
this difference should not have any substantial impact
on the calculated proton affinities.
The changes in the ground state electronic energies on protonation,
$\Delta E^{elec}_{\rm Rn}$, are ordered as:
DMol $<$ DGauss(LSD) $\approx$ deMon(LSD) $<$
DGauss(NLSD) $\approx$ deMon(NLSD) $<$ MP2(6d) $\approx$ MP2(5d)
$<$ RHF(6d) $\approx$ RHF(5d).
This finding suggests that the LSD approximation
generally underestimates electronic energy changes on protonation,
and gradient (NLSD) corrections bring the results significantly closer
to the MP2 values.
The DH5D and DH6D basis sets are far from optimal for calculations
of proton affinities since they do not include diffuse s and p-type
functions\rf{13}. It is dramatically pronounced in the calculated energies
for anions. However, they were chosen here for comparison with DFT
calculations.
Inclusion of diffuse functions is absolutely required\rf{2,12,32}
for meaningful calculations of electronic energy for anions by {\em ab initio}
methods..
To illustrate the effect of diffuse functions, we performed
RHF and MP2 calculations for formate anion,
and the corresponding neutral acid with DH6D basis set
augmented with diffuse s-type and p-type functions
on oxygen atoms. As expected, this led to a substantial lowering of
electronic energy for the formate anion, while the effect was much
smaller for the neutral acid, i.e., the $\Delta E^{elec}_{\rm Rn}$ was
significantly smaller for the basis sets including diffuse functions.
Surprisingly, the DFT methods seem less sensitive to the lack of diffuse
functions. This apparent insensitivity of the DFT calculations to the
inclusion of diffuse functions may be explained by the fact that
the underlying quantity in the DFT calculation is the charge density.
Diffuse functions do not contribute substantially to the occupied
molecular orbitals, and the charge density is the sum of one-particle
densities derived from the corresponding occupied orbitals.
This reasoning is indirectly supported by our {\em ab initio}
results. The HF energy for the formate, is substantially less affected
by adding diffuse functions (0.0110 hartree) than
the MP2 energy (0.0222 hartree). The energies of formic acid
calculated with and without diffuse functions differ by 0.0032 and 0.0090,
for HF and MP2; respectively.
Table VI collects zero point energies ($ZPE$) derived from calculated
vibrational frequencies by eq 2, and the temperature-dependent portion
of vibrational enthalpies, $E^{\prime vib}$,
computed from eq 3 for the molecules studied.
All values of $E^{\prime vib}$ are small for molecules of this
size, in our case on the order of 1 kcal/mol or less. Consequently, the
differences between $E^{\prime vib}$ values for the
parent and protonated molecules are
negligibly small compared to the experimental error in proton affinity
measurements.
For larger molecules $\Delta E^{\prime vib}$ can also be safely
omitted in proton affinity calculations since the largest contributions
to $E^{\prime vib}$ come from the lowest frequency vibrations
(e.g., torsions around single bonds) and these are for the most part not
substantially affected by protonation.
\\
\\
\centerline{\large\bf Table VI}
\\
\\
The value of $ZPE$ is proportional to the sum of vibrational frequencies,
and therefore it is primarily affected by larger frequences
(e.g., bond stretching and bending vibrations).
The differences in $ZPE$ for protonated and parent molecule
is therefore not small and in our case approaches 10 kcal/mol. Protonation
results in the formation of a new covalent X---H bond and also affects
geometry and strength of vicinal bonds. For that reason, $\Delta ZPE$
must be accounted for in proton affinity calculations.
Also, systematic over- or under-estimation of frequencies will
bias the value of $ZPE$ in the same direction.
Moreover, it should be noted that $ZPE$ in our case is calculated
within the harmonic approximation. Experimental frequencies available
in the literature for our series of molecules are not
corrected for anharmonicity, and therefore,
direct comparisons with experimental values of $ZPE$ were
not attempted.
The magnitudes of calculated $ZPE$'s in our series of molecules
are ordered as:
DMol $\approx$ DGauss(LSD) $\approx$ deMon(LSD) $\approx$ deMon(NLSD)
$<$ MP2(6d) $\approx$ MP2(5d) $<$ RHF(6d) $\approx$ RHF(5d). RHF
frequences are usually larger than experimental
(even if measured frequencies are converted to harmonic ones)
and therefore the RHF calculated $ZPE$'s are too large\rf{2}.
The MP2 calculated frequencies are in much better agreement with
experiment, though in general they are also slightly higher
than measured harmonic frequencies\rf{2}. This
effect is clearly visible in stretches along polar X---H bonds.
Since LSD calculated
zero point energies are slighlty smaller than MP2 calculated ones,
they exhibit the plausible trend.
RHF and MP2 calculations with the DH6D$^{(+)}$ basis set for formic acid,
and the corresponding anion, produced $ZPE$'s very similar to
the DH6D basis set. The major contribution to error in proton affinities
of anions calculated with basis sets lacking diffuse functions is therefore
from the electronic energies.
Proton affinities were calculated from eq 4. The values of calculated
gas phase proton affinities together with the experimental values
are collected in Table VII. The errors in experimental values are roughly
2 kcal/mol for affinities within 167--204 kcal/mol range and the error
gets much larger outside this range due to the lack of adequate reference
proton affinities\rf{34}.
\\
\\
\centerline{\large\bf Table III}
\\
\\
The RHF(6d) and RHF(5d) calculated proton affinities are practically
identical. In our series of molecules they are all too large compared
to experimental values. Overestimation of proton affinities at the
RHF level is well illustrated by other authors for similar
molecules\rf{22,2,3,35,36}.
This trend in RHF calculations is mainly due to
the overestimated $\Delta E^{elec}_{\rm Rn}$, however, $\Delta ZPE$ is
usually also exaggerated since
the calculated frequency of the newly formed X---H bond
is too large. The MP2 results are practically
identical for MP2(6d) and MP2(5d) cases, and are in substantially
better agreement with experimental values. As was mentioned above, adding
diffuse functions on oxygen dramatically improves
the calculated proton affinity for the formate anion.
The proton affinities calculated within the
LSD approximation are substantially smaller than
the experimental values. For our series of molecules, LSD underestimates
$\Delta E^{elec}_{\rm Rn}$.
Proton affinities calculated with DGauss(NLSD) and deMon(NLSD)
are in excellent agreement with the experimental values.
They are even slightly better
($\sigma_{DGauss(NLSD)} = 2.4$ kcal/mol, $\sigma_{deMol(NLSD)} = 2.3$ kcal/mol)
then those from MP2 calculations
($\sigma_{MP2(6d)} = 2.6$ kcal/mol when
the result obtained with diffuse basis functions was used for formate)
and much better
than the ones from RHF calculations ($\sigma_{RHF(6d)} = 4.4$ kcal/mol).
It seems that for this particular series of molecules it does not matter
if gradient corrections were applied perturbationally
to the LSD energy (at LSD optimized
geometry) or if they were introduced in a self-consistent
way during geometry optimization. This may, however, be fortuitous due to
the fact that C---X bonds are slightly shorter and X---H and C---H
slightly longer than for the corresponding NLSD calculations.
It may provide for some cancellation of errors which results in similar
overall values of affinities.
\section*{Concluding Remarks}
We have shown that the DFT methodology is a good candidate for routine
calculations of proton affinities. However, it is evident that the
nonlocal gradient corrections have to be used to correctly estimate the
change in electronic energy of protonation. The results obtained with
Becke-Perdew corrections applied perturbationally or in a self-consistent
manner are of MP2 quality but require much less computation. Moreover,
since DFT methods scale formally with molecular size as $N^3$ (compared
to $N^5$ for MP2 {\em ab initio\/} approach) they can be applied to much larger
molecules.
\vskip 0.3in.
\noindent
{\bf Acknowledgment:} Authors are grateful for the Ohio Supercomputer Center
grant of computer time which made these calculations possible. JKL, RAH
and DDM appreciate the support of National Institutes of General Medical
Sciences, GM 29358.
\section*{References}
\begin{list}{}{\setlength{\leftmargin}{1em} \setlength{\labelwidth}{1em}
\setlength{\itemindent}{-1em} \setlength{\itemsep}{0in}}
\item[(1) ]
Dixon D.A.; Lias, S.G. In {\em Molecular Structure and Energetics, Vol. 2,
Physical Measurements\/}; Liebman, J. F.; Greenberg, A., Eds.,
VCH, Deerfield Beach, FL, 1987.
\item[(2) ] Hehre, W.J.; Radom, L.; Schleyer, P.v.R.; Pople, J.A.
{\em Ab Initio
Molecular Orbital Theory\/}; John Willey and Sons, New York, 1986.
\item[(3) ] Ozment, J. L.; Schmiedekamp, A. M. {\em Int. J. Quantum Chem.}
{\bf 1992}, {\em 43} 783.
\item[(4) ]
Parr, R.G.; Yang W. {\em Density Functional Theory of Atoms and Molecules\/};
Oxford University, New York, 1989.
\item[(5) ]
{\em Density Functional Methods in Chemistry}; Labanowski, J. K.;
Andzelm, J.W., Eds., Springer Verlag, New York, 1991.
\item[(6) ] Ziegler, T. {\em Chem. Rev.\/} {\bf 1991}, {\em 91\/}, 651.
\item[(7) ] we used the following constants:
$N = 6.0221367 \times 10^{23}$ molecules/mol,
1 J $\equiv$ 4.184 cal; R = 1.9872156 cal mol$^{-1}$ K$^{-1}$;
h = 6.6260755$\times 10^{-34}$ J s;
1 cm$^{-1}$ = 0.0028591439 kcal mol$^{-1}$;
1 hartree = 627.50955 kcal mol$^{-1}$;
1 bohr = 0.529177249 \AA .
\item[(8) ] Gaussian90, Revision H; Frisch, M. J.; Head-Gordon, M.;
Trucks, G. W.; Foresman, J. B.; Schlegel, H. B.;
Raghavachari, K.; Robb, M.; Binkley, J. S.; Gonzalez, C.;
Defrees, D. J.; Fox, D. J.; Whiteside, R. A.; Seeger, R.;
Melius, C. F.; Baker, J.; Martin, R. L.; Kahn, L. R.;
Stewart, J. J. P.; Topiol, S.; Pople, J. A.
Gaussian Inc., Pittsburgh PA, 1990.
\item[(9) ] Stanton, J. F.; Gauss, J.; Watts, J. D.; Lauderdale, W. J.;
Bartlett, R. J. {\em ACES2 quantum package}, Quantum Theory Project of the
University of Florida at Gainesville.
\item[(10) ] Dunning, T.H.; Hay, P.J. In {\em Modern
Theoretical Chemistry vol. 3\/}; Schaefer III, H.F., Ed.,
Plenum Press, New York, 1977; pp 1--28.
\item[(11) ] These basis sets differ
from the D95H** set supplied with Gaussian90 in the values of exponents at
polarization functions.
\item[(12) ] Clark, T.; Chandrasekhar, J.; Spitznagel, G.W.,
Schleyer P.v.R. {\em J. Comp. Chem.\/} {\bf 1983}, {\em 4\/}, 294.
\item[(13) ] Del Bene, J. E. {\em J. Phys. Chem.\/} {\bf 1993}, {\em 97}, 107.
\item[(14) ] Andzelm, J. In {\em Density Functional Methods
in Chemistry\/}; Labanowski, J.K.; Andzelm, J.W., Eds.,
Springer-Verlag, New York, 1991; pp 155--174.
\item[(15) ] DMol: Local Density Functional Program, Ver. 2.0; Biosym
Technologies, San Diego, CA.
\item[(16) ] St-Amant, A., Ph.D. thesis, Universit\'e de Montr\`eal, 1992;
St-Amant, A.; Salahub, D. R. {\em Chem. Phys. Lett.}
{\bf 1990}, {\em 169}, 387.
\item[(17) ] Vosco, S.H.; Wilk, L., Nuisar, M. {\em Can. J. Phys.} {\bf 1980},
{\em 58\/}, 1200.
\item[(18) ] von Barth, U.; Hedin, L. {\em J. Phys. C\/} {\bf 1972},
{\em 5\/}, 1629.
\item[(19) ] Huzinaga, S.; Andzelm, J.; Klobukowski, M.;
Radzio, E.; Sakai, Y.; Tatewaki, H. {\em Gaussian Basis Sets for Molecular
Calculations\/}; Elsevier, Amsterdam, 1984.
\item[(20) ] Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E.;
{\em Can. J. Chem.} {\bf 1992}, {\em 70}, 560.
\item[(21) ] Sybyl Molecular Modeling Software, Ver. 5.3, Tripos Associates,
St. Louis, MO.
\item[(22) ] von Nagy, E.I.; Kimura, K.
{\em J. Phys. Chem.\/} {\bf 1990}, {\em 94\/}, 8041.
\item[(23) ] Marcoccia, J.F.; Csizmadia, I.G.; Peterson, M.R.; Poirier,
R.A. {\em Gazz. Chim. Ital.\/} {\bf 1990}, {\em 120\/}, 77.
\item[(24) ] Harmony, M. D.; Laurie, V. W.; Kuczkowski, R. L.;
Schwendeman, R. H.; Ramsay, D. A.; Lovas, F. J.; Lafferty, W. J.;
Maki, A. G.; {\em J. Phys. Chem. Ref. Data} {\bf 1979}, {\em 8}, 619.
\item[(25) ] Callomon, J. H.; Hirota, E.; Kuchitsu, K.;
Lafferty, W. J.; Maki, A. G.; Pote, C. S. {\em Structure Data of Free
Polyatomic Molecules}; Hellwege, K.-H.; Hellwege, A. M., Eds.;
Landolt-Bornstein, Group II, New Series, Vol. 7, 1976.
\item[(26) ] Sasada, Y.; Takano, M.; Satoh, T.
{\em J. Molec. Spectr.\/} {\bf 1971}, {\em 38\/}, 33;
Sasada, Y. J. {\em Mol. Spectr.} {\bf 1988}, {\em 190}, 93.
\item[(27) ] Blake, J.; Tirado-Rives, J.; {\em MindTool for SunOS}, ver.
March 92. Available from anonymous ftp on rani.chem.yale.edu and
www.ccl.net.
\item[(28) ] Hill, R.A.; Labanowski, J.K.; Heisterberg, D.J.;
Miller, D.D. In {\em Density Functional Methods
in Chemistry\/}; Labanowski, J.K.; Andzelm, J.W., Eds.,
Springer-Verlag, New York, 1991; pp 357--372.
\item[(29) ] Becke, A.D. {\em Phys. Rev. A\/} {\bf 1988}, {\em 38\/}, 3098.
\item[(30) ] Perdew, J.P. {\em Phys. Rev. B\/} {\bf 1986}, {\em 33\/}, 8822;
erratum {\em ibid.\/} {\bf 1986}, {\em 38\/}, 7406.
\item[(31) ] Sim, F.; St-Amant, A.; Papai, I.; Salahub, D.R.
{\em J. Am. Chem. Soc.\/} {\bf 1992}, {\em 114\/}, 4391.
\item[(32) ] Masamura, M. {\em Theor. Chim. Acta\/}
{\bf 1989}, {\em 75\/}, 433.
\item[(33) ] McQuarrie, D. A.; {\em Statistical Mechanics}
Harper \& Row, New York, 1976.
\item[(34) ] Lias, S.G.; Liebman, J.F.; Levin, R.D. {\em J. Chem. Phys.
Ref. Data} {\bf 1984}, {\em 13\/}, 695; and updates available from the
NIST Web site: {\tt http://webbook.nist.gov/}.
\item[(35) ] Del Bene, J.E.; Frish, M.J.; Raghavachari, K.; Pople, J.A.
{\em J. Phys. Chem.} {\bf 1982}, {\em 86\/}, 1529.
\item[(36) ] Siggel, M.R.F.; Thomas, T.D.; Saethre, L.J.
{\em J. Am. Chem. Soc.\/} {\bf 1988}, {\em 110\/}, 91.
\item[(37) ] Callomon, J. H.; Hirota, E.; Iijima, T.; Kuchitsu, K.;
Lafferty, C. S.; Pote, C. S. {\em Structure Data of Free
Polyatomic Molecules}; Hellwege, K.-H.; Hellwege, A. M., Eds.;
Landolt-Bornstein, Group II, New Series, Vol. 15, 1986.
\item[(38) ] Benston O.J., Ewbank J.D., Paul D.W., Klimkowskim V.J., Schaefer L.
{\em Appl. Spectr.} {\bf 1984}, {\em 38}, 204.
\item[(39) ] D. Davies R. W., Robiette A.G., Gerry M.C.L., Bjarnov E., Winnewisser G.,
{\em J. Mol. Spectr.} {\bf 1980}, {\em 81}, 93.
\item[(40) ] Salahub, D.R.; Fournier, R.; M\l{}ynarski, P,;
Papai, I.; St-Amant, A.; Ushio, J. In {\em Density Functional Methods
in Chemistry\/}; Labanowski, J.K.; Andzelm, J.W., Eds.,
Springer-Verlag, New York, 1991; pp 77--100.
\item[(41) ] Bartmes, J.E.; McIver Jr., R.T. In {\sl Gas Phase Ion Chemistry,
vol 2\/}; Bowers, M.T., Ed., Academic Press, New York, 1979; pp 87--121.
\item[(42) ] Oie, T,; Topol, I.A.; Burt, S.K.; {\em J. Phys. Chem.}
{\bf 1994}, {\em 98}, 1121.
\item[(43) ] Chen, H.; Krasowski, M.; Fitzgerald G.;
{\em J. Chem. Phys.} {\bf 1993}, {\em 98}, 8710.
\end{list}
\end{document}