We present an origin-invariant approach to compute the full optical rotation tensor (Buckingham/Dunn tensor) in the length dipole gauge without recourse to London atomic orbitals, called LG(OI). The LG(OI) approach is simpler and less computationally demanding than the more common length gauge (LG)-London and modified velocity gauge (MVG) approaches, and it can be used with any approximate wave function or density functional method. We report an implementation at the coupled cluster with single and double excitations level (CCSD), for which we present the first simulations of the origin-invariant Buckingham/Dunn tensor in the LG. We compare LG(OI) and MVG results on a series of 22 organic molecules, showing good linear correlation between the approaches, although for small tensor elements, they provide values of opposite sign. We also attempt to decouple the effects of electron correlation and basis set incompleteness on the choice of gauge for specific rotation calculations on simple test systems. The simulations show a smooth convergence of the LG(OI) and MVG results with the basis set size toward the complete basis set limit. However, these preliminary results indicate that CCSD may not be close to a complete description of the electron correlation effects on this property even for small molecules and that basis set incompleteness may be a less important cause of discrepancy between choices of gauge than electron correlation incompleteness.

Optically active compounds are able to rotate the plane of polarization of impinging light, a phenomenon called optical rotation (OR).1 For isotropic systems (e.g., a gas or solution phase sample), only a spatially averaged OR can be measured, often expressed as a normalized quantity known as specific rotation. However, for oriented systems such as chiral crystals, one can measure the OR in a specific direction.2,3 Such measurements can only be performed on crystals because of the limited intensity of the signal, and the preparation of the sample is also difficult because the optical measurement requires a very smooth surface to distinguish the optical rotation signal (circular birefringence) from more intense signals (linear birefringence); thus, only a limited amount of experimental data are available.2,3 Simulations can, in principle, provide a direct comparison with experimental measurements by evaluating the Buckingham/Dunn optical activity tensor.4 However, only few attempts have been presented of simulations in crystals,5,6 and no general-purpose method is available. Furthermore, there are open theoretical questions that need to be addressed for the formulation of the Buckingham/Dunn in approximate methods, as discussed below. The development of accurate theoretical methods for the evaluation of this tensor could be used for comparisons with experiments to obtain structure–property relationships, assign the absolute configuration of the compound, and study the effect of intermolecular interactions on the optical activity.

Accurate quantum mechanical (QM) methods based on density functional theory (DFT) and coupled cluster (CC) theory have been developed for the calculation of chiroptical properties of molecules in isotropic media7–22 using response theory to evaluate the appropriate OR tensor.17,21–26 However, given the steep computational scaling of electronic QM calculations, all of these methods only provide an approximate solution to the Schrödinger equation using an incomplete basis expansion for the electron density. Thus, the numerical results depend on the choice of gauge for the electric dipole and quadrupole operators.17,27–29 Two typical choices are the length gauge (LG), which provides an origin-dependent tensor, and the velocity gauge (VG), with which the OR tensor is origin-invariant but has an unphysical static limit.19,26 For the LG and variational methods, such as Hartree–Fock (HF) and DFT, the origin-dependence issue is resolved using London orbitals, also known as gauge including atomic orbitals (GIAOs).8,21,30–32 However, GIAOs cannot be utilized with standard CC methods, where the reference orbitals are not reoptimized because orbital relaxation is neglected in the linear response (LR) equations to avoid unphysical poles in the LR function due to the reference wave function.23,24 Alternative formulations of CC methods with orbital optimization have been proposed,33,34 but they are not commonly used because they are computationally intensive and often have convergence issues. On the other hand, the modified VG (MVG)19,26 recipe requires the explicit evaluation and removal of the unphysical static limit, and so far, it has been the preferred approach for OR calculations at the CC level.17,18,35,36

We recently proposed a strategy to overcome the origin-dependence issue of LG calculations of specific rotation without the complication of London orbitals, an approach we called LG(OI).37 This is based on a transformation of the electric dipole–magnetic dipole polarizability tensor using the singular value decomposition (SVD) eigenvectors of the mixed-gauge electric dipole–electric dipole polarizability tensor. The LG(OI) approach is simpler than the LG-GIAO approach and faster than the MVG approach, but it shares with the latter the applicability to any approximate method. Therefore, we were able to present the first origin-invariant LG simulations of specific rotation with standard CC methods.

In this work, we extend the LG(OI) approach to the calculation of the full Buckingham/Dunn tensor and present the first origin-invariant LG-CC simulations of the full OR tensor at the single and double excitation level (CCSD). Although this approach can be applied only to individual molecules at this point, this is an important step toward the development of the corresponding method for solids using periodic boundary conditions, which, in turn, can be compared directly to experimental measurements. Here, we compare the CCSD-LG(OI) and MVG approaches on a test set of 22 organic molecules to assess whether they provide physically equivalent results. We also explore the effect of basis set and electron correlation incompleteness. We decouple these two sources of approximation using a model system where CCSD is exact and another system where CCSD is not exact, but we can use fairly large basis sets and extrapolate to the complete basis set (CBS) limit. This paper is organized as follows: the theory derivation is presented in Sec. II, details of the calculations are reported in Sec. III, the results of numerical simulations are discussed in Sec. IV, and concluding remarks are summarized in Sec. V.

In this section, we present the theory for the evaluation of the full OR tensor with the LG(OI) approach.37 In order to do that, it is useful to briefly review the equations for the Buckingham/Dunn optical activity tensor B, which is defined as4,25,38

(1)

where ϵ is the Levi–Cività operator, β is the electric dipole–magnetic dipole polarizability tensor, and A is the electric dipole–electric quadrupole polarizability tensor,

(2)
(3)

We use atomic units throughout this paper, except when otherwise specified. Since there has been some confusion in the literature for the use of the G′ symbol to indicate the electric dipole–magnetic dipole polarizability with or without the ω−1 factor, we use the notation β = −ω−1G′.25 The multipole operators in Eqs. (2) and (3) are expressed in the length dipole gauge as

(4)

which are, respectively, the electric dipole, the magnetic dipole, and the traceless electric quadrupole operators, respectively, with the position r and gradient ∇ operators implicitly summed over all the electrons of the molecule. The Greek indices denote Cartesian coordinates, ω is the frequency of the incident electromagnetic radiation, while |ψj⟩ and ωj are the jth excited state wave function and excitation frequency, respectively. These definitions are valid for non-resonant optical activity (ωjω) calculations; resonant optical activity is discussed in greater detail elsewhere.1,4,39,40

For isotropic media, the observed optical rotation is commonly reported as a normalized quantity in units of deg [dm (g/ml)]−1, known as specific rotation,

(5)

where β and ω are given in atomic units, is the reduced Planck’s constant (J s), NA is Avogadro’s number, c is the speed of light (m/s), me is the electron rest mass (kg), and M is the molecular mass (amu). The OR parameter β(ω) for isotropic media is

(6)

and we used the fact that the A contribution to B is traceless. Therefore, the A tensor is unnecessary in OR simulations of molecules in isotropic media. For an oriented system, the OR parameter β(ω) is given by the average of the components of the B tensor in Eq. (1) that are perpendicular to the direction of the beam.1 By defining a chiroptical response tensor B,25 

(7)

the OR parameter for a light beam in the direction of the unit vector n is

(8)

The tensors β and A in Eqs. (2) and (3) are both origin dependent according to

(9)
(10)

where a sum over common indices is implied in the last term on the right-hand side of both equations, O is a particular choice of origin of the coordinate system and

(11)

is a displaced origin, δ is the Kronecker delta tensor, and α is the electric dipole–electric dipole polarizability tensor expressed with two gauge representations for the dipole operator: the superscript (R, R) in Eq. (10) indicates the length gauge representation for both occurrences of the dipole operator [shown in Eq. (4)], while the superscript (R, P) in Eq. (9) indicates a mixed representation with the length gauge for one occurrence of the dipole operator and the velocity gauge for the other one,

(12)

where p is the momentum operator including an implicit summation over the number of electrons in the molecule, and there is a factor of ω−1 to be included to the mixed-gauge LR function to account for the change in representation of the electric dipole operator.19 In an exact calculation, the individual origin dependent terms for the two tensors β and A perfectly cancel out when they are combined in B in Eq. (1) such that the latter is origin invariant. However, this is no longer the case in approximate calculations, and the individual elements of B as well as its trace are origin dependent.

There are typically two ways to overcome the origin-dependence issue in approximate calculations. For the length gauge and variational methods (such as HF and DFT), one can use London orbitals8,21,30–32 (also known as GIAOs), but a such 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.23 An approach that works with any approximate method is to express the electric multipoles in the velocity gauge [see Eq. (12)] and calculate β and A tensors accordingly.19,26 The downside of this approach is that it has an unphysical static limit that needs to be evaluated explicitly and subtracted out, a procedure known as modified velocity gauge (MVG).19,26 Note that the LG-GIAO and MVG approaches provide different numerical values of the B tensor and specific rotation for approximate methods and that MVG is computationally more expensive than the LG-GIAO approach.

However, it is possible to obtain an origin-invariant version of the B tensor in the length gauge for any approximate method without using GIAOs. The approach is based on the same transformation we suggested for the specific rotation in Ref. 37, which we called LG(OI). First, we need a form of the A tensor that has a qualitatively similar origin dependence as the β vector in Eq. (9), i.e., dependent on the mixed-gauge electric dipole–electric dipole polarizability α(R,P). This is accomplished by using the following form of the quadrupole operator:

(13)

which is also used in the MVG approach to form the A tensor (albeit in conjunction with the velocity electric dipole).26 With this choice of quadrupole operator,

(14)

such that the A tensor uses a mixed-gauge representation: length for the dipole operator and velocity for the quadrupole operator. This form of the A tensor preserves the expression of the B tensor in Eq. (1) but leads to the following expression of its origin dependence:

(15)

which is now compatible with that of the β tensor in Eq. (9). However, the use of the mixed-gauge form of the A tensor is not sufficient to lead to an origin-invariant B tensor because the α(R,P) tensor is not symmetric. Therefore, the contributions of the β and A(R,P) tensors to the B tensor do not perfectly cancel out. To achieve that, we diagonalize the α(R,P) tensor with a singular value decomposition (SVD),

(16)

where αD(R,P) is diagonal and U and V are unitary transformations. We then apply the inverse transformation to the β and A(R,P) tensors,

(17)
(18)

Using the β̃ and Ã(R,P) tensors in Eq. (1) makes the B tensor fully origin invariant without recourse to GIAOs. The transformations in Eqs. (17) and (18) are sensitive to the phase of the eigenvectors and their order in the unitary matrices. Therefore, it is important to avoid inadvertently changing the handedness of the coordinate system with the transformation. One can argue that the use of the quadrupole operator in Eq. (13) makes the overall approach not purely “length gauge” but rather “mixed gauge.” However, we prefer to maintain the name LG(OI) because the isotropic specific rotation is still computed in the LG and to avoid introducing more acronyms for essentially the same theoretical approach.

We know that the LG and VG approaches are only equivalent for an exact calculation and that approximations coming from electron correlation and basis set incompleteness in practical simulations lead to different numerical results with the two gauge choices. However, the effects of the two types of incompleteness are difficult to disentangle. Nevertheless, it would be desirable to have a criterion to compare the quality of two different calculations on the same system, at least within the same choice of gauge. For the velocity gauge, a criterion could be the relative magnitude of the unphysical static limit compared to that of the MVG result. For the LG approach, a criterion could be based on the lack of origin invariance of the specific rotation, which is related to the degree of symmetry of the mixed-gauge polarizability tensor α(R,P). This degree of symmetry can be expressed as

(19)

where ‖·‖F represents the Frobenius norm of a matrix and the subscript A indicates the anti-symmetric part of the α(R,P) tensor. The systematic improvement of a chosen model chemistry should correspond to an increase in Δs toward 1.

A legitimate question about the LG(OI) B tensor is what orientation of the molecule corresponds to this tensor values. In an exact calculation, α(R,P) is symmetric and UV in Eq. (16). Therefore, the transformed B tensor corresponds to an orientation where the α(R,P) tensor is diagonal. For approximate calculations, the closest we can get to this orientation is to use the eigenvectors that diagonalize the symmetric part of the α(R,P) tensor, which we collect in the unitary matrix W. As we approach the exact solution, U, VW. Once the molecule is reoriented according to the rotation defined by the W matrix, we can analytically calculate the origin displacement vector that makes the diagonal elements of the LG β tensor equal to the corresponding elements of the LG(OI) β̃ tensor so that they provide the same value of specific rotation. Using Eq. (9), the components of this particular displacement vector d are the solution of the following linear system of equations:

(20)

The superscript S serves as a reminder that the molecule is oriented according to the W rotation matrix, and we neglected the (R, P) superscript on the polarizability tensor elements for clarity. The LG B tensor obtained with this choice of molecular coordinates is the closest to the LG(OI) result and can be considered the best LG result obtainable for a particular approximate model chemistry.

In terms of computational cost, the LG(OI) approach is more efficient than both the LG-GIAO and MVG approaches because no extra terms due to GIAOs are necessary and there is no static limit to evaluate explicitly. All tensors can be computed using standard linear response theory, where a perturbed density is evaluated due to length gauge electric dipole perturbation and then contracted with the appropriate multipole integrals to obtain β, A(R,P), and α(R,P).17,21–25 After these tensors are available, the LG(OI) transformations in Eqs. (16)(18) are trivial. For the implementation at the CC level, one would need to use the non-symmetric form of the LR function, which requires the evaluation of the perturbed T and Λ amplitudes for one perturbation.24,41 In that case, the LG(OI) method is completely equivalent in cost to a standard LG calculation and less demanding than an MVG calculation. However, in our current implementation in GAUSSIAN42 at the CCSD level, only the symmetric form of the LR function is available, which requires the evaluation of the perturbed T amplitudes for all perturbations.24,41 This route makes the LG(OI) and MVG approaches equivalent in cost. Nevertheless, this is an issue of the implementation and not of the method per se, and it will be resolved in future versions of the software.

All calculations were performed with a development version of the GAUSSIAN suite of programs.42 The tensors were computed using standard LR theory at the CCSD level with frozen core orbitals for the correlation energy evaluation and frozen orbitals for the LR function. The aug-cc-pVDZ43,44 basis set was used for the 22 test molecules, taken from Ref. 26. For molecules 14, 21, and 22, geometries were taken from Ref. 37, but reoriented as described in Sec. II. For the remaining molecules, we simply used the geometries from Ref. 26 and rotated the MVG B tensor directly using the W rotation matrix. For the two-electron model system, the aug-cc-pVXZ43–45 and aug-mcc-pVYZ46,47 basis sets were used, with X = 2–6 and Y = 7, 8 (basis sets with X = 7, 8 were not available on the basis set exchange for H).48–50 For (−)-hydrogen peroxide, a series of Dunning and Pople51–55 basis sets are used as detailed in Sec. IV. Results are presented with the LG, LG(OI), and MVG approaches at the sodium D line, i.e. 589.3 nm, because it is far from resonance for all systems.

We start this section with a comparison of specific rotation calculations between the LG(OI) and MVG approaches to investigate the numerical differences induced by the use of approximate electronic structure methods. We use [α]D initially because it is a more convenient quantity than the whole OR tensor for this analysis. As discussed in Sec. II, the difference between the LG and VG results is due to approximations in the level of theory used to compute the electron density (i.e., electron correlation incompleteness) and to the incompleteness of practical basis sets. These two factors are difficult to disentangle in calculations on typical molecules. Therefore, we first use a toy two-electron model, H42+, for which CCSD provides the exact answer for a non-relativistic molecular Hamiltonian within the Born–Oppenheimer approximation (we also neglect vibrational contributions). With this model, we can directly study the effect of the basis set on the specific rotation. We built the model by using the H2O2 geometry37 and replacing the O’s with H’s and by making two somewhat random choices of origin: O1 is the origin where the LG and LG(OI) approaches provide the same result for the specific rotation of H2O2 with the aug-cc-pVDZ basis set, and it is located somewhere inside the molecule, while O2 is displaced by −5 Å in every Cartesian direction from O1. Obviously, only the LG results are affected by the choice of origin, and O1/O2 represent two points with no special meaning for this system.

In Fig. 1, we report the change in [α]D for H42+ with the increasing size of the basis set for all choices of gauge; this figure also reports the total CCSD energy and the Δs factor defined in Eq. (19). All numerical values are reported in Tables S1 and S2 of the supplementary material. The figure shows a smooth convergence of the specific rotation with the basis set size for all choices of gauge, except for VG. For the latter, there is an oscillatory behavior that starts from a value with the wrong sign for ζ = 2 [−290 deg dm−1 (g/ml)−1]; the unphysical static limit is still >1 deg with the ζ = 5 basis set. The LG(OI) and MVG values are always within 1 deg from each other and are essentially converged with the ζ = 5 basis set. This indicates that basis set incompleteness is not a major cause of discrepancy between choices of gauge if the origin-dependence issue of the LG approach and the unphysical static limit issue of the VG approach are properly addressed. The LG approach converges to the LG(OI) value for the ζ = 5 basis set even with a significantly displaced geometry (O2). The convergence of the specific rotation follows the convergence of the energy [compare panels (a) and (b) of the figure], which supports that similar extrapolation formulas can be used to estimate the complete basis set (CBS) limit.36 The Δs factor also increases until reaching a plateau with the ζ = 5 basis set consistently with the convergence of [α]D.

FIG. 1.

Data computed for H42+ with all choices of gauge and the aug-cc-pVXZ (X = 2–6) and aug-mcc-pVYZ (Y = 7, 8) basis sets at CCSD level; for LG, data are presented at two choices of the origin, O1 and O2. Panel (a) Specific rotation at the sodium D line in deg dm−1 (g/ml)−1. Panel (b) Combined plot of the CCSD energy (ECC in a.u.) and Δs factor [Eq. (19)].

FIG. 1.

Data computed for H42+ with all choices of gauge and the aug-cc-pVXZ (X = 2–6) and aug-mcc-pVYZ (Y = 7, 8) basis sets at CCSD level; for LG, data are presented at two choices of the origin, O1 and O2. Panel (a) Specific rotation at the sodium D line in deg dm−1 (g/ml)−1. Panel (b) Combined plot of the CCSD energy (ECC in a.u.) and Δs factor [Eq. (19)].

Close modal

A more realistic but still small system to investigate the effect of electron correlation incompleteness on the specific rotation is (−)-hydrogen peroxide, whose geometry is taken from Ref. 37. Figure 2 reports the [α]D values computed with the LG(OI) and MVG approaches with Dunning basis sets of increasing angular momentum functions up to ζ = 5 and estimates of the CBS limit; the latter are computed with a two-point inverse-power extrapolation formula from the ζ = 4 and ζ = 5 values with exponent n = 5, as suggested in Ref. 36. Table I reports the values in Fig. 2 together with [α]D and Δs computed with the VG, MVG, and LG(OI) approaches for the same Dunning basis sets43–47 and a series of Pople basis sets51–55 of increasing size. From a comparison of the data in Fig. 2 and Table I, the specific rotation smoothly decreases in magnitude for both the MVG and LG(OI) approaches as the size of the Dunning basis sets increases. Simultaneously, the VG unphysical static limit decreases in magnitude, while the Δs factor increases, as one would expect. However, the two choices of gauge do not converge to the same CBS limit, and the difference (∼40 deg) is twice as large as the magnitude of the LG(OI) [α]D and two thirds the magnitude of the MVG [α]D. The MVG–VG difference for the CBS estimate is even larger (∼50 deg). These results indicate that CCSD is not quite converged in terms of electron correlation for this molecule.

FIG. 2.

Specific rotation [deg dm−1 (g/ml)−1] at the sodium D line for H2O2 computed with the LG(OI) and MVG approaches at the CCSD/aug-cc-pVζZ level. The CBS values are obtained with a two-point extrapolation from the ζ = 4 and ζ = 5 values (see text for details).

FIG. 2.

Specific rotation [deg dm−1 (g/ml)−1] at the sodium D line for H2O2 computed with the LG(OI) and MVG approaches at the CCSD/aug-cc-pVζZ level. The CBS values are obtained with a two-point extrapolation from the ζ = 4 and ζ = 5 values (see text for details).

Close modal
TABLE I.

Values of [α]D (deg dm−1 (g/ml)−1) and Δs for H2O2 computed with the VG, MVG, and LG(OI) approaches at the CCSD level and a series of Dunning and Pople basis sets.

[α]DΔs
MVGVGLG(OI)LG(OI)
CBS(5Z-4Z) −56 −5 −19 ⋯ 
aug-cc-pV5Z −61 −32 −22 0.9989 
aug-cc-pVQZ −70 −86 −30 0.9988 
aug-cc-pVTZ −92 −272 −45 0.9986 
aug-cc-pVDZ −140 −899 −80 0.9983 
6-311++G(3df,3pd) −127 −756 −57 0.9995 
6-311++G(2d,2p) −169 −2677 −98 0.9992 
6-311++G(d,p) −141 −3756 −65 0.9996 
6-31++G(d,p) −138 −2942 −39 0.9988 
[α]DΔs
MVGVGLG(OI)LG(OI)
CBS(5Z-4Z) −56 −5 −19 ⋯ 
aug-cc-pV5Z −61 −32 −22 0.9989 
aug-cc-pVQZ −70 −86 −30 0.9988 
aug-cc-pVTZ −92 −272 −45 0.9986 
aug-cc-pVDZ −140 −899 −80 0.9983 
6-311++G(3df,3pd) −127 −756 −57 0.9995 
6-311++G(2d,2p) −169 −2677 −98 0.9992 
6-311++G(d,p) −141 −3756 −65 0.9996 
6-31++G(d,p) −138 −2942 −39 0.9988 

On the other hand, the Pople basis sets provide the wrong trend of [α]D with the basis set size, as the property increases in magnitude until the largest basis set is used. This poor performance is accompanied by very large values of the static limit for the VG approach, especially for the smaller sets in the series. In contrast, the Δs factors are as large as or larger than those for the more accurate Dunning basis sets. In other words, the calculations with the Pople basis sets suffer from a smaller origin dependence of the LG results compared to those with the Dunning basis sets even if the latter provide overall more reliable results based on the trend with the basis set size.

We now report calculations of the full B tensor with the LG(OI) and MVG approaches for the test set of 22 organic molecules we compiled in Ref. 26. First, we use molecules 14, 21, and 22 from the original LG(OI) paper37 to demonstrate that the LG(OI) B tensor is indeed origin invariant. These results are shown in Tables IIIV, where the B tensor values and the corresponding contributions from the β̃ and à tensors are reported at the same geometries of Ref. 37, i.e., with the origin in the center of mass and displaced by −1000 Å in every Cartesian direction. These data confirm that while the β̃ and à tensors are origin dependent, they combine to make a B tensor that is origin invariant (within numerical precision).

TABLE II.

Values of the β̃, Ã, and B tensors (a.u.) computed with the LG(OI) approach for molecule 14 at the CCSD/aug-cc-pVDZ level. The origin is located at the center of mass (O) and displaced by −1000 Å in every direction (O′); the geometries are taken from Ref. 37. The differences between the O and O′ values for the B tensors elements are due to the numerical precision of the calculations, as discussed in Ref. 37.

OO′
β̃ÃBβ̃ÃB
xx 4.7267 −1.5875 3.1392 4.7281 −1.5872 3.1408 
yy −7.2375 1.5638 −5.6737 −7.2345 1.5610 −5.6735 
zz −3.5626 0.0238 −3.5388 −3.5590 0.0262 −3.5328 
xy 0.0060 −0.0550 −0.0491 −1647.0100 1646.9600 −0.0479 
xz −1.2853 0.4307 −0.8546 3902.4600 −3903.3200 −0.8512 
yz 1.2498 −0.2015 1.0483 −1530.1500 1531.2000 1.0516 
OO′
β̃ÃBβ̃ÃB
xx 4.7267 −1.5875 3.1392 4.7281 −1.5872 3.1408 
yy −7.2375 1.5638 −5.6737 −7.2345 1.5610 −5.6735 
zz −3.5626 0.0238 −3.5388 −3.5590 0.0262 −3.5328 
xy 0.0060 −0.0550 −0.0491 −1647.0100 1646.9600 −0.0479 
xz −1.2853 0.4307 −0.8546 3902.4600 −3903.3200 −0.8512 
yz 1.2498 −0.2015 1.0483 −1530.1500 1531.2000 1.0516 
TABLE III.

Values of the β̃, Ã, and B tensors (a.u.) computed with the LG(OI) approach for molecule 21 at the CCSD/aug-cc-pVDZ level. The origin is located at the center of mass (O) and displaced by −1000 Å in every direction (O′); the geometries are taken from Ref. 37. The differences between the O and O′ values for the B tensors elements are due to the numerical precision of the calculations, as discussed in Ref. 37.

OO′
β̃ÃBβ̃ÃB
xx −0.4291 0.0059 −0.4232 −0.4286 0.0054 −0.4232 
yy −0.2843 0.3912 0.1070 −0.2846 0.3915 0.1070 
zz 0.5021 −0.3971 0.1049 0.5023 −0.3969 0.1054 
xy 0.8011 0.0202 0.8213 −1349.0100 1349.8300 0.8212 
xz 0.0000 0.0000 0.0000 1777.0900 −1777.0900 0.0002 
yz 0.0000 0.0000 0.0000 −111.0810 111.0810 0.0003 
OO′
β̃ÃBβ̃ÃB
xx −0.4291 0.0059 −0.4232 −0.4286 0.0054 −0.4232 
yy −0.2843 0.3912 0.1070 −0.2846 0.3915 0.1070 
zz 0.5021 −0.3971 0.1049 0.5023 −0.3969 0.1054 
xy 0.8011 0.0202 0.8213 −1349.0100 1349.8300 0.8212 
xz 0.0000 0.0000 0.0000 1777.0900 −1777.0900 0.0002 
yz 0.0000 0.0000 0.0000 −111.0810 111.0810 0.0003 
TABLE IV.

Values of the β̃, Ã, and B tensors (a.u.) computed with the LG(OI) approach for molecule 22 at the CCSD/aug-cc-pVDZ level. The origin is located at the center of mass (O) and displaced by −1000 Å in every direction (O′); the geometries are taken from Ref. 37. The differences between the O and O′ values for the B tensors elements are due to the numerical precision of the calculations, as discussed in Ref. 37.

OO′
β̃ÃBβ̃ÃB
xx −1.3295 −1.4304 −2.7599 −1.331 4 −1.428 5 −2.7599 
yy −0.3650 0.0034 −0.3616 −0.365 5 0.003 9 −0.3616 
zz 2.1746 1.4270 3.6015 2.172 2 1.424 5 3.5967 
xy 6.3784 −1.0017 5.3767 10 863.600 0 −10 858.200 0 5.3767 
xz 0.0000 0.0000 0.0000 −81.456 7 81.454 0 −0.0028 
yz 0.0000 0.0000 0.0000 −12 407.700 0 12 407.700 0 −0.0020 
OO′
β̃ÃBβ̃ÃB
xx −1.3295 −1.4304 −2.7599 −1.331 4 −1.428 5 −2.7599 
yy −0.3650 0.0034 −0.3616 −0.365 5 0.003 9 −0.3616 
zz 2.1746 1.4270 3.6015 2.172 2 1.424 5 3.5967 
xy 6.3784 −1.0017 5.3767 10 863.600 0 −10 858.200 0 5.3767 
xz 0.0000 0.0000 0.0000 −81.456 7 81.454 0 −0.0028 
yz 0.0000 0.0000 0.0000 −12 407.700 0 12 407.700 0 −0.0020 

For the same three molecules, Tables VVII report the LG, LG(OI), and MVG results computed with the same geometries but rotated and translated as discussed in Sec. II (the geometries at these orientations are reported in Tables S3–S5 of the supplementary material). These data confirm that the orientation obtained from the W rotation and the translation in Eq. (20) provides very good agreement between the LG and LG(OI) results, and it is the right orientation to compare the LG(OI) and MVG B tensors. Therefore, as mentioned in Sec. III, for the remaining molecules in the test set, we use the geometries from Ref. 26 for the LG(OI) calculations, and we simply rotate the MVG tensors presented in that work with the W matrix for the comparison with the LG(OI) results. Note that the LG(OI) and MVG approaches provide tensor elements of the same sign in almost all cases for the three molecules (see Tables VVII) The only exceptions are element xy for molecule 14 and element yy for molecule 21. For molecules 14 and 22, MVG tends to provide elements that are smaller in magnitude than those obtained with LG(OI), while the trend is opposite for 21. However, the % Diff values in Tables VVII show no clear trend across the various elements, and the differences in specific rotation (i.e., in the trace of the B tensor) are determined by different cancellation effects between the diagonal tensor elements.

TABLE V.

Values of [α]D [deg dm−1 (g/ml)−1] and of the B tensor elements (a.u.) computed with the LG, LG(OI), and MVG approaches for molecule 14 at the CCSD/aug-cc-pVDZ level. The last column (% Diff) reports the % difference between the LG(OI) and MVG results.

LGLG (OI)MVG% Diff
[α]D −724 −724 −549 −24.1 
xx 3.1338 3.1392 2.4325 −22.5 
yy −5.6727 −5.6737 −4.5313 −20.1 
zz −3.5344 −3.5388 −2.5111 −29.0 
xy −0.0496 −0.0491 0.0574 −216.9 
xz −0.8620 −0.8546 −0.4874 −43.0 
yz 1.0431 1.0483 0.4843 −53.8 
LGLG (OI)MVG% Diff
[α]D −724 −724 −549 −24.1 
xx 3.1338 3.1392 2.4325 −22.5 
yy −5.6727 −5.6737 −4.5313 −20.1 
zz −3.5344 −3.5388 −2.5111 −29.0 
xy −0.0496 −0.0491 0.0574 −216.9 
xz −0.8620 −0.8546 −0.4874 −43.0 
yz 1.0431 1.0483 0.4843 −53.8 
TABLE VI.

Values of [α]D [deg dm−1 (g/ml)−1] and of the B tensor elements (a.u.) computed with the LG, LG(OI), and MVG approaches for molecule 21 at the CCSD/aug-cc-pVDZ level. The last column (% Diff) reports the % difference between the LG(OI) and MVG results.

LGLG (OI)MVG% Diff
[α]D −80 −80 −140 74.3 
xx −0.4202 −0.4232 −0.5348 26.4 
yy 0.1096 0.1070 −0.0184 −117.2 
zz 0.0993 0.1049 0.1849 76.2 
xy 0.8200 0.8213 0.8813 7.3 
xz 0.0000 0.0000 0.0000 0.0 
yz 0.0000 0.0000 0.0000 0.0 
LGLG (OI)MVG% Diff
[α]D −80 −80 −140 74.3 
xx −0.4202 −0.4232 −0.5348 26.4 
yy 0.1096 0.1070 −0.0184 −117.2 
zz 0.0993 0.1049 0.1849 76.2 
xy 0.8200 0.8213 0.8813 7.3 
xz 0.0000 0.0000 0.0000 0.0 
yz 0.0000 0.0000 0.0000 0.0 
TABLE VII.

Values of [α]D [deg dm−1 (g/ml)−1] and of the B tensor elements (a.u.) computed with the LG, LG(OI), and MVG approaches for molecule 22 at the CCSD/aug-cc-pVDZ level. The last column (% Diff) reports the % difference between the LG(OI) and MVG results.

LGLG(OI)MVG% Diff
[α]D 91 91 133 46.4 
xx −2.7852 −2.7599 −2.6363 −4.5 
yy −0.4039 −0.3616 −0.1883 −47.9 
zz 3.6692 3.6015 3.5272 −2.1 
xy 5.3826 5.3767 5.1156 −4.9 
xz −0.0007 0.0000 0.0000 0.0 
yz −0.0076 0.0000 0.0000 0.0 
LGLG(OI)MVG% Diff
[α]D 91 91 133 46.4 
xx −2.7852 −2.7599 −2.6363 −4.5 
yy −0.4039 −0.3616 −0.1883 −47.9 
zz 3.6692 3.6015 3.5272 −2.1 
xy 5.3826 5.3767 5.1156 −4.9 
xz −0.0007 0.0000 0.0000 0.0 
yz −0.0076 0.0000 0.0000 0.0 

The B tensor elements and the trace for the entire set of 22 molecules are compared in the correlation plots in Fig. 3. In these plots, the tensor elements and the trace are divided by the molecular mass to produce an intensive quantity, which allows us to compare molecules of various sizes on equal footing. We performed linear fits on the data both with and without molecule 14 [i.e., the notorious (S)-(−)-norbornenone], as in Ref. 26, and the slope and R2 parameters for these fits are plotted in Fig. 4. It is clear from Figs. 3 and 4 that there is a strong linear correlation between the MVG and LG(OI) B tensor components, with the R2 of each fit approaching 1.0. This strong correlation persists regardless of whether molecule 14 is included in the fit; however, this molecule has an outsized influence on the predicted slope. While the fit with molecule 14 suggests that MVG B components consistently underestimate the corresponding LG(OI) tensor elements by 10%–20%, the fit excluding 14 predicts a consistent overestimation relative to LG(OI) for all components. In this case, the overestimation is smaller in magnitude (<5%), except for the yy (∼10%) and yz (15%) elements. The MVG trace is consistently smaller than that with LG(OI), but again removing molecule 14 from the fit reduces this underestimation to only about 5%. The R2 parameters are all >0.95, and all but the yy component increase when molecule 14 is excluded from the fitting. The MVG yy elements for molecules 21 and 22 significantly underestimate their LG(OI) counterparts, with the component for 22 even changing sign. This underestimation leads to a slightly larger R2 in the fit with 14, since the slope is <1.0 in this case. The large effect of including/excluding molecule 14 suggests that this level of theory and basis set are somewhat far from convergence for this compound.

FIG. 3.

MVG vs LG(OI) correlation plots for the mass-normalized elements and trace (a.u. × 103) of the B tensor computed at the CCSD/aug-cc-pVDZ level for the molecules from Ref. 26. The linear fits were performed with (cyan) and without (yellow) molecule 14. Molecule 14 is denoted by the unfilled marker in each of the plots.

FIG. 3.

MVG vs LG(OI) correlation plots for the mass-normalized elements and trace (a.u. × 103) of the B tensor computed at the CCSD/aug-cc-pVDZ level for the molecules from Ref. 26. The linear fits were performed with (cyan) and without (yellow) molecule 14. Molecule 14 is denoted by the unfilled marker in each of the plots.

Close modal
FIG. 4.

Slope and R2 correlation parameters of the linear fits of B components including (left) and excluding (right) molecule 14.

FIG. 4.

Slope and R2 correlation parameters of the linear fits of B components including (left) and excluding (right) molecule 14.

Close modal

Despite the general accuracy of the fits, there are still individual cases for which the MVG and LG(OI) B tensor elements differ qualitatively: molecules 7 (xx), 8 (yz), 14 (xy), 15 (xz), and 21 (yy) all have one MVG B element that differs in sign from the LG(OI) B tensor. These sign changes are not specific to any particular Cartesian component and tend to occur in cases where the absolute value of the LG(OI) B elements are small (<0.10 a.u). These small values lead to relatively small squared errors, so the R2 parameters of the fits remain close to 1.0. Nevertheless, it is not obvious whether these deviations stem from incompleteness of the basis or insufficient treatment of correlation; this will be explored further in a future study.

In this work, we present the theory for the origin-invariant formulation of the full OR tensor using the length-dipole gauge without recourse to London atomic orbitals (also known as GIAOs). This approach, LG(OI), is based on the use of the mixed-gauge electric quadrupole operator [see Eq. (13)] to evaluate the electric dipole–electric quadrupole tensor A(R,P) and the subsequent transformation of this tensor using the SVD eigenvectors of the mixed-gauge electric dipole–electric dipole polarizability tensor α(R,P) [see Eq. (15)]. The LG(OI) approach is considerably less complex to implement than the common LG-GIAO approach used in HF and DFT, and it is currently the only option for origin-invariant calculations with standard LR-CC methods in the length gauge. Simultaneously, the LG(OI) approach is, in principle, computationally more efficient than the MVG approach because there is no unphysical static limit.

In Sec. IV, we attempt to decouple the effect of electron correlation and basis set incompleteness in the evaluation of molecular specific rotation. Using a two-electron model system for which CCSD is exact, we show that the LG(OI) and MVG approaches are equivalent even for incomplete basis sets and they smoothly converge to the CBS limit. The standard (origin-dependent) LG method also converges smoothly to the CBS limit independently of the choice of origin, but this is not the case for the VG approach where the unphysical static limit shows an erratic behavior with different basis sets (see Fig. 1). The H2O2 example in Fig. 2 shows that the LG(OI) and MVG approaches smoothly converge toward the CBS limit when using an approximate electronic structure method (CCSD). However, the two approaches converge to a significantly different CBS limit, and the VG unphysical static limit is also comparatively large (see Table I), indicating that CCSD is surprisingly far from convergence in the description of electron correlation even for such a small molecule. A comparison with the results for the two-electron model system suggests that electron correlation incompleteness may be a more important cause of discrepancy between choices of gauge than basis set incompleteness. We are currently performing a thorough basis set study on multiple test molecules to explore this behavior further. Furthermore, a comparison between results obtained with Dunning and Pople basis sets in Table I indicates that while the symmetry factor Δs [see Eq. (19)] systematically increases toward one with well-balanced basis sets (e.g., Dunning’s), its value cannot be easily used as a metric for the numerical convergence of the calculations between basis sets of different type. In fact, Δs is larger when using Pople basis sets than with Dunning basis sets even if the latter provide more reliable results.

The results for the full OR tensor B computed with LG(OI) and MVG in Tables VVII indicate that the two approaches tend to provide qualitatively similar results. This is further supported by Figs. 3 and 4, which show a clear linear correlation between the mass-scaled MVG and LG(OI) tensor elements for a wider range of molecules. However, there remain unexplained discrepancies between the two gauges. In some cases, like the xy element for molecule 14 in Table II and the yy element for molecule 21 in Table III, the LG(OI) and MVG values have opposite signs. Large component differences, like those seen for the yy element of molecule 14 in Fig. 3, can have a pronounced effect on the obtained linear fits for these components, even when the sign is consistent. It is unclear at this point whether the deviations of these components are solely due to the small basis set used for these calculations (aug-cc-pVDZ) or if the level of theory plays a role. We will further investigate the effect of the basis set in a follow-up study.

In summary, the LG(OI) approach expresses the Buckingham/Dunn tensor in an origin-invariant formulation that is uniform across any level of theory, and it can be used for the evaluation of the optical rotation in oriented systems such as chiral crystals, providing a direct comparison with experimental measurements.

See the supplementary material for the numerical values of the quantities plotted in Fig. 1 of the main text (Tables S1 and S2), the Cartesian coordinates of molecules 14, 21, and 22 corresponding to the tensor values in Tables VVII of the main text (Tables S3–S5), and the LG(OI) and MVG B tensors for the remaining molecules (Tables S6–S24).

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

The data that support the findings of this study are available within the article and its supplementary material.

1.
L. D.
Barron
,
Molecular Light Scattering and Optical Activity
, 2nd ed. (
Cambridge University Press
,
2004
).
2.
W.
Kaminsky
,
K.
Claborn
, and
B.
Kahr
,
Chem. Soc. Rev.
33
,
514
(
2004
).
3.
J. H.
Freudenthal
,
E.
Hollis
, and
B.
Kahr
,
Chirality
21
,
E20
(
2009
).
4.
A. D.
Buckingham
and
M. B.
Dunn
,
J. Chem. Soc. A
1971
,
1988
.
5.
H.
Zhong
,
Z. H.
Levine
,
D. C.
Allan
, and
J. W.
Wilkins
,
Phys. Rev. Lett.
69
,
379
(
1992
).
6.
H.
Zhong
,
Z. H.
Levine
,
D. C.
Allan
, and
J. W.
Wilkins
,
Phys. Rev. B
48
,
1384
(
1993
).
7.
P. L.
Polavarapu
,
Mol. Phys.
91
,
551
(
1997
).
8.
J. R.
Cheeseman
,
M. J.
Frisch
,
F. J.
Devlin
, and
P. J.
Stephens
,
J. Phys. Chem. A
104
,
1039
(
2000
).
9.
S.
Grimme
,
Chem. Phys. Lett.
339
,
380
(
2001
).
10.
P. L.
Polavarapu
,
Chirality
14
,
768
(
2002
).
11.
J.
Autschbach
,
S.
Patchkovskii
,
T.
Ziegler
,
S. J. a.
van Gisbergen
, and
E.
Jan Baerends
,
J. Chem. Phys.
117
,
581
(
2002
).
12.
K.
Ruud
and
T.
Helgaker
,
Chem. Phys. Lett.
352
,
533
(
2002
).
13.
P. J.
Stephens
,
D. M.
McCann
,
J. R.
Cheeseman
, and
M. J.
Frisch
,
Chirality
17
,
52
(
2005
).
14.
S.
Grimme
,
A.
Bahlmann
, and
G.
Haufe
,
Chirality
14
,
793
(
2002
).
15.
K.
Ruud
,
P. J.
Stephens
,
F. J.
Devlin
,
P. R.
Taylor
,
J. R.
Cheeseman
, and
M. J.
Frisch
,
Chem. Phys. Lett.
373
,
606
(
2003
).
16.
M. C.
Tam
,
N. J.
Russ
, and
T. D.
Crawford
,
J. Chem. Phys.
121
,
3550
(
2004
).
17.
T. D.
Crawford
,
Theor. Chem. Acc.
115
,
227
(
2006
).
18.
T. D.
Crawford
and
P. J.
Stephens
,
J. Phys. Chem. A
112
,
1339
(
2008
).
19.
T. B.
Pedersen
,
H.
Koch
,
L.
Boman
, and
A. M. J.
Sánchez de Merás
,
Chem. Phys. Lett.
393
,
319
(
2004
).
20.
T. D.
Crawford
,
L. S.
Owens
,
M. C.
Tam
,
P. R.
Schreiner
, and
H.
Koch
,
J. Am. Chem. Soc.
127
,
1368
(
2005
).
21.
M.
Krykunov
and
J.
Autschbach
,
J. Chem. Phys.
123
,
114103
(
2005
).
22.
23.
H.
Koch
and
P.
Jørgensen
,
J. Chem. Phys.
93
,
3333
(
1990
).
24.
O.
Christiansen
,
P.
Jørgensen
, and
C.
Hättig
,
Int. J. Quantum Chem.
68
,
1
(
1998
).
25.
26.
K.
Zhang
,
T.
Balduf
, and
M.
Caricato
,
Chirality
33
,
303
(
2021
).
27.
S.
Pelloni
and
P.
Lazzeretti
,
J. Chem. Phys.
140
,
074105
(
2014
).
28.
P.
Lazzeretti
,
Int. J. Quantum Chem.
114
,
1364
(
2014
).
29.
M. C.
Caputo
,
S.
Pelloni
, and
P.
Lazzeretti
,
Int. J. Quantum Chem.
115
,
900
(
2015
).
31.
32.
T.
Helgaker
,
K.
Ruud
,
K. L.
Bak
,
P.
Jørgensen
, and
J.
Olsen
,
Faraday Discuss.
99
,
165
(
1994
).
33.
T. B.
Pedersen
,
B.
Fernández
, and
H.
Koch
,
J. Chem. Phys.
114
,
6983
(
2001
).
34.
G. D.
Lindh
,
T. J.
MacH
, and
T. D.
Crawford
,
Chem. Phys.
401
,
125
(
2012
).
35.
T. J.
Mach
and
T. D.
Crawford
,
J. Phys. Chem. A
115
,
10045
(
2011
).
36.
S.
Haghdani
,
P.-O.
Åstrand
, and
H.
Koch
,
J. Chem. Theory Comput.
12
,
535
(
2016
).
37.
M.
Caricato
,
J. Chem. Phys.
153
,
151101
(
2020
).
38.
T.
Balduf
and
M.
Caricato
,
J. Phys. Chem. C
123
,
4329
(
2019
).
39.
P.
Norman
,
K.
Ruud
, and
T.
Helgaker
,
J. Chem. Phys.
120
,
5027
(
2004
).
40.
M.
Krykunov
and
J.
Autschbach
,
J. Chem. Phys.
125
,
034102
(
2006
).
41.
J. F.
Stanton
and
J.
Gauss
,
Int. Rev. Phys. Chem.
19
,
61
(
2000
).
42.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
G. A.
Petersson
,
H.
Nakatsuji
,
X.
Li
,
M.
Caricato
,
A. V.
Marenich
,
J.
Bloino
,
B. G.
Janesko
,
R.
Gomperts
,
B.
Mennucci
,
H. P.
Hratchian
,
J. V.
Ortiz
,
A. F.
Izmaylov
,
J. L.
Sonnenberg
,
D.
Williams-Young
,
F.
Ding
,
F.
Lipparini
,
F.
Egidi
,
J.
Goings
,
B.
Peng
,
A.
Petrone
,
T.
Henderson
,
D.
Ranasinghe
,
V. G.
Zakrzewski
,
J.
Gao
,
N.
Rega
,
G.
Zheng
,
W.
Liang
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
K.
Throssell
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M. J.
Bearpark
,
J. J.
Heyd
,
E. N.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
T. A.
Keith
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A. P.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
J. M.
Millam
,
M.
Klene
,
C.
Adamo
,
R.
Cammi
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
O.
Farkas
,
J. B.
Foresman
, and
D. J.
Fox
, Gaussian Development Version Revision j.11,
2020
.
43.
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
90
,
1007
(
1989
).
44.
R. A.
Kendall
,
T. H.
Dunning
, and
R. J.
Harrison
,
J. Chem. Phys.
96
,
6796
(
1992
).
45.
K. A.
Peterson
,
D. E.
Woon
, and
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
100
,
7410
(
1994
).
46.
S. L.
Mielke
,
B. C.
Garrett
, and
K. A.
Peterson
,
J. Chem. Phys.
116
,
4142
(
2002
).
47.
S. L.
Mielke
,
D. W.
Schwenke
, and
K. A.
Peterson
,
J. Chem. Phys.
122
,
224313
(
2005
).
48.
B. P.
Pritchard
,
D.
Altarawy
,
B.
Didier
,
T. D.
Gibson
, and
T. L.
Windus
,
J. Chem. Inf. Model.
59
,
4814
(
2019
).
50.
K. L.
Schuchardt
,
B. T.
Didier
,
T.
Elsethagen
,
L.
Sun
,
V.
Gurumoorthi
,
J.
Chase
,
J.
Li
, and
T. L.
Windus
,
J. Chem. Inf. Model.
47
,
1045
(
2007
).
51.
R.
Ditchfield
,
W. J.
Hehre
, and
J. A.
Pople
,
J. Chem. Phys.
54
,
724
(
1971
).
52.
P. C.
Hariharan
and
J. A.
Pople
,
Theor. Chim. Acta
28
,
213
(
1973
).
53.
M. J.
Frisch
,
J. A.
Pople
, and
J. S.
Binkley
,
J. Chem. Phys.
80
,
3265
(
1984
).
54.
R.
Krishnan
,
J. S.
Binkley
,
R.
Seeger
, and
J. A.
Pople
,
J. Chem. Phys.
72
,
650
(
1980
).
55.
T.
Clark
,
J.
Chandrasekhar
,
G. W.
Spitznagel
, and
P. V. R.
Schleyer
,
J. Comput. Chem.
4
,
294
(
1983
).

Supplementary Material