Quasi-relativistic approach to analytical gradients of parity violating potentials
Sascha A Brück, Nityananda Sahu, Konstantin Gaul, Robert Berger
QQuasi-relativistic approach to analytical gradients of parity violatingpotentials a) Sascha A Brück, Nityananda Sahu, Konstantin Gaul, and Robert Berger
1, 2, 3 Frankfurt Institute for Advanced Studies, Ruth-Moufang-Straße 1, 60438 Frankfurt am Main,Germany Fachbereich Chemie, Philipps–Universität Marburg, Hans-Meerwein-Straße 4, 35032 Marburg,Germany Clemens-Schöpf-Institut, Technische Universität Darmstadt, Alarich-Weiss-Straße 4, 64287, Darmstadt,Germany (Dated: 22 February 2021)
An analytic gradient approach for the computation of derivatives of parity-violating (PV) potentials with respect todisplacements of the nuclei in chiral molecules is described and implemented within a quasirelativistic mean-fieldframework. Calculated PV potential gradients are utilised for estimating PV frequency splittings between enantiomersin 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 reportedtheoretical values. The influence of non-separable anharmonic effects (multi-mode effects) on the vibrational frequencyshifts, which are readily accessible with the present analytic derivative approach, are estimated for the C–F stretchingfundamental of all four molecules and computed for each of the fundamentals in CHBrClF and CHAtFI. Multi-modeeffects are found to be significant, in particular for the C–F stretching modes, being for some modes and cases of similarsize as the single-mode contribution.
I. INTRODUCTION
Early after the proof of parity violation (PV) was providedin 1957 by the famous experiment of β − disintegration of Conuclei , which realised a proposal by the theoreticians Leeand Yang , it was suggested by Yamagata that PV weak in-teractions can induce a tiny energy difference between a chiralmolecule and its non-identical mirror-image. Diverse ex-perimental schemes have been proposed to detect PV effectsin chiral molecules, ranging from Mössbauer spectroscopy to vibrational spectroscopy, electron paramagnetic res-onance spectroscopy, rotational spectroscopy andnuclear magnetic resonance spectroscopy to time-dependent approaches and quantum-beat experiments .We refer the reader for more details and an extended overviewto a collections of reviews on this subject . The mostaccurate experimental attempts reported so far were in thehigh-resolution infrared spectroscopy of bromochlorofluo-romethane (CHBrClF). An upper bound of ∆ ν PV / ν ≈ − was obtained for the relative PV difference of the C-Fstretching frequency between ( R )- and ( S )-enantiomers. Thetheoretically predicted value for this relative frequencysplitting is of the order of 10 − , however. Later, with animproved experimental set-up, a measurement of the samecompound with a resolution of 5 × − has been reported in2002 and with a new set-up, it may be hoped that a precisionof 10 − can be reached. As nuclear spin-independentelectroweak PV effects in chiral compounds scale approxi-mately with nuclear charge Z to the power of five (in the pres-ence of spin-orbit coupling), compounds containing heavymetal nuclei could be of greater experimental value than the a) Parts of this work were reported in preliminary form in S. Brück, Diplomathesis, University of Frankfurt, 2011. originally used organic molecules. Theoretical searches in thisdirection have been already made on molecules containing forinstance bismuth, rhenium, mercury and astatine.
In general, a measurement of the PV energy as a differencebetween electronic energies of separated enantiomers wouldbe very difficult since it comes on top of the rest mass energyof the molecule, but the scheme proposed by Quack , forinstance, circumvents this by directly measuring the PV en-ergy via the PV induced time-dependent interconversion be-tween states of opposite parity . An alternative for detect-ing molecular parity violation is to measure the difference be-tween PV energy differences arising due to the PV potential[ E PV ( (cid:126) q ) ] in the vibrational and rotational transitions in chiralcompounds, with (cid:126) q denoting the vector of dimensionless re-duced normal coordinates. Often also V PV is used as symbolfor the PV potential. In first order perturbation theory, E PV ( (cid:126) q ) gives rise to a PV shift in the n th vibrational energy level of agiven enantiomer ( R or S ) according to E R , Sn , PV ≈ (cid:104) Ψ R , Sn | E PV ( (cid:126) q ) | Ψ R , Sn (cid:105) , (1)where | Ψ R , Sn (cid:105) denote the n th vibrational state of the R - and S -enantiomer, which are obtained by solving the parity-conserving (PC) rovibrational Schrödinger equation for eachenantiomer. Since, E PV ( (cid:126) q ) is parity odd, E Rn , PV and E Sn , PV val-ues are numerically of equal magnitude but have an oppositesign. Thus, the corresponding PV energy difference betweenthe n th vibrational levels of the two enantiomers is ∆ E vib , PV = E Sn , PV − E Rn , PV ≈ (cid:104) Ψ Sn | E PV ( (cid:126) q ) | Ψ Sn (cid:105) . (2)The relative change in vibrational ( ∆ E vib , PV ) and rotational( ∆ E rot , PV ) energies between right- and left-handed moleculesis expected to scale to same order of magnitude as that of therelative change in electronic ( ∆ E el , PV ) energy. ∆ E el , PV E el ≈ ∆ E vib , PV E vib ≈ ∆ E rot , PV E rot (3) a r X i v : . [ phy s i c s . c h e m - ph ] F e b The introduction of PV within electroweak quantum chem-istry not only affects the energy of an enantiomer, but also itsequilibrium structure if the PC vibrational potential gets modi-fied by the PV potential. In general, the PV potential inducesa minute change in the equilibrium structure of a moleculecompared to the equilibrium structure without PV contribu-tion, if the PV energy gradient ( (cid:126) ∇ E PV ) does not vanish there.Thus, the weak interaction leads to a change of the equilib-rium structure of a chiral molecule, which is different for bothenantiomers due to PV effects.This effect could in principle be measured by microwavespectroscopy, since it leads to a shift of the rotationalconstants. In practice, the change of the structure due to (cid:126) ∇ E PV at the minimum of the parity conserving potentialcan be estimated with the help of the vibrational Hessian F MW , of the PC potential, given in mass-weighted Carte-sian displacement coordinates. From this 3 N nuc × N nuc -dimensional vibrational Hessian, with N nuc being the numberof nuclei in the system, contributions of infinitesimal trans-lational and rotational displacements can be eliminated withthe 3 N nuc × ( N nuc − ) dimensional projection matrix A byforming ( A T F MW A ) − , so that the corresponding Hessian ininternal Cartesian displacement coordinates results, which canbe inverted. When M is a ( N nuc × N nuc ) diagonal matrixcontaining the masses of the various atoms, the final displace-ments ( δ PV (cid:126) R ) of the atoms from their equilibrium position dueto (cid:126) ∇ E PV read δ PV (cid:126) R = − M − / A ( A T F MW A ) − A T M − / (cid:126) ∇ E PV . (4)Coordinates and displacements are transformed subsequentlyto the principal axis system of the PC equilibrium structureand the splittings of the diagonal elements ( ∆ I ) of the momentof inertia tensor due to PV effects are approximated to first or-der in δ PV (cid:126) R . For the xx component of ∆ I we have for instance: ∆ I x = m A ( y A δ PV y A + z A δ PV z A ) (5)and analogous for the other diagonal elements. Then, thechanges in diagonal elements are used to estimate the split-tings of rotational constants ( ∆ X R ) between two enantiomerswithin the approximate expression ∆ X R X R ≈ − ∆ I X I X (6)with X R being one of the rotational constants ( A , B or C )and I X being the corresponding eigenvalues of the momentof inertia tensor. From Eq. 4 the importance of the gradi-ent of the parity violating potential with respect to displace-ments of the nuclei becomes particularly evident. The numer-ical calculation of this term, however, is tedious, in particu-lar within a (quasi)relativistic electronic structure framework,which is why we present in this work an approach for calculat-ing this gradient analytically within a quasirelativistic mean-field framework.Another promising experiment for detecting parity viola-tion in chiral molecules is vibrational spectroscopy. Most ofthe effort has been put into the measurement of vibrational frequency shifts due to PV effects. Therelative vibrational frequency shift for a transition from state n to m is determined by taking the difference of the vibrationallyaveraged PV potentials of the enantiomers and dividing by thecorresponding transition energy h ν mn ∆ ν mn , PV ν mn = ( (cid:104) Ψ m | E PV ( (cid:126) q ) | Ψ m (cid:105) − (cid:104) Ψ n | E PV ( (cid:126) q ) | Ψ n (cid:105) ) h ν mn (7)where the factor two enters into the equation since the dif-ference of the R and S enantiomer is twice ( cf. Eq. 2 ) thedifference between one enantiomer and the parity conserv-ing case, where no shift occurs. The vibrationally averagedpotentials depend on the multi-dimensional PV energy sur-face and not only on a single point energy at the equilibriumstructure. In addition, the complete rovibrational wavefunc-tion would be needed to compute the expectation value. Asthis problem is extremely complex even for relatively smallmolecules, the PV effects of electronic, vibrational and rota-tional degrees of freedom are typically separated in a first step.Still the computational effort is in general too high, so as a sec-ond step, the multi-dimensional problem is usually split intoone-dimensional problems, where the movements along thenormal coordinates are treated separately as if they were inde-pendent of each other.
This can be augmented by addingcontributions from the potential depending on a smaller num-ber of modes order by order.
In practice, the PC poten-tial ( V BO ) and PV potential [ E PV ( (cid:126) q ) ] can be evaluated at one-dimensional cuts along the dimensionless reduced normal co-ordinates (cid:126) q . This approach has already been used to esti-mate the vibrational frequency splitting of the C-F stretchingmode in chiral halogenated methane derivatives.
A different approach comes from perturbation theory, where the PC and PV potentials are expanded in a Taylor se-ries. In this, the influence of multi-mode effects are estimatedby the calculations of the derivatives of the PV potentials withrespect to all normal coordinates. The contributions of lowestnon-vanishing order to the n th vibrational energy levels of anormal mode ν r within the perturbative treatment are (cid:104) Ψ n r | E PV ( (cid:126) q ) | Ψ n r (cid:105) ≈ E PV ( (cid:126) R )+ (cid:18) n r + (cid:19) (cid:18) ∂ E PV ∂ q r − ∑ s h ω s ∂ V BO ∂ q r ∂ q s ∂ E PV ∂ q s (cid:19) (8)where E PV ( (cid:126) R ) is the PV energy at the equilibrium structureand ω s denotes the harmonic vibrational angular frequencyof normal mode ν s . The terms ∂ E PV ∂ q s and ∂ E PV ∂ q r are the first(gradient) and second derivatives of the PV energy along nor-mal coordinates ( q s , q r ), respectively. Ignoring all the r (cid:54) = s terms in the second sum of Eq. 8 gives the one-dimensional(1D) perturbative estimate of the vibrational energy for nor-mal mode ν r as (cid:104) Ψ n r | E PV ( q r ) | Ψ n r (cid:105) = E PV ( (cid:126) R )+ (cid:18) n r + (cid:19) (cid:18) ∂ E PV ∂ q r − h ω r ∂ V BO ∂ q r ∂ E PV ∂ q r (cid:19) . (9)The drawback of this method is that the one-dimensionalterms of the PC and PV potentials are not included to infiniteorder. But the advantage is that, once the Cartesian gradientof the PV energy at the equilibrium structure is known, thetwo-dimensional coupling term can be included without toomuch effort, because only the semidiagonal cubic force con-stants are needed, whereas the gradient of PV energy alongall the normal mode can conveniently be determined by pro-jecting the Cartesian PV energy gradient onto displacementsalong to the normal coordinates. Hence, contributions fromother modes ( cf. Eq. 10) in addition to the single-mode termsdescribed by Eq. 9 can easily be accounted for by consideringfirst order multi-mode (MM) effects on the vibrational energylevels as shown in Eq. 8 (cid:104) Ψ n r | E PV ( (cid:126) q ) | Ψ n r (cid:105) MM = − (cid:18) n r + (cid:19) ∑ s (cid:54) = r h ω s ∂ V BO ∂ q r ∂ q s ∂ E PV ∂ q s . (10)We define with these MM terms the perturbative approxima-tion of the vibrational energy shifts with 2D coupling termsas (cid:104) Ψ n r | E PV ( (cid:126) q ) | Ψ n r (cid:105) = (cid:104) Ψ n r | E PV ( (cid:126) q ) | Ψ n r (cid:105) + (cid:104) Ψ n r | E PV ( (cid:126) q ) | Ψ n r (cid:105) MM . (11)Thus, for the calculation of the influence of PV on the vi-brational and rotational spectra, which is one of the maingoals of the current manuscript, the knowledge of the gra-dient of the PV potential along all the normal modes is es-sential. The present article describes an analytical approachfor calculating the gradient of the PV energy ( (cid:126) ∇ E PV ) at HFand LDA level of theory (see the following section II). Thesevalues obtained from this analytic gradient approach, withcorresponding computational details being described in sec-tion III, are then utilised in section IV for determining therelative shifts in rotational constants and vibrational frequen-cies, due to PV effects in chiral methane derivatives. Besidesthis, the influence of non-separable anharmonic effects (multi-mode effects) on PV induced vibrational frequency shifts forall the vibrational normal modes is calculated for CHBrClFand CHAtFI molecules. II. THEORY
For the description of the vibrational movement of nuclei inthe Born-Oppenheimer approximation, it is often sufficient toevaluate the electronic structure for a single fixed arrangementof the nuclei described by the coordinates (cid:126) R and to treat thedisplacement of nuclei perturbatively (see also Refs. 56–58).The electroweak parity-violating potential is very smallcompared to the PC potential due to the appearance of theFermi coupling constant which is 2 . × − E h a inatomic units. Therefore, we can treat the parity-violating elec-troweak Hamiltonian ˆ H PV as a small addition to the parityconserving molecular Hamiltonian ˆ H :ˆ H = ˆ H + λ PV ˆ H PV , (12)where we have introduced λ PV as a formal perturbation pa-rameter. In a mean-field approach the leading order parity violatingcontribution E PV ( (cid:126) R + (cid:126) η , λ PV ) = E ( (cid:126) R + (cid:126) η , λ PV ) − E ( (cid:126) R ) tothe gradient of the variational energy with respect to nucleardisplacements is: (cid:126) ∇ (cid:126) η E PV ( (cid:126) R + (cid:126) η , λ PV ) (cid:12)(cid:12)(cid:12) (cid:126) η = (cid:126) ≈ λ PV ∂ (cid:126) ∇ (cid:126) η E ( (cid:126) R + (cid:126) η , λ PV ) ∂ λ PV (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126) η = (cid:126) λ PV = (13)For details on the specific case of electroweak parity violation,see Appendix A, whereas for the general case of variationalperturbation theory, see e.g. Ref. 59.In this paper, we want to focus on the molecu-lar Hamiltonian approximated within a (quasi-relativistic)two-component zeroth-order regular approximation (ZORA)framework on the level of Generalised Hartree-Fock(GHF) or Generalised Kohn-Sham (GKS) density functionaltheory (DFT) (see Ref. 62 and 63):ˆ H = ˆ H ZORA = V GHF , GKS + ˆ V nuc ( (cid:126) r ) + c (cid:126) σ · ˆ (cid:126) p ω ( (cid:126) r ) (cid:126) σ · ˆ (cid:126) p (cid:124) (cid:123)(cid:122) (cid:125) ˆ h ZORA (14)Here, V GHF , GKS is the effective electron repulsion potentialwithin the GHF or GKS framework, ˆ (cid:126) p is the linear mo-mentum operator, (cid:126) σ is the vector of Pauli spin matrices and ω ( (cid:126) r ) = c − ˜ V ( (cid:126) r ) is the ZORA factor with the ZORA modelpotential ˜ V as proposed by van Wüllen to alleviate the gaugedependence of ZORA. The ZORA one-electron operator hasan electron spin-independent and an electron spin-dependentcontribution:ˆ h ZORA = ˆ V nuc + c ˆ (cid:126) p · (cid:16) ω ( (cid:126) r ) ˆ (cid:126) p (cid:17) + c ı (cid:126) σ · ˆ (cid:126) p × (cid:16) ω ( (cid:126) r ) ˆ (cid:126) p (cid:17) (15)Within ZORA, the nuclear-spin independent parity-violating electroweak one-electron Hamiltonian appears as ˆ h PV = G F √ N nuc ∑ A = (cid:110) ω ( (cid:126) r ) Q W , A ρ A ( (cid:126) r ) , c (cid:126) σ · ˆ (cid:126) p (cid:111) + , (16)where { A , B } + = AB + BA is the anti-commutator, Q W , A is theweak charge of nucleus A and ρ A is the normalised nucleardensity distribution. This operator has only electron spin-dependent contributions.We expand the ZORA two-component HF or KS molec-ular orbitals (MOs) φ i in a linear combination of real one-component basis functions χ µ and complex two-componentcoefficients (cid:126) C µ i = (cid:32) C ( α ) µ i C ( β ) µ i (cid:33) as φ i = ∑ µ (cid:126) C µ i χ µ (17)In this two-component framework, we can define four com-plex one-component density matrices ( κ = , , , D ( κ ) µν = N orb ∑ i = n i (cid:126) C † µ i σ κ (cid:126) C ν i (cid:124) (cid:123)(cid:122) (cid:125) D ( κ ) i µν , (18)where the 0 th component of the Pauli spin matrices is the 2 × D = ∑ κ = ( σ κ ) ∗ ⊗ D ( κ ) . (19)We can write the expectation value of the one electronZORA operator as h ZORA = ∑ µν R (cid:40) ∑ κ = D ( κ ) µν h ( κ ) ZORA , µν (cid:41) , (20) and the energy contribution of the parity-violating elec-troweak potential as h PV = ∑ µν R (cid:40) ∑ κ = D ( κ ) µν h ( κ ) PV , µν (cid:41) . (21)In the following, the effective potential V GHF , GKS will berepresented in the space of basis functions in terms of the ma-trix of contracted two-electron integrals G . In the general caseof hybrid DFT, G ( κ ) is constructed as G ( κ ) µν ( D ) = ∑ ρσ (cid:20) δ k D ( κ ) ρσ ( µν | ρσ ) − a X D ( κ ) ρσ ( µσ | ρν ) + a DFT (cid:68) χ µ (cid:12)(cid:12)(cid:12) V ( κ ) XC ( D ) (cid:12)(cid:12)(cid:12) χ ν (cid:69)(cid:21) , (22)where the Mulliken notation for two elec-tron integrals is employed: ( µν | ρσ ) = (cid:82)(cid:82) d r d r χ µ ( (cid:126) r ) χ ρ ( (cid:126) r ) | (cid:126) r − (cid:126) r | χ ν ( (cid:126) r ) χ σ ( (cid:126) r ) . In case ofpure DFT (non-hybrid) we have a X = a X = a DFT = (cid:68) χ µ (cid:12)(cid:12)(cid:12) V ( κ ) XC (cid:12)(cid:12)(cid:12) χ ν (cid:69) are always real. In this paper, we restrictthe discussion to the spin-unpolarised local density approx-imation (LDA), in which the exchange-correlation potentialhas the form V ( ) XC , LDA = δ F XC , LDA [ ρ e ( (cid:126) r ; D ( ) )] δ ρ e ( (cid:126) r ; D ( ) ) , (23) with F XC , LDA being the LDA density functional and the elec-tronic number density function being ρ e ( (cid:126) r ; D ( ) ) = ∑ µν R (cid:110) D ( ) µν (cid:111) χ µ ( (cid:126) r ) χ ν ( (cid:126) r ) . (24)For the general form of the exchange-correlation potential seeAppendix B.The expectation value of the electron repulsion potential is V GHF , GKS = R (cid:40) ∑ κ = D ( κ ) µν G ( κ ) µν ( D ) (cid:41) . (25)The total ZORA energy in presence of the full perturbationˆ H PV , which implies λ PV =
1, is given by E ∞ = h ZORA + V GHF , GKS + h PV = R (cid:40) ∑ µν (cid:34) ∑ κ = D ∞ , ( κ ) µν h ( κ ) ZORA , µν + ∑ κ = D ∞ , ( κ ) µν G ( κ ) µν ( D ) + ∑ κ = D ∞ , ( κ ) µν h ( κ ) PV , µν (cid:35)(cid:41) , (26)where, ∞ refers to density matrices obtained variationally inthe presence of the perturbation, which is analogous to an in-finite order perturbation theory treatment of ˆ H PV (for details see Eq. (A3) in the Appendix), provided that this converges.Thus, the gradient with respect to nuclear displacements canbe written as (cid:126) ∇ (cid:126) η E ∞ = R (cid:40) ∑ µν (cid:34) ∑ κ = (cid:16) (cid:126) ∇ (cid:126) η D ∞ , ( κ ) µν (cid:17) h ( κ ) ZORA , µν + D ∞ , ( κ ) µν (cid:16) (cid:126) ∇ (cid:126) η h ( κ ) ZORA , µν (cid:17) + ∑ κ = (cid:16) (cid:126) ∇ (cid:126) η D ∞ , ( κ ) µν (cid:17) G ( κ ) µν ( D ∞ ) + D ∞ , ( κ ) µν (cid:16) (cid:126) ∇ (cid:126) η G ( κ ) µν ( D ∞ ) (cid:17) + ∑ κ = (cid:16) (cid:126) ∇ (cid:126) η D ∞ , ( κ ) µν (cid:17) h ( κ ) PV , µν + D ∞ , ( κ ) µν (cid:16) (cid:126) ∇ (cid:126) η h ( κ ) PV , µν (cid:17)(cid:35)(cid:41) (27)For a first order property that does not depend on the basisfunctions such as the parity-violating potential, this expres- sion can be simplified by using the orthonormality conditionof the HF equations and with Eq. (13) we receive the leadingorder parity-violating energy gradient (see Appendix C): (cid:126) ∇ (cid:126) η E PV ≈ R (cid:40) ∑ µν (cid:34) ∑ κ = D (cid:48) ( κ ) µν (cid:16) (cid:126) ∇ (cid:126) η h ( κ ) ZORA , µν (cid:17) + ∑ κ = D ( κ ) µν (cid:126) G ( κ ) grad , µν ( D , D (cid:48) ) + ∑ κ = D ( κ ) µν (cid:16) (cid:126) ∇ (cid:126) η h ( κ ) PV , µν (cid:17) − W (cid:48) ( ) µν (cid:16) (cid:126) ∇ (cid:126) η S µν (cid:17)(cid:35)(cid:41) . (28)where we have introduced the energy weighted density matrix (EWDM) W as W ( κ ) = N orb ∑ i = n i ε i D ( κ ) i , (29)and the matrix (cid:126) G ( κ ) grad of contracted gradients of two-electron integrals with elements: (cid:126) G ( κ ) grad , µν ( D , D (cid:48) ) = ∑ ρσ (cid:20) δ k D (cid:48) ( κ ) ρσ (cid:126) ∇ (cid:126) η ( µν | ρσ ) − a X D (cid:48) ( κ ) ρσ (cid:126) ∇ (cid:126) η ( µσ | ρν ) + a DFT (cid:68) (cid:126) ∇ (cid:126) η χ µ (cid:12)(cid:12)(cid:12) ˆ V (cid:48) ( κ ) XC ( D , D (cid:48) ) (cid:12)(cid:12)(cid:12) χ ν (cid:69)(cid:21) . (30)Note, that in the DFT case, derivatives of the exchange-correlation potential have to be computed for the calculationof the perturbed two-electron gradients, which in the spin-unpolarised LDA case are V (cid:48) XC , LDA ( D ( ) , D (cid:48) ( ) ) = δ V XC , LDA ( D ( ) ) δ ρ ( (cid:126) r ; D ( ) ) ρ ( (cid:126) r ; D (cid:48) ( ) ) . (31)For other functionals see Appendix B. D is the SCF den-sity matrix of the unperturbed system, i.e. received with theHamiltonian ˆ H only, and D (cid:48) and W (cid:48) are the perturbed densitymatrix and EWDM of first order in λ PV .The nuclear displacement gradients of the full unperturbedZORA Hamiltonian, i.e. the gradients of one-electron ZORAintegrals (cid:126) ∇ (cid:126) η h ( κ ) ZORA , µν , as well as gradients of two-electronintegrals and the exchange-correlation potentials needed for (cid:126) G ( κ ) grad , µν ( D ) have been implemented in Ref. 65. What remainsto be calculated are the perturbed density matrices D (cid:48) ( κ ) , theperturbed spin-independent EWDM W (cid:48) ( ) and the gradient ofthe PV integrals (cid:126) ∇ (cid:126) η h ( κ ) PV , µν . We describe the scheme for com-puting the latter in the following subsection II A, before wediscuss in subsection II B the linear response scheme that isused for computation of perturbed density matrices. A. Gradient of the PV integrals
The matrix elements of the one electron PV operator in ba-sis set representation is h ( κ ) PV , µν = (cid:68) χ µ , B (cid:12)(cid:12)(cid:12) ˆ h ( κ ) PV (cid:12)(cid:12)(cid:12) χ ν , C (cid:69) = ι ¯ h G F √ N nuc ∑ A = (cid:68) χ µ , B (cid:12)(cid:12)(cid:12) ˆ h ( κ ) PV , A (cid:12)(cid:12)(cid:12) χ ν , C (cid:69) = ι ¯ h G F √ N nuc ∑ A = (cid:0)(cid:10) χ µ , B (cid:12)(cid:12) ω ( (cid:126) r ) Q W , A ρ nuc , A ( (cid:126) r ) (cid:12)(cid:12) ∂ κ χ ν , C (cid:11) − (cid:10) ∂ κ χ µ , B (cid:12)(cid:12) ω ( (cid:126) r ) Q W , A ρ nuc , A ( (cid:126) r ) (cid:12)(cid:12) χ ν , C (cid:11)(cid:1) (32)with ∂ κ = ∂∂ x κ denoting a component of the electronic fourderivative, with the four vector being x ( , , , ) = ( 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 ba-sis functions depend on the nuclear coordinates, the geometrygradient of the PV integrals is (cid:126) ∇ (cid:126) η h ( κ ) PV , A , µν = (cid:68) (cid:126) ∇ (cid:126) η χ µ , B (cid:12)(cid:12)(cid:12) ˆ h ( κ ) PV , A (cid:12)(cid:12)(cid:12) χ ν , C (cid:69) + (cid:68) χ µ , B (cid:12)(cid:12)(cid:12) (cid:126) ∇ (cid:126) η ˆ h ( κ ) PV , A (cid:12)(cid:12)(cid:12) χ ν , C (cid:69) + (cid:68) χ µ , B (cid:12)(cid:12)(cid:12) ˆ h ( κ ) PV , A (cid:12)(cid:12)(cid:12) (cid:126) ∇ (cid:126) η χ ν , C (cid:69) (33)where, the second term is the Hellmann-Feynman term andthe first and third term are the basis function contributions.The integrals of the PV operator are purely imaginary and dueto Hermiticity of the PV operator, the basis function contribu-tions are connected by (cid:68) (cid:126) ∇ (cid:126) η χ µ , B (cid:12)(cid:12)(cid:12) ˆ h ( κ ) PV , A (cid:12)(cid:12)(cid:12) χ ν , C (cid:69) = − (cid:68) χ ν , C (cid:12)(cid:12)(cid:12) ˆ h ( κ ) PV , A (cid:12)(cid:12)(cid:12) (cid:126) ∇ (cid:126) η χ µ , B (cid:69) . (34)We split the nuclear displacements (cid:126) η into the separate con-tributions from different nuclei A . The basis functions de-pend on (cid:126) r − (cid:126) R A . Therefore, only the subset of basis func-tions centered at A contributes to the corresponding inte-grals. Furthermore, as we deal with Gaussian basis functions,derivatives with respect to nuclear coordinates can be repre- sented by derivatives with respect to electronic coordinates as ∂ A , κ χ B = − δ AB ∂ κ χ B , where ∂ A , κ denotes a component of thefour derivative of nucleus A . For the basis function contribu-tion to the PV gradient integrals, we arrive at (cid:68) χ µ , B (cid:12)(cid:12)(cid:12) ˆ h ( κ ) PV , A (cid:12)(cid:12)(cid:12) (cid:126) ∇ D χ ν , C (cid:69) = ı ¯ h G F √ δ DC N nuc ∑ A = (cid:16) − (cid:68) χ µ , B (cid:12)(cid:12)(cid:12) ω ( (cid:126) r ) Q W , A ρ nuc , A ( (cid:126) r ) (cid:12)(cid:12)(cid:12) (cid:126) ∇ ∂ κ χ ν , C (cid:69) + (cid:68) ∂ κ χ µ , B (cid:12)(cid:12)(cid:12) ω ( (cid:126) r ) Q W , A ρ A ( (cid:126) r ) (cid:12)(cid:12)(cid:12) (cid:126) ∇ χ ν , C (cid:69)(cid:17) . (35)For a Gaussian shaped nuclear density distribution, the Hellmann-Feynman term reads (cid:68) χ µ , B (cid:12)(cid:12)(cid:12) (cid:126) ∇ D ˆ h ( κ ) PV (cid:12)(cid:12)(cid:12) χ ν , C (cid:69) = ı ¯ h G F √ (cid:32) − (cid:42) χ µ , B (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ( (cid:126) ∇ D ω ( (cid:126) r )) ∑ A Q W , A ρ A ( (cid:126) r ) + ω ( (cid:126) r ) Q W , D ( (cid:126) ∇ ρ D ( (cid:126) r )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂ κ χ ν , C (cid:43) + (cid:42) ∂ κ χ µ , B (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ( − (cid:126) ∇ D ω ( (cid:126) r )) ∑ A Q W , A ρ A ( (cid:126) r ) + ω ( (cid:126) r ) Q W , D ( (cid:126) ∇ ρ D ( (cid:126) r )) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) χ ν , C (cid:43)(cid:33) , (36)with the gradient of a Gaussian normalised nuclear densitydistribution being (cid:126) ∇ ρ A ( (cid:126) r ) = − ζ A (cid:18) ζ A π (cid:19) / ( (cid:126) r − (cid:126) R A ) exp (cid:26) − ζ A (cid:12)(cid:12)(cid:12) (cid:126) r − (cid:126) R A (cid:12)(cid:12)(cid:12) (cid:27) . (37)The geometry gradient of the ZORA factor (cid:126) ∇ A ω is discussedin Ref. 62. The construction of ZORA, PV operators and itsderivatives is done by means of numerical integration.The above integrals are evaluated numerically on a grid. Inthe present work, an atom centered grids is used employingthe Treuter-Ahlrichs version of a Becke grid, where thepartitioning is done by weight functions w i for each grid point i . As these weight functions depend on the nuclear coordi-nates they give an additional contribution to the PV gradientwhich is calculated as N grid ∑ i (cid:16) χ µ , B ( (cid:126) r i ) ˆ h ( κ ) PV , A ( (cid:126) r i ) χ ν , C ( (cid:126) r i ) (cid:17) (cid:126) ∇ (cid:126) η w i ( (cid:126) R ) (38)The grid points move when the nuclei are slightly displaced.To account for this effect, we employ that translational in-variance of the molecule holds and, therewith, the net forceon the molecule has to be zero. We, therefore, subtract thenet force that results from the numerical integration proce-dure from the numerically integrated gradient contribution ( (cid:126) ∇ D h ( κ ) PV , A , µν ) num : (cid:126) ∇ D h ( κ ) PV , A , µν = ( (cid:126) ∇ D h ( κ ) PV , A , µν ) num − ∑ B ( (cid:126) ∇ B h ( κ ) PV , A , µν ) num . (39) B. Linear response computation of perturbed densitymatrices
In first order the perturbed density matrix can be writtenin terms of the unoccupied-occupied block T uo of an anti- Hermitian transformation matrix T (see Appendix D 1 for de-tails): D (cid:48) ( κ ) µν ( T uo ) = occ ∑ i unocc ∑ a (cid:104) (cid:126) C † µ a σ ( κ ) (cid:126) C ν i T ∗ ai + (cid:126) C † ν i σ ( κ ) (cid:126) C µ a T ai (cid:105) . (40)For how to compute perturbed density matrices to arbitraryorder see Ref. 68.We introduce the MO transformed PV operator as H PV , MO , i j = ∑ κ = ∑ µν (cid:126) C † µ i σ κ (cid:126) C ν j h ( κ ) PV , µν . (41)If we would assume that the first order perturbed wave func-tion would be calculated by transformation of the orbital co-efficients as (cid:126) C (cid:48) µ i ≈ unocc ∑ a (cid:126) C µ a ( ε a − ε i ) − H PV , MO , ai (42)then the perturbed density matrix would simply calculated as D (cid:48) ( κ ) µν ( ∆ − ◦ H uoPV , MO ) with ◦ being the Hadamard product andthe matrices ∆ , ∆ − having the elements ∆ ai = ( ε a − ε i ) and ∆ − ai = ( ε a − ε i ) − , respectively. Here a and i are defined inthe space of unoccupied and occupied orbitals, respectively.However, this does not account for the response of the or-bitals to the perturbation but correspond to a simple sum overstates approach in which it is assumed that the electronic Hes-sian is diagonal and therefore the perturbed SCF equations areuncoupled. Within HF and KS, however, the electronic Hes-sian is not diagonal as the two-electron matrix G is a functionof the orbitals and one has to solve the coupled perturbed HF(CPHF) or coupled perturbed KS (CPKS) equations (see Ap-pendix D 2 and for a detailed derivation e.g. Ref.69). To solvethe response equations, we use the reduced form (see Refs. 70and 71 and Appendix D 2 for details): ∑ b j A b j T b j + B b j T ∗ b j = − H uoPV , MO (43)where indices b run over unoccupied orbitals and indices j over occupied orbitals. The elements of the electronic Hessianare A ai , b j = ( ε a − ε i ) δ ab δ i j + ˜ G ai , jb (44) B ai , b j = ˜ G ai , b j , (45)with ˜ G ai , b j being an element of the four-index MO trans-formed two-electron tensor:˜ G ai , b j = ∑ κ = ∑ µν D ( κ ) ai µν G ( κ ) µν ( D b j ) , (46)with the transition density matrix D ai having the elements D ( κ ) ai , µν = (cid:126) C † µ a σ κ (cid:126) C ν i . (47)Eq. (43) is solved within a preconditioned conjugate gradientalgorithm. Thereby, as first guess, we employ trial vectors ˜ T that represent the uncoupled solutions:˜ T ( ) = − H uoPV , MO ◦ ∆ − (48) R ( ) = ˜ T ( ) ◦ ∆ (49) P ( ) = ˜ T ( ) . (50)In an iterative procedure in each step i from the trial vector˜ T ( i − ) , we construct the perturbed density matrices ( D (cid:48) ( κ ) ) ( i ) following Eq. (40) and calculate the contracted electronicHessian as˜ E [ ] , ( i ) = ˜ T ( i − ) ◦ ∆ + G uoMO ( D (cid:48) ( ˜ T ( i − ) )) , (51)with the two-index MO transformed two-electron matrix: G MO , ai ( D (cid:48) ) = ∑ κ = ∑ µν (cid:126) C † µ a σ κ (cid:126) C ν i G ( κ ) µν ( D (cid:48) ) , (52)From this we update α ( i ) = r ( i − ) Tr (cid:16)(cid:0) ˜ T ( i − ) (cid:1) † ˜ E [ ] , ( i ) (cid:17) (53) R ( i ) = R ( i − ) − α ( i ) ˜ E [ ] , ( i ) (54) P ( i ) = − R ( i ) ◦ ∆ − (55) r ( i ) = (cid:12)(cid:12)(cid:12)(cid:12) Tr (cid:18)(cid:16) R ( i ) (cid:17) † P ( i ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) (56)For all calculations in this paper, the Fletcher-Reeves weightof the precondition was employed, which is calculated as β FR = r ( i ) r ( i − ) and the algorithm is followed until convergenceof the norm (cid:113) r ( i ) r ( ) to a given threshold. The trial matrix is updated as˜ T ( i ) = P ( i ) + β ˜ T ( i − ) . (57)Upon convergence, the solution is received as the sum over allweighted trial matrices from the N iter iterations: T uo = N iter ∑ i = α ( i ) ˜ T ( i − ) (58)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 W (cid:48) ( κ ) µν ( T uo ) = occ ∑ i unocc ∑ a (cid:104) (cid:126) C † µ a σ ( κ ) (cid:126) C ν i (cid:0) ε i T ∗ ai + ( F (cid:48) ai ) ∗ (cid:1) + (cid:126) C † ν i σ ( κ ) (cid:126) C µ a (cid:0) ε i T ai + F (cid:48) ai (cid:1)(cid:105) , (59)with the MO transformed Fock matrix being F (cid:48) = H uoPV , MO + G uoMO ( D (cid:48) ( T uo )) . (60)In case of DFT contributions from the exchange correlationpotential to the perturbed G matrix have to be calculated viaderivatives of the exchange correlation potential as shown asin 31 for the case of spin unpolarised LDA; for a more generaldiscussion of exchange-correlation functionals see AppendixB and Ref. 70). III. COMPUTATIONAL DETAILS
In the pilot implementation of the PV energy gradient, un-perturbed LCAO coefficients ( C ), orbital energies ( ε ), PV op-erators as well as two electron integrals ( G ) were computedwith a modified version of the Turbomole programpackage. Therein, the conjugate gradient algorithm forsolving the linear response equations [Eqs. (48-58)] was im-plemented within MATLAB and the modified Turbomoleprogram was called to compute the perturbed G matrix in AObasis following Eq. (22). The resulting perturbed density ma-trices [Eqs. (40) and (59)] were used in a modified version ofthe gradient implementation of Ref. 65 to assemble the ana-lytic gradient of the parity-violating potential.The calculation of the two-electron part (HF case andCoulomb contribution in LDA) of the perturbed contractedtwo-electron matrix G is carried out via contraction of the twoelectron tensor ( µν | ρσ ) with the two-particle density matrix Γ ( κ , λ ) which can be constructed from the one-particle densitymatrices as Γ ( κ , λ ) µνρσ ( D ( κ ) , D ( λ ) ) = a κ , λ C D ( κ ) µν D ( λ ) ρσ − a κ , λ X D ( κ ) µσ D ( λ ) ρν (61)Note, that for the employed ZORA operator without two-electron spin-orbit or spin-spin coupling terms, only two-electron densities with κ = λ are required at the HF or DFTlevel. Furthermore, we introduced the scaling factors for thedirect Coulomb contribution a κ , λ C , which is in all present cal-culations set as a κ , λ C = δ κ δ λ and a more flexible scaling pa-rameter for the exchange contribution a κ , λ X , which is in thepresent implementation set to be constant a κ , λ X = a X .For the calculation of the two-electron contribution to thePV energy gradient, the perturbed two particle density matri-ces are needed: Γ (cid:48) ( κ , λ ) µνρσ = Γ ( κ , λ ) µνρσ ( D (cid:48) ( κ ) , D ( λ ) ) = a κ , λ C D (cid:48) ( κ ) µν D ( λ ) ρσ − a κ , λ X D (cid:48) ( κ ) µσ D ( λ ) ρν (62) Whereas in the present implementation Eq. (62) is em-ployed directly, in the pilot implementation the perturbed two-electron density matrices were calculated for practical reasonsvia Γ (cid:48) ( κ , λ ) (cid:39) (cid:104) Γ ( κ , λ ) (cid:16) D ( κ ) + D (cid:48) ( κ ) , D ( λ ) + D (cid:48) ( λ ) (cid:17) − Γ ( κ , λ ) (cid:16) D ( κ ) , D ( λ ) (cid:17) − Γ ( κ , λ ) (cid:16) D (cid:48) ( κ ) , D (cid:48) ( λ ) (cid:17)(cid:105) , (63)where a constant, numerical scaling factor of2 √ / [ α G F / ( E h a )] ≈ . × was used to scale upartificially the very small numerical values in the perturbeddensity matrices to avoid substractive cancellation in Eq. 63.The present implementation of (cid:126) ∇ E PV as defined in Eq. (28)as well as the response equations [Eqs. (48-58)] were includedin a modified version of the Turbomole program. Theresults from the current production-level implementation areidentical to those of the pilot implementation. The two differ-ent implementations provided an internal test of our results.The implementation proceeds as follows:1. For a given molecular structure, unperturbed LCAO co-efficients and orbital energies are received from a SCFcomputation at the two-component ZORA level withthe modified version of Turbomole andwritten on disk.2. Integrals of the PV operator [Eq. (32)] are computedand from this the initial guess for the response equa-tions is formed [Eq. (41)] with the program describedin Ref. 77.3. The linear response equations are solved in an itera-tive manner following Eqs. (48-58) within the programof Ref. 77 and the perturbed density matrices and per-turbed energy weighted density matrix are written ondisk.4. Within a modified version of the program described inRef. 65 the gradient integrals of the PV operator [Eq.(33-38)], of the ZORA operator (cid:126) ∇ (cid:126) η h ( κ ) ZORA , µν and of thetwo-electron integrals (cid:126) ∇ (cid:126) η ( µν | ρσ ) are computed andcombined with the unperturbed and perturbed densitymatrices following Eq. (28). For two-electron contribu-tions the perturbed two electron matrix is formed [Eq.(62)] and contracted with the gradient of two-electronintegrals (cid:126) ∇ (cid:126) η ( µν | ρσ ) .In the following, the scheme described above for comput-ing (cid:126) ∇ E PV is used to study electroweak parity-violating ef- fects in halogenated methane derivatives; CHBrClF, CHClFI,CHBrFI and CHAtFI. As mentioned in the Introduc-tion, their vibrational spectra, in particular for CHBrClFhave been extensively studied theoretically andexperimentally to detect PV effects. A high reso-lution of 5 × − has been achieved for the C-F stretchingmode of CHBrClF with CO laser spectroscopy, which isnevertheless about three or four orders of magnitude largerthan the theoretical predictions of the size of the effectfor the molecule under investigation. Therefore, the vi-brational and rotational spectra of heavier homologues willbe studied herein, as has been done previously for the vi-brational frequency shifts in a one-dimensional anharmonicapproximation. Whereas in the previous work single-modeanharmonic effects were included variationally by solving theone-dimensional anharmonic vibrational Schrödinger equa-tion, we will use in the present work vibrational perturba-tion theory to account for anharmonic effects. Most impor-tantly, we can include now also multi-mode contributions toparity-violating frequency shifts efficiently, which were ne-glected in essentially all previous studies on parity-violatingfrequency shifts in chiral molecules with the notable exceptionof Ref. 53, where CDBrClF was studied in a four-dimensionalanharmonic model that span the C–F stretching as well as C-Dstretching and the two C-D bending modes.For all the methane derivatives mentioned above, we use thesame molecular structures, electronic structure methods andbasis sets from earlier work to allow for direct comparison.As in Ref. 40, the S -enantiomers of a given chiral compoundis considered when parity-violating potentials, and in the cur-rent work also gradients of the parity-violating potentials, arereported. Splittings of a property A between enantiomers aregiven as ∆ A = A S − A R . The equilibrium structures, harmonicvibrational frequencies and potential energy surfaces (PES)were computed on the CCSD(T) level with the cc-pVDZ ba-sis set for the first to third row elements and quasi-relativisticStuttgart pseudopotentials in the neutral atom reference sys-tem together with energy optimised valence basis sets for Br,I and At. Equilibriums structures, harmonic vibrationalfrequencies, corresponding normal coordinates and displacedstructures as reported in Ref. 40 were reused in the presentwork. The PES have been determined by the SURF module of MOLPRO. . The double zeta basis is expected to be-have comparatively poorly, but is kept herein to speed up thecalculations and to obtain results that are directly compara-ble to previous work. An improved description of the anhar-monic PES, in particular for the promising astatine containingmethane derivative will be left for a later study.In computations of the PV operator a Weinberg parameterof sin θ W = . Q W , A = ( − θ W ) Z A − N A , with Z A beingthe number of protons and N A being the number of neutronsin nucleus A . In all two-component ZORA calculations aneven tempered basis set, with the exception of an uncontractedaug-cc-pVDZ basis set for hydrogen, has been used in orderto compare the results to previous works. The parametersof the even tempered series are α i = γβ N − iN , i = , . . . , N , with N = γ = .
02 and α = α − and α − have been used, while in thecase of d functions, α − has been chosen for elements ofthe second and third row of the periodic table of elements, α − for the fourth and fifth row and α − for the sixthrow, which is the only one to contain f functions as well withexponents α − .For LDA, Dirac exchange and VWN5 correlation po-tentials have been used. For LDA a standard DFT integra-tion grid was used, whereas matrix elements of the ZORAand PV operators were computed on a very dense grid. MOshave been converged until the relative change of the SCF en-ergy and spin-orbit energy dropped below 10 − E h and 10 − E h respectively. The threshold for negelection of gradients oftwo-electron integrals was set 10 − E h .Following the regular spectroscopic notation, the normalmodes of the methane derivatives are named in descendingorder of the frequency values. IV. RESULTS AND DISCUSSION
We first make a comparison of the directional derivatives ofthe PV energy along the C-F stretching mode ( ν ) as obtainedfrom the analytical PV energy gradients with those of the nu-merical ones from an earlier study. The PV energy gradientsare utilised for estimating the shifts of the rotational constantscorresponding to the equilibrium structures. The vibrationalfrequency shifts in vibrational transitions for the C-F stretch-ing mode are calculated within a perturbative treatment usingthe PV energy gradients. At the end, the multi-mode effects inthe 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 char-acteristic C-F stretching mode ( ν ) which is amenable to high-resolution CO laser spectroscopy. The directional deriva-tive of the PV energy along the C-F stretching normal co- TABLE I. Analytical directional derivative ( (cid:126) ∇ E PV ( (cid:126) R ) · ˆ u q in 10 − cm − ) along the dimensionless reduced normal coordinates corre-sponding to the C-F stretching mode ( ν ) of the chiral halogenatedmethane derivatives as computed at the equilibrium structure.molecules (cid:126) ∇ E PV ( (cid:126) R ) · ˆ u q HF LDACHBrClF analytical 0 . . a . ( ) . ( ) CHClFI analytical 4 . . a . ( ) . ( ) CHBrFI analytical 7 . . a . ( ) . ( ) CHAtFI analytical 60 . − . a . ( ) − . ( ) a The numerical derivative was obtained in Ref. 40 from a polynomial fit ofa one-dimensional cut through the parity-violating potential along thedimensionless reduced normal coordinates q . ordinate ( q ) is obtained herein by projecting the analyticalCartesian PV gradient (cid:126) ∇ E PV onto the displacement vector ˆ u q corresponding to a unit shift along the dimensionless reducednormal coordinate q . Table I depicts the HF and LDA leveldirectional derivative along C-F stretching mode at the equi-librium structure. See the Supplementary information for thecorresponding Cartesian HF and LDA level (cid:126) ∇ E PV for each ofthe equilibrium structures.In reference to Table I, we note that a direct comparison ofthe analytical and numerical directional derivative of the par-ity violating potential is partially hampered by the fact that theanalytical gradient is calculated at the equilibrium structureof the molecule, whereas the numerical directional derivativewas determined in Ref. 40 from the linear term of a polyno-mial fit of the parity-violating potential as calculated alonga one-dimensional cut along the C-F stretching normal co-ordinate in the range q = − . . .
3. The two become bettercomparable, if one computes also the analytical directionalderivatives at several points along the same normal mode dis-placement and performs a polynomial fit. Therefore, as a nextstep, the coordinate dependence of the PV gradient is exam-ined, checking the higher order terms as well, by fitting the (cid:126) ∇ E PV · ˆ u q r to the third order polynomial. (cid:126) ∇ E PV ( q r ) · ˆ u q r ≈ ∑ i = a i q ir (64)where ˆ u q r is the unit vector along the normal coordinate q r and a i corresponds to i ! ∂ i (cid:126) ∇ E PV · ˆ u qr ∂ q ir (cid:12)(cid:12)(cid:12)(cid:12) q r = in a Taylor series ex-pansion of (cid:126) ∇ E PV · ˆ u q r . Thus, a is similar to the value of the (cid:126) ∇ E PV · ˆ u q r at the equilibrium structure along the normal mode ν r . The other fitting coefficients a , a and a respectively rep-resent the first, second and third order derivatives of (cid:126) ∇ E PV · ˆ u q r TABLE II. Fitting coefficients of the HF PV energy ( E PV ) and theHF PV energy gradients ( (cid:126) ∇ E PV ) along the C-F stretching mode ( ν )of the chiral halogenated methane derivatives in 10 − cm − .molecules E PVa (cid:126) ∇ E PV CHBrClF b − . ( ) b . ( ) . ( ) a b − . ( ) − . ( ) a /2b . ( ) . ( ) a /4CHClFI b − . ( ) b . ( ) . ( ) a b − . ( ) − . ( ) a /2b . ( ) . ( ) a /3b . ( ) . ( ) a /4CHBrFI b − . ( ) b . ( ) . ( ) a b − . ( ) − . ( ) a /2b . ( ) . ( ) a /3b . ( ) . ( ) a /4CHAtFI b − . ( ) b . ( ) . ( ) a b − . ( ) − . ( ) a /2b . ( ) . ( ) a /3b − . ( ) − . ( ) a /4 a Reference with respect to the dimensionless reduced normal coordinate q r for ν r mode.For this purpose, structures, which are displaced along thenormal coordinates of a vibrational mode are employed togenerate a one-dimensional cut through the PV potential en-ergy surface for q that varies from − + q = ν ). Figures 1–4 display the HF and LDA level (cid:126) ∇ E PV ( q ) · ˆ u q along the dimensionless reduced normal coordinates cor-responding to the C-F stretching mode of the methane deriva-tives. The fitting coefficients need to be adapted to comparethe structure dependence of the PV energy E PV ( (cid:126) q r ) with thedirectional derivative of the PV energy (cid:126) ∇ E PV · ˆ u q along thedimensionless reduced normal coordinates q correspondingto the C-F stretching mode, because the former can be fittedto a polynomial expansion, too, namely E PV ( q r ) ≈ ∑ i = b i q ir , (65)where b i corresponds to i ! ∂ i E PV ( q r ) ∂ q ir (cid:12)(cid:12)(cid:12) q r = in a Taylor seriesexpansion of E PV . Thus, one has to contrast E PV ( (cid:126) q r ) with (cid:82) (cid:126) ∇ E PV ( q r ) · ˆ u q r d q r and hence the values of a , a / a / a / b , b , b and b respec-tively. The term b gives the PV energy at the equilibriumposition. FIG. 1. HF and LDA level PV energy gradient (fitted to ∑ i = a i q ir )along the C-F stretching mode of CHBrClF. See Table II for thefitting coefficients ( a i ) due to HF PV energy gradient. For LDAPV energy gradient, the fitting coefficients (in 10 − cm − ) are a = . ( ) , a = − . ( ) , a = − . ( ) and a = − . ( ) respectively.FIG. 2. HF and LDA level PV energy gradient (fitted to ∑ i = a i q ir )along the C-F stretching mode of CHClFI. See Table II for the fittingcoefficients ( a i ) due to HF PV energy gradient. For LDA PV energygradient, the fitting coefficients (in 10 − cm − ) are a = . ( ) , a = − . ( ) , a = − . ( ) and a = − . ( ) respec-tively. The comparison of these fitting coefficients for PV energy(taken from Ref. 40) and directional derivative of the PV en-ergy along the C-F stretching mode for the HF level is pre-sented in Table II. To test the reliability of the analytical (cid:126) ∇ E PV ,it has also been compared to the numerical (cid:126) ∇ E PV · ˆ u q of theearlier work. Both numerical and analytical (cid:126) ∇ E PV · ˆ u q val-ues for all the systems agree well. The numerical values, in-cluding the averaged numerical and averaged analytical gra-dient match within the achieved precision, with the numbersin parentheses denoting the error due to the fit procedure. Thecorresponding comparisons of these fitting PV energy and PVenergy gradient coefficients along the C-F stretching mode atthe LDA level are given in the Supporting Information.For the determination of the second derivatives ∂ E PV ∂ q r (sameas a in Table II), which are needed for calculating the vi-brational energy levels as shown in Eqs. 8 and 9, one canexpect the availability of the analytical gradient to simplifythe calculations, because fewer points are needed from thePV energy surface. As a consistency check, a linear fit ofthe analytically calculated directional derivatives computed at q = + . q = q = − .
125 has been opted for in all of1
FIG. 3. HF and LDA level PV energy gradient (fitted to ∑ i = a i q ir )along the C-F stretching mode of CHBrFI. See Table II for the fittingcoefficients ( a i ) due to HF PV energy gradient. For LDA PV energygradient, the fitting coefficients (in 10 − cm − ) are a = . ( ) , a = − . ( ) , a = − . ( ) and a = − . ( ) respec-tively.FIG. 4. HF and LDA level PV energy gradient (fitted to ∑ i = a i q ir )along the C-F stretching mode of CHAtFI. See Table II for the fittingcoefficients ( a i ) due to HF PV energy gradient. For LDA PV energygradient, the fitting coefficients (in 10 − cm − ) are a = − . ( ) , a = − . ( ) , a = . ( ) and a = − . ( ) respectively. the molecules, where the fitting values ( a ) are found to be ingood agreement with the values presented in Table II. Theselinear fit values are reported in the supplementary material. B. Shifts of the rotational constants
The values shown in Tables I and II reveal that the equilib-rium structure of these halogenated methane derivatives pos-sess non-zero (cid:126) ∇ E PV . Due to these non-zero (cid:126) ∇ E PV values, thePV potential can induces a minute change in the equilibriumstructure. This leads to a shift of the rotational constants,which could in principle be measured by microwave spec-troscopy. As already been discussed in the Introduction, thechange of the structure and hence the shifts in the rotationalconstants due to the existence of non-zero (cid:126) ∇ E PV at the mini-mum of the parity conserving potential is calculated with thehelp of the vibrational Hessian F (see Eqs. 4 and 6). Thechange of the inertia tensor ( ∆ I x ) in the principal axis systemand subsequently the change of the rotational constants ( ∆ X R )can be approximated assuming that the result depends linearlyon the displacements. Shifts in the rotational constants are, in general, sensitive tothe uncertainty in the PV energy gradients (cid:126) ∇ E PV , so that onlythree digits are estimated to be reliable ( cf. Table III). Thenumbers obtained on the LDA levels are found to be aboutone quarter of the HF results for the CHBrClF molecule. Thecalculated relative effect in the microwave spectrum ( ∆ X / X ≈ − ) of CHBrClF, although it agrees well with the earlierfinding by Quack and Stohner , is still far from the currentexperimental resolution. Similarly, the corresponding relativeshifts in CHClFI and CHBrFI are about one order of magni-tude larger in absolute value than for CHBrClF, but still far be-low the present experimental resolution (see Table III). Con-sidering the current scenario, the PV induced shifts of the ro-tational constants in CHBrClF, CHClFI and CHBrFI are notexpected to be measurable with present experimental setups.On the other hand, due to the heavy mass of astatine, a largershift in the rotational constants is expected for CHAtFI. Forthis heavier At analogue, shifts are found to be about two orthree-orders of magnitude higher in absolute value than thelighter molecules considered herein.Another comment as to the accuracy of these estimates ofrotational energy shifts is in order: Herein, we have followedthe most simple approach and computed only the PV shiftat the equilibrium structure as it would arise from the PV gra-dient contribution. An improved treatment would include es-timates of PV induced shifts of vibrationally averaged rota-tional constants as well as their influence on Coriolis couplingterms and centrifugal distortion constants. C. Vibrational transitions for C-F stretching mode
For the vibrational spectrum, the results from the separableanharmonic adiabatic approximation (SAAA) used in ear-lier work have been compared to the one-dimensional (1D)perturbative approach (Eq. 9) of the present work. The influ-ence of multi-mode coupling terms to the vibrational energylevels is also accounted for after adding the perturbative multi-mode contribution (Eq. 10) to the single-mode results. TableIV presents the vibrationally averaged HF and LDA parity vi-olating potential energy levels up to 3 rd vibrational state forthe C-F stretching mode ( ν ) of the S -enantiomer of CHBr-ClF, CHClFI, CHBrFI and CHAtFI molecules along with acomparison to SAAA results. For the lowest transition (from n = n =
1) in CHBr-ClF, the full 1D and the perturbative 1D treatment match toabout < cf. Table IV). A similar influence of non-separable anharmonic effects on PV has been reported forsimilar molecules. Nonetheless, the conclusion to be drawnis that the effect of the 2D treatment of the parity conserv-ing potential is more important than the higher order terms ofthe 1D parity-violating potential. When comparing the differ-ent electronic structure methods, one finds significant differ-ences, but the values for vibrational frequency shifts in CHBr-2
TABLE III. Relative shift of the rotational constants in chiral halogenated methane derivatives.HF LDA ∆ A / A ∆ B / B ∆ C / C ∆ A / A ∆ B / B ∆ C / C CHBrClF 2 . × − . × − . × − . × − . × − . × − CHClFI 1 . × − . × − . × − . × − . × − . × − CHBrFI 1 . × − . × − . × − . × − − . × − . × − CHAtFI 6 . × − . × − . × − . × − . × − . × − TABLE IV. Vibrationally averaged HF and LDA parity-violating potential E Sn , PV for energy levels n in 10 − cm − for the C-F stretchingmode ( ν ) of the (S)-enantiomer of CHBrClF, CHClFI, CHBrFI and CHAtFI.HF LDAMolecule n Full 1D a Perturbed 1D Perturbed 2D 2D effects (%) b Full 1D a Perturbed 1D Perturbed 2D 2D effects (%) b CHBrClF 0 − . − . − . − .
80 0 . . . . − . − . − . − .
47 0 . . . . − . − . − . − .
28 0 . . . . − . − . − . − .
14 0 . . . . ← . . . .
32 0 . . . . − . − . − . − .
11 5 . . . . − . − . − . − .
37 5 . . . . − . − . − . − .
64 6 . . . . − . − . − . − .
93 6 . . . . ← . . . .
06 0 . . . . − . − . − . .
32 20 . . . − . − . − . − . .
98 21 . . . − . − . − . − . .
68 22 . . . − . − . − . − . .
42 23 . . . − . ← .
962 0 . . − .
79 1 . . . − . − . − . − . − .
96 938 .
77 938 .
061 966 .
961 3 . − . − . − . − .
88 887 .
70 883 .
325 970 .
026 9 . − . − . − . − .
69 840 .
46 828 .
589 973 .
091 17 . − . − . − . − .
41 797 .
14 773 .
853 976 .
155 26 . ← − . − . − . − . − . − .
736 3 . − . a Reference = (
Perturbed2D − Perturbed1D ) / Perturbed1D
ClF are in the same order of magnitude, deviating less than30%, with the same sign, the latter of which is not the casefor the PV energy shifts. When compared with absolute val-ues, the 2D contributions are found to increase for the highervibrational levels. When these calculated PV frequency shifts(Table IV) are multiplied by two and then divided by the corre-sponding vibrational transition frequencies, they produce thedimensionless parity-violating relative vibrational frequencysplittings between the C-F stretching fundamental of S- andR-enantiomers. The calculated relative vibrational frequencysplittings corresponding to the lowest transition for the C-Fstretching vibration of CHBrClF molecule is found to be ingood agreement with the earlier theoretical values.
For CHClFI, these parity-violating 1 ← cf . Table IV). PV transition frequencyshifts in the C-F stretching fundamental of CHBrFI due to full1D and perturbed 1D are roughly twice to that of the valuesfor CHClFI, whereas perturbed 2D effects in CHBrFI mea-sures about 75% higher than the 2D contributions in CHClFImolecule. Akin to CHBrClF, the difference between the full1D and perturbed 1D approach is again negligible, whereas3 FIG. 5. HF parity-violating energy along the dimensionless reducednormal coordinates of CHBrClF for all the normal modes.
2D contributions are significant.A similar observation is also noticed for CHAtFI, whichshows the largest absolute value of the PV vibrational fre-quency difference in the C-F stretching fundamental amongthe four halogenated methane derivatives discussed herein.This large shift is expected and mainly attributed to the pres-ence of the heavier astatine nucleus. The full 1D and per-turbative 1D treatment do not differ much. When taking the2D terms into account, the values change significantly withthe multi-mode contributions causing a reduction in abso-lute value compare to the perturbative 1D estimate of the fre-quency shift. This is expected to result from the fact, that thecouplings of the C-F stretching mode ν to the C-H bendingmode ν , where the frequencies are close ( ∆ ω ≈ ν , where there is a factor of aboutthree ( ≈ n = n =
1) are altered.Thus, the 2D perturbative terms are found to be importantin all cases, which increase gradually for higher vibrationallevels. Thus, the multi-mode effects due to the perturbative2D approach increases the fundamental vibrational frequencyshift for the CHBrClF, CHClFI molecules whereas the shiftsin CHBrFI molecule are lower in absolute value as comparedto the perturbative 1D treatment counterparts. The heaviestAt derivative shows the maximum frequency shift after theinclusion of the perturbative 2D effects at both HF and LDAlevel.
D. Multimode effects in CHBrClF and CHAtFI
The influence of multi-mode effects from non-separable anharmonic effects have been estimated for CHBr-ClF and the heaviest CHAtFI molecules after calculating thePV energy gradients along all the normal modes. Figures 5and 6, respectively, display the variation of the HF and LDAPV energies along the dimensionless reduced normal coordi-nates of all the normal modes of CHBrClF molecule. In bothcases, normal coordinates q , q , q and q show steep slopesfor the PV energy with respect to displacements along therespective dimensionless reduced normal coordinates. Thismay be naively expected, because the bending motion can al- FIG. 6. LDA parity-violating energy along the dimensionless re-duced normal coordinates of CHBrClF for all the normal modes. ter the chiral structure of the molecule to a larger extent thana stretching vibration. In view of this, the vibrationally av-eraged PV potential has been computed along all the normalmodes of CHBrClF and the heavier CHAtFI analogue. Akinto the observation in CHBrClF molecule, the bending vibra-tional modes in CHAtFI exhibit larger deviation in the PV en-ergies as compared to other normal modes and hence, largervibrational splittings and an enhancement of the perturbative2D effects can be expected due to these modes.Now, the vibrationally averaged HF and LDA PV gradientsare utilised for calculating the 1D (Eq. 9) and 2D (Eq. 8)perturbative treatment for the vibrational energy levels up to3 rd vibrational states. For this purpose, in addition to the PVenergy gradient (cid:126) ∇ E PV , which is available along all the normalmodes, the second derivative ( a in Table II which is ∼ ∂ E PV ∂ q r in Eqs. 9 and 8) of the PV potential is needed. But it is notnecessary to compute a full profile, because with the avail-ability of the analytic gradient, the second derivative ( a ) canbe computed from the linear fit of only a few energy gradi-ent points close to the equilibrium structure ( q = − .
125 and q = + . q = . n = n =
1) for all the vibrationalmodes are estimated from Eq. 8. The q , q , q and q bending normal coordinates (although the corresponding fun-damentals give rise to the least intense peaks in the infraredspectrum) contribute significantly to the multi-mode effectsin the various fundamental transitions. Figure 7 displays thecontributions towards the multi-mode effects for each of thevibrational modes of CHBrClF and CHAtFI molecules at HFand LDA level of theory. Contributions from ν (in red) and ν (in orange) are the largest and they contribute to almostevery fundamental, whereas q and q contribute moderatelyin CHBrClF. The normal coordinate of the mode with themost intense fundamental in the vibrational spectrum, the C-F stretching fundamental ( ν ), and other normal coordinateshardly contribute to multi-mode effects in other fundamentals.The observations are slightly different in CHAtFI. The contri-butions from q (At-C-F bending in red) and q (F-C-I bend-ing in orange) normal coordinates are significant. In additionto this, the C-At stretching vibration ( q ) has shown a substan-tial contribution to the multi-mode effects. This is expected,4 FIG. 7. Contributions (in cm − ) to multi-mode effects from each of the normal coordinates to other fundamental vibrational transitions (1 ← S )-enantiomer of CHBrClF and CHAtFI molecules at HF and LDA level of theory.FIG. 8. 1D and 2D perturbative treatment of the vibrationally av-eraged HF and LDA parity-violating potential E Sn , PV leading to theshown PV frequency shifts (in cm − ) for the fundamental transition(1 ←
0) in all the normal modes of the ( S )-enantiomer of CHBrClFand CHAtFI molecules. as astatine is the heaviest nucleus of the molecule and, dueto the steep Z -scaling of parity-violating effects, with Z beingthe nuclear charge, should also have a pronounced influenceon the parity violating shifts. The C-F stretching fundamen- tal ( ν ) gets a large contributions from the H-C-F bending ( q in olive) normal coordinate, most probably due to the strongcoupling between themselves since their frequencies are veryclose. Otherwise, the contributions from other vibrations arenot that significant.The total observed shifts associated with the fundamentaltransition (from n = n =
1) for all the vibrational modesof the ( S )-enantiomer of these two chiral derivatives from the1D and 2D perturbative treatment are displayed in Figure 8.It should be noted that one needs to multiply these values by2 and then divide by the vibrational frequencies of the cor-responding fundamentals to get the relative PV vibrationalfrequency splittings between the C–F stretching fundamentalof the S- and R-enantiomers of CHBrClF and CHAtFI. Thecorresponding PV shifts up to the 3 rd vibrational energy levelas well as the PV shifts in the fundamental transitions (from n = n =
1) are included in the supporting material.For CHBrClF, the vibrational frequency shifts from the 1Dperturbative treatment evaluated at HF level are larger in mag-nitude for ν , ν , ν and ν than that of the ν C-F stretch-ing mode. The steeper slope of the curves due to the largevariation in the PV energies (which increases the numericalvalue of (cid:126) ∇ E PV ) along q , q , q and q in Figure 5 favoursthese large vibrational shifts in the bending modes. On theother hand, the 2D perturbation results in large shifts for mostmodes as compared to the C-F stretching mode ( ν ). This isexpected since the parity-violating multi-mode effects are sig-nificant in other bending modes and stretching modes otherthan the C-F stretching vibration (See Figure 7). From thecomparison of 1D and 2D effects in frequency shifts within in-dividual modes, it is found that the 2D effects increase the fre-quency shifts only in ν , ν , ν and ν , whereas other funda-5mentals follow the opposite trend ( cf. Figure 8) where a can-cellation of the 1D contribution is caused due to the 2D treat-ment. Total calculated PV shifts in the fundamentals transi-tions from both of the perturbative 1D and 2D treatments aligneither towards positive or negative directions. No alternationof the sign is observed even after adding multi-mode contri-butions to 1D perturbative ones. Overall, the estimated vibra-tional shift in all the fundamental transitions matches to thesame order of magnitude of previous theoretical values.
The case for CHAtFI is slightly complicated. The 2D termsare important for all modes where ν , ν and ν fundamen-tals show frequency shifts that are larger in absolute value ascompared to the C-F stretching ( ν ) mode ( cf. Figure 8). Mul-timode effects increase the frequency shifts by a factor of ∼ ν when compared to the perturbative 1D corrections.The correction due to multi-mode effects for the most intenseC-F stretching mode reduces the absolute value as comparedto the perturbative 1D result, i.e. the perturbative 2D treatmentdecreases the energy gap between the ground ( n =
0) and 1 st ( n =
1) vibrational state. Even, as mentioned above (see TableIV), at LDA level the sign is altered ( cf.
Figure 8). A substan-tial increase in the ν fundamental at LDA level is seen afterincluding the multi-mode effects which is about 215% of theperturbative 1D value (See Figure 8). Although the funda-mentals ν and ν have relatively larger frequency shifts forthe vibrational transition from n = n =
1, the values de-crease when multi-mode effects are introduced. Multimodecontributions in ν increases the frequency shifts in magni-tude as compared to 1D effects. It may be noted that the asso-ciated transition frequencies of these fundamentals are not inthe range of the CO laser used in the previous experimen-tal set up. But, in spite of their lower intensity in the infraredspectrum, this might not be a severe limitation due to new de-velopments in contemporary laser technology with quantumcascade laser being available in a broader frequency range.We should emphasise, finally, that we reused in our presentwork the equilibrium structures and harmonic vibrationalforce fields from a previous study to allow for direct com-parison of results, Whereas a sophisticated CCSD(T) ap-proach was used in Ref. 40, the basis sets were only of limitedquality, which impacts not only on the equilibrium structuresand harmonic force constants obtained, but also on the qual-ity of molecular properties as was noted in our recent study of nuclear electric quadrupole coupling constants in CHBrClFand CHClFI. For CHBrClF and its deuterated isotopomer, im-proved anharmonic ab initio force fields have been reportedfor instance in Ref. 89. In particular for the compound withthe radioactive halogen, CHAtFI, for which no experimen-tal information is available and for which very pronouncedmulti-mode effects on the C-F stretching fundamental werecomputed herein, a study with an improved description of theconventional parity-conserving effects is indicated. This be-comes even more important by virtue of the potential role ofheavy elemental chiral molecules in the search for dark mat-ter candidates and recent advances in laser spectroscopyof radioactive molecules with short-lived nuclei. V. CONCLUDING REMARKS
In the present article, the gradient of the molecular parity-violating (PV) potential has been derived and implementedwithin a quasirelativistic framework at the level of HF andLDA, and is applied to rotational and vibrational spectroscopyof polyatomic molecules. A systematic study of frequencyshifts in the rotational and vibrational spectra of chiral poly-halomethanes has been carried out.A one-dimensional perturbative treatment has been com-pared to a previously calculated full anharmonic one-dimensional approach and found to be in good agreement. Forthe hypothetical astatine compound, the predicted frequencysplitting is in the right order of magnitude to be measuredwhen treating the problem in the one-dimensional approxima-tion. The results form the multi-mode effects in CHBrClF andCHAtFI reflects the importance of nonseparable anharmoniceffects and suggests that the C–F stretching mode is not ide-ally suited for a measurement of electroweak PV effects, sinceother bending and stretching modes contribute significantlyand can not be neglected. In view of this, choosing a propervibrational transition is equally important as choosing a suit-able molecule for conducting PV experiments and calculatingthe PV effects on the vibrational transition of chiral deriva-tives. Multimode effects are crucial for getting more insightsabout the vibrational transitions as well as rotational spectrain chiral compounds containing heavier atoms.The present analytic derivative approach within a quasi-relativistic mean-field scheme leads to an important simpli-fication for routine calculation of PV frequency shift in rota-tional spectra of chiral molecules and of multi-mode contribu-tions to the PV frequency shifts in rovibrational spectroscopy.The implementation of the analytic (cid:126) ∇ E PV at other DFT levelssuch as GGA (e.g. BLYP) and GGA-based hybrid-functionals(e.g. B3LYP) is currently underway, which will increase itsapplicability to diverse chiral systems containing heavier nu-clei, in particular transition metals. Also, such a developmenthelps in measuring the functional dependencies and influenceof electron correlation effects for the theoretical estimation ofthe PV effects in the rotational and vibrational transitions ofchiral compounds. This will provide valuable information forfuture experiments aiming at the detection of molecular parityviolation. ACKNOWLEDGMENTS
The authors are particularly thankfull to Yunlong Xiao, So-phie Nahrwold, Timur Isaev and Christoph van Wüllen forstimulating discussions. Financial support by the DeutscheForschungsgemeinschaft (DFG, German Research Founda-tion) – Projektnummer 328961117 – SFB 1319 ELCH andthe VolkswagenFoundation is gratefully acknowledged. Thecenter for scientific computing (CSC) Frankfurt is thanked forcomputer time.6
Appendix A: Variational Perturbation Theory
Following Refs. 56–58 the total energy E ( (cid:126) R ) = (cid:68) Ψ (cid:12)(cid:12)(cid:12) ˆ H ( (cid:126) R ) (cid:12)(cid:12)(cid:12) Ψ (cid:69) parametrically depends on the nuclear coor-dinates. Within a mean-field approach, the energy satisfies (cid:126) ∇ (cid:126) t E ( (cid:126) R ) = (cid:126) ∇ (cid:126) t (cid:68) Ψ ( (cid:126) t ) (cid:12)(cid:12)(cid:12) ˆ H ( (cid:126) R ) (cid:12)(cid:12)(cid:12) Ψ ( (cid:126) t ) (cid:69) = (cid:126) , (A1)where (cid:126) t is the vector of orbital rotations and (cid:126) ∇ (cid:126) t is the gradientwith respect to (cid:126) t . The variational energy of slightly displacednuclei E ( (cid:126) R + (cid:126) η ) can be expanded in a Taylor series in thenuclear displacements (cid:126) η as (see e.g. Ref. 59): E ( (cid:126) R + (cid:126) η ) = E ( (cid:126) R ) + (cid:126) ∇ (cid:126) η E ( (cid:126) R + (cid:126) η ) (cid:12)(cid:12)(cid:12) (cid:126) η = (cid:126) · (cid:126) η + (cid:126) η T · (cid:20) (cid:126) ∇ (cid:126) η ⊗ (cid:126) ∇ (cid:126) η E ( (cid:126) R + (cid:126) η ) (cid:12)(cid:12)(cid:12) (cid:126) η = (cid:126) (cid:21) · (cid:126) η + . . . , (A2)where (cid:126) ∇ (cid:126) η is the gradient with respect to the nuclear displace-ment vector.Thus, the total variational energy is a function of λ PV aswell: E ∞ = E ( (cid:126) R + (cid:126) η , λ PV )= E ( (cid:126) R + (cid:126) η ) + ∂ E ( (cid:126) R + (cid:126) η , λ PV ) ∂ λ PV (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ PV = λ PV + ∂ E ( (cid:126) R + (cid:126) η , λ PV ) ∂ λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ PV = λ + . . . . (A3)In combination with Eq. (A2) the parity violating contribution E PV ( (cid:126) R + (cid:126) η , λ PV ) = E ( (cid:126) R + (cid:126) η , λ PV ) − E ( (cid:126) R ) to the gradient ofthe energy with respect to nuclear displacements is: (cid:126) ∇ (cid:126) η E PV ( (cid:126) R + (cid:126) η , λ PV ) (cid:12)(cid:12)(cid:12) (cid:126) η = (cid:126) = ∂ (cid:126) ∇ (cid:126) η E ( (cid:126) R + (cid:126) η , λ PV ) ∂ λ PV (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126) η = (cid:126) λ PV = λ PV + ∂ (cid:126) ∇ (cid:126) η E ( (cid:126) R + (cid:126) η , λ PV ) ∂ λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126) η = (cid:126) λ PV = λ + . . . ≈ ∂ (cid:126) ∇ (cid:126) η E ( (cid:126) R + (cid:126) η , λ PV ) ∂ λ PV (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:126) η = (cid:126) λ PV = λ PV (A4) Appendix B: Spin polarised density functionals
We employ the approximation of non-collinear DFT. Then,we can write a exchange-correlation potential operator ˆ V ( κ ) XC as(see also Ref. 65)ˆ V ( κ ) XC = ˆ V ( κ ) XC , LDA + (cid:126) ˆ V ( κ ) XC , GGA · (cid:126) ∇ , (B1)where (cid:126) ˆ V ( κ ) XC , GGA = δ F XC , GGA δ (cid:126) ∇ ρ ( κ ) e ( (cid:126) r ; D ( κ ) ) and for the explicit definitionsof (cid:126) ˆ V XC see e.g. Refs. 63 and 65.We define the spin or number density function ρ ( κ ) e as ρ ( κ ) e ( (cid:126) r ; D ( κ ) ) = ∑ µν R (cid:110) D ( κ ) µν (cid:111) χ µ ( (cid:126) r ) χ ν ( (cid:126) r ) . (B2)Note, that this definition of the spin density coincides withwhat is usually called in non-collinear DFT the magnetiza-tion density and should not be confused with the length of themagnetization density which is what is called spin density innon-collinear DFT.We can write the exchange-correlation energy as E XC = ∑ κ ∑ µν ℜ D ( κ ) µν (cid:68) χ µ (cid:12)(cid:12)(cid:12) ˆ V ( κ ) XC (cid:12)(cid:12)(cid:12) χ ν (cid:69) (B3)and the gradient with respect to nuclear displacements is (cid:126) ∇ (cid:126) η E XC = ∑ κ ∑ µν ℜ D ( κ ) µν (cid:68) (cid:126) ∇ (cid:126) η χ µ (cid:12)(cid:12)(cid:12) ˆ V ( κ ) XC (cid:12)(cid:12)(cid:12) χ ν (cid:69) . (B4)In presence of a perturbing operator such as the parity-violating potential we have to compute the first order per-turbed exchange-correlation potential operator asˆ V (cid:48) ( κ ) XC = ∑ λ δ ˆ V ( κ ) XC δ ρ ( λ ) e ( (cid:126) r ; D ( λ ) ) ρ ( λ ) e ( (cid:126) r ; D (cid:48) ( λ ) )+ ∑ λ δ ˆ V ( κ ) XC δ (cid:126) ∇ ρ ( λ ) e ( (cid:126) r ; D ( λ ) ) (cid:126) ∇ ρ ( λ ) e ( (cid:126) r ; D (cid:48) ( λ ) ) , (B5)where ρ ( κ ) e ( (cid:126) r ; D (cid:48) ( κ ) ) is the perturbed density function of firstorder in λ PV . Appendix C: Nuclear displacement gradient of a first orderproperty
We can express terms that are proportional to gradients ofthe density matrix in terms of the Fock matrix F ( κ ) ( D ∞ ) = h ( κ ) ZORA + h ( κ ) PV + G ( κ ) ( D ∞ ) (cid:126) ∇ (cid:126) η E ∞ = R (cid:40) ∑ µν (cid:34) ∑ κ = D ∞ , ( κ ) µν (cid:16) (cid:126) ∇ (cid:126) η h ( κ ) ZORA , µν (cid:17) + ∑ κ = D ∞ , ( κ ) µν (cid:126) G ( κ ) grad , µν ( D ∞ ) + ∑ κ = D ∞ , ( κ ) µν (cid:16) (cid:126) ∇ (cid:126) η h ( κ ) PV , µν (cid:17) + ∑ κ = (cid:16) (cid:126) ∇ (cid:126) η D ∞ , ( κ ) µν (cid:17) F ( κ ) µν ( D ∞ ) (cid:35)(cid:41) , (C1)7where we have defined the matrix (cid:126) G ( κ ) grad of contracted gradients of two-electron integrals with elements: (cid:126) G ( κ ) grad , µν ( D ) = ∑ ρσ (cid:20) δ k D ( κ ) ρσ (cid:126) ∇ (cid:126) η ( µν | ρσ ) − a X D ( κ ) ρσ (cid:126) ∇ (cid:126) η ( µσ | ρν ) + a DFT (cid:68) (cid:126) ∇ (cid:126) η χ µ (cid:12)(cid:12)(cid:12) V ( κ ) XC ( D ) (cid:12)(cid:12)(cid:12) χ ν (cid:69)(cid:21) . (C2)The canonical SCF equations are FC = SC ε ; C † SC = , (C3)where S = σ ⊗ S is the overlap matrix constructed from theone component overlap matrix S of the basis functions and F = ∑ κ = σ κ ⊗ F ( κ ) ( D ) . The matrix of coefficients is chosenas C = (cid:18) C ( α ) C ( β ) (cid:19) . (C4)Exploiting the Hermiticity of F and ε we find from this ( (cid:126) ∇ (cid:126) η C † ) FC + C † F ( (cid:126) ∇ (cid:126) η C ) =( (cid:126) ∇ (cid:126) η C † ) SC ε + ε C † S ( (cid:126) ∇ (cid:126) η C ) (C5) Furthermore, the gradient of the orthonormality condition inEq. C3 gives ( (cid:126) ∇ (cid:126) η C † ) SC + C † S ( (cid:126) ∇ (cid:126) η C ) = − C † ( (cid:126) ∇ (cid:126) η S ) C (C6)Therewith, the contribution to the energy gradient from theFock matrix in Eq. (C1) can be expressed via the orthonor-mality condition as ∑ µν ∑ κ = (cid:16) (cid:126) ∇ (cid:126) η D ∞ , ( κ ) µν (cid:17) F ( κ ) µν ( D ∞ ) = ∑ µν ∑ i n i ∑ κ = (cid:104) ( (cid:126) ∇ (cid:126) η (cid:126) C ∞ µ i ) † σ ( κ ) F ( κ ) µν ( D ∞ ) (cid:126) C ∞ ν i + ( (cid:126) C ∞ µ i ) † σ ( κ ) F ( κ ) µν ( D ∞ )( (cid:126) ∇ (cid:126) η (cid:126) C ∞ ν i ) (cid:105) = ∑ µν ∑ i n i ε ∞ i (cid:104) ( (cid:126) ∇ (cid:126) η (cid:126) C ∞ µ i ) † S µν (cid:126) C ∞ ν i + ( (cid:126) C ∞ µ i ) † S µν ( (cid:126) ∇ (cid:126) η (cid:126) C ∞ ν i ) (cid:105) = − ∑ µν ∑ i n i ε ∞ i ( (cid:126) C ∞ µ i ) † ( (cid:126) ∇ (cid:126) η S µν ) (cid:126) C ∞ ν i = − ∑ µν W ∞ , ( ) µν (cid:126) ∇ (cid:126) η S µν (C7)where we introduced the energy weighted density matrix(EWDM) W as W ( κ ) = N orb ∑ i = n i ε i D ( κ ) i . (C8) Appendix D: Linear response equations1. First order perturbed density matrix
The infinite order coefficients can be expressed in terms oforbital rotations, expressed via the anti-Hermitian matrix T ∞ as C ∞ = C e T ∞ (D1)were C is the initial guess of orthonormal orbitals and ( C ∞ ) † = e − T ∞ C †0 . From the definition of the density ma-trix [Eq. (18)] we see D ∞ = (cid:0) C ∞ N ( C ∞ ) † (cid:1) T = ( C ∞ ) ∗ N ( C ∞ ) T , where N is a diagonal matrix containing the occupation num-bers n i . Thus, the perturbed density matrix of infinite order D ∞ can be written in terms of orbital rotation matrix T ∞ as D ∞ = C ∗ e − T T ∞ N e T T ∞ C T . (D2)Moreover, we can write D ∞ in terms of the unperturbed orbitalcoefficients C , received from orbital rotations (cid:126) t in absenceof the perturbation with the corresponding rotation matrix T and an the rotation matrix T as: D ∞ = C ∗ e − T T e − T T N e T T e T T C T = C ∗ e − T T N e T T C T . (D3)When assuming N to be successively sorted for occupied andunoccupied orbitals we can write the transformation matrix T in blocked form as T = (cid:18) T oo − T †uo T uo T uu (cid:19) , (D4)where, o means occupied and u means unoccupied. In firstorder the perturbed density matrix can be written in terms of8the unoccupied-occupied block T uo of T : D (cid:48) = C ∗ (cid:0) NT T − T T N (cid:1) C T = C ∗ (cid:0) NT T + T ∗ N (cid:1) C T = C ∗ (cid:18) T T uo T ∗ uo (cid:19) C T . (D5)Thus, we arrive at D (cid:48) ( κ ) µν ( T uo ) = occ ∑ i unocc ∑ a (cid:104) (cid:126) C † µ a σ ( κ ) (cid:126) C ν i T ∗ ai + (cid:126) C † ν i σ ( κ ) (cid:126) C µ a T ai (cid:105) . (D6)
2. Coupled perturbed HF/KS equations
The coupled perturbed HF (CPHF) or coupled perturbedKS (CPKS) equations read (see Ref. 69). FC (cid:48) + F (cid:48) C = S (cid:48) C ε + SC (cid:48) ε + SC ε (cid:48) (D7)where, a prime denotes matrices linear in λ PV and unprimedmatrices are unperturbed ( λ PV = F (cid:48) is for perturbation independent basis functions as wehave in the case of the perturbation due to the PV potential: F (cid:48) µν = ∑ κ σ κ h ( κ ) PV , µν + σ κ G ( κ ) µν ( D (cid:48) ( T uo )) (D8)where, we introduced the self consistent transformation ma-trix of the orbital coefficients T for which C (cid:48) = CT . FromEq. (D5) we see that only the unoccupied-occupied rotationschange the orbitals in first order, such that we can focus on theof unoccupied-occupied block of Eq (D7). Thus, for pertur-bation independent basis functions the CPHF/CPKS equationsreduce to the coupled equations ε u T uo + C †u F (cid:48) C o = T uo ε o (D9) ε o T ou + C †o F (cid:48) C u = T ou ε u , (D10)with the unoccupied (u) and occupied (o) orbital sub-blocksof C and ε . From this, one arrives at the linear response equa-tions (see e.g. Refs. 70 and 71) ∑ b j (cid:18) A b j B b j B ∗ b j A ∗ b j (cid:19) (cid:18) T b j hT ∗ b j (cid:19) = − (cid:18) H uoPV , MO ( H uoPV , MO ) ∗ (cid:19) , (D11)where the index b runs over unoccpuied orbitals and the index j runs over occupied orbitals. The Hermiticity factor h resultsfrom the orthonormality condition which for perturbation in-dependent basis functions yields C † SCT = − T † C † SC ⇔ T † = − T , which is in accordance with Eq. (D2). For ob-servables such as H PV we have h = + C. S. Wu, E. Ambler, R. W. Hayward, D. D. Hoppes, and R. P. Hudson,“Experimental test of parity conservation in beta decay,” Phys. Rev. ,1413–1415 (1957). T. D. Lee and C. N. Yang, “Question of parity conservation in weak inter-actions,” Phys. Rev. , 254–258 (1956). Y. Yamagata, “A hypothesis for the asymmetric appearance of biomoleculeson earth,” J. Theor. Biol. , 495–498 (1966). É. Gajzágó and G. Marx, “Energy difference of mirror molecules,” AtomkiKözl. Suppl. , 177–184 (1974). V. S. Letokhov, “On difference of energy levels of left and right moleculesdue to weak interactions,” Phys. Lett. A , 275–276 (1975). B. Y. Zel’dovich, D. B. Saakyan, and I. I. Sobel’man, “Energy differencebetween right-hand and left-hand molecules, due to parity nonconservationin weak interactions of electrons with nuclei,” JETP Lett. , 94–97 (1977). B. Y. Zel’dovich, D. B. Saakyan, and I. I. Sobel’man, “Energy differencebetween right-hand and left-hand molecules, due to parity nonconservationin weak interactions of electrons with nuclei,” Pis’ma Zhurnal Eksperimen-talnoi I Teoreticheskoi Fiziki , 106–109 (1977). R. A. Hegstrom, D. W. Rein, and P. G. H. Sandars, “Calculation of theparity nonconserving energy difference between mirror-image molecules,”J. Chem. Phys. , 2329–2341 (1980). I. B. Khriplovich, “On the energy difference between optical isomers re-sulting from parity nonconservation,” Sov. Phys. JETP , 177–183 (1980). V. G. Gorshkov, M. G. Kozlov, and L. N. Labzovsky, “ P -odd effects inpolyatomic molecules,” Zhurnal Eksperimentalnoi I Teoreticheskoi Fiziki , 1807–1819 (1982). V. G. Gorshkov, M. G. Kozlov, and L. N. Labzowsky, “ P -odd effects inpolyatomic molecules,” Sov. Phys. JETP , 1042–1048 (1982). M. Quack, “On the measurement of the parity violating energy differencebetween enantiomers,” Chem. Phys. Lett. , 147–153 (1986). R. N. Compton and R. M. Pagni, “The chirality of biomolecules,” AdvancesIn Atomic, Molecular, And Optical Physics, Vol 48 , 219–261 (2002). O. N. Kompanets, A. R. Kukudzhanov, V. S. Letokhov, and L. L.Gervits, “Narrow resonances of saturated absorption of the asymmetricalmolecule CHFClBr and the possibility of weak current detection in molec-ular physics,” Opt. Commun. , 414–416 (1976). E. Arimondo, P. Glorieux, and T. Oka, “Observation of inverted in-frared lamb dips in separated optical isomers,” Opt. Commun. , 369–372(1977). A. Bauder, A. Beil, D. Luckhaus, F. Müller, and M. Quack, “Com-bined high resolution infrared and microwave study of bromochlorofluo-romethane,” J. Chem. Phys. , 7558–7570 (1997). C. Daussy, T. Marrel, A. Amy-Klein, C. T. Nguyen, C. J. Bordé, andC. Chardonnet, “Limit on the parity nonconserving energy difference be-tween the enantiomers of a chiral molecule by laser spectroscopy,” Phys.Rev. Lett. , 1554–1557 (1999). R. A. Harris and L. Stodolski, “The effect of the parity violating electron-nucleus interaction on the spin-spin coupling Hamiltonian of chiralmolecules,” J. Chem. Phys. , 3862–3863 (1980). M. Schnell and J. Küpper, “Tailored molecular samples for precision spec-troscopy experiments,” Faraday Disc. , 33–49 (2011). A. L. Barra, J. B. Robert, and L. Wiesenfeld, “Parity non-conservationand NMR observables. Calculation of Tl resonance frequency differencesin enantiomers,” Phys. Lett. A , 443–447 (1986). A. L. Barra, J. B. Robert, and L. Wiesenfeld, “Possible observation ofparity nonconservation by high-resolution NMR,” Europhys. Lett. , 217–222 (1988). A. L. Barra and J. B. Robert, “Parity non-conservation and NMR parame-ters,” Mol. Phys. , 875–886 (1996). J. Eills, J. W. Blanchard, L. Bougas, M. G. Kozlov, A. Pines, and D. Bud-ker, “Measuring molecular parity nonconservation using nuclear-magnetic-resonance spectroscopy,” Phys. Rev. A , 042119 (2017). R. A. Harris and L. Stodolski, “Quantum beats in optical activity and weakinteractions,” Phys. Lett. B , 313–317 (1978). R. A. Harris and L. Stodolski, “On the time dependence of optical activity,”J. Chem. Phys. , 2145–2155 (1981). R. Berger, “Molecular parity violation in electronically excited states,”Phys. Chem. Chem. Phys. , 12–17 (2003). M. Quack, “Structure and dynamics of chiral molecules,” Angew. Chem.Int. Ed. , 571–586 (1989). M. Quack, “How important is parity violation for molecular and biomolec-ular chirality?” Angew. Chem. Int. Ed. , 4618–4630 (2002). R. Berger, “Parity-violation effects in molecules,” in
Relativistic ElectronicStructure Theory, Part: 2, Applications , edited by P. Schwerdtfeger (Else-vier, Netherlands, 2004) Chap. 4, pp. 188–288. J. Crassous, C. Chardonnet, T. Saue, and P. Schwerdtfeger, “Recent ex-perimental and theoretical developments towards the observation of parityviolation (pv) effects in molecules by spectroscopy,” Org. Biomol. Chem. , 2218–2224 (2005). M. Quack, J. Stohner, and M. Willeke, “High-resolution spectroscopicstudies and theory of parity violation in chiral molecules,” Annu. Rev. Phys.Chem. , 741–769 (2008). P. Schwerdtfeger, “The search for parity violation in chiral molecules,”in
Computational Spectroscopy: Methods, Experiments and Applications ,edited by J. Grunenberg (Wiley, Netherlands, 2010) Chap. 7, pp. 201–221. R. Berger and J. Stohner, “Parity violation,” Wiley In-terdiscip. Rev.-Comput. Mol. Sci. , e1396 (2019),https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1396. M. Quack and J. Stohner, “Influence of parity violating weak nuclear po-tentials on vibrational and rotational frequencies in chiral molecules,” Phys.Rev. Lett. , 3807–3810 (2000). J. K. Laerdahl and P. Schwerdtfeger and H. M. Quiney, “Theoretical analy-sis of parity-violating energy differences between the enantiomers of chiralmolecules,” Phys. Rev. Lett , 3811–3814 (2000). R. G. Viglione, R. Zanasi, P. Lazzeretti, and A. Ligabue, “Theoretical deter-mination of parity-violating vibrational frequency differences between theenantiomers of the CHFClBr molecule,” Phys. Rev. A , 052516 (2000). M. Quack and J. Stohner, “How do parity violating weak nuclear inter-actions influence rovibrational frequencies in chiral molecules?” Z. Phys.Chem. , 675–703 (2000). P. Schwerdtfeger, J. K. Laerdahl, and C. Chardonnet, “Calculation ofparity-violating effects for the C-F stretching mode of chiral methyl flu-orides,” Phys. Rev. A , 042508 (2002). M. Quack and J. Stohner, “Parity violation in chiral molecules,” Chimia ,530–538 (2005). R. Berger and J. L. Stuber, “Electroweak interactions in chiral molecules:Two-component density functional theory study of vibrational frequencyshifts in polyhalomethanes,” Mol. Phys. , 41–49 (2007). C. Thierfelder, G. Rauhut, and P. Schwerdtfeger, “Relativistic coupled-cluster study of the parity-violation energy shift of chfclbr,” Phys. Rev. A , 032513 (2010). M. Ziskind, C. Daussy, T. Marrel, and C. Chardonnet, “Improved sensitiv-ity in the search for a parity-violating energy difference in the vibrationalspectrum of the enantiomers of CHFClBr,” Eur. Phys. J. D , 219–225(2002). B. Darquie, C. Stoeffler, A. Shelkovnikov, C. Daussy, A. Amy-Klein,C. Chardonnet, S. Zrig, L. Guy, J. Crassous, P. Soulard, P. Asselin, T. R.Huet, P. Schwerdtfeger, R. Bast, and T. Saue, “Progress Toward the FirstObservation of Parity Violation in Chiral Molecules by High-ResolutionLaser Spectroscopy,” Chirality , 870–884 (2010). A. Cournol, M. Manceau, M. Pierens, L. Lecordier, D. B. A. Tran, R. San-tagata, B. Argence, A. Goncharov, O. Lopez, M. Abgrall, Y. L. Coq, R. L.Targat, H. A. Martinez, W. K. Lee, D. Xu, P. E. Pottie, R. J. Hendricks, T. E.Wall, J. M. Bieniewska, B. E. Sauer, M. R. Tarbutt, A. Amy-Klein, S. K.Tokunaga, and B. Darquié, “A new experiment to test parity symmetry incold chiral molecules using vibrational spectroscopy,” Quantum Electron. , 288–292 (2019). F. Faglioni and P. Lazzeretti, “Parity violation effect on vibrational spectra,”Phys. Rev. A , 032101 (2003). P. Schwerdtfeger, J. Gierlich, and T. Bollwein, “Large parity-violation ef-fects in heavy-metal-containing chiral compounds,” Angew. Chem. Int. Ed. , 1293–1296 (2003). R. Bast and P. Schwerdtfeger, “Parity-violation effects in the C-F stretchingmode of heavy-atom methyl fluorides,” Phys. Rev. Lett. , 023001 (2003). A. Messiah,
Quantenmechanik , Vol. 1 (Walter de Gruyter, Berlin, 1976). T. Marrel, M. Ziskind, C. Daussy, and C. Chardonnet, “High precisionrovibrational and hyperfine analysis of the ν = , 195–209 (2001). J. K. Laerdahl, R. Wesendrup, and P. Schwerdtfeger, “Parity-violatinginteractions and biochemical homochirality,” ChemPhysChem , 60–62(2000). M. Quack and J. Stohner, “Molecular chirality and the fundamental symme-tries of physics: Influence of parity violation on rovibrational frequenciesand thermodynamic properties,” Chirality , 745–753 (2001). P. Schwerdtfeger, T. Saue, J. N. P. van Stralen, and L. Visscher, “Relativis-tic second-order many-body and density-functional theory for the parity-violation contribution to the C-F stretching mode in CHFClBr,” Phys. Rev.A , 012103 (2005). M. Quack and J. Stohner, “Combined multidimensional anharmonic andparity violating effects in cdbrclf,” J. Chem. Phys. , 11228–11240(2003). G. Rauhut, “Efficient calculation of potential energy surfaces for the gen-eration of vibrational wave functions,” J. Chem. Phys. , 9313–9322(2004), https://doi.org/10.1063/1.1804174. A. D. Buckingham and W. Urland, “Isotope effects onmolecular properties,” Chem. Rev. , 113–117 (1975),https://doi.org/10.1021/cr60293a005. D. Jayatilaka, P. E. Maslen, R. D. Amos, and N. C. Handy,“Higher analytic derivatives,” Mol. Phys. , 271–291 (1992),https://doi.org/10.1080/00268979200100221. R. Bast, U. Ekström, B. Gao, T. Helgaker, K. Ruud, and A. J. Thorvaldsen,“The ab initio calculation of molecular electric, magnetic and geometricproperties,” Phys. Chem. Chem. Phys. , 2627–2651 (2011). T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen, andK. Ruud, “Recent advances in wave function-based methods of molecular-property calculations,” Chem. Rev. , 543–631 (2012), pMID: 22236047,https://doi.org/10.1021/cr2002239. H. Sellers, “Variational energy derivatives and perturba-tion theory,” Int. J. Quantum Chem. , 271–277 (1988),https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.560330403. C. Chang, M. Pelissier, and P. Durand, “Regular two-component Pauli-likeeffective Hamiltonians in Dirac theory,” Phys. Scr. , 394–404 (1986). E. van Lenthe, J. G. Snijders, and E.-J. Baerends, “The zero-order regularapproximation for relativistic effects: The effect of the spin-orbit couplingin closed shell molecules,” J. Chem. Phys. , 6505 (1996). C. van Wüllen, “Molecular density functional calculations in the regularrelativistic approximation: Method, application to coinage metal diatomics,hydrides, fluorides and chlorides, and comparison with first-order relativis-tic calculations,” J. Chem. Phys. , 392–399 (1998). C. van Wüllen, “A Quasirelativistic Two-component Density Functionaland Hartree-Fock Program,” Z. Phys. Chem , 413–426 (2010). R. Berger, N. Langermann, and C. van Wüllen, “Zeroth order regularapproximation approach to molecular parity violation,” Phys. Rev. A ,042105 (2005). C. van Wüllen and N. Langermann, “Gradients for two-componentquasirelativistic methods. application to dihalogenides of element 116,” J.Chem. Phys. , 114106 (2007), https://doi.org/10.1063/1.2711197. O. Treutler and R. Ahlrichs, “Efficient molecular numerical integrationschemes,” J. Chem. Phys. , 346 (1995). A. D. Becke, “Density-functional exchange-energy approximation withcorrect asymptotic-behavior,” Phys. Rev. A , 3098–3100 (1988). M. Ringholm, D. Jonsson, and K. Ruud, “A general, recursive, and open-ended response code,” J. Comput. Chem. , 622–633 (2014). J. Olsen, D. L. Yeager, and P. Jørgensen, “Triplet excitation properties inlarge scale multiconfiguration linear response calculations,” J. Chem. Phys. , 381–388 (1989). P. Sałek, O. Vahtras, T. Helgaker, and H. Ågren, “Density-functional the-ory of linear and nonlinear time-dependent molecular properties,” J. Chem.Phys. , 9630–9645 (2002), https://doi.org/10.1063/1.1516805. T. Saue and H. J. A. Jensen, “Linear response at the 4-component rela-tivistic level: Application to the frequency-dependent dipole polarizabil-ities of the coinage metal dimers,” J. Chem. Phys. , 522–536 (2003),https://doi.org/10.1063/1.1522407. S. Nahrwold and R. Berger, “Zeroth order regular approximation approachto parity violating nuclear magnetic resonance shielding tensors,” J. Chem.Phys. , 214101 (2009). T. A. Isaev and R. Berger, “Electron correlation and nuclear charge de-pendence of parity-violating properties in open-shell diatomic molecules,”Phys. Rev. A , 062515 (2012). M. Häser and R. Ahlrichs, “Improvements on the direct SCF method,” J.Comput. Chem. , 104–111 (1989). R. Ahlrichs, M. Bär, M. Häser, H. Horn, and C. Kölmel, “Electronic struc-ture calculations on workstation computers: The program system turbo-mole,” Chem. Phys. Lett. , 165–169 (1989). MATLAB, (The MathWorks Inc., Natick, Mas-sachusetts, 2018). K. Gaul and R. Berger, “Toolbox approach for quasi-relativistic calcula-tion of molecular properties for precision tests of fundamental physics,” J.Chem. Phys. , 044101 (2020), arXiv:1907.10432 [physics.chem-ph]. A. Bergner, M. Dolg, W. Küchle, H. Stoll, and H. Preuß, “
Ab initio energy-adjusted pseudopotentials for elements of groups 13-17,” Mol. Phys. ,1431–1441 (1993). W. Küchle, M. Dolg, H. Stoll, and H. Preuss, “
Ab initio pseudopotentialsfor Hg through Rn I. parameter sets and atomic calculations,” Mol. Phys. , 1245–1263 (1991). H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, “notitle,” WIREs Comput. Mol. Sci. , 242–253 (2012). Q. Ma and H.-J. Werner, “Explicitly correlated local coupled-cluster meth-ods using pair natural orbitals,” WIREs Comput. Mol. Sci , e1371 (2018),https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1371. H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßel-mann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, T. F. Miller,A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut, and M. Sibaev,“The molpro quantum chemistry package,” J. Chem. Phys. , 144107(2020), https://doi.org/10.1063/5.0005081. J. K. Laerdahl and P. Schwerdtfeger, “Fully relativistic ab initio calculationof the energies of chiral molecules including parity-violating weak interac-tions,” Phys. Rev. A , 4439–4453 (1999). P. A. M. Dirac, “Note on exchange phenomena in the Thomas atom,” Proc.Cambridge Phil. Soc. , 376–385 (1930). S. H. Vosko, L. Wilk, and M. Nuisar, “Accurate spin-dependent electronliquid correlation energies for local spin density calculations: A criticalanalysis,” Can. J. Phys. , 1200–1211 (1980). K. Gaul, M. G. Kozlov, T. A. Isaev, and R. Berger, “Chiral molecules assensitive probes for direct detection of P -odd cosmic fields,” Phys. Rev.Lett. , 123004 (2020), arXiv:2005.02429 [hep-ph]. K. Gaul, M. G. Kozlov, T. A. Isaev, and R. Berger, “Parity nonconservinginteractions of electrons in chiral molecules with cosmic fields,” Phys. Rev.A , 032816 (2020), arXiv:2005.03938 [physics.chem-ph]. K. Gaul and R. Berger, “Quasi-relativistic study of nuclear elec-tric quadrupole coupling constants in chiral molecules contain-ing heavy elements,” Molecular Physics , e1797199 (2020),https://doi.org/10.1080/00268976.2020.1797199. G. Rauhut, V. Barone, and P. Schwerdtfeger, “Vibrational analyses forCHFClBr and CDFClBr based on high level ab initio calculations,” J.Chem. Phys. , 054308 (2006), https://doi.org/10.1063/1.2236112. R. F. Garcia Ruiz, R. Berger, J. Billowes, C. L. Binnersley, M. L. Bissell,A. A. Breier, A. J. Brinson, K. Chrysalidis, T. E. Cocolios, B. S. Cooper,K. T. Flanagan, T. F. Giesen, R. P. de Groote, S. Franchoo, F. P. Gustafsson,T. A. Isaev, Á. Koszorús, G. Neyens, H. A. Perrett, C. M. Ricketts, S. Rothe,L. Schweikhard, A. R. Vernon, K. D. A. Wendt, F. Wienholtz, S. G. Wilkins,and X. F. Yang, “Spectroscopy of short-lived radioactive molecules,” Nature581