Origin invariant optical rotation in the length dipole gauge without London atomic orbitals

Optical rotation (OR) is one of the manifestations of the optical activity of chiral molecules, where the plane of polarization of linearly polarized light passing through a sample with an enantiomeric excess of a chiral compound is rotated by a certain angle.¹ The magnitude of the rotation is the same for the two enantiomers but the direction is opposite. The OR is usually reported as a normalized quantity called specific rotation, [α]_ω or [α]_λ, where ω and λ are the frequency and wavelength of the impinging light. Measurement of the OR can help in the assignment of the absolute configuration of a sample when they are aided by accurate simulations of this molecular property. This is fundamental for practical applications in biochemistry or in the pharmaceutical industry, because of the homochiral nature of all biologically relevant molecules and supramolecular systems.

Over the past few decades, accurate quantum mechanical methods based on density functional theory (DFT) and coupled cluster (CC) theory have been developed for the calculation of chiroptical properties.^2–17 These methods use response theory to evaluate the perturbed electron density from the external perturbation and compute the appropriate OR tensor. Given that all of these methods represent an approximate solution to the Schrödinger equation with a finite basis expansion for the electron density, the numerical results depend on the choice of gauge for the electric dipole operator.^12,18 Two typical choices are the length and velocity gauges, which have advantages and drawbacks. Length gauge (LG) calculations may converge faster toward the complete basis set limit,^3,19 but they provide a result that is origin dependent due to the properties of the magnetic dipole operator. Velocity gauge (VG) calculations are intrinsically origin independent, but they include an unphysical static contribution.¹⁴ The issue of the origin dependence of the LG approach is alleviated (but not eliminated) in DFT by using London orbitals, also known as gauge including atomic orbitals (GIAOs).^16,20,21 The issue of the static contribution in the VG approach is addressed by computing such a contribution explicitly and subtracting it out, an approach called modified velocity gauge (MVG).¹⁴ The LG approach is typically favored in DFT calculations because the extra computational cost introduced by the GIAO terms is significantly smaller than that required for the evaluation of the static contribution in the MVG approach. However, the origin dependence of the LG approach cannot be easily eliminated in CC calculations because the response function is typically computed without the orbital response in order to avoid unphysical poles in the response function.²² This frozen-orbital approach makes the GIAOs ineffectual. Alternative formulations of CC methods with orbital optimization have been proposed,^23,24 but they are not commonly used because they are computationally intensive and often have convergence issues. Therefore, the vast majority of CC calculations of OR use the MVG approach.

In this communication, we introduce a simple approach to eliminate the origin dependence in LG calculations without recourse to GIAOs or complicated formulations of CC methods. Although this approach can be applied to any approximate method, here, we focus on CC because of the current lack of production-level alternatives. This origin invariant LG method [LG(OI)] is based on the singular value decomposition of the electric dipole–electric dipole polarizability tensor using a mixed formulation of the dipole operator. This step provides the unitary transformation necessary to make the diagonal elements and the trace of the OR tensor in the length gauge origin invariant. The computational cost of this procedure is essentially equivalent to that of MVG for a single frequency calculation.

To the best of our knowledge, this communication reports the first origin invariant LG calculations of OR with conventional CC methods. Given the sensitivity of CC methods to the choice of basis set and the possible better dependence on the basis set size of the LG approach compared to MVG, this work provides the means to compute this important molecular property more accurately and efficiently than possible before. This paper is organized as follows: the derivation of the equations is presented in Sec. II, preliminary calculations are discussed in Sec. III, and concluding remarks are summarized in Sec. IV.

The OR is related to the mixed electric dipole–magnetic dipole polarizability tensor,

where μ is the electric dipole operator, m is the magnetic dipole operator, with Greek indices denoting Cartesian coordinates, ω is the frequency of the incident electromagnetic radiation, and |ψ_j⟩ and ω_j0 are the jth excited state wave function and excitation frequency, respectively. We use atomic units throughout the paper, except when otherwise specified. These definitions are valid for non-resonant optical activity (ω_j0 ≉ ω) calculations; resonant optical activity is discussed in greater detail elsewhere.^1,25,26 For molecules in a homogeneous phase, gas or solution, only a spatially averaged OR can be measured, which is proportional to the trace of the β tensor in Eq. (1). The OR is commonly reported in terms of specific rotation {deg [dm (g/ml)]⁻¹}, given by

where β and ω are given in atomic units, ℏ is the reduced Planck’s constant (J s), N_A is Avogadro’s number, c is the speed of light (m/s), m_e is the electron rest mass (kg), and M is the molecular mass (amu). Technically, the full OR tensor (the Buckingham/Dunn B tensor) includes contributions from the electric dipole–electric quadrupole polarizability tensor.^25,27,28 However, the latter does not contribute to [α]_ω in Eq. (2) because the quadrupole operator is traceless, and thus, this tensor will be neglected here. Nevertheless, such contribution cannot be ignored for spatially oriented systems.²⁸ 

Following the notation of Lazzeretti and co-workers,^18,29,30 the origin dependence of the β tensor can be expressed as

where O is a particular choice of origin of the coordinate system and

is a displaced origin, ϵ is the Levi-Civita tensor, and α is the electric dipole–electric dipole polarizability tensor. For exact wave functions, the second term on the right-hand side of Eq. (3) does not contribute to the trace of β because of the symmetry properties of the ϵ and α tensors, and therefore, Tr[β(O)] = Tr[β(O′)]. However, in approximate calculations, this may not be the case, depending on the choice of representation for the electric dipole operator. In particular, using the length gauge, Eq. (3) becomes¹⁸ 

where the superscripts indicate the representation for the electric and magnetic dipole operators,

where r and p are the position and momentum operators, and the sums run over the number of electrons. The equivalence between the R and P definitions of the electric dipole operator for exact wave functions comes from the commutation relation,

where H is the Hamiltonian operator and i is the imaginary unit. Using the hypervirial relationships, Eq. (7) leads to²⁹ 

Therefore, one can choose any equivalent definition of the electric dipole operator in Eq. (3) and obtain the same numerical result. However, the electric dipole–electric dipole polarizability tensor using the mixed representation of the electric dipole operator

is not a symmetric tensor for approximate wave functions because the relations in Eqs. (7) and (8) do not hold. Thus, β^(R,L) is origin dependent for approximate methods. As mentioned in the Introduction, this problem can be partially overcome in approximate DFT methods by employing GIAOs, but such a strategy is not feasible with conventional CC methods because the response of the molecular orbitals is neglected to avoid unphysical poles in the response function.²² Therefore, CC-LG calculations remain origin dependent.

Such limitations can be overcome by extending a procedure suggested by Pelloni and Lazzeretti in Ref. 18. The OR computed in the velocity gauge is origin invariant because it is expressed as

Since α^(P,P) is symmetric even for approximate methods, the contribution from the second term on the right-hand side of Eq. (10) to the trace of β^(P,L) disappears. Although Tr[β^(P,L)] is origin invariant, the diagonal elements are still origin dependent. Pelloni and Lazzeretti propose to render the individual diagonal elements origin invariant by re-orienting the molecule along the principal axis of the α^(P,P) tensor.¹⁸ In practice, it is sufficient to evaluate the unitary transformation that diagonalized α^(P,P) and directly apply such transformation to β^(P,L) in Eq. (10).

A similar procedure can be used to make the diagonal elements of the β^(R,L) tensor in Eq. (5) origin invariant. One can perform a singular value decomposition of α^(R,P),

such that $αD(R,P)$ is diagonal, and U and V^† are unitary transformations. By applying the inverse transformation to the β^(R,L) tensor in Eq. (5),

We obtain an OR tensor with origin invariant diagonal elements and trace

This approach requires the evaluation of α^(R,P) in Eq. (9), but is valid with any approximate method and does not require the use of GIAOs. We will refer to the approach to evaluate the β^(R,L) tensor in Eq. (12) as “LG(OI),” for origin invariant.

Note that Pedersen et al. have defined an origin-dependent vector,¹⁴ 

such that Eq. (5) becomes

The Δ vector has been used as a “diagnostic” of the quality of LG results.¹³ However, Δ can only say “how quickly” the origin dependence of the OR grows by displacing the origin from one point to another, see Eq. (15), but it does not provide any information about the error for a specific choice of origin. The LG(OI) approach solves this issue because the transformation in Eq. (12) renders Δ null and all choices of origin are equivalent, see Eq. (13).

In production calculations, the sum-over-state formulas for the various polarizability tensors in Eqs. (1) and (9) are impractical, and linear response theory is used instead. For CC, the linear response function is^12,22

where Φ⁰ is the reference wave function, typically the Hartree–Fock (HF) determinant, and T/Λ are the CC excitation/de-excitation operators.^31,32 The overbar indicates the similarity transformation,

and the subscript N indicates the normal product form of an operator,^31,32

$T−ωX$ is the CC T operator perturbed by X at frequency −ω, and analogously for perturbation Y at frequency ω. The operator C^±ω takes the complex conjugate of the entire expression and adds it to the original one, and the operator P(XY) swaps the perturbed operators,

The expression in Eq. (16) corresponds to the various polarizability tensors depending on the choice of X and Y perturbations,

The evaluation of the linear response function in Eq. (16) requires the solution of four independent linear systems of equations to evaluate the perturbed T amplitudes for the X/Y perturbations at the +/−ω frequencies. Since the dipole operators have three Cartesian components, we need to compute 12 sets of amplitudes to evaluate the β^(R,L) tensor and 6 more for the α^(R,P) tensor, for a total of 18 linear systems of equations per ω. This is the same number of amplitude sets that are necessary for the MVG approach since we need to evaluate the amplitudes corresponding to the static term (ω = 0). Therefore, the MVG and LG(OI) approaches are computationally equivalent for a single frequency calculation. For a multi-frequency calculation, the MVG is less demanding because the static term needs to be evaluated only once. Nevertheless, the LG approach may converge more quickly to the basis set limit than MVG,¹⁹ thus providing more accurate results for a finite basis set.

The LG(OI) approach is less computationally demanding than the MVG and LG-GIAO approaches for variational methods like HF and DFT because one would only need to evaluate the perturbed density for the length gauge electric dipole, without extra GIAO terms, then contract it with the magnetic dipole integrals to obtain β^(R,L) and with the velocity gauge electric dipole integrals to obtain α^(R,P). The LG(OI) approach as presented here cannot be used to compute the full origin invariant B tensor because the off-diagonal elements of $β̃(R,L)$ are still origin dependent; however, evaluating the full B tensor is uncommon.

We present preliminary calculations of specific rotation using the CC-LG(OI) method on three prototypical chiral molecules: (−)-hydrogen peroxide (1), (S)-(+)-2,3-pentadiene (2), and (S)-(−)-norbornenone (3). The origin of the coordinate system is set in the center of mass of the molecule (O ≡ O_CM) because this is a well-defined choice. A displaced origin is also considered as in Eq. (4), with a displacement vector d = d(1,1,1)^† and d = −1000 Å. The geometries of 2 and 3 were taken from Ref. 33 and reoriented so that the origin is set in O_CM. The geometry of 1 was optimized at the B3LYP/aug-cc-pVTZ^34–36 level for consistency with the other two molecules. All geometries are available in the supplementary material. The OR calculations were performed with the CC with single and double excitations (CCSD) method and the aug-cc-pVDZ^35,36 basis set at three wavelengths: 355 nm, 589.3 nm, and 633 nm. The first and last wavelengths are used in gas phase measurements of specific rotation,^33,37 while the middle one is the sodium D-line that is typically used in solution phase measurements. The calculations used the following thresholds for integrals, energy, and wave function: 10⁻¹² a.u. for the two-electron integral accuracy, 10⁻⁸ convergence in the root-mean-square (rms) and 10⁻⁶ convergence in the maximum element (Max) of the HF density matrix, respectively; the thresholds for the CCSD part are: 10⁻⁶ convergence of the rms and Max of the T and Λ amplitudes and 10⁻⁸ a.u. convergence in the CC correlation energy, 10⁻¹⁰ and 10⁻⁹ convergence in the rms and Max of the residuals in the linear response equations, respectively.³⁸ All calculations were performed with a development version of the GAUSSIAN suite of programs, and the thresholds above are the default values in GAUSSIAN.³⁹ 

The specific rotations calculated with the conventional LG and LG(OI) approaches with the origin set at O and O′ for the test molecules are reported in Table I. The corresponding individual diagonal elements are reported in Tables S4–S6 of the supplementary material. The results with the LG approach and the O′ origin are clearly nonsensical, while the LG(OI) results are origin invariant, as expected. The differences between the O and O′ results for the LG(OI) approach, most evident at 355 nm, are only due to the numerical precision of the calculations. We verified this by repeating the calculation for molecule 1 with tighter thresholds: two orders of magnitude tighter for the integrals, four orders of magnitude for the HF equations, and three orders of magnitude for the CCSD equations. The [α]_λ values for 1 with the tighter thresholds are reported in Table S7 of the supplementary material: the values computed with the origin at O are essentially identical to those in Table I, and the results with the origin at O′ now perfectly match those with the origin at O. It is likely that this is the case also for the other two molecules and that the default thresholds are appropriate when the origin is placed at or near the center of mass of the molecule.

A comparison between the LG and LG(OI) results with the origin set at O indicates that this is a reasonable choice for molecules 1 and 3, but for 2, there is a significant residual origin dependence of the order of 18%–25%. This indicates that the center of mass may not always be an appropriate choice of origin of the coordinate system for the conventional CC-LG approach.

In this communication, we present a simple and straightforward way to perform origin invariant specific rotation calculations in the length gauge with approximate methods. The approach, LG(OI), can be used with any approximate method and does not require the evaluation of the molecular orbital response or GIAOs. Instead, one needs to evaluate the electric dipole–electric dipole polarizability with a mixed dipole representation, α^(R,P) in Eq. (9). A singular value decomposition of the α^(R,P) tensor provides the unitary transformations to remove the origin dependence of the diagonal elements of the OR tensor, β^(R,L) in Eq. (10). The transformed $β̃(R,L)$ tensor in Eq. (12) provides origin invariant [α]_ω.

We implemented the LG(OI) approach with the CCSD method and tested it on three chiral molecules. We compare the results obtained with the LG and LG(OI) approaches by putting the origin of the coordinate system at the center of mass of the molecule and at a position displaced by 1000 Å in every Cartesian direction. The results, reported in Table I, show that indeed the LG(OI) results are origin invariant (within the numerical accuracy of the calculation). Additionally, locating the origin at the center of mass may not always be an appropriate choice for conventional LG calculations, as molecule 2 shows a difference of as much as 25% from the origin invariant [α]_ω.

The LG(OI) approach affords the ability to compute origin invariant OR for chiral molecules without resorting to GIAOs, which are not useful for conventional CC. It is now possible to investigate more accurately the basis set dependence of the LG approach in computing this molecular property at the CC level and possibly determine a well-defined location for the origin of the coordinate system such that the conventional LG approach (which does not require the evaluation of the α^(R,P) tensor) can be directly applied. These efforts are currently under way in our lab. Given the expected better convergence of the LG results toward the basis set limit compared to the MVG approach,¹⁹ CC-LG(OI) calculations may provide very accurate estimates for the specific rotation of optically active molecules.

See the supplementary material for the optimized geometries of the test molecules, tables with the diagonal elements of the β tensors, and a table for the specific rotation of molecule 1 computed using tighter thresholds for the calculation.

The author gratefully acknowledges the support from the National Science Foundation through Grant No. CHE-1650942.

The data that support the findings of this study are available within the supplementary material, which contains the optimized geometries of the test molecules, tables with the diagonal elements of the β tensors, and a table for the specific rotation of molecule 1 computed using tighter thresholds for the calculation.