Electronic circular dichroism spectra using the algebraic diagrammatic construction schemes of the polarization propagator up to third order

Electronic circular dichroism (ECD) is a widely used linear optical effect to gain information about the stereochemistry of molecular systems.¹ Different enantiomers of a chiral molecule preferentially interact with one of the two circularly polarized components of plane polarized light. This results in a rotation of the plane of polarization, as observed in optical rotation measurements, and differential absorption, as detected by ECD spectroscopy. As a spectroscopic tool, ECD has found use in a wide range of scientific fields ranging from biochemistry, where it is commonly employed to detect the secondary structures of proteins,² to pharmaceutical synthesis, where the enantiomeric purity is of importance.³ Especially important in ECD spectroscopy is the sign of the measured signal as enantiomers produce mirror-like ECD spectra. This allows for quick and straightforward quality control of synthesized chiral molecules once accurate theoretical predictions have been carried out. Furthermore, achiral systems interacting with chiral systems, e.g., a solvent arranging itself around a solute in a structurally chiral manner, can give rise to a CD signal, which is a process referred to as induced circular dichroism (ICD). Thus, ICD provides information not only on the absolute configuration of the chiral system but also on the orientation of its achiral (e.g., solvent) counterpart.⁴ 

For the prediction of ECD spectra, many theoretical methods have been developed over the last decades with methods suited for different molecular sizes and environments. Progress has been made, and highly accurate ECD spectra can be simulated and compared directly with experiment.¹ For proteins and macro-molecules, semi-empirical methods have been shown to accurately reproduce ECD spectra of low-lying excitations.^5,6 Similarly, for molecules on the order of hundreds of atoms, time-dependent Hartree–Fock (TD-HF)^7–9 and time-dependent density functional theory (TD-DFT)^10–14 have proven successful. For ECD spectra of small chiral molecules, on the order of tens of atoms, configuration interaction (CI)¹⁵ and linear response based on a complete active space self-consistent field (CAS-SCF) reference¹⁶ can be applied. Furthermore, the well-known equation of motion coupled cluster (EOM-CC)^17,18 and the similar linear response coupled cluster (LR-CC) approach^19–22 have been applied to the calculation of ECD spectra and proven to be highly accurate.

The algebraic diagrammatic construction (ADC)²³ scheme for the polarization propagator is an ab initio method, which shares many features with CC excited-state methods. ADC requires the solution of a Hermitian eigenvalue equation and provides a hierarchy of approximations. In comparison with CC methods, this has the advantage of not requiring the solution of both the left- and right-hand side equations. It has been employed to predict several linear and non-linear properties, e.g., two-photon absorption (TPA)²⁴ and resonant inelastic x-ray scattering (RIXS).²⁵ Furthermore, static properties such as static polarizabilities and C₆ dispersion coefficients^26,27 have recently been investigated. The computational effort required for ADC(n) schemes scales as n⁴ for ADC(1), n⁵ for ADC(2), and n⁶ for ADC(2)-x and ADC(3/2) with n being the number of basis functions. For comparison, the corresponding CC schemes, CC2, coupled-cluster singles and doubles (CCSD), and CC3, formally scale as n⁵, n⁶, and n⁷, respectively, and, moreover, require an iterative solution for the corresponding ground state as well.

Here, we report the first implementation of ADC-based methods to simulate chiro-optical properties such as ECD.

We choose methyloxirane and methylthiirane molecules as well as their dimethylated derivatives to investigate the performance of the ADC methods for the calculation of ECD spectra. For these molecules, highly resolved gas phase spectra exist for comparison. Furthermore, we include the inherently chiral hydrogen peroxide H₂O₂ and dihydrogen disulfide H₂S₂ for their small size and well studied chiro-optical properties. Beside these relatively small molecules, we evaluate the ECD spectra of camphor, norcamphor, and fenchone using a third order ADC scheme in comparison with experimental gas phase spectra.^28–32 The latter molecules have previously been investigated using TD-DFT but not yet using a third order ab initio method. We also include epinephrine, a common medical drug and hormone,³³ for which there are aqueous ECD solution spectra reported.³⁴ The molecular structures of the investigated systems are shown in Fig. 1.

A complication in the calculation of ECD spectra is the lack of gauge invariance of approximate solutions of the Schrödinger equation, which can lead to origin dependent rotatory strengths in the length gauge. This is circumvented in one of the two ways: (i) using magnetic-field dependent “gauge-including” atomic orbitals (GIAOs, also referred to as London orbitals)^11,35–42 or alternatively (ii) by a gauge transformation to the velocity gauge, which is formally origin invariant.^19,43–47 The latter strategy is referred to as the velocity gauge formalism or gradient form of the rotatory strength, and this method will be used throughout this work. For the sake of completeness, the length gauge formulation of the rotatory strength is also briefly presented.

Experimentally, an ECD spectrum is recorded by measuring the difference in absorption between the left and right circularly polarized light,

where ϵ_L(ω) and ϵ_R(ω) are the absorption coefficients, at incident frequency ω, for the left and right circularly polarized light, respectively. The measured signal may then be related to the underlying rotatory strength, governing the excitation by integration over a frequency range [ω₁, ω₂],⁴⁸ as

where the photon energy ω is in eV and Δϵ(ω) is given in L mol⁻¹ cm⁻¹. The rotatory strength, R, is in the usually reported units of 10⁻⁴⁰ esu cm erg/G or 10⁻⁴⁰ c.g.s. units.⁴⁹ The crucial assumption of Eq. (2) is that the integrated peak results from one electronic transition only. This is not necessarily the case and not known a priori, which emphasizes the need for robust theoretical models assisting in the interpretation of the experimental ECD spectra.

The fundamental property governing the preferential absorption of circularly polarized light is the rotatory strength, R_n0. It is given by the well-known Rosenfeld equation⁵⁰ as the scalar product of the electric and magnetic transition dipole moments between the ground state |Ψ₀$$ and an excited state |Ψ_n$$⁠,

where $μ^$ and $m^$ are the electric and magnetic dipole operators, respectively. From the Rosenfeld equation, we can see that the electric transition dipole moment and the magnetic transition dipole moment must be non-orthogonal for a transition to be active in ECD. Using the electric dipole operator, $μ^=−er^$⁠, and the magnetic dipole operator, $m^=ieℏ2mecr^×∇^$⁠,⁴³ within Eq. (3) yields the rotatory strength in its length gauge formulation,

where e is the elementary charge, m_e is the electron mass, c is the speed of light, and ℏ is the reduced Planck constant. In the sum-over-states expression, this form of the rotatory strength sums up to zero, which results in zero rotatory strengths far away from resonance,⁵¹ but lacks gauge-origin independence.⁵² Hence, a translation of the origin of the position operator alters the calculated rotatory strength. As mentioned previously, we circumvent the origin dependence of the rotatory strengths in the length gauge by employing the velocity gauge. By using the fact that the position operator commutes with every term in the electronic molecular Hamiltonian, except the kinetic energy operator,^47,53

and applying this to Eq. (4) yields the rotatory strength in the velocity (gradient) form

The corresponding oscillator strength of the same transition is computed in the length gauge as⁴³ 

which in analogy to the rotatory strength can be transformed to yield the oscillator strength in the velocity gauge,

The two forms of the oscillator strength are origin independent, and they are usually computed together as a measure of basis sets completeness.

In the context of the algebraic diagrammatic construction scheme for the polarization propagator, the excitation energies and transition moments required in order to evaluate Eq. (6) are obtained by finding the eigenvectors and eigenvalues of the Hermitian ADC matrix M, which is a representation of the electronic Hamiltonian $Ĥ$ in the so-called intermediate states $|Ψ̃I$ shifted by the exact ground-state energy, E₀.^54,55 The ADC matrix is connected to the diagonal representation of the shifted Hamiltonian in exact states,

by unitary transformation

with the eigenvector matrix X. Ground-to-excited state transition moments for a given operator $Ô$ can be evaluated by contracting an eigenvector xⁿ with the vector of so-called modified transition moments $F(Ô)$⁠,

which correspond to the transition moments between the ground state and the intermediate states,

To efficiently evaluate transition moments for different operators, the modified transition moments are not constructed explicitly, but rather the transition density ρ⁰ⁿ is calculated,

where $cp†$ and c_q are creation and annihilation operators in second quantization. The transition moments are then evaluated as

with the property integrals $opq=ϕp|Ô|ϕq$ in the molecular orbital basis {ϕ_p}. By consistently constructing the intermediate states up to a given order n of perturbation theory using the Møller–Plesset Hamiltonian partitioning, one arrives at the corresponding ADC(n) scheme. Using the eigenvalues and eigenvectors of the third-order ADC matrix with the second-order modified transition moments is referred to as ADC(3/2). This mixed third and second order scheme is used throughout this paper and called ADC(3) for brevity.^56,57

Full-ADC is formally gauge invariant. However, this is not true for truncated ADC schemes. The remaining difference between the length and velocity forms of the oscillator and rotatory strengths occurs in the highest-order contributions of their perturbational expansion or, in other words, in $O(n)$ for ADC(n). Consequently, the gauge-variance of ADC(3/2) is of the same order as ADC(2), i.e., the differences occur in the $O(2)$ contributions. For further theoretical details regarding the length and velocity gauge behavior of ADC(n) methods, see Ref. 58.

To simulate the ECD spectra, the computed rotatory strengths are converted from a.u. to 10⁻⁴⁰ c.g.s. using the conversion factor 471.44 and subsequently broadened by a Lorentzian function,

where Δϵ_n is the peak intensity of the nth transition given in L mol⁻¹ cm⁻¹, R_n0 is the rotatory strength in 10⁻⁴⁰ c.g.s., γ_n is the Lorentzian broadening factor, and ω and ω_n0 are the incident optical frequency and excitation energy in eV. A broadening factor corresponding to full width at half maximum of γ_n = 1000 cm⁻¹ is used throughout this work.⁵⁹ 

Calculations of oscillator and rotatory strengths in both length and velocity gauges using ADC were performed using a development version of Q-Chem 5.2.⁶⁰ For the CC calculations of the same properties, Dalton2018 was used.⁶¹ Dunning basis sets were employed throughout the paper.^62–63 The origin for the position operator was set to the hetero-atom for the oxiranes, thiiranes, and camphor and its derivatives and to one of the oxygen/sulfur atoms in the symmetrical H₂O₂ and H₂S₂ molecules. Geometries of methyloxirane, dimethyloxirane, methylthiirane, and dimethylthiirane were optimized using MP2/aug-cc-pVTZ and those of H₂O₂ and H₂S₂ at the MP2/d-aug-cc-pVTZ level. For camphor, norcamphor, and fenchone, DFT/CAM-B3LYP⁶⁴ and aug-cc-pVTZ were used for the geometry optimization.

Methyloxirane has been widely used for the evaluation of quantum chemical methods for the simulation of ECD spectroscopy, primarily due to its small size as well as its available gas phase circular dichroism spectrum.^65,66 Furthermore, there is extensive literature available, both experimental and theoretical work, covering electronic structure,⁶⁷ excited state properties,⁶⁸ and two-photon absorption ECD.⁶⁹ CC methods have been proven successful in reproducing the gas phase ECD spectrum.⁷⁰ The inclusion of molecular dynamics has addressed the water solution spectrum,⁷¹ and the use of a polarizable embedding model has also addressed the spectrum in various solvents.^72–74 

It has been shown that the inclusion of diffuse functions is important to ensure correct signs of the optical rotatory strength of methyloxirane.⁷⁵ Pople basis sets seem inadequate even with diffuse functions, while Dunning basis sets have been proven more capable.¹¹ Therefore, we performed an extensive basis set study for ADC(2), ADC(3), CC2, and CCSD using Dunning basis sets and comparing to experimentally measured values of methyloxirane. The excitation energies and oscillator and rotatory strengths calculated at the ADC(2), ADC(3), CC2, and CCSD levels of theory for the nine lowest excited states of (R)-methyloxirane and (R)-methylthiirane are collected in Table I (see Tables S1–S4 in the supplementary material for all excited states computed). The simulated ECD and UV/VIS spectra at the same levels of theory of the 15 energetically lowest excited states are shown in Fig. 2.

Starting with ADC(2) and CC2, the smallest basis of the cc-pVXZ series (X: D, T, Q), cc-pVDZ, yields nearly identical excitation energies for ADC(2) and CC2 with the first vertically excited state calculated at 8.24 eV at the ADC(2) level and 8.35 eV at the CC2 level, which should be compared with the experimental value of 7.08 eV.⁶⁵ Increasing the basis set size from double- to quadruple-zeta lowers the excitation energy to 7.51 eV and 7.54 eV for ADC(2) and CC2, respectively. However, the rotatory strength of the first transition is largely unaffected and is roughly overestimated by a factor of 2 compared to experiment. However, the reported value in Ref. 65 is based on Eq. (2) with the limits of integration determined by a low-level CI calculation possibly impairing the analysis. In general, careful analysis is always needed when comparing computed and experimental rotatory strengths.

Turning to the singly augmented, aug-cc-pVXZ series (X: D, T), aug-cc-pVDZ yields similar excitation energies for the first vertically excited state at the ADC(2) and CC2 level with 6.23 eV and 6.37 eV, respectively. Increasing the basis from double- to triple-zeta improves the excitation energies to 6.52 eV and 6.55 eV for ADC(2) and CC2, respectively. In general, the rotatory strengths are much more sensitive to augmentation of the basis set than the excitation energies. Even with aug-cc-pVDZ, the rotatory strength of the first transition is now within 20% of experiment and within 15% for aug-cc-pVTZ, which has been overestimated by roughly 100% using non-augmented basis sets.

Moving to the doubly augmented, d-aug-cc-pVXZ series (X: D, T), the rotatory strength of the first excited state is within 10% of experiment for both ADC(2) and CC2 using d-aug-cc-pVDZ, but the excitation energy (<0.04 eV) is hardly affected. Going to d-aug-cc-pVTZ, the rotatory strength is largely unaffected, while the excitation energy is increased, as was the case for the singly augmented aug-cc-pVTZ basis. However, double augmentation has a substantial impact on the rotatory strengths of higher-lying excited states. For cc-pVXZ and aug-cc-pVXZ basis sets, the second vertically excited state has an incorrect positive sign of the rotatory strength at ADC(2) and CC2 levels of theory but is correctly described as negative using d-aug-cc-pVXZ basis sets. Higher order methods such as ADC(3), CCSD, or CC3 reproduce the correct sign regardless of which basis set is employed. This illustrates the importance of diffuse functions as well as the degree of electron correlation included for a correct assignment of higher excited states.

Overall, the differences between the length and velocity forms of the rotatory strength are relatively small for the first electronic transition, i.e., on the order of ∼10–20% for cc-pVXZ, ∼8–10% for aug-cc-pVXZ, and ∼7–10% for d-aug-cc-pVXZ at the ADC(2) and CC2 level of theory. The difference is on the same order of magnitude for energetically higher-lying excited states with minor exceptions (to see the effect on a shift in the gauge-origin on the length gauge rotatory strength, refer Fig. S5 of the supplementary material). The difference between the length and velocity forms of the oscillator strength is very similar to the one observed for the rotatory strengths.

A comparison of the performance of the individual methods in the calculation of the ECD spectrum of methyloxirane when the d-aug-cc-pVDZ basis set is used reveals ADC(2) and CC2 to yield very similar spectra, as obtained from the first 15 excitations calculated (Table S3 in the supplementary material). All spectra are shifted for the first state to correspond to the experimental value of 7.10 eV. The shifts required are +0.92 eV for ADC(2) and +0.78 eV for CC2, respectively. Two prominent features of the experimental spectrum, A and B, are well reproduced by ADC(2) and CC2 (Fig. 2). Furthermore, at these levels of theory, A–C arise from oxygen lone-pair excitations to Rydberg orbitals. According to our calculations, A is a convolution of two states of negative amplitude and B is a convolution of three states with an overall positive amplitude. Finally, C corresponds to a convolution of three states and is incorrectly predicted as positive by both methods.

The spectra calculated at the ADC(3) and CCSD level of theory are again shifted by −0.74 eV and −0.02 eV, respectively, for the first vertical excitation energy to match the experimental value. At these levels of theory, in addition to A and B, C is also correctly predicted as negative albeit red-shifted by ∼0.5 eV for ADC(3).

Methylthiirane is the sulfur containing analog of methyloxirane, and again a high-quality experimental gas phase ECD spectrum,^65,66 as well as experimental vibrational circular dichroism (VCD)⁷⁶ and Raman optical activity (ROA)⁷⁷ spectra, is available. The ECD spectrum and the optical rotatory dispersion have also been addressed in theoretical studies.^78,79 The experimental ECD spectrum of methylthiirane possesses three distinct features, A, B, and C in analogy to methyloxirane. Using the d-aug-cc-pVDZ basis set, ADC(2) and CC2 yield very similar transition energies and rotatory strengths for the energetically lowest 15 computed states (Table S4 in the supplementary material). All excited states correspond to electronic transitions from sulfur lone-pair orbitals into Rydberg orbitals. ADC(2) and CC2 yield ECD spectra correctly predicting the signs of A, B, and C. The excitation energies are red-shifted by roughly 0.4 eV in comparison to experiment (Fig. 2). The first peak arises from a weakly absorbing state with a negative rotatory strength, which has also been experimentally observed but was not shown in the experimental spectrum as it was barely visible. The first significant feature of the experimental spectrum, A, is a convolution of two peaks, which are correctly predicted as negative, but the amplitude is underestimated by a factor of 2 at these levels of theory. B and C correspond to convolutions of multiple states, the signs of which are correctly predicted using ADC(2) and CC2.

Proceeding to ADC(3), the excitation energies require no shift and transition amplitudes agree very well with experiment. The two electronic states constituting the spectral feature A, a₁ and a₂, lie within ∼0.02 eV of the experimental spectrum at the ADC(3) level, but the rotatory strengths of these states appear to be swapped compared to experiment. However, it cannot be excluded that these arise from vibrational effects. The improvement is of equal quality for CCSD, although slightly red-shifted compared to experiment. B consists of four vibronic peaks in the experimental spectrum. At the ADC(3) level, two excited electronic states are computed, exhibiting the correct excitation energies and signs of their rotatory strengths. C is a convolution of several peaks in the experimental spectrum and is incorrectly described as negative by both ADC(3) and CCSD.

In addition, for dimethyloxirane and dimethylthiirane, high-quality experimental gas phase ECD spectra are available^65,66 and show features similar to the ones of methyloxirane and methylthiirane. The excitation energies and oscillator and rotatory strengths calculated at the ADC(2), ADC(3), CC2, and CCSD levels of theory for the nine lowest excited states of (R,R)-dimethyloxirane and (R,R)-dimethylthiirane are collected in Table II (see Tables S5–S8 in the supplementary material for all excited states computed). The simulated ECD and UV/VIS spectra at the same levels of theory of the 15 energetically lowest excited states are shown in Fig. 3.

At the ADC(2) and CC2 levels of theory using the d-aug-cc-pVDZ basis set, both methods produce very similar spectra for the energetically lowest 15 states computed with the ADC(2) excitation energies red-shifted by ∼0.15 eV–0.20 eV compared to CC2. All spectra are shifted for the first state to correspond to the experimental value of 7.0 eV. The shifts required are +0.92 eV and +0.78 eV for ADC(2) and CC2, respectively. At these levels of theory, the vertical excited states arise exclusively from lone-pair orbital excitations from the oxygen atom into Rydberg orbitals. Furthermore, the spectral feature A (Fig. 3) arises solely from the lowest transition, the excitation energy of which is red-shifted by ∼0.8 eV from experiment, reproduced as 6.08 eV and 6.22 eV for ADC(2) and CC2, respectively. The rotatory strength, however, is in reasonable agreement with experiment, although it is roughly overestimated by 40% for both ADC(2) and CC2. Furthermore, ADC(2) and CC2 show B arising solely from the second excited state for which the excitation energy is calculated at 6.38 eV and 6.53 eV for ADC(2) and CC2, respectively, which is again red-shifted by ∼0.9 eV compared to experiment. The rotatory strength is overestimated by an order of magnitude. The band C at these levels of theory arises from the third and fourth excited states, which in comparison with experiment are red-shifted by roughly 1 eV for both ADC(2) and CC2. D arises from the fifth and sixth excited states, which are again too low in energy by roughly 1.1 eV.

The spectra of the higher order methods ADC(3) and CCSD are also shifted to correspond to experiment. The shifts required are −0.81 eV and −0.07 eV for ADC(3) and CCSD, respectively. ADC(3) again yields higher excitation energies, as was already the case for methyloxirane. The CCSD excitation energies on the other hand agree to within ∼0.10 eV with experiment. At the ADC(3) level, the spectral feature A arises from two states with negative rotatory strengths, separated by 0.10 eV. These peaks appear in agreement with the experimental ECD spectrum. We note that other studies have assigned these peaks to vibronic structures originating from up to two electronic states.⁸¹ Furthermore, these two peaks are not reproduced in the CCSD spectrum, which, alike ADC(2) and CC2, shows A arising from a single excited electronic state, whose width closely agrees with the experimental one. At the ADC(3) level of theory, B is a convolution of the third and fourth excited state, the first sharply positive and latter negative. This results in a relatively weak positive feature, in reasonably good agreement with experiment. A similar convolution is observed at the CCSD level, however, of the second and third excited state. The larger energetic separation between the two states in CCSD results in a more intense positive B-band feature, whereas the negative component merges with the negative C band. C is again a convolution of many states at the ADC(3) level resulting in a broad positive band, not present in experiment, followed by a blue-shifted negative peak. CCSD predicts the overall negative band of C, which is originating from the fourth and fifth excited states (with a shoulder due to the third one), in better agreement with experiment. D is correctly predicted as positive and arising from the tenth excited state using ADC(3) and a convolution of many states using CCSD.

The three bands, A–C, of dimethylthiirane are again investigated using the d-aug-cc-pVDZ basis set and are characterized by lone-pair excitations into Rydberg orbitals by ADC(2) and CC2. The excitation energy of the first excited state is calculated as 5.03 eV and 5.08 eV by ADC(2) and CC2, respectively, which are blue-shifted by roughly 0.5 eV, compared to the reported experimental value of 4.5 eV. The first rotatory strength exhibits an error of more than an order of magnitude, 5.7 (×10⁻⁴⁰ c.g.s.) and 4.3 (×10⁻⁴⁰ c.g.s.) for ADC(2) and CC2, respectively, compared to the experimental value of 0.1 (×10⁻⁴⁰ c.g.s.).⁸⁰ However, this is a weakly active state and as such experimental limitations may also play a role.

The first distinct feature A in the experimental spectrum is a convolution of the second to fifth excited states at these levels of theory. The second feature, B, is a convolution of the sixth to eighth excited state, and the third band, C, arises from higher-lying excited states above the eighth state. It is interesting to recognize that the calculated spectra at ADC(2) and CC2 levels could be mistaken as that of red-shifted (S,S)-dimethylthiirane; however, this is not the case for the higher order methods.

Using ADC(3) and CCSD with the d-aug-cc-pVDZ basis set, the lowest excitation energy remains mostly unaltered at 4.97 eV and 5.04 eV, respectively. Unlike ADC(2), ADC(3) yields a rotatory strength for the first weak peak of similar magnitude to that of the experiment, deviating by ∼10%, however, still overestimated by an order of magnitude for CCSD. Here, 14 excited states are considered, and for the first and second band feature A and B, ADC(3) results agree well with experiment, although the rotatory strength of A is still overestimated by a factor of 2. Significantly, the last band feature, C, is correctly predicted as negative at the ADC(3) level resulting from a convolution of several excited states but incorrectly predicted as positive using CCSD. We note, however, that more excited states may be needed at the CCSD level to reproduce this spectral region.

The chiro-optical properties of H₂O₂ and H₂S₂ have been the subject of benchmark calculations due to their small size and complex electronic structure.^14,82,83 The dependence of the optical rotatory strength on the dihedral angle of H₂S₂ has been studied at the CIS level. It was discovered that H₂S₂ has degenerate excited states at a dihedral angle of ∼90°, which are poorly described using a HF wavefunction. Later, the performance of various density functionals was benchmarked against MRCI and CC2 calculations.⁸⁴ As can be expected from a sulfur containing chromophore, low-lying excitations are characterized by Rydberg to a certain extent.⁸⁴ Fewer studies have been carried out on H₂O₂, and for neither molecule, the ECD spectra have been simulated using a third order ab initio method. The excitation energies and oscillator and rotatory strengths calculated at the ADC(2), ADC(3), CC2, CCSD, and CC3 levels of theory for the nine lowest excited states of H₂O₂ and H₂S₂ are collected in Table III (see Tables S9–S12 in the supplementary material for all excited states computed). The simulated ECD and UV/VIS spectra at the same levels of theory of the 15 energetically lowest excited states are shown in Fig. 4.