ARTICLE pubs.acs.org/JPCA
Which DFT Functional Performs Well in the Calculation of Methylcobalamin? Comparison of the B3LYP and BP86 Functionals and Evaluation of the Impact of Empirical Dispersion Correction Hajime Hirao* Division of Chemistry and Biological Chemistry, School of Physical and Mathematical Sciences, Nanyang Technological University, 21 Nanyang Link, Singapore 637371
bS Supporting Information ABSTRACT: Controversy remains regarding the suitable density functionals for the calculation of vitamin B12 systems that contain cobalt. To identify the optimum functionals, geometry optimization calculations were performed on a full-size model of methylcobalamin (MeCbl) using the B3LYP, B3LYP-D, BP86, and BP86-D methods in conjunction with the 6-31G* basis set. Single-point energy evaluations were also performed with the 6-311+G(2d,p) basis set. Consistent with previous studies, the BP86-optimized geometry showed fairly good agreement with the experimental geometry. Various factors that may influence the homolytic bond dissociation energy (BDE) of the CoC bond of MeCbl were systematically evaluated with these methods. Our analysis demonstrated that dispersion was the largest correction term that influenced the magnitude of BDE. Previous studies have shown that B3LYP significantly underestimates BDE, whereas BP86 gives BDE values that are fairly close to the experimental values (3637 kcal/mol). The same trend in the relative magnitudes of the BDEs was observed in the present calculations. However, BP86 underestimated the BDE for a full model of MeCbl. When the amount of HartreeFock exchange in the B3LYP functional was reduced to 15% and the dispersion correction was made (i.e., B3LYP*-D), the calculated BDE was in good accord with experimental values. B3P86-D also performed well. A detailed analysis was undertaken to determine which atoms in cobalamin have large dispersion interactions with a methyl fragment of MeCbl.
1. INTRODUCTION Methylcobalamin (MeCbl) is a form of vitamin B12.14 In its isolated state, MeCbl has a coordination sphere in which the central cobalt ion is coordinated to a substituted corrin macrocycle, dimethylbenzimidazole (DMB) of a peripheral substituent tethered to ring D of the macrocycle, and a methyl group (Scheme 1).5 Several biologically important methyltransferase enzymes, such as methionine synthase (MetH), acetyl-CoA synthase, and methanol:coenzyme M methyltransferase,2b are known to utilize MeCbl or its analogue as a cofactor for their catalytic reactions. These reactions involve methyl transfer steps, in which MeCbl acts as a methyl donor. For example, one of the two major steps of the MetH reaction involves the transfer of a methyl cation from MeCbl to homocysteine to produce methionine. This process is essentially regarded as an SN2 reaction, and the CoC bond cleavage is therefore viewed as heterolysis. This differs from the reaction of another form of vitamin B12, adenosylcobalamin (AdoCbl), in which the CoC bond is homolytically cleaved to trigger subsequent reactions.3 The homolytic cleavage of the CoC bond of MeCbl has also attracted considerable attention, because it poses significant challenges for computational studies. In fact, much effort has r 2011 American Chemical Society
been devoted to the calculation of the bond dissociation energy (BDE) for the following process: MeCbl f • CH3 þ Cbl
ð1Þ
in which the CoC bond of MeCbl is homolytically cleaved to generate a methyl radical (•CH3) and pentacoordinated cobalamin, denoted simply as Cbl. Upon homolytic cleavage, the formal oxidation state of Co changes from +3 to +2. A number of theoretical calculations, mostly based on density functional theory (DFT) methods, have been performed to compute BDE.612 Carefully designed experiments performed by Martin and Finke13 and Hung and Grabowski14 found the BDE values to be 37 ( 3 and 36 ( 4 kcal/mol, respectively, and these values have been used as references for the assessment of various computational methods. Interestingly, the importance of the homolytic nature of the CoC cleavage in the reaction of MetH has recently been pointed out by Alfonso-Prieto et al.15a and Kumar et al.15b They proposed the possible involvement of Received: June 6, 2011 Revised: July 15, 2011 Published: August 02, 2011 9308
dx.doi.org/10.1021/jp2052807 | J. Phys. Chem. A 2011, 115, 9308–9313
The Journal of Physical Chemistry A Scheme 1. Methylcobalamin
CoC homolysis of MeCbl in the methionine formation step, which is assisted by the preceding electron transfer from the deprotonated homocysteine to a corrin π* orbital. Although this enzymatic CoC bond cleavage may not be as simple as shown in eq 1, the homolytic character of the cleavage is now gaining increasing biological relevance. In the majority of cases, past computational studies of the CoC bond cleavage have been performed with truncated models, typically using a corrin ring without any substituent. Simplified axial ligands, such as imidazole or benzimidazole have also been used to substitute for the DMB ligand. Thus, calculations have been made by omitting the acetamide, propanamide, and methyl substituents of corrin, as well as the linker connecting the corrin and DMB parts (Scheme 1). Although such simplification has been unavoidable in practice to make feasible computations of the large cofactor, its impact on the computed CoC BDE is unclear. To reduce this kind of uncertainty, the use of a full-size model of MeCbl is desirable. To the best of our knowledge, only Rovira and Biarnes have performed a full quantum mechanical calculation of the CoC BDE for the entire MeCbl system using the CarParrinello molecular dynamics method,11 combined with the B86-exchange/ P86-correlation functionals, a plane wave basis set with a 70 Ry kinetic energy cutoff, and the MartinsTroullier pseudopotentials. They obtained geometrical parameters for MeCbl that well accord with experimental values, and also obtained a BDE value of 31.7 kcal/mol. D€olker et al. used the IMOMM method to perform calculations on a full model of MeCbl, applying the BP86/[LANL2DZ-P,6-31G(d)] and MM3 molecular mechanics methods in a hybrid manner.12 They obtained a BDE value of 36.0 kcal/mol. A key outcome of previous computational studies is that hybrid functionals, such as B3LYP tend to underestimate significantly the BDE of MeCbl. In contrast, generalized gradient approximation (GGA) functionals, such as BP86, have been shown to produce larger BDE values, which are in excellent
ARTICLE
agreement with the experimental BDEs. The smaller BDE obtained with hybrid functionals was attributed to the poor ability of the exact HartreeFock (HF) exchange (which is included in B3LYP) to deal with the distinctly different electronic states of the fragments in eq 1 in a balanced way.7 Systematic errors may arise because MeCbl on the left-hand side has a closed-shell character, whereas Cbl on the right-hand side has an open-shell character. This rationalization appears to explain quite well the relative magnitudes of the CoC BDEs calculated with hybrid and GGA functionals. Nevertheless, the tendency for B3LYP to underestimate BDE does not necessarily guarantee the accuracy of BP86. The possibility remains that the agreement of the BP86 and experimental BDEs may have been achieved accidentally, as a result of using truncated models or omitting key factors that could significantly influence the BDE. To identify the optimum functionals for the study of B12 cofactors, a systematic study is obviously required that uses a full-size model of MeCbl and includes various factors involved in the bond cleavage process. It is noteworthy that the good performance of BP86 has recently been questioned by Siegbahn et al.10 Applying the B3LYP*-D method to a small model of MeCbl, they pointed out the importance of the dispersion effect (which is not adequately evaluated by conventional DFT methods) in the calculation of the BDEs of alkylcobalamins. They showed that the dispersion effect is as large as 11.3 kcal/mol for MeCbl, accounting to some extent for the underestimated part of the BDE calculated with hybrid functionals.10a Importantly, the magnitude of the dispersion correction term was significantly larger than that of any other correction terms that have ever been examined.7 When used with the B3LYP* functional,16 the dispersion-corrected DFT method gave a BDE value of 32.0 kcal/mol for MeCbl, which is a significant improvement on the B3LYP value (16.2 kcal/mol). Chen et al. argued that, if the dispersion correction is made on the BP86-calculated BDE, the agreement of the BP86 value with the experimental value may no longer be good.10b In fact, the importance of the dispersion effect has recently been highlighted for other transition-metal complexes or enzymes.1719 Lonsdale et al. demonstrated that the inclusion of dispersion effects in the B3LYP-D method significantly stabilized the complex between P450cam compound I and the camphor substrate, relative to their separate states.17 They also showed that the barrier to H-abstraction is lowered by the inclusion of the dispersion effect. Interestingly, the dispersion effect significantly changed the geometry of the H-abstraction transition state of P450cam, in a manner that enhanced the interaction between camphor and the porphyrin ligand. Although the dispersion effect seems to be less significant when structurally similar states are compared,20 recent studies indicate that when the physical quantity of interest is associated with the difference between structurally different states, its impact can be substantial. These previous studies prompted us to reexamine the DFT and DFT-D methods of calculating MeCbl, for which no consensus has been reached regarding the optimum DFT functionals that provide the most reasonable results. For future reference, we wanted to determine which functional really performs well in the calculation of MeCbl. In this study, we focused particularly on the B3LYP and BP86 functionals, which are two representative functionals that have been extensively used, especially for the study of vitamin B12 systems. We used a full-size model of MeCbl to eliminate the potential error that may arise from the use of truncated models. The simplest approximation to BDE should 9309
dx.doi.org/10.1021/jp2052807 |J. Phys. Chem. A 2011, 115, 9308–9313
The Journal of Physical Chemistry A
ARTICLE
the energy difference of the fragments in eq 1. However, experimentally obtained BDE values will also contain other factors.7 Therefore, we decided to evaluate and compare several factors that may influence BDE. As far as we are aware, the effect of the dispersion correction on the geometry of MeCbl has not been investigated. Understanding this effect was another goal of this study.
2. METHODS The accurate crystal structure of MeCbl, as determined by Randaccio et al. using synchrotron diffraction data, was used as the initial structure for our DFT analysis.5c All DFT calculations were performed with Gaussian 09.21 DFT/6-31G* and DFT-D/ 6-31G* geometry optimizations were performed in the gas phase with the B3LYP and BP86 functionals;22,23 i.e., B3LYP, B3LYPD, BP86, and BP86-D were used for this purpose. Grimme’s protocol was used for the DFT-D calculations.24 The 6-311 +G(2d,p) basis set was used to improve energy evaluations of the geometries optimized with the 6-31G* basis set. For the sake of simplicity, we henceforth call the 6-31G* and 6-311+G(2d,p) basis sets B1 and B2, respectively. The numbers of basis functions for MeCbl were 1587 (B1) and 3076 (B2). In the DFT-D method,24 the dispersion-corrected total energy (EDFT-D) is expressed as EDFT-D ¼ EDFT þ Edisp
N1
N
ij
∑ ∑ 66 fdamp ðRijÞ i ¼ 1 j ¼ i þ 1 Rij C
where C6 coefficients relate to atomic ionization potentials and static dipole polarizabilities.24 The damping function is given by 1 1 þ
edðRij =Rr 1Þ
ð5Þ
where d is a constant (=20.0) and Rr is the sum of atomic vdW radii (R0). The same scale factor s6 (=1.05) was derived for B3LYP and BP86.24 To assess the impact of the amount of HF exchange in the B3LYP functional and the correlation functional on the relative magnitude of BDE, B3LYP*-D and B3P86-D were also examined. In the former, the HF-exchange of B3LYP was reduced to 15%. Because the current version of Gaussian 09 does not allow us to perform geometry optimization with these methods, single-point calculations were performed on top of the BP86-optimized geometries (i.e., B3LYP*-D//BP86 and B3P86-D//BP86). Moreover, the s6 value has not been derived specifically for B3LYP* and B3P86, so we assumed that the s6 values for these functionals are identical to that for B3LYP and BP86 (i.e., 1.05). The individual terms in eq 3 were also inspected using our originally developed code.
ð6Þ
The first term on the right-hand side of eq 6 is the energy difference of the fragments obtained with B1: ΔEB1 ¼ EB1 ðCblÞ þ EB1 ð• CH3 Þ EB1 ðMeCblÞ
ð7Þ
where the dispersion effect is excluded from the energy terms. The other terms are correction terms that were evaluated for the B1-optimized geometries. Δdisp is the dispersion correction for a given geometry and originates from the Edisp term in eq 2. The ΔZPE,B1 term is for zero-point vibrational energy (ZPE) correction, which is obtained from vibrational frequency calculations of the fragments. ΔH,B1 is an enthalpy correction term defined as ΔH, B1 ¼ ΔHB1 ðΔEB1 þ Δdisp þ ΔZPE, B1 Þ
ð8Þ
where ΔHB1 is the enthalpy difference of the fragments in eq 1 but contains the dispersion effect in DFT-D calculations. The temperature and pressure were set at 298.15 K and 1 atm, respectively. The basis-set effect, ΔBS,X, is ΔBS, X ¼ ΔEX ΔEB1
ð9Þ
which has a nonzero value only in calculating BDEB2. The correction for the basis-set-superposition error, ΔBSSE,X, is calculated as
ð3Þ
where s6 is a scale factor, C6ij is a dispersion coefficient, Rij is the distance between atoms i and j, and fdamp is a damping function that is necessary to avoid near-singularities for small Rij. The factor s6 is introduced because different functionals give different behaviors of the intermolecular potential. C6ij is defined as a geometric mean of C6 coefficients parametrized for atoms: qffiffiffiffiffiffiffiffiffiffi ij j C6 ¼ Ci6 C6 ð4Þ
fdamp ðRij Þ ¼
BDEX ¼ ΔEB1 þ Δdisp þ ΔZPE, B1 þ ΔH, B1 þ ΔBS, X þ ΔBSSE, X þ ΔsolvN, X þ ΔREL, X
ð2Þ
where EDFT is DFT energy and Edisp is dispersion energy. Edisp was calculated as Edisp ¼ s6
The geometries of all fragments in eq 1 (MeCbl, •CH3, and Cbl) were fully optimized. Bond dissociation energy BDEX, evaluated with basis set X (X d B1 or B2), were assumed to be defined as
ΔBSSE, X ¼ ΔEX, CP ΔEX
ð10Þ
where ΔEX,CP is the counterpoise (CP) corrected ΔEX, which is obtained by applying the CP method to the B1-optimized geometry of MeCbl.25 ΔsolvN,X is a solvation-correction term. Because the experiments were performed in ethylene glycol (denoted “solv1”, ε = 40.245)13 and in water (denoted “solv2”, ε = 78.3553),14 the effects of these solvents on BDE were evaluated by applying the default self-consistent reaction field (SCRF) method implemented in Gaussian 09.26 Finally, the relativistic effect (ΔREL,X), which sometimes plays a critical role in transition-metal systems,18 was evaluated by DouglasKroll Hess second-order scalar relativistic calculations.27 VMD was used to visualize the molecules.28
3. RESULTS AND DISCUSSION 3.1. Geometry of MeCbl. Table 1 summarizes the six key bond distances associated with the central cobalt ion for the optimized geometries of MeCbl. Overall, fairly good agreement between the theoretical and experimental geometries was observed. The best agreement was obtained with the BP86 functional, for which the average deviation of the six distances from the experimental ones was only 0.009 Å. In fact, the better performance of BP86 than B3LYP in predicting geometries is widely recognized not only for alkylcobalamin6,7 but also for other transition-metal complexes.29 The dispersion correction of BP86 (BP86-D) slightly shortened all key distances, and the average deviation became 0.015 Å larger than that for BP86. It follows that the dispersion correction of BP86 did not improve the accuracy of the geometry. The CoNDMB distance obtained with B3LYP is known to be too long.6,7 In the present calculation 9310
dx.doi.org/10.1021/jp2052807 |J. Phys. Chem. A 2011, 115, 9308–9313
The Journal of Physical Chemistry A
ARTICLE
Table 1. Comparison of Key Distances (Å) in Optimized Geometries
a
r(CoCMe)
r(CoNDMB)
r(CoNA)
r(CoNB)
r(CoNC)
r(CoND)
|Δr|ava
B3LYP
1.957
2.238
1.890
1.938
1.930
1.894
0.027
B3LYP-D
1.957
2.102
1.882
1.923
1.921
1.880
0.016
BP86
1.966
2.181
1.870
1.932
1.921
1.875
0.009
BP86-D
1.965
2.072
1.862
1.916
1.915
1.857
0.024
exptlb
1.979
2.162
1.877
1.922
1.918
1.874
Average deviation from the experimental values. b Taken from ref 5c.
Table 2. Comparison of Computed and Experimental BDEs ethylene glycol solvent method
water solvent
B1
B2
B1
B2
B3LYP B3LYP-D
20.0 35.3
16.8 31.9
20.0 35.2
16.7 31.8
BP86
32.6
30.4
32.6
30.3
BP86-D
46.8
44.1
46.7
44.0
B3LYP*-D//BP86
40.5
37.7
40.5
37.7
B3P86-D//BP86
41.3
38.7
41.2
exptl
37 ( 3a
38.6 36 ( 4b
a
Obtained by thermolysis in ethylene glycol.13 b Obtained by photoacoustic calorimetry in water.14
Figure 1. Superposition of the experimental (blue) and BP86-optimized (yellow) geometries. Only heavy atoms are shown.
based on a full model of MeCbl, the distance was calculated to be 2.238 Å, which is 0.076 Å longer than the experimental value. This was the largest of all the deviations considered for the distances involving Co. B3LYP-D reduced this deviation, although the extent of the reduction of the CoNDMB distance was slightly too large. Therefore, the resultant CoNDMB distance (2.102 Å) was somewhat shorter than the experimental value (2.162 Å). The distances between Co and equatorial nitrogen atoms (NAND) were longer than the experimental values when B3LYP was used, but this trend was mitigated with the use of B3LYP-D. As a consequence, the average deviation was reduced from 0.027 Å (B3LYP) to 0.016 Å (B3LYP-D) by including the dispersion term. In contrast to the case of BP86, the addition of the dispersion term improved the B3LYP geometry. Figure 1 shows the superposition of the experimental and BP86-optimized geometries. It can be seen that the atoms in the two geometries overlap one another very well, although some clear geometric deviations can be seen at the corrin-DMB linker and the propanamide substituent on ring C. It is likely that these deviations originate from interactions with neighboring molecules in the crystal, whereas this effect is absent in our isolated MeCbl model. 3.2. Bond Dissociation Energy. Table 2 summarizes the theoretical BDEs of MeCbl obtained with several different methods, together with the experimental values. Their energy
components (see eq 6) are listed in Table 3. Comparison of the values obtained for the ethylene glycol solution and the water solution shows that the solvent effect was very small, which is consistent with similar experimental BDEs obtained in different solvents. It was found that B1 gives BDEs that were a few kcal/ mol larger than the B2 BDEs. When we consider the functional dependence of BDE, B3LYP is seen to provide too small BDE values, which is consistent with previous studies. As noted by Siegbahn et al., the addition of the dispersion-correction term has a significant impact on the BDE value. In fact, the B3LYP-D/B2 value is about 15 kcal/mol larger than the B3LYP/B2 value. Within the B3LYP-D/B2 BDE, the dispersion correction term amounts to 13.3 kcal/mol, which has the largest magnitude of all the correction terms considered (Table 3). Although the BDE values obtained with B3LYP-D/B1 (35.235.3 kcal/mol) showed relatively good agreement with the experimental values, the value for B3LYP-D/B2//B1 was 31.831.9 kcal/mol, which is ∼5 kcal/mol smaller than the experimental values. The BP86calculated BDE value was also smaller than the experimental value and even smaller than the B3LYP-D value. This is not necessarily consistent with previous studies based on truncated models, which showed good agreement between the BP86 and experimental BDEs. The smaller BDE obtained here with BP86 may be partly attributable to the use of a full-size model of MeCbl. Another reason for the small BDE obtained with BP86 is the large contribution of the ZPE effect, which reduced BDE by ∼5 kcal/mol. This effect has not always been taken into account in previous studies. Jensen and Ryde reported that other pure GGA functionals, such as BLYP and BPW91 give smaller BDE values than BP86. From this observation and our current results, it seems likely that pure GGA functionals tend to underestimate the BDE of MeCbl, which may partly be attributable to their inability to describe the dispersion effect. However, the addition of the dispersion term to BP86 led to an overestimation of BDE (44.044.1 kcal/mol at the BP86-D/B2//B1 level). 9311
dx.doi.org/10.1021/jp2052807 |J. Phys. Chem. A 2011, 115, 9308–9313
The Journal of Physical Chemistry A
ARTICLE
Table 3. Summary of the Components of BDE (kcal/mol) ΔEB1
Δdisp
ΔZPE,B1
ΔH,B1
ΔBS,B2
ΔBSSE,B1
ΔBSSE,B2
Δsolv1,B1
Δsolv2,B1
Δsolv1,B2
Δsolv2,B2
Δrel,B1
Δrel,B2
B3LYP
28.4
0.0
5.5
1.5
7.3
6.0
0.7
1.2
1.2
1.1
1.2
2.7
1.5
B3LYP-D
29.8
13.3
5.3
1.4
7.0
5.7
0.7
1.2
1.2
1.1
1.2
3.0
1.6
BP86
40.7
0.0
5.1
1.4
6.0
6.1
0.8
1.0
1.0
0.9
0.9
2.6
1.1
BP86-D
41.5
13.2
5.6
1.5
5.9
5.8
0.8
1.0
1.1
1.0
1.0
2.9
1.1
B3LYP*-D//BP86
33.2
15.5
5.1
1.4
6.8
6.1
0.8
1.0
1.1
1.0
1.0
2.6
1.3
B3P86-D//BP86
33.2
15.5
5.1
1.4
5.9
5.4
0.8
1.1
1.2
1.1
1.1
2.7
1.4
Although the dispersion correction of B3LYP (i.e., B3LYP-D) improved the BDE, the values (31.831.9 kcal/mol with B2) were still smaller than the experimental values. On the assumption that the empirical dispersion energy is accurate, this result implies that there is room for improvement in the EDFT part of eq 2. Therefore, we examined another two methods, B3LYP*-D and B3P86-D, using the BP86-optimized geometry. Siegbahn et al. suggested that the amount of HF exchange in B3LYP should be reduced from 20% to 15% (B3LYP*). Applying the B3LYP*-D method to a simplified MeCbl model, they obtained a BDE value of 32 kcal/mol. Interestingly, the application of this method to a full model of MeCbl yielded 40.5 and 37.7 kcal/mol with B1 and B2, respectively (Table 2). These values agree well with the experimental values, with B2 showing better agreement. Jensen and Ryde found that the P86 correlation functional tended to give larger BDE values than the LYP correlation functional. This prompted us to think that using B3P86-D, without altering the HF exchange, may also improve the BDE value. Indeed, the BDE value calculated with this method also showed good agreement with the experimental values (41.241.3 kcal/mol with B1 and 38.638.7 kcal/mol with B2). Closer inspection of the correction terms of BDE (Table 3) revealed that Δdisp (13.215.5 kcal/mol for DFT-D), ΔZPE,B1 (5.6 to 5.1), ΔBS,B2 (7.3 to 5.9), and ΔBSSE,B1 (6.1 to 5.4) were large in magnitude, whereas the other contributions were relatively small. The enthalpy (ΔH,B1) and relativistic (ΔREL,B1) corrections had positive values, whereas ΔBSSE,B2, ΔsolvN,B1, and ΔsolvN,B2 had negative values. The large ΔBS,B2 seems to be largely because of the large BSSE (∼6 kcal/mol) with B1, compared with the BSSE of
