Evaluation of Coupled Perturbed and Density Functional Methods of Computing the Parity-Violating Energy Difference between Enantiomers

Abstract

We present new coupled-perturbed Hartree-Fock (CPHF) and density functional theory (DFT) computations of the parity-violating energy difference (PVED) between enantiomers for H₂O₂ and H₂S₂. Our DFT PVED computations are the first for H₂S₂ and the first with the new HCTH and OLYP functionals. Like other “second generation” PVED computations, our results are an order of magnitude larger than the original “first generation” uncoupled-perturbed Hartree-Fock computations of Mason and Tranter. We offer an explanation for the dramatically larger size in terms of cancellation of contributions of opposing signs, which also explains the basis set sensitivity of the PVED, and its conformational hypersensitivity (addressed in the following paper). This paper also serves as a review of the different types of “second generation” PVED computations: we set our work in context, comparing our results with those of four other groups, and noting the good agreement between results obtained by very different methods. DFT PVEDs tend to be somewhat inflated compared to the CPHF values, but this is not a problem when only sign and order of magnitude are required. Our results with the new OLYP functional are less inflated than those with other functionals, and OLYP is also more efficient computationally. We therefore conclude that DFT computation offers a promising approach for low-cost extension to larger biosystems, especially polymers. The following two papers extend to terrestrial and extra-terrestrial amino acids respectively, and later work will extend to polymers.

Introduction

The “PVED” is the tiny parity-violating energy difference between enantiomers produced by the weak neutral current mediated by the Z⁰ boson. The existence of the PVED was first recognized by Yamagata (1966), who proposed that this small energy difference may have determinately selected life’s use of the natural L-amino acids rather than their “unnatural” D mirror images.

The PVED can be calculated using the theoretical framework framework (“Theoretical Framework of PVED Calculations”) originally developed by Hegstrom et al. (1979). Following the first calculations by Hegstrom et al. (1980), Mason, Tranter and MacDermott (Mason and Tranter 1984, 1985; Tranter 1985, 1986a, b; MacDermott 1989; MacDermott 1995, 2000) carried out a “first generation” of ab initio computations of the PVED of biorelevant molecules at the Uncoupled-Perturbed Hartree-Fock level, which we briefly review in ““First Generation” PVED Computations”. They concluded that many natural enantiomers, including the L-amino acids and the D-nucleotide units of right-hand B DNA, were more stable than their “unnatural” mirror images by about 10^− 20 hartree (or 10^− 17 kT at room temperature).

More recently, several more groups entered the field with a “second generation” of PVED computations using a sometimes bewildering variety of more sophisticated methods, all of which yielded PVEDs an order of magnitude larger than those of the first generation. We review these methods in ““Second Generation” PVED Computations” and “Comparison with other Second Generation Computations”, clarifying the differences between them.

In this paper, we present new coupled-perturbed and density functional computations of the PVED. Use of density functional theory makes these computations more efficient and less expensive than other second generation computations, and thus better suited to extension to larger biomolecules and polymers. “New Coupled-Perturbed Hartree-Fock Computations of the PVES” and “New Density Functional Computations of the PVES” describe our methods and their application to the test molecules H₂O₂ and H₂S₂, and in “Comparison with other Second Generation Computations” we compare our PVED results for these molecules with other second generation computations by other groups.

Theoretical Framework of PVED Calculations

We now summarize the theory of the PVED, so as to clarify the approximations used in our approach.

The PVED between enantiomers arises predominantly from weak neutral current (WNC) interactions between electrons and nuclei; the effect of electron-electron WNC interactions is estimated to be about two orders of magnitude smaller, and is therefore neglected. The electron-nucleus WNC interactions impart a parity-violating energy shift (PVES) of + E _pv to one enantiomer and an equal and opposite PVES of − E _pv to the other enantiomer, giving overall a parity-violating energy difference (PVED) of ΔE _pv  = 2E _pv .

The Hamiltonian for the electron-nucleus WNC interaction (Hegstrom et al. 1979; Mason and Tranter 1984) is

$$ \hat H_{pv} = \left({\frac{G_F}{2\sqrt{2}}}\right) \left({\frac{1}{2m_ec}}\right) \sum_{i}\sum_{A} Q_A^c \left\{\mathbf{\sigma}_i\cdot\mathbf{p}_i, \delta^3(\mathbf{r}_i - \mathbf{r}_A)\right\}_+ $$

(1)

or (since ${\mathbf{s}} = \hbar {\mathbf{\sigma}}/2$)

$$ \hat H_{pv} = \Gamma \sum_{i}\sum_{A} Q_A^c \left\{\mathbf{s}_i\cdot\mathbf{p}_i, \delta^3(\mathbf{r}_i - \mathbf{r}_A)\right\}_+ $$

(2)

where

$$ \Gamma = {\frac{G_F}{2\sqrt{2}m_ec\hbar}} $$

(3)

and s _i , p _i are the spin and momentum of electron i, the summations are over all electrons i and nuclei A, the anti-commutator ensures Hermiticity, and the delta function expresses the contact nature of the weak interaction. The electroweak coupling constant G _F has the value 1.435 ×10^− 62 Jm³. The “electroweak charge” $Q_A^c$ of nucleus A is given by

$$ Q_A^c = Z_{\!A} \left(1 - 4 \sin^2\Theta_W\right) - N_A $$

(4)

where Z _A and N _A are respectively the atomic number and neutron number of nucleus A, and Θ _W is the Weinberg angle. In common with most other computations, we use the theoretical value $\Theta_W = 30^\circ$ or $\sin^2\Theta_W = 0.25$, giving $Q_A^c = -N_A$.

It is convenient to work in atomic units, in which

$$ \hbar = m_e = e = 4\pi\epsilon_0=1, \quad\quad c = 1/\alpha $$

(5)

(where α is the fine structure constant), so that G _F becomes 2.222×10^− 14 atomic units, Γ is now given by

$$ \Gamma = {\frac{G_F\alpha}{2\sqrt{2}}} = 5.73395 \times 10^{-17} {\mathrm{atomic ~units}} $$

(6)

and energies are expressed in hartrees, where

$${\mathrm{E_h}} = 1~ {\mathrm{hartree}} = 4.35975 \times 10^{-18} {\mathrm{J}} \equiv 2625.50\enspace {\mathrm {kJ mol}}^{-1} $$

(7)

A hartree is thus about 1000 kT at room temperature.

The PVES is then given by

$$ E_{pv} = \big\langle\Psi_0|\hat H_{pv}|\Psi_0\big\rangle $$

(8)

where Ψ₀ is the multi-electron molecular wavefunction. But because $\hat H_{pv}$ is pure imaginary and E _pv must be real, E _pv is identically zero unless a spin-orbit coupling correction (Hegstrom et al. 1979) is used to give

$$ E_{pv} = \left\langle\Psi_0^{\mathrm{corr}}|\hat H_{pv}|\Psi_0^{\mathrm{corr}}\right\rangle $$

(9)

where the corrected wavefunction is

$$ \left|\Psi_0^{\mathrm{corr}}\right\rangle = |\Psi_0\rangle + \sum_{F \ne 0} {\frac{\big\langle\Psi_F|\hat H_{so}|\Psi_0\big\rangle}{E_0 - E_F}}|\Psi_F\rangle $$

(10)

and the spin-orbit coupling Hamiltonian is

$$ \hat H_{so} = \sum_{i}\sum_{A}\xi_{iA}~\mathbf{l}_{iA}\cdot\mathbf{s}_i $$

(11)

where l _iA  = r _iA  ∧ p _i is the orbital angular momentum of electron i with respect to nucleus A. The spin-orbit coupling constant is given by

$$ \xi_{iA} = {\frac{Z^{\mathrm{eff}}_Ae^2}{8\pi\epsilon_0m_e^2c^2r_{iA}^3}} $$

(12)

(see e.g. Atkins and Friedmann 2005) so in atomic units $\hat H_{so}$ is given by

$$ \hat H_{so} = {\frac{\alpha^2}{2}}\sum_{i}\sum_{A}Z_{\!A}^{\mathrm{eff}}{{\mathbf{l}_{iA}\cdot\mathbf{s}_i}\over{r_{iA}^3}} \ $$

(13)

Two-electron spin-other-orbit interactions are omitted in our approach. This is a valid approximation only because the effective nuclear charge $Z_{\!A}^{\mathrm{eff}}$, determined from Slater’s rules, is used instead of the bare nuclear charge Z _A .

The PVES then becomes

$$ E_{pv} = 2 \sum_{F \ne 0} {\mathrm{Re}} \left\{ {\frac{\big\langle\Psi_0|\hat H_{pv}|\Psi_F\big\rangle\big\langle\Psi_F|\hat H_{so}|\Psi_0\big\rangle}{E_0 - E_F}} \right\} \label{eq14} $$

(14)

Most computations—except the CI computations of Quack and co-workers (Bakasov et al. 1998; Bakasov and Quack 1999; Berger and Quack 2000)—approximate both the ground state Ψ₀ and the excited state Ψ _F as single determinants, with Ψ _F restricted to single excited relative to Ψ₀. By writing the multi-electron wavefunctions $|\Psi_0\rangle$ and $|\Psi_F\rangle$ as simple products of one-electron molecular spin-orbitals $|\psi_j m_s\rangle$, and then factorizing into spin $|m_s\rangle$ and orbital $|\psi_j\rangle$ parts and summing over spins, the spin-separated form of Eq. 14, suitable for restricted Hartree-Fock computations, then becomes

$$ E_{pv} = \sum_{j \in \mathrm{occ}}\sum_{k\in\mathrm{unocc}}{\frac{\big\langle\psi_j\big|\hat V_{\!pv}\big|\psi_k\big\rangle \cdot\big\langle\psi_k\big|\hat V_{\!so}\big|\psi_j\big\rangle} {E_0 - E_{F(j\rightarrow k)}}}\label{eq15} $$

(15)

where ψ _j and ψ _k are spatial MOs that are respectively occupied and unoccupied in Ψ₀, E _F(j→k) is the energy of the state Ψ _F formed by single-electron promotion from ψ _j to ψ _k , and the spin-free one-electron operators are given by

$$ \hat V_{\!pv}(i) = {{G_F\alpha}\over{2\sqrt{2}}}\sum_{A} Q_W^c(A)\left\{\mathbf{p}_i, \delta^3(\mathbf{r}_{iA})\right\}_+\label{eq16} $$

(16)

$$ \hat V_{\!so}(i) = {\frac{\alpha^2}{2}}\sum_{A}Z_{\!A}^{\mathrm{eff}}{\frac{\mathbf{l}_{iA}}{r_{iA}^3}}\label{eq17} $$

(17)

The MOs {ψ _a } are then expanded in terms of an AO basis set {ϕ _μ },

$$ \psi_a = \sum_{\mu}c_{\mu a}\phi_{\mu} \label{eq18} $$

(18)

where c _μa are the MO coefficients. Computation of the PVES then reduces to determining the matrix elements of $\hat V_{\!pv}$ and $\hat V_{\!so}$ in the AO basis. The single centre approximation (Mason and Tranter 1984) is then invoked, meaning that these integrals are taken to be non-zero only between atomic orbitals on the same atom. The integrals are then converted to the MO basis for use in Eq. 15.

“First Generation” PVED Computations

The earliest computations of Mason and Tranter (1984) approximated the energy difference between the multi-electron ground state Ψ₀ and the excited state Ψ _F(j→k) by the energy difference between the MOs ψ _j and ψ _k :

$$ E_0 - E_{F(j\rightarrow k)} \approx \varepsilon_j - \varepsilon_k\label{eq19} $$

(19)

This is the Uncoupled-Perturbed Hartree-Fock (UPHF) approximation, which assumes that the energies of the other orbitals do not change when an electron is excited from ψ _j to ψ _k . First generation UPHF computations seemed to show that the “natural” enantiomers are indeed more stable in most cases, e. g.  several L-amino acids (Mason and Tranter 1984, 1985; Tranter 1985) (in the α-helix and β-sheet conformations, and in what was then thought to be the solution conformation), D-glyceraldehyde (Tranter 1986b), β-D-deoxyribose (Tranter et al. 1992), and short fragments of B DNA and A RNA (MacDermott 1995), were all found to be PVED-stabilized by about 10^− 20 hartree per molecule (or per monomer unit within a polymer fragment). Much larger PVEDs (10^− 19 to 10^− 17 hartree per monomer unit) were obtained for possible primitive replicators in which the normal phosphate links of DNA are modified by thiosubstitution (MacDermott et al. 1989). Kondepudi (1987) devised an amplification mechanism, based on autocatalysis and enantiomeric antagonism, which could amplify a PVED of 10^− 20 hartree (10^− 17 kT at room temperature) to homochirality in only 10⁴ years—a very short time on an evolutionary timescale. The much larger PVEDs of thiosubstituted DNA analogues (which may be relevant to the possible origin of life in the sulphurous environment of deep-sea volcanic vents) should therefore be amplifiable on a much shorter timescale.

“Second Generation” PVED Computations

The two main deficiencies of the first generation computations are (a) their unsatisfactory treatment of the energy denominator E ₀ − E _F , which is set equal to ε _j  − ε _k , thus neglecting any perturbative coupling, and (b) the use of single determinants for the excited state Ψ _F , thus neglecting electron correlation effects.

Since 1996, several groups have produced a “second generation” of PVED computations that seek to remedy these deficiencies with a variety of more sophisticated methods. Lazzeretti and co-workers introduced perturbative coupling with their “RPA” method (Lazzeretti and Zanasi 1997; Zanasi and Lazzeretti 1998; Zanasi et al. 1999), which is essentially equivalent to our Coupled-Perturbed Hartree-Fock (CPHF) method presented in “New Coupled-Perturbed Hartree-Fock Computations of the PVES”. Quack and co-workers used an improved treatment of the excited states Ψ _F —the Configuration Interaction Singles Restricted Hartree-Fock (CIS-RHF) method—to introduce electron correlation (Bakasov et al. 1998; Bakasov and Quack 1999). They later used a more general Multiconfiguration Linear Response (MC-LR) method (Berger and Quack 2000) that encompasses the effects of both CI electron correlation and perturbative coupling. Schwerdtfeger and co-workers carried out fully relativistic Dirac Hartree-Fock (DHF) computations (Laerdahl and Schwerdtfeger 1999; Wesendrup et al. 2003). Helgaker and co-workers (Hennum et al. 2002) carried out the first ever density functional theory (DFT) computations of the PVED. Their study focuses on H₂O₂ only, and uses mainly the B3LYP functional, but in this paper we extend for the first time to H₂S₂ and to the newer HCTH and OLYP functionals.

All four groups of second generation calculations—“RPA”, CI/MC-LR, DHF and DFT—produce PVEDs that are consistently an order of magnitude larger than the corresponding first generation calculations (although almost always of the same sign), and also in remarkably close agreement with each other despite being based on very different approaches.

In the following sections we describe our own “second generation” improvements to the earlier “first generation” UPHF treatment: in “New Coupled-Perturbed Hartree-Fock Computations of the PVES” we use Coupled-Perturbed Hartree-Fock (CPHF) to improve treatment of the energy denominator E ₀ − E _F(j→k) in expression (15) for the PVES, and in “New Density Functional Computations of the PVES” we use Density Functional Theory (DFT) to recover the effects of electron correlation.

New Coupled-Perturbed Hartree-Fock Computations of the PVES

The effect of CPHF on the expression for the PVES is best seen by analogy with UPHF theory. Expression (15) is first written in the form

$$ E_{pv} = \sum_{j \in {\mathrm{occ}}} \sum_{k\in{\mathrm{unocc}}} U_{jk}\cdot\big\langle\psi_k|\hat V_{\!so}|\psi_j\big\rangle \label{eq20} $$

(20)

where, in UPHF theory, U _jk is given by

$$ U_{jk}^{\mathrm {UPHF}} = {\frac{\big\langle\psi_j|\hat V_{\!pv}|\psi_k\big\rangle}{\varepsilon_j - \varepsilon_k}}\label{eq21} $$

(21)

Equivalently, U may be considered as a solution to the matrix equation

$$ AU = B \label{eq22} $$

(22)

where

$$ B_{jk} = \big\langle\psi_j|\hat V_{\!pv}|\psi_k\big\rangle \label{eq23} $$

(23)

$$ A_{jk}^{\mathrm{UPHF}} = \varepsilon_j - \varepsilon_k \label{eq24} $$

(24)

and

$$ U = A^{-1}B \label{eq25} $$

(25)

The uncoupled nature of the UPHF approximation is evident from the fact that $U_{jk}^{\mathrm{UPHF}}$ may be determined with reference only to the energy (Eq. 24) and integral (Eq. 23) of the j→k transition: each $U_{jk}^{\mathrm{UPHF}}$ element may be determined independently for a given j and k. But in CPHF, Eq. 22 is solved as a matrix equation, i.e. a set of simultaneous equations (in contrast to UPHF, in which Eq. 22 may be separated into a set of independent equations). In CPHF, the form of the A matrix depends on whether the perturbation matrix B contains a real or imaginary operator. In this case, $\hat V_{\!pv}$ is imaginary, and the A matrix is given (Amos et al. 1987) by

$$ A_{jk}^{\mathrm{CPHF, Im}} = \sum_{l\in{\mathrm{occ}}}\sum_{m\in{\mathrm{unocc}}}\{[\varepsilon_j - \varepsilon_k]\delta_{jl}\delta_{km} - [(jm|kl) - (jl|km)]\} \label{eq26} $$

(26)

where the following shorthand is used for the two-electron integrals

$$ \left(pq|rs\right) = \int d{\mathbf{r}}_1d{\mathbf{r}}_2\psi_p^*({\mathbf{r}}_1)\psi_q({\mathbf{r}}_1)r_{12}^{-1}\psi_r^*({\mathbf{r}}_2)\psi_s({\mathbf{r}}_2) \label{eq27} $$

(27)

The perturbative coupling that is central to CPHF theory is evident from the fact that to calculate $U_{jk}^{\mathrm{CPHF}}$, a knowledge of all occupied MOs {ψ _l }, and all unoccupied MOs {ψ _m } is required, whereas in UPHF only the occupied MO ψ _j and the unoccupied MO ψ _k are required to calculate $U_{jk}^{\mathrm{UPHF}}$.

We have implemented both UPHF and CPHF computations of the PVED by adding new codes to the Cambridge Analytic Derivatives Package (CADPAC) (Amos et al. 1998). CADPAC already contains routines to carry out SCF computations and to solve the CPHF equations (Eq. 25) for a supplied set of perturbative matrix elements B, so computation of the PVES from Eq. 20 reduces to evaluating the matrix elements of $\hat V_{\!pv}$ and $\hat V_{\!so}$. As in earlier computations, the single centre approximation (see ““First Generation” PVED Computations”) is used for these matrix elements. The matrix elements of $\hat V_{\!pv}$ and $\hat V_{\!so}$ are computed in the AO basis set {ϕ _μ } and then converted to the MO basis {ψ _a }. The only non-zero matrix elements in $\hat V_{\!pv}$ in the AO basis are those between s and p orbitals (due to the delta function), and the only non-zero matrix elements of $\hat V_{\!so}$ in a gaussian basis of s, p, d and f-type orbitals are those between two p-orbitals, two d-orbitals, two f-orbitals, or a p-orbital and an f-orbital. The value of the effective nuclear charge Z ^eff used for the matrix element $\langle\phi_\mu|\hat V_{\!so}|\phi_\nu\rangle$ is $\sqrt{Z^{\mathrm{eff}}(\mu)Z^{\mathrm{eff}}(\nu)}$ where Z ^eff(μ), Z ^eff(ν) are the effective nuclear charges, determined from Slater’s rules, for electrons in the atomic orbitals ϕ _μ , ϕ _ν respectively.

The PVES involves a scalar product of matrix elements of the operators $\hat V_{\!pv}$ and $\hat V_{\!so}$. It is therefore possible, and sometimes useful, to write E _pv as the sum of “x, y and z components” (as in Bakasov and Quack (1999)), although the values of these individual components are merely functions of the orientation of the molecule with respect to the x, y and z axes in the input file.

The geometries of the test molecules H₂O₂ and H₂S₂ are the same as those used by Mason and Tranter (1984) and in subsequent second generation computations by other groups (Lazzeretti and Zanasi 1997; Bakasov et al. 1998): for H₂O₂, the O–O and H–O bondlengths are respectively 149.0pm and 97.0pm, and the H–O–O bond angle is 100.0°; for H₂S₂ the S–S and H–S bondlengths are respectively 205.5pm and 135.2pm, and the H–S–S bond angle is 92.0°. The S-enantiomers correspond to a positive H–O–O–H or H–S–S–H torsion angle ϕ (0 < ϕ < π), and the R enantiomers to a negative torsion angle (0 > ϕ > − π). The sign convention for the torsion angle is that of Eliel and Wilen (1994): in a molecule A–B–C–D, the torsion angle ϕ between the planes A–B–C and B–C–D is positive if, looking along the B–C axis from the B side, the turn from A to D is clockwise. Care must be taken because some ab initio packages (including, unfortunately, CADPAC) use the opposite convention (that ϕ is positive if the turn from A to D is anticlockwise); GAUSSIAN, however, uses the more standard convention described here. We use the standard (Eliel and Wilen 1994) convention throughout, taking care to reverse the sign of all torsion angles in the CADPAC input files.

Figure 1 shows our UPHF and CPHF results for the PVES of H₂O₂ as a function of torsion angle, using a 6-31G basis set. The most significant feature is that the CPHF results are approximately one order of magnitude larger than the UPHF results.

Fig. 1

The PVES of H₂O₂ (in units of 10^− 20 hartree) as a function of torsion angle ϕ at UPHF and CPHF levels of theory with a 6-31G basis set

The reasons behind the order of magnitude increase in CPHF PVES over UPHF PVES may not be immediately obvious, since introduction of perturbative coupling via Eq. 26 adds only a small correction to the uncoupled energy denominator (24). The explanation lies in the fact that the PVES is a sum of contributions from different MOs, different nuclei, and different q components for q = 1, 2, 3 (or x, y, z). The individual contributions are generally much larger than the overall PVES, which results from near cancellation of many large terms of opposing sign. Modest changes in individual contributions can therefore lead to much larger changes in the overall PVES. This is illustrated in Tables 1 and 2, which compare UPHF and CPHF results for the x, y and z-components of the PVES of H₂O₂ and H₂S₂ respectively: each q-component is increased on going from UPHF to CPHF, but the overall PVES is increased much more.

Table 1 The three q-components of the PVES of H₂O₂ (ϕ = 120°), computed at UPHF and CPHF levels of theory with a 6-31G basis set

Table 2 The three q-components of the PVES of H₂S₂ (ϕ = 60°), computed at UPHF and CPHF levels of theory with 3-21G, 4-31G and 6-31G basis sets

From the results for H₂S₂ in Table 2 it is clear that the CPHF PVES is dependent on the choice of basis set, and that smaller basis sets such as 3-21G and 4-31G are unsatisfactory. The first columns of Tables 3 and 4 show the CPHF PVES of H₂O₂ (ϕ = 120°) and H₂S₂ (ϕ = 60°) respectively for a larger range of basis sets. While there is some variation in the computed E _pv for different basis sets, it is clear that for H₂O₂ (ϕ = 120°), $E_{pv} \approx 35 \times 10^{-20}$ hartree (E_h) for all the larger basis sets (6-31G and larger); significant variation from this value is found only for the smaller basis sets (3-21G, 4-31G, etc.). The unreliability of smaller basis sets for PVES computations stems from the contact nature of the weak interaction (note the delta function in $\hat V_{\!pv}$): this requires good modelling of the MOs at the nucleus, which is not possible with basis sets that have too few s-type Gaussian primitives. PVES computations also require good modelling of the valence electrons, since these sample the chirality of the molecule, and we see that some increase in PVES occurs on going from 6-31G to 6-311G, DZ or TZ, although the effect is less pronounced for H₂S₂ than for H₂O₂, and is not as large as the increase in PVES achieved by improved modelling at the nuclei. It is clear that adding d-type polarization functions (as in 6-31G*, 6-31G** vs 6-31G, or DZP vs DZ) makes little difference: this is because (a) the matrix elements of $\hat V_{\!pv}$ between d-orbitals are zero, (b) in the matrix elements of $\hat V_{\!so}$ between d-orbitals, the Z ^eff value from Slater’s rules is very low, because d-orbitals are strongly shielded. The basis set dependence of the PVES, even among the larger basis sets (for which one might expect the basis set limit to be approached) may be a further consequence of E _pv being a sum of almost cancelling large terms of opposite sign. For many molecular observables, such as spectroscopic constants, this level of basis set dependence might be considered unsatisfactory. But for computations of the PVES with reference to the possible electroweak origin of biomolecular chirality, where sign and order of magnitude are more important than any claims of supposed chemical accuracy, the basis set dependence in Tables 3 and 4 for 6-31G and larger may be judged to be entirely satisfactory.

Table 3 PVES of H₂O₂ (ϕ = 120°) computed for various basis sets using CPHF, and DFT with the well-known BLYP and B3LYP functionals and the newer HCTH functional (Hamprecht et al. 1998); CPKS is used for the hybrid functional B3LYP, but UPKS is used for the non-hybrid BLYP and HCTH functionals

Table 4 PVES of H₂S₂ (ϕ = 60°) computed with various basis sets using CPHF, and DFT with the well-known B3LYP functional and the newer HCTH (Hamprecht et al. 1998) and OLYP (Handy and Cohen 2001; Lee et al. 1988) functionals; CPKS is used for the hybrid functional B3LYP, but UPKS is used for the non-hybrid HCTH and OLYP functionals

Figure 2 shows the ϕ-dependence of the CPHF PVES of H₂S₂ computed using 6-31G and some considerably larger basis sets. The close agreement demonstrates that there is probably little to be gained (for purposes of obtaining sign and order of magnitude estimate) from using basis sets larger than the 6-31G basis set that was used for most of the first generation UPHF computations (MacDermott 1995, 2000). However, this conclusion may or may not be general, and would need to be verified for each new class of molecules studied. In the following papers we show that the basis set dependence of the PVES of several amino acids follows a similar pattern to that seen here for H₂O₂ and H₂S₂, the results for 6-31G being little different from those using considerably larger basis sets.

Fig. 2

The CPHF PVES of H₂S₂ (in units of 10^− 20 hartree) as function of torsion angle ϕ with 6-31G, 6-31G** and DZ basis sets

Finally, it is noteworthy that although CPHF methods revise the magnitude of the PVES upwards by an order of magnitude as compared with UPHF, the sign of the PVES remains unchanged for the H₂O₂ and H₂S₂ test molecules studied here, and also the amino acids considered in the following paper.

New Density Functional Computations of the PVES

Our CPHF computations introduce perturbative coupling, but take no account of electron correlation. However, CADPAC has the ability to do density functional theory (DFT) computations by Kohn-Sham (KS) SCF as an alternative to regular Hartree-Fock (HF) SCF, and KS-SCF accounts for the effect of electron correlation, whereas HF-SCF does not. Our UPHF and CPHF computations use respectively Eqs. 15 and 20 to obtain the PVES using a wavefunction obtained from CADPAC’s HF-SCF procedure. We were therefore able to introduce electron correlation to both our uncoupled-perturbed and coupled-perturbed computations by using instead a wavefunction from CADPAC’s KS-SCF DFT procedure to give the PVES at an Uncoupled-Perturbed Kohn-Sham (UPKS) or Coupled-Perturbed Kohn-Sham (CPKS) level from Eqs. 15 and 20 respectively. However, the Kohn-Sham equations, unlike the Hartree-Fock equations, are not “coupled” unless one uses a “hybrid” exchange correlational functional (such as B3LYP or B3P86) that includes a semi-empirical mixing of HF exchange energy. In the absence of HF exchange energy, there can be no perturbative coupling, so CPKS is identical to UPKS for PVES computations with non-hybrid (GGA) functionals (such as LDA, BLYP, BP86, HCTH). UPHF thus has neither perturbative coupling nor any treatment of electron correlation; CPHF has perturbative coupling but no treatment of electron correlation; non-hybrid DFT (UPKS = CPKS) and UPKS hybrid-DFT have some treatment of electron correlation but no perturbative coupling; while CPKS hybrid-DFT has both electron correlation and perturbative coupling.

Table 5 shows the DFT PVES of H₂O₂ (ϕ = 120°), computed at both UPKS and CPKS levels with a range of functionals using a 6-31G basis set. The UPKS and CPKS results are identical with non-hybrid functionals, providing a useful check on the accuracy of our computer coding. The striking feature of the DFT results is that, like the CPHF results, they are an order of magnitude larger than the UPHF results. The DFT results are of the same sign and order of magnitude as the CPHF results. The non-hybrid DFT results (UKPS = CPKS) and the hybrid-CPKS results are slightly larger in magnitude than the corresponding CPHF results ($E_{pv} = 34.26 \times 10^{-20} {\mathrm{E_h}}$ for H₂O₂ (ϕ = 120°) at 6-31G CPHF-level). This is due to the well-known deficiency of DFT in underestimating perturbative transition energies (i.e. the orbital energy difference ε _j  − ε _k ). The notoriously poor LDA functional gives PVESs that are 30% larger than the CPHF PVESs, but the B3P86-CPKS and HCTH-(CPKS = UPKS) PVESs are only 14-18% larger than the CPHF PVES. This small overestimation is entirely acceptable within the sign and order of magnitude accuracy aimed for when computing the PVES of biomolecules in the context of a possible electroweak origin of chirality.

Table 5 DFT PVES of H₂O₂ (ϕ = 120°), computed at the uncoupled-perturbed and coupled-perturbed levels for a range of exchange-correlation functionals using the 6-31G basis set

Note that the UPKS hybrid-DFT results (using B3LYP or B3P86) are significantly smaller than the corresponding CPKS results. This is because hybrid functionals contain an admixture of HF exchange energy, and in a UPKS computation this corresponds to mixing in a fraction of the UPHF PVES (which underestimates E _pv by an order of magnitude), resulting in UP-hybrid-DFT results that are noticeably smaller, and in fact intermediate between the UPHF and CPHF values. It is thus clear that if a hybrid functional is used, it is essential to invoke perturbative coupling with a CPKS calculation. But if a non-hybrid functional is used, a UPKS calculation is sufficient—which is potentially a great advantage for larger molecules, as perturbative coupling is more computationally demanding (CADPAC, for example, can accommodate up to 250 basis functions in coupled-perturbed computations, but up to 1000 basis functions in uncoupled-perturbed computations).

Furthermore, we note that either electron correlation (through DFT) or perturbative coupling (through CPHF) results in an approximate order of magnitude increase in PVES over the UPHF results, but that treatment of both effects (through CP-hybrid-DFT) does not result in any extra increase in PVES. This phenomenon is also seen in the MC-LR results of Quack (see next section).

Tables 3 and 4 show the effects of varying basis set on the DFT PVES of H₂O₂ (ϕ = 120°) and H₂S₂ (ϕ = 60°) for different functionals; perturbative coupling is invoked for the hybrid functional B3LYP, but not for the non-hybrid BLYP, HCTH and OLYP functionals. The general trends in the variation of the DFT PVES with basis set parallel those described above for the CPHF PVES, but with the DFT PVES values consistently slightly larger than the CPHF values. The results from B3LYP, and more especially BLYP, show considerably more variation with basis set than the CPHF values; by contrast, the HCTH functional provides a narrower, more satisfactory distribution about the mean PVES for different basis sets, and the overestimation by DFT as compared with CPHF seems less serious for HCTH. A further advantage of HCTH over B3LYP is that it is a non-hybrid functional, so it is not necessary to invoke the computationally more demanding CPKS procedure: a less demanding UPKS computation gives the same result. Still better than HCTH (Hamprecht et al. 1998) is OLYP (Handy and Cohen 2001; Lee et al. 1988), a new non-hybrid functional which gives similar results to HCTH, but requires 3 to 4 times less CPU time. However, when basis sets with polarization are used (6-31G**, DZP, etc), some of the H₂S₂ DFT results (but not the H₂O₂ results) seem surprisingly overinflated (by as much as 70–80%) compared with the corresponding CPHF results. The cause of this overinflation is not immediately obvious, but is perhaps illustrative of the well-known tendency of DFT to give better results with some molecules, basis sets and functionals than with others. But we have already established that there is in general little advantage in adding polarization functions—and the H₂S₂ DFT results for basis sets without polarization functions (6-31G, DZ, etc) do not suffer from overinflation, and agree well with the CPHF results. However, it is clearly important to use DFT with caution: before proceeding with large numbers of computations for a particular class of molecules (e. g. in this case sulphur-containing molecules), the performance of different DFT functionals should be checked against the CPHF “gold standard” for the same molecules and basis sets.

Figure 3 compares the ϕ-dependence of the DFT PVES of H₂S₂ using the OLYP functional with that of the CPHF PVES for a 6-31G basis set: the OLYP functional overestimates the PVES as compared with CPHF by less than 15%. DFT-UPKS computation with the computationally efficient non-hybrid OLYP functional and a 6-31G or 6-311G basis set thus emerges as the most promising way forward for extension to larger molecules and polymers.

Fig. 3

Comparison of the DFT PVES (UPKS with the OLYP functional) and the CPHF PVES for H₂S₂ (in units of 10^− 18 hartree) as a function of torsion angle ϕ using a 6-31G basis set

Comparison with other Second Generation Computations

We consider first the computations of Lazzeretti and Zanasi (1997) using the “RPA” technique, which is equivalent to CPHF. As has already been pointed out by Berger and Quack (2000), these RPA results are consistently too large because they omit the two-electron term in the spin-orbit coupling operator without compensating for this by using an effective nuclear charge (such as that from Slater’s rules) instead of the bare nuclear charge. In an s,p-basis, the only matrix elements of $\hat V_{\!so}$ are those between p-orbitals, and Slater’s rules show that screening reduces the nuclear charge seen by an oxygen p-orbital from 8 to 4.55; this suggests that the RPA results of Lazzeretti and Zanasi (1997) should be scaled by a factor of approximately 4.55/8, or 1/1.76 for H₂O₂, and a corresponding factor of 1/1.35 for H₂S₂, to enable proper comparison between their RPA results and our CPHF results. Table 6 shows that when this scaling factor is included, their RPA results are in excellent agreement with our CPHF results.

Table 6 RPA PVES of H₂O₂ and H₂S₂ of Lazzeretti and Zanasi (1997) compared with our CPHF results

We turn now to the relativistic Dirac Hartree-Fock (DHF) computations of Laerdahl and Schwerdtfeger (1999). These authors tabulate the parity-violating matrix elements on each nucleus (note that no spin-orbit correction is needed in the relativistic formulation), so we have converted these matrix elements to total PVES values for comparison with our CPHF values in Table 7. There is evident sign and order of magnitude correspondence between their DHF values and our CPHF values, and furthermore their DHF PVES varies with dihedral angle ϕ in a similar manner to our CPHF values.

Table 7 DHF PVES of H₂O₂ and H₂S₂ of Laerdahl and Schwerdtfeger (1999) compared with our CPHF results

Turning now to the computations of Quack and co-workers, it is important to note, for the purposes of comparing results, that in their original paper (Bakasov et al. 1998) the PVES values computed using their CIS formalism must be multiplied by a factor of 2 due to computational artifacts (Bakasov et al. 1999). Table 8 compares our results (MacDHC: MacDermott, Hyde and Cohen) with the results of Berger and Quack (B&Q), obtained using a variety of methods (Berger and Quack 2000). In computing the spin-orbit coupling matrix elements, Berger and Quack (2000) used Z ^eff = 0.66Z for oxygen, while we use Z ^eff = 4.55 from Slater’s rules; for proper comparison, therefore, their results should be multipled by a factor (4.55/8)/0.66 = 0.86. When this is done, the agreement between their results (B&Q x 0.86) and our results (MacDHC) is good. Note that their RPA results include perturbative coupling but not electron correlation, and so should be compared to our CPHF results; likewise their CIS-RHF and our UPKS results include electron correlation but not perturbative coupling, while their various MC-LR results (CIS-LR, CISDT-LR, RASCIS-LR, RASCI(full)-LR, CASSCR-LR and RASSCF-LR) and our CPKS results include both perturbative coupling and electron correlation. The (computationally highly expensive) CI and MC-LR techniques thus provide PVES results that are not vastly different from the (computationally more efficient) CPHF and DFT techniques demonstrated in this paper.

Table 8 H₂O₂ PVES results of Berger and Quack (B&Q) at dihedral angles ϕ = 30°, 120°, compared with our results (MacDHC) using various computational methods with a 6-31G basis set

We now consider the DFT computations of Helgaker and co-workers (Hennum et al. 2002), which focus on H₂O₂ only and use mainly the B3LYP functional. They use DFT linear-response theory, which should encompass the perturbative coupling that would be necessary with B3LYP, which is a hybrid functional. Our B3LYP H₂O₂ results are in good agreement with theirs where we have used the same basis sets. They use an effective nuclear charge of 0.66Z for all atoms and orbitals, so in Table 9 we have multiplied their results (HHK) by 0.86 (because for oxygen p-orbitals we use an effective nuclear charge of 4.55, and (4.55/8)/0.66 = 0.86); when multiplied by this factor, the HHK results are in good agreement with ours (MacDHC). Although the factor 0.86 is appropriate for comparing their results with ours in an s,p-basis, their use of the same screening factor, 0.66Z, for all atoms and orbitals inflates their values as compared to ours when using basis sets with polarization functions (because we use a different, much lower Slater factor for d-orbitals). Interestingly, Hennum et al. (2002) obtain even larger results using B3LYP with some extra-large basis sets, especially those augmented with steep p functions. This is almost certainly a real effect (rather than the typical DFT overestimation of the PVED), because their results agree with those obtained by Laerdahl and Schwerdtfeger (1999) using a totally different method—Dirac Hartree-Fock—with similar steep p functions. This demonstrates the importance of obtaining a good representation of the wavefunction close to the nucleus when good numerical accuracy is required (since the weak force is a contact interaction). However, the 6-31G results of Hennum et al. (2002) actually recover 60% of the value obtained with the larger basis sets augmented with steep p functions, demonstrating that 6-31G remains adequate where only the sign and order of magnitude is required, as in origin of biomolecular chirality applications.

Table 9 H₂O₂ DFT PVES results of Hennum et al. (2002) (HHK) using the B3LYP functional, compared with our DFT results (MacDHC) using the B3LYP, HCTH and OLYP functionals, all at a fixed dihedral angle of ϕ = 45°

It is clear that very similar results are obtained using either pertubative coupling techniques (such as CPHF/RPA) or electron correlation techniques (such as CI or non-hybrid DFT), and that the PVES values do not change greatly with both perturbative coupling and electron correlation (as in MC-LR or CP-hybrid-DFT). There is therefore no extra advantage in using the computationally more demanding MC-LR and CP-hybrid-DFT techniques, and we believe the most efficient way forward for larger molecules lies in UPKS computations with non-hybrid functionals such as HCTH.

Density Functional vs Wavefunction Methods

Of the available techniques mentioned above—density functional theory (DFT) and the three types of wavefunction methods (CPHF/RPA, CI, and Dirac-HF)—DFT is definitely the cheapest computationally, and therefore the most readily extendable to larger systems. In DFT, the CPU time rises only as the third or fourth power of the number of electrons in the molecule, which is closer to the scaling of ordinary HF-SCF than to the more unfavourable scaling of the CI methods to which DFT is an alternative.

We therefore conclude that DFT is ideal for computing the PVED of large biomolecules, especially polymers. The tendency of DFT to overestimate the PVED (by up to 70%, depending on the molecule, the functional and the basis set, but more typically by 15–30%, and by less than 15% with the OLYP functional) is not important where only sign and order of magnitude is required, as in origin of life applications. But for comparison with possible future experimental measurement of the PVED (Quack 2002; MacDermott and Hegstrom 2004), as in fundamental physics applications, a Dirac Hartree-Fock approach as in Laerdahl and Schwerdtfeger (1999) (with steep s and p functions and perhaps other improvements to the wavefunction close to the nucleus) is essential: density-functional theory would not be sufficiently accurate for this purpose, and relativistic methods not only allow proper treatment of heavier nuclei, but are also free from inaccuracies arising from approximate treatments of the spin-orbit coupling matrix elements (since no spin-orbit correction to the wavefunction is needed to get a non-zero PVES in the Dirac treatment).

Conclusion

In this paper we presented coupled-perturbed and density functional computations of the PVED for H₂O₂ and H₂S₂, using a range of functionals and basis sets. Our computations are in good agreement with results from other groups using a variety of alternative approaches.

Of the currently available techniques, density functional theory (DFT) is the cheapest computationally, and therefore the most readily extendable to large biomolecules and polymers. Most importantly, in DFT computations there is no need to invoke perturbative coupling if a non-hybrid GGA functional is used, because coupled-perturbed Kohn-Sham and uncoupled-perturbed Kohn-Sham computations give identical results for non-hybrid functionals. Being able to do an uncoupled-perturbed DFT computation is a great advantage, because coupled-perturbed computations become rapidly more expensive for large molecules or larger basis sets. When using a hybrid functional, however (as in Hennum et al. (2002)), it is essential to invoke the more expensive coupled-perturbed DFT approach.

Of the available non-hybrid exchange-correlation functionals, the new HCTH (Hamprecht et al. 1998) and OLYP (Handy and Cohen 2001; Lee et al. 1988) functionals show the least basis set variation, and OLYP is also very fast computationally. Our results also show that little advantage accrues from using a basis set larger than 6-31G when only the sign and order of magnitude of the PVES is required.

Our overall conclusion is that for accurate determination of the PVED of small molecules, coupled-perturbed Hartree-Fock (or better still Dirac Hartree-Fock) is preferable to DFT where possible, because of DFT’s tendency to overestimate the PVED. But for cheap sign and order of magnitude estimation of the PVED of large biomolecules and polymers, DFT—specifically uncoupled-perturbed Kohn-Sham using a 6-31G basis set with a non-hybrid functional such as OLYP—is the most efficient and realistic way forward.