Does DFT-D estimate accurate energies for the binding of ligands to metal complexes?

We have studied the homolytic dissociation of a methyl radical from a model of methyl cobalamin. For this reaction, density functional theory with an atom-pairwise dispersion correction (DFT-D) gives a dispersion contribution to the bond dissociation energy (BDE) of 22–51 kJ/mol depending on the functional, i.e. much more than common estimates for the total dispersion interaction energy of the methyl group in typical solvents. We show that this large energy correction results from many rather small (0–2 kJ/mol) interactions that arise between the ligand and the metal and the other ligands when a short metal–ligand bond is formed. The energy terms result mostly from atom pairs connected by two or three bonds, i.e. terms that normally are ignored or scaled down at the molecular mechanics level, and have large contributions from r –8 terms. The added dispersion energy diminishes the variation in the calculated BDE observed among various generalised-gradient approximation functionals, whereas a gap still persists to the results of hybrid functionals. Model calculations at the local MP2 and CCSD (second-order perturbation theory and coupled cluster theory with single and double excitations) levels are in agreement with the dispersion interactions estimated by DFT-D (23–29 kJ/mol). However, both the DFT-D and the wavefunction-based results include middle-range correlation effects that vary much among different DFT methods owing to their different density-based description also in the short-range regime. Therefore, it is not meaningful to discuss which DFT method gives the most accurate estimate of the dispersion contribution to the BDE. Moreover, for a balanced treatment of dispersion during the binding reaction in solution, the dispersion energy of the ligand and the unbound complex with the surroundings needs also to be considered, which decreases the net dispersion contribution to binding by ~20 kJ/mol.


Introduction
 , CO, O 2 , or NO to metal corrins or porphyrins 23,24 or P t Bu 3 to Pd(0) complexes. 25 However, the binding of ligands to metals involve the formation of a short metal-ligand bond; therefore, it becomes difficult to define and quantify dispersion which is normally considered as a long-range correlation effect. Moreover, short-range (static) correlation effects become important, which makes results sensitive to the choice of the approximate exchange-correlation functional in DFT.
Considering the importance of metal-ligand binding in biochemical and organometallic reactions and the size of this effect, it is of highest interest to assess the applicability of the DFT-D approach to the formation of metal-ligand bonds and to decide whether it gives accurate estimates of the dispersion energy in such situations. Moreover, we will discuss how such reactions are best studied and how they should be coupled to continuum-solvation methods to make the reactions balanced.

Model systems
In this paper, we study the methyl cobalamin (related to coenzyme B12) system shown in Figure 1. In particular, we study the homolytic dissociation reaction: Co III CorImMe  Co II CorIm + Me  (1), where Cor is a corrin ring, Im is imidazole, and Me is a methyl group (CH3). This involves the homolytic cleavage of a Co-C bond, giving rise to Co II and a methyl radical.
In some cases, we have also studied a simpler model system, in which the corrin ring is replaced by four ammonia groups and Im by a fifth ammonia group: Co III (NH 3 ) 5 (Me)  Co II (NH 3 ) 5 + Me  (2) For this reaction, we assume that all reactants have C s symmetry (which affects the reaction energy by less than 1 kJ/mol), whereas for the reaction in Eqn. 1, no symmetry restraints were imposed.

DFT calculations
All DFT calculations were performed with the Turbomole 6.1 software 26 and the def2-SV(P) basis sets. 27 Most calculations were performed with the Becke-1988- Perdew-1986 functional (BP86), 28,29 but several other functionals were also tested: PBE, 30 TPSS, 31 BLYP, 28,32 B97, 33 TPSSH, 34 B3LYP, 35 PBE0, 36 and BHLYP. 37 Calculations with semi-local functionals were sped up by expansion of the Coulomb interactions in auxiliary basis sets, the resolution-of-identity (RI) approximation. 38,39 The structures were optimised until the change in energy between two iterations was below 2.6 J/mol (10 -6 a.u.) and the norm of the maximum norm of the gradients was below 10 -3 a.u. All open-shell systems were treated with the unrestricted formalism.
DFT-D2 8 calculations were performed with Turbomole, whereas single-point DFT-D3 10 calculations were performed with the stand-alone dftd3 program. 17 Single-point calculations with the polarised continuum model (PCM) 40,41 were performed with the Gaussian-03 software, 42 using the integral-equation formalism (IEF-PCM 43 ) and the UAKS radii (united atom topological model for Kohn-Sham theory). The calculations were performed at the BP86/def2-SV(P) level of theory (but in variance to the electrostatic solvation energy, the non-polar solvation energy terms, cavitation, dispersion, and repulsion, are independent of the level of theory or basis set). As a solvent, we used either water (with built in parameters) or ethylene glycol with the following parameters: dielectric constant = 40.245, but 2.05 at infinite frequency, radius = 2.51 Å, and density 0.0108 Å -3 . However, the latter solvent is not thoroughly parametrised and it was not possible to specify all the parameters needed for PCM. Moreover, the molecule is less spherical than water (making the solvent radius ill-defined) and it contains several atoms, so it is less well-defined what the dispersion and repulsion parameters should be for the continuum solvent (many popular MM water models, e.g. TIP3P, 44 make van der Waals interactions with only the oxygen atom). Therefore, only the results with water are discussed in the results section.

LMP2 and LCCSD calculations
Local correlation calculations, namely local second-order Møller-Plesset perturbation theory (LMP2) 45 and local coupled-cluster singles and doubles (LCCSD) 46 were carried out on the Co(NH3)5(Me) 2+ complex in its closed-shell singlet state. The def2-TZVPP basis set 27 was used throughout. Density-fitting (DF) approximations 47,48 were used for both the Hartree-Fock and the correlation part, with the def2-TZVPP/JKFIT 49 and def2-TZVPP/MP2FIT 50 auxiliary basis sets, respectively. Localized Pipek-Mezey orbitals 51 were used and the orbital domains were determined according to a natural population analysis occupation threshold of TNPA=0.03. 52 All calculations were carried out with a development version of Molpro 2010.2. 53

Result and Discussion
In this paper we study the homolytic cleavage of the Co(III)-C bond for the standard CoCorImMe model of methyl cobalamin, shown in Figure 1, i.e. the reaction in Eqn. 1. This reaction has been thoroughly studied with a variety of DFT methods and it has been shown that the homolytic Co-C bond dissociation energy (BDE; i.e. the energy of the reaction in Eqn. 1) strongly depends on the functional used. 54,55,56,57,58,59 In particular, dispersionuncorrected semi-local functionals like BP86 and BLYP give BDEs (140-160 kJ/mol) that are closer to the experimental result (15513 kJ/mol in ethylene glycol 60,61 ) than do hybrid functionals like B3LYP and B3P86 (100-120 kJ/mol). 56 The B1LYP functional gives an even lower BDE of 85 kJ/mol, illustrating that the BDE strongly depends on the amount of exact exchange in the functional. The effects of improved basis, treatment of the basis-set superposition error, relativity, solvation, zero-point vibrational energy, and thermal corrections have been considered, but do not change the results qualitatively. 56 However, Siegbahn and coworkers recently showed that with the DFT-D2 method in combination with the B3LYP functional, the BDE increases by 47 kJ/mol, 23 which would bring the results of hybrid methods closer to the experimental result. At a first glance, such a results seems spurious, considering that the estimated total dispersion energy of a methyl radical in solvents like cyclohexane, heptane, benzene, and water is 17-20 kJ/mol (i.e. the total dispersion interaction of a methyl radical with the surrounding solvent to infinity according to PCM calculations).
In this paper, we will analyse and explain how the dispersion energy can be so large in this complex and examine whether this is reasonable or not by comparison to wave functionbased calculations on smaller models. Moreover, we will discuss how balanced reaction energies can be obtained for reactions of this type in solution and enzymes.

DFT-D2 results
First, we calculated the BDE energy for the reaction in Eqn. 1 at the BP86/def2-SV(P) level. The result without any dispersion correction was 158 kJ/mol, in good agreement with previous results (156 kJ/mol with a different basis 56 ). With a single-point DFT-D2 correction, the BDE increased by 47 kJ/mol, i.e. the same value observed by Siegbahn et al. (at the B3LYP level). 23 With full geometry optimisation using BP86 and DFT-D2, the BDE remained at 205 kJ/mol. The dispersion energy can be easily understood, because it is a simple sum of r -6 terms for all pairs of atoms in the molecule: where r ij is the distance between the atoms, C ij is the dispersion coefficient for the two atoms i and j, s is a functional-dependent scaling factor, and f dmp is a damping function. This is very similar to the dispersion part of the standard Lennard-Jones term used in most molecular mechanics (MM) force fields. The only difference is that in MM methods, this term is omitted or scaled down for atoms that are connected by 1-3 bonds. Such information is not available in QM calculations; instead the damping function is included, which gives a similar effect.
Damping is necessary because of the singularity in the dispersion potential, which results from a break-down of the perturbation expansion on which Eqn. 3 is based (for a recent discussion of the effect of the damping function see ref. 62 ). Thus, DFT-D can be seen as a simple QM/MM method. To understand the large effect of dispersion on the BDE, we divided the dispersion energy into contributions from the various parts of the complex, viz. the methyl group (Me), Co, Im, and the Cor ring. The results in Table 1 show that the main contribution comes from interactions between the methyl group and the Cor ring (35 kJ/mol), but there are also large contributions from the interactions between the methyl group and Co (7 kJ/mol), as well as from the interactions between the imidazole group and the Cor ring (5 kJ/mol). The latter is a relaxation effect (i.e. it depends on the change in geometry between the Co(III) and Co(II) states) and may at first seem counter-intuitive, considering that the Co-NIm bond length is actually longer in CoCorImMe (2.17 Å) than in CoCorIm (2.13 Å). However, this is more than compensated by the fact that the Co ion resides essentially in the corrin plane in CoCorImMe, whereas it moves 0.15 Å out of it towards the imidazole ligand in CoCorIm. Thereby, the imidazole group actually comes closer to the corrin plane in the CoCorImMe complex. The other contributions are 1 kJ/mol or less.
Further division of the major terms shows that half of the Me-Cor interactions stem from C in Me (17 kJ/mol), whereas the other half is nearly equally divided between the three H atoms of Me (6 kJ/mol each). On the other hand, all the Me-Co interaction are attributed to the H atoms, giving 2 kJ/mol each. These three interactions are also the largest individual terms involving the methyl group. Thus, the individual Me-Cor interactions are small (up to -1.3 kJ/mol) but numerous. As can be seen in Table 2, they are dominated by interactions with the corrin NCor and C CCA atoms (i.e. the C atoms next to the N Cor atoms; 26 kJ/mol).
Going back to the MM perspective, we can conclude that all contributions from Me-Me (i.e. 1-3 interactions between the three H Me atoms in the methyl group) and Me-Co interactions, as well as the C Me -N Im/Cor and N Im -N Cor (C Me is the C atom in the methyl group and N Im the N atom in imidazole) interactions represent 1-3 interactions (interactions between atoms separated by two bonds) that normally are ignored in MM calculations. They sum up to a BDE contribution of 11 kJ/mol, i.e. 24% of the total dispersion contribution to the BDE. Moreover, the H Me -N Im/Cor , C Me -C ICA , C Me -C CCA , N Im -C CCA , and C ICA -N Cor (C ICA is the two C atoms in imidazole next to N Im ) interactions represent 1-4 interactions between atoms separated by three bonds, which normally are scaled down in MM. Those interactions sum up to 17 kJ/mol, i.e. 36% of the total dispersion contribution to the BDE. Therefore, 60% of the dispersion energy stems from interactions that are not included in a typical force field. However, the remaining contribution of 19 kJ/mol is still large and should not be problematic because it is long-ranged. Moreover, the scale factor for the 1-4 interactions is typically quite close to 1, e.g. 0.83 in the Amber force fields, 63 indicating that about 14 kJ/mol of the 1-4 interaction energy may also be reasonable (i.e. 33 kJ/mol in total). This analysis allows us to understand the source of the large dispersion energy for the binding of a methyl group to cobalamin in DFT-D. The Co-C bond is short, 1.98 Å, which brings the methyl group in close contact to many atoms in the metal complex. In particular, the methyl group is close to the many atoms in the corrin ring, with four N Cor atoms within ~2.7 Å, eight C CCA atoms within 3.3-3.7 Å, and three C meta atoms with distances of 3.7-4.0 Å. All these inter-atomic distances are smaller than typical van der Waals contacts, e.g. what can be estimated for a methyl radical in solution. The dispersion effect is further enhanced by the fact that the geometry of complex hardly changes by the binding of the methyl group (i.e. that the binding site is vacant before binding), so that there is little change in the dispersion energy of the other groups upon binding (for complexes with more flexible coordination sphere, it is likely that the ligands come closer when one ligand dissociates, giving a larger dispersion energy in the dissociated complex). On the contrary, the only significant contribution from the other ligands in this octahedral Co complex is an increase in the BDE energy, owing to the movement of the Co ion towards the imidazole ligand in the five-coordinate state.

DFT-D3 calculations
Next, we performed single-point energy calculations with the more recent DFT-D3 method. 10 It is more accurate, in particular for metal-containing systems, but also more involved, including both r -6 and r -8 terms, of which the latter can be expected to be important for such short-range interactions observed in the CoCorImMe complex. The DFT-D3 dispersion correction to the BDE for the BP86 functional is 6 kJ/mol smaller than the corresponding DFT-D2 correction, 41 kJ/mol. Interestingly, it is dominated by the r -8 term, 27 kJ/mol, which is almost twice as large as the r -6 term, 14 kJ/mol. This indicates that the dispersion is dominated by short-or medium-range correlations for which the definition of the dispersion energy becomes less clear.
Next, we calculated the dispersion correction to the BDE with eight additional density functionals. As mentioned above, previous studies have shown that the calculated BDE depends strongly on the functional used, and especially the amount of exact exchange. 56 The DFT-D3 dispersion correction also depends on the functional. It would be satisfactory if the variation in the calculated BDEs decreases when the dispersion correction is added. This would indicate that the missing dispersion or middle-range correlation, described by the r -8 term, is the main cause of the difference between the functionals.
As can be seen in Table 3, this is indeed the case in the group of semi-local generalisedgradient approximation (GGA) functionals. Without the dispersion correction, the BDE varies between 128 kJ/mol (B97) and 173 kJ/mol (PBE). With D3-correction the span is reduced from 45 kJ/mol to 20 kJ/mol, and the calculated BDEs follow the order PBE (199 kJ/mol) > BP86 > TPSS > BLYP > B97 (179 kJ/mol). Thus, there is a good anti-correlation between the BDE without dispersion and the dispersion correction (r 2 = 0.81).
However, if the hybrid functionals are included in the comparison, the correlation completely disappears (r 2 = 0.01). As mentioned above, the calculated BDE is lower for all the hybrid functionals except TPSSH and follows roughly the amount of non-local Fock exchange (50% in BHLYP, 25% in PBE0, 20% in B3LYP, and 10% in TPSS), with the only exception of PBE0 and B3LYP, where the effect of the underlying functional reverses the rather small difference in the amount of Fock exchange. However, the DFT-D3 dispersion correction follows another order: B97 (51 kJ/mol) > BP86 > BLYP > B3LYP > TPSS > TPSSH > BHLYP > PBE0 > PBE (22 kJ/mol), showing little difference between the pure and hybrid functionals or the amount of Fock exchange. Therefore, the gap of the BDE between the GGA and hybrid functionals is not reduced by the D3-correction, and for the hybrids, the range of BDEs for the functionals is not improved at the DFT-D3 level, as opposed to the case of the GGAs.

Wavefunction calculations
Now, we come to the most important question, viz. to what extent we can trust the DFT-D3 results quantitatively. Unfortunately, this is hard to answer conclusively for such a relatively large and complicated system. For normal non-bonded systems, the dispersion energy can be estimated as the difference in Hartree-Fock (HF) and MP2 (or even better coupled-cluster theory) interaction energies. However, for transition-metal complexes, HF gives very poor results (for our CoCorImMe model, the BDE at the HF level is -178 kJ/mol, for example 56 ) and the performance of MP2 can also be rather poor, typically worse than DFT, owing to the presence of extensive static correlation and orbital-relaxation effects. Therefore, more advanced (e.g. multiconfigurational) levels of theory are needed, 64 but they are normally too expensive to use for such complicated complexes and the choice of the active space is problematic. Likewise, it is questionable if other methods to calculate various contributions to the interaction energy, e.g. the symmetry adapted perturbation theory, 65 work properly for complexes involving the formation of a short metal-ligand bond.
Therefore, we decided to use local correlation methods. First, we employed the DF-LMP2 method to calculate the dispersion contribution to the BDE energy. In order to avoid dealing with strong static correlation effects and to reduce the system size, we made use of a simpler model, viz. Co(NH 3 ) 5 Me 2+ . Local correlation methods, besides having a lower computational cost than their canonical counterparts, also allow for a decomposition of the correlation energy: Because both occupied and virtual spaces are local, one can divide the excitations into different classes according to the orbitals involved. 66 If one considers dispersion interactions, the natural physical picture will be that of an excitation involving two electrons, each in an occupied orbital of a different interacting fragment, being excited within the same fragment. A more detailed discussion of the different excitation classes can be found in Ref. 66. In the problem at hand, it is difficult to exactly partition the occupied space into two distinct fragments because the methyl group shares a bond with the metal centre (the virtual space is not problematic because each virtual orbital is tagged to a single atomic centre). We circumvented the problem through the use of two separate calculations that should give welldefined upper and lower bounds to the dispersion energy. In a first calculation, the Co-C orbital is considered to belong exclusively to the methyl fragment. All dispersion-type excitations are included, except those involving the bonding orbital and the Co orbitals. This value gives an upper bound to the dispersion interactions, because the orbital does not belong exclusively to the methyl. In a second calculation, the bonding orbital is neglected and only dispersion interactions involving the three C-H orbitals are included. This should give a lower bound and uses a similar definition to that used in a previous study of bridged dimers. 67 With the def2-TZVPP basis set, we obtained dispersion energies of 25-29 kJ/mol. This is slightly lower than what is estimated by DFT-D2 and DFT-D3 for the BDE reaction in Eqn. 2 for the BP86 functional, 38 and 30 kJ/mol, respectively (Table 4). Next, we repeated the same calculations with the DF-LCCSD method. We obtained slightly lower dispersion energies of 23-26 kJ/mol. This is expected, because MP2 normally overestimates dispersion, whereas CCSD underestimates it (owing to the lack of higher-order excitations). On the other hand, the basis set used is rather small, which often cancels the MP2 overestimation. Therefore, we can estimate that the true dispersion for this complex is around 27 kJ/mol. However, as can be seen in Table 4 and in analogy with the results for the full model in Table 3, the DFT-D correction to the BDE strongly depends on the density functional used (15-46 kJ/mol for DFT-D3 and 21-45 kJ/mol for DFT-D2). As noted above, a significant portion of the dispersion effect is relatively short-ranged. Because the electronic part of the various density functionals employed performs differently in this regime (show a rather different repulsive behaviour), the DFT-D correction also has to correct for deficiencies of the density-based treatment. This view is also supported by the large contribution of the r -8 term to the BDE. Based on these results, it seems clear that the short-range part of the DFT-D dispersion energy is a model-dependent quantity that should not be overinterpreted. Only the (non-empirical) long-range part, which is more or less correct by construction, can be directly compared to long-range wavefunction-based results. 10 In passing, we note that a similar ambiguity exists in the interpretation of the MP2 or CCSD(T) calculations, which are based on an uncorrelated reference point (Hartree-Fock) and attribute all correlation effects as being of dispersive nature. This is seen for example in DF-LMP2 calculations on the full model system without the central Co ion (with an overall charge of -2 and maintaining the complex geometry). This removes the coordination of the methyl, providing a well-defined separation between the different orbital spaces. The energy separation scheme used above shows that only 46% of the correlation effect (of the interactions involving the methyl group) are actually of dispersive nature. With these aspects of the problem in mind, we consider the agreement between the employed wave function and DFT methods as being quite satisfactory.

Reactions in solution
Finally, we point out that any reasonable electronic-structure method (e.g. DFT-D or CCSD(T)) takes into account only the new dispersion energy between the complex (CoCorIm in our case) and the ligand (Me), but not any interactions with the surroundings before and after binding. This is completely valid if the reaction takes place in a vacuum. However, most chemical reactions take place in solution and then this assumption becomes problematic, because the treatment of dispersion will be unbalanced, leading to a too large BDE (in the complex, dispersion interactions between CoCorIm and Me are included, but before binding these groups will have dispersive interactions with the solvent instead, which are ignored). Although this quenching effect of isolated molecule dispersion seems to be trivial, it is often overlooked in applications and hence considered in our example in more detail.
In QM calculations of biological reactions, it is common to estimate the effect of the surrounding solvent from some sort of continuum-solvation model, e.g. the polarised continuum model (PCM) 68,69 or the related COSMO method, 70 by solving the Poisson-Boltzmann equation, 71,72 or by the generalised Born approach. 73,74 These methods primarily calculate the electrostatic part of the solvation energy. 75 However, to obtain accurate estimates of the total solvation energy, also the non-polar part of the solvation energy needs to be included, viz. the cost of forming a cavity in the solvent and the dispersion and repulsion energy between the solvent and the solute. 75 Unfortunately, it is seldom specified whether the non-polar part is included in the solvation energy or not. For example, in the calculations of Siegbahn et al., 23 only the polar part was included (Siegbahn, pers. commun.), making the treatment of dispersion unbalanced when comparisons to solution-phase experimental data are made.
For a reaction in a homogeneous solvent, the dispersion energy for all reactants can be calculated by the PCM method. This would make the treatment of dispersion balanced in the reaction and the dispersion energy is treated at a similar level of approximation as the DFT-D energy, i.e. by MM r -6 terms. 76 In particular, the PCM dispersion term is a pure energy term (not a free energy), calculated as the volume integral of the dispersion energy of the solute (treated as MM atoms) over the space not occupied by the solute and assuming a constant density of the solvent. This is probably the best approximation possible without any knowledge of the detailed solvent structure around the solute.
The PCM repulsion energy is calculated in the same way, i.e. as volume integrals over the Lennard-Jones r -12 terms of the QM atoms, treated in a MM manner and it is also a pure energy term. 76 It should be included if we aim at reproducing the reaction in solution, including all solvent effects.
The last PCM term, the cavitation energy, is more problematic, because it is a free-energy term, including also the entropy cost of forming the cavity and reorganising the solvent around it. It is fitted to reproduce the total solvation free energy of non-polar compounds using a function of the radii of each of the solute atoms to the powers of 1, 2, and 3 (i.e. including both area and volume terms). 77 On the other hand, it should be noted that neither of the four solvation terms (electrostatic, cavitation, dispersion, and repulsion) is experimentally observable alone, only the total solvation free energy.
For an enzyme reaction, the use of the non-polar PCM energy terms is less evident. In particular, it becomes questionable to include the cavitation energy if the reaction takes place in an active site that is buried inside the enzyme and is dehydrated before ligand binding. Then, the cavitation cost is paid already when the protein folds and it does not change during the reaction. However, the other two terms remain valid, although it is no longer evident what solvent should be used, because the surroundings is no longer homogeneous. A much better approach is to use QM/MM calculations, in which these terms are included in an atomistic (but not always in a thermodynamically averaged) way. However, if the ligand is taken from solution, it is once again evident that the best thing is to include all PCM terms for the pure solvent.
Going back to our Co-C bond cleavage in Eqn. 1, we obtain the PCM energies listed in Table 5. It can be seen that the electrostatic (polar) contribution to the BDE energy is only -2 kJ/mol. This is in agreement with previous results. 56 Thus, the polar solvation energy has little effect on the BDE.
On the other hand, the dispersion contribution to the BDE is -22 kJ/mol. Thus, if we make a balanced treatment for all three reactants, the dispersion energy decreases to 25 and 19 kJ/mol for DFT-D2 and DFT-D3 at the BP86 level, respectively. This shows that it is mandatory to consider the dispersion of all reactants, e.g. by using PCM solvation. The dispersion term comes almost entirely from the methyl group, whereas the contributions from the metal complex essentially cancel between the CoCorImMe and CoCorIm states. This also applies to the repulsion term, which is smaller than the dispersion, 6 kJ/mol, and has the opposite sign, as expected.
On the other hand, the cavitation energy is large and positive, 20 kJ/mol. The repulsion and cavity terms together more than counteract the dispersion term so that the free energy of the reaction in Eqn. 1 actually is disfavoured by solvation effects, but only by 2 kJ/mol. The reason for this is that the methyl radical is non-polar and has a low solvation in water -PCM predicts a hydration free energy of 7 kJ/mol, close to that of methane (8 kJ/mol 78 ). However, this does not necessarily affect the BDE, which is the enthalpy (not free energy) of the reaction in Eqn. 1. Unfortunately, PCM does not provide any estimate of the division of the cavitation energy in enthalpy and entropy terms, so we cannot obtain any estimate of the effect of the total solvation energy on the BDE. Therefore, it is not yet possible to conclusively decide which functional gives the best estimate of the BDE for methyl cobalamin.
A solution to this dilemma could be to consider instead the free energy of the Co-C dissociation. However, then we will encounter another major theoretical problem, viz. the change in translational and rotational entropy during ligand binding. The entropy in reactions treated by QM methods is normally estimated from the harmonic frequencies calculated with QM, using standard statistical mechanics methods and a rigid-rotor ideal-gas harmonicoscillator approximation. 79 This means that the translational entropy is taken from the Sackur-Tetrode equation, giving a penalty of ~50 kJ/mol when a ligand binds to a host. This fact has been much discussed, especially in the ligand-binding community. 80,81,82,83,84 There is a consensus that this is an overestimate, owing to the ideal-gas assumption (a molecule has less freedom to translate and rotate in solution than in gas phase), but there is still no agreement by how much. It has been suggested that this energy should halved or reduced by ~32 kJ/mol, 82 but these estimates are very uncertain and it is not known how general they are. Therefore, it has to be accepted that theoretical estimates of binding free energies have an uncertainty of ~30 kJ/mol.

Conclusions
In this paper, we examine the recent suggestion 23,24,25 that the binding of small ligands to metal complexes involves a dispersion component of 30-60 kJ/mol, i.e. much more than the total dispersion energy of the same ligand in solution. As a typical example, we have chosen the homolytic cleavage of the Co-C bond in the methyl cobalamin model CoCorImMe.
We show that the large dispersion energy is caused by the formation of a short metalligand bond that brings the ligand in close contact with many atoms of the other metal ligands, with distances that are shorter than for the free ligand in solution. Thus, the large dispersion energy is the sum of many rather small terms (<2.3 kJ/mol). However, 24% of the energy comes from atoms that are two bonds apart, i.e. interactions that normally are ignored at the MM level, and 36% of the energy comes from atoms that are three bonds apart, i.e. interactions that normally are scaled down at the MM level. This shows that these energies are sensitive to the damping function and that they mix in middle-range correlation effects. This is supported by the fact that the dispersion energy is dominated by the r -8 term in the DFT-D3 approach and that the DFT-D energies calculated with different functionals vary much, e.g. from 22 to 51 kJ/mol for the CoCorImMe BDE.
Nonetheless, DF-LMP2 and DF-LCCSD calculations indicate that the DFT-D methods provide quite accurate estimates of the dispersion energy for this type of complexes, e.g. 23-29 kJ/mol for the reaction in Eqn. 2, compared to 15-46 kJ/mol for DFT-D3 with nine different functionals. Because short-range dispersion is a model-dependent and not very welldefined quantity, we recommend that results of DFT-D, as well as of local wavefunction calculations, should not be overinterpreted in this context. Moreover, we also show that if the reaction does not take place in a vacuum, a direct application dispersion-including electronic structure methods alone will overestimate the effect of dispersion for the binding reaction because the dispersion energy of the unbound ligand and complex is ignored. We show that these energy terms can be estimated by the PCM method and that they will reduce the net dispersion energy by 22 kJ/mol.
Strictly speaking, all PCM solvation energies should be included for a reaction in a homogeneous solvent (i.e. also the electrostatic, repulsion, and cavitation parts). However, it should be noted that the electrostatic and cavitation parts are free energies, containing entropic effects also. For enzyme reactions, it is better to use QM/MM methods to obtain atomistic detail and balanced reactions (but still PCM for the free ligand). Finally, we have also pointed out the problem of estimating the translational entropy of the free ligand, which may introduce an uncertainty of ~30 kJ/mol.
In conclusion, we show that DFT-D3 reasonably accounts for changes of dispersion-type inter-and intra-molecular correlation energies when ligands bind to metal complexes, but the reaction needs to be balanced by inter-molecular dispersion terms also for the unbound reactants in solution.  Table 2. Contributions to the dominant Me-Cor BDE dispersion term from the various atoms in the methyl group (CMe and the three HMe atoms) and the corrin ring (the four NCCA atoms, the eight CCCA atoms next to the N atoms, the three meta Cmeta atoms that connect the fiverings, the remaining eight CB atoms, the two HA atoms bound to CA atoms, the three HM atoms bound to the CM atoms, and the 16 HB atoms bound to the CB atoms). The atom names are shown in Figure 1. The results are based on single-point DFT-D2 energy calculations on the BP86/def2-SV(P) structure of CoCorImMe.     1. The structure of Co III CorImMe optimised at the BP86/def2-SV(P) level with atom names marked.

TOC Graphics
With local correlated and continuum-solvation methods we discuss the large DFT-D dispersion contributions to metal binding. Dispersion Energy (kJ/mol)