Ab Initio Benchmark Calculations on Ca(II ... - ACS Publications

ARTICLE pubs.acs.org/JPCA

Ab Initio Benchmark Calculations on Ca(II) Complexes and Assessment of Density Functional Theory Methodologies Dimas Su
arez,† Víctor M. Ray
on,‡ Natalia Díaz,† and Hayd
ee Vald
es*,†,§ † ‡

Departamento de Química Física y Analítica, Facultad de Química Universidad de Oviedo, 33007 Oviedo, Spain Departamento de Química Física y Química Inorg
anica, Facultad de Ciencias Universidad de Valladolid, 47005 Valladolid, Spain

bS Supporting Information ABSTRACT: A set of benchmark results for the geometries, binding energies, and protonation affinities of 24 complexes of small organic ligands with Ca(II) is provided. The chosen level of theory is CCSD(T)/CBS obtained by means of a composite procedure. The performance of four density functionals, namely, PW91, PBE, B3LYP, and TPSS and several Pople-type basis sets, namely, 6-31G(d), 6-31+G(d), 6-31+G(2d,p) and 6-311+G(d) have been assessed. Additionally, the nature of the metal ligand bonding has been analyzed by means of the Symmetry Adapted Perturbation Theory (SAPT). We have found that the B3LYP hybrid functional, in conjunction with either the polarized double-ζ 6-31+G(2d,p) or the triple-ζ 6-311+ G(d) basis sets, yields the closest results compared to the benchmark data. The SAPT analysis stresses the importance of induction effects in the binding of these complexes and suggests that consideration of classical electrostatic contributions alone may not be reliable enough for the prediction of relative binding energies for Ca(II) complexes.

’ INTRODUCTION Calcium ion is one of the metals most commonly found in biological systems, playing important roles in diverse biochemical processes.1 Inside eukaryotic cells, Ca(II) concentration levels regulate muscle contraction, secretion, glycolysis, gluconeogenesis and ion transport.2 It also participates as a signaling mediator in cell cycle regulation, proliferation and apoptosis.3 Outside cells, Ca(II) plays crucial roles in many matrix
matrix, cellmatrix and cell
cell contacts, in the mineralization of bones, in blood clotting, and in receptor structures.4 In all these cases, binding of calcium ions to proteins seems to affect their structural stability or to modulate substrate binding and the catalytic activity of enzymes. According to the analysis of the Cambridge Structural Database (CSD), the calcium ion is able to interact with six, seven or eight ligands. The metal
ligand bonding is commonly considered to be predominantly ionic, so the ligand disposition is mainly determined by the possibilities of packing donor atoms around the Ca(II) cation.5
7 The calcium complexes show a strong preference to bind to oxygen atoms as the donors, whereas complexes with nitrogen or sulfur donors are very rare.8 In addition, carboxylate complexes with Ca(II) show several types of geometrical behavior (syn-monodentate, antimonodentate, and symmetrical bidentate) showing a greater tendency toward bidentate carboxylate binding. Interestingly, it has been shown that the geometry found in all the Ca(II)-containing protein structures in the Protein Data Bank determined at resolution e1.6 Å, agrees with that predicted from the CSD for ligands that are analogues of amino acid side chains in proteins.9 Binding of Ca(II) to proteins is thus a topic of great interest for which deeper understanding requires the proper study of the r 2011 American Chemical Society

calcium ligand abilities. A plausible way to achieve such knowledge is by means of computational chemistry. Then, following the same philosophy than in our previous paper about the Zn(II) ligand abilities,10 we have selected a set of complexes modeling the Ca(II)
ligand interactions and proceeded similarly. In other words, we have first performed high-level quantum chemical calculations on the selected complexes in order to obtain a reliable energetic, structural, and electronic description of the Ca(II)
ligand binding abilities. Second, we have used the data obtained by means of the high-level quantum chemical calculations as benchmark data to assess the performance of various density functional theory (DFT) functionals. And third, we have analyzed the nature of the ligand
metal bond by means of an energy decomposition method aiming to obtain relevant information on the relative importance of the different electrostatic, induction and dispersion contributions to the interaction energy. The assessment of the performance of the DFT functionals deserves a special comment. Nowadays, DFT has become the most popular correlated methodology for the study of biomolecules due to its fast computer performance and its capability for describing both noncovalent interactions11,12 and chemical reactions13 as well as for providing molecular properties and descriptors.14 However, it should not be forgotten that the versatility of a particular functional has to be carefully tested in advance15 because many functionals have not been parametrized to fit data Special Issue: Pavel Hobza Festschrift Received: May 31, 2011 Revised: August 11, 2011 Published: September 28, 2011 11331

dx.doi.org/10.1021/jp205101z | J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A for inorganic or organometallic molecules and their applicability outside the scope for which they were developed is not always guaranteed. Hence, assessing the performance of the DFT functionals against ab initio benchmark data is always recommendable. Ab initio calculations are time-consuming, and consequently, the size of the system to be studied is considerably limited especially when the accuracy of coupled-cluster (CC) methods combined with infinite (extrapolated to the complete basis set limit, CBS) basis sets is required. For this reason, the performance of the DFT functionals is typically assessed against ab initio benchmark data on small model systems (for Ca(II) complexes, see for example, refs 16
22). According to this philosophy and as commented before, various DFT functionals are here assessed against ab initio benchmark calculations for a conveniently chosen set of Ca(II)
ligand complexes mimicking interactions that potentially exist in proteins. In particular, the following ligands were chosen: water, methanol, formic acid, acetic acid, formaldehyde, acetone, formamide, acetamide, Nmethyl acetamide (NMA), methyl acetate, ammonia, methylamine, methanimine, 1H-imidazole, benzene, hydrogen sulfide, methanethiol, hydroxide, methanolate, acetate, imidazolate, hydrosulfide, methanethiolate, formate, and methyl phosphate. Note that the selected set of ligands is mainly populated by molecules with O donor atoms in different chemical environments (sp2 O, sp3 O, neutral, ionic, ...) as calcium ions in biological systems show a clear preference for binding to carboxylate/ carbonyl groups and water molecules. However, we have also considered other ligands with N and S donor as well as benzene because it is important to properly describe the whole spectrum of possible Ca(II)
ligand interactions. The third goal of this study is the analysis of the nature of the metal
ligand bond in these model complexes. Polarization interactions (also referred to as induction) are of utmost importance in the properties and dynamics of many biologically relevant systems.23 However, translation of these polarization effects into a particular model function (force field) to be efficiently used in atomistic molecular simulations is quite challenging and a matter of current debate.24
29 Most nonpolarizable force fields describe nonbonded interactions by means of a pairwise additive function that typically includes an electrostatic and a Lennard
Jones contribution. Regarding the electrostatic term, each atom is assigned an effective partial (fixed) charge that is adjusted to model, in an average way, polarization effects in the condensed phase (we note that because polarization interactions are highly nonadditive they cannot be properly described by a pairwise additive function). We would like to point out that this adjustment is not performed for metal ions because most standard force fields retain their original integer (bare) charges [i.e., +2 for Ca(II)]. The Lennard
Jones parameters are even more empirical than electrostatic charges because they are frequently used to compensate for inaccuracies in the force field.30 This machinery yields reasonable results whenever induction effects are not dominant in the system. For those cases where polarization effects cannot be neglected polarizable force fields have to be used instead. As commented above, derivation of specific functions that include induction effects and that are suitable for biological simulations is a hot topic of current research. A second important point is to assess the magnitude and importance of polarization effects for a particular interaction, for example the interaction between a specific metal ion and nucleophilic ligands. To this end, quantum mechanical (QM) calculations have proven very useful and, in particular, techniques aimed to partition the interaction energy

ARTICLE

provide a rigorous way to estimate the importance of individual interaction energy terms (electrostatic, induction, dispersion, and exchange-repulsion).10,29,31
37 In this work, we will assess the importance of polarization effects in Ca(II) complexes by means of the Symmetry Adapted Perturbation Theory (SAPT). We hope that these results, in conjunction with those provided by other research groups and partitioning schemes, will contribute to the improvement and refinement of polarizable molecular force fields. Ultimately, we hope that all this information may gradually help to settle the physical bases governing metal binding in proteins.

’ COMPUTATIONAL METHODS Ab Initio Calculations. Ab initio benchmark calculations of the monoligand Ca(II) complexes studied in this work were carried out assuming the frozen core approximation in all the correlated calculations treating the Ca(II) ions likewise third period elements (i.e., core orbitals for Ca(II) ions include only 1s, 2s, and 2p). We systematically used the Dunning’s basis set augmented with diffuse functions [aug-cc-pVnZ (n = 3
5)]38,39 for all the atoms except for the Ca(II) ion for which we used the cc-pVnZ basis sets developed by Koput and Peterson to describe both its 4s valence and 3s3p core spaces.40 Nevertheless and for the sake of simplicity, herein we will only use the notation aug-ccpVnZ [n = 3
5)]. The basis sets for calcium were taken from the EMSL Basis Set Exchange.41,42 Molecular geometries and harmonic frequencies were carried out at the MP2/aug-cc-pVTZ level of theory and electronic energies were refined by performing CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVTZ single-point energy calculations (coupled cluster single and double excitation augmented with a noniterative treatment of triple excitations).43 Subsequently, electronic energies were extrapolated toward the CBS (complete basis set) limit from MP2/aug-cc-pVnZ (n = 4, 5)//MP2/aug-cc-pVTZ energies using two different extrapolation formulas:44

En ¼ ECBS þ An
3

ð1Þ 2

En ¼ ECBS þ Ann
1 þ Be
ðn
1Þ

ð2Þ

where n is the cardinal number of the basis set and ECBS, A, and B are fitting parameters, with ECBS being the resulting estimate of the CBS limit. The average ECBS value obtained from these two expressions has been reported as a conservative estimate of the actual CBS limit.44 Herein, extrapolations using eqs 1 and 2 with the aug-cc-pVnZ (n = 4, 5) basis sets were used systematically on the MP2 correlation energies. The HF energies were not extrapolated, and the 5Z HF values were simply taken as the most accurate estimates of the HF limits. Although accurate, the calculation of the CCSD(T) CBS limit here employed is in general impractical and, thus, we decided to approximate the corresponding CCSD(T) CBS values by means of the following “composite” formula: ECCSDðTÞ CBS ≈ECCSDðTÞ=aug-cc-pVTZ þ ðEMP2=CBS
EMP2=aug-cc-pVTZ Þ

ð3Þ Additionally, the sensitivity of the computed binding energies and molecular properties with respect to basis set effects and/or the method employed for geometry optimization were assessed by carrying out further ab initio calculations for some of the Ca(II)
L complexes (L = OH
and H2O). For these small 11332

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A complexes, molecular geometries were reoptimized at the CCSD/aug-cc-pVTZ level. Besides, the CBS limit of the CCSD(T) energies was also estimated by means of single-point CCSD(T)/aug-cc-pVnZ (n = 4, 5) calculations using eqs 1 and 2. For the sake of completeness, CCSD(T)/aug-cc-pVTZ//CCSD/ aug-cc-pVTZ energies were also calculated. By means of the test calculations on the [Ca(H2O)]2+ and [Ca(OH)]+ complexes, we also examined whether the same CBS limit is consistently attained from both electronic energies that included the counterpoise (CP) correction45 for Basis Set Superposition Error (BSSE) and uncorrected energies. Finally, the validity of the frozen core approximation was assessed by carrying out full electron correlation calculations on the Ca2+ water/hydroxide systems. In this case, molecular geometries were optimized at the CCSD/aug-cc-pVTZ level and followed by single-point CCSD(T)/aug-cc-pVnZ (n = 3
5). It must be noticed that in the full electron calculations, we used the following core
valence basis sets for the Ca2+ ions: CVT+2dZ, CVQ+2dZ, and CV5Z, which have been developed and recommended by Martin and coworkers46 for including deep core correlation effects in Ca. Density Functional Calculations. In this work, we aim to assess the performance of four of the most commonly used and widely available in quantum chemistry program packages “generalpurpose” functionals. Specifically, we selected two functionals based on the generalized gradient approximation (GGA), PW9147 and PBE,48 one hybrid functional, B3LYP,49
52 and one metaGGA, TPSS53 which, unlike the formers, is a nonempirical functional. In conjunction with these four functionals, we tested several basis set of double-ζ and triple-ζ quality also equally accessible. Consequently, we selected a “6-31G” double-ζ and a “6-311G” triple-ζ basis sets. For the main group elements these basis sets consists of the Pople’s 6-31G54 and 6-311G55 basis. The only exception is sulfur for which the 6-311G set from McLean and Chandler56 was used. Currently, there are two double-ζ sets, usually referred to as “6-31G”, for calcium. The first one was initially developed by Rassolov et al. in 199857 and later improved by the same authors in 2001, including the 3d shell in the valence space.58 Its contraction [5s,4p,2d] or [66631,6631,31] has the same structure as the 6-31G set for the first two rows of the periodic table. The second available “6-31G” double-ζ basis does not actually have the expected “6-31G” structure because its contraction is [5s,3p,1d] or [63311,531,3] and has been developed as part of the Gaussian-2 (G2) theory in 1997.59 We would like to point out that in the Gaussian03 program package the default 6-31G basis is the old one from Rassolov et al. (year 1998), which does not include the 3d shell in the valence space. We carried out a small set of test calculations to assess the performance of these two ‘6-31G’ basis sets. The results, available as Supporting Information (see Tables S1 and S2), show that both basis sets yield similar results. We have selected the 6-31G set from G2 theory for the present study although the conclusions obtained should be valid for the improved Rassolov’s 6-31G basis set (year 2001) as well. In the following, we will refer to the “G2” 6-31G basis simply as 6-31G. Regarding the triple-ζ basis for the calcium atom, we selected the [8s,7p,1d] or [62111111, 3311111, 3] basis set available in G03 as the default 6-311G basis for this element.59 Similarly, this set has been developed as part of the Gaussian-2 (G2) theory. Both the 6-31G and the 6-311G basis sets were augmented with sp diffuse and d polarization functions. The following combinations have been tested in the current study: 6-31G(d), 6-31+G(d), 6-31+G(2d,p), and 6-311+G(d). The exponent of

ARTICLE

the diffuse primitive for calcium (0.0071) as well as the exponents of the d polarization functions
0.216 (6-31G) and 0.260 (6-311G) have been taken from ref 59. Cartesian functions were used with the 6-31G set, whereas spherical functions were used with the 6-311G. For the sake of completeness, geometries and binding energies were also computed by means of the MP2 method combined with these basis sets. The geometries of the tested Ca(II) complexes were fully optimized with the above-described functionals and basis sets. Analytical frequency calculations were carried out to check that the optimized structures were minima on the corresponding potential energy surfaces. The default integration grid in G03 (75 radial shells and 302 angular points) was used for all the DFT calculations. For test purposes, we also carried out calculations on the whole set of complexes using a larger (“Ultrafine”) integration grid (99 radial shells and 590 angular points). These calculations were performed using the B3LYP functional in conjunction with the aug-cc-pVTZ basis set. All of the wave function and density functional calculations carried out in this work were performed with the Gaussian0360 program package. SAPT Analyses. The nature of the Ca(II)
ligand interaction was analyzed by means of the Density Functional Theory
Symmetry Adapted Perturbation Theory (DFT-SAPT) developed by Hesselmann and Jansen61
63 and independently by Szalewicz and co-workers64,65 [in the latter case, the theory is usually referred to as SAPT(DFT)]. DFT-SAPT interaction energy is expressed as a sum of physically meaningful contributions, namely, electrostatic, induction, dispersion, and their corresponding exchange counterparts ð1Þ

ð1Þ

ð2Þ

ð2Þ

ð2Þ

ð2Þ

EDFT-SAPT ¼ Eelst þ Eexch þ Eind þ Eexch-ind þ Edisp þ Eexch-disp int

ð4Þ where the superscript refers to the order in which these terms appear in the perturbation expansion. The molecular orbitals and orbital energies resulting from a DFT calculation on the fragments are used to compute E(1) elst whereas Frequency Dependent (2) Dipole Polarizabilities (FDDS) are used to obtain E(2) ind and Edisp. These three contributions are potentially exact, meaning that they would be exact if the exact exchange-correlation (xc) potential was used. On the other hand, the exchange terms, which cannot be expressed in terms of fragments’ properties, are computed using the first- and second-order density matrices, which are only approximations to the true matrices. Thus, exchange contributions are not potentially exact. In any case, it has been shown that DFT-SAPT provides indeed exchange terms, which are in good agreement with those computed using SAPT. In this work we will group the DFT-SAPT contributions in terms arising at first- and second-order: ð1Þ

ð1Þ

ð2Þ

ð2Þ

Eð1Þ ¼ Eelst þ Eexch

ð5Þ ð2Þ

ð2Þ

Eð2Þ ¼ Eind þ Eexch-ind þ Edisp þ Eexch-disp

ð6Þ

In this way, we split the interaction energy in contributions associated to the interaction of frozen (unrelaxed) wave functions of the monomers (E(1)) and those contributions associated to the relaxation (through polarization and charge transfer) of fragments’ charge densities (E(2)). We finally point out that thirdand higher-order induction and exchange induction effects can be included at the HF (noncorrelated) level computing the 11333

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Geometrical arrangement of the studied Ca(II) complexes with neutral ligands along with the most important geometrical parameters, as computed at the MP2/aug-cc-pVTZ level of theory. For water, froze-core CCSD/aug-cc-pVTZ and full electron CCSD/aug-cc-pVTZ data are given in italics and in squared brackets, respectively. All distances are given in Å.

difference between the SAPT noncorrelated first- and secondorder energy and the supermolecular HF interaction energy: ð10Þ

ð10Þ

ð20Þ

ð20Þ

HF δEHF int ¼ Eint
Eelst
Eexch
Eind, r
Eexch-ind, r

ð7Þ

The final DFT-SAPT interaction energy reads, therefore EDFT-SAPT ¼ Eð1Þ þ Eð2Þ þ δEHF int int

ð8Þ

Further details of the DFT-SAPT [and the equivalent SAPT(DFT)] method can be found in the original references. All the DFT-SAPT calculations were performed on the MP2/aug-ccpVTZ geometries using an aug-cc-pVTZ dimer centered basis set (DCBS). The PBE0 xc potential was used throughout. In order to correct the wrong asymptotic behavior of the PBE exchange potential the correction approach of Gr€uning et al. was used.66 This approach requires the calculation for each fragment of the difference between the (exact) vertical ionization potential (IP) and the negative of the highest occupied molecular orbital (HOMO). The vertical IPs were taken from the NIST Chemistry webbook,67 whereas the HOMO energies were computed at the PBE0/aug-cc-pVTZ level. All the DFT-SAPT calculations were carried out with the MOLPRO suite of programs.68

’ RESULTS AND DISCUSSION Ab Initio Calculations. The optimized structures for the monoligand complexes formed between Ca(II) and the neutral ligands (L = water, methanol, formic acid, acetic acid,

formaldehyde, acetone, formamide, acetamide, N-methyl acetamide (NMA), methyl acetate, ammonia, methylamine, methanimine, 1H-imidazole, benzene, hydrogen sulfide, and methanethiol) and anionic ligands (L = hydroxide, methanolate, acetate, imidazolate, hydrosulfide, methanethiolate, formate, and methyl phosphate), respectively, are collected in Figures 1 and 2. For each complex shown in Figures 1 and 2, the equilibrium Ca
L distance involving the Ca atom and the donor atom(s) of the ligand, the Natural Population Analysis (NPA) charge of the Ca atom, and the harmonic frequency of the stretching normal mode associated with the Ca
L bond are also reported. As mentioned in Computational Methods, the consistency of the computational scheme employed for computing the binding energies was tested in the case of the complexes of Ca(II) with water and hydroxide by using more elaborated computational schemes as well (see Table 1). We will analyze first these validation calculations, and then we will present the benchmark data for all the Ca(II) complexes. Validation Calculations on the Water and Hydroxide Complexes. Test calculations performed on the Ca(II)
L (L = water and hydroxide) complexes that assess the adequacy of the various choices made in the computational scheme defined in eqs 1
3 (i.e., “composite” approximation, MP2 geometries, valence correlation, CBS extrapolation, ...) are collected in Table 1. Such calculations have been done aiming to ensure that the benchmark calculations here performed are accurate and could replace missing experimental data. 11334

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A

ARTICLE

Figure 2. Geometrical arrangement of the studied Ca(II) complexes with anionic ligands along with the most important geometrical parameters as computed at the MP2/aug-cc-pVTZ level of theory. For hydroxide, froze-core CCSD/aug-cc-pVTZ and full electron CCSD/aug-cc-pVTZ data are given in italics and in squared brackets, respectively. All distances are given in Å.

Table 1. Binding Energies (in kcal/mol) for the Ca(II) Complexes with Water and Hydroxide Obtained with Different Ab Initio Computational Protocols ligand

MP2/aug-cc-pVTZ

water hydroxide


57.14
338.36

water hydroxide


54.75
334.49

CCSD(T)/aug-cc-pVTZ

MP2/CBSa

composite methodb


57.49
341.81


57.33
342.68

Frozen Core, MP2/aug-cc-pVTZ Geometries
56.97
339.22

Frozen Core, MP2/aug-cc-pVTZ Geometries, CP-Corrected Energies
54.32
334.83


56.60
340.40


56.18
340.74

CCSD(T)/CBSa

ligand

CCSD/aug-cc-pVTZ

CCSD(T)/aug-cc-pVTZ

water hydroxide


56.91
339.29

water hydroxide

Full Electron, CCSD/aug-cc-pVTZ Geometries, Martin’s Basis Sets for Ca
58.55
55.80
339.27
339.55

Frozen Core, CCSD/aug-cc-pVTZ Geometries
56.98
339.23


57.23
342.83
56.83
343.10

Full Electron, CCSD/aug-cc-pVTZ Geometries, Martin’s Basis Sets for Ca, CP-Corrected Energies water hydroxide


57.94
337.86


55.16
338.03


56.77
342.91

a

Obtained from CBS extrapolation of the correlation energy based on eqs 1 and 2 and using the HF/aug-cc-pV5Z energies. b Using an additive combination of electronic energies.

A first conclusion that can be drawn from the analysis of the data collected in Table 1 is that binding energies (De) obtained

using the MP2-optimized geometries and the Dunning et al. basis sets combined with the “composite method” (
57.33 and 11335

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A

ARTICLE

Table 2. Binding Energies (in kcal/mol) for the Monoligand Ca(II) Complexes Obtained with the MP2 and CCSD(T) Methodsa ligand

MP2/aug-cc-pVTZ De

MP2/CBSb De

CCSD(T)/ aug-cc-pVTZ De

composite methodc De

composite methodd D0

water


57.14


56.97


57.49


57.33


55.55

methanol


67.79


67.58


67.63


67.42


66.29

formic acid


81.19


81.62


81.01


81.44


79.88

acetic acid


83.91


84.29


83.49


83.87


83.38

formaldehyde


67.22


67.58


66.96


67.32


65.74

acetone


89.11


89.45


88.79


89.13


88.25

formamide


96.01


96.58


96.06


96.64


94.36

acetamide N-methyl acetamide


103.73
109.38


104.29
110.00


103.67
109.11


104.23
109.73


102.12
107.86

methyl acetate


89.64


90.29


89.30


89.95


88.95

ammonia


66.00


65.86


65.86


65.72


63.56

methylamine


73.43


73.32


72.82


72.71


71.33

methanimine


73.65


73.43


72.43


72.21


70.41

1H-imidazol


97.98


98.05


95.62


95.69


94.47

benzene


86.88


86.02


86.22


85.36


84.03

hydrogen sulfide methanetiol


45.90
57.82


45.61
57.56


46.67
58.41


46.38
58.15


44.96
57.31
340.66

hydroxide


338.36


339.22


341.81


342.68

metoxi


337.03


336.45


339.87


339.29


335.56

formate


315.35


316.13


317.55


318.33


316.18

acetate


323.67


324.20


325.49


326.02


324.15

imidazolate


279.65


279.86


277.91


278.12


277.04

hydrosulfide


285.91


285.68


290.24


290.01


288.88

methanetiolate


292.57


292.05


295.88


295.37


293.73

a

Molecular geometries were optimized at the MP2/aug-cc-pVTZ level. b Obtained from CBS extrapolation of the MP2 correlation energy based on eqs 1 and 2 and using the HF/aug-cc-pV5Z energies. c Using an additive combination of electronic energies (CCSD(T)/aug-cc-pVTZ + MP2/CBS
MP2/aug-cc-pVTZ). d Using an additive combination of electronic energies and including MP2/aug-cc-pVTZ ZPVE energies.


342.68 for L = H2O and OH
, respectively) are nearly identical to those obtained at the CCSD(T)/CBS//CCSD/aug-cc-pVTZ level of theory (
57.23 and
342.83 kcal/mol), showing, thus, that it seems reasonable to rely on the “composite method” data when the CCSD(T)/CBS calculations are not feasible. The relevance of deep core
valence correlation effects was evaluated by carrying out full electron CCSD calculations on the Ca(II)
L (L = water and hydroxide) complexes and using the aug-cc-pVnZ basis sets for the nonmetal atoms with the larger basis sets for the Ca(II) ion developed by Martin and co-workers,46 which contain about 30% more Gaussian type functions than the cc-pVnZ basis sets of Koput and Peterson40 in order to take into account both outer-core valence and deep-core-valence effects in calcium. We found that the CCSD(T)/CBS binding energies computed within the frozen core approximation (
57.23/
342.83 kcal/ mol for H2O/OH
) deviate little from the corresponding De using the all-electrons basis set, namely,
0.4 and 0.3 kcal/mol for H2O and OH
, respectively. Thus, it may be reasonable to estimate that deep core
valence correlation effects may have an impact lower than 1.0 kcal/mol in absolute value for the monoligand Ca2+ complexes. We have additionally extrapolated the CP corrected binding energies to the basis set limit. The “composite method” De energies obtained from CP-corrected electronic energies (
56.18 and
340.74 kcal/mol) are about 1
2 kcal/mol above the CPuncorrected values (
57.33 and
342.68 kcal/mol). Notice, however, that both the CP-corrected and CP-uncorrected energies using the Martin et al. basis sets are very similar to each

other (
56.77/
342.91 and
56.83/
343.10 kcal/mol, respectively). Most importantly, we note that the “composite method” energies using the CP-uncorrected data are closer to the Martin et al. results than the “composite” CP-corrected values. The difference, as commented above, is smaller than 0.5 kcal/mol. As mentioned in Computational Methods, we also assessed the reliability of the CBS extrapolation scheme by computing the binding energy change when going from MP2/aug-cc-pVQZ to MP2/aug-cc-pV5Z and from MP2/aug-cc-pV5Z to MP2/CBS. In the first case, the differences are
0.25 and
1.09 kcal/mol for H2O and OH
, respectively, whereas in the second case such differences are 0.22 kcal/mol (H2O) and
0.20 kcal/mol (OH
). When the core
valence correlation is taken into account, these differences become
0.25 kcal/mol (H2O) and
0.80 kcal/mol (OH
) for the CCSD(T)/aug-cc-pVQZ to CCSD(T)/aug-ccpV5Z and 0.21 kcal/mol (H2O) and
0.09 kcal/mol (OH
), for the CCSD(T)/aug-cc-pV5Z to CCSD(T)/CBS. In summary, the largest factors that in principle could affect the benchmark “composite” energies are the deep core
valence correlation effects and BSSE. Nevertheless, and according to the different De energies reported in Table 1, the combined effects of the two factors should not be larger than 1
2 kcal/mol. Consequently, carrying out the benchmark calculations within the frozen core approximation without including the CP correction represents a reasonable compromise between cost and accuracy, which in turn should be close to chemical accuracy (∼1 kcal/mol). 11336

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A

ARTICLE

Table 3. Protonation Energies (PEs in kcal/mol) for Various Acid
Base Pairs Both in Their Isolated Form and in the Corresponding Ca(II) Complexesa MP2/CBSb

composite methodc


397.15


393.13


397.02 (
389.01)


353.60
388.27


355.99
391.13


353.38
388.11


355.77 (
349.97)
390.98 (
381.09)

CH3SH/CH3S



360.25


362.65


360.18


362.58 (
356.17)

HCOOH/HCO2



344.37


347.41


344.12


347.17 (
338.79)

CH3COOH/CH3CO2



352.52


355.27


352.31


355.05 (
346.44)

1H-imidazol/1H-imidazolate


352.91


356.78


352.89


356.76 (
348.01)

Acid
Base Pair

MP2/aug-cc-pVTZ

CCSD(T)/aug-cc-pVTZ

H2O/OH



393.26

H2S/HS
CH3OH/CH3O


A/B

Ca
A/Ca
B


H2O/OH


112.04


114.90


108.82


111.68 (
103.90)

H2S/HS



113.59


115.92


109.81


112.14 (
106.05)

CH3OH/CH3O
CH3SH/CH3S



119.02
125.50


122.26
128.16


115.87
122.71


119.11 (
111.82)
125.36 (
119.75)

HCOOH/HCO2



110.21


112.91


107.57


110.27 (
102.48)

CH3COOH/CH3CO2



112.76


115.36


110.30


112.90 (
105.67)

1H-imidazol/1H-imidazolate


171.24


174.97


170.60


174.33 (
165.44)

a

Molecular geometries were obtained at the MP2/aug-cc-pVTZ level. b Obtained from CBS extrapolation of the MP2 correlation energy based on eqs 1) and (2 and using the HF/aug-cc-pV5Z energies. c Using an additive combination of electronic energies (CCSD(T)/aug-cc-pVTZ+ MP2/CBS
MP2/ aug-cc-pVTZ). Values in parentheses include the MP2/aug-cc-pVTZ ZPVE energies.

Regarding geometries, we have observed that small changes on the equilibrium geometries of the Ca(II) water/hydroxide complexes affect little the binding energies. Thus, the structures optimized at the frozen-core or full-electron CCSD/aug-cc-pVTZ level have slightly larger Ca
O distances, by 0.01
0.04 Å, than the reference MP2 structures, but very similar De values are obtained (see Table 1). Additionally, NPA charges of the Ca atom computed using the MP2 and CCSD electronic densities are similar and the MP2 and CCSD harmonic stretching frequencies for the Ca
ligand bond match each other closely (see Figures 1 (Ca(II) 3 3 3 H2O) and 2 (Ca(II) 3 3 3 OH
)). Ab Initio Benchmark Data. Table 2 contains ab initio binding energies for the full set of Ca(II) complexes here studied calculated at various levels of theory. The “composite” De values (5th column) provide the benchmark data. Notice however that the MP2/aug-cc-pVTZ level (2nd column) shows a good performance as compared with the “composite” energy calculations. Indeed, binding energies computed on MP2/aug-cc-pVTZ geometries are quite similar regardless of the correlated method (CCSD(T) vs MP2) or the basis set used to perform the corresponding single-point energy calculations. Also true is that including energy corrections by using higher correlated methods and larger basis sets does not alter significantly the MP2/aug-ccpVTZ De values. In fact, those two improvements frequently “work” in opposite directions and turn out canceling each other to some extent implying that the average difference between the MP2/aug-cc-pVTZ calculations and the more accurate “composite” energy values is relatively small: the mean unsigned percentage error (MU%E) for neutral ligands is 0.7% (largest error for imidazole, 2.4%) and for anionic ligands 0.9% (largest error for hydrosulfide, 1.4%). Then, all these results suggest that predictions obtained at the MP2/aug-cc-pVTZ level of theory for Ca(II) complexes should be reliable enough. Certainly, this is a very satisfactory conclusion since both CCSD(T)/aug-cc-pVTZ and MP2/aug-cc-pVnZ (n = 4, 5) calculations are significantly

more computationally expensive. Table 2 also shows that most of the improvement when going from MP2/aug-cc-pVTZ to the composite values comes from the basis set extension. Thus, the MU%Es between the MP2/CBS and the composite values are merely 0.5% for neutral ligands (largest error for benzene, 1%) and 0.2% for anionic ligands (largest error for hydroxide, 0.3%). Indeed the mean unsigned errors (MUEs) are in all cases below 1 kcal/mol. For the neutral ligands, the optimized structures correspond to monodentate Ca(II) complexes characterized by a single interatomic Ca
X distance of 2.12, 2.31, and 2.78 Å for X = O, N, and S, respectively (averages of the corresponding Ca
X distances). It turns out that, for a giving donor atom, the binding energy correlates well with the Ca
X distance (e.g., the linear correlation coefficient, r, between the Ca
X distances and the composite De values amounts to 0.813 and 0.972 for X = O and N, respectively). NPA charges suggest a significant charge transfer when Ca(II) binds to neutral ligands: within the 0.08
0.2 e
range. The amount of charge transfer depends on the electronegativity and hybridization of the donor X atom and we have found that it is relatively well correlated with the De energies (r ∼
0.882 (for X = O) and r ∼
0.844 (for X = N)). This suggests that charge-transfer interactions could also contribute to the total binding of these complexes. We also point out that the frequency of the stretching motion of the Ca
X bond is uncorrelated with the strength of the metal
ligand bond expressed in terms of binding energies. The strong interaction between the Ca(II) cation and the anionic ligands is reflected in the magnitude of the computed binding energies (from
338 to
280 kcal/mol), which is accompanied by a typical shortening of the Ca
X distances within a range of 0.40
0.1 Å with respect to the neutral ligands and an important charge transfer from ligand to Ca(II) (0.21
0.12 e
). For the sake of completeness, we have also derived the protonation energies (PE) of some of the ligands studied in this 11337

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A

ARTICLE

Table 4. Mean Unsigned Percentage Errors (MU%Es) in the Metal
Ligand Bond Distances of the Ca(II) Complexes with Respect to the MP2/aug-cc-pVTZ Reference Values basis set

PW91

PBE

B3LYP

TPSS

MP2

6-31G(d)

1.4

1.2

0.9

1.2

1.1

6-31+G(d)

0.7

0.6

0.4

0.6

1.7

6-31+G(2d,p)

1.3

1.2

0.8

1.1

1.2

6-311+G(d)

0.4

0.4

0.6

0.4

2.0

aug-cc-pVTZ

1.3

1.2

0.8

1.1

work at the different correlated levels of theory. According to the PE values reported in Table 3, the nature of the correlated method (CCSD(T) vs MP2) has a more pronounced effect (∼2 kcal/mol) on the computed PEs than in the case of the Ca
L binding energies, whereas the CBS limit of the MP2 PEs is very similar to the MP2/aug-cc-pVTZ values. When comparing the PEs of the isolated A/B pairs with those of their Ca-bound counterparts, we see that the Ca(II) cation decreases the PEs by an amount ranging within
285,
183 kcal/mol (composite energy values including ZPVE) for the different acid
base pairs, respectively. Thus, the actual impact on the acid
base strength of the Ca-ligands is quite variable, with the water/hydroxide pair being the most affected one. DFT Calculations. The performance of the test DFT functionals will be discussed in terms of Mean Unsigned Percentage Errors (MU%Es). We point out that in our previous study on monoligand Zn(II) complexes, we mainly discussed mean signed errors (MSEs). The reason was that for the zinc complexes all tested functionals showed a systematic bias with respect to the reference values (either underestimating or overestimating them) for a given property and basis set. For the calcium complexes, we have not observed such systematic deviations and we have therefore chosen MU%Es to assess the performance of the density functionals. We also note that, for the sake of simplicity, part of the results of the present study has been collected in the tables presented in the Supporting Information. DFT Molecular Geometries of the Ca
L complexes. For the majority of the complexes analyzed in this work DFT methods, when used in conjunction with Pople-type basis sets, predict the same geometry for the lowest energy structure than the MP2/aug-cc-pVTZ reference level. There are, however, few exceptions that deserve to be commented. The complex Ca(II)
CH3COOH has a C
O
Ca angle of 104° (see Figure 1) at the MP2/aug-cc-pVTZ level of theory. With the Pople-type basis sets, particularly with the 6-31+G(2d,p) basis, this angle tends to open up to 180°. Indeed, with Pople-type basis sets, both isomers (C
O
Ca ∼ 104° and C
O
Ca ∼ 180°) are true minima on the potential energy surface except for the 6-31+ G(2d,p) basis for which the first isomer does not exist. We have nevertheless considered the MP2/aug-cc-pVTZ-like isomer (C
O
Ca ∼ 104°) for this complex, except with the 6-31+ G(2d,p) basis set. The Ca(II) complexes with N-methyl acetamide and methyl acetate show one imaginary frequency with some functional and basis sets (the values range between 30 and 70 cm
1). Appearance of this imaginary frequency is due to the rotation of methyl groups and has a negligible effect in the binding energies. We therefore considered for these complexes the MP2/aug-cc-pVTZ geometries shown in Figure 1. The Ca(II) complex with methanethiolate has C3v symmetry at the MP2/aug-cc-pVTZ level, but Cs symmetry with the Pople-type

Table 5. Mean Unsigned Percentage Errors (MU%Es) in the Binding Energies of the Ca(II) Complexes with Neutral and Anionic Ligands basis set

PW91

PBE

B3LYP

TPSS

MP2

Neutral Ligands 6-31G(d)

11.1

9.8

6.4

7.9

4.1

6-31+G(d)

6.5

5.4

2.8

4.3

4.0

6-31+G(2d,p)

5.5

4.5

2.1

3.4

4.6

6-311+G(d) aug-cc-pVTZ

5.1 6.7

4.2 5.9

2.3 3.4

3.2 4.6

5.3

6-31G(d)

9.1

8.6

6.5

7.5

2.2

6-31+G(d)

2.8

2.5

0.9

2.1

2.8

6-31+G(2d,p)

2.8

2.4

0.8

2.0

2.6

6-311+G(d)

2.0

1.8

0.9

1.4

3.3

aug-cc-pVTZ

3.5

3.2

1.5

2.7

Anionic Ligands

basis sets (also with the MP2 method). The C
S
Ca angle varies between 130 and 170°, with B3LYP and MP2 yielding the smaller values. We have considered these Cs symmetry structures in this work. Finally, and in agreement with what we found in our previous study on Zn(II) complexes, every functional in conjunction with Pople basis sets predicts C1 symmetry for the imidazolate complexes. This is in disagreement with the MP2/ aug-cc-pVTZ results, which predict that the planar Cs structure is the true minimum on the potential energy surface. At the B3LYP/aug-cc-pVTZ level of theory, the complex has again Cs symmetry, so it seems that the disagreement between the DFT and the MP2/aug-cc-pVTZ geometries is due to the Pople-type basis sets. For comparison purposes, we have considered a Cs symmetry for this complex. Table 4 collects the MU%E in the metal
ligand bond distances of all the calcium complexes here studied. For the present analysis, MP2/aug-cc-pVTZ geometries have been taken as the reference ones. We point out that the MU%Es have been computed for the whole set of complexes including neutral and anionic ligands since similar errors were obtained for both of them. All the functional we have tested give, in general, metal
ligand bond distances that are shorter than the MP2/aug-ccpVTZ reference values. With the 6-311+G(d) basis set, the distances may be overestimated for some complexes, but this is due to the fact that this basis yields distances very close to the reference values. The DFT/aug-cc-pVTZ distances are also systematically shorter than the reference. Thus, the general trend is that DFT functionals give too short Ca(II)
ligand bond distances. This result is interesting because the Zn(II) complexes show the opposite behavior.10 Table 4 shows, however, that the errors for the Ca(II) complexes are not very large, the worst results are obtained at the PW91/6-31G(d) level, and their MU% E is 1.4%. Focusing on the different basis sets, the following behavior can be observed: MU%E(6-31G(d)) > MU%E(6-31+ G(2d,p)) > MU%E(6-31+G(d)) > MU%E(6-311+G(d)). The performance of the MP2 method with these basis sets is totally different since (a) MP2 overestimates the metal
ligand bond distances with Pople basis sets, (b) MP2 metal
ligand bond distances show larger deviation with respect to the reference ones than those obtained with the various DFT functionals, and (c) comparison of the MU%Es obtained with the four Pople basis sets show the opposite trend to that found with the DFT 11338

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A

ARTICLE

Table 6. Mean Unsigned Percentage Errors (MU%Es) in the Counterpoise (CP) Corrected Binding Energies of the Ca(II) Complexes with Neutral and Anionic Ligands basis set

PW91

PBE

B3LYP

TPSS

MP2

Neutrals 6-31G(d)

8.6

7.4

4.3

5.8

6.0

6-31+G(d)

5.6

4.6

2.2

3.5

6.7

6-31+G(2d,p)

5.0

4.0

1.8

3.0

7.0

6-311+G(d) aug-cc-pVTZ

4.2 3.8

3.4 1.8

1.9 2.7

2.4 2.1

7.5

6-31G(d)

7.2

6.8

4.7

6.0

1.9

6-31+G(d)

2.4

2.1

0.8

1.7

4.0

6-31+G(2d,p)

2.3

2.0

0.6

1.6

3.6

6-311+G(d)

1.6

1.4

0.7

1.0

4.5

aug-cc-pVTZ

2.3

0.9

2.0

1.5

Anionic

functionals, namely, MU%E(6-31G(d)) < MU%E(6-31+G(2d,p)) < MU%E(6-31+G(d)) < MU%E(6-311+G(d)). Table 4 also shows that the DFT/aug-cc-pVTZ results are comparable to those obtained with the 6-31G(d) basis but worse than those obtained with the larger 6-31+G(d), 6-31+G(2d,p), and 6-311+ G(d) basis sets. The MP2 and DFT/aug-cc-pVTZ results clearly point to an error cancellation between the underestimation of the metal
ligand distances by the tested density functionals and their overestimation by Pople-type basis sets. Let us finally point out that geometry optimizations carried out using a larger (“Ultrafine”) integration grid at the B3LYP/aug-cc-pVTZ level differ very little from the “medium” grid results: the largest difference in the metal
ligand bond distances is 0.003 Å (see Table S3 in the Supporting Information). DFT Binding Energies. Tables 5 and 6 collect the MU%Es of the DFT binding energies, corrected and noncorrected from the BSSE, with respect to the reference composite ones for the different Ca(II)
L complexes here studied. Notice that we will discuss the DFT binding energies of the neutral and anionic ligands separately because their binding energies for complexation are very different in magnitude. In our previous study on Zn(II) complexes,10 we noticed that all tested functionals overshot the composite CCSD(T)/CBS energies and we related this fact to the overestimation of polarizabilities by the density functionals.69 Indeed, for the zinc complexes we showed that induction forces (collecting charge transfer and polarization) played a relevant role in the metal
ligand bonding interactions. As we will show here, the case of the calcium complexes is analogous, that is, the chosen functionals tend to overestimate the interaction energies. However, the overestimation is not so large and, in some particular cases, even underestimations are occasionally observed. We will see later that, in fact, induction effects in Ca(II) complexes are much smaller than in Zn(II) complexes. We point out in passing that the larger binding energies found for the Zn(II) complexes with respect to Ca(II) nicely correlate with the higher ionization energy of zinc that should induce a larger polarization and charge transfer from the bases to the metal cations. Taking into account the importance of induction effects in the binding of these complexes (recall also the correlation between the NPA charges and De energies commented above), this correlation is indeed what should be expected and suggests that consideration of electrostatic

contributions alone may not be reliable enough for the prediction of relative binding energies. This suggestion will be address in further detail below. The results of the complexes with neutral ligands will be discussed first (see Table 5). In general, all the tested density functionals overestimate the binding energies. Only for the 6-311+ G(d) basis set the results are close enough to the reference values [CCSD(T)/CBS] that some underestimations are observed. Table 5 shows that, for the Pople-type basis sets, the MU%E ranges from 2.1 to 11.1%. The largest errors are found with the 6-31G(d) basis set. Inclusion of diffuse and further polarization functions, in particular, 6-31+G(2d,p) and 6-311+G(d), diminishes the difference between DFT and the reference data to 2.1
5.5%. Because these two bases provide a better agreement with the reference data than aug-cc-pVTZ, it seems that some error compensation is taking place when using Pople-type sets. As a matter of fact, the MP2 energies do not show any improvement when increasing the basis set size. Regarding the assessment of the tested functionals, from the analysis of Table 5, it seems clear that both B3LYP and TPSS functionals, particularly the first one, show the best performance. Because, as anticipated above and discussed in further detail below, induction forces are the largest contributors to the interaction energy on these complexes, we believe that the more or less accurate predictions obtained by different density functionals should be related to a correct description of this particular intermolecular force. As expected, the complexes with anionic ligands, as compared to the neutral ligands, give larger absolute errors. However, relative errors are smaller than those previously discussed for the neutral ligands (see Table 5). For instance, the B3LYP functional gives MU%E = 1.5% with the aug-cc-pVTZ basis set (for neutral ligands MU%E is 3.4%) and MU%E = 0.8% when the 6-31+ G(2d,p) basis set is used (for neutral ligands is 2.1%). We observe again that Pople-type basis sets yield systematically better results than the aug-cc-pVTZ basis and that the B3LYP and TPSS functionals perform the best. We would like to point out that binding energies computed with a larger (“Ultrafine”) integration grid at the B3LYP/aug-cc-pVTZ level are within 0.1 kcal/ mol from the “medium” grid results except for two complexes, namely, Ca(II)
acetic acid (difference: 0.34 kcal/mol) and Ca(II)
methyl acetate (0.15 kcal/mol; see Table S4). Thus, larger integration grids are in principle not required for the calculation of binding energies in these Ca(II) complexes. Table 6 collects the MU%Es corresponding to the binding energies corrected for the BSSE. For both neutral and anionic ligands, comparison of the data shown in Tables 5 and 6 show that the CP-corrected DFT energies are in general in better agreement with the reference values than the uncorrected ones. This is not surprising because the tested functional overestimates the binding energies with errors in most cases larger than the BSSE itself. However, the improvement obtained after correcting the binding energies from the BSSE is not very large. As a matter of fact, the average BSSE value at the B3LYP/6-311+G(d) level for the set of neutral ligands is not very large, 0.59 kcal/mol (the largest value: 1.56 kcal/mol for water), and for the anionic ligands, 0.93 kcal/mol (largest value: 2.58 kcal/mol for hydroxide). As expected, the BSSE is larger at the MP2 level of theory: for the 6-311+G(d) basis set we found an average value for neutrals of 2.2 kcal/mol (largest: 3.91 kcal/mol for benzene) and for the anionic ligands 3.95 kcal/mol (largest: 4.84 kcal/mol for hydrosulfide). Because Pople-type basis sets underestimate the binding energies, the corrected MP2 energies show worst 11339

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A

ARTICLE

Table 7. Mean Unsigned Percentage Errors (MU%Es) in the Protonation Energies of the Isolated and Complexed Anionic Ligands (L = OH
, CH3O
, HCO2
, CH3CO2
, Im
, HS
, CH3S
) basis set

PW91

PBE

B3LYP

TPSS

MP2

Isolated 6-31G(d)

3.4

3.3

3.6

3.7

3.0

6-31+G(d)

1.8

1.8

1.3

1.0

1.8

6-31+G(2d,p)

0.8

0.8

0.3

0.3

0.6

6-311+G(d)

1.9

1.9

1.4

1.0

1.7

aug-cc-pVTZ

0.8

0.9

0.4

0.4

6-31G(d) 6-31+G(d)

6.5 9.0

6.4 8.9

2.3 4.8

3.8 6.0

3.4 1.3

6-31+G(2d,p)

5.9

5.7

2.2

3.0

2.4

6-311+G(d)

7.4

7.3

3.3

4.3

1.2

aug-cc-pVTZ

7.1

7.0

2.9

4.2

Complexed

agreement with the reference than the uncorrected values. On the other hand, the tested density functionals in conjunction with the aug-cc-pVTZ basis yield small BSSEs: the B3LYP average for neutral ligands is 0.20 kcal/mol (largest: 0.37 kcal/mol for benzene) and for anionic ligands 0.40 kcal/mol (largest: 0.65 kcal/mol for hydroxide). Thus, taking into account the relatively small effect of the BSSE on the binding energies, on the one hand, and the computational effort needed for the calculation of the CP-corrected energies (which involves three additional calculations: both fragments in the dimer centered basis set and the ligand in its complexed geometry with its own basis set) on the other, it seems that raw binding energies provide a better quality/cost ratio than the CP-corrected ones. It is finally worth stressing that the CP method does not correct for basis set incompleteness errors apart from the BSSE itself. Thus, the 6-31G(d) binding energies in Table 6 show again the worst agreement with the reference. In principle, computation of relative energies for ligand exchange processes can benefit from cancellation of “systematic” errors in the computed binding energies (De). Hence, we also obtained the linear correlation coefficients (r) of DFT results versus the ab initio reference results in order to ascertain more clearly the differences and similarities in the systematic deviations of the computed energies (discriminating between the neutral and the anionic ligands; see Table S5 in the Supporting Information). Thus, we found that the DFT binding energies of neutral ligands show quite a good correlation (r > 0.99) with respect to the reference data for all the combinations of DFT functionals and basis sets used. For the anionic ligands, the linear regression fits of the binding energies tend to be less satisfactory, especially for the small 6-31G(d) basis set that has r values ∼0.96
0.97. However, when the 6-31+G(2d,p) and aug-ccpVTZ basis sets are used, the trends observed in the DFT binding energies of the neutral and anionic ligands match quite closely the trends observed in the ab initio reference energies, resulting in similar r coefficients. For example, the B3LYP/aug-cc-pVTZ De values for the neutral and anionic ligands have the same r value (0.9967). Therefore, we expect that the relative energies for processes involving the rupture of bonds between Ca(II) and either neutral or anionic ligands should be determined quite

accurately using a DFT method combined with a flexible basis set. DFT Protonation Energies. Table 7 collects the MU%Es in the protonation energies (PEs) of the isolated acid
base pairs as well as for the complexed anionic ligands. A first observation that can be made from the analysis of these data is that the former (isolated acid
base pairs) show significantly smaller MU%E than the latter. Indeed, unless for the 6-31G(d) basis set, MU%Es are below 2% and typically closer to 1%. The best performance of all the functional tested is obtained, as in the case of the binding energies, with the B3LYP and TPSS functionals, for which the MU%Es are quite close to zero (0.3% each), particularly in conjunction with the 6-31+G(2d,p) basis set. For the complexed anionic ligands, results collected in Table 7 show MU%Es that are clearly larger than those for the isolated ligands. Moreover, the computed protonation energies display similar or even slightly larger MU%Es with respect to the reference than the binding energies (compare to Table 5). We also note that for PEs the MP2 values are in better agreement with respect to the reference results than the DFT ones. As a general conclusion, it can be said then that one should be careful when using complexed anionic ligands PEs, even those obtained with the best levels of theory, in order to find out whether a ligand is protonated or deprotonated when coordinated to the Ca(II) cation. DFT-SAPT Analyses. We would finally like to comment briefly on the nature of the interaction of the Ca(II) complexes with neutral ligands. Figure 3 shows the E(1) and E(2) contributions (see Computational Methods for the definitions) as well as the total DFT-SAPT interaction energies for the 17 Ca(II) complexes with neutral ligands. A complete set of numerical DFTSAPT data has been collected in Tables S6 and S7 in the Supporting Information. Let us first comment briefly on two general conclusions that can be drawn from the DFT-SAPT results. First, it is encouraging to see that the DFT-SAPT interaction energies are in very good agreement with the composite ones. The average percentage error for the seventeen Ca(II) complexes is a mere 2.9% being the largest deviation 6.5% (5.6 kcal/mol) in Ca(II) 3 3 3 Bz. Thus, we conclude that DFTSAPT provide for these systems fairly reliable interaction energies in spite of their relatively large magnitude. Second, the δEHF int term is not very large for these calcium complexes being in all cases smaller than 5% of the total DFT-SAPT interaction energy (except Ca(II) 3 3 3 Bz, 12.4%). We point out that for the Zn(II) complexes δEHF int represented 10
20% of the total SAPT2 energy, showing that third-order contributions are far more important for these systems in agreement with their larger binding energies. The error in the DFT-SAPT interaction energies without inclusion of the δEHF int term is 4.0%, not much larger than the error obtained including this term. Indeed, for 8 out of the 17 complexes, inclusion of δEHF int slightly worsens the agreement with the reference. Explicit calculation of third order corrections, which is feasible within DFT-SAPT but without response, is much more computationally demanding and, taking into account the relatively small magnitude of the δEHF int term, has not been attempted in this work. Figure 3 shows that for the whole set of complexes the secondorder terms are larger than the first-order ones. Since dispersion is almost negligible for these systems (see Table S6), the DFTSAPT results highlight the importance of charge transfer and polarization contributions to the metal
ligand bonding in Ca(II) complexes. For the sulfur-containing ligands, E(1) is indeed 11340

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A

ARTICLE

Figure 3. First- and second-order contributions (kcal/mol) to the interaction energy of the 17 Ca(II) complexes with neutral ligands, as obtained at the DFT-SAPT/aug-cc-pVTZ level of theory.

Figure 4. First- and second-order contributions (kcal/mol) to the interaction energy of 12 Ca(II) and Zn(II) complexes, as obtained at the DFTSAPT/aug-cc-pVTZ level of theory.

very small and most of the binding energy comes from the second-order terms. In the particular case of Ca(II) 3 3 3 Bz, the first-order contributions are almost negligible, 0.41 kcal/mol, and the stability of the interaction comes again from the relaxation of the charge densities. This result highlights once more the importance of polarization effects when an aromatic ring interacts with an ionic partner.70
73 The Ca(II) binding energies to formaldehyde and ammonia have been analyzed before by Corral et al.74 Interestingly, the relative binding energies of these two complexes (larger for formaldehyde) seem to be reversed with respect to their protonation affinities (larger for ammonia). The analysis of Corral et al. lead them to conclude that the main contribution to the binding preference of Ca(II) appear to stem from polarization effects. The SAPT results shown in Figure 3 (see also Table S6) show indeed that the electrostatic interaction is larger for the ammonia complex but that the induction contributions (second, E(2) plus higher order, δEHF int ) reverse the final interaction energies. We therefore agree with Corral et al. in that classical electrostatics alone does not seem to be able to reproduce the Ca(II) binding preference for these two complexes. Interestingly, the binding energies of ammonia and formaldehyde

to Zn(II) are in agreement with their relative proton affinities.10 SAPT shows (see Table 11 in ref 10) that, for these complexes, the electrostatic contribution favors the ammonia complex by 20 kcal/mol, whereas the induction forces favor formaldehyde by about 5 kcal/mol. This is the same qualitative trend found for the Ca(II) complexes. However, the electrostatic preference for ammonia in the Zn(II) complex is so large that it overcomes the preference for formaldehyde in terms of the induction forces and the final binding energy is larger for the former. These two examples stress again the important role played by induction forces and therefore suggest that models based on classical electrostatic interactions may not be reliable enough for the prediction of relative ligand binding energies. We would finally like to point out that Gresh and Garmer studied some of the complexes considered in this work using a different partitioning scheme, the Restricted Variational Space (RVS) procedure, and arrived to similar conclusions about the role of the polarization effects.75,76 We also point out that SAPT results for the complexes with water, formaldehyde, and benzene have also been recently published.32,77 It is worthwhile to compare the Ca(II) results to those obtained in our previous work for Zn(II).10 In Figure 4 we compare 11341

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A the E(1) and E(2) contributions for 12 Ca(II) and Zn(II) complexes. Because the Zn(II) results had been obtained at the SAPT2/6-31+G(2d,p) level of theory,10 we have recalculated four representative complexes (namely, those with water, ammonia, benzene, and hydrogen sulfide) at the DFT-SAPT/augcc-pVTZ level just to confirm that the results obtained are very similar (this is not indeed an unexpected result since we already showed10 that for Zn(II) complexes the 6-31+G(2d,p) basis yielded results in good agreement to those obtained with aug-ccpVTZ). Figure 4 clearly shows that, whereas the magnitude of the first-order term is similar for Ca(II) and Zn(II) complexes, the second-order contributions are about 2.5
3 times larger for the later. As above-mentioned, the dissociation energies of the Ca(II) complexes are about half of the corresponding Zn(II) ones. We related this fact to the higher ionization energy of zinc, which should induce a larger charge transfer from the bases to the metal cation. Figure 4 and Table S7 show that this is indeed the case: it is the much larger polarization of the ligands’ charge densities in the Zn(II) complexes that yields larger interaction energies for these compounds. We would like to stress again that the different roles played by the second-order forces in the interaction energy of Ca(II) and Zn(II) complexes has necessary implications for the description of these systems using molecular mechanics models. Indeed, it may well be that models based on classical electrostatics alone may not yield reliable relative binding energies, taking into account the importance of polarization effects. We would like to finish this section, stressing the importance of performing rigorous bonding analysis by means of methodologies like SAPT such that further developments in molecular mechanics and force fields can reflect the underlying physical forces that govern the interaction potentials in these complexes. We point out that efforts in this direction are already on the way (see, for example,refs 25 and 78).

’ CONCLUSIONS In this study we have carried out high-level ab initio calculations for a set of monoligand Ca(II) complexes to provide a reliable set of benchmark results on metal
ligand distances, binding energies, and protonation energies. Geometries have been obtained at the MP2/aug-cc-pVTZ level of theory, and a “composite” CCSD(T)/CBS scheme has been chosen as the reference method for the protonation and binding energies. We believe that these levels of theory should provide results reliable enough for benchmark purposes. We have, for example, estimated that the “composite” CCSD(T)/CBS binding energies should be already close to chemical accuracy (∼1 kcal/mol). An interesting first conclusion of this work is that MP2/aug-ccpVTZ binding energies are very close to the “composite” reference ones, the MU%E is smaller than 1%, meaning that they should be within 1 kcal/mol from the CCSD(T)/CBS results. We have additionally tested four widely used density functionals and several Pople-type basis sets (including an increasing number of polarization and diffuse functions). The B3LYP and TPSS functionals provide the closest agreement with the reference values, particularly the former. In conjunction with the 6-31+ G(2d,p) and 6-311+G(d) basis sets, both functionals yield MU% Es in the binding energies below 3%. Besides, we have shown that DFT binding energies yield a very good linear correlation when compared to the “composite” reference results. Therefore, the computation of accurate relative energies, which are important for the study of ligand exchange processes, should provide

ARTICLE

reliable results when using a DFT method combined with a flexible basis set. We have additionally estimated the relative weight of the different contributions to the metal
ligand bond using the Symmetry Adapted Perturbation Theory. This analysis should provide useful information for the assessment of current “nonpolarizable” molecular mechanics potentials and for the future development of polarizable force fields. Our results have shown that induction and electrostatic forces contribute similarly to the metal
ligand bond in the studied Ca(II) complexes, suggesting that an explicit treatment of polarization is advisible to properly describe the coordination structure of Ca(II) complexes.

’ ASSOCIATED CONTENT

bS

Supporting Information. Mean unsigned percentage errors in the binding energies and metal
ligand bond distances for the whole set of complexes, as obtained with Blaudeau’s and Rassolov’s 6-31G basis sets, metal
ligand bond distances, and binding energies at the B3LYP/aug-cc-pVTZ level with an “ultrafine” grid; linear correlation coefficients corresponding to the binding energies between BFT/MP2 and the benchmark results, detailed SAPT contributions for the Ca(II) complexes; and comparative SAPT results for some selected Ca(II) and Zn(II) complexes. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Present Addresses §

Departament de Fisicoquímica, Facultat de Farmacia, Universitat de Barcelona, 08028 Barcelona, Spain.

’ ACKNOWLEDGMENT This research has been supported by the Junta de Castilla y Le
on (Grant VA006B07). N.D. and D.S. are also grateful to the Spanish MEC for support via Grant CTQ2007-63266. This article is dedicated to the memory of our colleague and friend
lvarez Blanco. Dr. Miguel A ’ REFERENCES (1) Even€as, J.; Malmendal, A.; Fors
en, S. Curr. Opin. Chem. Biol. 1998, 2, 293. (2) Dudev, T.; Lim, C. Chem. Rev. 2003, 103, 773. (3) Berridge, M. J.; Bootman, M. D.; Lipp, P. Nature 1998, 395, 645. (4) Maurer, P.; Hohenester, E. Matrix Biol. 1997, 15, 569. (5) Harding, M. M. Acta Crystallogr., Sect. D 2000, 56, 857. (6) Rao, J. S.; Dinadayalane, T. C.; Leszcynsky, J.; Sastry, G. N. J. Phys. Chem. A 2008, 112, 12944. (7) Tunell, I.; Lim, C. Inorg. Chem. 2006, 45, 4811. (8) Harding, M. M. Acta Crystallogr., Sect. D 1999, 55, 1432. (9) Harding, M. M. Acta Crystallogr., Sect. D 2001, 57, 401. (10) Ray
on, V. M.; Vald
es, H.; Díaz, N.; Su
arez, D. J. Chem. Theory Comput. 2008, 4, 14. (11) Riley, K. E.; Pitonak, M.; Jurecka, P.; Hobza, P. Chem. Rev. 2010, 110, 5023. (12) Grimme, S. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 211. (13) Ramos, M. J.; Fernandes, P. A. Acc. Chem. Res. 2007, 41, 689. 11342

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343

The Journal of Physical Chemistry A (14) Geerlings, P.; De Proft, F.; Langenaeker, W. Chem. Rev. 2003, 103, 1793. (15) Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. J. Phys. Chem. A 2007, 111, 10439. (16) Alcamí, M.; Gonz
alez, A. I.; M
o, O.; Y
a~nez, M. Chem. Phys. Lett. 1999, 307, 244. (17) Corral, I.; M
o, O.; Y
a~nez, M.; Scott, A. P.; Radom, L. J. Phys. Chem. A 2003, 107, 10456. (18) da Costa, L. M.; Carneiro, J. W. d. M.; Coelho Paes, L. W.; Romeiro, G. A. J. Mol. Struct.: THEOCHEM 2009, 911, 46. (19) da Costa, L. M.; de Mesquita Carneiro, J. W.; Romeiro, G. A.; Paes, L. W. C. J. Mol. Model. 2011, 17, 243. (20) Haworth, N. L.; Sullivan, M. B.; Wilson, A. K.; Martin, J. M.; Radom, L. J. Phys. Chem. A 2005, 109, 9156. (21) Petrie, S.; Radom, L. Int. J. Mass Spectrom. 1999, 192, 173. (22) Viswanathan, B.; Barden, C. J.; Ban, F.; Boyd, R. J. Mol. Phys. 2005, 103, 337. (23) Barnes, P.; Finney, J. L.; Nicholas, J. D.; Quinn, E. E. Nature 1979, 282, 459. (24) Cieplak, P.; Dupradeau, F.-Y.; Duan, Y.; Wang, Y. J. Phys.: Condens. Matter 2009, 21, 333102. (25) Luque, F. J.; Dehez, F.; Chipot, C.; Orozco, M. WIREs Computational Molecular Science 2011, 1, 844. (26) Leontyev, I.; Stuchebrukhov, A. Phys. Chem. Chem. Phys. 2011, 13, 2613. (27) Lopes, P. E. M.; Roux, B.; MacKerell, A. D., Jr Theor. Chem. Acc. 2009, 124, 11. (28) Warshel, A.; Kato, M.; Pisliakov, A. V. J. Chem. Theory Comput 2007, 3, 2034. (29) Wu, J. C.; Piquemal, J.-P.; Chaudret, R.; Reinhardt, P.; Ren, P. J. Chem. Theory Comput. 2010, 6, 2059. (30) Project, E.; Nachliel, E.; Gutman, M. J. Comput. Chem. 2008, 29, 1163. (31) de Courcy, B.; Piquemal, J.-P.; Gresh, N. J. Chem. Theory Comput. 2008, 4, 1659. (32) Dehez, F.; Archambault, F.; Soteras Guti
errez, I.; Luque, F. J.; Chipot, C. Mol. Phys. 2008, 106, 1685. (33) Frison, G.; Ohanessian, G. J. Comput. Chem. 2008, 29, 416. (34) Kol
ar, M.; Berka, K.; Jurecka, P.; Hobza, P. ChemPhysChem 2010, 11, 2399. (35) Zgarbov
a, M.; Otyepka, M.; Sponer, J.; Hobza, P.; Jurecka, P. Phys. Chem. Chem. Phys. 2010, 12, 10476. (36) Piquemal, J.-P.; Chevreau, H.; Gresh, N. J. Chem. Theory Comput. 2007, 3, 824. (37) de Courcy, B.; Pedersen, L. G.; Parisel, O.; Gresh, N.; Silvi, B.; Pilm, J.; Piquemal, J.-P. J. Chem. Theory Comput. 2010, 6, 1048. (38) Kendall, R. A.; Dunning, J. T. H.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796. (39) Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 98, 1358. (40) Koput, J.; Peterson, K. A. J. Phys. Chem. A 2002, 106, 9595. (41) Feller, D. J. Comput. Chem. 1996, 17, 1571. (42) Schuchardt, K. L.; Didier, B. T.; Elsethagen, T.; Sun, L.; Gurumoorthi, V.; Chase, J.; Li, J.; Windus, T. L. J. Chem. Inf. Model. 2007, 47, 1045. (43) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479. (44) Peterson, K. A.; Puzzarini, C. Theor. Chem. Acc. 2005, 114, 283. (45) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (46) Iron, M. A.; Oren, M.; Martin, J. M. Mol. Phys. 2003, 101. (47) Perdew, J. P.; Burke, K.; Wang, Y. Phys. Rev. B 1996, 57, 14999. (48) Perdew, J. P.; Burke, K.; Enzerhof, M. Phys. Rev. Lett. 1996, 78, 1396. (49) Becke, A. D. J. Chem. Phys. 1986, 95, 7184. (50) Becke, A. D. J. Chem. Phys. 1988, 88, 2547. (51) Becke, A. D. J. Chem. Phys. 1988, 88, 1053. (52) Lee, C.; Yang, R. G.; Parr, R. G. Phys. Rev. B 1988, 37, 785. (53) Tao, J.; Perdew, J. P.; Starovenov, V. N.; Scuseria, G. E. Phys. Rev. Lett. 2003, 91, 146401.

ARTICLE

(54) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (55) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 4654. (56) McLean, A. D.; Chandler, G. S. J. Chem. Phys. 1980, 72, 5639. (57) Rassolov, V. A.; Pople, J. A.; Ratner, M. A.; Windus, T. L. J. Chem. Phys. 1998, 109, 1223. (58) Rassolov, V. A.; Ratner, M. A.; Pople, J. A.; Redfern, P. C.; Curtiss, L. A. J. Comput. Chem. 2001, 22, 976. (59) Bladeau, J.-P.; McGrath, M. P.; Curtiss, L. A.; Radom, L. J. Chem. Phys. 1997, 107, 5016. (60) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (61) Hesselmann, A.; Jansen, G. Chem. Phys. Lett. 2002, 357, 464. (62) Hesselmann, A.; Jansen, G. Chem. Phys. Lett. 2002, 362, 319. (63) Hesselmann, A.; Jansen, G. Chem. Phys. Lett. 2003, 367, 778. (64) Misquitta, A. J.; Szalewicz, K. Chem. Phys. Lett. 2002, 357, 301. (65) Misquitta, A. J.; Jeziorski, B.; Szalewicz, K. Phys. Rev. Lett. 2003, 91, 33201. (66) Gr€uning, M.; Gritsenko, O. V.; van Gisbergen, S. J. A.; Baerends, J. E. J. Chem. Phys. 2001, 114, 652. (67) Lias, S. G.; Levin, R. D.; Kafafi, S. A. Ion Energetics Data. In NIST Chemistry Webbook, NIST Standard Reference Database Number 69; Linstrom, P.J.; Mallard, W.G., Eds.; National Institute of Standards and Technology: Gaithersburg, MD, http://webbook.nist.gob. (68) Werner, H.-J.; Knowles, P. J.; Lindh, R.; Manby, F. R.; Sch€utz, M.; Celani, P.; Korona, T.; Rauhut, G.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Hampel, C.; Hetzer, G.; Lloyd, A. W.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; Palmieri, P.; Pitzer, R.; Schumann, U.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T. MOLPRO, a package of ab initio programs; 2006; see http://www.molpro.net. (69) Piquemal, J. P.; Marquez, A.; Parisel, O.; Giessner-Prettre, C. J. Comput. Chem. 2005, 26, 1052. (70) Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 4177. (71) Cubero, E.; Luque, F. J.; Orozco, M. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 5976.
ngy
an, J. G.; Soteras Guti
errez, I.; Luque, F. J.; (72) dehez, F.; A Schulten, K.; Chipot, C. J. Chem. Theory Comput. 2007, 3, 1914. (73) Gresh, N.; Pullman, B. Biochim. Biophys. Acta 1980, 625, 356. (74) Corral, I.; M
o, O.; Y
a~ nez, M.; Radom, L. J. Phys. Chem. A 2005, 109, 6735. (75) Garmer, D. R.; Gresh, N. J. Am. Chem. Soc. 1994, 116, 3556. (76) Gresh, N.; Garmer, D. J. J. Comput. Chem. 1996, 17, 1481. (77) Soteras, I.; Orozco, M.; Luque, F. J. Phys. Chem. Chem. Phys. 2008, 10, 2616. (78) Gresh, N.; Cisneros, G. A.; Darden, T. A.; Piquemal, J.-P. J. Chem. Theory Comput. 2007, 3, 1960.

11343

dx.doi.org/10.1021/jp205101z |J. Phys. Chem. A 2011, 115, 11331–11343