Quasi-relativistic study of nuclear electric quadrupole coupling constants in chiral molecules containing heavy elements

Abstract

A quasi-relativistic (two-component) zeroth-order regular approximation model potential approach is employed to calculate nuclear electric quadrupole coupling tensors of the chiral polyhalomethane derivatives CHBrClF and CHClFI. Each of these two chiral compounds contains two nuclei with electric quadrupole moments, which poses special challenges for an analysis of high-resolution rovibrational spectra and, thus, renders predictions of electric field gradients and resulting nuclear electric quadrupole coupling tensors in such molecules particularly valuable. The theoretical approach is studied with respect to the influence of relativistic effects, by comparison to scalar-relativistic and non-relativistic calculations, and electron correlation effects, which are included within different flavours of density functional theory. The quality of the calculations is assessed by comparison to results from experiments reported previously by Bauder et al. and Soulard et al., respectively. It is shown that the local density approximation exchange-correlation functional gives results that agree with experiment within 12 % for all nuclear quadrupole coupling tensor components of ⁷⁹Br and ¹²⁷I. This benchmark test of the quasi-relativistic approach indicates its applicability to support upcoming microwave experiments with chiral heavy-element containing main group compounds.

GRAPHICAL ABSTRACT

Keywords:

• Electric field gradients
• ZORA
• relativistic effects
• nuclear quadrupole moments
• DFT

1. Introduction

Nuclear electric quadrupole moments (NEQMs) are a window into nuclear structure [1]. Presently, the most accurate method for determining NEQMs relies on the fruitful interplay between atomic or molecular precision spectroscopy and electronic structure theory [2]. In atomic and molecular spectroscopy NEQM coupling tensors are determined, which are a product of the NEQM and the electric field gradient (EFG) tensor probed at the position of the nucleus. The latter is not directly accessible by experiment, but can be calculated with electronic structure methods. Once the NEQM of a nucleus is determined, the measurement of NEQM coupling tensors in high-resolution rovibrational spectra of molecules provides important information in particular on the electronic density distribution in the vicinity of the nuclei. Due to this property, relativistic effects can play a major role in the calculation of EFGs of heavy element containing compounds [3,4]. As several atomic and molecular properties of current fundamental interest, such as a possible permanent electric dipole moment in molecules and atoms caused by parity- and time-reversal violating electric dipole moments of the nucleons, nuclear Schiff moments, nuclear magnetic quadrupole moments or an electric dipole moment of the electron, depend on the behaviour of the electronic wavefunction near the nucleus and benefit from relativistic enhancement effects [5–7], NEQM coupling tensors can serve as tests for the accuracy of theoretical approaches for the description of symmetry violating effects in molecular systems.

Whereas for an accurate prediction of NEQM coupling tensors in molecules with light nuclei a perturbative description of relativistic effects gives excellent results (see for instance, the direct perturbation theory studies of Stopkowicz and Gauss in Refs. [8,9]), in heavier compounds relativistic effects are preferentially included self-consistently. Many different relativistic approaches, ranging from four-component Dirac–Hatree–Fock [3,4,10] and coupled cluster methods [11–15] to two- and one-component Douglas–Kroll–Hess (DKH) and zeroth-order regular approximation (ZORA) approaches [16–21] have been considered. It was shown in some of these previous studies that ZORA performs well for the description of these relativistic effects in heavy elemental compounds [16,18,19,21] as the EFG pronouncedly depends on the behaviour of the valence electrons near the nuclei, whereas contribution from the core orbitals themselves remain usually below 15 % [18].

Jürgen Gauss and colleagues have performed combined experimental and theoretical high-level studies of NEQM coupling tensors in rotational spectra of the methane derivatives CH${}_2$BrF [22,23] and CH${}_2$FI [24,25], which feature each a single quadrupolar halogen nucleus. They were also able to extract NEQM coupling tensors of rare isotopes such as deuterium in the corresponding deuterated derivatives.

In high-resolution rotational spectroscopy of chiral polyhalogenated molecules, additional challenges have to be overcome as different halogen nuclei contribute to the spectral splittings in experiment [26,27]. Furthermore, the full NEQM coupling tensors needs to be determined. This renders quantum chemical calculations of NEQM couplings in chiral molecules particularly valuable.

In this paper, we demonstrate the applicability of model potential ZORA-DFT methods [7,28] for heavy-element containing chiral molecules that are of special interest for a measurement of molecular parity violation (see e.g. Ref [29] for a recent review on molecular parity violation), i.e. the polyhalomethanes CHBrClF [26,30–37] and CHClFI [33,37,38]. Both molecules were studied intensely by experiment giving a good reference to validate our methods for such chiral systems. Chiral compounds that contain even heavier nuclei will be in the focus of future research.

2. Methodology and computational details

The static interaction of a nuclear electric quadrupole moment (NEQM) described by the symmetric second rank tensor ${\mathit{Q}}_A$ of nucleus A with the electron cloud of a molecule with $N_{\mathrm{elec}}$ electrons and $N_{\mathrm{nuc}}$ nuclei is described in a four-component relativistic electronic structure framework by the Hamiltonian (see e.g. Ref. [39]; we follow here the sign convention of Refs. [26,27]): (1) $\begin{array}{rl}{\hat{H}}_{\mathrm{Q},A}& ={\hat{H}}_{\mathrm{Q},\mathrm{elec},A}+{\hat{H}}_{\mathrm{Q},\mathrm{nuc},A}\\ {\hat{H}}_{\mathrm{Q},\mathrm{elec},A}& =\frac{1}{6}\sum _i^{N_{\mathrm{elec}}}\mathrm{Tr}({\mathit{Q}}_A\cdot {\mathit{V}}_{\mathrm{iA}}^{(2)}){\mathbf{1}}_{4\times 4,i}\\ {\hat{H}}_{\mathrm{Q},\mathrm{nuc},A}& =\frac{1}{6}\sum _{B\ne A}^{N_{\mathrm{nuc}}}\mathrm{Tr}({\mathit{Q}}_A\cdot {\mathit{V}}_{\mathrm{BA}}^{(2)}),\end{array}$(1) with ${\mathit{V}}_{\mathrm{iA}}^{(2)}$ denoting the negative of the EFG at nucleus A from electron i, ${\mathit{V}}_{\mathrm{BA}}^{(2)}$ denoting the negative of the EFG at nucleus A from nucleus B and $N_{\mathrm{elec}}$ being the number of electrons in the system.

The quadrupole moment tensor can be written in terms of the nuclear electric quadrupole moment $e{Q}_A$ as ${\mathit{Q}}_A=e{Q}_A{\mathit{T}}_A$, where ${\mathit{T}}_A$ is a dimensionless symmetric second rank tensor depending on the nuclear spin. With this, we can also define an effective quadrupole coupling tensor ${\mathcal{X}}_A$ for the NEQM entering the effective spin Hamiltonian as (2) $\begin{array}{r}{H}_{\mathrm{Q},A}=\frac{1}{6}\mathrm{Tr}({\mathit{T}}_A\cdot {\mathcal{X}}_A),\end{array}$(2) with elements of the NEQM coupling tensor being (see e.g. Ref. [40]) (3) $\begin{array}{rl}{\mathcal{X}}_{\mathrm{\alpha \beta },A}& =e{Q}_A[⟨{\mathrm{\Psi }}_{\mathrm{e}}\left|\sum _i^{N_{\mathrm{elec}}}{\mathbf{1}}_{4\times 4,i}{V}_{\mathrm{\alpha \beta },\mathrm{iA}}^{(2)}\right|{\mathrm{\Psi }}_{\mathrm{e}}⟩\\ & +\sum _B^{N_{\mathrm{nuc}}}{V}_{\mathrm{\alpha \beta },\mathrm{BA}}^{(2)}].\end{array}$(3) Here ${\mathrm{\Psi }}_{\mathrm{e}}$ is the multi-electron wave function. The first of the terms in square brackets in Equation (3(3) $\begin{array}{rl}{\mathcal{X}}_{\mathrm{\alpha \beta },A}& =e{Q}_A[⟨{\mathrm{\Psi }}_{\mathrm{e}}\left|\sum _i^{N_{\mathrm{elec}}}{\mathbf{1}}_{4\times 4,i}{V}_{\mathrm{\alpha \beta },\mathrm{iA}}^{(2)}\right|{\mathrm{\Psi }}_{\mathrm{e}}⟩\\ & +\sum _B^{N_{\mathrm{nuc}}}{V}_{\mathrm{\alpha \beta },\mathrm{BA}}^{(2)}].\end{array}$(3) ) is due to the contribution from the electrons in the system, whereas the second-term accounts for the electric field gradient caused by all other nuclei. For point-like particles, the negative of the EFG is in SI units (4) $\begin{array}{rl}{V}_{\mathrm{\alpha \beta },\mathrm{ij}}^{(2)}& ={\frac{q_i}{4\pi {ϵ}_0}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{r}_j^{\alpha }\mathrm{\partial }{r}_j^{\beta }}\frac{1}{|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j|}|}_{{\overrightarrow{r}}_j={\overrightarrow{R}}_j\ne {\overrightarrow{r}}_i}\\ & =\frac{q_i}{4\pi {ϵ}_0}\frac{3(r_i^{\alpha }-R_j^{\alpha })(r_i^{\beta }-R_j^{\beta })-{\delta }_{\mathrm{\alpha \beta }}{|{\overrightarrow{r}}_i-{\overrightarrow{R}}_j|}^2}{{|{\overrightarrow{r}}_i-{\overrightarrow{R}}_j|}^5},\end{array}$(4) with $q_i=-e$ for electrons and $q_B=e{Z}_B$ for nuclei, with $Z_B$ being the nuclear charge number of nucleus B. As the point-nucleus potential has a pole at ${\overrightarrow{r}}_i={\overrightarrow{R}}_j$ there is an additional term to its derivative in the distributional sense proportional to $\delta ({\overrightarrow{R}}_j)$. This term does not contribute to the NEQC tensor and, thus, is not included in the present calculations. Here, ${ϵ}_0$ is the electric constant, $\overrightarrow{r}=(x,y,z{)}^{\mathrm{T}}$ is the position vector and ${\delta }_{\mathrm{\alpha \beta }}$ is the Kronecker delta, with $r^{\alpha }$ and $r^{\beta }$ denoting one of the Cartesian components x, y, z. Effects of the finite nuclear size on the form of the EFG operator are negligibly small for the molecules studied herein, but may become more pronounced for molecules that contain heavier elements as we have shown for the electric field at the nucleus in another context [6].

Due to appearance of the $4\times 4$ identity matrix in the electronic expectation value in Equation (3(3) $\begin{array}{rl}{\mathcal{X}}_{\mathrm{\alpha \beta },A}& =e{Q}_A[⟨{\mathrm{\Psi }}_{\mathrm{e}}\left|\sum _i^{N_{\mathrm{elec}}}{\mathbf{1}}_{4\times 4,i}{V}_{\mathrm{\alpha \beta },\mathrm{iA}}^{(2)}\right|{\mathrm{\Psi }}_{\mathrm{e}}⟩\\ & +\sum _B^{N_{\mathrm{nuc}}}{V}_{\mathrm{\alpha \beta },\mathrm{BA}}^{(2)}].\end{array}$(3) ), both the large and the small component densities of the relativistic electronic wavefunction contribute within a four-component formulation of the EFG operator to the electronic part of Equation (4(4) $\begin{array}{rl}{V}_{\mathrm{\alpha \beta },\mathrm{ij}}^{(2)}& ={\frac{q_i}{4\pi {ϵ}_0}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{r}_j^{\alpha }\mathrm{\partial }{r}_j^{\beta }}\frac{1}{|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j|}|}_{{\overrightarrow{r}}_j={\overrightarrow{R}}_j\ne {\overrightarrow{r}}_i}\\ & =\frac{q_i}{4\pi {ϵ}_0}\frac{3(r_i^{\alpha }-R_j^{\alpha })(r_i^{\beta }-R_j^{\beta })-{\delta }_{\mathrm{\alpha \beta }}{|{\overrightarrow{r}}_i-{\overrightarrow{R}}_j|}^2}{{|{\overrightarrow{r}}_i-{\overrightarrow{R}}_j|}^5},\end{array}$(4) ). When the multi-electron wave function is approximated as a single determinant composed from $N_{\mathrm{orb}}$ four-component molecular orbitals ${\psi }_{\mathrm{e},i}$ with large (or upper) components ${\psi }_{\mathrm{e},i}^{\mathrm{L}}$ and small (or lower) components ${\psi }_{\mathrm{e},i}^{\mathrm{S}}$, the contribution to the electronic part reads as (5) $\begin{array}{rl}{\mathcal{X}}_{\mathrm{\alpha \beta },A}^{\mathrm{elec}}& =\sum _{i=1}^{N_{\mathrm{orb}}}{n}_{i}e{Q}_A[⟨{\psi }_{\mathrm{e},i}^{\mathrm{L}}\left|{\mathbf{1}}_{2\times 2}{V}_{\mathrm{\alpha \beta },1A}^{(2)}\right|{\psi }_{\mathrm{e},i}^{\mathrm{L}}⟩\\ & +⟨{\psi }_{\mathrm{e},i}^{\mathrm{S}}\left|{\mathbf{1}}_{2\times 2}{V}_{\mathrm{\alpha \beta },1A}^{(2)}\right|{\psi }_{\mathrm{e},i}^{\mathrm{S}}⟩],\end{array}$(5) with $n_i$ being the occupation number of molecular orbital i and with 1 serving as a dummy number to indicate the single electron occupying the molecular orbital.

In the present paper, we follow our toolbox approach for the calculations of two-component properties within ZORA [7]. For this purpose, we compute the ZORA Hartree–Fock (HF) or Kohn–Sham (KS) molecular orbitals by solving the ZORA equations self-consistently: (6) $\begin{array}{rl}& [{\mathit{V}}_{\mathrm{cGHF}/\mathrm{cGKS}}+\overrightarrow{\mathit{\sigma }}\cdot \hat{\overrightarrow{p}}\frac{c^2}{2m_{\mathrm{e}}{c}^2-\tilde{V}(\overrightarrow{r})}\overrightarrow{\mathit{\sigma }}\cdot \hat{\overrightarrow{p}}]{\psi }_{\mathrm{e},i}^{\mathrm{ZORA}}\\ & ={ϵ}_i{\psi }_{\mathrm{e},i}^{\mathrm{ZORA}},\end{array}$(6) where $\tilde{V}$ is the model potential as proposed by van Wüllen [28], which is employed in order to alleviate the gauge dependence of ZORA. and ${\mathit{V}}_{\mathrm{cGHF}/\mathrm{cGKS}}$ is the effective potential within the complex generalised HF (cGHF) or KS (cGKS) approximation.

In this implementation, the large component is approximated as the ZORA wave function ${\psi }^{\mathrm{L}}\approx {\psi }^{\mathrm{ZORA}}$ and terms involving the small component can be calculated by automatic generation of an approximate small component in the sense of ZORA: (7) $\begin{array}{r}{\psi }_{\mathrm{e},i}^{\mathrm{S}}\approx \frac{c}{2m_{\mathrm{e}}{c}^2-\tilde{V}(\overrightarrow{r})}\overrightarrow{\mathit{\sigma }}\cdot \hat{\overrightarrow{p}}{\psi }_{\mathrm{e},i}^{\mathrm{ZORA}},\end{array}$(7) In the sense of our toolbox approach [7], the electronic EFG operator can be represented by the four component density function ${\mathrm{\Gamma }}^{0,0,0}$, a tensor operator $\hat{\mathsf{t}}={\mathit{V}}_{\mathrm{iA}}^{(2)}$, a scalar operator $\hat{s}=e{Q}_A$ and no differential operator, where we have used the notation detailed in Ref. [7]. Within this notation, contributions to ${\mathrm{\Gamma }}^{0,0,0}$ from the approximate ZORA two component number density functions ${\rho }^{\mathrm{LL},0}$ and ${\rho }^{\mathrm{SS},0}$ have to be evaluated. The contribution of the nuclear-nuclear repulsion was directly evaluated as written in Equation (3(3) $\begin{array}{rl}{\mathcal{X}}_{\mathrm{\alpha \beta },A}& =e{Q}_A[⟨{\mathrm{\Psi }}_{\mathrm{e}}\left|\sum _i^{N_{\mathrm{elec}}}{\mathbf{1}}_{4\times 4,i}{V}_{\mathrm{\alpha \beta },\mathrm{iA}}^{(2)}\right|{\mathrm{\Psi }}_{\mathrm{e}}⟩\\ & +\sum _B^{N_{\mathrm{nuc}}}{V}_{\mathrm{\alpha \beta },\mathrm{BA}}^{(2)}].\end{array}$(3) ) and added subsequently to the electronic contribution to receive the NEQC tensor. For more details on our toolbox approach and the corresponding notation, we refer the reader to Ref. [7].

In the following, we will refer to property calculations that include contributions from both ${\rho }^{\mathrm{LL},0}$ and ${\rho }^{\mathrm{SS},0}$ and use the ZORA orbitals in renormalised form via $\frac{{\psi }_{\mathrm{e},i}^{\mathrm{ZORA}}}{\sqrt{1+\int {\mathrm{d}}^{3}r({\psi }_{\mathrm{e},i}^{\mathrm{S}}{)}^{†}{\psi }_{\mathrm{e},i}^{\mathrm{S}}}}$ as (LL+SS) and to those which consider contributions from only ${\rho }^{\mathrm{LL},0}$ and do not use renormalised orbitals as (LL). The neglect of the contributions from ${\rho }^{\mathrm{SS},0}$ and from renormalisation is commonly referred to as the picture change error, which is given herein by the difference between LL+SS and LL results. It shall be noted that, for the EFG operator, corrections stemming from renormalisation remain essentially negligible as mainly valence orbitals contribute to the expectation value. One commonly accounts for nuclear electric quadrupole moments within first-order perturbation theory only, so that contributions from electric quadrupole moments of the various nuclei remain additive.

Quasi-relativistic two-component (2c) calculations are performed in this work within ZORA at the level of cGHF or cGKS with a modified version [7,41–46] of the quantum chemistry program package Turbomole [47] to determine the unperturbed HF or KS molecular orbitals.

Structure parameters of CHBrClF and CHClFI, as well as electronic densities at the level of two-component ZORA (2cZORA) were employed as described in Ref. [37]. Electronic densities were calculated on the level of one component ZORA-cGHF and ZORA-cGKS (1cZORA) and on the non-relativistic level of restricted HF (RHF) and restricted KS (RKS) with the same basis set employed in Ref. [37]. Density functionals used in the present study are the local density approximation (LDA) [48–50] as well as the Lee, Yang and Parr correlation functional (LYP) [51] combined with a generalised gradient exchange functional by Becke (BLYP) [52] or the hybrid Becke three parameter exchange functional (B3LYP)[49,53–55].

A ZORA-model potential $\tilde{V}$ as proposed by van Wüllen [28] was employed with additional damping [56].

In all self-consistent calculations of the unperturbed electron densities a finite nucleus was used, described by a normalised spherical Gaussian nuclear density distribution ${\rho }_{\mathrm{nuc},A}(\overrightarrow{r})=({\zeta }_A^{3/2}/{\pi }^{3/2}){\mathrm{e}}^{-{\zeta }_A|\overrightarrow{r}-{\overrightarrow{r}}_A{|}^2}$, where ${\zeta }_A=(3/(2r_{\mathrm{nuc},A}^2))$ and the root mean square radius $r_{\mathrm{nuc},A}$ of nucleus A was used as suggested by Visscher and Dyall [57]. The mass numbers A were chosen to correspond to the isotopes ${}^1$H, ${}^{12}$C, ${}^{19}$F, ${}^{35}$Cl, ${}^{79}$Br, ${}^{127}$I.

All electronic expectation values of electric field gradients were calculated by numerical integration on an ultra fine grid with our ZORA property toolbox approach described in Ref. [7]. Resulting EFG tensors were multiplied by the corresponding value of the NEQM taken from Ref. [2] and rotated into the principal axis system of the moment of inertia tensor following the same axis convention as used in Ref. [26] and Ref. [27], respectively (see Figure 1).

Figure 1. Choice of principal axes in CHBrClF (left) and CHClFI (right) in centre of mass coordinates. Equilibrium structures are shown as optimised in Ref. [37] at the level of CCSD(T) with small basis sets. The relative length chosen here for the principal axes corresponds to the relative size of the rotational constant. The a axis in CHClFI is scaled (shortened) by a factor of 2/3 relative to the b and c axes for a better representation. Atoms are labelled with mass numbers of used isotopes for ${}^{12}{\mathrm{C}}^1{\mathrm{H}}^{79}\mathrm{B}{\mathrm{r}}^{35}\mathrm{C}{\mathrm{l}}^{19}\mathrm{F}$ and ${}^{12}{\mathrm{C}}^1{\mathrm{H}}^{35}\mathrm{C}{\mathrm{l}}^{19}{\mathrm{F}}^{127}\mathrm{I}$.

3. Results

The rotational constants of the molecular structures obtained in Ref. [37] at the level of CCSD(T) with a small basis set are $A_{\mathrm{e}}=6394$ MHz, $B_{\mathrm{e}}=1941$ MHz and $C_{\mathrm{e}}=1537$ MHz for CHBrClF and $A_{\mathrm{e}}=6230$ MHz, $B_{\mathrm{e}}=1405$ MHz and $C_{\mathrm{e}}=1175$ MHz for CHClFI. Compared to the experimental rotational constants of the vibrational ground state reported in Refs. [26,27] these are off by 1 to 5 %. As EFGs can be rather sensitive to a change in molecular structure, this might slightly affect the accuracy of our results (see e.g. Ref. [11]). In Tables 1 and 2, the results for CHBrClF are listed and those of CHClFI can be found in Tables 3 and 4.

3.1. CHBrClF

Relativistic effects are not important for $\mathcal{X}({}^{35}\mathrm{Cl})$ and spin-orbit coupling effects are negligible for the NEQM coupling of both Cl and Br ($<1\mathrm{\%}$). Among the relativistic corrections, scalar-relativistic effects are most important (about $92\mathrm{\%}$ to $94\mathrm{\%}$ of the total relativistic effects for Br, $95\mathrm{\%}$ to $100\mathrm{\%}$ for Cl). The so-called picture change error has an influence of about 20 % on the total relativistic corrections, whereas spin-orbit coupling contributes only around 6 % to 8 % of the total relativistic effects for Br, and at most 5 % for Cl.

As there are only comparatively light nuclei (Z<50) in CHBrClF, electron correlation effects play a more prominent role than relativistic effects. Our calculations show that HF overestimates the absolute values of $\mathcal{X}$ consistently by about 15 %.

In order to gauge the dependence on the molecular structure and the basis set and to allow for direct comparison to post-HF results obtained in Ref. [58] on the second-order Møller–Plesset perturbation theory (MP2) level in the non-relativistic limit, we calculated also NEQM coupling tensors for the structure that was energy-minimised in Ref. [58] at the level of MP2/6-311+G(2df,2pd). These additional calculations were performed at the HF level with the large even tempered basis set (see methodology and Ref. [37]) (Hobi-G) and with the 6-311+G(2df,2pd) basis set (Hobi-B+G). Furthermore, with the structure of Ref. [37] employed for all other calculations, the non-relativistic NEQM tensor was calculated at the level of HF/6-311+G(2df,2pd) (Hobi-B). These calculations show a weak dependence on the basis set ($<3\mathrm{\%}$ for $\mathcal{X}({}^{35}\mathrm{Cl})$ and $<1\mathrm{\%}$ $\mathcal{X}({}^{79}\mathrm{Br})$) and a mild dependence on the molecular structure ($4\mathrm{\%}$ to $6\mathrm{\%}$ for the diagonal elements of $\mathcal{X}({}^{35}\mathrm{Cl})$ and $\mathcal{X}({}^{79}\mathrm{Br})$). The pure correlation effects, estimated as the difference between MP2 [58] and HF results, are found to be on the order of $\sim 5\mathrm{\%}$. This leads to deviations of our HF result (with basis and structure from Ref. [37]) from the MP2 result of Ref. [58] of up to $20\mathrm{\%}$ (aa component of $\mathcal{X}({}^{35}\mathrm{Cl})$).

We tested density functionals of different flavour: hybrid, generalised gradient approximation (GGA) and LDA functionals. As to be expected, the largest deviation was found between LDA and HF, being on average of about 5 % (largest deviation is 11 % for the bb component of $\mathcal{X}({}^{35}\mathrm{Cl})$). Electron correlation described by DFT is found to lower absolute values, which is in line with the correlation effects of MP2. Thus, our study shows that LDA gives the best results among the functionals considered here (up to 14 % deviation for $\mathcal{X}({}^{35}\mathrm{Cl})$ and up to 8 % deviation of $\mathcal{X}({}^{79}\mathrm{Br})$ from experiment).

Assuming the structure dependence studied at the level of HF is the same for DFT methods, the NEQM coupling tensors may get even closer to experimental results if the molecular structure of Ref. [58] is employed. For example, a correction of the 2cZORA(LL+SS)-LDA values by the difference between the HF/Hobi-G data and the HF results would yield deviations below $5\mathrm{\%}$ from experiment for all components of $\mathcal{X}({}^{79}\mathrm{Br})$.

3.2. CHClFI

Whereas relativistic contributions do not play an important role for $\mathcal{X}({}^{35}\mathrm{Cl})$ as in CHBrClF, relativistic effects become as expected very important for $\mathcal{X}({}^{127}\mathrm{I})$ and exceed electron correlation effects.

Scalar-relativistic effects are most important giving for $\mathcal{X}({}^{127}\mathrm{I})>90\mathrm{\%}$ of the total relativistic effects. Spin-orbit coupling contributes at most about $8\mathrm{\%}$. Again the picture change error is very large, making a difference of about 20 % of the total relativistic effects. The small influence of spin-orbit coupling can be understood as there is effectively only one heavier atom, and this heavy atom is a main group element (see Ref. [59]).

Again the LDA functional performs best of all tested density functionals and results obtained are in excellent agreement with experiment for the diagonal ($<3\mathrm{\%}$ deviation) and in good agreement for the off-diagonal tensor components (absolute deviations are of the same size, but relative deviations are between $5\mathrm{\%}$ and $12\mathrm{\%}$). Correlation effects are on the order of $8\mathrm{\%}$ and thus, contrary to the situation in CHBrClF or the NEQM coupling for Cl, less important than relativistic effects, which amount to about 15 %.

3.3. Discussion

We find surprisingly good results for NEQM coupling tensors of halogen nuclei in chiral methane derivatives when compared to experimental values. We emphasise, however, that we purposely did not use the currently most sophisticated computational approaches to determine the equilibrium structures of these molecule but reused structures determined previously  [37] on the CCSD(T) level with a comparatively small basis set to obtain consistent and comparable results for the different molecular properties. Moreover, we did not include any vibrational corrections. Although the latter can be expected to widely cancel out as the sign of nuclear and electronic contributions is opposite, there are cases where vibrational corrections can have a considerable influence on the calculated NEQM coupling tensor [60]. However, at the non-relativistic HF level we could show for an improved molecular structure that the values of the tensor are altered at most on the level of $6\mathrm{\%}$ and we do not expect vibrational effects to be much larger. Furthermore, these structure effects suggested even better agreement with experiment if the molecular structure of Ref. [58] is employed. We note in passing that agreement with experiment may also be improved by use of a NEQM that is scaled by the difference between experiment and theory obtained from calculation of a series of similar compounds as done e.g. in Ref. [61] for the NEQMs of ${}^{79}$Br and ${}^{81}$Br for different methods.

With the methods used herein, larger molecules can be calculated easily, and comparison of CHBrClF and the slightly heavier CHClFI suggests that the performance of the present approaches should improve in molecules with even heavier nuclei for which relativistic effects dominate correlation effects. As usual, some caution has to be taken when using DFT for the estimation of molecular properties, and a range of functionals should be applied in order to allow for an estimate of the importance of electron correlation effects. Results on the performance of different flavours of functionals can vary for different compound classes. For example our results for chiral halomethane derivatives deviate from the findings for mercury halides in Ref. [15], where hybrid functionals were found to outperform LDA and GGA functionals. In other systems, in particular those that contain the transition metal atom Cu, the performance of a variety of density functionals and especially LDA may, however, turn out to be not sufficiently reliable on the whole [62,63]. However, our present study indicates the applicability of 2cZORA-LDA for heavy-element containing chiral main group compounds and motivates future studies of NEQM couplings in these systems with our approach.

4. Conclusion and outlook

In the chiral methane derivatives CHBrClF and CHClFI, which serve as prototypical molecular probes for violations of fundamental symmetries, we have studied nuclear electric quadrupole moment couplings at the level of ZORA-GHF and ZORA-GKS with different density functionals. Our results show good agreement between two-component ZORA-LDA calculations of NEQM couplings and reported experimental values, which is in line with previous studies of main group compounds with DFT and approximate two-component methods. Our methods are easily applicable to larger molecules containing heavier elements and we can expect from our analysis an even improved performance for compounds with heavier nuclei. In future, the method used herein may be applied to other related chiral molecules that are spectroscopically relevant with regard to fundamental aspects of chirality, such as chiral polyhalocubanes [64] or polyhaloadamantanes [65]. Furthermore, this benchmark test enables us to provide support for upcoming microwave experiments with chiral heavy-elemental main group molecules to which more sophisticated methods are not directly applicable. One can also envision our methods to be employed within high-accuracy studies based on composite approaches, in which few single-point calculations for NEQM couplings are performed with coupled cluster methods, whereas vibrational corrections are accounted for with the help of increments determined with the present quasi-relativistic ZORA-DFT methods.

Acknowledgments

This paper is dedicated to Jürgen Gauss on the occasion of his 60th birthday. The authors thank Pekka Pyykkö for enlightening discussions about nuclear quadrupole moments and Melanie Schnell, Daniel Obenchain for raising our interest in electric field gradient calculations of chiral molecules with heavy nuclei. The authors gratefully acknowledge discussion with and help from Jürgen Stohner regarding nuclear electric quadrupole coupling constants in CHBrClF. We are much obliged to Oliver Kreuz, Carsten Zülch and Steffen Giesen for their support in testing our methods. Computer time provided by the Center for Scientific Computing (CSC) Frankfurt is gratefully acknowledged.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was financially supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) project no. 32896117 – SFB 1319 ELCH.

Appendix

Table A1. NEQM coupling tensor $\mathcal{X}$ of the ${}^{35}$Cl nucleus in CHBrClF calculated at level of Hartree–Fock (HF) and DFT with exchange-correlation functionals B3LYP, BLYP and LDA without inclusion of relativistic effects (NR) or relativistic effects included via the one-component (1cZORA) or two-component (2cZORA) ZORA approach. ZORA calculations of $\mathcal{X}$ are shown with (LL+SS) and without (LL) inclusion of picture change effects determined by consideration of small-component integrals and renormalisation. The NEQM of ${}^{35}$Cl was used as $Q^{35}\mathrm{Cl}=-8.11{\mathrm{fm}}^2$ as proposed in Ref. [2]. Non-relativistic HF values are compared to those calculated with the molecular structure reported in Ref. [58] (Hobi-G), or the 6-311+G(2df,2pd) basis set used in Ref. [58] (Hobi-B), or both (Hobi-B+G).

Method		$\frac{{\mathcal{X}}_{\mathrm{aa}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{bb}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{cc}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{ab}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{ac}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{bc}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$
NR	HF/Hobi-G	−22.6	−8.42	31.1	−54.7	24.2	19.9
	HF/Hobi-B+G	−21.0	−7.94	28.9	−51.1	22.6	18.6
	HF/Hobi-B	−23.3	−7.94	31.2	−53.7	23.1	19.1
	HF	−24.2	−8.23	32.4	−55.9	24.0	20.0
	B3LYP	−23.8	−7.68	31.5	−53.5	23.4	19.4
	BLYP	−23.4	−7.38	30.8	−52.2	22.9	19.0
	5LDA	−23.2	−7.32	30.5	−51.1	22.6	18.6
1cZORA (LL)	HF	−24.6	−8.32	32.9	−56.7	24.3	20.3
	B3LYP	−24.2	−7.75	32.0	−54.2	23.6	19.7
	BLYP	−23.8	−7.44	31.3	−52.9	23.1	19.2
	LDA	−23.6	−7.37	31.0	−51.7	22.8	18.9
1cZORA (LL+SS)	HF	−24.5	−8.31	32.8	−56.5	24.3	20.2
	B3LYP	−24.1	−7.73	31.9	−54.0	23.6	19.6
	BLYP	−23.8	−7.42	31.2	−52.7	23.0	19.2
	LDA	−23.5	−7.36	30.9	−51.6	22.7	18.8
2cZORA (LL)	HF	−24.6	−8.32	32.9	−56.7	24.3	20.3
	B3LYP	−24.2	−7.75	32.0	−54.2	23.6	19.7
	BLYP	−23.8	−7.44	31.3	−52.8	23.1	19.2
	LDA	−23.6	−7.38	31.0	−51.7	22.8	18.9
2cZORA (LL+SS)	HF	−24.5	−8.31	32.9	−56.5	24.3	20.2
	B3LYP	−24.1	−7.74	31.9	−54.0	23.6	19.6
	BLYP	−23.8	−7.43	31.2	−52.7	23.0	19.2
	LDA	−23.5	−7.37	30.9	−51.6	22.8	18.8
Experiment [26]		−20.60	−8.54	29.14	−51	23	19
NR-MP2 [58]${}^a$		−19.91	−7.41	27.31	47.92	21.87	−17.54

^aNon-relativistic second-order Møller–Plesset perturbation theory calculations (NR-MP2) performed with the 6-311+G(2df,2pd) basis in Ref. [58], scaled by a factor of $\frac{-8.11{\mathrm{fm}}^2}{-8.17{\mathrm{fm}}^2}\approx 0.994$ to account for the different value of $Q_{{}^{35}\mathrm{Cl}}$ employed in Ref. [58]

Table A2. NEQM coupling tensor $\mathcal{X}$ of the ${}^{79}$Br nucleus in CHBrClF calculated at level of Hartree–Fock (HF) and DFT with exchange-correlation functionals B3LYP, BLYP and LDA without inclusion of relativistic effects (NR) or relativistic effects included via the one-component (1cZORA) or two-component (2cZORA) ZORA approach. ZORA calculations of $\mathcal{X}$ are shown with (LL+SS) and without (LL) inclusion of picture change effects determined by consideration of small-component integrals and renormalisation. The NEQM of ${}^{79}$Br was used as $Q_{{}^{79}\mathrm{Br}}=30.87\mathrm{f}{\mathrm{m}}^2$ as proposed in Ref. [2]. Non-relativistic HF values are compared to those calculated with the molecular structure reported in Ref. [58] (Hobi-G), or the 6-311+G(2df,2pd) basis set used in Ref. [58] (Hobi-B), or both (Hobi-B+G).

Method		$\frac{{\mathcal{X}}_{\mathrm{aa}}({}^{79}\mathrm{Br})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{bb}}({}^{79}\mathrm{Br})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{cc}}({}^{79}\mathrm{Br})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{ab}}({}^{79}\mathrm{Br})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{ac}}({}^{79}\mathrm{Br})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{bc}}({}^{79}\mathrm{Br})}{h\mathrm{MHz}}$
NR	HF/Hobi-G	459	−213	−246	−250	191	−57.1
	HF/Hobi-B+G	459	−213	−246	−251	191	−57.6
	HF/Hobi-B	484	−225	−259	−253	190	−56.2
	HF	483	−225	−258	−252	189	−55.8
	B3LYP	471	−219	−253	−243	188	−54.1
	BLYP	461	−214	−247	−236	185	−52.5
	LDA	455	−210	−245	−233	184	−51.6
1cZORA (LL)	HF	522	−243	−279	−272	205	−60.5
	B3LYP	507	−235	−272	−261	203	−58.3
	BLYP	495	−230	−265	−254	199	−56.5
	LDA	489	−226	−263	−250	197	−55.6
cZORA (LL+SS)	HF	514	−239	−275	−268	202	−59.5
	B3LYP	499	−232	−267	−257	200	−57.4
	BLYP	487	−226	−261	−250	196	−55.7
	LDA	481	−222	−259	−246	194	−54.8
2cZORA (LL)	HF	522	−243	−279	−272	205	−60.5
	B3LYP	507	−235	−272	−261	203	−58.3
	BLYP	494	−230	−265	−254	199	−56.5
	LDA	489	−226	−263	−250	197	−55.6
2cZORA (LL+SS)	HF	516	−240	−276	−269	203	−59.8
	B3LYP	501	−233	−269	−258	200	−57.7
	BLYP	489	−227	−262	−251	196	−55.9
	LDA	483	−223	−260	−247	195	−55.0
Experiment [26]		456.07	−212.12	−243.94	−249.3	199	−59
NR-MP2 [58]${}^a$		431.19	−199.83	−231.38	234.39	184.51	53.25

^aNon-relativistic second-order Møller–Plesset perturbation theory calculations (NR-MP2) performed with the 6-311+G(2df,2pd) basis in Ref. [58], scaled by a factor of $\frac{30.87\mathrm{f}{\mathrm{m}}^2}{31.3\mathrm{f}{\mathrm{m}}^2}\approx 0.986$ to account for the different value of $Q_{{}^{79}\mathrm{Br}}$ employed in Ref. [58]

Table A3. NEQM coupling tensor $\mathcal{X}$ of the ${}^{35}$Cl nucleus in CHClFI calculated at level of Hartree–Fock (HF) and DFT with exchange-correlation functionals B3LYP, BLYP and LDA without inclusion of relativistic effects (NR) or relativistic effects included via the one-component (1cZORA) or two-component (2cZORA) ZORA approach. ZORA calculations of $\mathcal{X}$ are shown with (LL+SS) and without (LL) inclusion of picture change effects determined by consideration of small-component integrals and renormalisation. The NEQM of ${}^{35}$Cl was used as $Q_{{}^{35}\mathrm{Cl}}=-8.112\mathrm{f}{\mathrm{m}}^2$ as proposed in Ref. [2].

Method		$\frac{{\mathcal{X}}_{\mathrm{aa}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{bb}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{cc}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{ab}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{ac}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{bc}}({}^{35}\mathrm{Cl})}{h\mathrm{MHz}}$
NR	HF	−17.3	−14.8	32.1	−56.0	−23.2	−21.2
	B3LYP	−17.4	−13.9	31.3	−53.8	−22.6	−20.5
	BLYP	−17.2	−13.5	30.7	−52.5	−22.1	−20.0
	LDA	−17.3	−13.2	30.5	−51.6	−21.9	−19.6
1cZORA (LL)	HF	−17.8	−14.9	32.7	−56.8	−23.3	−21.4
	B3LYP	−17.8	−14.0	31.8	−54.4	−22.7	−20.8
	BLYP	−17.7	−13.5	31.2	−53.1	−22.2	−20.3
	LDA	−17.7	−13.3	31.0	−52.1	−21.9	−19.9
1cZORA (LL+SS)	HF	−17.7	−14.9	32.5	−56.6	−23.2	−21.4
	B3LYP	−17.8	−13.9	31.7	−54.3	−22.6	−20.7
	BLYP	−17.6	−13.5	31.1	−53.0	−22.1	−20.2
	LDA	−17.6	−13.2	30.9	−52.0	−21.8	−19.9
2cZORA (LL)	HF	−17.7	−14.9	32.6	−56.6	−23.3	−21.4
	B3LYP	−17.7	−14.0	31.7	−54.3	−22.6	−20.7
	BLYP	−17.6	−13.5	31.1	−53.0	−22.1	−20.3
	LDA	−17.6	−13.3	30.9	−52.0	−21.9	−19.9
2cZORA (LL+SS)	HF	−17.6	−14.9	32.5	−56.5	−23.2	−21.4
	B3LYP	−17.7	−14.0	31.6	−54.2	−22.6	−20.7
	BLYP	−17.5	−13.5	31.0	−52.9	−22.1	−20.2
	LDA	−17.5	−13.3	30.8	−51.9	−21.9	−19.9
Experiment [27]		−13.40	−15.45	28.85	−51.4	−22.7	−21.1

Table A4. NEQM coupling tensor $\mathcal{X}$ of the ${}^{127}$I nucleus in CHClFI calculated at level of Hartree–Fock (HF) and DFT with exchange-correlation functionals B3LYP, BLYP and LDA without inclusion of relativistic effects (NR) or relativistic effects included via the one-component (1cZORA) or two-component (2cZORA) ZORA approach. ZORA calculations of $\mathcal{X}$ are shown with (LL+SS) and without (LL) inclusion of picture change effects determined by consideration of small-component integrals and renormalisation. The NEQM of ${}^{127}$I was used as $Q_{{}^{127}\mathrm{I}}=-68.822\mathrm{f}{\mathrm{m}}^2$ as proposed in Ref. [2].

Method		$\frac{{\mathcal{X}}_{\mathrm{aa}}({}^{127}\mathrm{I})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{bb}}({}^{127}\mathrm{I})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{cc}}({}^{127}\mathrm{I})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{ab}}({}^{127}\mathrm{I})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{ac}}({}^{127}\mathrm{I})}{h\mathrm{MHz}}$	$\frac{{\mathcal{X}}_{\mathrm{bc}}({}^{127}\mathrm{I})}{h\mathrm{MHz}}$
NR	HF	−1620	792	827	633	524	−118
	B3LYP	−1560	763	797	602	517	−112
	BLYP	−1520	742	775	583	506	−108
	LDA	−1500	732	771	576	505	−106
1cZORA (LL)	HF	−1950	956	998	764	634	−144
	B3LYP	−1860	908	949	718	617	−135
	BLYP	−1800	879	918	691	600	−129
	LDA	−1780	868	914	684	600	−127
1cZORA (LL+SS)	HF	−1880	920	961	736	611	−139
	B3LYP	−1790	875	915	692	594	−130
	BLYP	−1730	848	885	667	579	−124
	LDA	−1720	837	882	659	578	−123
2cZORA (LL)	HF	−1960	957	1000	766	636	−144
	B3LYP	−1850	903	944	715	614	−134
	BLYP	−1790	874	912	687	597	−128
	LDA	−1770	862	908	679	596	−126
2cZORA (LL+SS)	HF	−1900	932	973	745	618	−141
	B3LYP	−1800	881	920	696	598	−131
	BLYP	−1740	852	889	670	582	−125
	LDA	−1730	841	886	662	581	−123
Experiment [27]		−1704.58	838.58	866.0	699.77	622.7	−140.4