Derivation and implementation of the optical rotation tensor for chiral crystals

Measurements of optical rotation (OR) have long been used to characterize chiral systems, but it was only relatively recently that the first Hartree–Fock level simulations were performed to compute this property.¹ Great strides have been made toward more accurate, efficient, and predictive OR simulations since these initial calculations: extensions to post-Hartree-Fock levels of theory,^2–5 development of OR-optimized basis sets,^6–9 improved understanding of the role of solvation,^10–14 and approaches for decomposing the OR into more chemically intuitive contributions.^15–18 Despite these efforts, OR remains a very challenging property to compute accurately, and a clear structure–property relationship remains unknown.

Many of the current difficulties related to computing OR stem from comparisons made with solution phase experiments.¹⁹ Accurately replicating experimental OR in solution is not as simple as just determining the OR for the minimum energy structure; multiple configurations of the molecule can exist in solution, requiring a conformational search and subsequent Boltzmann averaging of the OR from each conformer. Not only does this require an extensive, costly search to ensure that all relevant conformers are found, but the energies of each conformer must also be determined to high accuracy, as any errors in the exponential energy-weighting could greatly distort the conformer populations and, thus, the average OR.²⁰ Compounding these issues is the need to model the surrounding solvent, which can influence the conformer populations and contribute to the OR itself through the formation of a chiral shell around a given conformer. Implicit electrostatic models of solvation account for some of these effects but neglect specific solvent–solute interactions and the induced chirality of the solvation shell.¹⁴ Incorporating some explicit solvent molecules has mixed results, as the decision of how many solvent molecules to include and where to place them is typically made in an ad hoc fashion. Including a full solvation shell would be less arbitrary, but requires molecular dynamics (MD) simulations to determine reasonable solvation shell configurations and the levels of the theory that can be used for these system sizes may not be sufficiently accurate.^{11–13,21–23} Larger solvation shells currently preclude the use of quantum mechanical (QM) calculations altogether, and there is no general approach that provides accurate results.

Comparisons with crystalline solids, particularly molecular crystals, would eliminate or greatly simplify a number of these issues. The rigidity and regularity of a crystal would simplify the “solvent” environment by eliminating the need for conformational averaging. The extended structure around a molecule can be determined from the crystal structure and can be readily simulated through the use of periodic boundary conditions (PBC). In principle, solid state experiments allow the entire OR tensor to be resolved, unlike solution or gas phase where only the isotropic OR can be measured. However, linear birefringence is 2–3 orders of magnitude larger than circular birefringence. Therefore, OR measurements are typically performed along special axes (called optic axes) where linear birefringence is zero, although the experimental determination of the directions corresponding to these axes is difficult.^24–28 Nonetheless, a direct comparison with experimental data for solids would provide complimentary information about this property that is not easily accessible for solvated systems.

Despite early attempts to determine the optical activity of periodic systems,^29–35 no general-purpose approach is available. Rerát and Kirtman recently presented an implementation of the electric dipole-magnetic dipole polarizability with PBC,³⁶ which is one of the tensors necessary to evaluate the OR in solids. Here, we present the first derivation of the full OR tensor at density functional theory with PBC (DFT-PBC) level, including both the electric dipole-magnetic dipole and electric dipole-electric quadrupole contributions. This implementation is based on the approach of Scuseria and co-workers to compute the PBC electric dipole-electric dipole polarizability tensor using an atomic orbital (AO) basis.^37–39 This work reports numerical tests that confirm the correct physical properties of the OR tensor (i.e., its origin invariance), a preliminary analysis of the dependence of the choice of approximate density functional, basis set, and choice of gauge, and a direct comparison with experimental data for the tartaric acid chiral crystal, (2R,3R)–C₄H₆O₆.

This paper is organized as follows: In Sec. II, expressions for the electric dipole-magnetic dipole and electric dipole-electric quadrupole polarizabilities are derived that satisfy periodic boundary conditions. Section III describes the computational procedure for the simulations. In Sec. IV, simulations on cluster models and crystalline systems are presented. Section V provides concluding remarks.

Optical activity is governed by the rank-2 Buckingham/Dunn B tensor,^40–42 

where ϵ is the Levi-Cività operator, β is the electric dipole-magnetic dipole polarizability tensor, and A is the electric dipole-electric quadrupole polarizability tensor. In the length gauge,

whereas in the velocity gauge,

where ⟨⟨·; ·⟩⟩ is the response function of the enclosed operators, and Θ is in traceless form.

Throughout, we use atomic units unless otherwise specified and denote tensors over Cartesian indices using bold font and matrices over an orbital basis with an underline. Since there has been some confusion in the literature regarding the use of the G′ symbol to indicate the electric dipole-magnetic dipole polarizability with or without the ω⁻¹ factor, we use the notation β = −ω⁻¹G′.⁴¹ 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 for β and A are valid for non-resonant optical activity (ω_j ≉ ω) calculations; resonant optical activity is discussed in greater detail elsewhere.^40,43–45

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

where ℏ 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). Another common set of units for isotropic OR is deg L dm⁻¹ mol⁻¹, which can be obtained by multiplying the specific rotation in Eq. (6) by a factor of M/1000. Note that Eq. (6) depends on ω both through the constant prefactor on the right-hand side as well through B. For an oriented system, the $13Tr(B)$ factor should be replaced by β_n, the average of the components of the B tensor that are perpendicular to the direction of the beam.⁴³ By defining a chiroptical response tensor $B$⁠,⁴¹ 

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

By defining this alternative chiroptical response tensor, one can easily compute the specific rotation for light traveling in an arbitrary direction relative to the oriented material.^41,43

The evaluation of polarizability tensors using the sum-over-states formulas in Eqs. (2)–(5) is impractical, and linear response theory is used instead.^45–50 This requires the evaluation of a perturbed density as the solution of coupled-perturbed Kohn–Sham (CPKS) equations for one of the multipole moments, which is then contracted with the integrals of the paired perturbation for a specific tensor. For PBC calculations, we use the CPKS implementation for the length gauge electric dipole μ^L from Scuseria and co-workers.³⁸ Therefore, extending β and A (and consequently $B$⁠) to periodic systems reduces to the problem of defining μ^V, m, and Θ operators that are translation invariant (i.e., consistent with periodic boundary conditions) that can be evaluated over an AO basis. The evaluation of μ^V is trivial, as it amounts to a Fourier transform to k-space of the real-space ∇_r integrals over the central and replicated cells, which already have the proper translational symmetry. The only complication is the extension of the periodic CPKS algorithm to handle anti-Hermitian matrices that are stored in lower-triangular form. For the derivation of the equations for m and Θ^V, we follow the approach used for μ^L.^38,51 This is a different approach than that used in Ref. 36 for m, although the final expression is the same. Once the PBC dipole integrals are available, it is sufficient to evaluate the perturbed density for the dipole operator P^μ (with either choice of gauge) and contract it with the magnetic dipole integrals to form β and with the electric quadrupole integrals to form A.

In real space, the matrix elements of the magnetic dipole for a periodic system can be written as

where μ, ν are atomic orbitals indices, D is a matrix with (in general) three translation vectors as columns, l is a vector indexing the cells along each translation vector, with 0 the center cell, and e^−iDl⋅p is the translation operator. In real space, the matrix in Eq. (9) is pure imaginary as for molecules. Due to the presence of the translation operator, the magnetic dipole matrix is not Hermitian, i.e., $mμν0l≠mνμ0−l$⁠,

The second step is the result of re^−iDl⋅p = e^−iDl⋅pr + [r, e^−iDl⋅p], where [·, ·] is the commutator. The third step evaluates the commutator in the second term by using the relationship [r, f(p)] = if′(p), where the prime indicates the function derivative with respect to p, which holds for any function of momentum f(p) expressed in atomic units.⁵² The last step uses that the magnetic dipole and momentum matrix elements are pure imaginary, so taking their complex conjugates is equivalent to negating them.

Expressing the magnetic dipole matrix in a Hermitian form will help us to collect terms that arise during the derivation and is convenient for storing/manipulating the matrix in our implementation. Defining

produces a Hermitian matrix by splitting the origin-dependent momentum term between the given matrix element and the corresponding adjoint matrix element.

We can now proceed to transform the magnetic dipole operator to a translation invariant form. This involves substituting the length operator r with a translation invariant form,^37,38,51

Applying this transformation to the magnetic dipole operator yields

The second step uses the identity exp(ik·r)p = p exp(ik·r) + [exp(ik·r), p]. The commutator is then evaluated in the third step by using [f(r), p] = if′(r).⁵² 

Using the transformed operator, we can now derive magnetic dipole matrix elements between Bloch functions,

with

where |qk⟩ is the qth Bloch function at point k in reciprocal space, u_qk(r) is the cell periodic part of that Bloch function, and $C̲$ is the matrix of crystal orbital (CO) coefficients for that k point.^37,38,53 To simplify the notation, we write the triple sum over cell indices as a single sum over the multi-index l,

The crystal orbital PBC magnetic dipole matrix is indicated as $m̲̂$ to distinguish it from the AO matrix $m̲$ and the Hermitian AO matrix $m̲̃$⁠. We now must express these reciprocal-space, CO matrix elements in terms of real-space, AO matrix elements so they can be evaluated.

We can start to evaluate the first term on the right-hand side of Eq. (14) by applying the operators sequentially to the cell periodic ket,

where Ξ is simply a collection of the terms produced by applying the momentum operator to the ket. Continuing to apply ∇_k× yields

We can evaluate the ket of the second term of Eq. (14) in the same way,

Note that Eqs. (19d)–(19f) and (20a)–(20c) are identical except for their sign, so these contributions will cancel out. Projecting the remaining terms in Eq. (19) on the cell periodic bra, we obtain an expression for the periodic magnetic dipole matrix,

where (21a)–(21c) are the corresponding terms to (19a)–(19c). $U̲$ is the vector of matrices that transforms the CO coefficients $C̲$ to their gradient with respect to k, $∇kC̲=C̲U̲$⁠, as defined in Ref. 37. We can make one further simplification by noting that

which we can utilize to transform the real-space magnetic dipole matrix $m̲$ to the Hermitian matrix $m̲̃$ in Eq. (11). The resulting expression for the periodic magnetic dipole is

Recently, Rerát and Kirtman independently published a derivation of the PBC magnetic dipole, along with some proof of concept calculations of the β tensor.³⁶ As Rerát and Kirtman note, ensuring that β has all real elements is equivalent to the constraint that $m̲̂$ is Hermitian. However, the second term in Eq. (23) is anti-Hermitian, and the third term is non-Hermitian. To ensure that the OR tensor is real, we explicitly symmetrize the magnetic dipole matrix. In practice, we only need to symmetrize the third term of Eq. (23), as the first term is already Hermitian and the second term, being anti-Hermitian, is canceled out by the symmetrization.

The derivation of the electric quadrupole for a periodic system proceeds in much the same way as for the magnetic dipole. The following will focus on the velocity form of the electric quadrupole operator, as this is the form used in both the modified velocity gauge (MVG)⁵⁴ and the origin-independent length gauge [LG(OI)].³ We also use the non-traceless form of the operator for this derivation to simplify the equations. In the velocity gauge, the real-space matrix elements of the electric quadrupole are written as

where ⊗ is the outer product, and ⊙ is the symmetric outer product. Again, the presence of the translation operator leads to a non-Hermitian matrix, which we can correct by reallocating the deviation from Hermiticity for a given matrix element evenly between itself and its corresponding adjoint element. We first determine the form of this deviation by moving the translation operator from the ket to the bra,

Knowing the form of the deviation, we can define a Hermitian quadrupole in the same manner as for the magnetic and electric dipoles,

We can now proceed to apply the transformation of Eq. (12) to the quadrupole operator in order to make it translation invariant,

The expression for the PBC quadrupole operator in Eq. (27) is very similar in form to that of the PBC magnetic dipole in Eq. (13), with the only difference being a change in sign and the substitution × → ⊙. Consequently, applying the PBC quadrupole operator to the cell periodic ket leads to similar contributions. Applying the first term of the transformed operator to the cell periodic ket, we obtain

where Γ just collects the intermediate terms produced by applying the momentum operator. Continuing to apply ∇_k⊙ yields

where δ is the Kronecker delta. The “extra” term (29g) [compare to Eq. (19)] appears since the outer product allows ∇_k to act on k in the second term of Eq. (28).

Applying to the ket of the second term proceeds similarly,

where these terms are equal and opposite in sign to Eqs. (29d)–(29g) and, thus, cancel out of the final expression. After applying the bra, the remaining terms simplify to

It is convenient to re-express the AO basis quadrupole $Q̲$ as the Hermitian matrix $Q̲̃$⁠. Noting that

the PBC electric quadrupole can be rewritten as

to incorporate $Q̲̃$⁠. To the best of our knowledge, this is the first derivation of the PBC electric quadrupole. Note the resemblance in form with Eq. (23), where the main difference is the use of a symmetric direct product rather than a cross product, along with a change in sign in each term. As was the case with Eq. (23), we need to ensure that $Q̲̂$ is Hermitian and we again achieve this via explicit symmetrization, with the second term in Eq. (33) canceling out due to being anti-Hermitian.

Once the perturbed density P^μ is evaluated from the periodic CPKS equations, the β and A tensors are computed as

and then combined as in Eqs. (1) and (7) to obtain the $B$ tensor. Note that both expressions in Eqs. (34) and (35) are missing a term that depends on the k-derivative of the perturbed density.⁵⁵ Based on the results in Sec. IV, this contribution seems numerically small. Although the tensors calculated within the MVG approach are origin invariant,⁵⁴ this is not directly the case with the LG approach. However, this issue is resolved by applying the LG(OI) transformation,^2,3 which avoids the need to employ gauge-including atomic orbitals.^48,56–58 This approach is based on the singular value decomposition (SVD) of the mixed-gauge electric dipole-electric dipole polarizability and the application of the inverse transformation to the β and A tensors.^2,3 The extra tensor can be readily computed by combining P^μ with the μ^V integrals and the SVD of a square matrix of dimension 3 is quick. Therefore, the cost of the LG(OI) transformation is negligible, and this approach is computationally more efficient than MVG because the evaluation of the tensors in the static limit is not required.

The optical rotation of chiral crystals can only be measured accurately along the optic axis because, in any other direction, the circular birefringence is obscured by the much larger linear birefringence.^24,59 The optic axis is defined as the direction perpendicular to the central circular section(s) of the ellipsoid called the indicatrix. The semi-axes of this ellipsoid are oriented along the principal axes of the dielectric tensor ɛ_r. For low-polar media, the dielectric tensor is related to the electric dipole-electric dipole polarizability tensor α (all in a.u.) as^24,59

where 1₃ is the 3 × 3 unit matrix, 4π is the dielectric permittivity of vacuum in a.u., and V is the volume of the unit cell also in a.u. Therefore, the principal axes of ɛ_r are the same as those for α. The length of the indicatrix semi-axes is equal to the refractive indices along the ɛ_r principal axes: $nα=(εr)αα$⁠, where the Greek letter indicates one of the principal axes. The common convention is to organize n_α in increasing order. By definition, the refractive indices in the directions perpendicular to the optic axis are equal to each other, while the value along the optic axis is different. Uniaxial crystals only have one optic axis and it coincides with one of the indicatrix semi-axes. The refractive index along the optic axis is called the extraordinary refractive index, n_e, while the two equal values in the directions perpendicular to the optic axis are called the ordinary refractive indices, n_o.

Biaxial crystals have two optic axes that lay on the plane defined by the x and z ɛ_r principal axes.^24,59 The z principal axis bisects the optic axes, and the angle 2θ between the principal axes can be found as

Therefore, a comparison with experimental data for a biaxial crystal requires two sequential transformations of the $B$ tensor. The first one reorients the tensor along the principal axes of ɛ_r,

where U_α is the matrix of the eigenvectors of the electric dipole-electric dipole polarizability tensor α, ordered in increasing order of the corresponding eigenvalues. The second transformation is a rotation around the ɛ_ry principal axis by the angle θ, so that the resulting z_θ axis is aligned with the optic axis,

where R_y(θ) is the rotation matrix. Finally, the z diagonal component of the transformed tensor, $Bzzoa$⁠, corresponds to the measured OR of the crystal along the optic axis. The OR for chiral crystals is usually expressed in deg/mm; the conversion from dm⁻¹ mol⁻¹ to deg/mm is obtained by multiplying the OR by $(100NAVdm)−1$⁠, when N_A is Avogadro’s number and V_dm is the unit cell volume in dm³. We emphasize two important details: (i) in approximate (i.e., virtually all) calculations, the optic axes are different for different choices of gauge because the numerical values of the α tensor depend on this choice; (ii) for LG(OI), the $B$ tensor is already oriented along the principal axis of the corresponding α tensor;^2,3 therefore, only the second transformation [in Eq. (39)] is necessary.

All calculations were performed with a development version of the GAUSSIAN suite of programs.⁶⁰ Two functionals were used: the pure functional SVWN5 (called LDA in Ref. 36) and the range-separated hybrid HSE06,⁶¹ which has exact exchange only in the short range. Both functionals allow faster simulations for periodic systems because the exact exchange is absent in the former functional and short-range in the latter one. A pruned (99 590) grid was used for both the energy and linear response part of the calculations for all systems. Several basis sets and k-point grids were used for different test systems, as detailed in Sec. IV. The OR tensors were computed with the MVG⁶² and LG(OI)^2,3 schemes.

The first system we consider is a chain of H₂O₂ molecules with a 1D periodicity along the direction of the O–O bond, with a single molecule in the unit cell and a lattice vector a = 4 Å, as in Ref. 36. The geometry of the unit cell is reported in Table S1 of the supplementary material. This is almost the same setup as in Ref. 36: the lattice vector is along the x axis, but now the origin is located at the center of the O–O bond, so that the C₂ rotation axis coincides with the y axis. A grid of 500 k points was used.

The results in Tables I and II demonstrate that the full OR tensor $B$ is origin invariant as for molecules. The results in the tables are obtained at the SVWN5/STO-3G level with two choices of gauge, MVG and LG(OI), respectively, at the original geometry and one translated by −100 Å in every Cartesian direction. For simplicity, we used the guess density built from the unperturbed COs (i.e., the SOS approach in Ref. 36) rather than the full relaxed density to avoid numerical noise from the solution of the coupled perturbed equations. The tables also report the β and A tensor contributions to $B$⁠, see Eqs. (2)–(5), which are individually origin dependent, but in opposite directions, such that these effects cancel out in $B$⁠. More specifically, for MVG in Table I, only the trace of β is origin invariant, but the individual tensor elements are not; the yy element is origin invariant because of the symmetry of this specific system, but it would not be in general. For LG(OI) in Table II, the diagonal elements of the β and A tensors are fixed by the LG(OI) transformation, and only the off-diagonal elements are origin-dependent. It is important to note that the elements of the β and A tensors are comparable in magnitude. Furthermore, each diagonal element of $B$ depends on the other two diagonal elements of β and A, see Eq. (7): for instance, $Bzz∝βxx+βyy$⁠. These details are not relevant for the isotropic OR of systems in gas and solution phase because this is proportional to $Tr(B)$⁠, and the A tensor is traceless, but they are fundamental for oriented systems. Therefore, only the $B$ tensor can be compared with experimental data.

The same model system is used to compare PBC and large molecular cluster calculations, which helps to verify that the implementation is correct. The SVWN5/cc-pVQZ model chemistry was employed as in Ref. 36, and the PBC calculation used a 1000 k-point mesh. We consider two molecular cluster sizes with 39 and 37 H₂O₂ units. The “unit cell” tensors are then obtained as the average of the 39-unit calculation and as the averaged difference between the 39- and 37-unit calculations. For this comparison, using a geometry with the origin in the center of the O–O bond is convenient because the unit cell in the PBC calculation matches the central unit in the cluster calculations. The results for the fully relaxed density with MVG and LG(OI) are shown in Figs. 1 and 2, respectively. The figures include the tensor elements and the spatially averaged trace (all in a.u.) for the full $B$ tensor as well as for the β and A contributions. The isotropic OR values (in units of dm⁻¹ mol⁻¹) are reported in Table III, and all the tensor elements are reported in Tables S2–S4 of the supplementary material.

The data in Figs. 1 and 2 show very good agreement between the molecular cluster and PBC calculations for all the elements of the $B$ tensor, irrespective of the choice of gauge. The largest difference is found for the xz element: 9%–11% for MVG and 15%–17% for LG(OI). The discrepancy is likely due simultaneously to the limited size of the cluster models and the lack of the ∇_kP^μ terms in the PBC expressions in Eqs. (34) and (35), respectively.⁵⁵ In this case, these effects are small but they need to be investigated further. The MVG and LG(OI) results also match each other for the isotropic OR (see the trace values in the figures and the results in Table III). This indicates that this basis set is close to completeness for this system and the calculations approach gauge invariance.⁶³ The PBC-MVG value in Table III is in very good agreement with the value reported in Ref. 36: 36.1 dm⁻¹ mol⁻¹ (the small discrepancy is due to the extreme sensitivity of this property to unavoidable differences in numerical thresholds between different implementations). The differences in the individual elements of the $B$ tensor between MVG and LG(OI) are due to the fact the SVD transformation in the latter corresponds to a rotation of the system;³ repeating the MVG calculation at this rotated geometry would lead to matching elements. More importantly, the molecular cluster and PBC calculations do not match each other for the β and A contributions to $B$⁠, as shown in panels (b) and (c) of Figs. 1 and 2. The yy and zz elements of the β contribution are much smaller in the PBC case compared to the molecular cluster case; simultaneously, the corresponding elements of the A contribution have the opposite sign. It is only in the full $B$ tensor that these contributions combine properly to provide the same picture.

We provide a preliminary comparison of the basis set and functional dependence of OR calculations between isolated molecules and periodic systems. Specifically, we compare the isolated H₂O₂ molecule at the geometry of the unit cell with 1D-periodic calculations as in the previous examples for two functionals: SVWN5 and HSE06, and three basis sets: aug-cc-pVTZ, cc-pVQZ, and aug-cc-pVQZ (in the following shortened to aTZ, QZ, and aQZ, respectively). A grid of 500 k points was used in the PBC calculations. To simplify the comparison, we focus on the isotropic OR [i.e., Tr$(B)$] rather than the whole tensor. The specific rotation values for all calculations are reported in Table S5 of the supplementary material. Here, we report relative % changes in OR calculated as

where X = MVG, LG(OI); Y = aTZ, QZ, aQZ; W = SVWN5, HSE06; and Z = Mol, PBC. In other words, Δ_QZ in Eq. (40) tests the change in basis set and choice of gauge within the same method (W) and system (Z), while Δ_HSE06 in Eq. (41) tests the change in functional within the same choice of gauge (X), basis set (Y), and system (Z).

The calculated values of Δ_QZ are plotted in Figs. 3(a) and 3(b). These results show that, for the molecular case, both functionals are close to basis set completeness and gauge invariance. However, the differences between aTZ and aQZ for any choice of the gauge are considerably smaller than with QZ. For instance, the aTZ-aQZ difference within the same gauge is 2%–3% and $<10$% between choices of gauge; on the other hand, the difference with either basis set and QZ is 15%–25% within the same gauge, and the MVG-LG(OI) difference for QZ is 16%–21%. This indicates that the use of diffuse functions is more important than increasing the basis set size from Tζ to Qζ. However, the situation is very different for the PBC calculations. Figures 3(a) and 3(b) show that QZ results are 2–2.5 times smaller than the aTZ ones, and 3–4 times smaller than the aQZ ones. The Δ_QZ values are smaller for HSE06/aTZ than for SVWN5, but the opposite trend manifests for aQZ. The difference between choices of the gauge is most significant for the aTZ basis set. The change in trend with respect to the molecular case is likely due to the spatial extension of diffuse functions, which may lead to strong interactions between repeated units in adjacent cells.

The difference between functionals (Δ_HSE06) is shown in Fig. 3(c). The SVWN5 results are significantly larger than those obtained with HSE06. For the molecular case, QZ values are twice as large with the pure functional than with the hybrid, and thrice as large with aTZ and aQZ. For the periodic calculations, this difference increases to a factor of three and four for QZ and aTZ, respectively, but it decreases to a factor of two for aQZ. The change from molecular to periodic cases is likely due to intermolecular interactions. However, the different behavior across basis sets may be another indication that diffuse functions may be problematic for PBC calculations. In any case, the choice of functional seems to have the largest impact on the calculated values, followed by the choice of basis set, and to a much smaller degree by the choice of gauge.

The last system that we consider is tartaric acid, (2R,3R)–C₄H₆O₆, for which experimental data are available.²⁵ We used the same geometry and basis set [6-31G(d)] as in Rérat and Kirtman,³⁶ but the calculations were performed with HSE06 rather than PBE0⁶⁴ because the former is computationally more efficient. The k point grid is 6 × 8 × 8. The geometry was also rotated so that the coordinate axes coincide with the principal axes of the symmetric part of the mixed-gauge electric dipole-electric dipole polarizability, which is the best representation of the coordinate system for the LG(OI) approach as explained in Ref. 3. The unit cell geometry is reported in Table S6 of the supplementary material.

Table IV reports a comparison of experimental and calculated data for the angle 2θ, the birefringence δ_αβ = n_α − n_β, and the optical rotation along the z axis, all of which are defined in Sec. II C. We report these values calculated along the ɛ_r principal axes (denoted by the superscript ɛ) and with the coordinate system rotated to have the z axis along the optic axis (denoted by the oa superscript), see discussion in Sec. II C. The full tensors and other intermediate quantities that were used to evaluate the data in the table are reported in Tables S7–S10 of the supplementary material. The agreement between theoretical and experimental values in Table IV is very good, considering the basis set utilized in the calculations. The error is 14%–16% for the angle 2θ and 19%–29% for the birefringence δ^ɛ, indicating that the electric dipole-electric dipole polarizability is reasonably accurate. The error for the OR is larger (about 50% if the latest experimental value of −12.3 ± 1.0 deg/mm is considered), but it is consistent with the error for the isotropic OR of molecules at this level of theory.^58,63 The choice of gauge is less important than that of the level of theory as the MVG and LG(OI) results are very close to each other. This is also confirmed by the polar plots representation of the OR tensors in Figs. 4 and 5. In these plots, the unit cell was rotated such that the coordinate system matches the orientation of the tensor, and z coincides with the optic axis. The difference in geometrical orientation between MVG and LG(OI) is small, and barely detectable in the figures. The polar plots are also virtually indistinguishable between figures. It is interesting to emphasize that the plots of the β and A contributions to $B$ were scaled down by a factor of 10 compared to the full $B$ plot, to fit within the same picture frame.