Quasi-relativistic approach to analytical gradients of parity violating potentials

An analytic gradient approach for the computation of derivatives of parity-violating (PV) potentials with respect to displacements of the nuclei in chiral molecules is described and implemented within a quasirelativistic mean-field framework. Calculated PV potential gradients are utilised for estimating PV frequency splittings between enantiomers in rotational and vibrational spectra of four chiral polyhalomethanes, i.e. CHBrClF, CHClFI, CHBrFI and CHAtFI. Values calculated within the single-mode approximation for the frequency shifts agree well with previously reported theoretical values. The influence of non-separable anharmonic effects (multi-mode effects) on the vibrational frequency shifts, which are readily accessible with the present analytic derivative approach, are estimated for the C--F stretching fundamental of all four molecules and computed for each of the fundamentals in CHBrClF and CHAtFI. Multi-mode effects are found to be significant, in particular for the C--F stretching modes, being for some modes and cases of similar size as the single-mode contribution.


I. INTRODUCTION
Early after the proof of parity violation (PV) was provided in 1957 by the famous experiment of β − disintegration of Co nuclei 1 , which realised a proposal by the theoreticians Lee and Yang 2 , it was suggested by Yamagata 3 that PV weak interactions can induce a tiny energy difference between a chiral molecule and its non-identical mirror-image. [4][5][6][7][8][9][10][11][12] Diverse experimental schemes have been proposed to detect PV effects in chiral molecules, ranging from Mössbauer spectroscopy 13 to vibrational spectroscopy, [14][15][16][17] electron paramagnetic resonance spectroscopy, 18 rotational spectroscopy 16,19 and nuclear magnetic resonance spectroscopy [20][21][22][23] to timedependent approaches and quantum-beat experiments 12,[24][25][26] . We refer the reader for more details and an extended overview to a collections of reviews on this subject [27][28][29][30][31][32][33] . The most accurate experimental attempts reported so far were in the high-resolution infrared spectroscopy of bromochlorofluoromethane (CHBrClF). An upper bound of ∆ν PV /ν ≈ 10 −13 was obtained for the relative PV difference of the C-F stretching frequency between (R)-and (S)-enantiomers. 17 The theoretically [34][35][36][37][38][39][40][41] predicted value for this relative frequency splitting is of the order of 10 −17 , however. Later, with an improved experimental set-up, a measurement of the same compound with a resolution of 5 × 10 −14 has been reported in 2002 42 and with a new set-up, it may be hoped that a precision of 10 −16 can be reached. 43,44 As nuclear spin-independent electroweak PV effects in chiral compounds scale approximately with nuclear charge Z to the power of five (in the presence of spin-orbit coupling), compounds containing heavy metal nuclei could be of greater experimental value than the a) Parts of this work were reported in preliminary form in S. Brück, Diploma thesis, University of Frankfurt, 2011. originally used organic molecules. Theoretical searches in this direction have been already made on molecules containing for instance bismuth, rhenium, mercury and astatine. 40,[45][46][47] In general, a measurement of the PV energy as a difference between electronic energies of separated enantiomers would be very difficult since it comes on top of the rest mass energy of the molecule, 5 but the scheme proposed by Quack 12 , for instance, circumvents this by directly measuring the PV energy via the PV induced time-dependent interconversion between states of opposite parity 27 . An alternative for detecting molecular parity violation is to measure the difference between PV energy differences arising due to the PV potential [E PV ( q)] in the vibrational and rotational transitions in chiral compounds, 5 with q denoting the vector of dimensionless reduced normal coordinates. Often also V PV is used as symbol for the PV potential. In first order perturbation theory, E PV ( q) gives rise to a PV shift in the n th vibrational energy level of a given enantiomer (R or S) according to 37 E R,S n,PV ≈ Ψ R,S n |E PV ( q)|Ψ R,S n , where |Ψ R,S n denote the n th vibrational state of the R-and S-enantiomer, which are obtained by solving the parityconserving (PC) rovibrational Schrödinger equation for each enantiomer. Since, E PV ( q) is parity odd, E R n,PV and E S n,PV values are numerically of equal magnitude but have an opposite sign. Thus, the corresponding PV energy difference between the n th vibrational levels of the two enantiomers is ∆E n,PV = E S n,PV − E R n,PV ≈ 2 Ψ S n |E PV ( q)|Ψ S n .
The relative change in vibrational (∆E vib,PV = ∆E m,PV − ∆E n,PV ) and rotational (∆E rot,PV ) transition energies between right-and left-handed molecules is expected to scale to same order of magnitude as that of the relative change in electronic (∆E el,PV ) transition energy. 48 The introduction of PV within electroweak quantum chemistry not only affects the energy of an enantiomer, but also its equilibrium structure if the PC vibrational potential gets modified by the PV potential. 37 In general, the PV potential induces a minute change in the equilibrium structure of a molecule compared to the equilibrium structure without PV contribution, if the PV energy gradient ( ∇E PV ) does not vanish there. Thus, the weak interaction leads to a change of the equilibrium structure of a chiral molecule, which is different for both enantiomers due to PV effects.
This effect could in principle be measured by microwave spectroscopy, since it leads to a shift of the rotational constants. 34 In practice, the change of the structure due to ∇E PV at the minimum of the parity conserving potential can be estimated with the help of the vibrational Hessian F MW , 37 of the PC potential, given in mass-weighted Cartesian displacement coordinates. From this 3N nuc × 3N nucdimensional vibrational Hessian, with N nuc being the number of nuclei in the system, contributions of infinitesimal translational and rotational displacements can be eliminated with the 3N nuc × (3N nuc − 6) dimensional projection matrix A by forming (A T F MW A) −1 , so that the corresponding Hessian in internal Cartesian displacement coordinates results, which can be inverted. When M is a (3N nuc × 3N nuc ) diagonal matrix containing the masses of the various atoms, the final displacements (δ PV R) of the atoms from their equilibrium position due to ∇E PV read δ PV R = −M −1/2 A(A T F MW A) −1 A T M −1/2 ∇E PV . (4) Coordinates and displacements are transformed subsequently to the principal axis system of the PC equilibrium structure and the splittings of the diagonal elements (∆I) of the moment of inertia tensor due to PV effects are approximated to first order in δ PV R. For the xx component of ∆I we have for instance: and analogous for the other diagonal elements. Then, the changes in diagonal elements are used to estimate the splittings of rotational constants (∆X R ) between two enantiomers within the approximate expression 37 with X R being one of the rotational constants (A, B or C) and I X being the corresponding eigenvalues of the moment of inertia tensor. From Eq. 4 the importance of the gradient of the parity violating potential with respect to displacements of the nuclei becomes particularly evident. The numerical calculation of this term, however, is tedious, in particular within a (quasi)relativistic electronic structure framework, which is why we present in this work an approach for calculating this gradient analytically within a quasirelativistic meanfield framework.
Another promising experiment for detecting parity violation in chiral molecules is vibrational spectroscopy. Most of the effort has been put into the measurement of vibrational frequency shifts due to PV effects. 14,16,17,34,36,38,40,42,[49][50][51][52] The relative vibrational frequency shift for a transition from state n to m is determined by taking the difference of the vibrationally averaged PV potentials of the enantiomers and dividing by the corresponding transition energy hν mn 37,40 ∆ν mn,PV ν mn = 2 ( Ψ m |E PV ( q)|Ψ m − Ψ n |E PV ( q)|Ψ n ) hν mn (7) where the factor two enters into the equation since the difference of the R and S enantiomer is twice (cf. Eq. 2 ) the difference between one enantiomer and the parity conserving case, where no shift occurs. 40 The vibrationally averaged potentials depend on the multi-dimensional PV energy surface and not only on a single point energy at the equilibrium structure. In addition, the complete rovibrational wavefunction would be needed to compute the expectation value. As this problem is extremely complex even for relatively small molecules, the PV effects of electronic, vibrational and rotational degrees of freedom are typically separated in a first step. Still the computational effort is in general too high, so as a second step, the multi-dimensional problem is usually split into one-dimensional problems, where the movements along the normal coordinates are treated separately as if they were independent of each other. 37,40 This can be augmented by adding contributions from the potential depending on a smaller number of modes order by order. 53,54 In practice, the PC potential (V BO ) and PV potential [E PV ( q)] can be evaluated at onedimensional cuts along the dimensionless reduced normal coordinates q. 37,40 This approach has already been used to estimate the vibrational frequency splitting of the C-F stretching mode in chiral halogenated methane derivatives. 34,35,37,38,40,52 A different approach comes from perturbation theory, 55 where the PC and PV potentials are expanded in a Taylor series. In this, the influence of multi-mode effects are estimated by the calculations of the derivatives of the PV potentials with respect to all normal coordinates. The contributions of lowest non-vanishing order to the n th vibrational energy levels of a normal mode ν r within the perturbative treatment are Ψ n r |E PV ( q)|Ψ n r ≈ E PV ( R 0 ) where E PV ( R 0 ) is the PV energy at the equilibrium structure and ω s denotes the harmonic vibrational angular frequency of normal mode ν s . The terms ∂ E PV ∂ q s and ∂ 2 E PV ∂ q 2 r are the first (gradient) and second partial derivatives of the PV energy along normal coordinates (q s , q r ), respectively. Ignoring all the r = s terms in the second sum of Eq. 8 gives the one-dimensional (1D) perturbative estimate of the vibrational energy for normal mode ν r as Ψ n r |E PV (q r )|Ψ n r 1D = E PV ( R 0 ) The drawback of this method is that the one-dimensional terms of the PC and PV potentials are not included to infinite order. But the advantage is that, once the Cartesian gradient of the PV energy at the equilibrium structure is known, the two-dimensional coupling term can be included without too much effort, because only the semidiagonal cubic force constants are needed, whereas the gradient of PV energy along all the normal mode can conveniently be determined by projecting the Cartesian PV energy gradient onto displacements along to the normal coordinates. Hence, contributions from other modes (cf. Eq. 10) in addition to the single-mode terms described by Eq. 9 can easily be accounted for by considering first order multi-mode (MM) effects on the vibrational energy levels as shown in Eq. 8 We define with these MM terms the perturbative approximation of the vibrational energy shifts with 2D coupling terms as Thus, for the calculation of the influence of PV on the vibrational and rotational spectra, which is one of the main goals of the current manuscript, the knowledge of the gradient of the PV potential along all the normal modes is essential. The present article describes an analytical approach for calculating the gradient of the PV energy ( ∇E PV ) at HF and LDA level of theory (see the following section II). These values obtained from this analytic gradient approach, with corresponding computational details being described in section III, are then utilised in section IV for determining the relative shifts in rotational constants and vibrational frequencies, due to PV effects in chiral methane derivatives. Besides this, the influence of non-separable anharmonic effects (multimode effects) on PV induced vibrational frequency shifts for all the vibrational normal modes is calculated for CHBrClF and CHAtFI molecules.

II. THEORY
For the description of the vibrational movement of nuclei in the Born-Oppenheimer approximation, it is often sufficient to evaluate the electronic structure for a single fixed arrangement of the nuclei described by the coordinates R 0 and to treat the displacement of nuclei perturbatively (see also Refs. [56][57][58].
The electroweak parity-violating potential is very small compared to the PC potential due to the appearance of the Fermi coupling constant which is 2.22249 × 10 −14 E h a 3 0 in atomic units. Therefore, we can treat the parity-violating electroweak HamiltonianĤ PV as a small addition to the parity conserving molecular HamiltonianĤ 0 : where we have introduced λ PV as a formal perturbation parameter.
In a mean-field approach the leading order parity violating contribution E PV ( R 0 + η, λ PV ) = E( R 0 + η, λ PV ) − E( R 0 ) to the gradient of the variational energy with respect to nuclear displacements is: For details on the specific case of electroweak parity violation, see Appendix A, whereas for the general case of variational perturbation theory, see e.g. Ref. 59.
In this paper, we want to focus on the molecular Hamiltonian approximated within a (quasi-relativistic) two-component zeroth-order regular approximation (ZORA) framework 60,61 on the level of Generalised Hartree-Fock (GHF) or Generalised Kohn-Sham (GKS) density functional theory (DFT) (see Ref. 62 and 63): Here, V GHF,GKS is the effective electron repulsion potential within the GHF or GKS framework,ˆ p is the linear momentum operator, σ is the vector of Pauli spin matrices and ω( r) = 1 2c 2 −Ṽ ( r) is the ZORA factor with the ZORA model potentialṼ as proposed by van Wüllen to alleviate the gauge dependence of ZORA. 62 The ZORA one-electron operator has an electron spin-independent and an electron spin-dependent contribution: Within ZORA, the nuclear-spin independent parityviolating electroweak one-electron Hamiltonian appears as 64 where {A, B} + = AB + BA is the anti-commutator, Q W,A is the weak charge of nucleus A and ρ A is the normalised nuclear density distribution. This operator has only electron spindependent contributions. We expand the ZORA two-component HF or KS molecular orbitals (MOs) φ i in a linear combination of real onecomponent basis functions χ µ and complex two-component In this two-component framework, we can define four complex one-component density matrices (κ = 0, 1, 2, 3): where the 0 th component of the Pauli spin matrices is the 2 × 2 identity matrix. The two-component density matrix can be written in terms of these one-component density matrices as We can write the expectation value of the one electron ZORA operator as and the energy contribution of the parity-violating electroweak potential as In the following, the effective potential V GHF,GKS will be represented in the space of basis functions in terms of the matrix of contracted two-electron integrals G. In the general case of hybrid DFT, G (κ) is constructed as where the Mulliken notation for two electron integrals is employed: In case of pure DFT (non-hybrid) we have a X = 0 and in case of pure HF we have a X = 1 and a DFT = 0.
We consider non-relativistic density functionals which do not depend on the current density, which are commonly employed even in relativistic electronic structure theory. In this case, the matrix elements of the exchange-correlation potential χ µ V (κ) XC χ ν are always real. In this paper, we restrict the discussion to the spin-unpolarised local density approximation (LDA), in which the exchange-correlation potential has the form with F XC,LDA being the LDA density functional and the electronic number density function being For the general form of the exchange-correlation potential see Appendix B. The expectation value of the electron repulsion potential is The total ZORA energy in presence of the full perturbation H PV , which implies λ PV = 1, is given by where, ∞ refers to density matrices obtained variationally in the presence of the perturbation, which is analogous to an infinite order perturbation theory treatment ofĤ PV (for details see Eq. (A3) in the Appendix), provided that this converges. Thus, the gradient with respect to nuclear displacements can be written as For a first order property that does not depend on the basis functions such as the parity-violating potential, this expres-sion can be simplified by using the orthonormality condition of the HF equations and with Eq. (13) we receive the leading order parity-violating energy gradient (see Appendix C): where we have introduced the energy weighted density matrix (EWDM) W as and the matrix G (κ) grad of contracted gradients of two-electron integrals with elements: Note, that in the DFT case, derivatives of the exchangecorrelation potential have to be computed for the calculation of the perturbed two-electron gradients, which in the spinunpolarised LDA case are V XC,LDA (D (0) , D (0) ) = δV XC,LDA (D (0) ) δ ρ( r; D (0) ) ρ( r; D (0) ) . (31) For other functionals see Appendix B. D is the SCF density matrix of the unperturbed system, i.e. received with the HamiltonianĤ 0 only, and D and W are the perturbed density matrix and EWDM of first order in λ PV . The nuclear displacement gradients of the full unperturbed ZORA Hamiltonian, i.e. the gradients of one-electron ZORA integrals ∇ η h (κ) ZORA,µν , as well as gradients of two-electron integrals and the exchange-correlation potentials needed for G (κ) grad,µν (D) have been implemented in Ref. 65. In this work, we have altered the previous one-electron gradient part to account also for a finite nucleus model when computing the derivatives of the one-electron integrals. What remains then to be obtained are the perturbed density matrices D (κ) , the perturbed spin-independent EWDM W (0) and the gradient of the PV integrals ∇ η h (κ) PV,µν . We describe the scheme for computing the latter in the following subsection II A, before we discuss in subsection II B the linear response scheme that is used for computation of perturbed density matrices.

A. Gradient of the PV integrals
The matrix elements of the one electron PV operator in basis set representation is with ∂ κ = ∂ ∂ x κ denoting a component of the electronic four derivative, with the four vector being x (0,1,2,3) = (ct, x, y, z) T . Indices A, B,C denote the nuclei at which the basis functions χ or nuclear density distribution ρ nuc are hooked. As the basis functions depend on the nuclear coordinates, the geometry gradient of the PV integrals is where, the second term is the Hellmann-Feynman term and the first and third term are the basis function contributions. The integrals of the PV operator are purely imaginary and due to Hermiticity of the PV operator, the basis function contributions are connected by We split the nuclear displacements η into the separate contributions from different nuclei A. The basis functions depend on r − R A . Therefore, only the subset of basis functions centered at A contributes to the corresponding integrals. Furthermore, as we deal with Gaussian basis functions, derivatives with respect to nuclear coordinates can be represented by derivatives with respect to electronic coordinates as ∂ A,κ χ B = −δ AB ∂ κ χ B , where ∂ A,κ denotes a component of the four derivative of nucleus A. For the basis function contribution to the PV gradient integrals, we arrive at For a Gaussian shaped nuclear density distribution, the Hellmann-Feynman term reads with the gradient of a Gaussian normalised nuclear density distribution being (37) The geometry gradient of the ZORA factor ∇ A ω is discussed in Ref. 62. The construction of ZORA, PV operators and its derivatives is done by means of numerical integration.
The above integrals are evaluated numerically on a grid. In the present work, an atom centered grids is used employing the Treuter-Ahlrichs 66 version of a Becke grid, 67 where the partitioning is done by weight functions w i for each grid point i. As these weight functions depend on the nuclear coordinates they give an additional contribution to the PV gradient which is calculated as The grid points move when the nuclei are slightly displaced. To account for this effect, we employ that translational invariance of the molecule holds and, therewith, the net force on the molecule has to be zero. We, therefore, subtract the net force that results from the numerical integration procedure from the numerically integrated gradient contribution PV,A,µν ) num :

B. Linear response computation of perturbed density matrices
In first order the perturbed density matrix can be written in terms of the unoccupied-occupied block T uo of an anti-Hermitian transformation matrix T (see Appendix D 1 for details): For how to compute perturbed density matrices to arbitrary order see Ref. 68.
We introduce the MO transformed PV operator as If we would assume that the first order perturbed wave function would be calculated by transformation of the orbital coefficients as then the perturbed density matrix would simply calculated as D However, this does not account for the response of the orbitals to the perturbation but correspond to a simple sum over states approach in which it is assumed that the electronic Hessian is diagonal and therefore the perturbed SCF equations are uncoupled. Within HF and KS, however, the electronic Hessian is not diagonal as the two-electron matrix G is a function of the orbitals and one has to solve the coupled perturbed HF (CPHF) or coupled perturbed KS (CPKS) equations (see Appendix D 2 and for a detailed derivation e.g. Ref. 69). To solve the response equations, we use the reduced form (see Refs. 70 and 71 and Appendix D 2 for details): where indices b run over unoccupied orbitals and indices j over occupied orbitals. The elements of the electronic Hessian are withG ai,b j being an element of the four-index MO transformed two-electron tensor: with the transition density matrix D ai having the elements Eq. (43) is solved iteratively within a preconditioned conjugate gradient algorithm. Thereby, as initial guess (0), we employ trial vectorsT that represent the uncoupled solutions: In an iterative procedure in each step i from the trial vector T (i−1) , we construct the perturbed density matrices (D (κ) ) (i) following Eq. (40) and calculate the contracted electronic Hessian as with the two-index MO transformed two-electron matrix: From this we update For all calculations in this paper, the Fletcher-Reeves weight of the precondition was employed, which is calculated as β FR = r (i) r (i−1) and the algorithm is followed until convergence of the norm r (i) r (0) to a given threshold. The trial matrix is updated as Upon convergence, the solution is received as the sum over all weighted trial matrices from the N iter iterations: From the solution vectors of the linear response equations T, new perturbed density matrices are calculated following Eq. (40) and the perturbed EWDM can be calculated as with the MO transformed Fock matrix being In case of DFT contributions from the exchange correlation potential to the perturbed G matrix have to be calculated via derivatives of the exchange correlation potential as shown as in 31 for the case of spin unpolarised LDA; for a more general discussion of exchange-correlation functionals see Appendix B and Ref. 70).
The calculation of the two-electron part (HF case and Coulomb contribution in LDA) of the perturbed contracted two-electron matrix G is carried out via contraction of the two electron tensor (µν|ρσ ) with the two-particle density matrix Γ (κ,λ ) which can be constructed from the one-particle density matrices as Note, that for the employed ZORA operator without two-electron spin-orbit or spin-spin coupling terms, only twoelectron densities with κ = λ are required at the HF or DFT level. Furthermore, we introduced the scaling factors for the direct Coulomb contribution a κ,λ C , which is in all present calculations set as a κ,λ C = δ 0κ δ 0λ and a more flexible scaling parameter for the exchange contribution a κ,λ X , which is in the present implementation set to be constant a κ,λ X = a X . For the calculation of the two-electron contribution to the PV energy gradient, the perturbed two particle density matrices are needed: Whereas in the present implementation Eq. (62) is employed directly, in the pilot implementation the perturbed twoelectron density matrices were calculated for practical reasons via where a constant, numerical scaling factor of 2 16 was used to scale up artificially the very small numerical values in the perturbed density matrices to avoid substractive cancellation in Eq. 63.
The present implementation of ∇E PV as defined in Eq. (28) as well as the response equations [Eqs. (48)(49)(50)(51)(52)(53)(54)(55)(56)(57)(58)] were included in a modified version 63,77 of the Turbomole program. 74,75 The results from the current production-level implementation are identical to those of the pilot implementation. The two different implementations provided an internal test of our results. The implementation proceeds as follows: 1. For a given molecular structure, unperturbed LCAO coefficients and orbital energies are received from a SCF computation at the two-component ZORA level with the modified version 62 ZORA,µν and of the two-electron integrals ∇ η (µν|ρσ ) are computed and combined with the unperturbed and perturbed density matrices following Eq. (28). For two-electron contributions the perturbed two electron matrix is formed [Eq. (62)] and contracted with the gradient of two-electron integrals ∇ η (µν|ρσ ).
In the following, the scheme described above for computing ∇E PV is used to study electroweak parity-violating effects in halogenated methane derivatives; CHBrClF, CHClFI, CHBrFI and CHAtFI. As mentioned in the Introduction, their vibrational spectra, in particular for CHBrClF have been extensively studied theoretically [34][35][36]38,40,51,52 and experimentally 14,16,17,42,49 to detect PV effects. A high resolution of 5×10 −14 has been achieved for the C-F stretching mode of CHBrClF with CO 2 laser spectroscopy, 42 which is nevertheless about three or four orders of magnitude larger than the theoretical predictions of the size of the effect for the molecule under investigation. Therefore, the vibrational and rotational spectra of heavier homologues will be studied herein, as has been done previously for the vibrational frequency shifts in a one-dimensional anharmonic approximation. 40 Whereas in the previous work single-mode anharmonic effects were included variationally by solving the one-dimensional anharmonic vibrational Schrödinger equation, we will use in the present work vibrational perturbation theory to account for anharmonic effects. Most importantly, we can include now also multi-mode contributions to parity-violating frequency shifts efficiently, which were neglected in essentially all previous studies on parity-violating frequency shifts in chiral molecules with the notable exceptions of Ref. 53, where CDBrClF was studied in a four-dimensional anharmonic model that span the C-F stretching as well as C-D stretching and the two C-D bending modes and of Ref. 78, where up to third order vibrational effects in chiral arsenic and lead compounds were studied.
For all the methane derivatives mentioned above, we use the same molecular structures, electronic structure methods and basis sets from earlier work 40 to allow for direct comparison. As in Ref. 40, the S-enantiomers of a given chiral compound is considered when parity-violating potentials, and in the current work also gradients of the parity-violating potentials, are reported. Splittings of a property A between enantiomers are given as ∆A = A S − A R . The equilibrium structures, harmonic vibrational frequencies and potential energy surfaces (PES) were computed on the CCSD(T) level with the cc-pVDZ basis set for the first to third row elements and quasi-relativistic Stuttgart pseudopotentials in the neutral atom reference system together with energy optimised valence basis sets for Br, I 79 and At. 80 Equilibriums structures, harmonic vibrational frequencies, corresponding normal coordinates and displaced structures as reported in Ref. 40 were reused in the present work. The PES have been determined by the SURF module 54 of MOLPRO. [81][82][83] . The double zeta basis is expected to behave comparatively poorly, but is kept herein to speed up the calculations and to obtain results that are directly comparable to previous work. An improved description of the anharmonic PES, in particular for the promising astatine containing methane derivative will be left for a later study.
In computations of the PV operator a Weinberg parameter of sin 2 θ W = 0.2319 has been used to determine the weak nu- For LDA, Dirac exchange 85 and VWN5 correlation 86 potentials have been used. For LDA a standard DFT integration grid was used, whereas matrix elements of the ZORA and PV operators were computed on a very dense grid. MOs have been converged until the change of the SCF energy and relative change of spin-orbit energy (except CHClFI at HF level) between two-successive iterations dropped below at least 10 −6 E h and 10 −12 respectively. MOs for CHClFI at HF level have been converged to less than 10 −15 for the relative spin-orbit energy change as the corresponding PV energy for the equilibrium structure was not converged to the desired accuracy with the 10 −12 criterion. In practice, the spin-orbit energy criterion was by far the more restrictive one, such that at the end of the iterative process, the change in SCF energy between two cycles typically dropped below 10 −9 E h . The threshold for negelection of gradients of two-electron integrals was set to 10 −15 E h a −1 0 . Following the regular spectroscopic notation, the normal modes of the methane derivatives are named in descending order of the frequency values.

IV. RESULTS AND DISCUSSION
We first make a comparison of the directional derivatives of the PV energy along the C-F stretching mode (ν 4 ) as obtained from the analytical PV energy gradients with those of the numerical ones from an earlier study. 40 The PV energy gradients are utilised for estimating the shifts of the rotational constants corresponding to the equilibrium structures. The vibrational frequency shifts in vibrational transitions for the C-F stretching mode are calculated within a perturbative treatment using the PV energy gradients. At the end, the multi-mode effects in the vibrational transitions for all the normal modes in CHBr-ClF and CHAtFI molecules are discussed.

A. PV energy gradients
All the chiral methane derivatives studied here show a characteristic C-F stretching mode (ν 4 ) which is amenable to highresolution CO 2 laser spectroscopy. The directional derivative of the PV energy along the C-F stretching normal coordinate (q 4 ) is obtained herein by projecting the analytical Cartesian PV gradient ∇E PV onto the Cartesian displacement vectorû q 4 corresponding to a unit shift along the dimensionless reduced normal coordinate q 4 . Table I depicts the HF and LDA level directional derivative along C-F stretching mode at the equilibrium structure. See the Supporting Information for the corresponding Cartesian HF and LDA level ∇E PV for each of the equilibrium structures.
In general, the analytically computed gradient is expected to be more accurate than a numerical counterpart obtained by a finite difference scheme, provided that the underlying self-consistent field solutions are well converged. We note in reference to Table I, that a direct comparison of the analytical and numerical directional derivative of the parity violating potential is partially hampered by the fact that the analytical gradient is calculated at the equilibrium structure of the molecule, whereas the numerical directional derivative was determined in Ref. 40 from the linear term of a polynomial fit of the parity-violating potential as calculated along a onedimensional cut along the C-F stretching normal coordinate q 4 in the range from −3 to 3. The two become better comparable, if one computes also the analytical directional derivatives at several points along the same normal mode displacement and performs a polynomial fit. Therefore, as a next step, the coordinate dependence of the PV gradient is examined, checking the higher order terms as well, by fitting the ∇E PV ·û q r to the third order polynomial.  (7) a Numerical derivatives obtained in Ref. 40 from a polynomial fit of a one-dimensional cut through the parity-violating potential along the dimensionless reduced normal coordinates q 4 . For the equilibrium structure of CHClFI, however, the HF value was slightly less tightly converged, so that this value was recomputed herein (see supplement) and the potential refitted, leading to only slightly improved fit values.
whereû q r is the Cartesian displacement vector corresponding to a unit displacement along the normal coordinate q r and a i would correspond to 1 in a Taylor series expansion of ∇E PV ·û q r . Thus, a 0 is similar to the value of the ∇E PV ·û q r at the equilibrium structure along the normal mode ν r . The other fitting coefficients a 1 , a 2 and a 3 are related to the first, second and third order derivatives, respectively, of ∇E PV ·û q r with respect to the dimensionless reduced normal coordinate q r for ν r mode. For this purpose, structures that are displaced along the normal coordinates of a vibrational mode are employed to generate a one-dimensional cut through the PV potential energy surface along q that varies from −3 to 3 through the equilibrium structure (q = 0). A total of 17 points is taken into consideration. Here, for each of the methane derivatives, we consider the one-dimensional cut for the C-F stretching mode (ν 4 ). Figures 1-4 display the HF and LDA level ∇E PV (q 4 ) ·û q 4 along the dimensionless reduced normal coordinates corresponding to the C-F stretching mode of the methane derivatives. The corresponding numerical values of ∇E PV (q 4 ) ·û q 4 are given in Tables S15-S18 in the Supporting Information file. The fitting coefficients need to be adapted to compare the structure dependence of the PV energy E PV ( q r ) with the directional derivative of the PV energy ∇E PV ·û q 4 along the dimensionless reduced normal coordinates q 4 corresponding to the C-F stretching mode, because the former can be fitted to a polynomial expansion, too, namely 0.000 60 (5) 0.000 60 (7)  FIG. 1. HF and LDA level PV energy gradient (fitted to ∑ 3 i=0 a i q i r ) along the C-F stretching mode of CHBrClF. See Table II for the fitting coefficients (a i ) due to HF PV energy gradient. For the fitting coefficients of LDA PV energy gradient, see Table S23 in the Supporting Information. where b i would correspond to 1 in a Taylor series expansion of E PV . Thus, one has to contrast E PV ( q r ) with ∇E PV (q r ) ·û q r dq r and hence the values of a 0 , a 1 /2, a 2 /3 and a 3 /4 are similar to the values of b 1 , b 2 , b 3 and b 4 respectively. The term b 0 corresponds in this approximation to the PV energy at the equilibrium structure.
The comparison of these fitting coefficients for PV energy (taken from Ref. 40) and directional derivative of the PV energy along the C-F stretching mode for the HF level is presented in Table II. Both numerical and analytical ∇E PV ·û q 4 values for all the systems are in good agreement, with the numbers in parentheses denoting the error due to the fit pro-FIG. 2. HF and LDA level PV energy gradient (fitted to ∑ 3 i=0 a i q i r ) along the C-F stretching mode of CHClFI. See Table II for the fitting coefficients (a i ) due to HF PV energy gradient. For the fitting coefficients of LDA PV energy gradient, see Table S23 in the Supporting Information.
FIG. 3. HF and LDA level PV energy gradient (fitted to ∑ 3 i=0 a i q i r ) along the C-F stretching mode of CHBrFI. See Table II for the fitting coefficients (a i ) due to HF PV energy gradient. For the fitting coefficients of LDA PV energy gradient, see Table S23 in the Supporting  Information. cedure, which involved one free parameter less for the directional derivative (eq. 64) than the PV potential (eq. 65). The corresponding comparisons of these fitting PV energy and PV energy gradient coefficients along the C-F stretching mode at the LDA level are given in the Supporting Information.
For the determination of the second partial derivatives FIG. 4. HF and LDA level PV energy gradient (fitted to ∑ 3 i=0 a i q i r ) along the C-F stretching mode of CHAtFI. See Table II for the fitting coefficients (a i ) due to HF PV energy gradient. For the fitting coefficients of LDA PV energy gradient, see Table S23 in the Supporting Information. Table II), which are needed for calculating the vibrational energy levels as shown in Eqs. 8 and 9, one can expect the availability of the analytical gradient to simplify the calculations, because fewer points are needed from the PV energy surface. As a consistency check, a linear fit of the analytically calculated directional derivatives computed at q = −0.125, q = 0 and q = 0.125 has been opted for in all of the molecules, where the fitting values (a 1 ) are found to be in good agreement with the values presented in Table II. These linear fit values for the C-F stretching mode are reported [Table S22] in the Supporting Information.

B. Shifts of the rotational constants
The values shown in Tables I and II reveal that the equilibrium structure of these halogenated methane derivatives possess non-zero ∇E PV . Due to these non-zero ∇E PV values, the PV potential can induces a minute change in the equilibrium structure. This leads to a shift of the rotational constants, which could in principle be measured by microwave spectroscopy. As already been discussed in the Introduction, the change of the structure and hence the shifts in the rotational constants due to the existence of non-zero ∇E PV at the minimum of the parity conserving potential is calculated with the help of the vibrational Hessian F (see Eqs. 4 and 6). The change of the inertia tensor (∆I x ) in the principal axis system and subsequently the change of the rotational constants (∆X R ) can be approximated assuming that the result depends linearly on the displacements. 37 Shifts in the rotational constants as reported in Table III are highly sensitive to the PV energy gradients ∇E PV . The numbers obtained on the LDA levels are found to be about one quarter of the HF results for the CHBrClF molecule. The calculated relative effect in the microwave spectrum (∆X/X ≈ 10 −17 ) of CHBrClF, although it agrees well with the earlier finding by Quack and Stohner 34 , is still far from the current experimental resolution. Similarly, the corresponding relative shifts in CHClFI and CHBrFI are about one order of magnitude larger in absolute value than for CHBrClF, but still far below the present experimental resolution (see Table III). Considering the current scenario, the PV induced shifts of the rotational constants in CHBrClF, CHClFI and CHBrFI are not expected to be measurable with present experimental setups. On the other hand, due to the heavy mass of astatine, a larger shift in the rotational constants is expected for CHAtFI. For this heavier At analogue, shifts are found to be about two or three-orders of magnitude higher in absolute value than the lighter molecules considered herein.
Another comment as to the accuracy of these estimates of rotational energy shifts is in order: Herein, we have followed the most simple approach 37 and computed only the PV shift at the equilibrium structure as it would arise from the PV gradient contribution. An improved treatment would include estimates of PV induced shifts of vibrationally averaged rotational constants as well as their influence on Coriolis coupling terms and centrifugal distortion constants.

C. Vibrational transitions for C-F stretching mode
For the vibrational spectrum, the results from the separable anharmonic adiabatic approximation (SAAA) 37 used in earlier work 40 have been compared to the one-dimensional (1D) perturbative approach (Eq. 9) of the present work. The influence of multi-mode coupling terms to the vibrational energy levels is also accounted for after adding the perturbative multimode contribution (Eq. 10) to the single-mode results. Table  IV presents the vibrationally averaged HF and LDA parity violating potential energy levels up to 3 rd vibrational state for the C-F stretching mode (ν 4 ) of the S-enantiomer of CHBr-ClF, CHClFI, CHBrFI and CHAtFI molecules along with a comparison to SAAA results. 40 For the lowest transition (from n 4 = 0 to n 4 = 1) in CHBr-ClF, the full 1D and the perturbative 1D treatment match to about < 1 %, which is far below the accuracy of the employed ZORA-cGKS and ZORA-cGHF methods and therefore negligible. The transition after perturbative inclusion of multi-mode effects up to 2D coupling terms, however, differs by about 50 % (cf . Table IV). A similar influence of nonseparable anharmonic effects on PV has been reported for similar molecules. 53 Nonetheless, the conclusion to be drawn is that the effect of the 2D treatment of the parity conserving potential is more important than the higher order terms of the 1D parity-violating potential. When comparing the different electronic structure methods, one finds significant differences, but the values for vibrational frequency shifts in CHBr-ClF are in the same order of magnitude, deviating less than 30%, with the same sign, the latter of which is not the case for the PV energy shifts. When compared with absolute values, the 2D contributions are found to increase for the higher vibrational levels. When these calculated PV frequency shifts (Table IV) are multiplied by two and then divided by the corresponding vibrational transition frequencies, they produce the dimensionless parity-violating relative vibrational frequency splittings between the C-F stretching fundamental of S-and R-enantiomers. The calculated relative vibrational frequency splittings corresponding to the lowest transition for the C-F stretching vibration of CHBrClF molecule is found to be in good agreement with the earlier theoretical values. [34][35][36]38,40,52 For CHClFI, these parity-violating 1 ← 0 transition frequency shifts are found to be about one order of magnitude higher than CHBrClF (See Table IV). Again, the shifts in vibrational energy levels from the full and perturbed 1D treatment match better with each other, whereas the values obtained after inclusion of 2D contributions differ by a larger margin as compared to perturbative 1D treatment. Substitution with the heavier element bromine to CHClFI in place of chlorine leads to an increase of the PV transition frequency shifts compared to CHClFI again, which is not quite an order of magnitude this time (cf. Table IV). PV transition frequency shifts in the C-F stretching fundamental of CHBrFI due to full 1D and perturbed 1D are roughly twice to that of the values for CHClFI, whereas perturbed 2D effects in CHBrFI measures about 75% higher than the 2D contributions in CHClFI molecule. Akin to CHBrClF, the difference between the full 1D and perturbed 1D approach is again negligible, whereas 2D contributions are significant.
A similar observation is also noticed for CHAtFI, which shows the largest absolute value of the PV vibrational frequency difference in the C-F stretching fundamental among the four halogenated methane derivatives discussed herein. This large shift is expected and mainly attributed to the presence of the heavier astatine nucleus. The full 1D and perturbative 1D treatment do not differ much. When taking the 2D terms into account, the values change significantly with the multi-mode contributions causing a reduction in absolute value compare to the perturbative 1D estimate of the frequency shift. This is expected to result from the fact, that the couplings of the C-F stretching mode ν 4 to the C-H bending mode ν 3 , where the frequencies are close (∆ω ≈ 1.5%), and to the C-H stretching mode ν 1 , where there is a factor of about three (≈ 2.9) between the fundamental frequencies are very strong. In the perturbative approach for LDA, even the sign and the order of magnitude of the effects for the fundamental transition (from n 4 = 0 to n 4 = 1) are altered.
Thus, the 2D perturbative terms are found to be important in all cases, which increase gradually for higher vibrational levels. Thus, the multi-mode effects due to the perturbative 2D approach increases the fundamental vibrational frequency shift for the CHBrClF, CHClFI molecules whereas the shifts in CHBrFI molecule are lower in absolute value as compared to the perturbative 1D treatment counterparts. The heaviest At derivative shows the maximum frequency shift after the inclusion of the perturbative 2D effects at both HF and LDA level.

D. Multimode effects in CHBrClF and CHAtFI
The influence of multi-mode effects 53,87,88 from nonseparable anharmonic effects have been estimated for CHBr-ClF and the heaviest CHAtFI molecules after calculating the PV energy gradients along all the normal modes. Figures 5  and 6, respectively, display the variation of the HF and LDA PV energies along the dimensionless reduced normal coordinates of all the normal modes of CHBrClF molecule. In both cases, normal coordinates q 8 , q 7 , q 3 and q 2 show steep slopes for the PV energy with respect to displacements along the respective dimensionless reduced normal coordinates. This may be naively expected, because the bending motion can alter the chiral structure of the molecule to a larger extent than a stretching vibration. In view of this, the vibrationally averaged PV potential has been computed along all the normal modes of CHBrClF and the heavier CHAtFI analogue. Akin to the observation in CHBrClF molecule, the bending vibrational modes in CHAtFI exhibit larger deviation in the PV energies as compared to other normal modes and hence, larger vibrational splittings and an enhancement of the perturbative 2D effects can be expected due to these modes. Now, the vibrationally averaged HF and LDA PV gradients are utilised for calculating the 1D (Eq. 9) and 2D (Eq. 8) perturbative treatment for the vibrational energy levels up to 3 rd vibrational states. For this purpose, in addition to the PV energy gradient ∇E PV , which is available along all the normal modes, the second derivative (a 1 in Table II which  necessary to compute a full profile, because with the availability of the analytic gradient, the second derivative (a 1 ) can be computed from the linear fit of only a few energy gradient points close to the equilibrium structure (q = −0.125 and q = 0.125). This term is derived from a linear fit from these two points with another at q = 0 (cf . Tables S19-S20 in the Supporting Information).
The multi-mode contributions for the fundamental vibrational transition (from n = 0 to n = 1) for all the vibrational modes are estimated from Eq. 8. The q 8 , q 7 , q 3 and q 2 bending normal coordinates (although the corresponding fundamentals give rise to the least intense peaks in the infrared spectrum) contribute significantly to the multi-mode effects in the various fundamental transitions. Figure 7 displays the contributions towards the multi-mode effects for each of the vibrational modes of CHBrClF and CHAtFI molecules at HF and LDA level of theory. Contributions from ν 8 (in red) and ν 7 (in orange) are the largest and they contribute to almost every fundamental, whereas q 6 and q 3 contribute moderately in CHBrClF. The normal coordinate of the mode with the most intense fundamental in the vibrational spectrum, the C-F stretching fundamental (ν 4 ), and other normal coordinates hardly contribute to multi-mode effects in other fundamentals. The observations are slightly different in CHAtFI. The contributions from q 8 (At-C-F bending in red) and q 7 (F-C-I bending in orange) normal coordinates are significant. In addition to this, the C-At stretching vibration (q 6 ) has shown a substantial contribution to the multi-mode effects. This is expected, as astatine is the heaviest nucleus of the molecule and, due to the steep Z-scaling of parity-violating effects, with Z being the nuclear charge, should also have a pronounced influence on the parity violating shifts. The C-F stretching fundamental (ν 4 ) gets a large contributions from the H-C-F bending (q 3 in olive) normal coordinate, most probably due to the strong coupling between themselves since their frequencies are very close. Otherwise, the contributions from other vibrations are not that significant.
The total observed shifts associated with the fundamental transition (from n = 0 to n = 1) for all the vibrational modes of the (S)-enantiomer of these two chiral derivatives from the 1D and 2D perturbative treatment are displayed in Figure 8. It should be noted that one needs to multiply these values by 2 and then divide by the vibrational frequencies of the corresponding fundamentals to get the relative PV vibrational frequency splittings between the C-F stretching fundamental of the S-and R-enantiomers of CHBrClF and CHAtFI. The corresponding PV shifts up to the 3 rd vibrational energy level as well as the PV shifts in the fundamental transitions (from n = 0 to n = 1) are included in Tables S25-S26 of the Supporting Information. For CHBrClF, the vibrational frequency shifts from the 1D perturbative treatment evaluated at HF level are larger in magnitude for ν 8 , ν 6 , ν 5 and ν 2 than that of the ν 4 C-F stretching mode. The steeper slope of the curves due to the large variation in the PV energies (which increases the numerical value of ∇E PV ) along q 8 , q 7 , q 3 and q 2 in Figure 5 favours these large vibrational shifts in the bending modes. On the other hand, the 2D perturbation results in large shifts for most modes as compared to the C-F stretching mode (ν 4 ). This is expected since the parity-violating multi-mode effects are significant in other bending modes and stretching modes other than the C-F stretching vibration (See Figure 7). From the comparison of 1D and 2D effects in frequency shifts within individual modes, it is found that the 2D effects increase the frequency shifts only in ν 9 , ν 7 , ν 4 and ν 3 , whereas other fundamentals follow the opposite trend (cf. Figure 8) where a cancellation of the 1D contribution is caused due to the 2D treatment. Total calculated PV shifts in the fundamentals transitions from both of the perturbative 1D and 2D treatments align either towards positive or negative directions. No alternation of the sign is observed even after adding multi-mode contributions to 1D perturbative ones. Overall, the estimated vibrational shift in all the fundamental transitions matches to the same order of magnitude of previous theoretical values. 34,36 The case for CHAtFI is slightly complicated. The 2D terms are important for all modes where ν 8 , ν 7 and ν 5 fundamentals show frequency shifts that are larger in absolute value as compared to the C-F stretching (ν 4 ) mode (cf. Figure 8). Multimode effects increase the frequency shifts by a factor of ∼ 1.5 in ν 7 when compared to the perturbative 1D corrections. The correction due to multi-mode effects for the most intense C-F stretching mode reduces the absolute value as compared to the perturbative 1D result, i.e. the perturbative 2D treatment decreases the energy gap between the ground (n = 0) and 1 st (n = 1) vibrational state. Even, as mentioned above (see Table  IV), at LDA level the sign is altered (cf. Figure 8). A substantial increase in the ν 2 fundamental at LDA level is seen after including the multi-mode effects which is about 215% of the perturbative 1D value (See Figure 8). Although the fundamentals ν 8 and ν 5 have relatively larger frequency shifts for the vibrational transition from n = 0 to n = 1, the values decrease when multi-mode effects are introduced. Multimode contributions in ν 7 increases the frequency shifts in magnitude as compared to 1D effects. It may be noted that the asso-ciated transition frequencies of these fundamentals are not in the range of the CO 2 laser used in the previous 42 experimental set up. But, in spite of their lower intensity in the infrared spectrum, this might not be a severe limitation due to new developments in contemporary laser technology with quantum cascade laser being available in a broader frequency range.
We should emphasise, finally, that we reused in our present work the equilibrium structures and harmonic vibrational force fields from a previous study 40 to allow for direct comparison of results, Whereas a sophisticated CCSD(T) approach was used in Ref. 40, the basis sets were only of limited quality, which impacts not only on the equilibrium structures and harmonic force constants obtained, but also on the quality of molecular properties as was noted in our recent study 89 of nuclear electric quadrupole coupling constants in CHBrClF and CHClFI. For CHBrClF and its deuterated isotopomer, improved anharmonic ab initio force fields have been reported for instance in Ref. 90. In particular for the compound with the radioactive halogen, CHAtFI, for which no experimental information is available and for which very pronounced multi-mode effects on the C-F stretching fundamental were computed herein, a study with an improved description of the conventional parity-conserving effects is indicated. This becomes even more important by virtue of the potential role of heavy elemental chiral molecules in the search for dark matter candidates 87,88 and recent advances in laser spectroscopy of radioactive molecules with short-lived nuclei. 91

V. CONCLUDING REMARKS
In the present article, the gradient of the molecular parityviolating (PV) potential has been derived and implemented within a quasirelativistic framework at the level of HF and LDA, and is applied to rotational and vibrational spectroscopy of polyatomic molecules. A systematic study of frequency shifts in the rotational and vibrational spectra of chiral polyhalomethanes has been carried out.
A one-dimensional perturbative treatment has been compared to a previously calculated full anharmonic onedimensional approach and found to be in good agreement. For the hypothetical astatine compound, the predicted frequency splitting is in the right order of magnitude to be measured when treating the problem in the one-dimensional approximation. The results form the multi-mode effects in CHBrClF and CHAtFI reflects the importance of nonseparable anharmonic effects and suggests that the C-F stretching mode is not ideally suited for a measurement of electroweak PV effects, since other bending and stretching modes contribute significantly and can not be neglected. In view of this, choosing a proper vibrational transition is equally important as choosing a suitable molecule for conducting PV experiments and calculating the PV effects on the vibrational transition of chiral derivatives. Multimode effects are crucial for getting more insights about the vibrational transitions as well as rotational spectra in chiral compounds containing heavier atoms.
The present analytic derivative approach within a quasirelativistic mean-field scheme leads to an important simpli-fication for routine calculation of PV frequency shift in rotational spectra of chiral molecules and of multi-mode contributions to the PV frequency shifts in rovibrational spectroscopy. The implementation of the analytic ∇E PV at other DFT levels such as GGA (e.g. BLYP) and GGA-based hybrid-functionals (e.g. B3LYP) is currently underway, which will increase its applicability to diverse chiral systems containing heavier nuclei, in particular transition metals. Also, such a development helps in measuring the functional dependencies and influence of electron correlation effects for the theoretical estimation of the PV effects in the rotational and vibrational transitions of chiral compounds. This will provide valuable information for future experiments aiming at the detection of molecular parity violation.

ACKNOWLEDGMENTS
The authors are particularly thankful to Yunlong Xiao, Sophie Nahrwold, Timur Isaev and Christoph van Wüllen for stimulating discussions. Financial support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Projektnummer 328961117 -SFB 1319 ELCH and the VolkswagenFoundation is gratefully acknowledged. The center for scientific computing (CSC) Frankfurt is thanked for computer time.

Appendix A: Variational Perturbation Theory
Following Refs. 56-58 the total energy E( R 0 ) = Ψ Ĥ ( R 0 ) Ψ parametrically depends on the nuclear coordinates. Within a mean-field approach, the energy satisfies where t is the vector of orbital rotations and ∇ t is the gradient with respect to t. The variational energy of slightly displaced nuclei E( R 0 + η) can be expanded in a Taylor series in the nuclear displacements η as (see e.g. Ref. 59): where ∇ η is the gradient with respect to the nuclear displacement vector. Thus, the total variational energy is a function of λ PV as well: In combination with Eq. (A2) the parity violating contribution E PV ( R 0 + η, λ PV ) = E( R 0 + η, λ PV )−E( R 0 ) to the gradient of the energy with respect to nuclear displacements is: Appendix B: Spin polarised density functionals We employ the approximation of non-collinear DFT. Then, we can write a exchange-correlation potential operatorV Note, that this definition of the spin density coincides with what is usually called in non-collinear DFT the magnetization density and should not be confused with the length of the magnetization density which is what is called spin density in non-collinear DFT. We can write the exchange-correlation energy as and the gradient with respect to nuclear displacements is In presence of a perturbing operator such as the parityviolating potential we have to compute the first order perturbed exchange-correlation potential operator aŝ We can express terms that are proportional to gradients of the density matrix in terms of the Fock matrix grad of contracted gradients of two-electron integrals with elements: The canonical SCF equations are where S = σ 0 ⊗ S 1c is the overlap matrix constructed from the one component overlap matrix S 1c of the basis functions and . The matrix of coefficients is chosen as Exploiting the Hermiticity of F and ε we find from this Furthermore, the gradient of the orthonormality condition in Eq. C3 gives Therewith, the contribution to the energy gradient from the Fock matrix in Eq. (C1) can be expressed via the orthonormality condition as where we introduced the energy weighted density matrix (EWDM) W as Appendix D: Linear response equations

First order perturbed density matrix
The infinite order coefficients can be expressed in terms of orbital rotations, expressed via the anti-Hermitian matrix T ∞ as were C 0 is the initial guess of orthonormal orbitals and (C ∞ ) † = e −T ∞ C † 0 . From the definition of the density matrix [Eq. (18) where N is a diagonal matrix containing the occupation numbers n i . Thus, the perturbed density matrix of infinite order D ∞ can be written in terms of orbital rotation matrix T ∞ as Moreover, we can write D ∞ in terms of the unperturbed orbital coefficients C, received from orbital rotations t 0 in absence of the perturbation with the corresponding rotation matrix T 0 and an the rotation matrix T as: When assuming N to be successively sorted for occupied and unoccupied orbitals we can write the transformation matrix T in blocked form as where, o means occupied and u means unoccupied. In first order the perturbed density matrix can be written in terms of the unoccupied-occupied block T uo of T: Thus, we arrive at (D6)
where, a prime denotes matrices linear in λ PV and unprimed matrices are unperturbed (λ PV = 0). The perturbed Fock matrix F is for perturbation independent basis functions as we have in the case of the perturbation due to the PV potential: where, we introduced the self consistent transformation matrix of the orbital coefficients T for which C = CT. From Eq. (D5) we see that only the unoccupied-occupied rotations change the orbitals in first order, such that we can focus on the of unoccupied-occupied block of Eq (D7). Thus, for perturbation independent basis functions the CPHF/CPKS equations reduce to the coupled equations with the unoccupied (u) and occupied (o) orbital sub-blocks of C and ε. From this, one arrives at the linear response equations (see e.g. Refs. 70 and 71) where the index b runs over unoccpuied orbitals and the index j runs over occupied orbitals. The Hermiticity factor h results from the orthonormality condition which for perturbation independent basis functions yields C † SCT = −T † C † SC ⇔ T † = −T, which is in accordance with Eq. (D2). For observables such as H PV we have h = +1.

Supporting Information for
Quasi-relativistic approach to analytical gradients of parity violating potentials Herein, we report analytically determined Cartesian gradient components of the parity violating potential as obtained on the Hartree-Fock (HF) and local density approximation (LDA) level at the equilibrium structures of the (S)-enantiomer of the various methane derivatives (see Tables S1-S4). They refer to the specific equilibrium structures and orientation as reported previously in the Supplementary Material to R. Berger and J. L. Stuber, Mol. Phys. 105, 41 (2007). Based on the calculations and verifications on the convergence of parity violating energies and parity violating energy gradients of few planar molecules 1 , we here give the values up to four digits after the decimal point. The parity violating energy at the equilibrium structures as obtained on the HF and LDA levels are given in table S5. Molecular orbitals (MOs) have been converged until the change of the SCF energy and relative change of spin-orbit energy between two successive iterations dropped below at least 10 −6 E h and 10 −12 respectively. MOs of CHClFI at the HF level, however, have been converged until the relative change of the spin-orbit energy was on the level of about 10 −15 as the corresponding PV energy for the equilibrium structure did not achieve the desired accuracy with the 10 −12 criterion. In practice, the spin-orbit energy criterion was by far the more restrictive one. At the end of the process, the change in SCF energy typically remained below 10 −9 E h .
In table S6, we repeat unscaled harmonic vibrational wavenumbers of all normal modes as computed in R. Berger and J. L. Stuber, Mol. Phys. 105, 41 (2007), but give for precision reasons one more digit. The Cartesian displacement coordinates of each of these modes, corresponding to a unit shift along the dimensionless reduced normal coordinates q i , are given in tables S8-S11. The displacements refer to the same equilibrium structure and orientation as given in the Supplementary Material to R. Berger and J. L. Stuber, Mol. Phys. 105, 41 (2007), where the displacements for q 4 were given. The various displacements are provided herein to clarify also their phase used in determining cuts through the parity violating potentials along the normal modes. Table S7 depicts the unscaled harmonic vibrational wavenumbers as computed from the anharmonic force field calculations at the CCSD(T) level of theory for the equilibrium structure of the (S)-enantiomer of CHBrClF, CHClFI, CHBrFI and CHAtFI molecules using the MOLPRO, which are utilized for calculating vibrationally averaged HF and LDA parity violating potential for energy levels.
In tables S12-S15, we give the semi-diagonal cubic force constants as computed with the SURF module of the program package Molpro. Signs have been adjusted where neccessary to reflect the same phases of the normal modes as reported in tables S8-S11.
In table S16, HF and LDA level parity violating energies along the dimensionless reduced normal coordinates q Tables S17-S20 depict the HF and LDA level PV energy gradients along the dimensionless reduced normal coordinates (q 4 in the range from −3.00 to +3.00) corresponding to the C-F stretching mode (ν 4 ) for the (S)-enantiomers of CHBrClF, CHClFI, CHBrFI and CHAtFI respectively. On the other hand, the HF and LDA level PV energy gradients along the dimensionless reduced normal coordinates q = −0.125, +0.00 and +0.125 corresponding to all (except C-F stretching mode ν 4 ) the normal modes of the (S)-enantiomer of CHBrClF and CHAtFI are given respectively in tables S21 and S22. The values are reported up to four digits after the decimal point.
In table S23, the fitting coefficients of the LDA PV energy and the LDA PV energy gradients along the C-F stretching mode (ν 4 ) of the methane derivatives are given. Whereas in table S24, the second order derivative terms of the HF and PV LDA energy are given. These are obtained from a linear fit of the PV energy gradients at q = −0.125, +0.00 and +0.125 along the C-F stretching mode.
In tables S25 and S26, the vibrationally averaged HF and LDA parity violating potential E S v,PV up to four lowest vibrational energy levels as obtianed from the 1-dimensional and 2-dimensional perturbative treatment are given for all (expect C-F stretching vibration ν 4 ) the normal modes of (S)-enantiomers of CHBrClF and CHAtFI molecules. The frequency shift for the fundamental transition for all the vibrational modes are calculated. The percentage of the multimode effects in each of the energy levels are also included.          Fig. 6 in the manuscript. S12 Table S17: HF and LDA level PV energy gradients (in 10 −13 cm −1 ) along the dimensionless reduced normal coordinates (q 4 in the range from −3.00 to +3.00) corresponding to the C-F stretching mode (ν 4 ) for the (S)-enantiomer of CHBrClF.  Fig. 1 in the manuscript. d See the red curve of Fig. 1 in the manuscript. Table S18: HF and LDA level PV energy gradients (in 10 −12 cm −1 ) along the dimensionless reduced normal coordinates (q 4 in the range from −3.00 to +3.00) corresponding to the C-F stretching mode (ν 4 ) for the (S)-enantiomer of CHClFI.  Fig. 2 in the manuscript. f See the red curve of Fig. 2 in the manuscript. S13 Table S19: HF and LDA level PV energy gradients (in 10 −12 cm −1 ) along the dimensionless reduced normal coordinates (q 4 in the range from −3.00 to +3.00) corresponding to the C-F stretching mode (ν 4 ) for the (S)-enantiomer of CHBrFI.       Table S25: Vibrationally averaged HF and LDA parity violating potential E S v,PV for energy levels n in 10 −12 cm −1 for all (expect C-F stretching vibration ν 4 ) the normal modes of (S)-enantiomer of CHBrClF.