Relativistic Four-Component DFT Calculations of Vibrational Frequencies

We investigate the effect of relativity on harmonic vibrational frequencies. Density functional theory (DFT) calculations using the four-component Dirac–Coulomb Hamiltonian have been performed for 15 hydrides (H2X, X = O, S, Se, Te, Po; XH3, X = N, P, As, Sb, Bi; and XH4, X = C, Si, Ge, Sn, Pb) as well as for HC≡CPbH3. The vibrational frequencies have been calculated using finite differences of the molecular energy with respect to geometrical distortions of the nuclei. The influences of the choice of basis set, exchange–correlation functional, and step length for the numerical differentiation on the calculated harmonic vibrational frequencies have been tested, and the method has been found to be numerically robust. Relativistic effects are noticeable for the heavier congeners H2Te and H2Po, SbH3 and BiH3, and SnH4 and PbH4 and are much more pronounced for the vibrational modes with higher frequencies. Spin–orbit effects constitute a very small fraction of the total relativistic effects, except for H2Te and H2Po. For HC≡CPbH3 we find that only the frequencies of the modes with large contributions from Pb displacements are significantly affected by relativity.


■ INTRODUCTION
For molecules containing heavy atoms, relativistic effects play a crucial role in their electronic structure and chemical bonding. 1 Relativistic effects are commonly separated into scalar relativistic effects, which are due to (among other contributions) the mass−velocity and Darwin corrections, and the effects due to the spin−orbit interaction. The former lead for instance to contraction of the inner-shell orbitals (the energies of core levels are lower than those for the nonrelativistic case), and the latter result in the spin−orbit splitting of molecular orbital energy levels. Furthermore, the contraction of the inner-shell orbitals in turn increases the screening of the nuclear charge for the outer-shell electrons, giving rise to an indirect effect that results in expansion of the valence orbitals. These relativistic effects affect the valence orbitals involved with chemical bonding and consequently the potential energy surfaces. 1 In most cases where potential energy surfaces are concerned, it is sufficient to account for scalar relativistic effects using for example effective core potentials, 2 but for systems where strong spin−orbit effects may be expected, it is important to have an apparatus to calculate total relativistic effects using the four-component Dirac−Coulomb (or Dirac−Coulomb−Breit) Hamiltonian. Many four-component calculations have been carried out for dissociation energies 1,3 and molecular gradients (first derivatives of the molecular energy with respect to distortions of the nuclei in the molecule) 4 as well as equilibrium geometries. 4,5 In this contribution, we present the results of calculations of the molecular Hessian (second derivatives of the molecular energy with respect to nuclear distortions) and harmonic vibrational frequencies with the Dirac−Coulomb Hamiltonian. Four-component methods for the analytic calculation of molecular Hessians (and thus also harmonic vibrational frequencies) are currently not available in any computational chemistry program package. Instead, our calculations have been carried out with an external driver to the existing program package DIRAC. 6 Test calculations have been performed for hydrides of elements from groups 14 (XH 4 , X = C, Si, Ge, Sn, Pb), 15 (XH 3 , X = N, P, As, Sb, Bi), and 16 (H 2 X, X = O, S, Se, Te, Po) and in addition for the acetylene derivative HCCPbH 3 . Vibrational frequencies have been computed with the use of both relativistic and nonrelativistic methods in order to study the importance of the relativistic effects. Such calculations have previously been reported for halogen diatomics 7 but not for polyatomic molecules.
There are well-established methods of calculating potential energy surfaces, including the molecular Hessian and vibra-tional frequencies, for molecular systems with substantial relativistic effects through the use of the zeroth-order regular approximation Hamiltonian, 8,9 other two-component Hamiltonians, 10−13 and relativistic effective core potentials. 2 In most cases these are sufficient for rendering the relativistic effects, apart from some systems with very strong spin−orbit coupling, such as some lanthanide compounds. 1,14 However, a fourcomponent protocol will be useful for benchmarking more approximate treatments of relativistic effects. We have recently also demonstrated that the geometry dependence of NMR spin−spin coupling constants depends more strongly on relativistic effects than the spin−spin coupling constants themselves. 15 This suggests that relativity may be important for zero-point vibrational (ZPV) corrections to NMR properties. For properties such as spin−spin coupling constants, a full relativistic treatment is necessary, and it is therefore important also to develop tools that allow vibrational frequencies to be calculated at the full four-component level of theory.
■ METHODS Numerical Derivatives. Our program works as an external driver to the DIRAC program package. 6 The method for calculating the molecular Hessian and thus also the harmonic vibrational frequencies is fully numerical. Computation of the Hessian is based on calculating the second derivatives of the molecular energy with respect to geometric distortions of the molecule using simple three-point formulas: 16 This means that the Hessian computation involves performing a number of energy calculations in which atoms are displaced from their original positions in all degrees of freedom. In the case of a molecule with N atoms, 18N 2 + 1 single-energy computations need to be run to determine the full Hessian. Once the Hessian is obtained, the vibrational frequencies are calculated by diagonalization of the Hessian in its massweighted form. When numerical differentiation is performed, it is important that an appropriate step length (h in the above equation) is used to ensure numerically accurate results. We performed test calculations of the vibrational frequencies for the water molecule with a number of different step lengths in the range of 10 −1 −10 −5 Å. The calculations turned out to be numerically stable for step lengths between 10 −2 −10 −4 Å. Similar test calculations for the H 2 Po molecule revealed that in the case of heavier atoms, numerical stability shifts in the direction of larger step lengths (5 × 10 −2 to 5 × 10 −3 Å). All of the results can be found in the Supporting Information. On this basis, in all of the subsequent calculations we used a step length of 10 −3 Å. The only exceptions were systems involving the heaviest atoms (Pb, Bi, and Po), for which we used a step length of 10 −2 Å. Using the differences between vibrational frequencies obtained with step lengths within the range of numerical stability, we were able to estimate the error bars to be about 3 cm −1 .
In order to test the numerical stability of the three-point formula, additional calculations were carried out using a fivepoint formula: 16 No significant differences were found in comparison with the three-point formula with the same step length.
Geometry optimizations were performed using the DIRAC program 6 at the same level of theory as for the Hesssian calculations carried out afterward in order to ensure that the molecular gradient was equal to zero (a condition for the harmonic approximation). The convergence treshold for the gradient was 10 −4 .
Single-Energy Calculations. The four-component Dirac−Kohn−Sham (DKS) energy calculations were carried out with the DIRAC program. 6 The uncontracted aug-cc-pVDZ basis set 17  Moreover, we carried out an investigation of the dependence of the results on the choice of exchange−correlation functional and basis set. In order to do so, four-component DKS calculations were performed using the PBE0 functional 26 (to be compared to B3LYP) and also two additional basis sets: (1) the uncontracted aug-cc-pVDZ basis set 17 on the hydrogen atoms and Dyall's uncontracted double-ζ basis set (dyall.v2z) 19,20 on all of the other atoms and (2) the uncontracted aug-cc-pVQZ basis set 17 on the hydrogen atoms and Dyall's uncontracted quadrupole-ζ basis set (dyall.v4z) 19,20 on all of the other atoms.
The convergence threshold for all of the single-energy calculations was 10 −6 .
Since an analytical method for calculating the molecular Hessian is implemented in the DALTON program, 27,28 some calculations were performed with this program for comparison. We note that DALTON allows only one-component nonrelativistic DFT calculations. All of the DALTON computations were run with the same uncontracted basis set and exchange− correlation functional as above and were carried out using the geometry optimized in DIRAC at the nonrelativistic level of theory (the same geometry as the one used for nonrelativistic numerical calculations of vibrational frequencies).

■ RESULTS AND DISCUSSION
In order to test the newly developed method for calculating harmonic vibrational frequencies, simple systems consisting of three, four, or five atoms have been chosen: • H 2 X where X = O, S, Se, Te, Po; • XH 3 where X = N, P, As, Sb, Bi; • XH 4 where X = C, Si, Sn, Pb.  Table 1. In most cases the frequencies obtained with PBE0 are larger, but the differences between the results obtained with PBE0 and B3LYP do not exceed 3% in any case. Taking this into consideration, it seems that in this case B3LYP and PBE0 would produce comparable results. The B3LYP functional has been chosen for the following calculations because of its good performance for vibrational properties reported in the literature. 29−31 The results of the four-component DFT calculations of vibrational frequencies for H 2 X systems, carried out with double-ζ, triple-ζ, and quadruple-ζ quality basis sets, can be found in Table 2. The differences between vibrational frequencies obtained with these three basis sets are almost negligible. The biggest differences occur between DZ and QZ for the H 2 O molecule, yet even in this case these differences are not larger than 1% of the values, being at most 31 cm −1 . In all other cases, the differences do not exceed 10 cm −1 . In light of the above findings, the triple-ζ-quality basis set appears to be an optimal compromise between accuracy and computational cost, and this basis set has been used in all of the following calculations.
Numerical versus Analytical Hessian. As numerical methods for calculating the molecular Hessian will inevitably have limitations on the numerical accuracy, we have tried to estimate these by comparing the numerical harmonic vibrational frequencies with the results obtained with the analytical nonrelativistic method implemented in the DALTON program. The comparison of the calculated harmonic vibrational frequencies can be found in Table 3. We obtained excellent agreement in case of the H 2 O, H 2 S and H 2 Se molecules.
Influence of Relativity on the Vibrational Frequencies. Harmonic vibrational frequencues calculated with the relativistic and nonrelativistic methods are summarized in Tables 3 − 5. As can be noted, in almost all cases the relativistic vibrational frequencies are smaller than the corresponding nonrelativistic values, that is, relativity decreases   the bond strength. In the case of the H 2 X systems, relativistic effects are significant for H 2 Te (2% for ω 1 and ω 2 ) and H 2 Po (10% for ω 1 and ω 2 , 3% for ω 3 ). Also in the case of the XH 3 molecules, relativistic effects are not negligible for the two heaviest congeners, 2% for ω 1 and ω 2 in SbH 3 and 8% for ω 1 , 5% for ω 2 , and 2% for ω 3 in BiH 3 . In the case of the XH 4 molecules, relativistic effects are negligible for all but one mode, ω 3 (E symmetry mode), for SnH 4 (2%) and PbH 4 (6%).
It should be stressed here that the precentage change in the values when calculated with relativistic and nonrelativistic methods varies for each vibrational mode.
In addition, for all of the molecules, relativistic fourcomponent calculations without spin−orbit effects have been performed (Table 4). These results show that in the case of XH 3 , all of the relativistic effects are in fact scalar relativistic effects, whereas in the case of H 2 X, spin−orbit effects play a B3LYP functional, aug-cc-pVTZ (on H) + dyall.v3z (on X) basis set. b No symmetry has been used, so frequencies of degenerate vibrations vary (by at most 2 cm −1 ). Arithmetic averages are given. c Fundamental vibrational frequencies are reported for the experimental data. For the calculated results, the following notation is used: rel, relativistic; nrel, nonrelativistic; no SO, no spin−orbit coupling. Comparison with Experimental Values. When comparing the results obtained with the experimental values, one should keep in mind that the diagonalization of the molecular Hessian gives us harmonic vibrational frequencies. Thus, anharmonicity is not taken into account, and this will lead to some difference between the results and the experimental values. In our case, the differences do not exceed 5%, in line with the expected magnitude of anharmonic corrections. 35 Vibrational Frequencies for HCCPbH 3 . To illustrate the usefulness of the presented method and to study the effects of relativity on vibrational frequencies for a more complex system where only some of the modes involve the heavy atom, we have calculated the vibrational frequencies for the acetylene derivative HCCPbH 3 with both relativistic and nonrelativistic approaches. The motivation for choosing this particular system was our previous work, 15 where we showed that for this molecule the relativistic effects on the derivatives of the indirect spin−spin coupling constants with respect to molecular geometry parameters tend to be more pronounced than the effects on the coupling constants themselves. The ZPV corrections calculated at the nonrelativistic level are therefore not necessarily reliable. The uncontracted aug-cc-pVDZ basis set 17 was used on the hydrogen and carbon atoms and Dyall's uncontracted triple-ζ basis set (dyall.v3z) 18−20 on the lead atom have been applied together with the B3LYP exchange−correlation functional. 21−24 The results are collected in Table 6.
In the case of the HCCPbH 3 molecule, it seems that only vibrations that involve the Pb atom are significantly affected by relativity. There is almost no difference between the relativistic and nonrelativistic values of vibrational frequencies for C−H stretching, C−C stretching, and C−C−H bending. This finding may be useful for future calculations of vibrational effects on molecular properties for large molecules, since it may allow for relativity to be taken into account only for selected localized modes. Similar findings were previously reported by Berger et al. 36 In addition to this, the relativistic effects are much more pronounced for deformation modes than for stretching modes. The difference between the vibrational frequencies calculated with relativistic and nonrelativistic methods for the C−C−Pb bend exceeds 20% of the value, whereas for the C−Pb stretch it is only little more than 5%.

■ CONCLUSIONS
We have presented a numerical method for calculating the molecular Hessian and harmonic vibrational frequencies with relativistic four-component DFT. Test calculations have been performed for hydrides of elements from groups 14, 15, and 16. We have achieved good agreement with an analytical nonrelativistic DFT method.
Relativistic effects become significant primarily for the hydrides containing atoms from the fifth and sixth rows of the periodic table and are much more pronounced for the vibrational modes with higher frequencies. Spin−orbit effects constitute a very small fraction of the relativistic effects on the whole, with the exception of H 2 Te and H 2 Po. Additional calculations for HCCPbH 3 show that only the frequencies of the modes with large contributions from Pb displacements are significantly affected by relativity.
This work is considered a stepping stone towards the development of a four-component relativistic numerical method for calculating ZPV corrections to NMR parameters (spin−spin coupling constants and shielding constants).
Vibrational frequencies for H 2 O and H 2 Po obtained with different step lengths and different optimized geometries (PDF)