- Email: [email protected]
Vibrational effects in the parity-violating contributions to the isotropic nuclear magnetic resonance chemical shift
Chemical Physics Letters 470 (2009) 166–171
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Vibrational effects in the parity-violating contributions to the isotropic nuclear magnetic resonance chemical shift Ville Weijo a,b,*, Mikkel Bo Hansen c, Ove Christiansen c, Pekka Manninen d a
Laboratory of Physics, Helsinki University of Technology, P.O. Box 1100 (Otakaari 1 M), 02015 Espoo, Finland Laboratory of Physical Chemistry, Department of Chemistry, P.O. Box 55 (A.I. Virtasen Aukio 1), University of Helsinki, FI-00014 Helsinki, Finland c The Lundbeck Foundation Center for Theoretical Chemistry and Center for Oxygen Microscopy and Imaging, Department of Chemistry, University of Aarhus, Langelandsgade 140, DK-8000 Århus C, Denmark d CSC-IT Center for Science Ltd., P.O. Box 405, FI-02101 Espoo, Finland b
a r t i c l e
i n f o
Article history: Received 5 May 2008 In final form 13 January 2009 Available online 19 January 2009
a b s t r a c t We investigate the effect of vibrational corrections to the parity-violating (PV) contributions to the nuclear magnetic resonance shielding constant for the CHFClBr molecule using a non-relativistic framework. Density-functional theory is used for the electronic structure part and vibrational configuration interaction for the vibrational part. Zero-point vibrational corrections are found to be of the order of a less than 10% with respect to the PV contributions calculated at the equilibrium geometry in relevant cases. Vibrational corrections are evaluated also for the fundamental vibrations and are found to be up to 42% of the total PV contribution in the heavy elements. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction The first formulations of the parity-violating (PV) theory of nuclear magnetic resonance (NMR) spectral parameters were presented by Barra et al. [1–3]. The first ab initio calculations of PV contributions to the isotropic NMR chemical shift were made for small molecules, such as H2O2 [4,5], more than a decade later. Advanced computational aspects, like electron correlation effects, basis sets, and special relativity have been explored in recent years [6–8] and all of those aspects have been found to have significant effect on the calculated values. In addition to NMR shielding constants, PV contributions to the NMR spin–spin coupling constants have been evaluated [9], also electron spin resonance spectroscopy has been considered as a possible method for observing PV effects in molecules [10]. As for now, all predicted effects for the spectroscopical methods mentioned above have been well below the experimental limits. The effects of molecular vibrations have, on the other hand, so far not been considered in magnetic resonance properties. Tools for molecular vibrational calculations have been developed in particular in recent years (see, e.g., Refs. [11–13] and references therein). Vibrational methods usually treat electronic and vibrational wave functions separately in the spirit of the adiabatic (Born–Oppenheimer) approximation. Computations then proceed by using electronic structure methods to generate points on the potential energy surface (PES) and the property surfaces. These * Corresponding author. Present address: Department of Chemistry, University of Tromsø, N-9037 Tromsø, Norway. E-mail address: [email protected] (V. Weijo). 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.01.022
points can then be used for constructing an intermediate representation of the PES and property surfaces as function of the relative nuclear positions. Subsequently, a wave function for the nuclear motion is obtained from the PES by solving approximately the Schrödinger equation for the motion of the nuclei. Then using this wave function and the generated property surfaces, the vibrational averaging can be carried out. Different variants of this strategy vary in their selections of PES generation options, vibrational wave function method, and vibrational basis sets. We note that the most widely used approach for calculating anharmonic vibrational averages is probably still the simplest vibrational perturbation approach, circumventing the explicit construction of vibrational wave functions. Explicitly constructing vibrational wavefunctions may be considered to be more robust in many ways but often also comes with an increased cost from typically requiring PES information over a wide range of structures. Moreover, the calculations of the vibrational wave functions can become computationally demanding. In this Letter, we will use explicit vibrational wave function methods to calculate vibrational corrections to the PV isotropic NMR shielding constant contributions of a CHFClBr molecule. The molecule has been used previously in calculations of PV shielding and nuclear spin–spin coupling contributions [8]. Also, the molecule has been considered for a candidate for observation of PV effects in vibrational spectroscopy (see, e.g., Refs. [14–18] for theory and Refs. [19–21] for experiments). We focus here on CHFClBr rather than H2O2, which has been used more commonly in PV NMR calculations [4–8]. Although CHFClBr is not experimentally well-suited for detection of PV effects in NMR due to the high-spin nuclei in it, CHFClBr is a better analogue to a hypothetical
167
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171
heavy-element-containing experimentally relevant molecule than H2O2. The latter is also significantly more difficult to address using automated computational tools due to the presence of a torsional HOOH mode giving rise to a double minimum potential. Vibrational corrections to the property in question vary significantly depending on the nucleus in question and the amount of vibrational excitation. Due to the relative size of the corrections, they cannot be neglected in high-accuracy calculations.
2. Computational details The potential and property surfaces used in this study has been calculated using density-functional theory with the Becke threeparameter exchange-correlation hybrid functional (B3LYP) and Dunning’s augmented correlation-consistent triple-zeta basis set (aug-cc-pVTZ). Vibrational calculations were made using MIDASCPP vibrational structure program [22] calling the DALTON electronic structure program [23] for energies and properties at given points. Calculated surfaces and structures are available from the authors on request. The nuclei used in vibrational and response calculations were 1 H, 13C, 19F, 35Cl and 79Br. The molecule is in its (R)-enantiomeric form. Geometries and potential/property surfaces are available from the authors on request. 2.1. Vibrational structure calculations The representation of the potential energy surfaces (PES) and the property surfaces is a central point in vibrational structure theory and still an issue of much contemporary research. We shall give a brief discussion of the strategy we use in this work, referring to a number of reviews and recent research papers for more detailed accounts [11–13,24,25]. The surfaces can be expanded in a series of n-mode potential terms,
V
ð1Þ
;V
ð2Þ
;V
ð3Þ
ð1Þ
;...;
(1)
K
nates, rather than the full Watson operator and hence we neglect vibration–rotation interactions in this study.
(2)
where V includes only one-mode terms, while V includes all terms with up to two-mode couplings and so on. If one does not truncate the series, one obtains the fully coupled potential. The potential energy surfaces are generated as a multi-dimensional high-degree polynomial fit to a set of grid points [24]. The outline of the procedure is as follow. First the molecular structure is optimized and a vibrational analysis carried out. The vibrational analysis provides harmonic frequencies and mass-weighted normal coordinates which are subsequently used in the MIDASCPP program to determine the boundaries of the grid surface as well as the positions of the coarse grid points (CGP). The grid cut-offs are determined from the classical turning points for each vibrational mode (in mass-weighted normal coordinates),
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2h 1 xtp ¼ vþ ; 2 x
Prior to the fitting a number of additional points may be generated by interpolation via a direct product natural N-cubic spline. We used a fine mesh of 12 per coarse grid point keeping only the inner most 90% to avoid border effects. The resulting points were fitted to a polynomial of degree 12. More commonly used method based on numerical differentiation and Taylor expansions of the potential and molecular property surfaces has been also used. The step-sizes used for the individual modes in the numerical differentiation procedure are scaled by the corresponding harmonic frequencies multiplied with the same constant factor. In this study a factor of 0.070 has been used. Surfaces generated this way are denoted mMtT, where m denote the maximum mode-coupling level in the energy/property potential and t the degree of the Taylor expansion. For example, 2M4T means that all terms up to fourth-order depending on one- and two-mode couplings are included. At this point we emphasize that PESs are obtained via the grid approach, whereas property surfaces are based on fourth-order Taylor expansions. The vibrational wave function is constructed in terms of Hartree products of excitations of distinguishable modes. The one-mode basis is expanded in a certain number of harmonic oscillator (HO) functions. The simplest such vibrational wave function consists of a single Hartree product, which is optimized using the vibrational self consistent field method (VSCF) to minimize the energy (see Refs. [11–13] and references therein). The VSCF method treats the interaction between the modes in an averaged fashion, and to alleviate this, the vibrational configuration interaction (VCI) method is used. In the VCI calculations we excite simultaneously up to two, three, or four modes (VCI[2], VCI[3], and VCI[4], respectively). In this work, we use a state-specific approach. Thus, all states have their own VSCF reference state, and the excitation space is generated relative to that VSCF state. The vibrational structure calculations employ a simple kinetic P @2 energy operator, i.e. 12 K @Q 2 in mass-weighted normal coordi-
ð2Þ
where v is a vibrational quantum number, in this study a value of v = 10 is used. The grid is then divided into a number of equal-sized domains in the one-mode, two-mode, etc. parts, thus determining the coarse grid points, which are subsequently calculated via an interface to an electronic structure program. The number of grid points needed in each dimension typically decreases with increasing mode-combination level, i.e., all one-dimensional cuts are determined from a large number of coarse grid points, all twodimensional cuts from a smaller number of points in each dimension, and so forth. A 641161 grid is based on 64 CGPs in each 1D cut and 162 CGPs in each 2D cut.
2.2. Property evaluation The leading-order parity-violating contribution to the isotropic NMR shielding tensor in the non-relativistic framework can be expressed as a linear response function (see, e.g., Ref. [9])
rPV K;s ¼
1 DD
cK
PV
OZ
hK; ; hB0 ;s
EE 0
;
ð3Þ
which combines the orbital Zeeman operator OZ
hB0 ;s ¼
1X ‘iO;s ; 2 i
ð4Þ
and the nuclear spin-dependent part of the total parity-violating operator
X GF a PV 2 hK; ¼ pffiffiffi kK ð1 4 sin hW Þ iri; ; dðriK Þ þ : 2 2 i
ð5Þ
Here, K refers to the nucleus in question and and s are Cartesian indices of the tensor. ‘iO is the angular momentum operator of electron i with respect to the gauge origin O, and riK is the distance vector between i and K. The SI-based atomic unit system is used throughout this work. The weak coupling constant is denoted GF = 2.22 1014 a.u. and the Weinberg angle is hW = 0.2320. kK is a nucleus-dependent factor that includes contributions from the nuclear structure and the nuclear anapole moment. The kK factors are not usually known experimentally or theoretically and are set to unity during the calculations. If accurate values will become
168
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171
available, results can be scaled accordingly. This issue also does not change the relative importance of vibrational and electronic contributions. There is also another leading-order contribution, but it does not contribute to the isotropic part of the tensor [9]. It can be shown that the combined leading-order term is gauge-independent at the complete basis set limit [5], and thus, although there exists a reference to the gauge origin in Eq. (4), the isotropic part should be free of problems in this respect. Vibrational contributions to molecular properties fall into two categories: (i) the motion of the nuclei may respond to an external perturbation, thus giving rise to pure vibrational terms as for example (hyper)-polarizabilities [13,26], and/or (ii) the property of interest may depend on the nuclear geometry, and hence one have to average the property over at least the vibrational ground state, giving rise to so called zero-point vibrational averages (ZPVA). In this Letter, we are dealing solely with issue (ii). Contribution (i) is zero in our concrete case, since the electronic expectation value of the orbital Zeeman operator 4 is zero for the closedshell molecule considered here. The ZPVA is defined as
hAiZPVA ¼ hUjAjUi;
ð6Þ
where jUi may be either a VSCF or a VCI state and A may be any property that is a function of nuclear coordinates and expanded as described above [24,27–29]. In this work, we concentrate on the vibrational correction, which is the difference between the vibrational average and the equilibrium geometry results. 3. Results and discussion Harmonic zero-point energy (ZPE) and vibrational frequencies of CHFClBr are listed in Table 1. Our values for harmonic frequencies of CHFClBr are mostly in line with the results of Rauhut et al. [30]. They made their calculations with the same aug-cc-pVTZ basis set and B3LYP functional, but they used the 12C isotope instead of 13C. After scaling the harmonic frequencies of stretching modes by a square root of their reduced mass ratios, our results agree with Rauhuts within 3 cm1, with the exception of the C–Br stretch. The difference can be explained by their use of a relativistic effective core potential (RECP) for Br. Although RECPs are often very useful for heavy elements like Br to take care of relativistic effects, the way they were implemented in the currently used program package made them unsuitable in this case, because PV contributions are very sensitive to the description of the electronic structure in the vicinity of the nucleus.
Table 1 Zero-point energies (ZPE) and fundamental excitation frequencies (in cm1) of 13 1 19 35 79 C H F Cl Br evaluated using the harmonic approximation and also results using different PESs are included (VCI[2] calculation with the 6 HO vibrational basis set). The assignment of the modes follow the notation of Ref. [30]. PES
Harm.a
2M4T
32181
64181
641161
64116141
Assignment
ZPE
4547.5 3143.3 1318.6 1211.2 1049.3 720.9 622.7 415.6 306.3 218.9
4453.3 3005.3 1285.9 1180.8 1023.0 708.0 611.1 411.4 303.8 215.4
4443.6 3006.6 1275.1 1169.5 1022.5 707.2 611.2 411.3 303.5 215.4
4444.7 3006.6 1275.1 1169.6 1022.5 707.2 611.1 411.3 303.5 215.4
4453.4 3012.2 1284.1 1181.5 1022.2 707.2 610.6 411.3 303.5 215.4
4453.0 3013.0 1283.9 1179.7 1022.0 707.3 610.0 411.1 303.2 215.4
CH str. CHF sciss. CHBr/Cl sciss. CF str. CCl str. CBr str. CFCl sciss. CFBr sciss. CClBr sciss.
m1 m2 m3 m4 m5 m6 m7 m8 m9 a
Results using the harmonic approximation.
3.1. Anharmonic effects to vibrational frequencies Effects of anharmonicity calculated using different PESs are also shown in Table 1. It has previously been demonstrated that the use of low order Taylor expanded PESs in vibrational calculations may provide results of low accuracy. We have therefore chosen to also use grid based PESs in the calculations of anharmonic energies and wavefunctions. The 64181 grid is saturated with respect to one-mode couplings. The differences between 64181 and 641161 grids can be approximately 10 cm1 for the high-energy modes, but again for the low-energy modes, the differences are less than 1 cm1. 641161 is most likely close to convergence with respect to the two-mode couplings, whereas using just 8 points in either dimension of the two-mode cuts is too crude an approximation. The distance between the points simply becomes too large. The three-mode couplings do not change the situation much, the difference between 641161 and 64116141 grids is 2 cm1 at most (in mode m3) and often much less. This suggest that the PES series of Eq. (1) is very rapidly converging, being close to convergence already at the V(2) level. The fast convergence with respect to mode-coupling level was illustrated with state-specific VCI[2] wave functions, which could be considered not flexible enough. However, the convergence of ZPE and fundamental frequencies with respect to the vibrational correlation treatment is also rather rapid, as can be seen from Table 2 reporting calculations using the 64116141 potential. Even at the VSCF level, the low-energy modes do not differ from VCI[4] results by more than a couple of cm1. With ZPE and modes m2 and m3, the similar difference is at most 10 cm1. The only exception to the very good convergence is the mode m1 where a difference of 50 cm1 is observed. VCI[3] results are within 2 cm1 from the VCI[2] results with an exception of mode m6, where a slightly erratic behavior is seen. There, the VCI[2] treatment first reduces the frequency by 6 cm1 compared to VSCF, and subsequently, VCI[3] restores the value close to the VSCF result again. The VCI[4] calculations lead to additional changes to VCI[3] of an order of a few tenths of a cm1 with a slightly larger effect for m1. All states are dominated by a single configuration meaning they are physically well described as a fundamental. The highest lying fundamental, m1, obtains somewhat larger contributions to its wave function from other components. A one-mode basis with six HO basis functions is also very close to its respective limit, and the difference between six HO and eight HO basis sets is at most 1 cm1 (in mode m1) and is often much smaller (values not shown here). Again, any experimental reference values for fundamental frequencies are not available for 13C1H19F35Cl79Br. However, if we compared the general features of a vibrationally correlated gridbased calculation to its harmonic counterpart, we would observe a similar trend of lowering of frequencies as can be found for
Table 2 Effects of vibrational correlation to the ZPE and fundamental frequencies of 13 1 19 35 79 C H F Cl Br (in cm1). The PES is from the 64116141 grid and the 6 HO basis set is used. State
VSCF
VCI[2]
VCI[3]
VCI[4]
ZPE
4460.5 2965.5 1289.3 1188.7 1024.2 709.0 615.9 412.0 303.8 215.5
4453.0 3013.0 1283.9 1179.7 1022.0 707.3 610.0 411.1 303.2 215.4
4452.4 3012.9 1282.0 1178.7 1021.5 705.1 617.2 410.8 302.8 215.2
4452.4 3011.7 1281.6 1178.5 1021.4 704.9 615.0 410.7 302.8 215.2
m1 m2 m3 m4 m5 m6 m7 m8 m9
169
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171 12
CHFClBr in vibrational perturbation theory calculations [30]. There, the computational errors are at most 50 cm1 compared to the experimental values, and we expect that our results achieve comparable accuracy, even though RECPs are not used. 3.2. Vibrationally averaged parity-violating contributions PV contributions to the NMR shielding at the equilibrium geometry for all the nuclei in CHFClBr are presented in the column ‘Eq.’ of Table 3. Similar PV contributions were also evaluated in Ref. [9] for the same system with the same functional but different basis set. There are a couple of differences in the results. Excluding the 79 Br contributions, signs of the contributions are reversed, which can be explained by different signs of the k-factors used in Ref. [9]. The contributions calculated in this study (excluding 79Br) are somewhat smaller in magnitude than their counterparts calculated in Ref. [9]. This behavior can be linked to basis set effects. It is known that general purpose basis sets lack the proper flexibility for core electrons (see, e.g., Refs. [7,8]) in this property. Augmenting aug-cc-pVTZ basis set with steep s and p-type Gaussians can result to large differences in PV contributions [7]. The drawback is the extra computational effort, which can be very large, especially if many computations are required as is in our case. Core orbitals, however, are only slightly affected by molecular environment and their effect to the relative magnitudes between vibrationally uncorrected and corrected PV contributions should be small. Basis set effects are also responsible for the sign difference in 79Br contributions (after including the difference in k factors). Table 3 also includes vibrational corrections to the equilibrium values using different levels of vibrational correlation methods and different basis sets. The PES is chosen to be the 64116141 grid, which was shown above to be sufficiently close to convergence in the vibrational frequency calculations. The property surface is based on the 2M4T Taylor expansion, because this level offers considerable computational savings over the grid based methods in time-consuming response calculations. The convergence of the ground state vibrationally corrections is very robust, as can be expected from the convergence of ground state energies in Table 2. Differences in various levels of VCI and HO basis sets are very small overall. The 6 HO basis set can be considered to be fully converged as the difference between the 4 HO and 6 HO results is of the order of per miles at most. Similar differences can be found between the VCI[3] and VCI[4] results, and thus the VCI[4] calculations with the 6 HO basis set is in practice converged. At this point we note that previous studies, see Ref. [17] and references therein, estimate the PV effects arising from molecular vibrational motion from uncoupled potential/property surfaces. In order to test the quality of this approximation, we tested several levels of both the PES and property surfaces. It was found that the PV effects are highly dependent on the quality of the PES, which in turn determines the accuracy of the vibrational wavefunction. For example, using an uncoupled anharmonic PES with any property
surface give vibrational corrections to the PV contribution which are very far from the VCI[4] results using a proper PES. On the other hand, with a proper vibrational wavefunction, e.g. VCI[3], obtained from a good PES the effect of mode–mode coupling in the property surfaces is more limited. The absolute deviations of VCI[3] using uncoupled property surfaces and a good PES from the most accurate results presented in Table 3 amounts to about 12%, 46%, 11%, 1%, and 1% for the 13C, 1H, 19F, 35Cl, and 79Br nuclei, respectively. The deviation for the 1H nucleus is rather large. This is, however, believed to be of minor importance, as the PV effect is unlikely to be measured for this nucleus due to its small value. From the above study it is therefore problematic to give a guideline for including mode–mode-coupling in the property surfaces. One should, however, make sure that the underlying vibrational wavefunction is satisfactory, which in turn is reflected by the quality of the PES entering the Hamiltonian. Vibrational corrections to the vibrationally excited states behave similarly with respect to basis sets and vibrational correlation. Vibrational corrections for all modes are shown in Table 4. The results were obtained using the VCI[4] method and 6 HO basis set. As compared to the equilibrium geometry reference values, vibrational corrections to the ground state (GS) are smaller. The largest correction is for 1H, where the correction is approximately 20%, while the smallest vibrational corrections are for 13C (of the order 1%). The vibrational corrections to 19F, 35Cl, and 79Br shielding contributions are 2.4%, 6.6%, and 8.4% of the equilibrium geometry values, respectively. Hence, the vibrational corrections do not systematically depend on the mass of the nucleus. They are, however found to be non-negligible. The vibrational corrections are generally significant for the considered vibrationally excited states, although some variations are observed. Most of the 13C corrections are equal or larger in magnitude than the GS correction, and depending on the mode, increase or decrease the magnitude of the total contribution. Correction from the mode m3 has the largest value of 1.23 1011 ppm, which is approximately 13% correction to the equilibrium value. Most of the 1H corrections increase the magnitude of the total contribution. The largest correction is close to be as large (mode m8) as the equilibrium geometry value, which is clearly substantial. In the case of 19 F corrections, the magnitude of the total PV contribution is reduced by less than 10% at most (mode m1). Also, most of the corrections are of a different sign than the equilibrium geometry result, thereby reducing the magnitude of the total PV contribution. Likewise in 35Cl corrections the magnitude of the total contribution is often reduced by the vibrational averaging of excited states. The magnitudes of the corrections are within a factor of four from the GS correction with a maximum at the mode m3, where the correction decreases the magnitude of the equilibrium value by 21%. Vibrational corrections in 79Br are mostly negative (increasing the magnitude) although four out of nine (excited states) are positive. Magnitudes are, again, roughly of the size of the GS correction, but the largest of the corrections (mode m4) is rather substantial and provides a 41% increase in magnitude to the equilibrium geometry values.
Table 3 Vibrationally averaged ground state corrections to the parityviolating contributions to the isotropic NMR shieldings in CHFClBr using different vibrational correlation methods and basis sets. Equilibrium geometry reference values (Eq.) are also given. All values in ppm. Nucleus
13
C H 19 F 35 Cl 79 Br 1
VCI[2]
VCI[3]
VCI[4]
Eq.
4 HO
6 HO
4 HO
6 HO
4 HO
6 HO
9.45 1013 3.06 1014 2.81 1012 1.76 1010 1.43 1010
9.25 1013 3.07 1014 2.86 1012 1.78 1010 1.45 1010
9.79 1013 3.11 1014 2.78 1012 1.76 1010 1.45 1010
9.60 1013 3.12 1014 2.83 1012 1.78 1010 1.47 1010
9.69 1013 3.11 1014 2.82 1012 1.78 1010 1.47 1010
9.68 1013 3.11 1014 2.82 1012 1.78 1010 1.47 1010
9.48 1011 1.51 1013 1.20 1010 2.68 109 1.76 109
170
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171
Table 4 Vibrational corrections to the parity-violating contributions to the NMR shielding constants of nuclei in CHFClBr at the ground state (GS) and vibrationally excited states. Equilibrium geometry reference values (Eq.) are also given. Vibrational calculations were made at the VCI[4] level with the 6 HO basis set. All values in ppm. Mode Eq. GS
m1 m2 m3 m4 m5 m6 m7 m8 m9
13
1
C
19
H
11
9.48 10 9.68 1013 5.26 1012 1.48 1012 1.23 1011 7.52 1012 7.18 1012 1.56 1012 7.83 1013 3.69 1012 3.92 1012
35
F
13
1.51 10 3.11 1014 2.74 1014 8.14 1014 1.65 1014 5.31 1014 1.01 1016 6.30 1014 5.35 1014 1.19 1013 8.84 1015
Connection between the magnitude of a vibrational correction and motions in a particular mode is not a clear one. For example, the largest vibrational correction to the 79Br PV contribution can be found from the mode m4, which consists mostly of stretching motion between C and F. Based on this, one could conclude that large variations in total molecular geometry would result in large vibrational corrections. Then, one could expect CH stretching mode (m1) to have even larger corrections, but in this case, the vibrational correction is two orders of magnitude smaller than that of the mode m4. Similar lack of systematic correlation between mode types and vibrational corrections can be found from all the other PV contributions, too. If one considers the case of a room temperature experiment, then the most relevant modes are the GS, m9, m8, the first overtone of m9, and perhaps also the mode m7. Here we exclude 1H corrections from the discussion, because the total 1H PV contribution is orders of magnitude smaller than that of the other nuclei. Vibrational corrections from the modes mentioned above (except the overtone) range from 10% to +16% of the equilibrium values. For order of magnitude estimates of PV effects in NMR spectrum these corrections are not too significant, but an experimental setup will suffer from line width broadening and associated loss of peak intensity. Also, if an accurate experiment is realized and, e.g., k factors are to be evaluated using quantum chemical calculations, variations in PV contributions of the low-lying modes impose a serious computational challenge. Experimentally relevant molecules can be larger than CHFClBr and most likely will contain even heavier elements. Accurate electronic structure calculations themselves on such systems are tricky and time consuming, but if vibrational structure calculations are included into the process, then the computational effort will truly be massive for high-accuracy. In general, the CHFClBr molecule is quite rigid. Somewhat larger relative changes in PV contributions due to the vibrational corrections could be expected for floppy molecules such as H2O2, where the PV contribution varies significantly near the equilibrium geometry as a function of the dihedral angle [4–8]. We believe that our numbers with the B3LYP functional and the aug-cc-pVTZ basis set for the electronic structure part and VCI[4] with 6 HO basis set for vibrational part provide good estimates of the relative magnitudes of vibrational corrections. The non-relativistic Hamiltonian used in our calculations, on the other hand, leaves room for errors. Relativistic effects have been shown to increase the PV contributions by orders of magnitudes for heavy elements [6,8], and in some cases, the spin–orbit interaction can even change the sign of the contribution. At the equilibrium geometry, however, the effect is not nearly as pronounced, and the relativistic effects only increase the magnitude of the PV contribution by a factor of 1.3 in third row elements [8]. Also, we expect that the property surface near the equilibrium geometry in the relativistic case is qualitatively similar to the non-relativistic case, at least in rigid molecules like
10
1.20 10 2.82 1012 1.10 1011 7.14 1012 6.30 1012 2.70 1012 6.70 1012 4.24 1012 6.55 1012 1.77 1012 2.58 1012
79
Cl
Br
09
2.68 10 1.78 1010 2.14 1010 1.82 1010 5.74 1010 5.49 1010 3.66 1011 1.85 1010 2.83 1010 3.10 1010 2.67 1010
1.76 1009 1.47 1010 4.34 1012 5.04 1011 5.64 1010 7.35 1010 2.81 1011 6.91 1011 1.82 1010 8.02 1011 2.76 1010
CHFClBr. Thus, we expect our findings to be fairly generalizable to other rigid molecules containing even heavier elements. 4. Conclusions We have computed vibrationally averaged parity-violating (PV) contributions to the nuclear magnetic resonance shielding constant at the vibrational ground and excited states of the CHFClBr molecule. Electronic structure calculations were made using the B3LYP functional and an all-electron aug-cc-pVTZ basis set. For vibrational calculations, we used vibrational self consistent field and configuration integration (VCI[n]) methods with harmonic oscillator (HO) eigenfunctions as a basis for modal excitations. Potential energy surfaces were generated automatically using grid-based methods. We found out that a reasonably sized three-mode coupled grid combined with (state-specific) the VCI[4] level of theory and the 6 HO basis set results in converged ground state energies and fundamental vibrational frequenices. In addition to pure vibrational frequencies, we evaluated the vibrational corrections to the PV contributions to the nuclear magnetic shielding constant for all nuclei in the molecule. Property surfaces for PV linear response functions were generated using the two-mode coupled fourth-order Taylor expansion. Vibrational corrections were calculated for the ground state (GS) and also for vibrationally excited fundamental states. The convergence of the PV contributions was roughly similar to the convergence of energies. In general, vibrational corrections to the PV contributions were smaller in magnitude than the PV contribution calculated at the equilibrium geometry. The ground state vibrational averaging increased the 1H, 13C, and 79Br contributions by as much as 20%, and decreased the 19F and 35Cl contributions by less than 10%. Vibrational corrections from the selected vibrationally excited states were mostly of the same sign as the GS correction and sometimes substantially larger in magnitude. For example, vibrational averaging of the mode m4 increases the vibrational correction of the 79Br PV contribution to be 41%, which is approximately five times larger than the GS correction. For the 1H nucleus, averaging over the m1 excited state resulted in contributions of the same magnitude as the equilibrium geometry result. Corrections from the lowest energy modes (m9–m7) are within a factor of two compared the GS correction (with an exception of 1H nucleus). In conclusion, vibrational averaging of the ground state changed the magnitude of the PV contributions by a 1–10% in relevant nuclei. Five-fold increases in magnitudes of the corrections could be seen in vibrationally excited states, but more modest differences were found in low-lying states. The results indicate that the vibrational corrections to the PV contributions are of order a few to 40% for lower lying states of fairly rigid molecules. A correction of as much as 79% was observed for the higher lying CH stretching
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171
mode. We therefore conclude that vibrational contributions to PV effects are non-negligible even at the 0 K limit where zero-point vibrational contributions should be taken into account even for rigid molecules such as CHFClBr. Acknowledgements One of the authors (V.W.) acknowledges financial support from the Jenny and Antti Wihuri Fund. O.C. acknowledges EUROHORCs for a EURYI award and support from the Danish Center for Scientific Computing (DCSC), the Danish national research foundation, and the Lundbeck Foundation. Computational resources were provided by CSC-IT Center for Science Ltd. (Espoo, Finland). References [1] [2] [3] [4] [5] [6] [7] [8]
A.L. Barra, J.B. Robert, L. Wiesenfeld, Phys. Lett. A 115 (1986) 443. A.L. Barra, J.B. Robert, L. Wiesenfeld, Europhys. Lett. 5 (1988) 217. A.L. Barra, J.B. Robert, Mol. Phys. 88 (1996) 875. G. Laubender, R. Berger, ChemPhysChem 4 (2003) 395. A. Soncini, F. Faglioni, P. Lazzeretti, Phys. Rev. A 68 (2003) 033402. R. Bast, P. Schwerdtfeger, T. Saue, J. Chem. Phys. 125 (2006) 064504. G. Laubender, R. Berger, Phys. Rev. A 74 (2006) 032105. V. Weijo, R. Bast, P. Manninen, T. Saue, J. Vaara, J. Chem. Phys. 126 (2007) 074107.
[9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]
171
V. Weijo, P. Manninen, J. Vaara, J. Chem. Phys. 123 (2005) 054501. V. Weijo, P. Manninen, Chem. Phys. Lett. 433 (2006) 37. J.M. Bowman, S. Carter, X.C. Huang, Int. Rev. Phys. Chem. 22 (2003) 533. R. Gerber, J. Jung, in: P. Jensen, P.R. Bunker (Eds.), Computational Molecular Spectroscopy, Wiley, Chichester, 2000, p. 365. O. Christiansen, Phys. Chem. Chem. Phys. 9 (2007) 2942. R.G. Viglione, R. Zanasi, P. Lazzeretti, A. Ligabue, Phys. Rev. A 62 (2000) 052516. J.K. Laerdahl, P. Schwerdtfeger, H.M. Quiney, Phys. Rev. Lett. 84 (2000) 3811. M. Quack, J. Stohner, Phys. Rev. Lett. 84 (2000) 3807. P. Schwerdtfeger, J.K. Laerdahl, C. Chardonnet, Phys. Rev. A 65 (2002) 042508. P. Schwerdtfeger, T. Saue, J.N.P. van Stralen, L. Visscher, Phys. Rev. A 71 (2005) 012103. C. Daussy, T. Marrel, A. Amy-Klein, C.T. Nguyen, C.J. Bordé, C. Chardonnet, Phys. Rev. Lett. 83 (1999) 1554. M. Ziskind, C. Daussy, T. Marrel, C. Chardonnet, Eur. Phys. J. D 20 (2002) 219. J. Crassous, C. Chardonnet, T. Saue, P. Schwerdtfeger, Org. Biomol. Chem. 3 (2005) 2218. MIDASCPP (Molecular interactions, dynamics and simulation chemistry program package in C++), 2007.. DALTON, A molecular electronic structure program, release 2.0, 2005. . D. Toffoli, J. Kongsted, O. Christiansen, J. Chem. Phys. 127 (2007) 204106. K. Yagi, S. Hirata, K. Hirao, Theor. Chem. Acc. 118 (2007) 681. D.M. Bishop, B. Kirtman, J. Chem. Phys. 95 (1991) 2646. J.K.G. Watson, Mol. Phys. 15 (1968) 479. S. Carter, S.J. Culik, J.M. Bowman, J. Chem. Phys. 107 (1997) 10458. J. Kongsted, O. Christiansen, J. Chem. Phys. 125 (2006) 124108. G. Rauhut, V. Barone, P. Schwerdtfeger, J. Chem. Phys. 125 (2006) 054308.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Vibrational effects in the parity-violating contributions to the isotropic nuclear magnetic resonance chemical shift Ville Weijo a,b,*, Mikkel Bo Hansen c, Ove Christiansen c, Pekka Manninen d a
Laboratory of Physics, Helsinki University of Technology, P.O. Box 1100 (Otakaari 1 M), 02015 Espoo, Finland Laboratory of Physical Chemistry, Department of Chemistry, P.O. Box 55 (A.I. Virtasen Aukio 1), University of Helsinki, FI-00014 Helsinki, Finland c The Lundbeck Foundation Center for Theoretical Chemistry and Center for Oxygen Microscopy and Imaging, Department of Chemistry, University of Aarhus, Langelandsgade 140, DK-8000 Århus C, Denmark d CSC-IT Center for Science Ltd., P.O. Box 405, FI-02101 Espoo, Finland b
a r t i c l e
i n f o
Article history: Received 5 May 2008 In final form 13 January 2009 Available online 19 January 2009
a b s t r a c t We investigate the effect of vibrational corrections to the parity-violating (PV) contributions to the nuclear magnetic resonance shielding constant for the CHFClBr molecule using a non-relativistic framework. Density-functional theory is used for the electronic structure part and vibrational configuration interaction for the vibrational part. Zero-point vibrational corrections are found to be of the order of a less than 10% with respect to the PV contributions calculated at the equilibrium geometry in relevant cases. Vibrational corrections are evaluated also for the fundamental vibrations and are found to be up to 42% of the total PV contribution in the heavy elements. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction The first formulations of the parity-violating (PV) theory of nuclear magnetic resonance (NMR) spectral parameters were presented by Barra et al. [1–3]. The first ab initio calculations of PV contributions to the isotropic NMR chemical shift were made for small molecules, such as H2O2 [4,5], more than a decade later. Advanced computational aspects, like electron correlation effects, basis sets, and special relativity have been explored in recent years [6–8] and all of those aspects have been found to have significant effect on the calculated values. In addition to NMR shielding constants, PV contributions to the NMR spin–spin coupling constants have been evaluated [9], also electron spin resonance spectroscopy has been considered as a possible method for observing PV effects in molecules [10]. As for now, all predicted effects for the spectroscopical methods mentioned above have been well below the experimental limits. The effects of molecular vibrations have, on the other hand, so far not been considered in magnetic resonance properties. Tools for molecular vibrational calculations have been developed in particular in recent years (see, e.g., Refs. [11–13] and references therein). Vibrational methods usually treat electronic and vibrational wave functions separately in the spirit of the adiabatic (Born–Oppenheimer) approximation. Computations then proceed by using electronic structure methods to generate points on the potential energy surface (PES) and the property surfaces. These * Corresponding author. Present address: Department of Chemistry, University of Tromsø, N-9037 Tromsø, Norway. E-mail address: [email protected] (V. Weijo). 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.01.022
points can then be used for constructing an intermediate representation of the PES and property surfaces as function of the relative nuclear positions. Subsequently, a wave function for the nuclear motion is obtained from the PES by solving approximately the Schrödinger equation for the motion of the nuclei. Then using this wave function and the generated property surfaces, the vibrational averaging can be carried out. Different variants of this strategy vary in their selections of PES generation options, vibrational wave function method, and vibrational basis sets. We note that the most widely used approach for calculating anharmonic vibrational averages is probably still the simplest vibrational perturbation approach, circumventing the explicit construction of vibrational wave functions. Explicitly constructing vibrational wavefunctions may be considered to be more robust in many ways but often also comes with an increased cost from typically requiring PES information over a wide range of structures. Moreover, the calculations of the vibrational wave functions can become computationally demanding. In this Letter, we will use explicit vibrational wave function methods to calculate vibrational corrections to the PV isotropic NMR shielding constant contributions of a CHFClBr molecule. The molecule has been used previously in calculations of PV shielding and nuclear spin–spin coupling contributions [8]. Also, the molecule has been considered for a candidate for observation of PV effects in vibrational spectroscopy (see, e.g., Refs. [14–18] for theory and Refs. [19–21] for experiments). We focus here on CHFClBr rather than H2O2, which has been used more commonly in PV NMR calculations [4–8]. Although CHFClBr is not experimentally well-suited for detection of PV effects in NMR due to the high-spin nuclei in it, CHFClBr is a better analogue to a hypothetical
167
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171
heavy-element-containing experimentally relevant molecule than H2O2. The latter is also significantly more difficult to address using automated computational tools due to the presence of a torsional HOOH mode giving rise to a double minimum potential. Vibrational corrections to the property in question vary significantly depending on the nucleus in question and the amount of vibrational excitation. Due to the relative size of the corrections, they cannot be neglected in high-accuracy calculations.
2. Computational details The potential and property surfaces used in this study has been calculated using density-functional theory with the Becke threeparameter exchange-correlation hybrid functional (B3LYP) and Dunning’s augmented correlation-consistent triple-zeta basis set (aug-cc-pVTZ). Vibrational calculations were made using MIDASCPP vibrational structure program [22] calling the DALTON electronic structure program [23] for energies and properties at given points. Calculated surfaces and structures are available from the authors on request. The nuclei used in vibrational and response calculations were 1 H, 13C, 19F, 35Cl and 79Br. The molecule is in its (R)-enantiomeric form. Geometries and potential/property surfaces are available from the authors on request. 2.1. Vibrational structure calculations The representation of the potential energy surfaces (PES) and the property surfaces is a central point in vibrational structure theory and still an issue of much contemporary research. We shall give a brief discussion of the strategy we use in this work, referring to a number of reviews and recent research papers for more detailed accounts [11–13,24,25]. The surfaces can be expanded in a series of n-mode potential terms,
V
ð1Þ
;V
ð2Þ
;V
ð3Þ
ð1Þ
;...;
(1)
K
nates, rather than the full Watson operator and hence we neglect vibration–rotation interactions in this study.
(2)
where V includes only one-mode terms, while V includes all terms with up to two-mode couplings and so on. If one does not truncate the series, one obtains the fully coupled potential. The potential energy surfaces are generated as a multi-dimensional high-degree polynomial fit to a set of grid points [24]. The outline of the procedure is as follow. First the molecular structure is optimized and a vibrational analysis carried out. The vibrational analysis provides harmonic frequencies and mass-weighted normal coordinates which are subsequently used in the MIDASCPP program to determine the boundaries of the grid surface as well as the positions of the coarse grid points (CGP). The grid cut-offs are determined from the classical turning points for each vibrational mode (in mass-weighted normal coordinates),
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2h 1 xtp ¼ vþ ; 2 x
Prior to the fitting a number of additional points may be generated by interpolation via a direct product natural N-cubic spline. We used a fine mesh of 12 per coarse grid point keeping only the inner most 90% to avoid border effects. The resulting points were fitted to a polynomial of degree 12. More commonly used method based on numerical differentiation and Taylor expansions of the potential and molecular property surfaces has been also used. The step-sizes used for the individual modes in the numerical differentiation procedure are scaled by the corresponding harmonic frequencies multiplied with the same constant factor. In this study a factor of 0.070 has been used. Surfaces generated this way are denoted mMtT, where m denote the maximum mode-coupling level in the energy/property potential and t the degree of the Taylor expansion. For example, 2M4T means that all terms up to fourth-order depending on one- and two-mode couplings are included. At this point we emphasize that PESs are obtained via the grid approach, whereas property surfaces are based on fourth-order Taylor expansions. The vibrational wave function is constructed in terms of Hartree products of excitations of distinguishable modes. The one-mode basis is expanded in a certain number of harmonic oscillator (HO) functions. The simplest such vibrational wave function consists of a single Hartree product, which is optimized using the vibrational self consistent field method (VSCF) to minimize the energy (see Refs. [11–13] and references therein). The VSCF method treats the interaction between the modes in an averaged fashion, and to alleviate this, the vibrational configuration interaction (VCI) method is used. In the VCI calculations we excite simultaneously up to two, three, or four modes (VCI[2], VCI[3], and VCI[4], respectively). In this work, we use a state-specific approach. Thus, all states have their own VSCF reference state, and the excitation space is generated relative to that VSCF state. The vibrational structure calculations employ a simple kinetic P @2 energy operator, i.e. 12 K @Q 2 in mass-weighted normal coordi-
ð2Þ
where v is a vibrational quantum number, in this study a value of v = 10 is used. The grid is then divided into a number of equal-sized domains in the one-mode, two-mode, etc. parts, thus determining the coarse grid points, which are subsequently calculated via an interface to an electronic structure program. The number of grid points needed in each dimension typically decreases with increasing mode-combination level, i.e., all one-dimensional cuts are determined from a large number of coarse grid points, all twodimensional cuts from a smaller number of points in each dimension, and so forth. A 641161 grid is based on 64 CGPs in each 1D cut and 162 CGPs in each 2D cut.
2.2. Property evaluation The leading-order parity-violating contribution to the isotropic NMR shielding tensor in the non-relativistic framework can be expressed as a linear response function (see, e.g., Ref. [9])
rPV K;s ¼
1 DD
cK
PV
OZ
hK; ; hB0 ;s
EE 0
;
ð3Þ
which combines the orbital Zeeman operator OZ
hB0 ;s ¼
1X ‘iO;s ; 2 i
ð4Þ
and the nuclear spin-dependent part of the total parity-violating operator
X GF a PV 2 hK; ¼ pffiffiffi kK ð1 4 sin hW Þ iri; ; dðriK Þ þ : 2 2 i
ð5Þ
Here, K refers to the nucleus in question and and s are Cartesian indices of the tensor. ‘iO is the angular momentum operator of electron i with respect to the gauge origin O, and riK is the distance vector between i and K. The SI-based atomic unit system is used throughout this work. The weak coupling constant is denoted GF = 2.22 1014 a.u. and the Weinberg angle is hW = 0.2320. kK is a nucleus-dependent factor that includes contributions from the nuclear structure and the nuclear anapole moment. The kK factors are not usually known experimentally or theoretically and are set to unity during the calculations. If accurate values will become
168
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171
available, results can be scaled accordingly. This issue also does not change the relative importance of vibrational and electronic contributions. There is also another leading-order contribution, but it does not contribute to the isotropic part of the tensor [9]. It can be shown that the combined leading-order term is gauge-independent at the complete basis set limit [5], and thus, although there exists a reference to the gauge origin in Eq. (4), the isotropic part should be free of problems in this respect. Vibrational contributions to molecular properties fall into two categories: (i) the motion of the nuclei may respond to an external perturbation, thus giving rise to pure vibrational terms as for example (hyper)-polarizabilities [13,26], and/or (ii) the property of interest may depend on the nuclear geometry, and hence one have to average the property over at least the vibrational ground state, giving rise to so called zero-point vibrational averages (ZPVA). In this Letter, we are dealing solely with issue (ii). Contribution (i) is zero in our concrete case, since the electronic expectation value of the orbital Zeeman operator 4 is zero for the closedshell molecule considered here. The ZPVA is defined as
hAiZPVA ¼ hUjAjUi;
ð6Þ
where jUi may be either a VSCF or a VCI state and A may be any property that is a function of nuclear coordinates and expanded as described above [24,27–29]. In this work, we concentrate on the vibrational correction, which is the difference between the vibrational average and the equilibrium geometry results. 3. Results and discussion Harmonic zero-point energy (ZPE) and vibrational frequencies of CHFClBr are listed in Table 1. Our values for harmonic frequencies of CHFClBr are mostly in line with the results of Rauhut et al. [30]. They made their calculations with the same aug-cc-pVTZ basis set and B3LYP functional, but they used the 12C isotope instead of 13C. After scaling the harmonic frequencies of stretching modes by a square root of their reduced mass ratios, our results agree with Rauhuts within 3 cm1, with the exception of the C–Br stretch. The difference can be explained by their use of a relativistic effective core potential (RECP) for Br. Although RECPs are often very useful for heavy elements like Br to take care of relativistic effects, the way they were implemented in the currently used program package made them unsuitable in this case, because PV contributions are very sensitive to the description of the electronic structure in the vicinity of the nucleus.
Table 1 Zero-point energies (ZPE) and fundamental excitation frequencies (in cm1) of 13 1 19 35 79 C H F Cl Br evaluated using the harmonic approximation and also results using different PESs are included (VCI[2] calculation with the 6 HO vibrational basis set). The assignment of the modes follow the notation of Ref. [30]. PES
Harm.a
2M4T
32181
64181
641161
64116141
Assignment
ZPE
4547.5 3143.3 1318.6 1211.2 1049.3 720.9 622.7 415.6 306.3 218.9
4453.3 3005.3 1285.9 1180.8 1023.0 708.0 611.1 411.4 303.8 215.4
4443.6 3006.6 1275.1 1169.5 1022.5 707.2 611.2 411.3 303.5 215.4
4444.7 3006.6 1275.1 1169.6 1022.5 707.2 611.1 411.3 303.5 215.4
4453.4 3012.2 1284.1 1181.5 1022.2 707.2 610.6 411.3 303.5 215.4
4453.0 3013.0 1283.9 1179.7 1022.0 707.3 610.0 411.1 303.2 215.4
CH str. CHF sciss. CHBr/Cl sciss. CF str. CCl str. CBr str. CFCl sciss. CFBr sciss. CClBr sciss.
m1 m2 m3 m4 m5 m6 m7 m8 m9 a
Results using the harmonic approximation.
3.1. Anharmonic effects to vibrational frequencies Effects of anharmonicity calculated using different PESs are also shown in Table 1. It has previously been demonstrated that the use of low order Taylor expanded PESs in vibrational calculations may provide results of low accuracy. We have therefore chosen to also use grid based PESs in the calculations of anharmonic energies and wavefunctions. The 64181 grid is saturated with respect to one-mode couplings. The differences between 64181 and 641161 grids can be approximately 10 cm1 for the high-energy modes, but again for the low-energy modes, the differences are less than 1 cm1. 641161 is most likely close to convergence with respect to the two-mode couplings, whereas using just 8 points in either dimension of the two-mode cuts is too crude an approximation. The distance between the points simply becomes too large. The three-mode couplings do not change the situation much, the difference between 641161 and 64116141 grids is 2 cm1 at most (in mode m3) and often much less. This suggest that the PES series of Eq. (1) is very rapidly converging, being close to convergence already at the V(2) level. The fast convergence with respect to mode-coupling level was illustrated with state-specific VCI[2] wave functions, which could be considered not flexible enough. However, the convergence of ZPE and fundamental frequencies with respect to the vibrational correlation treatment is also rather rapid, as can be seen from Table 2 reporting calculations using the 64116141 potential. Even at the VSCF level, the low-energy modes do not differ from VCI[4] results by more than a couple of cm1. With ZPE and modes m2 and m3, the similar difference is at most 10 cm1. The only exception to the very good convergence is the mode m1 where a difference of 50 cm1 is observed. VCI[3] results are within 2 cm1 from the VCI[2] results with an exception of mode m6, where a slightly erratic behavior is seen. There, the VCI[2] treatment first reduces the frequency by 6 cm1 compared to VSCF, and subsequently, VCI[3] restores the value close to the VSCF result again. The VCI[4] calculations lead to additional changes to VCI[3] of an order of a few tenths of a cm1 with a slightly larger effect for m1. All states are dominated by a single configuration meaning they are physically well described as a fundamental. The highest lying fundamental, m1, obtains somewhat larger contributions to its wave function from other components. A one-mode basis with six HO basis functions is also very close to its respective limit, and the difference between six HO and eight HO basis sets is at most 1 cm1 (in mode m1) and is often much smaller (values not shown here). Again, any experimental reference values for fundamental frequencies are not available for 13C1H19F35Cl79Br. However, if we compared the general features of a vibrationally correlated gridbased calculation to its harmonic counterpart, we would observe a similar trend of lowering of frequencies as can be found for
Table 2 Effects of vibrational correlation to the ZPE and fundamental frequencies of 13 1 19 35 79 C H F Cl Br (in cm1). The PES is from the 64116141 grid and the 6 HO basis set is used. State
VSCF
VCI[2]
VCI[3]
VCI[4]
ZPE
4460.5 2965.5 1289.3 1188.7 1024.2 709.0 615.9 412.0 303.8 215.5
4453.0 3013.0 1283.9 1179.7 1022.0 707.3 610.0 411.1 303.2 215.4
4452.4 3012.9 1282.0 1178.7 1021.5 705.1 617.2 410.8 302.8 215.2
4452.4 3011.7 1281.6 1178.5 1021.4 704.9 615.0 410.7 302.8 215.2
m1 m2 m3 m4 m5 m6 m7 m8 m9
169
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171 12
CHFClBr in vibrational perturbation theory calculations [30]. There, the computational errors are at most 50 cm1 compared to the experimental values, and we expect that our results achieve comparable accuracy, even though RECPs are not used. 3.2. Vibrationally averaged parity-violating contributions PV contributions to the NMR shielding at the equilibrium geometry for all the nuclei in CHFClBr are presented in the column ‘Eq.’ of Table 3. Similar PV contributions were also evaluated in Ref. [9] for the same system with the same functional but different basis set. There are a couple of differences in the results. Excluding the 79 Br contributions, signs of the contributions are reversed, which can be explained by different signs of the k-factors used in Ref. [9]. The contributions calculated in this study (excluding 79Br) are somewhat smaller in magnitude than their counterparts calculated in Ref. [9]. This behavior can be linked to basis set effects. It is known that general purpose basis sets lack the proper flexibility for core electrons (see, e.g., Refs. [7,8]) in this property. Augmenting aug-cc-pVTZ basis set with steep s and p-type Gaussians can result to large differences in PV contributions [7]. The drawback is the extra computational effort, which can be very large, especially if many computations are required as is in our case. Core orbitals, however, are only slightly affected by molecular environment and their effect to the relative magnitudes between vibrationally uncorrected and corrected PV contributions should be small. Basis set effects are also responsible for the sign difference in 79Br contributions (after including the difference in k factors). Table 3 also includes vibrational corrections to the equilibrium values using different levels of vibrational correlation methods and different basis sets. The PES is chosen to be the 64116141 grid, which was shown above to be sufficiently close to convergence in the vibrational frequency calculations. The property surface is based on the 2M4T Taylor expansion, because this level offers considerable computational savings over the grid based methods in time-consuming response calculations. The convergence of the ground state vibrationally corrections is very robust, as can be expected from the convergence of ground state energies in Table 2. Differences in various levels of VCI and HO basis sets are very small overall. The 6 HO basis set can be considered to be fully converged as the difference between the 4 HO and 6 HO results is of the order of per miles at most. Similar differences can be found between the VCI[3] and VCI[4] results, and thus the VCI[4] calculations with the 6 HO basis set is in practice converged. At this point we note that previous studies, see Ref. [17] and references therein, estimate the PV effects arising from molecular vibrational motion from uncoupled potential/property surfaces. In order to test the quality of this approximation, we tested several levels of both the PES and property surfaces. It was found that the PV effects are highly dependent on the quality of the PES, which in turn determines the accuracy of the vibrational wavefunction. For example, using an uncoupled anharmonic PES with any property
surface give vibrational corrections to the PV contribution which are very far from the VCI[4] results using a proper PES. On the other hand, with a proper vibrational wavefunction, e.g. VCI[3], obtained from a good PES the effect of mode–mode coupling in the property surfaces is more limited. The absolute deviations of VCI[3] using uncoupled property surfaces and a good PES from the most accurate results presented in Table 3 amounts to about 12%, 46%, 11%, 1%, and 1% for the 13C, 1H, 19F, 35Cl, and 79Br nuclei, respectively. The deviation for the 1H nucleus is rather large. This is, however, believed to be of minor importance, as the PV effect is unlikely to be measured for this nucleus due to its small value. From the above study it is therefore problematic to give a guideline for including mode–mode-coupling in the property surfaces. One should, however, make sure that the underlying vibrational wavefunction is satisfactory, which in turn is reflected by the quality of the PES entering the Hamiltonian. Vibrational corrections to the vibrationally excited states behave similarly with respect to basis sets and vibrational correlation. Vibrational corrections for all modes are shown in Table 4. The results were obtained using the VCI[4] method and 6 HO basis set. As compared to the equilibrium geometry reference values, vibrational corrections to the ground state (GS) are smaller. The largest correction is for 1H, where the correction is approximately 20%, while the smallest vibrational corrections are for 13C (of the order 1%). The vibrational corrections to 19F, 35Cl, and 79Br shielding contributions are 2.4%, 6.6%, and 8.4% of the equilibrium geometry values, respectively. Hence, the vibrational corrections do not systematically depend on the mass of the nucleus. They are, however found to be non-negligible. The vibrational corrections are generally significant for the considered vibrationally excited states, although some variations are observed. Most of the 13C corrections are equal or larger in magnitude than the GS correction, and depending on the mode, increase or decrease the magnitude of the total contribution. Correction from the mode m3 has the largest value of 1.23 1011 ppm, which is approximately 13% correction to the equilibrium value. Most of the 1H corrections increase the magnitude of the total contribution. The largest correction is close to be as large (mode m8) as the equilibrium geometry value, which is clearly substantial. In the case of 19 F corrections, the magnitude of the total PV contribution is reduced by less than 10% at most (mode m1). Also, most of the corrections are of a different sign than the equilibrium geometry result, thereby reducing the magnitude of the total PV contribution. Likewise in 35Cl corrections the magnitude of the total contribution is often reduced by the vibrational averaging of excited states. The magnitudes of the corrections are within a factor of four from the GS correction with a maximum at the mode m3, where the correction decreases the magnitude of the equilibrium value by 21%. Vibrational corrections in 79Br are mostly negative (increasing the magnitude) although four out of nine (excited states) are positive. Magnitudes are, again, roughly of the size of the GS correction, but the largest of the corrections (mode m4) is rather substantial and provides a 41% increase in magnitude to the equilibrium geometry values.
Table 3 Vibrationally averaged ground state corrections to the parityviolating contributions to the isotropic NMR shieldings in CHFClBr using different vibrational correlation methods and basis sets. Equilibrium geometry reference values (Eq.) are also given. All values in ppm. Nucleus
13
C H 19 F 35 Cl 79 Br 1
VCI[2]
VCI[3]
VCI[4]
Eq.
4 HO
6 HO
4 HO
6 HO
4 HO
6 HO
9.45 1013 3.06 1014 2.81 1012 1.76 1010 1.43 1010
9.25 1013 3.07 1014 2.86 1012 1.78 1010 1.45 1010
9.79 1013 3.11 1014 2.78 1012 1.76 1010 1.45 1010
9.60 1013 3.12 1014 2.83 1012 1.78 1010 1.47 1010
9.69 1013 3.11 1014 2.82 1012 1.78 1010 1.47 1010
9.68 1013 3.11 1014 2.82 1012 1.78 1010 1.47 1010
9.48 1011 1.51 1013 1.20 1010 2.68 109 1.76 109
170
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171
Table 4 Vibrational corrections to the parity-violating contributions to the NMR shielding constants of nuclei in CHFClBr at the ground state (GS) and vibrationally excited states. Equilibrium geometry reference values (Eq.) are also given. Vibrational calculations were made at the VCI[4] level with the 6 HO basis set. All values in ppm. Mode Eq. GS
m1 m2 m3 m4 m5 m6 m7 m8 m9
13
1
C
19
H
11
9.48 10 9.68 1013 5.26 1012 1.48 1012 1.23 1011 7.52 1012 7.18 1012 1.56 1012 7.83 1013 3.69 1012 3.92 1012
35
F
13
1.51 10 3.11 1014 2.74 1014 8.14 1014 1.65 1014 5.31 1014 1.01 1016 6.30 1014 5.35 1014 1.19 1013 8.84 1015
Connection between the magnitude of a vibrational correction and motions in a particular mode is not a clear one. For example, the largest vibrational correction to the 79Br PV contribution can be found from the mode m4, which consists mostly of stretching motion between C and F. Based on this, one could conclude that large variations in total molecular geometry would result in large vibrational corrections. Then, one could expect CH stretching mode (m1) to have even larger corrections, but in this case, the vibrational correction is two orders of magnitude smaller than that of the mode m4. Similar lack of systematic correlation between mode types and vibrational corrections can be found from all the other PV contributions, too. If one considers the case of a room temperature experiment, then the most relevant modes are the GS, m9, m8, the first overtone of m9, and perhaps also the mode m7. Here we exclude 1H corrections from the discussion, because the total 1H PV contribution is orders of magnitude smaller than that of the other nuclei. Vibrational corrections from the modes mentioned above (except the overtone) range from 10% to +16% of the equilibrium values. For order of magnitude estimates of PV effects in NMR spectrum these corrections are not too significant, but an experimental setup will suffer from line width broadening and associated loss of peak intensity. Also, if an accurate experiment is realized and, e.g., k factors are to be evaluated using quantum chemical calculations, variations in PV contributions of the low-lying modes impose a serious computational challenge. Experimentally relevant molecules can be larger than CHFClBr and most likely will contain even heavier elements. Accurate electronic structure calculations themselves on such systems are tricky and time consuming, but if vibrational structure calculations are included into the process, then the computational effort will truly be massive for high-accuracy. In general, the CHFClBr molecule is quite rigid. Somewhat larger relative changes in PV contributions due to the vibrational corrections could be expected for floppy molecules such as H2O2, where the PV contribution varies significantly near the equilibrium geometry as a function of the dihedral angle [4–8]. We believe that our numbers with the B3LYP functional and the aug-cc-pVTZ basis set for the electronic structure part and VCI[4] with 6 HO basis set for vibrational part provide good estimates of the relative magnitudes of vibrational corrections. The non-relativistic Hamiltonian used in our calculations, on the other hand, leaves room for errors. Relativistic effects have been shown to increase the PV contributions by orders of magnitudes for heavy elements [6,8], and in some cases, the spin–orbit interaction can even change the sign of the contribution. At the equilibrium geometry, however, the effect is not nearly as pronounced, and the relativistic effects only increase the magnitude of the PV contribution by a factor of 1.3 in third row elements [8]. Also, we expect that the property surface near the equilibrium geometry in the relativistic case is qualitatively similar to the non-relativistic case, at least in rigid molecules like
10
1.20 10 2.82 1012 1.10 1011 7.14 1012 6.30 1012 2.70 1012 6.70 1012 4.24 1012 6.55 1012 1.77 1012 2.58 1012
79
Cl
Br
09
2.68 10 1.78 1010 2.14 1010 1.82 1010 5.74 1010 5.49 1010 3.66 1011 1.85 1010 2.83 1010 3.10 1010 2.67 1010
1.76 1009 1.47 1010 4.34 1012 5.04 1011 5.64 1010 7.35 1010 2.81 1011 6.91 1011 1.82 1010 8.02 1011 2.76 1010
CHFClBr. Thus, we expect our findings to be fairly generalizable to other rigid molecules containing even heavier elements. 4. Conclusions We have computed vibrationally averaged parity-violating (PV) contributions to the nuclear magnetic resonance shielding constant at the vibrational ground and excited states of the CHFClBr molecule. Electronic structure calculations were made using the B3LYP functional and an all-electron aug-cc-pVTZ basis set. For vibrational calculations, we used vibrational self consistent field and configuration integration (VCI[n]) methods with harmonic oscillator (HO) eigenfunctions as a basis for modal excitations. Potential energy surfaces were generated automatically using grid-based methods. We found out that a reasonably sized three-mode coupled grid combined with (state-specific) the VCI[4] level of theory and the 6 HO basis set results in converged ground state energies and fundamental vibrational frequenices. In addition to pure vibrational frequencies, we evaluated the vibrational corrections to the PV contributions to the nuclear magnetic shielding constant for all nuclei in the molecule. Property surfaces for PV linear response functions were generated using the two-mode coupled fourth-order Taylor expansion. Vibrational corrections were calculated for the ground state (GS) and also for vibrationally excited fundamental states. The convergence of the PV contributions was roughly similar to the convergence of energies. In general, vibrational corrections to the PV contributions were smaller in magnitude than the PV contribution calculated at the equilibrium geometry. The ground state vibrational averaging increased the 1H, 13C, and 79Br contributions by as much as 20%, and decreased the 19F and 35Cl contributions by less than 10%. Vibrational corrections from the selected vibrationally excited states were mostly of the same sign as the GS correction and sometimes substantially larger in magnitude. For example, vibrational averaging of the mode m4 increases the vibrational correction of the 79Br PV contribution to be 41%, which is approximately five times larger than the GS correction. For the 1H nucleus, averaging over the m1 excited state resulted in contributions of the same magnitude as the equilibrium geometry result. Corrections from the lowest energy modes (m9–m7) are within a factor of two compared the GS correction (with an exception of 1H nucleus). In conclusion, vibrational averaging of the ground state changed the magnitude of the PV contributions by a 1–10% in relevant nuclei. Five-fold increases in magnitudes of the corrections could be seen in vibrationally excited states, but more modest differences were found in low-lying states. The results indicate that the vibrational corrections to the PV contributions are of order a few to 40% for lower lying states of fairly rigid molecules. A correction of as much as 79% was observed for the higher lying CH stretching
V. Weijo et al. / Chemical Physics Letters 470 (2009) 166–171
mode. We therefore conclude that vibrational contributions to PV effects are non-negligible even at the 0 K limit where zero-point vibrational contributions should be taken into account even for rigid molecules such as CHFClBr. Acknowledgements One of the authors (V.W.) acknowledges financial support from the Jenny and Antti Wihuri Fund. O.C. acknowledges EUROHORCs for a EURYI award and support from the Danish Center for Scientific Computing (DCSC), the Danish national research foundation, and the Lundbeck Foundation. Computational resources were provided by CSC-IT Center for Science Ltd. (Espoo, Finland). References [1] [2] [3] [4] [5] [6] [7] [8]
A.L. Barra, J.B. Robert, L. Wiesenfeld, Phys. Lett. A 115 (1986) 443. A.L. Barra, J.B. Robert, L. Wiesenfeld, Europhys. Lett. 5 (1988) 217. A.L. Barra, J.B. Robert, Mol. Phys. 88 (1996) 875. G. Laubender, R. Berger, ChemPhysChem 4 (2003) 395. A. Soncini, F. Faglioni, P. Lazzeretti, Phys. Rev. A 68 (2003) 033402. R. Bast, P. Schwerdtfeger, T. Saue, J. Chem. Phys. 125 (2006) 064504. G. Laubender, R. Berger, Phys. Rev. A 74 (2006) 032105. V. Weijo, R. Bast, P. Manninen, T. Saue, J. Vaara, J. Chem. Phys. 126 (2007) 074107.
[9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]
171
V. Weijo, P. Manninen, J. Vaara, J. Chem. Phys. 123 (2005) 054501. V. Weijo, P. Manninen, Chem. Phys. Lett. 433 (2006) 37. J.M. Bowman, S. Carter, X.C. Huang, Int. Rev. Phys. Chem. 22 (2003) 533. R. Gerber, J. Jung, in: P. Jensen, P.R. Bunker (Eds.), Computational Molecular Spectroscopy, Wiley, Chichester, 2000, p. 365. O. Christiansen, Phys. Chem. Chem. Phys. 9 (2007) 2942. R.G. Viglione, R. Zanasi, P. Lazzeretti, A. Ligabue, Phys. Rev. A 62 (2000) 052516. J.K. Laerdahl, P. Schwerdtfeger, H.M. Quiney, Phys. Rev. Lett. 84 (2000) 3811. M. Quack, J. Stohner, Phys. Rev. Lett. 84 (2000) 3807. P. Schwerdtfeger, J.K. Laerdahl, C. Chardonnet, Phys. Rev. A 65 (2002) 042508. P. Schwerdtfeger, T. Saue, J.N.P. van Stralen, L. Visscher, Phys. Rev. A 71 (2005) 012103. C. Daussy, T. Marrel, A. Amy-Klein, C.T. Nguyen, C.J. Bordé, C. Chardonnet, Phys. Rev. Lett. 83 (1999) 1554. M. Ziskind, C. Daussy, T. Marrel, C. Chardonnet, Eur. Phys. J. D 20 (2002) 219. J. Crassous, C. Chardonnet, T. Saue, P. Schwerdtfeger, Org. Biomol. Chem. 3 (2005) 2218. MIDASCPP (Molecular interactions, dynamics and simulation chemistry program package in C++), 2007.