The Problems: (see appendix for tabulated data)
1. Last problem set, we used AM1 and PM3 to calculate the heats of formation of the following molecules: F2, Cl2, ClF, and ClF3. Results were awful. We will now continue at the ab initio level. Since the first three are small molecules, we will attempt to reach basis set convergence with respect to various properties. To do that, we will use the correlation-consistent polarized valence double zeta (cc-pVDZ is the abbreviation used in the Gaussian94 keyword line), triple zeta (cc-pVTZ), and quadruple zeta (cc- pVQZ) basis sets of Dunning. The QZ set, besides being valence quadruple zeta, includes 2 sets of f functions and 3 sets of d functions on all atoms, i.e., it's big. The QZ calculations are so expensive that I have done them for you, and they are in the appendix--I have also filled in the CCSD(T)/cc- pVTZ energies there, so you only need to do up to a CCSD calculation to fill in the rest of the cc-pVTZ rows. Fill in all empty positions in the appendix tables. Note that all of the correlated calculations are done at the HF/cc-pVQZ geometries! You must use the right bond lengths from my calculations to get the right energies--bond lengths are: F2, 1.32749788 Å; Cl2, 1.97826267 Å; ClF, 1.58547516 Å.
To fill in the thermal correction block (because we want to compare to an enthalpy of reaction), perform a frequency calculation for each molecule at the HF/cc-pVDZ level. In the output file, take the value listed as Thermal correction to Enthalpy=.
The "best estimate" block at the bottom involves taking the best calculation that can be afforded with the cc-pVQZ basis, and correcting it for higher orders of correlation using results obtained with a smaller basis set, and finally adding the enthalpic correction. That is, Best estimate = CCSD/cc-pVQZ + [CCSD(T)/cc-pVTZ - CCSD/cc-pVTZ] + DDH. (The term in brackets thus analyzes the effects of triple excitations on the correlation energy using a basis set where the calculation can be afforded.)
Questions: How well converged are the RHF bond lengths by the cc-pVQZ level? Based on your results, estimate how much remaining error there might be relative to a "saturated" basis set. What about the convergence of the absolute energies? Are the correlated energies better converged or less well converged with respect to basis set size than the RHF energies? How does the best estimate for the heat of reaction to form ClF compare with the experimental value of -26.6 kcal/mol?
The bond lengths are very well converged. On going from DZ to TZ there are variations on the order of hundredths of an Ångstrom, on going from TZ to QZ this drops to thousandths of an Ångstrom. Presumably further doubling of the basis set size would have effects of only a few ten-thousandths, which is chemically irrelevant. The absolute energies are converging less quickly. Moreover, when one looks at the energy of reaction section, one finds that the convergence is not uniform for all three molecules. From TZ to QZ, there is a change of 2.3 kcal/mol in the reaction energy at the RHF level. The correlated energies are even less well converged (this is a well known feature of the correlation energy--it converges more slowly with respect to basis set than the RHF energy--the reason is that the correlated wave function is more diffuse than the RHF, so it is more important to have a very flexible basis set). Nevertheless, the best estimate of -27.3 is in quite good agreement with experiment (anything below about 1 kcal/mol is typically referred to as "chemical accuracy")
Other than using a still larger basis set, how could this calculation be further improved?
Probably the largest remaining source of error is the use of RHF geometries. Using correlated geometries would be better (e.g., MP2 or CCSD--other levels do not lend themselves very well to geometry optimization) This would be true for the energies and also possibly the thermal corrections.
Recall that neither AM1 nor PM3 could reproduce the correct geometry for ClF3. Optimize the geometry of ClF3 at the RHF/cc-pVDZ level and report it with all necessary degrees of freedom. Using all of your data so far, make a best guess at the heat of formation of this molecule, and explain your reasoning. Don't do any additional calculations to answer this--you can't afford it!
The molecule is indeed T shaped with an equatorial ClF bond of 1.591 Å and two axial ClF bonds of length 1.691 Å. The angle between the equatorial and axial bonds is 86.4deg.. The molecule has C2v symmetry. A best guess for heat of formation can be derived in the following manner. We can calculate the energy of reaction for F2 + ClF -> ClF3 at the RHF/cc-pVDZ, and that is +15.0 kcal/mol. That same level of theory for F2 + Cl2 -> 2ClF underestimated experiment by 2.4 kcal/mol (not exothermic enough), so we might correct our ClF3 estimate by that much, arriving at +12.6. Since the heats of formation are known for F2 and ClF (0.0 and -13.3 kcal/mol, respectively), DHf for ClF3 must be in the vicinity of -0.7 kcal/mol (note this is exactly in between the semiempirical results of Problem Set 2). A conservative error bar would probably be +/- 5 kcal/mol. Other reasonable ways to estimate the heat of formation may be imagined, and will be given appropriate credit in grading.
2. Recall the ammonia N-oxide example used in class (if you took notes, you've already got a starting Z-matrix, although you'll recall it wasn't a great initial guess). This molecule is a tautomer of hydroxylamine, NH2OH. The two tautomers are interconverted via a 1,2- proton shift. Calculate the relative stabilities of the two tautomers and the barrier to tautomerization at the levels indicated in the appendix table. Also, provide the bond distances rNO, rNH, and rOH for the transition state structure at all levels of theory (where H is the proton in flight). It will be more efficient if you first find all of the structures at the RHF level, and then use these to start calculations at other levels. Verify your RHF transition state structure with a frequency calculation.
Questions: What is the magnitude of the imaginary frequency in your RHF TS structure? Looking at the associated displacement vector, what atom is moving the most in the imaginary vibration (which, of course, corresponds to the reaction coordinate)? Does the correlation energy appear converged by the MP4 level with respect to tautomer energies and barrier height? How do the RHF and the DFT results compare to the correlated levels? If you were going to assign a best estimate with an error bar for the barrier, what would it be (i.e., pick the numbers in the sentence: "The barrier is ##.# +/- ##.# kcal/mol"). Explain your reasoning.
RHF imaginary frequency is 2016i (don't be confused by computer output, which flags this as an imaginary frequency by listing it as a negative number. It is not negative--the force constant is negative. Since the harmonic frequency is related to the square root of the force constant, it is imaginary). The vector indicates that the proton in flight is moving much, much more than the heavy atoms. This is to be expected based on their mass compared to the proton (not to mention the nature of the reaction coordinate). The relative tautomer energies vary over only 0.6 kcal/mol from MP2 to MP4, suggesting that reasonable convergence has been obtained. However, the reaction barrier varies by considerably more, 3.0 kcal/mol, indicating that it is not yet converged. This is typical for MPn methods--transition states typically have smaller HOMO-LUMO gaps, so they converge more slowly with respect to correlation energy. By comparison to MP2, RHF predicts a "tighter" transition state structure (rOH reduced by a full 0.11 Å) and a much higher barrier. It is typical of the RHF level to overestimate barrier heights in reactions where a pair of electrons is being transferred between different bonds. DFT (in this case the BLYP functional) has a looser structure than MP2, particularly in the NO bond length. It also predicts both the transition state structure and ammonia N-oxide to be significantly lower in energy than MPn. Based on these results, I would state that the barrier is about 52 +/- 5 kcal/mol. I picked a number between MP4 and DFT but closer to DFT because DFT energies converge more rapidly with respect to basis set size (and we're using a pretty small one) than MPn energies. 5 kcal/mol is the separation between BLYP and MP4, so it seems a reasonable guess at an error bar. However, it could be that errors associated with the 6-31G* basis are in excess of 5 kcal/mol. Considerable literature suggests that this will not be the case however--there always seem to be fortuitous cancellations of errors with the 6-31G* basis that cause it to give better results than one has a right to expect.