|Home | About | Journals | Submit | Contact Us | Français|
Heats of formation were calculated using coupled-cluster methods for a series of zinc complexes. The calculated values were evaluated against previously conducted computational studies using density functional methods as well as experimental values. Heats of formation for nine neutral ZnXn complexes [X = -Zn, -H, -O, -F2, -S, -Cl, -Cl2, -CH3, (-CH3)2] were determined at the CCSD and CCSD(T) levels using the 6–31G** and TZVP basis sets, as well as the LANL2DZ-6–31G** (LACVP**) and LANL2DZ-TZVP hybrid basis sets. The CCSD(T)/6–31G** level of theory was found to predict the heat of formation for the non-alkyl Zn complexes most accurately. The alkyl Zn species were problematic in that none of the methods that were tested accurately predicted the heat of formation for these complexes. For the seven non-alkyl species, the CCSD(T)/6–31G** level of theory was shown to predict the most accurate heat of formation values. In instances where experimental geometric parameters were available, these were most accurately predicted by the CCSD/6–31G** level of theory; going to CCSD(T) did not improve agreement with the experimental values.
Zinc complexes are critically important in biological systems, serving in both a structural and a catalytic capacity.1,2 Indeed, zinc trails only iron as the most ubiquitous transition metal in biological systems. Key complexes include, but are certainly not limited to, Human Carbonic Anhydrase,3–5 Carboxypeptidase,6 Alcohol Dehydrogenase7 and so-called “zinc fingers”8,9 which play structural roles in DNA recognition. Many of these systems have been studied extensively via X-ray crystallography and spectroscopic methods, including NMR.
The literature to date contains many computational studies on systems that contain zinc. In a 1991 study, Kaupp et al. probed the structures of ZnR2 complexes with 1,4-diaza-1,3-butadienes using pseudopotential calculations.10 A subsequent study by Kaupp and von Schnering probed the structures and binding energies of (ZnX2)2 dimers at the MP2 and HF levels of theory using pseudopotentials.11 Kabelac and Hobza examined the binding of Zn2+ with the nucleic acid bases adenine, guanine, cytosine and thymine at the MP2/TZVP level of theory.12 In recent work, Rayon and coworkers have probed binding energies and geometries of several ZnII complexes using MP2 and density functional methods for optimizations, along with single point calculations at the CCSD(T)/aug-cc-pVTZ level of theory.13
We have recently employed a variety of popular density functional methods spanning the GGA, meta-GGA, hybrid-GGA and meta-hybrid-GGA functional classes, a total of 12, in the calculation of both heat of formation (ΔHf) and ionization potential for transition metal complexes.14 This study utilized both the Pople style 6–31G** and the TZVP basis sets. Further efforts within the group have applied the same set of density functionals in conjunction with the pseudopotential based LANL2DZ15 basis set on Zn.
Although not as abundant as work with lower level methods, some literature precedent for the application of coupled-cluster methods to ΔHf calculations in small organic systems does exist. Dixon and coworkers have investigated iodine fluorides using CCSD(T) methods16 and applied similar methodology in the study of diazomethane.17 Feller and Franz used large basis set coupled cluster calculations to determine ΔHf values for furan and other related derivatives.18 Related work by Feller and coworkers includes the determination of ΔHf values for combustion products19, oxyfluorides20, and borohydrides21. Harding and coworkers evaluated the ΔHf value of vinyl chloride using high level ab initio calculations under the HEAT protocol with geometries optimized at the CCSD(T) level.22 Marenich and Boggs performed CCSD(T) calculations to investigate formyl and isoformyl radical species.23 Additionally, Lee has studied compounds of the form XONO2 (X = H, F, Cl) using coupled-cluster methods and determined a heat of formation value for FONO2 using isodesmic reactions.24 Martin and coworkers have probed thermodynamic properties of bromoalkanes25 and boron trifluoride26, including investigation of relativistic effects.
There is also some precedent for the use of these methods in ΔHf calculations on metal-containing systems. A popular approach is extrapolation towards the complete basis set (CBS) limit and has been applied to transition metal compounds from Sc-Zn. 27,28 Further work by Peterson and coworkers describes the application of the correlation consistent Composite Approach to the thermochemistry of transition metal systems.29 Sullivan et al. used CCSD(T) calculations to determine ΔHf values for oxide and hydroxide complexes of alkali and alkaline earth metals.30 Similarly, Nielsen and coworkers used CCSD(T) calculations in the study of tin—oxygen compounds.31 Finally, Lu and coworkers have recently published high level calculations on transition metal-ammonia complexes at the CCSD(T) level extrapolated to the CBS limit in order to accurately predict ionization potentials.32 While it may not be practical to apply the highly accurate CCSD(T)-CBS extrapolated model to larger transition metal containing systems, we desired to probe the efficacy of CCSD and CCSD(T) calculations with smaller basis sets in predicting the heats of formation of such systems, since to our knowledge this has not been examined closely. Moreover, as far as we are aware extensive coupled-cluster calculations including full geometry optimizations have not been used in the prediction of ΔHf values for Zn complexes.
We desired to expand our efforts towards the calculation of heat of formation values to include more powerful computational methods. Our current efforts are focused on the application of coupled-cluster methods toward this end. Higher level calculations such as CASPT2 were not employed, as these are not viable for larger Zn-containing systems and one of our primary reasons for these investigations is the uncovering of computational methodologies suitable for use in QM and QM/MM type calculations on biological systems incorporating one or more Zn centers. With this in mind we have focused on using relatively modest basis sets (e.g., 6–31G** and TZVP) combined with the CCSD (coupled-cluster with single and double excitations) and CCSD(T) (coupled-cluster with single and double and perturbative triple excitations) methods for the purpose of this study. CCSD(T) provides an excellent compromise between accuracy and computational cost and we wish to evaluate its performance in determining ΔHf values for ZnXn complexes. CCSD is less computationally intensive, and we evaluate this method as well because it may be more amenable to the study of larger Zn systems. In the future the effect of larger basis sets will be explored, but the present effort represents a balance between accuracy and computational efficacy.
We have chosen a series of zinc complexes to be the focus of our initial work with coupled cluster methods for ΔHf calculations. The nine chosen ZnXn compounds are Zn2, ZnH, ZnO, ZnF2, ZnS, ZnCl, ZnCl2, ZnCH3 and Zn(CH3)2 and these were selected based on availability of literature heats of formations. These compounds present a mix of both open and closed shell species, albeit with either singlet or doublet ground state multiplicities except in the cases of ZnS and ZnO (vide infra). Table 1 contains experimental heats of formation with errors and the calculated values using the “best” density functional method with both the 6–31G** and TZVP basis sets from the work of Riley and Merz.14 These “best” values represent the density functional providing the smallest deviation from the reported experimental values with each respective basis set. The errors associated with the literature values for this set are relatively small, with the noted exceptions of ZnCl and Zn(CH3)2 whose errors are 15.4 and 15.5% respectively relative to the reported values. One can readily identify the problem with the application of DFT methods in calculations on Zn complexes; there is no “universal method” providing optimal results. It is hoped that coupled-cluster methods can systematically reduce the overall error and either provide a more global approach to calculations on Zn complexes, or at least demonstrate that the best approach to heats of formation in these complexes is the identification of a complex-specific density functional.
All calculations were carried out on a SUN cluster featuring dual 2.5GHz Opteron nodes using the Gaussian 0334 suite of programs. All geometry optimizations incorporated standard gradient methods. For all single point calculations, the SCF=TIGHT keyword was used. The SCF=XQC keyword was applied in all instances, as SCF convergence was often problematic, especially for higher energy multiplets. CCSD and CCSD(T) calculations were run as implemented in Gaussian 03.35–40 Frequency calculations were conducted on all geometries (at the minimum energy multiplicity) to insure all calculated lowest energy structures resided at local minima on the potential energy surface. Where applicable, calculations were done at the UCCSD at UCCSD(T) levels. All other calculations are closed shell. The Pople type 6–31G** and triple-ζ quality TZVP basis sets were used as implemented in Gaussian 03.41,42 LACVP** calculations were run using the GEN keyword for the basis set. In these calculations, the LANL2DZ basis/pseudopotential was used for Zn and the 6–31G** basis set for the nonmetal atoms. A second set of calculations was run which applied the TZVP basis set to the nonmetals while retaining LANL2DZ on the Zn atom. T1 diagnostics were computed with Gaussian 03 at the CCSD/6–31G** and CCSD/TZVP levels of theory.43,44 This is a measure that identifies instances where multi-reference effects may be important. While multi-reference methods are beyond the intent of this investigation, but we have placed these values in the supporting information, as the results identify several compounds for which multi-reference approaches should be considered (ZnO, ZnS and ZnF2).
For all Zn species considered, we initially desired to optimize the 1, 3, 5 and 7 multiplicities for even electron species and the 2, 4, 6 and 8 multiplicities for odd electron species as done in our previous DFT work.14 This worked well for most CCSD calculations, although it was sometimes difficult to achieve SCF convergence for high energy multiplicities. CCSD(T) calculations failed for a large number of high energy multiplicities, although the CCSD ground state could always be converged for smaller complexes using CCSD(T) calculations. Zn(CH3)2 proved to have serious convergence problems in the geometry optimization procedure for the CCSD(T) calculations, which uses the numerical eigenvector following algorithm in Gaussian, and attempts at the CCSD(T) level were abandoned for this system at all multiplicities. At C1 symmetry, the lowest for this system, Zn(CH3)2 experienced difficulties with the number of variables as well. Enforcing D3h symmetry to decrease the variables considered did not provide any relief for convergence related problems.
Heats of formation (ΔHf) for all complexes were computed using the method outlined in the Gaussian white paper on Thermochemistry in the Gaussian 03 online manual.45 These calculations follow equation (1) that simplifies to equation (2), derived from the procedures outlined in the Gaussian white paper. The ‘M’ and ‘X’ designations in equation one correspond to the molecule and individual atoms, respectively. For convenience, equation (2) is provided in terms of the output provided by Gaussian 03.
ECORR is identified as the sum of electronic and thermal enthalpies from the output of the Gaussian frequency calculation (which includes thermal and ZPE corrections to the energy). EZn and Eatom are the energies of the Zn and nonmetal atoms at a given level of theory. The constant 31.17 (kcal/mol) in Equation 2 is the ΔHf (Zn, 298K) taken from the NIST chemistry WebBook46 and the respective ΔHf (atom, 298K) values for the nonmetals are found there as well. Finally, equation 3 was implemented in the calculation of root mean squared deviations (RMSD).
Summarized in Table 2 are calculated heats of formation for Zn complexes at the CCSD and CCSD(T) levels for all applications of the 6–31G** basis set (stand-alone, and as part of LACVP**). A plot of all calculated values versus the experimental ΔHf for an illustrative comparison is given in Figure 1. For each metal entry, most data points are grouped rather closely together. Significant deviations from the experimental values were found in calculated ΔHf values for several of the complexes investigated. Most ZnXn species at all levels of theory possessed a calculated ΔHf substantially higher than the experimental value. The exceptions were 1ZnF2 at both the CCSD/6–31G** and CCSD(T)/6–31G** levels of theory. These were slightly lower than the experimental value, at 1.0 and 2.9 kcal/mol respectively. All other complexes had overestimated ΔHf values at each method/basis set combination. ΔHf values for diatomic 1Zn2 are consistently calculated for both CCSD and CCSD(T) calculations with the 6–31G** and LACVP** basis sets, however the errors across the board are nearly 7 kcal/mol. For 2ZnH the calculated values are off by nearly the same amount for both coupled cluster methods with the LACVP** basis set/pseudopotential but with the 6–31G** basis set the calculated heats of formation are within 1.0 kcal/mol of the experimental value, which represents excellent agreement. The ΔHf values calculated for 1ZnS were significantly larger than the experimental values. The percentage errors in the calculated values ranged from 45.6% at CCSD(T)/6–31G** to 80.7% at CCSD/LACVP**. Very large errors were observed in the calculated ΔHf values for 1ZnCl2 at both the CCSD/LACVP** and CCSD(T)/LACVP** levels of theory (49.4 and 46.9% errors respectively).
The calculated ΔHf values for 2ZnCH3 were poor across all investigated levels of theory. The best calculated value was 61.4 kcal/mol at the CCSD(T) level of theory, which compared to the experimental value is 136.2% higher in energy. Both values calculated using the LACVP** basis set were significantly worse, with percentage errors approaching 200%. So while it may be argued that the CCSD(T)/6–31G** level of theory is the best level of theory presented for 2ZnCH3, it is quite clear that none of the these combinations were really appropriate for this species. In the cases where 1Zn(CH3)2 can be evaluated, the values are again poor. A full discussion of this complex is presented elsewhere (vide infra).
Overall, calculations using the 6–31G** basis set always perform better than their LACVP** counterparts. In 6 of 8 cases, CCSD(T)/6–31G** outperforms CCSD/6–31G** with the two exceptions being 2ZnH (by 0.1 kcal/mol) and 1ZnF2 (by 1.6 kcal/mol). CCSD(T)/LACVP** outperforms CCSD/LACVP** in all cases except for 2ZnH. The average error associated with CCSD(T)/6–31G** is -10.1 kcal/mol, a 2.5 kcal/mol improvement over the average for CCSD/6–31G**. The average error increases nearly twofold for both CCSD and CCSD(T) calculations using the LACVP** basis set.
Next we will focus on the geometries of species for which experimental data are available: 2ZnH, 1ZnF2, 1ZnCl2 and 1Zn(CH3)2. Table 3 contains a summary of Zn-X bond lengths for all complexes, with the aforementioned available literature values. For 2ZnH, the Zn-H bond length is calculated to within 0.001 A at the CCSD/6–31G** and CCSD(T)/6–31G** levels of theory. The deviation from experimental is significantly larger for both coupled-cluster methods using the LACVP** basis set; 0.057 A for CCSD and 0.058 A for CCSD(T). For 1ZnF2, the reverse trend is observed in that coupled-cluster methods incorporating the LACVP** basis set calculate Zn-F bond lengths to within 0.01 A, whereas methods utilizing strictly the Pople-style 6–31G** basis set arrive at equilibrium bond lengths 0.033 A lower than the experimental value. For 1ZnCl2, CCSD/6–31G** predicts the most accurate bond length, within 0.007 A of the experimental value (CCSD(T)/6–31G** is nearly as accurate, off by 0.008 A). Both methods with the LACVP** basis set are off by at least 0.04 A. Finally, the Zn-C bond distance in 1Zn(CH3)2 is calculated to within 0.002 A at the CCSD/6–31G** level of theory and within 0.063 A at CCSD/LACVP**. The Zn-X distances for which experimental data are not available are scattered at best amongst the levels of theory, but with no literature values available it is difficult to assess which methods are “correct” in their predictions. Three points are immediately obvious: (1) the 6–31G** basis set is preferred over LACVP** for these heat of formation predictions, (2) no obvious trend with respect to over- or underestimation of bond lengths and (3) applying the CCSD(T) level offers no significant improvement on calculated equilibrium bond lengths with a constant basis set; indeed, the geometries typically deviate more from the experimental value using the higher-level method.
In order to further probe the effect of basis set on these calculations, we ran all calculations using the TZVP basis set. Calculated heats of formation for ZnXn complexes at the CCSD and CCSD(T) levels for all applications of the TZVP basis set (stand-alone, and in conjunction with LANL2DZ) are summarized in Table 4 and Figure 3. All calculated ΔHf values were overestimated using these two basis sets. As with the 6–31G** basis sets, the two methyl Zn species had poorly predicted ΔHf values. The 1Zn2 and 2ZnH complexes were again consistently well calculated by all methods used. As with the 6–31G** and LACVP** basis sets, the calculated ΔHf values for ZnS were very poor, with nearly 100% errors at all four levels of theory. Large errors for ZnO were also observed in each instance. Errors in ZnS and ZnO are further discussed in the section dealing with the calculation of their ground state spin multiplicities. The ΔHf values for 1 ZnCl2 were best calculated at the CCSD/TZVP and CCSD(T)/TZVP levels of theory with percentage errors of 26.6 and 23.5% respectively, so clearly overall performance for calculations involving this dichloride were rather poor.
As was observed using the Pople-style basis set, TZVP and LANL2DZ-TZVP results for 2ZnCH3 were quite poor. All errors well surpassed 100% of the experimental value, ranging from a low of 162.7% for CCSD(T)/LANL2DZ-TZVP to a high of 175.4% at the CCSD/TZVP level of theory. Of eight methods utilized during the course of this study, not one does an adequate job of predicting the ΔHf value for 2ZnCH3.
The predicted bond lengths were compared against experimental values for 2ZnH, 1 ZnF2, 1ZnCl2 and 1Zn(CH3)2 (Table 5). The CCSD(T)/TZVP geometry for 2ZnH was closest to the experimental rZn–H, overestimating by 0.025 Å. The CCSD/TZVP level of theory best described the Zn-F bond length in 1ZnF2, off by only 0.006 A. Both the CCSD and CCSD(T) methods in conjunction with the TZVP basis set gave a Zn-Cl bond distance of 2.109, 0.037 A higher than the experimental rZnCl value in 1ZnCl2. 1Zn(CH3)2 could only be optimized at the CCSD/LANL2DZ-TZVP level of theory and rZn-C was overestimated by 0.056 A. For the three complexes that could be optimized at all levels of theory incorporating TZVP or LANL2DZ-TZVP, the TZVP basis set performed better although there was very little variation in the results for 1ZnF2. There was no clear separation between the CCSD and CCSD(T) methods, with CCSD better for 1ZnF2, CCSD(T) better for 2ZnH (both by slim margins) and both methods producing identical geometries for 1ZnCl2.
Overall, CCSD(T)/6–31G** and CCSD(T)/TZVP calculations best predict the heat of formation for the 7 non-alkyl Zn complexes. These methods are compared in Figure 3. The CCSD(T)/6–31G** level of theory is shown to slightly outperform the CCSD(T)/TZVP level in these ΔHf predictions. The difference in average errors between these two methods is 5.0 kcal/mol, with CCSD(T)/631G** averaging a −9.0 kcal/mol deviation and CCSD(T)/TZVP differing from the experimental by an average of −14.0 kcal/mol for the 7 non-alkyl complexes.
Generally, all levels of theory correctly predicted the appropriate spin ground states for the Zn species investigated. Open shell species were found to be ground state doublets, and closed shell species ground state singlets. There are two notable exceptions, and these are the cases of ZnO and ZnS. The expected ground state multiplicity, in so far as what species the heat of formation is determined for, is ambiguous. Experimental work indicates that the available heat of formation and Zn-X bond lengths are for the triplet.47,48 However, the most accurate calculations to date on these complexes predict a ground state singlet multiplicity.49,50 The triplet ground state was predicted for ZnO at the CCSD/LACVP**, CCSD(T)/TZVP and CCSD/LANL2DZ-TZVP levels of theory. The triplet was found to be the ground state of ZnS at only the CCSD/LACVP** and CCSD/LANL2DZ-TZVP levels of theory. The pseudopotential incorporating basis sets with the CCSD level of coupled-cluster theory most often arrive at the triplet ground state for both ZnO and ZnS. The lowest errors in calculated ΔHf values were observed using the TZVP and 6–31G** basis sets, which most frequently predict the singlet ground state multiplicities. For ZnO, theoretical studies support a triplet ground state for bond lengths in excess of 1.85A and a singlet for distances closer to the reported internuclear distance of 1.69 Å.50 These findings are in complete agreement with the predicted ground states in our study. Other investigations in our laboratory have turned up similar confusion in predicted ground state spins—in our LANL2DZ investigation,15 a singlet ground state was predicted for all GGA and meta-GGA functionals and triplet ground state for all hybrid-GGA and hybrid-meta-GGA functionals.
We decided to further our investigation by applying a Douglas-Kroll-Hess 2nd order relativistic correction (DKH) to calculations at the CCSD/6–31G** level of theory as implemented in Gaussian 03.51–55 This correction was applied during the course of both the geometry optimizations and frequency analyses. A comparison of CCSD/6–31G** ΔHf values with and without this correction is provided in table 6. Addition of the DKH correction actually results in a slight increase in the average heat of formation error. The error is decreased with this correction for 1Zn2, ZnO, and 1ZnF2 while it is increased for the remaining entries. It is worthwhile to point out that not only does the addition of a 2nd order relativistic correction improve the calculated ΔHf value for ZnO, it also results in the prediction of a triplet ground state for the species, as opposed to the singlet predicted by the standard CCSD/6–31G** calculation.
Table 7 contains results for CCSD(T)/6–31G** calculations including relativistic effects for six ZnXn complexes. The alkyl Zn species were omitted, and efforts to include this correction for 1Zn2 failed. The average error increases by 1.4 kcal/mol with the inclusion of this correction compared to the same set of six compounds using standard CCSD(T)/6–31G** calculations. Only the error for 1ZnF2 decrease with the 2nd order correction, while all other errors increase in magnitude. The most significant increase in error is associated with ZnO. There is also a ground state multiplicity change for this species, which is predicted to be a triplet with inclusion of the DKH correction whereas at the CCSD(T)/6–31G** level of theory it is predict to possess a ground state singlet multiplicity.
There are two viable conformations for 1Zn(CH3)2, specifically conformations which have pseudo-eclipsed hydrogen atoms, and pseudo-anti hydrogen atoms (Figure 4). Both conformations were examined during the course of this study. The pseudo-eclipsed conformation was found to be a minimum at both the CCSD/6–31G** and CCSD/TZVP levels of theory, while the pseudo-anti conformer was calculated to be a transition state as evidenced by the presence of one negative mode in the vibrational analysis. Still, the energy difference between the two conformations is quite small 0.03 kcal/mol at CCSD/6–31G** which implies virtually no barrier to free rotation of the methyl groups as this value is lower than kT. This is not surprising, as a bridging Zn atom separates these substituents. A plot of the unrelaxed potential energy surface for methyl rotation is displayed in Figure 5 (relative energy as a function of H-C-C-H twist angle).
The eclipsed conformation is perhaps favored over the gauche due to stabilizing hyperconjugative interactions with the d-orbitals of the bridging Zn atom and the C-H antibonding orbitals; the eclipsed conformation affords more of these interactions than does the gauche. Therefore the potential energy surface is essentially a 60° shift of that seen in ethane (where the pseudo-gauche conformation is most stable and the pseudo-eclipsed is a rotational transition state). We are currently investigating bulkier alkyl substituents (e.g. ZnEt2 and ZniPr2) to see if this trend of eclipsed conformational preference continues. Unfortunately, the size of these species will preclude any investigations using higher level methods such as CCSD(T), but we are confident that density functional methods and CCSD calculations will prove adequate in these efforts.
For nine ZnXn complexes, ΔHf values were calculated using the CCSD and CCSD(T) coupled cluster methods in conjunction with the 6–31G**, TZVP, LANL2DZ-6–31G** and LANL2DZ-TZVP basis sets. Generally, the 6–31G** basis set was found to outperform the other three and the CCSD(T)/6–31G** level of theory provided ΔHf values that compared most favorably to the experimental values, slightly outperforming CCSD(T)/TZVP calculations. For the four complexes for which experimental equilibrium bond lengths were available, the CCSD/6–31G** level of theory generally provided the best geometries, with CCSD(T)/6–31G** bond lengths actually deviating more from the experimental distances in most cases. The lone exception to this trend was 1 ZnF2, for which the CCSD/TZVP level of theory predicted the most accurate rZn-F value. CCSD(T) calculations could not be performed on 1Zn(CH3)2 due to the size of the system and CCSD calculations on this system provided poor heat of formation values. All levels of theory resulted in poor predicted heats of formation for 2ZnCH3.
In general, CCSD/6–31G** calculations seem appropriate for the prediction of heat of formation values for non alkylated species, with the CCSD(T) method providing slight improvement. The CCSD/6–31G** does a good job reproducing the experimental bond distances in these systems. The addition of a 2nd order Douglas-Kroll-Hess relativistic correction does not provide overall improvement in the calculated ΔHf values. Rather, this term offers modest improvement in some instances while actually resulting in an increased error for other entries. The fact that CCSD values are comparable is important, as this method is clearly more amenable to calculations on larger systems whereas the CCSD(T) calculations would quickly become too resource-intensive to be viable. We are currently extending our research to include complexes of other third row transition metals to determine the extent of applicability of coupled-cluster methods to these systems.
We thank the NIH (GM066859 and GM044974) for supporting this research. MNW wishes to thank the NIH for support in the form of an NRSA postdoctoral fellowship (F32GM079968). We thank the members of the Merz group for useful discussions during the course of this research. MNW specifically thanks Duane Williams for his numerous editorial comments.