«Jan M.L. Martin∗ Department of Organic Chemistry, Kimmelman Building, Room 262, Weizmann Institute of Science, 76100 Rehovot, Israel. Email: ...»
A deﬁnitive heat of vaporization of silicon through benchmark ab
initio calculations on SiF4
Jan M.L. Martin∗
Department of Organic Chemistry, Kimmelman Building, Room 262, Weizmann Institute of
Science, 76100 Rehovot, Israel. Email: firstname.lastname@example.org
Peter R. Taylor
San Diego Supercomputer Center and Department of Chemistry and Biochemistry MC0505,
University of California, San Diego, San Diego, CA 92093-0505, USA. Email: email@example.com
(Submitted to J. Phys. Chem. A February 1, 1999) Abstract In order to resolve a signiﬁcant uncertainty in the heat of vaporization of silicon — a fundamental parameter in gas-phase thermochemistry — ◦ ∆Hf,0 [Si(g)] has been determined from a thermochemical cycle involving the precisely known experimental heats of formation of SiF4 (g) and F(g) and a benchmark calculation of the total atomization energy (TAE0 ) of SiF4 using coupled-cluster methods. Basis sets up to [8s7p6d4f 2g1h] on Si and [7s6p5d4f 3g2h] on F have been employed, and extrapolations for residual basis set incompleteness applied. The contributions of inner-shell correlation (−0.08 kcal/mol), scalar relativistic eﬀects (−1.88 kcal/mol), atomic spinorbit splitting (−1.97 kcal/mol), and anharmonicity in the zero-point energy (+0.04 kcal/mol) have all been explicitly accounted for. Our benchmark ◦ TAE0 =565.89±0.22 kcal/mol leads to ∆Hf,0 [Si(g)]=107.15±0.38 kcal/mol ◦ (∆Hf,298 [Si(g)]=108.19±0.38 kcal/mol): between the JANAF/CODATA value of 106.5±1.9 kcal/mol and the revised value proposed by Grev and Schaefer [J. Chem. Phys. 1992, 97, 8389], 108.1±0.5 kcal/mol. The revision will be relevant for future computational studies on heats of formation of silicon compounds. Among standard computational thermochemistry methods, G2 and G3 theory exhibit large errors, while CBS-Q performs relatively well and the very recent W1 theory reproduces the present calibration result to
I. INTRODUCTION For three of the ﬁrst-and second-row elements, namely Be, B, and Si, the tabulated heats of formation of the atoms in the gas phase carry experimental uncertainties in excess of 1 kcal/mol. Aside from being propagated into uncertainties for experimental gas-phase thermochemical data for compounds involving these elements, they adversely aﬀect the accuracy of any directly computed heat of formation — be it ab initio or semiempirical — of any Be, B, or Si-containing compounds through the identity
Particularly given the importance of boron and silicon compounds, this is a rather unsatisfactory state of aﬀairs.
Recently we succeeded  in reducing the uncertainty for boron by almost an order of magnitude (from 3 kcal/mol to 0.4 kcal/mol) by means of a benchmark calculation of the total atomization energy (TAE0 ) of BF3 (g). By combining the latter with the experimentally precisely known  heat of formation of BF3, we were able to indirectly obtain the vaporization enthalphy of boron to high accuracy. It was thus shown that a 1977 experiment by Storms and Mueller , which was considered an outlier by the leading compilation of thermochemical tables , was in fact the correct value.
The heat of formation of Si(g) is given in the JANAF  as well as the CODATA  tables as 106.5±1.9 kcal/mol. Desai  reviewed the available data and recommended the JANAF/CODATA value, but with a reduced uncertainty of ±1.0 kcal/mol. Recently, Grev and Schaefer (GS)  found that their ab initio calculated TAE[SiH4 ], despite basis set incompleteness, was actually larger than the value derived from the experimental heats of formation of Si(g), H(g), and SiH4 (g). They suggested that the heat of vaporization of ◦ silicon be revised upwards to ∆Hf,0 [Si(g)]=108.07±0.50 kcal/mol, a suggestion supported by Ochterski et al. .
The calculations by GS neglected relativistic contributions, which were very recently considered by Collins and Grev (CG) . Using relativistic (Douglas-Kroll ) coupledcluster methods, these authors found that the TAE of SiH4 contains a relativistic contribution of −0.67 kcal/mol. Combined with the earlier calculations of GS, this yields ◦ ∆Hf,0 [Si(g)]=107.4±0.5 kcal/mol, within Desai’s reduced error bar. However, as discussed there , the experimental data for silane, SiH4, involve an ambiguity. The JANAF heat of formation of silane, 10.5±0.5 kcal/mol is in fact the Gunn and Green  measurement of 9.5 kcal/mol increased with a correction  of +1 kcal/mol for the phase transition Si(amorphous)→Si(cr), which was considered an artifact of the method of preparation by Gunn and Green. If one were to accept their argument, the GS and CG calculations on SiH4 ◦ would actually support the original JANAF/CODATA ∆Hf,0 [Si(g)].
No such ambiguities exist for tetraﬂuorosilane, SiF4, for which a very accurate experimental heat of formation has been determined  by direct combination of the pure elements in their respective standard states in a ﬂuorine bomb calorimeter. Johnson’s  heat of formation at 298.15 K, −386.18±0.11 kcal/mol, is slightly higher in absolute value and slightly more precise than the CODATA value of −386.0±0.2 kcal/mol, itself based on an earlier experiment from the same laboratory .
Clearly, if a benchmark quality (preferably ±0.3 kcal/mol or better) TAE[SiF4 (g)] could ◦ be calculated, then an unambiguous redetermination of ∆Hf,0 [Si(g)] would be possible. Our previous study on BF3 being at the limit of the then available computational hardware, a similar study on SiF4 — which contains an additional heavy atom and eight additional valence electrons, leading to an expected increase in CPU time and memory requirements by a factor of about 3.7 (see below) — could only be completed most recently, and is reported in the present contribution.
Most electronic structure calculations reported here were carried out using MOLPRO 97.3  running on SGI Octane and SGI Origin 2000 minisupercomputers at the Weizmann Institute of Science. The very largest calculation, a full-valence coupled-cluster calculation involving 620 basis functions, was carried out on the National Partnership for Advanced Computational Infrastructure CRAY T90 at the San Diego Supercomputer Center.
As in our previous study on BF3, all electron correlation calculations involved in determining the valence and inner-shell correlation contributions to TAE were carried out using the CCSD  and CCSD(T) [16,17] coupled-cluster methods. (For the energies of the constituent atoms, the deﬁnition of Ref.  for the open-shell CCSD(T) energy was employed.) Both the very low T1 diagnostic  of 0.012, and inspection of the largest coupled-cluster amplitudes, suggest a system essentially totally dominated by dynamical correlation. From experience it is known  that CCSD(T) yields results very close to the exact (full conﬁguration interaction) basis set correlation energy under such circumstances.
Basis set limits for the SCF and valence correlation limits were extrapolated (see below for details) from calculated results using the (A)VTZ+2d1f, (A)VQZ+2d1f, and (A)V5Z+2d1f basis sets. For silicon, those basis sets consist of the standard Dunning correlation consistent [20,21] cc-pVTZ, cc-pVQZ, and cc-pV5Z basis sets augmented with two high-exponent d and one high-exponent f functions with exponents obtained by progressively multiplying the highest exponent already present by a factor of 2.5. The addition of such ‘inner shell polarization functions’  has been shown [22–25] to be essential for smooth basis set convergence in second-row compounds, particularly those containing highly polar bonds such as SiF4 . (It should be recalled that inner shell polarization is a pure SCF eﬀect and bears little relationship to inner shell correlation. In the present case of SiF4, the contribution of the inner polarization functions to the SCF/(A)VTZ+2d1f TAE was found to be no less than 9.81 kcal/mol.) For ﬂuorine, the basis sets given correspond to Dunning (diﬀuse function)-augmented correlation consistent  aug-cc-pVTZ, aug-cc-pVQZ, and aug-cc-pV5Z basis sets — it was shown repeatedly (e.g. ) that the use of augmented basis sets on highly electronegative elements such as F in polar compounds is absolutely indispensable for accurate binding energies. The ﬁnal basis sets for SiF4 involve 235, 396, and 620 basis functions, respectively, for (A)VTZ+2d1f, (A)VQZ+2d1f, and (A)V5Z+2d1f.
The geometry of SiF4 was optimized by repeated parabolic interpolation at the CCSD(T)/cc-pVQZ+1 level, where the suﬃx ‘+1’ stands for the addition of a tight d func
derive re is available.) The inner-shell correlation contribution was determined by comparing the computed binding energies correlating all electrons except Si(1s), and correlating only valence electrons, using the MTsmall basis set . The latter is a variant of the Martin-Taylor core correlation basis set [31,32] in which the very tightest p, d, and f functions were deleted at no signiﬁcant loss in accuracy on the contributions to TAE.
The scalar relativistic contributions were obtained as expectation values of the ﬁrst-order Darwin and mass-velocity operators [33,34] at the ACPF (averaged coupled-pair functional ) level using the MTsmall basis set. All electrons were correlated in this calculation, and it should be noted that the MTsmall basis set is completely uncontracted and therefore ﬂexible enough in the s and p functions for this purpose. For the sake of illustration, this approach yields −0.67 kcal/mol for SiH4, identical to two decimal places with the more rigorous relativistic coupled-cluster value .
The contribution of atomic spin-orbit splitting derived from the experimental atomic ﬁne structures  of Si(3 P ) and F(2 S) is −1.968 kcal/mol. For comparison, we also carried out all-electron CASSCF/CI spin-orbit calculations  using the spdf part of a completely uncontracted aug-cc-pV5Z basis set, augmented with a single tight p, three tight d, and two tight f functions in even-tempered series with ratio 3.0. In this manner, we obtain a contribution of −1.940 kcal/mol. In short, to the accuracy relevant for this purpose it is immaterial whether the computed or the experimentally derived value is used.
The zero-point energy was obtained from the experimentally derived harmonic frequencies and anharmonicity constants of McDowell et al. . This leads to a value of 8.029 kcal/mol, whereas one would obtain 8.067 kcal/mol from one-half the sum of the harmonic
All relevant data are given in Table 1.
As expected, the SCF contribution of TAE converges quite rapidly. We have shown previously  that the SCF convergence behavior is best described by a geometric extrapolation A + B/C n of the type ﬁrst proposed by Feller , with extrapolation from the TAE contributions to be preferred over extrapolation from the constituent total energies. From the (A)VTZ+2d1f, (A)VQZ+2d1f, and (A)V5Z+2d1f results, i.e. Feller(TQ5), we obtain a basis set limit of 448.43 kcal/mol, 0.02 kcal/mol more than the SCF/(A)V5Z+2d1f result itself. An extrapolation from the (A)VDZ+2d, (A)VTZ+2d1f, and (A)VQZ+2d1f basis sets would have yielded 448.47 kcal/mol, an increment of 0.22 kcal/mol over the (A)VQZ+2d1f result.
Given the large number of valence electrons, connected triple excitations account for a rather small part of the binding energy: 9.61 kcal/mol at the CCSD(T)/(A)VQZ+2d1f level, compared to a CCSD valence correlation contribution of 114.85 kcal/mol and an SCF contribution of 448.25 kcal/mol. Since a CCSD(T)/(A)V5Z+2d1f calculation is beyond the limits particularly of memory and available CPU time for this system, this suggests an approach in which only the CCSD valence correlation contribution be obtained from the largest basis set, while the (T) contribution is obtained from an extrapolation on smaller basis sets. Indeed, Martin and de Oliveira (MdO) recently found in a systematic study  on a wide variety of ﬁrst-and second-row molecules that this essentially does not aﬀect the quality of the results, except when the (T) contribution is a dominant component to the binding energy. Helgaker and coworkers  previously noted the more rapid basis set convergence behavior of connected triple excitations as compared with the CCSD correlation energy.
The CCSD/(A)V5Z+2d1f calculation required over 3GB of memory, some 120 GB of disk space, and 43 hours of real time (81 hours of CPU time) running on 8 CPUs of the NPACI CRAY T90. (Close to 99% parallellism was achieved in the CCSD code simply by adapting it to use vendor-supplied parallel BLAS and LAPACK libraries.) To our knowledge, this is the largest coupled-cluster calculation ever carried out using a conventional algorithm.
We have considered two extrapolation formulas based on the asymptotic behavior of pair correlation energies [41,42], namely the 3-point extrapolation A + B/(l + 1/2)α due to Martin, and the 2-point extrapolation A + B/l3 formula due to Helgaker and coworkers .
(In both formulas, l stands for the maximum angular momentum present in the basis set.) MdO found  that both formulas tend to predict the same basis set limit if extrapolated from suﬃciently large basis sets, but that the limits predicted by the A + B/l3 formula are much more stable with respect to reduction of the sizes of the basis sets used in the extrapolation. This is at least in part related to the fact that the three-point extrapolation involves, of necessity, one value with an even smaller l than the two-point extrapolation.
As an illustration, let us consider the BF diatomic which was used to reﬁne the BF3 result . From the three-point A + B/(l + 1/2)α extrapolation applied to AVnZ (n=3,4,5) valence correlation contributions to De, we obtain 38.35 kcal/mol, compared to 38.76 kcal/mol for AVnZ (n=4,5,6). In contrast, a A + B/l3 extrapolation applied to AVnZ (n=Q,5) yields
38.78 kcal/mol, just like AVnZ (n=5,6) does; application to AVnZ (n=T,Q) results yields an overestimate of 39.08 kcal/mol.